JuliaControl · franckgaga · Sep 23, 2026 · Sep 20, 2026 · Sep 20, 2026 · Sep 20, 2026
diff --git a/Project.toml b/Project.toml
@@ -1,6 +1,6 @@
 name = "ModelPredictiveControl"
 uuid = "61f9bdb8-6ae4-484a-811f-bbf86720c31c"
-version = "2.13.2"
+version = "2.14.0"
 authors = ["Francis Gagnon"]
 
 [deps]

diff --git a/README.md b/README.md
@@ -82,7 +82,7 @@ for more detailed examples.
 
 ### 🎯 Model Predictive Control Features
 
-- 🏭️ **Plant Model**: Linear or nonlinear models exploiting multiple dispatch.
+- 🏭️ **Plant Model**: Linear, nonlinear ODE and nonlinear DAE models exploiting multiple dispatch.
 - ⛳️ **Objectives**: Tracking for inputs/outputs, move suppression, terminal costs, and economic costs.
 - ⏳️ **Horizons**: Distinct prediction/control horizons with custom move blocking.
 - 📸 **Linearization**: Auto-differentiation for exact Jacobians.

diff --git a/docs/src/internals/sim_model.md b/docs/src/internals/sim_model.md
@@ -22,6 +22,8 @@ ModelPredictiveControl.init_defectmat_dae
 ```@docs
 ModelPredictiveControl.f!
 ModelPredictiveControl.h!
+ModelPredictiveControl.fq_dae!
+ModelPredictiveControl.h_dae!
 ```
 
 ## Init State

diff --git a/docs/src/internals/state_estim.md b/docs/src/internals/state_estim.md
@@ -67,8 +67,6 @@ ModelPredictiveControl.get_nonlincon_oracle(::MovingHorizonEstimator, ::ModelPre
 ```@docs
 ModelPredictiveControl.f̂!
 ModelPredictiveControl.ĥ!
-ModelPredictiveControl.fq_dae!
-ModelPredictiveControl.ĥ_dae!
 ```
 
 ## Remove Operating Points

diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md
@@ -4,6 +4,10 @@
 Pages = ["nonlinmpc2.md"]
 ```
 
+!!! todo "Under Construction"
+    This tutorial is currently under construction. Only the modeling part is written for
+    now.
+
 ## Nonlinear Model (DAE)
 
 In this example, the goal is to control the pH of a solution in a continuously stirred

diff --git a/src/controller/linmpc.jl b/src/controller/linmpc.jl
@@ -1,5 +1,3 @@
-const DEFAULT_LINMPC_TRANSCRIPTION = SingleShooting()
-
 struct LinMPC{
     NT<:Real, 
     SE<:StateEstimator, 
@@ -241,7 +239,7 @@ function LinMPC(
     Wd = nothing,
     Wr = nothing,
     Cwt = DEFAULT_CWT,
-    transcription::ShootingMethod = DEFAULT_LINMPC_TRANSCRIPTION,
+    transcription::ShootingMethod = SingleShooting(),
     optim::JuMP.GenericModel = JuMP.Model(DEFAULT_QP_OPTIMIZER, add_bridges=true),
     kwargs...
 )
@@ -300,7 +298,7 @@ function LinMPC(
     Wd = nothing,
     Wr = nothing,
     Cwt  = DEFAULT_CWT,
-    transcription::ShootingMethod = DEFAULT_LINMPC_TRANSCRIPTION,
+    transcription::ShootingMethod = SingleShooting(),
     optim::JM = JuMP.Model(DEFAULT_QP_OPTIMIZER, add_bridges=true)
 ) where {NT<:Real, SE<:StateEstimator{NT}, JM<:JuMP.GenericModel}
     isa(estim.model, LinModel) || error(MSG_LINMODEL_ERR) 

diff --git a/src/controller/nonlinmpc.jl b/src/controller/nonlinmpc.jl
@@ -1,4 +1,3 @@
-const DEFAULT_NONLINMPC_TRANSCRIPTION = SingleShooting()
 const DEFAULT_NONLINMPC_HESSIAN = AutoSparse(
     AutoForwardDiff();
     sparsity_detector=TracerSparsityDetector(),
@@ -365,7 +364,7 @@ function NonLinMPC(
     gc ::Function = gc!,
     nc::Int = 0,
     p = model.p,
-    transcription::TranscriptionMethod = DEFAULT_NONLINMPC_TRANSCRIPTION,
+    transcription::TranscriptionMethod = SingleShooting(),
     optim::JuMP.GenericModel = JuMP.Model(DEFAULT_NLP_OPTIMIZER, add_bridges=false),
     gradient::AbstractADType = DEFAULT_GRADIENT,
     jacobian::AbstractADType = default_jacobian(transcription),
@@ -381,9 +380,6 @@ function NonLinMPC(
     )
 end
 
-default_estimator(model::SimModelODE; kwargs...) = UnscentedKalmanFilter(model; kwargs...)
-default_estimator(model::LinModel; kwargs...) = SteadyKalmanFilter(model; kwargs...)
-
 """
     NonLinMPC(estim::StateEstimator; <keyword arguments>)
 
@@ -441,7 +437,7 @@ function NonLinMPC(
     gc ::Function = gc!,
     nc = 0,
     p = estim.model.p,
-    transcription::TranscriptionMethod = DEFAULT_NONLINMPC_TRANSCRIPTION,
+    transcription::TranscriptionMethod = SingleShooting(),
     optim::JuMP.GenericModel = JuMP.Model(DEFAULT_NLP_OPTIMIZER, add_bridges=false),
     gradient::AbstractADType = DEFAULT_GRADIENT,
     jacobian::AbstractADType = default_jacobian(transcription),
@@ -467,6 +463,9 @@ function NonLinMPC(
     )
 end
 
+default_estimator(model::SimModelODE; kwargs...) = UnscentedKalmanFilter(model; kwargs...)
+default_estimator(model::LinModel; kwargs...) = SteadyKalmanFilter(model; kwargs...)
+
 """
     validate_JE(NT, JE) -> nothing
 

diff --git a/src/estimator/construct.jl b/src/estimator/construct.jl
@@ -3,13 +3,14 @@ struct StateEstimatorBuffer{NT<:Real}
     û ::Vector{NT}
     k̄ ::Vector{NT}
     x̂ ::Vector{NT}
-    a ::Vector{NT}
+    â ::Vector{NT}
     Z̃ ::Vector{NT}
     V̂ ::Vector{NT}
     Ŵ ::Vector{NT}
     X̂ ::Vector{NT}
-    A ::Vector{NT}
+    Â ::Vector{NT}
     Ŷ ::Vector{NT}
+    U ::Vector{NT}
     D ::Vector{NT}
     P̂ ::Matrix{NT}
     Q̂ ::Matrix{NT}
@@ -38,18 +39,19 @@ function StateEstimatorBuffer{NT}(
     transcription::TranscriptionMethod = SingleShooting()
 ) where NT <: Real
     nZ̃ = nε + get_nZ_mhe(transcription, He, nx̂, nk̄, nŵ, na)
-    nV̂, nŴ, nX̂, nA, nŶ, nD = nym*He, nŵ*He, nx̂*He, na*He, ny*He, nd*(He+1)
+    nV̂, nŴ, nX̂, nÂ, nŶ, nU, nD = nym*He, nŵ*He, nx̂*He, na*He, ny*He, nu*He, nd*(He+1)
     u  = Vector{NT}(undef, nu)
     û  = Vector{NT}(undef, nu)
     k̄  = Vector{NT}(undef, nk̄)
     x̂  = Vector{NT}(undef, nx̂)
-    a  = Vector{NT}(undef, na)
+    â  = Vector{NT}(undef, na)
     Z̃  = Vector{NT}(undef, nZ̃)
     V̂  = Vector{NT}(undef, nV̂)
     Ŵ  = Vector{NT}(undef, nŴ)
     X̂  = Vector{NT}(undef, nX̂)
-    A  = Vector{NT}(undef, nA)
+    Â  = Vector{NT}(undef, nÂ)
     Ŷ  = Vector{NT}(undef, nŶ)
+    U  = Vector{NT}(undef, nU)
     D  = Vector{NT}(undef, nD)
     P̂  = Matrix{NT}(undef, nx̂, nx̂)
     Q̂  = Matrix{NT}(undef, nx̂, nx̂)
@@ -60,7 +62,7 @@ function StateEstimatorBuffer{NT}(
     d  = Vector{NT}(undef, nd)
     empty = Vector{NT}(undef, 0)
     return StateEstimatorBuffer{NT}(
-        u, û, k̄, x̂, a, Z̃, V̂, Ŵ, X̂, A, Ŷ, D, P̂, Q̂, R̂, K̂, ym, ŷ, d, empty
+        u, û, k̄, x̂, â, Z̃, V̂, Ŵ, X̂, Â, Ŷ, U, D, P̂, Q̂, R̂, K̂, ym, ŷ, d, empty
     )
 end
 

diff --git a/src/estimator/execute.jl b/src/estimator/execute.jl
@@ -98,12 +98,12 @@ Same than [`f̂!`](@ref) for [`SimModel`](@ref) but without the `estim` argument
 """
 function f̂!(x̂0next, û0, k̄, model::SimModel, As, Cs_u, f̂op, x̂op, x̂0, u0, d0)
     # `@views` macro avoid copies with matrix slice operator e.g. [a:b]
-    @views xd, xs = x̂0[1:model.nx], x̂0[model.nx+1:end]
-    @views xdnext, xsnext = x̂0next[1:model.nx], x̂0next[model.nx+1:end]
-    mul!(û0, Cs_u, xs)      # ys_u = Cs_u*xs
+    @views x̂d, x̂s = x̂0[1:model.nx], x̂0[model.nx+1:end]
+    @views x̂dnext, x̂snext = x̂0next[1:model.nx], x̂0next[model.nx+1:end]
+    mul!(û0, Cs_u, x̂s)      # ys_u = Cs_u*xs
     û0 .+= u0               # û0 = u0 + ys_u  
-    f!(xdnext, k̄, model, xd, û0, d0, model.p)
-    mul!(xsnext, As, xs)
+    f!(x̂dnext, k̄, model, x̂d, û0, d0, model.p)
+    mul!(x̂snext, As, x̂s)
     x̂0next .+= f̂op .- x̂op
     return nothing
 end
@@ -135,42 +135,25 @@ Same than [`ĥ!`](@ref) for [`SimModel`](@ref) but without the `estim` argument
 """
 function ĥ!(ŷ0, model::SimModel, Cs_y::AbstractMatrix, x̂0, d0)
     # `@views` macro avoid copies with matrix slice operator e.g. [a:b]
-    @views xd, xs = x̂0[1:model.nx], x̂0[model.nx+1:end]
-    h!(ŷ0, model, xd, d0, model.p)  # y0 = h(xd, d0)
-    mul!(ŷ0, Cs_y, xs, 1, 1)        # ŷ0 = y0 + Cs_y*xs
+    @views x̂d, x̂s = x̂0[1:model.nx], x̂0[model.nx+1:end]
+    h!(ŷ0, model, x̂d, d0, model.p)  # y0 = h(xd, d0)
+    mul!(ŷ0, Cs_y, x̂s, 1, 1)        # ŷ0 = y0 + Cs_y*xs
     return nothing
 end
 
 """
-    fq_dae!(ẋ0, q0, model, x0, a0, u0, d0)
+    ĥ!(ŷ0, estim::StateEstimator, model::SimModel, x̂0, â0, d0)
 
-Call `model.fq!` for [`NonLinModelDAE`](@ref) or `model.f!` for [`NonLinModel`](@ref).
+Same than [`ĥ!`](@ref) for [`SimModel`](@ref), but call [`h_dae!`](@ref) internally.
 """
-function fq_dae!(ẋ0, q0, model::NonLinModelDAE, x0, a0, u0, d0)
-    return model.fq!(ẋ0, q0, x0, a0, u0, d0, model.p)
-end
-fq_dae!(ẋ0, _ , model::NonLinModel, x0, _ , u0, d0)= model.f!(ẋ0, x0, u0, d0, model.p)
-
-"""
-    ĥ_dae!(ŷ0, estim::StateEstimator, model::NonLinModelDAE, x0, a0, d0)
-
-Similar than [`ĥ!`](@ref) but with a algebraic variable `a0` for [`NonLinModelDAE`](@ref) .
-"""
-function ĥ_dae!(ŷ0, estim::StateEstimator, model::NonLinModelDAE, x̂0, a0, d0)
+function ĥ!(ŷ0, estim::StateEstimator, model::SimModel, x̂0, â0, d0)
     # `@views` macro avoid copies with matrix slice operator e.g. [a:b]
-    @views xd, xs = x̂0[1:model.nx], x̂0[model.nx+1:end]
-    model.h!(ŷ0, xd, a0, d0, model.p)
-    mul!(ŷ0, estim.Cs_y, xs, 1, 1)        # ŷ0 = y0 + Cs_y*xs
+    @views x̂d, x̂s = x̂0[1:model.nx], x̂0[model.nx+1:end]
+    h_dae!(ŷ0, model, x̂d, â0, d0)
+    mul!(ŷ0, estim.Cs_y, x̂s, 1, 1)        # ŷ0 = y0 + Cs_y*xs
     return nothing
 end
 
-"""
-    ĥ_dae!(ŷ0, estim::StateEstimator, model::SimModelODE, x0, _ , d0)
-
-Ignore the algebraic variable argument for other [`SimModelODE`](@ref) types.
-"""
-ĥ_dae!(ŷ0, estim::StateEstimator, model::SimModelODE, x̂0, _, d0) = ĥ!(ŷ0, estim, model, x̂0, d0)
-
 """
     disturbedinput!(Û0, estim::StateEstimator, x̂0, X̂0, U0) -> Û0
 

diff --git a/src/estimator/mhe.jl b/src/estimator/mhe.jl
@@ -3,7 +3,7 @@ include("mhe/execute.jl")
 include("mhe/transcription.jl")
 
 "Return estimation horizon He and slack variables length nε for `MovingHorizonEstimator`."
-get_other_dims(estim::MovingHorizonEstimator) = (estim.He, estim.nε)
+get_other_dims(estim::MovingHorizonEstimator) = (estim.He, estim.nε, get_na(estim.model))
 
 "Print optimizer and other information for `MovingHorizonEstimator`."
 function print_details(io::IO, estim::MovingHorizonEstimator)
@@ -27,11 +27,13 @@ print_backends(::IO, ::MovingHorizonEstimator, ::LinModel) = nothing
 function print_estim_dim(io::IO, estim::MovingHorizonEstimator, n; firstchars=nothing)
     nu, nd = estim.model.nu, estim.model.nd
     nx̂, nym, nyu = estim.nx̂, estim.nym, estim.nyu
+    na = get_na(estim.model)
     He, nε = estim.He, estim.nε
     niu, niym = sum(estim.nint_u), sum(estim.nint_ym)
     println(io, "  │ ├$(lpad(He, n)) estimation steps He")
     println(io, "  │ ├$(lpad(nu, n)) manipulated inputs u ($niu integrating states)")
     println(io, "  │ ├$(lpad(nx̂, n)) estimated states x̂")
+    println(io, "  │ ├$(lpad(na, n)) algebraic variables a")
     println(io, "  │ ├$(lpad(nym, n)) measured outputs ym ($niym integrating states)")
     println(io, "  │ ├$(lpad(nyu, n)) unmeasured outputs yu")
     print(io,   "  │ └$(lpad(nd, n)) measured disturbances d")