From f4e87f7af65b163a82c390731efcc983f4ff6605 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Mon, 3 Aug 2026 14:39:13 -0400 Subject: [PATCH 1/6] bench: removing useless unit tests on MHE There is no real computations involved in `updatestate!` with `direct=true`, or `preparestate!` with `direct=false`. These results are just visual noise, let's remove them. --- benchmark/2_bench_state_estim.jl | 28 ---------------------------- 1 file changed, 28 deletions(-) diff --git a/benchmark/2_bench_state_estim.jl b/benchmark/2_bench_state_estim.jl index 25fee68d2..a1c7a213a 100644 --- a/benchmark/2_bench_state_estim.jl +++ b/benchmark/2_bench_state_estim.jl @@ -156,28 +156,11 @@ UNIT_ESTIM["MovingHorizonEstimator"]["preparestate!"]["LinModel"]["Current form" 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"] = - @benchmarkable( - preparestate!($mhe_lin_pred, $y, $d), - samples=samples, evals=evals, seconds=seconds, - ) UNIT_ESTIM["MovingHorizonEstimator"]["updatestate!"]["LinModel"]["Prediction form"] = @benchmarkable( updatestate!($mhe_lin_pred, $u, $y, $d), @@ -195,17 +178,6 @@ 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"] = - @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), - samples=samples, evals=evals, seconds=seconds, - ) UNIT_ESTIM["MovingHorizonEstimator"]["updatestate!"]["NonLinModel"]["Prediction form"] = @benchmarkable( updatestate!($mhe_nonlin_pred, $u, $y, $d), From beb8fb6ae05ee2e8d7e30f6c71c9600b7ad53233 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Mon, 3 Aug 2026 14:52:52 -0400 Subject: [PATCH 2/6] bench: new unit tests for MHE with `MultipleShooting` --- benchmark/2_bench_state_estim.jl | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/benchmark/2_bench_state_estim.jl b/benchmark/2_bench_state_estim.jl index a1c7a213a..be597378e 100644 --- a/benchmark/2_bench_state_estim.jl +++ b/benchmark/2_bench_state_estim.jl @@ -146,9 +146,11 @@ 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"] = @@ -161,6 +163,11 @@ UNIT_ESTIM["MovingHorizonEstimator"]["preparestate!"]["LinModel"]["Constant arr. preparestate!($mhe_lin_skf, $y, $d), samples=samples, evals=evals, seconds=seconds, ) +UNIT_ESTIM["MovingHorizonEstimator"]["preparestate!"]["LinModel"]["MultipleShooting"] = + @benchmarkable( + preparestate!($mhe_lin_ms, $y, $d), + samples=samples, evals=evals, seconds=seconds, + ) UNIT_ESTIM["MovingHorizonEstimator"]["updatestate!"]["LinModel"]["Prediction form"] = @benchmarkable( updatestate!($mhe_lin_pred, $u, $y, $d), @@ -178,6 +185,11 @@ UNIT_ESTIM["MovingHorizonEstimator"]["preparestate!"]["NonLinModel"]["Current fo preparestate!($mhe_nonlin_curr, $y, $d), samples=samples, evals=evals, seconds=seconds, ) +UNIT_ESTIM["MovingHorizonEstimator"]["preparestate!"]["NonLinModel"]["MultipleShooting"] = + @benchmarkable( + preparestate!($mhe_nonlin_ms, $y, $d), + samples=samples, evals=evals, seconds=seconds, + ) UNIT_ESTIM["MovingHorizonEstimator"]["updatestate!"]["NonLinModel"]["Prediction form"] = @benchmarkable( updatestate!($mhe_nonlin_pred, $u, $y, $d), From 9eb3c8bd26d23442c4db033fbc0df943ab277aa8 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Mon, 3 Aug 2026 14:55:19 -0400 Subject: [PATCH 3/6] bench: MHE et MS case study on CSTR --- benchmark/2_bench_state_estim.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/benchmark/2_bench_state_estim.jl b/benchmark/2_bench_state_estim.jl index be597378e..80d9c6e13 100644 --- a/benchmark/2_bench_state_estim.jl +++ b/benchmark/2_bench_state_estim.jl @@ -238,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) @@ -261,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 From ab39b38e14765842bfa3300210a064ad3bc48349 Mon Sep 17 00:00:00 2001 From: franckgaga Date: Tue, 4 Aug 2026 08:16:28 -0400 Subject: [PATCH 4/6] doc: adding `transcription` to seconde MHE constructor --- src/estimator/mhe/construct.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/estimator/mhe/construct.jl b/src/estimator/mhe/construct.jl index dab45443f..bdb17ff27 100644 --- a/src/estimator/mhe/construct.jl +++ b/src/estimator/mhe/construct.jl @@ -572,6 +572,7 @@ end gc!=(_,_,_,_,_,_,_,_,_,_,_) -> nothing, gc=gc!, nc=0, + transcription=SingleShooting(), optim=default_optim_mhe(model, nc), gradient=AutoForwardDiff(), jacobian=AutoForwardDiff(), From e5ec04ff6816e70298edf66846b0dff836ca6d8f Mon Sep 17 00:00:00 2001 From: franckgaga Date: Tue, 4 Aug 2026 11:24:20 -0400 Subject: [PATCH 5/6] bench: new `MultipleShooting` test for the MHE on the inverted pendulum --- benchmark/2_bench_state_estim.jl | 79 ++++++++++++++++++++++---------- 1 file changed, 55 insertions(+), 24 deletions(-) diff --git a/benchmark/2_bench_state_estim.jl b/benchmark/2_bench_state_estim.jl index 80d9c6e13..ca5273d67 100644 --- a/benchmark/2_bench_state_estim.jl +++ b/benchmark/2_bench_state_estim.jl @@ -294,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 @@ -311,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 @@ -335,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 @@ -345,29 +367,38 @@ 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) 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), @@ -377,4 +408,4 @@ 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() - ) \ No newline at end of file + ) From b1e5e68124d5b7789a97efd31c1b19daee0df19d Mon Sep 17 00:00:00 2001 From: franckgaga Date: Tue, 4 Aug 2026 11:35:01 -0400 Subject: [PATCH 6/6] bench: `UnoSolver` on MHE with inverted pendulum Using the `funnelsqp` preset. It is very efficient compared to `Ipopt` (about 6 times faster) --- benchmark/2_bench_state_estim.jl | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/benchmark/2_bench_state_estim.jl b/benchmark/2_bench_state_estim.jl index ca5273d67..0829ca803 100644 --- a/benchmark/2_bench_state_estim.jl +++ b/benchmark/2_bench_state_estim.jl @@ -368,6 +368,25 @@ 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) +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( @@ -409,3 +428,13 @@ CASE_ESTIM["Pendulum"]["MovingHorizonEstimator"]["MadNLP"]["Prediction form (Hes 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() + )