PR from fork (#59)

Co-authored-by: corentinmartinon <corentin.martinon@free.fr>

Project.toml

@@ -8,10 +8,10 @@ ADNLPModels = "54578032-b7ea-4c30-94aa-7cbd1cce6c9a"
 CTBase = "54762871-cc72-4466-b8e8-f6c8b58076cd"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
-ADNLPModels = "0.5"
+ADNLPModels = "0.7"
 CTBase = "0.7"
 DocStringExtensions = "0.9"
-NLPModelsIpopt = "0.10"
-julia = "1.8"
+julia = "1.9"

src/CTDirect.jl

@@ -2,6 +2,7 @@ module CTDirect
 
 using CTBase
 using DocStringExtensions
+using Symbolics                   # for optimized auto diff
 using NLPModelsIpopt, ADNLPModels # nlp modeling and resolution
 
 # Other declarations

src/problem.jl

@@ -303,7 +303,7 @@ end
 
 
 # IPOPT constraints  (add bounds computation here at first call ?)
-function ipopt_constraint(xu, ctd)
+function ipopt_constraint!(c, xu, ctd)    
     """
     compute the constraints for the NLP : 
         - discretization of the dynamics via the trapeze method
@@ -327,7 +327,6 @@ function ipopt_constraint(xu, ctd)
     tf = get_final_time(xu, ctd)
     N = ctd.dim_NLP_steps
     h = (tf - t0) / N
-    c = zeros(eltype(xu), ctd.dim_NLP_constraints)
     v = get_variable(xu, ctd)
 
     # state equation
@@ -417,5 +416,5 @@ function ipopt_constraint(xu, ctd)
         index = index + 1
     end
 
-    return c
+    return c # needed even for inplace version, AD error otherwise oO
 end

src/solution.jl

@@ -11,7 +11,8 @@ function _OptimalControlSolution(ocp, ipopt_solution, ctd)
     ctd.NLP_constraints_violation = ipopt_solution.primal_feas
     ctd.NLP_iterations = ipopt_solution.iter
     ctd.NLP_solution = ipopt_solution.solution
-    ctd.NLP_sol_constraints = ipopt_constraint(ipopt_solution.solution, ctd)
+    ctd.NLP_sol_constraints = zeros(ctd.dim_NLP_constraints)
+    ipopt_constraint!(ctd.NLP_sol_constraints, ipopt_solution.solution, ctd)
 
     # parse NLP variables, constraints and multipliers 
     X, U, v, P, sol_control_constraints, sol_state_constraints, sol_mixed_constraints, sol_variable_constraints, mult_control_constraints, mult_state_constraints, mult_mixed_constraints, mult_variable_constraints, mult_state_box_lower, mult_state_box_upper, mult_control_box_lower, mult_control_box_upper, mult_variable_box_lower, mult_variable_box_upper = parse_ipopt_sol(ctd)

src/solve.jl

@@ -46,16 +46,14 @@ function solve(ocp::OptimalControlModel,
         xu0 = initial_guess(ctd)
         l_var, u_var = variables_bounds(ctd)
         lb, ub = constraints_bounds(ctd)
-        nlp = ADNLPModel(xu -> ipopt_objective(xu, ctd), xu0, l_var, u_var, xu -> ipopt_constraint(xu, ctd), lb, ub) 
+        nlp = ADNLPModel!(xu -> ipopt_objective(xu, ctd), xu0, l_var, u_var, (c, xu) -> ipopt_constraint!(c, xu, ctd), lb, ub, backend = :optimized)
     end
 
     # solve
     if :ipopt in method
         # https://github.com/JuliaSmoothOptimizers/NLPModelsIpopt.jl/blob/main/src/NLPModelsIpopt.jl#L119
-        # options of ipopt: https://coin-or.github.io/Ipopt/OPTIONS.html
         # callback: https://github.com/jump-dev/Ipopt.jl#solver-specific-callback
         # sb="yes": remove ipopt header +++ make that default
-        # solve by IPOPT: +++ later use more advanced call for callback use
         print_level = display ?  print_level : 0
         ipopt_solution = ipopt(nlp, print_level=print_level, mu_strategy=mu_strategy, sb="yes"; kwargs...)
     end

test/old/test_integrator.jl

@@ -22,7 +22,7 @@ ocp = prob.model
         @test ctd.has_lagrange_cost         == true
         @test ctd.has_mayer_cost            == false
         @test ctd.dim_NLP_state           == ocp.state_dimension + 1
-        @test ctd.mayer                   == nothing
+        @test ctd.mayer                   === nothing
         f_Mayer_test(t,x,u)=[ocp.dynamics(t,x,u);ocp.lagrange(t,x,u)]
         @test ctd.dynamics_lagrange_to_mayer(0,[0;2],1) == f_Mayer_test(0,[0;2],1)
         @test ctd.has_free_final_time     == false

test/suite/goddard_all_constraints.jl

@@ -50,7 +50,7 @@ function F1(x)
     return [ 0, Tmax/m, -b*Tmax ]
 end
 dynamics!(ocp, (x, u, v) -> F0(x) + u*F1(x) )
-sol1 = solve(ocp, grid_size=30, print_level=0, tol=1e-8)
+sol1 = solve(ocp, grid_size=100, print_level=0, tol=1e-8)
 @testset verbose = true showtiming = true ":goddard :max_rf :all_constraints_types" begin
     @test sol1.objective ≈ 1.0125 rtol=1e-2
 end

test/test_prof.jl

@@ -0,0 +1,65 @@
+#using BenchmarkTools
+#using Traceur
+using Profile
+#using PProf
+#using JET 
+
+using CTDirect
+using CTBase
+
+println("Test: profiling")
+
+# define OCP
+ocp = Model(variable=true)
+Cd = 310
+Tmax = 3.5
+β = 500
+b = 2
+r0 = 1
+v0 = 0
+vmax = 0.1
+m0 = 1
+mf = 0.6
+x0 = [ r0, v0, m0 ]
+state!(ocp, 3)
+control!(ocp, 1)
+variable!(ocp, 1)
+time!(ocp, 0, Index(1))
+constraint!(ocp, :initial, x0, :initial_constraint)
+constraint!(ocp, :final, Index(3), mf, :final_constraint)
+constraint!(ocp, :state, (x,v)->x[2], -Inf, vmax, :state_con_v_ub)
+constraint!(ocp, :control, (u,v)->u, -Inf, 1, :control_con_u_ub)
+constraint!(ocp, :mixed, (x,u,v)->x[3], mf, Inf, :mixed_con_m_lb)
+constraint!(ocp, :variable, v->v, -Inf, 10, :variable_con_tf_ubx)
+constraint!(ocp, :state, 1:2, [r0,v0], [r0+0.2, Inf], :state_box_rv)
+constraint!(ocp, :control, Index(1), 0, Inf, :control_box_lb)
+constraint!(ocp, :variable, Index(1), 0.01, Inf, :variable_box_tfmin)
+objective!(ocp, :mayer,  (x0, xf, v) -> xf[1], :max)
+function F0(x)
+    r, v, m = x
+    D = Cd * v^2 * exp(-β*(r - 1))
+    return [ v, -D/m - 1/r^2, 0 ]
+end
+function F1(x)
+    r, v, m = x
+    return [ 0, Tmax/m, -b*Tmax ]
+end
+dynamics!(ocp, (x, u, v) -> F0(x) + u*F1(x) )
+
+# Solver
+println("First run for compilation")
+@time sol = solve(ocp, grid_size=50, print_level=0, tol=1e-12)
+println("Second run for benchmark")
+@timev sol = solve(ocp, grid_size=50, print_level=5, tol=1e-12)
+
+#= basic benchmark: goddard with 50 steps (second run, 30% compilation time on first one)
+NLP stats: 
+205var, 154const eq, 154 const ineq 
+31570 nnz jac eq/ineq, 21115 nnz hess (non sparse matrices)
+954 nnz jac eq, 154 nnz jac ineq, 561 nnz hess (sparse)
+
+base:                               50iter / 89GB / 64s
+inplace constraints:                50iter / 77GB / 58s
+AD optimized backend*:              50iter / 387MB / 2.5s
+(*slightly worse without inplace)
+=#