## Moving Horizon Estimation

## Main.MovingHorizonEstimation History

Hide minor edits - Show changes to output

Changed line 205 from:

plt.plot(time,Tc_meas,'k-',~~linewidth~~=2)

to:

plt.plot(time,Tc_meas,'k-',lw=2)

Changed line 212 from:

plt.plot(time,UA_mhe_store,'r:',~~linewidth~~=2)

to:

plt.plot(time,UA_mhe_store,'r:',lw=2)

Changed line 219 from:

plt.plot(time,T_mhe_store,'b-',~~linewidth~~=2)

to:

plt.plot(time,T_mhe_store,'b-',lw=2)

Changed line 226 from:

plt.plot(time,Ca_mhe_store,'m-',~~linewidth~~=2)

to:

plt.plot(time,Ca_mhe_store,'m-',lw=2)

Deleted lines 195-200:

if i==0:

# remove '.value' overrides to use '.MEAS'

for xi in [Ca,T,Tc]:

xi.value=[None]*len(s.time)

Tc_mhe.value=[None]*len(m.time)

Added line 36:

# Moving Horizon Estimation

Changed lines 49-53 from:

#~~%%~~ Simulation

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

#~~1~~ step of simulation, discretization matches MHE

s = GEKKO(

#

to:

# Simulation

s = GEKKO(remote=False,name='cstr-sim')

# One step of simulation, discretization matches MHE

s = GEKKO(remote=False,name='cstr-sim')

# One step of simulation, discretization matches MHE

Changed line 55 from:

#Receive measurement from simulated control

to:

# Receive measurement from simulated control

Changed line 60 from:

#~~State variables to watch~~

to:

# Simulator variables

Changed lines 62-64 from:

T = s.SV(value=~~305~~,lb=250,ub=500,name='t')

#~~other parameters~~

#

to:

T = s.SV(value=335,lb=250,ub=500,name='t')

# Parameters

# Parameters

Changed lines 71-72 from:

k0 = s.Param(value=7.~~2*10**10~~)

UA = s.Param(value=~~5*10**4~~)

UA = s.Param(value=

to:

k0 = s.Param(value=7.2e10)

UA = s.Param(value=5e4)

UA = s.Param(value=5e4)

Changed line 76 from:

#Variables

to:

# Variables

Changed line 80 from:

#Rate equations

to:

# Rate equations

Changed line 83 from:

#CSTR equations

to:

# CSTR equations

Changed line 87 from:

#Options

to:

# Options

Changed lines 93-107 from:

#~~%%~~ MHE

#Model

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

#~~6~~ time points in horizon

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

#Parameter to Estimate

UA_mhe = m.~~FV~~(~~value=10*10**4~~,~~name='ua'~~)

~~UA_mhe.STATUS = 1 #estimate~~

UA_mhe~~.FSTATUS~~ = ~~0 #no measurements~~

#upper and lower bounds for optimizer

UA_mhe.~~LOWER~~ = ~~10000~~

#Model

m = GEKKO(

#

#Parameter to Estimate

UA_mhe = m

UA_mhe

#upper and lower bounds for optimizer

UA_mhe.

to:

# MHE

m = GEKKO(remote=False,name='cstr-mhe')

# 11 time points in horizon (0.1 each step)

m.time = np.linspace(0,1.0,11)

# Parameter to Estimate

UA_mhe = m.FV(value=1e4,name='ua')

UA_mhe.STATUS = 1 # estimate

UA_mhe.FSTATUS = 0 # no measurements

# Upper and lower bounds for optimizer

UA_mhe.LOWER = 30000

m = GEKKO(remote=False,name='cstr-mhe')

# 11 time points in horizon (0.1 each step)

m.time = np.linspace(0,1.0,11)

# Parameter to Estimate

UA_mhe = m.FV(value=1e4,name='ua')

UA_mhe.STATUS = 1 # estimate

UA_mhe.FSTATUS = 0 # no measurements

# Upper and lower bounds for optimizer

UA_mhe.LOWER = 30000

Changed line 107 from:

#~~Measurement input~~

to:

# Cooling Jacket Temperature

Changed lines 109-121 from:

Tc_mhe.STATUS = 0 #don't estimate

Tc_mhe.FSTATUS = 1 #receive measurement

#~~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~~

#State to watch

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

#~~Other parameters~~

Tc_mhe.FSTATUS = 1 #receive measurement

#

T_mhe =

T_mhe.

T_mhe.

#State to watch

Ca_mhe =

#

to:

Tc_mhe.STATUS = 0 # don't estimate

Tc_mhe.FSTATUS = 1 # receive measurement

# Reactor Temperature

T_mhe = m.CV(value=335,lb=250,ub=500,name='t')

T_mhe.FSTATUS = 1 # minimize error with measurement

T_mhe.MEAS_GAP = 0.1 # measurement deadband gap

# Reactor Concentration

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

# Parameters

Tc_mhe.FSTATUS = 1 # receive measurement

# Reactor Temperature

T_mhe = m.CV(value=335,lb=250,ub=500,name='t')

T_mhe.FSTATUS = 1 # minimize error with measurement

T_mhe.MEAS_GAP = 0.1 # measurement deadband gap

# Reactor Concentration

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

# Parameters

Changed line 127 from:

k0 = m.Param(value=7.~~2*10**10~~)

to:

k0 = m.Param(value=7.2e10)

Changed line 131 from:

#Equation variables(2 other DOF from CV and FV)

to:

# Equation variables (2 other DOF from CV and FV)

Changed line 135 from:

#Reaction equations

to:

# Reaction equations

Changed lines 138-144 from:

#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~~

#Global Tuning

m.options.IMODE = 5 #MHE

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)

#Global Tuning

m.options.IMODE = 5

to:

# 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

# Global Tuning

m.options.IMODE = 5 # MHE

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

# Global Tuning

m.options.IMODE = 5 # MHE

Changed lines 147-151 from:

m.options.SOLVER = 3 #IPOPT

#~~%% Loop~~

# number of cycles to run

#

# number of cycles

to:

m.options.SOLVER = 3 # IPOPT

# Cycles to run

# Cycles to run

Changed lines 157-158 from:

time ~~=~~ np.linspace(0,cycles*dt-dt,cycles) ~~# time points for plot~~

to:

# time points for plot

time = np.linspace(0,cycles*dt-dt,cycles)

time = np.linspace(0,cycles*dt-dt,cycles)

Deleted lines 166-167:

Changed lines 168-169 from:

to:

# Process

Changed lines 172-173 from:

s.solve()

# retrieve Ca and T~~measurements from the process~~

# retrieve Ca and T

to:

s.solve(disp=False)

# retrieve Ca and T measurements

# retrieve Ca and T measurements

Changed line 177 from:

to:

# Estimator

Changed lines 183-185 from:

to:

# Solve MHE

m.solve(disp=False)

m.solve(disp=False)

Changed lines 197-203 from:

to:

if i==0:

# remove '.value' overrides to use '.MEAS'

for xi in [Ca,T,Tc]:

xi.value=[None]*len(s.time)

Tc_mhe.value=[None]*len(m.time)

print('MHE: Ca (est)=' + str(Ca_mhe_store[i]) + \

# remove '.value' overrides to use '.MEAS'

for xi in [Ca,T,Tc]:

xi.value=[None]*len(s.time)

Tc_mhe.value=[None]*len(m.time)

print('MHE: Ca (est)=' + str(Ca_mhe_store[i]) + \

Changed line 205 from:

' UA (~~estimated~~)=' + str(UA_mhe_store[i]) + \

to:

' UA (est)=' + str(UA_mhe_store[i]) + \

Changed line 208 from:

#~~%%~~ plot results

to:

# plot results

Changed lines 48-49 from:

# Simulation

to:

#%% Simulation

Changed line 52 from:

#~~ ~~1 step of simulation, discretization matches MHE

to:

#1 step of simulation, discretization matches MHE

Changed line 55 from:

#~~ ~~Receive measurement from simulated control

to:

#Receive measurement from simulated control

Changed line 60 from:

#~~ ~~State variables to watch

to:

#State variables to watch

Changed line 64 from:

#~~ Other~~ parameters

to:

#other parameters

Changed line 76 from:

#~~ ~~Variables

to:

#Variables

Changed line 80 from:

#~~ ~~Rate equations

to:

#Rate equations

Changed line 83 from:

#~~ ~~CSTR equations

to:

#CSTR equations

Changed line 87 from:

#~~ ~~Options

to:

#Options

Changed lines 93-96 from:

# MHE

#~~ ~~Model

#

to:

#%% MHE

#Model

#Model

Changed line 99 from:

#~~ ~~6 time points in horizon

to:

#6 time points in horizon

Changed line 102 from:

#~~ ~~Parameter to Estimate

to:

#Parameter to Estimate

Changed line 106 from:

#~~ Upper~~ and lower bounds for optimizer

to:

#upper and lower bounds for optimizer

Changed line 110 from:

#~~ ~~Measurement input

to:

#Measurement input

Changed line 115 from:

#~~ ~~Measurement to match simulation with

to:

#Measurement to match simulation with

Changed line 121 from:

#~~ ~~State to watch

to:

#State to watch

Changed line 124 from:

#~~ ~~Other parameters

to:

#Other parameters

Changed line 135 from:

#~~ ~~Equation variables~~ ~~(2 other DOF from CV and FV)

to:

#Equation variables(2 other DOF from CV and FV)

Changed line 139 from:

#~~ ~~Reaction equations

to:

#Reaction equations

Changed line 142 from:

#~~ ~~CSTR equations

to:

#CSTR equations

Changed lines 146-147 from:

#~~ ~~Global Tuning

to:

#Global Tuning

Changed lines 153-154 from:

# Loop

to:

#%% Loop

Added lines 172-173:

Changed line 211 from:

# plot results

to:

#%% plot results

Changed lines 36-236 from:

to:

# get latest gekko packge with:

# pip install gekko

# or

# pip install gekko --upgrade

# to upgrade the version to the latest

# or from the Python script:

# import pip

# pip.main(['install','gekko'])

from gekko import GEKKO

import numpy as np

import matplotlib.pyplot as plt

# Simulation

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

# 1 step of simulation, discretization matches MHE

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

# Receive measurement from simulated control

Tc = s.MV(value=300,name='tc')

Tc.FSTATUS = 1 #receive measurement

Tc.STATUS = 0 #don't optimize

# 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')

# 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)

# Variables

k = s.Var()

rate = s.Var()

# Rate equations

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

s.Equation(rate==k*Ca)

# 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))

# Options

s.options.IMODE = 4 #dynamic simulation

s.options.NODES = 3

s.options.SOLVER = 3

# MHE

# Model

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

# 6 time points in horizon

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

# Parameter to Estimate

UA_mhe = m.FV(value=10*10**4,name='ua')

UA_mhe.STATUS = 1 #estimate

UA_mhe.FSTATUS = 0 #no measurements

# Upper and lower bounds for optimizer

UA_mhe.LOWER = 10000

UA_mhe.UPPER = 100000

# Measurement input

Tc_mhe = m.MV(value=300,name='tc')

Tc_mhe.STATUS = 0 #don't estimate

Tc_mhe.FSTATUS = 1 #receive measurement

# 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

# State to watch

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

# 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)

# Equation variables (2 other DOF from CV and FV)

k = m.Var()

rate = m.Var()

# Reaction equations

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

m.Equation(rate==k*Ca_mhe)

# 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

# Global Tuning

m.options.IMODE = 5 #MHE

m.options.EV_TYPE = 1

m.options.NODES = 3

m.options.SOLVER = 3 #IPOPT

# Loop

# number of cycles to run

cycles = 50

# 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

# 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')

# 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()

# pip install gekko

# or

# pip install gekko --upgrade

# to upgrade the version to the latest

# or from the Python script:

# import pip

# pip.main(['install','gekko'])

from gekko import GEKKO

import numpy as np

import matplotlib.pyplot as plt

# Simulation

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

# 1 step of simulation, discretization matches MHE

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

# Receive measurement from simulated control

Tc = s.MV(value=300,name='tc')

Tc.FSTATUS = 1 #receive measurement

Tc.STATUS = 0 #don't optimize

# 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')

# 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)

# Variables

k = s.Var()

rate = s.Var()

# Rate equations

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

s.Equation(rate==k*Ca)

# 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))

# Options

s.options.IMODE = 4 #dynamic simulation

s.options.NODES = 3

s.options.SOLVER = 3

# MHE

# Model

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

# 6 time points in horizon

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

# Parameter to Estimate

UA_mhe = m.FV(value=10*10**4,name='ua')

UA_mhe.STATUS = 1 #estimate

UA_mhe.FSTATUS = 0 #no measurements

# Upper and lower bounds for optimizer

UA_mhe.LOWER = 10000

UA_mhe.UPPER = 100000

# Measurement input

Tc_mhe = m.MV(value=300,name='tc')

Tc_mhe.STATUS = 0 #don't estimate

Tc_mhe.FSTATUS = 1 #receive measurement

# 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

# State to watch

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

# 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)

# Equation variables (2 other DOF from CV and FV)

k = m.Var()

rate = m.Var()

# Reaction equations

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

m.Equation(rate==k*Ca_mhe)

# 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

# Global Tuning

m.options.IMODE = 5 #MHE

m.options.EV_TYPE = 1

m.options.NODES = 3

m.options.SOLVER = 3 #IPOPT

# Loop

# number of cycles to run

cycles = 50

# 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

# 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')

# 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()

Added lines 32-38:

(:toggle hide gekko button show="Show GEKKO (Python) Code":)

(:div id=gekko:)

(:source lang=python:)

# add source

(: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.''

Added lines 23-24:

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

# 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 [[https://jbrwww.che.wisc.edu/tech-reports/twmcc-2002-03.pdf|Preprint]], [[https://pubs.acs.org/doi/abs/10.1021/ie034308l|Article]]

# 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 [[https://pubs.acs.org/doi/abs/10.1021/ie9018116|Article]]

# 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. [[https://apm.byu.edu/prism/uploads/Members/eaton_hedengren_OPTE_springer.pdf|Preprint]], [[https://link.springer.com/article/10.1007/s11081-015-9295-9|Article]]

!!!! References

# 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 [[https://jbrwww.che.wisc.edu/tech-reports/twmcc-2002-03.pdf|Preprint]], [[https://pubs.acs.org/doi/abs/10.1021/ie034308l|Article]]

# 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 [[https://pubs.acs.org/doi/abs/10.1021/ie9018116|Article]]

# 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. [[https://apm.byu.edu/prism/uploads/Members/eaton_hedengren_OPTE_springer.pdf|Preprint]], [[https://link.springer.com/article/10.1007/s11081-015-9295-9|Article]]

Deleted lines 7-8:

Attach:download.png [[Attach:mhe_simulink.zip|MHE Example in Simulink]]

Added lines 8-9:

Attach:download.png [[Attach:mhe_simulink.zip|MHE Example in Simulink]]

Changed lines 25-29 from:

Attach:download.png [[Attach:cstr_mhe_solution.zip|CSTR MHE Solution in Simulink]]

to:

Attach:download.png [[Attach:cstr_mhe_solution.zip|CSTR MHE Solution in Simulink]]

(:html:)

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

(:htmlend:)

(:html:)

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

(:htmlend:)

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.''

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 ~~''~~T~~'~~_c_~~'''. See a [[https://apmonitor.com/che436/index.php/Main/CaseStudyCSTR|related CSTR case study]] for details on the model.

Design an estimator to predict the concentration of

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 '''T'_c_''''. See a [[https://apmonitor.com/che436/index.php/Main/CaseStudyCSTR|related CSTR case study]] for details on the model.

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_''''. See a [[https://apmonitor.com/che436/index.php/Main/CaseStudyCSTR|related CSTR case study]] for details on the model.

Changed lines 23-24 from:

Attach:download.png [[Attach:cstr_mhe_solution_~~python~~.zip|CSTR MHE Solution in ~~Python~~]]Attach:download.png [[Attach:cstr_mhe_solution_~~matlab~~.zip|CSTR MHE Solution in ~~MATLAB~~]]

to:

Attach:download.png [[Attach:cstr_mhe_solution_matlab.zip|CSTR MHE Solution in MATLAB]]

Attach:download.png [[Attach:cstr_mhe_solution_python.zip|CSTR MHE Solution in Python]]

Attach:download.png [[Attach:cstr_mhe_solution_python.zip|CSTR MHE Solution in Python]]

Changed lines 23-24 from:

Attach:download.png [[Attach:cstr_mhe_solution.zip|CSTR MHE Solution in Simulink]]

to:

Attach:download.png [[Attach:cstr_mhe_solution_python.zip|CSTR MHE Solution in Python]]Attach:download.png [[Attach:cstr_mhe_solution_matlab.zip|CSTR MHE Solution in MATLAB]]

Attach:download.png [[Attach:cstr_mhe_solution.zip|CSTR MHE Solution in Simulink]]

Attach:download.png [[Attach:cstr_mhe_solution.zip|CSTR MHE Solution in Simulink]]

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 ''T'_c_'''. See a [[https://apmonitor.com/che436/index.php/Main/CaseStudyCSTR|related CSTR case study]] for details on the model.

Changed line 23 from:

Attach:download.png [[Attach:cstr_mhe_solution.zip|CSTR MHE Solution]]

to:

Attach:download.png [[Attach:cstr_mhe_solution.zip|CSTR MHE Solution in Simulink]]

Changed lines 21-23 from:

!!!! ~~Solution~~

to:

!!!! Solution

Attach:download.png [[Attach:cstr_mhe_solution.zip|CSTR MHE Solution]]

Attach:download.png [[Attach:cstr_mhe_solution.zip|CSTR MHE Solution]]

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 ''T'_c_'''.

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'''.

Added lines 18-19:

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

Changed line 13 from:

Attach:download.png [[Attach:cstr_~~control~~.zip|CSTR Source Files]]

to:

Attach:download.png [[Attach:cstr_estimation.zip|CSTR Source Files]]

Added lines 12-13:

Attach:download.png [[Attach:cstr_control.zip|CSTR Source Files]]

Added lines 9-17:

!!!! 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.''

Attach:cstr.png

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

'''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.''

Attach:cstr.png

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

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.

Added lines 4-5:

Dynamic models constructed with equations that describe physical phenomena may need to be tuned by adjusting parameters so that predicted outputs match with experimental data. Physical models are based on the underlying physical principles that govern the problem and result from expressions such as a force or momentum balance and may include quantities such as velocity, acceleration, and position. Other quantities of interest may include anything that changes with respect to time such as reactor composition, temperature, mole fraction, etc. Mathematical models likely contain both physical and experimental elements. This section shows how to reconcile experimental data with the physical model through parameter estimation.

Added lines 1-5:

(:title Moving Horizon Estimation:)

(:keywords Kalman filter, Simulink, moving horizon, time window, dynamic data, validation, estimation, differential, algebraic, tutorial:)

(:description Dynamic state and parameter estimation with Moving Horizon Estimation:)

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.

(:keywords Kalman filter, Simulink, moving horizon, time window, dynamic data, validation, estimation, differential, algebraic, tutorial:)

(:description Dynamic state and parameter estimation with Moving Horizon Estimation:)

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.