python - Gekko:获得的解决方案有问题

Question

我正在尝试解决 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-',LineWidth=2,label=r'$K$')

plt.legend(loc='best')
plt.xlabel('Time')
plt.ylabel('Value')
plt.show()

c和k的求解路径如下图:

这显然是不对的，因为 c 应该会随着时间的推移而增加，因为只关注之前给出的解决方案。

我不确定我哪里错了。

Answer 1

当前编写的最优控制问题是无界的。 c 的值将趋于无穷大以最大化函数。我将上限设置为 100在 c解算器去了那个界限。我重新制定了模型以反射(reflect)当前的问题陈述。这里有一些建议:

1. 使用m.integral()使模型更具可读性的函数。
2. 初始化c值不是 0 (默认)。您可能还想使用 c>0.01 设置下限这样m.log(c)如果求解器尝试一个值 <0 则不是未定义的.
3. 仅在 Gekko 方程中使用 Gekko 函数，例如 factor = m.exp(rate) .使用 factor = np.exp(rate)相反，除非它在可以对其进行评估的 Gekko 方程式中。
4. 附上精确解图，以便您可以比较精确解和数值解。
5. 使用m.options.NODES=3与 c=m.MV()和 c.STATUS=1以提高求解精度。默认值为 m.options.NODES=2那不那么准确。
6. 您可以使用 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--',linewidth=2,label=r'$C_{exact}$')
plt.ylabel('Value'); plt.legend(loc='best')

plt.subplot(3,1,2)
plt.plot(m.time,k.value,'b-',linewidth=2,label=r'$K$')
plt.legend(loc='best')

plt.subplot(3,1,3)
plt.plot(m.time,objective.value,'g:',linewidth=2,label=r'$obj$')
plt.plot(m.time,iobj.value,'k',linewidth=2,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,10,nt)
r = 0.33
i = 0.5

# Variables
c = m.MV(4,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,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--',linewidth=2,label=r'$C_{exact}$')
plt.ylabel('Value'); plt.legend(loc='best')

plt.subplot(3,1,2)
plt.plot(m.time,k.value,'b-',linewidth=2,label=r'$K$')
plt.legend(loc='best')

plt.subplot(3,1,3)
plt.plot(m.time,objective.value,'g:',linewidth=2,label=r'$obj$')
plt.plot(m.time,iobj.value,'k',linewidth=2,label=r'$\int obj$')
plt.legend(loc='best')
plt.xlabel('Time')
plt.show()