## Main.MovingHorizonEstimation History

Changed lines 48-49 from:
1. Simulation
to:
1. Simulation
Changed line 52 from:
1. 1 step of simulation, discretization matches MHE
to:
1. 1 step of simulation, discretization matches MHE
Changed line 55 from:
1. Receive measurement from simulated control
to:
1. Receive measurement from simulated control
Changed line 60 from:
1. State variables to watch
to:
1. State variables to watch
Changed line 64 from:
1. Other parameters
to:
1. other parameters
Changed line 76 from:
1. Variables
to:
1. Variables
Changed line 80 from:
1. Rate equations
to:
1. Rate equations
Changed line 83 from:
1. CSTR equations
to:
1. CSTR equations
Changed line 87 from:
1. Options
to:
1. Options
Changed lines 93-96 from:
1. MHE
2. Model
to:
1. MHE
2. Model
Changed line 99 from:
1. 6 time points in horizon
to:
1. 6 time points in horizon
Changed line 102 from:
1. Parameter to Estimate
to:
1. Parameter to Estimate
Changed line 106 from:
1. Upper and lower bounds for optimizer
to:
1. upper and lower bounds for optimizer
Changed line 110 from:
1. Measurement input
to:
1. Measurement input
Changed line 115 from:
1. Measurement to match simulation with
to:
1. Measurement to match simulation with
Changed line 121 from:
1. State to watch
to:
1. State to watch
Changed line 124 from:
1. Other parameters
to:
1. Other parameters
Changed line 135 from:
1. Equation variables (2 other DOF from CV and FV)
to:
1. Equation variables(2 other DOF from CV and FV)
Changed line 139 from:
1. Reaction equations
to:
1. Reaction equations
Changed line 142 from:
1. CSTR equations
to:
1. CSTR equations
Changed lines 146-147 from:
1. Global Tuning
to:
1. Global Tuning
Changed lines 153-154 from:
1. Loop
to:
1. Loop
Changed line 211 from:
1. plot results
to:
1. plot results
January 24, 2018, at 12:14 AM by 10.37.134.137 -
Changed lines 36-236 from:
to:
1. get latest gekko packge with:
2. pip install gekko
3. or
5. to upgrade the version to the latest
6. or from the Python script:
7. import pip
8. pip.main(['install','gekko'])

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

1. Simulation

s = GEKKO(name='cstr-sim')

1. 1 step of simulation, discretization matches MHE

s.time = np.linspace(0,.1,2)

1. Receive measurement from simulated control

Tc = s.MV(value=300,name='tc') Tc.FSTATUS = 1 #receive measurement Tc.STATUS = 0 #don't optimize

1. State variables to watch

Ca = s.SV(value=.7, ub=1, lb=0,name='ca') T = s.SV(value=305,lb=250,ub=500,name='t')

1. Other parameters

q = s.Param(value=100) V = s.Param(value=100) rho = s.Param(value=1000) Cp = s.Param(value=0.239) mdelH = s.Param(value=50000) ER = s.Param(value=8750) k0 = s.Param(value=7.2*10**10) UA = s.Param(value=5*10**4) Ca0 = s.Param(value=1) T0 = s.Param(value=350)

1. Variables

k = s.Var() rate = s.Var()

1. Rate equations

s.Equation(k==k0*s.exp(-ER/T)) s.Equation(rate==k*Ca)

1. CSTR equations

s.Equation(V* Ca.dt() == q*(Ca0-Ca)-V*rate) s.Equation(rho*Cp*V* T.dt() == q*rho*Cp*(T0-T) + V*mdelH*rate + UA*(Tc-T))

1. Options

s.options.IMODE = 4 #dynamic simulation s.options.NODES = 3 s.options.SOLVER = 3

1. MHE
2. Model

m = GEKKO(name='cstr-mhe')

1. 6 time points in horizon

m.time = np.linspace(0,.5,6)

1. Parameter to Estimate

UA_mhe = m.FV(value=10*10**4,name='ua') UA_mhe.STATUS = 1 #estimate UA_mhe.FSTATUS = 0 #no measurements

1. Upper and lower bounds for optimizer

UA_mhe.LOWER = 10000 UA_mhe.UPPER = 100000

1. Measurement input

Tc_mhe = m.MV(value=300,name='tc') Tc_mhe.STATUS = 0 #don't estimate Tc_mhe.FSTATUS = 1 #receive measurement

1. Measurement to match simulation with

T_mhe = m.CV(value=325 ,lb=250,ub=500,name='t') T_mhe.STATUS = 1 #minimize error between simulation and measurement T_mhe.FSTATUS = 1 #receive measurement T_mhe.MEAS_GAP = 0.1 #measurement deadband gap

1. State to watch

Ca_mhe = m.SV(value=0.8, ub=1, lb=0,name='ca')

1. Other parameters

q = m.Param(value=100) V = m.Param(value=100) rho = m.Param(value=1000) Cp = m.Param(value=0.239) mdelH = m.Param(value=50000) ER = m.Param(value=8750) k0 = m.Param(value=7.2*10**10) Ca0 = m.Param(value=1) T0 = m.Param(value=350)

1. Equation variables (2 other DOF from CV and FV)

k = m.Var() rate = m.Var()

1. Reaction equations

m.Equation(k==k0*m.exp(-ER/T_mhe)) m.Equation(rate==k*Ca_mhe)

1. CSTR equations

m.Equation(V* Ca_mhe.dt() == q*(Ca0-Ca_mhe)-V*rate) #mol balance m.Equation(rho*Cp*V* T_mhe.dt() == q*rho*Cp*(T0-T_mhe) + V*mdelH*rate + UA_mhe*(Tc_mhe-T_mhe)) #energy balance

1. Global Tuning

m.options.IMODE = 5 #MHE m.options.EV_TYPE = 1 m.options.NODES = 3 m.options.SOLVER = 3 #IPOPT

1. Loop
2. number of cycles to run

cycles = 50

1. step in the jacket cooling temperature at cycle 6

Tc_meas = np.empty(cycles) Tc_meas[0:15] = 280 Tc_meas[5:cycles] = 300 dt = 0.1 # min time = np.linspace(0,cycles*dt-dt,cycles) # time points for plot

1. allocate storage

Ca_meas = np.empty(cycles) T_meas = np.empty(cycles) UA_mhe_store = np.empty(cycles) Ca_mhe_store = np.empty(cycles) T_mhe_store = np.empty(cycles)

for i in range(cycles):

    ## Process
# input Tc (jacket cooling temperature)
Tc.MEAS = Tc_meas[i]
# simulate process model, 1 time step
s.solve()
# retrieve Ca and T measurements from the process
Ca_meas[i] = Ca.MODEL
T_meas[i] = T.MODEL

## Estimator
# input process measurements
# input Tc (jacket cooling temperature)
Tc_mhe.MEAS = Tc_meas[i]
# input T (reactor temperature)
T_mhe.MEAS = T_meas[i] #CV

# solve process model, 1 time step
m.solve()
# check if successful
if m.options.APPSTATUS == 1:
# retrieve solution
UA_mhe_store[i] = UA_mhe.NEWVAL
Ca_mhe_store[i] = Ca_mhe.MODEL
T_mhe_store[i] = T_mhe.MODEL
else:
# failed solution
UA_mhe_store[i] = 0
Ca_mhe_store[i] = 0
T_mhe_store[i] = 0

print('MHE results: Ca (estimated)=' + str(Ca_mhe_store[i]) +         ' Ca (actual)=' + str(Ca_meas[i]) +         ' UA (estimated)=' + str(UA_mhe_store[i]) +         ' UA (actual)=50000')

1. plot results

plt.figure() plt.subplot(411) plt.plot(time,Tc_meas,'k-',linewidth=2) plt.axis([0,time[-1],275,305]) plt.ylabel('Jacket T (K)') plt.legend('T_c')

plt.subplot(412) plt.plot([0,time[-1]],[50000,50000],'k--') plt.plot(time,UA_mhe_store,'r:',linewidth=2) plt.axis([0,time[-1],10000,100000]) plt.ylabel('UA') plt.legend(['Actual UA','Predicted UA'],loc=4)

plt.subplot(413) plt.plot(time,T_meas,'ro') plt.plot(time,T_mhe_store,'b-',linewidth=2) plt.axis([0,time[-1],300,340]) plt.ylabel('Reactor T (K)') plt.legend(['Measured T','Predicted T'],loc=4)

plt.subplot(414) plt.plot(time,Ca_meas,'go') plt.plot(time,Ca_mhe_store,'m-',linewidth=2) plt.axis([0,time[-1],.6,1]) plt.ylabel('Reactor C_a (mol/L)') plt.legend(['Measured C_a','Predicted C_a'],loc=4) plt.show()

(:toggle hide gekko button show="Show GEKKO (Python) Code":) (:div id=gekko:) (:source lang=python:)

(:sourceend:) (:divend:)

Changed lines 11-12 from:

Objective: Design an estimator to predict an unknown parameter and state variable. Develop a model of the reactor and implement the model with moving horizon estimation in MATLAB, Python, or Simulink. Estimated time: 3 hours.

to:

Objective: Design an estimator to predict an unknown parameter and state variable. Use a model of the reactor and implement the estimator to detect the current states (temperature and concentration) as well as the unmeasured heat transfer coefficient (U). Estimated time: 2-3 hours.

The estimator can be any type such as a Kalman filter, Extended Kalman filter, Unscented Kalman Filter (particle filter), or an observer that can detect the states (T and Ca) along with the unknown parameter (U). The following solutions demonstrate an implementation of Moving Horizon Estimation.

Changed lines 31-39 from:

(:htmlend:)

to:

(:htmlend:)

#### References

1. Haseltine, E.L., Rawlings, J.B., Critical Evaluation of Extended Kalman Filtering and Moving-Horizon Estimation, Industrial & Engineering Chemistry Research 2005 44 (8), 2451-2460, DOI: 10.1021/ie034308l Preprint, Article
2. Spivey, B.J., Hedengren, J.D., Edgar, T.F., Constrained Nonlinear Estimation for Industrial Process Fouling, Industrial & Engineering Chemistry Research 2010 49 (17), 7824-7831, DOI: 10.1021/ie9018116 Article
3. Hedengren, J. D., Eaton, A. N., Overview of Estimation Methods for Industrial Dynamic Systems, Optimization and Engineering, Springer, Vol 18 (1), 2017, pp. 155-178, DOI: 10.1007/s11081-015-9295-9. Preprint, Article
May 13, 2015, at 10:59 PM by 10.10.147.4 -
Deleted lines 7-8:
May 13, 2015, at 06:58 PM by 10.5.113.160 -
May 10, 2015, at 02:21 AM by 45.56.3.184 -
Changed lines 25-29 from:
to:

(:html:) <iframe width="560" height="315" src="https://www.youtube.com/embed/yQWgSByYjd8?rel=0" frameborder="0" allowfullscreen></iframe> (:htmlend:)

May 09, 2015, at 10:50 PM by 45.56.3.184 -
Changed line 11 from:

Objective: Design an estimator to predict an unknown parameter and state variable. Develop a model of the reactor and implement the model with moving horizon estimation in Simulink. Estimated time: 3 hours.

to:

Objective: Design an estimator to predict an unknown parameter and state variable. Develop a model of the reactor and implement the model with moving horizon estimation in MATLAB, Python, or Simulink. Estimated time: 3 hours.

May 09, 2015, at 10:49 PM by 45.56.3.184 -
Changed lines 17-19 from:

A reactor is used to convert a hazardous chemical A to an acceptable chemical B in waste stream before entering a nearby lake. This particular reactor is dynamically modeled as a Continuously Stirred Tank Reactor (CSTR) with a simplified kinetic mechanism that describes the conversion of reactant A to product B with an irreversible and exothermic reaction. It is desired to maintain the temperature at a constant setpoint that maximizes the destruction of A (highest possible temperature).

Design an estimator to predict the concentration of A leaving the reactor and the heat transfer coefficient UA from the measured reactor temperature T and jacket temperature Tc. See a related CSTR case study for details on the model.

to:

A reactor is used to convert a hazardous chemical A to an acceptable chemical B in waste stream before entering a nearby lake. This particular reactor is dynamically modeled as a Continuously Stirred Tank Reactor (CSTR) with a simplified kinetic mechanism that describes the conversion of reactant A to product B with an irreversible and exothermic reaction. It is desired to maintain the temperature at a constant setpoint that maximizes the destruction of A (highest possible temperature). First, however, an estimator must predict the concentration of A because there is no direct measurement of this quantity. The reaction kinetics and dynamic equations are well-known but there is a parameter in the model, the heat transfer coefficient UA, that is unknown.

Design an estimator to predict the concentration of A leaving the reactor and the heat transfer coefficient UA from the measured reactor temperature T and jacket temperature Tc. See a related CSTR case study for details on the model.

May 09, 2015, at 01:44 PM by 45.56.3.184 -
Changed lines 23-24 from:
to:
May 09, 2015, at 01:43 PM by 45.56.3.184 -
Changed lines 23-24 from:
to:
May 09, 2015, at 06:14 AM by 45.56.3.184 -
Changed line 19 from:

Design an estimator to predict the concentration of A leaving the reactor and the heat transfer coefficient UA from the measured reactor temperature T and jacket temperature Tc.

to:

Design an estimator to predict the concentration of A leaving the reactor and the heat transfer coefficient UA from the measured reactor temperature T and jacket temperature Tc. See a related CSTR case study for details on the model.

May 09, 2015, at 06:12 AM by 45.56.3.184 -
Changed line 23 from:
to:
May 09, 2015, at 06:11 AM by 45.56.3.184 -
Changed lines 21-23 from:

to:

#### Solution

May 09, 2015, at 12:42 AM by 10.5.113.160 -
Changed line 19 from:

Design an estimator to predict the concentration of A leaving the reactor and the heat transfer coefficient UA from the measured reactor temperature T and jacket temperature T'_c'.

to:

Design an estimator to predict the concentration of A leaving the reactor and the heat transfer coefficient UA from the measured reactor temperature T and jacket temperature Tc.

May 09, 2015, at 12:41 AM by 10.5.113.160 -
Changed line 19 from:

Design an estimator to predict the concentration of A leaving the reactor and the heat transfer coefficient UA from the measured temperature.

to:

Design an estimator to predict the concentration of A leaving the reactor and the heat transfer coefficient UA from the measured reactor temperature T and jacket temperature T'_c'.

May 09, 2015, at 12:40 AM by 10.5.113.160 -

Design an estimator to predict the concentration of A leaving the reactor and the heat transfer coefficient UA from the measured temperature.

May 09, 2015, at 12:36 AM by 10.5.113.160 -
Changed line 13 from:
to:
May 09, 2015, at 12:36 AM by 10.5.113.160 -
May 08, 2015, at 06:38 PM by 10.5.113.160 -

#### Exercise

Objective: Design an estimator to predict an unknown parameter and state variable. Develop a model of the reactor and implement the model with moving horizon estimation in Simulink. Estimated time: 3 hours.

A reactor is used to convert a hazardous chemical A to an acceptable chemical B in waste stream before entering a nearby lake. This particular reactor is dynamically modeled as a Continuously Stirred Tank Reactor (CSTR) with a simplified kinetic mechanism that describes the conversion of reactant A to product B with an irreversible and exothermic reaction. It is desired to maintain the temperature at a constant setpoint that maximizes the destruction of A (highest possible temperature).

#### Solution

April 30, 2015, at 10:56 PM by 45.56.3.184 -
Changed lines 7-8 from:

Moving Horizon Estimation (MHE) uses dynamic optimization and a backward time horizon of measurements to optimally adjust parameters and states. The data may include noise (random fluctuations), drift (gradual departure from true values), outliers (sudden and temporary departure from true values), or other inaccuracies. Nonlinear programming solvers are employed to numerically converge the dynamic optimization problem.

to:

Moving Horizon Estimation (MHE) uses dynamic optimization and a backward time horizon of measurements to optimally adjust parameters and states. The data may include noise (random fluctuations), drift (gradual departure from true values), outliers (sudden and temporary departure from true values), or other inaccuracies. Nonlinear programming solvers are employed to numerically converge the dynamic optimization problem.