Skip to content

Commit 2890943

Browse files
franckgagafranckgaga
committed
Merge branch 'main' into extension
2 parents 58c9604 + 3b14da7 commit 2890943

File tree

7 files changed

+115
-12
lines changed

7 files changed

+115
-12
lines changed

‎Project.toml‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
name = "ModelPredictiveControl"
22
uuid = "61f9bdb8-6ae4-484a-811f-bbf86720c31c"
3-
version = "1.13.1"
3+
version = "1.13.3"
44
authors = ["Francis Gagnon"]
55

66
[deps]
@@ -20,6 +20,7 @@ RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
2020
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
2121
SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
2222
SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35"
23+
StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
2324

2425
[weakdeps]
2526
LinearMPC = "82e1c212-e1a2-49d2-b26a-a31d6968e3bd"
@@ -30,7 +31,7 @@ LinearMPCext = "LinearMPC"
3031
[compat]
3132
ControlSystemsBase = "1.18.2"
3233
DAQP = "0.6, 0.7.1"
33-
DifferentiationInterface = "0.6.45, 0.7"
34+
DifferentiationInterface = "0.7.11"
3435
Documenter = "1"
3536
FiniteDiff = "2"
3637
ForwardDiff = "0.10, 1"
@@ -48,6 +49,7 @@ RecipesBase = "1"
4849
SparseArrays = "1.10"
4950
SparseConnectivityTracer = "0.6.13, 1"
5051
SparseMatrixColorings = "0.4.14"
52+
StableRNGs = "1.0.4"
5153
TestItemRunner = "1"
5254
TestItems = "1"
5355
julia = "1.10"

‎README.md‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -51,8 +51,7 @@ Our goal is controlling the first output $y_1$, but the second one $y_2$ should
5151
35:
5252

5353
```julia
54-
mhe = MovingHorizonEstimator(model, He=10)
55-
mpc = LinMPC(mhe, Mwt=[1, 0], Nwt=[0.1])
54+
mpc = LinMPC(model, Mwt=[1, 0], Nwt=[0.1])
5655
mpc = setconstraint!(mpc, ymax=[Inf, 35])
5756
```
5857

‎benchmark/3_bench_predictive_control.jl‎

Lines changed: 69 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -308,6 +308,12 @@ nmpc_ipopt_tc = NonLinMPC(estim; Hp, Hc, Mwt, Nwt, Cwt, optim, transcription)
308308
nmpc_ipopt_tc = setconstraint!(nmpc_ipopt_tc; umin, umax)
309309
JuMP.unset_time_limit_sec(nmpc_ipopt_tc.optim)
310310

311+
optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false)
312+
transcription, hessian = TrapezoidalCollocation(), true
313+
nmpc_ipopt_tc_hess = NonLinMPC(estim; Hp, Hc, Mwt, Nwt, Cwt, optim, transcription, hessian)
314+
nmpc_ipopt_tc_hess = setconstraint!(nmpc_ipopt_tc_hess; umin, umax)
315+
JuMP.unset_time_limit_sec(nmpc_ipopt_tc_hess.optim)
316+
311317
optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false)
312318
transcription = TrapezoidalCollocation(f_threads=true)
313319
nmpc_ipopt_tct = NonLinMPC(estim; Hp, Hc, Mwt, Nwt, Cwt, optim, transcription)
@@ -371,6 +377,11 @@ CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["Ipopt"]["TrapezoidalCollocatio
371377
sim!($nmpc_ipopt_tc, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false),
372378
samples=samples, evals=evals, seconds=seconds
373379
)
380+
CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["Ipopt"]["TrapezoidalCollocation (Hessian)"] =
381+
@benchmarkable(
382+
sim!($nmpc_ipopt_tc_hess, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false),
383+
samples=samples, evals=evals, seconds=seconds
384+
)
374385
CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["Ipopt"]["TrapezoidalCollocation (threaded)"] =
375386
@benchmarkable(
376387
sim!($nmpc_ipopt_tct, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0, progress=false),
@@ -406,18 +417,36 @@ empc_ipopt_ss = NonLinMPC(estim2; Hp, Hc, Nwt, Mwt=Mwt2, Cwt, JE, Ewt, optim, tr
406417
empc_ipopt_ss = setconstraint!(empc_ipopt_ss; umin, umax)
407418
JuMP.unset_time_limit_sec(empc_ipopt_ss.optim)
408419

420+
optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false)
421+
transcription, hessian = SingleShooting(), true
422+
empc_ipopt_ss_hess = NonLinMPC(estim2; Hp, Hc, Nwt, Mwt=Mwt2, Cwt, JE, Ewt, optim, transcription, hessian, p)
423+
empc_ipopt_ss_hess = setconstraint!(empc_ipopt_ss_hess; umin, umax)
424+
JuMP.unset_time_limit_sec(empc_ipopt_ss_hess.optim)
425+
409426
optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false)
410427
transcription = MultipleShooting()
411428
empc_ipopt_ms = NonLinMPC(estim2; Hp, Hc, Nwt, Mwt=Mwt2, Cwt, JE, Ewt, optim, transcription, p)
412429
empc_ipopt_ms = setconstraint!(empc_ipopt_ms; umin, umax)
413430
JuMP.unset_time_limit_sec(empc_ipopt_ms.optim)
414431

432+
optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false)
433+
transcription, hessian = MultipleShooting(), true
434+
empc_ipopt_ms_hess = NonLinMPC(estim2; Hp, Hc, Nwt, Mwt=Mwt2, Cwt, JE, Ewt, optim, transcription, hessian, p)
435+
empc_ipopt_ms_hess = setconstraint!(empc_ipopt_ms_hess; umin, umax)
436+
JuMP.unset_time_limit_sec(empc_ipopt_ms_hess.optim)
437+
415438
optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false)
416439
transcription = TrapezoidalCollocation()
417440
empc_ipopt_tc = NonLinMPC(estim2; Hp, Hc, Nwt, Mwt=Mwt2, Cwt, JE, Ewt, optim, transcription, p)
418441
empc_ipopt_tc = setconstraint!(empc_ipopt_tc; umin, umax)
419442
JuMP.unset_time_limit_sec(empc_ipopt_tc.optim)
420443

444+
optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false)
445+
transcription, hessian = TrapezoidalCollocation(), true
446+
empc_ipopt_tc_hess = NonLinMPC(estim2; Hp, Hc, Nwt, Mwt=Mwt2, Cwt, JE, Ewt, optim, transcription, hessian, p)
447+
empc_ipopt_tc_hess = setconstraint!(empc_ipopt_tc_hess; umin, umax)
448+
JuMP.unset_time_limit_sec(empc_ipopt_tc_hess.optim)
449+
421450
optim = JuMP.Model(MadNLP.Optimizer, add_bridges=false)
422451
transcription = SingleShooting()
423452
empc_madnlp_ss = NonLinMPC(estim2; Hp, Hc, Nwt, Mwt=Mwt2, Cwt, JE, Ewt, optim, transcription, p)
@@ -432,16 +461,31 @@ CASE_MPC["Pendulum"]["NonLinMPC"]["Economic"]["Ipopt"]["SingleShooting"] =
432461
sim!($empc_ipopt_ss, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0, progress=false),
433462
samples=samples, evals=evals, seconds=seconds
434463
)
464+
CASE_MPC["Pendulum"]["NonLinMPC"]["Economic"]["Ipopt"]["SingleShooting (Hessian)"] =
465+
@benchmarkable(
466+
sim!($empc_ipopt_ss_hess, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0, progress=false),
467+
samples=samples, evals=evals, seconds=seconds
468+
)
435469
CASE_MPC["Pendulum"]["NonLinMPC"]["Economic"]["Ipopt"]["MultipleShooting"] =
436470
@benchmarkable(
437471
sim!($empc_ipopt_ms, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0, progress=false),
438472
samples=samples, evals=evals, seconds=seconds
439473
)
474+
CASE_MPC["Pendulum"]["NonLinMPC"]["Economic"]["Ipopt"]["MultipleShooting (Hessian)"] =
475+
@benchmarkable(
476+
sim!($empc_ipopt_ms_hess, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0, progress=false),
477+
samples=samples, evals=evals, seconds=seconds
478+
)
440479
CASE_MPC["Pendulum"]["NonLinMPC"]["Economic"]["Ipopt"]["TrapezoidalCollocation"] =
441480
@benchmarkable(
442481
sim!($empc_ipopt_tc, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0, progress=false),
443482
samples=samples, evals=evals, seconds=seconds
444483
)
484+
CASE_MPC["Pendulum"]["NonLinMPC"]["Economic"]["Ipopt"]["TrapezoidalCollocation (Hessian)"] =
485+
@benchmarkable(
486+
sim!($empc_ipopt_tc_hess, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0, progress=false),
487+
samples=samples, evals=evals, seconds=seconds
488+
)
445489
CASE_MPC["Pendulum"]["NonLinMPC"]["Economic"]["MadNLP"]["SingleShooting"] =
446490
@benchmarkable(
447491
sim!($empc_madnlp_ss, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0, progress=false),
@@ -480,6 +524,14 @@ nmpc2_ipopt_ms = NonLinMPC(estim2;
480524
nmpc2_ipopt_ms = setconstraint!(nmpc2_ipopt_ms; umin, umax)
481525
JuMP.unset_time_limit_sec(nmpc2_ipopt_ms.optim)
482526

527+
optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false)
528+
transcription, hessian = MultipleShooting(), true
529+
nmpc2_ipopt_ms_hess = NonLinMPC(estim2;
530+
Hp, Hc, Nwt=Nwt, Mwt=[0.5, 0], Cwt, gc!, nc, p=Pmax, optim, transcription, hessian
531+
)
532+
nmpc2_ipopt_ms_hess = setconstraint!(nmpc2_ipopt_ms_hess; umin, umax)
533+
JuMP.unset_time_limit_sec(nmpc2_ipopt_ms_hess.optim)
534+
483535
optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false)
484536
transcription = TrapezoidalCollocation()
485537
nmpc2_ipopt_tc = NonLinMPC(estim2;
@@ -488,9 +540,13 @@ nmpc2_ipopt_tc = NonLinMPC(estim2;
488540
nmpc2_ipopt_tc = setconstraint!(nmpc2_ipopt_tc; umin, umax)
489541
JuMP.unset_time_limit_sec(nmpc2_ipopt_tc.optim)
490542

491-
# TODO: test custom constraints with MadNLP and SingleShooting, see comment above.
492-
# TODO: test custom constraints with MadNLP and MultipleShooting, see comment above.
493-
# TODO: test custom constraints with MadNLP and TrapezoidalCollocation, see comment above.
543+
optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false)
544+
transcription, hessian = TrapezoidalCollocation(), true
545+
nmpc2_ipopt_tc_hess = NonLinMPC(estim2;
546+
Hp, Hc, Nwt=Nwt, Mwt=[0.5, 0], Cwt, gc!, nc, p=Pmax, optim, transcription, hessian
547+
)
548+
nmpc2_ipopt_tc = setconstraint!(nmpc2_ipopt_tc_hess; umin, umax)
549+
JuMP.unset_time_limit_sec(nmpc2_ipopt_tc_hess.optim)
494550

495551
samples, evals, seconds = 100, 1, 15*60
496552
CASE_MPC["Pendulum"]["NonLinMPC"]["Custom constraints"]["Ipopt"]["SingleShooting"] =
@@ -503,11 +559,21 @@ CASE_MPC["Pendulum"]["NonLinMPC"]["Custom constraints"]["Ipopt"]["MultipleShooti
503559
sim!($nmpc2_ipopt_ms, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0, progress=false),
504560
samples=samples, evals=evals, seconds=seconds
505561
)
562+
CASE_MPC["Pendulum"]["NonLinMPC"]["Custom constraints"]["Ipopt"]["MultipleShooting (Hessian)"] =
563+
@benchmarkable(
564+
sim!($nmpc2_ipopt_ms_hess, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0, progress=false),
565+
samples=samples, evals=evals, seconds=seconds
566+
)
506567
CASE_MPC["Pendulum"]["NonLinMPC"]["Custom constraints"]["Ipopt"]["TrapezoidalCollocation"] =
507568
@benchmarkable(
508569
sim!($nmpc2_ipopt_tc, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0, progress=false),
509570
samples=samples, evals=evals, seconds=seconds
510571
)
572+
CASE_MPC["Pendulum"]["NonLinMPC"]["Custom constraints"]["Ipopt"]["TrapezoidalCollocation (Hessian)"] =
573+
@benchmarkable(
574+
sim!($nmpc2_ipopt_tc_hess, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0, progress=false),
575+
samples=samples, evals=evals, seconds=seconds
576+
)
511577

512578
# ----------------- Case study: Pendulum successive linearization -------------------------
513579
linmodel = linearize(model, x=[0, 0], u=[0])

‎src/ModelPredictiveControl.jl‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@ module ModelPredictiveControl
33
using PrecompileTools
44
using LinearAlgebra, SparseArrays
55
using Random: randn
6+
using StableRNGs: StableRNG
67

78
using RecipesBase
89

@@ -14,6 +15,8 @@ using DifferentiationInterface: hessian!, value_gradient_and_hessian!, prepare_h
1415
using DifferentiationInterface: Constant, Cache
1516
using SparseConnectivityTracer: TracerSparsityDetector
1617
using SparseMatrixColorings: GreedyColoringAlgorithm, sparsity_pattern
18+
using SparseMatrixColorings: NaturalOrder, LargestFirst, SmallestLast
19+
using SparseMatrixColorings: IncidenceDegree, DynamicLargestFirst, RandomOrder
1720

1821
import ProgressLogging
1922

‎src/controller/nonlinmpc.jl‎

Lines changed: 14 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@ const DEFAULT_NONLINMPC_JACDENSE = AutoForwardDiff()
55
const DEFAULT_NONLINMPC_JACSPARSE = AutoSparse(
66
AutoForwardDiff();
77
sparsity_detector=TracerSparsityDetector(),
8-
coloring_algorithm=GreedyColoringAlgorithm(),
8+
coloring_algorithm=GreedyColoringAlgorithm(ALL_COLORING_ORDERS, postprocessing=true),
99
)
1010
const DEFAULT_NONLINMPC_HESSIAN = DEFAULT_NONLINMPC_JACSPARSE
1111

@@ -291,10 +291,21 @@ NonLinMPC controller with a sample time Ts = 10.0 s:
291291
AutoSparse(
292292
AutoForwardDiff();
293293
sparsity_detector = TracerSparsityDetector(),
294-
coloring_algorithm = GreedyColoringAlgorithm()
294+
coloring_algorithm = GreedyColoringAlgorithm(
295+
(
296+
NaturalOrder(),
297+
LargestFirst(),
298+
SmallestLast(),
299+
IncidenceDegree(),
300+
DynamicLargestFirst(),
301+
RandomOrder(StableRNG(0), 0)
302+
),
303+
postprocessing = true
304+
)
295305
)
296306
```
297-
This is also the sparse backend selected for the Hessian of the Lagrangian function if
307+
that is, it will test many coloring orders at preparation and keep the best. This is
308+
also the sparse backend selected for the Hessian of the Lagrangian function if
298309
`oracle=true` and `hessian=true`, which is the second exception. Second order
299310
derivatives are only supported with `oracle=true` option.
300311

‎src/estimator/mhe/construct.jl‎

Lines changed: 14 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@ const DEFAULT_NONLINMHE_JACOBIAN = AutoForwardDiff()
55
const DEFAULT_NONLINMHE_HESSIAN = AutoSparse(
66
AutoForwardDiff();
77
sparsity_detector=TracerSparsityDetector(),
8-
coloring_algorithm=GreedyColoringAlgorithm(),
8+
coloring_algorithm=GreedyColoringAlgorithm(ALL_COLORING_ORDERS, postprocessing=true),
99
)
1010

1111
@doc raw"""
@@ -382,9 +382,21 @@ MovingHorizonEstimator estimator with a sample time Ts = 10.0 s:
382382
AutoSparse(
383383
AutoForwardDiff();
384384
sparsity_detector = TracerSparsityDetector(),
385-
coloring_algorithm = GreedyColoringAlgorithm()
385+
coloring_algorithm = GreedyColoringAlgorithm(
386+
(
387+
NaturalOrder(),
388+
LargestFirst(),
389+
SmallestLast(),
390+
IncidenceDegree(),
391+
DynamicLargestFirst(),
392+
RandomOrder(StableRNG(0), 0)
393+
),
394+
postprocessing = true
395+
)
386396
)
387397
```
398+
that is, it will test many coloring orders at preparation and keep the best.
399+
388400
The slack variable ``ε`` relaxes the constraints if enabled, see [`setconstraint!`](@ref).
389401
It is disabled by default for the MHE (from `Cwt=Inf`) but it should be activated for
390402
problems with two or more types of bounds, to ensure feasibility (e.g. on the estimated

‎src/general.jl‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,16 @@ const DEFAULT_LWT = 0.0
66
const DEFAULT_CWT = 1e5
77
const DEFAULT_EWT = 0.0
88

9+
"All deterministic algorithms for matrix coloring order in `SparseMatrixColoring.jl`."
10+
const ALL_COLORING_ORDERS = (
11+
NaturalOrder(),
12+
LargestFirst(),
13+
SmallestLast(),
14+
IncidenceDegree(),
15+
DynamicLargestFirst(),
16+
RandomOrder(StableRNG(0), 0)
17+
)
18+
919
"Termination status that means 'no solution available'."
1020
const ERROR_STATUSES = (
1121
JuMP.INFEASIBLE, JuMP.DUAL_INFEASIBLE, JuMP.LOCALLY_INFEASIBLE,

0 commit comments

Comments
 (0)