GEKKO can't find solution for many if-statements in the problem

Viewed 48

I'm trying to solve an optimisation problem where one of constraints is a function which selects non-zero mass materials and calculates a property of selected materials. In general I have mix of materials and I'm trying to optimise that mix calculating different parameters. One of the parameters is sum of some materials property only for non-zero mass materials - basically many if statements. This is an example of such function.

from gekko import GEKKO
m = GEKKO()
m.options.SOLVER = 1

properties = [2, 4, 7, 3, 2, 9, 6, 2]
m1 = m.Var(0, lb=0, ub=1)
m2 = m.Var(0.5, lb=0, ub=1)
m3 = m.Var(0.1, lb=0, ub=1)
m4 = m.Var(0.2, lb=0, ub=1)
m5 = m.Var(0.4, lb=0, ub=1)
m6 = m.Var(0.1, lb=0, ub=1)
m7 = m.Var(0.3, lb=0, ub=1)
m8 = m.Var(0.1, lb=0, ub=1)

m.Equation(m1+m2+m3+m4 == 1.0)
m.Equation(m5 + m6 + m7 + m8 == 1.0)

vars = [m1, m2, m3, m4, m5, m6, m7, m8]

actual_score_materials = []
max_score_materials = []
for idx, material in enumerate(vars):
        is_nonzero = m.if3(-material, 1, 0)
        var = m.Intermediate(properties[idx] * is_nonzero)
        actual_score_materials.append(var)
        max_score_materials.append((properties[idx] + 10) * is_nonzero)

actual_score = m.sum(actual_score_materials)
max_score = m.sum(max_score_materials)

obj = m.Intermediate(actual_score / max_score / 0.8)

m.Minimize(obj)
m.solve()

The problem is that solver can't find any solution for such defined problem (I've tried set up many different starting points with no success). Is there a way in GEKKO to define such problems? Thanks

1 Answers

Try setting a lower bound on max_score.

max_score = m.Var(lb=0.01)
m.Equation(max_score == m.sum(max_score_materials))

Without this lower bound, the objective function is NaN, indicating that it is setting max_score at zero to maximize the objective. As a result, the optimizer gets stuck. At the optimal solution, max_score.value=[115].

from gekko import GEKKO
m = GEKKO()


properties = [2, 4, 7, 3, 2, 9, 6, 2]
m1 = m.Var(0, lb=0, ub=1)
m2 = m.Var(0.5, lb=0, ub=1)
m3 = m.Var(0.1, lb=0, ub=1)
m4 = m.Var(0.2, lb=0, ub=1)
m5 = m.Var(0.4, lb=0, ub=1)
m6 = m.Var(0.1, lb=0, ub=1)
m7 = m.Var(0.3, lb=0, ub=1)
m8 = m.Var(0.1, lb=0, ub=1)

m.Equation(m1+m2+m3+m4 == 1.0)
m.Equation(m5+m6+m7+m8 == 1.0)

vars = [m1, m2, m3, m4, m5, m6, m7, m8]

actual_score_materials = []
max_score_materials = []
for idx, material in enumerate(vars):
        is_nonzero = m.if3(-material, 1, 0)
        var = m.Intermediate(properties[idx] * is_nonzero)
        actual_score_materials.append(var)
        max_score_materials.append((properties[idx] + 10) * is_nonzero)

actual_score = m.sum(actual_score_materials)
max_score = m.Var(lb=0.01)
m.Equation(max_score == m.sum(max_score_materials))

obj = m.Intermediate(actual_score / max_score / 0.8)

m.Minimize(obj)

m.options.SOLVER = 1
m.solve()

Here is the solver output:

 --------- APM Model Size ------------
 Each time step contains
   Objects      :            2
   Constants    :            0
   Variables    :           59
   Intermediates:            9
   Connections  :           18
   Equations    :           53
   Residuals    :           44
 
 Number of state variables:             59
 Number of total equations: -           45
 Number of slack variables: -           16
 ---------------------------------------
 Degrees of freedom       :             -2
 
 * Warning: DOF <= 0
 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
Iter:     1 I:  0 Tm:      0.00 NLPi:    7 Dpth:    0 Lvs:    2 Obj:  3.83E-01 Gap:       NaN
--Integer Solution:   4.00E-01 Lowest Leaf:   3.83E-01 Gap:   1.72E-02
Iter:     2 I:  0 Tm:      0.00 NLPi:    3 Dpth:    1 Lvs:    1 Obj:  4.00E-01 Gap:  1.72E-02
--Integer Solution:   3.80E-01 Lowest Leaf:   3.80E-01 Gap:   0.00E+00
Iter:     3 I:  0 Tm:      0.00 NLPi:    3 Dpth:    1 Lvs:    1 Obj:  3.80E-01 Gap:  0.00E+00
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :   3.800000000046566E-002 sec
 Objective      :   0.380434782608696     
 Successful solution
 ---------------------------------------------------
Related