python - 如何修复 Python GEKKO 最优控制代码中的“未找到解决方案”错误
问题描述
我试图重现 K. Renee Fister 和 Jennifer Hughes Donnelly 于 2005 年撰写的论文“免疫疗法:最优控制理论方法”的图 1 中的结果。为此,我使用 Python 的 GEKKO 编写了一个数值最优控制求解器包裹。我使用了与论文中相同的初始条件、控制界限、参数值和模型方程。但是,当我运行代码时,出现以下错误:
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
我希望程序的输出提供两个数字:一个 ODE 动力学和一个最优控制解决方案的图。
我尝试过以多种方式更改代码:修改目标泛函、时间步数和更改最佳控制模式,但是,每次都得到相同的错误。下面是我正在使用的代码:
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()
此代码是通过修改此 Youtube 视频中提供的示例 GEKKO 代码得出的。任何解决此问题的帮助将不胜感激!
解决方案
当求解器找不到解决方案并报告“未找到解决方案”时,有一种故障排除方法来诊断问题。首先要做的是查看求解器输出m.solve(disp=True)
。求解器可能已经确定了一个不可行的问题,或者它达到了最大迭代次数而没有收敛到一个解决方案。
不可行的问题
如果求解器由于不可行的方程而失败,那么它会发现变量和方程的组合是不可解的。您可以尝试放宽变量界限或使用infeasibilities.txt
运行目录中的文件确定哪个方程不可行。从本地运行目录中检索infeasibilities.txt
文件,您可以使用m.open_folder()
when来查看该文件m=GEKKO(remote=False)
。
最大迭代限制
如果求解器达到默认迭代限制 ( m.options.MAX_ITER=250
),那么您可以尝试增加此限制,或者尝试以下策略。
- 通过设置
m.options.SOLVER=1
APOPT、m.options.SOLVER=2
BPOPT、m.options.SOLVER=3
IPOPT 或m.options.SOLVER=0
尝试所有可用的求解器来尝试不同的求解器。 - 首先通过求解变量数等于方程数的平方问题找到可行解。Gekko 有几个选项可以帮助解决这个问题,包括
m.options.COLDSTART=1
(为所有 FV 和 MV 设置 STATUS=0)或m.options.COLDSTART=2
(设置 STATUS=0 并执行块对角三角分解以找到可能的不可行方程)。 - 一旦找到可行的解决方案,尝试使用该解决方案作为初始猜测进行优化。
我将使用您在 YouTube 视频中引用的Luus 示例最优控制问题来演示此策略。这个问题在没有任何这些策略的情况下成功解决,但我将对其进行修改以展示如何使用这些方法。
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()
假设求解器不成功。您可以通过替换m.solve()
为以下内容进行初始化:
# 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()
初始化策略在文章中也有更详细的描述:
- Safdarnejad, SM, Hedengren, JD, Lewis, NR, Haseltine, E.,动态系统优化的初始化策略,计算机和化学工程,2015 年,第一卷。78,第 39-50 页,DOI:10.1016/j.compchemeng.2015.04.016。文章链接
推荐阅读
- ruby-on-rails - 将命令行参数读入 rb 文件
- python-3.x - Whoosh:从索引中检索文档编号和标题
- python - Tkinter get() 函数没有返回任何内容
- scala - Hikari 连接池无法从 sql 错误中恢复:当前事务被中止,命令被忽略,直到事务块结束
- r - 如何从 r 和 n 反算 P 值
- cart - BigCommerce - 在产品页面上显示“立即购买”按钮,购物车自动清空
- python - 在 python 中使用 str.format() 女巫类
- php - php dns_get_record ("something.com", DNS_TXT) 返回 false
- hy - 从 for 循环生成代码的宏
- html - Ionic 4 - 如何将模态背景更改为透明?