CSTR reactor

cstr.py

import numpy as np
import casadi as ca
from numpy import exp,pi



def ode():
  # Plant parameter values
  T0 = 350    # K
  c0 = 1      # kmol/m^3
  r = .219    # m
  k0 = 7.2e10 # min^-1
  E = 8750    # K
  U = 54.94   # kJ/(min m^2 K)
  rho = 1000  # kg/m^3
  Cp = .239   # kJ/(kg K)
  dH = -5e4   # kJ/kmol

  # Declare states
  c = ca.SX.sym('c') # Concentration of A
  T = ca.SX.sym('T') # Reactor temperature
  h = ca.SX.sym('h') # Level of the tank
  x = ca.vertcat((c, T, h))  # State vector

  # Disturbances
  F0 = ca.SX.sym('F0')
  d = F0

  # Declare controls
  Tc = ca.SX.sym('Tc')
  F = ca.SX.sym('F')
  u = ca.vertcat((Tc, F)) # Control vector

  # Reaction rate
  rate = k0*c*exp(-E/T)

  # Time derivatives
  c_dot = F0*(c0 - c)/(pi*r**2*h) - rate
  T_dot = F0*(T0 - T)/(pi*r**2*h) - dH/(rho*Cp)*rate + 2*U/(r*rho*Cp)*(Tc - T)
  h_dot = (F0 - F)/(pi*r**2)
  x_dot = ca.vertcat((c_dot, T_dot, h_dot))

  # NOTE: In current version of CasADi (2.3), having a "vertcat" as an input is only
  # allowed for SX, not for MX. In MX the problem has to be reformulated as in the
  # first two exercises. Support for this in MX will be added in future versions.

  # Continuous-time dynamics
  f_ct = ca.SXFunction([x, ca.vertcat((u,d))], [x_dot])
  f_ct.setOption("name","f")
  f_ct.init()
  return f_ct

cs = .878  # kmol/m^3
Ts = 324.5 # K
hs = .659  # m
xs = ca.DMatrix([cs, Ts, hs])
Tcs = 300  # K
Fs = .1    # 0.1 m^3/min
F0s = .1

def steady_state():
    # Steady-state operating conditions
    us = ca.DMatrix([Tcs, Fs])
    ds = ca.DMatrix(F0s)
    return (xs, us, ds)

def control_bounds():
    # Control bounds.
    _,us,_ = steady_state()
    delta = ca.DMatrix([.025*Tcs,.25*Fs])
    umax = us + delta
    umin = us - delta
    return umin, umax