changed: store P̃Δu, P̃u and Tu as sparse matrices

Author: franckgaga

src/controller/construct.jl

@@ -675,7 +675,7 @@ constraints:
 in which ``\mathbf{U_{min}}`` and ``\mathbf{U_{max}}`` vectors respectively contains
 ``\mathbf{u_{min}}`` and ``\mathbf{u_{max}}`` repeated ``H_p`` times.
 """
-function relaxU(Pu::Matrix{NT}, C_umin, C_umax, nϵ) where NT<:Real
+function relaxU(Pu::AbstractMatrix{NT}, C_umin, C_umax, nϵ) where NT<:Real
     if nϵ == 1 # Z̃ = [Z; ϵ]
         # ϵ impacts Z → U conversion for constraint calculations:
         A_Umin, A_Umax = -[Pu  C_umin], [Pu -C_umax] 
@@ -717,7 +717,7 @@ bound, which is more precise than a linear inequality constraint. However, it is
 convenient to treat it as a linear inequality constraint since the optimizer `OSQP.jl` does
 not support pure bounds on the decision variables.
 """
-function relaxΔU(PΔu::Matrix{NT}, C_Δumin, C_Δumax, ΔUmin, ΔUmax, nϵ) where NT<:Real
+function relaxΔU(PΔu::AbstractMatrix{NT}, C_Δumin, C_Δumax, ΔUmin, ΔUmax, nϵ) where NT<:Real
     nZ = size(PΔu, 2)
     if nϵ == 1 # Z̃ = [Z; ϵ]
         ΔŨmin, ΔŨmax = [ΔUmin; NT[0.0]], [ΔUmax; NT[Inf]] # 0 ≤ ϵ ≤ ∞
@@ -754,7 +754,7 @@ Denoting the decision variables augmented with the slack variable
 in which ``\mathbf{Y_{min}, Y_{max}}`` and ``\mathbf{Y_{op}}`` vectors respectively contains
 ``\mathbf{y_{min}, y_{max}}`` and ``\mathbf{y_{op}}`` repeated ``H_p`` times.
 """
-function relaxŶ(E::Matrix{NT}, C_ymin, C_ymax, nϵ) where NT<:Real
+function relaxŶ(E::AbstractMatrix{NT}, C_ymin, C_ymax, nϵ) where NT<:Real
     if nϵ == 1 # Z̃ = [Z; ϵ]
         if iszero(size(E, 1))
             # model is not a LinModel, thus Ŷ constraints are not linear:
@@ -792,7 +792,7 @@ the inequality constraints:
 \end{bmatrix}
 ```
 """
-function relaxterminal(ex̂::Matrix{NT}, c_x̂min, c_x̂max, nϵ) where {NT<:Real}
+function relaxterminal(ex̂::AbstractMatrix{NT}, c_x̂min, c_x̂max, nϵ) where {NT<:Real}
     if nϵ == 1 # Z̃ = [Z; ϵ]
         if iszero(size(ex̂, 1))
             # model is not a LinModel and transcription is a SingleShooting, thus terminal
@@ -821,7 +821,7 @@ It returns the ``\mathbf{Ẽŝ}`` matrix that appears in the defect equation
 \mathbf{A_ŝ Z̃} = - \mathbf{F_ŝ}
 ```
 """
-function augmentdefect(Eŝ::Matrix{NT}, nϵ) where NT<:Real
+function augmentdefect(Eŝ::AbstractMatrix{NT}, nϵ) where NT<:Real
     if nϵ == 1 # Z̃ = [Z; ϵ]
         Ẽŝ = [Eŝ zeros(NT, size(Eŝ, 1), 1)]
     else # Z̃ = Z (only hard constraints)

src/controller/explicitmpc.jl

@@ -15,9 +15,9 @@ struct ExplicitMPC{
     R̂u::Vector{NT}
     R̂y::Vector{NT}
     lastu0::Vector{NT}
-    P̃Δu::Matrix{NT}
-    P̃u ::Matrix{NT} 
-    Tu ::Matrix{NT}
+    P̃Δu::SparseMatrixCSC{NT, Int}
+    P̃u ::SparseMatrixCSC{NT, Int}
+    Tu ::SparseMatrixCSC{NT, Int}
     Tu_lastu0::Vector{NT}
     Ẽ::Matrix{NT}
     F::Vector{NT}

src/controller/linmpc.jl

@@ -24,9 +24,9 @@ struct LinMPC{
     R̂u::Vector{NT}
     R̂y::Vector{NT}
     lastu0::Vector{NT}
-    P̃Δu::Matrix{NT}
-    P̃u ::Matrix{NT} 
-    Tu ::Matrix{NT}
+    P̃Δu::SparseMatrixCSC{NT, Int}
+    P̃u ::SparseMatrixCSC{NT, Int}
+    Tu ::SparseMatrixCSC{NT, Int}
     Tu_lastu0::Vector{NT}
     Ẽ::Matrix{NT}
     F::Vector{NT}

src/controller/nonlinmpc.jl

@@ -40,9 +40,9 @@ struct NonLinMPC{
     R̂u::Vector{NT}
     R̂y::Vector{NT}
     lastu0::Vector{NT}
-    P̃Δu::Matrix{NT}
-    P̃u ::Matrix{NT}
-    Tu ::Matrix{NT}
+    P̃Δu::SparseMatrixCSC{NT, Int}
+    P̃u ::SparseMatrixCSC{NT, Int}
+    Tu ::SparseMatrixCSC{NT, Int}
     Tu_lastu0::Vector{NT}
     Ẽ::Matrix{NT}
     F::Vector{NT}

src/controller/transcription.jl

@@ -81,21 +81,22 @@ in which ``\mathbf{P_{Δu}}`` is defined in the Extended Help section.
     - ``\mathbf{P_{Δu}} = \mathbf{I}`` if `transcription` is a [`SingleShooting`](@ref)
     - ``\mathbf{P_{Δu}} = [\begin{smallmatrix}\mathbf{I} & \mathbf{0} \end{smallmatrix}]``
       if `transcription` is a [`MultipleShooting`](@ref)
+    The matrix is store as as `SparseMatrixCSC` to support both cases efficiently.
 """
 function init_ZtoΔU end
 
 function init_ZtoΔU(
     estim::StateEstimator{NT}, ::SingleShooting, _ , Hc
 ) where {NT<:Real}
-    PΔu = Matrix{NT}(I, estim.model.nu*Hc, estim.model.nu*Hc)
+    PΔu = sparse(Matrix{NT}(I, estim.model.nu*Hc, estim.model.nu*Hc))
     return PΔu
 end
 
 function init_ZtoΔU(
     estim::StateEstimator{NT}, ::MultipleShooting, Hp, Hc
 ) where {NT<:Real}
-    I_nu_Hc = Matrix{NT}(I, estim.model.nu*Hc, estim.model.nu*Hc)
-    PΔu = [I_nu_Hc zeros(NT, estim.model.nu*Hc, estim.nx̂*Hp)]
+    I_nu_Hc = sparse(Matrix{NT}(I, estim.model.nu*Hc, estim.model.nu*Hc))
+    PΔu = [I_nu_Hc spzeros(NT, estim.model.nu*Hc, estim.nx̂*Hp)]
     return PΔu
 end
 
@@ -144,30 +145,33 @@ The ``\mathbf{P_u}`` and ``\mathbf{T_u}`` matrices are defined in the Extended H
     - ``\mathbf{P_u} = \mathbf{P_u^†}`` if `transcription` is a [`SingleShooting`](@ref)
     - ``\mathbf{P_u} = [\begin{smallmatrix}\mathbf{P_u^†} & \mathbf{0} \end{smallmatrix}]``
       if `transcription` is a [`MultipleShooting`](@ref)
+    The conversion matrices are stored as `SparseMatrixCSC` since it was benchmarked that it
+    is generally more performant than normal dense matrices, even for small `nu`, `Hp` and 
+    `Hc` values. Using `Bool` element type and `BitMatrix` is also slower.
 """
 function init_ZtoU(
     estim::StateEstimator{NT}, transcription::TranscriptionMethod, Hp, Hc, nb
 ) where {NT<:Real}
     nu = estim.model.nu
-    # Pu and Tu are `Matrix{NT}`, conversion is faster than `Matrix{Bool}` or `BitMatrix`
-    I_nu = Matrix{NT}(I, nu, nu)
+    I_nu = sparse(Matrix{NT}(I, nu, nu))
     PuDagger = Matrix{NT}(undef, nu*Hp, nu*Hc)
     for i=1:Hc
         ni    = nb[i]
         Q_ni  = repeat(I_nu, ni, 1)
         iRows = (1:nu*ni) .+ @views nu*sum(nb[1:i-1])
-        PuDagger[iRows, :] = [repeat(Q_ni, 1, i) zeros(nu*ni, nu*(Hc-i))]
+        PuDagger[iRows, :] = [repeat(Q_ni, 1, i) spzeros(nu*ni, nu*(Hc-i))]
     end
+    PuDagger = sparse(PuDagger)
     Pu = init_PUmat(estim, transcription, Hp, Hc, PuDagger)
     Tu = repeat(I_nu, Hp)
     return Pu, Tu
 end
 
-
-
 init_PUmat( _ , ::SingleShooting, _ , _ , PuDagger) = PuDagger
-function init_PUmat(estim, ::MultipleShooting, Hp, _ , PuDagger)
-    return [PuDagger zeros(eltype(PuDagger), estim.model.nu*Hp, estim.nx̂*Hp)]
+function init_PUmat(
+    estim, ::MultipleShooting, Hp, _ , PuDagger::AbstractMatrix{NT}
+) where NT<:Real
+    return [PuDagger spzeros(NT, estim.model.nu*Hp, estim.nx̂*Hp)]
 end
 
 @doc raw"""