Direct shooting

Project.toml

@@ -1,6 +1,6 @@
 name = "CTDirect"
 uuid = "790bbbee-bee9-49ee-8912-a9de031322d5"
-version = "1.0.5"
+version = "1.0.6"
 authors = ["Pierre Martinon <pierrecmartinon@gmail.com>"]
 
 [workspace]
@@ -10,10 +10,12 @@ projects = ["test", "docs"]
 ADNLPModels = "54578032-b7ea-4c30-94aa-7cbd1cce6c9a"
 CTModels = "34c4fa32-2049-4079-8329-de33c2a22e2d"
 CTSolvers = "d3e8d392-8e4b-4d9b-8e92-d7d4e3650ef6"
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 ExaModels = "1037b233-b668-4ce9-9b63-f9f681f55dd2"
 SolverCore = "ff4d7338-4cf1-434d-91df-b86cb86fb843"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
 
 [compat]
 ADNLPModels = "0.8"
@@ -24,6 +26,7 @@ CTParser = "0.8"
 CTSolvers = "0.4"
 CUDA = "5"
 CommonSolve = "0.2"
+DifferentialEquations = "7.17.0"
 DocStringExtensions = "0.9"
 ExaModels = "0.9"
 MadNLP = "0.9"
@@ -32,6 +35,7 @@ NLPModels = "0.21"
 NLPModelsIpopt = "0.11"
 SolverCore = "0.3"
 SparseArrays = "1"
+SparseConnectivityTracer = "1.2.1"
 SplitApplyCombine = "1"
 Test = "1"
 julia = "1.10"

src/CTDirect.jl

@@ -7,165 +7,58 @@ import CTModels
 import CTSolvers, CTSolvers.Strategies, CTSolvers.Options
 import SolverCore
 import SparseArrays
-# using CTBase
-# using NLPModels
-
-# ----------------------------------------------------------------------
-# TYPES
-# ----------------------------------------------------------------------
-const AbstractModel = CTModels.AbstractModel
 
 # ---------------------------------------------------------------------------
 # Abstract discretizer type
 # ---------------------------------------------------------------------------
-"""
-$(TYPEDEF)
-
-Abstract type for optimal control problem discretization strategies.
-
-This type defines the interface for converting continuous optimal control problems
-into discrete optimization problems suitable for numerical solvers. Subtypes must
-implement the callable interface to perform the actual discretization.
-
-# Interface Requirements
-
-Subtypes must implement:
-- `(discretizer::AbstractDiscretizer)(ocp::AbstractModel)`: Perform discretization
-
-# Example
-```julia-repl
-julia> using CTDirect
-
-julia> MyDiscretizer <: AbstractDiscretizer end
-
-julia> # Implement the required interface
-julia> (d::MyDiscretizer)(ocp) = discretize_collocation(ocp)
-```
-
-See also: `Collocation`, `discretize`
-"""
+const AbstractModel = CTModels.AbstractModel
 abstract type AbstractDiscretizer <: Strategies.AbstractStrategy end
 
-"""
-$(TYPEDSIGNATURES)
-
-Discretize an optimal control problem using the specified discretizer.
-
-This function applies a discretization strategy to convert a continuous optimal
-control problem into a form suitable for numerical optimization.
-
-# Arguments
-- `ocp::AbstractModel`: The optimal control problem to discretize
-- `discretizer::AbstractDiscretizer`: The discretization strategy to apply
-
-# Returns
-- The discretized problem (type depends on the specific discretizer)
-
-# Example
-```julia-repl
-julia> using CTDirect
-
-julia> ocp = create_ocp()  # Create your OCP
-julia> discretized = discretize(ocp, Collocation())
-```
+__discretizer()::AbstractDiscretizer = Collocation()
 
-See also: `AbstractDiscretizer`, `Collocation`
-"""
-function discretize(
-    ocp::AbstractModel, 
-    discretizer::AbstractDiscretizer
-)
+function discretize(ocp::AbstractModel, discretizer::AbstractDiscretizer)
     return discretizer(ocp)
 end
-
-"""
-$(TYPEDSIGNATURES)
-
-Returns the default discretizer instance used by the convenience method.
-
-This function provides the default collocation discretizer that is used when
-no specific discretizer is provided to the `discretize` function.
-
-# Returns
-- `AbstractDiscretizer`: The default collocation discretizer instance
-
-# Notes
-- This function is internal and primarily used for providing default values
-- The default is currently `Collocation()` but may change in future versions
-
-See also: `discretize`, `Collocation`
-"""
-__discretizer()::AbstractDiscretizer = Collocation()
-
-"""
-$(TYPEDSIGNATURES)
-
-Discretize an optimal control problem using the default discretizer.
-
-This is a convenience method that uses the default collocation discretizer
-when no specific discretizer is provided.
-
-# Arguments
-- `ocp::AbstractModel`: The optimal control problem to discretize
-- `discretizer::AbstractDiscretizer`: The discretization strategy (default: Collocation())
-
-# Returns
-- The discretized problem using the specified or default discretizer
-
-# Example
-```julia-repl
-julia> using CTDirect
-
-julia> ocp = create_ocp()  # Create your OCP
-julia> discretized = discretize(ocp)  # Uses default Collocation
-julia> discretized_custom = discretize(ocp, MyCustomDiscretizer())
-```
-
-# Notes
-- The default discretizer is `Collocation()`
-- For custom discretization, provide a specific discretizer instance
-
-See also: `AbstractDiscretizer`, `Collocation`
-"""
-function discretize(
-    ocp::AbstractModel;
-    discretizer::AbstractDiscretizer=__discretizer(),
-)
+function discretize(ocp::AbstractModel; discretizer::AbstractDiscretizer=__discretizer())
     return discretize(ocp, discretizer)
 end
 
 # ---------------------------------------------------------------------------
-# Discretization schemes: see disc/
+# Discretization schemes: see ode/
 # ---------------------------------------------------------------------------
 """
 $(TYPEDEF)
 
-Abstract type representing a discretization strategy for an optimal
+Abstract type representing a discretization scheme strategy for an optimal
 control problem.  
 
-Concrete subtypes of `Discretization` define specific schemes for
+Concrete subtypes of `Scheme` define specific schemes for
 transforming a continuous-time problem into a discrete-time
 representation suitable for numerical solution.
 
 # Example
 
 ```julia-repl
-julia> struct MyDiscretization <: Discretization end
-MyDiscretization
+julia> struct MyScheme <: Scheme end
+MyScheme
 ```
 """
-abstract type Discretization end
+abstract type Scheme end
 
 
 # includes
+include("DOCP_data.jl")
+include("DOCP_variables.jl")
+include("DOCP_functions.jl")
+
+include("ode/common.jl")
+include("ode/euler.jl")
+include("ode/irk.jl")
+include("ode/midpoint.jl")
+include("ode/trapeze.jl")
+include("ode/variable.jl")
+
 include("collocation.jl")
-include("collocation_core.jl")
-include("collocation_variables.jl")
-include("collocation_functions.jl")
-include("disc/common.jl")
-include("disc/euler.jl")
-include("disc/irk.jl")
-include("disc/midpoint.jl")
-include("disc/trapeze.jl")
+include("direct_shooting.jl")
 
 end

src/collocation_core.jl → src/DOCP_data.jl

@@ -1,4 +1,4 @@
-# Discretized Optimal Control Problem DOCP
+# Data for iscretized Optimal Control Problem DOCP
 
 """
 $(TYPEDEF)
@@ -146,6 +146,7 @@ DOCPtime(10, [0.0, 0.1, …, 1.0], [0.0, 0.1, …, 1.0])
 """
 struct DOCPtime
     steps::Int
+    control_steps::Int
     normalized_grid::Vector{Float64}
     fixed_grid::Vector{Float64}
 end
@@ -172,7 +173,7 @@ julia> DOCPtime(ocp, 10, nothing)
 DOCPtime(10, [0.0, 0.1, …, 1.0], [0.0, 0.1, …, 1.0])
 ```
 """
-function DOCPtime(ocp::CTModels.Model, grid_size::Int, time_grid)
+function DOCPtime(ocp::CTModels.Model, grid_size::Int, control_steps::Int, time_grid)
 
     # 1. build/recover normalized time grid
     if time_grid === nothing
@@ -209,7 +210,7 @@ function DOCPtime(ocp::CTModels.Model, grid_size::Int, time_grid)
         NLP_fixed_time_grid = @. t0 + (NLP_normalized_time_grid * (tf - t0))
     end
 
-    return DOCPtime(dim_NLP_steps, NLP_normalized_time_grid, NLP_fixed_time_grid)
+    return DOCPtime(dim_NLP_steps, control_steps, NLP_normalized_time_grid, NLP_fixed_time_grid)
 end
 
 """
@@ -262,7 +263,7 @@ DOCP{...}(...)
 ```
 """
 mutable struct DOCP{
-    D<:CTDirect.Discretization, O<:CTModels.Model
+    D<:CTDirect.Scheme, O<:CTModels.Model
     }
 
     # discretization scheme
@@ -289,12 +290,7 @@ mutable struct DOCP{
     dim_NLP_constraints::Int
 
     # constructor
-    function DOCP(
-        ocp::CTModels.Model;
-        grid_size=__grid_size(),
-        time_grid=__time_grid(),
-        scheme=__scheme(),
-        )
+    function DOCP(ocp::CTModels.Model, grid_size::Int, control_steps::Int, scheme::Symbol, time_grid)
 
         # boolean flags
         flags = DOCPFlags(ocp)
@@ -303,10 +299,10 @@ mutable struct DOCP{
         dims = DOCPdims(ocp)
 
         # time grid
-        time = DOCPtime(ocp, grid_size, time_grid)
+        time = DOCPtime(ocp, grid_size, control_steps, time_grid)
 
         # discretization method 
-        disc_args = [time.steps, dims.NLP_x, dims.NLP_u, dims.NLP_v, dims.path_cons, dims.boundary_cons]
+        disc_args = [dims, time]
 
         if scheme == :trapeze
             discretization, dim_NLP_variables, dim_NLP_constraints = 
@@ -331,6 +327,10 @@ mutable struct DOCP{
             discretization, dim_NLP_variables, dim_NLP_constraints = 
             CTDirect.Gauss_Legendre_3(disc_args...)
 
+        elseif scheme == :variable
+            discretization, dim_NLP_variables, dim_NLP_constraints = 
+            CTDirect.VariableStepODE(disc_args...)
+
         else
             error(
                 "Unknown discretization method: ",
@@ -519,7 +519,7 @@ function build_OCP_solution(docp::DOCP, nlp_solution::SolverCore.AbstractExecuti
     # state and control variables (allocs seem needed here)
     N = docp.time.steps
     X = zeros(N + 1, docp.dims.NLP_x)
-    U = zeros(N + 1, docp.dims.NLP_u)
+    U = zeros(N*docp.time.control_steps + 1, docp.dims.NLP_u)
     v = zeros(docp.dims.NLP_v)
 
     # multipliers for box constraints
@@ -539,10 +539,24 @@ function build_OCP_solution(docp::DOCP, nlp_solution::SolverCore.AbstractExecuti
     end
 
     # retrieve state, control and optimization variables
+    T_state = T
     X[:] = getter(nlp_solution; val=:state)'
     U[:] = getter(nlp_solution; val=:control)'
     v[:] = getter(nlp_solution; val=:variable)
 
+    # NB. later, use function in control parametrization
+    # collocation case: T_control = T
+    T_control = zeros(N*docp.time.control_steps + 1)
+    k = 1
+    for i=1:N
+        h_i = T[i+1] - T[i]
+        for j=1:docp.time.control_steps
+            T_control[k] = T[i] + (j-1) * h_i / docp.time.control_steps
+            k += 1
+        end  
+    end
+    T_control[end] = T[end]
+
     # lower bounds multiplier
     if !is_empty(multipliers_L)
         mult_state_box_lower[:] = getter(nlp_solution; val=:state_l)'
@@ -557,15 +571,17 @@ function build_OCP_solution(docp::DOCP, nlp_solution::SolverCore.AbstractExecuti
         mult_variable_box_upper[:] = getter(nlp_solution; val=:variable_u)
     end
 
-    # costate and constraints multipliers
+    # costate
+    P = zeros(N, docp.dims.NLP_x)
+    P[:] = getter(nlp_solution; val=:costate)'
+    T_costate = T[1:end-1]
+
+    # constraints multipliers
     dpc = docp.dims.path_cons
     dbc = docp.dims.boundary_cons
-    P = zeros(N, docp.dims.NLP_x)
     mult_path_constraints = zeros(N + 1, dpc)
     mult_boundary_constraints = zeros(dbc)
-
-    # costate
-    P[:] = getter(nlp_solution; val=:costate)'
+    T_path = T
 
     # path constraints multipliers (+++ not yet implemented for exa)
     # NB. normalize wrt time step
@@ -580,9 +596,13 @@ function build_OCP_solution(docp::DOCP, nlp_solution::SolverCore.AbstractExecuti
     # boundary constraints multipliers (+++ not yet implemented for exa)
     isnothing(exa_getter) && (mult_boundary_constraints = getter(nlp_solution; val=:mult_boundary_constraints))
 
+    # use method with separate time grids
     return CTModels.build_solution(
         ocp_model(docp),
-        T,
+        T_state,
+        T_control,
+        T_costate,
+        T_path,
         X,
         U,
         v,