Part 3: A simple problem with 3 foods and 3 constraints
We document the optimization procedure with nloptr using a small example of 3 foods, and 3 constraints (energy, protein, ghge).
foods <-read.csv('data/foods.csv', sep =',')data.table::setDT(foods) # use data.table format# we only take 3 foodsfd <- foods[food %in%c('Bread', 'Vegetables', 'Red meat')]fd
We compute the constraints for 3 foods together. Note that the consumption for each of the 28 food groups are different, hence we can not use 3/28 times the total energy, protein…; we need to use a weighted average.
The cosntraints can be set to any positive number that make sense.
# original constraint for all foods (maximum)const_max_allfoods <-c(9314.3, 98.2, 85.8, 234.7, 39.2, 8.6, 3.8)# 3 foods contribution (maximum)const_max_3foods <-t(as.matrix(fd$intake)) %*%as.matrix(fd[, .(energy, protein, fat, carbs, sugar, alcohol, ghge)])const_max_3foods
energy protein fat carbs sugar alcohol ghge
[1,] 3099.047 38.6252 22.8452 86.7278 1.1238 0 1.8588
In the optimization problem, we want to have a range (lower, upper) for the metrics. We can assume that the lower range is 90% of the maximum.
\(x_1, x_2, x_3\) are the target diet (grams or other units) we want to find, these three are unknown.
\(X_1, X_2, X_3\) are the current diet for bread, vegetables and red meat. We have these values in the data: 175.4, 154.6, 117.6. In the function we use fd$intake to programmatically extract the values.
Now we define the inequality constraints. For this demo, we only include 3 metrics (energy, protein, ghge) for simplicity.
# select the metrics we want to keep c3 <- c3foods[, c('energy', 'protein', 'ghge')]c3
energy protein ghge
lwr 2789.142 34.76268 1.67292
upr 3099.047 38.62520 1.85880
The inequality constraints need to be reformulated for nloptr to work. More specifically, nloptr optimizes a function \(f(x)\) subject to a set of equality and inequality constraints. For the inequality constraints, they need to be the form of \(g(x) <= 0\).
We do not need to explicitly specify the non-negativity constraints for \(x\) inside the function. They are specified outside the function as a set of lower and upper bounds. You can set these range to be anything that make sense; for now we make it close to the current diet.
Initial values are the values you tell the program to start the search. This is also up to you; sometimes initial values can have a large impact on the final results when the optimization problem is difficult (e.g. has multiple local optima). We can start close to the values of the current diet.
# lower and upper bounds of x (3 foods)lb <-c(160, 140, 100)ub <-c(180, 160, 120)# Initial values# (try different ones!)x0 <-c(175, 150, 110)
Options should also be supplied to the program as a named list. These are specifications for the algorithm, tolerance and maximum evaluation. For more details please check the function manual.
Pay attention to the last line here: the size inside rep() need to match the number of inequality constraints. In this example it should be 6.
Now we run the optimization by calling nloptr from nloptr package. If you have loaded the package before, you don’t need the :: between the two.
We print the result and check the optimized values.
# run the algorithmres <- nloptr::nloptr(x0 = x0, # initial value for xeval_f = objective, # objective functionlb = lb, # lower bound for xub = ub, # upper bound for xeval_g_ineq = inequalconstr, # inequality constraintopts = opts # options )print(res)
Call:
nloptr::nloptr(x0 = x0, eval_f = objective, lb = lb, ub = ub,
eval_g_ineq = inequalconstr, opts = opts)
Minimization using NLopt version 2.7.1
NLopt solver status: 4 ( NLOPT_XTOL_REACHED: Optimization stopped because
xtol_rel or xtol_abs (above) was reached. )
Number of Iterations....: 35522
Termination conditions: xtol_rel: 1e-15 maxeval: 160000
Number of inequality constraints: 6
Number of equality constraints: 0
Optimal value of objective function: 9.22904150236166e-26
Optimal value of controls: 175.4 154.6 117.6
The optimal values are 175.4, 154.6, 117.6. These are the same as the current diet.
Note
This result is not surprising; our current diet on the three foods already satisfies the constraints, and it is the closest to the current diet (\((x_1 - X_1)^2 +(x_2 - X_2)^2 + (x_3 - X_3)^2 = 0\)) when \(x_1 = X_1, x_2 = X_2, x_3 = X_3\).
Reduce ghge to 80%
We can try to modify the constraint values to see if the diet will be different.
For example, we can limit ghge. The current upper and lower range are [1.673, 1.859], let us try to limit it to [1.338, 1.487], which is 80% of the original values.
Now we run the optimization again with the new constraint, inequalconstr_alt.
# run the algorithmres_alt <- nloptr::nloptr(x0 = x0, # initial value for xeval_f = objective, # objective functionlb = lb, # lower bound for xub = ub, # upper bound for xeval_g_ineq = inequalconstr_alt, # inequality constraint (NEW)opts = opts # options )print(res_alt)
Call:
nloptr::nloptr(x0 = x0, eval_f = objective, lb = lb, ub = ub,
eval_g_ineq = inequalconstr_alt, opts = opts)
Minimization using NLopt version 2.7.1
NLopt solver status: 5 ( NLOPT_MAXEVAL_REACHED: Optimization stopped because
maxeval (above) was reached. )
Number of Iterations....: 160000
Termination conditions: xtol_rel: 1e-15 maxeval: 160000
Number of inequality constraints: 6
Number of equality constraints: 0
Current value of objective function: 566.401428410334
Current value of controls: 168.806 140 100
The new results are 168.806, 140, 100.
Compare the two runs
Bread: 175.4 -> 168.8 (3.7% reduction)
Vegetables: 154.6 -> 140 (9.4% reduction)
Red meat: 117.6 -> 100 (15.0% reduction)
This makes sense, as red meat is the largest contributor for ghge.
What we have shown here is a very crude demonstration of how to construct the optimization problem, and a sanity check of the results. A lot more can be tried out.