BREAKING CHANGE: Move ode_solver from OptimalControlProgram to DynamicsList.add #959

EveCharbie · 2025-04-01T11:17:30Z

BREAKING CHANGE: Move ode_solver from OptimalControlProgram to DynamicsList.add

Summary:

This PR introduces a significant breaking change to the Bioptim API regarding the specification of the Ordinary Differential Equation (ODE) solver. The ode_solver argument has been removed from the OptimalControlProgram constructor and is now a required argument for each dynamics configuration added via DynamicsList.add().

Motivation:

Previously, a single ode_solver was specified at the OptimalControlProgram level, applying uniformly to all phases of a multiphase problem. This lacked flexibility for scenarios where different phases might benefit from different integration schemes (e.g., using a stiff solver for one phase and a standard Runge-Kutta for another).

Moving the ode_solver specification to the DynamicsList.add method:

Enhances Flexibility: Allows users to define a specific ODE solver for each phase independently.
Improves Logical Grouping: Associates the ODE solver directly with the dynamics function it integrates, making the configuration more intuitive.
Increases Consistency: Aligns the ode_solver setting with other phase-specific dynamics arguments like expand_dynamics and phase_dynamics.

⚠️ BREAKING CHANGE NOTICE FOR USERS ⚠️

This is a breaking change and will require modifications to existing Bioptim scripts.

If your code previously defined an OptimalControlProgram by passing the ode_solver argument, you will encounter an error with this new version. You must update your code as follows:

Remove the ode_solver argument from your OptimalControlProgram(...) call.
Add the ode_solver argument to each call where you define phase dynamics using dynamics.add(...).

Example Migration:

Let's illustrate with a snippet from a typical prepare_ocp function:

Code BEFORE this PR:

from bioptim import OptimalControlProgram, DynamicsList, DynamicsFcn, OdeSolver

# --- Inside prepare_ocp function ---
my_ode_solver = OdeSolver.RK4() # Or any other solver

# Dynamics defined WITHOUT the solver
dynamics = DynamicsList()
dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=True, phase=0)
# ... add more dynamics for other phases if needed
dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=True, phase=1)

# OCP defined WITH the global solver
ocp = OptimalControlProgram(
    bio_model=...,
    dynamics=dynamics,
    n_shooting=...,
    phase_time=...,
    objective_functions=...,
    constraints=...,
    ode_solver=my_ode_solver, # <--- ode_solver defined here for ALL phases
    control_type=...
)

Code AFTER this PR:

from bioptim import OptimalControlProgram, DynamicsList, DynamicsFcn, OdeSolver

# --- Inside prepare_ocp function ---
# Define the solver(s) needed
ode_solver_phase_0 = OdeSolver.RK4()
ode_solver_phase_1 = OdeSolver.RK4() # Can be the same or different
# Example: ode_solver_phase_1 = OdeSolver.COLLOCATION(polynomial_degree=4)

# Dynamics defined WITH the solver for EACH phase
dynamics = DynamicsList()
dynamics.add(
    DynamicsFcn.TORQUE_DRIVEN,
    expand_dynamics=True,
    phase=0,
    ode_solver=ode_solver_phase_0 # <--- ode_solver specified HERE for phase 0
)
# ... add more dynamics for other phases if needed
dynamics.add(
    DynamicsFcn.TORQUE_DRIVEN,
    expand_dynamics=True,
    phase=1,
    ode_solver=ode_solver_phase_1 # <--- ode_solver specified HERE for phase 1
)

# OCP defined WITHOUT the global solver argument
ocp = OptimalControlProgram(
    bio_model=...,
    dynamics=dynamics,
    n_shooting=...,
    phase_time=...,
    objective_functions=...,
    constraints=...,
    # ode_solver=... <--- REMOVED from here
    control_type=...
)

# Note: If using the same solver for all phases, you can still define it once
# and pass the same variable to each dynamics.add() call:
# common_ode_solver = OdeSolver.RK4()
# dynamics.add(..., ode_solver=common_ode_solver)
# dynamics.add(..., ode_solver=common_ode_solver)

This change is

EveCharbie · 2025-04-01T17:45:30Z

Only ACADOS is failing.
@pariterre @Ipuch could one of you guys review this offline please ?

Ipuch

Reviewed 7 of 72 files at r1, all commit messages.
Reviewable status: 7 of 72 files reviewed, 10 unresolved discussions

bioptim/__init__.py line 185 at r1 (raw file):

from .misc.__version__ import __version__
from .misc.enums import (
    Axis,

I'm not sure about this refactor...

Code quote:

from .dynamics.configure_problem import ConfigureProblem, DynamicsFcn, DynamicsList, Dynamics
from .dynamics.dynamics_evaluation import DynamicsEvaluation
from .dynamics.dynamics_functions import DynamicsFunctions
from .dynamics.fatigue.effort_perception import EffortPerception, TauEffortPerception
from .dynamics.fatigue.fatigue_dynamics import FatigueList
from .dynamics.fatigue.michaud_fatigue import MichaudFatigue, MichaudTauFatigue
from .dynamics.fatigue.xia_fatigue import XiaFatigue, XiaTauFatigue, XiaFatigueStabilized
from .dynamics.ode_solvers import OdeSolver, OdeSolverBase
from .gui.online_callback_server import PlottingServer
from .gui.plot import CustomPlot
from .interfaces import Solver
from .limits.constraints import ConstraintFcn, ConstraintList, Constraint, ParameterConstraintList
from .limits.fatigue_path_conditions import FatigueBounds, FatigueInitialGuess
from .limits.multinode_constraint import MultinodeConstraintFcn, MultinodeConstraintList, MultinodeConstraint
from .limits.multinode_objective import MultinodeObjectiveFcn, MultinodeObjectiveList, MultinodeObjective
from .limits.objective_functions import ObjectiveFcn, ObjectiveList, Objective, ParameterObjectiveList
from .limits.path_conditions import BoundsList, InitialGuessList, Bounds, InitialGuess
from .limits.penalty_controller import PenaltyController
from .limits.penalty_helpers import PenaltyHelpers
from .limits.phase_transition import PhaseTransitionFcn, PhaseTransitionList, PhaseTransition
from .misc.__version__ import __version__
from .misc.enums import (
    Axis,

bioptim/optimization/optimal_control_program.py line 523 at r1 (raw file):

            is_ode_solver = isinstance(dynamics[i_dyn].ode_solver, OdeSolverBase)
            if not is_ode_solver:
                raise RuntimeError("ode_solver should be built an instance of OdeSolver")

It should be a method of Dynamic

Code quote:

            if dyn.ode_solver is None:
                dynamics[i_dyn].ode_solver = self._set_default_ode_solver()

            is_ode_solver = isinstance(dynamics[i_dyn].ode_solver, OdeSolverBase)
            if not is_ode_solver:
                raise RuntimeError("ode_solver should be built an instance of OdeSolver")

bioptim/optimization/optimal_control_program.py line 580 at r1 (raw file):

        # Prepare path constraints and dynamics of the program
        NLP.add(self, "dynamics_type", dynamics, False)
        ode_solver = [dyn.ode_solver for dyn in dynamics]

Should it be method of DynamicsList ?

bioptim/optimization/optimal_control_program.py line 522 at r1 (raw file):

        )
        if not is_ode_solver and not is_list_ode_solver:
            raise RuntimeError("ode_solver should be built an instance of OdeSolver or a list of OdeSolver")

If the idea of this PR is preserved, we can put these checks in some utils functions to move import them wherever we want.

Code quote:

        if ode_solver is None:
            ode_solver = self._set_default_ode_solver()

        is_ode_solver = isinstance(ode_solver, OdeSolverBase)
        is_list_ode_solver = (
            all([isinstance(ode, OdeSolverBase) for ode in ode_solver])
            if isinstance(ode_solver, list) or isinstance(ode_solver, tuple)
            else False
        )
        if not is_ode_solver and not is_list_ode_solver:
            raise RuntimeError("ode_solver should be built an instance of OdeSolver or a list of OdeSolver")

bioptim/optimization/stochastic_optimal_control_program.py line 637 at r1 (raw file):

            "duplicate_starting_point=True"
            ")"
        )

You don't check it anymore ? somewhere else ? (I havent read the whole PR at this stage.)

Code quote:

def _check_has_no_ode_solver_defined(**kwargs):
    if "ode_solver" in kwargs:
        raise ValueError(
            "The ode_solver cannot be defined for a stochastic ocp. "
            "The value is chosen based on the type of problem solved:"
            "\n- TRAPEZOIDAL_EXPLICIT: OdeSolver.TRAPEZOIDAL() "
            "\n- TRAPEZOIDAL_IMPLICIT: OdeSolver.TRAPEZOIDAL() "
            "\n- COLLOCATION: "
            "OdeSolver.COLLOCATION("
            "method=problem_type.method, "
            "polynomial_degree=problem_type.polynomial_degree, "
            "duplicate_starting_point=True"
            ")"
        )

bioptim/optimization/variational_optimal_control_program.py line 74 at r1 (raw file):

                "ode_solver cannot be defined in VariationalOptimalControlProgram since the integration is"
                " done by the variational integrator."
            )

for memory

Code quote:

        if "ode_solver" in kwargs:
            raise ValueError(
                "ode_solver cannot be defined in VariationalOptimalControlProgram since the integration is"
                " done by the variational integrator."
            )

bioptim/optimization/non_linear_program.py line 289 at r1 (raw file):

            if not_direct_collocation and (x_init_all_point or a_init_all_point):
                raise ValueError("InterpolationType.ALL_POINTS must only be used with direct collocation")

thanks

Code quote:

        if x_init is not None or a_init is not None:
            not_direct_collocation = not self.ode_solver.is_direct_collocation
            x_init_all_point = x_init.type == InterpolationType.ALL_POINTS if x_init is not None else False
            a_init_all_point = a_init.type == InterpolationType.ALL_POINTS if a_init is not None else False

            if not_direct_collocation and (x_init_all_point or a_init_all_point):
                raise ValueError("InterpolationType.ALL_POINTS must only be used with direct collocation")

bioptim/optimization/stochastic_optimal_control_program.py line 91 at r1 (raw file):

        # Integrator
        for dyn in dynamics:
            dyn.ode_solver = self._set_default_ode_solver()

I would expect a:

dyn.check_has_no_ode_defined(...)

bioptim/optimization/non_linear_program.py line 149 at r1 (raw file):

        self.n_threads = None
        self.ns = None
        self.ode_solver = OdeSolver.RK4()

pretty dangerous now. I would be conservative and try something to set this attribute from self.dynamics to transfer the responsability to nlp.

bioptim/optimization/non_linear_program.py line 436 at r1 (raw file):

            return 1
        if self.ode_solver.is_direct_collocation:
            return self.dynamics[node_idx].shape_xall[1] - (1 if not self.ode_solver.duplicate_starting_point else 0)

Preserve the last version should be a priority, I think. But it might hide something more deep on why these methods are here ...

Code quote:

    def n_states_stepwise_steps(self, node_idx) -> int:
        """
        Parameters
        ----------
        node_idx: int
            The index of the node

        Returns
        -------
        The number of states
        """
        if node_idx >= self.ns:
            return 1
        if self.ode_solver.is_direct_collocation:
            return self.dynamics[node_idx].shape_xall[1] - (1 if not self.ode_solver.duplicate_starting_point else 0)

Ipuch · 2025-04-03T13:13:17Z

I really want to take the benefits of this refactor. Because this is a decision with a lot of impact. I reviewed only a part of it to see the main consequences. But the PR doesn't seem to justify why it's absolutely wanted. Don't work more on this one before we talk about it.

Ipuch · 2025-04-03T14:54:26Z

Disussing the best choice with EVE :

#Solution RETENUE

ode_solver = OdeSolverList()
ode_solver.add(OdeSolver.RK4, Dynamics(TORQUE_DRIVEN, contact_type, with_ ...))
ode_solver.add(OdeSolver.RK4, Dynamics(TORQUE_DRIVEN, contact_type, with_ ...))
ode_solver.integrate()

EveCharbie

Reviewable status: 7 of 72 files reviewed, 9 unresolved discussions (waiting on @Ipuch)

bioptim/__init__.py line 185 at r1 (raw file):