Thursday, February 8, 2024

Small non-convex MINLP: Pyomo vs GAMS

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

As Pyomo models are often somewhat difficult to read due to the large amount of "noise" (all kinds of syntactic stuff causing a low signal-to-noise ratio), it is a good idea to look at a mathematical model. 

The first constraint in Pyomo can be read as: \[\sum_{i=1}^{\color{darkred}x} 0.2\cdot\color{darkred}a \cdot i \ge 70\] A slight rewrite is: \[0.2\cdot\color{darkred}a \cdot (1+\dots+\color{darkred}x) \ge 70\] Summations with a variable upper bound are not that easy to handle with modeling tools. Often, we use a series of binary variables. Below I use a different technique. 

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}\]

We use here \[1+\dots+n = \frac{n\cdot(n+1)}{2}\] We replaced the summation with a variable upper bound to something that is much easier to handle inside a constraint. This is now easy to write down in GAMS:

variable 'objective';

integer variables

 'Number of batches'

 'Batch size'



x.lo = 1;

a.lo = 5;

a.up = 20;



  obj     'objective'

  limit   'minimum production'



obj.. z =e= a*x;

limit.. a*x*(x+1) =g= 2*max(70/0.2,90/0.8);


model m /all/;

option minlp=baron;

solve m minimizing z using minlp;

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

We solve with a global MINLP solver as the problem is non-convex. This gives:

                           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

This is a tiny model, so it does not take any time.


The original Pyomo model was a bit of a mess. There is a bit of a dichotomy between the Python language and math. And between what Python sees and what the underlying solvers see. What makes sense in Python, may not make sense in a Pyomo constraint. We can fix this easily by taking a step back from Python and look at the math. Sometimes Python code is just not at the right abstraction level for proper reasoning. As a result, new users may be totally lost even if they are proficient Python programmers. 

Models, whether implemented in a modeling language or a programming language, have two intended audiences: (1) the machine, or the solver if you want and (2) human readers such as colleagues, successors, or yourself. The latter group is arguably more important. But often, a readable model also has better changes to be correct and faster to solve. Just because it is easier to reason about it. This makes it the more important that model logic is written in a way, so that it is structure-revealing. Hopefully, for a well-written model, we don't need to re-engineer the math, like we did in this example. 


No comments:

Post a Comment