Skip to content

Commit 391bf73

Browse files
committed
changed: slack var for MHE begins the decision vector
1 parent 27b109e commit 391bf73

File tree

8 files changed

+406
-320
lines changed

8 files changed

+406
-320
lines changed

src/controller/construct.jl

Lines changed: 41 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,36 @@
1+
"Include all the data for the constraints of [`PredictiveController`](@ref)"
2+
struct ControllerConstraint{NT<:Real}
3+
ẽx̂ ::Matrix{NT}
4+
fx̂ ::Vector{NT}
5+
gx̂ ::Matrix{NT}
6+
jx̂ ::Matrix{NT}
7+
kx̂ ::Matrix{NT}
8+
vx̂ ::Matrix{NT}
9+
Umin ::Vector{NT}
10+
Umax ::Vector{NT}
11+
ΔŨmin ::Vector{NT}
12+
ΔŨmax ::Vector{NT}
13+
Ymin ::Vector{NT}
14+
Ymax ::Vector{NT}
15+
x̂min ::Vector{NT}
16+
x̂max ::Vector{NT}
17+
A_Umin ::Matrix{NT}
18+
A_Umax ::Matrix{NT}
19+
A_ΔŨmin ::Matrix{NT}
20+
A_ΔŨmax ::Matrix{NT}
21+
A_Ymin ::Matrix{NT}
22+
A_Ymax ::Matrix{NT}
23+
A_x̂min ::Matrix{NT}
24+
A_x̂max ::Matrix{NT}
25+
A ::Matrix{NT}
26+
b ::Vector{NT}
27+
i_b ::BitVector
28+
C_ymin ::Vector{NT}
29+
C_ymax ::Vector{NT}
30+
c_x̂min ::Vector{NT}
31+
c_x̂max ::Vector{NT}
32+
i_g ::BitVector
33+
end
134

235
@doc raw"""
336
setconstraint!(mpc::PredictiveController; <keyword arguments>) -> mpc
@@ -36,8 +69,8 @@ constraints are all soft by default. See Extended Help for time-varying constrai
3669
- `c_ymax = fill(1.0,ny)` : `ymax` softness weights ``\mathbf{c_{y_{max}}}``.
3770
- `c_x̂min = fill(1.0,nx̂)` : `x̂min` softness weights ``\mathbf{c_{x̂_{min}}}``.
3871
- `c_x̂max = fill(1.0,nx̂)` : `x̂max` softness weights ``\mathbf{c_{x̂_{max}}}``.
39-
- all the keyword arguments above but with a capital letter, except for the terminal
40-
constraints, e.g. `Ymax` or `C_Δumin` : for time-varying constraints (see Extended Help).
72+
- all the keyword arguments above but with a first capital letter, except for the terminal
73+
constraints, e.g. `Ymax` or `C_Δumin`: for time-varying constraints (see Extended Help).
4174
4275
# Examples
4376
```jldoctest
@@ -80,7 +113,7 @@ LinMPC controller with a sample time Ts = 4.0 s, OSQP optimizer, SteadyKalmanFil
80113
\mathbf{Y_{min} - C_{y_{min}}} ϵ ≤&&\ \mathbf{Ŷ} &≤ \mathbf{Y_{max} + C_{y_{max}}} ϵ
81114
\end{alignat*}
82115
```
83-
For this, use the same keyword arguments as above but with a capital letter:
116+
For this, use the same keyword arguments as above but with a first capital letter:
84117
- `Umin` / `Umax` / `C_umin` / `C_umax` : ``\mathbf{U}`` constraints `(nu*Hp,)`.
85118
- `ΔUmin` / `ΔUmax` / `C_Δumin` / `C_Δumax` : ``\mathbf{ΔU}`` constraints `(nu*Hc,)`.
86119
- `Ymin` / `Ymax` / `C_ymin` / `C_ymax` : ``\mathbf{Ŷ}`` constraints `(ny*Hp,)`.
@@ -142,7 +175,7 @@ function setconstraint!(
142175
isnothing(C_Δumax) && !isnothing(c_Δumax) && (C_Δumax = repeat(c_Δumax, Hc))
143176
isnothing(C_ymin) && !isnothing(c_ymin) && (C_ymin = repeat(c_ymin, Hp))
144177
isnothing(C_ymax) && !isnothing(c_ymax) && (C_ymax = repeat(c_ymax, Hp))
145-
if !all(isnothing.([C_umin, C_umax, C_Δumin, C_Δumax, C_ymin, C_ymax, c_x̂min, c_x̂max]))
178+
if !all(isnothing.((C_umin, C_umax, C_Δumin, C_Δumax, C_ymin, C_ymax, c_x̂min, c_x̂max)))
146179
!isinf(C) || throw(ArgumentError("Slack variable weight Cwt must be finite to set softness parameters"))
147180
notSolvedYet || error("Cannot set softness parameters after calling moveinput!")
148181
end
@@ -725,7 +758,7 @@ function relaxŶ(::LinModel{NT}, C, C_ymin, C_ymax, E) where {NT<:Real}
725758
if !isinf(C) # ΔŨ = [ΔU; ϵ]
726759
# ϵ impacts predicted output constraint calculations:
727760
A_Ymin, A_Ymax = -[E C_ymin], [E -C_ymax]
728-
# ϵ has no impact on output predictions
761+
# ϵ has no impact on output predictions:
729762
= [E zeros(NT, size(E, 1), 1)]
730763
else # ΔŨ = ΔU (only hard constraints)
731764
= E
@@ -764,9 +797,9 @@ the inequality constraints:
764797
"""
765798
function relaxterminal(::LinModel{NT}, C, c_x̂min, c_x̂max, ex̂) where {NT<:Real}
766799
if !isinf(C) # ΔŨ = [ΔU; ϵ]
767-
# ϵ impacts terminal constraint calculations:
800+
# ϵ impacts terminal state constraint calculations:
768801
A_x̂min, A_x̂max = -[ex̂ c_x̂min], [ex̂ -c_x̂max]
769-
# ϵ has no impact on terminal state predictions
802+
# ϵ has no impact on terminal state predictions:
770803
ẽx̂ = [ex̂ zeros(NT, size(ex̂, 1), 1)]
771804
else # ΔŨ = ΔU (only hard constraints)
772805
ẽx̂ = ex̂
@@ -834,6 +867,6 @@ function validate_weights(model, Hp, Hc, M_Hp, N_Hc, L_Hp, C, E=nothing)
834867
(!isdiag(N_Hc) || any(diag(N_Hc).<0)) && throw(ArgumentError("N_Hc should be a positive semidefinite diagonal matrix"))
835868
(!isdiag(L_Hp) || any(diag(L_Hp).<0)) && throw(ArgumentError("L_Hp should be a positive semidefinite diagonal matrix"))
836869
size(C) () && throw(ArgumentError("Cwt should be a real scalar"))
837-
C < 0 && throw(ArgumentError("Cwt weight should be ≥ 0"))
870+
C < 0 && throw(ArgumentError("Cwt weight should be ≥ 0"))
838871
!isnothing(E) && size(E) () && throw(ArgumentError("Ewt should be a real scalar"))
839872
end

src/controller/execute.jl

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -357,10 +357,9 @@ function optim_objective!(mpc::PredictiveController{NT}) where {NT<:Real}
357357
optim = mpc.optim
358358
model = mpc.estim.model
359359
ΔŨvar::Vector{VariableRef} = optim[:ΔŨvar]
360-
lastΔŨ = mpc.ΔŨ
361360
# initial ΔŨ (warm-start): [Δu_{k-1}(k); Δu_{k-1}(k+1); ... ; 0_{nu × 1}; ϵ_{k-1}]
362-
ϵ0 = !isinf(mpc.C) ? [lastΔŨ[end]] : empty(mpc.ΔŨ)
363-
ΔŨ0 = [lastΔŨ[(model.nu+1):(mpc.Hc*model.nu)]; zeros(NT, model.nu); ϵ0]
361+
ϵ0 = !isinf(mpc.C) ? mpc.ΔŨ[end] : empty(mpc.ΔŨ)
362+
ΔŨ0 = [mpc.ΔŨ[(model.nu+1):(mpc.Hc*model.nu)]; zeros(NT, model.nu); ϵ0]
364363
set_start_value.(ΔŨvar, ΔŨ0)
365364
set_objective_linear_coef!(mpc, ΔŨvar)
366365
try

src/controller/nonlinmpc.jl

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -131,7 +131,7 @@ This method uses the default state estimator :
131131
- `Lwt=fill(0.0,model.nu)` : main diagonal of ``\mathbf{L}`` weight matrix (vector).
132132
- `Cwt=1e5` : slack variable weight ``C`` (scalar), use `Cwt=Inf` for hard constraints only.
133133
- `Ewt=0.0` : economic costs weight ``E`` (scalar).
134-
- `JE=(_,_,_)->0.0` : economic function ``J_E(\mathbf{U}_E, \mathbf{Ŷ}_E, \mathbf{D̂}_E)``.
134+
- `JE=(_,_,_)->0` : economic function ``J_E(\mathbf{U}_E, \mathbf{Ŷ}_E, \mathbf{D̂}_E)``.
135135
- `M_Hp` / `N_Hc` / `L_Hp` : diagonal matrices ``\mathbf{M}_{H_p}, \mathbf{N}_{H_c},
136136
\mathbf{L}_{H_p}``, for time-varying weights (generated from `Mwt/Nwt/Lwt` args if omitted).
137137
- `optim=JuMP.Model(Ipopt.Optimizer)` : nonlinear optimizer used in the predictive
@@ -176,7 +176,7 @@ function NonLinMPC(
176176
Lwt = fill(DEFAULT_LWT, model.nu),
177177
Cwt = DEFAULT_CWT,
178178
Ewt = DEFAULT_EWT,
179-
JE::Function = (_,_,_) -> 0.0,
179+
JE::Function = (_,_,_) -> 0,
180180
M_Hp = nothing,
181181
N_Hc = nothing,
182182
L_Hp = nothing,
@@ -196,7 +196,7 @@ function NonLinMPC(
196196
Lwt = fill(DEFAULT_LWT, model.nu),
197197
Cwt = DEFAULT_CWT,
198198
Ewt = DEFAULT_EWT,
199-
JE::Function = (_,_,_) -> 0.0,
199+
JE::Function = (_,_,_) -> 0,
200200
M_Hp = nothing,
201201
N_Hc = nothing,
202202
L_Hp = nothing,
@@ -239,7 +239,7 @@ function NonLinMPC(
239239
Lwt = fill(DEFAULT_LWT, estim.model.nu),
240240
Cwt = DEFAULT_CWT,
241241
Ewt = DEFAULT_EWT,
242-
JE::JEFunc = (_,_,_) -> 0.0,
242+
JE::JEFunc = (_,_,_) -> 0,
243243
M_Hp = nothing,
244244
N_Hc = nothing,
245245
L_Hp = nothing,
@@ -402,7 +402,7 @@ The method mutates the `g` vector in argument and returns it.
402402
"""
403403
function con_nonlinprog!(g, mpc::NonLinMPC, ::SimModel, x̂end, Ŷ, ΔŨ)
404404
nx̂, nŶ = mpc.estim.nx̂, length(Ŷ)
405-
ϵ = !isinf(mpc.C) ? ΔŨ[end] : 0.0 # ϵ = 0.0 if Cwt=Inf (meaning: no relaxation)
405+
ϵ = isinf(mpc.C) ? 0 : ΔŨ[end] # ϵ = 0 if Cwt=Inf (meaning: no relaxation)
406406
for i in eachindex(g)
407407
mpc.con.i_g[i] || continue
408408
if i nŶ

src/estimator/mhe.jl

Lines changed: 4 additions & 144 deletions
Original file line numberDiff line numberDiff line change
@@ -1,142 +1,5 @@
1-
@doc raw"""
2-
Include all the data for the constraints of [`MovingHorizonEstimator`](@ref).
3-
4-
The bounds on the estimated state at arrival ``\mathbf{x̂}_k(k-N_k+1)`` is separated from
5-
the other state constraints ``\mathbf{x̂}_k(k-N_k+2), \mathbf{x̂}_k(k-N_k+3), ...`` since
6-
the former is always a linear inequality constraint (it's a decision variable). The fields
7-
`x̂min` and `x̂max` refer to the bounds at the arrival, and `X̂min` and `X̂max`, the others.
8-
"""
9-
struct EstimatorConstraint{NT<:Real}
10-
Ẽx̂ ::Matrix{NT}
11-
Fx̂ ::Vector{NT}
12-
Gx̂ ::Matrix{NT}
13-
Jx̂ ::Matrix{NT}
14-
x̂min ::Vector{NT}
15-
x̂max ::Vector{NT}
16-
X̂min ::Vector{NT}
17-
X̂max ::Vector{NT}
18-
Ŵmin ::Vector{NT}
19-
Ŵmax ::Vector{NT}
20-
V̂min ::Vector{NT}
21-
V̂max ::Vector{NT}
22-
A_x̂min ::Matrix{NT}
23-
A_x̂max ::Matrix{NT}
24-
A_X̂min ::Matrix{NT}
25-
A_X̂max ::Matrix{NT}
26-
A_Ŵmin ::Matrix{NT}
27-
A_Ŵmax ::Matrix{NT}
28-
A_V̂min ::Matrix{NT}
29-
A_V̂max ::Matrix{NT}
30-
A ::Matrix{NT}
31-
b ::Vector{NT}
32-
i_b ::BitVector
33-
i_g ::BitVector
34-
end
35-
36-
struct MovingHorizonEstimator{
37-
NT<:Real,
38-
SM<:SimModel,
39-
JM<:JuMP.GenericModel
40-
} <: StateEstimator{NT}
41-
model::SM
42-
# note: `NT` and the number type `JNT` in `JuMP.GenericModel{JNT}` can be
43-
# different since solvers that support non-Float64 are scarce.
44-
optim::JM
45-
con::EstimatorConstraint{NT}
46-
::Vector{NT}
47-
lastu0::Vector{NT}
48-
::Vector{NT}
49-
He::Int
50-
i_ym::Vector{Int}
51-
nx̂ ::Int
52-
nym::Int
53-
nyu::Int
54-
nxs::Int
55-
As ::Matrix{NT}
56-
Cs_u::Matrix{NT}
57-
Cs_y::Matrix{NT}
58-
nint_u ::Vector{Int}
59-
nint_ym::Vector{Int}
60-
::Matrix{NT}
61-
B̂u ::Matrix{NT}
62-
::Matrix{NT}
63-
B̂d ::Matrix{NT}
64-
D̂d ::Matrix{NT}
65-
::Matrix{NT}
66-
F ::Vector{NT}
67-
G ::Matrix{NT}
68-
J ::Matrix{NT}
69-
ẽx̄::Matrix{NT}
70-
fx̄::Vector{NT}
71-
::Hermitian{NT, Matrix{NT}}
72-
::Vector{NT}
73-
p::Vector{NT}
74-
P̂0::Hermitian{NT, Matrix{NT}}
75-
::Hermitian{NT, Matrix{NT}}
76-
::Hermitian{NT, Matrix{NT}}
77-
invP̄::Hermitian{NT, Matrix{NT}}
78-
invQ̂_He::Hermitian{NT, Matrix{NT}}
79-
invR̂_He::Hermitian{NT, Matrix{NT}}
80-
::Matrix{NT}
81-
::Union{Vector{NT}, Missing}
82-
Ym::Union{Vector{NT}, Missing}
83-
U ::Union{Vector{NT}, Missing}
84-
D ::Union{Vector{NT}, Missing}
85-
::Union{Vector{NT}, Missing}
86-
x̂arr_old::Vector{NT}
87-
P̂arr_old::Hermitian{NT, Matrix{NT}}
88-
Nk::Vector{Int}
89-
function MovingHorizonEstimator{NT, SM, JM}(
90-
model::SM, He, i_ym, nint_u, nint_ym, P̂0, Q̂, R̂, optim::JM
91-
) where {NT<:Real, SM<:SimModel{NT}, JM<:JuMP.GenericModel}
92-
nu, nd = model.nu, model.nd
93-
He < 1 && throw(ArgumentError("Estimation horizon He should be ≥ 1"))
94-
nym, nyu = validate_ym(model, i_ym)
95-
As, Cs_u, Cs_y, nint_u, nint_ym = init_estimstoch(model, i_ym, nint_u, nint_ym)
96-
nxs = size(As, 1)
97-
nx̂ = model.nx + nxs
98-
nŵ = nx̂
99-
Â, B̂u, Ĉ, B̂d, D̂d = augment_model(model, As, Cs_u, Cs_y)
100-
validate_kfcov(nym, nx̂, Q̂, R̂, P̂0)
101-
lastu0 = zeros(NT, model.nu)
102-
= [zeros(NT, model.nx); zeros(NT, nxs)]
103-
P̂0 = Hermitian(P̂0, :L)
104-
Q̂, R̂ = Hermitian(Q̂, :L), Hermitian(R̂, :L)
105-
invP̄ = Hermitian(inv(P̂0), :L)
106-
invQ̂_He = Hermitian(repeatdiag(inv(Q̂), He), :L)
107-
invR̂_He = Hermitian(repeatdiag(inv(R̂), He), :L)
108-
= zeros(NT, nx̂, nym)
109-
E, F, G, J, ex̄, fx̄, Ex̂, Fx̂, Gx̂, Jx̂ = init_predmat_mhe(
110-
model, He, i_ym, Â, B̂u, Ĉ, B̂d, D̂d
111-
)
112-
con, Ẽ, ẽx̄ = init_defaultcon_mhe(model, He, nx̂, nym, E, ex̄, Ex̂, Fx̂, Gx̂, Jx̂)
113-
nZ̃ = nx̂ + nŵ*He
114-
# dummy values, updated before optimization:
115-
H̃, q̃, p = Hermitian(zeros(NT, nZ̃, nZ̃), :L), zeros(NT, nZ̃), zeros(NT, 1)
116-
= zeros(NT, nZ̃)
117-
X̂, Ym = zeros(NT, nx̂*He), zeros(NT, nym*He)
118-
U, D, Ŵ = zeros(NT, nu*He), zeros(NT, nd*He), zeros(NT, nx̂*He)
119-
x̂arr_old = zeros(NT, nx̂)
120-
P̂arr_old = copy(P̂0)
121-
Nk = [0]
122-
estim = new{NT, SM, JM}(
123-
model, optim, con,
124-
Z̃, lastu0, x̂,
125-
He,
126-
i_ym, nx̂, nym, nyu, nxs,
127-
As, Cs_u, Cs_y, nint_u, nint_ym,
128-
Â, B̂u, Ĉ, B̂d, D̂d,
129-
Ẽ, F, G, J, ẽx̄, fx̄,
130-
H̃, q̃, p,
131-
P̂0, Q̂, R̂, invP̄, invQ̂_He, invR̂_He,
132-
M̂,
133-
X̂, Ym, U, D, Ŵ,
134-
x̂arr_old, P̂arr_old, Nk
135-
)
136-
init_optimization!(estim, model, optim)
137-
return estim
138-
end
139-
end
1+
include("mhe/construct.jl")
2+
include("mhe/execute.jl")
1403

1414
function Base.show(io::IO, estim::MovingHorizonEstimator)
1425
nu, nd = estim.model.nu, estim.model.nd
@@ -148,7 +11,7 @@ function Base.show(io::IO, estim::MovingHorizonEstimator)
14811
print_estim_dim(io, estim, n)
14912
end
15013

151-
"Print the overall dimensions of the state estimator `estim` with left padding `n`."
14+
"Print the overall dimensions of the MHE `estim` with left padding `n`."
15215
function print_estim_dim(io::IO, estim::MovingHorizonEstimator, n)
15316
nu, nd = estim.model.nu, estim.model.nd
15417
nx̂, nym, nyu = estim.nx̂, estim.nym, estim.nyu
@@ -159,7 +22,4 @@ function print_estim_dim(io::IO, estim::MovingHorizonEstimator, n)
15922
println(io, "$(lpad(nym, n)) measured outputs ym ($(sum(estim.nint_ym)) integrating states)")
16023
println(io, "$(lpad(nyu, n)) unmeasured outputs yu")
16124
print(io, "$(lpad(nd, n)) measured disturbances d")
162-
end
163-
164-
include("mhe/construct.jl")
165-
include("mhe/execute.jl")
25+
end

0 commit comments

Comments
 (0)