如何用GEKKO模拟带扰动的换热器温度控制过程？

问题核心 #

你遇到的关键矛盾是：热平衡方程里包含( \frac{dT_{in}}{dt} )项，但如果把( T_{in} )定义为Param类型，既不能直接调用dt()方法，也会因错误代入导致( T_{out} )出现异常阶跃。原热平衡方程整理后应为：

( \frac{dT_{out}}{dt} = \frac{m_s \lambda_s}{m_w C_p} + \frac{dT_{in}}{dt} )

解决思路 #

把( T_{in} )定义为Gekko的Var类型（而非Param），这样就能支持微分操作，再通过约束实现阶跃扰动。同时搭建控制器，通过调整蒸汽流量( m_s )让( T_{out} )稳定在设定值。

完整代码实现 #

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

# 初始化模型
m = GEKKO(remote=False)

# 时间配置：总时长600s，步长1s
m.time = np.linspace(0, 600, 601)

# 物理参数
Cp = 4.186e3       # 水比热，J/(kg·K)
lamb_s = 2085e3    # 6.5bar蒸汽潜热，J/kg
T_set = 110 + 273.15  # 出水温度设定值，K

# 核心变量定义
T_out = m.Var(value=100+273.15)  # 出水温度，初始100℃
m_s = m.MV(value=0.01)           # 蒸汽流量（操纵变量）
m_s.STATUS = 1                   # 允许控制器调整
m_s.LOWER = 0                    # 流量下限
m_s.UPPER = 0.1                  # 流量上限
m_s.DCOST = 0.1                  # 抑制频繁调整的惩罚项

# 扰动变量1：水流量阶跃
m_w = m.Param(value=1.0)
m_w.value[200:] = 1.2  # t=200s时，m_w从1.0kg/s跳至1.2kg/s

# 扰动变量2：进水温度阶跃（定义为Var以支持微分）
T_in = m.Var(value=100+273.15)
# t=300s时，T_in从100℃阶跃至90℃（263.15K）
m.fix(T_in, pos=300, val=263.15)

# 热平衡方程（严格匹配推导公式）
m.Equation(T_out.dt() == (m_s * lamb_s)/(m_w * Cp) + T_in.dt())

# 控制目标：T_out跟踪设定值
m.Obj((T_out - T_set)**2)

# 求解配置
m.options.IMODE = 6  # MPC控制模式
m.options.CV_TYPE = 2  # 平方误差目标
m.solve(disp=True)

# 结果可视化
plt.figure(figsize=(10,8))

plt.subplot(3,1,1)
plt.plot(m.time, [t-273.15 for t in T_out.value], label='T_out (℃)')
plt.plot(m.time, np.full(len(m.time), 110), 'r--', label='设定值')
plt.ylabel('温度')
plt.legend()
plt.grid(True)

plt.subplot(3,1,2)
plt.plot(m.time, m_w.value, label='m_w (kg/s)')
plt.plot(m.time, [t-273.15 for t in T_in.value], label='T_in (℃)')
plt.ylabel('扰动变量')
plt.legend()
plt.grid(True)

plt.subplot(3,1,3)
plt.plot(m.time, m_s.value, label='m_s (kg/s)')
plt.ylabel('蒸汽流量')
plt.xlabel('时间 (s)')
plt.legend()
plt.grid(True)

plt.show()

关键细节说明 #

1. T_in的Var定义：只有Var类型支持dt()微分操作，完美匹配热平衡方程的需求。
2. 阶跃实现方式：用m.fix(T_in, pos=300, val=263.15)在指定时间点固定T_in的值，实现精准阶跃；也可以用方程约束替代：

   step = m.Param(value=np.where(m.time>=300, 1, 0))
   m.Equation(T_in == (100+273.15) - 10*step)
   

3. 控制器配置：通过将m_s设为MV（操纵变量），结合平方误差目标，实现MPC自动控制，抵消扰动并维持T_out稳定。

内容的提问来源于stack exchange，提问作者Andre Furtado