diff --git a/Project.toml b/Project.toml index 945810f0..cdfe16e6 100644 --- a/Project.toml +++ b/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" diff --git a/src/CTDirect.jl b/src/CTDirect.jl index f88d7f07..23e3618e 100644 --- a/src/CTDirect.jl +++ b/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") diff --git a/src/DOCP_variables.jl b/src/DOCP_variables.jl index 267cc804..9fc0241b 100644 --- a/src/DOCP_variables.jl +++ b/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 diff --git a/src/collocation.jl b/src/collocation.jl index e6d8ff2b..ca1812c4 100644 --- a/src/collocation.jl +++ b/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) diff --git a/src/ode/common.jl b/src/ode/common.jl index 4db4deef..761abfbb 100644 --- a/src/ode/common.jl +++ b/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 diff --git a/test/ci/test_zero_control.jl b/test/ci/test_zero_control.jl new file mode 100644 index 00000000..52b2162b --- /dev/null +++ b/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 diff --git a/test/ci/test_zero_control_allocations.jl b/test/ci/test_zero_control_allocations.jl new file mode 100644 index 00000000..902b45b5 --- /dev/null +++ b/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 diff --git a/test/problems/autonomous_system.jl b/test/problems/autonomous_system.jl new file mode 100644 index 00000000..65f95247 --- /dev/null +++ b/test/problems/autonomous_system.jl @@ -0,0 +1,138 @@ +# Parameter estimation problems without control (dimension 0) +# Using CTModels directly since CTParser does not yet support zero control dimension +# CTModels must be loaded before including this file + +# Problem 1: Estimate initial condition of harmonic oscillator +function estimate_initial_condition() + pre = CTModels.PreModel() + CTModels.time!(pre; t0=0.0, tf=π/2) # quarter period + CTModels.state!(pre, 2) + CTModels.variable!(pre, 2) # v = [x₁(0), x₂(0)] to estimate + + # Dynamics: harmonic oscillator ẋ = [-x₂, x₁] + dynamics!(r, t, x, u, v) = r .= [-x[2], x[1]] + CTModels.dynamics!(pre, dynamics!) + + # Objective: minimize distance to target final state [0, 1] + # Solution should be v ≈ [1, 0] + CTModels.objective!(pre, :min; + mayer=(x0, xf, v) -> (xf[1] - 0.0)^2 + (xf[2] - 1.0)^2 + ) + + CTModels.definition!(pre, quote + v ∈ R², variable + t ∈ [0, π/2], time + x ∈ R², state + x(0) == v + ẋ(t) == [-x₂(t), x₁(t)] + (xf[1])^2 + (xf[2] - 1)^2 → min + end) + CTModels.time_dependence!(pre; autonomous=true) + + # Initial condition is the variable to estimate + f_initial(r, x0, xf, v) = r .= x0 .- v + CTModels.constraint!(pre, :boundary; + f=f_initial, + lb=[0.0, 0.0], + ub=[0.0, 0.0], + label=:initial + ) + + ocp = CTModels.build(pre) + return ((ocp=ocp, obj=nothing, name="estimate_initial", init=())) +end + +# Problem 2: Estimate parameter in dynamics (rotation rate) +function estimate_rotation_rate() + pre = CTModels.PreModel() + CTModels.time!(pre; t0=0.0, tf=1.0) + CTModels.state!(pre, 2) + CTModels.variable!(pre, 1) # v = [α] rotation rate to estimate + + # Dynamics: ẋ = α*[-x₂, x₁] (rotation with unknown rate) + dynamics!(r, t, x, u, v) = r .= v[1] .* [-x[2], x[1]] + CTModels.dynamics!(pre, dynamics!) + + # Objective: minimize distance to target final state [0, 1] + # and regularize parameter (solution should be α ≈ π/2) + CTModels.objective!(pre, :min; + mayer=(x0, xf, v) -> (xf[1] - 0.0)^2 + (xf[2] - 1.0)^2 + 0.01*v[1]^2 + ) + + CTModels.definition!(pre, quote + α ∈ R, variable + t ∈ [0, 1], time + x ∈ R², state + 0 ≤ α ≤ 10 + x(0) == [1, 0] + ẋ(t) == α * [-x₂(t), x₁(t)] + (xf[1])^2 + (xf[2] - 1)^2 + 0.01*α^2 → min + end) + CTModels.time_dependence!(pre; autonomous=false) + + # Known initial condition + f_initial(r, x0, xf, v) = r .= x0 .- [1.0, 0.0] + CTModels.constraint!(pre, :boundary; + f=f_initial, + lb=[0.0, 0.0], + ub=[0.0, 0.0], + label=:initial + ) + + # Bounds on parameter + CTModels.constraint!(pre, :variable; rg=1, lb=0.0, ub=10.0, label=:alpha_bounds) + + ocp = CTModels.build(pre) + return ((ocp=ocp, obj=nothing, name="estimate_parameter", init=())) +end + +# Problem 3: Least squares fit with path constraint +function least_squares_with_constraint() + pre = CTModels.PreModel() + CTModels.time!(pre; t0=0.0, tf=1.0) + CTModels.state!(pre, 2) + CTModels.variable!(pre, 2) # v = [x₁(0), x₂(0)] to estimate + + # Dynamics: harmonic oscillator + dynamics!(r, t, x, u, v) = r .= [-x[2], x[1]] + CTModels.dynamics!(pre, dynamics!) + + # Objective: minimize distance to measurements along trajectory + # Synthetic measurements at t=0.5: [0.7, 0.7] + CTModels.objective!(pre, :min; + lagrange=(t, x, u, v) -> (t - 0.5)^2 * ((x[1] - 0.7)^2 + (x[2] - 0.7)^2), + mayer=(x0, xf, v) -> 0.01 * (v[1]^2 + v[2]^2) + ) + + CTModels.definition!(pre, quote + v ∈ R², variable + t ∈ [0, 1], time + x ∈ R², state + x(0) == v + ẋ(t) == [-x₂(t), x₁(t)] + x₁(t)^2 + x₂(t)^2 ≤ 2 + ∫((t - 0.5)^2 * ((x₁(t) - 0.7)^2 + (x₂(t) - 0.7)^2)) + 0.01*(v[1]^2 + v[2]^2) → min + end) + CTModels.time_dependence!(pre; autonomous=true) + + # Initial condition is the variable + f_initial(r, x0, xf, v) = r .= x0 .- v + CTModels.constraint!(pre, :boundary; + f=f_initial, + lb=[0.0, 0.0], + ub=[0.0, 0.0], + label=:initial + ) + + # Path constraint: stay within radius + f_path(r, t, x, u, v) = r .= x[1]^2 + x[2]^2 + CTModels.constraint!(pre, :path; + f=f_path, + lb=-Inf, + ub=2.0, + label=:radius + ) + + ocp = CTModels.build(pre) + return ((ocp=ocp, obj=nothing, name="least_squares_constraint", init=())) +end diff --git a/test/runtests.jl b/test/runtests.jl index ec17c6d2..b14a541b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -27,4 +27,4 @@ const SHOWTIMING = true # CPU: run all scripts in subfolder ci/ include.(filter(contains(r".jl$"), readdir("./ci"; join=true))) end -end +end \ No newline at end of file