我正在尝试在 PYTHON 上使用 GEKKO 来控制 CSTR 罐的液位,同时操纵入口流量 q。我使用 pid 控制器尝试了同样的问题,它成功了。然而,在 GEKKO 上,高度并未跟踪其设定点。一旦我做了双重测试:在 200 的流速下,高度达到 800,当我将流速降低到 2 时,高度约为 0。但是,当我将 GEKKO 中的高度设定点设置为 700 或 800 时,流速并没有止步于200,而是无限地不断增加。
下面是我的代码:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from gekko import GEKKO
# Steady State Initial Condition
u2_ss=100.0
h_ss=399.83948182
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]
c1=5.0
Ac=10.0
# initial conditions
h0=399.83948182
q0=100.0
m.h= m.CV(value=h0)
m.q=m.MV(value=q0,lb=0,ub=100000)
m.Equation(m.h.dt()==(m.q-c1*m.h**0.5)/Ac)
#MV tuning
m.q.STATUS = 1
m.q.FSTATUS = 0
m.q.DMAX = 100
m.q.DMAXHI = 100
m.q.DMAXLO = -100
#CV tuning
DT = 0.5 # deadband
m.h.STATUS = 1
m.h.FSTATUS = 1
m.h.TR_INIT = 1
m.h.TAU = 1.0
m.options.CV_TYPE = 1
m.options.IMODE = 6
m.options.SOLVER = 3
#%% define CSTR model
def cstr(x,t,u2,Ac):
q=u2
c1=5.0
Ac=10.0
# States (2):
# the height of the tank (m)
h=x[0]
# Parameters:
# Calculate height derivative
dhdt=(q-c1*h**0.5)/Ac
if x[0]>=1000 and dhdt>0:
dh1dt = 0
# Return xdot:
xdot = np.zeros(1)
xdot[0]= dhdt
return xdot
# Time Interval (min)
t = np.linspace(0,20,100)
# 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:50] = 500.0
hsp[100:150]=700.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+1],Ac))
# retrieve measurements
h[i+1]= y[-1][0]
# insert measurement
m.h.MEAS=h[i+1]
# solve MPC
m.solve(disp=True)
m.h.SPHI = hsp[i+1] + DT
m.h.SPLO = hsp[i+1] - DT
# 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)