From a09b32fdb87c9420c298e3a1d32fa05786b5aea2 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Sun, 26 Oct 2025 18:17:48 -0400 Subject: [PATCH 01/19] wip: supporting exact hessian --- src/controller/nonlinmpc.jl | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/src/controller/nonlinmpc.jl b/src/controller/nonlinmpc.jl index 66473eebb..ac9930eca 100644 --- a/src/controller/nonlinmpc.jl +++ b/src/controller/nonlinmpc.jl @@ -68,7 +68,7 @@ struct NonLinMPC{ estim::SE, Hp, Hc, nb, weights::CW, JE::JEfunc, gc!::GCfunc, nc, p::PT, transcription::TM, optim::JM, - gradient::GB, jacobian::JB, oracle + gradient::GB, jacobian::JB, hessian::HB, oracle ) where { NT<:Real, SE<:StateEstimator, @@ -77,6 +77,7 @@ struct NonLinMPC{ JM<:JuMP.GenericModel, GB<:AbstractADType, JB<:AbstractADType, + HB<:Union{AbstractADType, Nothing}, PT<:Any, JEfunc<:Function, GCfunc<:Function, @@ -218,6 +219,8 @@ This controller allocates memory at each time step for the optimization. function, see [`DifferentiationInterface` doc](@extref DifferentiationInterface List). - `jacobian=default_jacobian(transcription)` : an `AbstractADType` backend for the Jacobian of the nonlinear constraints, see `gradient` above for the options (default in Extended Help). +- `hessian=false` : an `AbstractADType` backend for the Hessian of the Lagrangian, see + `gradient` above for the options (`false` to skip it and use `optim` approximation). - `oracle=JuMP.solver_name(optim)=="Ipopt"`: use the efficient [`VectorNonlinearOracle`](@extref MathOptInterface MathOptInterface.VectorNonlinearOracle) for the nonlinear constraints (not supported by most optimizers for now). - additional keyword arguments are passed to [`UnscentedKalmanFilter`](@ref) constructor @@ -314,6 +317,7 @@ function NonLinMPC( optim::JuMP.GenericModel = JuMP.Model(DEFAULT_NONLINMPC_OPTIMIZER, add_bridges=false), gradient::AbstractADType = DEFAULT_NONLINMPC_GRADIENT, jacobian::AbstractADType = default_jacobian(transcription), + hessian::Union{AbstractADType, Bool, Nothing} = false, oracle::Bool = JuMP.solver_name(optim)=="Ipopt", kwargs... ) @@ -321,7 +325,7 @@ function NonLinMPC( return NonLinMPC( estim; Hp, Hc, Mwt, Nwt, Lwt, Cwt, Ewt, JE, gc, nc, p, M_Hp, N_Hc, L_Hp, - transcription, optim, gradient, jacobian, oracle + transcription, optim, gradient, jacobian, hessian, oracle ) end @@ -379,6 +383,7 @@ function NonLinMPC( optim::JuMP.GenericModel = JuMP.Model(DEFAULT_NONLINMPC_OPTIMIZER, add_bridges=false), gradient::AbstractADType = DEFAULT_NONLINMPC_GRADIENT, jacobian::AbstractADType = default_jacobian(transcription), + hessian::Union{AbstractADType, Bool, Nothing} = false, oracle::Bool = JuMP.solver_name(optim)=="Ipopt" ) where { NT<:Real, @@ -394,15 +399,24 @@ function NonLinMPC( validate_JE(NT, JE) gc! = get_mutating_gc(NT, gc) weights = ControllerWeights(estim.model, Hp, Hc, M_Hp, N_Hc, L_Hp, Cwt, Ewt) + hessian = default_hessian(gradient, jacobian, hessian) return NonLinMPC{NT}( estim, Hp, Hc, nb, weights, JE, gc!, nc, p, - transcription, optim, gradient, jacobian, oracle + transcription, optim, gradient, jacobian, hessian, oracle ) end default_jacobian(::SingleShooting) = DEFAULT_NONLINMPC_JACDENSE default_jacobian(::TranscriptionMethod) = DEFAULT_NONLINMPC_JACSPARSE +function default_hessian(gradient, jacobian, hessian::Bool) + if hessian + return DEFAULT_NONLINMPC_HESSIAN + else + return nothing + end +end + """ validate_JE(NT, JE) -> nothing From 8592d487e683f8d6600b923d5de5a52205cd44c9 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Thu, 30 Oct 2025 17:29:22 -0400 Subject: [PATCH 02/19] wip: exact hessian --- src/ModelPredictiveControl.jl | 5 ++-- src/controller/nonlinmpc.jl | 49 +++++++++++++++++++++++++++-------- 2 files changed, 41 insertions(+), 13 deletions(-) diff --git a/src/ModelPredictiveControl.jl b/src/ModelPredictiveControl.jl index ec246c121..59184e050 100644 --- a/src/ModelPredictiveControl.jl +++ b/src/ModelPredictiveControl.jl @@ -7,8 +7,9 @@ using Random: randn using RecipesBase using DifferentiationInterface: ADTypes.AbstractADType, AutoForwardDiff, AutoSparse -using DifferentiationInterface: gradient!, jacobian!, prepare_gradient, prepare_jacobian -using DifferentiationInterface: value_and_gradient!, value_and_jacobian! +using DifferentiationInterface: gradient!, value_and_gradient!, prepare_gradient +using DifferentiationInterface: jacobian!, value_and_jacobian!, prepare_jacobian +using DifferentiationInterface: hessian!, value_gradient_and_hessian!, prepare_hessian using DifferentiationInterface: Constant, Cache using SparseConnectivityTracer: TracerSparsityDetector using SparseMatrixColorings: GreedyColoringAlgorithm, sparsity_pattern diff --git a/src/controller/nonlinmpc.jl b/src/controller/nonlinmpc.jl index ac9930eca..646fbcc1f 100644 --- a/src/controller/nonlinmpc.jl +++ b/src/controller/nonlinmpc.jl @@ -7,6 +7,7 @@ const DEFAULT_NONLINMPC_JACSPARSE = AutoSparse( sparsity_detector=TracerSparsityDetector(), coloring_algorithm=GreedyColoringAlgorithm(), ) +const DEFAULT_NONLINMPC_HESSIAN = DEFAULT_NONLINMPC_JACSPARSE struct NonLinMPC{ NT<:Real, @@ -15,7 +16,8 @@ struct NonLinMPC{ TM<:TranscriptionMethod, JM<:JuMP.GenericModel, GB<:AbstractADType, - JB<:AbstractADType, + JB<:AbstractADType, + HB<:Union{AbstractADType, Nothing}, PT<:Any, JEfunc<:Function, GCfunc<:Function @@ -28,6 +30,7 @@ struct NonLinMPC{ con::ControllerConstraint{NT, GCfunc} gradient::GB jacobian::JB + hessian::HB oracle::Bool Z̃::Vector{NT} ŷ::Vector{NT} @@ -117,9 +120,9 @@ struct NonLinMPC{ nZ̃ = get_nZ(estim, transcription, Hp, Hc) + nϵ Z̃ = zeros(NT, nZ̃) buffer = PredictiveControllerBuffer(estim, transcription, Hp, Hc, nϵ) - mpc = new{NT, SE, CW, TM, JM, GB, JB, PT, JEfunc, GCfunc}( + mpc = new{NT, SE, CW, TM, JM, GB, JB, HB, PT, JEfunc, GCfunc}( estim, transcription, optim, con, - gradient, jacobian, oracle, + gradient, jacobian, hessian, oracle, Z̃, ŷ, Hp, Hc, nϵ, nb, weights, @@ -619,7 +622,7 @@ function get_nonlinops(mpc::NonLinMPC, optim::JuMP.GenericModel{JNT}) where JNT< # ----------- common cache for all functions ---------------------------------------- model = mpc.estim.model transcription = mpc.transcription - grad, jac = mpc.gradient, mpc.jacobian + grad, jac, hess = mpc.gradient, mpc.jacobian, mpc.hessian nu, ny, nx̂, nϵ = model.nu, model.ny, mpc.estim.nx̂, mpc.nϵ nk = get_nk(model, transcription) Hp, Hc = mpc.Hp, mpc.Hc @@ -645,14 +648,28 @@ function get_nonlinops(mpc::NonLinMPC, optim::JuMP.GenericModel{JNT}) where JNT< gi .= @views g[i_g] return nothing end - Z̃_∇gi = fill(myNaN, nZ̃) # NaN to force update_predictions! at first call + function ℓ_gi(Z̃_λ_gi, ΔŨ, x̂0end, Ue, Ŷe, U0, Ŷ0, Û0, K0, X̂0, gc, geq, g, gi) + Z̃, λ = @views Z̃_λ_gi[begin:begin+nZ̃-1], Z̃_λ_gi[begin+nZ̃:end] + update_predictions!(ΔŨ, x̂0end, Ue, Ŷe, U0, Ŷ0, Û0, K0, X̂0, gc, g, geq, mpc, Z̃) + gi .= @views g[i_g] + return dot(λ, gi) + end + Z̃_∇gi = fill(myNaN, nZ̃) # NaN to force update at first call + Z̃_λ_gi = fill(myNaN, nZ̃ + ngi) ∇gi_context = ( Cache(ΔŨ), Cache(x̂0end), Cache(Ue), Cache(Ŷe), Cache(U0), Cache(Ŷ0), Cache(Û0), Cache(K0), Cache(X̂0), Cache(gc), Cache(geq), Cache(g) ) - ∇gi_prep = prepare_jacobian(gi!, gi, jac, Z̃_∇gi, ∇gi_context...; strict) - ∇gi = init_diffmat(JNT, jac, ∇gi_prep, nZ̃, ngi) + ∇gi_prep = prepare_jacobian(gi!, gi, jac, Z̃_∇gi, ∇gi_context...; strict) + ∇²gi_context = ( + Cache(ΔŨ), Cache(x̂0end), Cache(Ue), Cache(Ŷe), Cache(U0), Cache(Ŷ0), + Cache(Û0), Cache(K0), Cache(X̂0), + Cache(gc), Cache(geq), Cache(g), Cache(gi) + ) + ∇²gi_prep = prepare_hessian(ℓ_gi, hess, Z̃_λ_gi, ∇²gi_context...; strict) + ∇gi = init_diffmat(JNT, jac, ∇gi_prep, nZ̃, ngi) + ∇²ℓ_gi = init_diffmat(JNT, hess, ∇²gi_prep, nZ̃ + ngi, nZ̃ + ngi) function update_con!(gi, ∇gi, Z̃_∇gi, Z̃_arg) if isdifferent(Z̃_arg, Z̃_∇gi) Z̃_∇gi .= Z̃_arg @@ -668,23 +685,33 @@ function get_nonlinops(mpc::NonLinMPC, optim::JuMP.GenericModel{JNT}) where JNT< update_con!(gi, ∇gi, Z̃_∇gi, Z̃_arg) return diffmat2vec!(∇gi_arg, ∇gi) end + function ∇²gi_func!(∇²ℓ_arg, Z̃_arg, λ_arg) + Z̃_λ_gi[1:begin:begin+nZ̃-1] .= Z̃_arg + Z̃_λ_gi[[begin+nZ̃:end]] .= λ_arg + hessian!(ℓ_gi, ∇²ℓ_gi, ∇²gi_prep, hess, Z̃_λ_gi, ∇²gi_context) + return diffmat2vec!(∇²ℓ_arg, ∇²ℓ_gi) + end gi_min = fill(-myInf, ngi) gi_max = zeros(JNT, ngi) - ∇gi_structure = init_diffstructure(∇gi) + ∇gi_structure = init_diffstructure(∇gi) + ∇²gi_structure = init_diffstructure(∇²ℓ_gi) + display(∇²ℓ_gi) g_oracle = MOI.VectorNonlinearOracle(; dimension = nZ̃, l = gi_min, u = gi_max, eval_f = gi_func!, jacobian_structure = ∇gi_structure, - eval_jacobian = ∇gi_func! + eval_jacobian = ∇gi_func!, + hessian_lagrangian_structure = ∇²gi_structure, + eval_hessian_lagrangian = ∇²gi_func! ) # ------------- equality constraints : nonlinear oracle ------------------------------ function geq!(geq, Z̃, ΔŨ, x̂0end, Ue, Ŷe, U0, Ŷ0, Û0, K0, X̂0, gc, g) update_predictions!(ΔŨ, x̂0end, Ue, Ŷe, U0, Ŷ0, Û0, K0, X̂0, gc, g, geq, mpc, Z̃) return nothing end - Z̃_∇geq = fill(myNaN, nZ̃) # NaN to force update_predictions! at first call + Z̃_∇geq = fill(myNaN, nZ̃) # NaN to force update at first call ∇geq_context = ( Cache(ΔŨ), Cache(x̂0end), Cache(Ue), Cache(Ŷe), Cache(U0), Cache(Ŷ0), Cache(Û0), Cache(K0), Cache(X̂0), @@ -722,7 +749,7 @@ function get_nonlinops(mpc::NonLinMPC, optim::JuMP.GenericModel{JNT}) where JNT< update_predictions!(ΔŨ, x̂0end, Ue, Ŷe, U0, Ŷ0, Û0, K0, X̂0, gc, g, geq, mpc, Z̃) return obj_nonlinprog!(Ŷ0, U0, mpc, model, Ue, Ŷe, ΔŨ) end - Z̃_∇J = fill(myNaN, nZ̃) # NaN to force update_predictions! at first call + Z̃_∇J = fill(myNaN, nZ̃) # NaN to force update at first call ∇J_context = ( Cache(ΔŨ), Cache(x̂0end), Cache(Ue), Cache(Ŷe), Cache(U0), Cache(Ŷ0), Cache(Û0), Cache(K0), Cache(X̂0), From bbf48dd369cd6b995ebd60b59efb41159dc8bccc Mon Sep 17 00:00:00 2001 From: franckgaga Date: Thu, 30 Oct 2025 17:41:16 -0400 Subject: [PATCH 03/19] doc: some comments on hessian --- src/controller/nonlinmpc.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/controller/nonlinmpc.jl b/src/controller/nonlinmpc.jl index 646fbcc1f..1d0150042 100644 --- a/src/controller/nonlinmpc.jl +++ b/src/controller/nonlinmpc.jl @@ -706,6 +706,8 @@ function get_nonlinops(mpc::NonLinMPC, optim::JuMP.GenericModel{JNT}) where JNT< hessian_lagrangian_structure = ∇²gi_structure, eval_hessian_lagrangian = ∇²gi_func! ) + #TODO: verify if I must fill only upper/lower triangular part ? + #TODO: add Hessian for 1. Jfunc and 2. geq # ------------- equality constraints : nonlinear oracle ------------------------------ function geq!(geq, Z̃, ΔŨ, x̂0end, Ue, Ŷe, U0, Ŷ0, Û0, K0, X̂0, gc, g) update_predictions!(ΔŨ, x̂0end, Ue, Ŷe, U0, Ŷ0, Û0, K0, X̂0, gc, g, geq, mpc, Z̃) From fd52b787429a3c350e9782168d48bc8abd0ca7eb Mon Sep 17 00:00:00 2001 From: franckgaga Date: Fri, 31 Oct 2025 23:14:07 -0400 Subject: [PATCH 04/19] =?UTF-8?q?added:=20`hessian`=20now=20works=20well?= =?UTF-8?q?=20for=20`NonLinMPC`=20=F0=9F=8D=BE?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/controller/nonlinmpc.jl | 163 +++++++++++++++++++++++++----------- src/general.jl | 24 +++++- 2 files changed, 137 insertions(+), 50 deletions(-) diff --git a/src/controller/nonlinmpc.jl b/src/controller/nonlinmpc.jl index 1d0150042..a40f3abe9 100644 --- a/src/controller/nonlinmpc.jl +++ b/src/controller/nonlinmpc.jl @@ -642,34 +642,37 @@ function get_nonlinops(mpc::NonLinMPC, optim::JuMP.GenericModel{JNT}) where JNT< Û0::Vector{JNT}, X̂0::Vector{JNT} = zeros(JNT, nU), zeros(JNT, nX̂) gc::Vector{JNT}, g::Vector{JNT} = zeros(JNT, nc), zeros(JNT, ng) gi::Vector{JNT}, geq::Vector{JNT} = zeros(JNT, ngi), zeros(JNT, neq) + λi::Vector{JNT}, λeq::Vector{JNT} = zeros(JNT, ngi), zeros(JNT, neq) # -------------- inequality constraint: nonlinear oracle ----------------------------- function gi!(gi, Z̃, ΔŨ, x̂0end, Ue, Ŷe, U0, Ŷ0, Û0, K0, X̂0, gc, geq, g) update_predictions!(ΔŨ, x̂0end, Ue, Ŷe, U0, Ŷ0, Û0, K0, X̂0, gc, g, geq, mpc, Z̃) gi .= @views g[i_g] return nothing end - function ℓ_gi(Z̃_λ_gi, ΔŨ, x̂0end, Ue, Ŷe, U0, Ŷ0, Û0, K0, X̂0, gc, geq, g, gi) - Z̃, λ = @views Z̃_λ_gi[begin:begin+nZ̃-1], Z̃_λ_gi[begin+nZ̃:end] + function ℓ_gi(Z̃, λi, ΔŨ, x̂0end, Ue, Ŷe, U0, Ŷ0, Û0, K0, X̂0, gc, geq, g, gi) update_predictions!(ΔŨ, x̂0end, Ue, Ŷe, U0, Ŷ0, Û0, K0, X̂0, gc, g, geq, mpc, Z̃) gi .= @views g[i_g] - return dot(λ, gi) + return dot(λi, gi) end Z̃_∇gi = fill(myNaN, nZ̃) # NaN to force update at first call - Z̃_λ_gi = fill(myNaN, nZ̃ + ngi) ∇gi_context = ( Cache(ΔŨ), Cache(x̂0end), Cache(Ue), Cache(Ŷe), Cache(U0), Cache(Ŷ0), Cache(Û0), Cache(K0), Cache(X̂0), Cache(gc), Cache(geq), Cache(g) ) ∇gi_prep = prepare_jacobian(gi!, gi, jac, Z̃_∇gi, ∇gi_context...; strict) - ∇²gi_context = ( - Cache(ΔŨ), Cache(x̂0end), Cache(Ue), Cache(Ŷe), Cache(U0), Cache(Ŷ0), - Cache(Û0), Cache(K0), Cache(X̂0), - Cache(gc), Cache(geq), Cache(g), Cache(gi) - ) - ∇²gi_prep = prepare_hessian(ℓ_gi, hess, Z̃_λ_gi, ∇²gi_context...; strict) - ∇gi = init_diffmat(JNT, jac, ∇gi_prep, nZ̃, ngi) - ∇²ℓ_gi = init_diffmat(JNT, hess, ∇²gi_prep, nZ̃ + ngi, nZ̃ + ngi) + ∇gi = init_diffmat(JNT, jac, ∇gi_prep, nZ̃, ngi) + ∇gi_structure = init_diffstructure(∇gi) + if !isnothing(hess) + ∇²gi_context = ( + Cache(ΔŨ), Cache(x̂0end), Cache(Ue), Cache(Ŷe), Cache(U0), Cache(Ŷ0), + Cache(Û0), Cache(K0), Cache(X̂0), + Cache(gc), Cache(geq), Cache(g), Cache(gi) + ) + ∇²gi_prep = prepare_hessian(ℓ_gi, hess, Z̃_∇gi, Constant(λi), ∇²gi_context...; strict) + ∇²ℓ_gi = init_diffmat(JNT, hess, ∇²gi_prep, nZ̃, nZ̃) + ∇²gi_structure = lowertriangle_indices(init_diffstructure(∇²ℓ_gi)) + end function update_con!(gi, ∇gi, Z̃_∇gi, Z̃_arg) if isdifferent(Z̃_arg, Z̃_∇gi) Z̃_∇gi .= Z̃_arg @@ -682,20 +685,17 @@ function get_nonlinops(mpc::NonLinMPC, optim::JuMP.GenericModel{JNT}) where JNT< return gi_arg .= gi end function ∇gi_func!(∇gi_arg, Z̃_arg) - update_con!(gi, ∇gi, Z̃_∇gi, Z̃_arg) - return diffmat2vec!(∇gi_arg, ∇gi) + update_con!(gi, ∇gi, Z̃_∇gi, Z̃_arg) + return diffmat2vec!(∇gi_arg, ∇gi, ∇gi_structure) end function ∇²gi_func!(∇²ℓ_arg, Z̃_arg, λ_arg) - Z̃_λ_gi[1:begin:begin+nZ̃-1] .= Z̃_arg - Z̃_λ_gi[[begin+nZ̃:end]] .= λ_arg - hessian!(ℓ_gi, ∇²ℓ_gi, ∇²gi_prep, hess, Z̃_λ_gi, ∇²gi_context) - return diffmat2vec!(∇²ℓ_arg, ∇²ℓ_gi) + Z̃_∇gi .= Z̃_arg + λi .= λ_arg + hessian!(ℓ_gi, ∇²ℓ_gi, ∇²gi_prep, hess, Z̃_∇gi, Constant(λi), ∇²gi_context...) + return diffmat2vec!(∇²ℓ_arg, ∇²ℓ_gi, ∇²gi_structure) end gi_min = fill(-myInf, ngi) gi_max = zeros(JNT, ngi) - ∇gi_structure = init_diffstructure(∇gi) - ∇²gi_structure = init_diffstructure(∇²ℓ_gi) - display(∇²ℓ_gi) g_oracle = MOI.VectorNonlinearOracle(; dimension = nZ̃, l = gi_min, @@ -703,16 +703,18 @@ function get_nonlinops(mpc::NonLinMPC, optim::JuMP.GenericModel{JNT}) where JNT< eval_f = gi_func!, jacobian_structure = ∇gi_structure, eval_jacobian = ∇gi_func!, - hessian_lagrangian_structure = ∇²gi_structure, - eval_hessian_lagrangian = ∇²gi_func! + hessian_lagrangian_structure = isnothing(hess) ? Tuple{Int,Int}[] : ∇²gi_structure, + eval_hessian_lagrangian = isnothing(hess) ? nothing : ∇²gi_func! ) - #TODO: verify if I must fill only upper/lower triangular part ? - #TODO: add Hessian for 1. Jfunc and 2. geq # ------------- equality constraints : nonlinear oracle ------------------------------ function geq!(geq, Z̃, ΔŨ, x̂0end, Ue, Ŷe, U0, Ŷ0, Û0, K0, X̂0, gc, g) update_predictions!(ΔŨ, x̂0end, Ue, Ŷe, U0, Ŷ0, Û0, K0, X̂0, gc, g, geq, mpc, Z̃) return nothing end + function ℓ_geq(Z̃, λeq, ΔŨ, x̂0end, Ue, Ŷe, U0, Ŷ0, Û0, K0, X̂0, gc, geq, g) + update_predictions!(ΔŨ, x̂0end, Ue, Ŷe, U0, Ŷ0, Û0, K0, X̂0, gc, g, geq, mpc, Z̃) + return dot(λeq, geq) + end Z̃_∇geq = fill(myNaN, nZ̃) # NaN to force update at first call ∇geq_context = ( Cache(ΔŨ), Cache(x̂0end), Cache(Ue), Cache(Ŷe), Cache(U0), Cache(Ŷ0), @@ -720,7 +722,18 @@ function get_nonlinops(mpc::NonLinMPC, optim::JuMP.GenericModel{JNT}) where JNT< Cache(gc), Cache(g) ) ∇geq_prep = prepare_jacobian(geq!, geq, jac, Z̃_∇geq, ∇geq_context...; strict) - ∇geq = init_diffmat(JNT, jac, ∇geq_prep, nZ̃, neq) + ∇geq = init_diffmat(JNT, jac, ∇geq_prep, nZ̃, neq) + ∇geq_structure = init_diffstructure(∇geq) + if !isnothing(hess) + ∇²geq_context = ( + Cache(ΔŨ), Cache(x̂0end), Cache(Ue), Cache(Ŷe), Cache(U0), Cache(Ŷ0), + Cache(Û0), Cache(K0), Cache(X̂0), + Cache(gc), Cache(geq), Cache(g) + ) + ∇²geq_prep = prepare_hessian(ℓ_geq, hess, Z̃_∇geq, Constant(λeq), ∇²geq_context...; strict) + ∇²ℓ_geq = init_diffmat(JNT, hess, ∇²geq_prep, nZ̃, nZ̃) + ∇²geq_structure = lowertriangle_indices(init_diffstructure(∇²ℓ_geq)) + end function update_con_eq!(geq, ∇geq, Z̃_∇geq, Z̃_arg) if isdifferent(Z̃_arg, Z̃_∇geq) Z̃_∇geq .= Z̃_arg @@ -734,53 +747,109 @@ function get_nonlinops(mpc::NonLinMPC, optim::JuMP.GenericModel{JNT}) where JNT< end function ∇geq_func!(∇geq_arg, Z̃_arg) update_con_eq!(geq, ∇geq, Z̃_∇geq, Z̃_arg) - return diffmat2vec!(∇geq_arg, ∇geq) + return diffmat2vec!(∇geq_arg, ∇geq, ∇geq_structure) + end + function ∇²geq_func!(∇²ℓ_arg, Z̃_arg, λ_arg) + Z̃_∇geq .= Z̃_arg + λeq .= λ_arg + hessian!(ℓ_geq, ∇²ℓ_geq, ∇²geq_prep, hess, Z̃_∇geq, Constant(λeq), ∇²geq_context...) + return diffmat2vec!(∇²ℓ_arg, ∇²ℓ_geq, ∇²geq_structure) end geq_min = geq_max = zeros(JNT, neq) - ∇geq_structure = init_diffstructure(∇geq) geq_oracle = MOI.VectorNonlinearOracle(; dimension = nZ̃, l = geq_min, u = geq_max, eval_f = geq_func!, jacobian_structure = ∇geq_structure, - eval_jacobian = ∇geq_func! + eval_jacobian = ∇geq_func!, + hessian_lagrangian_structure = isnothing(hess) ? Tuple{Int,Int}[] : ∇²geq_structure, + eval_hessian_lagrangian = isnothing(hess) ? nothing : ∇²geq_func! ) # ------------- objective function: splatting syntax --------------------------------- function J!(Z̃, ΔŨ, x̂0end, Ue, Ŷe, U0, Ŷ0, Û0, K0, X̂0, gc, g, geq) update_predictions!(ΔŨ, x̂0end, Ue, Ŷe, U0, Ŷ0, Û0, K0, X̂0, gc, g, geq, mpc, Z̃) return obj_nonlinprog!(Ŷ0, U0, mpc, model, Ue, Ŷe, ΔŨ) end - Z̃_∇J = fill(myNaN, nZ̃) # NaN to force update at first call - ∇J_context = ( + Z̃_J = fill(myNaN, nZ̃) # NaN to force update at first call + J_context = ( Cache(ΔŨ), Cache(x̂0end), Cache(Ue), Cache(Ŷe), Cache(U0), Cache(Ŷ0), Cache(Û0), Cache(K0), Cache(X̂0), Cache(gc), Cache(g), Cache(geq), ) - ∇J_prep = prepare_gradient(J!, grad, Z̃_∇J, ∇J_context...; strict) - ∇J = Vector{JNT}(undef, nZ̃) - function update_objective!(J, ∇J, Z̃_∇J, Z̃_arg) - if isdifferent(Z̃_arg, Z̃_∇J) - Z̃_∇J .= Z̃_arg - J[], _ = value_and_gradient!(J!, ∇J, ∇J_prep, grad, Z̃_∇J, ∇J_context...) + ∇J_prep = prepare_gradient(J!, grad, Z̃_J, J_context...; strict) + ∇J = Vector{JNT}(undef, nZ̃) + if !isnothing(hess) + ∇²J_prep = prepare_hessian(J!, hess, Z̃_J, J_context...; strict) + ∇²J = init_diffmat(JNT, hess, ∇²J_prep, nZ̃, nZ̃) + end + update_objective! = if !isnothing(hess) + function (J, ∇J, ∇²J, Z̃_J, Z̃_arg) + if isdifferent(Z̃_arg, Z̃_J) + Z̃_J .= Z̃_arg + J[], _ = value_gradient_and_hessian!(J!, ∇J, ∇²J, hess, Z̃_J, J_context...) + end + end + else + update_objective! = function (J, ∇J, Z̃_∇J, Z̃_arg) + if isdifferent(Z̃_arg, Z̃_∇J) + Z̃_∇J .= Z̃_arg + J[], _ = value_and_gradient!(J!, ∇J, ∇J_prep, grad, Z̃_∇J, J_context...) + end + end + end + J_func = if !isnothing(hess) + function (Z̃_arg::Vararg{T, N}) where {N, T<:Real} + update_objective!(J, ∇J, ∇²J, Z̃_J, Z̃_arg) + return J[]::T + end + else + function (Z̃_arg::Vararg{T, N}) where {N, T<:Real} + update_objective!(J, ∇J, Z̃_J, Z̃_arg) + return J[]::T end - end - function J_func(Z̃_arg::Vararg{T, N}) where {N, T<:Real} - update_objective!(J, ∇J, Z̃_∇J, Z̃_arg) - return J[]::T end ∇J_func! = if nZ̃ == 1 # univariate syntax (see JuMP.@operator doc): + if !isnothing(hess) + function (Z̃_arg) + update_objective!(J, ∇J, ∇²J, Z̃_J, Z̃_arg) + return ∇J[] + end + else + function (Z̃_arg) + update_objective!(J, ∇J, Z̃_J, Z̃_arg) + return ∇J[] + end + end + else # multivariate syntax (see JuMP.@operator doc): + if !isnothing(hess) + function (∇J_arg::AbstractVector{T}, Z̃_arg::Vararg{T, N}) where {N, T<:Real} + update_objective!(J, ∇J, ∇²J, Z̃_J, Z̃_arg) + return ∇J_arg .= ∇J + end + else + function (∇J_arg::AbstractVector{T}, Z̃_arg::Vararg{T, N}) where {N, T<:Real} + update_objective!(J, ∇J, Z̃_J, Z̃_arg) + return ∇J_arg .= ∇J + end + end + end + ∇²J_func! = if nZ̃ == 1 # univariate syntax (see JuMP.@operator doc): function (Z̃_arg) - update_objective!(J, ∇J, Z̃_∇J, Z̃_arg) - return ∇J[] + update_objective!(J, ∇J, ∇²J, Z̃_J, Z̃_arg) + return ∇²J[] end else # multivariate syntax (see JuMP.@operator doc): - function (∇J_arg::AbstractVector{T}, Z̃_arg::Vararg{T, N}) where {N, T<:Real} - update_objective!(J, ∇J, Z̃_∇J, Z̃_arg) - return ∇J_arg .= ∇J + function (∇²J_arg::AbstractMatrix{T}, Z̃_arg::Vararg{T, N}) where {N, T<:Real} + update_objective!(J, ∇J, ∇²J, Z̃_J, Z̃_arg) + return fill_lowertriangle!(∇²J_arg, ∇²J) end end - J_op = JuMP.add_nonlinear_operator(optim, nZ̃, J_func, ∇J_func!, name=:J_op) + J_op = if !isnothing(hess) + JuMP.add_nonlinear_operator(optim, nZ̃, J_func, ∇J_func!, ∇²J_func!, name=:J_op) + else + JuMP.add_nonlinear_operator(optim, nZ̃, J_func, ∇J_func!, name=:J_op) + end return g_oracle, geq_oracle, J_op end diff --git a/src/general.jl b/src/general.jl index 7dec4709b..b5e1498ad 100644 --- a/src/general.jl +++ b/src/general.jl @@ -71,9 +71,27 @@ function init_diffstructure(A::AbstractSparseMatrix) end init_diffstructure(A::AbstractMatrix)= Tuple.(CartesianIndices(A))[:] -"Store the differentiation matrix `A` in the vector `v` as required by `JuMP.jl.`" -diffmat2vec!(v::AbstractVector, A::AbstractSparseMatrix) = v .= nonzeros(A) -diffmat2vec!(v::AbstractVector, A::AbstractMatrix) = v[:] = A +"Get the lower-triangular indices from the differentiation matrix structure." +function lowertriangle_indices(diffmat_struct::Vector{Tuple{Int, Int}}) + return [(i,j) for (i,j) in diffmat_struct if i ≥ j] +end + +"Fill the lower triangular part of A in-place with the corresponding part in B." +function fill_lowertriangle!(A::AbstractMatrix, B::AbstractMatrix) + for j in axes(A, 2), i in axes(A, 1) + (i ≥ j) && (A[i, j] = B[i, j]) + end + return A +end + +"Store the diff. matrix `A` in the vector `v` with list of nonzero indices `i_vec`" +function diffmat2vec!(v::AbstractVector, A::AbstractMatrix, i_vec::Vector{Tuple{Int, Int}}) + for i in eachindex(v) + i_A, j_A = i_vec[i] + v[i] = A[i_A, j_A] + end + return v +end backend_str(backend::AbstractADType) = string(nameof(typeof(backend))) function backend_str(backend::AutoSparse) From aefc4ffff399b1adf53616b5dc93296d3a500717 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Sat, 1 Nov 2025 00:50:19 -0400 Subject: [PATCH 05/19] added: validate hessian backend --- src/ModelPredictiveControl.jl | 3 ++- src/controller/nonlinmpc.jl | 51 ++++++++++++++++++++++++++--------- src/general.jl | 16 ++++++++--- 3 files changed, 52 insertions(+), 18 deletions(-) diff --git a/src/ModelPredictiveControl.jl b/src/ModelPredictiveControl.jl index 59184e050..423475c5f 100644 --- a/src/ModelPredictiveControl.jl +++ b/src/ModelPredictiveControl.jl @@ -6,7 +6,8 @@ using Random: randn using RecipesBase -using DifferentiationInterface: ADTypes.AbstractADType, AutoForwardDiff, AutoSparse +using DifferentiationInterface: ADTypes.AbstractADType, AutoForwardDiff +using DifferentiationInterface: AutoSparse, SecondOrder using DifferentiationInterface: gradient!, value_and_gradient!, prepare_gradient using DifferentiationInterface: jacobian!, value_and_jacobian!, prepare_jacobian using DifferentiationInterface: hessian!, value_gradient_and_hessian!, prepare_hessian diff --git a/src/controller/nonlinmpc.jl b/src/controller/nonlinmpc.jl index a40f3abe9..b26c8dfdb 100644 --- a/src/controller/nonlinmpc.jl +++ b/src/controller/nonlinmpc.jl @@ -223,7 +223,8 @@ This controller allocates memory at each time step for the optimization. - `jacobian=default_jacobian(transcription)` : an `AbstractADType` backend for the Jacobian of the nonlinear constraints, see `gradient` above for the options (default in Extended Help). - `hessian=false` : an `AbstractADType` backend for the Hessian of the Lagrangian, see - `gradient` above for the options (`false` to skip it and use `optim` approximation). + `gradient` above for the options. The default `false` skip it and use the quasi-Newton + method of `optim`, which is always the case if `oracle=false` (see Extended Help). - `oracle=JuMP.solver_name(optim)=="Ipopt"`: use the efficient [`VectorNonlinearOracle`](@extref MathOptInterface MathOptInterface.VectorNonlinearOracle) for the nonlinear constraints (not supported by most optimizers for now). - additional keyword arguments are passed to [`UnscentedKalmanFilter`](@ref) constructor @@ -282,9 +283,9 @@ NonLinMPC controller with a sample time Ts = 10.0 s: the `gc` argument (both `gc` and `gc!` accepts non-mutating and mutating functions). By default, the optimization relies on dense [`ForwardDiff`](@extref ForwardDiff) - automatic differentiation (AD) to compute the objective and constraint derivatives. One - exception: if `transcription` is not a [`SingleShooting`](@ref), the `jacobian` argument - defaults to this [sparse backend](@extref DifferentiationInterface AutoSparse-object): + automatic differentiation (AD) to compute the objective and constraint derivatives. Two + exceptions: if `transcription` is not a [`SingleShooting`](@ref), the `jacobian` + argument defaults to this [sparse backend](@extref DifferentiationInterface AutoSparse-object): ```julia AutoSparse( AutoForwardDiff(); @@ -292,6 +293,10 @@ NonLinMPC controller with a sample time Ts = 10.0 s: coloring_algorithm = GreedyColoringAlgorithm() ) ``` + This is also the sparse backend selected for the Hessian of the Lagrangian function if + `oracle=true` and `hessian=true`, which is the second exception. Second order + derivatives are only supported with `oracle=true` option. + Optimizers generally benefit from exact derivatives like AD. However, the [`NonLinModel`](@ref) state-space functions must be compatible with this feature. See [`JuMP` documentation](@extref JuMP Common-mistakes-when-writing-a-user-defined-operator) for common mistakes when writing these functions. @@ -402,7 +407,7 @@ function NonLinMPC( validate_JE(NT, JE) gc! = get_mutating_gc(NT, gc) weights = ControllerWeights(estim.model, Hp, Hc, M_Hp, N_Hc, L_Hp, Cwt, Ewt) - hessian = default_hessian(gradient, jacobian, hessian) + hessian = validate_hessian(hessian, gradient, oracle) return NonLinMPC{NT}( estim, Hp, Hc, nb, weights, JE, gc!, nc, p, transcription, optim, gradient, jacobian, hessian, oracle @@ -412,14 +417,6 @@ end default_jacobian(::SingleShooting) = DEFAULT_NONLINMPC_JACDENSE default_jacobian(::TranscriptionMethod) = DEFAULT_NONLINMPC_JACSPARSE -function default_hessian(gradient, jacobian, hessian::Bool) - if hessian - return DEFAULT_NONLINMPC_HESSIAN - else - return nothing - end -end - """ validate_JE(NT, JE) -> nothing @@ -513,6 +510,34 @@ function test_custom_functions(NT, model::SimModel, JE, gc!, nc, Uop, Yop, Dop, return nothing end +""" + validate_hessian(hessian, gradient, oracle) -> backend + +Validate `hessian` argument and return the differentiation backend. +""" +function validate_hessian(hessian, gradient, oracle) + if hessian == true + backend = DEFAULT_NONLINMPC_HESSIAN + elseif hessian == false || isnothing(hessian) + backend = nothing + else + backend = hessian + end + if oracle == false && !isnothing(backend) + error("Second order derivatives are only supported with oracle=true.") + end + if oracle == true && !isnothing(backend) + hess = dense_backend(backend) + grad = dense_backend(gradient) + if hess != grad + @info "The objective function gradient will be computed with the hessian "* + "backend ($(backend_str(hess)))\n instead of the one in gradient "* + "argument ($(backend_str(grad))) for efficiency." + end + end + return backend +end + """ addinfo!(info, mpc::NonLinMPC) -> info diff --git a/src/general.jl b/src/general.jl index b5e1498ad..6804bdcea 100644 --- a/src/general.jl +++ b/src/general.jl @@ -69,11 +69,11 @@ function init_diffstructure(A::AbstractSparseMatrix) I, J = findnz(A) return collect(zip(I, J)) end -init_diffstructure(A::AbstractMatrix)= Tuple.(CartesianIndices(A))[:] +init_diffstructure(A::AbstractMatrix) = Tuple.(CartesianIndices(A))[:] -"Get the lower-triangular indices from the differentiation matrix structure." -function lowertriangle_indices(diffmat_struct::Vector{Tuple{Int, Int}}) - return [(i,j) for (i,j) in diffmat_struct if i ≥ j] +"Get the lower-triangular indices from the differentiation matrix structure `i_vec`." +function lowertriangle_indices(i_vec::Vector{Tuple{Int, Int}}) + return [(i,j) for (i,j) in i_vec if i ≥ j] end "Fill the lower triangular part of A in-place with the corresponding part in B." @@ -100,6 +100,14 @@ function backend_str(backend::AutoSparse) " $(nameof(typeof(backend.coloring_algorithm))))" return str end +function backend_str(backend::SecondOrder) + str = "SecondOrder ($(nameof(typeof(backend.outer))),"* + " $(nameof(typeof(backend.inner))))" + return str +end +dense_backend(backend::AbstractADType) = backend +dense_backend(backend::AutoSparse) = backend.dense_ad +dense_backend(backend::SecondOrder) = backend.inner "Verify that x and y elements are different using `!==`." isdifferent(x, y) = any(xi !== yi for (xi, yi) in zip(x, y)) From bffa05ad934fd7dddb1c41b2ddbc2d0397a77391 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Sat, 1 Nov 2025 11:32:22 -0400 Subject: [PATCH 06/19] changed: two separate functions for objective and constraints --- src/controller/nonlinmpc.jl | 92 ++++++++++++++++++++++++++----------- 1 file changed, 64 insertions(+), 28 deletions(-) diff --git a/src/controller/nonlinmpc.jl b/src/controller/nonlinmpc.jl index b26c8dfdb..d9a69289d 100644 --- a/src/controller/nonlinmpc.jl +++ b/src/controller/nonlinmpc.jl @@ -589,8 +589,9 @@ function init_optimization!( end end if mpc.oracle - g_oracle, geq_oracle, J_op = get_nonlinops(mpc, optim) - optim[:J_op] = J_op + J_op = get_nonlinobj_op(mpc, optim) + g_oracle, geq_oracle = get_nonlincon_oracle(mpc, optim) + else J_func, ∇J_func!, g_funcs, ∇g_funcs!, geq_funcs, ∇geq_funcs! = get_optim_functions( mpc, optim @@ -616,7 +617,7 @@ Re-construct nonlinear constraints and add them to `mpc.optim`. """ function reset_nonlincon!(mpc::NonLinMPC) if mpc.oracle - g_oracle, geq_oracle = get_nonlinops(mpc, mpc.optim) + g_oracle, geq_oracle = get_nonlincon_oracle(mpc, mpc.optim) set_nonlincon!(mpc, mpc.optim, g_oracle, geq_oracle) else set_nonlincon_leg!(mpc, mpc.estim.model, mpc.transcription, mpc.optim) @@ -624,30 +625,22 @@ function reset_nonlincon!(mpc::NonLinMPC) end """ - get_nonlinops(mpc::NonLinMPC, optim) -> g_oracle, geq_oracle, J_op + get_nonlincon_oracle(mpc::NonLinMPC, optim) -> g_oracle, geq_oracle -Return the operators for the nonlinear optimization of `mpc` [`NonLinMPC`](@ref) controller. +Return the nonlinear constraint oracles for [`NonLinMPC`](@ref) `mpc`. Return `g_oracle` and `geq_oracle`, the inequality and equality [`VectorNonlinearOracle`](@extref MathOptInterface MathOptInterface.VectorNonlinearOracle) for the two respective constraints. Note that `g_oracle` only includes the non-`Inf` -inequality constraints, thus it must be re-constructed if they change. Also return `J_op`, -the [`NonlinearOperator`](@extref JuMP NonlinearOperator) for the objective function, based -on the splatting syntax. This method is really intricate and that's because of 3 elements: - -- These functions are used inside the nonlinear optimization, so they must be type-stable - and as efficient as possible. All the function outputs and derivatives are cached and - updated in-place if required to use the efficient [`value_and_jacobian!`](@extref DifferentiationInterface DifferentiationInterface.value_and_jacobian!). -- The splatting syntax for objective functions implies the use of `Vararg{T,N}` (see the [performance tip](@extref Julia Be-aware-of-when-Julia-avoids-specializing)) - and memoization to avoid redundant computations. This is already complex, but it's even - worse knowing that the automatic differentiation tools do not support splatting. -- The signature of gradient and hessian functions is not the same for univariate (`nZ̃ == 1`) - and multivariate (`nZ̃ > 1`) operators in `JuMP`. Both must be defined. +inequality constraints, thus it must be re-constructed if they change. This method is really +intricate because the oracles are used inside the nonlinear optimization, so they must be +type-stable and as efficient as possible. All the function outputs and derivatives are +ached and updated in-place if required to use the efficient [`value_and_jacobian!`](@extref DifferentiationInterface DifferentiationInterface.value_and_jacobian!). """ -function get_nonlinops(mpc::NonLinMPC, optim::JuMP.GenericModel{JNT}) where JNT<:Real +function get_nonlincon_oracle(mpc::NonLinMPC, ::JuMP.GenericModel{JNT}) where JNT<:Real # ----------- common cache for all functions ---------------------------------------- model = mpc.estim.model transcription = mpc.transcription - grad, jac, hess = mpc.gradient, mpc.jacobian, mpc.hessian + jac, hess = mpc.jacobian, mpc.hessian nu, ny, nx̂, nϵ = model.nu, model.ny, mpc.estim.nx̂, mpc.nϵ nk = get_nk(model, transcription) Hp, Hc = mpc.Hp, mpc.Hc @@ -658,7 +651,6 @@ function get_nonlinops(mpc::NonLinMPC, optim::JuMP.GenericModel{JNT}) where JNT< nΔŨ, nUe, nŶe = nu*Hc + nϵ, nU + nu, nŶ + ny strict = Val(true) myNaN, myInf = convert(JNT, NaN), convert(JNT, Inf) - J::Vector{JNT} = zeros(JNT, 1) ΔŨ::Vector{JNT} = zeros(JNT, nΔŨ) x̂0end::Vector{JNT} = zeros(JNT, nx̂) K0::Vector{JNT} = zeros(JNT, nK) @@ -667,7 +659,7 @@ function get_nonlinops(mpc::NonLinMPC, optim::JuMP.GenericModel{JNT}) where JNT< Û0::Vector{JNT}, X̂0::Vector{JNT} = zeros(JNT, nU), zeros(JNT, nX̂) gc::Vector{JNT}, g::Vector{JNT} = zeros(JNT, nc), zeros(JNT, ng) gi::Vector{JNT}, geq::Vector{JNT} = zeros(JNT, ngi), zeros(JNT, neq) - λi::Vector{JNT}, λeq::Vector{JNT} = zeros(JNT, ngi), zeros(JNT, neq) + λi::Vector{JNT}, λeq::Vector{JNT} = ones(JNT, ngi), ones(JNT, neq) # -------------- inequality constraint: nonlinear oracle ----------------------------- function gi!(gi, Z̃, ΔŨ, x̂0end, Ue, Ŷe, U0, Ŷ0, Û0, K0, X̂0, gc, geq, g) update_predictions!(ΔŨ, x̂0end, Ue, Ŷe, U0, Ŷ0, Û0, K0, X̂0, gc, g, geq, mpc, Z̃) @@ -694,7 +686,9 @@ function get_nonlinops(mpc::NonLinMPC, optim::JuMP.GenericModel{JNT}) where JNT< Cache(Û0), Cache(K0), Cache(X̂0), Cache(gc), Cache(geq), Cache(g), Cache(gi) ) - ∇²gi_prep = prepare_hessian(ℓ_gi, hess, Z̃_∇gi, Constant(λi), ∇²gi_context...; strict) + ∇²gi_prep = prepare_hessian( + ℓ_gi, hess, Z̃_∇gi, Constant(λi), ∇²gi_context...; strict + ) ∇²ℓ_gi = init_diffmat(JNT, hess, ∇²gi_prep, nZ̃, nZ̃) ∇²gi_structure = lowertriangle_indices(init_diffstructure(∇²ℓ_gi)) end @@ -755,7 +749,9 @@ function get_nonlinops(mpc::NonLinMPC, optim::JuMP.GenericModel{JNT}) where JNT< Cache(Û0), Cache(K0), Cache(X̂0), Cache(gc), Cache(geq), Cache(g) ) - ∇²geq_prep = prepare_hessian(ℓ_geq, hess, Z̃_∇geq, Constant(λeq), ∇²geq_context...; strict) + ∇²geq_prep = prepare_hessian( + ℓ_geq, hess, Z̃_∇geq, Constant(λeq), ∇²geq_context...; strict + ) ∇²ℓ_geq = init_diffmat(JNT, hess, ∇²geq_prep, nZ̃, nZ̃) ∇²geq_structure = lowertriangle_indices(init_diffstructure(∇²ℓ_geq)) end @@ -791,7 +787,47 @@ function get_nonlinops(mpc::NonLinMPC, optim::JuMP.GenericModel{JNT}) where JNT< hessian_lagrangian_structure = isnothing(hess) ? Tuple{Int,Int}[] : ∇²geq_structure, eval_hessian_lagrangian = isnothing(hess) ? nothing : ∇²geq_func! ) - # ------------- objective function: splatting syntax --------------------------------- + return g_oracle, geq_oracle +end + +""" + get_nonlinobj_op(mpc::NonLinMPC, optim::JuMP.GenericModel{JNT}) -> J_op + +Return the nonlinear operator for the objective function of `mpc` [`NonLinMPC`](@ref). + +It is based on the splatting syntax. This method is really intricate and that's because of: + +- These functions are used inside the nonlinear optimization, so they must be type-stable + and as efficient as possible. All the function outputs and derivatives are cached and + updated in-place if required to use the efficient [`value_and_gradient!`](@extref DifferentiationInterface DifferentiationInterface.value_and_jacobian!). +- The splatting syntax for objective functions implies the use of `Vararg{T,N}` (see the [performance tip](@extref Julia Be-aware-of-when-Julia-avoids-specializing)) + and memoization to avoid redundant computations. This is already complex, but it's even + worse knowing that the automatic differentiation tools do not support splatting. +- The signature of gradient and hessian functions is not the same for univariate (`nZ̃ == 1`) + and multivariate (`nZ̃ > 1`) operators in `JuMP`. Both must be defined. +""" +function get_nonlinobj_op(mpc::NonLinMPC, optim::JuMP.GenericModel{JNT}) where JNT<:Real + model = mpc.estim.model + transcription = mpc.transcription + grad, hess = mpc.gradient, mpc.hessian + 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 = length(mpc.con.i_g) + nc, neq = mpc.con.nc, mpc.con.neq + nZ̃, nU, nŶ, nX̂, nK = length(mpc.Z̃), Hp*nu, Hp*ny, Hp*nx̂, Hp*nk + nΔŨ, nUe, nŶe = nu*Hc + nϵ, nU + nu, nŶ + ny + strict = Val(true) + myNaN = convert(JNT, NaN) + J::Vector{JNT} = zeros(JNT, 1) + ΔŨ::Vector{JNT} = zeros(JNT, nΔŨ) + x̂0end::Vector{JNT} = zeros(JNT, nx̂) + K0::Vector{JNT} = zeros(JNT, nK) + Ue::Vector{JNT}, Ŷe::Vector{JNT} = zeros(JNT, nUe), zeros(JNT, nŶe) + U0::Vector{JNT}, Ŷ0::Vector{JNT} = zeros(JNT, nU), zeros(JNT, nŶ) + Û0::Vector{JNT}, X̂0::Vector{JNT} = zeros(JNT, nU), zeros(JNT, nX̂) + gc::Vector{JNT}, g::Vector{JNT} = zeros(JNT, nc), zeros(JNT, ng) + geq::Vector{JNT} = zeros(JNT, neq) function J!(Z̃, ΔŨ, x̂0end, Ue, Ŷe, U0, Ŷ0, Û0, K0, X̂0, gc, g, geq) update_predictions!(ΔŨ, x̂0end, Ue, Ŷe, U0, Ŷ0, Û0, K0, X̂0, gc, g, geq, mpc, Z̃) return obj_nonlinprog!(Ŷ0, U0, mpc, model, Ue, Ŷe, ΔŨ) @@ -870,12 +906,12 @@ function get_nonlinops(mpc::NonLinMPC, optim::JuMP.GenericModel{JNT}) where JNT< return fill_lowertriangle!(∇²J_arg, ∇²J) end end - J_op = if !isnothing(hess) - JuMP.add_nonlinear_operator(optim, nZ̃, J_func, ∇J_func!, ∇²J_func!, name=:J_op) + if !isnothing(hess) + @operator(optim, J_op, nZ̃, J_func, ∇J_func!, ∇²J_func!) else - JuMP.add_nonlinear_operator(optim, nZ̃, J_func, ∇J_func!, name=:J_op) + @operator(optim, J_op, nZ̃, J_func, ∇J_func!) end - return g_oracle, geq_oracle, J_op + return J_op end """ From 37d4e31da4e52165aadb9846d4c0d208ba68f7cf Mon Sep 17 00:00:00 2001 From: franckgaga Date: Sat, 1 Nov 2025 11:58:55 -0400 Subject: [PATCH 07/19] debug: new `diffmat2vec!` in mHE --- docs/src/internals/predictive_control.md | 3 ++- src/estimator/mhe/construct.jl | 4 ++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/docs/src/internals/predictive_control.md b/docs/src/internals/predictive_control.md index b57b4c0e9..4e8e74d17 100644 --- a/docs/src/internals/predictive_control.md +++ b/docs/src/internals/predictive_control.md @@ -24,7 +24,8 @@ ModelPredictiveControl.relaxterminal ModelPredictiveControl.init_quadprog ModelPredictiveControl.init_stochpred ModelPredictiveControl.init_matconstraint_mpc -ModelPredictiveControl.get_nonlinops(::NonLinMPC, ::ModelPredictiveControl.GenericModel) +ModelPredictiveControl.get_nonlinobj_op(::NonLinMPC, ::ModelPredictiveControl.GenericModel) +ModelPredictiveControl.get_nonlincon_oracle(::NonLinMPC, ::ModelPredictiveControl.GenericModel) ``` ## Update Quadratic Optimization diff --git a/src/estimator/mhe/construct.jl b/src/estimator/mhe/construct.jl index e8f63e4b5..40474f426 100644 --- a/src/estimator/mhe/construct.jl +++ b/src/estimator/mhe/construct.jl @@ -1372,6 +1372,7 @@ function get_nonlinops( ∇gi_prep = prepare_jacobian(gi!, gi, jac, Z̃_∇gi, ∇gi_context...; strict) estim.Nk[] = 0 ∇gi = init_diffmat(JNT, jac, ∇gi_prep, nZ̃, ngi) + ∇gi_structure = init_diffstructure(∇gi) function update_con!(gi, ∇gi, Z̃_∇gi, Z̃_arg) if isdifferent(Z̃_arg, Z̃_∇gi) Z̃_∇gi .= Z̃_arg @@ -1385,11 +1386,10 @@ function get_nonlinops( end function ∇gi_func!(∇gi_vec, Z̃_arg) update_con!(gi, ∇gi, Z̃_∇gi, Z̃_arg) - return diffmat2vec!(∇gi_vec, ∇gi) + return diffmat2vec!(∇gi_vec, ∇gi, ∇gi_structure) end gi_min = fill(-myInf, ngi) gi_max = zeros(JNT, ngi) - ∇gi_structure = init_diffstructure(∇gi) g_oracle = MOI.VectorNonlinearOracle(; dimension = nZ̃, l = gi_min, From 71c924600aff7ef78ad46a599bf950745b18c026 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Sat, 1 Nov 2025 13:20:21 -0400 Subject: [PATCH 08/19] test: new simple tests with `hessian=true` --- test/3_test_predictive_control.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/test/3_test_predictive_control.jl b/test/3_test_predictive_control.jl index a7f9edb28..e75dd002b 100644 --- a/test/3_test_predictive_control.jl +++ b/test/3_test_predictive_control.jl @@ -724,10 +724,12 @@ end @test size(nmpc17.con.Aeq, 1) == nmpc17.estim.nx̂*nmpc17.Hp nmpc18 = NonLinMPC(nonlinmodel, Hp=10, gradient=AutoFiniteDiff(), - jacobian=AutoFiniteDiff() + jacobian=AutoFiniteDiff(), + hessian=AutoFiniteDiff() ) @test nmpc18.gradient == AutoFiniteDiff() @test nmpc18.jacobian == AutoFiniteDiff() + @test nmpc18.hessian == AutoFiniteDiff() nonlinmodel2 = NonLinModel{Float32}(f, h, Ts, 2, 4, 2, 1, solver=nothing) nmpc15 = NonLinMPC(nonlinmodel2, Hp=15) @@ -742,6 +744,7 @@ end @test_throws ErrorException NonLinMPC(nonlinmodel, Hp=15, gc! = (_,_,_,_)->[0.0], nc=1) @test_throws ArgumentError NonLinMPC(nonlinmodel, transcription=TrapezoidalCollocation()) @test_throws ArgumentError NonLinMPC(nonlinmodel, transcription=TrapezoidalCollocation(2)) + @test_throws ErrorException NonLinMPC(linmodel1, oracle=false, hessian=AutoFiniteDiff()) @test_logs (:warn, Regex(".*")) NonLinMPC(nonlinmodel, Hp=15, JE=(Ue,_,_,_)->Ue) @test_logs (:warn, Regex(".*")) NonLinMPC(nonlinmodel, Hp=15, gc=(Ue,_,_,_,_)->Ue, nc=0) @@ -856,7 +859,8 @@ end nmpc10 = setconstraint!(NonLinMPC( nonlinmodel, Nwt=[0], Hp=100, Hc=1, gradient=AutoFiniteDiff(), - jacobian=AutoFiniteDiff()), + jacobian=AutoFiniteDiff(), + hessian=true), ymax=[100], ymin=[-100] ) preparestate!(nmpc10, [0], [0]) From 456eb6c9734755145f4045133f6e8056a10cfc10 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Sat, 1 Nov 2025 14:43:40 -0400 Subject: [PATCH 09/19] test: improve coverage with `SecondOrder` backends --- test/3_test_predictive_control.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/test/3_test_predictive_control.jl b/test/3_test_predictive_control.jl index e75dd002b..f4280d477 100644 --- a/test/3_test_predictive_control.jl +++ b/test/3_test_predictive_control.jl @@ -730,6 +730,11 @@ end @test nmpc18.gradient == AutoFiniteDiff() @test nmpc18.jacobian == AutoFiniteDiff() @test nmpc18.hessian == AutoFiniteDiff() + nmpc19 = NonLinMPC(nonlinmodel, Hc=1, Hp=10, Cwt=Inf, + hessian=SecondOrder(AutoForwardDiff(), AutoForardDiff()) + ) + @test nmpc19.hessian == SecondOrder(AutoForwardDiff(), AutoForardDiff()) + @test_nowarn repr(nmpc19) # printing SecondOrder backends, for coverage nonlinmodel2 = NonLinModel{Float32}(f, h, Ts, 2, 4, 2, 1, solver=nothing) nmpc15 = NonLinMPC(nonlinmodel2, Hp=15) From d97d30cdb039f36417f93fb47b2098e790da635c Mon Sep 17 00:00:00 2001 From: franckgaga Date: Sat, 1 Nov 2025 15:17:15 -0400 Subject: [PATCH 10/19] test: debug test --- test/3_test_predictive_control.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/3_test_predictive_control.jl b/test/3_test_predictive_control.jl index f4280d477..311ded449 100644 --- a/test/3_test_predictive_control.jl +++ b/test/3_test_predictive_control.jl @@ -731,7 +731,7 @@ end @test nmpc18.jacobian == AutoFiniteDiff() @test nmpc18.hessian == AutoFiniteDiff() nmpc19 = NonLinMPC(nonlinmodel, Hc=1, Hp=10, Cwt=Inf, - hessian=SecondOrder(AutoForwardDiff(), AutoForardDiff()) + hessian=SecondOrder(AutoForwardDiff(), AutoForwardDiff()) ) @test nmpc19.hessian == SecondOrder(AutoForwardDiff(), AutoForardDiff()) @test_nowarn repr(nmpc19) # printing SecondOrder backends, for coverage From 38436a58d1c95df3ef265a5c47f84f4c5cdc5e37 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Sat, 1 Nov 2025 15:46:50 -0400 Subject: [PATCH 11/19] test : idem --- test/3_test_predictive_control.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/3_test_predictive_control.jl b/test/3_test_predictive_control.jl index 311ded449..e230e16fb 100644 --- a/test/3_test_predictive_control.jl +++ b/test/3_test_predictive_control.jl @@ -733,7 +733,7 @@ end nmpc19 = NonLinMPC(nonlinmodel, Hc=1, Hp=10, Cwt=Inf, hessian=SecondOrder(AutoForwardDiff(), AutoForwardDiff()) ) - @test nmpc19.hessian == SecondOrder(AutoForwardDiff(), AutoForardDiff()) + @test nmpc19.hessian == SecondOrder(AutoForwardDiff(), AutoForwardDiff()) @test_nowarn repr(nmpc19) # printing SecondOrder backends, for coverage nonlinmodel2 = NonLinModel{Float32}(f, h, Ts, 2, 4, 2, 1, solver=nothing) From 35b7f03ee55096969a323eab7f41b59f44941a50 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Sat, 1 Nov 2025 16:30:00 -0400 Subject: [PATCH 12/19] bench: new `NonLinMPC` benchs with Hessian --- benchmark/3_bench_predictive_control.jl | 42 +++++++++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/benchmark/3_bench_predictive_control.jl b/benchmark/3_bench_predictive_control.jl index 8a624a654..20a82203a 100644 --- a/benchmark/3_bench_predictive_control.jl +++ b/benchmark/3_bench_predictive_control.jl @@ -46,10 +46,18 @@ nmpc_nonlin_ss = NonLinMPC( nonlinmodel, transcription=SingleShooting(), Mwt=[1, 1], Nwt=[0.1, 0.1], Lwt=[0.1, 0.1], Hp=10 ) +nmpc_nonlin_ss_hess = NonLinMPC( + nonlinmodel_c, transcription=SingleShooting(), hessian=true, + Mwt=[1], Nwt=[0.1], Lwt=[0.1], Hp=10 +) nmpc_nonlin_ms = NonLinMPC( nonlinmodel, transcription=MultipleShooting(), Mwt=[1, 1], Nwt=[0.1, 0.1], Lwt=[0.1, 0.1], Hp=10 ) +nmpc_nonlin_ms_hess = NonLinMPC( + nonlinmodel_c, transcription=MultipleShooting(), hessian=true, + Mwt=[1], Nwt=[0.1], Lwt=[0.1], Hp=10 +) nmpc_nonlin_tc = NonLinMPC( nonlinmodel_c, transcription=TrapezoidalCollocation(), Mwt=[1], Nwt=[0.1], Lwt=[0.1], Hp=10 @@ -74,12 +82,24 @@ UNIT_MPC["NonLinMPC"]["moveinput!"]["NonLinModel"]["SingleShooting"] = setup=preparestate!($nmpc_nonlin_ss, $y, $d), samples=samples, evals=evals, seconds=seconds ) +UNIT_MPC["NonLinMPC"]["moveinput!"]["NonLinModel"]["SingleShootingHessian"] = + @benchmarkable( + moveinput!($nmpc_nonlin_ss, $y, $d), + setup=preparestate!($nmpc_nonlin_ss_hess, $y, $d), + samples=samples, evals=evals, seconds=seconds + ) UNIT_MPC["NonLinMPC"]["moveinput!"]["NonLinModel"]["MultipleShooting"] = @benchmarkable( moveinput!($nmpc_nonlin_ms, $y, $d), setup=preparestate!($nmpc_nonlin_ms, $y, $d), samples=samples, evals=evals, seconds=seconds ) +UNIT_MPC["NonLinMPC"]["moveinput!"]["NonLinModel"]["MultipleShootingHessian"] = + @benchmarkable( + moveinput!($nmpc_nonlin_ms, $y, $d), + setup=preparestate!($nmpc_nonlin_ms_hess, $y, $d), + samples=samples, evals=evals, seconds=seconds + ) UNIT_MPC["NonLinMPC"]["moveinput!"]["NonLinModel"]["TrapezoidalCollocation"] = @benchmarkable( moveinput!($nmpc_nonlin_tc, $y_c, $d_c), @@ -258,12 +278,24 @@ nmpc_ipopt_ss = NonLinMPC(estim; Hp, Hc, Mwt, Nwt, Cwt, optim, transcription) nmpc_ipopt_ss = setconstraint!(nmpc_ipopt_ss; umin, umax) JuMP.unset_time_limit_sec(nmpc_ipopt_ss.optim) +optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false) +transcription, hessian = SingleShooting(), true +nmpc_ipopt_ss_hess = NonLinMPC(estim; Hp, Hc, Mwt, Nwt, Cwt, optim, transcription, hessian) +nmpc_ipopt_ss_hess = setconstraint!(nmpc_ipopt_ss_hess; umin, umax) +JuMP.unset_time_limit_sec(nmpc_ipopt_ss_hess.optim) + optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false) transcription = MultipleShooting() nmpc_ipopt_ms = NonLinMPC(estim; Hp, Hc, Mwt, Nwt, Cwt, optim, transcription) nmpc_ipopt_ms = setconstraint!(nmpc_ipopt_ms; umin, umax) JuMP.unset_time_limit_sec(nmpc_ipopt_ms.optim) +optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false) +transcription, hessian = MultipleShooting(), true +nmpc_ipopt_ms_hess = NonLinMPC(estim; Hp, Hc, Mwt, Nwt, Cwt, optim, transcription, hessian) +nmpc_ipopt_ms_hess = setconstraint!(nmpc_ipopt_ms_hess; umin, umax) +JuMP.unset_time_limit_sec(nmpc_nmpc_ipopt_ms_hess.optim) + optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false) transcription = MultipleShooting(f_threads=true) nmpc_ipopt_mst = NonLinMPC(estim; Hp, Hc, Mwt, Nwt, Cwt, optim, transcription) @@ -308,11 +340,21 @@ CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["Ipopt"]["SingleShooting"] = sim!($nmpc_ipopt_ss, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0), samples=samples, evals=evals, seconds=seconds ) +CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["Ipopt"]["SingleShooting (Hessian)"] = + @benchmarkable( + sim!($nmpc_ipopt_ss_hess, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0), + samples=samples, evals=evals, seconds=seconds + ) CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["Ipopt"]["MultipleShooting"] = @benchmarkable( sim!($nmpc_ipopt_ms, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0), samples=samples, evals=evals, seconds=seconds ) +CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["Ipopt"]["MultipleShooting (Hessian)"] = + @benchmarkable( + sim!($nmpc_ipopt_ms_hess, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0), + samples=samples, evals=evals, seconds=seconds + ) CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["Ipopt"]["MultipleShooting (threaded)"] = @benchmarkable( sim!($nmpc_ipopt_mst, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0), From 8ea2a06c5c6bfdb166c9aaedc1c13350e30f0a7b Mon Sep 17 00:00:00 2001 From: franckgaga Date: Sat, 1 Nov 2025 16:32:42 -0400 Subject: [PATCH 13/19] bench: debug bench --- benchmark/3_bench_predictive_control.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmark/3_bench_predictive_control.jl b/benchmark/3_bench_predictive_control.jl index 20a82203a..d833f383d 100644 --- a/benchmark/3_bench_predictive_control.jl +++ b/benchmark/3_bench_predictive_control.jl @@ -294,7 +294,7 @@ optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_b transcription, hessian = MultipleShooting(), true nmpc_ipopt_ms_hess = NonLinMPC(estim; Hp, Hc, Mwt, Nwt, Cwt, optim, transcription, hessian) nmpc_ipopt_ms_hess = setconstraint!(nmpc_ipopt_ms_hess; umin, umax) -JuMP.unset_time_limit_sec(nmpc_nmpc_ipopt_ms_hess.optim) +JuMP.unset_time_limit_sec(nmpc_ipopt_ms_hess.optim) optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false) transcription = MultipleShooting(f_threads=true) From c576842e8033207c4d8a0e5878effa0eaa2960fc Mon Sep 17 00:00:00 2001 From: franckgaga Date: Sat, 1 Nov 2025 16:54:38 -0400 Subject: [PATCH 14/19] added: pretty-print `hessian` backend for `NonLinMPC` --- src/controller/nonlinmpc.jl | 1 + src/general.jl | 1 + 2 files changed, 2 insertions(+) diff --git a/src/controller/nonlinmpc.jl b/src/controller/nonlinmpc.jl index d9a69289d..6c91fa4bc 100644 --- a/src/controller/nonlinmpc.jl +++ b/src/controller/nonlinmpc.jl @@ -981,4 +981,5 @@ end function print_backends(io::IO, mpc::NonLinMPC) println(io, "├ gradient: $(backend_str(mpc.gradient))") println(io, "├ jacobian: $(backend_str(mpc.jacobian))") + println(io, "├ hessian: $(backend_str(mpc.hessian))") end diff --git a/src/general.jl b/src/general.jl index 6804bdcea..b608bb415 100644 --- a/src/general.jl +++ b/src/general.jl @@ -94,6 +94,7 @@ function diffmat2vec!(v::AbstractVector, A::AbstractMatrix, i_vec::Vector{Tuple{ end backend_str(backend::AbstractADType) = string(nameof(typeof(backend))) +backend_str(backend::Nothing) = "nothing" function backend_str(backend::AutoSparse) str = "AutoSparse ($(nameof(typeof(backend.dense_ad))),"* " $(nameof(typeof(backend.sparsity_detector))),"* From 415b992e27b6e89bbaba102dd6d0bd62b00310ee Mon Sep 17 00:00:00 2001 From: franckgaga Date: Sat, 1 Nov 2025 17:04:18 -0400 Subject: [PATCH 15/19] bench: add `progress=false` to `sim!` calls --- benchmark/2_bench_state_estim.jl | 8 +++---- benchmark/3_bench_predictive_control.jl | 30 ++++++++++++------------- 2 files changed, 19 insertions(+), 19 deletions(-) diff --git a/benchmark/2_bench_state_estim.jl b/benchmark/2_bench_state_estim.jl index a05afe808..25dc87042 100644 --- a/benchmark/2_bench_state_estim.jl +++ b/benchmark/2_bench_state_estim.jl @@ -335,21 +335,21 @@ JuMP.set_attribute(mhe_pendulum_madnlp_pred.optim, "tol", 1e-7) samples, evals, seconds = 25, 1, 15*60 CASE_ESTIM["Pendulum"]["MovingHorizonEstimator"]["Ipopt"]["Current form"] = @benchmarkable( - sim!($mhe_pendulum_ipopt_curr, $N, $u; plant=$plant, x_0=$x_0, x̂_0=$x̂_0), + sim!($mhe_pendulum_ipopt_curr, $N, $u; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false), samples=samples, evals=evals, seconds=seconds ) CASE_ESTIM["Pendulum"]["MovingHorizonEstimator"]["Ipopt"]["Prediction form"] = @benchmarkable( - sim!($mhe_pendulum_ipopt_pred, $N, $u; plant=$plant, x_0=$x_0, x̂_0=$x̂_0), + sim!($mhe_pendulum_ipopt_pred, $N, $u; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false), samples=samples, evals=evals, seconds=seconds ) CASE_ESTIM["Pendulum"]["MovingHorizonEstimator"]["MadNLP"]["Current form"] = @benchmarkable( - sim!($mhe_pendulum_madnlp_curr, $N, $u; plant=$plant, x_0=$x_0, x̂_0=$x̂_0), + sim!($mhe_pendulum_madnlp_curr, $N, $u; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false), samples=samples, evals=evals, seconds=seconds ) CASE_ESTIM["Pendulum"]["MovingHorizonEstimator"]["MadNLP"]["Prediction form"] = @benchmarkable( - sim!($mhe_pendulum_madnlp_pred, $N, $u; plant=$plant, x_0=$x_0, x̂_0=$x̂_0), + sim!($mhe_pendulum_madnlp_pred, $N, $u; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false), samples=samples, evals=evals, seconds=seconds ) \ No newline at end of file diff --git a/benchmark/3_bench_predictive_control.jl b/benchmark/3_bench_predictive_control.jl index d833f383d..e157e8c49 100644 --- a/benchmark/3_bench_predictive_control.jl +++ b/benchmark/3_bench_predictive_control.jl @@ -337,42 +337,42 @@ JuMP.unset_time_limit_sec(nmpc_madnlp_ss.optim) samples, evals, seconds = 100, 1, 15*60 CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["Ipopt"]["SingleShooting"] = @benchmarkable( - sim!($nmpc_ipopt_ss, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0), + sim!($nmpc_ipopt_ss, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false), samples=samples, evals=evals, seconds=seconds ) CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["Ipopt"]["SingleShooting (Hessian)"] = @benchmarkable( - sim!($nmpc_ipopt_ss_hess, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0), + sim!($nmpc_ipopt_ss_hess, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false), samples=samples, evals=evals, seconds=seconds ) CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["Ipopt"]["MultipleShooting"] = @benchmarkable( - sim!($nmpc_ipopt_ms, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0), + sim!($nmpc_ipopt_ms, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false), samples=samples, evals=evals, seconds=seconds ) CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["Ipopt"]["MultipleShooting (Hessian)"] = @benchmarkable( - sim!($nmpc_ipopt_ms_hess, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0), + sim!($nmpc_ipopt_ms_hess, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false), samples=samples, evals=evals, seconds=seconds ) CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["Ipopt"]["MultipleShooting (threaded)"] = @benchmarkable( - sim!($nmpc_ipopt_mst, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0), + sim!($nmpc_ipopt_mst, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false), samples=samples, evals=evals, seconds=seconds ) CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["Ipopt"]["TrapezoidalCollocation"] = @benchmarkable( - sim!($nmpc_ipopt_tc, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0), + sim!($nmpc_ipopt_tc, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false), samples=samples, evals=evals, seconds=seconds ) CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["Ipopt"]["TrapezoidalCollocation (threaded)"] = @benchmarkable( - sim!($nmpc_ipopt_tct, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0), + sim!($nmpc_ipopt_tct, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false), samples=samples, evals=evals, seconds=seconds ) CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["MadNLP"]["SingleShooting"] = @benchmarkable( - sim!($nmpc_madnlp_ss, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0), + sim!($nmpc_madnlp_ss, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false), samples=samples, evals=evals, seconds=seconds ) @@ -418,22 +418,22 @@ JuMP.unset_time_limit_sec(empc_madnlp_ss.optim) samples, evals, seconds = 100, 1, 15*60 CASE_MPC["Pendulum"]["NonLinMPC"]["Economic"]["Ipopt"]["SingleShooting"] = @benchmarkable( - sim!($empc_ipopt_ss, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0), + sim!($empc_ipopt_ss, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0, progress=false), samples=samples, evals=evals, seconds=seconds ) CASE_MPC["Pendulum"]["NonLinMPC"]["Economic"]["Ipopt"]["MultipleShooting"] = @benchmarkable( - sim!($empc_ipopt_ms, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0), + sim!($empc_ipopt_ms, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0, progress=false), samples=samples, evals=evals, seconds=seconds ) CASE_MPC["Pendulum"]["NonLinMPC"]["Economic"]["Ipopt"]["TrapezoidalCollocation"] = @benchmarkable( - sim!($empc_ipopt_tc, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0), + sim!($empc_ipopt_tc, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0, progress=false), samples=samples, evals=evals, seconds=seconds ) CASE_MPC["Pendulum"]["NonLinMPC"]["Economic"]["MadNLP"]["SingleShooting"] = @benchmarkable( - sim!($empc_madnlp_ss, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0), + sim!($empc_madnlp_ss, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0, progress=false), samples=samples, evals=evals, seconds=seconds ) @@ -484,17 +484,17 @@ JuMP.unset_time_limit_sec(nmpc2_ipopt_tc.optim) samples, evals, seconds = 100, 1, 15*60 CASE_MPC["Pendulum"]["NonLinMPC"]["Custom constraints"]["Ipopt"]["SingleShooting"] = @benchmarkable( - sim!($nmpc2_ipopt_ss, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0), + sim!($nmpc2_ipopt_ss, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0), progress=false, samples=samples, evals=evals, seconds=seconds ) CASE_MPC["Pendulum"]["NonLinMPC"]["Custom constraints"]["Ipopt"]["MultipleShooting"] = @benchmarkable( - sim!($nmpc2_ipopt_ms, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0), + sim!($nmpc2_ipopt_ms, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0, progress=false), samples=samples, evals=evals, seconds=seconds ) CASE_MPC["Pendulum"]["NonLinMPC"]["Custom constraints"]["Ipopt"]["TrapezoidalCollocation"] = @benchmarkable( - sim!($nmpc2_ipopt_tc, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0), + sim!($nmpc2_ipopt_tc, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0, progress=false), samples=samples, evals=evals, seconds=seconds ) From 6a0639b20333bd292bb1461934aa591eaa9d3c46 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Sat, 1 Nov 2025 17:14:20 -0400 Subject: [PATCH 16/19] bench: debug benchmarks --- benchmark/3_bench_predictive_control.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/benchmark/3_bench_predictive_control.jl b/benchmark/3_bench_predictive_control.jl index e157e8c49..8bb3b0410 100644 --- a/benchmark/3_bench_predictive_control.jl +++ b/benchmark/3_bench_predictive_control.jl @@ -47,16 +47,16 @@ nmpc_nonlin_ss = NonLinMPC( Mwt=[1, 1], Nwt=[0.1, 0.1], Lwt=[0.1, 0.1], Hp=10 ) nmpc_nonlin_ss_hess = NonLinMPC( - nonlinmodel_c, transcription=SingleShooting(), hessian=true, - Mwt=[1], Nwt=[0.1], Lwt=[0.1], Hp=10 + nonlinmodel, transcription=SingleShooting(), hessian=true, + Mwt=[1, 1], Nwt=[0.1, 0.1], Lwt=[0.1, 0.1], Hp=10 ) nmpc_nonlin_ms = NonLinMPC( nonlinmodel, transcription=MultipleShooting(), Mwt=[1, 1], Nwt=[0.1, 0.1], Lwt=[0.1, 0.1], Hp=10 ) nmpc_nonlin_ms_hess = NonLinMPC( - nonlinmodel_c, transcription=MultipleShooting(), hessian=true, - Mwt=[1], Nwt=[0.1], Lwt=[0.1], Hp=10 + nonlinmodel, transcription=MultipleShooting(), hessian=true, + Mwt=[1, 1], Nwt=[0.1, 0.1], Lwt=[0.1, 0.1], Hp=10 ) nmpc_nonlin_tc = NonLinMPC( nonlinmodel_c, transcription=TrapezoidalCollocation(), From b332b83170788aeaf56b86c6096af354358b8ac5 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Sat, 1 Nov 2025 17:19:30 -0400 Subject: [PATCH 17/19] test: update `jldoctests` --- src/controller/nonlinmpc.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/controller/nonlinmpc.jl b/src/controller/nonlinmpc.jl index 6c91fa4bc..cabafde4c 100644 --- a/src/controller/nonlinmpc.jl +++ b/src/controller/nonlinmpc.jl @@ -242,6 +242,7 @@ NonLinMPC controller with a sample time Ts = 10.0 s: ├ transcription: MultipleShooting ├ gradient: AutoForwardDiff ├ jacobian: AutoSparse (AutoForwardDiff, TracerSparsityDetector, GreedyColoringAlgorithm) +├ hessian: nothing └ dimensions: ├ 20 prediction steps Hp ├ 10 control steps Hc @@ -355,10 +356,11 @@ julia> mpc = NonLinMPC(estim, Hp=20, Cwt=1e6) NonLinMPC controller with a sample time Ts = 10.0 s: ├ estimator: UnscentedKalmanFilter ├ model: NonLinModel -├ optimizer: Ipopt +├ optimizer: Ipopt ├ transcription: SingleShooting ├ gradient: AutoForwardDiff ├ jacobian: AutoForwardDiff +├ hessian: nothing └ dimensions: ├ 20 prediction steps Hp ├ 2 control steps Hc From eb92527b38a32028729198141d804377384b159c Mon Sep 17 00:00:00 2001 From: franckgaga Date: Sat, 1 Nov 2025 17:27:51 -0400 Subject: [PATCH 18/19] bench: debug bench --- benchmark/3_bench_predictive_control.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmark/3_bench_predictive_control.jl b/benchmark/3_bench_predictive_control.jl index 8bb3b0410..427f2e7e6 100644 --- a/benchmark/3_bench_predictive_control.jl +++ b/benchmark/3_bench_predictive_control.jl @@ -484,7 +484,7 @@ JuMP.unset_time_limit_sec(nmpc2_ipopt_tc.optim) samples, evals, seconds = 100, 1, 15*60 CASE_MPC["Pendulum"]["NonLinMPC"]["Custom constraints"]["Ipopt"]["SingleShooting"] = @benchmarkable( - sim!($nmpc2_ipopt_ss, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0), progress=false, + sim!($nmpc2_ipopt_ss, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0, progress=false), samples=samples, evals=evals, seconds=seconds ) CASE_MPC["Pendulum"]["NonLinMPC"]["Custom constraints"]["Ipopt"]["MultipleShooting"] = From df1bb376eb8fc85360c3807fe56e9d53dec0ff3b Mon Sep 17 00:00:00 2001 From: franckgaga Date: Sat, 1 Nov 2025 18:41:32 -0400 Subject: [PATCH 19/19] test: improve coverage --- test/3_test_predictive_control.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/3_test_predictive_control.jl b/test/3_test_predictive_control.jl index e230e16fb..db8ed3721 100644 --- a/test/3_test_predictive_control.jl +++ b/test/3_test_predictive_control.jl @@ -730,7 +730,8 @@ end @test nmpc18.gradient == AutoFiniteDiff() @test nmpc18.jacobian == AutoFiniteDiff() @test nmpc18.hessian == AutoFiniteDiff() - nmpc19 = NonLinMPC(nonlinmodel, Hc=1, Hp=10, Cwt=Inf, + nonlinmodel_simple = NonLinModel((x,u,_,_)->x+u,(x,_,_)->x, 1, 1, 1, 1, solver=nothing) + nmpc19 = NonLinMPC(nonlinmodel_simple, Hc=1, Hp=10, Cwt=Inf, hessian=SecondOrder(AutoForwardDiff(), AutoForwardDiff()) ) @test nmpc19.hessian == SecondOrder(AutoForwardDiff(), AutoForwardDiff())