Skip to content

Commit

Permalink
Merge pull request #405 from Ipuch/SolutionIntegrator
Browse files Browse the repository at this point in the history
Add the choice for an integrator of the solution
pariterre authored Nov 12, 2021
2 parents ef7c099 + e16441c commit 46c3f9d
Showing 8 changed files with 140 additions and 80 deletions.
12 changes: 11 additions & 1 deletion README.md
Original file line number Diff line number Diff line change
@@ -88,6 +88,7 @@ As a tour guide that uses this binder, you can watch the `bioptim` workshop that
- [InterpolationType](#enum-interpolationtype)
- [Shooting](#enum-shooting)
- [CostType](#enum-costtype)
- [SolutionIntegrator](#enum-solutionintegrator)

[Examples](#examples)
- [Run examples](#run-examples)
@@ -1383,7 +1384,6 @@ The accepted values are:
- SPLINE: Requires five columns. It performs a cubic spline to interpolate between the nodes.
- CUSTOM: User defined interpolation function


### Enum: Shooting
The type of integration to perform
- MULTIPLE: resets the state at each node
@@ -1396,6 +1396,16 @@ The type of cost
- CONSTRAINTS: The constraints
- ALL: All the previously described cost type

### Enum: SolutionIntegrator
The type of integrator used to integrate the solution of the optimal control problem
- DEFAULT: The default integrator initially chosen with [OdeSolver](#class-odesolver)
- SCIPY_RK23: The scipy integrator RK23
- SCIPY_RK45: The scipy integrator RK45
- SCIPY_DOP853: The scipy integrator DOP853
- SCIPY_BDF: The scipy integrator BDF
- SCIPY_LSODA: The scipy integrator LSODA


# Examples
In this section, you will find the description of all the examples implemented with bioptim. They are ordered in
separate files. Each subsection corresponds to the different files, dealing with different examples and topics.
12 changes: 11 additions & 1 deletion bioptim/__init__.py
Original file line number Diff line number Diff line change
@@ -165,7 +165,17 @@
from .limits.path_conditions import BoundsList, Bounds, InitialGuessList, InitialGuess, QAndQDotBounds
from .limits.fatigue_path_conditions import FatigueBounds, FatigueInitialGuess
from .limits.penalty_node import PenaltyNode, PenaltyNodeList
from .misc.enums import Axis, Node, InterpolationType, PlotType, ControlType, CostType, Shooting, VariableType
from .misc.enums import (
Axis,
Node,
InterpolationType,
PlotType,
ControlType,
CostType,
Shooting,
VariableType,
SolutionIntegrator,
)
from .misc.mapping import BiMappingList, BiMapping, Mapping
from .optimization.non_linear_program import NonLinearProgram
from .optimization.optimal_control_program import OptimalControlProgram
32 changes: 22 additions & 10 deletions bioptim/gui/plot.py
Original file line number Diff line number Diff line change
@@ -10,7 +10,7 @@
from casadi import Callback, nlpsol_out, nlpsol_n_out, Sparsity, DM

from ..limits.path_conditions import Bounds
from ..misc.enums import PlotType, ControlType, InterpolationType, Shooting
from ..misc.enums import PlotType, ControlType, InterpolationType, Shooting, SolutionIntegrator
from ..misc.mapping import Mapping
from ..optimization.solution import Solution

@@ -197,7 +197,7 @@ def __init__(
automatically_organize: bool = True,
show_bounds: bool = False,
shooting_type: Shooting = Shooting.MULTIPLE,
use_scipy_integrator: bool = False,
integrator: SolutionIntegrator = SolutionIntegrator.DEFAULT,
):
"""
Prepares the figures during the simulation
@@ -212,8 +212,9 @@ def __init__(
If the axes should fit the bounds (True) or the data (False)
shooting_type: Shooting
The type of integration method
use_scipy_integrator: bool
Use the scipy solve_ivp integrator for RungeKutta 45 instead of currently defined integrator
integrator: SolutionIntegrator
Use the ode defined by OCP or use a separate integrator provided by scipy
"""
for i in range(1, ocp.n_phases):
if len(ocp.nlp[0].states["q"]) != len(ocp.nlp[i].states["q"]):
@@ -235,7 +236,8 @@ def __init__(

self.t = []
self.t_integrated = []
self.use_scipy_integrator = use_scipy_integrator
self.integrator = integrator

if isinstance(self.ocp.original_phase_time, (int, float)):
self.tf = [self.ocp.original_phase_time]
else:
@@ -294,12 +296,14 @@ def __update_time_vector(self):
self.t_integrated = []
last_t = 0
for phase_idx, nlp in enumerate(self.ocp.nlp):
n_int_steps = nlp.ode_solver.steps_scipy if self.use_scipy_integrator else nlp.ode_solver.steps
n_int_steps = (
nlp.ode_solver.steps_scipy if self.integrator != SolutionIntegrator.DEFAULT else nlp.ode_solver.steps
)
dt_ns = self.tf[phase_idx] / nlp.ns
time_phase_integrated = []
last_t_int = copy(last_t)
for _ in range(nlp.ns):
if nlp.ode_solver.is_direct_collocation and not self.use_scipy_integrator:
if nlp.ode_solver.is_direct_collocation and self.integrator == SolutionIntegrator.DEFAULT:
time_phase_integrated.append(np.array(nlp.dynamics[0].step_time) * dt_ns + last_t_int)
else:
time_phase_integrated.append(np.linspace(last_t_int, last_t_int + dt_ns, n_int_steps + 1))
@@ -426,7 +430,11 @@ def legend_without_duplicate_labels(ax):
)
elif plot_type == PlotType.INTEGRATED:
plots_integrated = []
n_int_steps = nlp.ode_solver.steps_scipy if self.use_scipy_integrator else nlp.ode_solver.steps
n_int_steps = (
nlp.ode_solver.steps_scipy
if self.integrator != SolutionIntegrator.DEFAULT
else nlp.ode_solver.steps
)
zero = np.zeros(n_int_steps + 1)
color = self.plot_func[variable][i].color if self.plot_func[variable][i].color else "tab:brown"
for cmp in range(nlp.ns):
@@ -585,7 +593,7 @@ def update_data(self, v: dict):
continuous=False,
shooting_type=self.shooting_type,
keep_intermediate_points=True,
use_scipy_integrator=self.use_scipy_integrator,
integrator=self.integrator,
).states
data_controls = sol.controls
data_params = sol.parameters
@@ -598,7 +606,11 @@ def update_data(self, v: dict):
self.__update_xdata()

for i, nlp in enumerate(self.ocp.nlp):
step_size = nlp.ode_solver.steps_scipy + 1 if self.use_scipy_integrator else nlp.ode_solver.steps + 1
step_size = (
nlp.ode_solver.steps_scipy + 1
if self.integrator != SolutionIntegrator.DEFAULT
else nlp.ode_solver.steps + 1
)
n_elements = nlp.ns * step_size + 1

state = np.ndarray((0, n_elements))
13 changes: 13 additions & 0 deletions bioptim/misc/enums.py
Original file line number Diff line number Diff line change
@@ -103,3 +103,16 @@ class VariableType(Enum):

STATES = "states"
CONTROLS = "controls"


class SolutionIntegrator(Enum):
"""
Selection of integrator to use integrate function
"""

DEFAULT = "DEFAULT"
SCIPY_RK23 = "RK23"
SCIPY_RK45 = "RK45"
SCIPY_DOP853 = "DOP853"
SCIPY_BDF = "BDF"
SCIPY_LSODA = "LSODA"
10 changes: 5 additions & 5 deletions bioptim/optimization/optimal_control_program.py
Original file line number Diff line number Diff line change
@@ -28,7 +28,7 @@
from ..limits.penalty import PenaltyOption
from ..limits.objective_functions import ObjectiveFunction
from ..misc.__version__ import __version__
from ..misc.enums import ControlType, SolverType, Shooting, PlotType, CostType
from ..misc.enums import ControlType, SolverType, Shooting, PlotType, CostType, SolutionIntegrator
from ..misc.mapping import BiMappingList, Mapping
from ..misc.utils import check_version
from ..optimization.parameters import ParameterList, Parameter
@@ -749,7 +749,7 @@ def prepare_plots(
automatically_organize: bool = True,
show_bounds: bool = False,
shooting_type: Shooting = Shooting.MULTIPLE,
use_scipy_integrator: bool = False,
integrator: SolutionIntegrator = SolutionIntegrator.DEFAULT,
) -> PlotOcp:
"""
Create all the plots associated with the OCP
@@ -762,8 +762,8 @@ def prepare_plots(
If the ylim should fit the bounds
shooting_type: Shooting
What type of integration
use_scipy_integrator: bool
Use the scipy solve_ivp integrator for RungeKutta 45 instead of currently defined integrator
integrator: SolutionIntegrator
Use the ode defined by OCP or use a separate integrator provided by scipy
Returns
-------
@@ -775,7 +775,7 @@ def prepare_plots(
automatically_organize=automatically_organize,
show_bounds=show_bounds,
shooting_type=shooting_type,
use_scipy_integrator=use_scipy_integrator,
integrator=integrator,
)

def solve(
43 changes: 26 additions & 17 deletions bioptim/optimization/solution.py
Original file line number Diff line number Diff line change
@@ -9,7 +9,7 @@
from matplotlib import pyplot as plt

from ..limits.path_conditions import InitialGuess, InitialGuessList
from ..misc.enums import ControlType, CostType, Shooting, InterpolationType, SolverType
from ..misc.enums import ControlType, CostType, Shooting, InterpolationType, SolverType, SolutionIntegrator
from ..misc.utils import check_version
from ..optimization.non_linear_program import NonLinearProgram
from ..optimization.optimization_variable import OptimizationVariableList, OptimizationVariable
@@ -492,7 +492,7 @@ def integrate(
keep_intermediate_points: bool = False,
merge_phases: bool = False,
continuous: bool = True,
use_scipy_integrator: bool = False,
integrator: SolutionIntegrator = SolutionIntegrator.DEFAULT,
) -> Any:
"""
Integrate the states
@@ -509,8 +509,8 @@ def integrate(
continuous: bool
If the arrival value of a node should be discarded [True] or kept [False]. The value of an integrated
arrival node and the beginning of the next one are expected to be almost equal when the problem converged
use_scipy_integrator: bool
Ignore the dynamics defined by OCP and use a separate integrator provided by scipy
integrator: SolutionIntegrator
Use the ode defined by OCP or use a separate integrator provided by scipy
Returns
-------
@@ -535,7 +535,7 @@ def integrate(
"Shooting.SINGLE_CONTINUOUS and continuous=False cannot be used simultaneously it is a contradiction"
)

out = self.__perform_integration(shooting_type, keep_intermediate_points, continuous, use_scipy_integrator)
out = self.__perform_integration(shooting_type, keep_intermediate_points, continuous, integrator)

if merge_phases:
if continuous:
@@ -548,24 +548,25 @@ def integrate(
return out

def __perform_integration(
self, shooting_type: Shooting, keep_intermediate_points: bool, continuous: bool, use_scipy_integrator: bool
self, shooting_type: Shooting, keep_intermediate_points: bool, continuous: bool, integrator: SolutionIntegrator
):
n_direct_collocation = sum([nlp.ode_solver.is_direct_collocation for nlp in self.ocp.nlp])
if n_direct_collocation > 0 and not use_scipy_integrator:

if n_direct_collocation > 0 and integrator == SolutionIntegrator.DEFAULT:
if continuous:
raise RuntimeError(
"Integration with direct collocation must be not continuous if not use_scipy_integrator"
"Integration with direct collocation must be not continuous if a scipy integrator is used"
)

if shooting_type != Shooting.MULTIPLE:
raise RuntimeError(
"Integration with direct collocation must using shooting_type=Shooting.MULTIPLE "
"if not use_scipy_integrator"
"if a scipy integrator is not used"
)

# Copy the data
out = self.copy(skip_data=True)
out.recomputed_time_steps = use_scipy_integrator
out.recomputed_time_steps = integrator != SolutionIntegrator.DEFAULT
out._states = []
for _ in range(len(self._states)):
out._states.append({})
@@ -575,7 +576,7 @@ def __perform_integration(
for p, nlp in enumerate(self.ocp.nlp):
param_scaling = nlp.parameters.scaling
n_states = self._states[p]["all"].shape[0]
n_steps = nlp.ode_solver.steps_scipy if use_scipy_integrator else nlp.ode_solver.steps
n_steps = nlp.ode_solver.steps_scipy if integrator != SolutionIntegrator.DEFAULT else nlp.ode_solver.steps
if not continuous:
n_steps += 1
if keep_intermediate_points:
@@ -596,7 +597,11 @@ def __perform_integration(
)
x0 += np.array(val)[:, 0]
else:
col = slice(0, n_steps) if nlp.ode_solver.is_direct_collocation and not use_scipy_integrator else 0
col = (
slice(0, n_steps)
if nlp.ode_solver.is_direct_collocation and integrator == SolutionIntegrator.DEFAULT
else 0
)
x0 = self._states[p]["all"][:, col]

for s in range(self.ns[p]):
@@ -607,13 +612,17 @@ def __perform_integration(
else:
raise NotImplementedError(f"ControlType {nlp.control_type} " f"not yet implemented in integrating")

if use_scipy_integrator:
if integrator != SolutionIntegrator.DEFAULT:
t_init = sum(out.phase_time[:p]) / nlp.ns
t_end = sum(out.phase_time[: (p + 2)]) / nlp.ns
n_points = n_steps + 1 if continuous else n_steps
t_eval = np.linspace(t_init, t_end, n_points) if keep_intermediate_points else [t_init, t_end]
integrated = solve_ivp(
lambda t, x: np.array(nlp.dynamics_func(x, u, params))[:, 0], [t_init, t_end], x0, t_eval=t_eval
lambda t, x: np.array(nlp.dynamics_func(x, u, params))[:, 0],
[t_init, t_end],
x0,
t_eval=t_eval,
method=integrator.value,
).y

next_state_col = (
@@ -853,7 +862,7 @@ def graphs(
show_bounds: bool = False,
show_now: bool = True,
shooting_type: Shooting = Shooting.MULTIPLE,
use_scipy_integrator: bool = False,
integrator: SolutionIntegrator = SolutionIntegrator.DEFAULT,
):
"""
Show the graphs of the simulation
@@ -868,14 +877,14 @@ def graphs(
If the show method should be called. This is blocking
shooting_type: Shooting
The type of interpolation
use_scipy_integrator: bool
integrator: SolutionIntegrator
Use the scipy solve_ivp integrator for RungeKutta 45 instead of currently defined integrator
"""

if self.is_merged or self.is_interpolated or self.is_integrated:
raise NotImplementedError("It is not possible to graph a modified Solution yet")

plot_ocp = self.ocp.prepare_plots(automatically_organize, show_bounds, shooting_type, use_scipy_integrator)
plot_ocp = self.ocp.prepare_plots(automatically_organize, show_bounds, shooting_type, integrator)
plot_ocp.update_data(self.vector)
if show_now:
plt.show()
95 changes: 50 additions & 45 deletions tests/test_simulate.py

Large diffs are not rendered by default.

3 changes: 2 additions & 1 deletion tests/utils.py
Original file line number Diff line number Diff line change
@@ -18,6 +18,7 @@
InitialGuess,
Shooting,
Solver,
SolutionIntegrator,
)


@@ -117,7 +118,7 @@ def simulate(sol, decimal_value=7):
merge_phases=True,
shooting_type=Shooting.SINGLE_CONTINUOUS,
keep_intermediate_points=True,
use_scipy_integrator=False,
integrator=SolutionIntegrator.DEFAULT,
)

# Evaluate the final error of the single shooting integration versus the finale node

0 comments on commit 46c3f9d

Please sign in to comment.