The McDonald's Diet Problem A Case Study in Optimization Using AMPL

We now scale up to the whole McDonald's menu for which we have usable nutritional information -- 63 items in all. We return to the food an nutrition limits that we used in our first small-scale run, with the addition of recommended daily limits for far, saturated fat, cholesterol, and sodium. (It's also recommended that "calories from fat" be at most 30% of total calories, but that can't be fit so easily into the current model -- we'll return to it shortly.) Since we're providing many more possibilities than before, we expect that the optimization will now find a more economical combination of foods. Indeed the resulting diet is far cheaper than any we've seen previously -- only \$5.36, compared to upwards of \$14. But this diet fails the "laugh test" . . . you aren't going to consume 16 containers of honey and 4 of sweet-and-sour sauce, unless you also buy a lot of "Chicken McNuggets". Here the expansion of the data has brought up some more constraints that weren't anticipated in dealing with the smaller-scale data.

```   stretto% ampl

ampl: model diet1.mod;
ampl: data diet2.dat;

ampl data: solve;
MINOS 5.4: ignoring integrality of 63 variables
MINOS 5.4: optimal solution found.
15 iterations, objective 5.363036219
Objective = Total_Cost

ampl: option omit_zero_rows 1;

Cheeseburger   2.05861
"Sweet 'N Sour Sauce"   4.12272
Honey  16.1668
Cheerios   2.26957
'1% Lowfat Milk'   1.77758
'Orange Juice'   0.408178
;

ampl: display n_min,Diet.body,n_max;

:        n_min   Diet.body    n_max      :=
Cal       2000   2017.93     Infinity
CalFat       0    306.548    Infinity
Fat          0     32.7316        100
SatFat       0     12.9956         30
Chol         0    123.605         375
Sodium       0   3000            3000
Carbo      350    375             375
Protein     55     55        Infinity
VitA       100    100        Infinity
VitC       100    100        Infinity
Calcium    100    100        Infinity
Iron       100    100        Infinity
;
```

Before we fix up the model, we'll return to our previous tricks for making the solution more palatable: setting an upper limit of two of any food, and requiring an integer solution. We change the upper limits, switch the solver to "cplex", type solve, and ... and ... and then ... and then many minutes later, we find that the solver has reached some kind of "node limit" and has returned with a solution that doesn't seem to be necessarily optimal. It costs \$9.62, much more than the previous fractional solution, and is stil fairly ridiculous -- too many drinks, and inappropriate condiments (although you might try putting the hot mustard sauce on your hamburgers).

After conferring with an expert, we find a command for getting around the mysterious node limit. Then, after an even longer wait, the solver produces an optimal solution with a cost of only \$7.90, though with not much greater palatability.

This experience is not a fluke. It's generally much harder to solve for integer values than for continuous ones. And it takes more specialized technical knowledge to get the best performance out of integer solvers. You can quite easily formulate an integer program that is too big for your computer to solve, whereas its fairly hard to come up with a useful linear program that's too big for the computer you're likely to buy today.

```   ampl: let {j in FOOD} f_max[j] := 2;

ampl: option solver cplex;

ampl: solve;
CPLEX 2.1: node limit with integer solution.
Current node limit = 20000; objective 9.62
56219 simplex iterations
20000 branch-and-bound nodes
Objective = Total_Cost

Hamburger  2
'Hot Mustard Sauce'  2
"Sweet 'N Sour Sauce"  1
'Bacon Bits'  1
Cheerios  2
'Chocolate Shake'  2
'Strawberry Shake'  1
'Orange Juice'  1
'H-C Orange Drink (small)'  1
;

ampl: display n_min,Diet.body,n_max;

:        n_min Diet.body    n_max      :=
Cal       2000    2195     Infinity
CalFat       0     415     Infinity
Fat          0      46          100
SatFat       0      19           30
Chol         0     186          375
Sodium       0    2925         3000
Carbo      350     371          375
Protein     55      72     Infinity
VitA       100     138     Infinity
VitC       100     190     Infinity
Calcium    100     137     Infinity
Iron       100     110     Infinity
;

ampl: option cplex_options 'nodes 200000';

ampl: solve;
CPLEX 2.1: nodes 200000
CPLEX 2.1: optimal integer solution; objective 7.9
66614 simplex iterations
23906 branch-and-bound nodes
Objective = Total_Cost

Hamburger  2
"Sweet 'N Sour Sauce"  2
Honey  2
Cheerios  1
Wheaties  1
'Cinnamon Raisin Danish'  2
'Raspberry Danish'  1
'1% Lowfat Milk'  2
'Orange Juice'  1
;

ampl: display n_min,Diet.body,n_max;

:        n_min Diet.body    n_max      :=
Cal       2000    2440     Infinity
CalFat       0     765     Infinity
Fat          0      84          100
SatFat       0      28           30
Chol         0     235          375
Sodium       0    2920         3000
Carbo      350     357          375
Protein     55      63     Infinity
VitA       100     101     Infinity
VitC       100     171     Infinity
Calcium    100     110     Infinity
Iron       100     102     Infinity
;
```

Switching to a different solver (named "osl") that is available, we find that we get much more efficient performance. We also show here a way to get some output on the progress of the solver. It records a best integer solution so far ("Best") and a lower bound on the best possible ("Bound"); the difference between these ("Gap") must be forced to zero before the current best solution can be declared optimal.

```   ampl: option solver osl;
ampl: option osl_options 'bbdisplay 2 bbdispfreq 100 timing 1';

ampl: solve;
OSL 1.2:
bbdisplay 2
bbdispfreq 100
timing 1

The primal algorithm has been chosen
Switching to Devex pricing--Success rate

Nodes  Integer  Simplex
Searched   Soln's    Iters  Best            Bound           Gap
21        1       85  10.43           6.8560894       3.6e+00
100        1      335  10.43           7.0060121       3.4e+00
186        2      618  8.55            7.1694099       1.4e+00
200        2      671  8.55            7.2504978       1.3e+00
300        2     1004  8.55            7.2936297       1.3e+00
400        2     1247  8.55            7.3292927       1.2e+00
472        3     1350  7.9             7.3502945       5.5e-01
500        3     1408  7.9             7.3584184       5.4e-01
600        3     1678  7.9             7.4142911       4.9e-01
700        3     1965  7.9             7.5014493       4.0e-01
800        3     2183  7.9             7.5382228       3.6e-01
900        3     2443  7.9             7.5772978       3.2e-01
1000        3     2668  7.9             7.6049684       3.0e-01
1100        3     2880  7.9             7.6204369       2.8e-01
1200        3     3103  7.9             7.6628878       2.4e-01
1300        3     3299  7.9             7.6878921       2.1e-01
1400        3     3477  7.9             7.7403236       1.6e-01
1500        3     3630  7.9             7.7710816       1.3e-01
1600        3     3785  7.9             7.802429        9.8e-02
1700        3     3922  7.9             7.8329599       6.7e-02
1800        3     4066  7.9             7.8658743       3.4e-02
1900        3     4168  7.9             7.8967828       3.2e-03

Times (seconds):
Input =  0.0666667
Solve =  61.4
Output = 0.0166667

OSL 1.2: optimal solution; objective 7.9
4190 simplex iterations; 1917 branch-and-bound nodes
3 feasible integer solutions found
Objective = Total_Cost
```

But wait . . . A newer version of cplex is now available. It exhibits much more competitive performance, in terms of total solve time. The competitive situation is actually even more complicated, because integer solvers have many settings that can be adjusted to improve performance. There are supplements to the AMPL book that explain these settings for both the cplex and osl solvers. Clearly it's not a simple matter to say which solver is better in general.

```   ampl: option solver cplex3;
ampl: option cplex_options 'mipdisplay 2 mipinterval 100 timing 1';

ampl: solve;
CPLEX 3.0: mipdisplay 2
mipinterval 100
timing 1

MIP Presolve modified 3 coefficients.
Reduced MIP has 11 rows, 64 columns, and 490 nonzeros.

Nodes                                         Cuts/
Node  Left     Objective  IInf  Best Integer     Best Node   ItCnt

0     0        6.8489     6                      6.8489      42
*    48    19        9.2400     0        9.2400        6.8854     126
*    60    25        8.9200     0        8.9200        6.9015     150
*    75    30        8.8900     0        8.8900        6.9021     176
*    98    37        8.4900     0        8.4900        6.9073     224
100    37        6.9246     7        8.4900        6.9246     227
200    98        cutoff              8.4900        7.0537     494
300   149        8.4845     3        8.4900        7.1039     717
400   196        7.4372     5        8.4900        7.1533     964
500   245        7.4666     7        8.4900        7.1760    1155
600   292        7.5399     8        8.4900        7.2124    1389
700   349        8.3061     4        8.4900        7.2525    1606
800   403    infeasible              8.4900        7.2659    1864
900   455        cutoff              8.4900        7.2856    2101
1000   502        cutoff              8.4900        7.3024    2326
*  1046   417        8.0500     0        8.0500        7.3117    2410
1100   430        7.3250     6        8.0500        7.3250    2535
1200   465        7.9429     3        8.0500        7.3368    2758
1300   506        cutoff              8.0500        7.3507    2939
1400   542        8.0304     4        8.0500        7.3692    3140
1500   568        7.7988     4        8.0500        7.3840    3329
1600   593        7.9564     4        8.0500        7.4055    3493
1700   621        7.4184     5        8.0500        7.4184    3695
1800   657        7.5516     5        8.0500        7.4341    3902
1900   683        7.7420     4        8.0500        7.4427    4099
2000   701        cutoff              8.0500        7.4547    4299
2100   730        8.0493     3        8.0500        7.4693    4524
2200   750        7.4814     6        8.0500        7.4787    4696
2300   778        7.5702     4        8.0500        7.4951    4892
2400   800        7.5581     6        8.0500        7.5107    5054
2500   819        8.0398     4        8.0500        7.5205    5224
2600   841        8.0331     5        8.0500        7.5266    5411
2700   860        cutoff              8.0500        7.5392    5577
2800   873        7.5498     5        8.0500        7.5465    5744
2900   881        7.5577     6        8.0500        7.5574    5938
3000   909        cutoff              8.0500        7.5667    6143
3100   932        7.9943     4        8.0500        7.5788    6338
3200   946        7.6392     6        8.0500        7.5869    6503
3300   967        cutoff              8.0500        7.5953    6704
3400   985    infeasible              8.0500        7.6036    6873
3500   996        7.6647     5        8.0500        7.6182    7067
3600  1004        7.7547     5        8.0500        7.6255    7255
3700  1032        cutoff              8.0500        7.6346    7423
3800  1046        8.0264     2        8.0500        7.6417    7586
3900  1055        cutoff              8.0500        7.6497    7765
4000  1060        7.7646     5        8.0500        7.6558    7920
4100  1075        cutoff              8.0500        7.6636    8111
4200  1092        7.6711     5        8.0500        7.6711    8283
4300  1101        8.0181     4        8.0500        7.6777    8461
4400  1111        cutoff              8.0500        7.6844    8653
4500  1108        7.7624     5        8.0500        7.6921    8811
4600  1110        7.7006     5        8.0500        7.6988    8976
4700  1114        cutoff              8.0500        7.7075    9144
4800  1116        7.7153     6        8.0500        7.7142    9329
4900  1127        7.7196     5        8.0500        7.7196    9501
5000  1127        7.7786     6        8.0500        7.7285    9676
5100  1127        7.8828     6        8.0500        7.7337    9848
5200  1134        8.0072     7        8.0500        7.7394    9990
5300  1133        cutoff              8.0500        7.7465   10162
5400  1136        7.7547     5        8.0500        7.7518   10335
5500  1136        7.8516     5        8.0500        7.7584   10478
5600  1144        7.7627     4        8.0500        7.7627   10629
5700  1129        7.8148     3        8.0500        7.7717   10791
5800  1135        cutoff              8.0500        7.7776   10935
5900  1129        7.8555     3        8.0500        7.7826   11104
6000  1128        7.9558     3        8.0500        7.7865   11254
6100  1132        cutoff              8.0500        7.7942   11428
6200  1119        7.7981     5        8.0500        7.7981   11599
6300  1123        7.9140     4        8.0500        7.8027   11778
6400  1122        7.8069     5        8.0500        7.8069   11948
6500  1116        cutoff              8.0500        7.8118   12097
6600  1108        7.9152     3        8.0500        7.8158   12259
6700  1095        cutoff              8.0500        7.8224   12431
6800  1081        cutoff              8.0500        7.8271   12588
6900  1068        7.8368     6        8.0500        7.8339   12727
7000  1051    infeasible              8.0500        7.8381   12857
7100  1052        cutoff              8.0500        7.8427   12991
7200  1039        cutoff              8.0500        7.8501   13142
7300  1019        7.9690     4        8.0500        7.8562   13286
7400   999        7.8623     5        8.0500        7.8623   13416
7500   979        8.0289     5        8.0500        7.8680   13562
7600   948        cutoff              8.0500        7.8760   13709
7700   918        8.0000     4        8.0500        7.8816   13840
7800   894        cutoff              8.0500        7.8882   14008
7900   861        7.9014     3        8.0500        7.8949   14154
8000   833    infeasible              8.0500        7.9000   14282
*  8021     0        7.9000     0        7.9000                 14306

Times (seconds):
Input =  0.0166667
Solve =  54.8833
Output = 0

CPLEX 3.0: optimal integer solution; objective 7.9
14306 simplex iterations
8021 branch-and-bound nodes
Objective = Total_Cost
```