The MidpointVI class implements a variational integrator using a midpoint quadrature. The integrator works with any system, including those with constraints, forces, and kinematic configuration variables. The integrator implements full first and second derivatives of the discrete dynamics.
For information on variational integrators (including relevant notation used in this manual), see Discrete Dynamics and Variational Integrators.
The Midpoint Variational Integrator is defined by the approximations for \(L_d\) and \(f^\pm\):
where \(L\) and \(f\) are the continuous Lagrangian and generalized forcing on the system.
Create a new empty mechanical system. system is a valid System object that will be simulation.
solving the DEL equation to advance the integrator.
MidpointVI makes use of multithreading to speed up the calculations. num_threads sets the number of threads used by this integrator. If num_threads is None, the integrator will use the number of available processors reported by Python’s multiprocessing module.
Number of configuration variables in the system.
Number of dynamic configuration variables in the system.
Number of kinematic configuration variables in the system.
Number of input variables in the system.
Number of constraints in the system.
In continuous dynamics the System only calculates the second derivative, leaving the actual simulation to be done by separate numerical integrator. The discrete variational integrator, on the other hand, performs the actual simulation by finding the next state. Because of this, MidpointVI objects store two discrete states: one for the previous time step \(t_1\), and one for the new/current time step \(t_2\).
Also unlike the continuous System, these variables are rarely set directly. They are typically modified by the initialization and stepping methods.
State variables for the previous time. q1 is the entire configuration (dynamic and kinematic). p1 is discrete momemtum, which only has a dynamic component, not a kinematic component.
The input vector at \(t_1\).
The finite-differenced velocity of the kinematic configuration variables at \(t_2\) i.e. \(v_2=\frac{q_{k2}-q_{k1}}{t_2-t_1}\). If \(t_2=t_1\) the built-in Python constant None is returned.
The constraint force vector at \(t_1\). These are the constraint forces used for the system to move from \(t_1\) to \(t_2\).
Initialize the integrator from a known state (time, configuration, and momentum). This prepares the integrator to start simulating from time \(t_1\).
lambda1 can optionally specify the initial constraint vector. This serves at the initial guess for the simulation’s root solving algorithm. If you have a simulation that you are trying to re-simulate, but from a different starting time (e.g, you saved a forward simulation and now want to move backwards and calculate the linearization at each time step for a LQR problem.), it is a good idea to save lambda1 during the simulation and set it here. Otherwise, the root solver can find a slightly different solution which eventually diverges. If not specified, it defaults to the zero vector.
Initialize the integrator from two consecutive time and configuration pairs. This calculates \(p_1\) from the two pairs and initializes the integrator with the state (t1, q1, p1).
lambda1 is optional. See initialize_from_state().
Step the integrator forward to time t2 . This advances the time and solves the DEL equation. The current state will become the previous state (ie, \(t_2 \Rightarrow t_1\), \(q_2 \Rightarrow q_1\), \(p_2 \Rightarrow p_1\)). The solution will be saved as the new state, available through t2, q2, and p2. lambda will be updated with the new constraint force, and u1 will be updated with the value of u1.
lambda1 and q2 can be specified to seed the root solving algorthm. If they are None, the previous values will be used.
The method returns the number of root solver iterations needed to find the solution.
Raises a ConvergenceError exception if the root solver cannot find a solution after max_iterations.
Evaluate the DEL equation at the current states. For dynamically consistent states, this should be zero. Otherwise it is the remainder of the DEL.
The value is returned as a numpy array.
Calculate the first derivatives of \(q_2\) with respect to the previous configuration, the previous momentum, the input forces, or the kinematic inputs.
If both the parameters are None, the entire derivative is returned as a numpy array with derivatives across the rows.
If any parameters are specified, they must be appropriate objects. The function will return the specific information requested. For example, q2_dq1() will calculate the derivative of the new value of q with respect to the previous value of q1.
Calculating the derivative of \(q_2\) with respect to the i-th configuration variable \(\deriv[q_2]{q_1^i}\):
>>> q1 = system.configs[i]
>>> mvi.q2_dq1(q1=q1)
array([ 0.133, -0.017, 0.026, ..., -0.103, -0.017, 0.511])
Equivalently:
>>> mvi.q2_dq1()[:,i]
array([ 0.133, -0.017, 0.026, ..., -0.103, -0.017, 0.511])
Calculating the derivative of the j-th new configuration with respect to the i-th kinematic input:
>>> q = system.configs[j]
>>> k2 = system.kin_configs[i]
>>> mvi.q2_dk2(q, k2)
0.023027007157071705
# Or...
>>> mvi.q2_dk2()[j,i]
0.023027007157071705
# Or...
>>> mvi.q2_dk2()[q.index, k2.k_index]
0.023027007157071705
Calculate second derivatives of the new configuration.
If no parameters specified, the entire second derivative is returned as a numpy array. The returned arrays are indexed with the two derivatives variables as the first two dimensions and the new state as the last dimension.
To calculate the derivative of the k-th new configuration with respect to the i-th previous configuration and j previous momentum:
>>> q = system.configs[k]
>>> q1 = system.configs[i]
>>> p1 = system.configs[j]
>>> mvi.q2_dq1dp1(q, q1, p1)
6.4874251262289155e-06
# Or...
>>> result = mvi.q2_dq1dp1(q)
>>> result[i, j]
6.4874251262289155e-06
>>> result[q1.index, p1.index]
6.4874251262289155e-06
# Or....
>>> result = mvi.q2_dq1dp1()
>>> result[i, j, k]
6.4874251262289155e-06
>>> result[q1.index, p1.index, q.index]
6.4874251262289155e-06
Calculate the first and second derivatives of \(p_2\). These follow the same conventions as the derivatives of \(q_2\).
Calculate the first and second derivatives of \(\lambda_1\). These follow the same conventions as the derivatives of \(q_2\). The constraint dimension is ordered according to System.constraints.