Merge branch 'main' into Ipopt_VectorNonlinearOracle

Author: franckgaga

.github/dependabot.yml

@@ -0,0 +1,7 @@
+# https://docs.github.com/github/administering-a-repository/configuration-options-for-dependency-updates
+version: 2
+updates:
+  - package-ecosystem: "github-actions"
+    directory: "/" # Location of package manifests
+    schedule:
+      interval: "weekly"

.github/workflows/CI.yml

@@ -34,7 +34,7 @@ jobs:
       - name: Set JULIA_DEBUG environment variable if applicable
         if: ${{ runner.debug == '1' }}
         run: echo "JULIA_DEBUG=ModelPredictiveControl" >> $GITHUB_ENV
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v5
       - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}

.github/workflows/DocCleanup.yml

@@ -9,7 +9,7 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - name: Checkout gh-pages branch
-        uses: actions/checkout@v2
+        uses: actions/checkout@v5
         with:
           ref: gh-pages
 

.github/workflows/benchmark.yml

@@ -2,11 +2,18 @@ name: Benchmark
 on:
   pull_request_target:
     branches: [ main ]
+    types: [labeled, opened, synchronize, reopened]
+concurrency:
+  # Skip intermediate builds: always.
+  # Cancel intermediate builds: only if it is a pull request build.
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }}
 permissions:
   pull-requests: write    # needed to post comments
 jobs:
   bench:
     runs-on: ubuntu-latest
+    if: ${{ contains(github.event.pull_request.labels.*.name, 'benchmark') }}
     steps:
       - uses: MilesCranmer/AirspeedVelocity.jl@action-v1
         with:

.github/workflows/documentation.yml

@@ -13,8 +13,8 @@ jobs:
       contents: write
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v5
+      - uses: julia-actions/setup-julia@v2
         with:
           version: '1'
       - name: Install dependencies

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelPredictiveControl"
 uuid = "61f9bdb8-6ae4-484a-811f-bbf86720c31c"
 authors = ["Francis Gagnon"]
-version = "1.8.2"
+version = "1.11.0"
 
 [deps]
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
@@ -21,7 +21,7 @@ SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
 SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35"
 
 [compat]
-ControlSystemsBase = "1.9"
+ControlSystemsBase = "1.18.2"
 DAQP = "0.6, 0.7.1"
 DifferentiationInterface = "0.6.45, 0.7"
 Documenter = "1"

README.md

@@ -4,6 +4,7 @@
 [![codecov](https://codecov.io/gh/JuliaControl/ModelPredictiveControl.jl/branch/main/graph/badge.svg?token=K4V0L113M4)](https://codecov.io/gh/JuliaControl/ModelPredictiveControl.jl)
 [![doc-stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://JuliaControl.github.io/ModelPredictiveControl.jl/stable)
 [![doc-dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://JuliaControl.github.io/ModelPredictiveControl.jl/dev)
+[![arXiv](https://img.shields.io/badge/arXiv-2411.09764-b31b1b.svg)](https://arxiv.org/abs/2411.09764)
 
 An open source [model predictive control](https://en.wikipedia.org/wiki/Model_predictive_control)
 package for Julia.
@@ -115,6 +116,7 @@ for more detailed examples.
 - supported transcription methods of the optimization problem:
   - direct single shooting
   - direct multiple shooting
+  - trapezoidal collocation
 - additional information about the optimum to ease troubleshooting
 - real-time control loop features:
   - implementations that carefully limits the allocations

benchmark/0_bench_setup.jl

@@ -16,6 +16,11 @@ linmodel = setop!(LinModel(sys, Ts, i_d=[3]), uop=[10, 50], yop=[50, 30], dop=[5
 nonlinmodel = NonLinModel(f_lin!, h_lin!, Ts, 2, 4, 2, 1, p=linmodel, solver=nothing)
 nonlinmodel = setop!(nonlinmodel, uop=[10, 50], yop=[50, 30], dop=[5])
 u, d, y = [10, 50], [5], [50, 30]
+nonlinmodel_c = NonLinModel(
+    (ẋ,x,u,d,_) -> ẋ .= -0.001x .+ u .+ d , 
+    (y,x,d,_) -> y .= x .+ 0.001d, 500, 1, 1, 1, 1
+)
+u_c, d_c, y_c = [1], [0], [0]
 
 G = [ tf(1.90, [18, 1]) tf(1.90, [18, 1]);
       tf(-0.74,[8, 1])  tf(0.74, [8, 1]) ]

benchmark/3_bench_predictive_control.jl

@@ -50,6 +50,10 @@ nmpc_nonlin_ms = NonLinMPC(
     nonlinmodel, transcription=MultipleShooting(),
     Mwt=[1, 1], Nwt=[0.1, 0.1], Lwt=[0.1, 0.1], Hp=10    
 )
+nmpc_nonlin_tc = NonLinMPC(
+    nonlinmodel_c, transcription=TrapezoidalCollocation(),
+    Mwt=[1], Nwt=[0.1], Lwt=[0.1], Hp=10    
+)
 
 samples, evals, seconds = 10000, 1, 60
 UNIT_MPC["NonLinMPC"]["moveinput!"]["LinModel"]["SingleShooting"] =
@@ -76,6 +80,12 @@ UNIT_MPC["NonLinMPC"]["moveinput!"]["NonLinModel"]["MultipleShooting"] =
         setup=preparestate!($nmpc_nonlin_ms, $y, $d),
         samples=samples, evals=evals, seconds=seconds
     )
+UNIT_MPC["NonLinMPC"]["moveinput!"]["NonLinModel"]["TrapezoidalCollocation"] =
+    @benchmarkable(
+        moveinput!($nmpc_nonlin_tc, $y_c, $d_c),
+        setup=preparestate!($nmpc_nonlin_tc, $y_c, $d_c),
+        samples=samples, evals=evals, seconds=seconds
+    )
 
 ## ----------------------------------------------------------------------------------------
 ## ---------------------- CASE STUDIES ----------------------------------------------------
@@ -230,7 +240,6 @@ CASE_MPC["CSTR"]["LinMPC"]["With feedforward"]["Ipopt"]["MultipleShooting"] =
         samples=samples, evals=evals
     )
 
-
 # ----------------- Case study: Pendulum noneconomic -----------------------------
 model, p = pendulum_model, pendulum_p
 σQ = [0.1, 1.0]; σR=[5.0]; nint_u=[1]; σQint_u=[0.1]
@@ -255,13 +264,32 @@ nmpc_ipopt_ms = NonLinMPC(estim; Hp, Hc, Mwt, Nwt, Cwt, optim, transcription)
 nmpc_ipopt_ms = setconstraint!(nmpc_ipopt_ms; umin, umax)
 JuMP.unset_time_limit_sec(nmpc_ipopt_ms.optim)
 
+optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false)
+transcription = MultipleShooting(f_threads=true)
+nmpc_ipopt_mst = NonLinMPC(estim; Hp, Hc, Mwt, Nwt, Cwt, optim, transcription)
+nmpc_ipopt_mst = setconstraint!(nmpc_ipopt_mst; umin, umax)
+JuMP.unset_time_limit_sec(nmpc_ipopt_mst.optim)
+
+optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false)
+transcription = TrapezoidalCollocation()
+nmpc_ipopt_tc = NonLinMPC(estim; Hp, Hc, Mwt, Nwt, Cwt, optim, transcription)
+nmpc_ipopt_tc = setconstraint!(nmpc_ipopt_tc; umin, umax)
+JuMP.unset_time_limit_sec(nmpc_ipopt_tc.optim)
+
+optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false)
+transcription = TrapezoidalCollocation(f_threads=true)
+nmpc_ipopt_tct = NonLinMPC(estim; Hp, Hc, Mwt, Nwt, Cwt, optim, transcription)
+nmpc_ipopt_tct = setconstraint!(nmpc_ipopt_tct; umin, umax)
+JuMP.unset_time_limit_sec(nmpc_ipopt_tct.optim)
+
 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) 
 
-# TODO: does not work well with MadNLP and MultipleShooting, figure out why. Current theory: 
+# TODO: does not work well with MadNLP and MultipleShooting or TrapezoidalCollocation, 
+# figure out why. Current theory: 
 # MadNLP LBFGS approximation is less robust than Ipopt version. Re-test when exact Hessians
 # will be supported in ModelPredictiveControl.jl. The following attributes kinda work with 
 # the MadNLP LBFGS approximation but super slow (~1000 times slower than Ipopt):
@@ -285,6 +313,21 @@ CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["Ipopt"]["MultipleShooting"] =
         sim!($nmpc_ipopt_ms, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0),
         samples=samples, evals=evals, seconds=seconds
     )
+CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["Ipopt"]["MultipleShooting (threaded)"] =
+    @benchmarkable(
+        sim!($nmpc_ipopt_mst, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0),
+        samples=samples, evals=evals, seconds=seconds
+    )
+CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["Ipopt"]["TrapezoidalCollocation"] =
+    @benchmarkable(
+        sim!($nmpc_ipopt_tc, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0),
+        samples=samples, evals=evals, seconds=seconds
+    )
+CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["Ipopt"]["TrapezoidalCollocation (threaded)"] =
+    @benchmarkable(
+        sim!($nmpc_ipopt_tct, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0),
+        samples=samples, evals=evals, seconds=seconds
+    )
 CASE_MPC["Pendulum"]["NonLinMPC"]["Noneconomic"]["MadNLP"]["SingleShooting"] = 
     @benchmarkable(
         sim!($nmpc_madnlp_ss, $N, $ry; plant=$plant, x_0=$x_0, x̂_0=$x̂_0),
@@ -316,13 +359,19 @@ empc_ipopt_ms = NonLinMPC(estim2; Hp, Hc, Nwt, Mwt=Mwt2, Cwt, JE, Ewt, optim, tr
 empc_ipopt_ms = setconstraint!(empc_ipopt_ms; umin, umax)
 JuMP.unset_time_limit_sec(empc_ipopt_ms.optim)
 
+optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false)
+transcription = TrapezoidalCollocation()
+empc_ipopt_tc = NonLinMPC(estim2; Hp, Hc, Nwt, Mwt=Mwt2, Cwt, JE, Ewt, optim, transcription, p)
+empc_ipopt_tc = setconstraint!(empc_ipopt_tc; umin, umax)
+JuMP.unset_time_limit_sec(empc_ipopt_tc.optim)
+
 optim = JuMP.Model(MadNLP.Optimizer, add_bridges=false)
 transcription = SingleShooting()
 empc_madnlp_ss = NonLinMPC(estim2; Hp, Hc, Nwt, Mwt=Mwt2, Cwt, JE, Ewt, optim, transcription, p)
 empc_madnlp_ss = setconstraint!(empc_madnlp_ss; umin, umax)
 JuMP.unset_time_limit_sec(empc_madnlp_ss.optim)
 
-# TODO: test EMPC with MadNLP and MultipleShooting, see comment above.
+# TODO: test EMPC with MadNLP and MultipleShooting and TrapezoidalCollocation, see comment above.
 
 samples, evals, seconds = 100, 1, 15*60
 CASE_MPC["Pendulum"]["NonLinMPC"]["Economic"]["Ipopt"]["SingleShooting"] = 
@@ -335,6 +384,11 @@ CASE_MPC["Pendulum"]["NonLinMPC"]["Economic"]["Ipopt"]["MultipleShooting"] =
         sim!($empc_ipopt_ms, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0),
         samples=samples, evals=evals, seconds=seconds
     )
+CASE_MPC["Pendulum"]["NonLinMPC"]["Economic"]["Ipopt"]["TrapezoidalCollocation"] =
+    @benchmarkable(
+        sim!($empc_ipopt_tc, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0),
+        samples=samples, evals=evals, seconds=seconds
+    )
 CASE_MPC["Pendulum"]["NonLinMPC"]["Economic"]["MadNLP"]["SingleShooting"] = 
     @benchmarkable(
         sim!($empc_madnlp_ss, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0),
@@ -373,8 +427,17 @@ nmpc2_ipopt_ms = NonLinMPC(estim2;
 nmpc2_ipopt_ms = setconstraint!(nmpc2_ipopt_ms; umin, umax)
 JuMP.unset_time_limit_sec(nmpc2_ipopt_ms.optim)
 
+optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false)
+transcription = TrapezoidalCollocation()
+nmpc2_ipopt_tc = NonLinMPC(estim2; 
+    Hp, Hc, Nwt=Nwt, Mwt=[0.5, 0], Cwt, gc!, nc, p=Pmax, optim, transcription
+)
+nmpc2_ipopt_tc = setconstraint!(nmpc2_ipopt_tc; umin, umax)
+JuMP.unset_time_limit_sec(nmpc2_ipopt_tc.optim)
+
 # TODO: test custom constraints with MadNLP and SingleShooting, see comment above.
 # TODO: test custom constraints with MadNLP and MultipleShooting, see comment above.
+# TODO: test custom constraints with MadNLP and TrapezoidalCollocation, see comment above.
 
 samples, evals, seconds = 100, 1, 15*60
 CASE_MPC["Pendulum"]["NonLinMPC"]["Custom constraints"]["Ipopt"]["SingleShooting"] = 
@@ -387,6 +450,11 @@ CASE_MPC["Pendulum"]["NonLinMPC"]["Custom constraints"]["Ipopt"]["MultipleShooti
         sim!($nmpc2_ipopt_ms, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0),
         samples=samples, evals=evals, seconds=seconds
     )
+CASE_MPC["Pendulum"]["NonLinMPC"]["Custom constraints"]["Ipopt"]["TrapezoidalCollocation"] =
+    @benchmarkable(
+        sim!($nmpc2_ipopt_tc, $N, $ry; plant=$plant2, x_0=$x_0, x̂_0=$x̂_0),
+        samples=samples, evals=evals, seconds=seconds
+    )
 
 # ----------------- Case study: Pendulum successive linearization -------------------------
 linmodel = linearize(model, x=[0, 0], u=[0])

docs/Project.toml

@@ -16,5 +16,5 @@ Documenter = "1"
 JuMP = "1"
 LinearAlgebra = "1.10"
 Logging = "1.10"
-ModelingToolkit = "9.50 - 9.76"
+ModelingToolkit = "10"
 Plots = "1"