Merge pull request #200 from JuliaControl/cherry_pick

Cherry pick doc and general improvements commits from `hessian_objective`

Author: franckgaga

docs/Project.toml

@@ -16,5 +16,5 @@ Documenter = "1"
 JuMP = "1"
 LinearAlgebra = "1.10"
 Logging = "1.10"
-ModelingToolkit = "9.50"
+ModelingToolkit = "9.50 - 9.76"
 Plots = "1"

docs/src/manual/mtk.md

@@ -23,7 +23,8 @@ the last section.
     will work for all corner cases.
 
 !!! compat
-    The example relies on features and bugfixes of `ModelingToolkit.jl` v9.50.
+    The example works on `ModelingToolkit.jl` v9.50 to 9.76 (corresponding to the following
+    `[compat]` entry: `ModelingToolkit = "9.50 - 9.76"`).
 
 We first construct and instantiate the pendulum model:
 

src/controller/execute.jl

@@ -395,7 +395,7 @@ function optim_objective!(mpc::PredictiveController{NT}) where {NT<:Real}
     model, optim = mpc.estim.model, mpc.optim
     nu, Hc = model.nu, mpc.Hc 
     Z̃var::Vector{JuMP.VariableRef} = optim[:Z̃var]
-    Z̃0 = set_warmstart!(mpc, mpc.transcription, Z̃var)
+    Z̃s = set_warmstart!(mpc, mpc.transcription, Z̃var)
     set_objective_linear_coef!(mpc, Z̃var)
     try
         JuMP.optimize!(optim)
@@ -405,7 +405,7 @@ function optim_objective!(mpc::PredictiveController{NT}) where {NT<:Real}
             MOIU.reset_optimizer(optim)
             JuMP.optimize!(optim)
         else
-            rethrow(err)
+            rethrow()
         end
     end
     if !issolved(optim)
@@ -426,7 +426,7 @@ function optim_objective!(mpc::PredictiveController{NT}) where {NT<:Real}
         @debug info2debugstr(getinfo(mpc))
     end
     if iserror(optim)
-        mpc.Z̃ .= Z̃0
+        mpc.Z̃ .= Z̃s
     else
         mpc.Z̃ .= JuMP.value.(Z̃var)
     end
@@ -488,7 +488,7 @@ Call `periodsleep(mpc.estim.model)`.
 periodsleep(mpc::PredictiveController, busywait=false) = periodsleep(mpc.estim.model, busywait)
 
 """
-    setstate!(mpc::PredictiveController, x̂, P̂=nothing) -> mpc
+    setstate!(mpc::PredictiveController, x̂[, P̂]) -> mpc
 
 Call [`setstate!`](@ref) on `mpc.estim` [`StateEstimator`](@ref).
 """

src/controller/transcription.jl

@@ -972,14 +972,14 @@ where ``\mathbf{Δu}(k+j|k-1)`` is the input increment for time ``k+j`` computed
 last control period ``k-1``, and ``ϵ(k-1)``, the slack variable of the last control period.
 """
 function set_warmstart!(mpc::PredictiveController, transcription::SingleShooting, Z̃var)
-    nu, Hc, Z̃0 = mpc.estim.model.nu, mpc.Hc, mpc.buffer.Z̃
+    nu, Hc, Z̃s = mpc.estim.model.nu, mpc.Hc, mpc.buffer.Z̃
     # --- input increments ΔU ---
-    Z̃0[1:(Hc*nu-nu)] .= @views mpc.Z̃[nu+1:Hc*nu]
-    Z̃0[(Hc*nu-nu+1):(Hc*nu)] .= 0
+    Z̃s[1:(Hc*nu-nu)] .= @views mpc.Z̃[nu+1:Hc*nu]
+    Z̃s[(Hc*nu-nu+1):(Hc*nu)] .= 0
     # --- slack variable ϵ ---
-    mpc.nϵ == 1 && (Z̃0[end] = mpc.Z̃[end])
-    JuMP.set_start_value.(Z̃var, Z̃0)
-    return Z̃0
+    mpc.nϵ == 1 && (Z̃s[end] = mpc.Z̃[end])
+    JuMP.set_start_value.(Z̃var, Z̃s)
+    return Z̃s
 end
 
 @doc raw"""
@@ -1009,17 +1009,17 @@ last control period ``k-1``, expressed as a deviation from the operating point
 ``\mathbf{x̂_{op}}``.
 """
 function set_warmstart!(mpc::PredictiveController, transcription::MultipleShooting, Z̃var)
-    nu, nx̂, Hp, Hc, Z̃0 = mpc.estim.model.nu, mpc.estim.nx̂, mpc.Hp, mpc.Hc, mpc.buffer.Z̃
+    nu, nx̂, Hp, Hc, Z̃s = mpc.estim.model.nu, mpc.estim.nx̂, mpc.Hp, mpc.Hc, mpc.buffer.Z̃
     # --- input increments ΔU ---
-    Z̃0[1:(Hc*nu-nu)] .= @views mpc.Z̃[nu+1:Hc*nu]
-    Z̃0[(Hc*nu-nu+1):(Hc*nu)] .= 0
+    Z̃s[1:(Hc*nu-nu)] .= @views mpc.Z̃[nu+1:Hc*nu]
+    Z̃s[(Hc*nu-nu+1):(Hc*nu)] .= 0
     # --- predicted states X̂0 ---
-    Z̃0[(Hc*nu+1):(Hc*nu+Hp*nx̂-nx̂)]       .= @views mpc.Z̃[(Hc*nu+nx̂+1):(Hc*nu+Hp*nx̂)]
-    Z̃0[(Hc*nu+Hp*nx̂-nx̂+1):(Hc*nu+Hp*nx̂)] .= @views mpc.Z̃[(Hc*nu+Hp*nx̂-nx̂+1):(Hc*nu+Hp*nx̂)]
+    Z̃s[(Hc*nu+1):(Hc*nu+Hp*nx̂-nx̂)]       .= @views mpc.Z̃[(Hc*nu+nx̂+1):(Hc*nu+Hp*nx̂)]
+    Z̃s[(Hc*nu+Hp*nx̂-nx̂+1):(Hc*nu+Hp*nx̂)] .= @views mpc.Z̃[(Hc*nu+Hp*nx̂-nx̂+1):(Hc*nu+Hp*nx̂)]
     # --- slack variable ϵ ---
-    mpc.nϵ == 1 && (Z̃0[end] = mpc.Z̃[end])
-    JuMP.set_start_value.(Z̃var, Z̃0)
-    return Z̃0
+    mpc.nϵ == 1 && (Z̃s[end] = mpc.Z̃[end])
+    JuMP.set_start_value.(Z̃var, Z̃s)
+    return Z̃s
 end
 
 getΔŨ!(ΔŨ, mpc::PredictiveController, ::SingleShooting, Z̃) = (ΔŨ .= Z̃) # since mpc.P̃Δu = I

src/estimator/execute.jl

@@ -330,7 +330,7 @@ function validate_args(estim::StateEstimator, ym, d, u=nothing)
 end
 
 """
-    setstate!(estim::StateEstimator, x̂, P̂=nothing) -> estim
+    setstate!(estim::StateEstimator, x̂[, P̂]) -> estim
 
 Set `estim.x̂0` to `x̂ - estim.x̂op` from the argument `x̂`, and `estim.P̂` to `P̂` if applicable. 
 

src/estimator/mhe/execute.jl

@@ -364,7 +364,7 @@ Optimize objective of `estim` [`MovingHorizonEstimator`](@ref) and return the so
 
 If supported by `estim.optim`, it warm-starts the solver at:
 ```math
-\mathbf{Z̃} = 
+\mathbf{Z̃_s} = 
 \begin{bmatrix}
     ϵ_{k-1}                         \\
     \mathbf{x̂}_{k-1}(k-N_k+p)       \\ 
@@ -391,12 +391,12 @@ function optim_objective!(estim::MovingHorizonEstimator{NT}) where NT<:Real
     X̂0  = Vector{NT}(undef, nx̂*Nk)
     û0, ŷ0, x̄, k0 = buffer.û, buffer.ŷ, buffer.x̂, buffer.k
     ϵ_0 = estim.nϵ ≠ 0 ? estim.Z̃[begin] : empty(estim.Z̃)
-    Z̃_0 = [ϵ_0; estim.x̂0arr_old; estim.Ŵ]
-    V̂, X̂0 = predict!(V̂, X̂0, û0, k0, ŷ0, estim, model, Z̃_0)
-    J_0 = obj_nonlinprog!(x̄, estim, model, V̂, Z̃_0)
-    # initial Z̃0 with Ŵ=0 if objective or constraint function not finite :
-    isfinite(J_0) || (Z̃_0 = [ϵ_0; estim.x̂0arr_old; zeros(NT, nŵ*estim.He)])
-    JuMP.set_start_value.(Z̃var, Z̃_0)
+    Z̃s = [ϵ_0; estim.x̂0arr_old; estim.Ŵ]
+    V̂, X̂0 = predict!(V̂, X̂0, û0, k0, ŷ0, estim, model, Z̃s)
+    J_0 = obj_nonlinprog!(x̄, estim, model, V̂, Z̃s)
+    # warm-start Z̃s with Ŵ=0 if objective or constraint function not finite :
+    isfinite(J_0) || (Z̃s = [ϵ_0; estim.x̂0arr_old; zeros(NT, nŵ*estim.He)])
+    JuMP.set_start_value.(Z̃var, Z̃s)
     # ------- solve optimization problem --------------
     try
         JuMP.optimize!(optim)
@@ -428,7 +428,7 @@ function optim_objective!(estim::MovingHorizonEstimator{NT}) where NT<:Real
         @debug info2debugstr(getinfo(estim))
     end
     if iserror(optim)
-        estim.Z̃ .= Z̃_0
+        estim.Z̃ .= Z̃s
     else
         estim.Z̃ .= JuMP.value.(Z̃var)
     end

src/general.jl

@@ -51,7 +51,7 @@ function limit_solve_time(optim::GenericModel, Ts)
         if isa(err, MOI.UnsupportedAttribute{MOI.TimeLimitSec})
             @warn "Solving time limit is not supported by the optimizer."
         else
-            rethrow(err)
+            rethrow()
         end
     end
 end

test/2_test_state_estim.jl

@@ -1364,7 +1364,7 @@ end
     X̂_mhe = zeros(4, 6)
     X̂_kf  = zeros(4, 6)
     for i in 1:6
-        y = [50,31] #+ randn(2)
+        y = [50,31] + randn(2)
         x̂_mhe = preparestate!(mhe, y, [25])
         x̂_kf  = preparestate!(kf,  y, [25])
         X̂_mhe[:,i] = x̂_mhe