From e02f51230204722d95fb8cdea4f5aaea873eb9cc Mon Sep 17 00:00:00 2001 From: franckgaga Date: Wed, 23 Sep 2026 08:57:45 -0400 Subject: [PATCH 01/13] doc: minor correction --- src/estimator/mhe/construct.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/estimator/mhe/construct.jl b/src/estimator/mhe/construct.jl index b0ade0423..492b79f3d 100644 --- a/src/estimator/mhe/construct.jl +++ b/src/estimator/mhe/construct.jl @@ -259,7 +259,7 @@ end Construct a moving horizon estimator (MHE) based on `model`. -It supports ([`LinModel`](@ref), [`NonLinModel`](@ref), [`NonLinModelDAE`](@ref)) and +It supports [`LinModel`](@ref), [`NonLinModel`](@ref), [`NonLinModelDAE`](@ref) and constraints on the estimates. Additionally, `model` is not linearized like the [`ExtendedKalmanFilter`](@ref), and the probability distribution is not approximated like the [`UnscentedKalmanFilter`](@ref). The computational costs are drastically higher, From 2b324be0a66e530c05f913e4849317d5b2d1c8f1 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Wed, 23 Sep 2026 09:22:31 -0400 Subject: [PATCH 02/13] test: new MHE construction with DEA and default transcription --- test/2_test_state_estim.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/2_test_state_estim.jl b/test/2_test_state_estim.jl index d17376c63..6cb418212 100644 --- a/test/2_test_state_estim.jl +++ b/test/2_test_state_estim.jl @@ -1059,6 +1059,9 @@ end mhe2 = MovingHorizonEstimator(dae; He=3, transcription, direct=false) @test mhe2.direct == false + mhe3 = MovingHorizonEstimator(dae, He=3, hessian=true) + @test mhe.transcription isa OrthogonalCollocation + @test_throws ArgumentError MovingHorizonEstimator(dae, He=3, transcription=SingleShooting()) @test_throws ArgumentError MovingHorizonEstimator(dae, He=3, transcription=MultipleShooting()) end From 6f780ebf0b5bc28a0f0e804970e33dd703f947e8 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Wed, 23 Sep 2026 09:27:53 -0400 Subject: [PATCH 03/13] test: debug test --- test/2_test_state_estim.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/2_test_state_estim.jl b/test/2_test_state_estim.jl index 6cb418212..620f254e0 100644 --- a/test/2_test_state_estim.jl +++ b/test/2_test_state_estim.jl @@ -1060,7 +1060,7 @@ end @test mhe2.direct == false mhe3 = MovingHorizonEstimator(dae, He=3, hessian=true) - @test mhe.transcription isa OrthogonalCollocation + @test mhe3.transcription isa OrthogonalCollocation @test_throws ArgumentError MovingHorizonEstimator(dae, He=3, transcription=SingleShooting()) @test_throws ArgumentError MovingHorizonEstimator(dae, He=3, transcription=MultipleShooting()) From 1f4906da0dfe8be9dc27a86d08d54b06481be037 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Wed, 23 Sep 2026 12:56:48 -0400 Subject: [PATCH 04/13] doc: wip pH neutralization example --- docs/src/manual/nonlinmpc2.md | 155 ++++++++++++++++++++++++++-------- 1 file changed, 119 insertions(+), 36 deletions(-) diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md index 63a724903..e5fb5884d 100644 --- a/docs/src/manual/nonlinmpc2.md +++ b/docs/src/manual/nonlinmpc2.md @@ -5,8 +5,8 @@ Pages = ["nonlinmpc2.md"] ``` !!! todo "Under Construction" - This tutorial is currently under construction. Only the modeling part is written for - now. + This tutorial is currently under construction. Only the modeling and the estimation + parts are written for now. ## Nonlinear Model (DAE) @@ -89,8 +89,9 @@ The pH is computed with: 0 = a_H^3 + (c_B + K_a) a_H^2 + \big(K_a(c_A + c_B) + K_w \big) a_H - K_w K_a ``` We could extract the positive real root of this expression inside the output function - `h!` to transform the system to an ODE, effectively avoiding the increased complexity - of DAEs. The tutorial will still treat the system as a DAE to illustrate its API. + `h!` to transform the system to an ODE, effectively avoiding the complexity of DAEs. + When possible, plant model should be constructed with the specialized [`NonLinModel`](@ref) + for ODEs. This tutorial will still treat the system as a DAE to illustrate its API. ### Mass Balance @@ -117,17 +118,17 @@ because of the weir, the following relations compute the outflow terms: \end{aligned} ``` -The code is: +The [`NonLinModelDAE`](@ref) constructor expects that the state dynamics and the algebraic +equation is combined into a single `fq!(ẋ, res, x, a, u, d, p) -> nothing` function that +modifies both `ẋ` and `res` arguments in-place (an out-of-place option is also available), +with the state dynamics and the residual of the algebraic equation, respectively: ```@example 1 using ModelPredictiveControl -V = 1000.0 # reactor volume [L] -c_Ain = 0.1 # feed concentration of weak acid [mol/L] -c_Bin = 0.1 # feed concentration of strong base [mol/L] -Kw = 1.0e-14 # water dissociation constant [mol^2/L^2] -Ka = 1.75e-5 # acid dissociation constant [mol/L] - +calc_ċ_A(c_Ain, q_Ain, c_Aout, q_out) = (60/V)*(q_Ain * c_Ain - q_out * c_Aout) +calc_ċ_B(c_Bin, q_Bin, c_Bout, q_out) = (60/V)*(q_Bin * c_Bin - q_out * c_Bout) +calc_res(a_H, c_A, c_B, Kw, Ka) = a_H + c_B - (Kw / a_H) - (Ka * c_A / (Ka + a_H)) function fq!(ẋ, res, x, a, u, d, p) c_Ain, c_Bin, Kw, Ka, V = p q_Ain, q_Bin = d[1], u[1] # [L/min], [L/min] @@ -136,58 +137,140 @@ function fq!(ẋ, res, x, a, u, d, p) q_out = q_Ain + q_Bin # [L/min] c_Aout = c_A # [mol/L] c_Bout = c_B # [mol/L] - ẋ[1] = (60/V)*(q_Ain * c_Ain - q_out * c_Aout) - ẋ[2] = (60/V)*(q_Bin * c_Bin - q_out * c_Bout) - res[1] = a_H + c_B - (Kw / a_H) - (Ka * c_A / (Ka + a_H)) + ẋ[1] = calc_ċ_A(c_Ain, q_Ain, c_Aout, q_out) + ẋ[2] = calc_ċ_B(c_Bin, q_Bin, c_Bout, q_out) + res[1] = calc_res(a_H, c_A, c_B, Kw, Ka) return nothing end +``` + +A similar in-place function is expected for the model output: -function h!(y, _, a, _ , _ ) +```@example 1 +calc_pH(a_H) = -log10(a_H) +function h!(y, _ , a, _ , _ ) a_H = a[1] pH = try - -log10(a_H) + calc_pH(a_H) catch myerror myerror isa DomainError ? NaN : rethrow() end y[1] = pH return nothing end +``` -Ts = 0.5 # Sample time [h] +By default, Julia throws a `DomainError` if `log10` is called with a negative number. The +`try` blocks is necessary to let the optimizer explores undefined domains by a returning a +`NaN` value in such cases. For similar reasons, providing an initial guess for the algebraic +variable `as_0` is crucial here to prioritize positive ``a_H`` concentration and ensure a +defined ``\mathrm{pH}`` value: + +```@example 1 +V = 1000.0 # reactor volume [L] +c_Ain = 0.1 # feed concentration of weak acid [mol/L] +c_Bin = 0.1 # feed concentration of strong base [mol/L] +Kw = 1.0e-14 # water dissociation constant [mol^2/L^2] +Ka = 1.75e-5 # acid dissociation constant [mol/L] + +Ts = 0.5 # Sample time [h] nu, nx, na, ny, nd = 1, 2, 1, 1, 1 p = [c_Ain, c_Bin, Kw, Ka, V] +vu, vd = [raw"$q_{Bin}$ (L/min)"], [raw"$q_{Ain}$ (L/min)"] +vx, vy = [raw"$c_A$ (mol/L)", raw"$c_B$ (mol/L)"], [raw"$\mathrm{pH}$"] + model = NonLinModelDAE(fq!, h!, Ts, nu, nx, na, ny, nd; p, as_0=[1e-5]) -vu, vd = ["\$q_B\$ (L/min)"], ["\$q_A\$ (L/min)"] -vx, vy = ["\$c_A\$ (mol/L)", "\$c_B\$ (mol/L)"], ["\$\\mathrm{pH}\$"] model = setname!(model, u=vu, x=vx, y=vy, d=vd) +``` -u = [10.0] -d = [10.0] -x_0 = [0.051, 0.049] -N = 61 -Y_data, U_data, D_data, X_data = zeros(ny, N), zeros(nu, N), zeros(nd, N), zeros(nx, N) -x = x_0 -let x=x, u=u, d=d - setstate!(model, x) +By default, an [`OrthogonalCollocation`](@ref) with 3 collocation points transcribes the +state dynamics and the algebraic equations into an optimization problem. A simple open-loop +simulation of `model` with: + +1. a bump on the base flow rate ``\mathbf{u} = q_{Bin}`` +2. a bump on the acid flow rare ``\mathbf{d} = q_{Ain}`` +3. a bump on the acid feed concentration ``c_{Ain}`` (as an unmeasured disturbance) + +validates that our DAE is well-posed: + +```@example 1 +function simDAE(model, N; x_0) + ny, ny, nd, nx = model.ny, model.ny, model.nd, model.nx + Y_data, U_data, D_data, X_data = zeros(ny, N), zeros(nu, N), zeros(nd, N), zeros(nx, N) + C_Ain_0 = model.p[1] + setstate!(model, x_0) + x = x_0 for i=1:N - d = i ≤ 2N÷3 ? [10.0] : [9.8] + u = i ≤ (1N÷4) ? [10.0] : [9.7] + d = i ≤ (2N÷4) ? [10.0] : [9.8] + c_Ain = i ≤ (3N÷4) ? C_Ain_0 : (C_Ain_0 + 0.05) + model.p[1] = c_Ain y = model(d) - u = i ≤ N÷3 ? [10.0] : [9.7] Y_data[:, i] = y U_data[:, i] = u D_data[:, i] = d X_data[:, i] = x x = updatestate!(model, u, d) end + model.p[1] = C_Ain_0 + return SimResult(model, U_data, Y_data, D_data; X_data) end -res = SimResult(model, U_data, Y_data, D_data; X_data) +x_0 = [0.0505, 0.0495] +N = 101 +res = simDAE(model, N; x_0) +``` +We plot the results by modifying the x-axis label to substitute the default time units +to hours: + +```@example 1 using Plots -#theme(:default) -theme(:dark) -default(fontfamily="Computer Modern"); scalefontsizes(1.1) -#p = plot(res, plotx=true, plotd=false, xlabel="Time (h)") -#xlabel!(p[3], "") -p = plot(res, plotd=true, xlabel="Time (h)") +plot(res, plotu=true, plotd=true, xlabel="Time (h)") +savefig("plot1_DAEpH.svg"); nothing # hide +``` + +![plot1_DAEpH](plot1_DAEpH.svg) + +## Adaptive Moving Horizon Estimation + +The default settings of the [`MovingHorizonEstimator`](@ref) assume that the measured +output is disturbed by a random-walk stochastic process (the pH). This is generally enough +to estimate the unmeasured disturbances in steady-state (the acid feed concentration). +To improve the interpretability of the results and the estimation performances, we can +instead disable the default stochastic model and construct an adaptive estimator. We first +need to augment the dynamics with our estimated parameter, the acid feed concentration +``c_{Ain}``: + +```@example 1 +calc_ċ_Ain( _ ) = 0 +function f̂q!(ẋ, res, x, a, u, d, p̂) + c_Bin, Kw, Ka, V = p̂ + q_Ain, q_Bin = d[1], u[1] + c_A, c_B = x[1], x[2] + c_Ain = x[3] + a_H = a[1] + q_out = q_Ain + q_Bin + c_Aout = c_A + c_Bout = c_B + ẋ[1] = calc_ċ_A(c_Ain, q_Ain, c_Aout, q_out) + ẋ[2] = calc_ċ_B(c_Bin, q_Bin, c_Bout, q_out) + ẋ[3] = calc_ċ_Ain(c_Ain) + res[1] = calc_res(a_H, c_A, c_B, Kw, Ka) + return nothing +end +ĥ!(y, x, a, d, p̂) = h!(y, x, a, d, p) +p̂ = [c_Bin, Kw, Ka, V] +nx̂ = 3 +vx̂ = [vx; raw"$c_{Ain}$ (mol/L)"] +model_aug = NonLinModelDAE(f̂q!, ĥ!, Ts, nu, nx̂, na, ny, nd; p, as_0=[1e-5]) +model_aug = setname!(model_aug, u=vu, x=vx̂, y=vy, d=vd) +``` + +Since `calc_ċ_Ain` returns `0`, the ``c_{Ain}`` parameter is assumed to be time-invariant. +More precisely, this feed property is assumed to be disturbed by a random-walk, instead of +the measured output directly. + +```@example 1 +mhe = MovingHorizonEstimator(model_aug, He=5, nint_ym=0) ``` From 3c0d90fb5caad9fc261ae55fdac7d8594f9d3d05 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Wed, 23 Sep 2026 18:02:22 -0400 Subject: [PATCH 05/13] doc: wip MHE on pH neutralization --- docs/src/manual/nonlinmpc2.md | 187 ++++++++++++++++++++++------------ 1 file changed, 121 insertions(+), 66 deletions(-) diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md index e5fb5884d..60619a78d 100644 --- a/docs/src/manual/nonlinmpc2.md +++ b/docs/src/manual/nonlinmpc2.md @@ -58,7 +58,7 @@ K_a &= \frac{[\mathrm{H}^+][\mathrm{Ac}^-]}{[\mathrm{H}\mathrm{Ac}]} \end{aligned} ``` -We respectively denote the algebraic variable and the two states with: +The following notation highlights three key concentrations in the model: ```math \begin{aligned} @@ -77,12 +77,6 @@ balance leads the algebraic equation: 0 = a_H + c_B - \frac{K_w}{a_H} - \frac{K_a c_A}{K_a + a_H} ``` -The pH is computed with: - -```math -\mathrm{pH} = -10 \log_{10}(a_H) -``` - !!! details "Reduction to an ODE" The algebraic equation can be further manipulated to produce this cubic expression: ```math @@ -93,6 +87,12 @@ The pH is computed with: When possible, plant model should be constructed with the specialized [`NonLinModel`](@ref) for ODEs. This tutorial will still treat the system as a DAE to illustrate its API. +The pH is: + +```math +\mathrm{pH} = -10 \log_{10}(a_H) ⟹ a_H = 10^{-\mathrm{pH}} +``` + ### Mass Balance Applying a mass balance on the weak acid ``A`` and the strong base ``B`` invariants leads @@ -118,14 +118,31 @@ because of the weir, the following relations compute the outflow terms: \end{aligned} ``` +### Model Construction + +The state and the algebraic vectors are respectively defined as: + +```math +\begin{aligned} + \mathbf{x} &= \begin{bmatrix} c_A \\ c_B \end{bmatrix} \\ + \mathbf{a} &= \mathrm{pH} +\end{aligned} +``` + +Alternatively, defining the algebraic vector as ``\mathbf{a} = a_H`` is a valid realization, +but a vastly inferior choice numerically since ``a_H`` spans around 14 orders of magnitude +(``10^{-1}`` to ``10^{-14}``) while the pH is bounded between roughly 1 to 14. Moreover, +it avoids a `log10` call that is undefined for negative values. + The [`NonLinModelDAE`](@ref) constructor expects that the state dynamics and the algebraic -equation is combined into a single `fq!(ẋ, res, x, a, u, d, p) -> nothing` function that +equation are combined into a single `fq!(ẋ, res, x, a, u, d, p) -> nothing` function that modifies both `ẋ` and `res` arguments in-place (an out-of-place option is also available), with the state dynamics and the residual of the algebraic equation, respectively: ```@example 1 using ModelPredictiveControl +calc_a_H(pH) = 10.0^(-pH) calc_ċ_A(c_Ain, q_Ain, c_Aout, q_out) = (60/V)*(q_Ain * c_Ain - q_out * c_Aout) calc_ċ_B(c_Bin, q_Bin, c_Bout, q_out) = (60/V)*(q_Bin * c_Bin - q_out * c_Bout) calc_res(a_H, c_A, c_B, Kw, Ka) = a_H + c_B - (Kw / a_H) - (Ka * c_A / (Ka + a_H)) @@ -133,10 +150,11 @@ function fq!(ẋ, res, x, a, u, d, p) c_Ain, c_Bin, Kw, Ka, V = p q_Ain, q_Bin = d[1], u[1] # [L/min], [L/min] c_A, c_B = x[1], x[2] # [mol/L], [mol/L] - a_H = a[1] # [mol/L] - q_out = q_Ain + q_Bin # [L/min] + pH = a[1] # [-] + q_out = q_Ain + q_Bin # [L/min] c_Aout = c_A # [mol/L] c_Bout = c_B # [mol/L] + a_H = calc_a_H(pH) ẋ[1] = calc_ċ_A(c_Ain, q_Ain, c_Aout, q_out) ẋ[2] = calc_ċ_B(c_Bin, q_Bin, c_Bout, q_out) res[1] = calc_res(a_H, c_A, c_B, Kw, Ka) @@ -147,31 +165,18 @@ end A similar in-place function is expected for the model output: ```@example 1 -calc_pH(a_H) = -log10(a_H) -function h!(y, _ , a, _ , _ ) - a_H = a[1] - pH = try - calc_pH(a_H) - catch myerror - myerror isa DomainError ? NaN : rethrow() - end - y[1] = pH - return nothing -end +h!(y, _ , a , _ , _ ) = (y .= a; nothing) ``` -By default, Julia throws a `DomainError` if `log10` is called with a negative number. The -`try` blocks is necessary to let the optimizer explores undefined domains by a returning a -`NaN` value in such cases. For similar reasons, providing an initial guess for the algebraic -variable `as_0` is crucial here to prioritize positive ``a_H`` concentration and ensure a -defined ``\mathrm{pH}`` value: +Providing an initial guess for the algebraic variable `as_0` is important to prioritize +positive pH solution: ```@example 1 -V = 1000.0 # reactor volume [L] c_Ain = 0.1 # feed concentration of weak acid [mol/L] c_Bin = 0.1 # feed concentration of strong base [mol/L] Kw = 1.0e-14 # water dissociation constant [mol^2/L^2] Ka = 1.75e-5 # acid dissociation constant [mol/L] +V = 1000.0 # reactor volume [L] Ts = 0.5 # Sample time [h] nu, nx, na, ny, nd = 1, 2, 1, 1, 1 @@ -180,45 +185,45 @@ p = [c_Ain, c_Bin, Kw, Ka, V] vu, vd = [raw"$q_{Bin}$ (L/min)"], [raw"$q_{Ain}$ (L/min)"] vx, vy = [raw"$c_A$ (mol/L)", raw"$c_B$ (mol/L)"], [raw"$\mathrm{pH}$"] -model = NonLinModelDAE(fq!, h!, Ts, nu, nx, na, ny, nd; p, as_0=[1e-5]) -model = setname!(model, u=vu, x=vx, y=vy, d=vd) +plant = NonLinModelDAE(fq!, h!, Ts, nu, nx, na, ny, nd; p=p, as_0=[7]) +plant = setname!(plant, u=vu, x=vx, y=vy, d=vd) ``` By default, an [`OrthogonalCollocation`](@ref) with 3 collocation points transcribes the state dynamics and the algebraic equations into an optimization problem. A simple open-loop -simulation of `model` with: +simulation of `plant` with: 1. a bump on the base flow rate ``\mathbf{u} = q_{Bin}`` 2. a bump on the acid flow rare ``\mathbf{d} = q_{Ain}`` -3. a bump on the acid feed concentration ``c_{Ain}`` (as an unmeasured disturbance) +3. a bump on the acid feed concentration ``c_{Ain}`` (an unmeasured disturbance) validates that our DAE is well-posed: ```@example 1 -function simDAE(model, N; x_0) - ny, ny, nd, nx = model.ny, model.ny, model.nd, model.nx +function simDAE(plant, N; x_0) + ny, ny, nd, nx = plant.ny, plant.ny, plant.nd, plant.nx Y_data, U_data, D_data, X_data = zeros(ny, N), zeros(nu, N), zeros(nd, N), zeros(nx, N) - C_Ain_0 = model.p[1] - setstate!(model, x_0) + c_Ain_0 = plant.p[1] + setstate!(plant, x_0) x = x_0 for i=1:N - u = i ≤ (1N÷4) ? [10.0] : [9.7] - d = i ≤ (2N÷4) ? [10.0] : [9.8] - c_Ain = i ≤ (3N÷4) ? C_Ain_0 : (C_Ain_0 + 0.05) - model.p[1] = c_Ain - y = model(d) + u = i ≤ (1N÷4) ? [10.0] : [11.0] + d = i ≤ (2N÷4) ? [10.0] : [12.0] + c_Ain = i ≤ (3N÷4) ? c_Ain_0 : (c_Ain_0 - 0.01) + plant.p[1] = c_Ain + y = plant(d) Y_data[:, i] = y U_data[:, i] = u D_data[:, i] = d X_data[:, i] = x - x = updatestate!(model, u, d) + x = updatestate!(plant, u, d) end - model.p[1] = C_Ain_0 - return SimResult(model, U_data, Y_data, D_data; X_data) + plant.p[1] = c_Ain_0 + return SimResult(plant, U_data, Y_data, D_data; X_data) end x_0 = [0.0505, 0.0495] N = 101 -res = simDAE(model, N; x_0) +res = simDAE(plant, N; x_0) ``` We plot the results by modifying the x-axis label to substitute the default time units @@ -236,41 +241,91 @@ savefig("plot1_DAEpH.svg"); nothing # hide The default settings of the [`MovingHorizonEstimator`](@ref) assume that the measured output is disturbed by a random-walk stochastic process (the pH). This is generally enough -to estimate the unmeasured disturbances in steady-state (the acid feed concentration). -To improve the interpretability of the results and the estimation performances, we can -instead disable the default stochastic model and construct an adaptive estimator. We first -need to augment the dynamics with our estimated parameter, the acid feed concentration -``c_{Ain}``: +to estimate the unmeasured disturbances in steady-state (the acid feed concentration, in +this case study). To improve the interpretability of the results and the estimation +performances, we can instead disable the default stochastic model and construct an adaptive +estimator. We first need to augment the dynamics with our estimated parameter, the acid feed +concentration ``c_{Ain}``: ```@example 1 calc_ċ_Ain( _ ) = 0 -function f̂q!(ẋ, res, x, a, u, d, p̂) +function f̂q!(dx̂, res, x̂, a, u, d, p̂) c_Bin, Kw, Ka, V = p̂ q_Ain, q_Bin = d[1], u[1] - c_A, c_B = x[1], x[2] - c_Ain = x[3] - a_H = a[1] - q_out = q_Ain + q_Bin + c_A, c_B = x̂[1], x̂[2] + c_Ain = x̂[3] + pH = a[1] + q_out = q_Ain + q_Bin c_Aout = c_A - c_Bout = c_B - ẋ[1] = calc_ċ_A(c_Ain, q_Ain, c_Aout, q_out) - ẋ[2] = calc_ċ_B(c_Bin, q_Bin, c_Bout, q_out) - ẋ[3] = calc_ċ_Ain(c_Ain) + c_Bout = c_B + a_H = calc_a_H(pH) + dx̂[1] = calc_ċ_A(c_Ain, q_Ain, c_Aout, q_out) + dx̂[2] = calc_ċ_B(c_Bin, q_Bin, c_Bout, q_out) + dx̂[3] = calc_ċ_Ain(c_Ain) res[1] = calc_res(a_H, c_A, c_B, Kw, Ka) return nothing end -ĥ!(y, x, a, d, p̂) = h!(y, x, a, d, p) +ĥ!(y, x̂, a, d, p̂) = h!(y, x̂, a, d, p̂) p̂ = [c_Bin, Kw, Ka, V] -nx̂ = 3 +nx̂ = nx + 1 vx̂ = [vx; raw"$c_{Ain}$ (mol/L)"] -model_aug = NonLinModelDAE(f̂q!, ĥ!, Ts, nu, nx̂, na, ny, nd; p, as_0=[1e-5]) -model_aug = setname!(model_aug, u=vu, x=vx̂, y=vy, d=vd) +model = NonLinModelDAE(f̂q!, ĥ!, Ts, nu, nx̂, na, ny, nd; p=p̂, as_0=[7]) +model = setname!(model, u=vu, x=vx̂, y=vy, d=vd) +``` + +Since `calc_ċ_Ain` always returns `0`, the ``c_{Ain}`` parameter is assumed to be +time-invariant. More precisely, this concentration of the acid feed is assumed to be +disturbed by a random-walk, instead of the measured output. Among all the settings of the +[`MovingHorizonEstimator`](@ref), a proper tuning of the covariance matrices through `σQ` +`σR` and `σP_0`, a past horizon `He` long enough to see the whole dynamics, and an exact +Hessian matrix can improve the stability on a highly nonlinear and stiff plant model like +here. We can also bound the three estimated states to positive values since they are +concentration in mol/L: + +```@example 1 +nint_ym=0; nint_u=0; # disable the default stochastic model +He = 10; hessian=true +σQ = [0.01, 0.01, 0.1]; σR=[0.5]; σP_0 = 0.1*σQ +mhe = MovingHorizonEstimator(model; nint_ym, nint_u, He, hessian, σQ, σR, σP_0) +using JuMP; unset_time_limit_sec(mhe.optim) # no wall time limit at optimization +mhe = setconstraint!(mhe, x̂min=[0, 0, 0]) ``` -Since `calc_ċ_Ain` returns `0`, the ``c_{Ain}`` parameter is assumed to be time-invariant. -More precisely, this feed property is assumed to be disturbed by a random-walk, instead of -the measured output directly. +The state constraints are shown in round bracket next to the decision variables. There are +33 of them (3 states × 10 datapoints in the pasts + 3 others for the arrival estimate). The +arrival covariance ``\mathbf{P̄}`` is constant by default for [`NonLinModelDAE`](@ref). A +proper tuning of `σP_0` and `He` reduces the impact of this approximation. We can now +reproduce the last simulated scenario and see how `mhe` performs: ```@example 1 -mhe = MovingHorizonEstimator(model_aug, He=5, nint_ym=0) +function simMHE(mhe, plant, N; x_0, x̂_0) + ny, ny, nd, nx, nx̂ = plant.ny, plant.ny, plant.nd, plant.nx, mhe.nx̂ + Y_data, U_data, D_data, X_data = zeros(ny, N), zeros(nu, N), zeros(nd, N), zeros(nx, N) + Ŷ_data, X̂_data = zeros(ny, N), zeros(nx̂, N) + c_Ain_0 = plant.p[1] + initstate!(mhe, [7], [10], [10]) + setstate!(plant, x_0); setstate!(mhe, x̂_0) + x = x_0 + for i=1:N + u = i ≤ (1N÷4) ? [10.0] : [11.0] + d = i ≤ (2N÷4) ? [10.0] : [12.0] + c_Ain = i ≤ (3N÷4) ? c_Ain_0 : (c_Ain_0 - 0.01) + plant.p[1] = c_Ain + y = plant(d) + x̂ = preparestate!(mhe, y, d) + ŷ = mhe(d) + Y_data[:, i] = y + U_data[:, i] = u + D_data[:, i] = d + X_data[:, i] = x + Ŷ_data[:, i] = ŷ + X̂_data[:, i] = x̂ + x = updatestate!(plant, u, d) + end + plant.p[1] = c_Ain_0 + return SimResult(mhe, U_data, Y_data, D_data; plant, X_data, X̂_data, Ŷ_data) +end +x̂_0=[x_0; c_Ain] +res = simMHE(mhe, plant, 100; x_0, x̂_0) +plot(res, plotd=false, plotxwithx̂=true) ``` From bbca02e578bbef1b4b682c060c99365de11b753b Mon Sep 17 00:00:00 2001 From: franckgaga Date: Wed, 23 Sep 2026 18:10:14 -0400 Subject: [PATCH 06/13] doc: change x-label units --- docs/src/manual/nonlinmpc2.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md index 60619a78d..f14bb78f0 100644 --- a/docs/src/manual/nonlinmpc2.md +++ b/docs/src/manual/nonlinmpc2.md @@ -327,5 +327,5 @@ function simMHE(mhe, plant, N; x_0, x̂_0) end x̂_0=[x_0; c_Ain] res = simMHE(mhe, plant, 100; x_0, x̂_0) -plot(res, plotd=false, plotxwithx̂=true) +p = plot(res, plotd=false, plotxwithx̂=true, xlabel="Time (h)") ``` From c43e2d39130e8d59e07c677fdaa838747a164732 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Wed, 23 Sep 2026 18:40:18 -0400 Subject: [PATCH 07/13] doc: wip mhe tuning --- docs/src/manual/nonlinmpc2.md | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md index f14bb78f0..70e6eb87c 100644 --- a/docs/src/manual/nonlinmpc2.md +++ b/docs/src/manual/nonlinmpc2.md @@ -284,8 +284,8 @@ concentration in mol/L: ```@example 1 nint_ym=0; nint_u=0; # disable the default stochastic model -He = 10; hessian=true -σQ = [0.01, 0.01, 0.1]; σR=[0.5]; σP_0 = 0.1*σQ +He = 10; hessian = true +σQ = [0.01, 0.01, 0.005]; σR=[0.1]; σP_0 = 100*[0.001, 0.001, 1e-6] mhe = MovingHorizonEstimator(model; nint_ym, nint_u, He, hessian, σQ, σR, σP_0) using JuMP; unset_time_limit_sec(mhe.optim) # no wall time limit at optimization mhe = setconstraint!(mhe, x̂min=[0, 0, 0]) @@ -295,9 +295,10 @@ The state constraints are shown in round bracket next to the decision variables. 33 of them (3 states × 10 datapoints in the pasts + 3 others for the arrival estimate). The arrival covariance ``\mathbf{P̄}`` is constant by default for [`NonLinModelDAE`](@ref). A proper tuning of `σP_0` and `He` reduces the impact of this approximation. We can now -reproduce the last simulated scenario and see how `mhe` performs: +reproduce the last simulated scenario and see how `mhe` performs under pH measurement noise: ```@example 1 +using Random function simMHE(mhe, plant, N; x_0, x̂_0) ny, ny, nd, nx, nx̂ = plant.ny, plant.ny, plant.nd, plant.nx, mhe.nx̂ Y_data, U_data, D_data, X_data = zeros(ny, N), zeros(nu, N), zeros(nd, N), zeros(nx, N) @@ -311,7 +312,7 @@ function simMHE(mhe, plant, N; x_0, x̂_0) d = i ≤ (2N÷4) ? [10.0] : [12.0] c_Ain = i ≤ (3N÷4) ? c_Ain_0 : (c_Ain_0 - 0.01) plant.p[1] = c_Ain - y = plant(d) + y = plant(d) + 0.1*randn(1) x̂ = preparestate!(mhe, y, d) ŷ = mhe(d) Y_data[:, i] = y @@ -321,6 +322,7 @@ function simMHE(mhe, plant, N; x_0, x̂_0) Ŷ_data[:, i] = ŷ X̂_data[:, i] = x̂ x = updatestate!(plant, u, d) + x̂ = updatestate!(mhe, y, u, d) end plant.p[1] = c_Ain_0 return SimResult(mhe, U_data, Y_data, D_data; plant, X_data, X̂_data, Ŷ_data) From 0b81f507b878415d2518dcfa9c8c873ea00d74f7 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Thu, 24 Sep 2026 15:34:16 -0400 Subject: [PATCH 08/13] =?UTF-8?q?doc:=20statisfied=20by=20the=20MHE=20tuni?= =?UTF-8?q?ng=20(not=20easy=20=F0=9F=98=AE=E2=80=8D=F0=9F=92=A8)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- docs/src/manual/nonlinmpc2.md | 72 ++++++++++++++++++++--------------- 1 file changed, 41 insertions(+), 31 deletions(-) diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md index 70e6eb87c..9dda6f16f 100644 --- a/docs/src/manual/nonlinmpc2.md +++ b/docs/src/manual/nonlinmpc2.md @@ -168,8 +168,8 @@ A similar in-place function is expected for the model output: h!(y, _ , a , _ , _ ) = (y .= a; nothing) ``` -Providing an initial guess for the algebraic variable `as_0` is important to prioritize -positive pH solution: +Providing an initial guess for the state `xs_0` and algebraic variable `as_0` is important +for DAEs, to prioritize positive pH and concentration solution, *inter alia*: ```@example 1 c_Ain = 0.1 # feed concentration of weak acid [mol/L] @@ -185,13 +185,17 @@ p = [c_Ain, c_Bin, Kw, Ka, V] vu, vd = [raw"$q_{Bin}$ (L/min)"], [raw"$q_{Ain}$ (L/min)"] vx, vy = [raw"$c_A$ (mol/L)", raw"$c_B$ (mol/L)"], [raw"$\mathrm{pH}$"] -plant = NonLinModelDAE(fq!, h!, Ts, nu, nx, na, ny, nd; p=p, as_0=[7]) +transcription = TrapezoidalCollocation() +xs_0, as_0 = [0.025, 0.025], [7] + +plant = NonLinModelDAE(fq!, h!, Ts, nu, nx, na, ny, nd; p=p, xs_0, as_0, transcription) plant = setname!(plant, u=vu, x=vx, y=vy, d=vd) ``` -By default, an [`OrthogonalCollocation`](@ref) with 3 collocation points transcribes the -state dynamics and the algebraic equations into an optimization problem. A simple open-loop -simulation of `plant` with: +We use a [`TrapezoidalCollocation`](@ref) transcription instead of the default +[`OrthogonalCollocation`](@ref), since it is less computationnaly expensive and its accuracy +and stability is good enough for this case study. A simple open-loop simulation of `plant` +with: 1. a bump on the base flow rate ``\mathbf{u} = q_{Bin}`` 2. a bump on the acid flow rare ``\mathbf{d} = q_{Ain}`` @@ -222,7 +226,7 @@ function simDAE(plant, N; x_0) return SimResult(plant, U_data, Y_data, D_data; X_data) end x_0 = [0.0505, 0.0495] -N = 101 +N = 81 res = simDAE(plant, N; x_0) ``` @@ -269,7 +273,8 @@ ĥ!(y, x̂, a, d, p̂) = h!(y, x̂, a, d, p̂) p̂ = [c_Bin, Kw, Ka, V] nx̂ = nx + 1 vx̂ = [vx; raw"$c_{Ain}$ (mol/L)"] -model = NonLinModelDAE(f̂q!, ĥ!, Ts, nu, nx̂, na, ny, nd; p=p̂, as_0=[7]) +x̂s_0 = [xs_0; 0.1] +model = NonLinModelDAE(f̂q!, ĥ!, Ts, nu, nx̂, na, ny, nd; p=p̂, xs_0=x̂s_0, as_0, transcription) model = setname!(model, u=vu, x=vx̂, y=vy, d=vd) ``` @@ -284,50 +289,55 @@ concentration in mol/L: ```@example 1 nint_ym=0; nint_u=0; # disable the default stochastic model -He = 10; hessian = true -σQ = [0.01, 0.01, 0.005]; σR=[0.1]; σP_0 = 100*[0.001, 0.001, 1e-6] +He = 8; hessian = true +σQ = [0.0015, 0.0015, 2e-4]; σR=[0.05]; σP_0 = [0.05, 0.05, 1e-3] #5e-4] mhe = MovingHorizonEstimator(model; nint_ym, nint_u, He, hessian, σQ, σR, σP_0) -using JuMP; unset_time_limit_sec(mhe.optim) # no wall time limit at optimization mhe = setconstraint!(mhe, x̂min=[0, 0, 0]) +using JuMP; unset_time_limit_sec(mhe.optim) # no wall time limit during optimization ``` -The state constraints are shown in round bracket next to the decision variables. There are -33 of them (3 states × 10 datapoints in the pasts + 3 others for the arrival estimate). The -arrival covariance ``\mathbf{P̄}`` is constant by default for [`NonLinModelDAE`](@ref). A -proper tuning of `σP_0` and `He` reduces the impact of this approximation. We can now -reproduce the last simulated scenario and see how `mhe` performs under pH measurement noise: +The state constraints are shown in round brackets next to the decision variables. There are +27 of them (3 states × 8 datapoints in the pasts + 3 arrival estimates). The arrival +covariance ``\mathbf{P̄}`` is constant by default for [`NonLinModelDAE`](@ref). A proper +tuning of `σP_0` and `He` reduces the impact of this approximation. We can now reproduce the +last simulated scenario and see how `mhe` performs under pH measurement noise: ```@example 1 using Random function simMHE(mhe, plant, N; x_0, x̂_0) ny, ny, nd, nx, nx̂ = plant.ny, plant.ny, plant.nd, plant.nx, mhe.nx̂ - Y_data, U_data, D_data, X_data = zeros(ny, N), zeros(nu, N), zeros(nd, N), zeros(nx, N) + Y_data, U_data, D_data = zeros(ny, N), zeros(nu, N), zeros(nd, N) + X_data = zeros(nx+1, N) # to store also the actual c_Ain value Ŷ_data, X̂_data = zeros(ny, N), zeros(nx̂, N) c_Ain_0 = plant.p[1] - initstate!(mhe, [7], [10], [10]) setstate!(plant, x_0); setstate!(mhe, x̂_0) + initstate!(mhe, [7], [10], [10]) x = x_0 for i=1:N u = i ≤ (1N÷4) ? [10.0] : [11.0] d = i ≤ (2N÷4) ? [10.0] : [12.0] c_Ain = i ≤ (3N÷4) ? c_Ain_0 : (c_Ain_0 - 0.01) plant.p[1] = c_Ain - y = plant(d) + 0.1*randn(1) - x̂ = preparestate!(mhe, y, d) - ŷ = mhe(d) - Y_data[:, i] = y - U_data[:, i] = u - D_data[:, i] = d - X_data[:, i] = x - Ŷ_data[:, i] = ŷ - X̂_data[:, i] = x̂ + y = evaloutput(plant, d) + ym = y + 0.05*randn(1) + dm = d + 0.10*randn(1) + x̂ = preparestate!(mhe, ym, dm) + ŷ = evaloutput(mhe, dm) + Y_data[:, i] = ym + U_data[:, i] = u + D_data[:, i] = dm + X_data[1:2, i] = x + X_data[3, i] = c_Ain + Ŷ_data[:, i] = ŷ + X̂_data[:, i] = x̂ x = updatestate!(plant, u, d) - x̂ = updatestate!(mhe, y, u, d) + x̂ = updatestate!(mhe, ym, u, dm) end plant.p[1] = c_Ain_0 return SimResult(mhe, U_data, Y_data, D_data; plant, X_data, X̂_data, Ŷ_data) end -x̂_0=[x_0; c_Ain] -res = simMHE(mhe, plant, 100; x_0, x̂_0) -p = plot(res, plotd=false, plotxwithx̂=true, xlabel="Time (h)") +x̂_0 = [0.025, 0.025, c_Ain] +res = simMHE(mhe, plant, N; x_0, x̂_0)# N; x_0, x̂_0) +# T = @elapsed +p = plot(res, plotd=false, plotu=false, plotxwithx̂=true, plotx̂min=false, xlabel="Time (h)") ``` From 71ac51c7d4ff824a0a6b9c8d4d199eae1fb463f7 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Thu, 24 Sep 2026 16:10:23 -0400 Subject: [PATCH 09/13] doc: almost finished MHE on pH neutralization process --- docs/src/manual/nonlinmpc2.md | 28 ++++++++++++++++++++++------ 1 file changed, 22 insertions(+), 6 deletions(-) diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md index 9dda6f16f..fb970fa96 100644 --- a/docs/src/manual/nonlinmpc2.md +++ b/docs/src/manual/nonlinmpc2.md @@ -290,10 +290,10 @@ concentration in mol/L: ```@example 1 nint_ym=0; nint_u=0; # disable the default stochastic model He = 8; hessian = true -σQ = [0.0015, 0.0015, 2e-4]; σR=[0.05]; σP_0 = [0.05, 0.05, 1e-3] #5e-4] +σQ = [0.0015, 0.0015, 2e-4]; σR=[0.05]; σP_0 = [0.05, 0.05, 5e-4] mhe = MovingHorizonEstimator(model; nint_ym, nint_u, He, hessian, σQ, σR, σP_0) -mhe = setconstraint!(mhe, x̂min=[0, 0, 0]) using JuMP; unset_time_limit_sec(mhe.optim) # no wall time limit during optimization +mhe = setconstraint!(mhe, x̂min=[0, 0, 0]) ``` The state constraints are shown in round brackets next to the decision variables. There are @@ -307,11 +307,10 @@ using Random function simMHE(mhe, plant, N; x_0, x̂_0) ny, ny, nd, nx, nx̂ = plant.ny, plant.ny, plant.nd, plant.nx, mhe.nx̂ Y_data, U_data, D_data = zeros(ny, N), zeros(nu, N), zeros(nd, N) - X_data = zeros(nx+1, N) # to store also the actual c_Ain value + X_data = zeros(nx+1, N) # nx+1 to store the actual c_Ain value in the last row Ŷ_data, X̂_data = zeros(ny, N), zeros(nx̂, N) c_Ain_0 = plant.p[1] setstate!(plant, x_0); setstate!(mhe, x̂_0) - initstate!(mhe, [7], [10], [10]) x = x_0 for i=1:N u = i ≤ (1N÷4) ? [10.0] : [11.0] @@ -337,7 +336,24 @@ function simMHE(mhe, plant, N; x_0, x̂_0) return SimResult(mhe, U_data, Y_data, D_data; plant, X_data, X̂_data, Ŷ_data) end x̂_0 = [0.025, 0.025, c_Ain] -res = simMHE(mhe, plant, N; x_0, x̂_0)# N; x_0, x̂_0) -# T = @elapsed +res = simMHE(mhe, plant, N; x_0, x̂_0) p = plot(res, plotd=false, plotu=false, plotxwithx̂=true, plotx̂min=false, xlabel="Time (h)") +xlabel!(p[2], ""); xlabel!(p[3], "") # remove xlabel on c_A and c_B plots +savefig(p, "plot2_DAEpH.svg"); nothing # hide ``` + +![plot2_DAEpH](plot2_DAEpH.svg) + +The estimated acid concentration does not perfectly converge towards the actual value, +but it is a well-known issue of adaptive estimation and control. A persistent excitation on +``\mathbf{u}`` like an additive dither signal would presumably improve the estimation +performances. With a sampling time of 30 min, the solving of the optimization problem is +obviously fast enough for realtime execution and application to closed-loop control: + +```@example +T = @elapsed simMHE(mhe, plant, N; x_0, x̂_0) +println("Total optimization and simulation time for $N time steps: $T s") +``` + +Perhaps more importantly, the fast simulations ease the tuning of the estimation horizon and +covariance matrices, for iterative and trial-and-error approaches. From cbc543657e12a981d49329d998b64b0c3600a177 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Thu, 24 Sep 2026 16:19:15 -0400 Subject: [PATCH 10/13] doc: debug markdown block --- docs/src/manual/nonlinmpc2.md | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md index fb970fa96..f5c04b90d 100644 --- a/docs/src/manual/nonlinmpc2.md +++ b/docs/src/manual/nonlinmpc2.md @@ -298,9 +298,10 @@ mhe = setconstraint!(mhe, x̂min=[0, 0, 0]) The state constraints are shown in round brackets next to the decision variables. There are 27 of them (3 states × 8 datapoints in the pasts + 3 arrival estimates). The arrival -covariance ``\mathbf{P̄}`` is constant by default for [`NonLinModelDAE`](@ref). A proper -tuning of `σP_0` and `He` reduces the impact of this approximation. We can now reproduce the -last simulated scenario and see how `mhe` performs under pH measurement noise: +covariance ``\mathbf{P̄}`` is constant by default for [`NonLinModelDAE`](@ref), specified by +`σP_0` argument. A proper tuning of `σP_0` and `He` reduces the impact of the constant +arrival approximation. We can now reproduce the last simulated scenario and see how `mhe` +performs under pH and flow rate measurement noise: ```@example 1 using Random @@ -350,7 +351,7 @@ but it is a well-known issue of adaptive estimation and control. A persistent ex performances. With a sampling time of 30 min, the solving of the optimization problem is obviously fast enough for realtime execution and application to closed-loop control: -```@example +```@example 1 T = @elapsed simMHE(mhe, plant, N; x_0, x̂_0) println("Total optimization and simulation time for $N time steps: $T s") ``` From fa9c1eff554235754a6c312b581fcd6da2d7aee7 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Thu, 24 Sep 2026 16:34:34 -0400 Subject: [PATCH 11/13] doc: minor correction in pH neutralization --- docs/src/manual/nonlinmpc2.md | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md index f5c04b90d..8c3146faf 100644 --- a/docs/src/manual/nonlinmpc2.md +++ b/docs/src/manual/nonlinmpc2.md @@ -281,11 +281,12 @@ model = setname!(model, u=vu, x=vx̂, y=vy, d=vd) Since `calc_ċ_Ain` always returns `0`, the ``c_{Ain}`` parameter is assumed to be time-invariant. More precisely, this concentration of the acid feed is assumed to be disturbed by a random-walk, instead of the measured output. Among all the settings of the -[`MovingHorizonEstimator`](@ref), a proper tuning of the covariance matrices through `σQ` -`σR` and `σP_0`, a past horizon `He` long enough to see the whole dynamics, and an exact -Hessian matrix can improve the stability on a highly nonlinear and stiff plant model like -here. We can also bound the three estimated states to positive values since they are -concentration in mol/L: +[`MovingHorizonEstimator`](@ref), a proper tuning of the covariance matrices through `σQ`, +`σR` and `σP_0`, and a past horizon `He` long enough to see the main dynamics can improve +the stability on a highly nonlinear and stiff plant model like here. An exact Hessian matrix +with `hessian=true` also helps for DAEs, since the dynamics are encoded in the nonlinear +equality constraints. We can also bound the three estimated states to positive values since +they are concentration in mol/L: ```@example 1 nint_ym=0; nint_u=0; # disable the default stochastic model @@ -345,11 +346,12 @@ savefig(p, "plot2_DAEpH.svg"); nothing # hide ![plot2_DAEpH](plot2_DAEpH.svg) -The estimated acid concentration does not perfectly converge towards the actual value, -but it is a well-known issue of adaptive estimation and control. A persistent excitation on -``\mathbf{u}`` like an additive dither signal would presumably improve the estimation -performances. With a sampling time of 30 min, the solving of the optimization problem is -obviously fast enough for realtime execution and application to closed-loop control: +The estimated acid feed concentration ``c_{Ain}`` does not perfectly converge towards the +actual value, but it is a well-known issue of adaptive estimation and control. A persistent +excitation on ``\mathbf{u}`` like an additive dither signal would presumably improve the +estimation performances. With a sampling time of 30 min, the solving of the optimization +problem is obviously fast enough for realtime execution and application to closed-loop +control: ```@example 1 T = @elapsed simMHE(mhe, plant, N; x_0, x̂_0) From 72035ac2a93a1089f116e7ce77980afdc9ff0043 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Thu, 24 Sep 2026 16:46:37 -0400 Subject: [PATCH 12/13] doc: remove outputs when not necessary --- docs/src/manual/nonlinmpc2.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md index 8c3146faf..95854dafa 100644 --- a/docs/src/manual/nonlinmpc2.md +++ b/docs/src/manual/nonlinmpc2.md @@ -159,13 +159,13 @@ function fq!(ẋ, res, x, a, u, d, p) ẋ[2] = calc_ċ_B(c_Bin, q_Bin, c_Bout, q_out) res[1] = calc_res(a_H, c_A, c_B, Kw, Ka) return nothing -end +end; ``` A similar in-place function is expected for the model output: ```@example 1 -h!(y, _ , a , _ , _ ) = (y .= a; nothing) +h!(y, _ , a , _ , _ ) = (y .= a; nothing); ``` Providing an initial guess for the state `xs_0` and algebraic variable `as_0` is important @@ -313,6 +313,7 @@ function simMHE(mhe, plant, N; x_0, x̂_0) Ŷ_data, X̂_data = zeros(ny, N), zeros(nx̂, N) c_Ain_0 = plant.p[1] setstate!(plant, x_0); setstate!(mhe, x̂_0) + initstate!(mhe, [10], [7], [10]) x = x_0 for i=1:N u = i ≤ (1N÷4) ? [10.0] : [11.0] From f8e54a43dc48c9c96dfb9d110b20b3924edd457d Mon Sep 17 00:00:00 2001 From: franckgaga Date: Thu, 24 Sep 2026 16:57:41 -0400 Subject: [PATCH 13/13] doc: idem --- docs/src/manual/nonlinmpc2.md | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md index 95854dafa..e87912306 100644 --- a/docs/src/manual/nonlinmpc2.md +++ b/docs/src/manual/nonlinmpc2.md @@ -159,13 +159,15 @@ function fq!(ẋ, res, x, a, u, d, p) ẋ[2] = calc_ċ_B(c_Bin, q_Bin, c_Bout, q_out) res[1] = calc_res(a_H, c_A, c_B, Kw, Ka) return nothing -end; +end +nothing # hide ``` A similar in-place function is expected for the model output: ```@example 1 -h!(y, _ , a , _ , _ ) = (y .= a; nothing); +h!(y, _ , a , _ , _ ) = (y .= a; nothing) +nothing # hide ``` Providing an initial guess for the state `xs_0` and algebraic variable `as_0` is important