changed: support multiple shooting in the input conversion matrices

Author: franckgaga

src/controller/construct.jl

@@ -447,20 +447,20 @@ end
 
 
 @doc raw"""
-    init_ΔUtoU(model, Hp, Hc) -> S, T
+    init_ZtoU(estim, Hp, Hc, transcription) -> S, T
 
-Init manipulated input increments to inputs conversion matrices.
+Init decision variables to inputs over ``H_p`` conversion matrices.
 
 The conversion from the input increments ``\mathbf{ΔU}`` to manipulated inputs over ``H_p`` 
 are calculated by:
 ```math
-\mathbf{U} = \mathbf{S} \mathbf{ΔU} + \mathbf{T} \mathbf{u}(k-1)
+\mathbf{U} = \mathbf{S} \mathbf{Z} + \mathbf{T} \mathbf{u}(k-1)
 ```
 The ``\mathbf{S}`` and ``\mathbf{T}`` matrices are defined in the Extended Help section.
 
 # Extended Help
 !!! details "Extended Help"
-    The ``\mathbf{U}`` vector and the two conversion matrices are defined as:
+    The ``\mathbf{U}`` vector and the conversion matrices are defined as:
     ```math
     \mathbf{U} = \begin{bmatrix}
         \mathbf{u}(k + 0)                                           \\
@@ -469,7 +469,7 @@ The ``\mathbf{S}`` and ``\mathbf{T}`` matrices are defined in the Extended Help
         \mathbf{u}(k + H_c - 1)                                     \\
         \vdots                                                      \\
         \mathbf{u}(k + H_p - 1)                                     \end{bmatrix} , \quad
-    \mathbf{S} = \begin{bmatrix}
+    \mathbf{S^†} = \begin{bmatrix}
         \mathbf{I}  & \mathbf{0}    & \cdots    & \mathbf{0}        \\
         \mathbf{I}  & \mathbf{I}    & \cdots    & \mathbf{0}        \\
         \vdots      & \vdots        & \ddots    & \vdots            \\
@@ -484,12 +484,24 @@ The ``\mathbf{S}`` and ``\mathbf{T}`` matrices are defined in the Extended Help
         \vdots                                                      \\
         \mathbf{I}                                                  \end{bmatrix}
     ```
+    and, depending on the transcription method, we have:
+    - ``\mathbf{S} = \mathbf{S^†}`` if `transcription == :singleshooting`
+    - ``\mathbf{S} = [\begin{smallmatrix}\mathbf{S^†} \mathbf{0} \end{smallmatrix}]`` if 
+      `transcription == :multipleshooting` 
 """
-function init_ΔUtoU(model::SimModel{NT}, Hp, Hc) where {NT<:Real}
+function init_ZtoU(estim::StateEstimator{NT}, Hp, Hc, transcription) where {NT<:Real}
+    model = estim.model
     # S and T are `Matrix{NT}` since conversion is faster than `Matrix{Bool}` or `BitMatrix`
     I_nu = Matrix{NT}(I, model.nu, model.nu)
     S_Hc = LowerTriangular(repeat(I_nu, Hc, Hc))
-    S = [S_Hc; repeat(I_nu, Hp - Hc, Hc)]
+    if transcription == :singleshooting
+        O = zeros(NT, model.nu*Hp, 0)
+    elseif transcription == :multipleshooting
+        O = zeros(NT, model.nu*Hp, estim.nx̂*Hp)
+    else
+        throw(ArgumentError("transcription method not recognized"))
+    end
+    S = hcat([S_Hc; repeat(I_nu, Hp - Hc, Hc)], O)
     T = repeat(I_nu, Hp)
     return S, T
 end

src/controller/nonlinmpc.jl

@@ -9,6 +9,7 @@ struct NonLinMPC{
     GCfunc<:Function
 } <: PredictiveController{NT}
     estim::SE
+    transcription::Symbol
     # note: `NT` and the number type `JNT` in `JuMP.GenericModel{JNT}` can be
     # different since solvers that support non-Float64 are scarce.
     optim::JM
@@ -47,7 +48,8 @@ struct NonLinMPC{
     buffer::PredictiveControllerBuffer{NT}
     function NonLinMPC{NT}(
         estim::SE, 
-        Hp, Hc, M_Hp, N_Hc, L_Hp, Cwt, Ewt, JE::JEfunc, gc!::GCfunc, nc, p::P, optim::JM
+        Hp, Hc, M_Hp, N_Hc, L_Hp, Cwt, Ewt, JE::JEfunc, gc!::GCfunc, nc, p::P, 
+        transcription, optim::JM
     ) where {
             NT<:Real, 
             SE<:StateEstimator, 
@@ -62,7 +64,7 @@ struct NonLinMPC{
         weights = ControllerWeights{NT}(model, Hp, Hc, M_Hp, N_Hc, L_Hp, Cwt, Ewt)
         # dummy vals (updated just before optimization):
         R̂y, R̂u, T_lastu = zeros(NT, ny*Hp), zeros(NT, nu*Hp), zeros(NT, nu*Hp)
-        S, T = init_ΔUtoU(model, Hp, Hc)
+        S, T = init_ZtoU(estim, Hp, Hc, transcription)
         E, G, J, K, V, B, ex̂, gx̂, jx̂, kx̂, vx̂, bx̂ = init_predmat(estim, model, Hp, Hc)
         # dummy vals (updated just before optimization):
         F, fx̂  = zeros(NT, ny*Hp), zeros(NT, nx̂)
@@ -81,7 +83,7 @@ struct NonLinMPC{
         ΔŨ = zeros(NT, nΔŨ)
         buffer = PredictiveControllerBuffer{NT}(nu, ny, nd, Hp, Hc, nϵ)
         mpc = new{NT, SE, JM, P, JEfunc, GCfunc}(
-            estim, optim, con,
+            estim, transcription, optim, con,
             ΔŨ, ŷ,
             Hp, Hc, nϵ,
             weights,
@@ -122,9 +124,8 @@ subject to [`setconstraint!`](@ref) bounds, and the custom inequality constraint
 ```
 and the slack ``ϵ`` and the ``\mathbf{Z}`` vector as the decision variables:
 
-- ``\mathbf{Z} = \mathbf{ΔŨ}`` if the transcription method is `:singleshooting`
-- ``\mathbf{Z} = [\begin{smallmatrix} \mathbf{ΔŨ}' & \mathbf{X̂}' \end{smallmatrix}]'`` if the transcription method is `:multipleshoting`
-- ``\mathbf{Z} = [\begin{smallmatrix} \mathbf{ΔŨ}' & \mathbf{X̂}' & \mathbf{Ĉ}' \end{smallmatrix}]'`` if the transcription method is `:directcollocation`
+- ``\mathbf{Z} = \mathbf{ΔU}`` if the transcription method is `:singleshooting`
+- ``\mathbf{Z} = [\begin{smallmatrix} \mathbf{ΔU}' & \mathbf{X̂}' \end{smallmatrix}]'`` if the transcription method is `:multipleshoting`
 
 The economic function ``J_E`` can penalizes solutions with high economic costs. Setting all
 the weights to 0 except ``E``  creates a pure economic model predictive controller (EMPC). 
@@ -176,7 +177,7 @@ This controller allocates memory at each time step for the optimization.
 - `nc=0` : number of custom inequality constraints.
 - `p=model.p` : ``J_E`` and ``\mathbf{g_c}`` functions parameter ``\mathbf{p}`` (any type).
 - `transcription=:singleshooting` : transcription method for the nonlinear optimization 
-    problem, can be `:singleshooting`, `:multipleshoting` or `:directcollocation`.
+    problem, can be `:singleshooting` or `:multipleshoting`.
 - `optim=JuMP.Model(Ipopt.Optimizer)` : nonlinear optimizer used in the predictive
    controller, provided as a [`JuMP.Model`](https://jump.dev/JuMP.jl/stable/api/JuMP/#JuMP.Model)
    (default to [`Ipopt`](https://github.com/jump-dev/Ipopt.jl) optimizer).
@@ -254,13 +255,15 @@ function NonLinMPC(
     gc ::Function = gc!,
     nc::Int = 0,
     p = model.p,
+    transcription = :singleshooting,
     optim::JuMP.GenericModel = JuMP.Model(DEFAULT_NONLINMPC_OPTIMIZER, add_bridges=false),
     kwargs...
 )
     estim = UnscentedKalmanFilter(model; kwargs...)
     return NonLinMPC(
         estim; 
-        Hp, Hc, Mwt, Nwt, Lwt, Cwt, Ewt, JE, gc, nc, p, M_Hp, N_Hc, L_Hp, optim
+        Hp, Hc, Mwt, Nwt, Lwt, Cwt, Ewt, JE, gc, nc, p, M_Hp, N_Hc, L_Hp, 
+        transcription, optim
     )
 end
 
@@ -281,13 +284,15 @@ function NonLinMPC(
     gc ::Function = gc!,
     nc::Int = 0,
     p = model.p,
+    transcription = :singleshooting,
     optim::JuMP.GenericModel = JuMP.Model(DEFAULT_NONLINMPC_OPTIMIZER, add_bridges=false),
     kwargs...
 )
     estim = SteadyKalmanFilter(model; kwargs...)
     return NonLinMPC(
         estim; 
-        Hp, Hc, Mwt, Nwt, Lwt, Cwt, Ewt, JE, gc, nc, p, M_Hp, N_Hc, L_Hp, optim
+        Hp, Hc, Mwt, Nwt, Lwt, Cwt, Ewt, JE, gc, nc, p, M_Hp, N_Hc, L_Hp, 
+        transcription, optim
     )
 end
 
@@ -332,6 +337,7 @@ function NonLinMPC(
     gc ::Function = gc!,
     nc = 0,
     p::P = estim.model.p,
+    transcription = :singleshooting,
     optim::JM = JuMP.Model(DEFAULT_NONLINMPC_OPTIMIZER, add_bridges=false),
 ) where {
     NT<:Real, 
@@ -347,7 +353,8 @@ function NonLinMPC(
     validate_JE(NT, JE)
     gc! = get_mutating_gc(NT, gc)
     return NonLinMPC{NT}(
-        estim, Hp, Hc, M_Hp, N_Hc, L_Hp, Cwt, Ewt, JE, gc!, nc, p, optim
+        estim, Hp, Hc, M_Hp, N_Hc, L_Hp, Cwt, Ewt, JE, gc!, nc, p, 
+        transcription, optim
     )
 end