2

我正在尝试在 PYTHON 上使用 GEKKO 来控制 CSTR 罐的液位，同时操纵入口流量 q。我使用 pid 控制器尝试了同样的问题，它成功了。然而，在 GEKKO 上，高度并未跟踪其设定点。一旦我做了双重测试：在 200 的流速下，高度达到 800，当我将流速降低到 2 时，高度约为 0。但是，当我将 GEKKO 中的高度设定点设置为 700 或 800 时，流速并没有止步于200，而是无限地不断增加。另外，我试过把 qout=5.0, Ac=30 和 h0=50，也没有用。

下面是我的代码：


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

m.q=m.MV(value=q0,lb=0,ub=100)
m.h= m.CV(value=h0)
m.qout=m.Param(value=5)

m.Equation(m.h.dt()==(m.q- m.qout)/Ac)


#MV tuning
m.q.STATUS = 1
m.q.FSTATUS = 0
m.q.DMAX = 100
m.q.DMAXHI = 5  
m.q.DMAXLO = -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.0)/Ac

    # Return xdot:
    xdot = np.zeros(1)
    xdot[0]= dhdt
    return xdot


# Time Interval (min)
t = np.linspace(0,6,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] = 55.0
hsp[100:150]=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+1],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)
4

1 回答 1

2

问题在于您对模拟器的输入y = odeint(cstr,x0,ts,args=(u2[i+1],Ac))。它应该使用前一个循环中的值：y = odeint(cstr,x0,ts,args=(u2[i],Ac)). 这是一个更新的脚本。

仿真结果

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

m.q=m.MV(value=q0,lb=0,ub=100)
m.h= m.CV(value=h0)
m.qout=m.Param(value=5)

m.Equation(Ac * m.h.dt()==m.q- m.qout)


#MV tuning
m.q.STATUS = 1
m.q.FSTATUS = 0
m.q.DMAX = 100
m.q.DMAXHI = 5  
m.q.DMAXLO = -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,6,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] = 55.0
hsp[100:150]=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=False)


    # retrieve new q value

    u2[i+1] = m.q.NEWVAL
    # update initial conditions
    x0[0]= h[i+1]

    #%% Plot the results
    if i%10==0:
        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)
于 2019-10-03T14:35:50.977 回答