From 27bbc3969438f37c31b53400bdfcaa1c9db357fa Mon Sep 17 00:00:00 2001 From: Charbie Date: Tue, 1 Apr 2025 15:16:55 +0400 Subject: [PATCH 01/20] moved ode_solver --- bioptim/__init__.py | 21 -------------- bioptim/examples/acados/cube.py | 5 ++-- bioptim/examples/custom_model/main.py | 2 +- .../fatigue/static_arm_with_fatigue.py | 2 +- .../getting_started/custom_constraint.py | 5 ++-- .../getting_started/custom_dynamics.py | 3 +- .../getting_started/custom_initial_guess.py | 6 ++-- .../getting_started/custom_objectives.py | 4 +-- .../getting_started/custom_parameters.py | 3 +- .../custom_phase_transitions.py | 9 +++--- .../example_continuity_as_objective.py | 5 ++-- .../example_cyclic_movement.py | 3 +- .../example_external_forces.py | 2 +- .../example_inequality_constraint.py | 3 +- .../example_joints_acceleration_driven.py | 3 +- .../example_multinode_constraints.py | 7 ++--- .../example_multinode_objective.py | 3 +- .../getting_started/example_multiphase.py | 7 ++--- .../example_parameter_scaling.py | 5 ++-- .../example_variable_scaling.py | 3 +- bioptim/examples/getting_started/pendulum.py | 3 +- .../pendulum_constrained_states_controls.py | 3 +- .../holonomic_constraints/two_pendulums.py | 3 +- .../two_pendulums_algebraic.py | 2 +- .../muscle_activations_tracker.py | 2 +- .../muscle_excitations_tracker.py | 2 +- .../examples/muscle_driven_ocp/static_arm.py | 2 +- .../static_arm_with_contact.py | 3 +- ...nequality_constraint_muscle_excitations.py | 2 +- .../muscle_activations_contacts_tracker.py | 3 +- .../multiphase_time_constraint.py | 7 ++--- .../pendulum_min_time_Mayer.py | 3 +- .../optimal_time_ocp/time_constraint.py | 3 +- bioptim/examples/sqp_method/pendulum.py | 3 +- .../obstacle_avoidance_direct_collocation.py | 12 ++++---- .../rockit_matrix_lyapunov.py | 7 ++--- .../symmetry_by_constraint.py | 3 +- .../symmetry_by_mapping.py | 3 +- .../example_pendulum_time_dependent.py | 5 ++-- .../torque_driven_ocp/example_quaternions.py | 3 +- .../torque_driven_ocp/example_soft_contact.py | 4 +-- .../maximize_predicted_height_CoM.py | 3 +- .../ocp_mass_with_ligament.py | 2 +- .../pendulum_with_passive_torque.py | 2 +- bioptim/examples/torque_driven_ocp/slider.py | 7 ++--- .../torque_activation_driven.py | 2 +- .../track_markers_2D_pendulum.py | 3 +- .../track_markers_with_torque_actuators.py | 7 ++--- .../torque_driven_ocp/trampo_quaternions.py | 3 +- bioptim/examples/track/optimal_estimation.py | 6 ++-- .../examples/track/track_marker_on_segment.py | 3 +- bioptim/examples/track/track_segment_on_rt.py | 3 +- bioptim/gui/serializable_class.py | 3 -- bioptim/optimization/non_linear_program.py | 18 ++++++------ .../optimization/optimal_control_program.py | 28 +++++++------------ .../receding_horizon_optimization.py | 1 - .../stochastic_optimal_control_program.py | 24 ++++------------ .../variational_optimal_control_program.py | 7 +---- 58 files changed, 113 insertions(+), 188 deletions(-) diff --git a/bioptim/__init__.py b/bioptim/__init__.py index 059f79cc4..dab7e4c66 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/examples/acados/cube.py b/bioptim/examples/acados/cube.py index cf80559c7..30c32d60f 100644 --- a/bioptim/examples/acados/cube.py +++ b/bioptim/examples/acados/cube.py @@ -24,7 +24,9 @@ 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 +45,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..7f54968a9 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 2e53a5a30..b960a6a9b 100644 --- a/bioptim/examples/getting_started/custom_dynamics.py +++ b/bioptim/examples/getting_started/custom_dynamics.py @@ -161,6 +161,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, ) @@ -168,6 +169,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, ) @@ -201,7 +203,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 9227dce53..5a8d452e1 100644 --- a/bioptim/examples/getting_started/custom_initial_guess.py +++ b/bioptim/examples/getting_started/custom_initial_guess.py @@ -126,7 +126,10 @@ 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() @@ -248,7 +251,6 @@ def prepare_ocp( u_init=u_init, objective_functions=objective_functions, constraints=constraints, - ode_solver=ode_solver, x_scaling=x_scaling, xdot_scaling=xdot_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..a0ff0178c 100644 --- a/bioptim/examples/getting_started/custom_objectives.py +++ b/bioptim/examples/getting_started/custom_objectives.py @@ -120,7 +120,8 @@ 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 +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/getting_started/custom_parameters.py b/bioptim/examples/getting_started/custom_parameters.py index 538a95d74..6826d6df9 100644 --- a/bioptim/examples/getting_started/custom_parameters.py +++ b/bioptim/examples/getting_started/custom_parameters.py @@ -221,7 +221,7 @@ 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 +249,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..f9ffa8070 100644 --- a/bioptim/examples/getting_started/custom_phase_transitions.py +++ b/bioptim/examples/getting_started/custom_phase_transitions.py @@ -115,10 +115,10 @@ 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 +201,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..b09192bba 100644 --- a/bioptim/examples/getting_started/example_cyclic_movement.py +++ b/bioptim/examples/getting_started/example_cyclic_movement.py @@ -75,7 +75,7 @@ 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 +115,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 3d5d664d8..5dff2d3b0 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 76f3963f3..21ae55620 100644 --- a/bioptim/examples/getting_started/example_inequality_constraint.py +++ b/bioptim/examples/getting_started/example_inequality_constraint.py @@ -93,7 +93,7 @@ def prepare_ocp( # Dynamics dynamics = DynamicsList() dynamics.add( - DynamicsFcn.TORQUE_DRIVEN, with_contact=True, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, with_contact=True, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics ) # Constraints @@ -150,7 +150,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..a01e00108 100644 --- a/bioptim/examples/getting_started/example_multinode_constraints.py +++ b/bioptim/examples/getting_started/example_multinode_constraints.py @@ -116,9 +116,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, 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 +197,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..83078bd59 100644 --- a/bioptim/examples/getting_started/example_multinode_objective.py +++ b/bioptim/examples/getting_started/example_multinode_objective.py @@ -90,7 +90,7 @@ 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 +115,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..94b2b5be4 100644 --- a/bioptim/examples/getting_started/example_multiphase.py +++ b/bioptim/examples/getting_started/example_multiphase.py @@ -120,9 +120,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, 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() @@ -162,7 +162,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..1863484e5 100644 --- a/bioptim/examples/getting_started/example_parameter_scaling.py +++ b/bioptim/examples/getting_started/example_parameter_scaling.py @@ -80,7 +80,7 @@ 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 +104,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 +206,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 +242,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..a8020cd31 100644 --- a/bioptim/examples/getting_started/example_variable_scaling.py +++ b/bioptim/examples/getting_started/example_variable_scaling.py @@ -72,7 +72,7 @@ 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 +107,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..036cb7b7d 100644 --- a/bioptim/examples/getting_started/pendulum.py +++ b/bioptim/examples/getting_started/pendulum.py @@ -79,7 +79,7 @@ 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 +114,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..e59ca982e 100644 --- a/bioptim/examples/getting_started/pendulum_constrained_states_controls.py +++ b/bioptim/examples/getting_started/pendulum_constrained_states_controls.py @@ -85,7 +85,7 @@ 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 +123,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 1b043c786..289a67960 100644 --- a/bioptim/examples/holonomic_constraints/two_pendulums_algebraic.py +++ b/bioptim/examples/holonomic_constraints/two_pendulums_algebraic.py @@ -143,6 +143,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, ) @@ -218,7 +219,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 54132e248..d44e82acf 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, ) diff --git a/bioptim/examples/muscle_driven_ocp/muscle_excitations_tracker.py b/bioptim/examples/muscle_driven_ocp/muscle_excitations_tracker.py index c83afebcb..76f3b06d6 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, ) 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 3faf69108..57cd2ce5f 100644 --- a/bioptim/examples/muscle_driven_ocp/static_arm_with_contact.py +++ b/bioptim/examples/muscle_driven_ocp/static_arm_with_contact.py @@ -65,7 +65,7 @@ def prepare_ocp( # Dynamics dynamics = DynamicsList() - dynamics.add(DynamicsFcn.MUSCLE_DRIVEN, with_residual_torque=True, with_contact=True) + dynamics.add(DynamicsFcn.MUSCLE_DRIVEN, with_residual_torque=True, with_contact=True, ode_solver=ode_solver) # raise RuntimeError("This example is broken, since contact dynamics with muscle is not implemented") # Path constraint @@ -100,7 +100,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 8b09516c1..aad732339 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 @@ -50,6 +50,7 @@ def prepare_ocp(biorbd_model_path, phase_time, n_shooting, min_bound, ode_solver with_excitations=True, with_residual_torque=True, with_contact=True, + ode_solver=ode_solver, expand_dynamics=expand_dynamics, ) @@ -113,7 +114,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 223a91c3f..32140003e 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 @@ -55,7 +55,7 @@ def prepare_ocp( # Dynamics dynamics = DynamicsList() - dynamics.add(DynamicsFcn.MUSCLE_DRIVEN, with_residual_torque=True, with_contact=True) + dynamics.add(DynamicsFcn.MUSCLE_DRIVEN, with_residual_torque=True, with_contact=True, ode_solver=ode_solver) # Path constraint q_at_first_node = [0, 0, -0.75, 0.75] @@ -91,7 +91,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..359d1054c 100644 --- a/bioptim/examples/optimal_time_ocp/multiphase_time_constraint.py +++ b/bioptim/examples/optimal_time_ocp/multiphase_time_constraint.py @@ -97,10 +97,10 @@ 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 +160,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..a23cd1743 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,7 @@ 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 +111,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..846b0e486 100644 --- a/bioptim/examples/optimal_time_ocp/time_constraint.py +++ b/bioptim/examples/optimal_time_ocp/time_constraint.py @@ -75,7 +75,7 @@ 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 +104,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..020e7de3b 100644 --- a/bioptim/examples/sqp_method/pendulum.py +++ b/bioptim/examples/sqp_method/pendulum.py @@ -68,7 +68,7 @@ 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 +93,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 69cdde787..f898fdea9 100644 --- a/bioptim/examples/stochastic_optimal_control/obstacle_avoidance_direct_collocation.py +++ b/bioptim/examples/stochastic_optimal_control/obstacle_avoidance_direct_collocation.py @@ -553,6 +553,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( @@ -564,14 +569,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, @@ -587,7 +588,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 a75440a22..c0ffd8379 100644 --- a/bioptim/examples/stochastic_optimal_control/rockit_matrix_lyapunov.py +++ b/bioptim/examples/stochastic_optimal_control/rockit_matrix_lyapunov.py @@ -194,6 +194,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( @@ -206,10 +208,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, @@ -219,7 +219,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..f9f1275ba 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,7 @@ 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 +109,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..5a24dfee1 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,7 @@ 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 +145,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 b38f6a2b4..9e0137a7c 100644 --- a/bioptim/examples/torque_driven_ocp/example_pendulum_time_dependent.py +++ b/bioptim/examples/torque_driven_ocp/example_pendulum_time_dependent.py @@ -151,7 +151,9 @@ 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 @@ -187,7 +189,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..1db624e24 100644 --- a/bioptim/examples/torque_driven_ocp/example_quaternions.py +++ b/bioptim/examples/torque_driven_ocp/example_quaternions.py @@ -207,7 +207,7 @@ 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 +252,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 018d16f66..1d625be28 100644 --- a/bioptim/examples/torque_driven_ocp/example_soft_contact.py +++ b/bioptim/examples/torque_driven_ocp/example_soft_contact.py @@ -52,6 +52,7 @@ def prepare_single_shooting( dynamics = Dynamics( DynamicsFcn.TORQUE_DRIVEN, soft_contacts_dynamics=SoftContactDynamics.ODE, + ode_solver=ode_solver, ) return OptimalControlProgram( @@ -59,7 +60,6 @@ def prepare_single_shooting( dynamics, n_shooting, final_time, - ode_solver=ode_solver, use_sx=use_sx, n_threads=n_threads, ) @@ -150,6 +150,7 @@ def prepare_ocp( dynamics = Dynamics( DynamicsFcn.TORQUE_DRIVEN, soft_contacts_dynamics=SoftContactDynamics.ODE, + ode_solver=ode_solver, phase_dynamics=phase_dynamics, ) @@ -191,7 +192,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 32b3c431d..b6a0684c6 100644 --- a/bioptim/examples/torque_driven_ocp/maximize_predicted_height_CoM.py +++ b/bioptim/examples/torque_driven_ocp/maximize_predicted_height_CoM.py @@ -96,6 +96,7 @@ def prepare_ocp( if use_actuators: dynamics.add( DynamicsFcn.TORQUE_ACTIVATIONS_DRIVEN, + ode_solver=ode_solver, with_contact=True, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics, @@ -104,6 +105,7 @@ def prepare_ocp( dynamics.add( DynamicsFcn.TORQUE_DRIVEN, with_contact=True, + ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics, ) @@ -155,7 +157,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..050fd2eda 100644 --- a/bioptim/examples/torque_driven_ocp/slider.py +++ b/bioptim/examples/torque_driven_ocp/slider.py @@ -76,9 +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, 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 +114,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..9a8ae2d67 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,7 @@ 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 +144,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..eabcc57f6 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,14 @@ 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 +123,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..fb08287a1 100644 --- a/bioptim/examples/torque_driven_ocp/trampo_quaternions.py +++ b/bioptim/examples/torque_driven_ocp/trampo_quaternions.py @@ -90,7 +90,7 @@ 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 +134,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..aeb9ba8cb 100644 --- a/bioptim/examples/track/track_marker_on_segment.py +++ b/bioptim/examples/track/track_marker_on_segment.py @@ -79,7 +79,7 @@ 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 +126,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..f9b103ab0 100644 --- a/bioptim/examples/track/track_segment_on_rt.py +++ b/bioptim/examples/track/track_segment_on_rt.py @@ -69,7 +69,7 @@ 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 +98,6 @@ def prepare_ocp( u_bounds=u_bounds, objective_functions=objective_functions, constraints=constraints, - ode_solver=ode_solver, ) diff --git a/bioptim/gui/serializable_class.py b/bioptim/gui/serializable_class.py index 7d4879b2d..132fdb59c 100644 --- a/bioptim/gui/serializable_class.py +++ b/bioptim/gui/serializable_class.py @@ -489,7 +489,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 +503,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 +518,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/optimization/non_linear_program.py b/bioptim/optimization/non_linear_program.py index f9320e59e..034b4a36c 100644 --- a/bioptim/optimization/non_linear_program.py +++ b/bioptim/optimization/non_linear_program.py @@ -280,13 +280,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: + if x_init is not None or a_init is not None: not_direct_collocation = not self.ode_solver.is_direct_collocation - init_all_point = x_init.type == InterpolationType.ALL_POINTS + 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, @@ -419,21 +419,23 @@ def n_states_decision_steps(self, node_idx) -> int: return 1 return self.dynamics[node_idx].shape_xf[1] + (1 if self.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.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 5325ea91a..11c7dd03e 100644 --- a/bioptim/optimization/optimal_control_program.py +++ b/bioptim/optimization/optimal_control_program.py @@ -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, @@ -216,8 +215,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 @@ -289,7 +286,6 @@ def __init__( parameter_init, parameter_constraints, parameter_objectives, - ode_solver, use_sx, bio_model, plot_mappings, @@ -432,7 +428,6 @@ def _check_arguments_and_build_nlp( parameter_init, parameter_constraints, parameter_objectives, - ode_solver, use_sx, bio_model, plot_mappings, @@ -509,18 +504,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") @@ -531,6 +514,14 @@ def _check_arguments_and_build_nlp( if not isinstance(dynamics, DynamicsList): raise ValueError("dynamics must be of type DynamicsList or Dynamics") + for i_dyn, dyn in enumerate(dynamics): + if dyn.ode_solver is None: + dynamics[i_dyn].ode_solver = self._set_default_ode_solver() + + is_ode_solver = isinstance(dynamics[i_dyn].ode_solver, OdeSolverBase) + if not is_ode_solver: + raise RuntimeError("ode_solver should be built an instance of OdeSolver") + # Type of CasADi graph self.cx = SX if use_sx else MX @@ -586,7 +577,8 @@ 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) + ode_solver = [dyn.ode_solver for dyn in dynamics] + NLP.add(self, "ode_solver", ode_solver, False) NLP.add(self, "control_type", control_type, True) # Prepare the variable mappings 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/stochastic_optimal_control_program.py b/bioptim/optimization/stochastic_optimal_control_program.py index 96a2c74c4..09bb82712 100644 --- a/bioptim/optimization/stochastic_optimal_control_program.py +++ b/bioptim/optimization/stochastic_optimal_control_program.py @@ -6,7 +6,7 @@ from .non_linear_program import NonLinearProgram as NLP from .optimization_vector import OptimizationVectorHelper from ..dynamics.configure_problem import DynamicsList, Dynamics -from ..dynamics.ode_solvers import OdeSolver +from ..dynamics.ode_solvers import OdeSolver, OdeSolverBase from ..limits.constraints import ( ConstraintFcn, ConstraintList, @@ -74,7 +74,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 @@ -87,6 +86,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( @@ -145,7 +148,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, @@ -621,22 +623,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 c3933cf5f..b0297e121 100644 --- a/bioptim/optimization/variational_optimal_control_program.py +++ b/bioptim/optimization/variational_optimal_control_program.py @@ -67,12 +67,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 +83,7 @@ def __init__( self.configure_torque_driven, expand_dynamics=expand, skip_continuity=True, + ode_solver=None, ) if qdot_bounds is None or not isinstance(qdot_bounds, BoundsList): From 915e47e611399722f0def635b45c7936bc59f55e Mon Sep 17 00:00:00 2001 From: Charbie Date: Tue, 1 Apr 2025 15:17:06 +0400 Subject: [PATCH 02/20] blacked --- bioptim/examples/acados/cube.py | 4 +--- .../getting_started/custom_constraint.py | 6 ++--- .../getting_started/custom_initial_guess.py | 7 +++--- .../getting_started/custom_objectives.py | 5 ++-- .../getting_started/custom_parameters.py | 4 +++- .../custom_phase_transitions.py | 16 +++++++++---- .../example_cyclic_movement.py | 4 +++- .../example_inequality_constraint.py | 6 ++++- .../example_multinode_constraints.py | 12 +++++++--- .../example_multinode_objective.py | 4 +++- .../getting_started/example_multiphase.py | 12 +++++++--- .../example_parameter_scaling.py | 4 +++- .../example_variable_scaling.py | 4 +++- bioptim/examples/getting_started/pendulum.py | 4 +++- .../pendulum_constrained_states_controls.py | 4 +++- .../multiphase_time_constraint.py | 24 ++++++++++++++++--- .../pendulum_min_time_Mayer.py | 4 +++- .../optimal_time_ocp/time_constraint.py | 4 +++- bioptim/examples/sqp_method/pendulum.py | 4 +++- .../symmetry_by_constraint.py | 4 +++- .../symmetry_by_mapping.py | 4 +++- .../example_pendulum_time_dependent.py | 4 +++- .../torque_driven_ocp/example_quaternions.py | 5 +++- bioptim/examples/torque_driven_ocp/slider.py | 12 +++++++--- .../track_markers_2D_pendulum.py | 4 +++- .../track_markers_with_torque_actuators.py | 19 ++++++++++++--- .../torque_driven_ocp/trampo_quaternions.py | 4 +++- .../examples/track/track_marker_on_segment.py | 4 +++- bioptim/examples/track/track_segment_on_rt.py | 4 +++- 29 files changed, 146 insertions(+), 50 deletions(-) diff --git a/bioptim/examples/acados/cube.py b/bioptim/examples/acados/cube.py index 30c32d60f..234310656 100644 --- a/bioptim/examples/acados/cube.py +++ b/bioptim/examples/acados/cube.py @@ -24,9 +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, - ode_solver=ode_solver, - expand_dynamics=expand_dynamics) + dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics) # Path constraint x_bounds = BoundsList() diff --git a/bioptim/examples/getting_started/custom_constraint.py b/bioptim/examples/getting_started/custom_constraint.py index 7f54968a9..3eef913ec 100644 --- a/bioptim/examples/getting_started/custom_constraint.py +++ b/bioptim/examples/getting_started/custom_constraint.py @@ -101,9 +101,9 @@ def prepare_ocp( objective_functions = Objective(ObjectiveFcn.Lagrange.MINIMIZE_CONTROL, key="tau", weight=100) # Dynamics - dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, - ode_solver=ode_solver, - 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() diff --git a/bioptim/examples/getting_started/custom_initial_guess.py b/bioptim/examples/getting_started/custom_initial_guess.py index 5a8d452e1..c98ad0c3a 100644 --- a/bioptim/examples/getting_started/custom_initial_guess.py +++ b/bioptim/examples/getting_started/custom_initial_guess.py @@ -126,10 +126,9 @@ def prepare_ocp( objective_functions = Objective(ObjectiveFcn.Lagrange.MINIMIZE_CONTROL, key="tau", weight=100) # Dynamics - dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, - ode_solver=ode_solver, - 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() diff --git a/bioptim/examples/getting_started/custom_objectives.py b/bioptim/examples/getting_started/custom_objectives.py index a0ff0178c..812d72de8 100644 --- a/bioptim/examples/getting_started/custom_objectives.py +++ b/bioptim/examples/getting_started/custom_objectives.py @@ -120,8 +120,9 @@ def prepare_ocp( ) # Dynamics - dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, - ode_solver=ode_solver, 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() diff --git a/bioptim/examples/getting_started/custom_parameters.py b/bioptim/examples/getting_started/custom_parameters.py index 6826d6df9..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, ode_solver=ode_solver, 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() diff --git a/bioptim/examples/getting_started/custom_phase_transitions.py b/bioptim/examples/getting_started/custom_phase_transitions.py index f9ffa8070..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, 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) + 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() diff --git a/bioptim/examples/getting_started/example_cyclic_movement.py b/bioptim/examples/getting_started/example_cyclic_movement.py index b09192bba..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, ode_solver=ode_solver, 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() diff --git a/bioptim/examples/getting_started/example_inequality_constraint.py b/bioptim/examples/getting_started/example_inequality_constraint.py index 21ae55620..5458340ec 100644 --- a/bioptim/examples/getting_started/example_inequality_constraint.py +++ b/bioptim/examples/getting_started/example_inequality_constraint.py @@ -93,7 +93,11 @@ def prepare_ocp( # Dynamics dynamics = DynamicsList() dynamics.add( - DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, with_contact=True, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + DynamicsFcn.TORQUE_DRIVEN, + ode_solver=ode_solver, + with_contact=True, + expand_dynamics=expand_dynamics, + phase_dynamics=phase_dynamics, ) # Constraints diff --git a/bioptim/examples/getting_started/example_multinode_constraints.py b/bioptim/examples/getting_started/example_multinode_constraints.py index a01e00108..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, 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 + ) + 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() diff --git a/bioptim/examples/getting_started/example_multinode_objective.py b/bioptim/examples/getting_started/example_multinode_objective.py index 83078bd59..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, ode_solver=ode_solver, 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() diff --git a/bioptim/examples/getting_started/example_multiphase.py b/bioptim/examples/getting_started/example_multiphase.py index 94b2b5be4..e0199e0af 100644 --- a/bioptim/examples/getting_started/example_multiphase.py +++ b/bioptim/examples/getting_started/example_multiphase.py @@ -120,9 +120,15 @@ def prepare_ocp( # Dynamics dynamics = DynamicsList() - 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 + ) + 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() diff --git a/bioptim/examples/getting_started/example_parameter_scaling.py b/bioptim/examples/getting_started/example_parameter_scaling.py index 1863484e5..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, ode_solver=ode_solver, 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() diff --git a/bioptim/examples/getting_started/example_variable_scaling.py b/bioptim/examples/getting_started/example_variable_scaling.py index a8020cd31..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, ode_solver=ode_solver, 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() diff --git a/bioptim/examples/getting_started/pendulum.py b/bioptim/examples/getting_started/pendulum.py index 036cb7b7d..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, ode_solver=ode_solver, 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() diff --git a/bioptim/examples/getting_started/pendulum_constrained_states_controls.py b/bioptim/examples/getting_started/pendulum_constrained_states_controls.py index e59ca982e..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, ode_solver=ode_solver, 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() diff --git a/bioptim/examples/optimal_time_ocp/multiphase_time_constraint.py b/bioptim/examples/optimal_time_ocp/multiphase_time_constraint.py index 359d1054c..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, ode_solver=ode_solver, 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, 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) + 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() 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 a23cd1743..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, 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 + ) # Path constraint x_bounds = BoundsList() diff --git a/bioptim/examples/optimal_time_ocp/time_constraint.py b/bioptim/examples/optimal_time_ocp/time_constraint.py index 846b0e486..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, ode_solver=ode_solver, 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) diff --git a/bioptim/examples/sqp_method/pendulum.py b/bioptim/examples/sqp_method/pendulum.py index 020e7de3b..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, ode_solver=ode_solver, 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() 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 f9f1275ba..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, 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() 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 5a24dfee1..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, 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() 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 9e0137a7c..3742907da 100644 --- a/bioptim/examples/torque_driven_ocp/example_pendulum_time_dependent.py +++ b/bioptim/examples/torque_driven_ocp/example_pendulum_time_dependent.py @@ -153,7 +153,9 @@ def prepare_ocp( dynamics.add( custom_configure, dynamic_function=time_dependent_dynamic, - ode_solver=ode_solver, expand_dynamics=expand, phase_dynamics=phase_dynamics + ode_solver=ode_solver, + expand_dynamics=expand, + phase_dynamics=phase_dynamics, ) # Path constraint diff --git a/bioptim/examples/torque_driven_ocp/example_quaternions.py b/bioptim/examples/torque_driven_ocp/example_quaternions.py index 1db624e24..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, ode_solver=ode_solver, 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 diff --git a/bioptim/examples/torque_driven_ocp/slider.py b/bioptim/examples/torque_driven_ocp/slider.py index 050fd2eda..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, 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 + ) + 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() 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 9a8ae2d67..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, 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 + ) # Path constraint x_bounds = BoundsList() 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 eabcc57f6..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, ode_solver=ode_solver, 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, 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, + ) else: raise ValueError("actuator_type is 1 (torque activations) or 2 (torque max constraints)") else: - 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() diff --git a/bioptim/examples/torque_driven_ocp/trampo_quaternions.py b/bioptim/examples/torque_driven_ocp/trampo_quaternions.py index fb08287a1..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, 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 + ) # Define control path constraint n_tau = bio_model.nb_tau # bio_model.nb_tau diff --git a/bioptim/examples/track/track_marker_on_segment.py b/bioptim/examples/track/track_marker_on_segment.py index aeb9ba8cb..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, 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 if constr: diff --git a/bioptim/examples/track/track_segment_on_rt.py b/bioptim/examples/track/track_segment_on_rt.py index f9b103ab0..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, 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() From fe34bdf08bd53b6316e0e2f86988816d3f04ac7c Mon Sep 17 00:00:00 2001 From: Charbie Date: Tue, 1 Apr 2025 16:05:04 +0400 Subject: [PATCH 03/20] forgot some ode_solvers --- bioptim/dynamics/configure_problem.py | 6 +++++- bioptim/examples/getting_started/example_multiphase.py | 8 +++++--- tests/shard1/test_controltype_none.py | 2 +- tests/shard1/test_dynamics.py | 3 +-- tests/shard2/test_cost_function_integration.py | 3 +-- tests/shard2/test_global_minimize_marker_velocity.py | 2 +- tests/shard3/test_graph.py | 6 ++---- tests/shard3/test_multinode_constraints.py | 7 +++---- tests/shard4/test_update_bounds_and_init.py | 8 +++----- tests/shard4/test_variable_time.py | 7 +++---- tests/shard5/test_helper_socp.py | 10 ---------- tests/shard6/test_dt_dependent.py | 2 +- tests/shard6/test_time_dependent_custom_model.py | 2 +- tests/shard6/test_time_dependent_problems.py | 2 +- 14 files changed, 28 insertions(+), 40 deletions(-) diff --git a/bioptim/dynamics/configure_problem.py b/bioptim/dynamics/configure_problem.py index bdbc74849..fec2e3b7a 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 ( @@ -1955,6 +1955,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, **extra_parameters: Any, ): @@ -1976,6 +1977,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. """ @@ -2002,6 +2005,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 diff --git a/bioptim/examples/getting_started/example_multiphase.py b/bioptim/examples/getting_started/example_multiphase.py index e0199e0af..cc8f1dc5a 100644 --- a/bioptim/examples/getting_started/example_multiphase.py +++ b/bioptim/examples/getting_started/example_multiphase.py @@ -119,15 +119,17 @@ def prepare_ocp( ) # Dynamics + if not isinstance(ode_solver, list): + ode_solver = [ode_solver] * 3 dynamics = DynamicsList() dynamics.add( - DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + 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, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + 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, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver[2], expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics ) # Constraints diff --git a/tests/shard1/test_controltype_none.py b/tests/shard1/test_controltype_none.py index 0afc6a40d..50036a390 100644 --- a/tests/shard1/test_controltype_none.py +++ b/tests/shard1/test_controltype_none.py @@ -191,6 +191,7 @@ def prepare_ocp( dynamic_function=custom_model.custom_dynamics, phase=i, expand_dynamics=True, + ode_solver=ode_solver, phase_dynamics=phase_dynamics, ) @@ -224,7 +225,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 eb9868c65..021b2f08b 100644 --- a/tests/shard1/test_dynamics.py +++ b/tests/shard1/test_dynamics.py @@ -1231,7 +1231,7 @@ def test_with_contact_error(dynamics_fcn, phase_dynamics): objective_functions = ObjectiveList() # Dynamics - dynamics = Dynamics(dynamics_fcn, with_contact=True, expand_dynamics=True, phase_dynamics=phase_dynamics) + dynamics = Dynamics(dynamics_fcn, with_contact=True, ode_solver=OdeSolver.RK4(), expand_dynamics=True, phase_dynamics=phase_dynamics) # Path constraint x_bounds = BoundsList() @@ -1253,5 +1253,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/shard2/test_cost_function_integration.py b/tests/shard2/test_cost_function_integration.py index c6d52d418..e44ab6d5a 100644 --- a/tests/shard2/test_cost_function_integration.py +++ b/tests/shard2/test_cost_function_integration.py @@ -87,7 +87,7 @@ 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 +112,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 d2a1f46ee..711b1a537 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/shard3/test_graph.py b/tests/shard3/test_graph.py index cde4014f4..258923a94 100644 --- a/tests/shard3/test_graph.py +++ b/tests/shard3/test_graph.py @@ -322,7 +322,7 @@ 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 +350,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 +416,7 @@ 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 +440,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_multinode_constraints.py b/tests/shard3/test_multinode_constraints.py index 565ea5829..831155b2c 100644 --- a/tests/shard3/test_multinode_constraints.py +++ b/tests/shard3/test_multinode_constraints.py @@ -30,9 +30,9 @@ 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 +82,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/shard4/test_update_bounds_and_init.py b/tests/shard4/test_update_bounds_and_init.py index 094727b00..9614b32d3 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, ) diff --git a/tests/shard4/test_variable_time.py b/tests/shard4/test_variable_time.py index 810538699..b6aa0a21f 100644 --- a/tests/shard4/test_variable_time.py +++ b/tests/shard4/test_variable_time.py @@ -54,9 +54,9 @@ 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 +134,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_helper_socp.py b/tests/shard5/test_helper_socp.py index bc0eda5d5..9982aff8e 100644 --- a/tests/shard5/test_helper_socp.py +++ b/tests/shard5/test_helper_socp.py @@ -26,16 +26,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_dt_dependent.py b/tests/shard6/test_dt_dependent.py index daed7a46d..92fb4909b 100644 --- a/tests/shard6/test_dt_dependent.py +++ b/tests/shard6/test_dt_dependent.py @@ -95,6 +95,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, ) @@ -131,7 +132,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 97724d80c..84cfb713a 100644 --- a/tests/shard6/test_time_dependent_custom_model.py +++ b/tests/shard6/test_time_dependent_custom_model.py @@ -296,6 +296,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, ) @@ -376,7 +377,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 b85c773d1..b43354116 100644 --- a/tests/shard6/test_time_dependent_problems.py +++ b/tests/shard6/test_time_dependent_problems.py @@ -164,6 +164,7 @@ def prepare_ocp( custom_configure, dynamic_function=time_dynamic, phase=i, + ode_solver=ode_solver, expand_dynamics=expand, phase_dynamics=phase_dynamics, ) @@ -210,7 +211,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, ) From 3044959f15c79ed7b9cf9292f451aaa5cdbf9dc4 Mon Sep 17 00:00:00 2001 From: Charbie Date: Tue, 1 Apr 2025 16:05:19 +0400 Subject: [PATCH 04/20] blacked --- .../getting_started/example_multiphase.py | 15 ++++++++++++--- tests/shard1/test_dynamics.py | 4 +++- tests/shard2/test_cost_function_integration.py | 4 +++- tests/shard3/test_graph.py | 8 ++++++-- tests/shard3/test_multinode_constraints.py | 12 +++++++++--- tests/shard4/test_update_bounds_and_init.py | 2 +- tests/shard4/test_variable_time.py | 12 +++++++++--- 7 files changed, 43 insertions(+), 14 deletions(-) diff --git a/bioptim/examples/getting_started/example_multiphase.py b/bioptim/examples/getting_started/example_multiphase.py index cc8f1dc5a..d7f0f1dc6 100644 --- a/bioptim/examples/getting_started/example_multiphase.py +++ b/bioptim/examples/getting_started/example_multiphase.py @@ -123,13 +123,22 @@ def prepare_ocp( ode_solver = [ode_solver] * 3 dynamics = DynamicsList() dynamics.add( - DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver[0], expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics + 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 + 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 + DynamicsFcn.TORQUE_DRIVEN, + ode_solver=ode_solver[2], + expand_dynamics=expand_dynamics, + phase_dynamics=phase_dynamics, ) # Constraints diff --git a/tests/shard1/test_dynamics.py b/tests/shard1/test_dynamics.py index 021b2f08b..1f43d1ff7 100644 --- a/tests/shard1/test_dynamics.py +++ b/tests/shard1/test_dynamics.py @@ -1231,7 +1231,9 @@ def test_with_contact_error(dynamics_fcn, phase_dynamics): objective_functions = ObjectiveList() # Dynamics - dynamics = Dynamics(dynamics_fcn, with_contact=True, ode_solver=OdeSolver.RK4(), expand_dynamics=True, phase_dynamics=phase_dynamics) + dynamics = Dynamics( + dynamics_fcn, with_contact=True, ode_solver=OdeSolver.RK4(), expand_dynamics=True, phase_dynamics=phase_dynamics + ) # Path constraint x_bounds = BoundsList() diff --git a/tests/shard2/test_cost_function_integration.py b/tests/shard2/test_cost_function_integration.py index e44ab6d5a..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, ode_solver=ode_solver, 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() diff --git a/tests/shard3/test_graph.py b/tests/shard3/test_graph.py index 258923a94..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, ode_solver=ode_solver, 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() @@ -416,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, ode_solver=ode_solver, 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() diff --git a/tests/shard3/test_multinode_constraints.py b/tests/shard3/test_multinode_constraints.py index 831155b2c..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, 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) + 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 diff --git a/tests/shard4/test_update_bounds_and_init.py b/tests/shard4/test_update_bounds_and_init.py index 9614b32d3..03ccee621 100644 --- a/tests/shard4/test_update_bounds_and_init.py +++ b/tests/shard4/test_update_bounds_and_init.py @@ -981,7 +981,7 @@ def test_update_noised_initial_guess_collocation(interpolation, phase_dynamics): ode_solver = OdeSolver.COLLOCATION(polynomial_degree=1) dynamics = DynamicsList() - dynamics.add(DynamicsFcn.TORQUE_DRIVEN, ode_solver=ode_solver,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 diff --git a/tests/shard4/test_variable_time.py b/tests/shard4/test_variable_time.py index b6aa0a21f..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, 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) + 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() From c5b6d2b44c62f822254d2c22f5a4701ae278f8f6 Mon Sep 17 00:00:00 2001 From: Charbie Date: Tue, 1 Apr 2025 16:09:47 +0400 Subject: [PATCH 05/20] . --- tests/shard5/test_helper_socp.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/shard5/test_helper_socp.py b/tests/shard5/test_helper_socp.py index 9982aff8e..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, ) From 3bf44e1a32433668c534037635b160e2dd2af289 Mon Sep 17 00:00:00 2001 From: eve-mac Date: Tue, 1 Apr 2025 18:24:50 +0400 Subject: [PATCH 06/20] fixing tests --- tests/shard1/test_prepare_all_examples.py | 2 +- tests/shard4/test_update_bounds_and_init.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/shard1/test_prepare_all_examples.py b/tests/shard1/test_prepare_all_examples.py index 02d36edf8..4b36281dc 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/shard4/test_update_bounds_and_init.py b/tests/shard4/test_update_bounds_and_init.py index 03ccee621..2ed9b8c42 100644 --- a/tests/shard4/test_update_bounds_and_init.py +++ b/tests/shard4/test_update_bounds_and_init.py @@ -1035,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) From a7e2872bee99b208f45078b87d22d95cc1303f5e Mon Sep 17 00:00:00 2001 From: Charbie Date: Mon, 7 Apr 2025 10:17:53 +0400 Subject: [PATCH 07/20] nlp.ode_solver -> nlp.dynamics_type.ode_solver --- README.md | 1 - bioptim/dynamics/configure_new_variable.py | 6 ++--- bioptim/dynamics/configure_problem.py | 6 +++-- bioptim/dynamics/dynamics_functions.py | 2 +- bioptim/dynamics/ode_solver_base.py | 10 ++++---- .../muscle_activations_tracker.py | 2 +- .../muscle_excitations_tracker.py | 2 +- bioptim/gui/graph.py | 4 +-- bioptim/gui/plot.py | 4 +-- bioptim/gui/serializable_class.py | 1 - bioptim/limits/penalty.py | 4 +-- bioptim/limits/penalty_controller.py | 2 +- bioptim/limits/penalty_option.py | 4 +-- bioptim/optimization/non_linear_program.py | 15 +++++------ .../optimization/optimal_control_program.py | 25 +++---------------- bioptim/optimization/solution/solution.py | 12 ++++----- .../stochastic_optimal_control_program.py | 2 +- tests/shard2/test_global_nmpc_final.py | 4 +-- .../test_global_stochastic_collocation.py | 2 +- tests/utils.py | 4 +-- 20 files changed, 46 insertions(+), 66 deletions(-) diff --git a/README.md b/README.md index a003bde57..c0a6b60d9 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, diff --git a/bioptim/dynamics/configure_new_variable.py b/bioptim/dynamics/configure_new_variable.py index 30bd04966..bc6320410 100644 --- a/bioptim/dynamics/configure_new_variable.py +++ b/bioptim/dynamics/configure_new_variable.py @@ -220,7 +220,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( @@ -280,7 +280,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.xdot_scaling[self.name].scaling) self.nlp.states_dot.append( @@ -295,7 +295,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 fec2e3b7a..8e4c867e0 100644 --- a/bioptim/dynamics/configure_problem.py +++ b/bioptim/dynamics/configure_problem.py @@ -1341,7 +1341,7 @@ 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 @@ -1997,6 +1997,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 @@ -2032,7 +2035,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 cf1339c78..dff2b271e 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 not with_contact and fatigue is None: qddot = DynamicsFunctions.get(nlp.states_dot["qddot"], nlp.states_dot.scaled.cx) tau_id = DynamicsFunctions.inverse_dynamics(nlp, q, qdot, qddot, with_contact, external_forces) diff --git a/bioptim/dynamics/ode_solver_base.py b/bioptim/dynamics/ode_solver_base.py index d783a4f06..36b3f9b22 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,26 @@ 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/examples/muscle_driven_ocp/muscle_activations_tracker.py b/bioptim/examples/muscle_driven_ocp/muscle_activations_tracker.py index d44e82acf..d5c472d32 100644 --- a/bioptim/examples/muscle_driven_ocp/muscle_activations_tracker.py +++ b/bioptim/examples/muscle_driven_ocp/muscle_activations_tracker.py @@ -369,7 +369,7 @@ 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 76f3b06d6..55e457922 100644 --- a/bioptim/examples/muscle_driven_ocp/muscle_excitations_tracker.py +++ b/bioptim/examples/muscle_driven_ocp/muscle_excitations_tracker.py @@ -372,7 +372,7 @@ 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/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 132fdb59c..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 diff --git a/bioptim/limits/penalty.py b/bioptim/limits/penalty.py index 68de6df6b..2fb684408 100644 --- a/bioptim/limits/penalty.py +++ b/bioptim/limits/penalty.py @@ -920,7 +920,7 @@ def minimize_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( @@ -1187,7 +1187,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: cx = horzcat(*([controller.states.cx_start] + controller.states.cx_intermediates_list)) integrated = controller.integrate( 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 034b4a36c..41bf1baf8 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,10 @@ 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): @@ -281,7 +278,7 @@ def _update_plot_bounds(self, keys, bounds, suffix): def update_init(self, x_init, u_init, a_init): if x_init is not None or a_init is not None: - not_direct_collocation = not self.ode_solver.is_direct_collocation + 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 @@ -417,7 +414,7 @@ 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, ode_solver: OdeSolver = None) -> int: """ @@ -431,7 +428,7 @@ def n_states_stepwise_steps(self, node_idx: int, ode_solver: OdeSolver = None) - ------- The number of states """ - ode_solver = ode_solver if ode_solver is not None else self.ode_solver + ode_solver = ode_solver if ode_solver is not None else self.dynamics_type.ode_solver if node_idx >= self.ns: return 1 if ode_solver.is_direct_collocation: diff --git a/bioptim/optimization/optimal_control_program.py b/bioptim/optimization/optimal_control_program.py index 11c7dd03e..205f6c646 100644 --- a/bioptim/optimization/optimal_control_program.py +++ b/bioptim/optimization/optimal_control_program.py @@ -10,7 +10,6 @@ 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 ..gui.check_conditioning import check_conditioning from ..gui.graph import OcpToConsole, OcpToGraph from ..gui.ipopt_output_plot import SaveIterationsInfo @@ -514,14 +513,6 @@ def _check_arguments_and_build_nlp( if not isinstance(dynamics, DynamicsList): raise ValueError("dynamics must be of type DynamicsList or Dynamics") - for i_dyn, dyn in enumerate(dynamics): - if dyn.ode_solver is None: - dynamics[i_dyn].ode_solver = self._set_default_ode_solver() - - is_ode_solver = isinstance(dynamics[i_dyn].ode_solver, OdeSolverBase) - if not is_ode_solver: - raise RuntimeError("ode_solver should be built an instance of OdeSolver") - # Type of CasADi graph self.cx = SX if use_sx else MX @@ -577,8 +568,6 @@ def _check_arguments_and_build_nlp( # Prepare path constraints and dynamics of the program NLP.add(self, "dynamics_type", dynamics, False) - ode_solver = [dyn.ode_solver for dyn in dynamics] - NLP.add(self, "ode_solver", ode_solver, False) NLP.add(self, "control_type", control_type, True) # Prepare the variable mappings @@ -615,7 +604,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, " @@ -861,7 +850,7 @@ 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, @@ -1434,7 +1423,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: @@ -1681,16 +1670,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/solution/solution.py b/bioptim/optimization/solution/solution.py index d46e18f61..1488de495 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,7 @@ 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 09bb82712..1ac90f0a3 100644 --- a/bioptim/optimization/stochastic_optimal_control_program.py +++ b/bioptim/optimization/stochastic_optimal_control_program.py @@ -6,7 +6,7 @@ from .non_linear_program import NonLinearProgram as NLP from .optimization_vector import OptimizationVectorHelper from ..dynamics.configure_problem import DynamicsList, Dynamics -from ..dynamics.ode_solvers import OdeSolver, OdeSolverBase +from ..dynamics.ode_solvers import OdeSolver from ..limits.constraints import ( ConstraintFcn, ConstraintList, 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/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/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, " From 8ed6ac281dcd233043514bc682ec7ed4e71c804d Mon Sep 17 00:00:00 2001 From: Charbie Date: Mon, 7 Apr 2025 10:18:06 +0400 Subject: [PATCH 08/20] blacked --- bioptim/dynamics/configure_problem.py | 6 +++++- bioptim/dynamics/ode_solver_base.py | 10 ++++++++-- .../muscle_driven_ocp/muscle_activations_tracker.py | 6 +++++- .../muscle_driven_ocp/muscle_excitations_tracker.py | 6 +++++- bioptim/optimization/non_linear_program.py | 12 +++++++++--- bioptim/optimization/optimal_control_program.py | 5 ++++- 6 files changed, 36 insertions(+), 9 deletions(-) diff --git a/bioptim/dynamics/configure_problem.py b/bioptim/dynamics/configure_problem.py index 8e4c867e0..5547be747 100644 --- a/bioptim/dynamics/configure_problem.py +++ b/bioptim/dynamics/configure_problem.py @@ -1341,7 +1341,11 @@ def configure_integrated_value( # TODO: compute values at collocation points # but for now only cx_start can be used - n_cx = nlp.dynamics_type.ode_solver.n_cx - 1 if isinstance(nlp.dynamics_type.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 diff --git a/bioptim/dynamics/ode_solver_base.py b/bioptim/dynamics/ode_solver_base.py index 36b3f9b22..c271de288 100644 --- a/bioptim/dynamics/ode_solver_base.py +++ b/bioptim/dynamics/ode_solver_base.py @@ -249,14 +249,20 @@ def prepare_dynamic_integrator(self, ocp, nlp): dynamics = dynamics * nlp.ns else: for node_index in range(1, nlp.ns): - dynamics.append(nlp.dynamics_type.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.dynamics_type.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 diff --git a/bioptim/examples/muscle_driven_ocp/muscle_activations_tracker.py b/bioptim/examples/muscle_driven_ocp/muscle_activations_tracker.py index d5c472d32..c409b4638 100644 --- a/bioptim/examples/muscle_driven_ocp/muscle_activations_tracker.py +++ b/bioptim/examples/muscle_driven_ocp/muscle_activations_tracker.py @@ -369,7 +369,11 @@ def main(): markers[:, :, i] = bio_model.markers()(q[:, i]) plt.figure("Markers") - 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 + 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 55e457922..8b04d2037 100644 --- a/bioptim/examples/muscle_driven_ocp/muscle_excitations_tracker.py +++ b/bioptim/examples/muscle_driven_ocp/muscle_excitations_tracker.py @@ -372,7 +372,11 @@ def main(): markers[:, :, i] = horzcat(*bio_model.markers()(q[:, i])) plt.figure("Markers") - 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 + 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/optimization/non_linear_program.py b/bioptim/optimization/non_linear_program.py index 41bf1baf8..491418251 100644 --- a/bioptim/optimization/non_linear_program.py +++ b/bioptim/optimization/non_linear_program.py @@ -207,9 +207,13 @@ def initialize(self, cx: MX | SX | Callable = None): 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.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.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.dynamics_type.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): @@ -414,7 +418,9 @@ 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.dynamics_type.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, ode_solver: OdeSolver = None) -> int: """ diff --git a/bioptim/optimization/optimal_control_program.py b/bioptim/optimization/optimal_control_program.py index 205f6c646..80b40d6f5 100644 --- a/bioptim/optimization/optimal_control_program.py +++ b/bioptim/optimization/optimal_control_program.py @@ -850,7 +850,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.dynamics_type.ode_solver.is_direct_collocation and nlp.dynamics_type.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, From 01f3cd1dc3025e924d5b5c1bc5214b9928af083f Mon Sep 17 00:00:00 2001 From: Charbie Date: Mon, 7 Apr 2025 10:20:01 +0400 Subject: [PATCH 09/20] reblack --- bioptim/optimization/solution/solution.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/bioptim/optimization/solution/solution.py b/bioptim/optimization/solution/solution.py index 1488de495..290cdda97 100644 --- a/bioptim/optimization/solution/solution.py +++ b/bioptim/optimization/solution/solution.py @@ -731,7 +731,9 @@ def _prepare_integrate(self, integrator: SolutionIntegrator): " Shooting.SINGLE, Shooting.MULTIPLE, or Shooting.SINGLE_DISCONTINUOUS_PHASE" ) - has_trapezoidal = sum([isinstance(nlp.dynamics_type.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, " From 27f19b7e1947e7e2ef3142630fe215b93ccb4085 Mon Sep 17 00:00:00 2001 From: Charbie Date: Mon, 7 Apr 2025 11:18:03 +0400 Subject: [PATCH 10/20] blacked --- bioptim/dynamics/ode_solvers.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/bioptim/dynamics/ode_solvers.py b/bioptim/dynamics/ode_solvers.py index 9c360f23c..3c4b426b6 100644 --- a/bioptim/dynamics/ode_solvers.py +++ b/bioptim/dynamics/ode_solvers.py @@ -165,8 +165,8 @@ def x_ode(self, nlp): def p_ode(self, nlp): if nlp.control_type in ( - ControlType.CONSTANT, - ControlType.CONSTANT_WITH_LAST_NODE, + ControlType.CONSTANT, + ControlType.CONSTANT_WITH_LAST_NODE, ): return nlp.controls.scaled.cx_start elif nlp.control_type == ControlType.LINEAR_CONTINUOUS: From f825c0d3745ba6afd1a6ab0b749b67abada83b60 Mon Sep 17 00:00:00 2001 From: Charbie Date: Mon, 7 Apr 2025 13:48:40 +0400 Subject: [PATCH 11/20] dynamics_type declaration --- tests/shard1/test_dynamics.py | 149 ++++++++++++++++------------ tests/shard3/test_ligaments.py | 23 +++-- tests/shard3/test_passive_torque.py | 54 ++++++---- 3 files changed, 132 insertions(+), 94 deletions(-) diff --git a/tests/shard1/test_dynamics.py b/tests/shard1/test_dynamics.py index 1f43d1ff7..17122ee96 100644 --- a/tests/shard1/test_dynamics.py +++ b/tests/shard1/test_dynamics.py @@ -121,6 +121,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, + with_contact=with_contact, + expand_dynamics=True, + phase_dynamics=phase_dynamics, + numerical_data_timeseries=numerical_time_series, + ) nlp.cx = cx nlp.time_cx = cx.sym("time", 1, 1) @@ -139,13 +146,7 @@ def test_torque_driven(with_contact, with_external_force, cx, phase_dynamics): NonLinearProgram.add( ocp, "dynamics_type", - Dynamics( - DynamicsFcn.TORQUE_DRIVEN, - with_contact=with_contact, - 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)] @@ -220,6 +221,14 @@ def test_torque_driven_soft_contacts_dynamics(with_contact, cx, implicit_contact nlp.model = BiorbdModel( TestUtils.bioptim_folder() + "/examples/getting_started/models/2segments_4dof_2contacts.bioMod" ) + nlp.dynamics_type = Dynamics( + DynamicsFcn.TORQUE_DRIVEN, + with_contact=with_contact, + soft_contacts_dynamics=implicit_contact, + expand_dynamics=True, + phase_dynamics=phase_dynamics, + ) + nlp.ns = N_SHOOTING nlp.cx = cx nlp.time_cx = cx.sym("time", 1, 1) @@ -239,13 +248,7 @@ def test_torque_driven_soft_contacts_dynamics(with_contact, cx, implicit_contact NonLinearProgram.add( ocp, "dynamics_type", - Dynamics( - DynamicsFcn.TORQUE_DRIVEN, - with_contact=with_contact, - soft_contacts_dynamics=implicit_contact, - expand_dynamics=True, - phase_dynamics=phase_dynamics, - ), + nlp.dynamics_type, False, ) @@ -303,6 +306,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, + with_contact=with_contact, + 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) @@ -320,13 +331,7 @@ def test_torque_derivative_driven(with_contact, with_external_force, cx, phase_d NonLinearProgram.add( ocp, "dynamics_type", - Dynamics( - DynamicsFcn.TORQUE_DERIVATIVE_DRIVEN, - with_contact=with_contact, - expand_dynamics=True, - phase_dynamics=phase_dynamics, - numerical_data_timeseries=numerical_timeseries, - ), + nlp.dynamics_type, False, ) @@ -445,6 +450,14 @@ def test_torque_derivative_driven_soft_contacts_dynamics(with_contact, cx, impli nlp.model = BiorbdModel( TestUtils.bioptim_folder() + "/examples/getting_started/models/2segments_4dof_2contacts.bioMod" ) + nlp.dynamics_type = Dynamics( + DynamicsFcn.TORQUE_DERIVATIVE_DRIVEN, + with_contact=with_contact, + soft_contacts_dynamics=implicit_contact, + expand_dynamics=True, + phase_dynamics=phase_dynamics, + ) + nlp.ns = N_SHOOTING nlp.cx = cx nlp.time_cx = cx.sym("time", 1, 1) @@ -463,13 +476,7 @@ def test_torque_derivative_driven_soft_contacts_dynamics(with_contact, cx, impli NonLinearProgram.add( ocp, "dynamics_type", - Dynamics( - DynamicsFcn.TORQUE_DERIVATIVE_DRIVEN, - with_contact=with_contact, - soft_contacts_dynamics=implicit_contact, - expand_dynamics=True, - phase_dynamics=phase_dynamics, - ), + nlp.dynamics_type, False, ) @@ -544,6 +551,8 @@ 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 @@ -555,7 +564,7 @@ def test_soft_contacts_dynamics_errors(dynamics, phase_dynamics): NonLinearProgram.add( ocp, "dynamics_type", - Dynamics(dynamics, soft_contacts_dynamics=True, expand_dynamics=True, phase_dynamics=phase_dynamics), + nlp.dynamics_type, False, ) phase_index = [i for i in range(ocp.n_phases)] @@ -589,6 +598,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, + with_contact=with_contact, + 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) @@ -606,13 +623,7 @@ def test_torque_activation_driven(with_contact, with_external_force, cx, phase_d NonLinearProgram.add( ocp, "dynamics_type", - Dynamics( - DynamicsFcn.TORQUE_ACTIVATIONS_DRIVEN, - with_contact=with_contact, - expand_dynamics=True, - phase_dynamics=phase_dynamics, - numerical_data_timeseries=numerical_timeseries, - ), + nlp.dynamics_type, False, ) phase_index = [i for i in range(ocp.n_phases)] @@ -713,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) @@ -729,13 +748,7 @@ def test_torque_activation_driven_with_residual_torque( NonLinearProgram.add( ocp, "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.dynamics_type, False, ) phase_index = [i for i in range(ocp.n_phases)] @@ -826,6 +839,8 @@ 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) @@ -844,7 +859,7 @@ def test_torque_driven_free_floating_base(cx, phase_dynamics): NonLinearProgram.add( ocp, "dynamics_type", - Dynamics(DynamicsFcn.TORQUE_DRIVEN_FREE_FLOATING_BASE, expand_dynamics=True, phase_dynamics=phase_dynamics), + nlp.dynamics_type, False, ) phase_index = [i for i in range(ocp.n_phases)] @@ -897,6 +912,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, + with_contact=with_contact, + expand_dynamics=True, + phase_dynamics=phase_dynamics, + numerical_data_timeseries=numerical_timeseries, + ) nlp.cx = cx nlp.time_cx = cx.sym("time", 1, 1) @@ -916,15 +940,7 @@ def test_muscle_driven(with_excitations, with_contact, with_residual_torque, wit NonLinearProgram.add( ocp, "dynamics_type", - Dynamics( - DynamicsFcn.MUSCLE_DRIVEN, - with_residual_torque=with_residual_torque, - with_excitations=with_excitations, - with_contact=with_contact, - expand_dynamics=True, - phase_dynamics=phase_dynamics, - numerical_data_timeseries=numerical_timeseries, - ), + nlp.dynamics_type, False, ) phase_index = [i for i in range(ocp.n_phases)] @@ -1071,6 +1087,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,11 +1112,7 @@ def test_joints_acceleration_driven(cx, phase_dynamics): NonLinearProgram.add( ocp, "dynamics_type", - Dynamics( - DynamicsFcn.JOINTS_ACCELERATION_DRIVEN, - expand_dynamics=True, - phase_dynamics=phase_dynamics, - ), + nlp.dynamics_type, False, ) np.random.seed(42) @@ -1148,6 +1165,14 @@ def configure(ocp, nlp, with_contact=None, numerical_data_timeseries=None): nlp.model = BiorbdModel( TestUtils.bioptim_folder() + "/examples/getting_started/models/2segments_4dof_2contacts.bioMod" ) + nlp.dynamics_type = Dynamics( + configure, + dynamic_function=custom_dynamic, + with_contact=with_contact, + expand_dynamics=True, + phase_dynamics=phase_dynamics, + ) + nlp.ns = N_SHOOTING nlp.cx = MX nlp.time_cx = nlp.cx.sym("time", 1, 1) @@ -1165,13 +1190,7 @@ def configure(ocp, nlp, with_contact=None, numerical_data_timeseries=None): NonLinearProgram.add( ocp, "dynamics_type", - Dynamics( - configure, - dynamic_function=custom_dynamic, - with_contact=with_contact, - expand_dynamics=True, - phase_dynamics=phase_dynamics, - ), + nlp.dynamics_type, False, ) phase_index = [i for i in range(ocp.n_phases)] diff --git a/tests/shard3/test_ligaments.py b/tests/shard3/test_ligaments.py index 9fe8c161e..9f32bcf33 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) @@ -58,10 +60,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)] @@ -102,6 +105,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) @@ -120,7 +125,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, ) @@ -162,6 +167,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) @@ -178,7 +185,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)] @@ -220,6 +227,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) @@ -237,10 +249,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_passive_torque.py b/tests/shard3/test_passive_torque.py index c679ab83c..3c428a1f0 100644 --- a/tests/shard3/test_passive_torque.py +++ b/tests/shard3/test_passive_torque.py @@ -42,11 +42,16 @@ 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) nlp.dt = cx.sym("dt", 1, 1) - nlp.initialize(cx) nlp.x_scaling = VariableScalingList() nlp.xdot_scaling = VariableScalingList() nlp.u_scaling = VariableScalingList() @@ -56,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)] @@ -103,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) @@ -121,11 +129,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, ) @@ -194,6 +198,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) @@ -206,15 +217,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)] @@ -306,6 +313,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) @@ -320,14 +333,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)] From ad19ef9fd8c55e32def1c5da72cdfe918caac2a4 Mon Sep 17 00:00:00 2001 From: Charbie Date: Mon, 7 Apr 2025 15:31:20 +0400 Subject: [PATCH 12/20] fake ode_solver for variational --- bioptim/dynamics/integrator.py | 7 +++++++ bioptim/dynamics/ode_solvers.py | 10 ++++++++++ bioptim/optimization/optimal_control_program.py | 1 + .../variational_optimal_control_program.py | 3 ++- tests/shard3/test_passive_torque.py | 1 + 5 files changed, 21 insertions(+), 1 deletion(-) diff --git a/bioptim/dynamics/integrator.py b/bioptim/dynamics/integrator.py index 399abbd9e..3dd4d5b2d 100644 --- a/bioptim/dynamics/integrator.py +++ b/bioptim/dynamics/integrator.py @@ -785,3 +785,10 @@ class CVODES(Integrator): Class for CVODES integrators """ + + +class VARIATIONAL(RK4): + """ + Fake class for variational integrator. + The real work is done in VariationalOptimalControlProgram. + """ diff --git a/bioptim/dynamics/ode_solvers.py b/bioptim/dynamics/ode_solvers.py index 3c4b426b6..0ad0934f5 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. + """ + + @property + def integrator(self): + return integrator.VARIATIONAL + + class TRAPEZOIDAL(OdeSolverBase): """ A trapezoidal ode solver diff --git a/bioptim/optimization/optimal_control_program.py b/bioptim/optimization/optimal_control_program.py index 80b40d6f5..b1da0b657 100644 --- a/bioptim/optimization/optimal_control_program.py +++ b/bioptim/optimization/optimal_control_program.py @@ -10,6 +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 from ..gui.check_conditioning import check_conditioning from ..gui.graph import OcpToConsole, OcpToGraph from ..gui.ipopt_output_plot import SaveIterationsInfo diff --git a/bioptim/optimization/variational_optimal_control_program.py b/bioptim/optimization/variational_optimal_control_program.py index b0297e121..96bd847c3 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 @@ -83,7 +84,7 @@ def __init__( self.configure_torque_driven, expand_dynamics=expand, skip_continuity=True, - ode_solver=None, + 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/shard3/test_passive_torque.py b/tests/shard3/test_passive_torque.py index 3c428a1f0..083cc6bcb 100644 --- a/tests/shard3/test_passive_torque.py +++ b/tests/shard3/test_passive_torque.py @@ -52,6 +52,7 @@ def test_torque_driven_with_passive_torque(with_passive_torque, cx, phase_dynami nlp.cx = cx nlp.time_cx = cx.sym("time", 1, 1) nlp.dt = cx.sym("dt", 1, 1) + nlp.initialize(cx) nlp.x_scaling = VariableScalingList() nlp.xdot_scaling = VariableScalingList() nlp.u_scaling = VariableScalingList() From b967964c7fbe5c68a7212268e28e65b9620529c4 Mon Sep 17 00:00:00 2001 From: Charbie Date: Mon, 7 Apr 2025 15:38:35 +0400 Subject: [PATCH 13/20] tests --- tests/shard6/test_configure_problem.py | 40 ++++++++++++-------------- 1 file changed, 19 insertions(+), 21 deletions(-) diff --git a/tests/shard6/test_configure_problem.py b/tests/shard6/test_configure_problem.py index 38674582a..c1920c049 100644 --- a/tests/shard6/test_configure_problem.py +++ b/tests/shard6/test_configure_problem.py @@ -13,12 +13,21 @@ 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.xdot_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 +37,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,13 +56,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.xdot_scaling = VariableScalingList() - nlp.u_scaling = VariableScalingList() - nlp.a_scaling = VariableScalingList() ocp = OptimalControlProgram(nlp, use_sx=(cx == SX)) nlp.model = BiorbdModel( @@ -108,13 +120,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.xdot_scaling = VariableScalingList() - nlp.u_scaling = VariableScalingList() - nlp.a_scaling = VariableScalingList() ocp = OptimalControlProgram(nlp, use_sx=(cx == SX)) nlp.model = BiorbdModel( @@ -133,13 +138,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.xdot_scaling = VariableScalingList() - nlp.u_scaling = VariableScalingList() - nlp.a_scaling = VariableScalingList() ocp = OptimalControlProgram(nlp, use_sx=(cx == SX)) nlp.model = BiorbdModel( From 2547d7cf46df20e3ae1480fd58ed72dc245e5e7a Mon Sep 17 00:00:00 2001 From: Charbie Date: Mon, 7 Apr 2025 15:38:44 +0400 Subject: [PATCH 14/20] blacked --- bioptim/dynamics/ode_solvers.py | 1 - tests/shard1/test_dynamics.py | 116 ++++++++++++++-------------- tests/shard3/test_ligaments.py | 6 +- tests/shard3/test_passive_torque.py | 34 ++++---- 4 files changed, 80 insertions(+), 77 deletions(-) diff --git a/bioptim/dynamics/ode_solvers.py b/bioptim/dynamics/ode_solvers.py index 0ad0934f5..d4c78d6f3 100644 --- a/bioptim/dynamics/ode_solvers.py +++ b/bioptim/dynamics/ode_solvers.py @@ -58,7 +58,6 @@ class VARIATIONAL(RK): def integrator(self): return integrator.VARIATIONAL - class TRAPEZOIDAL(OdeSolverBase): """ A trapezoidal ode solver diff --git a/tests/shard1/test_dynamics.py b/tests/shard1/test_dynamics.py index 17122ee96..1cd3503cb 100644 --- a/tests/shard1/test_dynamics.py +++ b/tests/shard1/test_dynamics.py @@ -122,12 +122,12 @@ def test_torque_driven(with_contact, with_external_force, cx, phase_dynamics): external_force_set=external_forces, ) nlp.dynamics_type = Dynamics( - DynamicsFcn.TORQUE_DRIVEN, - with_contact=with_contact, - expand_dynamics=True, - phase_dynamics=phase_dynamics, - numerical_data_timeseries=numerical_time_series, - ) + DynamicsFcn.TORQUE_DRIVEN, + with_contact=with_contact, + expand_dynamics=True, + phase_dynamics=phase_dynamics, + numerical_data_timeseries=numerical_time_series, + ) nlp.cx = cx nlp.time_cx = cx.sym("time", 1, 1) @@ -222,12 +222,12 @@ def test_torque_driven_soft_contacts_dynamics(with_contact, cx, implicit_contact TestUtils.bioptim_folder() + "/examples/getting_started/models/2segments_4dof_2contacts.bioMod" ) nlp.dynamics_type = Dynamics( - DynamicsFcn.TORQUE_DRIVEN, - with_contact=with_contact, - soft_contacts_dynamics=implicit_contact, - expand_dynamics=True, - phase_dynamics=phase_dynamics, - ) + DynamicsFcn.TORQUE_DRIVEN, + with_contact=with_contact, + soft_contacts_dynamics=implicit_contact, + expand_dynamics=True, + phase_dynamics=phase_dynamics, + ) nlp.ns = N_SHOOTING nlp.cx = cx @@ -307,12 +307,12 @@ def test_torque_derivative_driven(with_contact, with_external_force, cx, phase_d external_force_set=external_forces, ) nlp.dynamics_type = Dynamics( - DynamicsFcn.TORQUE_DERIVATIVE_DRIVEN, - with_contact=with_contact, - expand_dynamics=True, - phase_dynamics=phase_dynamics, - numerical_data_timeseries=numerical_timeseries, - ) + DynamicsFcn.TORQUE_DERIVATIVE_DRIVEN, + with_contact=with_contact, + expand_dynamics=True, + phase_dynamics=phase_dynamics, + numerical_data_timeseries=numerical_timeseries, + ) nlp.cx = cx nlp.time_cx = cx.sym("time", 1, 1) @@ -451,12 +451,12 @@ def test_torque_derivative_driven_soft_contacts_dynamics(with_contact, cx, impli TestUtils.bioptim_folder() + "/examples/getting_started/models/2segments_4dof_2contacts.bioMod" ) nlp.dynamics_type = Dynamics( - DynamicsFcn.TORQUE_DERIVATIVE_DRIVEN, - with_contact=with_contact, - soft_contacts_dynamics=implicit_contact, - expand_dynamics=True, - phase_dynamics=phase_dynamics, - ) + DynamicsFcn.TORQUE_DERIVATIVE_DRIVEN, + with_contact=with_contact, + soft_contacts_dynamics=implicit_contact, + expand_dynamics=True, + phase_dynamics=phase_dynamics, + ) nlp.ns = N_SHOOTING nlp.cx = cx @@ -551,7 +551,9 @@ 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.dynamics_type = Dynamics( + dynamics, soft_contacts_dynamics=True, expand_dynamics=True, phase_dynamics=phase_dynamics + ) nlp.ns = N_SHOOTING nlp.cx = MX @@ -599,12 +601,12 @@ def test_torque_activation_driven(with_contact, with_external_force, cx, phase_d external_force_set=external_forces, ) nlp.dynamics_type = Dynamics( - DynamicsFcn.TORQUE_ACTIVATIONS_DRIVEN, - with_contact=with_contact, - expand_dynamics=True, - phase_dynamics=phase_dynamics, - numerical_data_timeseries=numerical_timeseries, - ) + DynamicsFcn.TORQUE_ACTIVATIONS_DRIVEN, + with_contact=with_contact, + expand_dynamics=True, + phase_dynamics=phase_dynamics, + numerical_data_timeseries=numerical_timeseries, + ) nlp.cx = cx nlp.time_cx = cx.sym("time", 1, 1) @@ -725,12 +727,12 @@ def test_torque_activation_driven_with_residual_torque( 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, - ) + 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) @@ -839,7 +841,9 @@ 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.dynamics_type = Dynamics( + DynamicsFcn.TORQUE_DRIVEN_FREE_FLOATING_BASE, expand_dynamics=True, phase_dynamics=phase_dynamics + ) nlp.ns = N_SHOOTING nlp.cx = cx @@ -913,14 +917,14 @@ def test_muscle_driven(with_excitations, with_contact, with_residual_torque, wit external_force_set=external_forces, ) nlp.dynamics_type = Dynamics( - DynamicsFcn.MUSCLE_DRIVEN, - with_residual_torque=with_residual_torque, - with_excitations=with_excitations, - with_contact=with_contact, - expand_dynamics=True, - phase_dynamics=phase_dynamics, - numerical_data_timeseries=numerical_timeseries, - ) + DynamicsFcn.MUSCLE_DRIVEN, + with_residual_torque=with_residual_torque, + with_excitations=with_excitations, + with_contact=with_contact, + expand_dynamics=True, + phase_dynamics=phase_dynamics, + numerical_data_timeseries=numerical_timeseries, + ) nlp.cx = cx nlp.time_cx = cx.sym("time", 1, 1) @@ -1088,10 +1092,10 @@ def test_joints_acceleration_driven(cx, phase_dynamics): 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, - ) + DynamicsFcn.JOINTS_ACCELERATION_DRIVEN, + expand_dynamics=True, + phase_dynamics=phase_dynamics, + ) nlp.ns = N_SHOOTING nlp.cx = cx @@ -1166,12 +1170,12 @@ def configure(ocp, nlp, with_contact=None, numerical_data_timeseries=None): TestUtils.bioptim_folder() + "/examples/getting_started/models/2segments_4dof_2contacts.bioMod" ) nlp.dynamics_type = Dynamics( - configure, - dynamic_function=custom_dynamic, - with_contact=with_contact, - expand_dynamics=True, - phase_dynamics=phase_dynamics, - ) + configure, + dynamic_function=custom_dynamic, + with_contact=with_contact, + expand_dynamics=True, + phase_dynamics=phase_dynamics, + ) nlp.ns = N_SHOOTING nlp.cx = MX diff --git a/tests/shard3/test_ligaments.py b/tests/shard3/test_ligaments.py index 9f32bcf33..188d5c90f 100644 --- a/tests/shard3/test_ligaments.py +++ b/tests/shard3/test_ligaments.py @@ -228,9 +228,9 @@ def test_muscle_driven_with_ligament(with_ligament, cx, phase_dynamics): TestUtils.bioptim_folder() + "/examples/muscle_driven_ocp/models/arm26_with_ligament.bioMod" ) nlp.dynamics_type = Dynamics( - DynamicsFcn.MUSCLE_DRIVEN, - with_ligament=with_ligament, - ) + DynamicsFcn.MUSCLE_DRIVEN, + with_ligament=with_ligament, + ) nlp.ns = 5 nlp.cx = cx diff --git a/tests/shard3/test_passive_torque.py b/tests/shard3/test_passive_torque.py index 083cc6bcb..b3d2626b4 100644 --- a/tests/shard3/test_passive_torque.py +++ b/tests/shard3/test_passive_torque.py @@ -43,10 +43,10 @@ def test_torque_driven_with_passive_torque(with_passive_torque, cx, phase_dynami 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, - ) + DynamicsFcn.TORQUE_DRIVEN, + with_passive_torque=with_passive_torque, + phase_dynamics=phase_dynamics, + ) nlp.ns = 5 nlp.cx = cx @@ -107,10 +107,10 @@ def test_torque_derivative_driven_with_passive_torque(with_passive_torque, cx, p 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, - ) + DynamicsFcn.TORQUE_DERIVATIVE_DRIVEN, + with_passive_torque=with_passive_torque, + phase_dynamics=phase_dynamics, + ) nlp.ns = 5 nlp.cx = cx @@ -200,11 +200,11 @@ def test_torque_activation_driven_with_passive_torque(with_passive_torque, with_ 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, - ) + 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 @@ -315,10 +315,10 @@ def test_muscle_driven_with_passive_torque(with_passive_torque, cx, phase_dynami 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, - ) + DynamicsFcn.MUSCLE_DRIVEN, + with_passive_torque=with_passive_torque, + phase_dynamics=phase_dynamics, + ) nlp.ns = 5 nlp.cx = cx From 163046332d3316ce64c38cf3e47d0dad4527c207 Mon Sep 17 00:00:00 2001 From: Charbie Date: Tue, 8 Apr 2025 11:47:54 +0400 Subject: [PATCH 15/20] bad merge --- bioptim/dynamics/configure_new_variable.py | 27 ++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/bioptim/dynamics/configure_new_variable.py b/bioptim/dynamics/configure_new_variable.py index f34b909d1..431b7db67 100644 --- a/bioptim/dynamics/configure_new_variable.py +++ b/bioptim/dynamics/configure_new_variable.py @@ -290,6 +290,33 @@ def _declare_cx_and_plot(self): ), ) + if self.as_algebraic_states: + 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.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( + self.name, + cx, + cx_scaled, + self.nlp.variable_mappings[self.name], + node_index, + ) + if not self.skip_plot: + self.nlp.plot[f"{self.name}_algebraic"] = CustomPlot( + lambda t0, phases_dt, node_idx, x, u, p, a, d: ( + a[self.nlp.algebraic_states.key_index(self.name), :] + if a.any() + else np.ndarray((cx[0].shape[0], 1)) * np.nan + ), + plot_type=PlotType.INTEGRATED, + axes_idx=self.axes_idx, + legend=self.legend, + combine_to=self.combine_name, + ) + def _manage_fatigue_to_new_variable( name: str, From 655105bb6f2cde095baf1c8ca06d2f7b2d0e4c4b Mon Sep 17 00:00:00 2001 From: Charbie Date: Wed, 9 Apr 2025 08:54:56 +0400 Subject: [PATCH 16/20] blacked --- .../static_arm_with_contact.py | 7 +++- .../muscle_activations_contacts_tracker.py | 7 +++- tests/shard1/test_dynamics.py | 42 ++++++++++--------- 3 files changed, 35 insertions(+), 21 deletions(-) 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 6ffc3be75..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], ode_solver=ode_solver) + 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 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 3bb22e7ed..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], ode_solver=ode_solver) + 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] diff --git a/tests/shard1/test_dynamics.py b/tests/shard1/test_dynamics.py index 81124dea9..246934a15 100644 --- a/tests/shard1/test_dynamics.py +++ b/tests/shard1/test_dynamics.py @@ -145,7 +145,7 @@ 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", @@ -246,14 +246,14 @@ 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_type = Dynamics( - DynamicsFcn.TORQUE_DRIVEN, - contact_type=[ContactType.RIGID_EXPLICIT] if with_contact else (), - soft_contacts_dynamics=implicit_contact, - expand_dynamics=True, - phase_dynamics=phase_dynamics, - ) + DynamicsFcn.TORQUE_DRIVEN, + contact_type=[ContactType.RIGID_EXPLICIT] if with_contact else (), + soft_contacts_dynamics=implicit_contact, + expand_dynamics=True, + phase_dynamics=phase_dynamics, + ) NonLinearProgram.add( ocp, "dynamics_type", @@ -344,11 +344,11 @@ 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 - + NonLinearProgram.add( ocp, "dynamics_type", - nlp.dynamics_type + nlp.dynamics_type, False, ) @@ -489,7 +489,7 @@ def test_torque_derivative_driven_soft_contacts_dynamics(contact_type, cx, phase ocp = OptimalControlProgram(nlp, use_sx=(cx == SX)) nlp.control_type = ControlType.CONSTANT - + NonLinearProgram.add( ocp, "dynamics_type", @@ -617,14 +617,14 @@ 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_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, - ) + 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, + ) NonLinearProgram.add( ocp, "dynamics_type", @@ -1260,7 +1260,11 @@ def test_with_contact_error(dynamics_fcn, phase_dynamics): # Dynamics dynamics = Dynamics( - dynamics_fcn, contact_type=[ContactType.RIGID_EXPLICIT], ode_solver=OdeSolver.RK4(), 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 From 4a5378fe19b7f20263e2ef65212f591d02badf10 Mon Sep 17 00:00:00 2001 From: Charbie Date: Wed, 9 Apr 2025 09:09:10 +0400 Subject: [PATCH 17/20] doc --- bioptim/dynamics/integrator.py | 1 + bioptim/dynamics/ode_solvers.py | 1 + 2 files changed, 2 insertions(+) diff --git a/bioptim/dynamics/integrator.py b/bioptim/dynamics/integrator.py index 3dd4d5b2d..33cdd4786 100644 --- a/bioptim/dynamics/integrator.py +++ b/bioptim/dynamics/integrator.py @@ -791,4 +791,5 @@ 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_solvers.py b/bioptim/dynamics/ode_solvers.py index d4c78d6f3..13000dc01 100644 --- a/bioptim/dynamics/ode_solvers.py +++ b/bioptim/dynamics/ode_solvers.py @@ -52,6 +52,7 @@ def integrator(self): 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 From 9ea234c3de138600e3aa6ac5dbd97815c6e3eab7 Mon Sep 17 00:00:00 2001 From: Charbie Date: Wed, 9 Apr 2025 09:13:32 +0400 Subject: [PATCH 18/20] fixed tests after manual merge --- tests/shard1/test_dynamics.py | 13 +++---------- 1 file changed, 3 insertions(+), 10 deletions(-) diff --git a/tests/shard1/test_dynamics.py b/tests/shard1/test_dynamics.py index 246934a15..86833088d 100644 --- a/tests/shard1/test_dynamics.py +++ b/tests/shard1/test_dynamics.py @@ -226,8 +226,8 @@ def test_torque_driven_soft_contacts_dynamics(contact_type, cx, phase_dynamics): ) nlp.dynamics_type = Dynamics( DynamicsFcn.TORQUE_DRIVEN, - with_contact=with_contact, - soft_contacts_dynamics=implicit_contact, + contact_type=contact_type, + soft_contacts_dynamics=True if ContactType.SOFT_IMPLICIT in contact_type else False, expand_dynamics=True, phase_dynamics=phase_dynamics, ) @@ -247,13 +247,6 @@ 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_type = Dynamics( - DynamicsFcn.TORQUE_DRIVEN, - contact_type=[ContactType.RIGID_EXPLICIT] if with_contact else (), - soft_contacts_dynamics=implicit_contact, - expand_dynamics=True, - phase_dynamics=phase_dynamics, - ) NonLinearProgram.add( ocp, "dynamics_type", @@ -1177,7 +1170,7 @@ def configure(ocp, nlp, with_contact=None, numerical_data_timeseries=None, conta nlp.dynamics_type = Dynamics( configure, dynamic_function=custom_dynamic, - contact_type=[ContactType.RIGID_EXPLICIT] if with_contact else (), + contact_type=contact_type, expand_dynamics=True, phase_dynamics=phase_dynamics, ) From 3f68250ccf7f1eff77b2ab000153c003a7937a17 Mon Sep 17 00:00:00 2001 From: Charbie Date: Wed, 9 Apr 2025 11:20:15 +0400 Subject: [PATCH 19/20] forgot one --- tests/shard1/test_dynamics.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/shard1/test_dynamics.py b/tests/shard1/test_dynamics.py index 86833088d..6165a2365 100644 --- a/tests/shard1/test_dynamics.py +++ b/tests/shard1/test_dynamics.py @@ -462,8 +462,8 @@ def test_torque_derivative_driven_soft_contacts_dynamics(contact_type, cx, phase ) nlp.dynamics_type = Dynamics( DynamicsFcn.TORQUE_DERIVATIVE_DRIVEN, - contact_type=[ContactType.RIGID_EXPLICIT] if with_contact else (), - soft_contacts_dynamics=implicit_contact, + contact_type=contact_type, + soft_contacts_dynamics=True if ContactType.SOFT_IMPLICIT in contact_type else False, expand_dynamics=True, phase_dynamics=phase_dynamics, ) From 3308d888591392ad9cbc47c01dd4964c1742a813 Mon Sep 17 00:00:00 2001 From: Charbie Date: Wed, 9 Apr 2025 18:33:19 +0400 Subject: [PATCH 20/20] Doc: readme update --- README.md | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index e99b62be5..59a3b5432 100644 --- a/README.md +++ b/README.md @@ -733,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*). @@ -948,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`. @@ -956,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