Skip to content

Commit 10dde49

Browse files
committed
added: setmodel! for MHE
1 parent 89200e5 commit 10dde49

File tree

7 files changed

+104
-95
lines changed

7 files changed

+104
-95
lines changed

src/controller/execute.jl

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -564,12 +564,12 @@ function setmodel_controller!(mpc::PredictiveController, model::LinModel, x̂op_
564564
mpc.Yop[(1+ny*i):(ny+ny*i)] .= model.yop
565565
mpc.Dop[(1+nd*i):(nd+nd*i)] .= model.dop
566566
end
567-
con.U0min .-= mpc.Uop # convert U0 to U with the new operating point
568-
con.U0max .-= mpc.Uop # convert U0 to U with the new operating point
569-
con.Y0min .-= mpc.Yop # convert Y0 to Y with the new operating point
570-
con.Y0max .-= mpc.Yop # convert Y0 to Y with the new operating point
571-
con.x̂0min .-= estim.x̂op # convert x̂0 to with the new operating point
572-
con.x̂0max .-= estim.x̂op # convert x̂0 to with the new operating point
567+
con.U0min .-= mpc.Uop # convert U to U0 with the new operating point
568+
con.U0max .-= mpc.Uop # convert U to U0 with the new operating point
569+
con.Y0min .-= mpc.Yop # convert Y to Y0 with the new operating point
570+
con.Y0max .-= mpc.Yop # convert Y to Y0 with the new operating point
571+
con.x̂0min .-= estim.x̂op # convert to x̂0 with the new operating point
572+
con.x̂0max .-= estim.x̂op # convert to x̂0 with the new operating point
573573
con.A_Ymin .= A_Ymin
574574
con.A_Ymax .= A_Ymax
575575
con.A_x̂min .= A_x̂min

src/estimator/execute.jl

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -273,6 +273,9 @@ julia> setmodel!(kf, LinModel(ss(0.42, 0.5, 1, 0, 4.0))); kf.model.A
273273
"""
274274
function setmodel!(estim::StateEstimator, model::LinModel)
275275
validate_model(estim.model, model)
276+
uop_old = copy(estim.model.uop)
277+
yop_old = copy(estim.model.yop)
278+
dop_old = copy(estim.model.dop)
276279
# --- update model matrices and its operating points ---
277280
estim.model.A .= model.A
278281
estim.model.Bu .= model.Bu
@@ -286,7 +289,7 @@ function setmodel!(estim::StateEstimator, model::LinModel)
286289
estim.model.dop .= model.dop
287290
estim.model.xop .= model.xop
288291
estim.model.fop .= model.fop
289-
setmodel_estimator!(estim, model)
292+
setmodel_estimator!(estim, model, uop_old, yop_old, dop_old)
290293
return estim
291294
end
292295

@@ -301,7 +304,7 @@ end
301304
validate_model(::SimModel, ::SimModel) = error("Only LinModel is supported for setmodel!")
302305

303306
"Update the augmented model matrices of `estim` by default."
304-
function setmodel_estimator!(estim::StateEstimator, model::LinModel)
307+
function setmodel_estimator!(estim::StateEstimator, model::LinModel, _ , _ , _)
305308
As, Cs_u, Cs_y = estim.As, estim.Cs_u, estim.Cs_y
306309
Â, B̂u, Ĉ, B̂d, D̂d, x̂op, f̂op = augment_model(model, As, Cs_u, Cs_y, verify_obsv=false)
307310
# --- update augmented state-space matrices ---

src/estimator/internal_model.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -203,7 +203,7 @@ function init_internalmodel(As, Bs, Cs, Ds)
203203
end
204204

205205
"Update similar values for [`InternalModel`](@ref) estimator."
206-
function setmodel_estimator!(estim::InternalModel, model::LinModel)
206+
function setmodel_estimator!(estim::InternalModel, model::LinModel, _ , _ , _)
207207
Â, B̂u, Ĉ, B̂d, D̂d, x̂op, f̂op = matrices_internalmodel(model)
208208
# --- update augmented state-space matrices ---
209209
estim. .=
@@ -215,7 +215,7 @@ function setmodel_estimator!(estim::InternalModel, model::LinModel)
215215
estim.x̂0 .+= estim.x̂op # convert x̂0 to x̂ with the old operating point
216216
estim.x̂op .= x̂op
217217
estim.f̂op .= f̂op
218-
estim.x̂0 .-= estim.x̂op # convert x̂0 to with the new operating point
218+
estim.x̂0 .-= estim.x̂op # convert to x̂0 with the new operating point
219219
end
220220

221221
@doc raw"""

src/estimator/kalman.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -165,7 +165,7 @@ function SteadyKalmanFilter(model::SM, i_ym, nint_u, nint_ym, Q̂, R̂) where {N
165165
end
166166

167167
"Throw an error if `setmodel!` is called on a SteadyKalmanFilter"
168-
function setmodel_estimator!(estim::SteadyKalmanFilter, model::LinModel)
168+
function setmodel_estimator!(::SteadyKalmanFilter, ::LinModel, _ , _ , _)
169169
error("SteadyKalmanFilter does not support setmodel! (use KalmanFilter instead)")
170170
end
171171

src/estimator/luenberger.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -122,4 +122,4 @@ function update_estimate!(estim::Luenberger, u, ym, d=empty(estim.x̂0))
122122
end
123123

124124
"Throw an error if `setmodel!` is called on `Luenberger` observer."
125-
setmodel_estimator!(estim::Luenberger, ::LinModel) = error("Luenberger does not support setmodel!")
125+
setmodel_estimator!(::Luenberger,::LinModel,_,_,_) = error("Luenberger does not support setmodel!")

src/estimator/mhe/construct.jl

Lines changed: 9 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -92,7 +92,6 @@ struct MovingHorizonEstimator{
9292
invQ̂_He::Hermitian{NT, Matrix{NT}}
9393
invR̂_He::Hermitian{NT, Matrix{NT}}
9494
C::NT
95-
::Matrix{NT}
9695
X̂op::Vector{NT}
9796
X̂0 ::Vector{NT}
9897
Y0m::Vector{NT}
@@ -121,10 +120,11 @@ struct MovingHorizonEstimator{
121120
invP̄ = Hermitian(inv(P̂_0), :L)
122121
invQ̂_He = Hermitian(repeatdiag(inv(Q̂), He), :L)
123122
invR̂_He = Hermitian(repeatdiag(inv(R̂), He), :L)
124-
= zeros(NT, nx̂, nym)
125-
E, F, G, J, B, ex̄, fx̄, Ex̂, Fx̂, Gx̂, Jx̂, Bx̂ = init_predmat_mhe(
123+
E, G, J, B, ex̄, Ex̂, Gx̂, Jx̂, Bx̂ = init_predmat_mhe(
126124
model, He, i_ym, Â, B̂u, Ĉ, B̂d, D̂d, x̂op, f̂op
127125
)
126+
# dummy values (updated just before optimization):
127+
F, fx̄, Fx̂ = zeros(NT, nym*He), zeros(NT, nx̂), zeros(NT, nx̂*He)
128128
con, Ẽ, ẽx̄ = init_defaultcon_mhe(model, He, Cwt, nx̂, nym, E, ex̄, Ex̂, Fx̂, Gx̂, Jx̂, Bx̂)
129129
nZ̃ = size(Ẽ, 2)
130130
# dummy values, updated before optimization:
@@ -147,9 +147,7 @@ struct MovingHorizonEstimator{
147147
Ẽ, F, G, J, B, ẽx̄, fx̄,
148148
H̃, q̃, p,
149149
P̂_0, Q̂, R̂, invP̄, invQ̂_He, invR̂_He, Cwt,
150-
M̂,
151-
X̂op,
152-
X̂0, Y0m, U0, D0, Ŵ,
150+
X̂op, X̂0, Y0m, U0, D0, Ŵ,
153151
x̂0arr_old, P̂arr_old, Nk
154152
)
155153
init_optimization!(estim, model, optim)
@@ -558,8 +556,8 @@ function init_matconstraint_mhe(::LinModel{NT},
558556
if isempty(args)
559557
A = zeros(NT, length(i_b), 0)
560558
else
561-
A_x̂min, A_x̂max, A_X̂min, A_X̂max, A_Ŵmin, A_Ŵmax, A_V̂min, A_V̂max = args
562-
A = [A_x̂min; A_x̂max; A_X̂min; A_X̂max; A_Ŵmin; A_Ŵmax; A_V̂min; A_V̂max]
559+
A_x̃min, A_x̃max, A_X̂min, A_X̂max, A_Ŵmin, A_Ŵmax, A_V̂min, A_V̂max = args
560+
A = [A_x̃min; A_x̃max; A_X̂min; A_X̂max; A_Ŵmin; A_Ŵmax; A_V̂min; A_V̂max]
563561
end
564562
return i_b, i_g, A
565563
end
@@ -811,7 +809,7 @@ end
811809
@doc raw"""
812810
init_predmat_mhe(
813811
model::LinModel, He, i_ym, Â, B̂u, Ĉ, B̂d, D̂d, x̂op, f̂op
814-
) -> E, F, G, J, ex̄, fx̄, Ex̂, Fx̂, Gx̂, Jx̂
812+
) -> E, G, J, B, ex̄, Ex̂, Gx̂, Jx̂, Bx̂
815813
816814
Construct the MHE prediction matrices for [`LinModel`](@ref) `model`.
817815
@@ -987,11 +985,7 @@ function init_predmat_mhe(
987985
coef_Bx̂[iRow,:] = getpower(Âpow_csum, j)
988986
end
989987
Bx̂ = coef_Bx̂*f̂_op_n_x̂op
990-
# --- F vectors ---
991-
F = zeros(NT, nym*He) # dummy F vector (updated just before optimization)
992-
fx̄ = zeros(NT, nx̂) # dummy fx̄ vector (updated just before optimization)
993-
Fx̂ = zeros(NT, nx̂*He) # dummy Fx̂ vector (updated just before optimization)
994-
return E, F, G, J, B, ex̄, fx̄, Ex̂, Fx̂, Gx̂, Jx̂, Bx̂
988+
return E, G, J, B, ex̄, Ex̂, Gx̂, Jx̂, Bx̂
995989
end
996990

997991
"Return empty matrices if `model` is not a [`LinModel`](@ref), except for `ex̄`."
@@ -1009,10 +1003,7 @@ function init_predmat_mhe(
10091003
Jx̂ = zeros(NT, 0, model.nd*He)
10101004
B = zeros(NT, nym*He)
10111005
Bx̂ = zeros(NT, nx̂*He)
1012-
F = zeros(NT, nym*He)
1013-
fx̄ = zeros(NT, nx̂)
1014-
Fx̂ = zeros(NT, nx̂*He)
1015-
return E, F, G, J, B, ex̄, fx̄, Ex̂, Fx̂, Gx̂, Jx̂, Bx̂
1006+
return E, G, J, B, ex̄, Ex̂, Gx̂, Jx̂, Bx̂
10161007
end
10171008

10181009
"""

src/estimator/mhe/execute.jl

Lines changed: 80 additions & 65 deletions
Original file line numberDiff line numberDiff line change
@@ -461,80 +461,95 @@ end
461461

462462

463463
"Update the augmented model matrices of `estim` by default."
464-
function setmodel_estimator!(estim::MovingHorizonEstimator, model::LinModel)
465-
As, Cs_u, Cs_y = estim.As, estim.Cs_u, estim. Cs_y
466-
Â, B̂u, Ĉ, B̂d, D̂d = augment_model(model, As, Cs_u, Cs_y, verify_obsv=false)
464+
function setmodel_estimator!(
465+
estim::MovingHorizonEstimator, model::LinModel, uop_old, yop_old, dop_old
466+
)
467+
con = estim.con
468+
nx̂, nym, nu, nd, He = estim.nx̂, estim.nym, model.nu, model.nd, estim.He
469+
= isinf(estim.C) ? 0 : 1
470+
As, Cs_u, Cs_y = estim.As, estim.Cs_u, estim.Cs_y
471+
Â, B̂u, Ĉ, B̂d, D̂d, x̂op, f̂op = augment_model(model, As, Cs_u, Cs_y, verify_obsv=false)
472+
# --- update augmented state-space matrices ---
467473
estim. .=
468474
estim.B̂u .= B̂u
469475
estim.Ĉ .=
470476
estim.B̂d .= B̂d
471477
estim.D̂d .= D̂d
472-
# TODO: re-construct the MHE prediction matrices here:
473-
error("setmodel! for MovingHorizonEstimator is not implemented yet.")
474-
475-
476-
477-
478-
479-
480-
481-
482-
483-
484-
485-
# LINMPC:
478+
# --- update state estimate and its operating points ---
479+
x̂op_old = copy(estim.x̂op)
480+
X̂op_old = copy(estim.X̂op)
481+
estim.x̂0 .+= estim.x̂op # convert x̂0 to x̂ with the old operating point
482+
estim.x̂op .= x̂op
483+
estim.f̂op .= f̂op
484+
estim.x̂0 .-= estim.x̂op # convert x̂ to x̂0 with the new operating point
486485
# --- predictions matrices ---
487-
E, G, J, K, V, B, ex̂, gx̂, jx̂, kx̂, vx̂, bx̂ = init_predmat(estim, model, Hp, Hc)
488-
A_Ymin, A_Ymax, Ẽ = relaxŶ(model, mpc.C, con.C_ymin, con.C_ymax, E)
489-
A_x̂min, A_x̂max, ẽx̂ = relaxterminal(model, mpc.C, con.c_x̂min, con.c_x̂max, ex̂)
490-
mpc.Ẽ .=
491-
mpc.G .= G
492-
mpc.J .= J
493-
mpc.K .= K
494-
mpc.V .= V
495-
mpc.B .= B
486+
E, G, J, B, _ , Ex̂, Gx̂, Jx̂, Bx̂ = init_predmat_mhe(
487+
model, He, estim.i_ym, Â, B̂u, Ĉ, B̂d, D̂d, x̂op, f̂op
488+
)
489+
A_X̂min, A_X̂max, Ẽx̂ = relaxX̂(model, estim.C, con.C_x̂min, con.C_x̂max, Ex̂)
490+
A_V̂min, A_V̂max, Ẽ = relaxV̂(model, estim.C, con.C_v̂min, con.C_v̂max, E)
491+
estim. .=
492+
estim.G .= G
493+
estim.J .= J
494+
estim.B .= B
496495
# --- linear inequality constraints ---
497-
con.ẽx̂ .= ẽx̂
498-
con.gx̂ .= gx̂
499-
con.jx̂ .= jx̂
500-
con.kx̂ .= kx̂
501-
con.vx̂ .= vx̂
502-
con.bx̂ .= bx̂
503-
con.U0min .+= mpc.Uop # convert U0 to U with the old operating point
504-
con.U0max .+= mpc.Uop # convert U0 to U with the old operating point
505-
con.Y0min .+= mpc.Yop # convert Y0 to Y with the old operating point
506-
con.Y0max .+= mpc.Yop # convert Y0 to Y with the old operating point
507-
con.x̂0min .+= x̂op_old # convert x̂0 to x̂ with the old operating point
508-
con.x̂0max .+= x̂op_old # convert x̂0 to x̂ with the old operating point
509-
# --- operating points ---
510-
for i in 0:Hp-1
511-
mpc.Uop[(1+nu*i):(nu+nu*i)] .= model.uop
512-
mpc.Yop[(1+ny*i):(ny+ny*i)] .= model.yop
513-
mpc.Dop[(1+nd*i):(nd+nd*i)] .= model.dop
496+
con.Ẽx̂ .= Ẽx̂
497+
con.Gx̂ .= Gx̂
498+
con.Jx̂ .= Jx̂
499+
con.Bx̂ .= Bx̂
500+
# convert x̃0 to x̃ with the old operating point:
501+
con.x̃0min[end-nx̂+1:end] .-= x̂op_old
502+
con.x̃0max[end-nx̂+1:end] .-= x̂op_old
503+
# convert X̂0 to X̂ with the old operating point:
504+
con.X̂0min .-= X̂op_old
505+
con.X̂0max .-= X̂op_old
506+
for i in 0:He-1
507+
estim.X̂op[(1+nx̂*i):(nx̂+nx̂*i)] .= estim.x̂op
514508
end
515-
con.U0min .-= mpc.Uop # convert U0 to U with the new operating point
516-
con.U0max .-= mpc.Uop # convert U0 to U with the new operating point
517-
con.Y0min .-= mpc.Yop # convert Y0 to Y with the new operating point
518-
con.Y0max .-= mpc.Yop # convert Y0 to Y with the new operating point
519-
con.x̂0min .-= estim.x̂op # convert x̂0 to x̂ with the new operating point
520-
con.x̂0max .-= estim.x̂op # convert x̂0 to x̂ with the new operating point
521-
con.A_Ymin .= A_Ymin
522-
con.A_Ymax .= A_Ymax
523-
con.A_x̂min .= A_x̂min
524-
con.A_x̂max .= A_x̂max
525-
nUandΔŨ = length(con.U0min) + length(con.U0max) + length(con.ΔŨmin) + length(con.ΔŨmax)
526-
con.A[nUandΔŨ+1:end, :] = [con.A_Ymin; con.A_Ymax; con.A_x̂min; con.A_x̂max]
509+
# convert x̃ to x̃0 with the new operating point:
510+
con.x̃0min[end-nx̂+1:end] .+= estim.x̂op
511+
con.x̃0max[end-nx̂+1:end] .+= estim.x̂op
512+
# convert X̂ to X̂0 with the new operating point:
513+
con.X̂0min .+= estim.X̂op
514+
con.X̂0max .+= estim.X̂op
515+
con.A
516+
con.A_X̂min .= A_X̂min
517+
con.A_X̂max .= A_X̂max
518+
con.A_V̂min .= A_V̂min
519+
con.A_V̂max .= A_V̂max
520+
nx̃ = length(con.x̃0min) + length(con.x̃0max)
521+
nX̂ = length(con.X̂0min) + length(con.X̂0max)
522+
con.A[nx̃+1:nx̃+nX̂,:] .= [con.A_X̂min; con.A_X̂max]
523+
nŴ = length(con.Ŵmin) + length(con.Ŵmax)
524+
con.A[nx̃+nX̂+nŴ+1:end,:] .= [con.A_V̂min; con.A_V̂max]
527525
A = con.A[con.i_b, :]
528526
b = con.b[con.i_b]
529-
ΔŨvar::Vector{JuMP.VariableRef} = optim[:ΔŨvar]
530-
JuMP.delete(optim, optim[:linconstraint])
531-
JuMP.unregister(optim, :linconstraint)
532-
@constraint(optim, linconstraint, A*ΔŨvar .≤ b)
533-
# --- quadratic programming Hessian matrix ---
534-
= init_quadprog(model, mpc.Ẽ, mpc.S̃, mpc.M_Hp, mpc.Ñ_Hc, mpc.L_Hp)
535-
mpc.H̃ .=
536-
set_objective_hessian!(mpc, ΔŨvar)
537-
538-
527+
Z̃var::Vector{JuMP.VariableRef} = estim.optim[:Z̃var]
528+
JuMP.delete(estim.optim, estim.optim[:linconstraint])
529+
JuMP.unregister(estim.optim, :linconstraint)
530+
@constraint(estim.optim, linconstraint, A*Z̃var .≤ b)
531+
# --- data windows ---
532+
for i in 1:He
533+
# convert x̂0 to x̂ with the old operating point:
534+
estim.X̂0[(1+nx̂*(i-1)):(nx̂*i)] .-= x̂op_old
535+
# convert ym0 to ym with the old operating point:
536+
estim.Y0m[(1+nym*(i-1)):(nym*i)] .-= @views yop_old[estim.i_ym]
537+
# convert u0 to u with the old operating point:
538+
estim.U0[(1+nu*(i-1)):(nu*i)] .-= uop_old
539+
# convert d0 to d with the old operating point:
540+
estim.D0[(1+nd*(i-1)):(nd*i)] .-= dop_old
541+
# convert x̂ to x̂0 with the new operating point:
542+
estim.X̂0[(1+nx̂*(i-1)):(nx̂*i)] .+= x̂op
543+
# convert ym to y0m with the new operating point:
544+
estim.Y0m[(1+nym*(i-1)):(nym*i)] .+= @views yop_old[estim.i_ym]
545+
# convert u to u0 with the new operating point:
546+
estim.U0[(1+nu*(i-1)):(nu*i)] .+= uop_old
547+
# convert d to d0 with the new operating point:
548+
estim.D0[(1+nd*(i-1)):(nd*i)] .+= dop_old
549+
end
550+
estim.Z̃[nϵ+1:+nx̂] .-= x̂op_old
551+
estim.Z̃[nϵ+1:+nx̂] .+= x̂op
552+
estim.x̂0arr_old .-= x̂op_old
553+
estim.x̂0arr_old .+= x̂op
539554
return nothing
540555
end

0 commit comments

Comments
 (0)