Infer state-space f and h function types to eleminate the remaining abstract fields in SimModel object.

Author: franckgaga

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelPredictiveControl"
 uuid = "61f9bdb8-6ae4-484a-811f-bbf86720c31c"
 authors = ["Francis Gagnon"]
-version = "0.3.2"
+version = "0.3.3"
 
 [deps]
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"

example/juMPC.jl

@@ -130,7 +130,7 @@ function test_mpc(model, mpc)
 end
 
 @time u_data, y_data, r_data, d_data = test_mpc(linModel4, mpc)
-@profview u_data, y_data, r_data, d_data = test_mpc(linModel4, mpc)
+#@profview u_data, y_data, r_data, d_data = test_mpc(linModel4, mpc)
 #=
 using PlotThemes, Plots
 #theme(:default)

src/estimator/internal_model.jl

@@ -63,7 +63,7 @@ sum both values : ``\mathbf{ŷ = ŷ_d + ŷ_s}``.
 # Examples
 ```jldoctest
 julia> estim = InternalModel(LinModel([tf(3, [30, 1]); tf(-2, [5, 1])], 0.5), i_ym=[2])
-InternalModel{LinModel} estimator with a sample time Ts = 0.5 s and:
+InternalModel estimator with a sample time Ts = 0.5 s, LinModel and:
  1 manipulated inputs u
  2 states x̂
  1 measured outputs ym

src/estimator/kalman.jl

@@ -99,7 +99,7 @@ ones, for ``\mathbf{Ĉ^u, D̂_d^u}``).
 julia> model = LinModel([tf(3, [30, 1]); tf(-2, [5, 1])], 0.5);
 
 julia> estim = SteadyKalmanFilter(model, i_ym=[2], σR=[1], σQ_int=[0.01])
-SteadyKalmanFilter estimator with a sample time Ts = 0.5 s and:
+SteadyKalmanFilter estimator with a sample time Ts = 0.5 s, LinModel and:
  1 manipulated inputs u
  3 states x̂
  1 measured outputs ym
@@ -244,7 +244,7 @@ with ``\mathbf{P̂}_{-1}(0) = \mathrm{diag}\{ \mathbf{P}(0), \mathbf{P_{int}}(0)
 julia> model = LinModel([tf(3, [30, 1]); tf(-2, [5, 1])], 0.5);
 
 julia> estim = KalmanFilter(model, i_ym=[2], σR=[1], σP0=[100, 100], σQ_int=[0.01])
-KalmanFilter estimator with a sample time Ts = 0.5 s and:
+KalmanFilter estimator with a sample time Ts = 0.5 s, LinModel and:
  1 manipulated inputs u
  3 states x̂
  1 measured outputs ym
@@ -407,7 +407,7 @@ unmeasured ones, for ``\mathbf{ĥ^u}``).
 julia> model = NonLinModel((x,u,_)->0.1x+u, (x,_)->2x, 10, 1, 1, 1);
 
 julia> estim = UnscentedKalmanFilter(model, σR=[1], nint_ym=[2], σP0_int=[1, 1])
-UnscentedKalmanFilter{NonLinModel} estimator with a sample time Ts = 10.0 s and:
+UnscentedKalmanFilter estimator with a sample time Ts = 10.0 s, NonLinModel and:
  1 manipulated inputs u
  3 states x̂
  1 measured outputs ym

src/predictive_control.jl

@@ -416,7 +416,7 @@ default arguments.
 julia> model = NonLinModel((x,u,_)->0.5x+u, (x,_)->2x, 10, 1, 1, 1);
 
 julia> mpc = NonLinMPC(model, Hp=20, Hc=1, Cwt=1e6)
-NonLinMPC controller with a sample time Ts = 10.0 s, UnscentedKalmanFilter{NonLinModel} estimator and:
+NonLinMPC controller with a sample time Ts = 10.0 s, UnscentedKalmanFilter estimator and:
  1 manipulated inputs u
  2 states x̂
  1 measured outputs ym
@@ -444,7 +444,7 @@ julia> model = NonLinModel((x,u,_)->0.5x+u, (x,_)->2x, 10, 1, 1, 1);
 julia> estim = UnscentedKalmanFilter(model, σQ_int=[0.05]);
 
 julia> mpc = NonLinMPC(estim, Hp=20, Hc=1, Cwt=1e6)
-NonLinMPC controller with a sample time Ts = 10.0 s, UnscentedKalmanFilter{NonLinModel} estimator and:
+NonLinMPC controller with a sample time Ts = 10.0 s, UnscentedKalmanFilter estimator and:
  1 manipulated inputs u
  2 states x̂
  1 measured outputs ym
@@ -1232,8 +1232,8 @@ repeatdiag(A, n::Int) = kron(I(n), A)
 
 
 function Base.show(io::IO, mpc::PredictiveController)
-    println(io, "$(typeof(mpc).name.name) controller with a sample time "*
-                "Ts = $(mpc.estim.model.Ts) s, $(typeof(mpc).parameters[1]) estimator and:")
+    println(io, "$(typeof(mpc).name.name) controller with a sample time Ts = "*
+                "$(mpc.estim.model.Ts) s, $(typeof(mpc.estim).name.name) estimator and:")
     println(io, " $(mpc.estim.model.nu) manipulated inputs u")
     println(io, " $(mpc.estim.nx̂) states x̂")
     println(io, " $(mpc.estim.nym) measured outputs ym")

src/sim_model.jl

@@ -1,3 +1,5 @@
+const IntRangeOrVector = Union{UnitRange{Int}, Vector{Int}}
+
 @doc raw"""
 Abstract supertype of [`LinModel`](@ref) and [`NonLinModel`](@ref) types.
 
@@ -18,7 +20,7 @@ julia> y = model()
 """
 abstract type SimModel end
 
-struct LinModel <: SimModel
+struct LinModel{F<:Function, H<:Function} <: SimModel
     A   ::Matrix{Float64}
     Bu  ::Matrix{Float64}
     C   ::Matrix{Float64}
@@ -35,20 +37,15 @@ struct LinModel <: SimModel
     uop::Vector{Float64}
     yop::Vector{Float64}
     dop::Vector{Float64}
-    function LinModel(A, Bu, C, Bd, Dd, Ts, nu, nx, ny, nd)
+    function LinModel{F, H}(
+        A, Bu, C, Bd, Dd, f::F, h::H, Ts, nu, nx, ny, nd
+    ) where{F<:Function, H<:Function}
         size(A)  == (nx,nx) || error("A size must be $((nx,nx))")
         size(Bu) == (nx,nu) || error("Bu size must be $((nx,nu))")
         size(C)  == (ny,nx) || error("C size must be $((ny,nx))")
         size(Bd) == (nx,nd) || error("Bd size must be $((nx,nd))")
         size(Dd) == (ny,nd) || error("Dd size must be $((ny,nd))")
         Ts > 0 || error("Sampling time Ts must be positive")
-        # the `let` block captures and fixes A, Bu, Bd, C, Dd values (faster computations):
-        f = let A=A, Bu=Bu, Bd=Bd 
-            (x,u,d) -> A*x + Bu*u + Bd*d
-        end
-        h = let C=C, Dd=Dd
-            (x,d) -> C*x + Dd*d
-        end
         uop = zeros(nu)
         yop = zeros(ny)
         dop = zeros(nd)
@@ -57,8 +54,12 @@ struct LinModel <: SimModel
     end
 end
 
-
-const IntRangeOrVector = Union{UnitRange{Int}, Vector{Int}}
+"Infer the type of state-space `f` and `h` functions and construct the linear model."
+function LinModel_ssfunc(
+    A, Bu, C, Bd, Dd, f::F, h::H, Ts, nu, nx, ny, nd
+) where {F<:Function, H<:Function}
+    return LinModel{F, H}(A, Bu, C, Bd, Dd, f, h, Ts, nu, nx, ny, nd)
+end
 
 @doc raw"""
     LinModel(sys::StateSpace[, Ts]; i_u=1:size(sys,2), i_d=Int[])
@@ -159,7 +160,14 @@ function LinModel(
     Bd  = sys_dis.B[:,nu+1:end]
     C   = sys_dis.C;
     Dd  = sys_dis.D[:,nu+1:end]
-    return LinModel(A, Bu, C, Bd, Dd, Ts, nu, nx, ny, nd)
+    # the `let` block captures and fixes A, Bu, Bd, C, Dd values (faster computations):
+    f = let A=A, Bu=Bu, Bd=Bd 
+        (x,u,d) -> A*x + Bu*u + Bd*d
+    end
+    h = let C=C, Dd=Dd
+        (x,d) -> C*x + Dd*d
+    end
+    return LinModel_ssfunc(A, Bu, C, Bd, Dd, f, h, Ts, nu, nx, ny, nd)
 end
 
 @doc raw"""
@@ -208,7 +216,6 @@ function LinModel(sys::DelayLtiSystem, Ts::Real; kwargs...)
     return LinModel(sys_dis, Ts; kwargs...)
 end
 
-
 @doc raw"""
     steadystate(model::LinModel, u, d=Float64[])
 
@@ -226,6 +233,28 @@ function steadystate(model::LinModel, u, d=Float64[])
     return pinv(I - model.A)*(model.Bu*(u - model.uop) + model.Bd*(d - model.dop))
 end
 
+struct NonLinModel{F<:Function, H<:Function} <: SimModel
+    x::Vector{Float64}
+    f::F
+    h::H
+    Ts::Float64
+    nu::Int
+    nx::Int
+    ny::Int
+    nd::Int
+    uop::Vector{Float64}
+    yop::Vector{Float64}
+    dop::Vector{Float64}
+    function NonLinModel{F,H}(f::F, h::H, Ts, nu, nx, ny, nd) where {F<:Function,H<:Function}
+        Ts > 0 || error("Sampling time Ts must be positive")
+        validate_fcts(f, h)
+        uop = zeros(nu)
+        yop = zeros(ny)
+        dop = zeros(nd)
+        x = zeros(nx)
+        return new(x, f, h, Ts, nu, nx, ny, nd, uop, yop, dop)
+    end
+end
 
 @doc raw"""
     NonLinModel(f::Function, h::Function, Ts, nu, nx, ny, nd=0)
@@ -260,27 +289,10 @@ Discrete-time nonlinear model with a sample time Ts = 10.0 s and:
  0 measured disturbances d
 ```
 """
-struct NonLinModel <: SimModel
-    x::Vector{Float64}
-    f::Function
-    h::Function
-    Ts::Float64
-    nu::Int
-    nx::Int
-    ny::Int
-    nd::Int
-    uop::Vector{Float64}
-    yop::Vector{Float64}
-    dop::Vector{Float64}
-    function NonLinModel(f, h, Ts::Real, nu::Int, nx::Int, ny::Int, nd::Int = 0)
-        Ts > 0 || error("Sampling time Ts must be positive")
-        validate_fcts(f, h)
-        uop = zeros(nu)
-        yop = zeros(ny)
-        dop = zeros(nd)
-        x = zeros(nx)
-        return new(x, f, h, Ts, nu, nx, ny, nd, uop, yop, dop)
-    end
+function NonLinModel(
+    f::F, h::H, Ts::Real, nu::Int, nx::Int, ny::Int, nd::Int=0
+) where {F<:Function, H<:Function}
+    return NonLinModel{F, H}(f, h, Ts, nu, nx, ny, nd)
 end
 
 function validate_fcts(f, h)

src/state_estim.jl

@@ -34,8 +34,8 @@ end
 
 
 function Base.show(io::IO, estim::StateEstimator)
-    println(io, "$(typeof(estim)) estimator with "*
-                "a sample time Ts = $(estim.model.Ts) s and:")
+    println(io, "$(typeof(estim).name.name) estimator with a sample time "*
+                "Ts = $(estim.model.Ts) s, $(typeof(estim.model).name.name) and:")
     println(io, " $(estim.model.nu) manipulated inputs u")
     println(io, " $(estim.nx̂) states x̂")
     println(io, " $(estim.nym) measured outputs ym")