diff --git a/README.md b/README.md index 50a031f30..59a3b5432 100644 --- a/README.md +++ b/README.md @@ -707,7 +707,6 @@ OptimalControlProgram( objective_functions: [Objective, ObjectiveList], constraints: [Constraint, ConstraintList], parameters: ParameterList, - ode_solver: OdeSolver, control_type: [ControlType, list], all_generalized_mapping: BiMapping, q_mapping: BiMapping, @@ -734,8 +733,7 @@ In the case of a multiphase optimization, one model per phase should be passed i `objective_functions` is the objective function set of the ocp (see The objective functions section). `constraints` is the constraint set of the ocp (see The constraints section). `parameters` is the parameter set of the ocp (see The parameters section). -It is a list (one element for each phase) of np.ndarray of shape (6, i, n), where the 6 components are [Mx, My, Mz, Fx, Fy, Fz], for the ith force platform (defined by the externalforceindex) for each node n. -`ode_solver` is the ode solver used to solve the dynamic equations. +It is a list (one element for each phase) of np.ndarray of shape (6, i, n), where the 6 components are [Mx, My, Mz, Fx, Fy, Fz], for the ith force platform (defined by the externalforceindex) for each node n. `control_type` is the type of discretization of the controls (usually CONSTANT) (see ControlType section). `all_generalized_mapping` is used to reduce the number of degrees of freedom by linking them (see The mappings section). This one applies the same mapping to the generalized coordinates (*q*), velocities (*qdot*), and forces (*tau*). @@ -949,7 +947,7 @@ The `DynamicsFcn` is the one presented in the corresponding section below. #### The options The full signature of Dynamics is as follows: ```python -Dynamics(dynamics_type, configure: Callable, dynamic_function: Callable, phase: int) +Dynamics(dynamics_type, configure: Callable, dynamic_function: Callable, phase: int, ode_solver: OdeSolver, contact_type: list[ContactType], numerical_timeseries: dict[str, np.ndarray]) ``` The `dynamics_type` is the selected `DynamicsFcn`. It automatically defines both `configure` and `dynamic_function`. @@ -957,6 +955,9 @@ If a function is sent instead, this function is interpreted as `configure` and t If one is interested in changing the behavior of a particular `DynamicsFcn`, they can refer to the Custom dynamics functions right below. The `phase` is the index of the phase the dynamics applies to. +The `ode_solver` is the ode to use to "integrate" the dynamics function. +The `contact_type` is a list of contact types to consider in the dynamics equations. +The `numerical_timeseries` is a list of numerical values (one per node) to use in the dynamics. For example, it can be used to define experimental ground reaction forces. The `add()` method of `DynamicsList` usually takes care of this, but it can be useful when declaring the dynamics out of order. #### Custom dynamic functions diff --git a/bioptim/__init__.py b/bioptim/__init__.py index 92814908a..88abefc9e 100644 --- a/bioptim/__init__.py +++ b/bioptim/__init__.py @@ -161,46 +161,25 @@ """ from .dynamics.configure_problem import ConfigureProblem, DynamicsFcn, DynamicsList, Dynamics -from .dynamics.configure_problem import ConfigureProblem, DynamicsFcn, DynamicsList, Dynamics -from .dynamics.dynamics_evaluation import DynamicsEvaluation from .dynamics.dynamics_evaluation import DynamicsEvaluation from .dynamics.dynamics_functions import DynamicsFunctions -from .dynamics.dynamics_functions import DynamicsFunctions from .dynamics.fatigue.effort_perception import EffortPerception, TauEffortPerception -from .dynamics.fatigue.effort_perception import EffortPerception, TauEffortPerception -from .dynamics.fatigue.fatigue_dynamics import FatigueList from .dynamics.fatigue.fatigue_dynamics import FatigueList from .dynamics.fatigue.michaud_fatigue import MichaudFatigue, MichaudTauFatigue -from .dynamics.fatigue.michaud_fatigue import MichaudFatigue, MichaudTauFatigue from .dynamics.fatigue.xia_fatigue import XiaFatigue, XiaTauFatigue, XiaFatigueStabilized -from .dynamics.fatigue.xia_fatigue import XiaFatigue, XiaTauFatigue, XiaFatigueStabilized -from .dynamics.ode_solvers import OdeSolver, OdeSolverBase from .dynamics.ode_solvers import OdeSolver, OdeSolverBase from .gui.online_callback_server import PlottingServer -from .gui.online_callback_server import PlottingServer from .gui.plot import CustomPlot -from .gui.plot import CustomPlot -from .interfaces import Solver from .interfaces import Solver from .limits.constraints import ConstraintFcn, ConstraintList, Constraint, ParameterConstraintList -from .limits.constraints import ConstraintFcn, ConstraintList, Constraint, ParameterConstraintList -from .limits.fatigue_path_conditions import FatigueBounds, FatigueInitialGuess from .limits.fatigue_path_conditions import FatigueBounds, FatigueInitialGuess from .limits.multinode_constraint import MultinodeConstraintFcn, MultinodeConstraintList, MultinodeConstraint -from .limits.multinode_constraint import MultinodeConstraintFcn, MultinodeConstraintList, MultinodeConstraint from .limits.multinode_objective import MultinodeObjectiveFcn, MultinodeObjectiveList, MultinodeObjective -from .limits.multinode_objective import MultinodeObjectiveFcn, MultinodeObjectiveList, MultinodeObjective -from .limits.objective_functions import ObjectiveFcn, ObjectiveList, Objective, ParameterObjectiveList from .limits.objective_functions import ObjectiveFcn, ObjectiveList, Objective, ParameterObjectiveList from .limits.path_conditions import BoundsList, InitialGuessList, Bounds, InitialGuess -from .limits.path_conditions import BoundsList, InitialGuessList, Bounds, InitialGuess -from .limits.penalty_controller import PenaltyController from .limits.penalty_controller import PenaltyController from .limits.penalty_helpers import PenaltyHelpers -from .limits.penalty_helpers import PenaltyHelpers from .limits.phase_transition import PhaseTransitionFcn, PhaseTransitionList, PhaseTransition -from .limits.phase_transition import PhaseTransitionFcn, PhaseTransitionList, PhaseTransition -from .misc.__version__ import __version__ from .misc.__version__ import __version__ from .misc.enums import ( Axis, diff --git a/bioptim/dynamics/configure_new_variable.py b/bioptim/dynamics/configure_new_variable.py index ec62b3cca..431b7db67 100644 --- a/bioptim/dynamics/configure_new_variable.py +++ b/bioptim/dynamics/configure_new_variable.py @@ -214,7 +214,7 @@ def _declare_cx_and_plot(self): for node_index in range( self.nlp.n_states_nodes if self.nlp.phase_dynamics == PhaseDynamics.ONE_PER_NODE else 1 ): - n_cx = self.nlp.ode_solver.n_required_cx + 2 + n_cx = self.nlp.dynamics_type.ode_solver.n_required_cx + 2 cx_scaled = self.define_cx_scaled(n_col=n_cx, node_index=node_index) cx = self.define_cx_unscaled(cx_scaled, self.nlp.x_scaling[self.name].scaling) self.nlp.states.append( @@ -246,7 +246,7 @@ def _declare_cx_and_plot(self): for node_index in range( self.nlp.n_states_nodes if self.nlp.phase_dynamics == PhaseDynamics.ONE_PER_NODE else 1 ): - n_cx = self.nlp.ode_solver.n_required_cx + 2 + n_cx = self.nlp.dynamics_type.ode_solver.n_required_cx + 2 cx_scaled = self.define_cx_scaled(n_col=n_cx, node_index=node_index) cx = self.define_cx_unscaled(cx_scaled, np.ones_like(self.nlp.x_scaling[self.name].scaling)) self.nlp.states_dot.append( @@ -294,7 +294,7 @@ def _declare_cx_and_plot(self): for node_index in range( self.nlp.n_states_nodes if self.nlp.phase_dynamics == PhaseDynamics.ONE_PER_NODE else 1 ): - n_cx = self.nlp.ode_solver.n_required_cx + 2 + n_cx = self.nlp.dynamics_type.ode_solver.n_required_cx + 2 cx_scaled = self.define_cx_scaled(n_col=n_cx, node_index=node_index) cx = self.define_cx_unscaled(cx_scaled, self.nlp.a_scaling[self.name].scaling) self.nlp.algebraic_states.append( diff --git a/bioptim/dynamics/configure_problem.py b/bioptim/dynamics/configure_problem.py index 061320741..d61eaf57f 100644 --- a/bioptim/dynamics/configure_problem.py +++ b/bioptim/dynamics/configure_problem.py @@ -6,7 +6,7 @@ from .configure_new_variable import NewVariableConfiguration from .dynamics_functions import DynamicsFunctions from .fatigue.fatigue_dynamics import FatigueList -from .ode_solvers import OdeSolver +from .ode_solvers import OdeSolver, OdeSolverBase from ..gui.plot import CustomPlot from ..limits.constraints import ImplicitConstraintFcn from ..misc.enums import ( @@ -1356,7 +1356,11 @@ def configure_integrated_value( # TODO: compute values at collocation points # but for now only cx_start can be used - n_cx = nlp.ode_solver.n_cx - 1 if isinstance(nlp.ode_solver, OdeSolver.COLLOCATION) else 3 + n_cx = ( + nlp.dynamics_type.ode_solver.n_cx - 1 + if isinstance(nlp.dynamics_type.ode_solver, OdeSolver.COLLOCATION) + else 3 + ) if n_cx < 3: n_cx = 3 @@ -1961,6 +1965,7 @@ def __init__( skip_continuity: bool = False, state_continuity_weight: float | None = None, phase_dynamics: PhaseDynamics = PhaseDynamics.SHARED_DURING_THE_PHASE, + ode_solver: OdeSolver | OdeSolverBase = OdeSolver.RK4(), numerical_data_timeseries: dict[str, np.ndarray] = None, contact_type: list[ContactType] | tuple[ContactType] = (), **extra_parameters: Any, @@ -1983,6 +1988,8 @@ def __init__( otherwise it is an objective phase_dynamics: PhaseDynamics If the dynamics should be shared between the nodes or not + ode_solver: OdeSolver + The integrator to use to integrate this dynamics. numerical_data_timeseries: dict[str, np.ndarray] The numerical timeseries at each node. ex: the experimental external forces data should go here. contact_type: list[ContactType] | tuple[ContactType] @@ -2003,6 +2010,9 @@ def __init__( dynamic_function = extra_parameters["dynamic_function"] del extra_parameters["dynamic_function"] + if not isinstance(ode_solver, OdeSolverBase): + raise RuntimeError("ode_solver should be built an instance of OdeSolver") + super(Dynamics, self).__init__(type=dynamics_type, **extra_parameters) self.dynamic_function = dynamic_function self.configure = configure @@ -2011,6 +2021,7 @@ def __init__( self.skip_continuity = skip_continuity self.state_continuity_weight = state_continuity_weight self.phase_dynamics = phase_dynamics + self.ode_solver = ode_solver self.numerical_data_timeseries = numerical_data_timeseries self.contact_type = contact_type @@ -2038,7 +2049,6 @@ def add(self, dynamics_type: Callable | Dynamics | DynamicsFcn, **extra_paramete extra_parameters: dict Any parameters to pass to Dynamics """ - if isinstance(dynamics_type, Dynamics): self.copy(dynamics_type) diff --git a/bioptim/dynamics/dynamics_functions.py b/bioptim/dynamics/dynamics_functions.py index 592a0a7df..d8925ceec 100644 --- a/bioptim/dynamics/dynamics_functions.py +++ b/bioptim/dynamics/dynamics_functions.py @@ -154,7 +154,7 @@ def torque_driven( defects = None # TODO: contacts and fatigue to be handled with implicit dynamics - if nlp.ode_solver.defects_type == DefectType.IMPLICIT: + if nlp.dynamics_type.ode_solver.defects_type == DefectType.IMPLICIT: if len(contact_type) == 0 and fatigue is None: qddot = DynamicsFunctions.get(nlp.states_dot["qdot"], nlp.states_dot.scaled.cx) tau_id = DynamicsFunctions.inverse_dynamics(nlp, q, qdot, qddot, contact_type, external_forces) diff --git a/bioptim/dynamics/integrator.py b/bioptim/dynamics/integrator.py index 399abbd9e..33cdd4786 100644 --- a/bioptim/dynamics/integrator.py +++ b/bioptim/dynamics/integrator.py @@ -785,3 +785,11 @@ class CVODES(Integrator): Class for CVODES integrators """ + + +class VARIATIONAL(RK4): + """ + Fake class for variational integrator. + The real work is done in VariationalOptimalControlProgram. + TODO: The implementation of the variational integrator could be moved here (see issue #962). + """ diff --git a/bioptim/dynamics/ode_solver_base.py b/bioptim/dynamics/ode_solver_base.py index d783a4f06..c271de288 100644 --- a/bioptim/dynamics/ode_solver_base.py +++ b/bioptim/dynamics/ode_solver_base.py @@ -229,7 +229,7 @@ def initialize_integrator( "implicit_ode": nlp.implicit_dynamics_func, } - return nlp.ode_solver.integrator(ode, ode_opt) + return nlp.dynamics_type.ode_solver.integrator(ode, ode_opt) def prepare_dynamic_integrator(self, ocp, nlp): """ @@ -244,26 +244,32 @@ def prepare_dynamic_integrator(self, ocp, nlp): """ # Primary dynamics - dynamics = [nlp.ode_solver.initialize_integrator(ocp, nlp, dynamics_index=0, node_index=0)] + dynamics = [nlp.dynamics_type.ode_solver.initialize_integrator(ocp, nlp, dynamics_index=0, node_index=0)] if nlp.phase_dynamics == PhaseDynamics.SHARED_DURING_THE_PHASE: dynamics = dynamics * nlp.ns else: for node_index in range(1, nlp.ns): - dynamics.append(nlp.ode_solver.initialize_integrator(ocp, nlp, dynamics_index=0, node_index=node_index)) + dynamics.append( + nlp.dynamics_type.ode_solver.initialize_integrator( + ocp, nlp, dynamics_index=0, node_index=node_index + ) + ) nlp.dynamics = dynamics # Extra dynamics extra_dynamics = [] for i in range(len(nlp.extra_dynamics_func)): extra_dynamics += [ - nlp.ode_solver.initialize_integrator(ocp, nlp, dynamics_index=i, node_index=0, is_extra_dynamics=True) + nlp.dynamics_type.ode_solver.initialize_integrator( + ocp, nlp, dynamics_index=i, node_index=0, is_extra_dynamics=True + ) ] if nlp.phase_dynamics == PhaseDynamics.SHARED_DURING_THE_PHASE: extra_dynamics = extra_dynamics * nlp.ns else: for node_index in range(1, nlp.ns): extra_dynamics += [ - nlp.ode_solver.initialize_integrator( + nlp.dynamics_type.ode_solver.initialize_integrator( ocp, nlp, dynamics_index=i, node_index=node_index, is_extra_dynamics=True ) ] diff --git a/bioptim/dynamics/ode_solvers.py b/bioptim/dynamics/ode_solvers.py index 3c4b426b6..13000dc01 100644 --- a/bioptim/dynamics/ode_solvers.py +++ b/bioptim/dynamics/ode_solvers.py @@ -49,6 +49,16 @@ class RK8(RK): def integrator(self): return integrator.RK8 + class VARIATIONAL(RK): + """ + This is a fake enum to be able to use the variational integrator. + TODO: The implementation of the variational integrator could be moved here (see issue #962). + """ + + @property + def integrator(self): + return integrator.VARIATIONAL + class TRAPEZOIDAL(OdeSolverBase): """ A trapezoidal ode solver diff --git a/bioptim/examples/acados/cube.py b/bioptim/examples/acados/cube.py index cf80559c7..234310656 100644 --- a/bioptim/examples/acados/cube.py +++ b/bioptim/examples/acados/cube.py @@ -24,7 +24,7 @@ def prepare_ocp(biorbd_model_path, n_shooting, tf, ode_solver=OdeSolver.RK4(), u bio_model = BiorbdModel(biorbd_model_path) # Dynamics - dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics) + dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics) # Path constraint x_bounds = BoundsList() @@ -43,7 +43,6 @@ def prepare_ocp(biorbd_model_path, n_shooting, tf, ode_solver=OdeSolver.RK4(), u tf, x_bounds=x_bounds, u_bounds=u_bounds, - ode_solver=ode_solver, use_sx=use_sx, ) diff --git a/bioptim/examples/custom_model/main.py b/bioptim/examples/custom_model/main.py index 55452e65e..01d1eb18d 100644 --- a/bioptim/examples/custom_model/main.py +++ b/bioptim/examples/custom_model/main.py @@ -73,6 +73,7 @@ def prepare_ocp( dynamics.add( configure_dynamics, dynamic_function=dynamics, + ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics, ) @@ -111,7 +112,6 @@ def prepare_ocp( x_init=x_init, u_init=u_init, objective_functions=objective_functions, - ode_solver=ode_solver, use_sx=False, n_threads=n_threads, ) diff --git a/bioptim/examples/fatigue/static_arm_with_fatigue.py b/bioptim/examples/fatigue/static_arm_with_fatigue.py index 3aea8839c..9ceb21a3c 100644 --- a/bioptim/examples/fatigue/static_arm_with_fatigue.py +++ b/bioptim/examples/fatigue/static_arm_with_fatigue.py @@ -156,6 +156,7 @@ def prepare_ocp( expand_dynamics=expand_dynamics, fatigue=fatigue_dynamics, with_residual_torque=torque_level > 0, + ode_solver=ode_solver, phase_dynamics=phase_dynamics, ) @@ -208,7 +209,6 @@ def prepare_ocp( u_init=u_init, objective_functions=objective_functions, constraints=constraint, - ode_solver=ode_solver, use_sx=False, n_threads=n_threads, ) diff --git a/bioptim/examples/getting_started/custom_constraint.py b/bioptim/examples/getting_started/custom_constraint.py index 25520f97e..3eef913ec 100644 --- a/bioptim/examples/getting_started/custom_constraint.py +++ b/bioptim/examples/getting_started/custom_constraint.py @@ -101,7 +101,9 @@ def prepare_ocp( objective_functions = Objective(ObjectiveFcn.Lagrange.MINIMIZE_CONTROL, key="tau", weight=100) # Dynamics - dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics = Dynamics( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) # Constraints constraints = ConstraintList() @@ -129,7 +131,6 @@ def prepare_ocp( u_bounds=u_bounds, objective_functions=objective_functions, constraints=constraints, - ode_solver=ode_solver, ) diff --git a/bioptim/examples/getting_started/custom_dynamics.py b/bioptim/examples/getting_started/custom_dynamics.py index e02134b47..a2d656f8a 100644 --- a/bioptim/examples/getting_started/custom_dynamics.py +++ b/bioptim/examples/getting_started/custom_dynamics.py @@ -166,6 +166,7 @@ def prepare_ocp( custom_configure, dynamic_function=custom_dynamics, my_additional_factor=1, + ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics, ) @@ -173,6 +174,7 @@ def prepare_ocp( dynamics.add( DynamicsFcn.TORQUE_DRIVEN, dynamic_function=custom_dynamics, + ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics, ) @@ -206,7 +208,6 @@ def prepare_ocp( u_bounds=u_bounds, objective_functions=objective_functions, constraints=constraints, - ode_solver=ode_solver, use_sx=use_sx, ) diff --git a/bioptim/examples/getting_started/custom_initial_guess.py b/bioptim/examples/getting_started/custom_initial_guess.py index 5216b19e6..e83c84a6f 100644 --- a/bioptim/examples/getting_started/custom_initial_guess.py +++ b/bioptim/examples/getting_started/custom_initial_guess.py @@ -126,7 +126,9 @@ def prepare_ocp( objective_functions = Objective(ObjectiveFcn.Lagrange.MINIMIZE_CONTROL, key="tau", weight=100) # Dynamics - dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics = Dynamics( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) # Constraints constraints = ConstraintList() @@ -244,7 +246,6 @@ def prepare_ocp( u_init=u_init, objective_functions=objective_functions, constraints=constraints, - ode_solver=ode_solver, x_scaling=x_scaling, u_scaling=u_scaling, ) diff --git a/bioptim/examples/getting_started/custom_objectives.py b/bioptim/examples/getting_started/custom_objectives.py index ef585a2b2..812d72de8 100644 --- a/bioptim/examples/getting_started/custom_objectives.py +++ b/bioptim/examples/getting_started/custom_objectives.py @@ -120,7 +120,9 @@ def prepare_ocp( ) # Dynamics - dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics = Dynamics( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) # Path constraint x_bounds = BoundsList() @@ -145,7 +147,6 @@ def prepare_ocp( x_bounds=x_bounds, u_bounds=u_bounds, objective_functions=objective_functions, - ode_solver=ode_solver, ) diff --git a/bioptim/examples/getting_started/custom_parameters.py b/bioptim/examples/getting_started/custom_parameters.py index 538a95d74..c0aa78c72 100644 --- a/bioptim/examples/getting_started/custom_parameters.py +++ b/bioptim/examples/getting_started/custom_parameters.py @@ -221,7 +221,9 @@ def prepare_ocp( objective_functions.add(ObjectiveFcn.Lagrange.MINIMIZE_STATE, key="q", weight=1) # Dynamics - dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics = Dynamics( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) # Path constraint x_bounds = BoundsList() @@ -249,7 +251,6 @@ def prepare_ocp( parameter_objectives=parameter_objectives, parameter_bounds=parameter_bounds, parameter_init=parameter_init, - ode_solver=ode_solver, use_sx=use_sx, ) diff --git a/bioptim/examples/getting_started/custom_phase_transitions.py b/bioptim/examples/getting_started/custom_phase_transitions.py index fe53fe0ed..1dc94826f 100644 --- a/bioptim/examples/getting_started/custom_phase_transitions.py +++ b/bioptim/examples/getting_started/custom_phase_transitions.py @@ -115,10 +115,18 @@ def prepare_ocp( # Dynamics dynamics = DynamicsList() - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) # Constraints constraints = ConstraintList() @@ -201,7 +209,6 @@ def prepare_ocp( x_init=x_init, objective_functions=objective_functions, constraints=constraints, - ode_solver=ode_solver, phase_transitions=phase_transitions, ) diff --git a/bioptim/examples/getting_started/example_continuity_as_objective.py b/bioptim/examples/getting_started/example_continuity_as_objective.py index 2768a0bf8..22e35da6f 100644 --- a/bioptim/examples/getting_started/example_continuity_as_objective.py +++ b/bioptim/examples/getting_started/example_continuity_as_objective.py @@ -111,6 +111,7 @@ def prepare_ocp_first_pass( dynamics = Dynamics( DynamicsFcn.TORQUE_DRIVEN, state_continuity_weight=state_continuity_weight, + ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics, ) @@ -162,7 +163,6 @@ def prepare_ocp_first_pass( u_bounds=u_bounds, objective_functions=objective_functions, constraints=constraints, - ode_solver=ode_solver, use_sx=use_sx, n_threads=n_threads, ) @@ -211,7 +211,7 @@ def prepare_ocp_second_pass( objective_functions.add(ObjectiveFcn.Mayer.MINIMIZE_TIME, weight=100 / n_shooting) # Dynamics - dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN) + dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver) # Path constraint x_bounds = BoundsList() @@ -267,7 +267,6 @@ def prepare_ocp_second_pass( u_bounds=u_bounds, objective_functions=objective_functions, constraints=constraints, - ode_solver=ode_solver, use_sx=use_sx, n_threads=n_threads, ) diff --git a/bioptim/examples/getting_started/example_cyclic_movement.py b/bioptim/examples/getting_started/example_cyclic_movement.py index a2ca1a7c8..4c6ee47c4 100644 --- a/bioptim/examples/getting_started/example_cyclic_movement.py +++ b/bioptim/examples/getting_started/example_cyclic_movement.py @@ -75,7 +75,9 @@ def prepare_ocp( objective_functions = Objective(ObjectiveFcn.Lagrange.MINIMIZE_CONTROL, key="tau", weight=100) # Dynamics - dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics = Dynamics( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) # Constraints constraints = ConstraintList() @@ -115,7 +117,6 @@ def prepare_ocp( u_bounds=u_bounds, objective_functions=objective_functions, constraints=constraints, - ode_solver=ode_solver, phase_transitions=phase_transitions, ) diff --git a/bioptim/examples/getting_started/example_external_forces.py b/bioptim/examples/getting_started/example_external_forces.py index 073fd01e9..cb82c516b 100644 --- a/bioptim/examples/getting_started/example_external_forces.py +++ b/bioptim/examples/getting_started/example_external_forces.py @@ -85,6 +85,7 @@ def prepare_ocp( dynamics = DynamicsList() dynamics.add( DynamicsFcn.TORQUE_DRIVEN, + ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics, numerical_data_timeseries=numerical_time_series, @@ -116,7 +117,6 @@ def prepare_ocp( u_bounds=u_bounds, objective_functions=objective_functions, constraints=constraints, - ode_solver=ode_solver, n_threads=n_threads, use_sx=use_sx, ) diff --git a/bioptim/examples/getting_started/example_inequality_constraint.py b/bioptim/examples/getting_started/example_inequality_constraint.py index 1d3698b8c..3c41339d9 100644 --- a/bioptim/examples/getting_started/example_inequality_constraint.py +++ b/bioptim/examples/getting_started/example_inequality_constraint.py @@ -96,6 +96,7 @@ def prepare_ocp( dynamics.add( DynamicsFcn.TORQUE_DRIVEN, contact_type=[ContactType.RIGID_EXPLICIT], + ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics, ) @@ -154,7 +155,6 @@ def prepare_ocp( objective_functions=objective_functions, constraints=constraints, variable_mappings=dof_mapping, - ode_solver=ode_solver, ) diff --git a/bioptim/examples/getting_started/example_joints_acceleration_driven.py b/bioptim/examples/getting_started/example_joints_acceleration_driven.py index 33afe1f24..b764f1cfa 100644 --- a/bioptim/examples/getting_started/example_joints_acceleration_driven.py +++ b/bioptim/examples/getting_started/example_joints_acceleration_driven.py @@ -62,7 +62,7 @@ def prepare_ocp( objective_functions = Objective(ObjectiveFcn.Lagrange.MINIMIZE_CONTROL, key="qddot_joints") # Dynamics - dynamics = Dynamics(DynamicsFcn.JOINTS_ACCELERATION_DRIVEN) + dynamics = Dynamics(DynamicsFcn.JOINTS_ACCELERATION_DRIVEN, ode_solver=ode_solver) # Path constraint x_bounds = BoundsList() @@ -100,7 +100,6 @@ def prepare_ocp( x_bounds=x_bounds, u_bounds=u_bounds, objective_functions=objective_functions, - ode_solver=ode_solver, use_sx=use_sx, n_threads=n_threads, ) diff --git a/bioptim/examples/getting_started/example_multinode_constraints.py b/bioptim/examples/getting_started/example_multinode_constraints.py index f3be926e4..cbf918fcc 100644 --- a/bioptim/examples/getting_started/example_multinode_constraints.py +++ b/bioptim/examples/getting_started/example_multinode_constraints.py @@ -116,9 +116,15 @@ def prepare_ocp( # Dynamics dynamics = DynamicsList() - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) # Constraints constraints = ConstraintList() @@ -197,7 +203,6 @@ def prepare_ocp( constraints=constraints, multinode_constraints=multinode_constraints, multinode_objectives=multinode_objectives, - ode_solver=ode_solver, ) diff --git a/bioptim/examples/getting_started/example_multinode_objective.py b/bioptim/examples/getting_started/example_multinode_objective.py index 2e1a15fde..5f5eccdd7 100644 --- a/bioptim/examples/getting_started/example_multinode_objective.py +++ b/bioptim/examples/getting_started/example_multinode_objective.py @@ -90,7 +90,9 @@ def prepare_ocp( ) # Dynamics - dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics = Dynamics( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) # Path constraint x_bounds = BoundsList() @@ -115,7 +117,6 @@ def prepare_ocp( x_bounds=x_bounds, u_bounds=u_bounds, multinode_objectives=multinode_objectives, - ode_solver=ode_solver, use_sx=use_sx, n_threads=n_threads, # This has to be set to 1 by definition. ) diff --git a/bioptim/examples/getting_started/example_multiphase.py b/bioptim/examples/getting_started/example_multiphase.py index cb841718e..d7f0f1dc6 100644 --- a/bioptim/examples/getting_started/example_multiphase.py +++ b/bioptim/examples/getting_started/example_multiphase.py @@ -119,10 +119,27 @@ def prepare_ocp( ) # Dynamics + if not isinstance(ode_solver, list): + ode_solver = [ode_solver] * 3 dynamics = DynamicsList() - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, + ode_solver=ode_solver[0], + expand_dynamics=expand_dynamics, + phase_dynamics=phase_dynamics, + ) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, + ode_solver=ode_solver[1], + expand_dynamics=expand_dynamics, + phase_dynamics=phase_dynamics, + ) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, + ode_solver=ode_solver[2], + expand_dynamics=expand_dynamics, + phase_dynamics=phase_dynamics, + ) # Constraints constraints = ConstraintList() @@ -162,7 +179,6 @@ def prepare_ocp( objective_functions=objective_functions, constraints=constraints, multinode_objectives=multinode_objective, - ode_solver=ode_solver, control_type=control_type, ) diff --git a/bioptim/examples/getting_started/example_parameter_scaling.py b/bioptim/examples/getting_started/example_parameter_scaling.py index 958092797..c4435df8f 100644 --- a/bioptim/examples/getting_started/example_parameter_scaling.py +++ b/bioptim/examples/getting_started/example_parameter_scaling.py @@ -80,7 +80,9 @@ def generate_dat_to_track( objective_functions.add(ObjectiveFcn.Lagrange.MINIMIZE_STATE, key="q", weight=1) # Dynamics - dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics = Dynamics( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) # Path constraint x_bounds = BoundsList() @@ -104,7 +106,6 @@ def generate_dat_to_track( x_bounds=x_bounds, u_bounds=u_bounds, objective_functions=objective_functions, - ode_solver=ode_solver, use_sx=use_sx, ) @@ -207,6 +208,7 @@ def prepare_ocp( dynamics = Dynamics( DynamicsFcn.TORQUE_DRIVEN, state_continuity_weight=100, + ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics, ) @@ -242,7 +244,6 @@ def prepare_ocp( parameter_objectives=parameter_objectives, parameter_bounds=parameter_bounds, parameter_init=parameter_init, - ode_solver=ode_solver, use_sx=use_sx, ) diff --git a/bioptim/examples/getting_started/example_variable_scaling.py b/bioptim/examples/getting_started/example_variable_scaling.py index 10a0ec85a..48a1fbb63 100644 --- a/bioptim/examples/getting_started/example_variable_scaling.py +++ b/bioptim/examples/getting_started/example_variable_scaling.py @@ -72,7 +72,9 @@ def prepare_ocp( objective_functions = Objective(ObjectiveFcn.Lagrange.MINIMIZE_CONTROL, key="tau") # Dynamics - dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics = Dynamics( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) # Path constraint x_bounds = BoundsList() @@ -107,7 +109,6 @@ def prepare_ocp( x_scaling=x_scaling, u_scaling=u_scaling, objective_functions=objective_functions, - ode_solver=ode_solver, use_sx=use_sx, n_threads=n_threads, ) diff --git a/bioptim/examples/getting_started/pendulum.py b/bioptim/examples/getting_started/pendulum.py index 2e619e267..c094e0f74 100644 --- a/bioptim/examples/getting_started/pendulum.py +++ b/bioptim/examples/getting_started/pendulum.py @@ -79,7 +79,9 @@ def prepare_ocp( objective_functions = Objective(ObjectiveFcn.Lagrange.MINIMIZE_CONTROL, key="tau") # Dynamics - dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics = Dynamics( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) # Path bounds x_bounds = BoundsList() @@ -114,7 +116,6 @@ def prepare_ocp( x_bounds=x_bounds, u_bounds=u_bounds, objective_functions=objective_functions, - ode_solver=ode_solver, use_sx=use_sx, n_threads=n_threads, control_type=control_type, diff --git a/bioptim/examples/getting_started/pendulum_constrained_states_controls.py b/bioptim/examples/getting_started/pendulum_constrained_states_controls.py index 28e053a7e..32255ee00 100644 --- a/bioptim/examples/getting_started/pendulum_constrained_states_controls.py +++ b/bioptim/examples/getting_started/pendulum_constrained_states_controls.py @@ -85,7 +85,9 @@ def prepare_ocp( objective_functions = Objective(ObjectiveFcn.Lagrange.MINIMIZE_CONTROL, key="tau") # Dynamics - dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics = Dynamics( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) # Path state constraints constraints = ConstraintList() @@ -123,7 +125,6 @@ def prepare_ocp( final_time, constraints=constraints, objective_functions=objective_functions, - ode_solver=ode_solver, use_sx=use_sx, n_threads=n_threads, control_type=control_type, diff --git a/bioptim/examples/holonomic_constraints/two_pendulums.py b/bioptim/examples/holonomic_constraints/two_pendulums.py index 43d858dd5..a3363aeeb 100644 --- a/bioptim/examples/holonomic_constraints/two_pendulums.py +++ b/bioptim/examples/holonomic_constraints/two_pendulums.py @@ -133,7 +133,7 @@ def prepare_ocp( # Dynamics dynamics = DynamicsList() - dynamics.add(DynamicsFcn.HOLONOMIC_TORQUE_DRIVEN, expand_dynamics=expand_dynamics) + dynamics.add(DynamicsFcn.HOLONOMIC_TORQUE_DRIVEN, ode_solver=OdeSolver.RK4(), expand_dynamics=expand_dynamics) # Path Constraints constraints = ConstraintList() @@ -176,7 +176,6 @@ def prepare_ocp( dynamics, n_shooting, final_time, - ode_solver=OdeSolver.RK4(), x_bounds=x_bounds, u_bounds=u_bounds, x_init=x_init, diff --git a/bioptim/examples/holonomic_constraints/two_pendulums_algebraic.py b/bioptim/examples/holonomic_constraints/two_pendulums_algebraic.py index db81f783b..5fbea5916 100644 --- a/bioptim/examples/holonomic_constraints/two_pendulums_algebraic.py +++ b/bioptim/examples/holonomic_constraints/two_pendulums_algebraic.py @@ -144,6 +144,7 @@ def prepare_ocp( # dynamics.add(DynamicsFcn.HOLONOMIC_TORQUE_DRIVEN, expand_dynamics=expand_dynamics) dynamics.add( configure_holonomic_torque_driven, + ode_solver=ode_solver, dynamic_function=holonomic_torque_driven_with_qv, expand_dynamics=expand_dynamics, ) @@ -219,7 +220,6 @@ def prepare_ocp( a_init=a_init, objective_functions=objective_functions, variable_mappings=tau_variable_bimapping, - ode_solver=ode_solver, constraints=constraints, ), bio_model, diff --git a/bioptim/examples/muscle_driven_ocp/muscle_activations_tracker.py b/bioptim/examples/muscle_driven_ocp/muscle_activations_tracker.py index 39f0cdcb8..c5c5518db 100644 --- a/bioptim/examples/muscle_driven_ocp/muscle_activations_tracker.py +++ b/bioptim/examples/muscle_driven_ocp/muscle_activations_tracker.py @@ -284,6 +284,7 @@ def prepare_ocp( dynamics = DynamicsList() dynamics.add( DynamicsFcn.MUSCLE_DRIVEN, + ode_solver=ode_solver, with_residual_torque=use_residual_torque, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics, @@ -315,7 +316,6 @@ def prepare_ocp( u_bounds=u_bounds, u_init=u_init, objective_functions=objective_functions, - ode_solver=ode_solver, n_threads=n_threads, ) @@ -369,7 +369,11 @@ def main(): markers[:, :, i] = bio_model.markers()(q[:, i]) plt.figure("Markers") - n_steps_ode = ocp.nlp[0].ode_solver.steps + 1 if ocp.nlp[0].ode_solver.is_direct_collocation else 1 + n_steps_ode = ( + ocp.nlp[0].dynamics_type.ode_solver.steps + 1 + if ocp.nlp[0].dynamics_type.ode_solver.is_direct_collocation + else 1 + ) for i in range(markers.shape[1]): plt.plot( np.linspace(0, final_time, n_shooting_points + 1), diff --git a/bioptim/examples/muscle_driven_ocp/muscle_excitations_tracker.py b/bioptim/examples/muscle_driven_ocp/muscle_excitations_tracker.py index d4e6473d6..2dcfa2585 100644 --- a/bioptim/examples/muscle_driven_ocp/muscle_excitations_tracker.py +++ b/bioptim/examples/muscle_driven_ocp/muscle_excitations_tracker.py @@ -287,6 +287,7 @@ def prepare_ocp( DynamicsFcn.MUSCLE_DRIVEN, with_excitations=True, with_residual_torque=use_residual_torque, + ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics, ) @@ -322,7 +323,6 @@ def prepare_ocp( x_bounds=x_bounds, u_bounds=u_bounds, objective_functions=objective_functions, - ode_solver=ode_solver, ) @@ -372,7 +372,11 @@ def main(): markers[:, :, i] = horzcat(*bio_model.markers()(q[:, i])) plt.figure("Markers") - n_steps_ode = ocp.nlp[0].ode_solver.steps + 1 if ocp.nlp[0].ode_solver.is_direct_collocation else 1 + n_steps_ode = ( + ocp.nlp[0].dynamics_type.ode_solver.steps + 1 + if ocp.nlp[0].dynamics_type.ode_solver.is_direct_collocation + else 1 + ) for i in range(markers.shape[1]): plt.plot(np.linspace(0, 2, n_shooting_points + 1), markers_ref[:, i, :].T, "k") plt.plot(np.linspace(0, 2, n_shooting_points * n_steps_ode + 1), markers[:, i, :].T, "r--") diff --git a/bioptim/examples/muscle_driven_ocp/static_arm.py b/bioptim/examples/muscle_driven_ocp/static_arm.py index 016368d31..9e8426243 100644 --- a/bioptim/examples/muscle_driven_ocp/static_arm.py +++ b/bioptim/examples/muscle_driven_ocp/static_arm.py @@ -86,6 +86,7 @@ def prepare_ocp( dynamics.add( DynamicsFcn.MUSCLE_DRIVEN, with_residual_torque=True, + ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics, ) @@ -122,7 +123,6 @@ def prepare_ocp( x_init=x_init, u_init=u_init, objective_functions=objective_functions, - ode_solver=ode_solver, control_type=control_type, n_threads=n_threads, ) diff --git a/bioptim/examples/muscle_driven_ocp/static_arm_with_contact.py b/bioptim/examples/muscle_driven_ocp/static_arm_with_contact.py index 41ecd8e06..926ea1927 100644 --- a/bioptim/examples/muscle_driven_ocp/static_arm_with_contact.py +++ b/bioptim/examples/muscle_driven_ocp/static_arm_with_contact.py @@ -66,7 +66,12 @@ def prepare_ocp( # Dynamics dynamics = DynamicsList() - dynamics.add(DynamicsFcn.MUSCLE_DRIVEN, with_residual_torque=True, contact_type=[ContactType.RIGID_EXPLICIT]) + dynamics.add( + DynamicsFcn.MUSCLE_DRIVEN, + with_residual_torque=True, + contact_type=[ContactType.RIGID_EXPLICIT], + ode_solver=ode_solver, + ) # raise RuntimeError("This example is broken, since contact dynamics with muscle is not implemented") # Path constraint @@ -101,7 +106,6 @@ def prepare_ocp( x_bounds, u_bounds, objective_functions, - ode_solver=ode_solver, ) diff --git a/bioptim/examples/muscle_driven_with_contact/contact_forces_inequality_constraint_muscle_excitations.py b/bioptim/examples/muscle_driven_with_contact/contact_forces_inequality_constraint_muscle_excitations.py index f10c7183b..c1dd2caf9 100644 --- a/bioptim/examples/muscle_driven_with_contact/contact_forces_inequality_constraint_muscle_excitations.py +++ b/bioptim/examples/muscle_driven_with_contact/contact_forces_inequality_constraint_muscle_excitations.py @@ -51,6 +51,7 @@ def prepare_ocp(biorbd_model_path, phase_time, n_shooting, min_bound, ode_solver with_excitations=True, with_residual_torque=True, contact_type=[ContactType.RIGID_EXPLICIT], + ode_solver=ode_solver, expand_dynamics=expand_dynamics, ) @@ -114,7 +115,6 @@ def prepare_ocp(biorbd_model_path, phase_time, n_shooting, min_bound, ode_solver objective_functions=objective_functions, constraints=constraints, variable_mappings=dof_mapping, - ode_solver=ode_solver, ) diff --git a/bioptim/examples/muscle_driven_with_contact/muscle_activations_contacts_tracker.py b/bioptim/examples/muscle_driven_with_contact/muscle_activations_contacts_tracker.py index b0ff1a99f..638ee36af 100644 --- a/bioptim/examples/muscle_driven_with_contact/muscle_activations_contacts_tracker.py +++ b/bioptim/examples/muscle_driven_with_contact/muscle_activations_contacts_tracker.py @@ -56,7 +56,12 @@ def prepare_ocp( # Dynamics dynamics = DynamicsList() - dynamics.add(DynamicsFcn.MUSCLE_DRIVEN, with_residual_torque=True, contact_type=[ContactType.RIGID_EXPLICIT]) + dynamics.add( + DynamicsFcn.MUSCLE_DRIVEN, + with_residual_torque=True, + contact_type=[ContactType.RIGID_EXPLICIT], + ode_solver=ode_solver, + ) # Path constraint q_at_first_node = [0, 0, -0.75, 0.75] @@ -92,7 +97,6 @@ def prepare_ocp( x_bounds=x_bounds, u_bounds=u_bounds, objective_functions=objective_functions, - ode_solver=ode_solver, ) diff --git a/bioptim/examples/optimal_time_ocp/multiphase_time_constraint.py b/bioptim/examples/optimal_time_ocp/multiphase_time_constraint.py index 90ff36c1e..a222dcf64 100644 --- a/bioptim/examples/optimal_time_ocp/multiphase_time_constraint.py +++ b/bioptim/examples/optimal_time_ocp/multiphase_time_constraint.py @@ -97,10 +97,28 @@ def prepare_ocp( # Dynamics dynamics = DynamicsList() - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, phase=0, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, + phase=0, + ode_solver=ode_solver, + expand_dynamics=expand_dynamics, + phase_dynamics=phase_dynamics, + ) if n_phases == 3: - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, phase=1, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, phase=2, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, + phase=1, + ode_solver=ode_solver, + expand_dynamics=expand_dynamics, + phase_dynamics=phase_dynamics, + ) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, + phase=2, + ode_solver=ode_solver, + expand_dynamics=expand_dynamics, + phase_dynamics=phase_dynamics, + ) # Constraints constraints = ConstraintList() @@ -160,7 +178,6 @@ def prepare_ocp( u_bounds=u_bounds, objective_functions=objective_functions, constraints=constraints, - ode_solver=ode_solver, time_phase_mapping=time_phase_mapping, ) diff --git a/bioptim/examples/optimal_time_ocp/pendulum_min_time_Mayer.py b/bioptim/examples/optimal_time_ocp/pendulum_min_time_Mayer.py index 7c3a8c547..4bfb39130 100644 --- a/bioptim/examples/optimal_time_ocp/pendulum_min_time_Mayer.py +++ b/bioptim/examples/optimal_time_ocp/pendulum_min_time_Mayer.py @@ -86,7 +86,9 @@ def prepare_ocp( # Dynamics dynamics = DynamicsList() - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) # Path constraint x_bounds = BoundsList() @@ -111,7 +113,6 @@ def prepare_ocp( x_bounds=x_bounds, u_bounds=u_bounds, objective_functions=objective_functions, - ode_solver=ode_solver, control_type=control_type, ) diff --git a/bioptim/examples/optimal_time_ocp/time_constraint.py b/bioptim/examples/optimal_time_ocp/time_constraint.py index 162597a98..0a1cb7a14 100644 --- a/bioptim/examples/optimal_time_ocp/time_constraint.py +++ b/bioptim/examples/optimal_time_ocp/time_constraint.py @@ -75,7 +75,9 @@ def prepare_ocp( objective_functions = Objective(ObjectiveFcn.Lagrange.MINIMIZE_CONTROL, key="tau") # Dynamics - dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics = Dynamics( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) # Constraints constraints = Constraint(ConstraintFcn.TIME_CONSTRAINT, node=Node.END, min_bound=time_min, max_bound=time_max) @@ -104,7 +106,6 @@ def prepare_ocp( u_bounds=u_bounds, objective_functions=objective_functions, constraints=constraints, - ode_solver=ode_solver, ) diff --git a/bioptim/examples/sqp_method/pendulum.py b/bioptim/examples/sqp_method/pendulum.py index f928b5602..555159545 100644 --- a/bioptim/examples/sqp_method/pendulum.py +++ b/bioptim/examples/sqp_method/pendulum.py @@ -68,7 +68,9 @@ def prepare_ocp( objective_functions = Objective(ObjectiveFcn.Lagrange.MINIMIZE_CONTROL, key="tau") # Dynamics - dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics = Dynamics( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) # Path constraint x_bounds = BoundsList() @@ -93,7 +95,6 @@ def prepare_ocp( x_bounds=x_bounds, u_bounds=u_bounds, objective_functions=objective_functions, - ode_solver=ode_solver, use_sx=use_sx, n_threads=n_threads, ) diff --git a/bioptim/examples/stochastic_optimal_control/obstacle_avoidance_direct_collocation.py b/bioptim/examples/stochastic_optimal_control/obstacle_avoidance_direct_collocation.py index ac9df164b..361b3117e 100644 --- a/bioptim/examples/stochastic_optimal_control/obstacle_avoidance_direct_collocation.py +++ b/bioptim/examples/stochastic_optimal_control/obstacle_avoidance_direct_collocation.py @@ -557,6 +557,11 @@ def prepare_socp( ) else: + ode_solver = OdeSolver.COLLOCATION( + polynomial_degree=socp_type.polynomial_degree, + method=socp_type.method, + duplicate_starting_point=True, + ) dynamics.add( configure_optimal_control_problem, dynamic_function=lambda time, states, controls, parameters, algebraic_states, numerical_timeseries, nlp, with_noise: bio_model.dynamics( @@ -568,14 +573,10 @@ def prepare_socp( nlp, with_noise=with_noise, ), + ode_solver=ode_solver, phase_dynamics=phase_dynamics, expand_dynamics=expand_dynamics, ) - ode_solver = OdeSolver.COLLOCATION( - polynomial_degree=socp_type.polynomial_degree, - method=socp_type.method, - duplicate_starting_point=True, - ) return OptimalControlProgram( bio_model, @@ -591,7 +592,6 @@ def prepare_socp( control_type=ControlType.CONSTANT, n_threads=6, phase_transitions=phase_transitions, - ode_solver=ode_solver, ) diff --git a/bioptim/examples/stochastic_optimal_control/rockit_matrix_lyapunov.py b/bioptim/examples/stochastic_optimal_control/rockit_matrix_lyapunov.py index e89a087f8..f5bca9daa 100644 --- a/bioptim/examples/stochastic_optimal_control/rockit_matrix_lyapunov.py +++ b/bioptim/examples/stochastic_optimal_control/rockit_matrix_lyapunov.py @@ -196,6 +196,8 @@ def prepare_socp( ) else: + ode_solver = OdeSolver.COLLOCATION(polynomial_degree=socp_type.polynomial_degree, method=socp_type.method) + dynamics.add( configure_optimal_control_problem, dynamic_function=lambda time, states, controls, parameters, algebraic_states, nlp, with_noise: bio_model.dynamics( @@ -208,10 +210,8 @@ def prepare_socp( ), phase_dynamics=phase_dynamics, expand_dynamics=expand_dynamics, + ode_solver=ode_solver, ) - - ode_solver = OdeSolver.COLLOCATION(polynomial_degree=socp_type.polynomial_degree, method=socp_type.method) - return OptimalControlProgram( bio_model, dynamics, @@ -221,7 +221,6 @@ def prepare_socp( u_bounds=u_bounds, objective_functions=objective_functions, constraints=constraints, - ode_solver=ode_solver, n_threads=6, ) diff --git a/bioptim/examples/symmetrical_torque_driven_ocp/symmetry_by_constraint.py b/bioptim/examples/symmetrical_torque_driven_ocp/symmetry_by_constraint.py index 02b9d4ee5..ca1ed7b38 100644 --- a/bioptim/examples/symmetrical_torque_driven_ocp/symmetry_by_constraint.py +++ b/bioptim/examples/symmetrical_torque_driven_ocp/symmetry_by_constraint.py @@ -76,7 +76,9 @@ def prepare_ocp( # Dynamics dynamics = DynamicsList() - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) # Constraints constraints = ConstraintList() @@ -109,7 +111,6 @@ def prepare_ocp( u_bounds=u_bounds, objective_functions=objective_functions, constraints=constraints, - ode_solver=ode_solver, ) diff --git a/bioptim/examples/symmetrical_torque_driven_ocp/symmetry_by_mapping.py b/bioptim/examples/symmetrical_torque_driven_ocp/symmetry_by_mapping.py index 34b7dbe33..194ef915a 100644 --- a/bioptim/examples/symmetrical_torque_driven_ocp/symmetry_by_mapping.py +++ b/bioptim/examples/symmetrical_torque_driven_ocp/symmetry_by_mapping.py @@ -113,7 +113,9 @@ def prepare_ocp( # Dynamics dynamics = DynamicsList() - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) # Constraints constraints = ConstraintList() @@ -145,7 +147,6 @@ def prepare_ocp( u_bounds=u_bounds, objective_functions=objective_functions, constraints=constraints, - ode_solver=ode_solver, variable_mappings=dof_mappings, ) diff --git a/bioptim/examples/torque_driven_ocp/example_pendulum_time_dependent.py b/bioptim/examples/torque_driven_ocp/example_pendulum_time_dependent.py index b4ad2b739..3a166ff81 100644 --- a/bioptim/examples/torque_driven_ocp/example_pendulum_time_dependent.py +++ b/bioptim/examples/torque_driven_ocp/example_pendulum_time_dependent.py @@ -157,7 +157,11 @@ def prepare_ocp( dynamics = DynamicsList() expand = not isinstance(ode_solver, OdeSolver.IRK) dynamics.add( - custom_configure, dynamic_function=time_dependent_dynamic, expand_dynamics=expand, phase_dynamics=phase_dynamics + custom_configure, + dynamic_function=time_dependent_dynamic, + ode_solver=ode_solver, + expand_dynamics=expand, + phase_dynamics=phase_dynamics, ) # Path constraint @@ -193,7 +197,6 @@ def prepare_ocp( x_bounds=x_bounds, u_bounds=u_bounds, objective_functions=objective_functions, - ode_solver=ode_solver, use_sx=use_sx, control_type=control_type, ) diff --git a/bioptim/examples/torque_driven_ocp/example_quaternions.py b/bioptim/examples/torque_driven_ocp/example_quaternions.py index 2978704da..b18235823 100644 --- a/bioptim/examples/torque_driven_ocp/example_quaternions.py +++ b/bioptim/examples/torque_driven_ocp/example_quaternions.py @@ -207,7 +207,10 @@ def prepare_ocp( # Dynamics dynamics = DynamicsList() dynamics.add( - DynamicsFcn.TORQUE_DRIVEN_FREE_FLOATING_BASE, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + DynamicsFcn.TORQUE_DRIVEN_FREE_FLOATING_BASE, + ode_solver=ode_solver, + expand_dynamics=expand_dynamics, + phase_dynamics=phase_dynamics, ) # Define control path constraint @@ -252,7 +255,6 @@ def prepare_ocp( x_init=x_init, objective_functions=objective_functions, constraints=constraints, - ode_solver=ode_solver, ) diff --git a/bioptim/examples/torque_driven_ocp/example_soft_contact.py b/bioptim/examples/torque_driven_ocp/example_soft_contact.py index 19410bd13..be275d533 100644 --- a/bioptim/examples/torque_driven_ocp/example_soft_contact.py +++ b/bioptim/examples/torque_driven_ocp/example_soft_contact.py @@ -54,6 +54,7 @@ def prepare_single_shooting( DynamicsFcn.TORQUE_DRIVEN, soft_contacts_dynamics=SoftContactDynamics.ODE, contact_type=[ContactType.SOFT_EXPLICIT], + ode_solver=ode_solver, ) return OptimalControlProgram( @@ -61,7 +62,6 @@ def prepare_single_shooting( dynamics, n_shooting, final_time, - ode_solver=ode_solver, use_sx=use_sx, n_threads=n_threads, ) @@ -153,6 +153,7 @@ def prepare_ocp( DynamicsFcn.TORQUE_DRIVEN, soft_contacts_dynamics=SoftContactDynamics.ODE, contact_type=[ContactType.SOFT_EXPLICIT], + ode_solver=ode_solver, phase_dynamics=phase_dynamics, ) @@ -194,7 +195,6 @@ def prepare_ocp( x_init=x_init, objective_functions=objective_functions, constraints=constraints, - ode_solver=ode_solver, use_sx=use_sx, n_threads=n_threads, ) diff --git a/bioptim/examples/torque_driven_ocp/maximize_predicted_height_CoM.py b/bioptim/examples/torque_driven_ocp/maximize_predicted_height_CoM.py index 8a40438d4..02622d433 100644 --- a/bioptim/examples/torque_driven_ocp/maximize_predicted_height_CoM.py +++ b/bioptim/examples/torque_driven_ocp/maximize_predicted_height_CoM.py @@ -98,6 +98,7 @@ def prepare_ocp( dynamics.add( DynamicsFcn.TORQUE_ACTIVATIONS_DRIVEN, contact_type=[ContactType.RIGID_EXPLICIT], + ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics, ) @@ -105,6 +106,7 @@ def prepare_ocp( dynamics.add( DynamicsFcn.TORQUE_DRIVEN, contact_type=[ContactType.RIGID_EXPLICIT], + ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics, ) @@ -156,7 +158,6 @@ def prepare_ocp( objective_functions=objective_functions, constraints=constraints, variable_mappings=dof_mapping, - ode_solver=ode_solver, ) diff --git a/bioptim/examples/torque_driven_ocp/ocp_mass_with_ligament.py b/bioptim/examples/torque_driven_ocp/ocp_mass_with_ligament.py index 70cf69f47..b11dd9936 100644 --- a/bioptim/examples/torque_driven_ocp/ocp_mass_with_ligament.py +++ b/bioptim/examples/torque_driven_ocp/ocp_mass_with_ligament.py @@ -70,6 +70,7 @@ def prepare_ocp( dynamics.add( DynamicsFcn.TORQUE_DRIVEN, with_ligament=True, + ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics, ) @@ -97,7 +98,6 @@ def prepare_ocp( u_bounds=u_bounds, u_init=u_init, objective_functions=objective_functions, - ode_solver=ode_solver, n_threads=n_threads, use_sx=use_sx, ) diff --git a/bioptim/examples/torque_driven_ocp/pendulum_with_passive_torque.py b/bioptim/examples/torque_driven_ocp/pendulum_with_passive_torque.py index b7245ee3d..fd1887534 100644 --- a/bioptim/examples/torque_driven_ocp/pendulum_with_passive_torque.py +++ b/bioptim/examples/torque_driven_ocp/pendulum_with_passive_torque.py @@ -70,6 +70,7 @@ def prepare_ocp( dynamics = Dynamics( DynamicsFcn.TORQUE_DRIVEN, with_passive_torque=with_passive_torque, + ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics, ) @@ -98,7 +99,6 @@ def prepare_ocp( x_bounds=x_bounds, u_bounds=u_bounds, objective_functions=objective_functions, - ode_solver=ode_solver, ) diff --git a/bioptim/examples/torque_driven_ocp/slider.py b/bioptim/examples/torque_driven_ocp/slider.py index 02935836f..e919a3806 100644 --- a/bioptim/examples/torque_driven_ocp/slider.py +++ b/bioptim/examples/torque_driven_ocp/slider.py @@ -76,9 +76,15 @@ def prepare_ocp( # Dynamics dynamics = DynamicsList() - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) # Constraints constraints = ConstraintList() @@ -114,7 +120,6 @@ def prepare_ocp( u_bounds=u_bounds, objective_functions=objective_functions, constraints=constraints, - ode_solver=ode_solver, control_type=control_type, ) diff --git a/bioptim/examples/torque_driven_ocp/torque_activation_driven.py b/bioptim/examples/torque_driven_ocp/torque_activation_driven.py index 1b726eeb6..579a7527b 100644 --- a/bioptim/examples/torque_driven_ocp/torque_activation_driven.py +++ b/bioptim/examples/torque_driven_ocp/torque_activation_driven.py @@ -70,6 +70,7 @@ def prepare_ocp( dynamics.add( DynamicsFcn.TORQUE_ACTIVATIONS_DRIVEN, with_residual_torque=True, + ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics, ) @@ -108,7 +109,6 @@ def prepare_ocp( u_bounds=u_bounds, x_init=x_init, objective_functions=objective_functions, - ode_solver=ode_solver, ) diff --git a/bioptim/examples/torque_driven_ocp/track_markers_2D_pendulum.py b/bioptim/examples/torque_driven_ocp/track_markers_2D_pendulum.py index bb4853d87..31f863204 100644 --- a/bioptim/examples/torque_driven_ocp/track_markers_2D_pendulum.py +++ b/bioptim/examples/torque_driven_ocp/track_markers_2D_pendulum.py @@ -119,7 +119,9 @@ def prepare_ocp( # Dynamics dynamics = DynamicsList() - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) # Path constraint x_bounds = BoundsList() @@ -144,7 +146,6 @@ def prepare_ocp( x_bounds=x_bounds, u_bounds=u_bounds, objective_functions=objective_functions, - ode_solver=ode_solver, ) diff --git a/bioptim/examples/torque_driven_ocp/track_markers_with_torque_actuators.py b/bioptim/examples/torque_driven_ocp/track_markers_with_torque_actuators.py index 71fb3b2d0..e833c0e34 100644 --- a/bioptim/examples/torque_driven_ocp/track_markers_with_torque_actuators.py +++ b/bioptim/examples/torque_driven_ocp/track_markers_with_torque_actuators.py @@ -85,14 +85,27 @@ def prepare_ocp( if actuator_type: if actuator_type == 1: dynamics.add( - DynamicsFcn.TORQUE_ACTIVATIONS_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + DynamicsFcn.TORQUE_ACTIVATIONS_DRIVEN, + ode_solver=ode_solver, + expand_dynamics=expand_dynamics, + phase_dynamics=phase_dynamics, ) elif actuator_type == 2: - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, + ode_solver=ode_solver, + expand_dynamics=expand_dynamics, + phase_dynamics=phase_dynamics, + ) else: raise ValueError("actuator_type is 1 (torque activations) or 2 (torque max constraints)") else: - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, + ode_solver=ode_solver, + expand_dynamics=expand_dynamics, + phase_dynamics=phase_dynamics, + ) # Constraints constraints = ConstraintList() @@ -123,7 +136,6 @@ def prepare_ocp( u_bounds=u_bounds, objective_functions=objective_functions, constraints=constraints, - ode_solver=ode_solver, ) diff --git a/bioptim/examples/torque_driven_ocp/trampo_quaternions.py b/bioptim/examples/torque_driven_ocp/trampo_quaternions.py index 75a6450d4..6cc4490fe 100644 --- a/bioptim/examples/torque_driven_ocp/trampo_quaternions.py +++ b/bioptim/examples/torque_driven_ocp/trampo_quaternions.py @@ -90,7 +90,9 @@ def prepare_ocp( # Dynamics dynamics = DynamicsList() - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) # Define control path constraint n_tau = bio_model.nb_tau # bio_model.nb_tau @@ -134,7 +136,6 @@ def prepare_ocp( u_bounds=u_bounds, x_init=x_init, objective_functions=objective_functions, - ode_solver=ode_solver, ) diff --git a/bioptim/examples/track/optimal_estimation.py b/bioptim/examples/track/optimal_estimation.py index 88e6973f7..96672943c 100644 --- a/bioptim/examples/track/optimal_estimation.py +++ b/bioptim/examples/track/optimal_estimation.py @@ -105,7 +105,7 @@ def prepare_ocp_to_track( # Dynamics dynamics = DynamicsList() - dynamics.add(DynamicsFcn.TORQUE_DRIVEN) + dynamics.add(DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver) # Path constraint x_bounds = BoundsList() @@ -128,7 +128,6 @@ def prepare_ocp_to_track( u_bounds=u_bounds, objective_functions=objective_functions, variable_mappings=variable_mappings, - ode_solver=ode_solver, ) @@ -170,7 +169,7 @@ def prepare_optimal_estimation( # Dynamics dynamics = DynamicsList() - dynamics.add(DynamicsFcn.TORQUE_DRIVEN) + dynamics.add(DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver) # Path constraint x_bounds = BoundsList() @@ -191,7 +190,6 @@ def prepare_optimal_estimation( x_bounds=x_bounds, u_bounds=u_bounds, objective_functions=objective_functions, - ode_solver=ode_solver, ) diff --git a/bioptim/examples/track/track_marker_on_segment.py b/bioptim/examples/track/track_marker_on_segment.py index fe03ee1cf..36a073c0b 100644 --- a/bioptim/examples/track/track_marker_on_segment.py +++ b/bioptim/examples/track/track_marker_on_segment.py @@ -79,7 +79,9 @@ def prepare_ocp( # Dynamics dynamics = DynamicsList() - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) # Constraints if constr: @@ -126,7 +128,6 @@ def prepare_ocp( x_init=x_init, objective_functions=objective_functions, constraints=constraints, - ode_solver=ode_solver, use_sx=use_sx, ) diff --git a/bioptim/examples/track/track_segment_on_rt.py b/bioptim/examples/track/track_segment_on_rt.py index 2dfad0559..db7fcb641 100644 --- a/bioptim/examples/track/track_segment_on_rt.py +++ b/bioptim/examples/track/track_segment_on_rt.py @@ -69,7 +69,9 @@ def prepare_ocp( # Dynamics dynamics = DynamicsList() - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + ) # Constraints constraints = ConstraintList() @@ -98,7 +100,6 @@ def prepare_ocp( u_bounds=u_bounds, objective_functions=objective_functions, constraints=constraints, - ode_solver=ode_solver, ) diff --git a/bioptim/gui/graph.py b/bioptim/gui/graph.py index b771e8406..de0d4feff 100644 --- a/bioptim/gui/graph.py +++ b/bioptim/gui/graph.py @@ -325,7 +325,7 @@ def print(self): print(f"PHASE DURATION IS OPTIMIZED") print(f"SHOOTING NODES : {self.ocp.nlp[phase_idx].ns}") print(f"DYNAMICS: {self.ocp.nlp[phase_idx].dynamics_type.type.name}") - print(f"ODE: {self.ocp.nlp[phase_idx].ode_solver.integrator.__name__}") + print(f"ODE: {self.ocp.nlp[phase_idx].dynamics_type.ode_solver.integrator.__name__}") print(f"**********") print("") @@ -579,7 +579,7 @@ def _draw_nlp_node(self, g, phase_idx: Int): node_str += f"Phase duration: optimized
" node_str += f"Shooting nodes: {self.ocp.nlp[phase_idx].ns}
" node_str += f"Dynamics: {self.ocp.nlp[phase_idx].dynamics_type.type.name}
" - node_str += f"ODE: {self.ocp.nlp[phase_idx].ode_solver.integrator.__name__}
" + node_str += f"ODE: {self.ocp.nlp[phase_idx].dynamics_type.ode_solver.integrator.__name__}
" node_str += f"Control type: {self.ocp.nlp[phase_idx].control_type.name}" g.node(f"nlp_node_{phase_idx}", f"""<{node_str}>""") diff --git a/bioptim/gui/plot.py b/bioptim/gui/plot.py index cbf0c6c55..6614d0bd4 100644 --- a/bioptim/gui/plot.py +++ b/bioptim/gui/plot.py @@ -514,8 +514,8 @@ def _set_bounds_for_axis(self, nlp, variable, ctr, mapping_to_first_index, y_min def _get_custom_bounds(self, nlp, variable, ctr, mapping_to_first_index): """Get custom bounds for a variable""" repeat = 1 - if isinstance(nlp.ode_solver, OdeSolver.COLLOCATION): - repeat = nlp.ode_solver.polynomial_degree + 1 + if isinstance(nlp.dynamics_type.ode_solver, OdeSolver.COLLOCATION): + repeat = nlp.dynamics_type.ode_solver.polynomial_degree + 1 nlp.plot[variable].bounds.check_and_adjust_dimensions(len(mapping_to_first_index), nlp.ns) idx = mapping_to_first_index.index(ctr) diff --git a/bioptim/gui/serializable_class.py b/bioptim/gui/serializable_class.py index 7d4879b2d..d5f17c293 100644 --- a/bioptim/gui/serializable_class.py +++ b/bioptim/gui/serializable_class.py @@ -472,7 +472,6 @@ def __init__( self.algebraic_states = algebraic_states self.parameters = parameters self.numerical_timeseries = numerical_timeseries - self.ode_solver = ode_solver self.plot = plot @classmethod @@ -489,7 +488,6 @@ def from_nlp(cls, nlp): algebraic_states=OptimizationVariableContainerSerializable.from_container(nlp.algebraic_states), parameters=OptimizationVariableContainerSerializable.from_container(nlp.parameters), numerical_timeseries=OptimizationVariableContainerSerializable.from_container(nlp.numerical_timeseries), - ode_solver=OdeSolverSerializable.from_ode_solver(nlp.ode_solver), plot={key: CustomPlotSerializable.from_custom_plot(nlp.plot[key]) for key in nlp.plot}, ) @@ -504,7 +502,6 @@ def serialize(self): "algebraic_states": self.algebraic_states.serialize(), "parameters": self.parameters.serialize(), "numerical_timeseries": self.numerical_timeseries.serialize(), - "ode_solver": self.ode_solver.serialize(), "plot": {key: plot.serialize() for key, plot in self.plot.items()}, } @@ -520,7 +517,6 @@ def deserialize(cls, data: AnyDict): algebraic_states=OptimizationVariableContainerSerializable.deserialize(data["algebraic_states"]), parameters=OptimizationVariableContainerSerializable.deserialize(data["parameters"]), numerical_timeseries=OptimizationVariableContainerSerializable.deserialize(data["numerical_timeseries"]), - ode_solver=OdeSolverSerializable.deserialize(data["ode_solver"]), plot={key: CustomPlotSerializable.deserialize(plot) for key, plot in data["plot"].items()}, ) diff --git a/bioptim/limits/penalty.py b/bioptim/limits/penalty.py index f1cb6884f..cd6db7a96 100644 --- a/bioptim/limits/penalty.py +++ b/bioptim/limits/penalty.py @@ -934,7 +934,7 @@ def minimize_explicit_rigid_contact_forces_end_of_interval( t_span = controller.t_span.cx cx = ( horzcat(*([controller.states.cx_start] + controller.states.cx_intermediates_list)) - if controller.get_nlp.ode_solver.is_direct_collocation + if controller.get_nlp.dynamics_type.ode_solver.is_direct_collocation else controller.states.cx_start ) end_of_interval_states = controller.integrate( @@ -1200,7 +1200,7 @@ def state_continuity(penalty: PenaltyOption, controller: PenaltyController | lis t_span = controller.t_span.cx continuity = controller.states.cx_end - if controller.get_nlp.ode_solver.is_direct_collocation: + if controller.get_nlp.dynamics_type.ode_solver.is_direct_collocation: states_cx = horzcat(*([controller.states.cx_start] + controller.states.cx_intermediates_list)) algebraic_states_cx = horzcat( *([controller.algebraic_states.cx_start] + controller.algebraic_states.cx_intermediates_list) diff --git a/bioptim/limits/penalty_controller.py b/bioptim/limits/penalty_controller.py index c6ff8607c..e697a8d7e 100644 --- a/bioptim/limits/penalty_controller.py +++ b/bioptim/limits/penalty_controller.py @@ -102,7 +102,7 @@ def control_type(self) -> ControlType: @property def ode_solver(self) -> OdeSolver: - return self._nlp.ode_solver + return self._nlp.dynamics_type.ode_solver @property def phase_idx(self) -> int: diff --git a/bioptim/limits/penalty_option.py b/bioptim/limits/penalty_option.py index 7665214a2..a188db343 100644 --- a/bioptim/limits/penalty_option.py +++ b/bioptim/limits/penalty_option.py @@ -366,7 +366,7 @@ def _set_control_types(self, controllers: list[PenaltyController]): self.control_types = control_types def _set_subnodes_are_decision_states(self, controllers: list[PenaltyController]): - subnodes_are_decision_states = [c.get_nlp.ode_solver.is_direct_collocation for c in controllers] + subnodes_are_decision_states = [c.get_nlp.dynamics_type.ode_solver.is_direct_collocation for c in controllers] if self.subnodes_are_decision_states: # If it was already set (e.g. for multinode), we want to make sure it is consistent if self.subnodes_are_decision_states != subnodes_are_decision_states: @@ -613,7 +613,7 @@ def _check_sanity_of_penalty_interactions(self, controller): if self.derivative and self.explicit_derivative: raise ValueError("derivative and explicit_derivative cannot be true simultaneously") - if controller.get_nlp.ode_solver.is_direct_collocation and ( + if controller.get_nlp.dynamics_type.ode_solver.is_direct_collocation and ( controller.get_nlp.phase_dynamics == PhaseDynamics.SHARED_DURING_THE_PHASE and len(self.node_idx) > 1 and controller.ns + 1 in self.node_idx diff --git a/bioptim/optimization/non_linear_program.py b/bioptim/optimization/non_linear_program.py index aaf20f156..d2a269fff 100644 --- a/bioptim/optimization/non_linear_program.py +++ b/bioptim/optimization/non_linear_program.py @@ -61,8 +61,6 @@ class NonLinearProgram: The number of thread to use ns: int The number of shooting points - ode_solver: OdeSolverBase - The chosen ode solver parameters: ParameterContainer Reference to the optimized parameters in the underlying ocp par_dynamics: casadi.Function @@ -146,7 +144,6 @@ def __init__(self, phase_dynamics: PhaseDynamics, use_sx: bool): self.model: BioModel | StochasticBioModel | HolonomicBioModel | VariationalBioModel | None = None self.n_threads = None self.ns = None - self.ode_solver = OdeSolver.RK4() self.par_dynamics = None self.phase_idx = None self.phase_mapping = None @@ -209,10 +206,14 @@ def initialize(self, cx: MX | SX | Callable = None): self.algebraic_states.initialize_from_shooting(n_shooting=self.ns + 1, cx=self.cx) self.integrated_values.initialize_from_shooting(n_shooting=self.ns + 1, cx=self.cx) - self.states.initialize_intermediate_cx(n_shooting=self.ns + 1, n_cx=self.ode_solver.n_required_cx) - self.states_dot.initialize_intermediate_cx(n_shooting=self.ns + 1, n_cx=self.ode_solver.n_required_cx) + self.states.initialize_intermediate_cx(n_shooting=self.ns + 1, n_cx=self.dynamics_type.ode_solver.n_required_cx) + self.states_dot.initialize_intermediate_cx( + n_shooting=self.ns + 1, n_cx=self.dynamics_type.ode_solver.n_required_cx + ) self.controls.initialize_intermediate_cx(n_shooting=self.ns + 1, n_cx=1) - self.algebraic_states.initialize_intermediate_cx(n_shooting=self.ns + 1, n_cx=self.ode_solver.n_required_cx) + self.algebraic_states.initialize_intermediate_cx( + n_shooting=self.ns + 1, n_cx=self.dynamics_type.ode_solver.n_required_cx + ) self.integrated_values.initialize_intermediate_cx(n_shooting=self.ns + 1, n_cx=1) def update_bounds(self, x_bounds, u_bounds, a_bounds): @@ -280,13 +281,13 @@ def _update_plot_bounds(self, keys, bounds, suffix): def update_init(self, x_init, u_init, a_init): - if x_init is not None: - not_direct_collocation = not self.ode_solver.is_direct_collocation - init_all_point = x_init.type == InterpolationType.ALL_POINTS + if x_init is not None or a_init is not None: + not_direct_collocation = not self.dynamics_type.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 init_all_point: + 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") - # TODO @ipuch in PR #907, add algebraic states to the error message self._update_init( init=x_init, @@ -417,23 +418,27 @@ def n_states_decision_steps(self, node_idx) -> int: """ if node_idx >= self.ns: return 1 - return self.dynamics[node_idx].shape_xf[1] + (1 if self.ode_solver.duplicate_starting_point else 0) + return self.dynamics[node_idx].shape_xf[1] + ( + 1 if self.dynamics_type.ode_solver.duplicate_starting_point else 0 + ) - def n_states_stepwise_steps(self, node_idx) -> int: + def n_states_stepwise_steps(self, node_idx: int, ode_solver: OdeSolver = None) -> int: """ Parameters ---------- node_idx: int The index of the node - + ode_solver: OdeSolver + The ode solver to use (it is useful for reintegration of COLLOCATION solutions) Returns ------- The number of states """ + ode_solver = ode_solver if ode_solver is not None else self.dynamics_type.ode_solver 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) + if ode_solver.is_direct_collocation: + return self.dynamics[node_idx].shape_xall[1] - (1 if not ode_solver.duplicate_starting_point else 0) else: return self.dynamics[node_idx].shape_xall[1] diff --git a/bioptim/optimization/optimal_control_program.py b/bioptim/optimization/optimal_control_program.py index 5dd6a270a..a39b9f8a5 100644 --- a/bioptim/optimization/optimal_control_program.py +++ b/bioptim/optimization/optimal_control_program.py @@ -10,7 +10,7 @@ from .non_linear_program import NonLinearProgram as NLP from .optimization_vector import OptimizationVectorHelper from ..dynamics.configure_problem import DynamicsList, Dynamics, ConfigureProblem -from ..dynamics.ode_solvers import OdeSolver, OdeSolverBase +from ..dynamics.ode_solvers import OdeSolver from ..gui.check_conditioning import check_conditioning from ..gui.graph import OcpToConsole, OcpToGraph from ..gui.ipopt_output_plot import SaveIterationsInfo @@ -155,7 +155,6 @@ def __init__( parameter_init: InitialGuessList = None, parameter_objectives: ParameterObjectiveList = None, parameter_constraints: ParameterConstraintList = None, - ode_solver: list | OdeSolverBase | OdeSolver = None, control_type: ControlType | list = ControlType.CONSTANT, variable_mappings: BiMappingList = None, time_phase_mapping: BiMapping = None, @@ -213,8 +212,6 @@ def __init__( All the parameter objectives to optimize of the program parameter_constraints: ParameterConstraintList All the parameter constraints of the program - ode_solver: OdeSolverBase - The solver for the ordinary differential equations control_type: ControlType The type of controls for each phase variable_mappings: BiMappingList @@ -284,7 +281,6 @@ def __init__( parameter_init, parameter_constraints, parameter_objectives, - ode_solver, use_sx, bio_model, plot_mappings, @@ -423,7 +419,6 @@ def _check_arguments_and_build_nlp( parameter_init, parameter_constraints, parameter_objectives, - ode_solver, use_sx, bio_model, plot_mappings, @@ -500,18 +495,6 @@ def _check_arguments_and_build_nlp( elif not isinstance(parameter_constraints, ParameterConstraintList): raise RuntimeError("constraints should be built from an Constraint or ConstraintList") - 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") - if not isinstance(use_sx, bool): raise RuntimeError("use_sx should be a bool") @@ -577,7 +560,6 @@ def _check_arguments_and_build_nlp( # Prepare path constraints and dynamics of the program NLP.add(self, "dynamics_type", dynamics, False) - NLP.add(self, "ode_solver", ode_solver, True) NLP.add(self, "control_type", control_type, True) # Prepare the variable mappings @@ -614,7 +596,7 @@ def _prepare_dynamics(self): self.nlp[i].parameters = self.parameters # This should be remove when phase parameters will be implemented self.nlp[i].numerical_data_timeseries = self.nlp[i].dynamics_type.numerical_data_timeseries ConfigureProblem.initialize(self, self.nlp[i]) - self.nlp[i].ode_solver.prepare_dynamic_integrator(self, self.nlp[i]) + self.nlp[i].dynamics_type.ode_solver.prepare_dynamic_integrator(self, self.nlp[i]) if (isinstance(self.nlp[i].model, VariationalBiorbdModel)) and self.nlp[i].algebraic_states.shape > 0: raise NotImplementedError( "Algebraic states were not tested with variational integrators. If you come across this error, " @@ -860,7 +842,10 @@ def _declare_inner_phase_continuity(self, nlp: NLP) -> None: ConstraintFcn.STATE_CONTINUITY, node=Node.ALL_SHOOTING, penalty_type=PenaltyType.INTERNAL ) penalty.add_or_replace_to_penalty_pool(self, nlp) - if nlp.ode_solver.is_direct_collocation and nlp.ode_solver.duplicate_starting_point: + if ( + nlp.dynamics_type.ode_solver.is_direct_collocation + and nlp.dynamics_type.ode_solver.duplicate_starting_point + ): penalty = Constraint( ConstraintFcn.FIRST_COLLOCATION_HELPER_EQUALS_STATE, node=Node.ALL_SHOOTING, @@ -1433,7 +1418,7 @@ def set_warm_start(self, sol: Solution): for i in range(self.n_phases): x_interp = ( InterpolationType.EACH_FRAME - if self.nlp[i].ode_solver.is_direct_shooting + if self.nlp[i].dynamics_type.ode_solver.is_direct_shooting else InterpolationType.ALL_POINTS ) if self.n_phases == 1: @@ -1680,16 +1665,10 @@ def node_time(self, phase_idx: int, node_idx: int): return previous_phase_time + self.nlp[phase_idx].dt * node_idx - def _set_default_ode_solver(self): - """ - Set the default ode solver to RK4 - """ - return OdeSolver.RK4() - def _set_nlp_is_stochastic(self): """ Set the is_stochastic variable to False - because it's not relevant for traditional OCP, + because it's not relevant for traditional OCP,_ only relevant for StochasticOptimalControlProgram Note diff --git a/bioptim/optimization/receding_horizon_optimization.py b/bioptim/optimization/receding_horizon_optimization.py index ee980920e..c9a61992a 100644 --- a/bioptim/optimization/receding_horizon_optimization.py +++ b/bioptim/optimization/receding_horizon_optimization.py @@ -867,7 +867,6 @@ def _initialize_one_cycle(self, dt: float, states: np.ndarray, controls: np.ndar bio_model=model_class(**model_initializer), dynamics=self.nlp[0].dynamics_type, objective_functions=deepcopy(self.common_objective_functions), - ode_solver=self.nlp[0].ode_solver, n_shooting=self.cycle_len, phase_time=self.cycle_len * dt, x_init=x_init, diff --git a/bioptim/optimization/solution/solution.py b/bioptim/optimization/solution/solution.py index d46e18f61..290cdda97 100644 --- a/bioptim/optimization/solution/solution.py +++ b/bioptim/optimization/solution/solution.py @@ -266,8 +266,8 @@ def from_initial_guess(cls, ocp, sol: list): # For states for p, ss in enumerate(sol_states): repeat = 1 - if isinstance(ocp.nlp[p].ode_solver, OdeSolver.COLLOCATION): - repeat = ocp.nlp[p].ode_solver.polynomial_degree + 1 + if isinstance(ocp.nlp[p].dynamics_type.ode_solver, OdeSolver.COLLOCATION): + repeat = ocp.nlp[p].dynamics_type.ode_solver.polynomial_degree + 1 for key in ss.keys(): ns = ocp.nlp[p].ns + 1 if ss[key].init.type != InterpolationType.EACH_FRAME else ocp.nlp[p].ns ss[key].init.check_and_adjust_dimensions(len(ocp.nlp[p].states[key]), ns, "states") @@ -455,8 +455,8 @@ def _process_time_vector( times.append([t[[0, -1]] for t in times_tp[nlp.phase_idx][:-1]]) else: if time_alignment == TimeAlignment.STATES: - if nlp.ode_solver.is_direct_collocation: - if nlp.ode_solver.duplicate_starting_point: + if nlp.dynamics_type.ode_solver.is_direct_collocation: + if nlp.dynamics_type.ode_solver.duplicate_starting_point: times.append( [t if t.shape == (1, 1) else vertcat(t[0], t[:-1]) for t in times_tp[nlp.phase_idx]] ) @@ -722,7 +722,7 @@ def _prepare_integrate(self, integrator: SolutionIntegrator): The integrator to use for the integration """ - has_direct_collocation = sum([nlp.ode_solver.is_direct_collocation for nlp in self.ocp.nlp]) > 0 + has_direct_collocation = sum([nlp.dynamics_type.ode_solver.is_direct_collocation for nlp in self.ocp.nlp]) > 0 if has_direct_collocation and integrator == SolutionIntegrator.OCP: raise ValueError( "When the ode_solver of the Optimal Control Problem is OdeSolver.COLLOCATION, " @@ -731,7 +731,9 @@ def _prepare_integrate(self, integrator: SolutionIntegrator): " Shooting.SINGLE, Shooting.MULTIPLE, or Shooting.SINGLE_DISCONTINUOUS_PHASE" ) - has_trapezoidal = sum([isinstance(nlp.ode_solver, OdeSolver.TRAPEZOIDAL) for nlp in self.ocp.nlp]) > 0 + has_trapezoidal = ( + sum([isinstance(nlp.dynamics_type.ode_solver, OdeSolver.TRAPEZOIDAL) for nlp in self.ocp.nlp]) > 0 + ) if has_trapezoidal and integrator == SolutionIntegrator.OCP: raise ValueError( "When the ode_solver of the Optimal Control Problem is OdeSolver.TRAPEZOIDAL, " diff --git a/bioptim/optimization/stochastic_optimal_control_program.py b/bioptim/optimization/stochastic_optimal_control_program.py index 9b4c75dec..0967cb6db 100644 --- a/bioptim/optimization/stochastic_optimal_control_program.py +++ b/bioptim/optimization/stochastic_optimal_control_program.py @@ -73,7 +73,6 @@ def __init__( **kwargs, ): _check_multi_threading_and_problem_type(problem_type, **kwargs) - _check_has_no_ode_solver_defined(**kwargs) _check_has_no_phase_dynamics_shared_during_the_phase(problem_type, **kwargs) self.problem_type = problem_type @@ -86,6 +85,10 @@ def __init__( if parameter_init is None: parameter_init = InitialGuessList() + # Integrator + for dyn in dynamics: + dyn.ode_solver = self._set_default_ode_solver() + if "motor_noise" not in parameters.keys(): n_motor_noise = bio_model.motor_noise_magnitude.shape[0] parameters.add( @@ -144,7 +147,6 @@ def __init__( parameter_init=parameter_init, parameter_objectives=parameter_objectives, parameter_constraints=parameter_constraints, - ode_solver=None, control_type=control_type, variable_mappings=variable_mappings, time_phase_mapping=time_phase_mapping, @@ -619,22 +621,6 @@ def _check_multi_threading_and_problem_type(problem_type, **kwargs): ) -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" - ")" - ) - - def _check_has_no_phase_dynamics_shared_during_the_phase(problem_type, **kwargs): if not isinstance(problem_type, SocpType.COLLOCATION): if "phase_dynamics" in kwargs: diff --git a/bioptim/optimization/variational_optimal_control_program.py b/bioptim/optimization/variational_optimal_control_program.py index a65d9b24a..e8000937b 100644 --- a/bioptim/optimization/variational_optimal_control_program.py +++ b/bioptim/optimization/variational_optimal_control_program.py @@ -8,6 +8,7 @@ from .optimal_control_program import OptimalControlProgram from ..dynamics.configure_problem import ConfigureProblem, DynamicsList from ..dynamics.dynamics_evaluation import DynamicsEvaluation +from ..dynamics.ode_solvers import OdeSolver from ..limits.constraints import ParameterConstraintList from ..limits.multinode_constraint import MultinodeConstraintList from ..limits.objective_functions import ParameterObjectiveList @@ -67,12 +68,6 @@ def __init__( " define it." ) - if "ode_solver" in kwargs: - raise ValueError( - "ode_solver cannot be defined in VariationalOptimalControlProgram since the integration is" - " done by the variational integrator." - ) - if "x_init" in kwargs or "x_bounds" in kwargs: raise ValueError( "In VariationalOptimalControlProgram q_init and q_bounds must be used instead of x_init and x_bounds " @@ -89,6 +84,7 @@ def __init__( self.configure_torque_driven, expand_dynamics=expand, skip_continuity=True, + ode_solver=OdeSolver.VARIATIONAL(), # This is a fake ode_solver to be able to use the variational integrator ) if qdot_bounds is None or not isinstance(qdot_bounds, BoundsList): diff --git a/tests/shard1/test_controltype_none.py b/tests/shard1/test_controltype_none.py index e5c9fcaac..e8cb062e1 100644 --- a/tests/shard1/test_controltype_none.py +++ b/tests/shard1/test_controltype_none.py @@ -188,6 +188,7 @@ def prepare_ocp( dynamic_function=custom_model.custom_dynamics, phase=i, expand_dynamics=True, + ode_solver=ode_solver, phase_dynamics=phase_dynamics, ) @@ -221,7 +222,6 @@ def prepare_ocp( x_bounds=x_bounds, objective_functions=objective_functions, constraints=constraints, - ode_solver=ode_solver, use_sx=use_sx, ) diff --git a/tests/shard1/test_dynamics.py b/tests/shard1/test_dynamics.py index 81ce546d6..6165a2365 100644 --- a/tests/shard1/test_dynamics.py +++ b/tests/shard1/test_dynamics.py @@ -124,6 +124,13 @@ def test_torque_driven(with_contact, with_external_force, cx, phase_dynamics): TestUtils.bioptim_folder() + "/examples/getting_started/models/2segments_4dof_2contacts.bioMod", external_force_set=external_forces, ) + nlp.dynamics_type = Dynamics( + DynamicsFcn.TORQUE_DRIVEN, + contact_type=[ContactType.RIGID_EXPLICIT] if with_contact else (), + expand_dynamics=True, + phase_dynamics=phase_dynamics, + numerical_data_timeseries=numerical_time_series, + ) nlp.cx = cx nlp.time_cx = cx.sym("time", 1, 1) @@ -138,16 +145,11 @@ def test_torque_driven(with_contact, with_external_force, cx, phase_dynamics): ocp = OptimalControlProgram(nlp, use_sx=(cx == SX)) nlp.control_type = ControlType.CONSTANT + NonLinearProgram.add( ocp, "dynamics_type", - Dynamics( - DynamicsFcn.TORQUE_DRIVEN, - contact_type=[ContactType.RIGID_EXPLICIT] if with_contact else (), - expand_dynamics=True, - phase_dynamics=phase_dynamics, - numerical_data_timeseries=numerical_time_series, - ), + nlp.dynamics_type, False, ) phase_index = [i for i in range(ocp.n_phases)] @@ -222,6 +224,14 @@ def test_torque_driven_soft_contacts_dynamics(contact_type, cx, phase_dynamics): TestUtils.bioptim_folder() + "/examples/muscle_driven_with_contact/models/2segments_4dof_2soft_contacts_1muscle.bioMod" ) + nlp.dynamics_type = Dynamics( + DynamicsFcn.TORQUE_DRIVEN, + contact_type=contact_type, + soft_contacts_dynamics=True if ContactType.SOFT_IMPLICIT in contact_type else False, + expand_dynamics=True, + phase_dynamics=phase_dynamics, + ) + nlp.ns = N_SHOOTING nlp.cx = cx nlp.time_cx = cx.sym("time", 1, 1) @@ -237,18 +247,10 @@ def test_torque_driven_soft_contacts_dynamics(contact_type, cx, phase_dynamics): ocp = OptimalControlProgram(nlp, use_sx=(cx == SX)) nlp.control_type = ControlType.CONSTANT - nlp.dynamics = Dynamics( - DynamicsFcn.TORQUE_DRIVEN, - contact_type=contact_type, - soft_contacts_dynamics=True if ContactType.SOFT_IMPLICIT in contact_type else False, - expand_dynamics=True, - phase_dynamics=phase_dynamics, - ) - NonLinearProgram.add( ocp, "dynamics_type", - nlp.dynamics, + nlp.dynamics_type, False, ) @@ -315,6 +317,14 @@ def test_torque_derivative_driven(with_contact, with_external_force, cx, phase_d TestUtils.bioptim_folder() + "/examples/getting_started/models/2segments_4dof_2contacts.bioMod", external_force_set=external_forces, ) + nlp.dynamics_type = Dynamics( + DynamicsFcn.TORQUE_DERIVATIVE_DRIVEN, + contact_type=[ContactType.RIGID_EXPLICIT] if with_contact else (), + expand_dynamics=True, + phase_dynamics=phase_dynamics, + numerical_data_timeseries=numerical_timeseries, + ) + nlp.cx = cx nlp.time_cx = cx.sym("time", 1, 1) nlp.dt = cx.sym("dt", 1, 1) @@ -328,17 +338,10 @@ def test_torque_derivative_driven(with_contact, with_external_force, cx, phase_d ocp = OptimalControlProgram(nlp, use_sx=(cx == SX)) nlp.control_type = ControlType.CONSTANT - nlp.dynamics = Dynamics( - DynamicsFcn.TORQUE_DERIVATIVE_DRIVEN, - contact_type=[ContactType.RIGID_EXPLICIT] if with_contact else [], - expand_dynamics=True, - phase_dynamics=phase_dynamics, - numerical_data_timeseries=numerical_timeseries, - ) NonLinearProgram.add( ocp, "dynamics_type", - nlp.dynamics, + nlp.dynamics_type, False, ) @@ -457,6 +460,14 @@ def test_torque_derivative_driven_soft_contacts_dynamics(contact_type, cx, phase TestUtils.bioptim_folder() + "/examples/muscle_driven_with_contact/models/2segments_4dof_2soft_contacts_1muscle.bioMod", ) + nlp.dynamics_type = Dynamics( + DynamicsFcn.TORQUE_DERIVATIVE_DRIVEN, + contact_type=contact_type, + soft_contacts_dynamics=True if ContactType.SOFT_IMPLICIT in contact_type else False, + expand_dynamics=True, + phase_dynamics=phase_dynamics, + ) + nlp.ns = N_SHOOTING nlp.cx = cx nlp.time_cx = cx.sym("time", 1, 1) @@ -472,17 +483,10 @@ def test_torque_derivative_driven_soft_contacts_dynamics(contact_type, cx, phase ocp = OptimalControlProgram(nlp, use_sx=(cx == SX)) nlp.control_type = ControlType.CONSTANT - nlp.dynamics = Dynamics( - DynamicsFcn.TORQUE_DERIVATIVE_DRIVEN, - contact_type=contact_type, - soft_contacts_dynamics=True if ContactType.SOFT_IMPLICIT in contact_type else False, - expand_dynamics=True, - phase_dynamics=phase_dynamics, - ) NonLinearProgram.add( ocp, "dynamics_type", - nlp.dynamics, + nlp.dynamics_type, False, ) @@ -533,6 +537,10 @@ def test_soft_contacts_dynamics_errors(dynamics, phase_dynamics): nlp.model = BiorbdModel( TestUtils.bioptim_folder() + "/examples/getting_started/models/2segments_4dof_2contacts.bioMod" ) + nlp.dynamics_type = Dynamics( + dynamics, soft_contacts_dynamics=True, expand_dynamics=True, phase_dynamics=phase_dynamics + ) + nlp.ns = N_SHOOTING nlp.cx = MX @@ -542,13 +550,10 @@ def test_soft_contacts_dynamics_errors(dynamics, phase_dynamics): ocp = OptimalControlProgram(nlp, use_sx=True) nlp.control_type = ControlType.CONSTANT - nlp.dynamics = ( - Dynamics(dynamics, soft_contacts_dynamics=True, expand_dynamics=True, phase_dynamics=phase_dynamics), - ) NonLinearProgram.add( ocp, "dynamics_type", - nlp.dynamics, + nlp.dynamics_type, False, ) phase_index = [i for i in range(ocp.n_phases)] @@ -584,6 +589,14 @@ def test_torque_activation_driven(with_contact, with_external_force, cx, phase_d TestUtils.bioptim_folder() + "/examples/getting_started/models/2segments_4dof_2contacts.bioMod", external_force_set=external_forces, ) + nlp.dynamics_type = Dynamics( + DynamicsFcn.TORQUE_ACTIVATIONS_DRIVEN, + contact_type=[ContactType.RIGID_EXPLICIT] if with_contact else (), + expand_dynamics=True, + phase_dynamics=phase_dynamics, + numerical_data_timeseries=numerical_timeseries, + ) + nlp.cx = cx nlp.time_cx = cx.sym("time", 1, 1) nlp.dt = cx.sym("dt", 1, 1) @@ -598,7 +611,7 @@ def test_torque_activation_driven(with_contact, with_external_force, cx, phase_d ocp = OptimalControlProgram(nlp, use_sx=(cx == SX)) nlp.control_type = ControlType.CONSTANT - nlp.dynamics = Dynamics( + nlp.dynamics_type = Dynamics( DynamicsFcn.TORQUE_ACTIVATIONS_DRIVEN, contact_type=[ContactType.RIGID_EXPLICIT] if with_contact else (), expand_dynamics=True, @@ -608,7 +621,7 @@ def test_torque_activation_driven(with_contact, with_external_force, cx, phase_d NonLinearProgram.add( ocp, "dynamics_type", - nlp.dynamics, + nlp.dynamics_type, False, ) phase_index = [i for i in range(ocp.n_phases)] @@ -711,6 +724,14 @@ def test_torque_activation_driven_with_residual_torque( TestUtils.bioptim_folder() + "/examples/torque_driven_ocp/models/2segments_2dof_2contacts.bioMod", external_force_set=external_forces, ) + nlp.dynamics_type = Dynamics( + DynamicsFcn.TORQUE_ACTIVATIONS_DRIVEN, + with_residual_torque=with_residual_torque, + expand_dynamics=True, + phase_dynamics=phase_dynamics, + numerical_data_timeseries=numerical_timeseries, + ) + nlp.cx = cx nlp.time_cx = cx.sym("time", 1, 1) nlp.dt = cx.sym("dt", 1, 1) @@ -724,17 +745,10 @@ def test_torque_activation_driven_with_residual_torque( ocp = OptimalControlProgram(nlp, use_sx=(cx == SX)) nlp.control_type = ControlType.CONSTANT - nlp.dynamics = Dynamics( - DynamicsFcn.TORQUE_ACTIVATIONS_DRIVEN, - with_residual_torque=with_residual_torque, - expand_dynamics=True, - phase_dynamics=phase_dynamics, - numerical_data_timeseries=numerical_timeseries, - ) NonLinearProgram.add( ocp, "dynamics_type", - nlp.dynamics, + nlp.dynamics_type, False, ) phase_index = [i for i in range(ocp.n_phases)] @@ -825,6 +839,10 @@ def test_torque_driven_free_floating_base(cx, phase_dynamics): nlp.model = BiorbdModel( TestUtils.bioptim_folder() + "/examples/getting_started/models/2segments_4dof_2contacts.bioMod" ) + nlp.dynamics_type = Dynamics( + DynamicsFcn.TORQUE_DRIVEN_FREE_FLOATING_BASE, expand_dynamics=True, phase_dynamics=phase_dynamics + ) + nlp.ns = N_SHOOTING nlp.cx = cx nlp.time_cx = cx.sym("time", 1, 1) @@ -840,13 +858,10 @@ def test_torque_driven_free_floating_base(cx, phase_dynamics): ocp = OptimalControlProgram(nlp, use_sx=(cx == SX)) nlp.control_type = ControlType.CONSTANT - nlp.dynamics = ( - Dynamics(DynamicsFcn.TORQUE_DRIVEN_FREE_FLOATING_BASE, expand_dynamics=True, phase_dynamics=phase_dynamics), - ) NonLinearProgram.add( ocp, "dynamics_type", - nlp.dynamics, + nlp.dynamics_type, False, ) phase_index = [i for i in range(ocp.n_phases)] @@ -900,6 +915,15 @@ def test_muscle_driven(with_excitations, with_contact, with_residual_torque, wit TestUtils.bioptim_folder() + "/examples/muscle_driven_ocp/models/arm26_with_contact.bioMod", external_force_set=external_forces, ) + nlp.dynamics_type = Dynamics( + DynamicsFcn.MUSCLE_DRIVEN, + with_residual_torque=with_residual_torque, + with_excitations=with_excitations, + contact_type=[ContactType.RIGID_EXPLICIT] if with_contact else (), + expand_dynamics=True, + phase_dynamics=phase_dynamics, + numerical_data_timeseries=numerical_timeseries, + ) nlp.cx = cx nlp.time_cx = cx.sym("time", 1, 1) @@ -916,19 +940,10 @@ def test_muscle_driven(with_excitations, with_contact, with_residual_torque, wit ocp = OptimalControlProgram(nlp, use_sx=(cx == SX)) nlp.control_type = ControlType.CONSTANT - nlp.dynamics = Dynamics( - DynamicsFcn.MUSCLE_DRIVEN, - with_residual_torque=with_residual_torque, - with_excitations=with_excitations, - contact_type=[ContactType.RIGID_EXPLICIT] if with_contact else (), - expand_dynamics=True, - phase_dynamics=phase_dynamics, - numerical_data_timeseries=numerical_timeseries, - ) NonLinearProgram.add( ocp, "dynamics_type", - nlp.dynamics, + nlp.dynamics_type, False, ) phase_index = [i for i in range(ocp.n_phases)] @@ -1075,6 +1090,11 @@ def test_joints_acceleration_driven(cx, phase_dynamics): # Prepare the program nlp = NonLinearProgram(phase_dynamics=phase_dynamics, use_sx=(cx == SX)) nlp.model = BiorbdModel(TestUtils.bioptim_folder() + "/examples/getting_started/models/double_pendulum.bioMod") + nlp.dynamics_type = Dynamics( + DynamicsFcn.JOINTS_ACCELERATION_DRIVEN, + expand_dynamics=True, + phase_dynamics=phase_dynamics, + ) nlp.ns = N_SHOOTING nlp.cx = cx @@ -1091,15 +1111,10 @@ def test_joints_acceleration_driven(cx, phase_dynamics): ocp = OptimalControlProgram(nlp, use_sx=(cx == SX)) nlp.control_type = ControlType.CONSTANT - nlp.dynamics = Dynamics( - DynamicsFcn.JOINTS_ACCELERATION_DRIVEN, - expand_dynamics=True, - phase_dynamics=phase_dynamics, - ) NonLinearProgram.add( ocp, "dynamics_type", - nlp.dynamics, + nlp.dynamics_type, False, ) np.random.seed(42) @@ -1152,6 +1167,14 @@ def configure(ocp, nlp, with_contact=None, numerical_data_timeseries=None, conta nlp.model = BiorbdModel( TestUtils.bioptim_folder() + "/examples/getting_started/models/2segments_4dof_2contacts.bioMod" ) + nlp.dynamics_type = Dynamics( + configure, + dynamic_function=custom_dynamic, + contact_type=contact_type, + expand_dynamics=True, + phase_dynamics=phase_dynamics, + ) + nlp.ns = N_SHOOTING nlp.cx = MX nlp.time_cx = nlp.cx.sym("time", 1, 1) @@ -1166,17 +1189,10 @@ def configure(ocp, nlp, with_contact=None, numerical_data_timeseries=None, conta ocp = OptimalControlProgram(nlp, use_sx=False) nlp.control_type = ControlType.CONSTANT - nlp.dynamics = Dynamics( - configure, - dynamic_function=custom_dynamic, - contact_type=contact_type, - expand_dynamics=True, - phase_dynamics=phase_dynamics, - ) NonLinearProgram.add( ocp, "dynamics_type", - nlp.dynamics, + nlp.dynamics_type, False, ) phase_index = [i for i in range(ocp.n_phases)] @@ -1237,7 +1253,11 @@ def test_with_contact_error(dynamics_fcn, phase_dynamics): # Dynamics dynamics = Dynamics( - dynamics_fcn, contact_type=[ContactType.RIGID_EXPLICIT], expand_dynamics=True, phase_dynamics=phase_dynamics + dynamics_fcn, + contact_type=[ContactType.RIGID_EXPLICIT], + ode_solver=OdeSolver.RK4(), + expand_dynamics=True, + phase_dynamics=phase_dynamics, ) # Path constraint @@ -1262,5 +1282,4 @@ def test_with_contact_error(dynamics_fcn, phase_dynamics): x_bounds=x_bounds, u_bounds=u_bounds, objective_functions=objective_functions, - ode_solver=OdeSolver.RK4(), ) diff --git a/tests/shard1/test_prepare_all_examples.py b/tests/shard1/test_prepare_all_examples.py index 1190d32e2..e2cb3d395 100644 --- a/tests/shard1/test_prepare_all_examples.py +++ b/tests/shard1/test_prepare_all_examples.py @@ -304,7 +304,7 @@ def test__getting_started__example_multiphase_different_ode_solvers(): with pytest.raises( RuntimeError, - match="ode_solver should be built an instance of OdeSolver or a list of OdeSolver", + match="ode_solver should be built an instance of OdeSolver", ): ocp_module.prepare_ocp( biorbd_model_path=bioptim_folder + "/models/cube.bioMod", diff --git a/tests/shard2/test_cost_function_integration.py b/tests/shard2/test_cost_function_integration.py index c6d52d418..182d8fb35 100644 --- a/tests/shard2/test_cost_function_integration.py +++ b/tests/shard2/test_cost_function_integration.py @@ -87,7 +87,9 @@ def prepare_ocp( raise ValueError("Wrong objective") # Dynamics - dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=True, phase_dynamics=phase_dynamics) + dynamics = Dynamics( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=True, phase_dynamics=phase_dynamics + ) # Path constraint x_bounds = BoundsList() @@ -112,7 +114,6 @@ def prepare_ocp( x_bounds=x_bounds, u_bounds=u_bounds, objective_functions=objective_functions, - ode_solver=ode_solver, use_sx=True, n_threads=1, control_type=control_type, diff --git a/tests/shard2/test_global_minimize_marker_velocity.py b/tests/shard2/test_global_minimize_marker_velocity.py index 7859e1ea5..56c05df40 100644 --- a/tests/shard2/test_global_minimize_marker_velocity.py +++ b/tests/shard2/test_global_minimize_marker_velocity.py @@ -105,6 +105,7 @@ def prepare_ocp( DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=not isinstance(ode_solver, OdeSolver.IRK), phase_dynamics=phase_dynamics, + ode_solver=ode_solver, ) # Path constraint @@ -133,7 +134,6 @@ def prepare_ocp( x_init=x_init, objective_functions=objective_functions, control_type=control_type, - ode_solver=ode_solver, ) diff --git a/tests/shard2/test_global_nmpc_final.py b/tests/shard2/test_global_nmpc_final.py index 367d3cd56..ccc40182a 100644 --- a/tests/shard2/test_global_nmpc_final.py +++ b/tests/shard2/test_global_nmpc_final.py @@ -62,7 +62,7 @@ def update_functions(_nmpc, cycle_idx, _sol): npt.assert_almost_equal(tau[:, -1], np.array([-0.00992505, 5.19414727, 2.34022319]), decimal=4) # check time - n_steps = nmpc.nlp[0].ode_solver.n_integration_steps + n_steps = nmpc.nlp[0].dynamics_type.ode_solver.n_integration_steps time = sol[0].stepwise_time(to_merge=SolutionMerge.NODES) assert time.shape == (n_cycles_total * cycle_len * (n_steps + 1) + 1, 1) assert time[0] == 0 @@ -197,7 +197,7 @@ def update_functions(_nmpc, cycle_idx, _sol): ) # check time - n_steps = nmpc.nlp[0].ode_solver.n_integration_steps + n_steps = nmpc.nlp[0].dynamics_type.ode_solver.n_integration_steps time = sol[0].stepwise_time(to_merge=SolutionMerge.NODES) assert time.shape == (n_cycles_total * cycle_len * (n_steps + 1) + 1, 1) assert time[0] == 0 diff --git a/tests/shard3/test_graph.py b/tests/shard3/test_graph.py index cde4014f4..31fd98824 100644 --- a/tests/shard3/test_graph.py +++ b/tests/shard3/test_graph.py @@ -322,7 +322,9 @@ def prepare_ocp_parameters( objective_functions.add(ObjectiveFcn.Mayer.MINIMIZE_TIME, weight=10) # Dynamics - dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=True, phase_dynamics=phase_dynamics) + dynamics = Dynamics( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=True, phase_dynamics=phase_dynamics + ) # Path constraint x_bounds = BoundsList() @@ -350,7 +352,6 @@ def prepare_ocp_parameters( parameter_objectives=parameter_objectives, parameter_bounds=parameter_bounds, parameter_init=parameter_init, - ode_solver=ode_solver, use_sx=use_sx, ) @@ -417,7 +418,9 @@ def prepare_ocp_custom_objectives( objective_functions.add(ObjectiveFcn.Mayer.MINIMIZE_MARKERS, list_index=7, index=[1, 2], target=target) # Dynamics - dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=True, phase_dynamics=phase_dynamics) + dynamics = Dynamics( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=True, phase_dynamics=phase_dynamics + ) # Path constraint x_bounds = BoundsList() @@ -441,7 +444,6 @@ def prepare_ocp_custom_objectives( x_bounds=x_bounds, u_bounds=u_bounds, objective_functions=objective_functions, - ode_solver=ode_solver, ) diff --git a/tests/shard3/test_ligaments.py b/tests/shard3/test_ligaments.py index 553d21e17..dbe0c941d 100644 --- a/tests/shard3/test_ligaments.py +++ b/tests/shard3/test_ligaments.py @@ -44,6 +44,8 @@ def test_torque_driven_with_ligament(with_ligament, cx, phase_dynamics): nlp.model = BiorbdModel( TestUtils.bioptim_folder() + "/examples/torque_driven_ocp/models/mass_point_with_ligament.bioMod" ) + nlp.dynamics_type = Dynamics(DynamicsFcn.TORQUE_DRIVEN, with_ligament=with_ligament) + nlp.ns = 5 nlp.cx = cx nlp.time_cx = cx.sym("time", 1, 1) @@ -57,10 +59,11 @@ def test_torque_driven_with_ligament(with_ligament, cx, phase_dynamics): nlp.u_bounds = np.zeros((nlp.model.nb_q, 1)) ocp = OptimalControlProgram(nlp, use_sx=(cx == SX)) nlp.control_type = ControlType.CONSTANT + NonLinearProgram.add( ocp, "dynamics_type", - Dynamics(DynamicsFcn.TORQUE_DRIVEN, with_ligament=with_ligament), + nlp.dynamics_type, False, ) phase_index = [i for i in range(ocp.n_phases)] @@ -101,6 +104,8 @@ def test_torque_derivative_driven_with_ligament(with_ligament, cx, phase_dynamic nlp.model = BiorbdModel( TestUtils.bioptim_folder() + "/examples/torque_driven_ocp/models/mass_point_with_ligament.bioMod" ) + nlp.dynamics_type = Dynamics(DynamicsFcn.TORQUE_DERIVATIVE_DRIVEN, with_ligament=with_ligament) + nlp.ns = 5 nlp.cx = cx nlp.time_cx = cx.sym("time", 1, 1) @@ -118,7 +123,7 @@ def test_torque_derivative_driven_with_ligament(with_ligament, cx, phase_dynamic NonLinearProgram.add( ocp, "dynamics_type", - Dynamics(DynamicsFcn.TORQUE_DERIVATIVE_DRIVEN, with_ligament=with_ligament), + nlp.dynamics_type, False, ) @@ -160,6 +165,8 @@ def test_torque_activation_driven_with_ligament(with_ligament, cx, phase_dynamic nlp.model = BiorbdModel( TestUtils.bioptim_folder() + "/examples/torque_driven_ocp/models/mass_point_with_ligament.bioMod" ) + nlp.dynamics_type = Dynamics(DynamicsFcn.TORQUE_ACTIVATIONS_DRIVEN, with_ligament=with_ligament) + nlp.ns = 5 nlp.cx = cx nlp.time_cx = cx.sym("time", 1, 1) @@ -175,7 +182,7 @@ def test_torque_activation_driven_with_ligament(with_ligament, cx, phase_dynamic NonLinearProgram.add( ocp, "dynamics_type", - Dynamics(DynamicsFcn.TORQUE_ACTIVATIONS_DRIVEN, with_ligament=with_ligament), + nlp.dynamics_type, False, ) phase_index = [i for i in range(ocp.n_phases)] @@ -217,6 +224,11 @@ def test_muscle_driven_with_ligament(with_ligament, cx, phase_dynamics): nlp.model = BiorbdModel( TestUtils.bioptim_folder() + "/examples/muscle_driven_ocp/models/arm26_with_ligament.bioMod" ) + nlp.dynamics_type = Dynamics( + DynamicsFcn.MUSCLE_DRIVEN, + with_ligament=with_ligament, + ) + nlp.ns = 5 nlp.cx = cx nlp.time_cx = cx.sym("time", 1, 1) @@ -233,10 +245,7 @@ def test_muscle_driven_with_ligament(with_ligament, cx, phase_dynamics): NonLinearProgram.add( ocp, "dynamics_type", - Dynamics( - DynamicsFcn.MUSCLE_DRIVEN, - with_ligament=with_ligament, - ), + nlp.dynamics_type, False, ) phase_index = [i for i in range(ocp.n_phases)] diff --git a/tests/shard3/test_multinode_constraints.py b/tests/shard3/test_multinode_constraints.py index 565ea5829..cb2c3a1b7 100644 --- a/tests/shard3/test_multinode_constraints.py +++ b/tests/shard3/test_multinode_constraints.py @@ -30,9 +30,15 @@ def prepare_ocp(biorbd_model_path, phase_1, phase_2, phase_dynamics) -> OptimalC # Dynamics dynamics = DynamicsList() - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=True, phase_dynamics=phase_dynamics) - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=True, phase_dynamics=phase_dynamics) - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=True, phase_dynamics=phase_dynamics) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=OdeSolver.RK4(), expand_dynamics=True, phase_dynamics=phase_dynamics + ) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=OdeSolver.RK4(), expand_dynamics=True, phase_dynamics=phase_dynamics + ) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=OdeSolver.RK4(), expand_dynamics=True, phase_dynamics=phase_dynamics + ) multinode_constraints = MultinodeConstraintList() # hard constraint @@ -82,7 +88,6 @@ def prepare_ocp(biorbd_model_path, phase_1, phase_2, phase_dynamics) -> OptimalC u_bounds=u_bounds, objective_functions=objective_functions, multinode_constraints=multinode_constraints, - ode_solver=OdeSolver.RK4(), ) diff --git a/tests/shard3/test_passive_torque.py b/tests/shard3/test_passive_torque.py index 15975a410..45bf04748 100644 --- a/tests/shard3/test_passive_torque.py +++ b/tests/shard3/test_passive_torque.py @@ -42,6 +42,12 @@ def test_torque_driven_with_passive_torque(with_passive_torque, cx, phase_dynami nlp.model = BiorbdModel( TestUtils.bioptim_folder() + "/examples/getting_started/models/2segments_4dof_2contacts.bioMod" ) + nlp.dynamics_type = Dynamics( + DynamicsFcn.TORQUE_DRIVEN, + with_passive_torque=with_passive_torque, + phase_dynamics=phase_dynamics, + ) + nlp.ns = 5 nlp.cx = cx nlp.time_cx = cx.sym("time", 1, 1) @@ -55,14 +61,11 @@ def test_torque_driven_with_passive_torque(with_passive_torque, cx, phase_dynami nlp.u_bounds = np.zeros((nlp.model.nb_q, 1)) ocp = OptimalControlProgram(nlp, use_sx=(cx == SX)) nlp.control_type = ControlType.CONSTANT + NonLinearProgram.add( ocp, "dynamics_type", - Dynamics( - DynamicsFcn.TORQUE_DRIVEN, - with_passive_torque=with_passive_torque, - phase_dynamics=phase_dynamics, - ), + nlp.dynamics_type, False, ) phase_index = [i for i in range(ocp.n_phases)] @@ -102,6 +105,12 @@ def test_torque_derivative_driven_with_passive_torque(with_passive_torque, cx, p nlp.model = BiorbdModel( TestUtils.bioptim_folder() + "/examples/getting_started/models/2segments_4dof_2contacts.bioMod" ) + nlp.dynamics_type = Dynamics( + DynamicsFcn.TORQUE_DERIVATIVE_DRIVEN, + with_passive_torque=with_passive_torque, + phase_dynamics=phase_dynamics, + ) + nlp.ns = 5 nlp.cx = cx nlp.time_cx = cx.sym("time", 1, 1) @@ -119,11 +128,7 @@ def test_torque_derivative_driven_with_passive_torque(with_passive_torque, cx, p NonLinearProgram.add( ocp, "dynamics_type", - Dynamics( - DynamicsFcn.TORQUE_DERIVATIVE_DRIVEN, - with_passive_torque=with_passive_torque, - phase_dynamics=phase_dynamics, - ), + nlp.dynamics_type, False, ) @@ -192,6 +197,13 @@ def test_torque_activation_driven_with_passive_torque(with_passive_torque, with_ nlp.model = BiorbdModel( TestUtils.bioptim_folder() + "/examples/getting_started/models/2segments_4dof_2contacts.bioMod" ) + nlp.dynamics_type = Dynamics( + DynamicsFcn.TORQUE_ACTIVATIONS_DRIVEN, + with_passive_torque=with_passive_torque, + with_residual_torque=with_residual_torque, + phase_dynamics=phase_dynamics, + ) + nlp.ns = 5 nlp.cx = cx nlp.time_cx = cx.sym("time", 1, 1) @@ -203,15 +215,11 @@ def test_torque_activation_driven_with_passive_torque(with_passive_torque, with_ nlp.u_bounds = np.zeros((nlp.model.nb_q, 1)) ocp = OptimalControlProgram(nlp, use_sx=(cx == SX)) nlp.control_type = ControlType.CONSTANT + NonLinearProgram.add( ocp, "dynamics_type", - Dynamics( - DynamicsFcn.TORQUE_ACTIVATIONS_DRIVEN, - with_passive_torque=with_passive_torque, - with_residual_torque=with_residual_torque, - phase_dynamics=phase_dynamics, - ), + nlp.dynamics_type, False, ) phase_index = [i for i in range(ocp.n_phases)] @@ -303,6 +311,12 @@ def test_muscle_driven_with_passive_torque(with_passive_torque, cx, phase_dynami # Prepare the program nlp = NonLinearProgram(phase_dynamics=phase_dynamics, use_sx=(cx == SX)) nlp.model = BiorbdModel(TestUtils.bioptim_folder() + "/examples/muscle_driven_ocp/models/arm26_with_contact.bioMod") + nlp.dynamics_type = Dynamics( + DynamicsFcn.MUSCLE_DRIVEN, + with_passive_torque=with_passive_torque, + phase_dynamics=phase_dynamics, + ) + nlp.ns = 5 nlp.cx = cx nlp.time_cx = cx.sym("time", 1, 1) @@ -316,14 +330,11 @@ def test_muscle_driven_with_passive_torque(with_passive_torque, cx, phase_dynami ocp = OptimalControlProgram(nlp, use_sx=(cx == SX)) nlp.control_type = ControlType.CONSTANT + NonLinearProgram.add( ocp, "dynamics_type", - Dynamics( - DynamicsFcn.MUSCLE_DRIVEN, - with_passive_torque=with_passive_torque, - phase_dynamics=phase_dynamics, - ), + nlp.dynamics_type, False, ) phase_index = [i for i in range(ocp.n_phases)] diff --git a/tests/shard4/test_update_bounds_and_init.py b/tests/shard4/test_update_bounds_and_init.py index 094727b00..2ed9b8c42 100644 --- a/tests/shard4/test_update_bounds_and_init.py +++ b/tests/shard4/test_update_bounds_and_init.py @@ -273,7 +273,7 @@ def test_update_noised_init_rk4(interpolation, phase_dynamics): phase_time = 1.0 dynamics = DynamicsList() - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, phase_dynamics=phase_dynamics) + dynamics.add(DynamicsFcn.TORQUE_DRIVEN, ode_solver=OdeSolver.RK4(), phase_dynamics=phase_dynamics) x_init = InitialGuessList() x_init["q"] = [0] * bio_model.nb_q @@ -285,7 +285,6 @@ def test_update_noised_init_rk4(interpolation, phase_dynamics): dynamics, n_shooting=ns, phase_time=phase_time, - ode_solver=OdeSolver.RK4(), x_init=x_init, u_init=u_init, ) @@ -979,10 +978,10 @@ def test_update_noised_initial_guess_collocation(interpolation, phase_dynamics): ntau = bio_model.nb_tau ns = 3 phase_time = 1.0 - solver = OdeSolver.COLLOCATION(polynomial_degree=1) + ode_solver = OdeSolver.COLLOCATION(polynomial_degree=1) dynamics = DynamicsList() - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, phase_dynamics=phase_dynamics) + dynamics.add(DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, phase_dynamics=phase_dynamics) x_init = InitialGuessList() x_init["q"] = [0] * bio_model.nb_q @@ -994,7 +993,6 @@ def test_update_noised_initial_guess_collocation(interpolation, phase_dynamics): dynamics, n_shooting=ns, phase_time=phase_time, - ode_solver=solver, x_init=x_init, u_init=u_init, ) @@ -1037,9 +1035,9 @@ def test_update_noised_initial_guess_collocation(interpolation, phase_dynamics): x.add("qdot", x_init[nq:, :], interpolation=interpolation) u.add("tau", np.zeros((ntau, ns)), interpolation=interpolation) elif interpolation == InterpolationType.ALL_POINTS: - x_init = np.zeros((nq * 2, ns * (solver.polynomial_degree + 1) + 1)) + x_init = np.zeros((nq * 2, ns * (ode_solver.polynomial_degree + 1) + 1)) for i in range(nq * 2): - x_init[i, :] = np.linspace(0, 1, ns * (solver.polynomial_degree + 1) + 1) + x_init[i, :] = np.linspace(0, 1, ns * (ode_solver.polynomial_degree + 1) + 1) x.add("q", x_init[:nq, :], interpolation=interpolation) x.add("qdot", x_init[nq:, :], interpolation=interpolation) u.add("tau", np.zeros((ntau, ns)), interpolation=interpolation) diff --git a/tests/shard4/test_variable_time.py b/tests/shard4/test_variable_time.py index 810538699..bc6d64006 100644 --- a/tests/shard4/test_variable_time.py +++ b/tests/shard4/test_variable_time.py @@ -54,9 +54,15 @@ def prepare_ocp(phase_time_constraint, use_parameter, phase_dynamics): # Dynamics dynamics = DynamicsList() - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, phase=0, expand_dynamics=True, phase_dynamics=phase_dynamics) - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, phase=1, expand_dynamics=True, phase_dynamics=phase_dynamics) - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, phase=2, expand_dynamics=True, phase_dynamics=phase_dynamics) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, phase=0, expand_dynamics=True, phase_dynamics=phase_dynamics + ) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, phase=1, expand_dynamics=True, phase_dynamics=phase_dynamics + ) + dynamics.add( + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, phase=2, expand_dynamics=True, phase_dynamics=phase_dynamics + ) # Constraints constraints = ConstraintList() @@ -134,7 +140,6 @@ def my_parameter_function(bio_model, value, extra_value): u_bounds=u_bounds, objective_functions=objective_functions, constraints=constraints, - ode_solver=ode_solver, parameters=parameters, parameter_init=parameter_init, parameter_bounds=parameter_bounds, diff --git a/tests/shard5/test_global_stochastic_collocation.py b/tests/shard5/test_global_stochastic_collocation.py index 2675f0a46..7793835e0 100644 --- a/tests/shard5/test_global_stochastic_collocation.py +++ b/tests/shard5/test_global_stochastic_collocation.py @@ -135,7 +135,7 @@ def test_arm_reaching_torque_driven_collocations(use_sx: bool): duration = sol_socp.decision_time()[-1] dt = duration / n_shooting p_sol = vertcat(ocp.nlp[0].model.motor_noise_magnitude, ocp.nlp[0].model.sensory_noise_magnitude) - polynomial_degree = socp.nlp[0].ode_solver.polynomial_degree + polynomial_degree = socp.nlp[0].dynamics_type.ode_solver.polynomial_degree # Constraint values x_opt = vertcat(q_sol, qdot_sol) diff --git a/tests/shard5/test_helper_socp.py b/tests/shard5/test_helper_socp.py index bc0eda5d5..a164f88f5 100644 --- a/tests/shard5/test_helper_socp.py +++ b/tests/shard5/test_helper_socp.py @@ -3,7 +3,6 @@ from bioptim import SocpType, PhaseDynamics from bioptim.optimization.stochastic_optimal_control_program import ( _check_multi_threading_and_problem_type, - _check_has_no_ode_solver_defined, _check_has_no_phase_dynamics_shared_during_the_phase, ) @@ -26,16 +25,6 @@ def test_check_multi_threading_and_problem_type_collocation(): _check_multi_threading_and_problem_type(SocpType.COLLOCATION) -# Tests for _check_has_no_ode_solver_defined() -def test_check_has_no_ode_solver_defined_with_ode_solver(): - with pytest.raises(ValueError, match="The ode_solver cannot be defined for a stochastic ocp"): - _check_has_no_ode_solver_defined(ode_solver="some_solver") - - -def test_check_has_no_ode_solver_defined_without_ode_solver(): - _check_has_no_ode_solver_defined() - - # Tests for _check_has_no_phase_dynamics_shared_during_the_phase() def test_check_has_no_phase_dynamics_shared_during_the_phase_shared_dynamics(): with pytest.raises(ValueError, match="The dynamics cannot be SHARED_DURING_THE_PHASE"): diff --git a/tests/shard6/test_configure_problem.py b/tests/shard6/test_configure_problem.py index 131b1a11b..082465c8b 100644 --- a/tests/shard6/test_configure_problem.py +++ b/tests/shard6/test_configure_problem.py @@ -13,12 +13,20 @@ VariableScalingList, FatigueList, XiaFatigue, + Dynamics, + DynamicsFcn, ) from ..utils import TestUtils class OptimalControlProgram: def __init__(self, nlp, use_sx): + nlp.time_cx = nlp.cx.sym("time", 1, 1) + nlp.dt = nlp.cx.sym("dt", 1, 1) + nlp.x_scaling = VariableScalingList() + nlp.u_scaling = VariableScalingList() + nlp.a_scaling = VariableScalingList() + self.cx = nlp.cx self.phase_dynamics = PhaseDynamics.SHARED_DURING_THE_PHASE self.n_phases = 1 @@ -28,6 +36,16 @@ def __init__(self, nlp, use_sx): self.parameters.initialize(parameters_list) self.implicit_constraints = ConstraintList() self.n_threads = 1 + nlp.dynamics_type = Dynamics( + DynamicsFcn.TORQUE_DRIVEN, + ) + NonLinearProgram.add( + self, + "dynamics_type", + nlp.dynamics_type, + False, + ) + nlp.initialize(nlp.cx) @pytest.mark.parametrize("cx", [MX, SX]) @@ -37,12 +55,6 @@ def test_configures(cx): nlp = NonLinearProgram(phase_dynamics=PhaseDynamics.SHARED_DURING_THE_PHASE, use_sx=(cx == SX)) nlp.ns = 30 nlp.cx = MX - nlp.time_cx = nlp.cx.sym("time", 1, 1) - nlp.dt = nlp.cx.sym("dt", 1, 1) - nlp.initialize(nlp.cx) - nlp.x_scaling = VariableScalingList() - nlp.u_scaling = VariableScalingList() - nlp.a_scaling = VariableScalingList() ocp = OptimalControlProgram(nlp, use_sx=(cx == SX)) nlp.model = BiorbdModel( @@ -155,12 +167,6 @@ def test_configure_soft_contacts(cx): nlp = NonLinearProgram(phase_dynamics=PhaseDynamics.SHARED_DURING_THE_PHASE, use_sx=(cx == SX)) nlp.ns = 30 nlp.cx = MX - nlp.time_cx = nlp.cx.sym("time", 1, 1) - nlp.dt = nlp.cx.sym("dt", 1, 1) - nlp.initialize(nlp.cx) - nlp.x_scaling = VariableScalingList() - nlp.u_scaling = VariableScalingList() - nlp.a_scaling = VariableScalingList() ocp = OptimalControlProgram(nlp, use_sx=(cx == SX)) nlp.model = BiorbdModel( @@ -183,12 +189,6 @@ def test_configure_muscles(cx): nlp = NonLinearProgram(phase_dynamics=PhaseDynamics.SHARED_DURING_THE_PHASE, use_sx=(cx == SX)) nlp.ns = 30 nlp.cx = MX - nlp.time_cx = nlp.cx.sym("time", 1, 1) - nlp.dt = nlp.cx.sym("dt", 1, 1) - nlp.initialize(nlp.cx) - nlp.x_scaling = VariableScalingList() - nlp.u_scaling = VariableScalingList() - nlp.a_scaling = VariableScalingList() ocp = OptimalControlProgram(nlp, use_sx=(cx == SX)) nlp.model = BiorbdModel( diff --git a/tests/shard6/test_dt_dependent.py b/tests/shard6/test_dt_dependent.py index b5b964052..49f596f30 100644 --- a/tests/shard6/test_dt_dependent.py +++ b/tests/shard6/test_dt_dependent.py @@ -99,6 +99,7 @@ def prepare_ocp_state_as_time( custom_configure, dynamic_function=dynamics, phase=i, + ode_solver=OdeSolver.RK4(n_integration_steps=5), expand_dynamics=True, phase_dynamics=phase_dynamics, ) @@ -135,7 +136,6 @@ def prepare_ocp_state_as_time( objective_functions=objective_functions, control_type=control_type, use_sx=use_sx, - ode_solver=OdeSolver.RK4(n_integration_steps=5), ) diff --git a/tests/shard6/test_time_dependent_custom_model.py b/tests/shard6/test_time_dependent_custom_model.py index 9f59281ab..bd2b1a3b8 100644 --- a/tests/shard6/test_time_dependent_custom_model.py +++ b/tests/shard6/test_time_dependent_custom_model.py @@ -301,6 +301,7 @@ def prepare_ocp( expand_dynamics=True, expand_continuity=False, phase=i, + ode_solver=OdeSolver.RK4(n_integration_steps=1), phase_dynamics=PhaseDynamics.SHARED_DURING_THE_PHASE, ) @@ -381,7 +382,6 @@ def prepare_ocp( parameter_init=parameters_init, control_type=ControlType.CONSTANT, use_sx=use_sx, - ode_solver=OdeSolver.RK4(n_integration_steps=1), ) diff --git a/tests/shard6/test_time_dependent_problems.py b/tests/shard6/test_time_dependent_problems.py index 5a69cbb48..b1191acde 100644 --- a/tests/shard6/test_time_dependent_problems.py +++ b/tests/shard6/test_time_dependent_problems.py @@ -172,6 +172,7 @@ def prepare_ocp( custom_configure, dynamic_function=time_dynamic, phase=i, + ode_solver=ode_solver, expand_dynamics=expand, phase_dynamics=phase_dynamics, ) @@ -218,7 +219,6 @@ def prepare_ocp( x_bounds=x_bounds, u_bounds=u_bounds, objective_functions=objective_functions, - ode_solver=ode_solver, control_type=control_type, use_sx=use_sx, ) diff --git a/tests/utils.py b/tests/utils.py index 6f4eb98f6..1dc1ada1e 100644 --- a/tests/utils.py +++ b/tests/utils.py @@ -120,7 +120,7 @@ def assert_warm_start(ocp, sol, state_decimal=2, control_decimal=2, param_decima @staticmethod def simulate(sol: Solution, decimal_value=7): - if sum([nlp.ode_solver.is_direct_collocation for nlp in sol.ocp.nlp]): + if sum([nlp.dynamics_type.ode_solver.is_direct_collocation for nlp in sol.ocp.nlp]): with pytest.raises( ValueError, match="When the ode_solver of the Optimal Control Problem is OdeSolver.COLLOCATION, " @@ -134,7 +134,7 @@ def simulate(sol: Solution, decimal_value=7): ) return - if sum([isinstance(nlp.ode_solver, OdeSolver.TRAPEZOIDAL) for nlp in sol.ocp.nlp]): + if sum([isinstance(nlp.dynamics_type.ode_solver, OdeSolver.TRAPEZOIDAL) for nlp in sol.ocp.nlp]): with pytest.raises( ValueError, match="When the ode_solver of the Optimal Control Problem is OdeSolver.TRAPEZOIDAL, "