CODEWITHSUNDEEP
pythonnonlinear-optimizationgekkomixed-integer-programming
Hamiduddin Hamdan
Hamiduddin Hamdan

Reputation: 21

Is there any benefit of choosing to formulate constraints in a way or another in GEKKO?

I have an MINLP problem and let's say the continuous variable Q can only be 0 when the binary variable z is 0. Two ways to formulate this would be:

m.Equation(Q*(1-z) == 0) (1)

or

m.Equation(Q < z*10000) (2)

whereby 10000 would be the upper bound to the continuous variable Q. Does (1) or (2) have any benefits over the other?

I've used (1) in my model for heat exchanger network synthesis and I got a good solution pretty quickly. Using (2) takes around 100x longer and it gives worse solution than the one given using (1).

From what I can see, the MINLP relaxation to (1) would still require z to be exactly 1 for Q to take on non-zero values. Does this have any effect on how APOPT solve the problem?

Upvotes: 1

Views: 36

Answers (2)

John Hedengren
John Hedengren

Reputation: 14376

I agree with Erwin's response. One other suggestion is to try using Q*(1-z)<eps where eps is a small value such as 1e-3 to improve convergence. This helps the solver explore the solution space to find the optimal solution. This strategy is also commonly used in Mathematical Programs with Complementarity Constraints. Here is an example that demonstrates the approach.

Option 1: Nonlinear Constraint Q*(1-z)<eps

from gekko import GEKKO
m = GEKKO(remote=False)
eps = 1e-3
Q = m.Var(value=0, lb=0, ub=1000)
z = m.Var(value=0.5, lb=0, ub=1, integer=True)

# Objective: maximize Q
m.Maximize(Q)

# Nonlinear constraint
m.Equation(Q * (1-z) < eps)

# Solve the problem
m.options.SOLVER = 1  # APOPT solver
m.solve(disp=True)

print(f"Q = {Q.value[0]}, z = {z.value[0]}")

Solution summary

 Number of state variables:              3
 Number of total equations: -            1
 Number of slack variables: -            1
 ---------------------------------------
 Degrees of freedom       :              1
 
 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
Iter: 1 I:  0 Tm: 0.00 NLPi: 12 Dpth: 0 Lvs: 0 Obj: -1.00E+03 Gap: 0.00E+00
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :   8.199999982025474E-003 sec
 Objective      :   -1000.00000000000     
 Successful solution
 ---------------------------------------------------
 
Q = 1000.0, z = 1.0

Option 2: Big-M Formulation

from gekko import GEKKO
m = GEKKO(remote=False)
Q = m.Var(value=0, lb=0, ub=1000)
z = m.Var(value=0.5, lb=0, ub=1, integer=True)

m.Maximize(Q)

# Big-M constraint with M = 10000
m.Equation(Q <= 10000 * z)

# Solve the problem
m.options.SOLVER = 1  # APOPT solver
m.solve(disp=True)

print(f"Q = {Q.value[0]}, z = {z.value[0]}")

The solution requires 2 iterations.

 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
Iter: 1 I: 0 Tm: 0.00 NLPi: 3 Dpth: 0 Lvs: 2 Obj: -1.00E+03 Gap: NaN
--Integer Solution:  -1.00E+03 Lowest Leaf: -1.00E+03 Gap: 0.00E+00
Iter: 2 I: 0 Tm: 0.00 NLPi: 2 Dpth: 1 Lvs: 2 Obj: -1.00E+03 Gap: 0.00E+00
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :   9.599999990314245E-003 sec
 Objective      :   -1000.00000000000     
 Successful solution
 ---------------------------------------------------

Q = 1000.0, z = 1.0

Gekko automatically translates user-defined inequality constraints into an equality constraint plus a slack variable as required by some solvers.

The inequality constraint Q < 10000*z becomes 10000*z-Q + s ==0 and s>0.

 Number of state variables:              3
 Number of total equations: -            1
 Number of slack variables: -            1
 ---------------------------------------
 Degrees of freedom       :              1

For APOPT solver, use Q*(1-z)<eps unless a linear form is required. If a Mixed-Integer Linear Programming (MILP) solver is required, use a Big-M or indicator constraints.

Upvotes: 1

Erwin Kalvelagen
Erwin Kalvelagen

Reputation: 16782

Q(1-z) = 0 is non-linear and non-convex while Q <= 10000z is linear (and convex). The last one is much better as long as the big-M constant is small.

If you can't reduce the size of the big-M constant, consider using indicator constraints (using suitable modeling tools and solvers).

Upvotes: 1

Related Questions