In [1], the following Pyomo model (Python fragment) is presented:

model.x = Var(name="Number of batches", domain=NonNegativeIntegers, initialize=10) model.a = Var(name="Batch Size", domain=NonNegativeIntegers, bounds=(5,20)) # Objective function def total_production(model): return model.x * model.a model.total_production = Objective(rule=total_production, sense=minimize) # Constraints # Minimum production of the two output products def first_material_constraint_rule(model): return sum(0.2 * model.a * i for i in range(1, value(model.x)+1)) >= 70 model.first_material_constraint = Constraint(rule=first_material_constraint_rule) def second_material_constraint_rule(model): return sum(0.8 * model.a * i for i in range(1, value(model.x)+1)) >= 90 model.second_material_constraint = Constraint(rule=second_material_constraint_rule) # At least one production run def min_production_rule(model): return model.x >= 1 model.min_production = Constraint(rule=min_production_rule)

This is a little bit incomplete: we miss imports, creating the model and a solve. However, we can see some real problematic issues here. The main problem is the use of value(model.x) inside a constraint. This is almost never intended, as this evaluates the initial point and not the current value during optimization. In GAMS terms, this is using x.L inside a constraint instead of x. I have seen this Pyomo mistake several times, and I think this is not very well explained and emphasized in the documentation. Pyomo constraints are generated before the solver gets involved, and some constructs are evaluated during that phase, instead of inside the solver. A difficult concept. A good way to see how the model passed on to the solver looks like is to print the scalar version using model.pprint():

2 Var Declarations a : Size=1, Index=None Key : Lower : Value : Upper : Fixed : Stale : Domain None : 5 : None : 20 : False : True : NonNegativeIntegers x : Size=1, Index=None Key : Lower : Value : Upper : Fixed : Stale : Domain None : 0 : 10 : None : False : False : NonNegativeIntegers 1 Objective Declarations total_production : Size=1, Index=None, Active=True Key : Active : Sense : Expression None : True : minimize : x*a 3 Constraint Declarations first_material_constraint : Size=1, Index=None, Active=True Key : Lower : Body : Upper : Active None : 70.0 : 0.2*a + 0.4*a + 0.6000000000000001*a + 0.8*a + a + 1.2000000000000002*a + 1.4000000000000001*a + 1.6*a + 1.8*a + 2.0*a : +Inf : True min_production : Size=1, Index=None, Active=True Key : Lower : Body : Upper : Active None : 1.0 : x : +Inf : True second_material_constraint : Size=1, Index=None, Active=True Key : Lower : Body : Upper : Active None : 90.0 : 0.8*a + 1.6*a + 2.4000000000000004*a + 3.2*a + 4.0*a + 4.800000000000001*a + 5.6000000000000005*a + 6.4*a + 7.2*a + 8.0*a : +Inf : True 6 Declarations: x a total_production first_material_constraint second_material_constraint min_production

We can see that these summations in the constraints are based on the initial value \(\color{darkred}x^{\mathrm{init}} = 10\): ten terms are being generated. We also see that the \(\color{darkred}x\) variable itself has disappeared. Note that range(p,q)in Python translates to something like \(p,\dots,q-1\).

Original Pyomo Model |
---|

\[\begin{align} \min\> & \color{darkred}a\cdot \color{darkred}x \\ & 0.2\cdot \color{darkred}a\cdot (1+\dots+\color{darkred}x) \ge 70\\ & 0.8\cdot \color{darkred}a\cdot (1+\dots+\color{darkred}x) \ge 90 \\ & \color{darkred}a \in \{5,\dots,20\} \\ & \color{darkred}x \in \{1,2,\dots\} \end{align}\] |

The constraints look a bit strange. I have never seen production limits like that. It says the \(i^\mathrm{th}\) batch has a yield proportional to \(i\). I really suspect this model is just not properly formulated. The logic being buried in Python code may have contributed to this. But anyway, we can substantially simplify this:

Cleaned-up Model |
---|

\[\begin{align} \min\> & \color{darkred}a\cdot \color{darkred}x \\ & \color{darkred}a\cdot \color{darkred}x \cdot (\color{darkred}x+1) \ge 2\cdot \max\{70/0.2,90/0.8\}\\ & \color{darkred}a \in \{5,\dots,20\} \\ & \color{darkred}x \in \{1,2,\dots\} \end{align}\] |

variablez 'objective';

integer variables

x 'Number of batches'a 'Batch size'

;

x.lo = 1;

a.lo = 5;

a.up = 20;

Equations

obj 'objective'limit 'minimum production'

;

obj.. z =e= a*x;

limit.. a*x*(x+1) =g= 2*

max(70/0.2,90/0.8);

modelm /all/;

optionminlp=baron;

solvemminimizingzusingminlp;

When dealing with a real model, I would store the numbers in the limit equation in parameters. That would give them a (meaningful) name.

LOWER LEVEL UPPER MARGINAL ---- EQU obj . . . 1.0000 ---- EQU limit 700.0000 780.0000 +INF . obj objective limit minimum production LOWER LEVEL UPPER MARGINAL ---- VAR z -INF 60.0000 +INF . ---- VAR x 1.0000 12.0000 +INF 5.0000 ---- VAR a 5.0000 5.0000 20.0000 12.0000 z objective x Number of batches a Batch size

### Conclusion

### References

- Pyomo constraints to enforce integer variables, https://stackoverflow.com/questions/77958306/pyomo-constraints-to-enforce-integer-variables

## No comments:

## Post a Comment