CODEWITHSUNDEEP
pythongekko
Junho Park
Junho Park

Reputation: 997

Steady State parameter estimation for CSTR using GEKKO

I would like to fit the reaction constants (k0 and EoverR) using steady-state concentration data of a CSTR with respect to the reactor temperature. I used just a simple linear equation to generate the operation data to fit. (Ca_data = -1.5*T_reactor/100 + 4.2)

Because this is steady-state data, there is no time step (m.time) needed. Please give advice on how to convert below simulation code into an estimation code for 'Ca vs. T_reactor'.

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO

# Feed Temperature (K)
Tf = 350
# Feed Concentration (mol/m^3)
Caf = 1

# Steady State Initial Conditions for the States
Ca_ss = 1
T_ss = 304

#%% GEKKO
m = GEKKO(remote=True)
m.time = np.linspace(0, 25, 251)

# Volumetric Flowrate (m^3/sec)
q = 100
# 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 = 8700
# Pre-exponential factor (1/sec)
k0 = 3.2e15
# 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 = 280
T0 = 304
Ca0 = 1.0

T = m.MV(value=T_ss)
rA = m.Var(value=0)
Ca = m.CV(value=Ca_ss)

m.Equation(rA == k0*m.exp(-EoverR/T)*Ca)
m.Equation(Ca.dt() == q/V*(Caf - Ca) - rA)

m.options.IMODE = 1
m.options.SOLVER = 3

T_reactor = np.linspace(220, 260, 11)
Ca_results = np.zeros(np.size(T_reactor))
for i in range(np.size(T_reactor)):
    T.Value = T_reactor[i]
    m.solve(disp=True)
    Ca_results[i] = Ca[-1]

Ca_data = -1.5*T_reactor/100 + 4.2 # for generating the operation data

# Plot the results
plt.plot(T_reactor,Ca_data,'bo',linewidth=3)
plt.plot(T_reactor,Ca_results,'r-',linewidth=3)
plt.ylabel('Ca (mol/L)')
plt.xlabel('Temperature (K)')
plt.legend(['Reactor Concentration'],loc='best')
plt.show()

Upvotes: 3

Views: 340

Answers (1)

John Hedengren
John Hedengren

Reputation: 14346

There is a steady state estimation mode (IMODE=2) in Gekko for linear or nonlinear regression. Two examples are nonlinear regression and energy price regression. For the posted problem, here are a few recommendations:

  • Solve the regression problem with one solve, instead of in a loop. That way, the parameters that you select will fit over the entire range, not just at one point.
  • Identify the parameters that should be adjusted to minimize the error between the data and model predictions. These should be an m.FV() type for parameters that are one value and m.MV() for parameters that have a different value for each time point.
  • Set Ca.FSTATUS=1 to tell the solver that it should try to match the predictions of Ca to the data loaded in Ca.value.
  • Set kf.STATUS=1 to tell the solver that it is a parameter that should be adjusted to minimize the Ca error.
  • Optional: make kf the adjustable parameter instead of k0 directly to improve problem scaling. Large values (such as >1e10 or <-1e10) can create problems for the solver (without automatic scaling) because the gradients are small. I created the new parameter kf as an Intermediate to also remove an additional equation.

Regression Fit

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
m = GEKKO(remote=True)
Tf = 350
Caf = 1
q = 100
V = 100
rho = 1000
Cp = 0.239
mdelH = 5e4
EoverR = 8700
k0 = 3.2e15
UA = 5e4
T = m.MV()
Ca = m.CV()

# new parameter to estimate
kf = m.FV(1,lb=0.5,ub=2.0)
kf.STATUS = 1

rA = m.Intermediate(kf*k0*m.exp(-EoverR/T)*Ca)
m.Equation(Ca.dt() == q/V*(Caf - Ca) - rA)

m.options.IMODE = 2
m.options.SOLVER = 3

# generate data
T_reactor = np.linspace(220, 260, 11)
Ca_data = -1.5*T_reactor/100 + 4.2

# insert data
T.value = T_reactor
Ca.value = Ca_data
Ca.FSTATUS = 1 # fit Ca

m.solve()

print('kf = ' + str(kf.value[0]))
print('k = ' + str(kf.value[0]*k0))

# Plot the results
plt.plot(T_reactor,Ca_data,'bo',linewidth=3)
plt.plot(T_reactor,Ca.value,'r-',linewidth=3)
plt.ylabel('Ca (mol/L)')
plt.xlabel('Temperature (K)')
plt.legend(['Reactor Concentration','Regression Fit'],loc='best')
plt.show()

You can select any number of parameters to estimate to improve the fit. It doesn't need to be limited to just kf. Your post mentions that EoverR is another potential parameter to estimate but this may not improve the fit substantially because k0 and EoverR are co-linear. Both parameters can be increased or decreased and give nearly the same solution. As a word of caution, there needs to be significant temperature variation to estimate both.

Upvotes: 1

Related Questions