added: tests with Float32 type parameters

Author: franckgaga

src/controller/explicitmpc.jl

@@ -51,10 +51,10 @@ struct ExplicitMPC{NT<:Real, SE<:StateEstimator} <: PredictiveController{NT}
         P̃_chol = cholesky(P̃)
         Ks, Ps = init_stochpred(estim, Hp)
         # dummy vals (updated just before optimization):
-        d0, D̂0, D̂E = zeros(nd), zeros(nd*Hp), zeros(nd + nd*Hp)
+        d0, D̂0, D̂E = zeros(NT, nd), zeros(NT, nd*Hp), zeros(NT, nd + nd*Hp)
         Ŷop, Dop = repeat(model.yop, Hp), repeat(model.dop, Hp)
         nvar = size(Ẽ, 2)
-        ΔŨ = zeros(nvar)
+        ΔŨ = zeros(NT, nvar)
         mpc = new{NT, SE}(
             estim,
             ΔŨ, ŷ,
@@ -167,9 +167,9 @@ function ExplicitMPC(
 ) 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(repeat(Mwt, Hp)))
-    isnothing(N_Hc) && (N_Hc = Diagonal(repeat(Nwt, Hc)))
-    isnothing(L_Hp) && (L_Hp = Diagonal(repeat(Lwt, 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)))
     return ExplicitMPC{NT, SE}(estim, Hp, Hc, M_Hp, N_Hc, L_Hp)
 end
 

src/controller/linmpc.jl

@@ -1,8 +1,14 @@
 const DEFAULT_LINMPC_OPTIMIZER = OSQP.MathOptInterfaceOSQP.Optimizer
 
-struct LinMPC{NT<:Real, SE<:StateEstimator} <: PredictiveController{NT}
+struct LinMPC{
+    NT<:Real, 
+    SE<:StateEstimator, 
+    JM<:JuMP.GenericModel
+} <: PredictiveController{NT}
     estim::SE
-    optim::JuMP.Model
+    # note: `NT` and the number type `JNT` in `JuMP.GenericModel{JNT}` can be
+    # different since solvers that support non-Float64 are scarce.
+    optim::JM
     con::ControllerConstraint{NT}
     ΔŨ::Vector{NT}
     ŷ ::Vector{NT}
@@ -34,9 +40,9 @@ struct LinMPC{NT<:Real, SE<:StateEstimator} <: PredictiveController{NT}
     D̂E::Vector{NT}
     Ŷop::Vector{NT}
     Dop::Vector{NT}
-    function LinMPC{NT, SE}(
-        estim::SE, Hp, Hc, M_Hp, N_Hc, L_Hp, Cwt, optim
-    ) where {NT<:Real, SE<:StateEstimator}
+    function LinMPC{NT, SE, JM}(
+        estim::SE, Hp, Hc, M_Hp, N_Hc, L_Hp, Cwt, optim::JM
+    ) where {NT<:Real, SE<:StateEstimator, JM<:JuMP.GenericModel}
         model = estim.model
         nu, ny, nd = model.nu, model.ny, model.nd
         ŷ = copy(model.yop) # dummy vals (updated just before optimization)
@@ -56,7 +62,7 @@ struct LinMPC{NT<:Real, SE<:StateEstimator} <: PredictiveController{NT}
         Ŷop, Dop = repeat(model.yop, Hp), repeat(model.dop, Hp)
         nvar = size(Ẽ, 2)
         ΔŨ = zeros(NT, nvar)
-        mpc = new{NT, SE}(
+        mpc = new{NT, SE, JM}(
             estim, optim, con,
             ΔŨ, ŷ,
             Hp, Hc, 
@@ -68,7 +74,7 @@ struct LinMPC{NT<:Real, SE<:StateEstimator} <: PredictiveController{NT}
             d0, D̂0, D̂E,
             Ŷop, Dop,
         )
-        init_optimization!(mpc)
+        init_optimization!(mpc, optim)
         return mpc
     end
 end
@@ -166,7 +172,7 @@ function LinMPC(
     M_Hp = nothing,
     N_Hc = nothing,
     L_Hp = nothing,
-    optim::JuMP.Model = JuMP.Model(DEFAULT_LINMPC_OPTIMIZER, add_bridges=false),
+    optim::JuMP.GenericModel = JuMP.Model(DEFAULT_LINMPC_OPTIMIZER, add_bridges=false),
     kwargs...
 )
     estim = SteadyKalmanFilter(model; kwargs...)
@@ -207,24 +213,24 @@ function LinMPC(
     M_Hp = nothing,
     N_Hc = nothing,
     L_Hp = nothing,
-    optim::JuMP.Model = JuMP.Model(DEFAULT_LINMPC_OPTIMIZER, add_bridges=false),
-) where {NT<:Real, SE<:StateEstimator{NT}}
+    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(repeat(Mwt, Hp)))
-    isnothing(N_Hc) && (N_Hc = Diagonal(repeat(Nwt, Hc)))
-    isnothing(L_Hp) && (L_Hp = Diagonal(repeat(Lwt, Hp)))
-    return LinMPC{NT, SE}(estim, Hp, Hc, M_Hp, N_Hc, L_Hp, Cwt, optim)
+    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)))
+    return LinMPC{NT, SE, JM}(estim, Hp, Hc, M_Hp, N_Hc, L_Hp, Cwt, optim)
 end
 
 """
-    init_optimization!(mpc::LinMPC)
+    init_optimization!(mpc::LinMPC, optim::JuMP.GenericModel)
 
 Init the quadratic optimization for [`LinMPC`](@ref) controllers.
 """
-function init_optimization!(mpc::LinMPC)
+function init_optimization!(mpc::LinMPC, optim::JuMP.GenericModel)
     # --- variables and linear constraints ---
-    optim, con = mpc.optim, mpc.con
+    con = mpc.con
     nvar = length(mpc.ΔŨ)
     set_silent(optim)
     limit_solve_time(mpc)

src/controller/nonlinmpc.jl

@@ -1,8 +1,15 @@
 const DEFAULT_NONLINMPC_OPTIMIZER = optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes")
 
-struct NonLinMPC{NT<:Real, SE<:StateEstimator, JEfunc<:Function} <: PredictiveController{NT}
+struct NonLinMPC{
+    NT<:Real, 
+    SE<:StateEstimator, 
+    JM<:JuMP.GenericModel, 
+    JEfunc<:Function
+} <: PredictiveController{NT}
     estim::SE
-    optim::JuMP.Model
+    # note: `NT` and the number type `JNT` in `JuMP.GenericModel{JNT}` can be
+    # different since solvers that support non-Float64 are scarce.
+    optim::JM
     con::ControllerConstraint{NT}
     ΔŨ::Vector{NT}
     ŷ ::Vector{NT}
@@ -35,27 +42,28 @@ struct NonLinMPC{NT<:Real, SE<:StateEstimator, JEfunc<:Function} <: PredictiveCo
     D̂E::Vector{NT}
     Ŷop::Vector{NT}
     Dop::Vector{NT}
-    function NonLinMPC{NT, SE, JEFunc}(
-        estim::SE, Hp, Hc, M_Hp, N_Hc, L_Hp, Cwt, Ewt, JE::JEFunc, optim
-    ) where {NT<:Real, SE<:StateEstimator, JEFunc<:Function}
+    function NonLinMPC{NT, SE, JM, JEFunc}(
+        estim::SE, Hp, Hc, M_Hp, N_Hc, L_Hp, Cwt, Ewt, JE::JEFunc, optim::JM
+    ) where {NT<:Real, SE<:StateEstimator, JM<:JuMP.GenericModel, JEFunc<:Function}
         model = estim.model
         nu, ny, nd = model.nu, model.ny, model.nd
         ŷ = copy(model.yop) # dummy vals (updated just before optimization)
         validate_weights(model, Hp, Hc, M_Hp, N_Hc, L_Hp, Cwt, Ewt)
         M_Hp, N_Hc, L_Hp = float(M_Hp), float(N_Hc), float(L_Hp) # debug julia 1.6
-        R̂y, R̂u = zeros(ny*Hp), zeros(nu*Hp) # dummy vals (updated just before optimization)
+        # dummy vals (updated just before optimization):
+        R̂y, R̂u = zeros(NT, ny*Hp), zeros(NT, nu*Hp)
         noR̂u = iszero(L_Hp)
         S, T = init_ΔUtoU(nu, Hp, Hc)
         E, F, G, J, K, V, ex̂, fx̂, gx̂, jx̂, kx̂, vx̂ = init_predmat(estim, model, Hp, Hc)
         con, S̃, Ñ_Hc, Ẽ = init_defaultcon(estim, Hp, Hc, Cwt, S, N_Hc, E, ex̂, fx̂, gx̂, jx̂, kx̂, vx̂)
         P̃, q̃, p = init_quadprog(model, Ẽ, S̃, M_Hp, Ñ_Hc, L_Hp)
         Ks, Ps = init_stochpred(estim, Hp)
         # dummy vals (updated just before optimization):
-        d0, D̂0, D̂E = zeros(nd), zeros(nd*Hp), zeros(nd + nd*Hp)
+        d0, D̂0, D̂E = zeros(NT, nd), zeros(NT, nd*Hp), zeros(NT, nd + nd*Hp)
         Ŷop, Dop = repeat(model.yop, Hp), repeat(model.dop, Hp)
         nvar = size(Ẽ, 2)
-        ΔŨ = zeros(nvar)
-        mpc = new{NT, SE, JEFunc}(
+        ΔŨ = zeros(NT, nvar)
+        mpc = new{NT, SE, JM, JEFunc}(
             estim, optim, con,
             ΔŨ, ŷ,
             Hp, Hc, 
@@ -67,7 +75,7 @@ struct NonLinMPC{NT<:Real, SE<:StateEstimator, JEfunc<:Function} <: PredictiveCo
             d0, D̂0, D̂E,
             Ŷop, Dop,
         )
-        init_optimization!(mpc)
+        init_optimization!(mpc, optim)
         return mpc
     end
 end
@@ -168,7 +176,7 @@ function NonLinMPC(
     M_Hp = nothing,
     N_Hc = nothing,
     L_Hp = nothing,
-    optim::JuMP.Model = JuMP.Model(DEFAULT_NONLINMPC_OPTIMIZER, add_bridges=false),
+    optim::JuMP.GenericModel = JuMP.Model(DEFAULT_NONLINMPC_OPTIMIZER, add_bridges=false),
     kwargs...
 )
     estim = UnscentedKalmanFilter(model; kwargs...)
@@ -188,7 +196,7 @@ function NonLinMPC(
     M_Hp = nothing,
     N_Hc = nothing,
     L_Hp = nothing,
-    optim::JuMP.Model = JuMP.Model(DEFAULT_NONLINMPC_OPTIMIZER, add_bridges=false),
+    optim::JuMP.GenericModel = JuMP.Model(DEFAULT_NONLINMPC_OPTIMIZER, add_bridges=false),
     kwargs...
 )
     estim = SteadyKalmanFilter(model; kwargs...)
@@ -231,13 +239,13 @@ function NonLinMPC(
     M_Hp = nothing,
     N_Hc = nothing,
     L_Hp = nothing,
-    optim::JuMP.Model = JuMP.Model(DEFAULT_NONLINMPC_OPTIMIZER, add_bridges=false),
-) where {NT<:Real, SE<:StateEstimator{NT}, JEFunc<:Function}
+    optim::JM = JuMP.Model(DEFAULT_NONLINMPC_OPTIMIZER, add_bridges=false),
+) where {NT<:Real, SE<:StateEstimator{NT}, JM<:JuMP.GenericModel, JEFunc<:Function}
     Hp = default_Hp(estim.model, Hp)
-    isnothing(M_Hp) && (M_Hp = Diagonal(repeat(Mwt, Hp)))
-    isnothing(N_Hc) && (N_Hc = Diagonal(repeat(Nwt, Hc)))
-    isnothing(L_Hp) && (L_Hp = Diagonal(repeat(Lwt, Hp)))
-    return NonLinMPC{NT, SE, JEFunc}(estim, Hp, Hc, M_Hp, N_Hc, L_Hp, Cwt, Ewt, JE, optim)
+    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)))
+    return NonLinMPC{NT, SE, JM, JEFunc}(estim, Hp, Hc, M_Hp, N_Hc, L_Hp, Cwt, Ewt, JE, optim)
 end
 
 """
@@ -256,13 +264,13 @@ function addinfo!(info, mpc::NonLinMPC)
 end
 
 """
-    init_optimization!(mpc::NonLinMPC)
+    init_optimization!(mpc::NonLinMPC, optim::JuMP.GenericModel)
 
 Init the nonlinear optimization for [`NonLinMPC`](@ref) controllers.
 """
-function init_optimization!(mpc::NonLinMPC{NT, SE, JEFunc}) where {NT<:Real, SE, JEFunc}
+function init_optimization!(mpc::NonLinMPC, optim::JuMP.GenericModel{JNT}) where JNT<:Real
     # --- variables and linear constraints ---
-    optim, con = mpc.optim, mpc.con
+    con = mpc.con
     nvar = length(mpc.ΔŨ)
     set_silent(optim)
     limit_solve_time(mpc)
@@ -277,10 +285,10 @@ function init_optimization!(mpc::NonLinMPC{NT, SE, JEFunc}) where {NT<:Real, SE,
     # inspired from https://jump.dev/JuMP.jl/stable/tutorials/nonlinear/tips_and_tricks/#User-defined-operators-with-vector-outputs
     Jfunc, gfunc = let mpc=mpc, model=model, ng=ng, nvar=nvar , nŶ=Hp*ny, nx̂=nx̂
         last_ΔŨtup_float, last_ΔŨtup_dual = nothing, nothing
-        Ŷ_cache::DiffCache{Vector{NT}, Vector{NT}} = DiffCache(zeros(NT, nŶ), nvar + 3)
-        g_cache::DiffCache{Vector{NT}, Vector{NT}} = DiffCache(zeros(NT, ng), nvar + 3)
-        x̂_cache::DiffCache{Vector{NT}, Vector{NT}} = DiffCache(zeros(NT, nx̂), nvar + 3)
-        function Jfunc(ΔŨtup::NT...)
+        Ŷ_cache::DiffCache{Vector{JNT}, Vector{JNT}} = DiffCache(zeros(JNT, nŶ), nvar + 3)
+        g_cache::DiffCache{Vector{JNT}, Vector{JNT}} = DiffCache(zeros(JNT, ng), nvar + 3)
+        x̂_cache::DiffCache{Vector{JNT}, Vector{JNT}} = DiffCache(zeros(JNT, nx̂), nvar + 3)
+        function Jfunc(ΔŨtup::JNT...)
             Ŷ = get_tmp(Ŷ_cache, ΔŨtup[1])
             ΔŨ = collect(ΔŨtup)
             if ΔŨtup != last_ΔŨtup_float
@@ -304,7 +312,7 @@ function init_optimization!(mpc::NonLinMPC{NT, SE, JEFunc}) where {NT<:Real, SE,
             end
             return obj_nonlinprog(mpc, model, Ŷ, ΔŨ)
         end
-        function gfunc_i(i, ΔŨtup::NTuple{N, NT}) where N
+        function gfunc_i(i, ΔŨtup::NTuple{N, JNT}) where N
             g = get_tmp(g_cache, ΔŨtup[1])
             if ΔŨtup != last_ΔŨtup_float
                 x̂ = get_tmp(x̂_cache, ΔŨtup[1])

src/estimator/internal_model.jl

@@ -142,9 +142,10 @@ function matrices_internalmodel(model::LinModel)
     return Â, B̂u, Ĉ, B̂d, D̂d
 end
 "Return empty matrices if `model` is not a [`LinModel`](@ref)."
-function matrices_internalmodel(model::SimModel)
+function matrices_internalmodel(model::SimModel{NT}) where NT<:Real
     nu, nx, nd = model.nu, model.nx, model.nd
-    return zeros(0, nx), zeros(0, nu), zeros(0, nx), zeros(0, nd), zeros(0, nd)
+    Â, B̂u, Ĉ, B̂d, D̂d = zeros(NT,0,nx), zeros(NT,0,nu), zeros(NT,0,nx), zeros(NT,0,nd), zeros(NT,0,nd)
+    return Â, B̂u, Ĉ, B̂d, D̂d
 end
 
 @doc raw"""
@@ -216,14 +217,16 @@ The [`InternalModel`](@ref) updates the deterministic `x̂d` and stochastic `x̂
 This estimator does not augment the state vector, thus ``\mathbf{x̂ = x̂_d}``. See 
 [`init_internalmodel`](@ref) for details. 
 """
-function update_estimate!(estim::InternalModel, u, ym, d=empty(estim.x̂))
+function update_estimate!(
+    estim::InternalModel{NT, SM}, u, ym, d=empty(estim.x̂)
+) where {NT<:Real, SM}
     model = estim.model
     x̂d, x̂s = estim.x̂d, estim.x̂s
     # -------------- deterministic model ---------------------
     ŷd = h(model, x̂d, d)
     x̂d[:] = f(model, x̂d, u, d) # this also updates estim.x̂ (they are the same object)
     # --------------- stochastic model -----------------------
-    ŷs = zeros(model.ny)
+    ŷs = zeros(NT, model.ny)
     ŷs[estim.i_ym] = ym - ŷd[estim.i_ym]   # ŷs=0 for unmeasured outputs
     x̂s[:] = estim.Âs*x̂s + estim.B̂s*ŷs
     return x̂d
@@ -247,11 +250,11 @@ of the measured ``\mathbf{ŷ_s^m} = \mathbf{y^m} - \mathbf{ŷ_d^m}`` and unmea
 This estimator does not augment the state vector, thus ``\mathbf{x̂ = x̂_d}``. See
 [`init_internalmodel`](@ref) for details.
 """
-function init_estimate!(estim::InternalModel, model::LinModel, u, ym, d)
+function init_estimate!(estim::InternalModel, model::LinModel{NT}, u, ym, d) where NT<:Real
     x̂d, x̂s = estim.x̂d, estim.x̂s
     x̂d[:] = (I - model.A)\(model.Bu*u + model.Bd*d)
     ŷd = h(model, x̂d, d)
-    ŷs = zeros(model.ny)
+    ŷs = zeros(NT, model.ny)
     ŷs[estim.i_ym] = ym - ŷd[estim.i_ym]  # ŷs=0 for unmeasured outputs
     x̂s[:] = (I-estim.Âs)\estim.B̂s*ŷs
     return nothing

src/estimator/kalman.jl

@@ -540,7 +540,7 @@ noise, respectively.
      Kalman, H∞, and Nonlinear Approaches", John Wiley & Sons, p. 433–459, https://doi.org/10.1002/0470045345.ch14, 
      ISBN9780470045343.
 """
-function update_estimate!(estim::UnscentedKalmanFilter{NT, SM}, u, ym, d) where {NT<:Real, SM}
+function update_estimate!(estim::UnscentedKalmanFilter{NT}, u, ym, d) where NT<:Real
     x̂, P̂, Q̂, R̂, K̂ = estim.x̂, estim.P̂, estim.Q̂, estim.R̂, estim.K̂
     nym, nx̂, nσ = estim.nym, estim.nx̂, estim.nσ
     γ, m̂, Ŝ = estim.γ, estim.m̂, estim.Ŝ

src/estimator/luenberger.jl

@@ -35,8 +35,8 @@ struct Luenberger{NT<:Real, SM<:LinModel} <: StateEstimator{NT}
             error("Cannot compute the Luenberger gain K̂ with specified poles p̂.")
         end
         Ĉm, D̂dm = Ĉ[i_ym, :], D̂d[i_ym, :] # measured outputs ym only
-        lastu0 = zeros(model.nu)
-        x̂ = [zeros(model.nx); zeros(nxs)]
+        lastu0 = zeros(NT, model.nu)
+        x̂ = [zeros(NT, model.nx); zeros(NT, nxs)]
         return new{NT, SM}(
             model, 
             lastu0, x̂,

src/model/nonlinmodel.jl

@@ -14,7 +14,7 @@ struct NonLinModel{NT<:Real, F<:Function, H<:Function} <: SimModel{NT}
         f::F, h::H, Ts, nu, nx, ny, nd
     ) where {NT<:Real, F<:Function, H<:Function}
         Ts > 0 || error("Sampling time Ts must be positive")
-        validate_fcts(f, h)
+        validate_fcts(NT, f, h)
         uop = zeros(NT, nu)
         yop = zeros(NT, ny)
         dop = zeros(NT, nd)
@@ -74,24 +74,18 @@ function NonLinModel(
 end
 
 "Validate `f` and `h` function argument signatures."
-function validate_fcts(f, h)
+function validate_fcts(NT, f, h)
     fargsvalid1 = hasmethod(f,
-        Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}
+        Tuple{Vector{NT}, Vector{NT}, Vector{NT}}
     )
-    fargsvalid2 = hasmethod(f,
-        Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}}
-    )
-    if !fargsvalid1 && !fargsvalid2
-        error("state function has no method of type "*
-            "f(x::Vector{Float64}, u::Vector{Float64}, d::Vector{Float64}) or "*
-            "f(x::Vector{ComplexF64}, u::Vector{Float64}, d::Vector{Float64})")
+    if !fargsvalid1
+        error("state function has no method with type signature "*
+              "f(x::Vector{$(NT)}, u::Vector{$(NT)}, d::Vector{$(NT)})")
     end
-    hargsvalid1 = hasmethod(h,Tuple{Vector{Float64}, Vector{Float64}})
-    hargsvalid2 = hasmethod(h,Tuple{Vector{ComplexF64}, Vector{Float64}})
-    if !hargsvalid1 && !hargsvalid2
-        error("output function has no method of type "*
-            "h(x::Vector{Float64}, d::Vector{Float64}) or "*
-            "h(x::Vector{ComplexF64}, d::Vector{Float64})")
+    hargsvalid1 = hasmethod(h,Tuple{Vector{NT}, Vector{NT}})
+    if !hargsvalid1
+        error("output function has no method with type signature "*
+              "h(x::Vector{$(NT)}, d::Vector{$(NT)})")
     end
 end