Skip to content

bench: adding new MovingHorizonEstimator cases with MultipleShooting #418

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
franckgaga merged 7 commits into from
Aug 4, 2026
+102 −47
Merged

bench: adding new MovingHorizonEstimator cases with MultipleShooting #418

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
f4e87f7
bench: removing useless unit tests on MHE
franckgaga Aug 3, 2026
beb8fb6
bench: new unit tests for MHE with `MultipleShooting`
franckgaga Aug 3, 2026
9eb3c8b
bench: MHE et MS case study on CSTR
franckgaga Aug 3, 2026
9c6b44d
Merge branch 'main' into bench_ms_mhe
franckgaga Aug 3, 2026
ab39b38
doc: adding `transcription` to seconde MHE constructor
franckgaga Aug 4, 2026
e5ec04f
bench: new `MultipleShooting` test for the MHE on the inverted pendulum
franckgaga Aug 4, 2026
b1e5e68
bench: `UnoSolver` on MHE with inverted pendulum
franckgaga Aug 4, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
148 changes: 101 additions & 47 deletions benchmark/2_bench_state_estim.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -146,36 +146,26 @@ UNIT_ESTIM["ExtendedKalmanFilter"]["evaloutput"]["NonLinModel"] =

mhe_lin_curr = MovingHorizonEstimator(linmodel, He=10, direct=true)
mhe_lin_pred = MovingHorizonEstimator(linmodel, He=10, direct=false)
mhe_lin_skf = MovingHorizonEstimator(linmodel, He=10, covestim=SteadyKalmanFilter(linmodel))
mhe_lin_ms = MovingHorizonEstimator(linmodel, He=10, transcription=MultipleShooting())
mhe_lin_skf = MovingHorizonEstimator(linmodel, He=10, covestim=SteadyKalmanFilter(linmodel))
mhe_nonlin_curr = MovingHorizonEstimator(nonlinmodel, He=10, direct=true)
mhe_nonlin_pred = MovingHorizonEstimator(nonlinmodel, He=10, direct=false)
mhe_nonlin_ms = MovingHorizonEstimator(nonlinmodel, He=10, transcription=MultipleShooting())

samples, evals, seconds = 10000, 1, 60
UNIT_ESTIM["MovingHorizonEstimator"]["preparestate!"]["LinModel"]["Current form"] =
@benchmarkable(
preparestate!($mhe_lin_curr, $y, $d),
samples=samples, evals=evals, seconds=seconds,
)
UNIT_ESTIM["MovingHorizonEstimator"]["updatestate!"]["LinModel"]["Current form"] =
@benchmarkable(
updatestate!($mhe_lin_curr, $u, $y, $d),
setup=preparestate!($mhe_lin_curr, $y, $d),
samples=samples, evals=evals, seconds=seconds,
)
UNIT_ESTIM["MovingHorizonEstimator"]["preparestate!"]["LinModel"]["Constant arr. cov."] =
@benchmarkable(
preparestate!($mhe_lin_skf, $y, $d),
samples=samples, evals=evals, seconds=seconds,
)
UNIT_ESTIM["MovingHorizonEstimator"]["updatestate!"]["LinModel"]["Constant arr. cov."] =
@benchmarkable(
updatestate!($mhe_lin_skf, $u, $y, $d),
setup=preparestate!($mhe_lin_skf, $y, $d),
samples=samples, evals=evals, seconds=seconds,
)
UNIT_ESTIM["MovingHorizonEstimator"]["preparestate!"]["LinModel"]["Prediction form"] =
UNIT_ESTIM["MovingHorizonEstimator"]["preparestate!"]["LinModel"]["MultipleShooting"] =
@benchmarkable(
preparestate!($mhe_lin_pred, $y, $d),
preparestate!($mhe_lin_ms, $y, $d),
samples=samples, evals=evals, seconds=seconds,
)
UNIT_ESTIM["MovingHorizonEstimator"]["updatestate!"]["LinModel"]["Prediction form"] =
Expand All @@ -195,15 +185,9 @@ UNIT_ESTIM["MovingHorizonEstimator"]["preparestate!"]["NonLinModel"]["Current fo
preparestate!($mhe_nonlin_curr, $y, $d),
samples=samples, evals=evals, seconds=seconds,
)
UNIT_ESTIM["MovingHorizonEstimator"]["updatestate!"]["NonLinModel"]["Current form"] =
UNIT_ESTIM["MovingHorizonEstimator"]["preparestate!"]["NonLinModel"]["MultipleShooting"] =
@benchmarkable(
updatestate!($mhe_nonlin_curr, $u, $y, $d),
setup=preparestate!($mhe_nonlin_curr, $y, $d),
samples=samples, evals=evals, seconds=seconds,
)
UNIT_ESTIM["MovingHorizonEstimator"]["preparestate!"]["NonLinModel"]["Prediction form"] =
@benchmarkable(
preparestate!($mhe_nonlin_pred, $y, $d),
preparestate!($mhe_nonlin_ms, $y, $d),
samples=samples, evals=evals, seconds=seconds,
)
UNIT_ESTIM["MovingHorizonEstimator"]["updatestate!"]["NonLinModel"]["Prediction form"] =
Expand Down Expand Up @@ -254,6 +238,12 @@ mhe_cstr_osqp_curr = MovingHorizonEstimator(model; He, nint_u, σQint_u, optim,
mhe_cstr_osqp_curr = setconstraint!(mhe_cstr_osqp_curr; v̂min, v̂max)
JuMP.unset_time_limit_sec(mhe_cstr_osqp_curr.optim)

optim = JuMP.Model(OSQP.Optimizer, add_bridges=true)
transcription = MultipleShooting()
mhe_cstr_osqp_ms = MovingHorizonEstimator(model; He, nint_u, σQint_u, optim, transcription)
mhe_cstr_osqp_ms = setconstraint!(mhe_cstr_osqp_ms; v̂min, v̂max)
JuMP.unset_time_limit_sec(mhe_cstr_osqp_ms.optim)

optim = JuMP.Model(OSQP.Optimizer, add_bridges=true)
direct = false
mhe_cstr_osqp_pred = MovingHorizonEstimator(model; He, nint_u, σQint_u, optim, direct)
Expand All @@ -277,6 +267,10 @@ CASE_ESTIM["CSTR"]["MovingHorizonEstimator"]["OSQP"]["Current form"] =
@benchmarkable(test_mhe($mhe_cstr_osqp_curr, $plant);
samples=samples, evals=evals
)
CASE_ESTIM["CSTR"]["MovingHorizonEstimator"]["OSQP"]["MultipleShooting"] =
@benchmarkable(test_mhe($mhe_cstr_osqp_ms, $plant);
samples=samples, evals=evals
)
CASE_ESTIM["CSTR"]["MovingHorizonEstimator"]["OSQP"]["Prediction form"] =
@benchmarkable(test_mhe($mhe_cstr_osqp_pred, $plant);
samples=samples, evals=evals
Expand All @@ -300,15 +294,17 @@ N = 35;

x_0 = [0.1, 0.1]; x̂_0 = [0, 0, 0]; u = [0.5]

### The MHE without exact Hessians does not work well on the inverted pendulum, commenting them:
# optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false)
# direct = true
# mhe_pendulum_ipopt_curr = MovingHorizonEstimator(
# model; He, σQ, σR, nint_u, σQint_u, optim, direct
# )
# mhe_pendulum_ipopt_curr = setconstraint!(mhe_pendulum_ipopt_curr; v̂min, v̂max)
# JuMP.unset_time_limit_sec(mhe_pendulum_ipopt_curr.optim)

optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false)
direct = true
mhe_pendulum_ipopt_curr = MovingHorizonEstimator(
model; He, σQ, σR, nint_u, σQint_u, optim, direct
)
mhe_pendulum_ipopt_curr = setconstraint!(mhe_pendulum_ipopt_curr; v̂min, v̂max)
JuMP.unset_time_limit_sec(mhe_pendulum_ipopt_curr.optim)
JuMP.set_attribute(mhe_pendulum_ipopt_curr.optim, "tol", 1e-7)

hessian = true
mhe_pendulum_ipopt_currh = MovingHorizonEstimator(
model; He, σQ, σR, nint_u, σQint_u, optim, direct, hessian
Expand All @@ -317,21 +313,42 @@ mhe_pendulum_ipopt_currh = setconstraint!(mhe_pendulum_ipopt_currh; v̂min, v̂m
JuMP.unset_time_limit_sec(mhe_pendulum_ipopt_currh.optim)

optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false)
direct = false
mhe_pendulum_ipopt_pred = MovingHorizonEstimator(
model; He, σQ, σR, nint_u, σQint_u, optim, direct
direct = true
hessian = true
transcription = MultipleShooting()
mhe_pendulum_ipopt_currhms = MovingHorizonEstimator(
model; He, σQ, σR, nint_u, σQint_u, optim, direct, hessian, transcription
)
mhe_pendulum_ipopt_pred = setconstraint!(mhe_pendulum_ipopt_pred; v̂min, v̂max)
JuMP.unset_time_limit_sec(mhe_pendulum_ipopt_pred.optim)
JuMP.set_attribute(mhe_pendulum_ipopt_pred.optim, "tol", 1e-7)
mhe_pendulum_ipopt_currhms = setconstraint!(mhe_pendulum_ipopt_currhms; v̂min, v̂max)
JuMP.unset_time_limit_sec(mhe_pendulum_ipopt_currhms.optim)

# ## The MHE without exact Hessians does not work well on the inverted pendulum, commenting them:
# optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false)
# direct = false
# mhe_pendulum_ipopt_pred = MovingHorizonEstimator(
# model; He, σQ, σR, nint_u, σQint_u, optim, direct
# )
# mhe_pendulum_ipopt_pred = setconstraint!(mhe_pendulum_ipopt_pred; v̂min, v̂max)
# JuMP.unset_time_limit_sec(mhe_pendulum_ipopt_pred.optim)

optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false)
direct = false
hessian = true
mhe_pendulum_ipopt_predh = MovingHorizonEstimator(
model; He, σQ, σR, nint_u, σQint_u, optim, direct, hessian
)
mhe_pendulum_ipopt_predh = setconstraint!(mhe_pendulum_ipopt_predh; v̂min, v̂max)
JuMP.unset_time_limit_sec(mhe_pendulum_ipopt_predh.optim)
JuMP.set_attribute(mhe_pendulum_ipopt_predh.optim, "tol", 1e-7)

optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false)
direct = false
hessian = true
transcription = MultipleShooting()
mhe_pendulum_ipopt_predhms = MovingHorizonEstimator(
model; He, σQ, σR, nint_u, σQint_u, optim, direct, hessian, transcription
)
mhe_pendulum_ipopt_predhms = setconstraint!(mhe_pendulum_ipopt_predhms; v̂min, v̂max)
JuMP.unset_time_limit_sec(mhe_pendulum_ipopt_predhms.optim)

optim = JuMP.Model(MadNLP.Optimizer, add_bridges=false)
direct = true
Expand All @@ -341,7 +358,6 @@ mhe_pendulum_madnlp_currh = MovingHorizonEstimator(
)
mhe_pendulum_madnlp_currh = setconstraint!(mhe_pendulum_madnlp_currh; v̂min, v̂max)
JuMP.unset_time_limit_sec(mhe_pendulum_madnlp_currh.optim)
JuMP.set_attribute(mhe_pendulum_madnlp_currh.optim, "tol", 1e-7)

optim = JuMP.Model(MadNLP.Optimizer, add_bridges=false)
direct = false
Expand All @@ -351,29 +367,57 @@ mhe_pendulum_madnlp_predh = MovingHorizonEstimator(
)
mhe_pendulum_madnlp_pred = setconstraint!(mhe_pendulum_madnlp_predh; v̂min, v̂max)
JuMP.unset_time_limit_sec(mhe_pendulum_madnlp_predh.optim)
JuMP.set_attribute(mhe_pendulum_madnlp_predh.optim, "tol", 1e-7)

optim = JuMP.Model(()->UnoSolver.Optimizer(preset="funnelsqp"), add_bridges=false)
direct = true
hessian = true
mhe_pendulum_uno_currh = MovingHorizonEstimator(
model; He, σQ, σR, nint_u, σQint_u, optim, direct, hessian
)
mhe_pendulum_uno_currh = setconstraint!(mhe_pendulum_uno_currh; v̂min, v̂max)
JuMP.unset_time_limit_sec(mhe_pendulum_uno_currh.optim)

optim = JuMP.Model(()->UnoSolver.Optimizer(preset="funnelsqp"), add_bridges=false)
direct = true
hessian = true
transcription = MultipleShooting()
mhe_pendulum_uno_currhms = MovingHorizonEstimator(
model; He, σQ, σR, nint_u, σQint_u, optim, direct, hessian, transcription
)
mhe_pendulum_uno_currhms = setconstraint!(mhe_pendulum_uno_currhms; v̂min, v̂max)
JuMP.unset_time_limit_sec(mhe_pendulum_uno_currhms.optim)

samples, evals, seconds = 25, 1, 15*60
CASE_ESTIM["Pendulum"]["MovingHorizonEstimator"]["Ipopt"]["Current form"] =
@benchmarkable(
sim!($mhe_pendulum_ipopt_curr, $N, $u; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false),
samples=samples, evals=evals, seconds=seconds, setup=GC.gc()
)
# CASE_ESTIM["Pendulum"]["MovingHorizonEstimator"]["Ipopt"]["Current form"] =
# @benchmarkable(
# sim!($mhe_pendulum_ipopt_curr, $N, $u; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false),
# samples=samples, evals=evals, seconds=seconds, setup=GC.gc()
# )
CASE_ESTIM["Pendulum"]["MovingHorizonEstimator"]["Ipopt"]["Current form (Hessian)"] =
@benchmarkable(
sim!($mhe_pendulum_ipopt_currh, $N, $u; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false),
samples=samples, evals=evals, seconds=seconds, setup=GC.gc()
)
CASE_ESTIM["Pendulum"]["MovingHorizonEstimator"]["Ipopt"]["Prediction form"] =
CASE_ESTIM["Pendulum"]["MovingHorizonEstimator"]["Ipopt"]["Current form (Hessian, MultipleShooting)"] =
@benchmarkable(
sim!($mhe_pendulum_ipopt_pred, $N, $u; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false),
sim!($mhe_pendulum_ipopt_currhms, $N, $u; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false),
samples=samples, evals=evals, seconds=seconds, setup=GC.gc()
)
# CASE_ESTIM["Pendulum"]["MovingHorizonEstimator"]["Ipopt"]["Prediction form"] =
# @benchmarkable(
# sim!($mhe_pendulum_ipopt_pred, $N, $u; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false),
# samples=samples, evals=evals, seconds=seconds, setup=GC.gc()
# )
CASE_ESTIM["Pendulum"]["MovingHorizonEstimator"]["Ipopt"]["Prediction form (Hessian)"] =
@benchmarkable(
sim!($mhe_pendulum_ipopt_predh, $N, $u; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false),
samples=samples, evals=evals, seconds=seconds, setup=GC.gc()
)
CASE_ESTIM["Pendulum"]["MovingHorizonEstimator"]["Ipopt"]["Prediction form (Hessian, MultipleShooting)"] =
@benchmarkable(
sim!($mhe_pendulum_ipopt_predhms, $N, $u; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false),
samples=samples, evals=evals, seconds=seconds, setup=GC.gc()
)
CASE_ESTIM["Pendulum"]["MovingHorizonEstimator"]["MadNLP"]["Current form (Hessian)"] =
@benchmarkable(
sim!($mhe_pendulum_madnlp_currh, $N, $u; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false),
Expand All @@ -383,4 +427,14 @@ CASE_ESTIM["Pendulum"]["MovingHorizonEstimator"]["MadNLP"]["Prediction form (Hes
@benchmarkable(
sim!($mhe_pendulum_madnlp_predh, $N, $u; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false),
samples=samples, evals=evals, seconds=seconds, setup=GC.gc()
)
)
CASE_ESTIM["Pendulum"]["MovingHorizonEstimator"]["Uno"]["Current form (Hessian)"] =
@benchmarkable(
sim!($mhe_pendulum_uno_currh, $N, $u; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false),
samples=samples, evals=evals, seconds=seconds, setup=GC.gc()
)
CASE_ESTIM["Pendulum"]["MovingHorizonEstimator"]["Uno"]["Current form (Hessian, MultipleShooting)"] =
@benchmarkable(
sim!($mhe_pendulum_uno_currhms, $N, $u; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false),
samples=samples, evals=evals, seconds=seconds, setup=GC.gc()
)
1 change: 1 addition & 0 deletions src/estimator/mhe/construct.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -572,6 +572,7 @@ end
gc!=(_,_,_,_,_,_,_,_,_,_,_) -> nothing,
gc=gc!,
nc=0,
transcription=SingleShooting(),
optim=default_optim_mhe(model, nc),
gradient=AutoForwardDiff(),
jacobian=AutoForwardDiff(),
Expand Down
Loading