Moving Horizon Estimator¶
This example is meant as a quickstart with the Moving Horizon Estimator implemented in HILO-MPC. We will use a model of a chemical reaction. We will se how:
To define a model of a chemical reaction system
Use a Moving Horizon Estimator (MHE)
The example is taken from the book [1]. For a more in-depth introduction to MHE have a look at:
let’s get started now.
Introduction¶
The model describe a gas-phase chemical reaction with three reactants.
\[\begin{split}\begin{align} \dot{C_a} &= - (k_1 C_a - k_{-1} C_b C_c) \\ \dot{C_b} &= k_1 C_a - k_{-1} C_b C_c - 2 (k_2 C_b^2 - k_{-2} C_c)\\ \dot{C_c} &= k_1 C_a - k_{-1} C_b C_c + (k_2 C_b^2 - k_{-2} C_c) \end{align}\end{split}\]we measure the total pressure of the vessel
\[y = \frac{R T}{V} (C_a + C_b + C_c)\]
Model¶
[1]:
import sys # You don't need this if you installed HILO-MPC using pip!
sys.path.append('../../../') # You don't need this if you installed HILO-MPC using pip!
from hilo_mpc import Model, MHE
from bokeh.io import output_notebook, show
from bokeh.plotting import figure
from bokeh.layouts import gridplot
import numpy as np
import warnings
warnings.filterwarnings("ignore")
[2]:
# Create model
model = Model(plot_backend='bokeh')
x = model.set_dynamical_states(['Ca', 'Cb', 'Cc'], units=['mol/l', 'mol/l', 'mol/l'], short_description=['Ca', 'Cb', 'Cc'])
model.set_measurements(['P'], units=['atm'], short_description=['Pressure'])
# Unwrap states
Ca = x[0]
Cb = x[1]
Cc = x[2]
# Known Parameters
k1 = 0.5
k_1 = 0.05
k2 = 0.2
k_2 = 0.01
dt = 0.25
RT = 32.84 # L atm/ (mol)
dCa = - (k1 * Ca - k_1 * Cb * Cc)
dCb = k1 * Ca - k_1 * Cb * Cc - 2 * (k2 * Cb ** 2 - k_2 * Cc)
dCc = k1 * Ca - k_1 * Cb * Cc + 1 * (k2 * Cb ** 2 - k_2 * Cc)
model.set_measurement_equations(RT * (Ca + Cb + Cc))
model.set_dynamical_equations([dCa, dCb, dCc])
model.setup(dt=dt)
# Initial conditions
x0_real = [0.5, 0.05, 0]
x0_est = [1, 0, 4]
model.set_initial_conditions(x0=x0_real)
n_steps = 120
Moving Horizon Estimator¶
[3]:
# Setup the MHE
mhe = MHE(model)
mhe.quad_arrival_cost.add_states(weights=[1/(0.5**2), 1/(0.5**2), 1/(0.5**2)], guess=x0_est)
mhe.quad_stage_cost.add_measurements(weights=[1/(0.25**2)])
mhe.quad_stage_cost.add_state_noise(weights=[1/(0.001**2), 1/(0.001**2), 1/(0.001**2)])
mhe.set_box_constraints(x_lb=[0, 0, 0])
mhe.horizon = 20
mhe.setup(options={'print_level': 0})
# Run the simulation
for i in range(n_steps):
model.simulate()
mhe.add_measurements(y_meas=model.solution['y'][:, -2])
x_est, _ = mhe.estimate()
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************
Plotting¶
[5]:
output_notebook()
p_tot = []
p = figure(background_fill_color="#fafafa", width=300, height=300)
p.scatter(x=np.array(mhe.solution['t']).squeeze(), y=np.array(mhe.solution['Ca']).squeeze(), legend_label='Estimated')
p.line(x=np.array(model.solution['t']).squeeze(), y=np.array(model.solution['Ca']).squeeze(), legend_label='Real')
p_tot.append(p)
p = figure(background_fill_color="#fafafa", width=300, height=300)
p.scatter(x=np.array(mhe.solution['t']).squeeze(), y=np.array(mhe.solution['Cb']).squeeze(), legend_label='Estimated')
p.line(x=np.array(model.solution['t']).squeeze(), y=np.array(model.solution['Cb']).squeeze(), legend_label='Real')
p_tot.append(p)
p = figure(background_fill_color="#fafafa", width=300, height=300)
p.scatter(x=np.array(mhe.solution['t']).squeeze(), y=np.array(mhe.solution['Cc']).squeeze(), legend_label='Estimated')
p.line(x=np.array(model.solution['t']).squeeze(), y=np.array(model.solution['Cc']).squeeze(), legend_label='Real')
p_tot.append(p)
show(gridplot(p_tot, ncols=2))