用经济成本函数求最优控制收敛

Question

我一直在使用gekko来优化一个生物反应器，以示例12 (https://apmonitor.com/wiki/index.php/Main/GekkoPythonOptimization)为基础。

我的模型稍微复杂一些，包括6个状态、7个状态和2个操纵变量。当我运行它的时间小值(t ~20)，模拟是能够收敛(虽然需要一个良好的时间分辨率(dt < 0.1)。但是，当我试图延长时间(例如，t= 30)时，它的失败与以下错误相当一致：

EXIT: Converged to a point of local infeasibility. Problem may be infeasible

我尝试了以下几点：

1. 使用不同的m.options.SOLVER = 1,2,3的求解器
2. 将m.options.MAX_ITER提高到10000
3. 将m.options.NODES降至1(低阶解似有助于收敛)
4. 通过在声明D = m.MV(value=0.1,lb=0.0,ub=0.1)中指定一个值，向MVs提供一个合理的初始猜测。从一些不同的帖子来看，这似乎是有帮助的。

我不太清楚该如何解决这个问题。对于一个简化的模型(3个状态、5个参数和2个MVs)，即使简化模型的速率常数是整个模型的子集，gekko仍然能够很好地优化它(虽然当我尝试大t时它失败了)。

我的代码如下：

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

#Parameters and IC
full_params = [0.027,2.12e-9,7.13e-3,168,168,0.035,1e-3]
full_x0 = [5e6,0.0,0.0,0.0,1.25e5,0.0]
mu,k1,k2,k3,k33,k4, f= full_params

#Initialize model
m = GEKKO()

#Time discretization
n_steps = 201
m.time = np.linspace(0,20,n_steps)

#Define MVs
D = m.MV(value=0.1,lb=0.0,ub=0.1)
D.STATUS = 1
D.DCOST = 0.0
Tin = m.MV(value=1e7,lb=0.0,ub=1e7)
Tin.STATUS = 1
Tin.DCOST = 0.0

#Define States
T = m.Var(value=full_x0[0])
Id = m.Var(value=full_x0[1])
Is = m.Var(value=full_x0[2])
Ic = m.Var(value=full_x0[3])
Vs = m.Var(value=full_x0[4])
Vd = m.Var(value=full_x0[5])

#Define equations
m.Equation(T.dt() == mu*T -k1*(Vs+Vd)*T + D*(Tin-T))
m.Equation(Id.dt() == k1*Vd*T -(k1*Vs -mu)*Id -D*Id)
m.Equation(Is.dt() == k1*Vs*T -(k1*Vd + k2)*Is -D*Is)
m.Equation(Ic.dt() == k1*(Vs*Id + Vd*Is) -k2*Ic -D*Ic)
m.Equation(Vs.dt() == k3*Is - (k1*(T+Id+Is+Ic) + k4 + D)*Vs)
m.Equation(Vd.dt() == k33*Ic + f*k3*Is - (k1*(T+Id+Is+Ic) + k4 + D)*Vd)

#Define objective function
J = m.Var(value=0) # objective (profit)
Jf = m.FV() # final objective
Jf.STATUS = 1
m.Connection(Jf,J,pos2="end")
m.Equation(J.dt() == D*(Vs + Vd))
m.Obj(-Jf)
m.options.IMODE = 6  # optimal control
m.options.NODES = 1  # collocation nodes
m.options.SOLVER = 3
m.options.MAX_ITER = 10000

#Solve
m.solve()

为了清楚起见，模型方程如下：

​

​

如果能提供任何帮助，我将不胜感激，例如，如何实现每个https://apmonitor.com/do/index.php/Main/ModelInitialization参数的缩放。谢谢!

Answer 1

尝试增加最后一次的值，直到求解者无法再找到解决方案，比如使用tf=28 (成功)。解的图显示，在解几乎无法收敛的时候，Tin被调整为零。我添加了几个额外的目标形式，它们无助于收敛(请参阅客观方法#1和#2)。J、Vs、Vd的值很高，但求解器并不是无法管理的。考虑缩放的一种方法是将单元(如从kg/day转换为kg/s )作为基础。Gekko通过初始条件自动缩放变量。

​

​

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

#Parameters and IC
full_params = [0.027,2.12e-9,7.13e-3,168,168,0.035,1e-3]
full_x0 = [5e6,0.0,0.0,0.0,1.25e5,0.0]
mu,k1,k2,k3,k33,k4, f= full_params

#Initialize model
m = GEKKO()

#Time discretization
tf = 28
n_steps = tf*10+1
m.time = np.linspace(0,tf,n_steps)

#Define MVs
D = m.MV(value=0.1,lb=0.0,ub=0.1)
D.STATUS = 1
D.DCOST = 0.0
Tin = m.MV(value=1e7,lb=0,ub=1e7)
Tin.STATUS = 1
Tin.DCOST = 0.0

#Define States
T = m.Var(value=full_x0[0])
Id = m.Var(value=full_x0[1])
Is = m.Var(value=full_x0[2])
Ic = m.Var(value=full_x0[3])
Vs = m.Var(value=full_x0[4])
Vd = m.Var(value=full_x0[5])

#Define equations
m.Equation(T.dt() == mu*T -k1*(Vs+Vd)*T + D*(Tin-T))
m.Equation(Id.dt() == k1*Vd*T -(k1*Vs -mu)*Id -D*Id)
m.Equation(Is.dt() == k1*Vs*T -(k1*Vd + k2)*Is -D*Is)
m.Equation(Ic.dt() == k1*(Vs*Id + Vd*Is) -k2*Ic -D*Ic)
m.Equation(Vs.dt() == k3*Is - (k1*(T+Id+Is+Ic) + k4 + D)*Vs)
m.Equation(Vd.dt() == k33*Ic + f*k3*Is - (k1*(T+Id+Is+Ic) + k4 + D)*Vd)

# Original Objective
if True:
    J = m.Var(value=0) # objective (profit)
    Jf = m.FV() # final objective
    Jf.STATUS = 1
    m.Connection(Jf,J,pos2="end")
    m.Equation(J.dt() == D*(Vs + Vd))
    m.Obj(-Jf)

# Objective Method 1
if False:
    p=np.zeros_like(m.time); p[-1]=1
    final = m.Param(p)
    J = m.Var(value=0) # objective (profit)
    m.Equation(J.dt() == D*(Vs + Vd))
    m.Maximize(J*final)

# Objective Method 2
if False:
    m.Maximize(D*(Vs + Vd))

m.options.IMODE = 6  # optimal control
m.options.NODES = 2  # collocation nodes
m.options.SOLVER = 3
m.options.MAX_ITER = 10000

#Solve
m.solve()

plt.figure(figsize=(10,8))
plt.subplot(3,1,1)
plt.plot(m.time,Tin.value,'r.-',label='Tin')
plt.legend(); plt.grid()
plt.subplot(3,1,2)
plt.semilogy(m.time,T.value,label='T')
plt.semilogy(m.time,Id.value,label='Id')
plt.semilogy(m.time,Is.value,label='Is')
plt.semilogy(m.time,Ic.value,label='Ic')
plt.legend(); plt.grid()
plt.subplot(3,1,3)
plt.semilogy(m.time,Vs.value,label='Vs')
plt.semilogy(m.time,Vd.value,label='Vd')
plt.semilogy(m.time,J.value,label='Objective')
plt.legend(); plt.grid()
plt.show()

问题中是否存在某种类型的约束，最终有利于减少？这可能是tf=30不可行的原因。获得可行解的另一种方法是用m.options.TIME_STEP=20求解，并在时间步骤20的先验解等于值的前提下解决问题。

#Solve
m.solve()
m.options.TIME_SHIFT=20
m.solve()

这样，解决方案就会及时前进，在各个部分进行优化。这种策略被用来优化一个高空长耐力(HALE)无人机，并被称为后退地平线控制。

马丁，R.A.，盖茨，N.，宁，A.，赫登仁，J.D.，“在站台约束下高空太阳能飞机飞行轨迹的动态优化”，“制导、控制和动力学杂志”，2018年，doi: 10.2514/1.G003737。