feat: Add support for zero control dimension (NLP_u = 0)

Project.toml

@@ -10,7 +10,6 @@ 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"
@@ -26,7 +25,6 @@ CTParser = "0.8"
 CTSolvers = "0.4"
 CUDA = "5"
 CommonSolve = "0.2"
-DifferentialEquations = "7.17.0"
 DocStringExtensions = "0.9"
 ExaModels = "0.9"
 MadNLP = "0.9"

src/CTDirect.jl

@@ -56,7 +56,7 @@ include("ode/euler.jl")
 include("ode/irk.jl")
 include("ode/midpoint.jl")
 include("ode/trapeze.jl")
-include("ode/variable.jl")
+#include("ode/variable.jl")
 
 include("collocation.jl")
 include("direct_shooting.jl")

src/DOCP_variables.jl

@@ -40,9 +40,11 @@ function __variables_bounds!(docp::DOCP)
     for i in 1:(docp.time.steps + 1)
         set_state_at_time_step!(var_l, x_lb, docp, i)
         set_state_at_time_step!(var_u, x_ub, docp, i)
-        for j in 1:docp.time.control_steps
-            set_control_at_time_step!(var_l, u_lb, docp, i; j=j)
-            set_control_at_time_step!(var_u, u_ub, docp, i; j=j)
+        if docp.dims.NLP_u > 0
+            for j in 1:docp.time.control_steps
+                set_control_at_time_step!(var_l, u_lb, docp, i; j=j)
+                set_control_at_time_step!(var_u, u_ub, docp, i; j=j)
+            end
         end
     end
 
@@ -132,8 +134,10 @@ function __initial_guess(docp::DOCP, init::CTModels.InitialGuess)
     for i in 1:(docp.time.steps + 1)
         ti = time_grid[i]
         set_state_at_time_step!(NLP_X, init.state(ti), docp, i)
-        for j in 1:docp.time.control_steps
-            set_control_at_time_step!(NLP_X, init.control(ti), docp, i; j=j)
+        if docp.dims.NLP_u > 0
+            for j in 1:docp.time.control_steps
+                set_control_at_time_step!(NLP_X, init.control(ti), docp, i; j=j)
+            end
         end
     end
 

src/collocation.jl

@@ -191,10 +191,15 @@ function (discretizer::Collocation)(ocp::AbstractModel)
         # N + 1 states, N controls
         state = hcat([x0[(1 + i * (n + m)):(1 + i * (n + m) + n - 1)] 
         for i in 0:N]...)
-        control = hcat([x0[(n + 1 + i * (n + m)):(n + 1 + i * (n + m) + m - 1)] 
-        for i in 0:(N - 1)]...,)
-        # see with JB: pass indeed to grid_size only for euler(_b), trapeze and midpoint
-        control = [control control[:, end]] 
+        if m > 0
+            control = hcat([x0[(n + 1 + i * (n + m)):(n + 1 + i * (n + m) + m - 1)] 
+            for i in 0:(N - 1)]...,)
+            # see with JB: pass indeed to grid_size only for euler(_b), trapeze and midpoint
+            control = [control control[:, end]]
+        else
+            # zero control dimension: create empty matrix with correct type
+            control = similar(x0, 0, N + 1)
+        end
         variable = x0[(end - q + 1):end]
         init = (variable, state, control)
 

src/ode/common.jl

@@ -80,6 +80,9 @@ function getter(nlp_solution, docp::DOCP; val::Symbol)
 
     # control
     elseif val == :control || val == :control_l || val == :control_u
+        if docp.dims.NLP_u == 0
+            return similar(data, 0, N*docp.time.control_steps + 1)
+        end
         V = zeros(docp.dims.NLP_u, N*docp.time.control_steps + 1)
         k = 1
         for i in 1:N

test/ci/test_zero_control.jl

@@ -0,0 +1,99 @@
+# LATER: put this in test_all_ocp instead ?
+println("testing: zero control dimension - parameter estimation")
+
+# Load test problems
+if !isdefined(Main, :estimate_initial_condition)
+    include("../problems/autonomous_system.jl")
+end
+
+# Test 1: Estimate initial condition with all schemes
+@testset verbose = true showtiming = true ":param_estimation :all_schemes" begin
+    prob = estimate_initial_condition()
+
+    @testset ":euler" begin
+        sol = solve_problem(prob; scheme=:euler, grid_size=50, display=false)
+        @test CTModels.successful(sol)
+        @test sol.objective >= 0
+    end
+
+    @testset ":midpoint" begin
+        sol = solve_problem(prob; scheme=:midpoint, grid_size=50, display=false)
+        @test CTModels.successful(sol)
+        @test sol.objective >= 0
+    end
+
+    @testset ":trapeze" begin
+        sol = solve_problem(prob; scheme=:trapeze, grid_size=50, display=false)
+        @test CTModels.successful(sol)
+        @test sol.objective >= 0
+    end
+end
+
+# Test 2: Estimate parameter in dynamics
+@testset verbose = true showtiming = true ":param_estimation :rotation_rate" begin
+    prob = estimate_rotation_rate()
+    sol = solve_problem(prob; scheme=:midpoint, grid_size=50, display=false)
+    @test CTModels.successful(sol)
+    @test sol.objective >= 0
+    # Verify variable is properly retrieved
+    @test length(variable(sol)) == 1
+end
+
+# Test 3: Least squares with path constraint
+@testset verbose = true showtiming = true ":param_estimation :with_constraint" begin
+    prob = least_squares_with_constraint()
+    sol = solve_problem(prob; scheme=:midpoint, grid_size=50, display=false)
+    @test CTModels.successful(sol)
+    @test sol.objective >= 0
+end
+
+# Test 4: Verify solution dimensions (zero control)
+@testset verbose = true showtiming = true ":param_estimation :solution_dimensions" begin
+    prob = estimate_initial_condition()
+    sol = solve_problem(prob; scheme=:midpoint, grid_size=50, display=false)
+
+    T = time_grid(sol, :state)
+
+    # State must have 2 dimensions
+    @test length(state(sol)(0.5)) == 2
+
+    # Control must be empty
+    @test length(control(sol)(0.5)) == 0
+
+    # Verify control() doesn't crash
+    @test control(sol) isa Function
+    for t in T
+        u = control(sol)(t)
+        @test u isa AbstractVector
+        @test length(u) == 0
+    end
+end
+
+# Test 5: Initial guess with empty control
+@testset verbose = true showtiming = true ":param_estimation :initial_guess" begin
+    prob = estimate_initial_condition()
+
+    # Functional initial guess
+    x_init = t -> [cos(π*t/4), sin(π*t/4)]
+    u_init = t -> Float64[]  # empty vector
+
+    sol = solve_problem(prob; 
+        scheme=:midpoint, 
+        grid_size=50, 
+        init=(state=x_init, control=u_init),
+        display=false
+    )
+    @test CTModels.successful(sol)
+end
+
+# Test 6: ADNLP manual backend (sparsity patterns)
+@testset verbose = true showtiming = true ":param_estimation :adnlp_manual" begin
+    prob = estimate_initial_condition()
+    sol = solve_problem(prob; 
+        scheme=:midpoint, 
+        grid_size=20, 
+        adnlp_backend=:manual,
+        display=false
+    )
+    @test CTModels.successful(sol)
+end

test/ci/test_zero_control_allocations.jl

@@ -0,0 +1,159 @@
+println("testing: zero control allocations")
+
+using CTDirect
+using CTModels
+
+# Load test problems
+if !isdefined(Main, :estimate_initial_condition)
+    include("../problems/autonomous_system.jl")
+end
+
+# Create DOCP with zero control dimension
+@testset verbose = true showtiming = true ":zero_control_allocations" begin
+
+    @testset "DOCP creation" begin
+        prob = estimate_initial_condition()
+        ocp = prob.ocp
+
+        # Build DOCP directly like in get_docp
+        grid_size = 10
+        control_steps = 1
+        scheme = :midpoint
+        time_grid = nothing
+        docp = CTDirect.DOCP(ocp, grid_size, control_steps, scheme, time_grid)
+
+        # Verify NLP_u = 0
+        @test docp.dims.NLP_u == 0
+
+        # Verify dim_NLP_variables is consistent
+        # For midpoint: steps * (NLP_x + NLP_u) + NLP_x + NLP_v
+        # = 10 * (2 + 0) + 2 + 2 = 24 (2 variables for initial condition)
+        @test docp.dim_NLP_variables == 10 * 2 + 2 + 2
+
+        # Verify dim_NLP_constraints is consistent
+        @test docp.dim_NLP_constraints > 0
+    end
+
+    @testset "Control getters return empty views" begin
+        prob = estimate_initial_condition()
+        ocp = prob.ocp
+
+        # Build DOCP directly
+        grid_size = 10
+        control_steps = 1
+        scheme = :midpoint
+        time_grid = nothing
+        docp = CTDirect.DOCP(ocp, grid_size, control_steps, scheme, time_grid)
+
+        # Create dummy xu vector
+        xu = zeros(docp.dim_NLP_variables)
+
+        # get_OCP_control_at_time_step should return empty view
+        u = CTDirect.get_OCP_control_at_time_step(xu, docp, 1)
+        @test u isa AbstractVector
+        @test length(u) == 0
+        @test eltype(u) == Float64
+    end
+
+    @testset "Control setters are no-ops" begin
+        prob = estimate_initial_condition()
+        ocp = prob.ocp
+
+        # Build DOCP directly
+        grid_size = 10
+        control_steps = 1
+        scheme = :midpoint
+        time_grid = nothing
+        docp = CTDirect.DOCP(ocp, grid_size, control_steps, scheme, time_grid)
+
+        xu = zeros(docp.dim_NLP_variables)
+        xu_copy = copy(xu)
+
+        # set_control_at_time_step! should not modify anything
+        CTDirect.set_control_at_time_step!(xu, Float64[], docp, 1)
+        @test xu == xu_copy
+    end
+
+    @testset "All schemes support zero control" begin
+        prob = estimate_initial_condition()
+        ocp = prob.ocp
+
+        for scheme in [:euler, :midpoint, :trapeze]
+            docp = CTDirect.DOCP(ocp, 10, 1, scheme, nothing)
+            @test docp.dims.NLP_u == 0
+            @test docp.dim_NLP_variables > 0
+        end
+    end
+
+    @testset "Variables bounds with zero control" begin
+        prob = estimate_initial_condition()
+        ocp = prob.ocp
+
+        docp = CTDirect.DOCP(ocp, 10, 1, :midpoint, nothing)
+        CTDirect.__variables_bounds!(docp)
+
+        # Verify bounds exist and have correct dimensions
+        @test length(docp.bounds.var_l) == docp.dim_NLP_variables
+        @test length(docp.bounds.var_u) == docp.dim_NLP_variables
+        # Bounds can contain -Inf/Inf for unbounded variables
+        @test all(docp.bounds.var_l .<= docp.bounds.var_u)
+    end
+
+    @testset "Initial guess with zero control" begin
+        prob = estimate_initial_condition()
+        ocp = prob.ocp
+
+        docp = CTDirect.DOCP(ocp, 10, 1, :midpoint, nothing)
+        init = CTModels.build_initial_guess(ocp, ())
+
+        x0 = CTDirect.__initial_guess(docp, init)
+
+        @test length(x0) == docp.dim_NLP_variables
+        @test all(isfinite, x0)
+    end
+
+    @testset "Control getters at all time steps" begin
+        prob = estimate_initial_condition()
+        ocp = prob.ocp
+
+        docp = CTDirect.DOCP(ocp, 10, 1, :midpoint, nothing)
+        xu = zeros(docp.dim_NLP_variables)
+
+        # Test at all time steps
+        for i in 1:(docp.time.steps + 1)
+            u = CTDirect.get_OCP_control_at_time_step(xu, docp, i)
+            @test length(u) == 0
+            @test eltype(u) == Float64
+        end
+    end
+
+    @testset "Zero control with optimization variable" begin
+        prob = estimate_rotation_rate()
+        ocp = prob.ocp
+
+        docp = CTDirect.DOCP(ocp, 10, 1, :midpoint, nothing)
+
+        @test docp.dims.NLP_u == 0
+        @test docp.dims.NLP_v == 1  # one optimization variable
+        @test docp.dim_NLP_variables == 10 * 2 + 2 + 1  # +1 for variable
+    end
+
+    @testset "Sparsity patterns with zero control" begin
+        prob = estimate_initial_condition()
+        ocp = prob.ocp
+
+        docp = CTDirect.DOCP(ocp, 10, 1, :midpoint, nothing)
+
+        # Jacobian pattern should be constructible
+        Is, Js = CTDirect.DOCP_Jacobian_pattern(docp)
+        @test length(Is) == length(Js)
+        @test all(1 .<= Is .<= docp.dim_NLP_constraints)
+        @test all(1 .<= Js .<= docp.dim_NLP_variables)
+
+        # Hessian pattern should be constructible
+        Is_h, Js_h = CTDirect.DOCP_Hessian_pattern(docp)
+        @test length(Is_h) == length(Js_h)
+        @test all(1 .<= Is_h .<= docp.dim_NLP_variables)
+        @test all(1 .<= Js_h .<= docp.dim_NLP_variables)
+    end
+end