Search code examples
nonlinear-optimizationgekko

Dyn.Opt problem solution doesn't converge on GEKKO on varying no. of time pts or values of maniplatd var or colc. nodes

I'm solving a Dynamic Optimization problem on gekko but the solution doesn't converge if I vary the number of time points or initial/lb/ub values for manipulated variable or change the collocation nodes. What could be the issue?

Just changed the values mentioned above, and on some combinations, the solution doesn't converge at all, and the following error message appears:

raise Exception(response)

Exception:  @error: Solution Not Found

Code:

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

m = GEKKO()     
nt = 101        # no. of time steps
m.time = np.linspace(0,147,nt)

k17 = 0.0508
k26 = 1
k27 = 0.0577
k28 = 0.0104
k32 = 2
k33 = 2
k37 = 0.0016
k38 = 0.0107
k47 = 0.006
k48 = 0.072
k53 = 2
k57 = 0.0201
k58 = 0.082
k61 = 2
k62 = 2
k63 = 2
k72 = 2
k77 = 0.0133
k78 = 0.011
k87 = 0.081
k88 = 0.0148
kd = 0.06
w = 0.043

s1 = m.Var(value=0,lb=0) 
s2 = m.Var(value=3,lb=0)
s3 = m.Var(value=0.01,lb=0)
s4 = m.Var(value=0.46,lb=0)
s5 = m.Var(value=0.27,lb=0)
p1 = m.Var(value=0.51,lb=0)
p2 = m.Var(value=0.33,lb=0)
p3 = m.Var(value=0.3,lb=0)
p4 = m.Var(value=0.34e-4,lb=0)
p5 = m.Var(value=2.01,lb=0)
p6 = m.Var(value=0.05,lb=0)
X = m.Var(value=0.09,lb=0)
Xd = m.Var(value=0.02,lb=0)
V = m.Var(value=5,lb=0,ub=10)         
u = m.MV(value=0.05,lb=0,ub=0.1) # manipulated variable
u.STATUS = 1
u.DCOST = 0

f1 = (8.443e-4)*X*s1/(8.989e5 + s1)
f2 = (2.481e6)*X*s1*s3/((6.495e4 + s1)*(7.076e2 + s3))
f3 = (3.968e5)*X*s1*s3/((3.723e4 + s1)*(2.782e3 + s3))
f4 = (1.09e2)*X*s3/(0.019+s3)
f5 = 7.283*X*s4/(1.92e3 + s4)
f6 = (3.337e5)*X*s2*s5/((2.719e4 + s2)*(4.488e4 + s5))
f7 = (3.977e3)*X*s2/(9.324e3 + s2)
f8 = (6.697e-6)*X*s2/(0.537+s2)
f9 = (3.261e4)*X*s2/(6.683e5 + s2)

P = np.zeros(nt)
P[-1] = 1.0
final = m.Param(value=P)

# Equations
m.Equation(V.dt()==u)     
m.Equation(s1.dt()==-f1-f2-f3-k17*f7+(2-s1)*u/V)
m.Equation(s2.dt()==-f6-k27*f7-k28*f8-f9-s2*u/V)
m.Equation(s3.dt()==-k32*f2-k33*f3-f4+f6-k37*f7-k38*f8+f9-s3*u/V)
m.Equation(s4.dt()==-f5+f6-k47*f7-k48*f8-s4*u/V)
m.Equation(s5.dt()==k53*f3+f5-f6-k57*f7-k58*f8-s5*u/V)
m.Equation(p1.dt()==k61*f1+k62*f2+k63*f3-p1*u/V)
m.Equation(p2.dt()==k72*f2-k77*f7-k78*f8-p2*u/V)
m.Equation(p3.dt()==f4-k87*f7-k88*f8-p3*u/V)
m.Equation(p4.dt()==f8-p4*u/V)
m.Equation(p5.dt()==f7-p5*u/V)
m.Equation(p6.dt()==f5+f9-p6*u/V)
m.Equation(X.dt()==w*X-kd*X*Xd-X*u/V)
m.Equation(Xd.dt()==kd*X*Xd-Xd*u/V)

m.Obj(-final*V*p4) 
m.options.IMODE = 6
m.options.NODES = 4
#m.options.COLDSTART=2
#m.options.MAX_ITER=1000
m.solve(disp=True) 

p4_ = np.multiply(p4.value,1000)

plt.figure(1)
plt.subplot(2,1,1)
plt.plot(m.time,u.value,'r-')
plt.subplot(2,1,2)
plt.plot(m.time,p4_,'b--')
Solution

  • There are several reasons that the solver fails to reach a solution with Solution Not Found. Here are some of the reasons and things to try:

    1. If the solver reports No Feasible Solution then the bounds are too restrictive. Widen the variable bounds until a solution is obtained.

    2. If the solver reaches the maximum iterations, try changing the solver with m.options.SOLVER=1 or letting the solver continue for more iterations with m.options.MAX_ITER=500.

    3. If the solver stops with the message Unbounded Solution then try setting bounds on the MV values.

    4. Try initialization with no degrees of freedom by setting {FV}.STATUS=0 and {MV}.STATUS=0). The simulation should have the same number of equations and variables as a square problem. Switching to m.options.IMODE=7 is a sequential simulation and can also help with the initialization. Setting m.options.COLDSTART=1 or m.options.COLDSTART=2 can help find an initial solution that helps the optimization. Additional initialization resources are available in the Dynamic Optimization course.

    Please post a copy of your code for more specific suggestions. There are several Dynamic Optimization Benchmark problems posted to the Dynamic Optimization Course.

    Response to Edit

    Try setting a smaller time scale until the problem converges. Here is a successful solution with a final time of 0.08.

    m.time = np.linspace(0,0.08,11)

    The solver does not find a solution after this time point. It reports an infeasible solution. One cause of this may be variable p4 that makes it infeasible because of the lower bound of zero.

    p4_solution

    Here is the diagnostic code with the plots:

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

m = GEKKO()
m.time = np.linspace(0,0.08,11)
nt = len(m.time)

k17 = 0.0508
k26 = 1
k27 = 0.0577
k28 = 0.0104
k32 = 2
k33 = 2
k37 = 0.0016
k38 = 0.0107
k47 = 0.006
k48 = 0.072
k53 = 2
k57 = 0.0201
k58 = 0.082
k61 = 2
k62 = 2
k63 = 2
k72 = 2
k77 = 0.0133
k78 = 0.011
k87 = 0.081
k88 = 0.0148
kd = 0.06
w = 0.043

s1 = m.Var(value=0,lb=0) 
s2 = m.Var(value=3,lb=0)
s3 = m.Var(value=0.01,lb=0)
s4 = m.Var(value=0.46,lb=0)
s5 = m.Var(value=0.27,lb=0)
p1 = m.Var(value=0.51,lb=0)
p2 = m.Var(value=0.33,lb=0)
p3 = m.Var(value=0.3,lb=0)
p4 = m.Var(value=0.34e-4,lb=0)
p5 = m.Var(value=2.01,lb=0)
p6 = m.Var(value=0.05,lb=0)
X = m.Var(value=0.09,lb=0)
Xd = m.Var(value=0.02,lb=0)
V = m.Var(value=5,lb=0.01)         
u = m.MV(value=0.05,lb=0,ub=0.1) # manipulated variable
u.STATUS = 1
u.DCOST = 0

f1 = m.Intermediate((8.443e-4)*X*s1/(8.989e5 + s1))
f2 = m.Intermediate((2.481e6)*X*s1*s3/((6.495e4 + s1)*(7.076e2 + s3)))
f3 = m.Intermediate((3.968e5)*X*s1*s3/((3.723e4 + s1)*(2.782e3 + s3)))
f4 = m.Intermediate((1.09e2)*X*s3/(0.019+s3))
f5 = m.Intermediate(7.283*X*s4/(1.92e3 + s4))
f6 = m.Intermediate((3.337e5)*X*s2*s5/((2.719e4 + s2)*(4.488e4 + s5)))
f7 = m.Intermediate((3.977e3)*X*s2/(9.324e3 + s2))
f8 = m.Intermediate((6.697e-6)*X*s2/(0.537+s2))
f9 = m.Intermediate((3.261e4)*X*s2/(6.683e5 + s2))

P = np.zeros(nt)
P[-1] = 1.0
final = m.Param(value=P)

# Equations
m.Equation(V.dt()==u)     
m.Equation(s1.dt()==-f1-f2-f3-k17*f7+(2-s1)*u/V)
m.Equation(s2.dt()==-f6-k27*f7-k28*f8-f9-s2*u/V)
m.Equation(s3.dt()==-k32*f2-k33*f3-f4+f6-k37*f7-k38*f8+f9-s3*u/V)
m.Equation(s4.dt()==-f5+f6-k47*f7-k48*f8-s4*u/V)
m.Equation(s5.dt()==k53*f3+f5-f6-k57*f7-k58*f8-s5*u/V)
m.Equation(p1.dt()==k61*f1+k62*f2+k63*f3-p1*u/V)
m.Equation(p2.dt()==k72*f2-k77*f7-k78*f8-p2*u/V)
m.Equation(p3.dt()==f4-k87*f7-k88*f8-p3*u/V)
m.Equation(p4.dt()==f8-p4*u/V)
m.Equation(p5.dt()==f7-p5*u/V)
m.Equation(p6.dt()==f5+f9-p6*u/V)
m.Equation(X.dt()==w*X-kd*X*Xd-X*u/V)
m.Equation(Xd.dt()==kd*X*Xd-Xd*u/V)

m.Maximize(final*V*p4) 
m.options.IMODE = 6
m.options.NODES = 3
m.options.SOLVER = 1
m.options.COLDSTART=0
#m.options.MAX_ITER=1000
m.solve(disp=True) 

p4_ = np.multiply(p4.value,1000)

plt.figure(1)
plt.subplot(2,1,1)
plt.plot(m.time,u.value,'r-')
plt.subplot(2,1,2)
plt.plot(m.time,p4_,'b--')

plt.figure(2)
plt.plot(m.time,s1,label='s1')
plt.plot(m.time,s2,label='s2')
plt.plot(m.time,s3,label='s3')
plt.plot(m.time,s4,label='s4')
plt.plot(m.time,s5,label='s5')
plt.legend()

plt.figure(3)
plt.plot(m.time,p1,label='p1')
plt.plot(m.time,p2,label='p2')
plt.plot(m.time,p3,label='p3')
plt.plot(m.time,p4,label='p4')
plt.plot(m.time,p5,label='p5')
plt.plot(m.time,p6,label='p6')
plt.legend()

plt.figure(4)
plt.subplot(3,1,1)
plt.plot(m.time,s1,label='s1')
plt.legend()
plt.subplot(3,1,2)
plt.plot(m.time,s3,label='s3')
plt.legend()
plt.subplot(3,1,3)
plt.plot(m.time,p4,label='p4')
plt.legend()

plt.show()