4

I am trying to use GEKKO on PYTHON to control the level of a CSTR tank while manipulating the inlet flow q. I tried the same problem using a pid controller and it worked. However, on GEKKO the height is not tracking its setpoint. Once I did the doublet test: at a flow rate of 200, the height reached 800 and as I decreased the flowrate to 2, the height was about 0. However, when im putting the height setpoint in GEKKO as 700 or 800, the flowrate is not stopping at 200, it is continuously increasing indefinitely.

Below is my code:


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)
Mohamad Ibrahim
  • 315
  • 2
  • 8

1 Answers1

2

Mohamad,

If you want to simulate a tank level control, please start with this code with the volumetric flowrate of inlet(qin, MV) and outlet(qout, DV). 

m = GEKKO(remote=False)
m.time = [0,0.02,0.04,0.06,0.08,0.1,0.12,0.15,0.2]
h0=50
q0=10
Ac=30 # m2,Cross section Area

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

m.Equation(m.h.dt() == (m.qin - m.qout)/Ac)
Junho Park
  • 997
  • 4
  • 12
  • I will try what you suggested Junho and get back to you ! Thanks a lot – Mohamad Ibrahim Sep 28 '19 at 11:22
  • Dear Junho, thank you for your comment. I tried what you suggested,however, the measured height also did not track the setpoint. I will add a new post to show you what I wrote as a code, if u run it, you will figure out that the same problem is still there. – Mohamad Ibrahim Sep 28 '19 at 11:33