Gekko solver for optimizing function; Certain local minimum, but it will be found just on rare occasions

Question

from gekko import GEKKO 
import numpy as np
import matplotlib.pyplot as plt

m = GEKKO()  # initialize gekko

nt=200
m.time = np.linspace(0,1,nt)
# Variables
x1 = m.Var(value=20,lb=0,)
x2 = m.Var(value=0,lb=0,)
x3 = m.Var(value=0,lb=0,)

p = np.zeros(nt)  # mark final time point
p[-1] = 1.0

# optimize final time
tf = m.FV(value=0.7,lb=0.001,ub=100.0)
tf.STATUS = 1

# Equations
m.Equation(x1.dt() == -5*x1*tf)
m.Equation(x2.dt() == (5* x1 - 0.5*x2)*tf)
m.Equation(x3.dt() == 0.5*x2*tf)


#Minifunc
#m.Minimize(tf)
m.Maximize(x2)

#options
m.options.IMODE = 6
m.options.SOLVER = 3 #IPOPT optimizer
m.options.RTOL = 1E-8
m.options.OTOL = 1E-8
m.options.NODES = 5
m.solve(disp=True,debug=1)

print('Final Time: ' + str(tf.value[0]))
#Plotting stuff 
tm = np.linspace(0,tf.value[0],nt)

plt.figure(1)  # plot results
plt.plot(tm, x1.value, "k-", label=r"$x_1$")
plt.plot(tm, x2.value, "b-", label=r"$x_2$")
plt.plot(tm, x3.value, "r--", label=r"$x_3$")

plt.legend(loc="best")
plt.xlabel("Time")
plt.ylabel("Value")
plt.show()

I am trying to get used to Gekko as an optimization tool for some kinetic data. Therefore, I tried to implement a simple ode system based on an example of the Gekko website. Nevertheless, it never found the global maximum unless I force the derivative to be dx2/dt>0 to get the right maximum. This is somewhat unsatisfying, as I would expect a solver to find the maximum of component x2 just by asking to find the Maximum via m.Maximize(x2) for such a simple system. Is there any setting I did wrong or is my solution badly proposed? Thanks for helping.

The result: shows clearly that the maximum is not found...

Answer 1

The unconstrained objective (without the constraint dx2dt>=0) is 10416.8 while the objective with the constraint is 9236.4 (worse). Additional constraints never improve the optimization result but can make the objective worse for convex optimization problems such as this one.

m.Equation(x2.dt() >= 0)

Here is the code that shows the results of adding the constraint:

import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt

m = GEKKO()  # initialize gekko

nt=200
m.time = np.linspace(0,1,nt)
# Variables
x1 = m.Var(value=20,lb=0,)
x2 = m.Var(value=0,lb=0,)
x3 = m.Var(value=0,lb=0,)

p = np.zeros(nt)  # mark final time point
p[-1] = 1.0

# optimize final time
tf = m.FV(value=0.7,lb=0.001,ub=100.0)
tf.STATUS = 1

# Equations
m.Equation(x1.dt() == -5*x1*tf)
m.Equation(x2.dt() == (5* x1 - 0.5*x2)*tf)
m.Equation(x3.dt() == 0.5*x2*tf)

#Minifunc
#m.Minimize(tf)
m.Maximize(x2)

#options
m.options.IMODE = 6
m.options.SOLVER = 1 #APOPT optimizer
m.options.RTOL = 1E-8
m.options.OTOL = 1E-8
m.options.NODES = 5
m.solve(disp=True,debug=1)

plt.figure(1)
tm = np.linspace(0,tf.value[0],nt)
plt.plot(tm, x2.value, "b:", label=r"$x_2$ Optimal")

# add constraint
m.Equation(x2.dt() >= 0)
m.options.TIME_SHIFT=0
m.solve(disp=True,debug=1)

tm = np.linspace(0,tf.value[0],nt)
plt.plot(tm, x2.value, "r--", label=r"$x_2$ Constrained")
plt.legend(loc="best"); plt.xlabel("Time"); plt.ylabel("Value")

print('Final Time: ' + str(tf.value[0]))

plt.figure(2)  # plot results
plt.plot(tm, x1.value, "k-", label=r"$x_1$")
plt.plot(tm, x2.value, "b-", label=r"$x_2$")
plt.plot(tm, x3.value, "r--", label=r"$x_3$")

plt.legend(loc="best"); plt.xlabel("Time"); plt.ylabel("Value")

plt.show()

The current objective statement is to maximize the value of x2 everywhere.

m.Maximize(x2)

Even though there is a maximum part-way through the simulation, the final objective is the summation of all the x2 values. Make the following changes to maximize only the final x2 point.

p = np.zeros(nt)  # mark final time point
p[-1] = 1.0
final = m.Param(p)

m.Maximize(x2*final)

The optimal solution is now just the final point for x2 that is 15.485.