延迟调用应该放在gekko代码中的什么地方？
where should the delay call be placed inside a gekko code?
我正在尝试使用 GEKKO MPC 在操纵入口流量的同时控制水箱的液位。我想将 GEKKO 控制器建模为 FOPDT。我得到了我需要的所有参数，但我想使用 delay 函数来计算时间延迟。我不确定这个函数的确切位置，因为当我把它放在代码中时它给我一个错误。当我删除它（即没有时间延迟）时，代码工作正常但我想更现实并延迟时间。这是附加的代码：
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from gekko import GEKKO
# Steady State Initial Condition
u2_ss=10.0
h_ss=50.0
x0 = np.empty(1)
x0[0]= h_ss
#%% GEKKO nonlinear MPC
m = GEKKO(remote=False)
m.time = [0,0.02,0.04,0.06,0.08,0.1,0.12,0.15,0.2]
Ac=30.0
# initial conditions
h0=50.0
q0=10.0
Kp=93.48425357240352
taup=1010.8757590561246
thetap= 3
m.q=m.MV(value=q0,lb=0,ub=100)
m.h= m.CV(value=h0)
m.delay(m.q,m.h,thetap)
m.Equation(taup * m.h.dt()==m.q*Kp -m.h)
#MV tuning
m.q.STATUS = 1
m.q.FSTATUS = 0
m.q.DMAX = 100
#CV tuning
m.h.STATUS = 1
m.h.FSTATUS = 1
m.h.TR_INIT = 2
m.h.TAU = 1.0
m.h.SP = 55.0
m.options.CV_TYPE = 2
m.options.IMODE = 6
m.options.SOLVER = 3
#%% define CSTR model
def cstr(x,t,u2,Ac):
q=u2
Ac=30.0
# States (2):
# the height of the tank (m)
h=x[0]
# Parameters:
# Calculate height derivative
dhdt=(q-5)/Ac
# Return xdot:
xdot = np.zeros(1)
xdot[0]= dhdt
return xdot
# Time Interval (min)
t = np.linspace(0,20,400)
# Store results for plotting
hsp=np.ones(len(t))*h_ss
h=np.ones(len(t))*h_ss
u2 = np.ones(len(t)) * u2_ss
# Set point steps
hsp[0:100] = 55.0
hsp[100:]=70.0
# Create plot
plt.figure(figsize=(10,7))
plt.ion()
plt.show()
# Simulate CSTR
for i in range(len(t)-1):
# simulate one time period (0.05 sec each loop)
ts = [t[i],t[i+1]]
y = odeint(cstr,x0,ts,args=(u2[i],Ac))
# retrieve measurements
h[i+1]= y[-1][0]
# insert measurement
m.h.MEAS=h[i+1]
m.h.SP=hsp[i+1]
# solve MPC
m.solve(disp=True)
# retrieve new q value
u2[i+1] = m.q.NEWVAL
# update initial conditions
x0[0]= h[i+1]
#%% Plot the results
plt.clf()
plt.subplot(2,1,1)
plt.plot(t[0:i],u2[0:i],'b--',linewidth=3)
plt.ylabel('inlet flow')
plt.subplot(2,1,2)
plt.plot(t[0:i],hsp[0:i],'g--',linewidth=3,label=r'$h_{sp}$')
plt.plot(t[0:i],h[0:i],'k.-',linewidth=3,label=r'$h_{meas}$')
plt.xlabel('time')
plt.ylabel('tank level')
plt.legend(loc='best')
plt.draw()
plt.pause(0.01)
您的模型的问题是微分方程和延迟模型试图根据 m.q 的值求解 m.h 的值。两个方程不能同时满足。 delay object requires that m.h is delayed version of m.q from 3 cycles ago. The differential equation requires the solution of the linear differential equation。他们不会为 m.h 生成相同的答案，因此这会导致解决方案不可行，因为求解器会正确报告。
m.delay(m.q,m.h,thetap)
m.Equation(taup * m.h.dt()==m.q*Kp -m.h)
相反，您应该创建一个新变量，例如 m.qd 作为 m.q 的延迟版本。 m.dq 就是微分方程的输入。
m.qd=m.Var()
m.delay(m.q,m.qd,thetap)
m.Equation(taup * m.h.dt()==m.qd*Kp -m.h)
其他与问题无关的问题
您的申请还有其他一些问题。
- 模拟器和控制器之间的时间同步不正确。您应该为模拟器和控制器使用相同的循环时间。我将仿真时间更改为
t = np.linspace(0,20,201)，循环时间为 0.1 分钟。
- 延迟模型要求控制器具有统一的时间间隔，因为它是离散模型。我将控制器时间间隔更改为
m.time = np.linspace(0,2,21) 或循环时间为 0.1 分钟。
- 模拟器（用 ODEINT 解决）没有输入延迟，所以控制器和模拟器之间存在模型不匹配。这仍然没问题，因为模型不匹配是一个现实的场景，但您需要意识到会有一些基于模拟器反馈的纠正控制操作。控制器能够将液位驱动到设定点，但由于模型不匹配和平方误差 objective.
，MV 中存在抖动
为了改善抖音，我切换到m.options.CV_TYPE=1，设置SPHI和SPLO死区，用m.options.TR_OPEN=50打开初始轨迹，并添加了一个使用 m.q.DCOST 移动压制。这些具有实现类似性能但没有阀门颤动的效果。
这里是源代码：
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from gekko import GEKKO
# Steady State Initial Condition
u2_ss=10.0
h_ss=50.0
x0 = np.empty(1)
x0[0]= h_ss
#%% GEKKO nonlinear MPC
m = GEKKO(remote=False)
m.time = np.linspace(0,2,21)
Ac=30.0
# initial conditions
h0=50.0
q0=10.0
Kp=93.48425357240352
taup=1010.8757590561246
thetap= 3
m.q=m.MV(value=q0,lb=0,ub=100)
m.qd=m.Var(value=q0)
m.h= m.CV(value=h0)
m.delay(m.q,m.qd,thetap)
m.Equation(taup * m.h.dt()==m.qd*Kp -m.h)
#MV tuning
m.q.STATUS = 1
m.q.FSTATUS = 0
m.q.DMAX = 100
m.q.DCOST = 1
#CV tuning
m.h.STATUS = 1
m.h.FSTATUS = 1
m.h.TR_INIT = 1
m.h.TR_OPEN = 50
m.h.TAU = 0.5
m.options.CV_TYPE = 1
m.options.IMODE = 6
m.options.SOLVER = 3
#%% define CSTR model
def cstr(x,t,u2,Ac):
q=u2
Ac=30.0
# States (2):
# the height of the tank (m)
h=x[0]
# Parameters:
# Calculate height derivative
dhdt=(q-5)/Ac
# Return xdot:
xdot = np.zeros(1)
xdot[0]= dhdt
return xdot
# Time Interval (min)
t = np.linspace(0,20,201)
# Store results for plotting
hsp=np.ones(len(t))*h_ss
h=np.ones(len(t))*h_ss
u2 = np.ones(len(t)) * u2_ss
# Set point steps
hsp[0:100] = 55.0
hsp[100:] = 70.0
# Create plot
plt.figure(figsize=(10,7))
plt.ion()
plt.show()
# Simulate CSTR
for i in range(len(t)-1):
# simulate one time period (0.05 sec each loop)
ts = [t[i],t[i+1]]
y = odeint(cstr,x0,ts,args=(u2[i],Ac))
# retrieve measurements
h[i+1]= y[-1][0]
# insert measurement
m.h.MEAS=h[i+1]
# for CV_TYPE = 1
m.h.SPHI=hsp[i+1]+0.05
m.h.SPLO=hsp[i+1]-0.05
# for CV_TYPE = 2
m.h.SP=hsp[i+1]
# solve MPC
m.solve(disp=False)
# retrieve new q value
u2[i+1] = m.q.NEWVAL
# update initial conditions
x0[0]= h[i+1]
#%% Plot the results
plt.clf()
plt.subplot(2,1,1)
plt.plot(t[0:i],u2[0:i],'b--',linewidth=3)
plt.ylabel('inlet flow')
plt.subplot(2,1,2)
plt.plot(t[0:i],hsp[0:i],'g--',linewidth=3,label=r'$h_{sp}$')
plt.plot(t[0:i],h[0:i],'k.-',linewidth=3,label=r'$h_{meas}$')
plt.xlabel('time')
plt.ylabel('tank level')
plt.legend(loc='best')
plt.draw()
plt.pause(0.01)
我正在尝试使用 GEKKO MPC 在操纵入口流量的同时控制水箱的液位。我想将 GEKKO 控制器建模为 FOPDT。我得到了我需要的所有参数，但我想使用 delay 函数来计算时间延迟。我不确定这个函数的确切位置，因为当我把它放在代码中时它给我一个错误。当我删除它（即没有时间延迟）时，代码工作正常但我想更现实并延迟时间。这是附加的代码：
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from gekko import GEKKO
# Steady State Initial Condition
u2_ss=10.0
h_ss=50.0
x0 = np.empty(1)
x0[0]= h_ss
#%% GEKKO nonlinear MPC
m = GEKKO(remote=False)
m.time = [0,0.02,0.04,0.06,0.08,0.1,0.12,0.15,0.2]
Ac=30.0
# initial conditions
h0=50.0
q0=10.0
Kp=93.48425357240352
taup=1010.8757590561246
thetap= 3
m.q=m.MV(value=q0,lb=0,ub=100)
m.h= m.CV(value=h0)
m.delay(m.q,m.h,thetap)
m.Equation(taup * m.h.dt()==m.q*Kp -m.h)
#MV tuning
m.q.STATUS = 1
m.q.FSTATUS = 0
m.q.DMAX = 100
#CV tuning
m.h.STATUS = 1
m.h.FSTATUS = 1
m.h.TR_INIT = 2
m.h.TAU = 1.0
m.h.SP = 55.0
m.options.CV_TYPE = 2
m.options.IMODE = 6
m.options.SOLVER = 3
#%% define CSTR model
def cstr(x,t,u2,Ac):
q=u2
Ac=30.0
# States (2):
# the height of the tank (m)
h=x[0]
# Parameters:
# Calculate height derivative
dhdt=(q-5)/Ac
# Return xdot:
xdot = np.zeros(1)
xdot[0]= dhdt
return xdot
# Time Interval (min)
t = np.linspace(0,20,400)
# Store results for plotting
hsp=np.ones(len(t))*h_ss
h=np.ones(len(t))*h_ss
u2 = np.ones(len(t)) * u2_ss
# Set point steps
hsp[0:100] = 55.0
hsp[100:]=70.0
# Create plot
plt.figure(figsize=(10,7))
plt.ion()
plt.show()
# Simulate CSTR
for i in range(len(t)-1):
# simulate one time period (0.05 sec each loop)
ts = [t[i],t[i+1]]
y = odeint(cstr,x0,ts,args=(u2[i],Ac))
# retrieve measurements
h[i+1]= y[-1][0]
# insert measurement
m.h.MEAS=h[i+1]
m.h.SP=hsp[i+1]
# solve MPC
m.solve(disp=True)
# retrieve new q value
u2[i+1] = m.q.NEWVAL
# update initial conditions
x0[0]= h[i+1]
#%% Plot the results
plt.clf()
plt.subplot(2,1,1)
plt.plot(t[0:i],u2[0:i],'b--',linewidth=3)
plt.ylabel('inlet flow')
plt.subplot(2,1,2)
plt.plot(t[0:i],hsp[0:i],'g--',linewidth=3,label=r'$h_{sp}$')
plt.plot(t[0:i],h[0:i],'k.-',linewidth=3,label=r'$h_{meas}$')
plt.xlabel('time')
plt.ylabel('tank level')
plt.legend(loc='best')
plt.draw()
plt.pause(0.01)
您的模型的问题是微分方程和延迟模型试图根据 m.q 的值求解 m.h 的值。两个方程不能同时满足。 delay object requires that m.h is delayed version of m.q from 3 cycles ago. The differential equation requires the solution of the linear differential equation。他们不会为 m.h 生成相同的答案，因此这会导致解决方案不可行，因为求解器会正确报告。
m.delay(m.q,m.h,thetap)
m.Equation(taup * m.h.dt()==m.q*Kp -m.h)
相反，您应该创建一个新变量，例如 m.qd 作为 m.q 的延迟版本。 m.dq 就是微分方程的输入。
m.qd=m.Var()
m.delay(m.q,m.qd,thetap)
m.Equation(taup * m.h.dt()==m.qd*Kp -m.h)
其他与问题无关的问题
您的申请还有其他一些问题。
- 模拟器和控制器之间的时间同步不正确。您应该为模拟器和控制器使用相同的循环时间。我将仿真时间更改为
t = np.linspace(0,20,201)，循环时间为 0.1 分钟。 - 延迟模型要求控制器具有统一的时间间隔，因为它是离散模型。我将控制器时间间隔更改为
m.time = np.linspace(0,2,21)或循环时间为 0.1 分钟。 - 模拟器（用 ODEINT 解决）没有输入延迟，所以控制器和模拟器之间存在模型不匹配。这仍然没问题，因为模型不匹配是一个现实的场景，但您需要意识到会有一些基于模拟器反馈的纠正控制操作。控制器能够将液位驱动到设定点，但由于模型不匹配和平方误差 objective. ，MV 中存在抖动
为了改善抖音，我切换到m.options.CV_TYPE=1，设置SPHI和SPLO死区，用m.options.TR_OPEN=50打开初始轨迹，并添加了一个使用 m.q.DCOST 移动压制。这些具有实现类似性能但没有阀门颤动的效果。
这里是源代码：
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from gekko import GEKKO
# Steady State Initial Condition
u2_ss=10.0
h_ss=50.0
x0 = np.empty(1)
x0[0]= h_ss
#%% GEKKO nonlinear MPC
m = GEKKO(remote=False)
m.time = np.linspace(0,2,21)
Ac=30.0
# initial conditions
h0=50.0
q0=10.0
Kp=93.48425357240352
taup=1010.8757590561246
thetap= 3
m.q=m.MV(value=q0,lb=0,ub=100)
m.qd=m.Var(value=q0)
m.h= m.CV(value=h0)
m.delay(m.q,m.qd,thetap)
m.Equation(taup * m.h.dt()==m.qd*Kp -m.h)
#MV tuning
m.q.STATUS = 1
m.q.FSTATUS = 0
m.q.DMAX = 100
m.q.DCOST = 1
#CV tuning
m.h.STATUS = 1
m.h.FSTATUS = 1
m.h.TR_INIT = 1
m.h.TR_OPEN = 50
m.h.TAU = 0.5
m.options.CV_TYPE = 1
m.options.IMODE = 6
m.options.SOLVER = 3
#%% define CSTR model
def cstr(x,t,u2,Ac):
q=u2
Ac=30.0
# States (2):
# the height of the tank (m)
h=x[0]
# Parameters:
# Calculate height derivative
dhdt=(q-5)/Ac
# Return xdot:
xdot = np.zeros(1)
xdot[0]= dhdt
return xdot
# Time Interval (min)
t = np.linspace(0,20,201)
# Store results for plotting
hsp=np.ones(len(t))*h_ss
h=np.ones(len(t))*h_ss
u2 = np.ones(len(t)) * u2_ss
# Set point steps
hsp[0:100] = 55.0
hsp[100:] = 70.0
# Create plot
plt.figure(figsize=(10,7))
plt.ion()
plt.show()
# Simulate CSTR
for i in range(len(t)-1):
# simulate one time period (0.05 sec each loop)
ts = [t[i],t[i+1]]
y = odeint(cstr,x0,ts,args=(u2[i],Ac))
# retrieve measurements
h[i+1]= y[-1][0]
# insert measurement
m.h.MEAS=h[i+1]
# for CV_TYPE = 1
m.h.SPHI=hsp[i+1]+0.05
m.h.SPLO=hsp[i+1]-0.05
# for CV_TYPE = 2
m.h.SP=hsp[i+1]
# solve MPC
m.solve(disp=False)
# retrieve new q value
u2[i+1] = m.q.NEWVAL
# update initial conditions
x0[0]= h[i+1]
#%% Plot the results
plt.clf()
plt.subplot(2,1,1)
plt.plot(t[0:i],u2[0:i],'b--',linewidth=3)
plt.ylabel('inlet flow')
plt.subplot(2,1,2)
plt.plot(t[0:i],hsp[0:i],'g--',linewidth=3,label=r'$h_{sp}$')
plt.plot(t[0:i],h[0:i],'k.-',linewidth=3,label=r'$h_{meas}$')
plt.xlabel('time')
plt.ylabel('tank level')
plt.legend(loc='best')
plt.draw()
plt.pause(0.01)