Add the choice for an integrator of the solution

README.md

@@ -88,6 +88,7 @@ As a tour guide that uses this binder, you can watch the `bioptim` workshop that
   - [InterpolationType](#enum-interpolationtype)
   - [Shooting](#enum-shooting)
   - [CostType](#enum-costtype)
+  - [SolutionIntegrator](#enum-solutionintegrator)
 
 [Examples](#examples)
 - [Run examples](#run-examples)
@@ -1383,7 +1384,6 @@ The accepted values are:
 - SPLINE: Requires five columns. It performs a cubic spline to interpolate between the nodes.
 - CUSTOM: User defined interpolation function
 
-
 ### Enum: Shooting
 The type of integration to perform
 - MULTIPLE: resets the state at each node
@@ -1396,6 +1396,16 @@ The type of cost
 - CONSTRAINTS: The constraints
 - ALL: All the previously described cost type
 
+### Enum: SolutionIntegrator
+The type of integrator used to integrate the solution of the optimal control problem
+- DEFAULT: The default integrator initially chosen with [OdeSolver](#class-odesolver)
+- SCIPY_RK23: The scipy integrator RK23
+- SCIPY_RK45: The scipy integrator RK45
+- SCIPY_DOP853: The scipy integrator DOP853
+- SCIPY_BDF: The scipy integrator BDF
+- SCIPY_LSODA: The scipy integrator LSODA
+
+
 # Examples
 In this section, you will find the description of all the examples implemented with bioptim. They are ordered in 
 separate files. Each subsection corresponds to the different files, dealing with different examples and topics.

bioptim/__init__.py

@@ -165,7 +165,17 @@
 from .limits.path_conditions import BoundsList, Bounds, InitialGuessList, InitialGuess, QAndQDotBounds
 from .limits.fatigue_path_conditions import FatigueBounds, FatigueInitialGuess
 from .limits.penalty_node import PenaltyNode, PenaltyNodeList
-from .misc.enums import Axis, Node, InterpolationType, PlotType, ControlType, CostType, Shooting, VariableType
+from .misc.enums import (
+    Axis,
+    Node,
+    InterpolationType,
+    PlotType,
+    ControlType,
+    CostType,
+    Shooting,
+    VariableType,
+    SolutionIntegrator,
+)
 from .misc.mapping import BiMappingList, BiMapping, Mapping
 from .optimization.non_linear_program import NonLinearProgram
 from .optimization.optimal_control_program import OptimalControlProgram

bioptim/gui/plot.py

@@ -10,7 +10,7 @@
 from casadi import Callback, nlpsol_out, nlpsol_n_out, Sparsity, DM
 
 from ..limits.path_conditions import Bounds
-from ..misc.enums import PlotType, ControlType, InterpolationType, Shooting
+from ..misc.enums import PlotType, ControlType, InterpolationType, Shooting, SolutionIntegrator
 from ..misc.mapping import Mapping
 from ..optimization.solution import Solution
 
@@ -197,7 +197,7 @@ def __init__(
         automatically_organize: bool = True,
         show_bounds: bool = False,
         shooting_type: Shooting = Shooting.MULTIPLE,
-        use_scipy_integrator: bool = False,
+        integrator: SolutionIntegrator = SolutionIntegrator.DEFAULT,
     ):
         """
         Prepares the figures during the simulation
@@ -212,8 +212,9 @@ def __init__(
             If the axes should fit the bounds (True) or the data (False)
         shooting_type: Shooting
             The type of integration method
-        use_scipy_integrator: bool
-            Use the scipy solve_ivp integrator for RungeKutta 45 instead of currently defined integrator
+        integrator: SolutionIntegrator
+             Use the ode defined by OCP or use a separate integrator provided by scipy
+
         """
         for i in range(1, ocp.n_phases):
             if len(ocp.nlp[0].states["q"]) != len(ocp.nlp[i].states["q"]):
@@ -235,7 +236,8 @@ def __init__(
 
         self.t = []
         self.t_integrated = []
-        self.use_scipy_integrator = use_scipy_integrator
+        self.integrator = integrator
+
         if isinstance(self.ocp.original_phase_time, (int, float)):
             self.tf = [self.ocp.original_phase_time]
         else:
@@ -294,12 +296,14 @@ def __update_time_vector(self):
         self.t_integrated = []
         last_t = 0
         for phase_idx, nlp in enumerate(self.ocp.nlp):
-            n_int_steps = nlp.ode_solver.steps_scipy if self.use_scipy_integrator else nlp.ode_solver.steps
+            n_int_steps = (
+                nlp.ode_solver.steps_scipy if self.integrator != SolutionIntegrator.DEFAULT else nlp.ode_solver.steps
+            )
             dt_ns = self.tf[phase_idx] / nlp.ns
             time_phase_integrated = []
             last_t_int = copy(last_t)
             for _ in range(nlp.ns):
-                if nlp.ode_solver.is_direct_collocation and not self.use_scipy_integrator:
+                if nlp.ode_solver.is_direct_collocation and self.integrator == SolutionIntegrator.DEFAULT:
                     time_phase_integrated.append(np.array(nlp.dynamics[0].step_time) * dt_ns + last_t_int)
                 else:
                     time_phase_integrated.append(np.linspace(last_t_int, last_t_int + dt_ns, n_int_steps + 1))
@@ -426,7 +430,11 @@ def legend_without_duplicate_labels(ax):
                         )
                     elif plot_type == PlotType.INTEGRATED:
                         plots_integrated = []
-                        n_int_steps = nlp.ode_solver.steps_scipy if self.use_scipy_integrator else nlp.ode_solver.steps
+                        n_int_steps = (
+                            nlp.ode_solver.steps_scipy
+                            if self.integrator != SolutionIntegrator.DEFAULT
+                            else nlp.ode_solver.steps
+                        )
                         zero = np.zeros(n_int_steps + 1)
                         color = self.plot_func[variable][i].color if self.plot_func[variable][i].color else "tab:brown"
                         for cmp in range(nlp.ns):
@@ -585,7 +593,7 @@ def update_data(self, v: dict):
             continuous=False,
             shooting_type=self.shooting_type,
             keep_intermediate_points=True,
-            use_scipy_integrator=self.use_scipy_integrator,
+            integrator=self.integrator,
         ).states
         data_controls = sol.controls
         data_params = sol.parameters
@@ -598,7 +606,11 @@ def update_data(self, v: dict):
             self.__update_xdata()
 
         for i, nlp in enumerate(self.ocp.nlp):
-            step_size = nlp.ode_solver.steps_scipy + 1 if self.use_scipy_integrator else nlp.ode_solver.steps + 1
+            step_size = (
+                nlp.ode_solver.steps_scipy + 1
+                if self.integrator != SolutionIntegrator.DEFAULT
+                else nlp.ode_solver.steps + 1
+            )
             n_elements = nlp.ns * step_size + 1
 
             state = np.ndarray((0, n_elements))

bioptim/misc/enums.py

@@ -103,3 +103,16 @@ class VariableType(Enum):
 
     STATES = "states"
     CONTROLS = "controls"
+
+
+class SolutionIntegrator(Enum):
+    """
+    Selection of integrator to use integrate function
+    """
+
+    DEFAULT = "DEFAULT"
+    SCIPY_RK23 = "RK23"
+    SCIPY_RK45 = "RK45"
+    SCIPY_DOP853 = "DOP853"
+    SCIPY_BDF = "BDF"
+    SCIPY_LSODA = "LSODA"

bioptim/optimization/optimal_control_program.py

@@ -28,7 +28,7 @@
 from ..limits.penalty import PenaltyOption
 from ..limits.objective_functions import ObjectiveFunction
 from ..misc.__version__ import __version__
-from ..misc.enums import ControlType, SolverType, Shooting, PlotType, CostType
+from ..misc.enums import ControlType, SolverType, Shooting, PlotType, CostType, SolutionIntegrator
 from ..misc.mapping import BiMappingList, Mapping
 from ..misc.utils import check_version
 from ..optimization.parameters import ParameterList, Parameter
@@ -749,7 +749,7 @@ def prepare_plots(
         automatically_organize: bool = True,
         show_bounds: bool = False,
         shooting_type: Shooting = Shooting.MULTIPLE,
-        use_scipy_integrator: bool = False,
+        integrator: SolutionIntegrator = SolutionIntegrator.DEFAULT,
     ) -> PlotOcp:
         """
         Create all the plots associated with the OCP
@@ -762,8 +762,8 @@ def prepare_plots(
             If the ylim should fit the bounds
         shooting_type: Shooting
             What type of integration
-        use_scipy_integrator: bool
-            Use the scipy solve_ivp integrator for RungeKutta 45 instead of currently defined integrator
+        integrator: SolutionIntegrator
+            Use the ode defined by OCP or use a separate integrator provided by scipy
 
         Returns
         -------
@@ -775,7 +775,7 @@ def prepare_plots(
             automatically_organize=automatically_organize,
             show_bounds=show_bounds,
             shooting_type=shooting_type,
-            use_scipy_integrator=use_scipy_integrator,
+            integrator=integrator,
         )
 
     def solve(

bioptim/optimization/solution.py

@@ -9,7 +9,7 @@
 from matplotlib import pyplot as plt
 
 from ..limits.path_conditions import InitialGuess, InitialGuessList
-from ..misc.enums import ControlType, CostType, Shooting, InterpolationType, SolverType
+from ..misc.enums import ControlType, CostType, Shooting, InterpolationType, SolverType, SolutionIntegrator
 from ..misc.utils import check_version
 from ..optimization.non_linear_program import NonLinearProgram
 from ..optimization.optimization_variable import OptimizationVariableList, OptimizationVariable
@@ -492,7 +492,7 @@ def integrate(
         keep_intermediate_points: bool = False,
         merge_phases: bool = False,
         continuous: bool = True,
-        use_scipy_integrator: bool = False,
+        integrator: SolutionIntegrator = SolutionIntegrator.DEFAULT,
     ) -> Any:
         """
         Integrate the states
@@ -509,8 +509,8 @@ def integrate(
         continuous: bool
             If the arrival value of a node should be discarded [True] or kept [False]. The value of an integrated
             arrival node and the beginning of the next one are expected to be almost equal when the problem converged
-        use_scipy_integrator: bool
-            Ignore the dynamics defined by OCP and use a separate integrator provided by scipy
+        integrator: SolutionIntegrator
+            Use the ode defined by OCP or use a separate integrator provided by scipy
 
         Returns
         -------
@@ -535,7 +535,7 @@ def integrate(
                 "Shooting.SINGLE_CONTINUOUS and continuous=False cannot be used simultaneously it is a contradiction"
             )
 
-        out = self.__perform_integration(shooting_type, keep_intermediate_points, continuous, use_scipy_integrator)
+        out = self.__perform_integration(shooting_type, keep_intermediate_points, continuous, integrator)
 
         if merge_phases:
             if continuous:
@@ -548,24 +548,25 @@ def integrate(
         return out
 
     def __perform_integration(
-        self, shooting_type: Shooting, keep_intermediate_points: bool, continuous: bool, use_scipy_integrator: bool
+        self, shooting_type: Shooting, keep_intermediate_points: bool, continuous: bool, integrator: SolutionIntegrator
     ):
         n_direct_collocation = sum([nlp.ode_solver.is_direct_collocation for nlp in self.ocp.nlp])
-        if n_direct_collocation > 0 and not use_scipy_integrator:
+
+        if n_direct_collocation > 0 and integrator == SolutionIntegrator.DEFAULT:
             if continuous:
                 raise RuntimeError(
-                    "Integration with direct collocation must be not continuous if not use_scipy_integrator"
+                    "Integration with direct collocation must be not continuous if a scipy integrator is used"
                 )
 
             if shooting_type != Shooting.MULTIPLE:
                 raise RuntimeError(
                     "Integration with direct collocation must using shooting_type=Shooting.MULTIPLE "
-                    "if not use_scipy_integrator"
+                    "if a scipy integrator is not used"
                 )
 
         # Copy the data
         out = self.copy(skip_data=True)
-        out.recomputed_time_steps = use_scipy_integrator
+        out.recomputed_time_steps = integrator != SolutionIntegrator.DEFAULT
         out._states = []
         for _ in range(len(self._states)):
             out._states.append({})
@@ -575,7 +576,7 @@ def __perform_integration(
         for p, nlp in enumerate(self.ocp.nlp):
             param_scaling = nlp.parameters.scaling
             n_states = self._states[p]["all"].shape[0]
-            n_steps = nlp.ode_solver.steps_scipy if use_scipy_integrator else nlp.ode_solver.steps
+            n_steps = nlp.ode_solver.steps_scipy if integrator != SolutionIntegrator.DEFAULT else nlp.ode_solver.steps
             if not continuous:
                 n_steps += 1
             if keep_intermediate_points:
@@ -596,7 +597,11 @@ def __perform_integration(
                         )
                     x0 += np.array(val)[:, 0]
             else:
-                col = slice(0, n_steps) if nlp.ode_solver.is_direct_collocation and not use_scipy_integrator else 0
+                col = (
+                    slice(0, n_steps)
+                    if nlp.ode_solver.is_direct_collocation and integrator == SolutionIntegrator.DEFAULT
+                    else 0
+                )
                 x0 = self._states[p]["all"][:, col]
 
             for s in range(self.ns[p]):
@@ -607,13 +612,17 @@ def __perform_integration(
                 else:
                     raise NotImplementedError(f"ControlType {nlp.control_type} " f"not yet implemented in integrating")
 
-                if use_scipy_integrator:
+                if integrator != SolutionIntegrator.DEFAULT:
                     t_init = sum(out.phase_time[:p]) / nlp.ns
                     t_end = sum(out.phase_time[: (p + 2)]) / nlp.ns
                     n_points = n_steps + 1 if continuous else n_steps
                     t_eval = np.linspace(t_init, t_end, n_points) if keep_intermediate_points else [t_init, t_end]
                     integrated = solve_ivp(
-                        lambda t, x: np.array(nlp.dynamics_func(x, u, params))[:, 0], [t_init, t_end], x0, t_eval=t_eval
+                        lambda t, x: np.array(nlp.dynamics_func(x, u, params))[:, 0],
+                        [t_init, t_end],
+                        x0,
+                        t_eval=t_eval,
+                        method=integrator.value,
                     ).y
 
                     next_state_col = (
@@ -853,7 +862,7 @@ def graphs(
         show_bounds: bool = False,
         show_now: bool = True,
         shooting_type: Shooting = Shooting.MULTIPLE,
-        use_scipy_integrator: bool = False,
+        integrator: SolutionIntegrator = SolutionIntegrator.DEFAULT,
     ):
         """
         Show the graphs of the simulation
@@ -868,14 +877,14 @@ def graphs(
             If the show method should be called. This is blocking
         shooting_type: Shooting
             The type of interpolation
-        use_scipy_integrator: bool
+        integrator: SolutionIntegrator
             Use the scipy solve_ivp integrator for RungeKutta 45 instead of currently defined integrator
         """
 
         if self.is_merged or self.is_interpolated or self.is_integrated:
             raise NotImplementedError("It is not possible to graph a modified Solution yet")
 
-        plot_ocp = self.ocp.prepare_plots(automatically_organize, show_bounds, shooting_type, use_scipy_integrator)
+        plot_ocp = self.ocp.prepare_plots(automatically_organize, show_bounds, shooting_type, integrator)
         plot_ocp.update_data(self.vector)
         if show_now:
             plt.show()