diff --git a/docs/problems.jl b/docs/problems.jl
index ba61045c..fa7596a0 100644
--- a/docs/problems.jl
+++ b/docs/problems.jl
@@ -32,7 +32,7 @@ function generate_documentation(
     using OptimalControlProblems    # to access the Beam model
     using OptimalControl            # to import the OptimalControl model
     using NLPModelsIpopt            # to solve the model with Ipopt
-    using DataFrames                # to store data
+    import DataFrames: DataFrame    # to store data
     using NLPModels                 # to retrieve data from the NLP solution
     using Plots                     # to plot the trajectories
     using Plots.PlotMeasures        # for leftmargin, bottommargin
@@ -100,7 +100,7 @@ function generate_documentation(
     push!(data_pb,
         (
             Problem=:$PROBLEM,
-            Grid_Size=OptimalControlProblems.metadata[:$PROBLEM][:N],
+            Grid_Size=metadata[:$PROBLEM][:N],
             Variables=get_nvar(nlp_oc),
             Constraints=get_ncon(nlp_oc),
         )
@@ -125,8 +125,8 @@ function generate_documentation(
     Visualise states, costates, and controls for the OptimalControl solution:
 
     ```@example main
-    x_vars = OptimalControlProblems.metadata[:$PROBLEM][:state_name]
-    u_vars = OptimalControlProblems.metadata[:$PROBLEM][:control_name]
+    x_vars = metadata[:$PROBLEM][:state_name]
+    u_vars = metadata[:$PROBLEM][:control_name]
 
     n = length(x_vars) # number of states
     m = length(u_vars) # number of controls
@@ -319,7 +319,7 @@ function generate_documentation(
     ```
 
     ```@example main
-    v_vars = OptimalControlProblems.metadata[:$PROBLEM][:variable_name]
+    v_vars = metadata[:$PROBLEM][:variable_name]
 
     function L2_norm(T, X)
         # T and X are supposed to be one dimensional
diff --git a/docs/src/dev-add.md b/docs/src/dev-add.md
index eb31e7b8..54f2589d 100644
--- a/docs/src/dev-add.md
+++ b/docs/src/dev-add.md
@@ -21,13 +21,13 @@ new_problem_meta = OrderedDict(
 
     For more details about the metadata, see the [MetaData](@ref problems-introduction-metadata) section.
 
-**2.** Define the **OptimalControl** model of the problem in a separate file in the `ext/OptimalControlModels` directory.
+**2.** Define the DOCP **OptimalControl** model of the problem in a file named `new_problem.jl` in the `ext/OptimalControlModels` directory.
 
 ```julia
 """
-    Description of the new problem
+    Documentation of the method
 """
-function OptimalControlProblems.new_problem(::OptimalControlBackend; N::Int=steps_number_data(:new_problem))
+function OptimalControlProblems.new_problem(::OptimalControlBackend, description::Symbol...; N::Int=steps_number_data(:new_problem), kwargs...)
 
     # if tf is fixed
     tf = final_time_data(:new_problem)
@@ -42,30 +42,32 @@ function OptimalControlProblems.new_problem(::OptimalControlBackend; N::Int=step
     init = () 
 
     # DOCP and NLP
-    docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
+    docp = direct_transcription(ocp, description...; init=init, grid_size=N, disc_method=:trapeze, kwargs...)
 
     return docp
 
 end
 ```
 
-**3.** Define the **JuMP** model of the problem in a new file in the `ext/JuMPModels` directory.
+**3.** Define the NLP **JuMP** model of the problem in a file named `new_problem.jl` in the `ext/JuMPModels` directory.
 
 ```julia
 """
-    Description of the new problem
+    Documentation of the method
 """
-function OptimalControlProblems.new_problem(::JuMPBackend; N::Int=steps_number_data(:new_problem))
+function OptimalControlProblems.new_problem(::JuMPBackend, args...; N::Int=steps_number_data(:new_problem), kwargs...)
 
     # if tf is fixed
     tf = final_time_data(:new_problem)
 
     # model
-    nlp = JuMP.Model()
+    model = JuMP.Model(args...; kwargs...)
 
     # Define the problem here
     # ...
 
-    return nlp
+    return model
 end
 ```
+
+**3.** Describe the problem in a file named `new_problem.jl` in the `ext/Descriptions` directory.
\ No newline at end of file
diff --git a/docs/src/index.md b/docs/src/index.md
index dec51526..22c5c109 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -15,6 +15,15 @@ using Pkg
 Pkg.add("OptimalControlProblems")
 ```
 
+## Credits (not exhaustive!)
+
+- [Nico77310](https://github.com/Nico77310)
+- [0Yassine0](https://github.com/0Yassine0)
+- [frapac](https://github.com/frapac)
+- [BaptisteCbl](https://github.com/BaptisteCbl)
+- [COPS: Large-Scale Optimization Problems](https://www.mcs.anl.gov/~more/cops) and [COPSBenchmark.jl](github.com/MadNLP/COPSBenchmark.jl)
+- [BOCOP - A collection of examples](https://project.inria.fr/bocop/files/2017/05/Examples-BOCOP.pdf)
+
 ## Reproducibility
 
 ```@raw html
diff --git a/docs/src/problems-introduction.md b/docs/src/problems-introduction.md
index 049e2b1d..836a1754 100644
--- a/docs/src/problems-introduction.md
+++ b/docs/src/problems-introduction.md
@@ -81,15 +81,14 @@ Depth = 1
 For each problem, additional data is provided in the [MetaData](https://github.com/control-toolbox/OptimalControlProblems.jl/tree/main/ext/MetaData) directory:
 
 ```@docs; canonical=false
-OptimalControlProblems.metadata
+metadata
 ```
 
-To list all metadata, use `OptimalControlProblems.metadata`.  
-To access the metadata of a specific problem, for example `chain`, run:
+To list all metadata, use `metadata`. To access the metadata of a specific problem, for example `chain`, run:
 
 ```@example main
 using OptimalControlProblems
-OptimalControlProblems.metadata[:chain]
+metadata[:chain]
 ```
 
 ## Problems characteristics
@@ -113,7 +112,7 @@ We detail below the characteristics of the optimal control problems (OCPs) and t
 
 ```@example main
 using NLPModels                 # to get the number of variables and constraints
-using DataFrames
+import DataFrames: DataFrame    # to store data
 using OptimalControl
 
 data_ocp = DataFrame(           # to store data of the OCPs
@@ -183,7 +182,7 @@ for problem in problems()
     )
 
     #
-    N = OptimalControlProblems.metadata[problem][:N] # get default number of steps
+    N = metadata[problem][:N] # get default number of steps
 
     push!(data_nlp,
         (
diff --git a/docs/src/tutorial-get.md b/docs/src/tutorial-get.md
index 8d6a39f4..1fe76afd 100644
--- a/docs/src/tutorial-get.md
+++ b/docs/src/tutorial-get.md
@@ -24,16 +24,19 @@ The `nlp` model represents the nonlinear programming problem (NLP) obtained afte
 
 !!! note
 
-    You also have access to the DOCP model, which corresponds to the discretised optimal control problem. Roughly speaking, the DOCP model is the union of the NLP and OCP models. For more details, see this [tutorial](@extref Tutorials Discretization-and-NLP-problem) or the documentation of [`CTDirect.DOCP`](@extref).
+    You also have access to the DOCP model, which corresponds to the discretised optimal control problem. Roughly speaking, the DOCP model is the union of the NLP and OCP models. For more details, see this [tutorial](@extref Tutorials Discretization-and-NLP-problem) or the documentation of [`CTDirect.DOCP`](@extref). To get the OCP model:
 
-### OCP
+    ```julia
+    ocp = ocp_model(docp)
+    ```
 
-You also have access to the OCP model, which corresponds to the optimal control problem.
+!!! note
 
-```@example main_oc
-ocp = ocp_model(docp)
-nothing # hide
-```
+    You can pass any `description` and `kwargs` of [`CTDirect.direct_transcription`](@extref) to the `beam` problem or any other.
+
+    ```julia
+    docp = beam(OptimalControlBackend(), :madnlp; grid_size=100, disc_method=:euler)
+    ```
 
 ### Number of variables, constraints, and nonzeros
 
@@ -72,7 +75,7 @@ println("nnzh = ", nnzh)
 The number of steps $N$ is stored in the metadata:
 
 ```@example main_oc
-OptimalControlProblems.metadata[:beam][:N]
+metadata[:beam][:N]
 ```
 
 !!! note
@@ -108,4 +111,9 @@ nlp = beam(JuMPBackend())
 ```
 
 !!! note
-    For details on how to interact with the JuMP model, see the [JuMP documentation](https://jump.dev/JuMP.jl).
+    For details on how to interact with the JuMP model, see the [JuMP documentation](https://jump.dev/JuMP.jl). In particular, you can pass any arguments and keyword arguments of [`JuMP.Model`](@extref) to the `beam` problem or any other.
+
+    ```julia
+    using Ipopt
+    nlp = beam(JuMPBackend(), Ipopt.Optimizer; add_bridges=true)
+    ``` 
diff --git a/docs/src/tutorial-solve.md b/docs/src/tutorial-solve.md
index 84657786..dd5bf529 100644
--- a/docs/src/tutorial-solve.md
+++ b/docs/src/tutorial-solve.md
@@ -145,8 +145,8 @@ println("p(tf) = ", p(tf))
 We can add the state, costate, and control to the plot:
 
 ```@example main
-n = length(OptimalControlProblems.metadata[problem][:state_name])   # dimension of the state
-m = length(OptimalControlProblems.metadata[problem][:control_name]) # dimension of the control
+n = length(metadata[problem][:state_name])   # dimension of the state
+m = length(metadata[problem][:control_name]) # dimension of the control
 
 for i in 1:n # state
     plot!(plt[i], t, t -> x(t)[i]; color=2, linestyle=:dash, label=:none)
diff --git a/ext/JuMPModels.jl b/ext/JuMPModels.jl
index c33c4201..a96b95bd 100644
--- a/ext/JuMPModels.jl
+++ b/ext/JuMPModels.jl
@@ -24,7 +24,7 @@ Compute the discretised time grid for a given optimal control problem solved wit
 
 # Arguments
 
-- `problem::Symbol`: The name of the problem as defined in `OptimalControlProblems.metadata`.
+- `problem::Symbol`: The name of the problem as defined in `metadata`.
 - `model::JuMP.GenericModel`: The JuMP model containing the problem solution.
 
 # Returns
@@ -41,18 +41,18 @@ julia> tgrid = OptimalControlProblems.time_grid(:my_problem, model)
 function OptimalControlProblems.time_grid(problem::Symbol, model::JuMP.GenericModel)
 
     # get N
-    x_vars = OptimalControlProblems.metadata[problem][:state_name]
+    x_vars = metadata[problem][:state_name]
     x_jp_var = JuMP.value.(model[Symbol(x_vars[1])])
     N = length(x_jp_var) - 1
 
     ## time grid: we assume that t0 = 0
-    time_data, time_value_or_index = OptimalControlProblems.metadata[problem][:final_time]
+    time_data, time_value_or_index = metadata[problem][:final_time]
 
     t0 = 0
     tf = if time_data == :fixed
         time_value_or_index
     elseif time_data == :free
-        v_vars = OptimalControlProblems.metadata[problem][:variable_name]
+        v_vars = metadata[problem][:variable_name]
         value.(model[Symbol(v_vars[time_value_or_index])])
     else
         error("the final time must be :fixed or :free, not: ", time_data)
@@ -69,7 +69,7 @@ Extract and interpolate the state trajectory from a JuMP model of an optimal con
 
 # Arguments
 
-- `problem::Symbol`: The name of the problem as defined in `OptimalControlProblems.metadata`.
+- `problem::Symbol`: The name of the problem as defined in `metadata`.
 - `model::JuMP.GenericModel`: The JuMP model containing the problem solution.
 
 # Returns
@@ -92,7 +92,7 @@ function OptimalControlProblems.state(problem::Symbol, model::JuMP.GenericModel)
     N = length(T) - 1
 
     # get dimension
-    state_names = OptimalControlProblems.metadata[problem][:state_name]
+    state_names = metadata[problem][:state_name]
     dim_x = length(state_names)
 
     # get state from the model
@@ -120,7 +120,7 @@ Extract and interpolate the control trajectory from a JuMP model of an optimal c
 
 # Arguments
 
-- `problem::Symbol`: The name of the problem as defined in `OptimalControlProblems.metadata`.
+- `problem::Symbol`: The name of the problem as defined in `metadata`.
 - `model::JuMP.GenericModel`: The JuMP model containing the problem solution.
 
 # Returns
@@ -143,7 +143,7 @@ function OptimalControlProblems.control(problem::Symbol, model::JuMP.GenericMode
     N = length(T) - 1
 
     # get dimension
-    control_names = OptimalControlProblems.metadata[problem][:control_name]
+    control_names = metadata[problem][:control_name]
     dim_u = length(control_names)
 
     # get control from the model
@@ -171,7 +171,7 @@ Extract and interpolate the costate trajectory (dual variables associated with s
 
 # Arguments
 
-- `problem::Symbol`: The name of the problem as defined in `OptimalControlProblems.metadata`.
+- `problem::Symbol`: The name of the problem as defined in `metadata`.
 - `model::JuMP.GenericModel`: The JuMP model containing the problem solution.
 
 # Returns
@@ -194,7 +194,7 @@ function OptimalControlProblems.costate(problem::Symbol, model::JuMP.GenericMode
     N = length(T) - 1
 
     # get dimension
-    costate_names = OptimalControlProblems.metadata[problem][:costate_name]
+    costate_names = metadata[problem][:costate_name]
     dim_x = length(costate_names)
 
     # get costate from the model
@@ -226,7 +226,7 @@ Extract scalar or vector decision variables (such as final time when free) from
 
 # Arguments
 
-- `problem::Symbol`: The name of the problem as defined in `OptimalControlProblems.metadata`.
+- `problem::Symbol`: The name of the problem as defined in `metadata`.
 - `model::JuMP.GenericModel`: The JuMP model containing the problem solution.
 
 # Returns
@@ -244,7 +244,7 @@ julia> v = OptimalControlProblems.variable(:my_problem, model)
 ```
 """
 function OptimalControlProblems.variable(problem::Symbol, model::JuMP.GenericModel)
-    variable_names = OptimalControlProblems.metadata[problem][:variable_name]
+    variable_names = metadata[problem][:variable_name]
 
     if isnothing(variable_names)
         return nothing
diff --git a/ext/JuMPModels/beam.jl b/ext/JuMPModels/beam.jl
index 5333e2a3..4fcceca1 100644
--- a/ext/JuMPModels/beam.jl
+++ b/ext/JuMPModels/beam.jl
@@ -27,14 +27,14 @@ julia> model = OptimalControlProblems.beam(JuMPBackend(); N=100)
 
 - Problem formulation available at: https://github.com/control-toolbox/bocop/tree/main/bocop
 """
-function OptimalControlProblems.beam(::JuMPBackend; N::Int=steps_number_data(:beam))
+function OptimalControlProblems.beam(::JuMPBackend, args...; N::Int=steps_number_data(:beam), kwargs...)
 
     # parameters
     tf = final_time_data(:beam)
     step = tf / N # t0 = 0
 
     # model
-    model = JuMP.Model()
+    model = JuMP.Model(args...; kwargs...)
 
     # variables and initial guess
     @variables(
diff --git a/ext/JuMPModels/bioreactor.jl b/ext/JuMPModels/bioreactor.jl
index dc894225..363754f5 100644
--- a/ext/JuMPModels/bioreactor.jl
+++ b/ext/JuMPModels/bioreactor.jl
@@ -31,7 +31,7 @@ julia> model = OptimalControlProblems.bioreactor(JuMPBackend(); N=100)
 - [control-toolbox/bocop](https://github.com/control-toolbox/bocop/tree/main/bocop)
 """
 function OptimalControlProblems.bioreactor(
-    ::JuMPBackend; N::Int=steps_number_data(:bioreactor)
+    ::JuMPBackend, args...; N::Int=steps_number_data(:bioreactor), kwargs...
 )
 
     # parameters
@@ -46,7 +46,7 @@ function OptimalControlProblems.bioreactor(
     T = final_time_data(:bioreactor)
 
     # model
-    model = JuMP.Model()
+    model = JuMP.Model(args...; kwargs...)
 
     # variables and initial guess
     @variables(
diff --git a/ext/JuMPModels/cart_pendulum.jl b/ext/JuMPModels/cart_pendulum.jl
index cd848ce0..517b42df 100644
--- a/ext/JuMPModels/cart_pendulum.jl
+++ b/ext/JuMPModels/cart_pendulum.jl
@@ -28,7 +28,7 @@ julia> model = OptimalControlProblems.cart_pendulum(JuMPBackend(); N=200)
 - [Cart–Pendulum Optimal Control Problem](https://arxiv.org/pdf/2303.16746)
 """
 function OptimalControlProblems.cart_pendulum(
-    ::JuMPBackend; N::Int=steps_number_data(:cart_pendulum)
+    ::JuMPBackend, args...; N::Int=steps_number_data(:cart_pendulum), kwargs...
 )
 
     # parameters
@@ -42,7 +42,7 @@ function OptimalControlProblems.cart_pendulum(
     max_v = 2
 
     # model
-    model = JuMP.Model()
+    model = JuMP.Model(args...; kwargs...)
 
     # variables and initial guess
     @variables(
diff --git a/ext/JuMPModels/chain.jl b/ext/JuMPModels/chain.jl
index 0da52413..90921b5b 100644
--- a/ext/JuMPModels/chain.jl
+++ b/ext/JuMPModels/chain.jl
@@ -27,7 +27,7 @@ julia> model = OptimalControlProblems.chain(JuMPBackend(); N=300)
 
 - [COPS Benchmark Problems – Hanging Chain](https://www.mcs.anl.gov/~more/cops/)
 """
-function OptimalControlProblems.chain(::JuMPBackend; N::Int=steps_number_data(:chain))
+function OptimalControlProblems.chain(::JuMPBackend, args...; N::Int=steps_number_data(:chain), kwargs...)
 
     # parameters
     L = 4
@@ -39,7 +39,7 @@ function OptimalControlProblems.chain(::JuMPBackend; N::Int=steps_number_data(:c
     tmin = b > a ? 1 / 4 : 3 / 4
 
     # model
-    model = JuMP.Model()
+    model = JuMP.Model(args...; kwargs...)
 
     # time
     @expressions(
diff --git a/ext/JuMPModels/dielectrophoretic_particle.jl b/ext/JuMPModels/dielectrophoretic_particle.jl
index 42c649ad..cf98a71d 100644
--- a/ext/JuMPModels/dielectrophoretic_particle.jl
+++ b/ext/JuMPModels/dielectrophoretic_particle.jl
@@ -29,7 +29,7 @@ julia> model = OptimalControlProblems.dielectrophoretic_particle(JuMPBackend();
   IEEE Transactions on Automatic Control, 51(7), 1100–1114.
 """
 function OptimalControlProblems.dielectrophoretic_particle(
-    ::JuMPBackend; N::Int=steps_number_data(:dielectrophoretic_particle)
+    ::JuMPBackend, args...; N::Int=steps_number_data(:dielectrophoretic_particle), kwargs...
 )
 
     # parameters
@@ -39,7 +39,7 @@ function OptimalControlProblems.dielectrophoretic_particle(
     c = 1
 
     # model
-    model = JuMP.Model()
+    model = JuMP.Model(args...; kwargs...)
 
     # state, control and variable (final time)
     @variable(model, x[0:N], start = 1)
diff --git a/ext/JuMPModels/double_oscillator.jl b/ext/JuMPModels/double_oscillator.jl
index 3e72797a..c7c60c1c 100644
--- a/ext/JuMPModels/double_oscillator.jl
+++ b/ext/JuMPModels/double_oscillator.jl
@@ -29,7 +29,7 @@ julia> model = OptimalControlProblems.double_oscillator(JuMPBackend(); N=200)
   IFAC-PapersOnLine, 51(2), 49–54.
 """
 function OptimalControlProblems.double_oscillator(
-    ::JuMPBackend; N::Int=steps_number_data(:double_oscillator)
+    ::JuMPBackend, args...; N::Int=steps_number_data(:double_oscillator), kwargs...
 )
 
     # parameters
@@ -41,7 +41,7 @@ function OptimalControlProblems.double_oscillator(
     tf = final_time_data(:double_oscillator)
 
     # model
-    model = Model()
+    model = JuMP.Model(args...; kwargs...)
 
     # state, control and initial guess
     @variables(
diff --git a/ext/JuMPModels/ducted_fan.jl b/ext/JuMPModels/ducted_fan.jl
index 1878104f..9fbc16c1 100644
--- a/ext/JuMPModels/ducted_fan.jl
+++ b/ext/JuMPModels/ducted_fan.jl
@@ -29,7 +29,7 @@ julia> model = OptimalControlProblems.ducted_fan(JuMPBackend(); N=100)
   Optimal Control Applications and Methods, 30(6), 537–561. [GP2009]
 """
 function OptimalControlProblems.ducted_fan(
-    ::JuMPBackend; N::Int=steps_number_data(:ducted_fan)
+    ::JuMPBackend, args...; N::Int=steps_number_data(:ducted_fan), kwargs...
 )
 
     # parameters
@@ -40,7 +40,7 @@ function OptimalControlProblems.ducted_fan(
     μ = 1000
 
     # model
-    model = Model()
+    model = JuMP.Model(args...; kwargs...)
 
     # state, control, variable (final time) and initial guess
     @variable(model, x₁[0:N], start = 0.1)
diff --git a/ext/JuMPModels/electric_vehicle.jl b/ext/JuMPModels/electric_vehicle.jl
index 31611bca..18716527 100644
--- a/ext/JuMPModels/electric_vehicle.jl
+++ b/ext/JuMPModels/electric_vehicle.jl
@@ -28,7 +28,7 @@ julia> model = OptimalControlProblems.electric_vehicle(JuMPBackend(); N=100)
 - Petit, N., & Sciarretta, A. (2011). *Optimal drive of electric vehicles using an inversion-based trajectory generation approach.* IFAC Proceedings Volumes, 44(1), 14519–14526. [PS2011]
 """
 function OptimalControlProblems.electric_vehicle(
-    ::JuMPBackend; N::Int=steps_number_data(:electric_vehicle)
+    ::JuMPBackend, args...; N::Int=steps_number_data(:electric_vehicle), kwargs...
 )
 
     # parameters
@@ -42,7 +42,7 @@ function OptimalControlProblems.electric_vehicle(
     α0, α1, α2, α3 = (3, 0.4, -1, 0.1)
 
     # model
-    model = Model()
+    model = JuMP.Model(args...; kwargs...)
 
     # state, control and initial guess
     @variable(model, x[0:N], start = 0.1)
diff --git a/ext/JuMPModels/glider.jl b/ext/JuMPModels/glider.jl
index fbd67e56..ca471d2b 100644
--- a/ext/JuMPModels/glider.jl
+++ b/ext/JuMPModels/glider.jl
@@ -27,7 +27,7 @@ julia> model = OptimalControlProblems.glider(JuMPBackend(); N=100)
 
 - Hang Glider Problem formulation as in: https://www.mcs.anl.gov/~more/cops/
 """
-function OptimalControlProblems.glider(::JuMPBackend; N::Int=steps_number_data(:glider))
+function OptimalControlProblems.glider(::JuMPBackend, args...; N::Int=steps_number_data(:glider), kwargs...)
 
     # parameters
     x_0 = 0
@@ -49,7 +49,7 @@ function OptimalControlProblems.glider(::JuMPBackend; N::Int=steps_number_data(:
     cL_max = 1.4
 
     # model
-    model = JuMP.Model()
+    model = JuMP.Model(args...; kwargs...)
 
     # state, control, variable (final time) and initial guess
     @variables(
diff --git a/ext/JuMPModels/insurance.jl b/ext/JuMPModels/insurance.jl
index 5a50bb1c..33d8a49a 100644
--- a/ext/JuMPModels/insurance.jl
+++ b/ext/JuMPModels/insurance.jl
@@ -28,7 +28,7 @@ julia> model = OptimalControlProblems.insurance(JuMPBackend(); N=100)
 - Problem formulation available at: https://github.com/control-toolbox/bocop/tree/main/bocop
 """
 function OptimalControlProblems.insurance(
-    ::JuMPBackend; N::Int=steps_number_data(:insurance)
+    ::JuMPBackend, args...; N::Int=steps_number_data(:insurance), kwargs...
 )
 
     # parameters
@@ -43,7 +43,7 @@ function OptimalControlProblems.insurance(
     α = 4
 
     # model
-    model = JuMP.Model()
+    model = JuMP.Model(args...; kwargs...)
 
     # state, control and initial guess
     @variables(
diff --git a/ext/JuMPModels/jackson.jl b/ext/JuMPModels/jackson.jl
index b95bb4a4..61199ae2 100644
--- a/ext/JuMPModels/jackson.jl
+++ b/ext/JuMPModels/jackson.jl
@@ -28,7 +28,7 @@ julia> model = OptimalControlProblems.jackson(JuMPBackend(); N=100)
 
 - Problem formulation available at: https://github.com/control-toolbox/bocop/tree/main/bocop
 """
-function OptimalControlProblems.jackson(::JuMPBackend; N::Int=steps_number_data(:jackson))
+function OptimalControlProblems.jackson(::JuMPBackend, args...; N::Int=steps_number_data(:jackson), kwargs...)
 
     # parameters
     tf = final_time_data(:jackson)
@@ -37,7 +37,7 @@ function OptimalControlProblems.jackson(::JuMPBackend; N::Int=steps_number_data(
     k3 = 1
 
     # model
-    model = JuMP.Model()
+    model = JuMP.Model(args...; kwargs...)
 
     @variables(
         model,
diff --git a/ext/JuMPModels/moonlander.jl b/ext/JuMPModels/moonlander.jl
index 6496e0f7..be89759a 100644
--- a/ext/JuMPModels/moonlander.jl
+++ b/ext/JuMPModels/moonlander.jl
@@ -29,7 +29,7 @@ julia> model = OptimalControlProblems.moonlander(JuMPBackend(); N=100)
 - Problem formulation available at: https://arxiv.org/pdf/2303.16746
 """
 function OptimalControlProblems.moonlander(
-    ::JuMPBackend; N::Int=steps_number_data(:moonlander)
+    ::JuMPBackend, args...; N::Int=steps_number_data(:moonlander), kwargs...
 )
 
     # parameters
@@ -41,7 +41,7 @@ function OptimalControlProblems.moonlander(
     max_thrust = 2g
 
     # define the problem
-    model = JuMP.Model()
+    model = JuMP.Model(args...; kwargs...)
 
     # state, control and final time variables
     @variables(
diff --git a/ext/JuMPModels/quadrotor.jl b/ext/JuMPModels/quadrotor.jl
index c6970041..4740e835 100644
--- a/ext/JuMPModels/quadrotor.jl
+++ b/ext/JuMPModels/quadrotor.jl
@@ -28,7 +28,7 @@ julia> model = OptimalControlProblems.quadrotor(JuMPBackend(); N=20)
 - Problem formulation available at: https://arxiv.org/pdf/2303.16746
 """
 function OptimalControlProblems.quadrotor(
-    ::JuMPBackend; N::Int=steps_number_data(:quadrotor)
+    ::JuMPBackend, args...; N::Int=steps_number_data(:quadrotor), kwargs...
 )
 
     # parameters
@@ -44,7 +44,7 @@ function OptimalControlProblems.quadrotor(
     vf = [0, 0, 0]
 
     # model
-    model = JuMP.Model()
+    model = JuMP.Model(args...; kwargs...)
 
     # 
     @variables(
diff --git a/ext/JuMPModels/robbins.jl b/ext/JuMPModels/robbins.jl
index b7cca767..37eacdac 100644
--- a/ext/JuMPModels/robbins.jl
+++ b/ext/JuMPModels/robbins.jl
@@ -27,7 +27,7 @@ julia> model = OptimalControlProblems.robbins(JuMPBackend(); N=100)
 
 - Problem formulation available at: https://github.com/control-toolbox/bocop/tree/main/bocop
 """
-function OptimalControlProblems.robbins(::JuMPBackend; N::Int=steps_number_data(:robbins))
+function OptimalControlProblems.robbins(::JuMPBackend, args...; N::Int=steps_number_data(:robbins), kwargs...)
 
     # parameters
     tf = final_time_data(:robbins)
@@ -39,7 +39,7 @@ function OptimalControlProblems.robbins(::JuMPBackend; N::Int=steps_number_data(
     step = tf / N
 
     # model
-    model = JuMP.Model()
+    model = JuMP.Model(args...; kwargs...)
 
     # state, control and initial guess
     @variables(
diff --git a/ext/JuMPModels/robot.jl b/ext/JuMPModels/robot.jl
index bc9afb51..4cf359c7 100644
--- a/ext/JuMPModels/robot.jl
+++ b/ext/JuMPModels/robot.jl
@@ -27,7 +27,7 @@ julia> model = OptimalControlProblems.robot(JuMPBackend(); N=100)
 
 - Problem formulation available at: https://github.com/MadNLP/COPSBenchmark.jl/blob/main/src/robot.jl
 """
-function OptimalControlProblems.robot(::JuMPBackend; N::Int=steps_number_data(:robot))
+function OptimalControlProblems.robot(::JuMPBackend, args...; N::Int=steps_number_data(:robot), kwargs...)
 
     # parameters
 
@@ -45,7 +45,7 @@ function OptimalControlProblems.robot(::JuMPBackend; N::Int=steps_number_data(:r
     θf = 2π/3
 
     # model
-    model = JuMP.Model()
+    model = JuMP.Model(args...; kwargs...)
 
     # state, control, variable (final time) and initial guess
     @variables(
diff --git a/ext/JuMPModels/rocket.jl b/ext/JuMPModels/rocket.jl
index d7aeaa3e..8798f4cd 100644
--- a/ext/JuMPModels/rocket.jl
+++ b/ext/JuMPModels/rocket.jl
@@ -27,7 +27,7 @@ julia> model = OptimalControlProblems.rocket(JuMPBackend(); N=200)
 
 - Problem formulation available at: https://github.com/MadNLP/COPSBenchmark.jl/blob/main/src/rocket.jl
 """
-function OptimalControlProblems.rocket(::JuMPBackend; N::Int=steps_number_data(:rocket))
+function OptimalControlProblems.rocket(::JuMPBackend, args...; N::Int=steps_number_data(:rocket), kwargs...)
 
     # parameters
     h0 = 1
@@ -44,7 +44,7 @@ function OptimalControlProblems.rocket(::JuMPBackend; N::Int=steps_number_data(:
     Tmax = Tc * m0 * g0
 
     # model
-    model = JuMP.Model()
+    model = JuMP.Model(args...; kwargs...)
 
     # state, control, variable (final time) and initial guess
     @variables(
diff --git a/ext/JuMPModels/space_shuttle.jl b/ext/JuMPModels/space_shuttle.jl
index bbf34c48..e5ead74c 100644
--- a/ext/JuMPModels/space_shuttle.jl
+++ b/ext/JuMPModels/space_shuttle.jl
@@ -29,7 +29,7 @@ julia> model = OptimalControlProblems.space_shuttle(JuMPBackend(); N=200)
 - Problem formulation and tutorial available at: https://jump.dev/JuMP.jl/stable/tutorials/nonlinear/space_shuttle_reentry_trajectory/
 """
 function OptimalControlProblems.space_shuttle(
-    ::JuMPBackend; N::Int=steps_number_data(:space_shuttle)
+    ::JuMPBackend, args...; N::Int=steps_number_data(:space_shuttle), kwargs...
 )
 
     ## Global variables
@@ -76,7 +76,7 @@ function OptimalControlProblems.space_shuttle(
     γ_t = deg2rad(-5)  # flight path angle (rad)
 
     # model
-    model = JuMP.Model()
+    model = JuMP.Model(args...; kwargs...)
 
     # state, control and variable (final time)
     @variables(
diff --git a/ext/JuMPModels/steering.jl b/ext/JuMPModels/steering.jl
index 30caa880..ad53a363 100644
--- a/ext/JuMPModels/steering.jl
+++ b/ext/JuMPModels/steering.jl
@@ -27,7 +27,7 @@ julia> model = OptimalControlProblems.steering(JuMPBackend(); N=200)
 
 - Problem formulation available at: https://github.com/MadNLP/COPSBenchmark.jl/blob/main/src/steering.jl
 """
-function OptimalControlProblems.steering(::JuMPBackend; N::Int=steps_number_data(:steering))
+function OptimalControlProblems.steering(::JuMPBackend, args...; N::Int=steps_number_data(:steering), kwargs...)
 
     # parameters
     a = 100
@@ -49,7 +49,7 @@ function OptimalControlProblems.steering(::JuMPBackend; N::Int=steps_number_data
     end
 
     # model
-    model = JuMP.Model()
+    model = JuMP.Model(args...; kwargs...)
 
     @variable(model, u_min <= u[i = 1:(N + 1)] <= u_max, start = 0)   # control
     @variable(model, x1[i = 1:(N + 1)], start = gen_x0(i, 1))           # state x1
diff --git a/ext/JuMPModels/truck_trailer.jl b/ext/JuMPModels/truck_trailer.jl
index 0bee1ff6..a9bce34f 100644
--- a/ext/JuMPModels/truck_trailer.jl
+++ b/ext/JuMPModels/truck_trailer.jl
@@ -28,7 +28,7 @@ julia> model = OptimalControlProblems.truck_trailer(JuMPBackend(); N=100)
 - Problem formulation available at: https://arxiv.org/pdf/2303.16746
 """
 function OptimalControlProblems.truck_trailer(
-    ::JuMPBackend; N::Int=steps_number_data(:truck_trailer)
+    ::JuMPBackend, args...; N::Int=steps_number_data(:truck_trailer), kwargs...
 )
 
     # parameters
@@ -55,7 +55,7 @@ function OptimalControlProblems.truck_trailer(
     θ0_tf = π / 2
 
     # model
-    model = JuMP.Model()
+    model = JuMP.Model(args...; kwargs...)
 
     # state, control, variable (final time) and initial guess
     @variables(
diff --git a/ext/JuMPModels/vanderpol.jl b/ext/JuMPModels/vanderpol.jl
index 44857fdd..8e4eeda6 100644
--- a/ext/JuMPModels/vanderpol.jl
+++ b/ext/JuMPModels/vanderpol.jl
@@ -27,7 +27,7 @@ julia> model = OptimalControlProblems.vanderpol(JuMPBackend(); N=100)
 - Problem formulation available at: https://github.com/control-toolbox/bocop/tree/main/bocop
 """
 function OptimalControlProblems.vanderpol(
-    ::JuMPBackend; N::Int=steps_number_data(:vanderpol)
+    ::JuMPBackend, args...; N::Int=steps_number_data(:vanderpol), kwargs...
 )
 
     # parameters
@@ -36,7 +36,7 @@ function OptimalControlProblems.vanderpol(
     ε = 1
 
     # model
-    model = JuMP.Model()
+    model = JuMP.Model(args...; kwargs...)
 
     # state, control and initial guess
     @variables(
diff --git a/ext/OptimalControlModels/beam.jl b/ext/OptimalControlModels/beam.jl
index 64d39c79..ed1760bd 100644
--- a/ext/OptimalControlModels/beam.jl
+++ b/ext/OptimalControlModels/beam.jl
@@ -28,7 +28,7 @@ julia> docp = OptimalControlProblems.beam(OptimalControlBackend(); N=100);
 - BOCOP repository: https://github.com/control-toolbox/bocop/tree/main/bocop
 """
 function OptimalControlProblems.beam(
-    ::OptimalControlBackend; N::Int=steps_number_data(:beam)
+    ::OptimalControlBackend, description::Symbol...; N::Int=steps_number_data(:beam), kwargs...
 )
 
     #
@@ -51,7 +51,7 @@ function OptimalControlProblems.beam(
     init = (state=[0.05, 0.1], control=0.1)
 
     # DOCP and NLP
-    docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
+    docp = direct_transcription(ocp, description...; init=init, grid_size=N, disc_method=:trapeze, kwargs...)
 
     return docp
 end
diff --git a/ext/OptimalControlModels/bioreactor.jl b/ext/OptimalControlModels/bioreactor.jl
index 8ff74fb0..36e4d502 100644
--- a/ext/OptimalControlModels/bioreactor.jl
+++ b/ext/OptimalControlModels/bioreactor.jl
@@ -28,7 +28,7 @@ julia> docp = OptimalControlProblems.bioreactor(OptimalControlBackend(); N=100);
 - BOCOP repository: https://github.com/control-toolbox/bocop/tree/main/bocop
 """
 function OptimalControlProblems.bioreactor(
-    ::OptimalControlBackend; N::Int=steps_number_data(:bioreactor)
+    ::OptimalControlBackend, description::Symbol...; N::Int=steps_number_data(:bioreactor), kwargs...
 )
 
     # METHANE PROBLEM
@@ -89,7 +89,7 @@ function OptimalControlProblems.bioreactor(
     init = (state=[50, 50, 50], control=0.5)
 
     # NLPModel + DOCP
-    docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
+    docp = direct_transcription(ocp, description...; init=init, grid_size=N, disc_method=:trapeze, kwargs...)
 
     return docp
 end
diff --git a/ext/OptimalControlModels/cart_pendulum.jl b/ext/OptimalControlModels/cart_pendulum.jl
index 6cba2f59..4c5e6f09 100644
--- a/ext/OptimalControlModels/cart_pendulum.jl
+++ b/ext/OptimalControlModels/cart_pendulum.jl
@@ -28,7 +28,7 @@ julia> docp = OptimalControlProblems.cart_pendulum(OptimalControlBackend(); N=10
 - Formulation inspired by OptimalControl approach for swing-up control problems.
 """
 function OptimalControlProblems.cart_pendulum(
-    ::OptimalControlBackend; N::Int=steps_number_data(:cart_pendulum)
+    ::OptimalControlBackend, description::Symbol...; N::Int=steps_number_data(:cart_pendulum), kwargs...
 )
 
     # parameters
@@ -102,7 +102,7 @@ function OptimalControlProblems.cart_pendulum(
     init = (state=xinit, control=uinit, variable=varinit)
 
     # NLPModel + DOCP
-    docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
+    docp = direct_transcription(ocp, description...; init=init, grid_size=N, disc_method=:trapeze, kwargs...)
 
     return docp
 end
diff --git a/ext/OptimalControlModels/chain.jl b/ext/OptimalControlModels/chain.jl
index 4d320693..53ef6a55 100644
--- a/ext/OptimalControlModels/chain.jl
+++ b/ext/OptimalControlModels/chain.jl
@@ -29,7 +29,7 @@ julia> docp = OptimalControlProblems.chain(OptimalControlBackend(); N=100);
 - Original problem source: [BOCOP repository](https://github.com/control-toolbox/bocop/tree/main/bocop)
 """
 function OptimalControlProblems.chain(
-    ::OptimalControlBackend; N::Int=steps_number_data(:chain)
+    ::OptimalControlBackend, description::Symbol...; N::Int=steps_number_data(:chain), kwargs...
 )
 
     # parameters
@@ -80,7 +80,7 @@ function OptimalControlProblems.chain(
     init = (state=xinit, control=uinit)
 
     # NLPModel + DOCP
-    docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
+    docp = direct_transcription(ocp, description...; init=init, grid_size=N, disc_method=:trapeze, kwargs...)
 
     return docp
 end
diff --git a/ext/OptimalControlModels/dielectrophoretic_particle.jl b/ext/OptimalControlModels/dielectrophoretic_particle.jl
index 1c1c2460..4b313d5c 100644
--- a/ext/OptimalControlModels/dielectrophoretic_particle.jl
+++ b/ext/OptimalControlModels/dielectrophoretic_particle.jl
@@ -29,7 +29,7 @@ julia> docp = OptimalControlProblems.dielectrophoretic_particle(OptimalControlBa
 - Formulation inspired by OptimalControl approach to time-optimal trajectory problems.
 """
 function OptimalControlProblems.dielectrophoretic_particle(
-    ::OptimalControlBackend; N::Int=steps_number_data(:dielectrophoretic_particle)
+    ::OptimalControlBackend, description::Symbol...; N::Int=steps_number_data(:dielectrophoretic_particle), kwargs...
 )
 
     # parameters
@@ -64,7 +64,7 @@ function OptimalControlProblems.dielectrophoretic_particle(
     init = (state=[1, 1], control=0.1, variable=5)
 
     # DOCP and NLP
-    docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
+    docp = direct_transcription(ocp, description...; init=init, grid_size=N, disc_method=:trapeze, kwargs...)
 
     return docp
 end
diff --git a/ext/OptimalControlModels/double_oscillator.jl b/ext/OptimalControlModels/double_oscillator.jl
index 8c9dde35..45e7edcf 100644
--- a/ext/OptimalControlModels/double_oscillator.jl
+++ b/ext/OptimalControlModels/double_oscillator.jl
@@ -29,7 +29,7 @@ julia> docp = OptimalControlProblems.double_oscillator(OptimalControlBackend();
 - Formulation follows OptimalControl approach to mechanical oscillator trajectory optimisation.
 """
 function OptimalControlProblems.double_oscillator(
-    ::OptimalControlBackend; N::Int=steps_number_data(:double_oscillator)
+    ::OptimalControlBackend, description::Symbol...; N::Int=steps_number_data(:double_oscillator), kwargs...
 )
 
     # parameters
@@ -75,7 +75,7 @@ function OptimalControlProblems.double_oscillator(
     init = (state=xinit, control=uinit)
 
     # DOCP and NLP
-    docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
+    docp = direct_transcription(ocp, description...; init=init, grid_size=N, disc_method=:trapeze, kwargs...)
 
     return docp
 end
diff --git a/ext/OptimalControlModels/ducted_fan.jl b/ext/OptimalControlModels/ducted_fan.jl
index 60ba5f98..69908c0c 100644
--- a/ext/OptimalControlModels/ducted_fan.jl
+++ b/ext/OptimalControlModels/ducted_fan.jl
@@ -29,7 +29,7 @@ julia> docp = OptimalControlProblems.ducted_fan(OptimalControlBackend(); N=250);
 - Problem instance follows OptimalControl formulation for ducted fan trajectory optimisation.
 """
 function OptimalControlProblems.ducted_fan(
-    ::OptimalControlBackend; N::Int=steps_number_data(:ducted_fan)
+    ::OptimalControlBackend, description::Symbol...; N::Int=steps_number_data(:ducted_fan), kwargs...
 )
 
     # parameters
@@ -96,7 +96,7 @@ function OptimalControlProblems.ducted_fan(
     init = (state=xinit, control=uinit, variable=varinit)
 
     # DOCP and NLP
-    docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
+    docp = direct_transcription(ocp, description...; init=init, grid_size=N, disc_method=:trapeze, kwargs...)
 
     return docp
 end
diff --git a/ext/OptimalControlModels/electric_vehicle.jl b/ext/OptimalControlModels/electric_vehicle.jl
index 2fd3bc7f..03c8885c 100644
--- a/ext/OptimalControlModels/electric_vehicle.jl
+++ b/ext/OptimalControlModels/electric_vehicle.jl
@@ -29,7 +29,7 @@ julia> docp = OptimalControlProblems.electric_vehicle(OptimalControlBackend(); N
 - Problem instance follows OptimalControl formulation for electric vehicle trajectory optimisation.
 """
 function OptimalControlProblems.electric_vehicle(
-    ::OptimalControlBackend; N::Int=steps_number_data(:electric_vehicle)
+    ::OptimalControlBackend, description::Symbol...; N::Int=steps_number_data(:electric_vehicle), kwargs...
 )
 
     # parameters
@@ -67,7 +67,7 @@ function OptimalControlProblems.electric_vehicle(
     init = (state=yinit, control=uinit)
 
     # DOCP and NLP
-    docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
+    docp = direct_transcription(ocp, description...; init=init, grid_size=N, disc_method=:trapeze, kwargs...)
 
     return docp
 end
diff --git a/ext/OptimalControlModels/glider.jl b/ext/OptimalControlModels/glider.jl
index 3dac7d2b..bf28bc08 100644
--- a/ext/OptimalControlModels/glider.jl
+++ b/ext/OptimalControlModels/glider.jl
@@ -29,7 +29,7 @@ julia> docp = OptimalControlProblems.glider(OptimalControlBackend(); N=500);
 - Problem inspired by glider dynamics with thermal updraft and lift modelling.
 """
 function OptimalControlProblems.glider(
-    ::OptimalControlBackend; N::Int=steps_number_data(:glider)
+    ::OptimalControlBackend, description::Symbol...; N::Int=steps_number_data(:glider), kwargs...
 )
 
     # parameters
@@ -107,7 +107,7 @@ function OptimalControlProblems.glider(
     init = (state=xinit, control=uinit, variable=tfinit)
 
     # DOCP and NLP
-    docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
+    docp = direct_transcription(ocp, description...; init=init, grid_size=N, disc_method=:trapeze, kwargs...)
 
     return docp
 end
diff --git a/ext/OptimalControlModels/insurance.jl b/ext/OptimalControlModels/insurance.jl
index 8182dfea..d9cd0fab 100644
--- a/ext/OptimalControlModels/insurance.jl
+++ b/ext/OptimalControlModels/insurance.jl
@@ -28,7 +28,7 @@ julia> docp = OptimalControlProblems.insurance(OptimalControlBackend(); N=500);
 - Problem formulation available at [Bocop repository](https://github.com/control-toolbox/bocop/tree/main/bocop)
 """
 function OptimalControlProblems.insurance(
-    ::OptimalControlBackend; N::Int=steps_number_data(:insurance)
+    ::OptimalControlBackend, description::Symbol...; N::Int=steps_number_data(:insurance), kwargs...
 )
 
     # parameters
@@ -100,7 +100,7 @@ function OptimalControlProblems.insurance(
     init = (state=xinit, control=uinit, variable=varinit)
 
     # DOCP and NLP
-    docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
+    docp = direct_transcription(ocp, description...; init=init, grid_size=N, disc_method=:trapeze, kwargs...)
 
     return docp
 end
diff --git a/ext/OptimalControlModels/jackson.jl b/ext/OptimalControlModels/jackson.jl
index d7957a9b..7f6a38f0 100644
--- a/ext/OptimalControlModels/jackson.jl
+++ b/ext/OptimalControlModels/jackson.jl
@@ -28,7 +28,7 @@ julia> docp = OptimalControlProblems.jackson(OptimalControlBackend(); N=500);
 - Problem formulation available at [Bocop repository](https://github.com/control-toolbox/bocop/tree/main/bocop)
 """
 function OptimalControlProblems.jackson(
-    ::OptimalControlBackend; N::Int=steps_number_data(:jackson)
+    ::OptimalControlBackend, description::Symbol...; N::Int=steps_number_data(:jackson), kwargs...
 )
 
     # parameters
@@ -66,7 +66,7 @@ function OptimalControlProblems.jackson(
     init = (state=xinit, control=uinit)
 
     # DOCP and NLP
-    docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
+    docp = direct_transcription(ocp, description...; init=init, grid_size=N, disc_method=:trapeze, kwargs...)
 
     return docp
 end
diff --git a/ext/OptimalControlModels/moonlander.jl b/ext/OptimalControlModels/moonlander.jl
index c7b9852b..3e1c5515 100644
--- a/ext/OptimalControlModels/moonlander.jl
+++ b/ext/OptimalControlModels/moonlander.jl
@@ -24,7 +24,7 @@ julia> docp = OptimalControlProblems.moonlander(OptimalControlBackend(); N=500);
 ```
 """
 function OptimalControlProblems.moonlander(
-    ::OptimalControlBackend; N::Int=steps_number_data(:moonlander)
+    ::OptimalControlBackend, description::Symbol...; N::Int=steps_number_data(:moonlander), kwargs...
 )
 
     # parameters
@@ -97,7 +97,7 @@ function OptimalControlProblems.moonlander(
     init = (state=xinit, control=uinit, variable=varinit)
 
     # DOCP and NLP
-    docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
+    docp = direct_transcription(ocp, description...; init=init, grid_size=N, disc_method=:trapeze, kwargs...)
 
     return docp
 end
diff --git a/ext/OptimalControlModels/quadrotor.jl b/ext/OptimalControlModels/quadrotor.jl
index 804ecd77..74fe21ff 100644
--- a/ext/OptimalControlModels/quadrotor.jl
+++ b/ext/OptimalControlModels/quadrotor.jl
@@ -24,7 +24,7 @@ julia> docp = OptimalControlProblems.quadrotor(OptimalControlBackend(); N=50);
 ```
 """
 function OptimalControlProblems.quadrotor(
-    ::OptimalControlBackend; N::Int=steps_number_data(:quadrotor)
+    ::OptimalControlBackend, description::Symbol...; N::Int=steps_number_data(:quadrotor), kwargs...
 )
 
     # parameters
@@ -112,7 +112,7 @@ function OptimalControlProblems.quadrotor(
     init = (state=xinit, control=uinit, variable=varinit)
 
     # DOCP and NLP
-    docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
+    docp = direct_transcription(ocp, description...; init=init, grid_size=N, disc_method=:trapeze, kwargs...)
 
     return docp
 end
diff --git a/ext/OptimalControlModels/robbins.jl b/ext/OptimalControlModels/robbins.jl
index 71ae0b0d..3b35c1f3 100644
--- a/ext/OptimalControlModels/robbins.jl
+++ b/ext/OptimalControlModels/robbins.jl
@@ -25,7 +25,7 @@ julia> docp = OptimalControlProblems.robbins(OptimalControlBackend(); N=500);
 ```
 """
 function OptimalControlProblems.robbins(
-    ::OptimalControlBackend; N::Int=steps_number_data(:robbins)
+    ::OptimalControlBackend, description::Symbol...; N::Int=steps_number_data(:robbins), kwargs...
 )
 
     # parameters
@@ -56,7 +56,7 @@ function OptimalControlProblems.robbins(
     init = (state=xinit, control=uinit)
 
     # DOCP and NLP
-    docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
+    docp = direct_transcription(ocp, description...; init=init, grid_size=N, disc_method=:trapeze, kwargs...)
 
     return docp
 end
diff --git a/ext/OptimalControlModels/robot.jl b/ext/OptimalControlModels/robot.jl
index 6044a0da..0fc69a69 100644
--- a/ext/OptimalControlModels/robot.jl
+++ b/ext/OptimalControlModels/robot.jl
@@ -24,7 +24,7 @@ julia> docp = OptimalControlProblems.robot(OptimalControlBackend(); N=250);
 ```
 """
 function OptimalControlProblems.robot(
-    ::OptimalControlBackend; N::Int=steps_number_data(:robot)
+    ::OptimalControlBackend, description::Symbol...; N::Int=steps_number_data(:robot), kwargs...
 )
 
     # parameters
@@ -94,7 +94,7 @@ function OptimalControlProblems.robot(
     init = (state=xinit, control=uinit, variable=tf)
 
     # DOCP and NLP
-    docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
+    docp = direct_transcription(ocp, description...; init=init, grid_size=N, disc_method=:trapeze, kwargs...)
 
     return docp
 end
diff --git a/ext/OptimalControlModels/rocket.jl b/ext/OptimalControlModels/rocket.jl
index f4d9dea7..3d11ba28 100644
--- a/ext/OptimalControlModels/rocket.jl
+++ b/ext/OptimalControlModels/rocket.jl
@@ -24,7 +24,7 @@ julia> docp = OptimalControlProblems.rocket(OptimalControlBackend(); N=500);
 ```
 """
 function OptimalControlProblems.rocket(
-    ::OptimalControlBackend; N::Int=steps_number_data(:rocket)
+    ::OptimalControlBackend, description::Symbol...; N::Int=steps_number_data(:rocket), kwargs...
 )
 
     # parameters
@@ -87,7 +87,7 @@ function OptimalControlProblems.rocket(
     init = (time=time_vec, state=xinit, control=Tmax/2, variable=1)
 
     # DOCP and NLP
-    docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
+    docp = direct_transcription(ocp, description...; init=init, grid_size=N, disc_method=:trapeze, kwargs...)
 
     return docp
 end
diff --git a/ext/OptimalControlModels/space_shuttle.jl b/ext/OptimalControlModels/space_shuttle.jl
index 46361a03..44d62d41 100644
--- a/ext/OptimalControlModels/space_shuttle.jl
+++ b/ext/OptimalControlModels/space_shuttle.jl
@@ -25,7 +25,7 @@ julia> docp = OptimalControlProblems.space_shuttle(OptimalControlBackend(); N=50
 ```
 """
 function OptimalControlProblems.space_shuttle(
-    ::OptimalControlBackend; N::Int=steps_number_data(:space_shuttle)
+    ::OptimalControlBackend, description::Symbol...; N::Int=steps_number_data(:space_shuttle), kwargs...
 )
 
     ## Global variables
@@ -158,7 +158,7 @@ function OptimalControlProblems.space_shuttle(
     init = (state=x_init, control=[α_s, β_s], variable=[tf_init])
 
     # DOCP and NLP
-    docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
+    docp = direct_transcription(ocp, description...; init=init, grid_size=N, disc_method=:trapeze, kwargs...)
 
     return docp
 end
diff --git a/ext/OptimalControlModels/steering.jl b/ext/OptimalControlModels/steering.jl
index 26254bb5..3a8d885f 100644
--- a/ext/OptimalControlModels/steering.jl
+++ b/ext/OptimalControlModels/steering.jl
@@ -24,7 +24,7 @@ julia> docp = OptimalControlProblems.steering(OptimalControlBackend(); N=500);
 ```
 """
 function OptimalControlProblems.steering(
-    ::OptimalControlBackend; N::Int=steps_number_data(:steering)
+    ::OptimalControlBackend, description::Symbol...; N::Int=steps_number_data(:steering), kwargs...
 )
 
     # parameters
@@ -70,7 +70,7 @@ function OptimalControlProblems.steering(
     init = (state=xinit, control=0, variable=1)
 
     # DOCP and NLP
-    docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
+    docp = direct_transcription(ocp, description...; init=init, grid_size=N, disc_method=:trapeze, kwargs...)
 
     return docp
 end
diff --git a/ext/OptimalControlModels/truck_trailer.jl b/ext/OptimalControlModels/truck_trailer.jl
index a8f8bca8..c33a8a79 100644
--- a/ext/OptimalControlModels/truck_trailer.jl
+++ b/ext/OptimalControlModels/truck_trailer.jl
@@ -24,7 +24,7 @@ julia> docp = OptimalControlProblems.truck_trailer(OptimalControlBackend(); N=20
 ```
 """
 function OptimalControlProblems.truck_trailer(
-    ::OptimalControlBackend; N::Int=steps_number_data(:truck_trailer)
+    ::OptimalControlBackend, description::Symbol...; N::Int=steps_number_data(:truck_trailer), kwargs...
 )
 
     # parameters
@@ -130,7 +130,7 @@ function OptimalControlProblems.truck_trailer(
     init = (state=xinit, control=uinit, variable=varinit)
 
     # DOCP and NLP
-    docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
+    docp = direct_transcription(ocp, description...; init=init, grid_size=N, disc_method=:trapeze, kwargs...)
 
     return docp
 end
diff --git a/ext/OptimalControlModels/vanderpol.jl b/ext/OptimalControlModels/vanderpol.jl
index 8c4d791a..64fc79a4 100644
--- a/ext/OptimalControlModels/vanderpol.jl
+++ b/ext/OptimalControlModels/vanderpol.jl
@@ -24,7 +24,7 @@ julia> docp = OptimalControlProblems.vanderpol(OptimalControlBackend(); N=500);
 ```
 """
 function OptimalControlProblems.vanderpol(
-    ::OptimalControlBackend; N::Int=steps_number_data(:vanderpol)
+    ::OptimalControlBackend, description::Symbol...; N::Int=steps_number_data(:vanderpol), kwargs...
 )
 
     # parameters
@@ -51,7 +51,7 @@ function OptimalControlProblems.vanderpol(
     init = (state=xinit, control=uinit)
 
     # DOCP and NLP
-    docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
+    docp = direct_transcription(ocp, description...; init=init, grid_size=N, disc_method=:trapeze, kwargs...)
 
     return docp
 end
diff --git a/src/OptimalControlProblems.jl b/src/OptimalControlProblems.jl
index 03540fad..d7e9f714 100644
--- a/src/OptimalControlProblems.jl
+++ b/src/OptimalControlProblems.jl
@@ -72,7 +72,7 @@ for file in files
     """
 
     code = quote
-        @doc $doc function $problem(model_backend::T; N::Int=0) where {T<:AbstractModelBackend}
+        @doc $doc function $problem(model_backend::T, args...; kwargs...) where {T<:AbstractModelBackend}
             throw(CTBase.ExtensionError(weakdeps[T]))
         end
         export $problem
@@ -110,7 +110,7 @@ const types = [
 ]
 
 """
-OptimalControlProblems.metadata::Dict()
+metadata::Dict()
 
 Dictionary containing metadata for all available optimal control problems.
 
@@ -130,7 +130,7 @@ The following keys are valid:
 # Example
 
 ```julia-repl
-julia> OptimalControlProblems.metadata[:my_problem][:name]
+julia> metadata[:my_problem][:name]
 "My Problem"
 ```
 """
diff --git a/test/Project.toml b/test/Project.toml
index dcb02017..b71fe07e 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,6 +1,7 @@
 [deps]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 CTBase = "54762871-cc72-4466-b8e8-f6c8b58076cd"
+CTDirect = "790bbbee-bee9-49ee-8912-a9de031322d5"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
diff --git a/test/figures/init/beam.pdf b/test/figures/init/beam.pdf
index c6eefb8a..124652da 100644
Binary files a/test/figures/init/beam.pdf and b/test/figures/init/beam.pdf differ
diff --git a/test/figures/solution/beam.pdf b/test/figures/solution/beam.pdf
index 95dd1f84..45c1bd9d 100644
Binary files a/test/figures/solution/beam.pdf and b/test/figures/solution/beam.pdf differ
diff --git a/test/runtests.jl b/test/runtests.jl
index 15c38240..f08c5c8d 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,5 +1,6 @@
 using Aqua
 using CTBase
+using CTDirect
 using Ipopt
 using JuMP
 using NLPModelsIpopt
@@ -49,6 +50,7 @@ const VERBOSE = true # print or not details during tests
 @testset "OptimalControlProblems tests" verbose=VERBOSE showtiming=true begin
     for name in (
         #:aqua,
+        :kwargs,
         :JuMP,                  # convergence tests for JuMP models
         :OptimalControl,        # convergence tests for OptimalControl models
         :init,                  # comparison between OptimalControl and JuMP: init
diff --git a/test/test_JuMP.jl b/test/test_JuMP.jl
index be8a29f6..c53c8dbe 100644
--- a/test/test_JuMP.jl
+++ b/test/test_JuMP.jl
@@ -2,7 +2,7 @@
 function test_JuMP()
     for f in LIST_OF_PROBLEMS
         @testset "$(f)" verbose=VERBOSE begin
-            N = OptimalControlProblems.metadata[f][:N]
+            N = metadata[f][:N]
 
             # do we keep or remove the problem from the list
             keep_problem = true
diff --git a/test/test_OptimalControl.jl b/test/test_OptimalControl.jl
index a77bd684..b7532b57 100644
--- a/test/test_OptimalControl.jl
+++ b/test/test_OptimalControl.jl
@@ -11,7 +11,7 @@ function test_OptimalControl()
 
     for f in LIST_OF_PROBLEMS
         @testset "$(f)" verbose=VERBOSE begin
-            N = OptimalControlProblems.metadata[f][:N]
+            N = metadata[f][:N]
 
             # do we keep or remove the problem from the list
             keep_problem = true
diff --git a/test/test_kwargs.jl b/test/test_kwargs.jl
new file mode 100644
index 00000000..00dbafdc
--- /dev/null
+++ b/test/test_kwargs.jl
@@ -0,0 +1,28 @@
+function test_kwargs()
+
+    N = 2
+    scheme = :euler
+    solver_backend = :madnlp
+    optimiser = Ipopt.Optimizer
+
+    for f in LIST_OF_PROBLEMS
+
+        @testset "$(string(f)) (kwargs)" verbose=VERBOSE begin
+
+            DEBUG && println("─ ", string(f), " (kwargs)")
+
+            # OptimalControl model
+            docp = OptimalControlProblems.eval(f)(OptimalControlBackend(), solver_backend; grid_size=N, disc_method=scheme)
+            @test docp isa CTDirect.DOCP
+            @test docp.time.steps == N
+            @test docp.discretization isa CTDirect.Euler
+
+            # JuMP model
+            nlp = OptimalControlProblems.eval(f)(JuMPBackend(), optimiser; add_bridges=true)
+            @test solver_name(nlp) == "Ipopt"
+
+        end
+
+    end
+
+end
\ No newline at end of file
diff --git a/test/test_quick.jl b/test/test_quick.jl
index 5741e429..e70f066a 100644
--- a/test/test_quick.jl
+++ b/test/test_quick.jl
@@ -19,7 +19,7 @@ function test_quick()
     max_r_err = -Inf # relative error max
 
     for f in LIST_OF_PROBLEMS
-        N = OptimalControlProblems.metadata[f][:N]
+        N = metadata[f][:N]
 
         @testset "$(string(f)) (objective)" verbose=VERBOSE begin
             DEBUG && println("\n", "┌─ ", string(f))
diff --git a/test/utils.jl b/test/utils.jl
index 9a34000d..f46bcc68 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -108,11 +108,11 @@ function comparison(; max_iter, test_name)
 
     # we loop over the problems
     for f in LIST_OF_PROBLEMS
-        N = OptimalControlProblems.metadata[f][:N] # get default N
-        x_vars = OptimalControlProblems.metadata[f][:state_name]
-        p_vars = OptimalControlProblems.metadata[f][:costate_name]
-        u_vars = OptimalControlProblems.metadata[f][:control_name]
-        v_vars = OptimalControlProblems.metadata[f][:variable_name]
+        N = metadata[f][:N] # get default N
+        x_vars = metadata[f][:state_name]
+        p_vars = metadata[f][:costate_name]
+        u_vars = metadata[f][:control_name]
+        v_vars = metadata[f][:variable_name]
 
         @testset "$(string(f)) ($(string(test_name)))" verbose=VERBOSE begin
             DEBUG && println("\n", "┌─ ", string(f), " (", string(test_name), ")")