rebase

Author: vyudu

Project.toml

@@ -77,7 +77,7 @@ LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800"
 
 [extensions]
 MTKBifurcationKitExt = "BifurcationKit"
-MTKCasADiDynamicOptExt = ["CasADi", "DataInterpolations", "DiffEqDevTools"]
+MTKCasADiDynamicOptExt = ["CasADi", "DataInterpolations"]
 MTKChainRulesCoreExt = "ChainRulesCore"
 MTKDeepDiffsExt = "DeepDiffs"
 MTKFMIExt = "FMI"

ext/MTKCasADiDynamicOptExt.jl

@@ -1,10 +1,10 @@
 module MTKCasADiDynamicOptExt
 using ModelingToolkit
 using CasADi
-using DiffEqDevTools, DiffEqBase
+using DiffEqBase
 using DataInterpolations
 using UnPack
-const MTK = MOdelingToolkit
+const MTK = ModelingToolkit
 
 struct CasADiDynamicOptProblem{uType, tType, isinplace, P, F, K} <:
        AbstractDynamicOptProblem{uType, tType, isinplace}
@@ -221,17 +221,22 @@ function add_solve_constraints!(prob, tableau; is_free_t)
     end
 end
 
-is_explicit(tableau) = tableau isa DiffEqDevTools.ExplicitRKTableau
 
 """
     solve(prob::CasADiDynamicOptProblem, casadi_solver, ode_solver; plugin_options, solver_options) 
 
 `plugin_options` and `solver_options` get propagated to the Opti object in CasADi.
 """
-function DiffEqBase.solve(prob::CasADiDynamicOptProblem, solver::Union{String, Symbol}, ode_solver::Union{String, Symbol}; plugin_options::Dict = Dict(), solver_options::Dict = Dict(), silent = false)
+function DiffEqBase.solve(prob::CasADiDynamicOptProblem, solver::Union{String, Symbol}, ode_solver::Symbol = :Default; plugin_options::Dict = Dict(), solver_options::Dict = Dict(), silent = false)
     model = prob.model
     opti = model.opti
 
+    if ode_solver == :Default
+        tableau = MTK.constructDefault()
+    else
+        tableau_getter = Symbol(:construct, ode_solver)
+        tableau = @eval Main.tableau_getter()
+    end
     solver!(opti, solver, plugin_options, solver_options)
     add_casadi_solve_constraints!(prob, tableau)
     solver!(cmodel, "$solver", plugin_options, solver_options)

src/systems/optimal_control_interface.jl

@@ -18,6 +18,7 @@ end
 
 function JuMPDynamicOptProblem end
 function InfiniteOptDynamicOptProblem end
+function CasADiDynamicOptProblem end
 
 function warn_overdetermined(sys, u0map)
     constraintsys = get_constraintsystem(sys)
@@ -27,6 +28,25 @@ function warn_overdetermined(sys, u0map)
     end
 end
 
+"""
+Default ODE Tableau: RadauIIA5
+"""
+function constructDefault(T::Type = Float64)
+    sq6 = sqrt(6)
+    A = [11 // 45-7sq6 / 360 37 // 225-169sq6 / 1800 -2 // 225+sq6 / 75
+         37 // 225+169sq6 / 1800 11 // 45+7sq6 / 360 -2 // 225-sq6 / 75
+         4 // 9-sq6 / 36 4 // 9+sq6 / 36 1//9]
+    c = [2 // 5 - sq6 / 10; 2 / 5 + sq6 / 10; 1]
+    α = [4 // 9 - sq6 / 36; 4 // 9 + sq6 / 36; 1 // 9]
+    A = map(T, A)
+    α = map(T, α)
+    c = map(T, c)
+    
+    DiffEqBase.ImplicitRKTableau(A, c, α, 5)
+end
+
+is_explicit(tableau) = tableau isa DiffEqBase.ExplicitRKTableau
+
 """
 Generate the control function f(x, u, p, t) from the ODESystem. 
 Input variables are automatically inferred but can be manually specified.