如何解决Gekko:获得的解决方案有问题
我正在尝试解决以下 GEKKO 中的最优控制问题:
我先验地知道 c
的路径是:
其中参数值为:r = 0.33,i = 0.5,K(0) = 10 和 T = 10。
我用 Python 编写了以下代码:
import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
m = GEKKO(remote=True)
nt = 101; m.time = np.linspace(0,10,nt)
r = 0.33
i = 0.5
# Variables
c = m.Var()
k = m.Var(value=10)
objective = m.Var()
rate = [-r*t/10 for t in range(0,101)]
factor = m.exp(rate)
p = np.zeros(nt)
p[-1] = 1.0
final = m.Param(value=p)
disc = m.Param(value=factor)
# Equations
m.Equation(k.dt() == i*k - c)
m.Equation(objective.dt() == disc*m.log(c))
# Objective Function
m.Maximize(final*objective)
m.options.IMODE = 6
m.solve()
plt.figure(1)
plt.plot(m.time,c.value,'k:',linewidth=2,label=r'$C$')
plt.plot(m.time,k.value,'b-',label=r'$K$')
plt.legend(loc='best')
plt.xlabel('Time')
plt.ylabel('Value')
plt.show()
这显然是不对的,因为 c
应该随着时间的推移而增加,而不仅仅是盯着事先给出的解决方案。
我不确定我错在哪里。
解决方法
当前编写的最优控制问题是无界的。 c
的值将趋于无穷以最大化函数。我在 100
上设置了 c
的上限,求解器就到了那个界限。我重新制定了模型以反映当前的问题陈述。以下是一些建议:
- 使用
m.integral()
函数使模型更具可读性。 - 将
c
初始化为0
以外的值(默认)。您可能还想使用c>0.01
设置下限,以便在求解器尝试使用值m.log(c)
时不会未定义<0
。 - 仅在 Gekko 方程中使用 Gekko 函数,例如
factor = m.exp(rate)
。使用factor = np.exp(rate)
代替,除非它在可以计算的 Gekko 方程中。 - 包括精确解的图,以便您可以比较精确解和数值解。
- 将
m.options.NODES=3
与c=m.MV()
和c.STATUS=1
结合使用可提高解的准确性。默认值m.options.NODES=2
不太准确。 - 您可以使用
m.free_initial(c)
释放初始条件来计算c
的初始值。
import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
m = GEKKO(remote=True)
nt = 101; m.time = np.linspace(0,10,nt)
r = 0.33
i = 0.5
# Variables
c = m.MV(4,lb=0.01,ub=100); c.STATUS=1
#m.free_initial(c)
k = m.Var(value=10)
objective = m.Var(0)
t = m.Param(m.time)
m.Equation(objective==m.exp(-r*t)*m.log(c))
# just to include on the plot
iobj = m.Intermediate(m.integral(objective))
p = np.zeros(nt)
p[-1] = 1.0
final = m.Param(value=p)
# Equations
m.Equation(k.dt() == i*k - c)
# Objective Function
m.Maximize(final*m.integral(objective))
m.options.IMODE = 6
m.solve()
plt.figure(1)
plt.subplot(3,1,1)
plt.plot(m.time,c.value,'k:',linewidth=2,label=r'$C_{gekko}$')
C_sol = r*10*np.exp((i-r)*m.time)/(1-np.exp(-r*10))
plt.plot(m.time,C_sol,'r--',label=r'$C_{exact}$')
plt.ylabel('Value'); plt.legend(loc='best')
plt.subplot(3,2)
plt.plot(m.time,k.value,'b-',label=r'$K$')
plt.legend(loc='best')
plt.subplot(3,3)
plt.plot(m.time,objective.value,'g:',label=r'$obj$')
plt.plot(m.time,iobj.value,'k',label=r'$\int obj$')
plt.legend(loc='best')
plt.xlabel('Time')
plt.show()
是否缺少此问题的其他信息?
编辑:添加了额外的约束 k>0
。
按照评论中的建议添加了额外的约束。由于最后一个 c
值似乎不会影响解决方案,因此最后的结果与精确解略有不同。
import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
m = GEKKO(remote=True)
nt = 101; m.time = np.linspace(0,lb=0.001,ub=100); c.STATUS=1; c.DCOST=1e-6
m.free_initial(c)
k = m.Var(value=10,lb=0)
objective = m.Var(0)
t = m.Param(m.time)
m.Equation(objective==m.exp(-r*t)*m.log(c))
# just to include on the plot
iobj = m.Intermediate(m.integral(objective))
p = np.zeros(nt)
p[-1] = 1.0
final = m.Param(value=p)
# Equations
m.Equation(k.dt() == i*k - c)
# Objective Function
m.Maximize(final*m.integral(objective))
m.options.IMODE = 6
m.options.NODES = 3
m.solve()
plt.figure(1)
plt.subplot(3,label=r'$\int obj$')
plt.legend(loc='best')
plt.xlabel('Time')
plt.show()
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。