How to Fix "Solution Not Found" Error in Python GEKKO Optimal Control Code

Question

I'm trying to reproduce the result in Figure 1 of the paper "Immunotherapy: An Optimal Control Theory Approach" by K. Renee Fister and Jennifer Hughes Donnelly, 2005. To do this, I have written a numerical optimal control solver using Python's GEKKO package. I have used the same initial conditions, control bounds, parameter values, and model equations as in the paper. However, when I run the code, I get the following error:

Traceback (most recent call last):
  File "xxxxxx", line 45, in <module>
    m.solve(disp=False) #solve
  File "xxxx", line 783, in solve
    raise Exception(response)
Exception:  @error: Solution Not Found

I expect the output of the program to provide two figures: one of the ODE dynamics and a plot of the optimal control solution.

I have tried changing the code in many ways: modifying the objective functional, number of time steps, and changing the optimal control mode, however, I get the same error each time. Below is the code I am using:

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

m = GEKKO()
nt = 1010
m.time = np.linspace(0,350,nt)

# Variables
X = m.Var(value=1)
Y = m.Var(value=1)
Z = m.Var(value=1)
OF = m.Var(value=0)
v = m.Var(value=0,lb=0,ub=1) #Control is initially 0 with a lower bound of 0 and an upper bound of 1

p = np.zeros(nt) #mark final time point
p[-1] = 1.0 #all zeros except the end, which is 1
final = m.Param(value=p) #final depends on integration limits

#Parameter Values
c = .000085
mu2 = .03
p1 = .1245
a = 1
r2 = .18
mu3 = 10
p2 = 5
g1 = 20000000 #2e7
g2 = 100000 #1e5
g3 = 1000 #1e3
b = 1*10**(-9)
s2 = 100000000
B = 100000000

# Equations
m.Equation(X.dt() == c*Y-mu2*X+(p1*X*Z)/(g1+Z))
m.Equation(Y.dt() == r2*Y*(1-b*Y)-(a*X*Y)/(g2+Y))
m.Equation(Z.dt() == (p2*X*Y)/(g3+Y)-mu3*Z+v*s2)
m.Equation(OF.dt() == X-Y+Z-B*v)

m.Obj(-OF*final)

m.options.IMODE = 6 #optimal control mode
m.solve(disp=False) #solve

plt.figure(figsize=(4,3)) #plot results
plt.subplot(2,1,1)
plt.plot(m.time,X.value,'k-',label=r'$S$')
plt.plot(m.time,Y.value,'b-',label=r'$R$')
plt.plot(m.time,Z.value,'g-',label=r'$E$')
plt.plot(m.time,OF.value,'r-',label=r'$OF$')
plt.legend()
plt.ylabel('CV')
plt.subplot(2,1,2)
plt.plot(m.time,v.value,'g--',label=r'$v$')
plt.legend()
plt.xlabel('Time')
plt.ylabel('Value')
plt.show()

This code was derived by modifying the example GEKKO code that was provided in this Youtube video. Any help resolving this issue would be much appreciated!

Answer 1

When the solver fails to find a solution and reports "Solution Not Found", there is a troubleshooting method to diagnose the problem. The first thing to do is to look at the solver output with m.solve(disp=True). The solver may have identified either an infeasible problem or it reached the maximum number of iterations without converging to a solution.

Infeasible Problem

If the solver failed because of infeasible equations then it found that the combination of variables and equations is not solvable. You can try to relax the variable bounds or identify which equation is infeasible with the infeasibilities.txt file in the run directory. Retrieve the infeasibilities.txt file from the local run directory that you can view with m.open_folder() when m=GEKKO(remote=False).

Maximum Iteration Limit

If the solver reached the default iteration limit (m.options.MAX_ITER=250) then you can either try to increase this limit or else try the strategies below.

• Try a different solver by setting m.options.SOLVER=1 for APOPT, m.options.SOLVER=2 for BPOPT, m.options.SOLVER=3 for IPOPT, or m.options.SOLVER=0 to try all the available solvers.
• Find a feasible solution first by solving a square problem where the number of variables is equal to the number of equations. Gekko a couple options to help with this including m.options.COLDSTART=1 (sets STATUS=0 for all FVs and MVs) or m.options.COLDSTART=2 (sets STATUS=0 and performs block diagonal triangular decomposition to find possible infeasible equations).
• Once a feasible solution is found, try optimizing with this solution as the initial guess.

I'll use the Luus example optimal control problem that you referenced with the YouTube video to demonstrate this strategy. This problem solves successfully without any of these strategies but I'll modify it to show how to use these methods.

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
m = GEKKO(remote=False)
nt = 101; m.time = np.linspace(0,2,nt)
x = m.Var(value=1)
u = m.MV(value=0,lb=-1,ub=1); u.STATUS=1
p = np.zeros(nt); p[-1] = 1.0; final = m.Param(value=p)
m.Equation(x.dt()==u)
m.Minimize(m.integral(0.5*x**2)*final)
m.options.IMODE = 6; m.solve()

plt.figure(1)
plt.plot(m.time,x.value,'k-',lw=2,label='x')
plt.plot(m.time,u.value,'r--',lw=2,label='u')
plt.legend(); plt.xlabel('Time'); plt.ylabel('Value')
plt.show()

Suppose that the solver were not successful. You could initialize by replacing m.solve() with the following:

# initialize to get a feasible solution
m.options.COLDSTART=2
m.solve()

# optimize, preserving the initial conditions (TIME_SHIFT=0)
m.options.TIME_SHIFT=0
m.options.COLDSTART=0
m.solve()

The intialization strategies are also described in more detail in the article:

• Safdarnejad, S.M., Hedengren, J.D., Lewis, N.R., Haseltine, E., Initialization Strategies for Optimization of Dynamic Systems, Computers and Chemical Engineering, 2015, Vol. 78, pp. 39-50, DOI: 10.1016/j.compchemeng.2015.04.016. article link