fixed the custom dynamics ?

bioptim/dynamics/configure_problem.py

@@ -134,7 +134,10 @@ def initialize(ocp, nlp):
         nlp.dynamics_type.type(ocp, nlp, **nlp.dynamics_type.extra_parameters)
 
     @staticmethod
-    def custom(ocp, nlp, **extra_params):
+    def custom(ocp,
+               nlp,
+               dynamics_constants_used_at_each_nodes: dict[list] = {},
+               **extra_params):
         """
         Call the user-defined dynamics configuration function
 
@@ -146,7 +149,7 @@ def custom(ocp, nlp, **extra_params):
             A reference to the phase
         """
 
-        nlp.dynamics_type.configure(ocp, nlp, **extra_params)
+        nlp.dynamics_type.configure(ocp, nlp, dynamics_constants_used_at_each_nodes, **extra_params)
 
     @staticmethod
     def torque_driven(
@@ -955,6 +958,7 @@ def configure_dynamics_function(ocp, nlp, dyn_func, **extra_params):
                 nlp.controls.scaled.mx_reduced,
                 nlp.parameters.scaled.mx_reduced,
                 nlp.algebraic_states.scaled.mx_reduced,
+                nlp.dynamics_constants.mx,
                 nlp,
                 **extra_params,
             )

bioptim/limits/penalty_helpers.py

@@ -135,22 +135,14 @@ def dynamics_constants(penalty, index, get_dynamics_constants: Callable, is_cons
             d = []
             phases, _, subnodes = _get_multinode_indices(penalty, is_constructing_penalty)
             for phase, sub in zip(phases, subnodes):
-                d += [_reshape_to_vector(get_dynamics_constants(phase, node, sub))]
-            return _vertcat(d)
+                d_tp = get_dynamics_constants(phase, node, sub)
+                if d_tp.shape != (0, 0):
+                    d += [_reshape_to_vector(get_dynamics_constants(phase, node, sub))]
+            return _vertcat(d) if len(d) > 0 else DM()
 
         d = get_dynamics_constants(penalty.phase, node, 0)  # cx_start
-        try:
-            if d.shape != (0, 0):
-                d = _reshape_to_vector(d)
-        except:
-            print("ici")
-
-        # if is_constructing_penalty:
-        #     d = _reshape_to_vector(get_dynamics_constants(penalty.phase, node, 0))  # cx_start
-        # else:
-        #     d0 = _reshape_to_vector(get_dynamics_constants(penalty.phase, node, 0))
-        #     d1 = _reshape_to_vector(get_dynamics_constants(penalty.phase, node, -1))
-        #     d = _vertcat([d0, d1])
+        if d.shape != (0, 0):
+            d = _reshape_to_vector(d)
 
         return d
 

bioptim/optimization/optimal_control_program.py

@@ -1620,7 +1620,7 @@ def _define_dynamics_constants(self, dynamics):
                     )
 
                     dynamics_constants[-1].append(
-                        name=key,
+                        name=f"{key}_{i_component}",
                         cx=[cx, cx, cx],
                         mx=mx,
                         bimapping=BiMapping(

bioptim/optimization/solution/solution.py

@@ -846,6 +846,7 @@ def _states_for_phase_integration(
         -------
         The states to integrate
         """
+        from ...interfaces.interface_utils import _get_dynamics_constants
 
         # In the case of multiple shootings, we don't need to do anything special
         if shooting_type == Shooting.MULTIPLE:
@@ -869,15 +870,21 @@ def _states_for_phase_integration(
             0,
             lambda p, n, sn: decision_controls[p][n][:, sn] if n < len(decision_controls[p]) else np.ndarray((0, 1)),
         )
-        s = PenaltyHelpers.states(
+        a = PenaltyHelpers.states(
             penalty,
             0,
             lambda p, n, sn: (
                 decision_algebraic_states[p][n][:, sn] if n < len(decision_algebraic_states[p]) else np.ndarray((0, 1))
             ),
         )
+        d_tp = PenaltyHelpers.dynamics_constants(
+            penalty,
+            0,
+            lambda p, n, sn: _get_dynamics_constants(self.ocp, p, n, sn)
+        )
+        d = np.array([]) if d_tp.shape == (0, 0) else np.array(d_tp)
 
-        dx = penalty.function[-1](t0, dt, x, u, params, s)
+        dx = penalty.function[-1](t0, dt, x, u, params, a, d)
         if dx.shape[0] != decision_states[phase_idx][0].shape[0]:
             raise RuntimeError(
                 f"Phase transition must have the same number of states ({dx.shape[0]}) "
@@ -897,6 +904,7 @@ def _integrate_stepwise(self) -> None:
         dict
             The integrated data structure similar in structure to the original _decision_states
         """
+        from ...interfaces.interface_utils import _get_dynamics_constants
 
         params = self._parameters.to_dict(to_merge=SolutionMerge.KEYS, scaled=True)[0][0]
         t_spans = self.t_span(time_alignment=TimeAlignment.CONTROLS)
@@ -908,6 +916,14 @@ def _integrate_stepwise(self) -> None:
 
         unscaled: list = [None] * len(self.ocp.nlp)
         for p, nlp in enumerate(self.ocp.nlp):
+            d = []
+            for n_idx in range(nlp.ns + 1):
+                d_tp = _get_dynamics_constants(self.ocp, p, n_idx, 0)
+                if d_tp.shape == (0, 0):
+                    d += [np.array([])]
+                else:
+                    d += [np.array(d_tp)]
+
             integrated_sol = solve_ivp_interface(
                 shooting_type=Shooting.MULTIPLE,
                 nlp=nlp,
@@ -916,6 +932,7 @@ def _integrate_stepwise(self) -> None:
                 u=u[p],
                 a=a[p],
                 p=params,
+                d=d,
                 method=SolutionIntegrator.OCP,
             )
 
@@ -1114,11 +1131,11 @@ def animate(
         self._check_models_comes_from_same_super_class()
 
         all_bioviz = []
-        for i, d in enumerate(data_to_animate):
+        for i, data in enumerate(data_to_animate):
             all_bioviz.append(
                 self.ocp.nlp[i].model.animate(
                     self.ocp,
-                    solution=d,
+                    solution=data,
                     show_now=show_now,
                     tracked_markers=tracked_markers,
                     **kwargs,