added: SimModelBuffer object to reduce the allocations

default `d` argument is now `model.buffer.empty`
(restore the error handling when no d is provided)

Author: franckgaga

src/controller/execute.jl

@@ -1,18 +1,15 @@
 @doc raw"""
-    initstate!(mpc::PredictiveController, u, ym, d=mpc.estim.model.dop) -> x̂
+    initstate!(mpc::PredictiveController, u, ym, d=[]) -> x̂
 
 Init the states of `mpc.estim` [`StateEstimator`](@ref) and warm start `mpc.ΔŨ` at zero.
 """
-function initstate!(mpc::PredictiveController, u, ym, d=mpc.estim.model.dop)
+function initstate!(mpc::PredictiveController, u, ym, d=mpc.estim.model.buffer.empty)
     mpc.ΔŨ .= 0
     return initstate!(mpc.estim, u, ym, d)
 end
 
 @doc raw"""
-    moveinput!(
-        mpc::PredictiveController, ry=mpc.estim.model.yop, d=mpc.estim.model.dop; 
-        <keyword arguments>
-    ) -> u
+    moveinput!(mpc::PredictiveController, ry=mpc.estim.model.yop, d=[]; <keyword args>) -> u
 
 Compute the optimal manipulated input value `u` for the current control period.
 
@@ -32,7 +29,7 @@ See also [`LinMPC`](@ref), [`ExplicitMPC`](@ref), [`NonLinMPC`](@ref).
 
 - `mpc::PredictiveController` : solve optimization problem of `mpc`.
 - `ry=mpc.estim.model.yop` : current output setpoints ``\mathbf{r_y}(k)``.
-- `d=mpc.estim.model.dop` : current measured disturbances ``\mathbf{d}(k)``.
+- `d=[]` : current measured disturbances ``\mathbf{d}(k)``.
 - `D̂=repeat(d, mpc.Hp)` or *`Dhat`* : predicted measured disturbances ``\mathbf{D̂}``, constant
    in the future by default or ``\mathbf{d̂}(k+j)=\mathbf{d}(k)`` for ``j=1`` to ``H_p``.
 - `R̂y=repeat(ry, mpc.Hp)` or *`Rhaty`* : predicted output setpoints ``\mathbf{R̂_y}``, constant
@@ -54,7 +51,7 @@ julia> ry = [5]; u = moveinput!(mpc, ry); round.(u, digits=3)
 function moveinput!(
     mpc::PredictiveController, 
     ry::Vector = mpc.estim.model.yop, 
-    d ::Vector = mpc.estim.model.dop;
+    d ::Vector = mpc.estim.model.buffer.empty;
     Dhat ::Vector = repeat(d,  mpc.Hp),
     Rhaty::Vector = repeat(ry, mpc.Hp),
     Rhatu::Vector = mpc.Uop,
@@ -508,11 +505,11 @@ end
 set_objective_linear_coef!(::PredictiveController, _ ) = nothing
 
 """
-    updatestate!(mpc::PredictiveController, u, ym, d=mpc.estim.model.dop) -> x̂
+    updatestate!(mpc::PredictiveController, u, ym, d=[]) -> x̂
 
 Call [`updatestate!`](@ref) on `mpc.estim` [`StateEstimator`](@ref).
 """
-function updatestate!(mpc::PredictiveController, u, ym, d=mpc.estim.model.dop)
+function updatestate!(mpc::PredictiveController, u, ym, d=mpc.estim.model.buffer.empty)
     return updatestate!(mpc.estim, u, ym, d)
 end
 updatestate!(::PredictiveController, _ ) = throw(ArgumentError("missing measured outputs ym"))

src/estimator/execute.jl

@@ -79,7 +79,7 @@ end
 
 
 @doc raw"""
-    initstate!(estim::StateEstimator, u, ym, d=estim.model.dop) -> x̂
+    initstate!(estim::StateEstimator, u, ym, d=[]) -> x̂
 
 Init `estim.x̂0` states from current inputs `u`, measured outputs `ym` and disturbances `d`.
 
@@ -111,7 +111,7 @@ julia> evaloutput(estim) ≈ y
 true
 ```
 """
-function initstate!(estim::StateEstimator, u, ym, d=estim.model.dop)
+function initstate!(estim::StateEstimator, u, ym, d=estim.model.buffer.empty)
     # --- validate arguments ---
     validate_args(estim, u, ym, d)
     # --- init state estimate ----
@@ -161,7 +161,7 @@ Left `estim.x̂0` estimate unchanged if `model` is not a [`LinModel`](@ref).
 init_estimate!(::StateEstimator, ::SimModel, _ , _ , _ ) = nothing
 
 @doc raw"""
-    evaloutput(estim::StateEstimator, d=estim.model.dop) -> ŷ
+    evaloutput(estim::StateEstimator, d=[]) -> ŷ
 
 Evaluate `StateEstimator` outputs `ŷ` from `estim.x̂0` states and disturbances `d`.
 
@@ -176,7 +176,7 @@ julia> ŷ = evaloutput(kf)
  20.0
 ```
 """
-function evaloutput(estim::StateEstimator{NT}, d=estim.model.dop) where NT <: Real
+function evaloutput(estim::StateEstimator{NT}, d=estim.model.buffer.empty) where NT <: Real
     validate_args(estim.model, d)
     ŷ0 = Vector{NT}(undef, estim.model.ny)
     d0 = d - estim.model.dop
@@ -187,10 +187,10 @@ function evaloutput(estim::StateEstimator{NT}, d=estim.model.dop) where NT <: Re
 end
 
 "Functor allowing callable `StateEstimator` object as an alias for `evaloutput`."
-(estim::StateEstimator)(d=estim.model.dop) = evaloutput(estim, d)
+(estim::StateEstimator)(d=estim.model.buffer.empty) = evaloutput(estim, d)
 
 @doc raw"""
-    updatestate!(estim::StateEstimator, u, ym, d=estim.model.dop) -> x̂
+    updatestate!(estim::StateEstimator, u, ym, d=[]) -> x̂
 
 Update `estim.x̂0` estimate with current inputs `u`, measured outputs `ym` and dist. `d`. 
 
@@ -207,7 +207,7 @@ julia> x̂ = updatestate!(kf, [1], [0]) # x̂[2] is the integrator state (nint_y
  0.0
 ```
 """
-function updatestate!(estim::StateEstimator, u, ym, d=estim.model.dop)
+function updatestate!(estim::StateEstimator, u, ym, d=[])
     validate_args(estim, u, ym, d)
     u0, ym0, d0 = remove_op!(estim, u, ym, d)
     x̂0next = update_estimate!(estim, u0, ym0, d0)

src/model/linmodel.jl

@@ -19,6 +19,7 @@ struct LinModel{NT<:Real} <: SimModel{NT}
     yname::Vector{String}
     dname::Vector{String}
     xname::Vector{String}
+    buffer::SimModelBuffer{NT}
     function LinModel{NT}(A, Bu, C, Bd, Dd, Ts) where {NT<:Real}
         A, Bu = to_mat(A, 1, 1), to_mat(Bu, 1, 1)
         nu, nx = size(Bu, 2), size(A, 2)
@@ -44,13 +45,15 @@ struct LinModel{NT<:Real} <: SimModel{NT}
         dname = ["\$d_{$i}\$" for i in 1:nd]
         xname = ["\$x_{$i}\$" for i in 1:nx]
         x0  = zeros(NT, nx)
+        buffer = SimModelBuffer{NT}(nx, nu, ny, nd)
         return new{NT}(
             A, Bu, C, Bd, Dd, 
             x0, 
             Ts, 
             nu, nx, ny, nd, 
             uop, yop, dop, xop, fop,
-            uname, yname, dname, xname
+            uname, yname, dname, xname,
+            buffer
         )
     end
 end

src/model/nonlinmodel.jl

@@ -17,6 +17,7 @@ struct NonLinModel{NT<:Real, F<:Function, H<:Function, DS<:DiffSolver} <: SimMod
     yname::Vector{String}
     dname::Vector{String}
     xname::Vector{String}
+    buffer::SimModelBuffer{NT}
     function NonLinModel{NT, F, H, DS}(
         f!::F, h!::H, solver::DS, Ts, nu, nx, ny, nd
     ) where {NT<:Real, F<:Function, H<:Function, DS<:DiffSolver}
@@ -31,13 +32,15 @@ struct NonLinModel{NT<:Real, F<:Function, H<:Function, DS<:DiffSolver} <: SimMod
         dname = ["\$d_{$i}\$" for i in 1:nd]
         xname = ["\$x_{$i}\$" for i in 1:nx]
         x0  = zeros(NT, nx)
+        buffer = SimModelBuffer{NT}(nx, nu, ny, nd)
         return new{NT, F, H, DS}(
             x0, 
             f!, h!, 
             solver, Ts, 
             nu, nx, ny, nd, 
             uop, yop, dop, xop, fop,
-            uname, yname, dname, xname
+            uname, yname, dname, xname,
+            buffer
         )
     end
 end

src/predictive_control.jl

@@ -3,7 +3,7 @@ Abstract supertype of all predictive controllers.
 
 ---
 
-    (mpc::PredictiveController)(ry, d=mpc.estim.model.dop; kwargs...) -> u
+    (mpc::PredictiveController)(ry, d=[]; kwargs...) -> u
 
 Functor allowing callable `PredictiveController` object as an alias for [`moveinput!`](@ref).
 
@@ -42,7 +42,7 @@ end
 "Functor allowing callable `PredictiveController` object as an alias for `moveinput!`."
 function (mpc::PredictiveController)(
     ry::Vector = mpc.estim.model.yop, 
-    d ::Vector = mpc.estim.model.dop;
+    d ::Vector = mpc.estim.model.buffer.empty;
     kwargs...
 )
     return moveinput!(mpc, ry, d; kwargs...)

src/sim_model.jl

@@ -5,7 +5,7 @@ Abstract supertype of [`LinModel`](@ref) and [`NonLinModel`](@ref) types.
 
 ---
 
-    (model::SimModel)(d=model.dop) -> y
+    (model::SimModel)(d=[]) -> y
 
 Functor allowing callable `SimModel` object as an alias for [`evaloutput`](@ref).
 
@@ -20,6 +20,29 @@ julia> y = model()
 """
 abstract type SimModel{NT<:Real} end
 
+struct SimModelBuffer{NT<:Real}
+    u::Vector{NT}
+    x::Vector{NT}
+    y::Vector{NT}
+    d::Vector{NT}
+    empty::Vector{NT}
+end
+
+@doc raw"""
+    SimModelBuffer(nu, nx, ny, nd) -> SimModelBuffer{NT}
+
+Create a buffer for `SimModel` objects for inputs, states, outputs, and disturbances.
+"""
+function SimModelBuffer{NT}(nu, nx, ny, nd) where NT <: Real
+    u = Vector{NT}(undef, nu)
+    x = Vector{NT}(undef, nx)
+    y = Vector{NT}(undef, ny)
+    d = Vector{NT}(undef, nd)
+    empty = Vector{NT}(undef, 0)
+    return SimModelBuffer{NT}(u, x, y, d, empty)
+end
+
+
 @doc raw"""
     setop!(model; uop=nothing, yop=nothing, dop=nothing, xop=nothing, fop=nothing) -> model
 
@@ -160,7 +183,7 @@ end
 detailstr(model::SimModel) = ""
 
 @doc raw"""
-    initstate!(model::SimModel, u, d=model.dop) -> x
+    initstate!(model::SimModel, u, d=[]) -> x
 
 Init `model.x0` with manipulated inputs `u` and measured disturbances `d` steady-state.
 
@@ -182,7 +205,7 @@ julia> x ≈ updatestate!(model, u)
 true
 ```
 """
-function initstate!(model::SimModel, u, d=model.dop)
+function initstate!(model::SimModel, u, d=model.buffer.empty)
     validate_args(model::SimModel, d, u)
     u0, d0 = u - model.uop, d - model.dop
     steadystate!(model, u0, d0)
@@ -191,7 +214,7 @@ function initstate!(model::SimModel, u, d=model.dop)
 end
 
 """
-    updatestate!(model::SimModel, u, d=model.dop) -> x
+    updatestate!(model::SimModel, u, d=[]) -> x
 
 Update `model.x0` states with current inputs `u` and measured disturbances `d`.
 
@@ -204,7 +227,7 @@ julia> x = updatestate!(model, [1])
  1.0
 ```
 """
-function updatestate!(model::SimModel{NT}, u, d=model.dop) where NT <: Real
+function updatestate!(model::SimModel{NT}, u, d=model.buffer.empty) where NT <: Real
     validate_args(model::SimModel, d, u)
     xnext0 = Vector{NT}(undef, model.nx)
     u0, d0 = u - model.uop, d - model.dop
@@ -217,7 +240,7 @@ function updatestate!(model::SimModel{NT}, u, d=model.dop) where NT <: Real
 end
 
 """
-    evaloutput(model::SimModel, d=model.dop) -> y
+    evaloutput(model::SimModel, d=[]) -> y
 
 Evaluate `SimModel` outputs `y` from `model.x0` states and measured disturbances `d`.
 
@@ -232,7 +255,7 @@ julia> y = evaloutput(model)
  20.0
 ```
 """
-function evaloutput(model::SimModel{NT}, d=model.dop) where NT <: Real
+function evaloutput(model::SimModel{NT}, d=model.buffer.empty) where NT <: Real
     validate_args(model, d)
     y0 = Vector{NT}(undef, model.ny)
     d0 = d - model.dop
@@ -262,7 +285,7 @@ to_mat(A::Real, dims...) = fill(A, dims)
 
 
 "Functor allowing callable `SimModel` object as an alias for `evaloutput`."
-(model::SimModel)(d=model.dop) = evaloutput(model::SimModel, d)
+(model::SimModel)(d=model.buffer.empty) = evaloutput(model::SimModel, d)
 
 include("model/linmodel.jl")
 include("model/solver.jl")

src/state_estim.jl

@@ -3,7 +3,7 @@ Abstract supertype of all state estimators.
 
 ---
 
-    (estim::StateEstimator)(d=estim.model.dop) -> ŷ
+    (estim::StateEstimator)(d=[]) -> ŷ
 
 Functor allowing callable `StateEstimator` object as an alias for [`evaloutput`](@ref).