如何在GEKKO中平滑MV函数以优化生物柴油浓度？

我正在研究一个类似《Model Reduction with Intermediates in GEKKO》的动态优化问题，目标是找到最优温度曲线以最大化终态生物柴油浓度（x4(tf)），并与多项式温度曲线对比。

我希望温度曲线平滑且可微，但仅通过T = m.MV(value=298,ub=338,lb=298)和T.STATUS = 1得到的结果不符合状态变量实际情况。我尝试调整T.DCOST、T.DMAX、MV_STEP_HOR等各类参数，仍未得到理想结果，甚至与示例效果相差甚远，请问是否有遗漏的调参参数？

我的代码 #

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

m = GEKKO()

tf = 100 # 最终时间
m.time = np.linspace(0, tf, tf) # 时间区间

# 温度变量
T = m.MV(value=298, ub=338, lb=298)
T.STATUS = 1
T.DCOST = 1e-5
T.DMAX = 0.7

# 状态变量
x1 = m.Var(value=0.3226)
x2 = m.Var(value=0)
x3 = m.Var(value=0)
x4 = m.Var(value=0)
x5 = m.Var(value=1.9356)
x6 = m.Var(value=0)

# 反应速率的Arrhenius方程参数
A1, b1 = 3.92e7, 6614.83
A2, b2 = 5.77e5, 4997.98
A3, b3 = 5.88e12, 9993.96
A4, b4 = 0.98e10, 7366.64
A5, b5 = 5.35e3, 3231.18
A6, b6 = 2.15e4, 4824.87
k1 = m.Intermediate(A1*m.exp(-b1/T))
k2 = m.Intermediate(A2*m.exp(-b2/T))
k3 = m.Intermediate(A3*m.exp(-b3/T))
k4 = m.Intermediate(A4*m.exp(-b4/T))
k5 = m.Intermediate(A5*m.exp(-b5/T))
k6 = m.Intermediate(A6*m.exp(-b6/T))

# 动态模型方程
m.Equation(x1.dt() == -k1*x1*x5 + k2*x2*x4)
m.Equation(x2.dt() == k1*x1*x5 - k2*x2*x4 - k3*x2*x5 + k4*x3*x4)
m.Equation(x3.dt() == k3*x2*x5 - k4*x3*x4 - k5*x3*x5 + k6*x6*x4)
m.Equation(x4.dt() == k1*x1*x5 - k2*x2*x4 + k3*x2*x5 - k4*x3*x4 + k5*x3*x5 - k6*x6*x4)
m.Equation(x5.dt() == -(k1*x1*x5 - k2*x2*x4 + k3*x2*x5 - k4*x3*x4 + k5*x3*x5 - k6*x6*x4))
m.Equation(x6.dt() == k5*x3*x5 - k6*x6*x4)

# 终态目标权重
p = np.zeros(tf)
p[-1] = 1.0
final = m.Param(value=p)

m.Maximize(x4*final)
m.options.IMODE = 6

m.solve(disp=False)

print(f'终态x4浓度: {x4.value[-1]:0.4f}')
print(f'终态温度: {T.value[-1]}')

# 绘制浓度与温度曲线
plt.figure(figsize=(7,4))
plt.subplot(2,1,1)
plt.plot(m.time, x4.value, label='生物柴油')
plt.xlabel('时间 (min)')
plt.ylabel('浓度 (mol/L)')
plt.legend(loc='best')
plt.grid(True)

plt.subplot(2,1,2)
plt.plot(m.time, T.value, label=f'GEKKO DCOST={T.DCOST}', c='r')
plt.xlabel('时间 (min)')
plt.ylabel('温度 (K)')
plt.legend(loc='best')
plt.grid(True)
plt.tight_layout()
plt.savefig('biodiesel.png',dpi=300)
plt.show()

运行结果 #

内容的提问来源于stack exchange，提问作者Luiz Miguel