使用Gekko的最优控制问题：“未找到解决方案”

Question

我需要用目标函数最小化解来解决一个最优控制问题。我的参考是APM示例，可在以下位置获得：https://www.youtube.com/watch?v=y26X-BSf8KU&list=PLLBUgWXdTBDjxcpH9hRuq-bsm_ti2UvoB&index=12

我的目标函数是Z=>；积分exp（-r*t）*C（t）.dt，其中C（t）=Ii.ai+Ij.aj。Ii和Ij是操纵变量，ai、aj和r是常数。目标是min Z

该问题受到以下约束：dKi/dt=Ii–deltai.Ki；dKj/dt=Ij–deltaj.Kj；Ki=（D-B.Kj-E）/A.Ki和Kj是变量，初始值Ki0和Kj0是已知的。D、 A、B、deltai、deltaj和E是常数系数

我开发了一个Gekko/python脚本，如下所示。但是，优化脚本没有实现一个解决方案（“异常：@error:solution not Found”）。我试图更改初始条件和初始猜测，但没有实现解决方案

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()

Answer 1

解算器表明您的问题是无限的

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

MV应该有边界吗

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

此外，有了这些边界，APOPT解算器报告问题不可行。似乎有3个方程可以解2个变量：Ki和Kj

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

如果你去掉最后一个等式，你就得到了一个成功的解。简言之，问题可能是无界的、过度指定的，或者两者兼而有之。以下是一个成功解决问题的版本：

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()

我还使用了新的m.integral()函数。您可能需要验证您的方程式，也可以尝试不同的解算器，看看是否有一个解算器提供了更有意义的信息

 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
                          -