Add the choice for an integrator of the solution

bioptim/__init__.py

@@ -165,7 +165,7 @@
 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 = False,
     ):
         """
         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
+        self.use_scipy_integrator = self.integrator.value is not None
         if isinstance(self.ocp.original_phase_time, (int, float)):
             self.tf = [self.ocp.original_phase_time]
         else:
@@ -585,7 +587,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

bioptim/misc/enums.py

@@ -103,3 +103,17 @@ class VariableType(Enum):
 
     STATES = "states"
     CONTROLS = "controls"
+
+
+class SolutionIntegrator(Enum):
+    """
+    Selection of integrator to use integrate function
+    """
+    DEFAULT = None
+    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,9 +548,12 @@ 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])
+
+        scipy_integrator = integrator.value
+        use_scipy_integrator = scipy_integrator is not None
+
         if n_direct_collocation > 0 and not use_scipy_integrator:
             if continuous:
                 raise RuntimeError(
@@ -613,8 +616,9 @@ def __perform_integration(
                     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
-                    ).y
+                        lambda t, x: np.array(nlp.dynamics_func(x, u, params))[:, 0], [t_init, t_end], x0,
+                        t_eval=t_eval,
+                        method=scipy_integrator, atol=1e-10, rtol=1e-10).y
 
                     next_state_col = (
                         (s + 1) * (nlp.ode_solver.steps + 1) if nlp.ode_solver.is_direct_collocation else s + 1
@@ -853,7 +857,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
@@ -875,7 +879,7 @@ def graphs(
         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()

tests/test_simulate.py

@@ -4,7 +4,7 @@
 import pytest
 
 import numpy as np
-from bioptim import Shooting, OdeSolver
+from bioptim import Shooting, OdeSolver, SolutionIntegrator
 
 
 def test_merge_phases_one_phase():
@@ -167,8 +167,8 @@ def test_interpolate_multiphases_merge_phase():
 
 
 @pytest.mark.parametrize("ode_solver", [OdeSolver.RK4, OdeSolver.COLLOCATION])
-@pytest.mark.parametrize("use_scipy", [True, False])
-def test_integrate(use_scipy, ode_solver):
+@pytest.mark.parametrize("integrator", [SolutionIntegrator.SCIPY_RK45, SolutionIntegrator.DEFAULT])
+def test_integrate(integrator, ode_solver):
     # Load pendulum
     from bioptim.examples.getting_started import pendulum as ocp_module
 
@@ -185,7 +185,7 @@ def test_integrate(use_scipy, ode_solver):
 
     sol = ocp.solve()
 
-    opts = {"shooting_type": Shooting.MULTIPLE, "keep_intermediate_points": False, "use_scipy_integrator": use_scipy}
+    opts = {"shooting_type": Shooting.MULTIPLE, "keep_intermediate_points": False, "integrator": integrator}
     with pytest.raises(
         ValueError,
         match="Shooting.MULTIPLE and keep_intermediate_points=False "
@@ -194,15 +194,15 @@ def test_integrate(use_scipy, ode_solver):
         _ = sol.integrate(**opts)
 
     opts["keep_intermediate_points"] = True
-    if ode_solver == OdeSolver.COLLOCATION and not use_scipy:
+    if ode_solver == OdeSolver.COLLOCATION and integrator.value is None:
         with pytest.raises(RuntimeError, match="Integration with direct collocation must be not continuous"):
             sol.integrate(**opts)
         return
 
     sol_integrated = sol.integrate(**opts)
     shapes = (4, 2, 2)
 
-    decimal = 5 if use_scipy else 8
+    decimal = 5 if integrator.value is not None else 8
     for i, key in enumerate(sol.states):
         np.testing.assert_almost_equal(
             sol_integrated.states[key][:, [0, -1]], sol.states[key][:, [0, -1]], decimal=decimal
@@ -241,7 +241,7 @@ def test_integrate_single_shoot(keep_intermediate_points, ode_solver):
 
     sol = ocp.solve()
 
-    opts = {"keep_intermediate_points": keep_intermediate_points, "use_scipy_integrator": False}
+    opts = {"keep_intermediate_points": keep_intermediate_points, "integrator": SolutionIntegrator.DEFAULT}
     if ode_solver == OdeSolver.COLLOCATION:
         with pytest.raises(RuntimeError, match="Integration with direct collocation must be not continuous"):
             sol.integrate(**opts)
@@ -296,7 +296,7 @@ def test_integrate_single_shoot_use_scipy(keep_intermediate_points, ode_solver):
 
     sol = ocp.solve()
 
-    opts = {"keep_intermediate_points": keep_intermediate_points, "use_scipy_integrator": True}
+    opts = {"keep_intermediate_points": keep_intermediate_points, "integrator": SolutionIntegrator.SCIPY_RK45}
 
     sol_integrated = sol.integrate(**opts)
     shapes = (4, 2, 2)
@@ -436,12 +436,13 @@ def test_integrate_single_shoot_use_scipy(keep_intermediate_points, ode_solver):
 @pytest.mark.parametrize("ode_solver", [OdeSolver.RK4, OdeSolver.COLLOCATION])
 @pytest.mark.parametrize("shooting", [Shooting.SINGLE_CONTINUOUS, Shooting.MULTIPLE, Shooting.SINGLE])
 @pytest.mark.parametrize("merge", [False, True])
-@pytest.mark.parametrize("use_scipy", [False, True])
-def test_integrate_non_continuous(shooting, merge, use_scipy, ode_solver):
+@pytest.mark.parametrize("integrator", [SolutionIntegrator.DEFAULT, SolutionIntegrator.SCIPY_RK45])
+def test_integrate_non_continuous(shooting, merge, integrator, ode_solver):
     # Load pendulum
     from bioptim.examples.getting_started import pendulum as ocp_module
 
     bioptim_folder = os.path.dirname(ocp_module.__file__)
+    use_scipy = integrator is not None
 
     n_shooting = 10
 
@@ -458,7 +459,7 @@ def test_integrate_non_continuous(shooting, merge, use_scipy, ode_solver):
         "shooting_type": shooting,
         "continuous": False,
         "keep_intermediate_points": False,
-        "use_scipy_integrator": use_scipy,
+        "integrator": integrator,
     }
     if shooting == Shooting.SINGLE_CONTINUOUS:
         with pytest.raises(
@@ -515,8 +516,8 @@ def test_integrate_non_continuous(shooting, merge, use_scipy, ode_solver):
 @pytest.mark.parametrize("ode_solver", [OdeSolver.RK4, OdeSolver.COLLOCATION])
 @pytest.mark.parametrize("shooting", [Shooting.SINGLE_CONTINUOUS, Shooting.MULTIPLE, Shooting.SINGLE])
 @pytest.mark.parametrize("keep_intermediate_points", [True, False])
-@pytest.mark.parametrize("use_scipy", [False, True])
-def test_integrate_multiphase(shooting, keep_intermediate_points, use_scipy, ode_solver):
+@pytest.mark.parametrize("integrator", [SolutionIntegrator.DEFAULT, SolutionIntegrator.SCIPY_RK45])
+def test_integrate_multiphase(shooting, keep_intermediate_points, integrator, ode_solver):
     # Load pendulum
     from bioptim.examples.getting_started import example_multiphase as ocp_module
 
@@ -527,11 +528,12 @@ def test_integrate_multiphase(shooting, keep_intermediate_points, use_scipy, ode
     sol = ocp.solve()
     n_shooting = [20, 30, 20]
 
+    use_scipy = integrator is not None
     opts = {
         "shooting_type": shooting,
         "continuous": False,
         "keep_intermediate_points": keep_intermediate_points,
-        "use_scipy_integrator": use_scipy,
+        "integrator": integrator,
     }
 
     if shooting == Shooting.SINGLE_CONTINUOUS:
@@ -603,8 +605,8 @@ def test_integrate_multiphase(shooting, keep_intermediate_points, use_scipy, ode
 @pytest.mark.parametrize("ode_solver", [OdeSolver.RK4, OdeSolver.COLLOCATION])
 @pytest.mark.parametrize("shooting", [Shooting.SINGLE_CONTINUOUS, Shooting.MULTIPLE, Shooting.SINGLE])
 @pytest.mark.parametrize("keep_intermediate_points", [True, False])
-@pytest.mark.parametrize("use_scipy", [False, True])
-def test_integrate_multiphase_merged(shooting, keep_intermediate_points, use_scipy, ode_solver):
+@pytest.mark.parametrize("integrator", [SolutionIntegrator.DEFAULT, SolutionIntegrator.SCIPY_RK45])
+def test_integrate_multiphase_merged(shooting, keep_intermediate_points, integrator, ode_solver):
     # Load pendulum
     from bioptim.examples.getting_started import example_multiphase as ocp_module
 
@@ -617,11 +619,12 @@ def test_integrate_multiphase_merged(shooting, keep_intermediate_points, use_sci
 
     sol = ocp.solve()
 
+    use_scipy= integrator.value is not None
     opts = {
         "shooting_type": shooting,
         "continuous": False,
         "keep_intermediate_points": keep_intermediate_points,
-        "use_scipy_integrator": use_scipy,
+        "integrator": integrator,
     }
 
     if shooting == Shooting.SINGLE_CONTINUOUS:
@@ -698,8 +701,8 @@ def test_integrate_multiphase_merged(shooting, keep_intermediate_points, use_sci
 
 @pytest.mark.parametrize("ode_solver", [OdeSolver.RK4, OdeSolver.COLLOCATION])
 @pytest.mark.parametrize("shooting", [Shooting.SINGLE_CONTINUOUS, Shooting.MULTIPLE, Shooting.SINGLE])
-@pytest.mark.parametrize("use_scipy", [False, True])
-def test_integrate_multiphase_non_continuous(shooting, use_scipy, ode_solver):
+@pytest.mark.parametrize("integrator", [SolutionIntegrator.DEFAULT, SolutionIntegrator.SCIPY_RK45])
+def test_integrate_multiphase_non_continuous(shooting, integrator, ode_solver):
     # Load pendulum
     from bioptim.examples.getting_started import example_multiphase as ocp_module
 
@@ -710,11 +713,12 @@ def test_integrate_multiphase_non_continuous(shooting, use_scipy, ode_solver):
     sol = ocp.solve()
     n_shooting = [20, 30, 20]
 
+    use_scipy = integrator.value is not None
     opts = {
         "shooting_type": shooting,
         "continuous": False,
         "keep_intermediate_points": True,
-        "use_scipy_integrator": use_scipy,
+        "integrator": integrator,
     }
 
     if shooting == Shooting.SINGLE_CONTINUOUS:
@@ -768,8 +772,8 @@ def test_integrate_multiphase_non_continuous(shooting, use_scipy, ode_solver):
 
 @pytest.mark.parametrize("ode_solver", [OdeSolver.RK4, OdeSolver.COLLOCATION])
 @pytest.mark.parametrize("shooting", [Shooting.SINGLE_CONTINUOUS, Shooting.MULTIPLE, Shooting.SINGLE])
-@pytest.mark.parametrize("use_scipy", [False, True])
-def test_integrate_multiphase_merged_non_continuous(shooting, use_scipy, ode_solver):
+@pytest.mark.parametrize("integrator", [SolutionIntegrator.DEFAULT, SolutionIntegrator.SCIPY_RK45])
+def test_integrate_multiphase_merged_non_continuous(shooting, integrator, ode_solver):
     # Load pendulum
     from bioptim.examples.getting_started import example_multiphase as ocp_module
 
@@ -779,11 +783,12 @@ def test_integrate_multiphase_merged_non_continuous(shooting, use_scipy, ode_sol
 
     sol = ocp.solve()
 
+    use_scipy = integrator.value is not None
     opts = {
         "shooting_type": shooting,
         "continuous": False,
         "keep_intermediate_points": False,
-        "use_scipy_integrator": use_scipy,
+        "integrator": integrator,
     }
     if shooting == Shooting.SINGLE_CONTINUOUS:
         with pytest.raises(