diff --git a/docs/make.jl b/docs/make.jl
index 38f6b41c..5f70dd4b 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -7,26 +7,26 @@ using CTDirect
 
 # to add docstrings from external packages
 Modules = [CTBase, CTFlows, CTDirect, OptimalControl]
-for Module ∈ Modules
+for Module in Modules
     isnothing(DocMeta.getdocmeta(Module, :DocTestSetup)) &&
-        DocMeta.setdocmeta!(Module, :DocTestSetup, :(using $Module); recursive = true)
+        DocMeta.setdocmeta!(Module, :DocTestSetup, :(using $Module); recursive=true)
 end
 
 repo_url = "github.com/control-toolbox/OptimalControl.jl"
 
-makedocs(
-    warnonly = [:cross_references, :autodocs_block],
-    sitename = "OptimalControl.jl",
-    format = Documenter.HTML(;
-        repolink = "https://" * repo_url,
-        prettyurls = false,
-        size_threshold_ignore = ["api-ctbase/types.md", "tutorial-plot.md"],
-        assets = [
+makedocs(;
+    warnonly=[:cross_references, :autodocs_block],
+    sitename="OptimalControl.jl",
+    format=Documenter.HTML(;
+        repolink="https://" * repo_url,
+        prettyurls=false,
+        size_threshold_ignore=["api-ctbase/types.md", "tutorial-plot.md"],
+        assets=[
             asset("https://control-toolbox.org/assets/css/documentation.css"),
             asset("https://control-toolbox.org/assets/js/documentation.js"),
         ],
     ),
-    pages = [
+    pages=[
         "Introduction" => "index.md",
         "Basic examples" => [
             "Energy minimisation" => "tutorial-double-integrator-energy.md",
@@ -64,4 +64,4 @@ makedocs(
     ],
 )
 
-deploydocs(; repo = repo_url * ".git", devbranch = "main")
+deploydocs(; repo=repo_url * ".git", devbranch="main")
diff --git a/src/OptimalControl.jl b/src/OptimalControl.jl
index 89ee0fdb..f9ece477 100644
--- a/src/OptimalControl.jl
+++ b/src/OptimalControl.jl
@@ -53,7 +53,8 @@ export OptimalControlModel, OptimalControlSolution
 export Autonomous, NonAutonomous
 export NonFixed, Fixed
 export Model, __OCPModel
-export variable!, time!, constraint!, dynamics!, objective!, state!, control!, remove_constraint!
+export variable!,
+    time!, constraint!, dynamics!, objective!, state!, control!, remove_constraint!
 export constraint
 export time_grid, control, state, variable, costate, objective
 export iterations, stopping, message, infos
@@ -65,7 +66,7 @@ export ParsingError
 # repl
 function __init__()
     isdefined(Base, :active_repl) && ct_repl()
-    nothing
+    return nothing
 end
 
 end
diff --git a/src/solve.jl b/src/solve.jl
index 4edfb471..d9317613 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -6,7 +6,7 @@ Return the list of available methods that can be used to solve the optimal contr
 function available_methods()
     # by order of preference: from top to bottom
     methods = ()
-    for method ∈ CTDirect.available_methods()
+    for method in CTDirect.available_methods()
         methods = add(methods, (:direct, method...))
     end
     return methods
@@ -84,14 +84,15 @@ julia> sol = solve(ocp, init=(state=[-0.5, 0.2], control=0.5))
     For more information on how to provide the initial guess, see the [tutorial on the initial guess](@ref tutorial-init).
 
 """
-function CommonSolve.solve(ocp::OptimalControlModel, description::Symbol...; kwargs...)::OptimalControlSolution
+function CommonSolve.solve(
+    ocp::OptimalControlModel, description::Symbol...; kwargs...
+)::OptimalControlSolution
 
     # get the full description
     method = getFullDescription(description, available_methods())
 
     # solve the problem
-    if :direct ∈ method 
+    if :direct ∈ method
         return CTDirect.direct_solve(ocp, clean(description)...; kwargs...)
     end
-
 end
diff --git a/test/problems/beam.jl b/test/problems/beam.jl
index ed095092..b4bd7e01 100644
--- a/test/problems/beam.jl
+++ b/test/problems/beam.jl
@@ -13,5 +13,5 @@ function beam()
         ∫(u(t)^2) → min
     end
 
-    return ((ocp = beam, obj = 8.898598, name = "beam", init = nothing))
+    return ((ocp=beam, obj=8.898598, name="beam", init=nothing))
 end
diff --git a/test/problems/bioreactor.jl b/test/problems/bioreactor.jl
index 29def3a4..34cebd09 100644
--- a/test/problems/bioreactor.jl
+++ b/test/problems/bioreactor.jl
@@ -55,7 +55,7 @@ function bioreactor_1day_periodic()
         ∫(mu2 * b(t) / (beta + c)) → max
     end
 
-    return ((ocp = bioreactor_1, obj = 0.614134, name = "bioreactor_1day_periodic", init = nothing))
+    return ((ocp=bioreactor_1, obj=0.614134, name="bioreactor_1day_periodic", init=nothing))
 end
 
 # N days (non periodic)
@@ -96,9 +96,6 @@ function bioreactor_Ndays()
     end
 
     return ((
-        ocp = bioreactor_N,
-        obj = 19.0745,
-        init = (state = [50, 50, 50],),
-        name = "bioreactor_Ndays",
+        ocp=bioreactor_N, obj=19.0745, init=(state=[50, 50, 50],), name="bioreactor_Ndays"
     ))
 end
diff --git a/test/problems/fuller.jl b/test/problems/fuller.jl
index d006c364..44c58427 100644
--- a/test/problems/fuller.jl
+++ b/test/problems/fuller.jl
@@ -11,5 +11,5 @@ function fuller()
         ∫(x₁(t)^2) → min
     end
 
-    return ((ocp = fuller, obj = 2.683944e-1, name = "fuller", init = nothing))
+    return ((ocp=fuller, obj=2.683944e-1, name="fuller", init=nothing))
 end
diff --git a/test/problems/goddard.jl b/test/problems/goddard.jl
index 5c1a41ce..ecfa6373 100644
--- a/test/problems/goddard.jl
+++ b/test/problems/goddard.jl
@@ -13,7 +13,7 @@ function F1(x, Tmax, b)
     return [0, Tmax / m, -b * Tmax]
 end
 
-function goddard(; vmax = 0.1, Tmax = 3.5, functional_constraints = false)
+function goddard(; vmax=0.1, Tmax=3.5, functional_constraints=false)
     # constants
     Cd = 310
     beta = 500
@@ -25,31 +25,33 @@ function goddard(; vmax = 0.1, Tmax = 3.5, functional_constraints = false)
     x0 = [r0, v0, m0]
 
     #ocp
-    goddard = Model(variable = true)
+    goddard = Model(; variable=true)
     state!(goddard, 3)
     control!(goddard, 1)
     variable!(goddard, 1)
-    time!(goddard, t0 = 0, indf = 1)
-    constraint!(goddard, :initial, lb = x0, ub = x0)
-    constraint!(goddard, :final, rg = 3, lb = mf, ub = Inf)
+    time!(goddard; t0=0, indf=1)
+    constraint!(goddard, :initial; lb=x0, ub=x0)
+    constraint!(goddard, :final; rg=3, lb=mf, ub=Inf)
     if functional_constraints
         # note: the equations do not handle r<1 well
         # without the box constraint on x, the default init (0.1) is not suitable
-        constraint!(goddard, :state, f = (x, v) -> x, lb = [r0, v0, mf], ub = [r0 + 0.2, vmax, m0])
-        constraint!(goddard, :control, f = (u, v) -> u, lb = 0, ub = 1)
+        constraint!(
+            goddard, :state; f=(x, v) -> x, lb=[r0, v0, mf], ub=[r0 + 0.2, vmax, m0]
+        )
+        constraint!(goddard, :control; f=(u, v) -> u, lb=0, ub=1)
     else
-        constraint!(goddard, :state, lb = [r0, v0, mf], ub = [r0 + 0.2, vmax, m0])
-        constraint!(goddard, :control, lb = 0, ub = 1)
+        constraint!(goddard, :state; lb=[r0, v0, mf], ub=[r0 + 0.2, vmax, m0])
+        constraint!(goddard, :control; lb=0, ub=1)
     end
-    constraint!(goddard, :variable, lb = 0.01, ub = Inf)
+    constraint!(goddard, :variable; lb=0.01, ub=Inf)
     objective!(goddard, :mayer, (x0, xf, v) -> xf[1], :max)
     dynamics!(goddard, (x, u, v) -> F0(x, Cd, beta) + u * F1(x, Tmax, b))
 
-    return ((ocp = goddard, obj = 1.01257, name = "goddard", init = (state = [1.01, 0.05, 0.8],)))
+    return ((ocp=goddard, obj=1.01257, name="goddard", init=(state=[1.01, 0.05, 0.8],)))
 end
 
 # abstratc definition
-function goddard_a(; vmax = 0.1, Tmax = 3.5)
+function goddard_a(; vmax=0.1, Tmax=3.5)
     # constants
     Cd = 310
     beta = 500
@@ -79,10 +81,5 @@ function goddard_a(; vmax = 0.1, Tmax = 3.5)
         r(tf) → max
     end
 
-    return ((
-        ocp = goddard_a,
-        obj = 1.01257,
-        name = "goddard_a",
-        init = (state = [1.01, 0.05, 0.8],),
-    ))
+    return ((ocp=goddard_a, obj=1.01257, name="goddard_a", init=(state=[1.01, 0.05, 0.8],)))
 end
diff --git a/test/problems/insurance.jl b/test/problems/insurance.jl
index 03398dc5..d9119d81 100644
--- a/test/problems/insurance.jl
+++ b/test/problems/insurance.jl
@@ -53,5 +53,5 @@ function insurance()
         ∫(U(t) * fx) → max
     end
 
-    return ((ocp = insurance, obj = 2.059511, name = "insurance", init = nothing))
+    return ((ocp=insurance, obj=2.059511, name="insurance", init=nothing))
 end
diff --git a/test/problems/jackson.jl b/test/problems/jackson.jl
index e36e81eb..85605373 100644
--- a/test/problems/jackson.jl
+++ b/test/problems/jackson.jl
@@ -22,5 +22,5 @@ function jackson()
         x[3](4) → max
     end
 
-    return ((ocp = jackson, obj = 0.192011, name = "jackson", init = nothing))
+    return ((ocp=jackson, obj=0.192011, name="jackson", init=nothing))
 end
diff --git a/test/problems/parametric.jl b/test/problems/parametric.jl
index 311dc0d5..71acb6dc 100644
--- a/test/problems/parametric.jl
+++ b/test/problems/parametric.jl
@@ -25,4 +25,4 @@ function myocp(ρ)
 end
 
 # return problem and objective
-return ((ocp = myocp(1.0), obj = -0.335936, name = "param"))
+return ((ocp=myocp(1.0), obj=-0.335936, name="param"))
diff --git a/test/problems/robbins.jl b/test/problems/robbins.jl
index 6c53191d..9bea842e 100644
--- a/test/problems/robbins.jl
+++ b/test/problems/robbins.jl
@@ -16,5 +16,5 @@ function robbins()
         ∫(alpha * x[1](t) + beta * x[1](t)^2 + gamma * u(t)^2) → min
     end
 
-    return ((ocp = robbins, obj = 20.0, name = "robbins", init = nothing))
+    return ((ocp=robbins, obj=20.0, name="robbins", init=nothing))
 end
diff --git a/test/problems/swimmer.jl b/test/problems/swimmer.jl
index e2b84a53..3283abd6 100644
--- a/test/problems/swimmer.jl
+++ b/test/problems/swimmer.jl
@@ -43,18 +43,20 @@ function swimmer()
             2 * cos(2 * (b1 - b3)) +
             12 * cos(2 * b1 - b3) +
             186 * cos(b3) +
-            37 * cos(2 * b3) - 6 * cos(b1 + b3) - 3 * cos(2 * (b1 + b3)) - 6 * cos(2 * b1 + b3) -
-            6 * cos(b1 + 2 * b3)
+            37 * cos(2 * b3) - 6 * cos(b1 + b3) - 3 * cos(2 * (b1 + b3)) -
+            6 * cos(2 * b1 + b3) - 6 * cos(b1 + 2 * b3)
 
         g11 =
             (
-                -42 * sin(b1 - th) - 2 * sin(2 * b1 - th) - 24 * sin(th) - 300 * sin(b1 + th) -
-                12 * sin(2 * b1 + th) - 6 * sin(b1 - th - 2 * b3) - sin(2 * b1 - th - 2 * b3) +
-                4 * sin(th - 2 * b3) - 12 * sin(b1 + th - 2 * b3) - sin(2 * b1 + th - 2 * b3) +
+                -42 * sin(b1 - th) - 2 * sin(2 * b1 - th) - 24 * sin(th) -
+                300 * sin(b1 + th) - 12 * sin(2 * b1 + th) - 6 * sin(b1 - th - 2 * b3) -
+                sin(2 * b1 - th - 2 * b3) + 4 * sin(th - 2 * b3) -
+                12 * sin(b1 + th - 2 * b3) - sin(2 * b1 + th - 2 * b3) +
                 18 * sin(b1 - th - b3) +
                 8 * sin(th - b3) - 54 * sin(b1 + th - b3) - 2 * sin(2 * b1 + th - b3) -
                 18 * sin(b1 - th + b3) - 38 * sin(th + b3) - 90 * sin(b1 + th + b3) -
-                6 * sin(b1 - th + 2 * b3) - 18 * sin(th + 2 * b3) - 30 * sin(b1 + th + 2 * b3)
+                6 * sin(b1 - th + 2 * b3) - 18 * sin(th + 2 * b3) -
+                30 * sin(b1 + th + 2 * b3)
             ) / (4 * aux)
 
         g12 =
@@ -62,8 +64,8 @@ function swimmer()
                 -42 * cos(b1 - th) - 2 * cos(2 * b1 - th) +
                 24 * cos(th) +
                 300 * cos(b1 + th) +
-                12 * cos(2 * b1 + th) - 6 * cos(b1 - th - 2 * b3) - cos(2 * b1 - th - 2 * b3) -
-                4 * cos(th - 2 * b3) +
+                12 * cos(2 * b1 + th) - 6 * cos(b1 - th - 2 * b3) -
+                cos(2 * b1 - th - 2 * b3) - 4 * cos(th - 2 * b3) +
                 12 * cos(b1 + th - 2 * b3) +
                 cos(2 * b1 + th - 2 * b3) +
                 18 * cos(b1 - th - b3) - 8 * cos(th - b3) +
@@ -82,7 +84,8 @@ function swimmer()
                 2 * cos(2 * b1) +
                 12 * cos(b1 - 2 * b3) +
                 30 * cos(b1 - b3) +
-                cos(2 * (b1 - b3)) - 4 * cos(2 * b3) - 6 * cos(b1 + b3) - 6 * cos(b1 + 2 * b3)
+                cos(2 * (b1 - b3)) - 4 * cos(2 * b3) - 6 * cos(b1 + b3) -
+                6 * cos(b1 + 2 * b3)
             ) / (2 * aux)
 
         g21 =
@@ -91,8 +94,9 @@ function swimmer()
                 4 * sin(2 * b1 - th) +
                 24 * sin(th) +
                 38 * sin(b1 + th) +
-                18 * sin(2 * b1 + th) - 2 * sin(b1 - th - 2 * b3) - sin(2 * b1 - th - 2 * b3) -
-                2 * sin(th - 2 * b3) - sin(2 * b1 + th - 2 * b3) - 54 * sin(b1 - th - b3) -
+                18 * sin(2 * b1 + th) - 2 * sin(b1 - th - 2 * b3) -
+                sin(2 * b1 - th - 2 * b3) - 2 * sin(th - 2 * b3) -
+                sin(2 * b1 + th - 2 * b3) - 54 * sin(b1 - th - b3) -
                 12 * sin(2 * b1 - th - b3) - 42 * sin(th - b3) + 18 * sin(b1 + th - b3) -
                 6 * sin(2 * b1 + th - b3) +
                 18 * sin(b1 - th + b3) +
@@ -106,10 +110,11 @@ function swimmer()
         g22 =
             (
                 8 * cos(b1 - th) + 4 * cos(2 * b1 - th) - 24 * cos(th) - 38 * cos(b1 + th) -
-                18 * cos(2 * b1 + th) - 2 * cos(b1 - th - 2 * b3) - cos(2 * b1 - th - 2 * b3) +
+                18 * cos(2 * b1 + th) - 2 * cos(b1 - th - 2 * b3) -
+                cos(2 * b1 - th - 2 * b3) +
                 2 * cos(th - 2 * b3) +
-                cos(2 * b1 + th - 2 * b3) - 54 * cos(b1 - th - b3) - 12 * cos(2 * b1 - th - b3) +
-                42 * cos(th - b3) - 18 * cos(b1 + th - b3) +
+                cos(2 * b1 + th - 2 * b3) - 54 * cos(b1 - th - b3) -
+                12 * cos(2 * b1 - th - b3) + 42 * cos(th - b3) - 18 * cos(b1 + th - b3) +
                 6 * cos(2 * b1 + th - b3) +
                 18 * cos(b1 - th + b3) +
                 6 * cos(2 * b1 - th + b3) - 300 * cos(th + b3) - 90 * cos(b1 + th + b3) -
@@ -131,5 +136,5 @@ function swimmer()
         #∫(u1^2 + u2^2) → min
     end
 
-    return ((ocp = swimmer, obj = 0.984273, name = "swimmer", init = nothing))
+    return ((ocp=swimmer, obj=0.984273, name="swimmer", init=nothing))
 end
diff --git a/test/problems/vanderpol.jl b/test/problems/vanderpol.jl
index 4f23beb3..6c99914b 100644
--- a/test/problems/vanderpol.jl
+++ b/test/problems/vanderpol.jl
@@ -9,9 +9,10 @@ function vanderpol()
         x ∈ R², state
         u ∈ R, control
         x(0) == [1, 0]
-        ẋ(t) == [x[2](t), epsilon * omega * (1 - x[1](t)^2) * x[2](t) - omega^2 * x[1](t) + u(t)]
+        ẋ(t) ==
+        [x[2](t), epsilon * omega * (1 - x[1](t)^2) * x[2](t) - omega^2 * x[1](t) + u(t)]
         ∫(0.5 * (x[1](t)^2 + x[2](t)^2 + u(t)^2)) → min
     end
 
-    return ((ocp = vanderpol, obj = 1.047921, name = "vanderpol", init = nothing))
+    return ((ocp=vanderpol, obj=1.047921, name="vanderpol", init=nothing))
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 0cd75a8c..fa13f25a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -12,7 +12,7 @@ include("problems/goddard.jl")
 
 #
 @testset verbose = true showtiming = true "Optimal control tests" begin
-    for name ∈ (
+    for name in (
         :abstract_ocp,
         :basic,
         :continuation,
diff --git a/test/test-nlp.jl b/test/test-nlp.jl
index 700df1e3..1e0199db 100644
--- a/test/test-nlp.jl
+++ b/test/test-nlp.jl
@@ -19,8 +19,8 @@ docp, nlp = direct_transcription(ocp)
 #println(nlp.meta.x0)
 
 using NLPModelsIpopt
-nlp_sol = ipopt(nlp; print_level = 5, mu_strategy = "adaptive", tol = 1e-8, sb = "yes")
-sol = OptimalControlSolution(docp, primal = nlp_sol.solution, dual = nlp_sol.multipliers)
+nlp_sol = ipopt(nlp; print_level=5, mu_strategy="adaptive", tol=1e-8, sb="yes")
+sol = OptimalControlSolution(docp; primal=nlp_sol.solution, dual=nlp_sol.multipliers)
 plot(sol)
 
 using MadNLP
@@ -30,10 +30,12 @@ set_initial_guess(docp, nlp, sol)
 #println(nlp.meta.x0)
 mad_nlp_sol = madnlp(nlp)
 
-docp, nlp = direct_transcription(ocp, init = sol)
+docp, nlp = direct_transcription(ocp; init=sol)
 mad_nlp_sol = madnlp(nlp)
 
-mad_sol = OptimalControlSolution(docp, primal = madnlp_sol.solution, dual = madnlp_sol.multipliers)
+mad_sol = OptimalControlSolution(
+    docp; primal=madnlp_sol.solution, dual=madnlp_sol.multipliers
+)
 plot(mad_sol)
 
 #=
diff --git a/test/test_abstract_ocp.jl b/test/test_abstract_ocp.jl
index 17d2aced..ecdfdfe2 100644
--- a/test/test_abstract_ocp.jl
+++ b/test/test_abstract_ocp.jl
@@ -16,13 +16,13 @@ function test_abstract_ocp()
     end
 
     @testset verbose = true showtiming = true ":double_integrator :min_tf :abstract" begin
-        sol1 = solve(ocp1, print_level = 0, tol = 1e-12)
+        sol1 = solve(ocp1; print_level=0, tol=1e-12)
         @test sol1.objective ≈ 2.0 rtol = 1e-2
     end
 
     @testset verbose = true showtiming = true ":goddard :max_rf :abstract :constr" begin
         ocp3 = goddard_a().ocp
-        sol3 = solve(ocp3, print_level = 0, tol = 1e-12)
+        sol3 = solve(ocp3; print_level=0, tol=1e-12)
         @test sol3.objective ≈ 1.0125 rtol = 1e-2
     end
 end
diff --git a/test/test_basic.jl b/test/test_basic.jl
index efff715f..cc44f999 100644
--- a/test/test_basic.jl
+++ b/test/test_basic.jl
@@ -9,6 +9,6 @@ function test_basic()
         ∫(0.5u(t)^2) → min
     end
 
-    sol = solve(ocp; display = false)
+    sol = solve(ocp; display=false)
     @test objective(sol) ≈ 6 atol = 1e-2
 end
diff --git a/test/test_continuation.jl b/test/test_continuation.jl
index ac7a6288..e4894e72 100644
--- a/test/test_continuation.jl
+++ b/test/test_continuation.jl
@@ -26,9 +26,9 @@ function test_continuation()
         @testset verbose = true showtiming = true ":parametric_ocp :warm_start" begin
             init = ()
             obj_list = []
-            for T = 1:5
+            for T in 1:5
                 ocp = ocp_T(T)
-                sol = solve(ocp, print_level = 0, init = init)
+                sol = solve(ocp; print_level=0, init=init)
                 init = sol
                 push!(obj_list, objective(sol))
             end
@@ -67,7 +67,7 @@ function test_continuation()
             obj_list = []
             for ρ in [0.1, 5, 10, 30, 100]
                 ocp = myocp(ρ)
-                sol = solve(ocp, print_level = 0, init = init)
+                sol = solve(ocp; print_level=0, init=init)
                 init = sol
                 push!(obj_list, objective(sol))
             end
@@ -77,14 +77,14 @@ function test_continuation()
 
     # parametric ocp definition
     if test3
-        sol0 = solve(goddard().ocp, print_level = 0)
+        sol0 = solve(goddard().ocp; print_level=0)
 
         @testset verbose = true showtiming = true ":global_variable :warm_start" begin
             sol = sol0
             Tmax_list = []
             obj_list = []
-            for Tmax = 3.5:-0.5:1
-                sol = solve(goddard(Tmax = Tmax).ocp, print_level = 0, init = sol)
+            for Tmax in 3.5:-0.5:1
+                sol = solve(goddard(; Tmax=Tmax).ocp; print_level=0, init=sol)
                 push!(Tmax_list, Tmax)
                 push!(obj_list, objective(sol))
             end
@@ -94,16 +94,16 @@ function test_continuation()
                 # plot obj(vmax)
                 pobj = plot(
                     Tmax_list,
-                    obj_list,
-                    label = "r(tf)",
-                    xlabel = "Maximal thrust (Tmax)",
-                    ylabel = "Maximal altitude r(tf)",
-                    seriestype = :scatter,
+                    obj_list;
+                    label="r(tf)",
+                    xlabel="Maximal thrust (Tmax)",
+                    ylabel="Maximal altitude r(tf)",
+                    seriestype=:scatter,
                 )
                 # plot multiple solutions
                 plot(sol0)
                 p = plot!(sol)
-                display(plot(pobj, p, layout = 2, reuse = false, size = (1000, 500)))
+                display(plot(pobj, p; layout=2, reuse=false, size=(1000, 500)))
             end
         end
     end
diff --git a/test/test_goddard_direct.jl b/test/test_goddard_direct.jl
index 056d9a9c..96f5241e 100644
--- a/test/test_goddard_direct.jl
+++ b/test/test_goddard_direct.jl
@@ -4,7 +4,7 @@ function test_goddard_direct()
     ocp, obj = Goddard()
 
     # initial guess (constant state and control functions)
-    init = (state = [1.01, 0.05, 0.8], control = 0.1, variable = 0.2)
-    sol = solve(ocp; grid_size = 10, display = false, init = init)
+    init = (state=[1.01, 0.05, 0.8], control=0.1, variable=0.2)
+    sol = solve(ocp; grid_size=10, display=false, init=init)
     @test objective(sol) ≈ obj atol = 5e-3
 end
diff --git a/test/test_goddard_indirect.jl b/test/test_goddard_indirect.jl
index 99cd238c..ca922026 100644
--- a/test/test_goddard_indirect.jl
+++ b/test/test_goddard_indirect.jl
@@ -65,7 +65,7 @@ function test_goddard_indirect()
         s[4] = H1(x1, p1)
         s[5] = H01(x1, p1)
         s[6] = g(x2)
-        s[7] = H0(xf, pf) # free tf
+        return s[7] = H0(xf, pf) # free tf
     end
 
     # tests
diff --git a/test/test_grid.jl b/test/test_grid.jl
index e58f6aad..79d4e6ea 100644
--- a/test/test_grid.jl
+++ b/test/test_grid.jl
@@ -4,66 +4,68 @@ function test_grid()
     ocp = Model()
     state!(ocp, 1)
     control!(ocp, 2)
-    time!(ocp, t0 = 0, tf = 1)
-    constraint!(ocp, :initial, lb = -1, ub = -1)
-    constraint!(ocp, :final, lb = 0, ub = 0)
-    constraint!(ocp, :control, lb = [0, 0], ub = [Inf, Inf])
+    time!(ocp; t0=0, tf=1)
+    constraint!(ocp, :initial; lb=-1, ub=-1)
+    constraint!(ocp, :final; lb=0, ub=0)
+    constraint!(ocp, :control; lb=[0, 0], ub=[Inf, Inf])
     dynamics!(ocp, (x, u) -> -x - u[1] + u[2])
     objective!(ocp, :lagrange, (x, u) -> (u[1] + u[2])^2)
 
     time_grid = [0, 0.1, 0.3, 0.6, 0.98, 0.99, 1]
-    sol = solve(ocp; time_grid = time_grid, display = false)
+    sol = solve(ocp; time_grid=time_grid, display=false)
     @test objective(sol) ≈ 0.309 rtol = 1e-2
 
     # 1. simple integrator min energy (dual control for test)
     ocp = Model()
     state!(ocp, 1)
     control!(ocp, 2)
-    time!(ocp, t0 = 0, tf = 1)
-    constraint!(ocp, :initial, lb = -1, ub = -1)
-    constraint!(ocp, :final, lb = 0, ub = 0)
-    constraint!(ocp, :control, lb = [0, 0], ub = [Inf, Inf])
+    time!(ocp; t0=0, tf=1)
+    constraint!(ocp, :initial; lb=-1, ub=-1)
+    constraint!(ocp, :final; lb=0, ub=0)
+    constraint!(ocp, :control; lb=[0, 0], ub=[Inf, Inf])
     dynamics!(ocp, (x, u) -> -x - u[1] + u[2])
     objective!(ocp, :lagrange, (x, u) -> (u[1] + u[2])^2)
-    sol0 = solve(ocp, print_level = 0)
+    sol0 = solve(ocp; print_level=0)
 
     # solve with explicit and non uniform time grid
     @testset verbose = true showtiming = true ":explicit_grid" begin
         time_grid = LinRange(0, 1, CTDirect.__grid_size() + 1)
-        sol = solve(ocp, time_grid = time_grid, print_level = 0)
+        sol = solve(ocp; time_grid=time_grid, print_level=0)
         @test objective(sol) == objective(sol0)
         @test iterations(sol) == iterations(sol0)
     end
 
     @testset verbose = true showtiming = true ":non_uniform_grid" begin
         time_grid = [0, 0.1, 0.3, 0.6, 0.98, 0.99, 1]
-        sol = solve(ocp, time_grid = time_grid, print_level = 0)
+        sol = solve(ocp; time_grid=time_grid, print_level=0)
         @test objective(sol) ≈ 0.309 rtol = 1e-2
     end
 
     # 2. integrator free times
-    ocp = Model(variable = true)
+    ocp = Model(; variable=true)
     state!(ocp, 2)
     control!(ocp, 1)
     variable!(ocp, 2)
-    time!(ocp, ind0 = 1, indf = 2)
-    constraint!(ocp, :initial, lb = [0, 0], ub = [0, 0])
-    constraint!(ocp, :final, lb = [1, 0], ub = [1, 0])
-    constraint!(ocp, :control, lb = -1, ub = 1)
-    constraint!(ocp, :variable, lb = [0.1, 0.1], ub = [10, 10])
-    constraint!(ocp, :variable, f = v -> v[2] - v[1], lb = 0.1, ub = Inf)
+    time!(ocp; ind0=1, indf=2)
+    constraint!(ocp, :initial; lb=[0, 0], ub=[0, 0])
+    constraint!(ocp, :final; lb=[1, 0], ub=[1, 0])
+    constraint!(ocp, :control; lb=-1, ub=1)
+    constraint!(ocp, :variable; lb=[0.1, 0.1], ub=[10, 10])
+    constraint!(ocp, :variable; f=v -> v[2] - v[1], lb=0.1, ub=Inf)
     dynamics!(ocp, (x, u, v) -> [x[2], u])
     objective!(ocp, :mayer, (x0, xf, v) -> v[1], :max)
-    sol0 = solve(ocp, print_level = 0)
+    sol0 = solve(ocp; print_level=0)
 
     @testset verbose = true showtiming = true ":explicit_grid" begin
-        sol = solve(ocp, time_grid = LinRange(0, 1, CTDirect.__grid_size() + 1), print_level = 0)
+        sol = solve(
+            ocp; time_grid=LinRange(0, 1, CTDirect.__grid_size() + 1), print_level=0
+        )
         @test objective(sol) == objective(sol0)
         @test iterations(sol) == iterations(sol0)
     end
 
     @testset verbose = true showtiming = true ":non_uniform_grid" begin
-        sol = solve(ocp, time_grid = [0, 0.1, 0.3, 0.5, 0.6, 0.8, 0.95, 1], print_level = 0)
+        sol = solve(ocp; time_grid=[0, 0.1, 0.3, 0.5, 0.6, 0.8, 0.95, 1], print_level=0)
         #plot(sol, show=true)
         @test objective(sol) ≈ 7.96 rtol = 1e-2
     end
@@ -82,16 +84,18 @@ function test_grid()
         ẋ(t) == [v(t), u(t)]
         ∫(u(t)^2) → min
     end
-    sol0 = solve(ocpT2, print_level = 0)
+    sol0 = solve(ocpT2; print_level=0)
 
     @testset verbose = true showtiming = true ":explicit_grid" begin
-        sol = solve(ocpT2, time_grid = LinRange(0, 1, CTDirect.__grid_size() + 1), print_level = 0)
+        sol = solve(
+            ocpT2; time_grid=LinRange(0, 1, CTDirect.__grid_size() + 1), print_level=0
+        )
         @test objective(sol) == objective(sol0)
         @test iterations(sol) == iterations(sol0)
     end
 
     @testset verbose = true showtiming = true ":non_uniform_grid" begin
-        sol = solve(ocpT2, time_grid = [0, 0.3, 1, 1.9, 2], print_level = 0)
+        sol = solve(ocpT2; time_grid=[0, 0.3, 1, 1.9, 2], print_level=0)
         @test objective(sol) ≈ 2.43 rtol = 1e-2
     end
 end
diff --git a/test/test_initial_guess.jl b/test/test_initial_guess.jl
index 5a2e308b..1fd1370b 100644
--- a/test/test_initial_guess.jl
+++ b/test/test_initial_guess.jl
@@ -28,7 +28,7 @@ function test_initial_guess()
         ẋ(t) == [x₂(t), u(t)]
         tf → min
     end
-    sol0 = solve(ocp, print_level = 0)
+    sol0 = solve(ocp; print_level=0)
 
     # constant initial guess
     x_const = [0.5, 0.2]
@@ -49,87 +49,75 @@ function test_initial_guess()
 
     # 1.a default initial guess
     @testset verbose = true showtiming = true ":default_init_no_arg" begin
-        sol = solve(ocp, print_level = 0, max_iter = maxiter)
+        sol = solve(ocp; print_level=0, max_iter=maxiter)
         @test(check_xf(sol, [0.1, 0.1]) && check_uf(sol, 0.1) && check_v(sol, 0.1))
     end
     @testset verbose = true showtiming = true ":default_init_()" begin
-        sol = solve(ocp, print_level = 0, init = (), max_iter = maxiter)
+        sol = solve(ocp; print_level=0, init=(), max_iter=maxiter)
         @test(check_xf(sol, [0.1, 0.1]) && check_uf(sol, 0.1) && check_v(sol, 0.1))
     end
     @testset verbose = true showtiming = true ":default_init_nothing" begin
-        sol = solve(ocp, print_level = 0, init = nothing, max_iter = maxiter)
+        sol = solve(ocp; print_level=0, init=nothing, max_iter=maxiter)
         @test(check_xf(sol, [0.1, 0.1]) && check_uf(sol, 0.1) && check_v(sol, 0.1))
     end
 
     # 1.b constant initial guess
     @testset verbose = true showtiming = true ":constant_x" begin
-        sol = solve(ocp, print_level = 0, init = (state = x_const,), max_iter = maxiter)
+        sol = solve(ocp; print_level=0, init=(state=x_const,), max_iter=maxiter)
         @test(check_xf(sol, x_const))
     end
     @testset verbose = true showtiming = true ":constant_u" begin
-        sol = solve(ocp, print_level = 0, init = (control = u_const,), max_iter = maxiter)
+        sol = solve(ocp; print_level=0, init=(control=u_const,), max_iter=maxiter)
         @test(check_uf(sol, u_const))
     end
     @testset verbose = true showtiming = true ":constant_v" begin
-        sol = solve(ocp, print_level = 0, init = (variable = v_const,), max_iter = maxiter)
+        sol = solve(ocp; print_level=0, init=(variable=v_const,), max_iter=maxiter)
         @test(check_v(sol, v_const))
     end
     @testset verbose = true showtiming = true ":constant_xu" begin
         sol = solve(
-            ocp,
-            print_level = 0,
-            init = (state = x_const, control = u_const),
-            max_iter = maxiter,
+            ocp; print_level=0, init=(state=x_const, control=u_const), max_iter=maxiter
         )
         @test(check_xf(sol, x_const) && check_uf(sol, u_const))
     end
     @testset verbose = true showtiming = true ":constant_xv" begin
         sol = solve(
-            ocp,
-            print_level = 0,
-            init = (state = x_const, variable = v_const),
-            max_iter = maxiter,
+            ocp; print_level=0, init=(state=x_const, variable=v_const), max_iter=maxiter
         )
         @test(check_xf(sol, x_const) && check_v(sol, v_const))
     end
     @testset verbose = true showtiming = true ":constant_uv" begin
         sol = solve(
-            ocp,
-            print_level = 0,
-            init = (control = u_const, variable = v_const),
-            max_iter = maxiter,
+            ocp; print_level=0, init=(control=u_const, variable=v_const), max_iter=maxiter
         )
         @test(check_uf(sol, u_const) && check_v(sol, v_const))
     end
     @testset verbose = true showtiming = true ":constant_xuv" begin
         sol = solve(
-            ocp,
-            print_level = 0,
-            init = (state = x_const, control = u_const, variable = v_const),
-            max_iter = maxiter,
+            ocp;
+            print_level=0,
+            init=(state=x_const, control=u_const, variable=v_const),
+            max_iter=maxiter,
         )
         @test(check_xf(sol, x_const) && check_uf(sol, u_const) && check_v(sol, v_const))
     end
 
     # 1. functional initial guess
     @testset verbose = true showtiming = true ":functional_x" begin
-        sol = solve(ocp, print_level = 0, init = (state = x_func,), max_iter = maxiter)
+        sol = solve(ocp; print_level=0, init=(state=x_func,), max_iter=maxiter)
         @test(check_xf(sol, x_func(time_grid(sol)[end])))
     end
     @testset verbose = true showtiming = true ":functional_u" begin
-        sol = solve(ocp, print_level = 0, init = (control = u_func,), max_iter = maxiter)
+        sol = solve(ocp; print_level=0, init=(control=u_func,), max_iter=maxiter)
         @test(check_uf(sol, u_func(time_grid(sol)[end])))
     end
     @testset verbose = true showtiming = true ":functional_xu" begin
         sol = solve(
-            ocp,
-            print_level = 0,
-            init = (state = x_func, control = u_func),
-            max_iter = maxiter,
+            ocp; print_level=0, init=(state=x_func, control=u_func), max_iter=maxiter
         )
         @test(
             check_xf(sol, x_func(time_grid(sol)[end])) &&
-            check_uf(sol, u_func(time_grid(sol)[end]))
+                check_uf(sol, u_func(time_grid(sol)[end]))
         )
     end
 
@@ -137,55 +125,61 @@ function test_initial_guess()
     t_vec = [0, 0.1, v_const]
     @testset verbose = true showtiming = true ":vector_txu :constant_v" begin
         sol = solve(
-            ocp,
-            print_level = 0,
-            init = (time = t_vec, state = x_vec, control = u_vec, variable = v_const),
-            max_iter = maxiter,
+            ocp;
+            print_level=0,
+            init=(time=t_vec, state=x_vec, control=u_vec, variable=v_const),
+            max_iter=maxiter,
+        )
+        @test(
+            check_xf(sol, x_vec[end]) && check_uf(sol, u_vec[end]) && check_v(sol, v_const)
         )
-        @test(check_xf(sol, x_vec[end]) && check_uf(sol, u_vec[end]) && check_v(sol, v_const))
     end
     t_matrix = [0 0.1 v_const]
     @testset verbose = true showtiming = true ":matrix_t :vector_xu :constant_v" begin
         sol = solve(
-            ocp,
-            print_level = 0,
-            init = (time = t_matrix, state = x_vec, control = u_vec, variable = v_const),
-            max_iter = maxiter,
+            ocp;
+            print_level=0,
+            init=(time=t_matrix, state=x_vec, control=u_vec, variable=v_const),
+            max_iter=maxiter,
+        )
+        @test(
+            check_xf(sol, x_vec[end]) && check_uf(sol, u_vec[end]) && check_v(sol, v_const)
         )
-        @test(check_xf(sol, x_vec[end]) && check_uf(sol, u_vec[end]) && check_v(sol, v_const))
     end
     @testset verbose = true showtiming = true ":matrix_x :vector_tu :constant_v" begin
         sol = solve(
-            ocp,
-            print_level = 0,
-            init = (time = t_vec, state = x_matrix, control = u_vec, variable = v_const),
-            max_iter = maxiter,
+            ocp;
+            print_level=0,
+            init=(time=t_vec, state=x_matrix, control=u_vec, variable=v_const),
+            max_iter=maxiter,
+        )
+        @test(
+            check_xf(sol, x_vec[end]) && check_uf(sol, u_vec[end]) && check_v(sol, v_const)
         )
-        @test(check_xf(sol, x_vec[end]) && check_uf(sol, u_vec[end]) && check_v(sol, v_const))
     end
 
     # 1.e mixed initial guess
     @testset verbose = true showtiming = true ":vector_tx :functional_u :constant_v" begin
         sol = solve(
-            ocp,
-            print_level = 0,
-            init = (time = t_vec, state = x_vec, control = u_func, variable = v_const),
-            max_iter = maxiter,
+            ocp;
+            print_level=0,
+            init=(time=t_vec, state=x_vec, control=u_func, variable=v_const),
+            max_iter=maxiter,
         )
         @test(
             check_xf(sol, x_vec[end]) &&
-            check_uf(sol, u_func(time_grid(sol)[end])) &&
-            check_v(sol, v_const)
+                check_uf(sol, u_func(time_grid(sol)[end])) &&
+                check_v(sol, v_const)
         )
     end
 
     # 1.f warm start
     @testset verbose = true showtiming = true ":warm_start" begin
-        sol = solve(ocp, print_level = 0, init = sol0, max_iter = maxiter)
+        sol = solve(ocp; print_level=0, init=sol0, max_iter=maxiter)
         @test(
             check_xf(sol, state(sol)(time_grid(sol)[end])) &&
-            check_uf(sol, control(sol)(time_grid(sol)[end])) &&
-            check_v(sol, variable(sol))
+                check_uf(sol, control(sol)(time_grid(sol)[end])) &&
+                check_v(sol, variable(sol))
         )
     end
 end
diff --git a/test/test_objective.jl b/test/test_objective.jl
index 79cce76a..bce9a71a 100644
--- a/test/test_objective.jl
+++ b/test/test_objective.jl
@@ -2,67 +2,67 @@
 function test_objective()
 
     # min tf
-    ocp = Model(variable = true)
+    ocp = Model(; variable=true)
     state!(ocp, 2)
     control!(ocp, 1)
     variable!(ocp, 1)
-    time!(ocp, t0 = 0, indf = 1)
-    constraint!(ocp, :initial, lb = [0, 0], ub = [0, 0])
-    constraint!(ocp, :final, lb = [1, 0], ub = [1, 0])
-    constraint!(ocp, :control, lb = -1, ub = 1)
-    constraint!(ocp, :variable, lb = 0.1, ub = 10)
+    time!(ocp; t0=0, indf=1)
+    constraint!(ocp, :initial; lb=[0, 0], ub=[0, 0])
+    constraint!(ocp, :final; lb=[1, 0], ub=[1, 0])
+    constraint!(ocp, :control; lb=-1, ub=1)
+    constraint!(ocp, :variable; lb=0.1, ub=10)
     dynamics!(ocp, (x, u, v) -> [x[2], u])
     objective!(ocp, :mayer, (x0, xf, v) -> v)
 
     @testset verbose = true showtiming = true ":min_tf :mayer" begin
-        sol = solve(ocp, print_level = 0, tol = 1e-12)
+        sol = solve(ocp; print_level=0, tol=1e-12)
         @test objective(sol) ≈ 2.0 rtol = 1e-2
     end
 
     # min tf (lagrange)
-    ocp = Model(variable = true)
+    ocp = Model(; variable=true)
     state!(ocp, 2)
     control!(ocp, 1)
     variable!(ocp, 1)
-    time!(ocp, t0 = 0, indf = 1)
-    constraint!(ocp, :initial, lb = [0, 0], ub = [0, 0])
-    constraint!(ocp, :final, lb = [1, 0], ub = [1, 0])
-    constraint!(ocp, :control, lb = -1, ub = 1)
-    constraint!(ocp, :variable, lb = 0.1, ub = 10)
+    time!(ocp; t0=0, indf=1)
+    constraint!(ocp, :initial; lb=[0, 0], ub=[0, 0])
+    constraint!(ocp, :final; lb=[1, 0], ub=[1, 0])
+    constraint!(ocp, :control; lb=-1, ub=1)
+    constraint!(ocp, :variable; lb=0.1, ub=10)
     dynamics!(ocp, (x, u, v) -> [x[2], u])
     objective!(ocp, :lagrange, (x, u, v) -> 1)
 
     @testset verbose = true showtiming = true ":min_tf :lagrange" begin
-        sol = solve(ocp, print_level = 0, tol = 1e-12)
+        sol = solve(ocp; print_level=0, tol=1e-12)
         @test objective(sol) ≈ 2.0 rtol = 1e-2
     end
 
     # max t0 (free t0 and tf)
-    ocp = Model(variable = true)
+    ocp = Model(; variable=true)
     state!(ocp, 2)
     control!(ocp, 1)
     variable!(ocp, 2)
-    time!(ocp, ind0 = 1, indf = 2)
-    constraint!(ocp, :initial, lb = [0, 0], ub = [0, 0])
-    constraint!(ocp, :final, lb = [1, 0], ub = [1, 0])
-    constraint!(ocp, :control, lb = -1, ub = 1)
-    constraint!(ocp, :variable, lb = [0.1, 0.1], ub = [10, 10])
-    constraint!(ocp, :variable, f = v -> v[2] - v[1], lb = 0.1, ub = Inf)
+    time!(ocp; ind0=1, indf=2)
+    constraint!(ocp, :initial; lb=[0, 0], ub=[0, 0])
+    constraint!(ocp, :final; lb=[1, 0], ub=[1, 0])
+    constraint!(ocp, :control; lb=-1, ub=1)
+    constraint!(ocp, :variable; lb=[0.1, 0.1], ub=[10, 10])
+    constraint!(ocp, :variable; f=v -> v[2] - v[1], lb=0.1, ub=Inf)
     dynamics!(ocp, (x, u, v) -> [x[2], u])
     objective!(ocp, :mayer, (x0, xf, v) -> v[1], :max)
 
     @testset verbose = true showtiming = true ":max_t0" begin
-        sol = solve(ocp, print_level = 0, tol = 1e-12)
+        sol = solve(ocp; print_level=0, tol=1e-12)
         @test objective(sol) ≈ 8.0 rtol = 1e-2
     end
 
     @testset verbose = true showtiming = true ":max_t0 :explicit_grid" begin
-        sol = solve(ocp, time_grid = LinRange(0, 1, 101), print_level = 0, tol = 1e-12)
+        sol = solve(ocp; time_grid=LinRange(0, 1, 101), print_level=0, tol=1e-12)
         @test objective(sol) ≈ 8.0 rtol = 1e-2
     end
 
     @testset verbose = true showtiming = true ":max_t0 :non_uniform_grid" begin
-        sol = solve(ocp, time_grid = [0, 0.1, 0.6, 0.95, 1], print_level = 0, tol = 1e-12)
+        sol = solve(ocp; time_grid=[0, 0.1, 0.6, 0.95, 1], print_level=0, tol=1e-12)
         @test objective(sol) ≈ 7.48 rtol = 1e-2
     end
 end