基于Gekko的MPC建模：为操纵变量添加过程动态遇阻

问题背景 #

我正在用Gekko开发模型预测控制（MPC）系统，需要给非线性过程的操纵变量（MV）添加动态特性，后续还要用公式计算增益。当前在为MV添加动态特性时遇到困难。

原始代码 #

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt  
import math 
import pandas as pd
import numpy as np

m = GEKKO(remote=False)
#m.open_folder()
t= np.linspace(0,500,100)
m.time=t

m.options.CV_TYPE = 2 # squared error
m.options.IMODE = 6  #MPC
m.options.SOLVER = 2

# Parameters
Tm = np.ones(len(t))*260
Tm[60:] = 265

Tr = m.Param(value=Tm)
Ti=m.Param(value=30)
HT=m.Param(value=230)

# Manipulated variable
F = m.Param(value=0.11, lb=0, ub=0.65)

j = m.MV(value=1, lb=0, ub=40)
j.STATUS = 1  
j.DCOST = 0.1 
j.DMAX = 5   

mv1 = m.MV(value=20, lb=18, ub=22)
mv1.STATUS = 1  
mv1.DCOST = 0.1 
mv1.DMAX = 0.1  
To = 8.6576*mv1 + 112.66

#MI
sp_m = np.ones(len(t))*4
sp_m[30:] = 10
sp_m[50:]=20
O_sp = m.Param(value=sp_m)

# Controlled Variable
P2 = m.Var(value=10) 

m.Obj((O_sp-O2)**2)

P1 = m.CV(value=94)
P1.STATUS = 1 
P1.SP = 95    
P1.TR_INIT = 1 
P1.TAU = 30   
P1.FDELAY=2  
P1.WSP=1   

# Process model
m.Equation(P1==(((To-(Ti+4))/(((16.15-0.019*(To+Ti)/2))*mv1))-0.06*F)*100)
m.Equation(P2==m.exp(0.792*m.log(j)-3.09*m.log(HT)+13.87*m.log(Tr)+3.06*m.log(1+F)-59.9))

m.options.IMODE = 6 # MPC
m.solve(disp=False)

import json
with open(m.path+'//results.json') as f:
    results = json.load(f)

plt.figure(figsize=(10,15))

plt.subplot(5,1,1)
plt.plot(m.time,mv1,'b-',label='MV Optimized')
plt.legend(loc='best')
plt.ylabel('MV1')

plt.subplot(5,1,2)
plt.plot(m.time,results['v2.tr'],'k-',label='Setpoint')
plt.plot(m.time,P1,'r--',label='CV Response')
plt.ylim([93, 96])
plt.ylabel('Conversion(%)')
plt.legend(loc='best')

plt.subplot(5,1,3)
plt.plot(m.time,results['p6'],'b-',label='MV Optimized')
plt.legend(loc='best')
plt.ylabel('MV2')

plt.subplot(5,1,4)
plt.plot(m.time,results['p8'],'k-',label='Setpoint')
plt.plot(m.time,results['v1'],'r--',label='CV Response')
plt.ylabel('P1(mu)')

plt.subplot(5,1,5)
plt.plot(m.time,Tm,'k-',label='Reactor temp(°C)')
plt.ylabel('T (°C)')
plt.xlabel('Time (min)')
plt.legend(loc='best')
plt.show()

需要修改的目标方程（为MV j 添加动态特性）：

P2==m.exp(0.792*m.log(j)-3.09*m.log(HT)+13.87 \ 
          *m.log(Tr)+3.06*m.log(1+F)-59.9))

尝试的无效方法 #

选项1（拉普拉斯逆变换） #

U1 =0.11/(1.8*(s**2)+2.88*s+1)*m.exp(-s) #transfer function
u1 = inverse_laplace_transform(U1,s,t)
P2==m.exp(0.792*m.log(j***u1**)-3.09*m.log(HT)\
          +13.87*m.log(Tr)+3.06*m.log(1+F)-59.9))

选项2（一阶加纯滞后模型） #

m.Equation(ft==k*(1-m.exp(-(t-tau_d)/tau))
m.Equation(P2==m.exp(0.792*m.log(j***ft**)-3.09*m.log(HT)\
            +13.87*m.log(Tr)+3.06*m.log(1+F)-59.9))

报错信息 #

Result: Equation without an equality (=) or inequality (> , <)
(((-(10.1010101010101-p1)))/(p2))
(((-(15.15151515151515-p1)))/(p2))
STOPPING...

解决方案 #

核心问题分析 #

1. Gekko不支持直接使用拉普拉斯变换或符号变量（如s），必须用微分方程形式描述动态特性
2. 尝试的代码存在语法错误（多余的***、未定义变量、方程书写不完整）
3. 动态特性需要通过新增状态变量实现，不能直接嵌入原方程

正确实现步骤 #

以给MV j 添加二阶加纯滞后动态特性为例（对应尝试的传递函数 U1 = 0.11/(1.8s²+2.88s+1) * e^(-s)）：

1. 定义动态特性的状态变量

将传递函数转换为微分方程，用m.delay()实现纯滞后，用状态变量表示二阶系统：

# 动态特性参数
K = 0.11
τ₁ = 1.8
τ₂ = 2.88
θ = 1  # 延迟时间

# 实现MV j的纯滞后
j_delayed = m.Var(value=j.value[0])
# 计算延迟对应的时间步数：总延迟时间 / 时间步长
delay_steps = int(θ / (m.time[1] - m.time[0]))
m.delay(j, j_delayed, delay_steps)

# 定义二阶系统的状态变量
x1 = m.Var(value=0)
x2 = m.Var(value=0)
j_dyn = m.Var(value=K * 0)  # 动态处理后的j值，用于P2方程

# 二阶系统微分方程
m.Equation(x1.dt() == x2)
m.Equation(x2.dt() == (j_delayed - x1 - τ₂*x2)/τ₁)
m.Equation(j_dyn == K * x1)

2. 修改P2方程，使用动态处理后的j_dyn

m.Equation(P2 == m.exp(0.792*m.log(j_dyn) - 3.09*m.log(HT) + 13.87*m.log(Tr) + 3.06*m.log(1+F) - 59.9))

3. 一阶加纯滞后模型简化实现（对应选项2思路）

如果需要一阶动态特性，可使用以下代码：

# 一阶加纯滞后参数
K = 1.0  # 增益，根据实际情况设置
tau = 5.0  # 时间常数
tau_d = 2.0  # 延迟时间

# 实现纯滞后
j_delayed = m.Var(value=j.value[0])
delay_steps = int(tau_d / (m.time[1]-m.time[0]))
m.delay(j, j_delayed, delay_steps)

# 一阶系统状态变量
j_dyn = m.Var(value=K * j.value[0])
m.Equation(j_dyn.dt() == (K*j_delayed - j_dyn)/tau)

# 修改P2方程
m.Equation(P2 == m.exp(0.792*m.log(j_dyn) - 3.09*m.log(HT) + 13.87*m.log(Tr) + 3.06*m.log(1+F) - 59.9))

注意事项 #

• 所有参数（如K、tau等）必须提前定义，避免未定义变量错误
• 检查方程语法：确保每个m.Equation()内有合法的等式/不等式，括号配对正确
• 动态特性的初始值要与MV初始值匹配，避免仿真初期的跳变

内容的提问来源于stack exchange，提问作者Adilton Lopes da Silva