Skip to content

Commit 3b9fe6d

Browse files
committed
added: trunc_cov function for MHE
The function will truncate the inverted covariance matrices only when Nk<He. There is also two methods, one for `AbstractMatrix` and one for `Diagonal`. It indices some allocations when Nk<He but at least `dot` in the objective function will be always computer on `Hermitian` or `Diagonal`.
1 parent 6fdbce5 commit 3b9fe6d

File tree

1 file changed

+28
-2
lines changed

1 file changed

+28
-2
lines changed

src/estimator/mhe/execute.jl

Lines changed: 28 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -266,7 +266,8 @@ function initpred!(estim::MovingHorizonEstimator, model::LinModel)
266266
# --- update H̃, q̃ and p vectors for quadratic optimization ---
267267
ẼZ̃ = @views [estim.ẽx̄[:, 1:nZ̃]; estim.Ẽ[1:nYm, 1:nZ̃]]
268268
FZ̃ = @views [estim.fx̄; estim.F[1:nYm]]
269-
invQ̂_Nk, invR̂_Nk = @views invQ̂_He[1:nŴ, 1:nŴ], invR̂_He[1:nYm, 1:nYm]
269+
invQ̂_Nk = trunc_cov(estim.cov.invQ̂_He, estim.nx̂, Nk, estim.He)
270+
invR̂_Nk = trunc_cov(estim.cov.invR̂_He, estim.nym, Nk, estim.He)
270271
M_Nk = [invP̄ zeros(nx̂, nYm); zeros(nYm, nx̂) invR̂_Nk]
271272
Ñ_Nk = [fill(C, nε, nε) zeros(nε, nx̂+nŴ); zeros(nx̂, nε+nx̂+nŴ); zeros(nŴ, nε+nx̂) invQ̂_Nk]
272273
M_Nk_ẼZ̃ = M_Nk*ẼZ̃
@@ -477,6 +478,28 @@ function set_warmstart_mhe!(V̂, X̂0, estim::MovingHorizonEstimator{NT}, Z̃var
477478
return Z̃s
478479
end
479480

481+
"Truncate the inverse covariance `invA_He` to the window size `Nk` if `Nk < He`."
482+
function trunc_cov(invA_He::Hermitian{<:Real, <:AbstractMatrix}, n, Nk, He)
483+
if Nk < He
484+
nA = Nk*n
485+
# avoid views since allocations only when Nk < He and we want type-stability:
486+
return Hermitian(invA_He[1:nA, 1:nA], :L)
487+
else
488+
return invA_He
489+
end
490+
end
491+
function trunc_cov(
492+
invA_He::Hermitian{NT, Diagonal{NT, Vector{NT}}}, n, Nk, He
493+
) where NT <:Real
494+
if Nk < He
495+
nA = Nk*n
496+
# avoid views since allocations only when Nk < He and we want type-stability:
497+
return Hermitian(Diagonal(invA_He.data.diag[1:nA]), :L)
498+
else
499+
return invA_He
500+
end
501+
end
502+
480503
"Correct the covariance estimate at arrival using `covestim` [`StateEstimator`](@ref)."
481504
function correct_cov!(estim::MovingHorizonEstimator)
482505
nym, nd = estim.nym, estim.model.nd
@@ -573,13 +596,16 @@ function obj_nonlinprog!(
573596
nε, Nk = estim.nε, estim.Nk[]
574597
nYm, nŴ, nx̂, invP̄ = Nk*estim.nym, Nk*estim.nx̂, estim.nx̂, estim.cov.invP̄
575598
nx̃ =+ nx̂
576-
invQ̂_Nk, invR̂_Nk = @views estim.cov.invQ̂_He[1:nŴ, 1:nŴ], estim.cov.invR̂_He[1:nYm, 1:nYm]
599+
invQ̂_Nk = trunc_cov(estim.cov.invQ̂_He, estim.nx̂, Nk, estim.He)
600+
invR̂_Nk = trunc_cov(estim.cov.invR̂_He, estim.nym, Nk, estim.He)
577601
x̂0arr, Ŵ, V̂ = @views Z̃[nx̃-nx̂+1:nx̃], Z̃[nx̃+1:nx̃+nŴ], V̂[1:nYm]
578602
x̄ .= estim.x̂0arr_old .- x̂0arr
579603
= 0 ? estim.C*Z̃[begin]^2 : zero(NT)
580604
return dot(x̄, invP̄, x̄) + dot(Ŵ, invQ̂_Nk, Ŵ) + dot(V̂, invR̂_Nk, V̂) +
581605
end
582606

607+
608+
583609
@doc raw"""
584610
predict_mhe!(V̂, X̂0, _, _, _, estim::MovingHorizonEstimator, model::LinModel, Z̃) -> V̂, X̂0
585611

0 commit comments

Comments
 (0)