Non-monotonic MDCEV forecasting

Michel Bierlaire, EPFL Fri Jul 25 2025, 17:27:35

Forecasting with a MDCEV model and the “non-monotonic utility” specification.

Example: non monotonic utility
Forecasting observation 0 / 2 [10 draws]
============ Comparison ===================
Brute force: {1: '221', 2: '111', 3: '158', 4: '10'} objective 196, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '221', 2: '111', 3: '158', 4: '10'} objective 196, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '27', 2: '402', 3: '67.1', 4: '4.33'} objective 191, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '27', 2: '402', 3: '67.1', 4: '4.33'} objective 191, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '31.2', 2: '322', 3: '142', 4: '4.99'} objective 180, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '31.2', 2: '322', 3: '142', 4: '4.99'} objective 180, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '2.5', 2: '482', 3: '12', 4: '3.13'} objective 279, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '2.5', 2: '482', 3: '12', 4: '3.13'} objective 279, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '30.1', 2: '413', 3: '50.4', 4: '6.12'} objective 210, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '30.1', 2: '413', 3: '50.4', 4: '6.11'} objective 210, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '19.6', 2: '439', 3: '34.2', 4: '6.72'} objective 223, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '19.6', 2: '439', 3: '34.2', 4: '6.72'} objective 223, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '6.09', 2: '463', 3: '26.9', 4: '4.44'} objective 309, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '6.09', 2: '463', 3: '26.9', 4: '4.44'} objective 309, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '87.1', 2: '163', 3: '239', 4: '10.8'} objective 182, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '87.1', 2: '163', 3: '239', 4: '10.8'} objective 182, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '18.6', 2: '397', 3: '78.2', 4: '5.94'} objective 241, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '18.6', 2: '397', 3: '78.2', 4: '5.94'} objective 241, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '13.9', 2: '409', 3: '72.3', 4: '4.49'} objective 217, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '13.9', 2: '409', 3: '72.3', 4: '4.49'} objective 217, constraint 500, choice set {1, 2, 3, 4}
Forecasting observation 1 / 2 [10 draws]
============ Comparison ===================
Brute force: {1: '19', 2: '400', 3: '74.5', 4: '6.43'} objective 226, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '19', 2: '400', 3: '74.5', 4: '6.42'} objective 226, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '31.8', 2: '242', 3: '220', 4: '6.17'} objective 173, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '31.8', 2: '242', 3: '220', 4: '6.17'} objective 173, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '6.71', 2: '56.2', 3: '433', 4: '4.45'} objective 261, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '6.71', 2: '56.1', 3: '433', 4: '4.46'} objective 261, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '149', 2: '322', 3: '25.4', 4: '4.23'} objective 230, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '149', 2: '322', 3: '25.4', 4: '4.23'} objective 230, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '4.01', 2: '421', 3: '71.2', 4: '3.68'} objective 225, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '4.01', 2: '421', 3: '71.2', 4: '3.68'} objective 225, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '9.85', 2: '76.2', 3: '411', 4: '3.41'} objective 260, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '9.85', 2: '76.2', 3: '411', 4: '3.41'} objective 260, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '12.9', 2: '416', 3: '62.1', 4: '8.71'} objective 230, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '12.9', 2: '416', 3: '62.1', 4: '8.71'} objective 230, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '10.3', 2: '307', 3: '167', 4: '15.8'} objective 198, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '10.3', 2: '307', 3: '167', 4: '15.8'} objective 198, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '51.3', 2: '425', 3: '16.2', 4: '7.79'} objective 248, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '51.3', 2: '425', 3: '16.2', 4: '7.78'} objective 248, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '14.7', 2: '255', 3: '224', 4: '6.09'} objective 173, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '14.7', 2: '255', 3: '224', 4: '6.1'} objective 173, constraint 500, choice set {1, 2, 3, 4}
Forecasting observation 0 / 2 [2000 draws]
Forecasting observation 1 / 2 [2000 draws]
Execution time for 2000 draws with brute force algorithm: 186 seconds
Forecasting observation 0 / 2 [2000 draws]
Forecasting observation 1 / 2 [2000 draws]
Execution time for 2000 draws with analytical algorithm: 4.07 seconds
                 1            2            3            4
count  2000.000000  2000.000000  2000.000000  2000.000000
mean     34.980702   328.944391   124.743144    11.331762
std      68.995633   139.551899   126.916481    31.995071
min       0.000000     0.000000     0.000000     0.000000
25%       6.475080   229.937979    30.649985     4.185199
50%      13.010862   375.453036    75.777429     5.980568
75%      26.390613   444.594704   177.272641     9.264476
max     497.565724   500.000000   500.000000   486.691690
                 1            2            3            4
count  2000.000000  2000.000000  2000.000000  2000.000000
mean     34.981379   328.948512   124.741180    11.328929
std      68.994759   139.551276   126.917956    31.978639
min       0.000000     0.000000     0.000000     0.000000
25%       6.477249   230.229019    30.649298     4.185054
50%      13.010914   375.448252    75.783575     5.979084
75%      26.389402   444.592920   177.240335     9.260488
max     497.564694   500.000000   500.000000   486.694703
                 1            2            3            4
count  2000.000000  2000.000000  2000.000000  2000.000000
mean     36.583961   320.455404   131.542141    11.418493
std      69.949585   141.228837   130.806341    32.359974
min       0.000000     0.000000     0.000000     0.000000
25%       6.976156   216.947521    33.133503     4.009581
50%      13.750494   366.381674    79.058918     5.887904
75%      31.763022   437.153893   194.313891     8.912869
max     499.580733   500.000000   499.403643   420.092523
                 1            2            3            4
count  2000.000000  2000.000000  2000.000000  2000.000000
mean     36.588174   320.454972   131.537215    11.419639
std      69.956395   141.232016   130.804352    32.373628
min       0.000000     0.000000     0.000000     0.000000
25%       6.979068   217.084955    33.136305     4.008129
50%      13.747908   366.385056    79.055459     5.887673
75%      31.760013   437.151264   194.311629     8.906076
max     499.580612   500.000000   499.403650   420.000183

import sys
import time

import numpy as np
import pandas as pd
from IPython.display import display

import biogeme.biogeme_logging as blog
from biogeme.database import Database
from biogeme.results_processing import EstimationResults
from non_monotonic_specification import the_non_monotonic
from process_data import database

logger = blog.get_screen_logger(level=blog.INFO)
logger.info('Example: non monotonic utility')

result_file = 'saved_results/non_monotonic.yaml'
try:
    results = EstimationResults.from_yaml_file(filename=result_file)
except FileNotFoundError as e:
    print(e)
    print(f'File {result_file} is missing.')
    sys.exit()

the_non_monotonic.estimation_results = results

# %
# We apply the model only on the first two rows of the database.
two_rows_of_database: Database = database.extract_rows([10, 11])
# %
budget_in_hours = 500

# %
# # Validation

# %
# As the implementation is still experimental, we compare the result obtained by the bruteforce algorithm and
# the analytical algorithm for a few draws.

# Note that minor discrepancies between the outcome of the two algorithms are likely to occur, due to numerical
# imprecision, inevitable in finite arithmetic.

# However, if there are major differences, it should be reported.

# %
number_of_draws = 10

# %
# We generate the draws
epsilons = [
    np.random.gumbel(
        loc=0, scale=1, size=(number_of_draws, the_non_monotonic.number_of_alternatives)
    )
    for _ in range(two_rows_of_database.num_rows())
]

# %
# We first compare the results obtained from the brute force and the analytical algorithms, for each draw.
the_non_monotonic.validate_forecast(
    database=two_rows_of_database, total_budget=budget_in_hours, epsilons=epsilons
)

# %
# # Forecasting
# We use a larger number of draws to obtain the forecast.

# %
number_of_draws = 2000

# %
# We generate the draws
epsilons = the_non_monotonic.generate_epsilons(
    number_of_observations=two_rows_of_database.num_rows(),
    number_of_draws=number_of_draws,
)
# %
# First, the brute force algorithm.
start_time = time.time()
optimal_consumptions_brute_force: list[pd.DataFrame] = the_non_monotonic.forecast(
    database=two_rows_of_database,
    total_budget=budget_in_hours,
    epsilons=epsilons,
    brute_force=True,
)
end_time = time.time()

# %
print(
    f'Execution time for {number_of_draws} draws with brute force algorithm: {end_time-start_time:.3g} seconds'
)

# %
# Then, the analytical algorithm.
start_time = time.time()
optimal_consumptions_analytical: list[pd.DataFrame] = the_non_monotonic.forecast(
    database=two_rows_of_database,
    total_budget=budget_in_hours,
    epsilons=epsilons,
    brute_force=False,
)
end_time = time.time()

# %
print(
    f'Execution time for {number_of_draws} draws with analytical algorithm: {end_time-start_time:.3g} seconds'
)

# %
# Results for the first observation, brute force method
display(optimal_consumptions_brute_force[0].describe())

# %
# Results for the first observation, analytical method
display(optimal_consumptions_analytical[0].describe())

# %
# Results for the second observation, brute force method
display(optimal_consumptions_brute_force[1].describe())

# %
# Results for the second observation, analytical method
display(optimal_consumptions_analytical[1].describe())