Optimal control problem using Gekko: “Solution Not Found”

Question

I need to solve an optimal control problem with an objective function minimization solution. The reference for me is the APM example available at: https://www.youtube.com/watch?v=y26X-BSf8KU&list=PLLBUgWXdTBDjxcpH9hRuq-bsm_ti2UvoB&index=12.

My objective function is Z => integration exp(-r*t)*C(t).dt, where C(t) = Ii.ai + Ij.aj. Ii and Ij are the manipulated variables and ai, aj and r are constants. The goal is min Z.

The problem is subject to the following constraints: dKi/dt = Ii – deltai.Ki; dKj/dt = Ij – deltaj.Kj; and Ki = (D-B.Kj-E)/A. Ki and Kj are variables and the initial values Ki0 and Kj0 are known. D, A, B, deltai, deltaj and E are constants coefficients.

I have developed a Gekko / python script as follow. However, the optimization script not achieves a solution (“Exception: @error: Solution Not Found”. I have tried to change the initial conditions and the initial guess but I haven´t achieve the solution.

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

#Create Gekko model
m=GEKKO()

#Time points
nt=101
m.time=np.linspace(0,10,nt)
t = m.Param(value=m.time)

#Constants
A=-0.0000159 
B=-0.0000506 
E=0.614 
deltai=0.031 
deltaj=0.1 
ai=10 
aj=27632 
r=0.085 

Ii=m.MV(value=100)
Ij=m.MV(value=100)
Ii.STATUS=1
Ii.COST=0
Ij.STATUS=1
Ij.COST=0

Ki=m.Var(value=15.548)
Kj=m.Var(value=0.932)

m.Equation(Ki.dt()==Ii-deltai*Ki)
m.Equation(Kj.dt()==Ij-deltaj*Kj)
m.Equation(Ki==(0.577-B*Kj-E)/A)

#Objective
Z=m.Var()

#Final objective
Zf=m.FV()
Zf.STATUS=1
m.Connection(Zf, Z, pos2='end')
m.Equation(Z.dt() == ((m.exp(-r*t)))*(Ii*ai+Ij*aj))

m.Obj(Zf)

#Options
m.options.IMODE=6
m.options.NODES=3
m.options.SOLVER=3
m.solve()

John Hedengren John Hedengren · Accepted Answer · 2020-12-16T22:51:43

The solver suggests that your problem is unbounded.

EXIT: Iterates diverging; problem might be unbounded.
 
 An error occured.
 The error code is            4
 
 
 ---------------------------------------------------
 Solver         :  IPOPT (v3.12)
 Solution time  :    1.44470000000729      sec
 Objective      :  -4.152133137209157E+021
 Unsuccessful with error code            0
 ---------------------------------------------------
 
 Creating file: infeasibilities.txt
 Use command apm_get(server,app,'infeasibilities.txt') to retrieve file
 @error: Solution Not Found
Traceback (most recent call last):
  File "C:\Users\johnh\Desktop\test.py", line 53, in <module>
    m.solve()
  File "C:\Python37\lib\site-packages\gekko\gekko.py", line 2174, in solve
    raise Exception(response)
Exception:  @error: Solution Not Found

Should the MVs have bounds?

Ii=m.MV(value=100,lb=0,ub=100)
Ij=m.MV(value=100,lb=0,ub=100)

Also, with these bounds the APOPT solver reports that the problem is infeasible. It appears that there are 3 equations for the solution of 2 variables: Ki and Kj.

m.Equation(Ki.dt()==Ii-deltai*Ki)
m.Equation(Kj.dt()==Ij-deltaj*Kj)
m.Equation(Ki*A==(0.577-B*Kj-E))

If you remove the last equation then you get a successful solution. In short, the problem may be unbounded, over-specified, or both. Here is a version that solves successfully:

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

m=GEKKO() #Create Gekko model

nt=101 #Time points
m.time=np.linspace(0,10,nt)
t = m.Param(value=m.time)
final = np.zeros(nt); final[-1] = 1
f = m.Param(final)

#Constants
A=-0.0000159; B=-0.0000506; E=0.614 
deltai=0.031; deltaj=0.1 
ai=10; aj=27632; r=0.085 

Ii=m.MV(value=100,lb=0,ub=100)
Ij=m.MV(value=100,lb=0,ub=100)
Ii.STATUS=1; Ii.COST=0
Ij.STATUS=1; Ij.COST=0

Ki=m.Var(value=15.548)
Kj=m.Var(value=0.932)

e = m.Intermediate(m.exp(-r*t))
m.Equation(Ki.dt()==Ii-deltai*Ki)
m.Equation(Kj.dt()==Ij-deltaj*Kj)
#m.Equation(Ki*A==(0.577-B*Kj-E))

#Minimize objective at final point
m.Minimize(f*m.integral(e*(Ii*ai+Ij*aj)))

#Options
m.options.IMODE=6; m.options.NODES=2
m.options.SOLVER=1
m.solve()

I also used the new m.integral() function. You may need to verify your equations and you can also try different solvers to see if one gives you more meaningful information.

 Each time step contains
   Objects      :            0
   Constants    :            0
   Variables    :            7
   Intermediates:            1
   Connections  :            0
   Equations    :            5
   Residuals    :            4
 
 Number of state variables:           1600
 Number of total equations: -         1400
 Number of slack variables: -            0
 ---------------------------------------
 Degrees of freedom       :            200
 
 ----------------------------------------------
 Dynamic Control with APOPT Solver
 ----------------------------------------------
 
 Iter    Objective  Convergence
    0  2.55775E+01  2.61383E+00
    1  1.99997E-03  9.61349E-09
    2  2.00000E-03  1.82031E-10
    3  2.00000E-03  2.22045E-16
    4  2.00000E-03  2.22045E-16
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :   0.498800000001211      sec
 Objective      :   1.999999949475750E-003
 Successful solution
 ---------------------------------------------------

Optimal control problem using Gekko: “Solution Not Found”

1 Answers