Part 4: Troubleshoot part 3

Previously we have coded the algorithm with nloptr function. A few issues have popped up: some diet do not fall within the constraints. Possible reasons:

overly progressive goal of reducing GHGE by 20% (0.8 times original)
numeric issues that do not treat each constraint equally (e.g. prioritizes energy, ignores GHGE which we care the most)
convergence of the current package

Steps taken so far (2023.08.27)

re-scale the values for different categories, run again
set a moderate reduction of GHGE, see if it works
identify what parameters to adjust
- different scaling to prioritize (so far used standard deviation)
- use soft rather hard constraints
- different values for lower/upper boundary of optimized diet
- …

library(data.table)
foods <- read.csv('data/foods.csv', sep = ',')
setDT(foods) # use data.table format

# start with 3 foods
fd <- foods[food %in% c('Bread', 'Vegetables', 'Red meat')]
fd

         food intake energy protein   fat carbs sugar alcohol  ghge
       <char>  <num>  <num>   <num> <num> <num> <num>   <num> <num>
1:      Bread  175.4 10.696   0.091 0.030 0.441 0.002       0 0.001
2: Vegetables  154.6  1.565   0.015 0.008 0.050 0.005       0 0.001
3:   Red meat  117.6  8.342   0.173 0.139 0.014 0.000       0 0.013

# split current diet (grams of intake) and contribution per gram
current_diet <- fd$intake
current_diet

[1] 175.4 154.6 117.6

contrib_pergram <- fd[, .(energy, protein, fat, carbs, sugar, alcohol, ghge)]
contrib_pergram

   energy protein   fat carbs sugar alcohol  ghge
    <num>   <num> <num> <num> <num>   <num> <num>
1: 10.696   0.091 0.030 0.441 0.002       0 0.001
2:  1.565   0.015 0.008 0.050 0.005       0 0.001
3:  8.342   0.173 0.139 0.014 0.000       0 0.013

Constraints for 3 foods

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. 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.

# 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(current_diet)) %*% as.matrix(contrib_pergram)
const_max_3foods

       energy protein     fat   carbs  sugar alcohol   ghge
[1,] 3099.047 38.6252 22.8452 86.7278 1.1238       0 1.8588

# exclude sugar, alcohol
const_max_3foods <- const_max_3foods[, c('energy', 'protein', 'fat', 'carbs', 'ghge')]

# range of target: 0.9-1 of the maximum
# this is one thing that we could adjust: hard threshould to soft
const_lwrupr <- rbind(const_max_3foods*0.9, const_max_3foods*1)
rownames(const_lwrupr) <- c('lwr', 'upr')
const_lwrupr <- data.table(const_lwrupr)
const_lwrupr

     energy  protein      fat    carbs    ghge
      <num>    <num>    <num>    <num>   <num>
1: 2789.142 34.76268 20.56068 78.05502 1.67292
2: 3099.047 38.62520 22.84520 86.72780 1.85880

For example, the total energy for the 3 foods together should be within [2789.14, 3099.05].

Standardize food contribution per gram

Since the range of constraints for 5 categories differ hugely, it could affect the numeric evaluation. We want them to be on comparative scales.

The current solution is to standardize the contribution in each category (e.g. energy) by its original value divided by the standard deviation. For example,

the contribution of energy per food becomes (2.25, 0.33, 1.76) rather than (10.69, 1.56, 8.34)
the contribution of ghge per food becomes (0.14, 0.14, 1.87) rather than (0.001, 0.001, 0.013)

Similarly, the upper and lower limit of the constraints also need to be re-scaled. After rescaling, the target will be on a range of hundreds, rather than 3000 vs 1.8.

# find sd for each category: energy, protein.. ghge
# divide by these coef
contrib_pergram <- contrib_pergram[, c('energy', 'protein', 'fat', 'carbs', 'ghge')]
contrib_pergram

   energy protein   fat carbs  ghge
    <num>   <num> <num> <num> <num>
1: 10.696   0.091 0.030 0.441 0.001
2:  1.565   0.015 0.008 0.050 0.001
3:  8.342   0.173 0.139 0.014 0.013

sd_coef <- apply(contrib_pergram, MARGIN = 2, sd)

contrib_pergram_std <- sweep(contrib_pergram, MARGIN = 2, 1/sd_coef, FUN = '*')
contrib_pergram_std

     energy   protein       fat     carbs      ghge
1 2.2562170 1.1516220 0.4276560 1.8621635 0.1443376
2 0.3301215 0.1898278 0.1140416 0.2111296 0.1443376
3 1.7596637 2.1893473 1.9814728 0.0591163 1.8763884

# standardize constraint
# test the previous constraint
const_lwrupr

     energy  protein      fat    carbs    ghge
      <num>    <num>    <num>    <num>   <num>
1: 2789.142 34.76268 20.56068 78.05502 1.67292
2: 3099.047 38.62520 22.84520 86.72780 1.85880

const_lwrupr_std <- sweep(const_lwrupr, MARGIN = 2, 1/sd_coef, FUN = '*')
const_lwrupr_std

    energy  protein      fat    carbs     ghge
1 588.3423 439.9282 293.0966 329.5946 241.4652
2 653.7137 488.8091 325.6629 366.2162 268.2947

In the future, different scaling factor can be applied; but it should be a positive number after scaling. Could try dividing the difference between max and min of this variable.

Run optimization

Now we are going to solve the optimization problem using nloptr. The details of the algorithm please refer to part 3.

The objective function remains the same for 3 foods. We want to minimize the difference between the current and optimized diet.

current_diet

[1] 175.4 154.6 117.6

objective <- function(x)
{
  return ( (x[1]- current_diet[1])^2 + 
             (x[2]- current_diet[2])^2 + 
             (x[3]- current_diet[3])^2)
}

We optimize using 5 categories (10 inequality constraints):

energy
protein
fat
carbs
ghge

We show two examples: one reduces ghge to 90% what the current diet produces; the other reduces ghge to 85%.

Example 1: reduce ghge to 90%

First we reduce the GHGE to 90%. (This is feasible, in contrast to the 80% reduction)

cstr <- copy(const_lwrupr_std) # this line is important: it keeps a copy of the original values
cstr$ghge <- cstr$ghge * 0.9 # reduce to 0.9
cstr

    energy  protein      fat    carbs     ghge
1 588.3423 439.9282 293.0966 329.5946 217.3187
2 653.7137 488.8091 325.6629 366.2162 241.4652

# use the standardized food contribution per gram 
# contrib_pergram_std

# define the inequality constraints
inequalconstr <- function (x) {
  
  cps <- contrib_pergram_std # call the values
  # cstr is used as the last item in each equation
  
  constr <- c(
    # energy
    - x[1]*cps$energy[1] - x[2]*cps$energy[2] - x[3]*cps$energy[3] + cstr$energy[1], # lower
    x[1]*cps$energy[1] + x[2]*cps$energy[2] + x[3]*cps$energy[3] - cstr$energy[2], # upper
    
    # protein
    - x[1]*cps$protein[1] - x[2]*cps$protein[2] - x[3]*cps$protein[3] + cstr$protein[1],
    x[1]*cps$protein[1] + x[2]*cps$protein[2] + x[3]*cps$protein[3] - cstr$protein[2],
    
    # fat
    - x[1]*cps$fat[1] - x[2]*cps$fat[2] - x[3]*cps$fat[3]+ cstr$fat[1],
    x[1]*cps$fat[1] + x[2]*cps$fat[2] + x[3]*cps$fat[3] - cstr$fat[2],

    # carbs
    - x[1]*cps$carbs[1] - x[2]*cps$carbs[2] - x[3]*cps$carbs[3]+ cstr$carbs[1],
    x[1]*cps$carbs[1] + x[2]*cps$carbs[2] + x[3]*cps$carbs[3] - cstr$carbs[2],
    
    # ghge
    - x[1]*cps$ghge[1] - x[2]*cps$ghge[2] - x[3]*cps$ghge[3]+ cstr$ghge[1],
    x[1]*cps$ghge[1] + x[2]*cps$ghge[2] + x[3]*cps$ghge[3] - cstr$ghge[2]
  )
  return (constr)
}

Set other parameters (this part is unchanged, apart from the number of constraints)

# 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)

opts <- list( "algorithm" = "NLOPT_GN_ISRES",
              "xtol_rel"= 1.0e-15,
              "maxeval"= 160000,
              "tol_constraints_ineq" = rep( 1.0e-10, 10 )) # match number of const

# run the algorithm
res <- nloptr::nloptr(
  x0          = x0,        # initial value for x
  eval_f      = objective, # objective function
  lb          = lb,        # lower bound for x
  ub          = ub,        # upper bound for x
  eval_g_ineq = inequalconstr, # inequality constraint
  opts        = 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: 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:  10 
Number of equality constraints:    0 
Current value of objective function:  202.054848606764 
Current value of controls: 174.3086 153.5102 103.4693

Now we print out the results in a more readable way

res_diet <- res$solution

# print with the current diet, and change percentage
# also print the boundary, in case it hit boundary

diet_result <- data.frame(
  name = c('Bread', 'Vegetables', 'Red meat'),
  current = current_diet, 
  new = res_diet, 
  percent_change = round((res_diet - current_diet)/current_diet, 3),
  lower_limit = lb, 
  upper_limit = ub
)
diet_result

        name current      new percent_change lower_limit upper_limit
1      Bread   175.4 174.3086         -0.006         160         180
2 Vegetables   154.6 153.5102         -0.007         140         160
3   Red meat   117.6 103.4693         -0.120         100         120

# verify whether it falls within 
output_newdiet <- t(as.matrix(res_diet)) %*% as.matrix(contrib_pergram_std)
# output_newdiet

# cstr
const_result <- t(rbind(output_newdiet, cstr))
colnames(const_result) <- c('new_diet','const_lwr', 'const_upr')
const_result <- data.table(const_result)
# conditions
const_result[, is_ok := 'Yes']
const_result[new_diet < const_lwr, is_ok := 'beyond lower']
const_result[new_diet > const_upr, is_ok := 'beyond upper']

# relative difference (since we rescaled the targets)
const_result[, relative_dev := 0]
const_result[is_ok == 'beyond lower', relative_dev := round((new_diet - const_lwr)/const_lwr, 3)]
const_result[is_ok == 'beyond upper', relative_dev := round((new_diet - const_upr)/const_upr, 3)]

# print out
const_result

Index: <is_ok>
   new_diet const_lwr const_upr  is_ok relative_dev
      <num>     <num>     <num> <char>        <num>
1: 626.0262  588.3423  653.7137    Yes            0
2: 456.4084  439.9282  488.8091    Yes            0
3: 297.0723  293.0966  325.6629    Yes            0
4: 363.1183  329.5946  366.2162    Yes            0
5: 241.4652  217.3187  241.4652    Yes            0

You should always check how much it actually deviates from the target constraints.

Example 2: reduce ghge to 85%

This is one example where the optimal constraints can not be reached. You’ll see that the boundaries that we set in the algorithm have been reached; suggesting that if we modify these parameters, the solution could still be found.

Be careful with HOW MUCH it deviates from the target!

cstr <- copy(const_lwrupr_std) # this line is important: it keeps a copy of the original values
cstr$ghge <- cstr$ghge * 0.85 # reduce to 0.85
cstr

    energy  protein      fat    carbs     ghge
1 588.3423 439.9282 293.0966 329.5946 205.2454
2 653.7137 488.8091 325.6629 366.2162 228.0505

# use the standardized food contribution per gram 
# contrib_pergram_std

# define the inequality constraints
inequalconstr <- function (x) {
  
  cps <- contrib_pergram_std # call the values
  # cstr is used as the last item in each equation
  
  constr <- c(
    # energy
    - x[1]*cps$energy[1] - x[2]*cps$energy[2] - x[3]*cps$energy[3] + cstr$energy[1], # lower
    x[1]*cps$energy[1] + x[2]*cps$energy[2] + x[3]*cps$energy[3] - cstr$energy[2], # upper
    
    # protein
    - x[1]*cps$protein[1] - x[2]*cps$protein[2] - x[3]*cps$protein[3] + cstr$protein[1],
    x[1]*cps$protein[1] + x[2]*cps$protein[2] + x[3]*cps$protein[3] - cstr$protein[2],
    
    # fat
    - x[1]*cps$fat[1] - x[2]*cps$fat[2] - x[3]*cps$fat[3]+ cstr$fat[1],
    x[1]*cps$fat[1] + x[2]*cps$fat[2] + x[3]*cps$fat[3] - cstr$fat[2],

    # carbs
    - x[1]*cps$carbs[1] - x[2]*cps$carbs[2] - x[3]*cps$carbs[3]+ cstr$carbs[1],
    x[1]*cps$carbs[1] + x[2]*cps$carbs[2] + x[3]*cps$carbs[3] - cstr$carbs[2],
    
    # ghge
    - x[1]*cps$ghge[1] - x[2]*cps$ghge[2] - x[3]*cps$ghge[3]+ cstr$ghge[1],
    x[1]*cps$ghge[1] + x[2]*cps$ghge[2] + x[3]*cps$ghge[3] - cstr$ghge[2]
  )
  return (constr)
}

# 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)

opts <- list( "algorithm" = "NLOPT_GN_ISRES",
              "xtol_rel"= 1.0e-15,
              "maxeval"= 160000,
              "tol_constraints_ineq" = rep( 1.0e-10, 10 )) # match number of const

# run the algorithm
res <- nloptr::nloptr(
  x0          = x0,        # initial value for x
  eval_f      = objective, # objective function
  lb          = lb,        # lower bound for x
  ub          = ub,        # upper bound for x
  eval_g_ineq = inequalconstr, # inequality constraint
  opts        = 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: 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:  10 
Number of equality constraints:    0 
Current value of objective function:  528.456412600945 
Current value of controls: 177.753 140 100

Now we check the results for the second problem.

res_diet <- res$solution

# print with the current diet, and change percentage
# also print the boundary, in case it hit boundary

diet_result <- data.frame(
  name = c('Bread', 'Vegetables', 'Red meat'),
  current = current_diet, 
  new = res_diet, 
  percent_change = round((res_diet - current_diet)/current_diet, 3),
  lower_limit = lb, 
  upper_limit = ub
)
diet_result

        name current     new percent_change lower_limit upper_limit
1      Bread   175.4 177.753          0.013         160         180
2 Vegetables   154.6 140.000         -0.094         140         160
3   Red meat   117.6 100.000         -0.150         100         120

# verify whether it falls within 
output_newdiet <- t(as.matrix(res_diet)) %*% as.matrix(contrib_pergram_std)
# output_newdiet

# cstr
const_result <- t(rbind(output_newdiet, cstr))
colnames(const_result) <- c('new_diet','const_lwr', 'const_upr')
const_result <- data.table(const_result)
# conditions
const_result[, is_ok := 'Yes']
const_result[new_diet < const_lwr, is_ok := 'beyond lower']
const_result[new_diet > const_upr, is_ok := 'beyond upper']

# relative difference (since we rescaled the targets)
const_result[, relative_dev := 0]
const_result[is_ok == 'beyond lower', relative_dev := round((new_diet - const_lwr)/const_lwr, 3)]
const_result[is_ok == 'beyond upper', relative_dev := round((new_diet - const_upr)/const_upr, 3)]

# print out
const_result

Index: <is_ok>
   new_diet const_lwr const_upr        is_ok relative_dev
      <num>     <num>     <num>       <char>        <num>
1: 623.2326  588.3423  653.7137          Yes        0.000
2: 450.2148  439.9282  488.8091          Yes        0.000
3: 290.1302  293.0966  325.6629 beyond lower       -0.010
4: 366.4749  329.5946  366.2162 beyond upper        0.001
5: 233.5025  205.2454  228.0505 beyond upper        0.024