Dyad

Appearance

LIBRARY
Examples.PHNeutralizationMPCExample.md

Examples.PHNeutralizationMPCExample ​

Closed-loop pH control of the neutralization reactor of the ModelPredictiveControl.jl DAE tutorial ("Nonlinear Design (DAE)"), with the DAE kept as such. The plant PHNeutralizationPlant is also the prediction model of the controller and the model of the estimator (dynamics = plant). ACADOSMPC manipulates the base inflow q_Bin to track a pH reference; it integrates the model by a 4-stage Gauss-Legendre IRK scheme (the tutorial transcribes it by trapezoidal collocation), and the pH is the algebraic variable of the model, which the cost weights through acados' algebraic output selection. The acid inflow q_Ain is a measured disturbance, which enters the controller as an online parameter (known_disturbance_inputs, its value on the p input). The acid feed concentration is an unmeasured disturbance, represented by the state variable d of the plant's integrating disturbance model.

In place of the adaptive moving-horizon estimator of the tutorial, a UKFObserver estimates the two concentrations and d, the latter driven by the process noise, with the tutorial's standard deviations (σQ = [0.0015, 0.0015, 2e-4], σR = 0.05, σP_0 = [0.05, 0.05, 5e-4]) and initial estimate [0.025, 0.025, 0], and outputs the estimated pH as the algebraic variable at its estimate. As in the tutorial, the pH is measured with a noise of standard deviation 0.05 and the acid inflow with one of 0.1 L/min, the sample time is half an hour, and the plant starts from c_A = 0.0505, c_B = 0.0495.

The scenario follows the three disturbances of the tutorial over 40 h, each at a quarter of the horizon: at 10 h the pH reference steps from r0 to r1 (where the tutorial bumps the base inflow in open loop), at 20 h the acid inflow steps from 10 to 12 L/min, and at 30 h the acid feed concentration drops from 0.1 to 0.09 mol/L. The increment penalty of ACADOSMPC is not available for a DAE prediction model, so the control is weighted by the small Q2, whose steady-state offset in the pH is of the order of Q2 u / (∂pH/∂u).

The solver starts every tick from the previous solution as it is (warm_start = Keep), not shifted by one stage. The cost weights the pH at the start of each shooting interval, and the terminal stage carries no algebraic variable, so the last control of a solution influences only the unweighted terminal state and is driven to its lower bound by Q2. A shift moves that control into the second-to-last interval, where the full Gauss-Newton steps that correct it cross the steep part of the titration curve and the integration diverges; with the shift, most ticks needed the retry of reset_on_failure. Verified in test/test_acados_tutorials.jl.

Usage ​

MPCComponents.Examples.PHNeutralizationMPCExample(r0=5.0, r1=5.5)

Parameters: ​

NameDescriptionUnitsDefault value
TsSample time [h]–0.5
NpPrediction horizon (intervals)–20
r0pH reference before the reference step–5.0
r1pH reference after the reference step–5.5

Variables ​

NameDescriptionUnits
pHpH of the plant–
pH_measMeasured pH–
pH_hatEstimated pH–
rpH reference–
c_ATotal concentration of the acid species [mol/L]–
c_BConcentration of the base [mol/L]–
c_AinAcid feed concentration [mol/L]–
c_A_hatEstimated total concentration of the acid species [mol/L]–
c_B_hatEstimated concentration of the base [mol/L]–
d_hatEstimated state of the disturbance model of the acid feed concentration [mol/L]–
c_Ain_hatEstimated acid feed concentration [mol/L]–
q_BinBase inflow, the control [L/min]–
q_Ain_measMeasured acid inflow [L/min]–
flagacados status–

Behavior ​

Behavior of this component cannot be rendered because it includes path variables.

Source ​

dyad
"""
Closed-loop pH control of the neutralization reactor of the ModelPredictiveControl.jl DAE tutorial
("Nonlinear Design (DAE)"), with the DAE kept as such. The plant `PHNeutralizationPlant` is also the
prediction model of the controller and the model of the estimator (`dynamics = plant`). `ACADOSMPC`
manipulates the base inflow `q_Bin` to track a pH reference; it integrates the model by a 4-stage
Gauss-Legendre IRK scheme (the tutorial transcribes it by trapezoidal collocation), and the pH is the
algebraic variable of the model, which the cost weights through acados' algebraic output selection.
The acid inflow `q_Ain` is a measured disturbance, which enters the controller as an online parameter
(`known_disturbance_inputs`, its value on the `p` input). The acid feed concentration is an unmeasured
disturbance, represented by the state variable `d` of the plant's integrating disturbance model.

In place of the adaptive moving-horizon estimator of the tutorial, a `UKFObserver` estimates the two
concentrations and `d`, the latter driven by the process noise, with the tutorial's standard
deviations (`σQ = [0.0015, 0.0015, 2e-4]`, `σR = 0.05`, `σP_0 = [0.05, 0.05, 5e-4]`) and initial
estimate `[0.025, 0.025, 0]`, and outputs the estimated pH as the algebraic variable at its estimate. As in
the tutorial, the pH is measured with a noise of standard deviation 0.05 and the acid inflow with one
of 0.1 L/min, the sample time is half an hour, and the plant starts from `c_A = 0.0505`,
`c_B = 0.0495`.

The scenario follows the three disturbances of the tutorial over 40 h, each at a quarter of the
horizon: at 10 h the pH reference steps from `r0` to `r1` (where the tutorial bumps the base inflow
in open loop), at 20 h the acid inflow steps from 10 to 12 L/min, and at 30 h the acid feed
concentration drops from 0.1 to 0.09 mol/L. The increment penalty of `ACADOSMPC` is not available for
a DAE prediction model, so the control is weighted by the small `Q2`, whose steady-state offset in
the pH is of the order of `Q2 u / (∂pH/∂u)`.

The solver starts every tick from the previous solution as it is (`warm_start = Keep`), not shifted by
one stage. The cost weights the pH at the start of each shooting interval, and the terminal stage
carries no algebraic variable, so the last control of a solution influences only the unweighted
terminal state and is driven to its lower bound by `Q2`. A shift moves that control into the
second-to-last interval, where the full Gauss-Newton steps that correct it cross the steep part of
the titration curve and the integration diverges; with the shift, most ticks needed the retry of
`reset_on_failure`. Verified in `test/test_acados_tutorials.jl`.
"""
example component PHNeutralizationMPCExample
  "Sample time [h]"
  structural parameter Ts::Real = 0.5
  "Prediction horizon (intervals)"
  structural parameter Np::Integer = 20
  "pH reference before the reference step"
  parameter r0::Real = 5.0
  "pH reference after the reference step"
  parameter r1::Real = 5.5
  periodicclock = DiscreteComponents.PeriodicClock(dt = Ts) {
    "Dyad": {"placement": {"diagram": {"x1": 150, "y1": 900, "x2": 250, "y2": 1000}}}
  }
  plant = MPCComponents.Examples.PHNeutralizationPlant(c_A0 = 0.0505, c_B0 = 0.0495) {
    "Dyad": {"placement": {"diagram": {"x1": 950, "y1": 250, "x2": 1050, "y2": 350}}}
  }
  "pH reference"
  ref = BlockComponents.Sources.Step(final offset = r0, final height = r1 - r0, start_time = 10.0) {"Dyad": {"placement": {"diagram": {"x1": 0, "y1": 50, "x2": 100, "y2": 150}}}}
  "Acid inflow, the measured disturbance [L/min]"
  q_Ain_src = BlockComponents.Sources.Step(offset = 10.0, height = 2.0, start_time = 20.0) {
    "Dyad": {"placement": {"diagram": {"x1": 750, "y1": 400, "x2": 850, "y2": 500}}}
  }
  "Deviation of the acid feed concentration, the unmeasured disturbance [mol/L]"
  c_Ain_src = BlockComponents.Sources.Step(offset = 0.0, height = -0.01, start_time = 30.0) {
    "Dyad": {"placement": {"diagram": {"x1": 750, "y1": 550, "x2": 850, "y2": 650}}}
  }
  "Input of the disturbance model of the plant, zero in simulation"
  zero = BlockComponents.Sources.Constant(k = 0.0) {
    "Dyad": {"placement": {"diagram": {"x1": 750, "y1": 700, "x2": 850, "y2": 800}}}
  }
  "Sampler of the pH"
  spH = DiscreteComponents.Sampler() {
    "Dyad": {"placement": {"diagram": {"x1": 1050, "y1": 750, "x2": 950, "y2": 850}}}
  }
  "Sampler of the acid inflow"
  sq = DiscreteComponents.Sampler() {
    "Dyad": {"placement": {"diagram": {"x1": 700, "y1": 400, "x2": 600, "y2": 500}}}
  }
  "Sampler of the pH reference"
  sr = DiscreteComponents.Sampler() {
    "Dyad": {"placement": {"diagram": {"x1": 150, "y1": 50, "x2": 250, "y2": 150}}}
  }
  "Measurement noise of the pH"
  noise_pH = DiscreteComponents.NormalNoise(sigma = 0.05, seed = 1) {
    "Dyad": {"placement": {"diagram": {"x1": 1050, "y1": 900, "x2": 950, "y2": 1000}}}
  }
  "Measurement noise of the acid inflow [L/min]"
  noise_q = DiscreteComponents.NormalNoise(sigma = 0.1, seed = 2) {
    "Dyad": {"placement": {"diagram": {"x1": 700, "y1": 550, "x2": 600, "y2": 650}}}
  }
  "Measured pH"
  add_pH = BlockComponents.Math.Add() {
    "Dyad": {"placement": {"diagram": {"x1": 850, "y1": 750, "x2": 750, "y2": 850}}}
  }
  "Measured acid inflow"
  add_q = BlockComponents.Math.Add() {
    "Dyad": {"placement": {"diagram": {"x1": 550, "y1": 400, "x2": 450, "y2": 500}}}
  }
  "Hold of the base inflow"
  zoh = DiscreteComponents.ZeroOrderHold() {
    "Dyad": {"placement": {"diagram": {"x1": 600, "y1": 200, "x2": 700, "y2": 300}}}
  }
  ukf = MPCComponents.UKFObserver(dynamics = plant, inputs = ["q_Bin_in", "q_Ain_in"], disturbance_inputs = ["dc_Ain_in", "w"], outputs = ["pH"], state_variables = ["c_A", "c_B", "d"], final Ts = Ts, R1 = diagonal([2.25e-6, 2.25e-6, 4e-8]), R2 = diagonal([0.0025]), x0 = [0.025, 0.025, 0.0], R0 = diagonal([0.0025, 0.0025, 2.5e-7])) {
    "Dyad": {"placement": {"diagram": {"x1": 500, "y1": 750, "x2": 400, "y2": 850}}}
  }
  mpc = MPCComponents.ACADOSMPC(dynamics = plant, inputs = ["q_Bin_in"], known_disturbance_inputs = ["q_Ain_in"], disturbance_inputs = ["dc_Ain_in", "w"], state_variables = ["c_A", "c_B", "d"], outputs = ["pH"], final Ts = Ts, final Np = Np, Q1 = diagonal([1.0]), Q2 = diagonal([1e-4]), umin = [0.0], umax = [30.0], u_guess = [10.0], x_scale = [0.05, 0.05, 0.01], u_scale = [10.0], integrator = MPCComponents.ACADOSIntegrator.IRK(), integrator_stages = 4, integrator_steps = 5, integrator_jac_reuse = false, warm_start = MPCComponents.ACADOSWarmStart.Keep(), reset_on_failure = true, x_guess = [0.0505, 0.0495, 0.0]) {
    "Dyad": {"placement": {"diagram": {"x1": 350, "y1": 200, "x2": 450, "y2": 300}}}
  }
  "pH of the plant"
  variable pH::Real
  "Measured pH"
  variable pH_meas::Real
  "Estimated pH"
  variable pH_hat::Real
  "pH reference"
  variable r::Real
  "Total concentration of the acid species [mol/L]"
  variable c_A::Real
  "Concentration of the base [mol/L]"
  variable c_B::Real
  "Acid feed concentration [mol/L]"
  variable c_Ain::Real
  "Estimated total concentration of the acid species [mol/L]"
  variable c_A_hat::Real
  "Estimated concentration of the base [mol/L]"
  variable c_B_hat::Real
  "Estimated state of the disturbance model of the acid feed concentration [mol/L]"
  variable d_hat::Real
  "Estimated acid feed concentration [mol/L]"
  variable c_Ain_hat::Real
  "Base inflow, the control [L/min]"
  variable q_Bin::Real
  "Measured acid inflow [L/min]"
  variable q_Ain_meas::Real
  "acados status"
  variable flag::Real
relations
  pH = plant.pH
  c_A = plant.c_A
  c_B = plant.c_B
  c_Ain = plant.c_Ain
  r = ref.y
  pH_meas = add_pH.y
  pH_hat = ukf.zhat[1]
  c_A_hat = ukf.xhat[1]
  c_B_hat = ukf.xhat[2]
  d_hat = ukf.xhat[3]
  c_Ain_hat = plant.c_Ain0 + d_hat
  q_Bin = mpc.u[1]
  q_Ain_meas = add_q.y
  flag = mpc.exitflag
  connect(plant.pH_out, spH.u) {}
  connect(spH.y, add_pH.u1) {}
  connect(noise_pH.y, add_pH.u2) {}
  connect(add_pH.y, ukf.y[1], periodicclock.y) {}
  connect(q_Ain_src.y, plant.q_Ain_in, sq.u) {}
  connect(sq.y, add_q.u1) {}
  connect(noise_q.y, add_q.u2) {}
  connect(add_q.y, ukf.u[2], mpc.p[1]) {}
  connect(c_Ain_src.y, plant.dc_Ain_in) {}
  connect(zero.y, plant.w) {}
  connect(ref.y, sr.u) {}
  connect(sr.y, mpc.r[1]) {}
  connect(ukf.xhat[1], mpc.x[1]) {}
  connect(ukf.xhat[2], mpc.x[2]) {}
  connect(ukf.xhat[3], mpc.x[3]) {}
  connect(mpc.u[1], zoh.u, ukf.u[1]) {}
  connect(zoh.y, plant.q_Bin_in) {}
metadata {"Dyad": {"tests": {"case1": {"stop": 40}}, "doc": {"behavior": false}}}
end
Flattened Source
dyad
"""
Closed-loop pH control of the neutralization reactor of the ModelPredictiveControl.jl DAE tutorial
("Nonlinear Design (DAE)"), with the DAE kept as such. The plant `PHNeutralizationPlant` is also the
prediction model of the controller and the model of the estimator (`dynamics = plant`). `ACADOSMPC`
manipulates the base inflow `q_Bin` to track a pH reference; it integrates the model by a 4-stage
Gauss-Legendre IRK scheme (the tutorial transcribes it by trapezoidal collocation), and the pH is the
algebraic variable of the model, which the cost weights through acados' algebraic output selection.
The acid inflow `q_Ain` is a measured disturbance, which enters the controller as an online parameter
(`known_disturbance_inputs`, its value on the `p` input). The acid feed concentration is an unmeasured
disturbance, represented by the state variable `d` of the plant's integrating disturbance model.

In place of the adaptive moving-horizon estimator of the tutorial, a `UKFObserver` estimates the two
concentrations and `d`, the latter driven by the process noise, with the tutorial's standard
deviations (`σQ = [0.0015, 0.0015, 2e-4]`, `σR = 0.05`, `σP_0 = [0.05, 0.05, 5e-4]`) and initial
estimate `[0.025, 0.025, 0]`, and outputs the estimated pH as the algebraic variable at its estimate. As in
the tutorial, the pH is measured with a noise of standard deviation 0.05 and the acid inflow with one
of 0.1 L/min, the sample time is half an hour, and the plant starts from `c_A = 0.0505`,
`c_B = 0.0495`.

The scenario follows the three disturbances of the tutorial over 40 h, each at a quarter of the
horizon: at 10 h the pH reference steps from `r0` to `r1` (where the tutorial bumps the base inflow
in open loop), at 20 h the acid inflow steps from 10 to 12 L/min, and at 30 h the acid feed
concentration drops from 0.1 to 0.09 mol/L. The increment penalty of `ACADOSMPC` is not available for
a DAE prediction model, so the control is weighted by the small `Q2`, whose steady-state offset in
the pH is of the order of `Q2 u / (∂pH/∂u)`.

The solver starts every tick from the previous solution as it is (`warm_start = Keep`), not shifted by
one stage. The cost weights the pH at the start of each shooting interval, and the terminal stage
carries no algebraic variable, so the last control of a solution influences only the unweighted
terminal state and is driven to its lower bound by `Q2`. A shift moves that control into the
second-to-last interval, where the full Gauss-Newton steps that correct it cross the steep part of
the titration curve and the integration diverges; with the shift, most ticks needed the retry of
`reset_on_failure`. Verified in `test/test_acados_tutorials.jl`.
"""
example component PHNeutralizationMPCExample
  "Sample time [h]"
  structural parameter Ts::Real = 0.5
  "Prediction horizon (intervals)"
  structural parameter Np::Integer = 20
  "pH reference before the reference step"
  parameter r0::Real = 5.0
  "pH reference after the reference step"
  parameter r1::Real = 5.5
  periodicclock = DiscreteComponents.PeriodicClock(dt = Ts) {
    "Dyad": {"placement": {"diagram": {"x1": 150, "y1": 900, "x2": 250, "y2": 1000}}}
  }
  plant = MPCComponents.Examples.PHNeutralizationPlant(c_A0 = 0.0505, c_B0 = 0.0495) {
    "Dyad": {"placement": {"diagram": {"x1": 950, "y1": 250, "x2": 1050, "y2": 350}}}
  }
  "pH reference"
  ref = BlockComponents.Sources.Step(final offset = r0, final height = r1 - r0, start_time = 10.0) {"Dyad": {"placement": {"diagram": {"x1": 0, "y1": 50, "x2": 100, "y2": 150}}}}
  "Acid inflow, the measured disturbance [L/min]"
  q_Ain_src = BlockComponents.Sources.Step(offset = 10.0, height = 2.0, start_time = 20.0) {
    "Dyad": {"placement": {"diagram": {"x1": 750, "y1": 400, "x2": 850, "y2": 500}}}
  }
  "Deviation of the acid feed concentration, the unmeasured disturbance [mol/L]"
  c_Ain_src = BlockComponents.Sources.Step(offset = 0.0, height = -0.01, start_time = 30.0) {
    "Dyad": {"placement": {"diagram": {"x1": 750, "y1": 550, "x2": 850, "y2": 650}}}
  }
  "Input of the disturbance model of the plant, zero in simulation"
  zero = BlockComponents.Sources.Constant(k = 0.0) {
    "Dyad": {"placement": {"diagram": {"x1": 750, "y1": 700, "x2": 850, "y2": 800}}}
  }
  "Sampler of the pH"
  spH = DiscreteComponents.Sampler() {
    "Dyad": {"placement": {"diagram": {"x1": 1050, "y1": 750, "x2": 950, "y2": 850}}}
  }
  "Sampler of the acid inflow"
  sq = DiscreteComponents.Sampler() {
    "Dyad": {"placement": {"diagram": {"x1": 700, "y1": 400, "x2": 600, "y2": 500}}}
  }
  "Sampler of the pH reference"
  sr = DiscreteComponents.Sampler() {
    "Dyad": {"placement": {"diagram": {"x1": 150, "y1": 50, "x2": 250, "y2": 150}}}
  }
  "Measurement noise of the pH"
  noise_pH = DiscreteComponents.NormalNoise(sigma = 0.05, seed = 1) {
    "Dyad": {"placement": {"diagram": {"x1": 1050, "y1": 900, "x2": 950, "y2": 1000}}}
  }
  "Measurement noise of the acid inflow [L/min]"
  noise_q = DiscreteComponents.NormalNoise(sigma = 0.1, seed = 2) {
    "Dyad": {"placement": {"diagram": {"x1": 700, "y1": 550, "x2": 600, "y2": 650}}}
  }
  "Measured pH"
  add_pH = BlockComponents.Math.Add() {
    "Dyad": {"placement": {"diagram": {"x1": 850, "y1": 750, "x2": 750, "y2": 850}}}
  }
  "Measured acid inflow"
  add_q = BlockComponents.Math.Add() {
    "Dyad": {"placement": {"diagram": {"x1": 550, "y1": 400, "x2": 450, "y2": 500}}}
  }
  "Hold of the base inflow"
  zoh = DiscreteComponents.ZeroOrderHold() {
    "Dyad": {"placement": {"diagram": {"x1": 600, "y1": 200, "x2": 700, "y2": 300}}}
  }
  ukf = MPCComponents.UKFObserver(dynamics = plant, inputs = ["q_Bin_in", "q_Ain_in"], disturbance_inputs = ["dc_Ain_in", "w"], outputs = ["pH"], state_variables = ["c_A", "c_B", "d"], final Ts = Ts, R1 = diagonal([2.25e-6, 2.25e-6, 4e-8]), R2 = diagonal([0.0025]), x0 = [0.025, 0.025, 0.0], R0 = diagonal([0.0025, 0.0025, 2.5e-7])) {
    "Dyad": {"placement": {"diagram": {"x1": 500, "y1": 750, "x2": 400, "y2": 850}}}
  }
  mpc = MPCComponents.ACADOSMPC(dynamics = plant, inputs = ["q_Bin_in"], known_disturbance_inputs = ["q_Ain_in"], disturbance_inputs = ["dc_Ain_in", "w"], state_variables = ["c_A", "c_B", "d"], outputs = ["pH"], final Ts = Ts, final Np = Np, Q1 = diagonal([1.0]), Q2 = diagonal([1e-4]), umin = [0.0], umax = [30.0], u_guess = [10.0], x_scale = [0.05, 0.05, 0.01], u_scale = [10.0], integrator = MPCComponents.ACADOSIntegrator.IRK(), integrator_stages = 4, integrator_steps = 5, integrator_jac_reuse = false, warm_start = MPCComponents.ACADOSWarmStart.Keep(), reset_on_failure = true, x_guess = [0.0505, 0.0495, 0.0]) {
    "Dyad": {"placement": {"diagram": {"x1": 350, "y1": 200, "x2": 450, "y2": 300}}}
  }
  "pH of the plant"
  variable pH::Real
  "Measured pH"
  variable pH_meas::Real
  "Estimated pH"
  variable pH_hat::Real
  "pH reference"
  variable r::Real
  "Total concentration of the acid species [mol/L]"
  variable c_A::Real
  "Concentration of the base [mol/L]"
  variable c_B::Real
  "Acid feed concentration [mol/L]"
  variable c_Ain::Real
  "Estimated total concentration of the acid species [mol/L]"
  variable c_A_hat::Real
  "Estimated concentration of the base [mol/L]"
  variable c_B_hat::Real
  "Estimated state of the disturbance model of the acid feed concentration [mol/L]"
  variable d_hat::Real
  "Estimated acid feed concentration [mol/L]"
  variable c_Ain_hat::Real
  "Base inflow, the control [L/min]"
  variable q_Bin::Real
  "Measured acid inflow [L/min]"
  variable q_Ain_meas::Real
  "acados status"
  variable flag::Real
relations
  pH = plant.pH
  c_A = plant.c_A
  c_B = plant.c_B
  c_Ain = plant.c_Ain
  r = ref.y
  pH_meas = add_pH.y
  pH_hat = ukf.zhat[1]
  c_A_hat = ukf.xhat[1]
  c_B_hat = ukf.xhat[2]
  d_hat = ukf.xhat[3]
  c_Ain_hat = plant.c_Ain0 + d_hat
  q_Bin = mpc.u[1]
  q_Ain_meas = add_q.y
  flag = mpc.exitflag
  connect(plant.pH_out, spH.u) {}
  connect(spH.y, add_pH.u1) {}
  connect(noise_pH.y, add_pH.u2) {}
  connect(add_pH.y, ukf.y[1], periodicclock.y) {}
  connect(q_Ain_src.y, plant.q_Ain_in, sq.u) {}
  connect(sq.y, add_q.u1) {}
  connect(noise_q.y, add_q.u2) {}
  connect(add_q.y, ukf.u[2], mpc.p[1]) {}
  connect(c_Ain_src.y, plant.dc_Ain_in) {}
  connect(zero.y, plant.w) {}
  connect(ref.y, sr.u) {}
  connect(sr.y, mpc.r[1]) {}
  connect(ukf.xhat[1], mpc.x[1]) {}
  connect(ukf.xhat[2], mpc.x[2]) {}
  connect(ukf.xhat[3], mpc.x[3]) {}
  connect(mpc.u[1], zoh.u, ukf.u[1]) {}
  connect(zoh.y, plant.q_Bin_in) {}
metadata {"Dyad": {"tests": {"case1": {"stop": 40}}, "doc": {"behavior": false}}}
end


Test Cases ​

Test Case case1 ​

  • Examples

  • Experiments

  • Analyses