Skip to content

Commit 44a5589

Browse files
authored
Merge pull request #201 from JuliaControl/move_block
added: move blocking feature in `LinMPC`, `ExplicitMPC` and `NonLinMPC`
2 parents 03cb818 + bd24cb7 commit 44a5589

18 files changed

+208
-90
lines changed

Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "ModelPredictiveControl"
22
uuid = "61f9bdb8-6ae4-484a-811f-bbf86720c31c"
33
authors = ["Francis Gagnon"]
4-
version = "1.7.0"
4+
version = "1.8.0"
55

66
[deps]
77
ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"

README.md

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -82,6 +82,7 @@ for more detailed examples.
8282
- input setpoint tracking
8383
- terminal costs
8484
- custom economic costs (economic model predictive control)
85+
- control horizon distinct from prediction horizon and custom move blocking
8586
- adaptive linear model predictive controller
8687
- manual model modification
8788
- automatic successive linearization of a nonlinear model

docs/src/internals/predictive_control.md

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,7 @@ The prediction methodology of this module is mainly based on Maciejowski textboo
1212
## Controller Construction
1313

1414
```@docs
15+
ModelPredictiveControl.move_blocking
1516
ModelPredictiveControl.init_ZtoΔU
1617
ModelPredictiveControl.init_ZtoU
1718
ModelPredictiveControl.init_predmat

docs/src/public/predictive_control.md

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -48,7 +48,8 @@ are shifted by one time step:
4848
in which ``\mathbf{U}`` and ``\mathbf{R̂_u}`` are both vectors of `nu*Hp` elements. Defining
4949
the manipulated input increment as ``\mathbf{Δu}(k+j) = \mathbf{u}(k+j) - \mathbf{u}(k+j-1)``,
5050
the control horizon ``H_c`` enforces that ``\mathbf{Δu}(k+j) = \mathbf{0}`` for ``j ≥ H_c``.
51-
For this reason, the vector that collects them is truncated up to ``k+H_c-1``:
51+
For this reason, the vector that collects them is truncated up to ``k+H_c-1`` (without
52+
any custom move blocking):
5253

5354
```math
5455
\mathbf{ΔU} =

src/controller/construct.jl

Lines changed: 66 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -456,6 +456,72 @@ end
456456
"Return `0` when model is not a [`LinModel`](@ref)."
457457
estimate_delays(::SimModel) = 0
458458

459+
460+
@doc raw"""
461+
move_blocking(Hp::Int, Hc::AbstractVector{Int}) -> nb
462+
463+
Get the move blocking vector `nb` from the `Hc` argument, and modify it to match `Hp`.
464+
465+
This feature is also known as manipulated variable blocking. The argument `Hc` is
466+
interpreted as the move blocking vector `nb`. It specifies the length of each step (or
467+
"block") in the ``\mathbf{ΔU}`` vector, to customize the pattern (in time steps, thus
468+
strictly positive integers):
469+
```math
470+
\mathbf{n_b} = \begin{bmatrix} n_1 & n_2 & \cdots & n_{H_c} \end{bmatrix}'
471+
```
472+
The vector that includes all the manipulated input increments ``\mathbf{Δu}`` is then
473+
defined as:
474+
```math
475+
\mathbf{ΔU} = \begin{bmatrix}
476+
\mathbf{Δu}(k + 0) \\[0.1em]
477+
\mathbf{Δu}(k + ∑_{i=1}^1 n_i) \\[0.1em]
478+
\mathbf{Δu}(k + ∑_{i=1}^2 n_i) \\[0.1em]
479+
\vdots \\[0.1em]
480+
\mathbf{Δu}(k + ∑_{i=1}^{H_c-1} n_i)
481+
\end{bmatrix}
482+
```
483+
The provided `nb` vector is modified to ensure `sum(nb) == Hp`:
484+
- If `sum(nb) < Hp`, a new element is pushed to `nb` with the value `Hp - sum(nb)`.
485+
- If `sum(nb) > Hp`, the intervals are truncated until `sum(nb) == Hp`. For example, if
486+
`Hp = 10` and `nb = [1, 2, 3, 6, 7]`, then `nb` is truncated to `[1, 2, 3, 4]`.
487+
"""
488+
function move_blocking(Hp_arg::Int, Hc_arg::AbstractVector{Int})
489+
Hp = Hp_arg
490+
nb = Hc_arg
491+
all(>(0), nb) || throw(ArgumentError("Move blocking vector must be strictly positive integers."))
492+
if sum(nb) < Hp
493+
newblock = [Hp - sum(nb)]
494+
nb = [nb; newblock]
495+
elseif sum(nb) > Hp
496+
nb = nb[begin:findfirst((Hp), cumsum(nb))]
497+
if sum(nb) > Hp
498+
# if the last block is too large, it is truncated to fit Hp:
499+
nb[end] = Hp - @views sum(nb[begin:end-1])
500+
end
501+
end
502+
return nb
503+
end
504+
505+
"""
506+
move_blocking(Hp::Int, Hc::Int) -> nb
507+
508+
Construct a move blocking vector `nb` that match the provided `Hp` and `Hc` integers.
509+
510+
The vector is filled with `1`s, except for the last element which is `Hp - Hc + 1`.
511+
"""
512+
function move_blocking(Hp_arg::Int, Hc_arg::Int)
513+
Hp, Hc = Hp_arg, Hc_arg
514+
nb = fill(1, Hc)
515+
if Hc > 0 # if Hc < 1, it will crash later with a clear error message
516+
nb[end] = Hp - Hc + 1
517+
end
518+
return nb
519+
end
520+
521+
"Get the actual control Horizon `Hc` (integer) from the move blocking vector `nb`."
522+
get_Hc(nb::AbstractVector{Int}) = length(nb)
523+
524+
459525
"""
460526
validate_args(mpc::PredictiveController, ry, d, D̂, R̂y, R̂u)
461527

src/controller/explicitmpc.jl

Lines changed: 11 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@ struct ExplicitMPC{
1010
Hp::Int
1111
Hc::Int
1212
::Int
13+
nb::Vector{Int}
1314
weights::CW
1415
R̂u::Vector{NT}
1516
R̂y::Vector{NT}
@@ -39,7 +40,7 @@ struct ExplicitMPC{
3940
Dop::Vector{NT}
4041
buffer::PredictiveControllerBuffer{NT}
4142
function ExplicitMPC{NT}(
42-
estim::SE, Hp, Hc, weights::CW
43+
estim::SE, Hp, Hc, nb, weights::CW
4344
) where {NT<:Real, SE<:StateEstimator, CW<:ControllerWeights}
4445
model = estim.model
4546
nu, ny, nd, nx̂ = model.nu, model.ny, model.nd, estim.nx̂
@@ -50,7 +51,7 @@ struct ExplicitMPC{
5051
lastu0 = zeros(NT, nu)
5152
transcription = SingleShooting() # explicit MPC only supports SingleShooting
5253
PΔu = init_ZtoΔU(estim, transcription, Hp, Hc)
53-
Pu, Tu = init_ZtoU(estim, transcription, Hp, Hc)
54+
Pu, Tu = init_ZtoU(estim, transcription, Hp, Hc, nb)
5455
E, G, J, K, V, B = init_predmat(model, estim, transcription, Hp, Hc)
5556
# dummy val (updated just before optimization):
5657
F = zeros(NT, ny*Hp)
@@ -70,7 +71,7 @@ struct ExplicitMPC{
7071
estim,
7172
transcription,
7273
Z̃, ŷ,
73-
Hp, Hc, nϵ,
74+
Hp, Hc, nϵ, nb,
7475
weights,
7576
R̂u, R̂y,
7677
lastu0,
@@ -129,12 +130,12 @@ ExplicitMPC controller with a sample time Ts = 4.0 s, SteadyKalmanFilter estimat
129130
function ExplicitMPC(
130131
model::LinModel;
131132
Hp::Int = default_Hp(model),
132-
Hc::Int = DEFAULT_HC,
133+
Hc::IntVectorOrInt = DEFAULT_HC,
133134
Mwt = fill(DEFAULT_MWT, model.ny),
134135
Nwt = fill(DEFAULT_NWT, model.nu),
135136
Lwt = fill(DEFAULT_LWT, model.nu),
136137
M_Hp = Diagonal(repeat(Mwt, Hp)),
137-
N_Hc = Diagonal(repeat(Nwt, Hc)),
138+
N_Hc = Diagonal(repeat(Nwt, get_Hc(move_blocking(Hp, Hc)))),
138139
L_Hp = Diagonal(repeat(Lwt, Hp)),
139140
kwargs...
140141
)
@@ -152,12 +153,12 @@ Use custom state estimator `estim` to construct `ExplicitMPC`.
152153
function ExplicitMPC(
153154
estim::SE;
154155
Hp::Int = default_Hp(estim.model),
155-
Hc::Int = DEFAULT_HC,
156+
Hc::IntVectorOrInt = DEFAULT_HC,
156157
Mwt = fill(DEFAULT_MWT, estim.model.ny),
157158
Nwt = fill(DEFAULT_NWT, estim.model.nu),
158159
Lwt = fill(DEFAULT_LWT, estim.model.nu),
159160
M_Hp = Diagonal(repeat(Mwt, Hp)),
160-
N_Hc = Diagonal(repeat(Nwt, Hc)),
161+
N_Hc = Diagonal(repeat(Nwt, get_Hc(move_blocking(Hp, Hc)))),
161162
L_Hp = Diagonal(repeat(Lwt, Hp)),
162163
) where {NT<:Real, SE<:StateEstimator{NT}}
163164
isa(estim.model, LinModel) || error(MSG_LINMODEL_ERR)
@@ -166,8 +167,10 @@ function ExplicitMPC(
166167
@warn("prediction horizon Hp ($Hp) ≤ estimated number of delays in model "*
167168
"($nk), the closed-loop system may be unstable or zero-gain (unresponsive)")
168169
end
170+
nb = move_blocking(Hp, Hc)
171+
Hc = get_Hc(nb)
169172
weights = ControllerWeights{NT}(estim.model, Hp, Hc, M_Hp, N_Hc, L_Hp)
170-
return ExplicitMPC{NT}(estim, Hp, Hc, weights)
173+
return ExplicitMPC{NT}(estim, Hp, Hc, nb, weights)
171174
end
172175

173176
setconstraint!(::ExplicitMPC; kwargs...) = error("ExplicitMPC does not support constraints.")

src/controller/linmpc.jl

Lines changed: 16 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,7 @@ struct LinMPC{
1919
Hp::Int
2020
Hc::Int
2121
::Int
22+
nb::Vector{Int}
2223
weights::CW
2324
R̂u::Vector{NT}
2425
R̂y::Vector{NT}
@@ -47,7 +48,7 @@ struct LinMPC{
4748
Dop::Vector{NT}
4849
buffer::PredictiveControllerBuffer{NT}
4950
function LinMPC{NT}(
50-
estim::SE, Hp, Hc, weights::CW,
51+
estim::SE, Hp, Hc, nb, weights::CW,
5152
transcription::TM, optim::JM
5253
) where {
5354
NT<:Real,
@@ -63,7 +64,7 @@ struct LinMPC{
6364
R̂y, R̂u, Tu_lastu0 = zeros(NT, ny*Hp), zeros(NT, nu*Hp), zeros(NT, nu*Hp)
6465
lastu0 = zeros(NT, nu)
6566
PΔu = init_ZtoΔU(estim, transcription, Hp, Hc)
66-
Pu, Tu = init_ZtoU(estim, transcription, Hp, Hc)
67+
Pu, Tu = init_ZtoU(estim, transcription, Hp, Hc, nb)
6768
E, G, J, K, V, B, ex̂, gx̂, jx̂, kx̂, vx̂, bx̂ = init_predmat(
6869
model, estim, transcription, Hp, Hc
6970
)
@@ -90,7 +91,7 @@ struct LinMPC{
9091
mpc = new{NT, SE, CW, TM, JM}(
9192
estim, transcription, optim, con,
9293
Z̃, ŷ,
93-
Hp, Hc, nϵ,
94+
Hp, Hc, nϵ, nb,
9495
weights,
9596
R̂u, R̂y,
9697
lastu0,
@@ -145,8 +146,9 @@ arguments. This controller allocates memory at each time step for the optimizati
145146
146147
# Arguments
147148
- `model::LinModel` : model used for controller predictions and state estimations.
148-
- `Hp=10+nk` : prediction horizon ``H_p``, `nk` is the number of delays in `model`.
149-
- `Hc=2` : control horizon ``H_c``.
149+
- `Hp::Int=10+nk` : prediction horizon ``H_p``, `nk` is the number of delays in `model`.
150+
- `Hc::Union{Int, Vector{Int}}=2` : control horizon ``H_c``, custom move blocking pattern is
151+
specified with a vector of integers (see [`move_blocking`](@ref) for details).
150152
- `Mwt=fill(1.0,model.ny)` : main diagonal of ``\mathbf{M}`` weight matrix (vector).
151153
- `Nwt=fill(0.1,model.nu)` : main diagonal of ``\mathbf{N}`` weight matrix (vector).
152154
- `Lwt=fill(0.0,model.nu)` : main diagonal of ``\mathbf{L}`` weight matrix (vector).
@@ -205,12 +207,12 @@ LinMPC controller with a sample time Ts = 4.0 s, OSQP optimizer, SteadyKalmanFil
205207
function LinMPC(
206208
model::LinModel;
207209
Hp::Int = default_Hp(model),
208-
Hc::Int = DEFAULT_HC,
210+
Hc::IntVectorOrInt = DEFAULT_HC,
209211
Mwt = fill(DEFAULT_MWT, model.ny),
210212
Nwt = fill(DEFAULT_NWT, model.nu),
211213
Lwt = fill(DEFAULT_LWT, model.nu),
212214
M_Hp = Diagonal(repeat(Mwt, Hp)),
213-
N_Hc = Diagonal(repeat(Nwt, Hc)),
215+
N_Hc = Diagonal(repeat(Nwt, get_Hc(move_blocking(Hp, Hc)))),
214216
L_Hp = Diagonal(repeat(Lwt, Hp)),
215217
Cwt = DEFAULT_CWT,
216218
transcription::TranscriptionMethod = DEFAULT_LINMPC_TRANSCRIPTION,
@@ -227,7 +229,8 @@ end
227229
228230
Use custom state estimator `estim` to construct `LinMPC`.
229231
230-
`estim.model` must be a [`LinModel`](@ref). Else, a [`NonLinMPC`](@ref) is required.
232+
`estim.model` must be a [`LinModel`](@ref). Else, a [`NonLinMPC`](@ref) is required. See
233+
[`ManualEstimator`](@ref) for linear controllers with nonlinear state estimation.
231234
232235
# Examples
233236
```jldoctest
@@ -248,12 +251,12 @@ LinMPC controller with a sample time Ts = 4.0 s, OSQP optimizer, KalmanFilter es
248251
function LinMPC(
249252
estim::SE;
250253
Hp::Int = default_Hp(estim.model),
251-
Hc::Int = DEFAULT_HC,
254+
Hc::IntVectorOrInt = DEFAULT_HC,
252255
Mwt = fill(DEFAULT_MWT, estim.model.ny),
253256
Nwt = fill(DEFAULT_NWT, estim.model.nu),
254257
Lwt = fill(DEFAULT_LWT, estim.model.nu),
255258
M_Hp = Diagonal(repeat(Mwt, Hp)),
256-
N_Hc = Diagonal(repeat(Nwt, Hc)),
259+
N_Hc = Diagonal(repeat(Nwt, get_Hc(move_blocking(Hp, Hc)))),
257260
L_Hp = Diagonal(repeat(Lwt, Hp)),
258261
Cwt = DEFAULT_CWT,
259262
transcription::TranscriptionMethod = DEFAULT_LINMPC_TRANSCRIPTION,
@@ -265,8 +268,10 @@ function LinMPC(
265268
@warn("prediction horizon Hp ($Hp) ≤ estimated number of delays in model "*
266269
"($nk), the closed-loop system may be unstable or zero-gain (unresponsive)")
267270
end
271+
nb = move_blocking(Hp, Hc)
272+
Hc = get_Hc(nb)
268273
weights = ControllerWeights{NT}(estim.model, Hp, Hc, M_Hp, N_Hc, L_Hp, Cwt)
269-
return LinMPC{NT}(estim, Hp, Hc, weights, transcription, optim)
274+
return LinMPC{NT}(estim, Hp, Hc, nb, weights, transcription, optim)
270275
end
271276

272277
"""

src/controller/nonlinmpc.jl

Lines changed: 15 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -33,6 +33,7 @@ struct NonLinMPC{
3333
Hp::Int
3434
Hc::Int
3535
::Int
36+
nb::Vector{Int}
3637
weights::CW
3738
JE::JEfunc
3839
p::PT
@@ -63,7 +64,7 @@ struct NonLinMPC{
6364
Dop::Vector{NT}
6465
buffer::PredictiveControllerBuffer{NT}
6566
function NonLinMPC{NT}(
66-
estim::SE, Hp, Hc, weights::CW,
67+
estim::SE, Hp, Hc, nb, weights::CW,
6768
JE::JEfunc, gc!::GCfunc, nc, p::PT,
6869
transcription::TM, optim::JM,
6970
gradient::GB, jacobian::JB
@@ -86,7 +87,7 @@ struct NonLinMPC{
8687
R̂y, R̂u, Tu_lastu0 = zeros(NT, ny*Hp), zeros(NT, nu*Hp), zeros(NT, nu*Hp)
8788
lastu0 = zeros(NT, nu)
8889
PΔu = init_ZtoΔU(estim, transcription, Hp, Hc)
89-
Pu, Tu = init_ZtoU(estim, transcription, Hp, Hc)
90+
Pu, Tu = init_ZtoU(estim, transcription, Hp, Hc, nb)
9091
E, G, J, K, V, B, ex̂, gx̂, jx̂, kx̂, vx̂, bx̂ = init_predmat(
9192
model, estim, transcription, Hp, Hc
9293
)
@@ -116,7 +117,7 @@ struct NonLinMPC{
116117
estim, transcription, optim, con,
117118
gradient, jacobian,
118119
Z̃, ŷ,
119-
Hp, Hc, nϵ,
120+
Hp, Hc, nϵ, nb,
120121
weights,
121122
JE, p,
122123
R̂u, R̂y,
@@ -188,8 +189,10 @@ This controller allocates memory at each time step for the optimization.
188189
189190
# Arguments
190191
- `model::SimModel` : model used for controller predictions and state estimations.
191-
- `Hp=nothing`: prediction horizon ``H_p``, must be specified for [`NonLinModel`](@ref).
192-
- `Hc=2` : control horizon ``H_c``.
192+
- `Hp::Int=10+nk` : prediction horizon ``H_p``, `nk` is the number of delays if `model` is a
193+
[`LinModel`](@ref) (must be specified otherwise).
194+
- `Hc::Union{Int, Vector{Int}}=2` : control horizon ``H_c``, custom move blocking pattern is
195+
specified with a vector of integers (see [`move_blocking`](@ref) for details).
193196
- `Mwt=fill(1.0,model.ny)` : main diagonal of ``\mathbf{M}`` weight matrix (vector).
194197
- `Nwt=fill(0.1,model.nu)` : main diagonal of ``\mathbf{N}`` weight matrix (vector).
195198
- `Lwt=fill(0.0,model.nu)` : main diagonal of ``\mathbf{L}`` weight matrix (vector).
@@ -281,12 +284,12 @@ NonLinMPC controller with a sample time Ts = 10.0 s, Ipopt optimizer, UnscentedK
281284
function NonLinMPC(
282285
model::SimModel;
283286
Hp::Int = default_Hp(model),
284-
Hc::Int = DEFAULT_HC,
287+
Hc::IntVectorOrInt = DEFAULT_HC,
285288
Mwt = fill(DEFAULT_MWT, model.ny),
286289
Nwt = fill(DEFAULT_NWT, model.nu),
287290
Lwt = fill(DEFAULT_LWT, model.nu),
288291
M_Hp = Diagonal(repeat(Mwt, Hp)),
289-
N_Hc = Diagonal(repeat(Nwt, Hc)),
292+
N_Hc = Diagonal(repeat(Nwt, get_Hc(move_blocking(Hp, Hc)))),
290293
L_Hp = Diagonal(repeat(Lwt, Hp)),
291294
Cwt = DEFAULT_CWT,
292295
Ewt = DEFAULT_EWT,
@@ -338,12 +341,12 @@ NonLinMPC controller with a sample time Ts = 10.0 s, Ipopt optimizer, UnscentedK
338341
function NonLinMPC(
339342
estim::SE;
340343
Hp::Int = default_Hp(estim.model),
341-
Hc::Int = DEFAULT_HC,
344+
Hc::IntVectorOrInt = DEFAULT_HC,
342345
Mwt = fill(DEFAULT_MWT, estim.model.ny),
343346
Nwt = fill(DEFAULT_NWT, estim.model.nu),
344347
Lwt = fill(DEFAULT_LWT, estim.model.nu),
345348
M_Hp = Diagonal(repeat(Mwt, Hp)),
346-
N_Hc = Diagonal(repeat(Nwt, Hc)),
349+
N_Hc = Diagonal(repeat(Nwt, get_Hc(move_blocking(Hp, Hc)))),
347350
L_Hp = Diagonal(repeat(Lwt, Hp)),
348351
Cwt = DEFAULT_CWT,
349352
Ewt = DEFAULT_EWT,
@@ -365,11 +368,13 @@ function NonLinMPC(
365368
@warn("prediction horizon Hp ($Hp) ≤ estimated number of delays in model "*
366369
"($nk), the closed-loop system may be unstable or zero-gain (unresponsive)")
367370
end
371+
nb = move_blocking(Hp, Hc)
372+
Hc = get_Hc(nb)
368373
validate_JE(NT, JE)
369374
gc! = get_mutating_gc(NT, gc)
370375
weights = ControllerWeights{NT}(estim.model, Hp, Hc, M_Hp, N_Hc, L_Hp, Cwt, Ewt)
371376
return NonLinMPC{NT}(
372-
estim, Hp, Hc, weights, JE, gc!, nc, p, transcription, optim, gradient, jacobian
377+
estim, Hp, Hc, nb, weights, JE, gc!, nc, p, transcription, optim, gradient, jacobian
373378
)
374379
end
375380

0 commit comments

Comments
 (0)