From b6a148efc153e833e50172dac2c67e9a4c6616f1 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Sat, 23 Aug 2025 13:08:00 -0400 Subject: [PATCH 01/19] WIP: `TrapezoidalMethod()` collocation --- docs/src/public/predictive_control.md | 6 +++++ src/ModelPredictiveControl.jl | 2 +- src/controller/transcription.jl | 36 +++++++++++++++++++++++++++ 3 files changed, 43 insertions(+), 1 deletion(-) diff --git a/docs/src/public/predictive_control.md b/docs/src/public/predictive_control.md index 439b585ef..2434e4ff5 100644 --- a/docs/src/public/predictive_control.md +++ b/docs/src/public/predictive_control.md @@ -109,3 +109,9 @@ SingleShooting ```@docs MultipleShooting ``` + +### TrapezoidalMethod + +```@docs +TrapezoidalMethod +``` diff --git a/src/ModelPredictiveControl.jl b/src/ModelPredictiveControl.jl index 38b7623bd..c5b003507 100644 --- a/src/ModelPredictiveControl.jl +++ b/src/ModelPredictiveControl.jl @@ -39,7 +39,7 @@ export MovingHorizonEstimator export ManualEstimator export default_nint, initstate! export PredictiveController, ExplicitMPC, LinMPC, NonLinMPC, setconstraint!, moveinput! -export TranscriptionMethod, SingleShooting, MultipleShooting +export TranscriptionMethod, SingleShooting, MultipleShooting, TrapezoidalMethod export SimResult, getinfo, sim! include("general.jl") diff --git a/src/controller/transcription.jl b/src/controller/transcription.jl index ba9f92103..b1fddc048 100644 --- a/src/controller/transcription.jl +++ b/src/controller/transcription.jl @@ -61,6 +61,42 @@ function get_nZ(estim::StateEstimator, ::MultipleShooting, Hp, Hc) return estim.model.nu*Hc + estim.nx̂*Hp end + +@doc raw""" + TrapezoidalMethod() + +An implicit trapezoidal transcription method (not yet implemented). + +This is presumably the simplest collocation method. It can handle moderately stiff systems +and is A-stable. However, it may not be as efficient as more advanced methods for highly +stiff systems. The decision variables are the same as for [`MultipleShooting`](@ref), hence +a similar algorithm complexity. + +# Extended Help + +!!! details "Extended Help" + The trapezoidal method estimates the defects with: + by: + ```math + \mathbf{Ŝ}(k) = \mathbf{Ẽ_ŝ Z̃} + \mathbf{K_ŝ x̂_0}(k) + 0.5\frac{T_s}\big(\mathbf{F̂}(k+1) + \mathbf{F̂}(k)\big) + ``` + where ``T_s`` is the sampling period, ``\mathbf{Ẽ}`` the matrix defined at + [`init_defectmat`](@ref), and ``\mathbf{F̂}(k+j)`` the stacked vector of system + + where ``\mathbf{f̂}(k+j) = \mathbf{f̂}\big(\mathbf{x̂}(k+j), \mathbf{u}(k+j), \mathbf{d}(k+j)\big)``. + This leads to the following defect constraints for ``j=0`` to ``H_p-1``: + ```math + \mathbf{ŝ}(k+j) = \mathbf{x̂}(k+j+1) - \mathbf{x̂}(k+j) - \frac{T_s}{2} \big( \mathbf{f̂}(k+j) + \mathbf{f̂}(k+j+1) \big) = 0 + ``` + which are added as equality constraints in the optimization problem. The initial state + ``\mathbf{x̂}(k)`` is given by the state estimator, and the future states + ``\mathbf{x̂}(k+j+1)`` are decision variables in the optimization problem. The method + requires evaluating the system dynamics at both the current and next time steps, which + can increase computational complexity compared to explicit methods like single shooting. +""" +struct TrapezoidalMethod <: TranscriptionMethod end + + @doc raw""" init_ZtoΔU(estim::StateEstimator, transcription::TranscriptionMethod, Hp, Hc) -> PΔu From 0689b5c6d7141cf3534aaba7850a9ab543131d4b Mon Sep 17 00:00:00 2001 From: franckgaga Date: Wed, 27 Aug 2025 15:07:21 -0400 Subject: [PATCH 02/19] changed: removed useless unpacking in transcriptions --- src/controller/transcription.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/controller/transcription.jl b/src/controller/transcription.jl index b1fddc048..4959a1644 100644 --- a/src/controller/transcription.jl +++ b/src/controller/transcription.jl @@ -1124,8 +1124,7 @@ function predict!( mpc::PredictiveController, model::NonLinModel, ::SingleShooting, U0, _ ) - nu, nx̂, ny, nd, nk = model.nu, mpc.estim.nx̂, model.ny, model.nd, model.nk - Hp, Hc = mpc.Hp, mpc.Hc + nu, nx̂, ny, nd, nk, Hp = model.nu, mpc.estim.nx̂, model.ny, model.nd, model.nk, mpc.Hp D̂0 = mpc.D̂0 x̂0 = @views mpc.estim.x̂0[1:nx̂] d0 = @views mpc.d0[1:nd] @@ -1209,7 +1208,7 @@ custom constraints are include in the `g` vector. function con_nonlinprog!( g, mpc::PredictiveController, ::NonLinModel, ::TranscriptionMethod, x̂0end, Ŷ0, gc, ϵ ) - nx̂, nŶ = length(x̂0end), length(Ŷ0) + nŶ = length(Ŷ0) for i in eachindex(g) mpc.con.i_g[i] || continue if i ≤ nŶ @@ -1276,7 +1275,7 @@ function con_nonlinprogeq!( geq, X̂0, Û0, K0, mpc::PredictiveController, model::NonLinModel, ::MultipleShooting, U0, Z̃ ) - nu, nx̂, ny, nd, nk = model.nu, mpc.estim.nx̂, model.ny, model.nd, model.nk + nu, nx̂, nd, nk = model.nu, mpc.estim.nx̂, model.nd, model.nk Hp, Hc = mpc.Hp, mpc.Hc nΔU, nX̂ = nu*Hc, nx̂*Hp D̂0 = mpc.D̂0 From df10988fb57478007e4ba7aac836b8cbbaab81c3 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Wed, 27 Aug 2025 16:21:35 -0400 Subject: [PATCH 03/19] changed: removed the hat for states in aug. model --- src/estimator/execute.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/estimator/execute.jl b/src/estimator/execute.jl index 790e4350a..0022b52d6 100644 --- a/src/estimator/execute.jl +++ b/src/estimator/execute.jl @@ -80,12 +80,12 @@ Same than [`f̂!`](@ref) for [`SimModel`](@ref) but without the `estim` argument """ function f̂!(x̂0next, û0, k0, model::SimModel, As, Cs_u, x̂0, u0, d0) # `@views` macro avoid copies with matrix slice operator e.g. [a:b] - @views x̂d, x̂s = x̂0[1:model.nx], x̂0[model.nx+1:end] - @views x̂d_next, x̂s_next = x̂0next[1:model.nx], x̂0next[model.nx+1:end] - mul!(û0, Cs_u, x̂s) # ŷs_u = Cs_u * x̂s - û0 .+= u0 - f!(x̂d_next, k0, model, x̂d, û0, d0, model.p) - mul!(x̂s_next, As, x̂s) + @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!(û0, Cs_u, xs) # ys_u = Cs_u*xs + û0 .+= u0 # û0 = u0 + ys_u + f!(xdnext, k0, model, xd, û0, d0, model.p) + mul!(xsnext, As, xs) return nothing end @@ -116,9 +116,9 @@ Same than [`ĥ!`](@ref) for [`SimModel`](@ref) but without the `estim` argument """ function ĥ!(ŷ0, model::SimModel, Cs_y, x̂0, d0) # `@views` macro avoid copies with matrix slice operator e.g. [a:b] - @views x̂d, x̂s = x̂0[1:model.nx], x̂0[model.nx+1:end] - h!(ŷ0, model, x̂d, d0, model.p) - mul!(ŷ0, Cs_y, x̂s, 1, 1) # ŷ0 = ŷ0 + Cs_y*x̂s + @views xd, xs = x̂0[1:model.nx], x̂0[model.nx+1:end] + h!(ŷ0, model, xd, d0, model.p) # y0 = h(xd, d0) + mul!(ŷ0, Cs_y, xs, 1, 1) # ŷ0 = y0 + Cs_y*xs return nothing end From 75c83cce2c3dd7872a9f5769b60f512240dd1a24 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Wed, 27 Aug 2025 18:36:22 -0400 Subject: [PATCH 04/19] changed: rename `TrapezoidalMethod` to `TrapezoidalCollocation` --- src/controller/transcription.jl | 70 ++++++++++++++++++++++++++++----- 1 file changed, 61 insertions(+), 9 deletions(-) diff --git a/src/controller/transcription.jl b/src/controller/transcription.jl index 4959a1644..0402aff94 100644 --- a/src/controller/transcription.jl +++ b/src/controller/transcription.jl @@ -63,7 +63,7 @@ end @doc raw""" - TrapezoidalMethod() + TrapezoidalCollocation() An implicit trapezoidal transcription method (not yet implemented). @@ -1284,19 +1284,71 @@ function con_nonlinprogeq!( d0 = @views mpc.d0[1:nd] #TODO: allow parallel for loop or threads? for j=1:Hp - u0 = @views U0[(1 + nu*(j-1)):(nu*j)] - û0 = @views Û0[(1 + nu*(j-1)):(nu*j)] - x̂0next = @views X̂0[(1 + nx̂*(j-1)):(nx̂*j)] - k0 = @views K0[(1 + nk*(j-1)):(nk*j)] + u0 = @views U0[(1 + nu*(j-1)):(nu*j)] + û0 = @views Û0[(1 + nu*(j-1)):(nu*j)] + k0 = @views K0[(1 + nk*(j-1)):(nk*j)] + d0next = @views D̂0[(1 + nd*(j-1)):(nd*j)] + x̂0next = @views X̂0[(1 + nx̂*(j-1)):(nx̂*j)] + x̂0next_Z̃ = @views X̂0_Z̃[(1 + nx̂*(j-1)):(nx̂*j)] + ŝnext = @views geq[(1 + nx̂*(j-1)):(nx̂*j)] f̂!(x̂0next, û0, k0, mpc.estim, model, x̂0, u0, d0) x̂0next .+= mpc.estim.f̂op .- mpc.estim.x̂op - x̂0next_Z̃ = @views X̂0_Z̃[(1 + nx̂*(j-1)):(nx̂*j)] - ŝnext = @views geq[(1 + nx̂*(j-1)):(nx̂*j)] - ŝnext .= x̂0next .- x̂0next_Z̃ + ŝnext .= x̂0next .- x̂0next_Z̃ x̂0 = x̂0next_Z̃ # using states in Z̃ for next iteration (allow parallel for) - d0 = @views D̂0[(1 + nd*(j-1)):(nd*j)] + d0 = d0next + end + return geq +end + + +function con_nonlinprogeq!( + geq, X̂0, Û0, K0, + mpc::PredictiveController, model::NonLinModel, ::TrapezoidalCollocation, U0, Z̃ +) + nu, nx̂, nd, nx = model.nu, mpc.estim.nx̂, model.nd, model.nx + Hp, Hc = mpc.Hp, mpc.Hc + nΔU, nX̂ = nu*Hc, nx̂*Hp + + Ts, p = model.Ts, model.p + As = mpc.estim.As + # K0 stores the state derivatives at each collocation point for all intervals + nk = model.nk # TODO: initialize K0 to the adequate size for a given collocation method + D̂0 = mpc.D̂0 + X̂0_Z̃ = @views Z̃[(nΔU+1):(nΔU+nX̂)] + x̂0 = @views mpc.estim.x̂0[1:nx̂] + d0 = @views mpc.d0[1:nd] + #TODO: allow parallel for loop or threads? + for j=1:Hp + u0 = @views U0[(1 + nu*(j-1)):(nu*j)] + û0 = @views Û0[(1 + nu*(j-1)):(nu*j)] + k0 = @views K0[(1 + nk*(j-1)):(nk*j)] + d0next = @views D̂0[(1 + nd*(j-1)):(nd*j)] + x̂0next = @views X̂0[(1 + nx̂*(j-1)):(nx̂*j)] + x̂0next_Z̃ = @views X̂0_Z̃[(1 + nx̂*(j-1)):(nx̂*j)] + ŝnext = @views geq[(1 + nx̂*(j-1)):(nx̂*j)] + xd, xs = @views x̂0[1:nx], x̂0[nx+1:end] + xdnext_Z̃, xsnext_Z̃ = @views x̂0next_Z̃[1:nx], x̂0next_Z̃[nx+1:end] + sdnext, ssnext = @views ŝnext[1:nx], ŝnext[nx+1:end] + + mul!(û0, Cs_u, xs) # ys_u = Cs_u*xs + û0 .+= u0 # û0 = u0 + ys_u + + ẋd, ẋdnext = @views k0[1:nx], k0[nx+1:2*nx] + # TODO: decide what to do with model.xop and model.fop + model.f!(ẋd, xd, û0, d0, p) + model.f!(ẋdnext, xdnext_Z̃, û0, d0next, p) # with ZOH on manipulated inputs u + + xsnext = @views x̂0next[nx+1:end] + mul!(xsnext, As, xs) + + sdnext .= @. xd - xdnext_Z̃ + (Ts/2)*(ẋd + ẋdnext) + ssnext .= @. xsnext - xsnext_Z̃ + + x̂0 = x̂0next_Z̃ # using states in Z̃ for next iteration (allow parallel for) + d0 = d0next end return geq end + con_nonlinprogeq!(geq,_,_,_,::PredictiveController,::NonLinModel,::SingleShooting, _,_)=geq con_nonlinprogeq!(geq,_,_,_,::PredictiveController,::LinModel,::TranscriptionMethod,_,_)=geq \ No newline at end of file From 7175b02a224986dc0830ce554958161d2a48f0f0 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Thu, 28 Aug 2025 08:24:41 -0400 Subject: [PATCH 05/19] debug: rename the transcription everywhere --- docs/src/public/predictive_control.md | 4 ++-- src/ModelPredictiveControl.jl | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/public/predictive_control.md b/docs/src/public/predictive_control.md index 2434e4ff5..1e73553de 100644 --- a/docs/src/public/predictive_control.md +++ b/docs/src/public/predictive_control.md @@ -110,8 +110,8 @@ SingleShooting MultipleShooting ``` -### TrapezoidalMethod +### TrapezoidalCollocation ```@docs -TrapezoidalMethod +TrapezoidalCollocation ``` diff --git a/src/ModelPredictiveControl.jl b/src/ModelPredictiveControl.jl index c5b003507..679768bfa 100644 --- a/src/ModelPredictiveControl.jl +++ b/src/ModelPredictiveControl.jl @@ -39,7 +39,7 @@ export MovingHorizonEstimator export ManualEstimator export default_nint, initstate! export PredictiveController, ExplicitMPC, LinMPC, NonLinMPC, setconstraint!, moveinput! -export TranscriptionMethod, SingleShooting, MultipleShooting, TrapezoidalMethod +export TranscriptionMethod, SingleShooting, MultipleShooting, TrapezoidalCollocation export SimResult, getinfo, sim! include("general.jl") From 08a09e32f9bebfa1f49656a7d3cc99ba75396b86 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Thu, 28 Aug 2025 08:28:05 -0400 Subject: [PATCH 06/19] debug: idem --- src/controller/transcription.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/controller/transcription.jl b/src/controller/transcription.jl index 0402aff94..695aa15b1 100644 --- a/src/controller/transcription.jl +++ b/src/controller/transcription.jl @@ -94,7 +94,7 @@ a similar algorithm complexity. requires evaluating the system dynamics at both the current and next time steps, which can increase computational complexity compared to explicit methods like single shooting. """ -struct TrapezoidalMethod <: TranscriptionMethod end +struct TrapezoidalCollocation <: TranscriptionMethod end @doc raw""" From 55e4aa9d67e04b8c55e00a29acedc96a8597d4fa Mon Sep 17 00:00:00 2001 From: franckgaga Date: Thu, 28 Aug 2025 10:17:21 -0400 Subject: [PATCH 07/19] added: `ShootingMethod` and `CollocationMethod` abstract types --- src/controller/construct.jl | 6 +- src/controller/linmpc.jl | 12 +-- src/controller/nonlinmpc.jl | 2 +- src/controller/transcription.jl | 185 ++++++++++++++++---------------- 4 files changed, 104 insertions(+), 101 deletions(-) diff --git a/src/controller/construct.jl b/src/controller/construct.jl index 711cb3040..08c251900 100644 --- a/src/controller/construct.jl +++ b/src/controller/construct.jl @@ -600,11 +600,11 @@ end ex̂, fx̂, gx̂, jx̂, kx̂, vx̂, bx̂, Eŝ, Fŝ, Gŝ, Jŝ, Kŝ, Vŝ, Bŝ, gc!=nothing, nc=0 - ) -> con, nϵ, P̃Δu, P̃u, Ẽ, Ẽŝ + ) -> con, nϵ, P̃Δu, P̃u, Ẽ Init `ControllerConstraint` struct with default parameters based on estimator `estim`. -Also return `P̃Δu`, `P̃u`, `Ẽ` and `Ẽŝ` matrices for the the augmented decision vector `Z̃`. +Also return `P̃Δu`, `P̃u` and `Ẽ` matrices for the the augmented decision vector `Z̃`. """ function init_defaultcon_mpc( estim::StateEstimator{NT}, @@ -660,7 +660,7 @@ function init_defaultcon_mpc( C_ymin , C_ymax , c_x̂min , c_x̂max , i_g, gc! , nc ) - return con, nϵ, P̃Δu, P̃u, Ẽ, Ẽŝ + return con, nϵ, P̃Δu, P̃u, Ẽ end "Repeat predictive controller constraints over prediction `Hp` and control `Hc` horizons." diff --git a/src/controller/linmpc.jl b/src/controller/linmpc.jl index edc906d0a..49f54cbc1 100644 --- a/src/controller/linmpc.jl +++ b/src/controller/linmpc.jl @@ -5,7 +5,7 @@ struct LinMPC{ NT<:Real, SE<:StateEstimator, CW<:ControllerWeights, - TM<:TranscriptionMethod, + TM<:ShootingMethod, JM<:JuMP.GenericModel } <: PredictiveController{NT} estim::SE @@ -54,7 +54,7 @@ struct LinMPC{ NT<:Real, SE<:StateEstimator, CW<:ControllerWeights, - TM<:TranscriptionMethod, + TM<:ShootingMethod, JM<:JuMP.GenericModel } model = estim.model @@ -71,7 +71,7 @@ struct LinMPC{ Eŝ, Gŝ, Jŝ, Kŝ, Vŝ, Bŝ = init_defectmat(model, estim, transcription, Hp, Hc) # dummy vals (updated just before optimization): F, fx̂, Fŝ = zeros(NT, ny*Hp), zeros(NT, nx̂), zeros(NT, nx̂*Hp) - con, nϵ, P̃Δu, P̃u, Ẽ, Ẽŝ = init_defaultcon_mpc( + con, nϵ, P̃Δu, P̃u, Ẽ = init_defaultcon_mpc( estim, weights, transcription, Hp, Hc, PΔu, Pu, E, @@ -156,7 +156,7 @@ arguments. This controller allocates memory at each time step for the optimizati - `N_Hc=Diagonal(repeat(Nwt,Hc))` : positive semidefinite symmetric matrix ``\mathbf{N}_{H_c}``. - `L_Hp=Diagonal(repeat(Lwt,Hp))` : positive semidefinite symmetric matrix ``\mathbf{L}_{H_p}``. - `Cwt=1e5` : slack variable weight ``C`` (scalar), use `Cwt=Inf` for hard constraints only. -- `transcription=SingleShooting()` : a [`TranscriptionMethod`](@ref) for the optimization. +- `transcription=SingleShooting()` : [`SingleShooting`](@ref) or [`MultipleShooting`](@ref). - `optim=JuMP.Model(OSQP.MathOptInterfaceOSQP.Optimizer)` : quadratic optimizer used in the predictive controller, provided as a [`JuMP.Model`](@extref) object (default to [`OSQP`](https://osqp.org/docs/parsers/jump.html) optimizer). @@ -215,7 +215,7 @@ function LinMPC( N_Hc = Diagonal(repeat(Nwt, get_Hc(move_blocking(Hp, Hc)))), L_Hp = Diagonal(repeat(Lwt, Hp)), Cwt = DEFAULT_CWT, - transcription::TranscriptionMethod = DEFAULT_LINMPC_TRANSCRIPTION, + transcription::ShootingMethod = DEFAULT_LINMPC_TRANSCRIPTION, optim::JuMP.GenericModel = JuMP.Model(DEFAULT_LINMPC_OPTIMIZER, add_bridges=false), kwargs... ) @@ -259,7 +259,7 @@ function LinMPC( N_Hc = Diagonal(repeat(Nwt, get_Hc(move_blocking(Hp, Hc)))), L_Hp = Diagonal(repeat(Lwt, Hp)), Cwt = DEFAULT_CWT, - transcription::TranscriptionMethod = DEFAULT_LINMPC_TRANSCRIPTION, + transcription::ShootingMethod = DEFAULT_LINMPC_TRANSCRIPTION, optim::JM = JuMP.Model(DEFAULT_LINMPC_OPTIMIZER, add_bridges=false), ) 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 index 5b5d2050a..8f8d4839d 100644 --- a/src/controller/nonlinmpc.jl +++ b/src/controller/nonlinmpc.jl @@ -94,7 +94,7 @@ struct NonLinMPC{ Eŝ, Gŝ, Jŝ, Kŝ, Vŝ, Bŝ = init_defectmat(model, estim, transcription, Hp, Hc) # dummy vals (updated just before optimization): F, fx̂, Fŝ = zeros(NT, ny*Hp), zeros(NT, nx̂), zeros(NT, nx̂*Hp) - con, nϵ, P̃Δu, P̃u, Ẽ, Ẽŝ = init_defaultcon_mpc( + con, nϵ, P̃Δu, P̃u, Ẽ = init_defaultcon_mpc( estim, weights, transcription, Hp, Hc, PΔu, Pu, E, diff --git a/src/controller/transcription.jl b/src/controller/transcription.jl index 695aa15b1..7ade0a714 100644 --- a/src/controller/transcription.jl +++ b/src/controller/transcription.jl @@ -1,9 +1,12 @@ """ Abstract supertype of all transcription methods of [`PredictiveController`](@ref). -The module currently supports [`SingleShooting`](@ref) and [`MultipleShooting`](@ref). +The module currently supports [`SingleShooting`](@ref), [`MultipleShooting`](@ref) and +[`TrapezoidalCollocation`](@ref) transcription methods. """ abstract type TranscriptionMethod end +abstract type ShootingMethod <: TranscriptionMethod end +abstract type CollocationMethod <: TranscriptionMethod end @doc raw""" SingleShooting() @@ -22,7 +25,7 @@ any custom move blocking): This method is generally more efficient for small control horizon ``H_c``, stable and mildly nonlinear plant model/constraints. """ -struct SingleShooting <: TranscriptionMethod end +struct SingleShooting <: ShootingMethod end @doc raw""" MultipleShooting() @@ -51,13 +54,13 @@ control horizon ``H_c``, unstable or highly nonlinear plant models/constraints. Sparse optimizers like `OSQP` or `Ipopt` and sparse Jacobian computations are recommended for this transcription method. """ -struct MultipleShooting <: TranscriptionMethod end +struct MultipleShooting <: ShootingMethod end "Get the number of elements in the optimization decision vector `Z`." function get_nZ(estim::StateEstimator, ::SingleShooting, Hp, Hc) return estim.model.nu*Hc end -function get_nZ(estim::StateEstimator, ::MultipleShooting, Hp, Hc) +function get_nZ(estim::StateEstimator, ::TranscriptionMethod, Hp, Hc) return estim.model.nu*Hc + estim.nx̂*Hp end @@ -94,7 +97,7 @@ a similar algorithm complexity. requires evaluating the system dynamics at both the current and next time steps, which can increase computational complexity compared to explicit methods like single shooting. """ -struct TrapezoidalCollocation <: TranscriptionMethod end +struct TrapezoidalCollocation <: CollocationMethod end @doc raw""" @@ -129,7 +132,7 @@ function init_ZtoΔU( end function init_ZtoΔU( - estim::StateEstimator{NT}, ::MultipleShooting, Hp, Hc + estim::StateEstimator{NT}, ::TranscriptionMethod, Hp, Hc ) where {NT<:Real} I_nu_Hc = sparse(Matrix{NT}(I, estim.model.nu*Hc, estim.model.nu*Hc)) PΔu = [I_nu_Hc spzeros(NT, estim.model.nu*Hc, estim.nx̂*Hp)] @@ -205,7 +208,7 @@ end init_PUmat( _ , ::SingleShooting, _ , _ , PuDagger) = PuDagger function init_PUmat( - estim, ::MultipleShooting, Hp, _ , PuDagger::AbstractMatrix{NT} + estim, ::TranscriptionMethod, Hp, _ , PuDagger::AbstractMatrix{NT} ) where NT<:Real return [PuDagger spzeros(NT, estim.model.nu*Hp, estim.nx̂*Hp)] end @@ -425,21 +428,13 @@ function init_predmat( return E, G, J, K, V, B, ex̂, gx̂, jx̂, kx̂, vx̂, bx̂ end -@doc raw""" - init_predmat(model::SimModel, estim, transcription::MultipleShooting, Hp, Hc) - -Return the terminal state matrices for [`SimModel`](@ref) and [`MultipleShooting`](@ref). - -The output prediction matrices are all empty matrices. The terminal state matrices are -given in the Extended Help section. +""" + init_predmat(model::NonLinModel, estim, transcription::SingleShooting, Hp, Hc) -# Extended Help -!!! details "Extended Help" - The terminal state matrices all appropriately sized zero matrices ``\mathbf{0}``, except - for ``\mathbf{e_x̂} = [\begin{smallmatrix}\mathbf{0} & \mathbf{I}\end{smallmatrix}]`` +Return empty matrices for [`SingleShooting`](@ref) of [`NonLinModel`](@ref) """ function init_predmat( - model::SimModel, estim::StateEstimator{NT}, transcription::MultipleShooting, Hp, Hc + model::NonLinModel, estim::StateEstimator{NT}, transcription::SingleShooting, Hp, Hc ) where {NT<:Real} nu, nx̂, nd = model.nu, estim.nx̂, model.nd nZ = get_nZ(estim, transcription, Hp, Hc) @@ -449,22 +444,25 @@ function init_predmat( K = zeros(NT, 0, nx̂) V = zeros(NT, 0, nu) B = zeros(NT, 0) - ex̂ = [zeros(NT, nx̂, Hc*nu + (Hp-1)*nx̂) I] - gx̂ = zeros(NT, nx̂, nd) - jx̂ = zeros(NT, nx̂, nd*Hp) - kx̂ = zeros(NT, nx̂, nx̂) - vx̂ = zeros(NT, nx̂, nu) - bx̂ = zeros(NT, nx̂) + ex̂, gx̂, jx̂, kx̂, vx̂, bx̂ = E, G, J, K, V, B return E, G, J, K, V, B, ex̂, gx̂, jx̂, kx̂, vx̂, bx̂ end -""" - init_predmat(model::SimModel, estim, transcription::TranscriptionMethod, Hp, Hc) +@doc raw""" + init_predmat(model::NonLinModel, estim, transcription::TranscriptionMethod, Hp, Hc) + +Return the terminal state matrices for [`NonLinModel`](@ref) and other [`TranscriptionMethod`](@ref). -Return empty matrices for all other cases. +The output prediction matrices are all empty matrices. The terminal state matrices are +given in the Extended Help section. + +# Extended Help +!!! details "Extended Help" + The terminal state matrices all appropriately sized zero matrices ``\mathbf{0}``, except + for ``\mathbf{e_x̂} = [\begin{smallmatrix}\mathbf{0} & \mathbf{I}\end{smallmatrix}]`` """ function init_predmat( - model::SimModel, estim::StateEstimator{NT}, transcription::TranscriptionMethod, Hp, Hc + model::NonLinModel, estim::StateEstimator{NT}, transcription::TranscriptionMethod, Hp, Hc ) where {NT<:Real} nu, nx̂, nd = model.nu, estim.nx̂, model.nd nZ = get_nZ(estim, transcription, Hp, Hc) @@ -474,7 +472,12 @@ function init_predmat( K = zeros(NT, 0, nx̂) V = zeros(NT, 0, nu) B = zeros(NT, 0) - ex̂, gx̂, jx̂, kx̂, vx̂, bx̂ = E, G, J, K, V, B + ex̂ = [zeros(NT, nx̂, Hc*nu + (Hp-1)*nx̂) I] + gx̂ = zeros(NT, nx̂, nd) + jx̂ = zeros(NT, nx̂, nd*Hp) + kx̂ = zeros(NT, nx̂, nx̂) + vx̂ = zeros(NT, nx̂, nu) + bx̂ = zeros(NT, nx̂) return E, G, J, K, V, B, ex̂, gx̂, jx̂, kx̂, vx̂, bx̂ end @@ -615,42 +618,42 @@ function init_matconstraint_mpc( return i_b, i_g, A, Aeq, neq end -"Init `i_b` without output constraints if `NonLinModel` and not `SingleShooting`." +"Init `i_b` without output & terminal constraints if `NonLinModel` and `SingleShooting`." function init_matconstraint_mpc( - ::NonLinModel{NT}, ::TranscriptionMethod, nc::Int, + ::NonLinModel{NT}, ::SingleShooting, nc::Int, i_Umin, i_Umax, i_ΔŨmin, i_ΔŨmax, i_Ymin, i_Ymax, i_x̂min, i_x̂max, args... ) where {NT<:Real} if isempty(args) A, Aeq, neq = nothing, nothing, nothing else - A_Umin, A_Umax, A_ΔŨmin, A_ΔŨmax, A_Ymin, A_Ymax, A_x̂min, A_x̂max, A_ŝ = args - A = [A_Umin; A_Umax; A_ΔŨmin; A_ΔŨmax; A_Ymin; A_Ymax; A_x̂min; A_x̂max] + A_Umin, A_Umax, A_ΔŨmin, A_ΔŨmax, _ , _ , _ , _ , A_ŝ = args + A = [A_Umin; A_Umax; A_ΔŨmin; A_ΔŨmax] Aeq = A_ŝ - nΔŨ, nZ̃ = size(A_ΔŨmin) - neq = nZ̃ - nΔŨ + neq = 0 end - i_b = [i_Umin; i_Umax; i_ΔŨmin; i_ΔŨmax; i_x̂min; i_x̂max] - i_g = [i_Ymin; i_Ymax; trues(nc)] + i_b = [i_Umin; i_Umax; i_ΔŨmin; i_ΔŨmax] + i_g = [i_Ymin; i_Ymax; i_x̂min; i_x̂max; trues(nc)] return i_b, i_g, A, Aeq, neq end -"Init `i_b` without output & terminal constraints if `NonLinModel` and `SingleShooting`." +"Init `i_b` without output constraints if `NonLinModel` and other `TranscriptionMethod`." function init_matconstraint_mpc( - ::NonLinModel{NT}, ::SingleShooting, nc::Int, + ::NonLinModel{NT}, ::TranscriptionMethod, nc::Int, i_Umin, i_Umax, i_ΔŨmin, i_ΔŨmax, i_Ymin, i_Ymax, i_x̂min, i_x̂max, args... ) where {NT<:Real} if isempty(args) A, Aeq, neq = nothing, nothing, nothing - else - A_Umin, A_Umax, A_ΔŨmin, A_ΔŨmax, A_Ymin, A_Ymax, A_x̂min, A_x̂max, A_ŝ = args - A = [A_Umin; A_Umax; A_ΔŨmin; A_ΔŨmax; A_Ymin; A_Ymax; A_x̂min; A_x̂max] + else + A_Umin, A_Umax, A_ΔŨmin, A_ΔŨmax, _ , _ , A_x̂min, A_x̂max, A_ŝ = args + A = [A_Umin; A_Umax; A_ΔŨmin; A_ΔŨmax; A_x̂min; A_x̂max] Aeq = A_ŝ - neq = 0 + nΔŨ, nZ̃ = size(A_ΔŨmin) + neq = nZ̃ - nΔŨ end - i_b = [i_Umin; i_Umax; i_ΔŨmin; i_ΔŨmax] - i_g = [i_Ymin; i_Ymax; i_x̂min; i_x̂max; trues(nc)] + i_b = [i_Umin; i_Umax; i_ΔŨmin; i_ΔŨmax; i_x̂min; i_x̂max] + i_g = [i_Ymin; i_Ymax; trues(nc)] return i_b, i_g, A, Aeq, neq end @@ -686,24 +689,21 @@ end """ init_nonlincon!( - mpc::PredictiveController, model::NonLinModel, transcription::MultipleShooting, + mpc::PredictiveController, model::NonLinModel, ::SingleShooting, gfuncs, ∇gfuncs!, geqfuncs, ∇geqfuncs! ) - -Init nonlinear constraints for [`NonLinModel`](@ref) and [`MultipleShooting`](@ref). -The nonlinear constraints are the output prediction `Ŷ` bounds, the custom inequality -constraints `gc` and all the nonlinear equality constraints `geq`. +Init nonlinear constraints for [`NonLinModel`](@ref) and [`SingleShooting`](@ref). + +The nonlinear constraints are the custom inequality constraints `gc`, the output +prediction `Ŷ` bounds and the terminal state `x̂end` bounds. """ function init_nonlincon!( - mpc::PredictiveController, ::NonLinModel, ::MultipleShooting, - gfuncs, ∇gfuncs!, - geqfuncs, ∇geqfuncs! -) + mpc::PredictiveController, ::NonLinModel, ::SingleShooting, gfuncs, ∇gfuncs!, _ , _ +) optim, con = mpc.optim, mpc.con ny, nx̂, Hp, nZ̃ = mpc.estim.model.ny, mpc.estim.nx̂, mpc.Hp, length(mpc.Z̃) - # --- nonlinear inequality constraints --- if length(con.i_g) ≠ 0 i_base = 0 for i in eachindex(con.Y0min) @@ -720,6 +720,20 @@ function init_nonlincon!( ) end i_base = 2Hp*ny + for i in eachindex(con.x̂0min) + name = Symbol("g_x̂0min_$i") + optim[name] = JuMP.add_nonlinear_operator( + optim, nZ̃, gfuncs[i_base+i], ∇gfuncs![i_base+i]; name + ) + end + i_base = 2Hp*ny + nx̂ + for i in eachindex(con.x̂0max) + name = Symbol("g_x̂0max_$i") + optim[name] = JuMP.add_nonlinear_operator( + optim, nZ̃, gfuncs[i_base+i], ∇gfuncs![i_base+i]; name + ) + end + i_base = 2Hp*ny + 2nx̂ for i in 1:con.nc name = Symbol("g_c_$i") optim[name] = JuMP.add_nonlinear_operator( @@ -727,37 +741,29 @@ function init_nonlincon!( ) end end - # --- nonlinear equality constraints --- - Z̃var = optim[:Z̃var] - for i in eachindex(geqfuncs) - name = Symbol("geq_$i") - geqfunc_i = optim[name] = JuMP.add_nonlinear_operator( - optim, nZ̃, geqfuncs[i], ∇geqfuncs![i]; name - ) - # set with @constrains here instead of set_nonlincon!, since the number of nonlinear - # equality constraints is known and constant (±Inf are impossible): - @constraint(optim, geqfunc_i(Z̃var...) == 0) - end return nothing end """ init_nonlincon!( - mpc::PredictiveController, model::NonLinModel, ::SingleShooting, + mpc::PredictiveController, model::NonLinModel, transcription::TranscriptionMethod, gfuncs, ∇gfuncs!, geqfuncs, ∇geqfuncs! ) + +Init nonlinear constraints for [`NonLinModel`](@ref) and other [`TranscriptionMethod`](@ref). -Init nonlinear constraints for [`NonLinModel`](@ref) and [`SingleShooting`](@ref). - -The nonlinear constraints are the custom inequality constraints `gc`, the output -prediction `Ŷ` bounds and the terminal state `x̂end` bounds. +The nonlinear constraints are the output prediction `Ŷ` bounds, the custom inequality +constraints `gc` and all the nonlinear equality constraints `geq`. """ function init_nonlincon!( - mpc::PredictiveController, ::NonLinModel, ::SingleShooting, gfuncs, ∇gfuncs!, _ , _ + mpc::PredictiveController, ::NonLinModel, ::TranscriptionMethod, + gfuncs, ∇gfuncs!, + geqfuncs, ∇geqfuncs! ) optim, con = mpc.optim, mpc.con ny, nx̂, Hp, nZ̃ = mpc.estim.model.ny, mpc.estim.nx̂, mpc.Hp, length(mpc.Z̃) + # --- nonlinear inequality constraints --- if length(con.i_g) ≠ 0 i_base = 0 for i in eachindex(con.Y0min) @@ -774,20 +780,6 @@ function init_nonlincon!( ) end i_base = 2Hp*ny - for i in eachindex(con.x̂0min) - name = Symbol("g_x̂0min_$i") - optim[name] = JuMP.add_nonlinear_operator( - optim, nZ̃, gfuncs[i_base+i], ∇gfuncs![i_base+i]; name - ) - end - i_base = 2Hp*ny + nx̂ - for i in eachindex(con.x̂0max) - name = Symbol("g_x̂0max_$i") - optim[name] = JuMP.add_nonlinear_operator( - optim, nZ̃, gfuncs[i_base+i], ∇gfuncs![i_base+i]; name - ) - end - i_base = 2Hp*ny + 2nx̂ for i in 1:con.nc name = Symbol("g_c_$i") optim[name] = JuMP.add_nonlinear_operator( @@ -795,6 +787,17 @@ function init_nonlincon!( ) end end + # --- nonlinear equality constraints --- + Z̃var = optim[:Z̃var] + for i in eachindex(geqfuncs) + name = Symbol("geq_$i") + geqfunc_i = optim[name] = JuMP.add_nonlinear_operator( + optim, nZ̃, geqfuncs[i], ∇geqfuncs![i]; name + ) + # set with @constrains here instead of set_nonlincon!, since the number of nonlinear + # equality constraints is known and constant (±Inf are impossible): + @constraint(optim, geqfunc_i(Z̃var...) == 0) + end return nothing end @@ -1148,17 +1151,17 @@ end @doc raw""" predict!( Ŷ0, x̂0end, _ , _ , _ , - mpc::PredictiveController, model::NonLinModel, transcription::MultipleShooting, + mpc::PredictiveController, model::NonLinModel, transcription::TranscriptionMethod, U0, Z̃ ) -> Ŷ0, x̂0end -Compute vectors if `model` is a [`NonLinModel`](@ref) and for [`MultipleShooting`](@ref). +Compute vectors if `model` is a [`NonLinModel`](@ref) and other [`TranscriptionMethod`](@ref). The method mutates `Ŷ0` and `x̂0end` arguments. """ function predict!( Ŷ0, x̂0end, _, _, _, - mpc::PredictiveController, model::NonLinModel, ::MultipleShooting, + mpc::PredictiveController, model::NonLinModel, ::TranscriptionMethod, U0, Z̃ ) nu, ny, nd, nx̂, Hp, Hc = model.nu, model.ny, model.nd, mpc.estim.nx̂, mpc.Hp, mpc.Hc @@ -1310,7 +1313,7 @@ function con_nonlinprogeq!( nΔU, nX̂ = nu*Hc, nx̂*Hp Ts, p = model.Ts, model.p - As = mpc.estim.As + As, Cs_u = mpc.estim.As, mpc.estim.Cs_u # K0 stores the state derivatives at each collocation point for all intervals nk = model.nk # TODO: initialize K0 to the adequate size for a given collocation method D̂0 = mpc.D̂0 From f9c520276fe940945e83ab02595ef62111608280 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Thu, 28 Aug 2025 10:38:22 -0400 Subject: [PATCH 08/19] added: validate `model` and `transcription` combo --- src/controller/explicitmpc.jl | 1 + src/controller/linmpc.jl | 1 + src/controller/nonlinmpc.jl | 1 + src/controller/transcription.jl | 14 +++++++++++++- 4 files changed, 16 insertions(+), 1 deletion(-) diff --git a/src/controller/explicitmpc.jl b/src/controller/explicitmpc.jl index 95621197e..49cc6b950 100644 --- a/src/controller/explicitmpc.jl +++ b/src/controller/explicitmpc.jl @@ -50,6 +50,7 @@ struct ExplicitMPC{ R̂y, R̂u, Tu_lastu0 = zeros(NT, ny*Hp), zeros(NT, nu*Hp), zeros(NT, nu*Hp) lastu0 = zeros(NT, nu) transcription = SingleShooting() # explicit MPC only supports SingleShooting + validate_transcription(model, transcription) PΔu = init_ZtoΔU(estim, transcription, Hp, Hc) Pu, Tu = init_ZtoU(estim, transcription, Hp, Hc, nb) E, G, J, K, V, B = init_predmat(model, estim, transcription, Hp, Hc) diff --git a/src/controller/linmpc.jl b/src/controller/linmpc.jl index 49f54cbc1..f6818daad 100644 --- a/src/controller/linmpc.jl +++ b/src/controller/linmpc.jl @@ -63,6 +63,7 @@ struct LinMPC{ # dummy vals (updated just before optimization): R̂y, R̂u, Tu_lastu0 = zeros(NT, ny*Hp), zeros(NT, nu*Hp), zeros(NT, nu*Hp) lastu0 = zeros(NT, nu) + validate_transcription(model, transcription) PΔu = init_ZtoΔU(estim, transcription, Hp, Hc) Pu, Tu = init_ZtoU(estim, transcription, Hp, Hc, nb) E, G, J, K, V, B, ex̂, gx̂, jx̂, kx̂, vx̂, bx̂ = init_predmat( diff --git a/src/controller/nonlinmpc.jl b/src/controller/nonlinmpc.jl index 8f8d4839d..4e601b9d4 100644 --- a/src/controller/nonlinmpc.jl +++ b/src/controller/nonlinmpc.jl @@ -86,6 +86,7 @@ struct NonLinMPC{ # dummy vals (updated just before optimization): R̂y, R̂u, Tu_lastu0 = zeros(NT, ny*Hp), zeros(NT, nu*Hp), zeros(NT, nu*Hp) lastu0 = zeros(NT, nu) + validate_transcription(model, transcription) PΔu = init_ZtoΔU(estim, transcription, Hp, Hc) Pu, Tu = init_ZtoU(estim, transcription, Hp, Hc, nb) E, G, J, K, V, B, ex̂, gx̂, jx̂, kx̂, vx̂, bx̂ = init_predmat( diff --git a/src/controller/transcription.jl b/src/controller/transcription.jl index 7ade0a714..7763a690b 100644 --- a/src/controller/transcription.jl +++ b/src/controller/transcription.jl @@ -100,6 +100,16 @@ a similar algorithm complexity. struct TrapezoidalCollocation <: CollocationMethod end +function validate_transcription(::LinModel, ::CollocationMethod) + throw(ArgumentError("Collocation methods are not supported for LinModel.")) + return nothing +end +function validate_transcription(::NonLinModel{<:Real, <:EmptySolver}, ::CollocationMethod) + throw(ArgumentError("Collocation methods require continuous-time NonLinModel.")) + return nothing +end +validate_transcription(::SimModel, ::TranscriptionMethod) = nothing + @doc raw""" init_ZtoΔU(estim::StateEstimator, transcription::TranscriptionMethod, Hp, Hc) -> PΔu @@ -206,7 +216,9 @@ function init_ZtoU( return Pu, Tu end -init_PUmat( _ , ::SingleShooting, _ , _ , PuDagger) = PuDagger +function init_PUmat(_,::SingleShooting,_,_,PuDagger::AbstractMatrix{NT}) where NT<:Real + return PuDagger +end function init_PUmat( estim, ::TranscriptionMethod, Hp, _ , PuDagger::AbstractMatrix{NT} ) where NT<:Real From c213281f3b6ae5d3af5bca37341438c531782178 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Thu, 28 Aug 2025 10:57:09 -0400 Subject: [PATCH 09/19] added: show `TranscriptionMethod` when printing `PredictiveController` --- src/predictive_control.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/predictive_control.jl b/src/predictive_control.jl index e88cab888..6c9c35f7b 100644 --- a/src/predictive_control.jl +++ b/src/predictive_control.jl @@ -33,6 +33,7 @@ function Base.show(io::IO, mpc::PredictiveController) n = maximum(ndigits.((Hp, Hc, nu, nx̂, nym, nyu, nd))) + 1 println(io, "$(typeof(mpc).name.name) controller with a sample time Ts = "* "$(mpc.estim.model.Ts) s, $(JuMP.solver_name(mpc.optim)) optimizer, "* + "$(typeof(mpc.transcription).name.name) transcription, "* "$(typeof(mpc.estim).name.name) estimator and:") println(io, "$(lpad(Hp, n)) prediction steps Hp") println(io, "$(lpad(Hc, n)) control steps Hc") From 781f43b8f8040e9f00d287fce802eedf75080fd9 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Thu, 28 Aug 2025 11:08:10 -0400 Subject: [PATCH 10/19] doc: `jldoctest` modifications --- src/controller/construct.jl | 2 +- src/controller/linmpc.jl | 4 ++-- src/controller/nonlinmpc.jl | 4 ++-- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/controller/construct.jl b/src/controller/construct.jl index 711cb3040..012ad5954 100644 --- a/src/controller/construct.jl +++ b/src/controller/construct.jl @@ -221,7 +221,7 @@ constraints are all soft by default. See Extended Help for time-varying constrai julia> mpc = LinMPC(setop!(LinModel(tf(3, [30, 1]), 4), uop=[50], yop=[25])); julia> mpc = setconstraint!(mpc, umin=[0], umax=[100], Δumin=[-10], Δumax=[+10]) -LinMPC controller with a sample time Ts = 4.0 s, OSQP optimizer, SteadyKalmanFilter estimator and: +LinMPC controller with a sample time Ts = 4.0 s, OSQP optimizer, SingleShooting transcription, SteadyKalmanFilter estimator and: 10 prediction steps Hp 2 control steps Hc 1 slack variable ϵ (control constraints) diff --git a/src/controller/linmpc.jl b/src/controller/linmpc.jl index edc906d0a..3b1fc1c95 100644 --- a/src/controller/linmpc.jl +++ b/src/controller/linmpc.jl @@ -167,7 +167,7 @@ arguments. This controller allocates memory at each time step for the optimizati julia> model = LinModel([tf(3, [30, 1]); tf(-2, [5, 1])], 4); julia> mpc = LinMPC(model, Mwt=[0, 1], Nwt=[0.5], Hp=30, Hc=1) -LinMPC controller with a sample time Ts = 4.0 s, OSQP optimizer, SteadyKalmanFilter estimator and: +LinMPC controller with a sample time Ts = 4.0 s, OSQP optimizer, SingleShooting transcription, SteadyKalmanFilter estimator and: 30 prediction steps Hp 1 control steps Hc 1 slack variable ϵ (control constraints) @@ -237,7 +237,7 @@ Use custom state estimator `estim` to construct `LinMPC`. julia> estim = KalmanFilter(LinModel([tf(3, [30, 1]); tf(-2, [5, 1])], 4), i_ym=[2]); julia> mpc = LinMPC(estim, Mwt=[0, 1], Nwt=[0.5], Hp=30, Hc=1) -LinMPC controller with a sample time Ts = 4.0 s, OSQP optimizer, KalmanFilter estimator and: +LinMPC controller with a sample time Ts = 4.0 s, OSQP optimizer, SingleShooting transcription, KalmanFilter estimator and: 30 prediction steps Hp 1 control steps Hc 1 slack variable ϵ (control constraints) diff --git a/src/controller/nonlinmpc.jl b/src/controller/nonlinmpc.jl index 5b5d2050a..7b950f88a 100644 --- a/src/controller/nonlinmpc.jl +++ b/src/controller/nonlinmpc.jl @@ -223,7 +223,7 @@ This controller allocates memory at each time step for the optimization. julia> model = NonLinModel((x,u,_,_)->0.5x+u, (x,_,_)->2x, 10.0, 1, 1, 1, solver=nothing); julia> mpc = NonLinMPC(model, Hp=20, Hc=1, Cwt=1e6) -NonLinMPC controller with a sample time Ts = 10.0 s, Ipopt optimizer, UnscentedKalmanFilter estimator and: +NonLinMPC controller with a sample time Ts = 10.0 s, Ipopt optimizer, SingleShooting transcription, UnscentedKalmanFilter estimator and: 20 prediction steps Hp 1 control steps Hc 1 slack variable ϵ (control constraints) @@ -327,7 +327,7 @@ julia> model = NonLinModel((x,u,_,_)->0.5x+u, (x,_,_)->2x, 10.0, 1, 1, 1, solver julia> estim = UnscentedKalmanFilter(model, σQint_ym=[0.05]); julia> mpc = NonLinMPC(estim, Hp=20, Hc=1, Cwt=1e6) -NonLinMPC controller with a sample time Ts = 10.0 s, Ipopt optimizer, UnscentedKalmanFilter estimator and: +NonLinMPC controller with a sample time Ts = 10.0 s, Ipopt optimizer, SingleShooting transcription, UnscentedKalmanFilter estimator and: 20 prediction steps Hp 1 control steps Hc 1 slack variable ϵ (control constraints) From c34531fd8f1b2d3c1a88d90db4b9961c5d525776 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Thu, 28 Aug 2025 11:18:12 -0400 Subject: [PATCH 11/19] doc: mention defaults in manual --- docs/src/manual/linmpc.md | 5 +++-- docs/src/manual/nonlinmpc.md | 6 ++++-- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/docs/src/manual/linmpc.md b/docs/src/manual/linmpc.md index 69edb189b..8d11aa9fa 100644 --- a/docs/src/manual/linmpc.md +++ b/docs/src/manual/linmpc.md @@ -78,8 +78,9 @@ mpc = setconstraint!(mpc, ymin=[48, -Inf]) in which `Hp` and `Hc` keyword arguments are respectively the predictive and control horizons, and `Mwt` and `Nwt`, the output setpoint tracking and move suppression weights. By -default, [`LinMPC`](@ref) controllers use [`OSQP`](https://osqp.org/) to solve the problem, -soft constraints on output predictions ``\mathbf{ŷ}`` to ensure feasibility, and a +default, [`LinMPC`](@ref) controllers use [`OSQP`](https://osqp.org/) and a direct +[`SingleShooting`](@ref) transcription method to solve the optimal control problem, soft +constraints on output predictions ``\mathbf{ŷ}`` to ensure feasibility, and a [`SteadyKalmanFilter`](@ref) to estimate the plant states[^1]. An attentive reader will also notice that the Kalman filter estimates two additional states compared to the plant model. These are the integrating states for the unmeasured plant disturbances, and they are diff --git a/docs/src/manual/nonlinmpc.md b/docs/src/manual/nonlinmpc.md index 1c5032b3b..8f8f8102f 100644 --- a/docs/src/manual/nonlinmpc.md +++ b/docs/src/manual/nonlinmpc.md @@ -119,8 +119,10 @@ umin, umax = [-1.5], [+1.5] nmpc = setconstraint!(nmpc; umin, umax) ``` -The option `Cwt=Inf` disables the slack variable `ϵ` for constraint softening. We test `mpc` -performance on `plant` by imposing an angular setpoint of 180° (inverted position): +The option `Cwt=Inf` disables the slack variable `ϵ` for constraint softening. By default, +[`NonLinMPC`](@ref) controllers use [`Ipopt`](https://coin-or.github.io/Ipopt/) and a direct +[`SingleShooting`](@ref) transcription method to solve the optimal control problem. We test +`mpc` performance on `plant` by imposing an angular setpoint of 180° (inverted position): ```@example man_nonlin using JuMP; unset_time_limit_sec(nmpc.optim) # hide From cf8a1cb4b0d6e2b6269cd11793ad93502ab2b766 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Thu, 28 Aug 2025 13:48:12 -0400 Subject: [PATCH 12/19] doc: details on operating points and continuous dynamics --- src/controller/transcription.jl | 90 +++++++++++++++++++-------------- 1 file changed, 53 insertions(+), 37 deletions(-) diff --git a/src/controller/transcription.jl b/src/controller/transcription.jl index 7763a690b..ac3bf74fd 100644 --- a/src/controller/transcription.jl +++ b/src/controller/transcription.jl @@ -56,32 +56,25 @@ for this transcription method. """ struct MultipleShooting <: ShootingMethod end -"Get the number of elements in the optimization decision vector `Z`." -function get_nZ(estim::StateEstimator, ::SingleShooting, Hp, Hc) - return estim.model.nu*Hc -end -function get_nZ(estim::StateEstimator, ::TranscriptionMethod, Hp, Hc) - return estim.model.nu*Hc + estim.nx̂*Hp -end - - @doc raw""" TrapezoidalCollocation() -An implicit trapezoidal transcription method (not yet implemented). +Construct an implicit trapezoidal [`TranscriptionMethod`](@ref). -This is presumably the simplest collocation method. It can handle moderately stiff systems -and is A-stable. However, it may not be as efficient as more advanced methods for highly -stiff systems. The decision variables are the same as for [`MultipleShooting`](@ref), hence -a similar algorithm complexity. +This is the simplest collocation method. It is only supported for continuous-time +[`NonLinModel`](@ref)s. It can handle moderately stiff systems and is A-stable. The decision +variables are the same as for [`MultipleShooting`](@ref), hence similar computational costs. +However, it may not be as efficient as more advanced collocation methods for highly stiff +systems. It currently assumes piecewise constant manipulated inputs (or zero-order hold) +between the samples, but linear interpolation will be added in the future. See Extended Help +for more details on the collocation/defect constraints. # Extended Help !!! details "Extended Help" - The trapezoidal method estimates the defects with: - by: + The implicit trapezoidal collocation estimates the defects with: ```math - \mathbf{Ŝ}(k) = \mathbf{Ẽ_ŝ Z̃} + \mathbf{K_ŝ x̂_0}(k) + 0.5\frac{T_s}\big(\mathbf{F̂}(k+1) + \mathbf{F̂}(k)\big) + \mathbf{s_d}(k) = \mathbf{x_0} \frac{T_s}{2}\big(\mathbf{F̂}(k+1) + \mathbf{F̂}(k)\big) ``` where ``T_s`` is the sampling period, ``\mathbf{Ẽ}`` the matrix defined at [`init_defectmat`](@ref), and ``\mathbf{F̂}(k+j)`` the stacked vector of system @@ -99,7 +92,6 @@ a similar algorithm complexity. """ struct TrapezoidalCollocation <: CollocationMethod end - function validate_transcription(::LinModel, ::CollocationMethod) throw(ArgumentError("Collocation methods are not supported for LinModel.")) return nothing @@ -110,6 +102,14 @@ function validate_transcription(::NonLinModel{<:Real, <:EmptySolver}, ::Collocat end validate_transcription(::SimModel, ::TranscriptionMethod) = nothing +"Get the number of elements in the optimization decision vector `Z`." +function get_nZ(estim::StateEstimator, ::SingleShooting, Hp, Hc) + return estim.model.nu*Hc +end +function get_nZ(estim::StateEstimator, ::TranscriptionMethod, Hp, Hc) + return estim.model.nu*Hc + estim.nx̂*Hp +end + @doc raw""" init_ZtoΔU(estim::StateEstimator, transcription::TranscriptionMethod, Hp, Hc) -> PΔu @@ -1015,11 +1015,10 @@ function linconstrainteq!(mpc::PredictiveController, model::LinModel, ::Multiple JuMP.set_normalized_rhs(linconeq, mpc.con.beq) return nothing end -linconstrainteq!(::PredictiveController, ::SimModel, ::SingleShooting) = nothing -linconstrainteq!(::PredictiveController, ::SimModel, ::MultipleShooting) = nothing +linconstrainteq!(::PredictiveController, ::SimModel, ::TranscriptionMethod) = nothing @doc raw""" - set_warmstart!(mpc::PredictiveController, transcription::SingleShooting, Z̃var) -> Z̃s + set_warmstart!(mpc::PredictiveController, ::SingleShooting, Z̃var) -> Z̃s Set and return the warm-start value of `Z̃var` for [`SingleShooting`](@ref) transcription. @@ -1050,9 +1049,9 @@ function set_warmstart!(mpc::PredictiveController, ::SingleShooting, Z̃var) end @doc raw""" - set_warmstart!(mpc::PredictiveController, transcription::MultipleShooting, Z̃var) -> Z̃s + set_warmstart!(mpc::PredictiveController, ::TranscriptionMethod, Z̃var) -> Z̃s -Set and return the warm-start value of `Z̃var` for [`MultipleShooting`](@ref) transcription. +Set and return the warm-start value of `Z̃var` for other [`TranscriptionMethod`](@ref). It warm-starts the solver at: ```math @@ -1075,7 +1074,7 @@ where ``\mathbf{x̂_0}(k+j|k-1)`` is the predicted state for time ``k+j`` comput last control period ``k-1``, expressed as a deviation from the operating point ``\mathbf{x̂_{op}}``. """ -function set_warmstart!(mpc::PredictiveController, ::MultipleShooting, Z̃var) +function set_warmstart!(mpc::PredictiveController, ::TranscriptionMethod, Z̃var) nu, nx̂, Hp, Hc, Z̃s = mpc.estim.model.nu, mpc.estim.nx̂, mpc.Hp, mpc.Hc, mpc.buffer.Z̃ # --- input increments ΔU --- Z̃s[1:(Hc*nu-nu)] .= @views mpc.Z̃[nu+1:Hc*nu] @@ -1307,6 +1306,7 @@ function con_nonlinprogeq!( x̂0next_Z̃ = @views X̂0_Z̃[(1 + nx̂*(j-1)):(nx̂*j)] ŝnext = @views geq[(1 + nx̂*(j-1)):(nx̂*j)] f̂!(x̂0next, û0, k0, mpc.estim, model, x̂0, u0, d0) + # handle operating points (but should be zeros for NonLinModel): x̂0next .+= mpc.estim.f̂op .- mpc.estim.x̂op ŝnext .= x̂0next .- x̂0next_Z̃ x̂0 = x̂0next_Z̃ # using states in Z̃ for next iteration (allow parallel for) @@ -1315,7 +1315,30 @@ function con_nonlinprogeq!( return geq end +@docs raw""" + con_nonlinprogeq!( + geq, X̂0, Û0, K0 + mpc::PredictiveController, model::NonLinModel, ::TrapezoidalCollocation, U0, Z̃ + ) + +Nonlinear equality constrains for [`NonLinModel`](@ref) and [`TrapezoidalCollocation`](@ref). +The method mutates the `geq`, `X̂0`, `Û0` and `K0` vectors in argument. The deterministic +and stochastic states are handled separately since collocation methods require +continuous-time state-space models, and the stochastic model of the unmeasured disturbances +is discrete-time. Also note that operating points in `model` are typically zeros for +[`NonLinModel`](@ref), but they are handled rigorously here if it's not the case. It should +be noted that linearization of continuous-time dynamics at non-equilibrium points leads to: +```math + \mathbf{ẋ_0}(t) ≈ \mathbf{A x_0}(t) + \mathbf{B_u u_0}(t) + \mathbf{B_d d_0}(t) +``` +as opposed to, for discrete-time models: +```math + \mathbf{x_0}(k+1) ≈ \mathbf{A x_0}(k) + \mathbf{B_u u_0}(k) + \mathbf{B_d d_0}(k) + + \mathbf{f_{op} - \mathbf{x_{op}} +``` +hence no need to add `model.fop` and subtract `model.xop` in the collocation equations. +""" function con_nonlinprogeq!( geq, X̂0, Û0, K0, mpc::PredictiveController, model::NonLinModel, ::TrapezoidalCollocation, U0, Z̃ @@ -1344,26 +1367,19 @@ function con_nonlinprogeq!( xd, xs = @views x̂0[1:nx], x̂0[nx+1:end] xdnext_Z̃, xsnext_Z̃ = @views x̂0next_Z̃[1:nx], x̂0next_Z̃[nx+1:end] sdnext, ssnext = @views ŝnext[1:nx], ŝnext[nx+1:end] - mul!(û0, Cs_u, xs) # ys_u = Cs_u*xs û0 .+= u0 # û0 = u0 + ys_u - - ẋd, ẋdnext = @views k0[1:nx], k0[nx+1:2*nx] - # TODO: decide what to do with model.xop and model.fop - model.f!(ẋd, xd, û0, d0, p) - model.f!(ẋdnext, xdnext_Z̃, û0, d0next, p) # with ZOH on manipulated inputs u - + ẋ0, ẋ0next = @views k0[1:nx], k0[nx+1:2*nx] + # no need to handle model.fop and model.xop, see docstring: + model.f!(ẋ0, xd, û0, d0, p) + model.f!(ẋ0next, xdnext_Z̃, û0, d0next, p) # assuming ZOH on manipulated inputs u xsnext = @views x̂0next[nx+1:end] mul!(xsnext, As, xs) - - sdnext .= @. xd - xdnext_Z̃ + (Ts/2)*(ẋd + ẋdnext) + sdnext .= @. xd - xdnext_Z̃ + (Ts/2)*(ẋ0 + ẋ0next) ssnext .= @. xsnext - xsnext_Z̃ - x̂0 = x̂0next_Z̃ # using states in Z̃ for next iteration (allow parallel for) d0 = d0next end return geq end - -con_nonlinprogeq!(geq,_,_,_,::PredictiveController,::NonLinModel,::SingleShooting, _,_)=geq -con_nonlinprogeq!(geq,_,_,_,::PredictiveController,::LinModel,::TranscriptionMethod,_,_)=geq \ No newline at end of file +con_nonlinprogeq!(geq,_,_,_,::PredictiveController,::SimModel,::TranscriptionMethod,_,_)=geq \ No newline at end of file From 1548b213e9d97dfd0971f97f87fa2ad71667800e Mon Sep 17 00:00:00 2001 From: franckgaga Date: Thu, 28 Aug 2025 13:49:24 -0400 Subject: [PATCH 13/19] doc: minor modification --- src/controller/transcription.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/controller/transcription.jl b/src/controller/transcription.jl index ac3bf74fd..02a5cd86b 100644 --- a/src/controller/transcription.jl +++ b/src/controller/transcription.jl @@ -1335,7 +1335,7 @@ be noted that linearization of continuous-time dynamics at non-equilibrium point as opposed to, for discrete-time models: ```math \mathbf{x_0}(k+1) ≈ \mathbf{A x_0}(k) + \mathbf{B_u u_0}(k) + \mathbf{B_d d_0}(k) - + \mathbf{f_{op} - \mathbf{x_{op}} + + \mathbf{f_{op}} - \mathbf{x_{op}} ``` hence no need to add `model.fop` and subtract `model.xop` in the collocation equations. """ From ab55d5c9bb9999a8a2d20dc3cf4365a45dcb8afb Mon Sep 17 00:00:00 2001 From: franckgaga Date: Thu, 28 Aug 2025 13:53:58 -0400 Subject: [PATCH 14/19] doc: debug --- src/controller/transcription.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/controller/transcription.jl b/src/controller/transcription.jl index 02a5cd86b..880cc7c8d 100644 --- a/src/controller/transcription.jl +++ b/src/controller/transcription.jl @@ -1315,7 +1315,7 @@ function con_nonlinprogeq!( return geq end -@docs raw""" +@doc raw""" con_nonlinprogeq!( geq, X̂0, Û0, K0 mpc::PredictiveController, model::NonLinModel, ::TrapezoidalCollocation, U0, Z̃ From 356e573a5d76365ce5f71de9a78791626ccafa6e Mon Sep 17 00:00:00 2001 From: franckgaga Date: Thu, 28 Aug 2025 14:56:02 -0400 Subject: [PATCH 15/19] doc: comment on the stochastic model with `TrapezoidalCollocation` --- src/controller/transcription.jl | 51 ++++++++++++++------------------- 1 file changed, 22 insertions(+), 29 deletions(-) diff --git a/src/controller/transcription.jl b/src/controller/transcription.jl index 880cc7c8d..a6fe8290a 100644 --- a/src/controller/transcription.jl +++ b/src/controller/transcription.jl @@ -22,8 +22,11 @@ any custom move blocking): \vdots \\ \mathbf{Δu}(k+H_c-1) \end{bmatrix} ``` -This method is generally more efficient for small control horizon ``H_c``, stable and mildly -nonlinear plant model/constraints. +This method computes the predictions by calling the augmented discrete-time model +recursively over the prediction horizon ``H_p`` in the objective function, or by updating +the linear coefficients of the quadratic optimization for [`LinModel`](@ref). It is +generally more efficient for small control horizon ``H_c``, stable and mildly nonlinear +plant model/constraints. """ struct SingleShooting <: ShootingMethod end @@ -48,7 +51,11 @@ operating point ``\mathbf{x̂_{op}}`` (see [`augment_model`](@ref)): where ``\mathbf{x̂}_i(k+j)`` is the state prediction for time ``k+j``, estimated by the observer at time ``i=k`` or ``i=k-1`` depending on its `direct` flag. Note that ``\mathbf{X̂_0 = X̂}`` if the operating point is zero, which is typically the case in practice -for [`NonLinModel`](@ref). This transcription method is generally more efficient for large +for [`NonLinModel`](@ref). + +This transcription computes the predictions by calling the augmented discrete-time model +in the equality constraint function recursively over ``H_p``, or by updating the linear +equality constraint vector for [`LinModel`](@ref). It is generally more efficient for large control horizon ``H_c``, unstable or highly nonlinear plant models/constraints. Sparse optimizers like `OSQP` or `Ipopt` and sparse Jacobian computations are recommended @@ -61,34 +68,20 @@ struct MultipleShooting <: ShootingMethod end Construct an implicit trapezoidal [`TranscriptionMethod`](@ref). -This is the simplest collocation method. It is only supported for continuous-time -[`NonLinModel`](@ref)s. It can handle moderately stiff systems and is A-stable. The decision -variables are the same as for [`MultipleShooting`](@ref), hence similar computational costs. -However, it may not be as efficient as more advanced collocation methods for highly stiff -systems. It currently assumes piecewise constant manipulated inputs (or zero-order hold) -between the samples, but linear interpolation will be added in the future. See Extended Help -for more details on the collocation/defect constraints. - -# Extended Help +This is the simplest collocation method. It supports continuous-time [`NonLinModel`](@ref)s +only. The decision variables are the same as for [`MultipleShooting`](@ref), hence similar +computational costs. It currently assumes piecewise constant manipulated inputs (or zero- +order hold) between the samples, but linear interpolation will be added soon. -!!! details "Extended Help" - The implicit trapezoidal collocation estimates the defects with: - ```math - \mathbf{s_d}(k) = \mathbf{x_0} \frac{T_s}{2}\big(\mathbf{F̂}(k+1) + \mathbf{F̂}(k)\big) - ``` - where ``T_s`` is the sampling period, ``\mathbf{Ẽ}`` the matrix defined at - [`init_defectmat`](@ref), and ``\mathbf{F̂}(k+j)`` the stacked vector of system +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 +moderately stiff systems and is A-stable. However, it may not be as efficient as more +advanced collocation methods for highly stiff systems. Note that the stochastic model of the +unmeasured disturbances is strictly discrete-time, it is thus transcribed separately using +[`MultipleShooting`](@ref). - where ``\mathbf{f̂}(k+j) = \mathbf{f̂}\big(\mathbf{x̂}(k+j), \mathbf{u}(k+j), \mathbf{d}(k+j)\big)``. - This leads to the following defect constraints for ``j=0`` to ``H_p-1``: - ```math - \mathbf{ŝ}(k+j) = \mathbf{x̂}(k+j+1) - \mathbf{x̂}(k+j) - \frac{T_s}{2} \big( \mathbf{f̂}(k+j) + \mathbf{f̂}(k+j+1) \big) = 0 - ``` - which are added as equality constraints in the optimization problem. The initial state - ``\mathbf{x̂}(k)`` is given by the state estimator, and the future states - ``\mathbf{x̂}(k+j+1)`` are decision variables in the optimization problem. The method - requires evaluating the system dynamics at both the current and next time steps, which - can increase computational complexity compared to explicit methods like single shooting. +Sparse optimizers like `Ipopt` and sparse Jacobian computations are recommended for this +transcription method. """ struct TrapezoidalCollocation <: CollocationMethod end From 4137880365860bf82f403db905607d15cd561011 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Thu, 28 Aug 2025 16:23:42 -0400 Subject: [PATCH 16/19] added: adapt `K0` size to number of collocation points --- src/controller/nonlinmpc.jl | 4 +++- src/controller/transcription.jl | 27 +++++++++++++++++++-------- 2 files changed, 22 insertions(+), 9 deletions(-) diff --git a/src/controller/nonlinmpc.jl b/src/controller/nonlinmpc.jl index 995cbcd57..76dad045a 100644 --- a/src/controller/nonlinmpc.jl +++ b/src/controller/nonlinmpc.jl @@ -562,8 +562,10 @@ Inspired from: [User-defined operators with vector outputs](@extref JuMP User-de function get_optim_functions(mpc::NonLinMPC, ::JuMP.GenericModel{JNT}) where JNT<:Real # ----------- common cache for Jfunc, gfuncs and geqfuncs ---------------------------- model = mpc.estim.model + transcription = mpc.transcription grad, jac = mpc.gradient, mpc.jacobian - nu, ny, nx̂, nϵ, nk = model.nu, model.ny, mpc.estim.nx̂, mpc.nϵ, model.nk + nu, ny, nx̂, nϵ = model.nu, model.ny, mpc.estim.nx̂, mpc.nϵ + nk = get_nk(model, transcription) Hp, Hc = mpc.Hp, mpc.Hc ng, nc, neq = length(mpc.con.i_g), mpc.con.nc, mpc.con.neq nZ̃, nU, nŶ, nX̂, nK = length(mpc.Z̃), Hp*nu, Hp*ny, Hp*nx̂, Hp*nk diff --git a/src/controller/transcription.jl b/src/controller/transcription.jl index a6fe8290a..6504c4d23 100644 --- a/src/controller/transcription.jl +++ b/src/controller/transcription.jl @@ -78,12 +78,18 @@ equality constraint function and by using the implicit trapezoidal rule. It can moderately stiff systems and is A-stable. However, it may not be as efficient as more advanced collocation methods for highly stiff systems. Note that the stochastic model of the unmeasured disturbances is strictly discrete-time, it is thus transcribed separately using -[`MultipleShooting`](@ref). +[`MultipleShooting`](@ref). Also note that the state Sparse optimizers like `Ipopt` and sparse Jacobian computations are recommended for this transcription method. """ -struct TrapezoidalCollocation <: CollocationMethod end +struct TrapezoidalCollocation <: CollocationMethod + nc::Int + function TrapezoidalCollocation() + nc = 2 # 2 collocation points per interval for trapezoidal rule + return new(nc) + end +end function validate_transcription(::LinModel, ::CollocationMethod) throw(ArgumentError("Collocation methods are not supported for LinModel.")) @@ -103,6 +109,10 @@ function get_nZ(estim::StateEstimator, ::TranscriptionMethod, Hp, Hc) return estim.model.nu*Hc + estim.nx̂*Hp end +"Get length of the `k` vector with all the solver intermediate steps or all the collocation pts." +get_nk(model::SimModel, ::ShootingMethod) = model.nk +get_nk(model::SimModel, transcription::CollocationMethod) = model.nx*transcription.nc + @doc raw""" init_ZtoΔU(estim::StateEstimator, transcription::TranscriptionMethod, Hp, Hc) -> PΔu @@ -1271,7 +1281,8 @@ end """ con_nonlinprogeq!( geq, X̂0, Û0, K0 - mpc::PredictiveController, model::NonLinModel, ::MultipleShooting, U0, Z̃ + mpc::PredictiveController, model::NonLinModel, transcription::MultipleShooting, + U0, Z̃ ) Nonlinear equality constrains for [`NonLinModel`](@ref) and [`MultipleShooting`](@ref). @@ -1311,7 +1322,8 @@ end @doc raw""" con_nonlinprogeq!( geq, X̂0, Û0, K0 - mpc::PredictiveController, model::NonLinModel, ::TrapezoidalCollocation, U0, Z̃ + mpc::PredictiveController, model::NonLinModel, transcription::TrapezoidalCollocation, + U0, Z̃ ) Nonlinear equality constrains for [`NonLinModel`](@ref) and [`TrapezoidalCollocation`](@ref). @@ -1334,16 +1346,15 @@ hence no need to add `model.fop` and subtract `model.xop` in the collocation equ """ function con_nonlinprogeq!( geq, X̂0, Û0, K0, - mpc::PredictiveController, model::NonLinModel, ::TrapezoidalCollocation, U0, Z̃ + mpc::PredictiveController, model::NonLinModel, transcription::TrapezoidalCollocation, + U0, Z̃ ) nu, nx̂, nd, nx = model.nu, mpc.estim.nx̂, model.nd, model.nx Hp, Hc = mpc.Hp, mpc.Hc nΔU, nX̂ = nu*Hc, nx̂*Hp - Ts, p = model.Ts, model.p As, Cs_u = mpc.estim.As, mpc.estim.Cs_u - # K0 stores the state derivatives at each collocation point for all intervals - nk = model.nk # TODO: initialize K0 to the adequate size for a given collocation method + nk = get_nk(model, transcription) D̂0 = mpc.D̂0 X̂0_Z̃ = @views Z̃[(nΔU+1):(nΔU+nX̂)] x̂0 = @views mpc.estim.x̂0[1:nx̂] From 29ed6637e5ac4d27165685b6ac35914a80ecd510 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Thu, 28 Aug 2025 16:35:57 -0400 Subject: [PATCH 17/19] doc: mention that the builtin estimator uses `solver` --- src/controller/transcription.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/controller/transcription.jl b/src/controller/transcription.jl index 6504c4d23..cf5caf3a1 100644 --- a/src/controller/transcription.jl +++ b/src/controller/transcription.jl @@ -78,7 +78,9 @@ equality constraint function and by using the implicit trapezoidal rule. It can moderately stiff systems and is A-stable. However, it may not be as efficient as more advanced collocation methods for highly stiff systems. Note that the stochastic model of the unmeasured disturbances is strictly discrete-time, it is thus transcribed separately using -[`MultipleShooting`](@ref). Also note that the state +[`MultipleShooting`](@ref). Also note that the built-in [`StateEstimator`](@ref) will still +use the `solver` provided at the construction of the [`NonLinModel`](@ref) to estimate the +plant states, not the trapezoidal rule (see `supersample` option for stiff systems). Sparse optimizers like `Ipopt` and sparse Jacobian computations are recommended for this transcription method. From d503153e4ffc5cd69615e23c1105dc622c05b2a0 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Thu, 28 Aug 2025 17:17:14 -0400 Subject: [PATCH 18/19] doc: move details on the stochastic part in Extended Help --- src/controller/transcription.jl | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/src/controller/transcription.jl b/src/controller/transcription.jl index cf5caf3a1..8f7f602b8 100644 --- a/src/controller/transcription.jl +++ b/src/controller/transcription.jl @@ -76,14 +76,20 @@ order hold) between the samples, but linear interpolation will be added soon. 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 moderately stiff systems and is A-stable. However, it may not be as efficient as more -advanced collocation methods for highly stiff systems. Note that the stochastic model of the -unmeasured disturbances is strictly discrete-time, it is thus transcribed separately using -[`MultipleShooting`](@ref). Also note that the built-in [`StateEstimator`](@ref) will still -use the `solver` provided at the construction of the [`NonLinModel`](@ref) to estimate the -plant states, not the trapezoidal rule (see `supersample` option for stiff systems). +advanced collocation methods for highly stiff systems. Note that the built-in [`StateEstimator`](@ref) +will still use the `solver` provided at the construction of the [`NonLinModel`](@ref) to +estimate the plant states, not the trapezoidal rule (see `supersample` option of +[`RungeKutta`](@ref) for stiff systems). See Extended Help for more details. Sparse optimizers like `Ipopt` and sparse Jacobian computations are recommended for this transcription method. + +# Extended Help +!!! details "Extended Help" + Note that the stochastic model of the unmeasured disturbances is strictly discrete-time, + as described in [`ModelPredictiveControl.init_estimstoch`](@ref). Collocation methods + require continuous-time dynamics. Because of this, the stochastic states are transcribed + separately using a [`MultipleShooting`](@ref) method. """ struct TrapezoidalCollocation <: CollocationMethod nc::Int From bf181aca88a0d67d39993361cdd9ca79bae2bac6 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Thu, 28 Aug 2025 17:25:40 -0400 Subject: [PATCH 19/19] bench: cancel current job on new commits --- .github/workflows/benchmark.yml | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/.github/workflows/benchmark.yml b/.github/workflows/benchmark.yml index 06eb8cb30..3c3f1dc1d 100644 --- a/.github/workflows/benchmark.yml +++ b/.github/workflows/benchmark.yml @@ -2,6 +2,11 @@ name: Benchmark on: pull_request_target: branches: [ main ] +concurrency: + # Skip intermediate builds: always. + # Cancel intermediate builds: only if it is a pull request build. + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} permissions: pull-requests: write # needed to post comments jobs: