Merge branch 'main' into extension

Author: franckgaga

.github/workflows/CI.yml

@@ -34,7 +34,7 @@ jobs:
       - name: Set JULIA_DEBUG environment variable if applicable
         if: ${{ runner.debug == '1' }}
         run: echo "JULIA_DEBUG=ModelPredictiveControl" >> $GITHUB_ENV
-      - uses: actions/checkout@v5
+      - uses: actions/checkout@v6
       - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}

.github/workflows/DocCleanup.yml

@@ -9,7 +9,7 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - name: Checkout gh-pages branch
-        uses: actions/checkout@v5
+        uses: actions/checkout@v6
         with:
           ref: gh-pages
 

.github/workflows/documentation.yml

@@ -13,7 +13,7 @@ jobs:
       contents: write
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v5
+      - uses: actions/checkout@v6
       - uses: julia-actions/setup-julia@v2
         with:
           version: '1'

Project.toml

@@ -1,6 +1,6 @@
 name = "ModelPredictiveControl"
 uuid = "61f9bdb8-6ae4-484a-811f-bbf86720c31c"
-version = "1.13.3"
+version = "1.14.1"
 authors = ["Francis Gagnon"]
 
 [deps]

benchmark/2_bench_state_estim.jl

@@ -172,6 +172,12 @@ UNIT_ESTIM["MovingHorizonEstimator"]["updatestate!"]["LinModel"]["Prediction for
         setup=preparestate!($mhe_lin_pred, $y, $d),
         samples=samples, evals=evals, seconds=seconds,
     )
+UNIT_ESTIM["MovingHorizonEstimator"]["getinfo!"]["LinModel"] =
+    @benchmarkable(
+        getinfo($mhe_lin_curr),
+        setup=preparestate!($mhe_lin_curr, $y, $d),
+        samples=samples, evals=evals, seconds=seconds,
+    )
 UNIT_ESTIM["MovingHorizonEstimator"]["preparestate!"]["NonLinModel"]["Current form"] =
     @benchmarkable(
         preparestate!($mhe_nonlin_curr, $y, $d),
@@ -193,7 +199,13 @@ UNIT_ESTIM["MovingHorizonEstimator"]["updatestate!"]["NonLinModel"]["Prediction
         updatestate!($mhe_nonlin_pred, $u, $y, $d),
         setup=preparestate!($mhe_nonlin_pred, $y, $d),
         samples=samples, evals=evals, seconds=seconds,
-    ) 
+    )
+UNIT_ESTIM["MovingHorizonEstimator"]["getinfo!"]["NonLinModel"] =
+    @benchmarkable(
+        getinfo($mhe_nonlin_curr),
+        setup=preparestate!($mhe_nonlin_curr, $y, $d),
+        samples=samples, evals=evals, seconds=seconds,
+    )
 
 ## ----------------------------------------------------------------------------------------
 ## ----------------- CASE STUDIES ---------------------------------------------------------

benchmark/3_bench_predictive_control.jl

@@ -25,6 +25,12 @@ UNIT_MPC["LinMPC"]["moveinput!"]["MultipleShooting"] =
         setup=preparestate!($linmpc_ms, $y, $d),
         samples=samples, evals=evals, seconds=seconds
     )
+UNIT_MPC["LinMPC"]["getinfo!"] =
+    @benchmarkable(
+        getinfo($linmpc_ss),
+        setup=(preparestate!($linmpc_ss, $y, $d); moveinput!($linmpc_ss, $y, $d)),
+        samples=samples, evals=evals, seconds=seconds
+    )
 
 empc = ExplicitMPC(linmodel, Mwt=[1, 1], Nwt=[0.1, 0.1], Lwt=[0.1, 0.1], Hp=10)
 
@@ -76,6 +82,12 @@ UNIT_MPC["NonLinMPC"]["moveinput!"]["LinModel"]["MultipleShooting"] =
         setup=preparestate!($nmpc_lin_ms, $y, $d),
         samples=samples, evals=evals, seconds=seconds
     )
+UNIT_MPC["NonLinMPC"]["getinfo!"]["LinModel"] =
+    @benchmarkable(
+        getinfo($nmpc_lin_ss),
+        setup=(preparestate!($nmpc_lin_ss, $y, $d); moveinput!($nmpc_lin_ss, $y, $d)),
+        samples=samples, evals=evals, seconds=seconds
+    )
 UNIT_MPC["NonLinMPC"]["moveinput!"]["NonLinModel"]["SingleShooting"] =
     @benchmarkable(
         moveinput!($nmpc_nonlin_ss, $y, $d),
@@ -106,6 +118,12 @@ UNIT_MPC["NonLinMPC"]["moveinput!"]["NonLinModel"]["TrapezoidalCollocation"] =
         setup=preparestate!($nmpc_nonlin_tc, $y_c, $d_c),
         samples=samples, evals=evals, seconds=seconds
     )
+UNIT_MPC["NonLinMPC"]["getinfo!"]["NonLinModel"] =
+    @benchmarkable(
+        getinfo($nmpc_nonlin_ss),
+        setup=(preparestate!($nmpc_nonlin_ss, $y, $d); moveinput!($nmpc_nonlin_ss, $y, $d)),
+        samples=samples, evals=evals, seconds=seconds
+    )
 
 ## ----------------------------------------------------------------------------------------
 ## ---------------------- CASE STUDIES ----------------------------------------------------

src/ModelPredictiveControl.jl

@@ -9,6 +9,9 @@ using RecipesBase
 
 using DifferentiationInterface: ADTypes.AbstractADType, AutoForwardDiff
 using DifferentiationInterface: AutoSparse, SecondOrder
+using DifferentiationInterface: gradient, jacobian, hessian
+using DifferentiationInterface: value_and_gradient, value_and_jacobian
+using DifferentiationInterface: value_gradient_and_hessian
 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

src/controller/execute.jl

@@ -105,8 +105,22 @@ The function should be called after calling [`moveinput!`](@ref). It returns the
 - `:d`   : current measured disturbance, ``\mathbf{d}(k)``
 
 For [`LinMPC`](@ref) and [`NonLinMPC`](@ref), the field `:sol` also contains the optimizer
-solution summary that can be printed. Lastly, the economical cost `:JE` and the custom
-nonlinear constraints `:gc` values at the optimum are also available for [`NonLinMPC`](@ref).
+solution summary that can be printed. Lastly, for [`NonLinMPC`](@ref), the following fields
+are also available:
+
+- `:JE`: economic cost value at the optimum, ``J_E``
+- `:gc`: custom nonlinear constraints values at the optimum, ``\mathbf{g_c}``
+- `:∇J` or *`:nablaJ`* : optimal gradient of the objective function, ``\mathbf{\nabla} J``
+- `:∇²J` or *`:nabla2J`* : optimal Hessian of the objective function, ``\mathbf{\nabla^2}J``
+- `:g` : optimal nonlinear inequality constraint values, ``\mathbf{g}``
+- `:∇g` or *`:nablag`* : optimal Jacobian of the inequality constraint, ``\mathbf{\nabla g}``
+- `:∇²ℓg` or *`:nabla2lg`* : optimal Hessian of the inequality Lagrangian, ``\mathbf{\nabla^2}\ell_{\mathbf{g}}``
+- `:geq` : optimal nonlinear equality constraint values, ``\mathbf{g_{eq}}``
+- `:∇geq` or *`:nablageq`* : optimal Jacobian of the equality constraint, ``\mathbf{\nabla g_{eq}}``
+- `:∇²ℓgeq` or *`:nabla2lgeq`* : optimal Hessian of the equality Lagrangian, ``\mathbf{\nabla^2}\ell_{\mathbf{g_{eq}}}``
+
+Note that retrieving optimal Hessians of Lagrangian are not fully supported yet. Their
+nonzero coefficients are random values for now.
 
 # Examples
 ```jldoctest

src/controller/nonlinmpc.jl

@@ -526,9 +526,11 @@ end
 """
     addinfo!(info, mpc::NonLinMPC) -> info
 
-For [`NonLinMPC`](@ref), add `:sol` and the optimal economic cost `:JE`.
+For [`NonLinMPC`](@ref), add `:sol`, the custom nonlinear objective `:JE`, the nonlinear
+constraint vectors and the various derivatives.
 """
 function addinfo!(info, mpc::NonLinMPC{NT}) where NT<:Real
+    # --- variables specific to NonLinMPC ---
     U, Ŷ, D̂, ŷ, d, ϵ = info[:U], info[:Ŷ], info[:D̂], info[:ŷ], info[:d], info[:ϵ]
     Ue = [U; U[(end - mpc.estim.model.nu + 1):end]]
     Ŷe = [ŷ; Ŷ]
@@ -539,6 +541,120 @@ function addinfo!(info, mpc::NonLinMPC{NT}) where NT<:Real
     info[:JE]  = JE 
     info[:gc] = LHS
     info[:sol] = JuMP.solution_summary(mpc.optim, verbose=true)
+    # --- objective derivatives ---
+    model, optim, con = mpc.estim.model, mpc.optim, mpc.con
+    transcription = mpc.transcription
+    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(con.i_g)
+    nc, neq = con.nc, con.neq
+    nU, nŶ, nX̂, nK = mpc.Hp*nu, Hp*ny, Hp*nx̂, Hp*nk
+    nΔŨ, nUe, nŶe = nu*Hc + nϵ, nU + nu, nŶ + ny  
+    ΔŨ     = zeros(NT, nΔŨ)
+    x̂0end  = zeros(NT, nx̂)
+    K0     = zeros(NT, nK)
+    Ue, Ŷe = zeros(NT, nUe), zeros(NT, nŶe)
+    U0, Ŷ0 = zeros(NT, nU),  zeros(NT, nŶ)
+    Û0, X̂0 = zeros(NT, nU),  zeros(NT, nX̂)
+    gc, g  = zeros(NT, nc),  zeros(NT, ng)
+    geq    = zeros(NT, neq)
+    J_cache = (
+        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),
+    )
+    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
+    if !isnothing(mpc.hessian)
+        _, ∇J, ∇²J = value_gradient_and_hessian(J!, mpc.hessian, mpc.Z̃, J_cache...)
+    else
+        ∇J, ∇²J = gradient(J!, mpc.gradient, mpc.Z̃, J_cache...), nothing
+    end
+    # --- inequality constraint derivatives ---
+    old_i_g = copy(con.i_g)
+    con.i_g .= 1 # temporarily set all constraints as finite so g is entirely computed
+    ∇g_cache = (
+        Cache(ΔŨ), Cache(x̂0end), Cache(Ue), Cache(Ŷe), Cache(U0), Cache(Ŷ0), 
+        Cache(Û0), Cache(K0), Cache(X̂0), 
+        Cache(gc), Cache(geq),
+    )
+    function g!(g, Z̃, ΔŨ, x̂0end, Ue, Ŷe, U0, Ŷ0, Û0, K0, X̂0, gc, geq)
+        update_predictions!(ΔŨ, x̂0end, Ue, Ŷe, U0, Ŷ0, Û0, K0, X̂0, gc, g, geq, mpc, Z̃)
+        return nothing
+    end
+    g, ∇g = value_and_jacobian(g!, g, mpc.jacobian, mpc.Z̃, ∇g_cache...)
+    if !isnothing(mpc.hessian) && any(old_i_g)
+        @warn(
+            "Retrieving optimal Hessian of the Lagrangian is not fully supported yet.\n"*
+            "Its nonzero coefficients are random values for now.", maxlog=1
+        )
+        function ℓ_g(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̃)
+            return dot(λ, g)
+        end
+        ∇²g_cache = (
+            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)
+        )
+        nonlincon = optim[:nonlinconstraint]
+        λ = JuMP.dual.(nonlincon) # FIXME: does not work for now
+        λ = rand(NT, ng)
+        ∇²ℓg = hessian(ℓ_g, mpc.hessian, mpc.Z̃, Constant(λ), ∇²g_cache...)
+    else
+        ∇²ℓg = nothing
+    end
+    con.i_g .= old_i_g # restore original finite/infinite constraint indices
+    # --- equality constraint derivatives ---
+    geq_cache = (
+        Cache(ΔŨ), Cache(x̂0end), Cache(Ue), Cache(Ŷe), Cache(U0), Cache(Ŷ0),
+        Cache(Û0), Cache(K0),   Cache(X̂0),
+        Cache(gc), Cache(g)
+    )
+    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
+    geq, ∇geq = value_and_jacobian(geq!, geq, mpc.jacobian, mpc.Z̃, geq_cache...)
+    if !isnothing(mpc.hessian) && con.neq > 0
+        @warn(
+            "Retrieving optimal Hessian of the Lagrangian is not fully supported yet.\n"*
+            "Its nonzero coefficients are random values for now.", maxlog=1
+        )
+        ∇²geq_cache = (
+            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)
+        )
+        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
+        nonlinconeq = optim[:nonlinconstrainteq]
+        λeq = JuMP.dual.(nonlinconeq) # FIXME: does not work for now
+        λeq = ones(NT, neq)
+        ∇²ℓgeq = hessian(ℓ_geq, mpc.hessian, mpc.Z̃, Constant(λeq), ∇²geq_cache...)
+    else
+        ∇²ℓgeq = nothing
+    end
+    info[:∇J] = ∇J
+    info[:∇²J] = ∇²J
+    info[:g] = g
+    info[:∇g] = ∇g
+    info[:∇²ℓg] = ∇²ℓg
+    info[:geq] = geq
+    info[:∇geq] = ∇geq
+    info[:∇²ℓgeq] = ∇²ℓgeq
+    # --- non-Unicode fields ---
+    info[:nablaJ] = ∇J
+    info[:nabla2J] = ∇²J
+    info[:nablag] = ∇g
+    info[:nabla2lg] = ∇²ℓg
+    info[:nablageq] = ∇geq
+    info[:nabla2lgeq] = ∇²ℓgeq
     return info
 end
 
@@ -942,10 +1058,10 @@ function set_nonlincon!(
         optim, JuMP.Vector{JuMP.VariableRef}, MOI.VectorNonlinearOracle{JNT}
     )
     map(con_ref -> JuMP.delete(optim, con_ref), nonlin_constraints)
-    optim[:g_oracle]   = g_oracle
-    optim[:geq_oracle] = geq_oracle
-    any(mpc.con.i_g) && @constraint(optim, Z̃var in g_oracle)
-    mpc.con.neq > 0  && @constraint(optim, Z̃var in geq_oracle)
+    JuMP.unregister(optim, :nonlinconstraint)
+    JuMP.unregister(optim, :nonlinconstrainteq)
+    any(mpc.con.i_g) && @constraint(optim, nonlinconstraint, Z̃var in g_oracle)
+    mpc.con.neq > 0  && @constraint(optim, nonlinconstrainteq, Z̃var in geq_oracle)
     return nothing
 end
 

src/estimator/mhe/construct.jl

@@ -1504,9 +1504,7 @@ function get_nonlincon_oracle(
         return dot(λi, gi)
     end
     Z̃_∇gi = fill(myNaN, nZ̃)      # NaN to force update_predictions! at first call
-    ∇gi_cache = (
-        Cache(V̂), Cache(X̂0), Cache(û0), Cache(k0), Cache(ŷ0), Cache(g)
-    )
+    ∇gi_cache = (Cache(V̂), Cache(X̂0), Cache(û0), Cache(k0), Cache(ŷ0), Cache(g))
     # temporarily "fill" the estimation window for the preparation of the gradient: 
     estim.Nk[] = He
     ∇gi_prep = prepare_jacobian(gi!, gi, jac, Z̃_∇gi, ∇gi_cache...; strict)
@@ -1577,7 +1575,7 @@ function set_nonlincon!(
         optim, JuMP.Vector{JuMP.VariableRef}, MOI.VectorNonlinearOracle{JNT}
     )
     map(con_ref -> JuMP.delete(optim, con_ref), nonlin_constraints)
-    optim[:g_oracle]   = g_oracle
-    any(estim.con.i_g) && @constraint(optim, Z̃var in g_oracle)
+    JuMP.unregister(optim, :nonlinconstraint)
+    any(estim.con.i_g) && @constraint(optim, nonlinconstraint, Z̃var in g_oracle)
     return nothing
 end