From 6f469ec8273200884336504a8922d42e25029fba Mon Sep 17 00:00:00 2001 From: franckgaga Date: Wed, 14 Oct 2026 10:16:54 -0400 Subject: [PATCH 1/6] changed: clearer notation for `NonLinModelDAE` sim --- src/model/nonlinmodeldae.jl | 46 ++++++++++++++++++------------------- 1 file changed, 22 insertions(+), 24 deletions(-) diff --git a/src/model/nonlinmodeldae.jl b/src/model/nonlinmodeldae.jl index 8b339a93d..a5b85bb1f 100644 --- a/src/model/nonlinmodeldae.jl +++ b/src/model/nonlinmodeldae.jl @@ -675,7 +675,7 @@ end """ update_predictions!(k̄, geq, model, Z) -TBW +Update `k̄` and `geq` in place for the open loop root solving of [`NonLinModelDAE`](@ref). """ function update_predictions!(k̄, geq, model, Z) x0, u0, d0 = model.x0_optim, model.u0_optim, model.d0_optim @@ -685,54 +685,52 @@ end """ con_nonlinprogeq!( - geq, k̄, model::NonLinModelDAE, ::TrapezoidalCollocation, x0, u0, d0, Z + geq, k̄dot, model::NonLinModelDAE, ::TrapezoidalCollocation, x0, u0, d0, Z ) -> geq -TBW +Nonlinear equality constraint for open loop root solving of [`NonLinModelDAE`](@ref). """ function con_nonlinprogeq!( - geq, k̄, model::NonLinModelDAE, ::TrapezoidalCollocation, x0, u0, d0, Z + geq, k̄dot, model::NonLinModelDAE, ::TrapezoidalCollocation, x0, u0, d0, Z ) nx, na = model.nx, model.na Ts = model.Ts - x0next_Z, aP_Z, a0_Z = @views Z[1:nx], Z[(nx+1):(nx+na)], Z[(nx+na+1):(nx+2na)] - sknext, q0, q1 = @views geq[1:nx], geq[(nx+1):(nx+na)], geq[(nx+na+1):(nx+2na)] - k̇0, k̇1 = @views k̄[1:nx], k̄[(nx+1):(2nx)] - model.fq!(k̇0, q0, x0, aP_Z, u0, d0, model.p) - model.fq!(k̇1, q1, x0next_Z, a0_Z, u0, d0, model.p) - sknext .= @. x0 - x0next_Z + 0.5*Ts*(k̇0 + k̇1) + x0next, a⁺, a0next = @views Z[1:nx], Z[(nx+1):(nx+na)], Z[(nx+na+1):(nx+2na)] + sknext, q⁺, q⁻next = @views geq[1:nx], geq[(nx+1):(nx+na)], geq[(nx+na+1):(nx+2na)] + k̇⁺, k̇⁻next = @views k̄dot[1:nx], k̄dot[(nx+1):(2nx)] + model.fq!(k̇⁺ , q⁺ , x0 , a⁺ , u0, d0, model.p) + model.fq!(k̇⁻next, q⁻next, x0next, a0next, u0, d0, model.p) + sknext .= @. x0 + 0.5*Ts*(k̇⁺ + k̇⁻next) - x0next return geq end """ con_nonlinprogeq!( - geq, k̄, model::NonLinModelDAE, transcription::OrthogonalCollocation, x0, u0, d0, Z + geq, k̄dot, model::NonLinModelDAE, transcription::OrthogonalCollocation, x0, u0, d0, Z ) -> geq -TBW +Same but using an [`OrthogonalCollocation`](@ref) transcription. """ function con_nonlinprogeq!( - geq, k̄, model::NonLinModelDAE, transcription::OrthogonalCollocation, x0, u0, d0, Z + geq, k̄dot, model::NonLinModelDAE, transcription::OrthogonalCollocation, x0, u0, d0, Z ) nx, na = model.nx, model.na Mo, no = model.Mo, transcription.no nk̄, nā = get_nk̄(model, transcription), get_nā(model, transcription) - aP_Z, k̄_Z, ā_Z = @views Z[(nx+1):(nx+na)], Z[(nx+na+1):(nx+na+nk̄)], Z[(nx+na+nk̄+1):end] - q0, sk̄, q̄ = @views geq[1:na], geq[(na+1):(na+nk̄)], geq[(na+nk̄+1):(na+nk̄+nā)] - @views model.fq!(k̄[1:nx], q0, x0, aP_Z, u0, d0, model.p) - Δk = k̄ + a⁺, k̄, ā = @views Z[(nx+1):(nx+na)], Z[(nx+na+1):(nx+na+nk̄)], Z[(nx+na+nk̄+1):end] + q⁺, sk̄, q̄ = @views geq[1:na], geq[(na+1):(na+nk̄)], geq[(na+nk̄+1):(na+nk̄+nā)] + @views model.fq!(k̄dot[1:nx], q⁺, x0, a⁺, u0, d0, model.p) + Δk = k̄dot for i=1:no - Δk[(1 + (i-1)*nx):(i*nx)] = @views k̄_Z[(1 + (i-1)*nx):(i*nx)] .- x0 + Δk[(1 + (i-1)*nx):(i*nx)] = @views k̄[(1 + (i-1)*nx):(i*nx)] .- x0 end mul!(sk̄, Mo, Δk) for i=1:no - k̇i = @views k̄[(1 + (i-1)*nx):(i*nx)] - qi = @views q̄[(1 + (i-1)*na):(i*na)] - ki_Z = @views k̄_Z[(1 + (i-1)*nx):(i*nx)] - ai_Z = @views ā_Z[(1 + (i-1)*na):(i*na)] - model.fq!(k̇i, qi, ki_Z, ai_Z, u0, d0, model.p) + k̇i, ki = @views k̄dot[(1 + (i-1)*nx):(i*nx)], k̄[(1 + (i-1)*nx):(i*nx)] + qi, ai = @views q̄[(1 + (i-1)*na):(i*na)], ā[(1 + (i-1)*na):(i*na)] + model.fq!(k̇i, qi, ki, ai, u0, d0, model.p) end - sk̄ .-= k̄ + sk̄ .-= k̄dot return geq end From cd960529f08a0b44cc67c13af01d57552601fd18 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Wed, 14 Oct 2026 10:27:55 -0400 Subject: [PATCH 2/6] doc: precision in `getinfo(::NonLinModelDAE)` --- src/model/nonlinmodeldae.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/model/nonlinmodeldae.jl b/src/model/nonlinmodeldae.jl index a5b85bb1f..c0bd589ee 100644 --- a/src/model/nonlinmodeldae.jl +++ b/src/model/nonlinmodeldae.jl @@ -877,10 +877,10 @@ The function should be called after calling [`updatestate!`](@ref) on `model` ob returns the dictionary `info` with the following fields: - `:xnext` : next state, ``\mathbf{x}(k+1)`` -- `:q` : current algebraic equation residuals `res`, ``\mathbf{q(x, a, u, d, p)}`` +- `:a` : right limit of current algebraic variable, ``\mathbf{a^+}(k)`` +- `:q` : current algebraic equation residuals `res`, ``\mathbf{q(x, a^+, u, d, p)}`` - `:y` : current output, ``\mathbf{y}(k)`` - `:x` : current state, ``\mathbf{x}(k)`` -- `:a` : current algebraic variable, ``\mathbf{a}(k)`` - `:u` : current manipulated input, ``\mathbf{u}(k)`` - `:d` : current measured disturbances, ``\mathbf{u}(k)`` From f67cf6446efc20a5a8adba057f85825b1e07b64b Mon Sep 17 00:00:00 2001 From: franckgaga Date: Wed, 14 Oct 2026 14:49:19 -0400 Subject: [PATCH 3/6] bench: re-enabling `UnoSolver` with `OrthogonalCollocation` The issue is solved here : https://github.com/cvanaret/Uno/issues/802 --- benchmark/3_bench_predictive_control.jl | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/benchmark/3_bench_predictive_control.jl b/benchmark/3_bench_predictive_control.jl index cf5f225ec..11ada2a58 100644 --- a/benchmark/3_bench_predictive_control.jl +++ b/benchmark/3_bench_predictive_control.jl @@ -359,12 +359,11 @@ nmpc_madnlp_ms_hess = NonLinMPC(estim; Hp, Hc, Mwt, Nwt, Cwt, optim, transcripti nmpc_madnlp_ms_hess = setconstraint!(nmpc_madnlp_ms_hess; umin, umax) JuMP.unset_time_limit_sec(nmpc_madnlp_ms_hess.optim) -# # commenting out this case study, it crashes on latest UnoSolver v0.3.4 -# optim = JuMP.Model(()->UnoSolver.Optimizer(preset="funnelsqp"), add_bridges=false) -# transcription, hessian = OrthogonalCollocation(), true -# nmpc_uno_oc_hess = NonLinMPC(estim; Hp, Hc, Mwt, Nwt, Cwt, optim, transcription, hessian) -# nmpc_uno_oc_hess = setconstraint!(nmpc_uno_oc_hess; umin, umax) -# JuMP.unset_time_limit_sec(nmpc_uno_oc_hess.optim) +optim = JuMP.Model(()->UnoSolver.Optimizer(preset="funnelsqp"), add_bridges=false) +transcription, hessian = OrthogonalCollocation(), true +nmpc_uno_oc_hess = NonLinMPC(estim; Hp, Hc, Mwt, Nwt, Cwt, optim, transcription, hessian) +nmpc_uno_oc_hess = setconstraint!(nmpc_uno_oc_hess; umin, umax) +JuMP.unset_time_limit_sec(nmpc_uno_oc_hess.optim) samples, evals, seconds = 100, 1, 15*60 CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["Ipopt"]["SingleShooting"] = @@ -432,11 +431,11 @@ CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["Uno"]["MultipleShooting (Hessi sim!($nmpc_uno_ms_hess, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false), samples=samples, evals=evals, seconds=seconds, setup=GC.gc() ) -# CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["Uno"]["OrthogonalCollocation (Hessian)"] = -# @benchmarkable( -# sim!($nmpc_uno_oc_hess, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false), -# samples=samples, evals=evals, seconds=seconds, setup=GC.gc() -# ) +CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["Uno"]["OrthogonalCollocation (Hessian)"] = + @benchmarkable( + sim!($nmpc_uno_oc_hess, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false), + samples=samples, evals=evals, seconds=seconds, setup=GC.gc() + ) # ----------------- Case study: Pendulum economic -------------------------------- model2, p = pendulum_model2, pendulum_p2 From bda81c1d3825bdaa0b099e4db99dacbc12738032 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Wed, 14 Oct 2026 15:48:40 -0400 Subject: [PATCH 4/6] doc: details on ZOH with `NonLinModelDAE` simulations --- src/model/nonlinmodeldae.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/model/nonlinmodeldae.jl b/src/model/nonlinmodeldae.jl index c0bd589ee..f8f2a7db4 100644 --- a/src/model/nonlinmodeldae.jl +++ b/src/model/nonlinmodeldae.jl @@ -240,6 +240,13 @@ NonLinModelDAE with a sample time Ts = 5.0 s: are no global direct transmissions from ``\mathbf{u}`` to ``\mathbf{y}`` with the functions ``\mathbf{q}`` and ``\mathbf{h}`` using [`SparseConnectivityTracer.jl`](@extref SparseConnectivityTracer.jl). + Note that open loop simulations discretize the manipulated inputs ``\mathbf{u}`` and + measured disturbances ``\mathbf{d}`` with a zero-order hold (piecewise constants), since + the values for the next time step cannot be fed to the simple [`updatestate!`](@ref) + method. It induces a small plant-model mismatch for the inputs if `h=1` in the + transcription method of the estimator/controller, and also for the disturbances, always + discretized as piecewise linear functions in the estimator/controller. + By default, a dense [`ForwardDiff`](@extref ForwardDiff) backend is used for the Jacobians of the nonlinear equality constraints. This is also the default backend for the Hessians if `hessian=true`. From 017833e1a89d77097a956aed2c9b03424989d4a5 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Wed, 14 Oct 2026 15:48:53 -0400 Subject: [PATCH 5/6] doc: minor addition --- src/transcription.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/transcription.jl b/src/transcription.jl index 3447c5561..deb8fbe7a 100644 --- a/src/transcription.jl +++ b/src/transcription.jl @@ -150,8 +150,8 @@ for descriptions of `f_threads` and `h_threads` keywords. The Extended Help deta the decision variables for [`NonLinModelDAE`](@ref). The `h` argument is `0` or `1`, for piecewise constant or linear manipulated inputs ``\mathbf{u}`` (`h=1` is slightly less expensive). Note that the various [`DiffSolver`](@ref) here assume zero-order hold, so `h=1` -will induce a plant-model mismatch if the plant is simulated with these solvers. Measured -disturbances ``\mathbf{d}`` are piecewise linear. +will induce a small plant-model mismatch if the plant is simulated with these solvers. +Measured disturbances ``\mathbf{d}`` are piecewise linear. This transcription computes the predictions by calling the continuous-time model in the equality constraint function and by using the implicit trapezoidal rule. It can handle From aee4b6ca7aee23c5e81bf73c6b48c5b206af5bef Mon Sep 17 00:00:00 2001 From: franckgaga Date: Wed, 14 Oct 2026 15:49:28 -0400 Subject: [PATCH 6/6] bump --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 6f212ff5c..022359586 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "ModelPredictiveControl" uuid = "61f9bdb8-6ae4-484a-811f-bbf86720c31c" -version = "2.14.3" +version = "2.14.4" authors = ["Francis Gagnon"] [deps]