added: NonLinModelDAE for differential and algebraic equations + sim! support

Project.toml

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

docs/make.jl

@@ -13,6 +13,7 @@ links = InterLinks(
     "JuMP" => "https://jump.dev/JuMP.jl/stable/objects.inv",
     "MathOptInterface" => "https://jump.dev/MathOptInterface.jl/stable/objects.inv",
     "DifferentiationInterface" => "https://juliadiff.org/DifferentiationInterface.jl/DifferentiationInterface/stable/objects.inv",
+    "SparseConnectivityTracer" => "https://adrianhill.de/SparseConnectivityTracer.jl/stable/objects.inv",
     "ForwardDiff" => "https://juliadiff.org/ForwardDiff.jl/stable/objects.inv",
     "LowLevelParticleFilters" => "https://baggepinnen.github.io/LowLevelParticleFilters.jl/stable/objects.inv",
     "LinearMPC" => "https://darnstrom.github.io/LinearMPC.jl/stable/objects.inv",

docs/src/internals/misc.md

@@ -4,6 +4,13 @@
 Pages = ["misc.md"]
 ```
 
+## Abstract Types
+
+```@docs
+ModelPredictiveControl.ShootingMethod
+ModelPredictiveControl.CollocationMethod
+```
+
 ## Orthogonal Collocation
 
 ```@docs

docs/src/internals/sim_model.md

@@ -4,15 +4,28 @@
 Pages = ["sim_model.md"]
 ```
 
+## Abstract Types
+
+```@docs
+ModelPredictiveControl.SimModelODE
+ModelPredictiveControl.SimModelDAE
+```
+
+## Model Construction
+
+```@docs
+ModelPredictiveControl.init_defectmat_dae
+```
+
 ## State-Space Functions
 
 ```@docs
 ModelPredictiveControl.f!
 ModelPredictiveControl.h!
 ```
 
-## Steady-State Calculation
+## Init State
 
 ```@docs
-ModelPredictiveControl.steadystate!
+ModelPredictiveControl.initstate_core!
 ```

docs/src/public/sim_model.md

@@ -37,6 +37,12 @@ LinModel
 NonLinModel
 ```
 
+## NonLinModelDAE
+
+```@docs
+NonLinModelDAE
+```
+
 ## Set Variable Names
 
 ```@docs

src/ModelPredictiveControl.jl

@@ -17,7 +17,7 @@ using DifferentiationInterface: gradient!, value_and_gradient!, prepare_gradient
 using DifferentiationInterface: jacobian!, value_and_jacobian!, prepare_jacobian 
 using DifferentiationInterface: hessian!, value_gradient_and_hessian!, prepare_hessian
 using DifferentiationInterface: Constant, Cache
-using SparseConnectivityTracer: TracerSparsityDetector
+using SparseConnectivityTracer: TracerSparsityDetector, jacobian_sparsity, jacobian_buffer
 using SparseMatrixColorings: GreedyColoringAlgorithm, sparsity_pattern
 using SparseMatrixColorings: NaturalOrder, LargestFirst, SmallestLast
 using SparseMatrixColorings: IncidenceDegree, DynamicLargestFirst, RandomOrder
@@ -41,7 +41,7 @@ import OSQP, Ipopt
 
 import FastGaussQuadrature
 
-export SimModel, LinModel, NonLinModel
+export SimModel, LinModel, NonLinModel, NonLinModelDAE
 export DiffSolver, RungeKutta, ForwardEuler
 export setop!, setname!
 export setstate!, setmodel!, preparestate!, updatestate!, evaloutput, linearize, linearize!
@@ -57,8 +57,8 @@ export TrapezoidalCollocation, OrthogonalCollocation
 export SimResult, getinfo, sim!
 
 include("general.jl")
-include("sim_model.jl")
 include("transcription.jl")
+include("sim_model.jl")
 include("state_estim.jl")
 include("predictive_control.jl")
 include("plot_sim.jl")

src/controller/construct.jl

@@ -94,7 +94,7 @@ end
 
 "Outer constructor to validate and convert weight matrices if necessary."
 function ControllerWeights(
-    model::SimModel{NT}, Hp, Hc, M_Hp, N_Hc, L_Hp, Cwt=Inf, Ewt=0
+    model::SimModelODE{NT}, Hp, Hc, M_Hp, N_Hc, L_Hp, Cwt=Inf, Ewt=0
 ) where {NT<:Real}
     validate_weights(model, Hp, Hc, M_Hp, N_Hc, L_Hp, Cwt, Ewt)
     M_Hp, N_Hc, L_Hp = NT.(M_Hp), NT.(N_Hc), NT.(L_Hp)
@@ -568,7 +568,7 @@ Estimate the default prediction horizon `Hp` for [`LinModel`](@ref).
 """
 default_Hp(model::LinModel) = DEFAULT_HP0 + estimate_delays(model)
 "Throw an error when model is not a [`LinModel`](@ref)."
-function default_Hp(::SimModel)
+function default_Hp(::SimModelODE)
     msg = "Prediction horizon Hp must be explicitly specified if model is not a LinModel."
     throw(ArgumentError(msg))
 end

src/controller/execute.jl

@@ -277,25 +277,25 @@ function initpred!(mpc::PredictiveController, model::LinModel, ry, d, lastu, D̂
 end
 
 @doc raw"""
-    initpred!(mpc::PredictiveController, model::SimModel, ry, d, lastu, D̂, R̂y, R̂u) -> nothing
+    initpred!(mpc::PredictiveController, model::SimModelODE, ry, d, lastu, D̂, R̂y, R̂u) -> nothing
 
 Init `lastu0, ŷ, F, d0, D̂0, D̂e, R̂y, R̂u` vectors when model is not a [`LinModel`](@ref).
 """
-function initpred!(mpc::PredictiveController, model::SimModel, ry, d, lastu, D̂, R̂y, R̂u)
+function initpred!(mpc::PredictiveController, model::SimModelODE, ry, d, lastu, D̂, R̂y, R̂u)
     initpred_common!(mpc, model, ry, d, lastu, D̂, R̂y, R̂u)
     return nothing
 end
 
 """
-    initpred_common!(mpc::PredictiveController, model::SimModel, ry, d, lastu, D̂, R̂y, R̂u) -> F
+    initpred_common!(mpc::PredictiveController, model::SimModelODE, ry, d, lastu, D̂, R̂y, R̂u) -> F
 
-Common computations of `initpred!` for all types of [`SimModel`](@ref).
+Common computations of `initpred!` for all types of [`SimModelODE`](@ref).
 
 Will also init `mpc.F` with 0 values, or with the stochastic predictions `Ŷs` if `mpc.estim`
 is an [`InternalModel`](@ref). The function returns `mpc.F`.
 """
 function initpred_common!(
-    mpc::PredictiveController, model::SimModel, ry, d, lastu, D̂, R̂y, R̂u
+    mpc::PredictiveController, model::SimModelODE, ry, d, lastu, D̂, R̂y, R̂u
 )
     mpc.lastu0 .= lastu .- model.uop
     mul!(mpc.Tu_lastu0, mpc.Tu, mpc.lastu0)
@@ -327,14 +327,14 @@ end
 predictstoch!(Ŷs, ::PredictiveController, ::StateEstimator) = (Ŷs .= 0; nothing)
 
 @doc raw"""
-    linconstraint_custom!(mpc::PredictiveController, model::SimModel)
+    linconstraint_custom!(mpc::PredictiveController, model::SimModelODE)
 
 Init the ``\mathbf{F_w}`` vector for the custom linear inequality constraints.
 
 See [`relaxW`](@ref) for the definition of the vector. The function does nothing if
 `mpc.con.nw < 1`.
 """
-function linconstraint_custom!(mpc::PredictiveController, model::SimModel)
+function linconstraint_custom!(mpc::PredictiveController, model::SimModelODE)
     mpc.con.nw < 1 && return nothing
     ny, nu, nd, buffer = model.ny, model.nu, model.nd, mpc.buffer
     Fw = mpc.con.Fw
@@ -363,7 +363,7 @@ function linconstraint_custom_outputs!(mpc::PredictiveController, model::LinMode
     return nothing
 end
 "Do nothing for other model types."
-linconstraint_custom_outputs!(::PredictiveController, ::SimModel) = nothing
+linconstraint_custom_outputs!(::PredictiveController, ::SimModelODE) = nothing
 
 """
     extended_vectors!(Ue, Ŷe, mpc::PredictiveController, U0, Ŷ0) -> Ue, Ŷe
@@ -506,7 +506,7 @@ end
 
 
 "By default, no need to update the objective function."
-set_objective_linear_coef!(::PredictiveController, ::SimModel, _) = nothing
+set_objective_linear_coef!(::PredictiveController, ::SimModelODE, _) = nothing
 
 "Update the linear coefficients of the quadratic objective with `mpc.q̃` for `LinModel`."
 function set_objective_linear_coef!(mpc::PredictiveController, ::LinModel, Z̃var)
@@ -790,7 +790,7 @@ function setmodel_controller!(mpc::PredictiveController, uop_old, x̂op_old)
 end
 
 "No need to set the objective Hessian by default (only needed for quadratic objective)."
-set_objective_hessian!(::PredictiveController, ::SimModel, _ ) = nothing
+set_objective_hessian!(::PredictiveController, ::SimModelODE, _ ) = nothing
 
 "Set the objective Hessian with `mpc.H̃` if the objective is quadratic."
 function set_objective_hessian!(mpc::PredictiveController, ::LinModel, Z̃var)