Merge pull request #277 from JuliaControl/mhe_hessian

added: `hessian` keyword argument for `MovingHorizonEstimator`

Author: franckgaga

Project.toml

@@ -1,6 +1,6 @@
 name = "ModelPredictiveControl"
 uuid = "61f9bdb8-6ae4-484a-811f-bbf86720c31c"
-version = "1.12.0"
+version = "1.13.0"
 authors = ["Francis Gagnon"]
 
 [deps]

benchmark/2_bench_state_estim.jl

@@ -305,6 +305,13 @@ 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
+)
+mhe_pendulum_ipopt_currh = setconstraint!(mhe_pendulum_ipopt_currh; v̂min, v̂max)
+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(
@@ -314,6 +321,14 @@ 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)
 
+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(MadNLP.Optimizer, add_bridges=false)
 direct = true
 mhe_pendulum_madnlp_curr = MovingHorizonEstimator(
@@ -338,11 +353,21 @@ CASE_ESTIM["Pendulum"]["MovingHorizonEstimator"]["Ipopt"]["Current form"] =
         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
     )
+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
+    )
 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
     )
+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
+    )
 CASE_ESTIM["Pendulum"]["MovingHorizonEstimator"]["MadNLP"]["Current form"] =
     @benchmarkable(
         sim!($mhe_pendulum_madnlp_curr, $N, $u; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false),

benchmark/3_bench_predictive_control.jl

@@ -318,7 +318,13 @@ optim = JuMP.Model(MadNLP.Optimizer, add_bridges=false)
 transcription = SingleShooting()
 nmpc_madnlp_ss = NonLinMPC(estim; Hp, Hc, Mwt, Nwt, Cwt, optim, transcription)
 nmpc_madnlp_ss = setconstraint!(nmpc_madnlp_ss; umin, umax)
-JuMP.unset_time_limit_sec(nmpc_madnlp_ss.optim) 
+JuMP.unset_time_limit_sec(nmpc_madnlp_ss.optim)
+
+optim = JuMP.Model(()->UnoSolver.Optimizer(preset="filtersqp"), add_bridges=false)
+transcription, hessian, oracle = MultipleShooting(), true, true
+nmpc_uno_ms_hess = NonLinMPC(estim; Hp, Hc, Mwt, Nwt, Cwt, optim, transcription, hessian, oracle)
+nmpc_uno_ms_hess = setconstraint!(nmpc_uno_ms_hess; umin, umax)
+JuMP.unset_time_limit_sec(nmpc_uno_ms_hess.optim)
 
 # TODO: does not work well with MadNLP and MultipleShooting or TrapezoidalCollocation, 
 # figure out why. Current theory: 
@@ -375,6 +381,11 @@ CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["MadNLP"]["SingleShooting"] =
         sim!($nmpc_madnlp_ss, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false),
         samples=samples, evals=evals, seconds=seconds
     )
+CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["Uno"]["MultipleShooting (Hessian)"] = 
+    @benchmarkable(
+        sim!($nmpc_uno_ms_hess, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false),
+        samples=samples, evals=evals, seconds=seconds
+    )    
 
 # ----------------- Case study: Pendulum economic --------------------------------
 model2, p = pendulum_model2, pendulum_p2

benchmark/Project.toml

@@ -13,3 +13,4 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
 SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35"
+UnoSolver = "1baa60ac-02f7-4b39-a7a8-2f4f58486b05"

benchmark/benchmarks.jl

@@ -1,6 +1,6 @@
 using BenchmarkTools
 using ModelPredictiveControl, ControlSystemsBase, LinearAlgebra
-using JuMP, OSQP, DAQP, Ipopt, MadNLP
+using JuMP, OSQP, DAQP, Ipopt, MadNLP, UnoSolver
 
 const SUITE = BenchmarkGroup(["ModelPredictiveControl"])
 

docs/src/internals/state_estim.md

@@ -18,7 +18,8 @@ ModelPredictiveControl.relaxX̂
 ModelPredictiveControl.relaxŴ
 ModelPredictiveControl.relaxV̂
 ModelPredictiveControl.init_matconstraint_mhe
-ModelPredictiveControl.get_nonlinops(::MovingHorizonEstimator, ::ModelPredictiveControl.GenericModel)
+ModelPredictiveControl.get_nonlinobj_op(::MovingHorizonEstimator, ::ModelPredictiveControl.GenericModel)
+ModelPredictiveControl.get_nonlincon_oracle(::MovingHorizonEstimator, ::ModelPredictiveControl.GenericModel)
 ```
 
 ## Augmented Model