Search code examples
pythonmathematical-optimizationgekko

Variable with conditionals

I'm working in a problem of optimization with Gekko and I have a complex equation with sums and many variables, but it is not a problem. The main difficulty is to write a variable that involves conditionals in its calculation, suppose this variable is called "X" and is calculated as mentioned in EQUATION, how could I implement this?

Equation

I've read Conditional Statements with Gekko, but I'm still in doubt.

Solution

  • Logical constraints can be implemented with m.if3() with binary switching variables (or m.if2() MPCC form).

    # conditionals
c1 = m.if3(a/b-10,1,0) # check if a/b>10
c2 = m.if3(10-a*b,1,0) # check if a*b<=10

    These switching variables can then be used to calculate X.

    # output
# when c1=0 (a/b>10), X=0
# when c1=1 (a/b<10) and c2=1 (a*b>10), X=a-b
# when c1=1 (a/b<10) and c2=0 (a*b<10), X=1
X = m.Var()
m.Equation(X==c1*c2*(a-b) + c1*(1-c2)*(1))

    results

    Here is a complete script for this simple test case in simulation where the value of a is a range of values.

    from gekko import GEKKO
m = GEKKO(remote=False)
a = m.Param([0.5,1,2,3,4,6,8,12,15,20,30,40,50])
b = 2
d = m.Intermediate(a/b); e=m.Intermediate(a*b)
# conditionals
c1 = m.if3(a/b-10,1,0) # check if a/b>10
c2 = m.if3(10-a*b,1,0) # check if a*b<=10
# output
# when c1=0 (a/b>10), X=0
# when c1=1 (a/b<10) and c2=1 (a*b>10), X=a-b
# when c1=1 (a/b<10) and c2=0 (a*b<10), X=1
X = m.Var()
m.Equation(X==c1*c2*(a-b) + c1*(1-c2)*(1))

# solve
m.options.IMODE=2
m.solve(disp=True)

print(a.value)
print(X.value)

# generate plot of results
import matplotlib.pyplot as plt
plt.subplot(4,1,1)
plt.plot(a.value,X.value,'b-o',label='X')
plt.legend(); plt.grid()
plt.subplot(4,1,2)
plt.plot(a.value,d.value,'k-o',label='a/b')
plt.plot([0.5,50],[10,10],'k--')
plt.legend(); plt.grid()
plt.subplot(4,1,3)
plt.plot(a.value,e.value,'r-o',label='a*b')
plt.plot([0.5,50],[10,10],'r--')
plt.legend(); plt.grid()
plt.subplot(4,1,4)
plt.plot(a.value,c1.value,'k:x',label='c1 (a/b<10)')
plt.plot(a.value,c2.value,'r:x',label='c2 (a*b<10)')
plt.legend(); plt.grid(); plt.xlabel('a')
plt.tight_layout(); plt.savefig('results.png',dpi=300)
plt.show()

    This can now be used in optimization to maximize, minimize, or reach a target value for X such as where the value of a is adjusted to maximize X:

    from gekko import GEKKO
m = GEKKO(remote=False)
a = m.Var(lb=0,ub=100)
b = 2
d = m.Intermediate(a/b); e=m.Intermediate(a*b)
# conditionals
c1 = m.if3(a/b-10,1,0) # check if a/b>10
c2 = m.if3(10-a*b,1,0) # check if a*b<=10
# output
# when c1=0 (a/b>10), X=0
# when c1=1 (a/b<10) and c2=1 (a*b>10), X=a-b
# when c1=1 (a/b<10) and c2=0 (a*b<10), X=1
X = m.Var()
m.Equation(X==c1*c2*(a-b) + c1*(1-c2)*(1))
m.Maximize(X)

# solve
m.options.IMODE=3
m.solve(disp=True)

print(f'a={a.value[0]}')
print(f'X={X.value[0]}')

    This gives the correct solution of:

    a=20.0
X=18.0

    Note that > and >= are equivalent for numerical optimization in Gekko. The solution of a=20.0 is at the switching condition between X=a-b and X=0. The simulation result shows X=0 at a=20 while the optimization result shows X=18 at a=20. There are default solver tolerances (m.options.RTOL for equations and m.options.RTOL for objective) that allow solutions within 1e-6 of the correct answer.