An AMPL model for the diet problem
1. Problem Description
Consider the problem of choosing prepared foods to meet certain nutritional requirements. Suppose that precooked dinners of the following kinds are available for the following prices per package:
BEEF |
Beef |
$3.19 |
CHK |
Chicken |
$2.59 |
FISH |
Fish |
$2.29 |
HAM |
Ham |
$2.89 |
MCH |
Macaroni & Cheese |
$1.89 |
MTL |
Meat Loaf |
$1.99 |
SPG |
Spaghetti |
$1.99 |
TUR |
Turkey |
$2.49 |
These dinners provide the following percentages, per package, of the minimum daily requirements for vitamins A, C, B1, and B2.
A |
C |
B1 |
B2 |
|
BEEF |
60% |
20% |
10% |
15% |
CHK |
8 |
0 |
20 |
20 |
FISH |
8 |
10 |
15 |
10 |
HAM |
40 |
40 |
35 |
10 |
MCH |
15 |
35 |
15 |
15 |
MTL |
70 |
30 |
15 |
15 |
SPG |
25 |
50 |
25 |
15 |
TUR |
60 |
20 |
15 |
10 |
The problem is to find the cheapest combination of packages that will meet a week’s requirements—that is, at least 700% of the daily requirement for each nutrient.
2. LP Formulation
Let us write XBEEF for the number of packages of beef dinner to be purchased, XCHK for the number of packages of beef dinner to be purchased, and so forth. Then the total cost of the diet will be:
Total cost =
3.19 XBEEF + 2.59 XCHK + 2.29 XFISH +2.89 XHAM +
1.89 XMCH + 1.99 XMTL + 1.99 XSPG + 2.49 XTUR
The total percentage of the vitamin A requirement is given by a similar formula, except that XBEEF , XCHK , and so forth are multiplied by the percentage per package instead of the cost per package:
Total percentage of vitamin A daily requirement met =
60 XBEEF + 8 XCHK + 8 XFISH + 40 XHAM +
15 XMCH + 70 XMTL + 25 XSPG + 60 XTUR
This amount needs to be greater than or equal to 700 percent. There is a similar formula for each of the other vitamins, and each of these also needs to be greater than or equal to 700 percent.
Putting these all together, we have the following linear program:
Minimize
3.19 XBEEF + 2.59 XCHK + 2.29 XFISH +2.89 XHAM +
1.89 XMCH + 1.99 XMTL + 1.99 XSPG + 2.49 XTUR
Subject to
60 XBEEF + 8 XCHK + 8 XFISH + 40 XHAM +
15 XMCH + 70 XMTL + 25 XSPG + 60 XTUR >= 700
20 XBEEF + 0 XCHK + 10 XFISH + 40 XHAM +
35 XMCH + 30 XMTL + 50 XSPG + 20 XTUR >= 700
10 XBEEF + 20 XCHK + 15XFISH + 35 XHAM +
15 XMCH + 15 XMTL + 25 XSPG + 15 XTUR >= 700
15 XBEEF + 20XCHK + 10 XFISH + 10 XHAM +
15 XMCH + 15XMTL + 15 XSPG + 10 XTUR >= 700
XBEEF, XCHK, XFISH, XHAM , XMCH, XMTL, XSPG , XTUR >= 0
Another way to describe this linear program is as follows:
Given: F, a set of foods
N, a set of nutrients
cj = cost per package of dinner j, for each j in F
pi,j = percentage of nutrient i per package of dinner j, for each i in N, j in F
lj = minimum weekly requirement of nutrient j, for each j in F
Decision variables: Xj = packages of dinner j to be purchased, for each j in F
Minimize:
Subject to:
3. AMPL Model
This model deals with two things: nutrients and foods. Thus we begin with an AMPL model by declaring sets of each:
set NUTR;
set FOOD;
Next we need to specify the numbers required by the model. Certainly a positive cost should be given for each food:
param cost {FOOD} > 0;
The cost of a particular food is written as, say, cost[BEEF], although references to specific values in an LP model are usually in terms of indices like j , rather than specific set members like BEEF.
To make the model somewhat more general than the LPs so far, we will also specify that for each food there are lower and upper limits on the number of packages in the diet:
param f_min {FOOD} >= 0;
param f_max {j in FOOD} >= f_min[j];
Notice that we need a dummy index j to run over FOOD in the declaration of f_max, in order to say that the maximum for each food must be greater than or equal to the corresponding minimum.
We will also find it convenient to specify similar lower and upper limits on the amount of each nutrient in the diet:
param n_min {NUTR} >= 0;
param n_max {i in NUTR} >= n_min[i];
Finally, for each combination of a nutrient and a food, we need a number that represents the amount of each nutrient i in a package of the food. Such a "product" of two sets is written by listing them both:
param amt {NUTR,FOOD} >= 0;
References to this parameter require two indices. For example, amt[i,j]is the amount of nutrient i in a package of j.
The decision variables for this model are the numbers of packages to buy of the different foods:
var Buy {j in FOOD} >= f_min[j], <= f_max[j];
The number of packages of some food j to be bought will be called Buy[j] ; in any acceptable solution it will have to lie between f_min[j] and f_max[j].
The total cost of buying a food j is the cost per package, cost[j], times the number of packages. Buy[j]. The objective to be minimized is the sum of this product over all foods j:
minimize total_cost: sum {j in FOOD} cost[j] * Buy[j];
Similarly, the amount of a nutrient i supplied by a food j is the nutrient per package, amt[i,j], times the number of package Buy[j]. The total amount of nutrient i supplied is the sum of this product over all food j:
sum {j in FOOD} amt[i,j] * Buy[j]
To complete the model, we need only specify that each such sum must lie between the appropriate bounds. Our constraint declaration begins
subject to diet {i in NUTR}:
to say that a constraint named diet[i] must be imposed for each member i of NUTR. The rest of the declaration gives the algebraic statement of the constraint for nutrient i: the variables must satisfy
n_min[i] <= sum {j in FOOD} amt[i,j] * Buy[j] <= n_max[i];
A "double inequality" like this is interpreted in the obvious way: the value of the sum in the middle must lie between n_min[i] and n_max[i]. The following is the complete model.
set NUTR;
set FOOD;
param cost {FOOD} > 0;
param f_min {FOOD} >= 0;
param f_max {j in FOOD} >= f_min[j];
param n_min {NUTR} >= 0;
param n_max {i in NUTR} >= n_min[i];
param amt {NUTR,FOOD} >= 0;
var Buy {j in FOOD} >= f_min[j], <= f_max[j];
minimize total_cost: sum {j in FOOD} cost[j] * Buy[j];
subject to diet {i in NUTR}:
n_min[i] <= sum {j in FOOD} amt[i,j] * Buy[j] <= n_max[i];
4. AMPL Data
By specifying appropriate data, we can solve any of the linear programs that correspond to the above model. Let’s begin by using the data from the beginning of this document. The following is the data in AMPL format.
set NUTR := A B1 B2 C ;
set FOOD := BEEF CHK FISH HAM MCH MTL SPG TUR ;
param: cost f_min f_max :=
BEEF 3.19 0 100
CHK 2.59 0 100
FISH 2.29 0 100
HAM 2.89 0 100
MCH 1.89 0 100
MTL 1.99 0 100
SPG 1.99 0 100
TUR 2.49 0 100 ;
param: n_min n_max :=
A 700 10000
C 700 10000
B1 700 10000
B2 700 10000 ;
param amt (tr):
A C B1 B2 :=
BEEF 60 20 10 15
CHK 8 0 20 20
FISH 8 10 15 10
HAM 40 40 35 10
MCH 15 35 15 15
MTL 70 30 15 15
SPG 25 50 25 15
TUR 60 20 15 10
;
The value of f_min and n_min are as given originally, while f_max and n_max are set, for the time being, to large values that won’t affect the optimal solution. In the table for amt, the notation(tr)indicates that we have "transposed" the table so the columns correspond to the first index (nutrients), and the rows to the second (foods). Alternatively, we could have changed the model to say
param amt {FOOD,NUTR}
in which case we would have had to write amt[j,i] in the constraint.
Suppose that model and data are stored in the files diet.mod and diet.dat respectively. Then AMPL is used as follows to read these files and to solve the resulting linear program:
ampl: model diet.mod;
ampl: data diet.dat;
ampl: solve;
CPLEX 2.0: optimal solution; objective 88.2
2 iterations (1 in phase I)
ampl: display Buy;
Buy [*] :=
BEEF 0
CHK 0
FISH 0
HAM 0
MCH 46.6667
MTL 0
SPG 0
TUR 0
;
Now suppose we want to make the following enhancements. To promote variety, the weekly diet must contain between 2 and 10 packages of each food. The amount of sodium and calories in each package is also given; total sodium must not exceed 40,000 mg, and total calories must be between 16,000 and 24,000. All of these changes can be made through a few modifications to the data. The following is the new AMPL data set after the changes.
set NUTR := A B1 B2 C NA CAL ;
set FOOD := BEEF CHK FISH HAM MCH MTL SPG TUR ;
param: cost f_min f_max :=
BEEF 3.19 2 10
CHK 2.59 2 10
FISH 2.29 2 10
HAM 2.89 2 10
MCH 1.89 2 10
MTL 1.99 2 10
SPG 1.99 2 10
TUR 2.49 2 10 ;
param: n_min n_max :=
A 700 20000
C 700 20000
B1 700 20000
B2 700 20000
NA 0 40000
CAL 16000 24000 ;
param amt (tr):
A C B1 B2 NA CAL :=
BEEF 60 20 10 15 938 295
CHK 8 0 20 20 2180 770
FISH 8 10 15 10 945 440
HAM 40 40 35 10 278 430
MCH 15 35 15 15 1182 315
MTL 70 30 15 15 896 400
SPG 25 50 25 15 1329 370
TUR 60 20 15 10 1397 450 ;
Putting this new data in file diet2.dat, we can run AMPL again:
ampl: model diet.mod;
ampl: data diet2.dat;
ampl: solve;
CPLEX 2.0: infeasible problem
7 iterations (7 in phase I)
The message infeasible tells us that we have constrained the diet too tightly; there is no way that all of the restriction can be satisfied.
AMPL let us examine a variety of values produced by a solver as it attempts to find a solution. We can look for the source of the infeasibility by displaying some values associated with this solution.
ampl: display diet.lb, diet.body, diet.ub;
: diet.lb diet.body diet.ub :=
A 700 1993.09 20000
B1 700 841.091 20000
B2 700 601.091 20000
C 700 1272.55 20000
CAL 16000 17222.9 24000
NA 0 40000 40000
;
The value for diet.lb and diet.ub are the "lower bounds" and "upper bounds" on the sum of the variables in the constraints diet[i] – in this case, just the values n_min[i] and n_max[i] ; diet.body is the current sum of the variables. We can see that the diet does not supply enough vitamin B2, while the amount of sodium(NA) has reached its upper bound. If we relax the sodium limit to 50,000 mg, a feasible solution becomes possible:
ampl: let n_max["NA"] := 50000; solve;
CPLEX 2.0: optimal solution; objective 118.0594032
11 iterations (7 in phase I)
ampl: display B
Buy [*] :=
BEEF 5.36061
CHK 2
FISH 2
HAM 10
MCH 10
MTL 10
SPG 9.30605
TUR 2
;
This is at least a start toward a palatable diet, although we have to spend $118.06, compared to $88.20 for the original, less restricted case. Clearly, it would be easy, now that the model is set up, to try many other possibilities. Once still disappointing aspect of the solution is the need to buy 5.36061 packages of beef, and 9.30365 of spaghetti. What if we can only buy whole packages? You might think that we could just round the optimal values to whole numbers, but it is not so easy to do so in a feasible way. Displaying the constraint bounds at the optimum.
ampl: display diet.lb, diet.body, diet.ub;
: diet.lb diet.body diet.ub :=
A 700 1956.29 20000
B1 700 1036.26 20000
B2 700 700 20000
C 700 1682.51 20000
CAL 16000 19794.6 24000
NA 0 50000 50000
;
we see that, for example, 6 packages of beef and 10 packages of spaghetti will violate the sodium limit, while 5 packages of beef and 9 packages of spaghetti provide insufficient vitamin B2. Even if we could find a nearby all-integer solution that satisfies the constraints, there would be no guarantee that it would be the least-cost all-integer solution. AMPL does provide for putting the integrality restriction directly into the declaration of the variables:
var Buy {j in FOOD} integer >= f_min[j], <= f_max[j];
This will only help, however, if you use a solver that can handle so-called integer programs. In general, integrality and other "discrete" restrictions make a model much harder to solve.