Прикрепленный код, который я написал: При запуске переменная, контролируемая уровнем, не отслеживает свою уставку.С другой стороны, контролируемая переменная температура очень хорошо отслеживает свою уставку.Я использую манипулирование температурой охлаждения и расходом на входе.Я пытаюсь контролировать уровень в баке, температуру и концентрацию.
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from gekko import GEKKO
# Steady State Initial Condition
u1_ss = 300.0
u2_ss=100.0
Ca_ss = 0.87725294
T_ss = 324.47544313
h_ss=75.82018806
# Feed Temperature (K)
Tf = 350
# Feed Concentration (mol/m^3)
Caf = 1
# Steady State Initial Conditions for the States
x0 = np.empty(2)
x0[0] = Ca_ss
x0[1] = T_ss
p0=np.empty(1)
p0[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=10
Ac=400.0
# Volume of CSTR (m^3)
V = 100
# Density of A-B Mixture (kg/m^3)
rho = 1000
# Heat capacity of A-B Mixture (J/kg-K)
Cp = 0.239
# Heat of reaction for A->B (J/mol)
mdelH = 5e4
# E - Activation energy in the Arrhenius Equation (J/mol)
# R - Universal Gas Constant = 8.31451 J/mol-K
EoverR = 8750
# Pre-exponential factor (1/sec)
k0 = 7.2e10
# U - Overall Heat Transfer Coefficient (W/m^2-K)
# A - Area - this value is specific for the U calculation (m^2)
UA = 5e4
# initial conditions
Tc0 = 300
T0 = 324.47544313
Ca0 = 0.87725294
h0=75.82018806
q0=100.0
tau = m.Const(value=0.5)
Kp = m.Const(value=1)
m.Tc = m.MV(value=Tc0,lb=250,ub=350)
m.T = m.CV(value=T_ss)
m.rA = m.Var(value=0)
m.Ca = m.CV(value=Ca_ss,lb=0,ub=1)
m.h=m.CV(value=h_ss)
m.q=m.MV(value=q0,lb=0,ub=1000)
m.Equation(m.rA == k0*m.exp(-EoverR/m.T)*m.Ca)
m.Equation(m.T.dt() == m.q/V*(Tf - m.T) \
+ mdelH/(rho*Cp)*m.rA \
+ UA/V/rho/Cp*(m.Tc-m.T))
m.Equation(m.Ca.dt() == (m.q)/V*(Caf - m.Ca) - m.rA)
m.Equation(m.h.dt()==(m.q-c1*pow(m.h,0.5))/Ac)
#MV tuning
m.Tc.STATUS = 1
m.Tc.FSTATUS = 0
m.Tc.DMAX = 100
m.Tc.DMAXHI = 20
m.Tc.DMAXLO = -100
m.q.STATUS = 1
m.q.FSTATUS = 0
m.q.DMAX = 10
#CV tuning
m.T.STATUS = 1
m.T.FSTATUS = 1
m.T.TR_INIT = 1
m.T.TAU = 1.0
DT = 0.5 # deadband
m.h.STATUS = 1
m.h.FSTATUS = 1
m.h.TR_INIT = 1
m.h.TAU = 1.0
m.Ca.STATUS = 1
m.Ca.FSTATUS = 0 # no measurement
m.Ca.TR_INIT = 0
m.options.CV_TYPE = 1
m.options.IMODE = 6
m.options.SOLVER = 3
# define CSTR model
def cstr(x,t,u1,u2,Tf,Caf):
# Inputs (3):
# Temperature of cooling jacket (K)
Tc = u1
q=u2
# Tf = Feed Temperature (K)
# Caf = Feed Concentration (mol/m^3)
# States (2):
# Concentration of A in CSTR (mol/m^3)
Ca = x[0]
# Temperature in CSTR (K)
T = x[1]
# Parameters:
# Volume of CSTR (m^3)
V = 100
# Density of A-B Mixture (kg/m^3)
rho = 1000
# Heat capacity of A-B Mixture (J/kg-K)
Cp = 0.239
# Heat of reaction for A->B (J/mol)
mdelH = 5e4
# E - Activation energy in the Arrhenius Equation (J/mol)
# R - Universal Gas Constant = 8.31451 J/mol-K
EoverR = 8750
# Pre-exponential factor (1/sec)
k0 = 7.2e10
# U - Overall Heat Transfer Coefficient (W/m^2-K)
# A - Area - this value is specific for the U calculation (m^2)
UA = 5e4
# reaction rate
rA = k0*np.exp(-EoverR/T)*Ca
# Calculate concentration derivative
dCadt = q/V*(Caf - Ca) - rA
# Calculate temperature derivative
dTdt = q/V*(Tf - T) \
+ mdelH/(rho*Cp)*rA \
+ UA/V/rho/Cp*(Tc-T)
# Return xdot:
xdot = np.zeros(2)
xdot[0] = dCadt
xdot[1] = dTdt
return xdot
def tank(p,t,u2,Ac):
q=u2
h=p[0]
dhdt=(q-c1*pow(h,0.5))/Ac
if p[0]>=300 and dhdt>0:
dhdt = 0
return dhdt
# Time Interval (min)
t = np.linspace(0,10,410)
# Store results for plotting
Ca = np.ones(len(t)) * Ca_ss
T = np.ones(len(t)) * T_ss
Tsp=np.ones(len(t))*T_ss
hsp=np.ones(len(t))*h_ss
h=np.ones(len(t))*h_ss
u1 = np.ones(len(t)) * u1_ss
u2 = np.ones(len(t)) * u2_ss
# Set point steps
Tsp[0:100] = 330.0
Tsp[100:200] = 350.0
hsp[200:300] = 150.0
hsp[300:] = 190.0
# Create plot
plt.figure(figsize=(10,7))
plt.ion()
plt.show()
# Simulate CSTR
for i in range(len(t)-1):
ts = [t[i],t[i+1]]
y = odeint(cstr,x0,ts,args=(u1[i+1],u2[i+1],Tf,Caf))
y1=odeint(tank,p0,ts,args=(u2[i+1],Ac))
Ca[i+1] = y[-1][0]
T[i+1] = y[-1][1]
h[i+1]=y1[-1][0]
# insert measurement
m.T.MEAS = T[i+1]
m.h.MEAS= h[i+1]
# solve MPC
m.solve(disp=True)
m.T.SPHI = Tsp[i+1] + DT
m.T.SPLO = Tsp[i+1] - DT
m.h.SPHI = hsp[i+1] + DT
m.h.SPLO = hsp[i+1] - DT
# retrieve new Tc value
u1[i+1] = m.Tc.NEWVAL
u2[i+1]= m.q.NEWVAL
# update initial conditions
x0[0] = Ca[i+1]
x0[1] = T[i+1]
p0[0]=h[i+1]
plt.clf()
# Plot the results
plt.subplot(5,1,1)
plt.plot(t[0:i],u1[0:i],'b--',linewidth=3)
plt.ylabel('Cooling T (K)')
plt.legend(['Jacket Temperature'],loc='best')
plt.subplot(5,1,2)
plt.plot(t[0:i],u2[0:i],'g--')
plt.xlabel('time')
plt.ylabel('flow in')
plt.subplot(5,1,3)
plt.plot(t[0:i],Ca[0:i],'r-',linewidth=3)
plt.ylabel('Ca (mol/L)')
plt.legend(['Reactor Concentration'],loc='best')
plt.subplot(5,1,4)
plt.plot(t[0:i],Tsp[0:i],'r-',linewidth=3,label=r'$T_{sp}$')
plt.plot(t[0:i],T[0:i],'k.-',linewidth=3,label=r'$T_{meas}$')
plt.ylabel('T (K)')
plt.xlabel('Time (min)')
plt.legend(loc='best')
plt.subplot(5,1,5)
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)