Goal Programming

In [1] a goal programming [2] model was stated, although the poster did not use that term.

The model is fairly simple:

• We have 6 different raw materials that need to be blended into a final product.
• The raw materials have 4 different ingredients.
• We want to blend one unit of final product that has some target values for some formulas. We have three of these goals. 

The data looks like:

----     17 SET b  raw materials

r1,    r2,    r3,    r4,    r5,    r6


----     17 SET i  ingredients

c,    s,    a,    f


----     17 SET t  targets

target1,    target2,    target3


----     17 PARAMETER data  raw materials data table

             c           s           a           f          lo          up

r1      51.290       4.160       0.970       0.380       0.300       0.370
r2      51.420       4.160       0.950       0.370       0.300       0.370
r3       6.880      63.360      13.580       3.060       0.070       0.130
r4      32.050       1.940                               0.005       0.015
r5       4.560      21.430       3.820      52.280       0.010       0.030
r6       0.190       7.450       4.580       0.420       0.110       0.130

The data table has the constituent ingredients of the raw materials and the bounds for raw material usage. The blending part of the model is simply: \[\begin{align}&\sum_r\color{darkred}x_r=1 \\&\color{darkred}y_i = \sum_r \color{darkred}x_r \cdot \color{darkblue}{\mathit{data}}_{r,i} \\ & \color{darkred}x_r \in [\color{darkblue}\ell_r,\color{darkblue}u_r] \end{align}\] Here the variable \(\color{darkred}x_r\) indicates the amount of raw material used and \(\color{darkred}y_i\) indicates the composition of the final product. 

The goals as stated in [1] are: \[\begin{align}\frac{100\cdot\color{darkred}y_c}{2.8\cdot\color{darkred}y_s +1.18\cdot \color{darkred}y_a+0.65\cdot \color{darkred}y_f} &  \approx 98.5 \\  \frac{\color{darkred}y_s}{\color{darkred}y_a+\color{darkred}y_f}    & \approx 2.5 \\   \frac{\color{darkred}y_a}{\color{darkred}y_f} &\approx 1.3 \end{align}\]

After adding slacks or deviations to these targets we can try to minimize them somehow. An obvious choice would be:

• Minimize the weighted sum of the squared deviations

This is by no means the only possible approach. But it is a sensible way to start analyzing this. A direct nonlinear formulation can look like:

 

Non-linear Model

\[\begin{align}\min & \sum_t \color{darkblue}w_t \cdot \color{darkred}\Delta_t^2 \\ &\frac{100\cdot\color{darkred}y_c}{2.8\cdot\color{darkred}y_s +1.18\cdot \color{darkred}y_a+0.65\cdot \color{darkred}y_f} = 98.5 + \color{darkred}\Delta_1 \\ & \frac{\color{darkred}y_s}{\color{darkred}y_a+\color{darkred}y_f} = 2.5 + \color{darkred}\Delta_2 \\ & \frac{\color{darkred}y_a}{\color{darkred}y_f} = 1.3 + \color{darkred}\Delta_3 \\ & \sum_r\color{darkred}x_r=1 \\&\color{darkred}y_i = \sum_r \color{darkred}x_r \cdot \color{darkblue}{\mathit{data}}_{r,i} \\ & \color{darkred}x_r \in [\color{darkblue}\ell_r,\color{darkblue}u_r] \end{align}\]

We can say a few things about this NLP model.

• This model has a quadratic objective and non-linear constraints.
• We can substitute out the variables \(\color{darkred}\Delta_t\). This has the advantage that only the objective is nonlinear. A model with a nonlinear objective and linear constraints is usually easier to solve. A good presolver may do this for us.
• We need to protect against division by zero. One ad-hoc approach is to use a nonzero initial value for \(\color{darkred}y_s,\color{darkred}y_a,\color{darkred}y_f\). A more interesting approach is using the fact that there are some inherent lower bounds on  \(\color{darkred}y_s,\color{darkred}y_a,\color{darkred}y_f\). For instance we can do: \[\begin{align}\min\>& \color{darkred}y_s\\ &\sum_r\color{darkred}x_r=1 \\ & \color{darkred}y_i = \sum_r \color{darkred}x_r \cdot \color{darkblue}{\mathit{data}}_{r,i} \\ & \color{darkred}x_r \in [\color{darkblue}\ell_r,\color{darkblue}u_r] \end{align} \] Similar for \(\color{darkred}y_a\) and \(\color{darkred}y_f\). When we do this, we find the following implied lower bounds:
  
  

  ----     55 PARAMETER lb  lowerbound on y(saf)
  
  s 10.105,    a  2.575,    f  1.176
  

  
  
  Adding these bounds to the corresponding variables will protect the division. NLP solvers will typically not evaluate non-linear functions with variables outside their bounds.
• We can solve this model with a local solver or with a global solver. If we want we can use the previous method to find good lower and upper bounds on all variables \(\color{darkred}y_i\).

A solution can look like:

----     89 VARIABLE x.L  amount of raw materials to use

r1 0.354,    r2 0.370,    r3 0.123,    r4 0.015,    r5 0.028,    r6 0.110


----     89 VARIABLE y.L  composition of output

c 38.654,    s 12.260,    a  2.977,    f  2.157


----     89 VARIABLE d.L  deviation of target

target1 -2.83830E-4,    target2      -0.112,    target3       0.081


----     89 VARIABLE z.L                   =        0.019  objective

We can linearize the constraints by observing that the target equations can be interpreted differently. We can write them as:  \[\begin{align} \color{darkred}y_c &  \approx 0.01 \cdot 98.5 \cdot \left( 2.8\cdot\color{darkred}y_s +1.18\cdot \color{darkred}y_a+0.65\cdot \color{darkred}y_f \right) \\  \color{darkred}y_s   & \approx 2.5 \cdot \left(\color{darkred}y_a+\color{darkred}y_f \right) \\   \color{darkred}y_a  &\approx 1.3 \cdot \color{darkred}y_f\end{align}\]

The resulting QP model (with a convex quadratic objective and linear constraints) can look like:

Quadratic Model

\[\begin{align}\min & \sum_t \color{darkblue}w_t \cdot \color{darkred}\Delta_t^2 \\ & \color{darkred}y_c  = 0.01 \cdot 98.5 \cdot \left( 2.8\cdot\color{darkred}y_s +1.18\cdot \color{darkred}y_a+0.65\cdot \color{darkred}y_f \right) + \color{darkred}\Delta_1 \\ & \color{darkred}y_s   = 2.5 \cdot \left(\color{darkred}y_a+\color{darkred}y_f \right) + \color{darkred}\Delta_2 \\ & \color{darkred}y_a  = 1.3 \cdot \color{darkred}y_f+ \color{darkred}\Delta_3 \\ & \sum_r\color{darkred}x_r=1 \\&\color{darkred}y_i = \sum_r \color{darkred}x_r \cdot \color{darkblue}{\mathit{data}}_{r,i} \\ & \color{darkred}x_r \in [\color{darkblue}\ell_r,\color{darkblue}u_r] \end{align}\]

Notice that we have redefined \( \color{darkred}\Delta_t\), so the weights \(\color{darkblue}w_t\) should be different. This is now a convex QP model, which is easier to solve.

The output is:

----    110 VARIABLE x.L  amount of raw materials to use

r1 0.355,    r2 0.370,    r3 0.125,    r4 0.015,    r5 0.025,    r6 0.110


----    110 VARIABLE y.L  composition of output

c 38.713,    s 12.310,    a  2.990,    f  2.013


----    110 VARIABLE d.L  deviation of target

target1 -0.003,    target2 -0.198,    target3  0.374


----    110 VARIABLE z.L                   =        0.179  objective

The result is close to the solution of the nonlinear model.

We can further linearize the model by introducing positive and negative slacks. This will allow us to form absolute deviations.

Linear Model

\[\begin{align}\min & \sum_t \color{darkblue}w_t \cdot \left( \color{darkred}E_t^+ + \color{darkred}E_t^- \right) \\ & \color{darkred}y_c  = 0.01 \cdot 98.5 \cdot \left( 2.8\cdot\color{darkred}y_s +1.18\cdot \color{darkred}y_a+0.65\cdot \color{darkred}y_f \right) + \color{darkred}E_1^+ - \color{darkred}E_1^- \\ & \color{darkred}y_s   = 2.5 \cdot \left(\color{darkred}y_a+\color{darkred}y_f \right) + \color{darkred}E_2^+ - \color{darkred}E_2^-\\ & \color{darkred}y_a  = 1.3 \cdot \color{darkred}y_f+ \color{darkred}E_3^+ -\color{darkred}E_3^-\\ & \sum_r\color{darkred}x_r=1 \\&\color{darkred}y_i = \sum_r \color{darkred}x_r \cdot \color{darkblue}{\mathit{data}}_{r,i} \\ & \color{darkred}x_r \in [\color{darkblue}\ell_r,\color{darkblue}u_r] \\ &  \color{darkred}E_t^+, \color{darkred}E_t^- \ge 0  \end{align}\]

This is now an LP that can be solved with any LP solver. The results are:

----    137 VARIABLE x.L  amount of raw materials to use

r1 0.358,    r2 0.370,    r3 0.130,    r4 0.015,    r5 0.016,    r6 0.111


----    137 VARIABLE y.L  composition of output

c 38.869,    s 12.459,    a  3.033,    f  1.536


----    137 VARIABLE e.L  pair of slacks

                  +

target2       1.036


----    137 VARIABLE z.L                   =        1.036  objective

One main difference between the solutions is that the quadratic model wants to distribute deviations over all targets. The linear model is reporting a corner solution, with many deviations equal to zero. This means one target is paying the price. This is typical behavior for linear models.

Conclusion

I have not discussed the lexicographic approach, where we have a hierarchy of goals. But even just using a weighted sum approach, we have quite some alternatives to consider.

References

1. How to solve multiple multivariate equation systems with constraints, https://stackoverflow.com/questions/72682801/how-to-solve-multiple-multivariate-equation-systems-with-constraints
2. Goal programming, https://en.wikipedia.org/wiki/Goal_programming

Appendix: GAMS model

$ontext     Goal programming example. $offtext set    r 'raw materials' /r1*r6/    i 'ingredients' /c,s,a,f/    t 'targets' /target1*target3/ ; table data(r,*) 'raw materials data table'       lo    up      c     s     a     f r1  0.300  0.370  51.29  4.16  0.97  0.38 r2  0.300  0.370  51.42  4.16  0.95  0.37 r3  0.070  0.130   6.88 63.36 13.58  3.06 r4  0.005  0.015  32.05  1.94 r5  0.010  0.030   4.56 21.43  3.82 52.28 r6  0.110  0.130   0.19  7.45  4.58  0.42 ; display r,i,t,data; *-------------------------------------------------------------------------- * establish lower bounds on y(s),y(a),y(f) *-------------------------------------------------------------------------- variables    x(r)  'amount of raw materials to use'    y(i)  'composition of output'    z     'objective' ; x.lo(r) = data(r,'lo'); x.up(r) = data(r,'up'); parameter s(i) 'select one y'; equations     obj1       'minimize one of y(s),y(a),y(f)'     blend(i)   'mixing'     unit       'one unit of output' ; obj1..     z =e= sum(i,s(i)*y(i)); blend(i).. y(i) =e= sum(r, x(r)*data(r,i)); unit..     sum(r, x(r)) =e= 1; model m1 /obj1,blend,unit/; set saf(i) /s,a,f/; parameter lb(saf) 'lowerbound on y(saf)'; loop(saf,    s(saf) = 1;    solve m1 minimizing z using lp;    lb(saf) = z.l;    s(saf) = 0; ); display lb; y.lo(saf) = lb(saf); *-------------------------------------------------------------------------- * nonlinear formulation *-------------------------------------------------------------------------- parameter w(t) 'weights'; w(t) = 1; variable    d(t)  'deviation of target' ; equations    nltarget1    nltarget2    nltarget3    obj2 ; nltarget1.. 100*y('c')/(2.8*y('s')+1.18*y('a')+0.65*y('f')) =e= 98.5 + d('target1'); nltarget2.. y('s')/(y('a')+y('f')) =e= 2.5 + d('target2'); nltarget3.. y('a')/y('f') =e= 1.3 + d('target3'); obj2.. z =e= sum(t,w(t)*sqr(d(t))); model m2 /obj2,blend,unit,nltarget1,nltarget2,nltarget3/; solve m2 minimizing z using nlp; *option nlp=baron,optcr=0; solve m2 minimizing z using nlp; display x.l,y.l,d.l,z.l; *-------------------------------------------------------------------------- * quadratic formulation *-------------------------------------------------------------------------- equations    lintarget1    lintarget2    lintarget3 ; lintarget1.. y('c') =e= 0.01*98.5*(2.8*y('s')+1.18*y('a')+0.65*y('f')) + d('target1'); lintarget2.. y('s') =e= 2.5*(y('a')+y('f')) + d('target2'); lintarget3.. y('a') =e= 1.3*y('f') + d('target3'); model m3 /obj2,blend,unit,lintarget1,lintarget2,lintarget3/; option qcp=cplex; solve m3 minimizing z using qcp; display x.l,y.l,d.l,z.l; *-------------------------------------------------------------------------- * linear formulation *-------------------------------------------------------------------------- set pm /'+','-'/; positive variable e(t,pm) 'pair of slacks'; equations    abstarget1    abstarget2    abstarget3    obj3 ; abstarget1.. y('c') =e= 0.01*98.5*(2.8*y('s')+1.18*y('a')+0.65*y('f')) + e('target1','+') - e('target1','-'); abstarget2.. y('s') =e= 2.5*(y('a')+y('f')) + e('target2','+') - e('target2','-'); abstarget3.. y('a') =e= 1.3*y('f') + e('target2','+') - e('target2','-'); obj3.. z =e= sum((t,pm),w(t)*e(t,pm)); model m4 /obj3,blend,unit,abstarget1,abstarget2,abstarget3/; solve m4 minimizing z using lp; display x.l,y.l,e.l,z.l;