system¶

The System class manages the simulation (integration) of a system whose equations are given by KanesMethod.

Many of the attributes are also properties, and can be directly modified.

Here is the procedure for using this class.

specify your options either via the constructor or via the attributes.

optionally, call generate_ode_function() if you want to customize how the ODE function is generated.

call integrate() to simulate your system.

The simplest usage of this class is as follows. First, we need a KanesMethod object on which we have already invoked kanes_equations():

>>> from sympy.physics.mechanics.models import n_link_pendulum_on_cart
>>> from pydy.system import System
>>> import numpy as np
>>> kane = n_link_pendulum_on_cart()
>>> times = np.linspace(0.0, 5.0, num=3)
>>> sys = System(kane, times=times)
>>> sys.integrate()
array([[0., 0., 0., 0.],
       [0., 0., 0., 0.],
       [0., 0., 0., 0.]])

In this case, we use defaults for the numerical values of the constants, specified quantities, initial conditions, etc. You probably won’t like these defaults. You can also specify such values via constructor keyword arguments or via the attributes:

>>> import sympy as sm
>>> sys = System(kane,
...              initial_conditions={kane.q[1]: 0.5},
...              times=times)
...
>>> g, l0, m0, m1 = list(sm.ordered(sys.constants_symbols))
>>> sys.constants = {m1: 5.0}
>>> sys.integrate()
array([[ 0.        ,  0.5       ,  0.        ,  0.        ],
       [-1.12276473,  4.19253522, -0.77003647,  1.86016638],
       [-1.00443253,  5.47085374, -0.4536987 , -0.7915558 ]])

To double-check the constants, specifieds, states and times in your problem, look at these properties:

>>> sys.coordinates
[q0(t), q1(t)]
>>> sys.speeds
[u0(t), u1(t)]
>>> sys.states
[q0(t), q1(t), u0(t), u1(t)]
>>> sys.constants_symbols
{g, l0, m0, m1}
>>> sys.specifieds_symbols
{F(t)}
>>> sys.times
array([0. , 2.5, 5. ])

You can also add additional equations to evaluate alongside the differential equations:

>>> k0 = sm.Symbol('k0')
>>> sys = System(kane,
...              initial_conditions={kane.q[1]: 0.5},
...              times=times,
...              outputs={k0: m0*sys.speeds[0]**2/2})
>>> x = sys.integrate()
>>> sys.outputs_symbols
[k0]
>>> sys.evaluate_outputs(x=x)
array([[0.00000000e+00],
       [8.92619991e-01],
       [7.01894807e-04]])

In the prior examples, the System generates the numerical ode function for you behind the scenes. If you want to customize how this function is generated, you must call generate_ode_function() on your own:

>>> rhs = sys.generate_ode_function(generator='cython')
>>> sys.evaluate_ode_function == rhs
True
>>> help(sys.evaluate_ode_function)
Help on function rhs in module pydy.codegen.ode_function_generators:

rhs(*args)
    Returns the derivatives of the states, i.e. numerically evaluates the right
    hand side of the first order differential equation.

    x' = f(x, t, r, p)

    Parameters
    ==========
    x : ndarray, shape(4,)
        The state vector is ordered as such:
            - q0(t)
            - q1(t)
            - u0(t)
            - u1(t)
    t : float
        The current time.
    r : dictionary; ndarray, shape(1,); function

        There are three options for this argument. (1) is more flexible but
        (2) and (3) are much more efficient.

        (1) A dictionary that maps the specified functions of time to floats,
        ndarrays, or functions that produce ndarrays. The keys can be a single
        specified symbolic function of time or a tuple of symbols. The total
        number of symbols must be equal to 1. If the value is a
        function it must be of the form g(x, t), where x is the current state
        vector ndarray and t is the current time float and it must return an
        ndarray of the correct shape. For example::

          r = {a: 1.0,
               (d, b) : np.array([1.0, 2.0]),
               (e, f) : lambda x, t: np.array(x[0], x[1]),
               c: lambda x, t: np.array(x[2])}

        (2) A ndarray with the specified values in the correct order and of the
        correct shape.

        (3) A function that must be of the form g(x, t), where x is the current
        state vector and t is the current time and it must return an ndarray of
        the correct shape.

        The specified inputs are, in order:
            - F(t)
    p : dictionary len(4) or ndarray shape(4,)
        Either a dictionary that maps the constants symbols to their numerical
        values or an array with the constants in the following order:
            - g
            - l0
            - m0
            - m1

    Returns
    =======
    dx : ndarray, shape(4,)
        The derivative of the state vector.
    y : ndarray, shape(1,)
        Values of the provided outputs.
            - y0(t)

>>> sys.integrate()
array([[ 0.        ,  0.5       ,  0.        ,  0.        ],
       [-0.31425675,  3.29123866, -1.33612873,  2.70246056],
       [-0.48148282,  5.77849021, -0.03746718, -0.08560791]])

class pydy.system.System(eom_method, constants=None, specifieds=None, ode_solver=None, initial_conditions=None, times=None, outputs=None, noncontributing_forces=None, constants_symbols=None, specifieds_symbols=None)[source]¶

Multibody dynamics system for simulation and numerical evaluation.

See the class’s attributes for a description of the arguments to this constructor.

The parameters to this constructor are all attributes of the System. With the exception of eom_method, these attributes can be modified directly at any future point.

Parameters:

eom_methodsympy.physics.mechanics.kane.KanesMethod: You must have called kanes_equations() before constructing this system.
constantsdict, optional (default: all 1.0): This dictionary maps SymPy Symbol objects to floats.
specifiedsdict, optional (default: all 0.0): This dictionary maps SymPy Functions of time objects, or tuples of them, to floats, NumPy arrays, or functions of the state and time.
ode_solverfunction, optional: This function computes the derivatives of the states. The default is scipy.integrate.odeint().
initial_conditionsdict, optional (default: all zero): This dictionary maps SymPy Functions of time objects to floats.
timesarray_like, shape(n,), optional: An array_like object, which contains time values over which equations are integrated. It has to be supplied before integrate() can be called.
outputsdictionary, optional: Maps functions of time or tuples of functions of time to expressions or iterables of expressions, respectively. In general, the expressions should be a function of the state, constants, and specifieds. Expressions that are linear in the functions of time and/or the time derivatives of the speeds are also supported, but not yet nonlinear functions of these variables.
noncontributing_forcesiterable of Functions of time, optional: If the eom_method includes noncontributig forces (Kane’s method), provide a list of variable names for these forces and they will be computed when evaluating the differential equations.
constants_symbolsiterable of Symbol, optional: If provided, the system’s equations will not be searched for the minimal set of constants. It is best to provide these for large system equations, as the search can be prohibitively long in duration.
specifieds_symbolsiterable of Functions of time, optional: If provided, the system’s equations will not be searched for the minimal set of specifieds. It is best to provide these for large system equations, as the search can be prohibitively long in duration.

property auxiliaries¶: Returns a list of the symbols representing the system’s auxiliary states which are the time integrals of any outputs that are linear functions of the time derivatices of the generalized speeds.

property constants¶: A dict that provides the numerical values for the constants in the problem (all non-dynamics symbols). Keys are the symbols for the constants, and values are floats. Constants that are not specified in this dict are given a default value of 1.0.

property constants_symbols¶: A set of the symbolic constants (not functions of time) in the system.

property constraints¶: A column matrix of configuration and nonholonomic constraints expressions, ordered as stored in KanesMethod.

property coordinates¶: Returns a list of the symbolic functions of time representing the system’s generalized coordinates.

property eom_method¶: This is a KanesMethod. The method used to generate the equations of motion. Read-only.

evaluate_constraints(x=None, t=None)[source]¶

Returns the values of the configuration and motion constraints at the initial condition or, alternatively, for the provided state vector.

Parameters:

xarray_like, shape(n,) or shape(m, n), optional: State vector of n states or a series of m state vectors.
tfloat or array_like, shape(m,), optional: Time or m time values.

Returns:

ndarray, shape(o,) or shape(m, o): Constraint vector of o constraints or a series of m constraint vectors.

Notes

To see the order of the state values use:

system = System(...)
system.states

or:

rhs = system.generate_ode_function()
help(rhs)

evaluate_holonomic_constraints(x=None, t=None)[source]¶

Returns the values of the configuration at the initial condition or, alternatively, for the provided state vector.

Parameters:

xarray_like, shape(n,) or shape(m, n), optional: State vector of n states or a series of m state vectors.
tfloat or array_like, shape(m,), optional: Time or m time values.

Returns:

ndarray, shape(o,) or shape(m, o): Constraint vector of o constraints or a series of m constraint vectors.

Notes

To see the order of the state values use:

system = System(...)
system.states

or:

rhs = system.generate_ode_function()
help(rhs)

evaluate_ode(x=None, t=None)[source]¶

Returns the right hand side of the differential equations. The default is to evaluate at the set initial_conditions at the first time value or with t=0 if times is not set. Pass in optional arguments to override using the initial state and time.

Parameters:

xarray_like, shape(n,) or shape(m, n), optional: State values at time t.
tfloat or array_like, shape(m,), optional: Time or m time values.

Returns:

xdndarray, shape(n,) or shape(m, n): Time derivative of the states at time t.

Notes

This method is present for convenience, it is not designed to be used where performance matters, use evaluate_ode_function directly when performance is needed.

To see the order of the state values use:

system = System(...)
system.states

or:

rhs = system.generate_ode_function()
help(rhs)

property evaluate_ode_function¶

A function generated by generate_ode_function() that computes the state derivatives:

xd = evaluate_ode_function(x, t, *args)

This function is used by the ode_solver.

To see the autogenerated docstring and expected arguments call help():

help(system.evaluate_ode_function)

evaluate_outputs(x=None, t=None)[source]¶

Returns an array of the evaluated outputs. The default is to evaluate at the initial conditions at the first time value. Pass in optional arguments to override the state or time.

Parameters:

xarray_like, shape(n,) or shape(m, n), optional: State values at time t.
tfloat or array_like, shape(m,), optional: Time or m time values.

Returns:

yndarray, shape(o,) or shape(m, o): o output values at time t.

Notes

This method is present for convenience, it is not designed to be used where performance matters, use evaluate_ode_function directly when performance is needed.

To see the order of the state values use:

system = System(...)
system.states

or:

rhs = system.generate_ode_function()
help(rhs)

evaluate_velocity_constraints(x=None, t=None)[source]¶

Returns the values of the velocity constraints at the initial condition or, alternatively, for the provided state vector.

Parameters:

xarray_like, shape(n,) or shape(m, n), optional: State vector of n states or a series of m state vectors.
tfloat or array_like, shape(m,), optional: Time or m time values.

Returns:

ndarray, shape(o,) or shape(m, o): Constraint vector of o constraints or a series of m constraint vectors.

Notes

To see the order of the state values use:

system = System(...)
system.states

or:

rhs = system.generate_ode_function()
help(rhs)

generate_ode_function(**kwargs)[source]¶

Returns a function generated from generate_ode_function() with the appropriate arguments and also sets the evaluate_ode_function attribute to the resulting function.

Parameters:

kwargs: All other kwargs are passed onto pydy.codegen.ode_function_generators.generate_ode_function(). Don’t specify the specifieds keyword argument though; the System class takes care of those.

Returns:

evaluate_ode_functionfunction: A function which evaluates the derivaties of the states.

Notes

If the Cython generator is selected and you have a custom ode_solver set, keyword argument force_c_contiguous will be automatically set to True. You can disable this by setting it to False but you must ensure ensure that ode solver only passes C contiguous arrays to the generated ode function. Forcing C contiguous arrays introduces a small performance penalty due to the necessity of copying arrays.

property initial_conditions¶: Initial conditions for all states (coordinates and speeds). Keys are the symbols for the coordinates and speeds, and values are floats. Coordinates or speeds that are not specified in this dict are given a default value of 0.0.

integrate(**solver_kwargs)[source]¶

Integrates the equations evaluate_ode_function using ode_solver.

It is necessary to have first generated an ode function. If you have not done so, we do so automatically by invoking generate_ode_function(). However, if you want to customize how this function is generated (e.g., change the generator to cython), you can call generate_ode_function() on your own (before calling integrate()).

Parameters:

**solver_kwargs: Optional arguments that are passed on to the ode_solver.

Returns:

x_historyndarray, shape(num_integrator_time_steps, num_states): The trajectory of states (coordinates and speeds) through the requested time interval. num_integrator_time_steps is either len(times) if len(times) > 2, or is determined by the ode_solver.

property noncontributing_forces¶: List of symbolic functions of time representing the noncontributing forces (force & torque measure numbers) associated with auxiliary speeds.

property num_auxiliaries¶: Returns the number of auxiliaries.

property num_constants¶: Returns the number of constants.

property num_constraints¶: Total number of configuration and nonholonomic constaints.

property num_coordinates¶: Returns the number of coordinates.

property num_holonomic_constraints¶: Number of configuration constraints.

property num_nonholonomic_constraints¶: Number of nonholonomic constraints.

property num_outputs¶: Returns the number of outputs.

property num_specifieds¶: Returns the number of specifieds.

property num_speeds¶: Returns the number of speeds.

property num_states¶: Returns the number of states.

property num_velocity_constraints¶: Number of motion constraints.

property ode_solver¶

A function that performs forward integration. It must have the same signature as scipy.integrate.odeint(), which is:

x_history = ode_solver(f, x0, t, args=f_args)

where f is a function f(x, t, *f_args), x0 are the initial conditions, x_history is the state time history, x is the state, t is the time, and args is a keyword argument takes arguments that are then passed to f. The default solver is scipy.integrate.odeint().

Examples

SciPy introduced a unified scipy.integrate.solve_ivp() API which can be used with PyDy. solve_ivp requires a function that has swapped first arguments and it returns a solution object where the trajectory is the transpose of what odeint outputs. You can make a custom ODE solver function to use solve_ivp like so:

>>> from pydy.models import multi_mass_spring_damper
>>> sys = multi_mass_spring_damper()
>>> sys.initial_conditions[sys.coordinates[0]] = 1.0
>>> sys.times = [1.0, 2.0, 3.0]
>>> from scipy.integrate import solve_ivp
>>> def custom_ode_solver(f, x0, ts, args=(), **kwargs):
...     return solve_ivp(lambda t, x: f(x, t, *args), ts[[0, -1]], x0,
...                      t_eval=ts, **kwargs).y.T
>>> sys.ode_solver = custom_ode_solver

This then allows one to easiliy change methods and settings following SciPy’s API:

>>> sys.integrate(method='LSODA', rtol=1e-10)
array([[ 1.00000000e+00, -5.67952532e-17],
       [ 6.59700039e-01, -5.33506568e-01],
       [ 1.50574778e-01, -4.19279930e-01]])
>>> sys.integrate(method='RK23', rtol=1e-12)
array([[ 1.        ,  0.        ],
       [ 0.65970115, -0.53350624],
       [ 0.15057689, -0.41928088]])

property outputs¶

Dictionary of functions of time or utple of functions of time mapped to SymPy expressions or iterables of expressions that represent extra functions of the state that should be evaluted alongside the ordinary differential equations. Acceptable key pairs for this dictionary take the following three forms:

A single function of time mapped to a function of the state:

outputs[p(t)] = k*x(t)

A tuple of functions of time mapped to functions of the state:

outputs[(f1(t), f2(t))] = (k*x(t), c*v(t))

A tuple of functions of time mapped to a system of linear equations in the functions and the time derivatives of the states:

outputs[(m1(t), m2(t))] = (m1(t) - 4*m2(t) + k*v(t).diff(t) + 2,
                           m1(t) + 3*m2(t) - omega(t).diff(t))

If equations of the last form are provided, this linear system will be numerically solved alongside the ordinary differential equations.

Notes

If your system has configuration or motion constraints, these will automatically be added to the outputs dictionary. If your system has noncontributing forces exposed and you provide names for those forces, these will automatically be added to the outputs dictionary.

set_dependent_initial_conditions(dep_vars=None, use_jac=False, **root_kwargs)[source]¶

Sets the initial conditions of the dependent coordinates and dependent speeds using the holonomic and nonholonomic constraints, respectively.

Parameters:

dep_varsiterable of Function()(t), optional: Dependent coordinates and speeds to solve for. The number of coordinates should be equal to the number of holonomic constraints. The number of speeds should be equal to the number of nonholonic constraints. If None, the dependent coordinates and speeds are those used in KanesMethod instantiation.
use_jacboolean, optional: If true the Jacobian of the constraint equations will be used to solve the constraint equations for the dependent states.
root_kwargs: Extra keyword arguments that are passed to scipy.optimize.root().

property specifieds¶

A dict that provides numerical values for the specified quantities in the problem (all dynamicsymbols that are not defined by the equations of motion). There are two possible formats. (1) is more flexible, but (2) is more efficient (by a factor of 3).

(1) Keys are the symbols for the specified quantities, or a tuple of symbols, and values are the floats, arrays of floats, or functions that generate the values. If a dictionary value is a function, it must have the same signature as f(x, t), the ode right-hand-side function (see the documentation for the ode_solver attribute). You needn’t provide values for all specified symbols. Those for which you do not give a value will default to 0.0.

(2) There are two keys: ‘symbols’ and ‘values’. The value for ‘symbols’ is an iterable of all the specified quantities in the order that you have provided them in ‘values’. Values is an ndarray, whose length is num_specifieds, or a function of x and t that returns an ndarray (also of length num_specifieds). NOTE: You must provide values for all specified symbols. In this case, we do not provide default values.

NOTE: If you switch formats with the same instance of System, you must call generate_ode_function() before calling integrate() again.

Examples

Here are examples for (1). Keys can be individual symbols, or a tuple of symbols. Length of a value must match the length of the corresponding key. Values can be functions that return iterables:

sys = System(km)
sys.specifieds = {(a, b, c): np.ones(3), d: lambda x, t: -3 * x[0]}
sys.specifieds = {(a, b, c): lambda x, t: np.ones(3)}

Here are examples for (2):

sys.specifieds = {'symbols': (a, b, c, d),
                  'values': np.ones(4)}
sys.specifieds = {'symbols': (a, b, c, d),
                  'values': lambda x, t: np.ones(4)}

property specifieds_symbols¶: A set of the dynamicsymbols you must specify.

property speeds¶: Returns a list of the symbolic functions of time representing the system’s generalized speeds.

property states¶: Returns a list of the symbolic functions of time representing the system’s states, i.e. generalized coordinates plus the generalized speeds. These are in the same order as used in integration (as passed into evaluate_ode_function) and match the order of the mass matrix and forcing vector.

property times¶: A 1D ndarray of monotonic time values over which the equations of motion are numerically integrated. Can be set with an array-like for a shape(n,) array.

system¶

PyDy

Navigation

Related Topics