doc: more consise setconstraint! and ExplicitMPC doc

Author: franckgaga

src/controller/construct.jl

@@ -52,25 +52,17 @@ for details on bounds and softness parameters ``\mathbf{c}``. The output and ter
 constraints are all soft by default. See Extended Help for time-varying constraints.
 
 # Arguments
-- `mpc::PredictiveController` : predictive controller to set constraints.
-- `umin  = fill(-Inf,nu)` : manipulated input lower bounds ``\mathbf{u_{min}}``.
-- `umax  = fill(+Inf,nu)` : manipulated input upper bounds ``\mathbf{u_{max}}``.
-- `Δumin = fill(-Inf,nu)` : manipulated input increment lower bounds ``\mathbf{Δu_{min}}``.
-- `Δumax = fill(+Inf,nu)` : manipulated input increment upper bounds ``\mathbf{Δu_{max}}``.
-- `ymin  = fill(-Inf,ny)` : predicted output lower bounds ``\mathbf{y_{min}}``.
-- `ymax  = fill(+Inf,ny)` : predicted output upper bounds ``\mathbf{y_{max}}``.
-- `x̂min  = fill(-Inf,nx̂)` : terminal constraint lower bounds ``\mathbf{x̂_{min}}``.
-- `x̂max  = fill(+Inf,nx̂)` : terminal constraint upper bounds ``\mathbf{x̂_{max}}``.
-- `c_umin  = fill(0.0,nu)` : `umin` softness weights ``\mathbf{c_{u_{min}}}``.
-- `c_umax  = fill(0.0,nu)` : `umax` softness weights ``\mathbf{c_{u_{max}}}``.
-- `c_Δumin = fill(0.0,nu)` : `Δumin` softness weights ``\mathbf{c_{Δu_{min}}}``.
-- `c_Δumax = fill(0.0,nu)` : `Δumax` softness weights ``\mathbf{c_{Δu_{max}}}``.
-- `c_ymin  = fill(1.0,ny)` : `ymin` softness weights ``\mathbf{c_{y_{min}}}``.
-- `c_ymax  = fill(1.0,ny)` : `ymax` softness weights ``\mathbf{c_{y_{max}}}``.
-- `c_x̂min  = fill(1.0,nx̂)` : `x̂min` softness weights ``\mathbf{c_{x̂_{min}}}``.
-- `c_x̂max  = fill(1.0,nx̂)` : `x̂max` softness weights ``\mathbf{c_{x̂_{max}}}``.
+- `mpc::PredictiveController` : predictive controller to set constraints
+- `umin=fill(-Inf,nu)` / `umax=fill(+Inf,nu)` : manipulated input bound ``\mathbf{u_{min/max}}``
+- `Δumin=fill(-Inf,nu)` / `Δumax=fill(+Inf,nu)` : manipulated input increment bound ``\mathbf{Δu_{min/max}}``
+- `ymin=fill(-Inf,ny)` / `ymax=fill(+Inf,ny)` : predicted output bound ``\mathbf{y_{min/max}}``
+- `x̂min=fill(-Inf,nx̂)` / `x̂max=fill(+Inf,nx̂)` : terminal constraint bound ``\mathbf{x̂_{min/max}}``
+- `c_umin=fill(0.0,nu)` / `c_umax=fill(0.0,nu)` : `umin` / `umax` softness weight ``\mathbf{c_{u_{min/max}}}``
+- `c_Δumin=fill(0.0,nu)` / `c_Δumax=fill(0.0,nu)` : `Δumin` / `Δumax` softness weight ``\mathbf{c_{Δu_{min/max}}}``
+- `c_ymin=fill(1.0,ny)` / `c_ymax=fill(1.0,ny)` : `ymin` / `ymax` softness weight ``\mathbf{c_{y_{min/max}}}``
+- `c_x̂min=fill(1.0,nx̂)` / `c_x̂max=fill(1.0,nx̂)` : `x̂min` / `x̂max` softness weight ``\mathbf{c_{x̂_{min/max}}}``
 - all the keyword arguments above but with a first capital letter, except for the terminal
-  constraints, e.g. `Ymax` or `C_Δumin`: for time-varying constraints (see Extended Help).
+  constraints, e.g. `Ymax` or `C_Δumin`: for time-varying constraints (see Extended Help)
 
 # Examples
 ```jldoctest
@@ -134,31 +126,7 @@ function setconstraint!(
     C_umax  = nothing, C_umin  = nothing,
     C_Δumax = nothing, C_Δumin = nothing,
     C_ymax  = nothing, C_ymin  = nothing,
-    # TODO:
-    # ------------ will be deleted in the future ---------------
-    ŷmin    = nothing, ŷmax    = nothing,
-    c_ŷmin  = nothing, c_ŷmax  = nothing,
-    # ----------------------------------------------------------
 )
-    # TODO:
-    # ----- these 4 `if`s will be deleted in the future --------
-    if !isnothing(ŷmin)
-        Base.depwarn("keyword arg ŷmin is deprecated, use ymin instead", :setconstraint!)
-        ymin = ŷmin
-    end
-    if !isnothing(ŷmax)
-        Base.depwarn("keyword arg ŷmax is deprecated, use ymax instead", :setconstraint!)
-        ymax = ŷmax
-    end
-    if !isnothing(c_ŷmin)
-        Base.depwarn("keyword arg ŷmin is deprecated, use ymin instead", :setconstraint!)
-        c_ymin = c_ŷmin
-    end
-    if !isnothing(c_ŷmax)
-        Base.depwarn("keyword arg ŷmax is deprecated, use ymax instead", :setconstraint!)
-        c_ymax = c_ŷmax
-    end
-    # ----------------------------------------------------------
     model, con, optim = mpc.estim.model, mpc.con, mpc.optim
     nu, ny, nx̂, Hp, Hc = model.nu, model.ny, mpc.estim.nx̂, mpc.Hp, mpc.Hc
     notSolvedYet = (termination_status(optim) == OPTIMIZE_NOT_CALLED)
@@ -339,11 +307,22 @@ end
 setnonlincon!(::PredictiveController, ::SimModel) = nothing
 
 """
-    default_Hp(model::LinModel, Hp)
+    default_Hp(model::LinModel)
+
+Estimate the default prediction horizon `Hp` for [`LinModel`](@ref).
+"""
+default_Hp(model::LinModel) = DEFAULT_HP0 + estimate_delays(model)
+"Throw an error when model is not a [`LinModel`](@ref)."
+function default_Hp(::SimModel)
+    throw(ArgumentError("Prediction horizon Hp must be explicitly specified if model is not a LinModel."))
+end
+
+"""
+    estimate_delays(model::LinModel)
 
-Estimate the default prediction horizon `Hp` with a security margin for [`LinModel`](@ref).
+Estimate the number of delays in `model` with a security margin.
 """
-function default_Hp(model::LinModel, Hp)
+function estimate_delays(model::LinModel)
     # TODO: also check for settling time (poles)
     # TODO: also check for non minimum phase systems (zeros)
     # TODO: replace sum with max delay between all the I/O
@@ -352,36 +331,10 @@ function default_Hp(model::LinModel, Hp)
     # atol=1e-3 to overestimate the number of delays : for closed-loop stability, it is
     # better to overestimate the default value of Hp, as a security margin.
     nk = sum(isapprox.(abs.(poles), 0.0, atol=1e-3)) # number of delays
-    if isnothing(Hp)
-        Hp = DEFAULT_HP0 + nk
-    end
-    if Hp ≤ nk
-        @warn("prediction horizon Hp ($Hp) ≤ estimated number of delays in model "*
-              "($nk), the closed-loop system may be unstable or zero-gain (unresponsive)")
-    end
-    return Hp
-end
-
-"""
-    default_Hp(model::SimModel, Hp)
-
-Throw an error if `isnothing(Hp)` when model is not a [`LinModel`](@ref).
-"""
-function default_Hp(::SimModel, Hp)
-    if isnothing(Hp)
-        # TODO:
-        # ------------ will be deleted in the future ------------------------------------
-        Base.depwarn("Hp=nothing is deprecated for NonLinModel, explicitly specify an "*
-                     "integer value", :NonLinMPC)
-        Hp = DEFAULT_HP0
-        # ------------- and replaced by this -------------------------------------------
-        # throw(ArgumentError("Prediction horizon Hp must be explicitly specified if "*
-        #                     "model is not a LinModel."))
-        # Hp = 0
-        # -----------------------------------------------------------------------------
-    end
-    return Hp
+    return nk
 end
+"Return `0` when model is not a [`LinModel`](@ref)."
+estimate_delays(::SimModel) = 0
 
 """
     validate_args(mpc::PredictiveController, ry, d, D̂, R̂y, R̂u)

src/controller/execute.jl

@@ -73,19 +73,19 @@ Get additional info about `mpc` [`PredictiveController`](@ref) optimum for troub
 The function should be called after calling [`moveinput!`](@ref). It returns the dictionary
 `info` with the following fields:
 
-- `:ΔU`  : optimal manipulated input increments over ``H_c``, ``\mathbf{ΔU}``.
-- `:ϵ`   : optimal slack variable, ``ϵ``.
-- `:J`   : objective value optimum, ``J``.
-- `:U`   : optimal manipulated inputs over ``H_p``, ``\mathbf{U}``.
-- `:u`   : current optimal manipulated input, ``\mathbf{u}(k)``.
-- `:d`   : current measured disturbance, ``\mathbf{d}(k)``.
-- `:D̂`   : predicted measured disturbances over ``H_p``, ``\mathbf{D̂}``.
-- `:ŷ`   : current estimated output, ``\mathbf{ŷ}(k)``.
-- `:Ŷ`   : optimal predicted outputs over ``H_p``, ``\mathbf{Ŷ}``.
-- `:x̂end`: optimal terminal states, ``\mathbf{x̂}_{k-1}(k+H_p)``.
-- `:Ŷs`  : predicted stochastic output over ``H_p`` of [`InternalModel`](@ref), ``\mathbf{Ŷ_s}``.
-- `:R̂y`  : predicted output setpoint over ``H_p``, ``\mathbf{R̂_y}``.
-- `:R̂u`  : predicted manipulated input setpoint over ``H_p``, ``\mathbf{R̂_u}``.
+- `:ΔU`  : optimal manipulated input increments over ``H_c``, ``\mathbf{ΔU}``
+- `:ϵ`   : optimal slack variable, ``ϵ``
+- `:J`   : objective value optimum, ``J``
+- `:U`   : optimal manipulated inputs over ``H_p``, ``\mathbf{U}``
+- `:u`   : current optimal manipulated input, ``\mathbf{u}(k)``
+- `:d`   : current measured disturbance, ``\mathbf{d}(k)``
+- `:D̂`   : predicted measured disturbances over ``H_p``, ``\mathbf{D̂}``
+- `:ŷ`   : current estimated output, ``\mathbf{ŷ}(k)``
+- `:Ŷ`   : optimal predicted outputs over ``H_p``, ``\mathbf{Ŷ}``
+- `:x̂end`: optimal terminal states, ``\mathbf{x̂}_{k-1}(k+H_p)``
+- `:Ŷs`  : predicted stochastic output over ``H_p`` of [`InternalModel`](@ref), ``\mathbf{Ŷ_s}``
+- `:R̂y`  : predicted output setpoint over ``H_p``, ``\mathbf{R̂_y}``
+- `:R̂u`  : predicted manipulated input setpoint over ``H_p``, ``\mathbf{R̂_u}``
 
 For [`LinMPC`](@ref) and [`NonLinMPC`](@ref), the field `:sol` also contains the optimizer
 solution summary that can be printed. Lastly, the optimal economic cost `:JE` is also

src/controller/explicitmpc.jl

@@ -40,7 +40,7 @@ struct ExplicitMPC{NT<:Real, SE<:StateEstimator} <: PredictiveController{NT}
         Cwt = Inf # no slack variable ϵ for ExplicitMPC
         Ewt = 0   # economic costs not supported for ExplicitMPC
         validate_weights(model, Hp, Hc, M_Hp, N_Hc, L_Hp, Cwt)
-        M_Hp, N_Hc, L_Hp = float(M_Hp), float(N_Hc), float(L_Hp) # debug julia 1.6
+        M_Hp, N_Hc, L_Hp = Diagonal{NT}(M_Hp), Diagonal{NT}(N_Hc), Diagonal{NT}(L_Hp) # debug julia 1.6
         # dummy vals (updated just before optimization):
         R̂y, R̂u = zeros(NT, ny*Hp), zeros(NT, nu*Hp)
         noR̂u = iszero(L_Hp)
@@ -88,19 +88,11 @@ The controller minimizes the following objective function at each discrete time
 
 See [`LinMPC`](@ref) for the variable definitions. This controller does not support
 constraints but the computational costs are extremely low (array division), therefore 
-suitable for applications that require small sample times. This method uses the default
-state estimator, a [`SteadyKalmanFilter`](@ref) with default arguments.
-
-# Arguments
-- `model::LinModel` : model used for controller predictions and state estimations.
-- `Hp=10+nk`: prediction horizon ``H_p``, `nk` is the number of delays in `model`.
-- `Hc=2` : control horizon ``H_c``.
-- `Mwt=fill(1.0,model.ny)` : main diagonal of ``\mathbf{M}`` weight matrix (vector).
-- `Nwt=fill(0.1,model.nu)` : main diagonal of ``\mathbf{N}`` weight matrix (vector).
-- `Lwt=fill(0.0,model.nu)` : main diagonal of ``\mathbf{L}`` weight matrix (vector).
-- `M_Hp` / `N_Hc` / `L_Hp` : diagonal matrices ``\mathbf{M}_{H_p}, \mathbf{N}_{H_c},
-  \mathbf{L}_{H_p}``, for time-varying weights (generated from `Mwt/Nwt/Lwt` args if omitted).
-- additional keyword arguments are passed to [`SteadyKalmanFilter`](@ref) constructor.
+suitable for applications that require small sample times. The keyword arguments are
+identical to [`LinMPC`](@ref), except for `Cwt` and `optim` which are not supported. 
+
+This method uses the default state estimator, a [`SteadyKalmanFilter`](@ref) with default
+arguments.
 
 # Examples
 ```jldoctest
@@ -120,14 +112,14 @@ ExplicitMPC controller with a sample time Ts = 4.0 s, SteadyKalmanFilter estimat
 """
 function ExplicitMPC(
     model::LinModel; 
-    Hp::Union{Int, Nothing} = nothing,
+    Hp::Int = default_Hp(model),
     Hc::Int = DEFAULT_HC,
     Mwt = fill(DEFAULT_MWT, model.ny),
     Nwt = fill(DEFAULT_NWT, model.nu),
     Lwt = fill(DEFAULT_LWT, model.nu),
-    M_Hp = nothing,
-    N_Hc = nothing,
-    L_Hp = nothing,
+    M_Hp = Diagonal(repeat(Mwt, Hp)),
+    N_Hc = Diagonal(repeat(Nwt, Hc)),
+    L_Hp = Diagonal(repeat(Lwt, Hp)),
     kwargs...
 ) 
     estim = SteadyKalmanFilter(model; kwargs...)
@@ -158,20 +150,21 @@ ExplicitMPC controller with a sample time Ts = 4.0 s, KalmanFilter estimator and
 """
 function ExplicitMPC(
     estim::SE;
-    Hp::Union{Int, Nothing} = nothing,
+    Hp::Int = default_Hp(estim.model),
     Hc::Int = DEFAULT_HC,
-    Mwt = fill(DEFAULT_MWT, estim.model.ny),
-    Nwt = fill(DEFAULT_NWT, estim.model.nu),
-    Lwt = fill(DEFAULT_LWT, estim.model.nu),
-    M_Hp = nothing,
-    N_Hc = nothing,
-    L_Hp = nothing
+    Mwt  = fill(DEFAULT_MWT, estim.model.ny),
+    Nwt  = fill(DEFAULT_NWT, estim.model.nu),
+    Lwt  = fill(DEFAULT_LWT, estim.model.nu),
+    M_Hp = Diagonal(repeat(Mwt, Hp)),
+    N_Hc = Diagonal(repeat(Nwt, Hc)),
+    L_Hp = Diagonal(repeat(Lwt, Hp)),
 ) where {NT<:Real, SE<:StateEstimator{NT}}
-    isa(estim.model, LinModel) || error("estim.model type must be LinModel") 
-    Hp = default_Hp(estim.model, Hp)
-    isnothing(M_Hp) && (M_Hp = Diagonal{NT}(repeat(Mwt, Hp)))
-    isnothing(N_Hc) && (N_Hc = Diagonal{NT}(repeat(Nwt, Hc)))
-    isnothing(L_Hp) && (L_Hp = Diagonal{NT}(repeat(Lwt, Hp)))
+    isa(estim.model, LinModel) || error("estim.model type must be a LinModel") 
+    nk = estimate_delays(estim.model)
+    if Hp ≤ nk
+        @warn("prediction horizon Hp ($Hp) ≤ estimated number of delays in model "*
+              "($nk), the closed-loop system may be unstable or zero-gain (unresponsive)")
+    end
     return ExplicitMPC{NT, SE}(estim, Hp, Hc, M_Hp, N_Hc, L_Hp)
 end
 

src/controller/linmpc.jl

@@ -48,7 +48,7 @@ struct LinMPC{
         ŷ = copy(model.yop) # dummy vals (updated just before optimization)
         Ewt = 0   # economic costs not supported for LinMPC
         validate_weights(model, Hp, Hc, M_Hp, N_Hc, L_Hp, Cwt)
-        M_Hp, N_Hc, L_Hp = float(M_Hp), float(N_Hc), float(L_Hp) # debug julia 1.6
+        M_Hp, N_Hc, L_Hp = Diagonal{NT}(M_Hp), Diagonal{NT}(N_Hc), Diagonal{NT}(L_Hp) # debug julia 1.6
         # dummy vals (updated just before optimization):
         R̂y, R̂u = zeros(NT, ny*Hp), zeros(NT, nu*Hp) 
         noR̂u = iszero(L_Hp)
@@ -93,7 +93,7 @@ The controller minimizes the following objective function at each discrete time
                       + C ϵ^2
 \end{aligned}
 ```
-in which the weight matrices are repeated ``H_p`` or ``H_c`` times:
+in which the weight matrices are repeated ``H_p`` or ``H_c`` times by default:
 ```math
 \begin{aligned}
     \mathbf{M}_{H_p} &= \text{diag}\mathbf{(M,M,...,M)}     \\
@@ -118,12 +118,14 @@ arguments.
 - `Mwt=fill(1.0,model.ny)` : main diagonal of ``\mathbf{M}`` weight matrix (vector).
 - `Nwt=fill(0.1,model.nu)` : main diagonal of ``\mathbf{N}`` weight matrix (vector).
 - `Lwt=fill(0.0,model.nu)` : main diagonal of ``\mathbf{L}`` weight matrix (vector).
+- `M_Hp=Diagonal(repeat(Mwt),Hp)` : diagonal weight matrix ``\mathbf{M}_{H_p}``.
+- `N_Hc=Diagonal(repeat(Nwt),Hc)` : diagonal weight matrix ``\mathbf{N}_{H_c}``.
+- `L_Hp=Diagonal(repeat(Lwt),Hp)` : diagonal weight matrix ``\mathbf{L}_{H_p}``.
 - `Cwt=1e5` : slack variable weight ``C`` (scalar), use `Cwt=Inf` for hard constraints only.
 - `optim=JuMP.Model(OSQP.MathOptInterfaceOSQP.Optimizer)` : quadratic optimizer used in
   the predictive controller, provided as a [`JuMP.Model`](https://jump.dev/JuMP.jl/stable/api/JuMP/#JuMP.Model)
   (default to [`OSQP`](https://osqp.org/docs/parsers/jump.html) optimizer).
-- `M_Hp` / `N_Hc` / `L_Hp` : diagonal matrices ``\mathbf{M}_{H_p}, \mathbf{N}_{H_c},
-  \mathbf{L}_{H_p}``, for time-varying weights (generated from `Mwt/Nwt/Lwt` args if omitted).
+
 - additional keyword arguments are passed to [`SteadyKalmanFilter`](@ref) constructor.
 
 # Examples
@@ -166,15 +168,15 @@ LinMPC controller with a sample time Ts = 4.0 s, OSQP optimizer, SteadyKalmanFil
 """
 function LinMPC(
     model::LinModel;
-    Hp::Union{Int, Nothing} = nothing,
+    Hp::Int = default_Hp(model),
     Hc::Int = DEFAULT_HC,
-    Mwt = fill(DEFAULT_MWT, model.ny),
-    Nwt = fill(DEFAULT_NWT, model.nu),
-    Lwt = fill(DEFAULT_LWT, model.nu),
+    Mwt  = fill(DEFAULT_MWT, model.ny),
+    Nwt  = fill(DEFAULT_NWT, model.nu),
+    Lwt  = fill(DEFAULT_LWT, model.nu),
+    M_Hp = Diagonal(repeat(Mwt, Hp)),
+    N_Hc = Diagonal(repeat(Nwt, Hc)),
+    L_Hp = Diagonal(repeat(Lwt, Hp)),
     Cwt = DEFAULT_CWT,
-    M_Hp = nothing,
-    N_Hc = nothing,
-    L_Hp = nothing,
     optim::JuMP.GenericModel = JuMP.Model(DEFAULT_LINMPC_OPTIMIZER, add_bridges=false),
     kwargs...
 )
@@ -207,22 +209,23 @@ LinMPC controller with a sample time Ts = 4.0 s, OSQP optimizer, KalmanFilter es
 """
 function LinMPC(
     estim::SE;
-    Hp::Union{Int, Nothing} = nothing,
+    Hp::Int = default_Hp(estim.model),
     Hc::Int = DEFAULT_HC,
-    Mwt = fill(DEFAULT_MWT, estim.model.ny),
-    Nwt = fill(DEFAULT_NWT, estim.model.nu),
-    Lwt = fill(DEFAULT_LWT, estim.model.nu),
-    Cwt = DEFAULT_CWT,
-    M_Hp = nothing,
-    N_Hc = nothing,
-    L_Hp = nothing,
+    Mwt  = fill(DEFAULT_MWT, estim.model.ny),
+    Nwt  = fill(DEFAULT_NWT, estim.model.nu),
+    Lwt  = fill(DEFAULT_LWT, estim.model.nu),
+    M_Hp = Diagonal(repeat(Mwt, Hp)),
+    N_Hc = Diagonal(repeat(Nwt, Hc)),
+    L_Hp = Diagonal(repeat(Lwt, Hp)),
+    Cwt  = DEFAULT_CWT,
     optim::JM = JuMP.Model(DEFAULT_LINMPC_OPTIMIZER, add_bridges=false),
 ) where {NT<:Real, SE<:StateEstimator{NT}, JM<:JuMP.GenericModel}
-    isa(estim.model, LinModel) || error("estim.model type must be LinModel") 
-    Hp = default_Hp(estim.model, Hp)
-    isnothing(M_Hp) && (M_Hp = Diagonal{NT}(repeat(Mwt, Hp)))
-    isnothing(N_Hc) && (N_Hc = Diagonal{NT}(repeat(Nwt, Hc)))
-    isnothing(L_Hp) && (L_Hp = Diagonal{NT}(repeat(Lwt, Hp)))
+    isa(estim.model, LinModel) || error("estim.model type must be a LinModel") 
+    nk = estimate_delays(estim.model)
+    if Hp ≤ nk
+        @warn("prediction horizon Hp ($Hp) ≤ estimated number of delays in model "*
+              "($nk), the closed-loop system may be unstable or zero-gain (unresponsive)")
+    end
     return LinMPC{NT, SE, JM}(estim, Hp, Hc, M_Hp, N_Hc, L_Hp, Cwt, optim)
 end