Tutorial

Simulation.sci

Introduction

This function, described in chapter 5 of the User Documentation, does all the simulation job: it is interesting to understand what it does and how, but you do not need to modify this function, so you can skip this page if you want. Here is how it is called
[T, Q, QDOT, Z, ACTUATIONSTATE, CONTACTSTATE] = Simulation(InitialTime, FinalTime, SamplingPeriod, InitialPosition, InitialVelocity)
The names of the input variables are explicit, with position and velocity being of course joint position and joint velocity. Outputs are:

• T: a row vector containing the sampling times of the simulation
• Q: a matrix containing the joint positions at each of these times (each column gives one position)
• QDOT: same as Q but with the joint velocities
• Z: a matrix similar to Q but containing the dynamical state of the actuators, in the case of the muscle model for example (see section 1.2 of the User Documentation)
• ACTUATIONSTATE: a matrix similar to Z containing the discrete (non-dynamical) state of the actuators
• CONTACTSTATE: a boolean matrix describing the state of the contacts at each time. More details on how these contacts are dealt with can be found here

The function

Initialization

The first part, besides the creation and the initialization of all needed variables, creates the main variable called x. It contains in the same column vector the joint position and the joint velocity x = [q, qdot]'.

Calculation

For each sampling time you defined, the following scheme is followed:

• The Scilab dynamical system solver is called with the following syntax:
  [x, DetectedEvents, WarmInfo1, WarmInfo2] = ode('root', x, t, NextSamplingTime, CompleteDynamics, NContactEvents+NActuationEvents, EventDetection);
  it solves the problem xdot = CompleteDynamics(t,x) between t and NextSamplingTime (with its own internal time steps) and stop if one component of the vector returned by the EventDetection(t,x) function changes of sign.
• Two cases are then possible:

  - either there was no event, that is EventDetection never crossed 0 between t and NextSamplingTime:

  a variable called WarmStart is set to true and the ode function will be called with the same parameters at the next time step (it is done by passing the optional parameters Warminfo1 and WarmInfo2 in the ode function)

  - or EventDetection did cross 0, meaning that an event occurred in the simulation, that needs to be treated explicitly:

  then ode stops at a time t' between t and NextSamplingTime, the stopping time is given in DetectedEvents and replaces the present t in the next call to ode. Functions detecting the changes in actuation and contact states are called, and they update the value of ContactState and ActuationState. Then ode is called again but the CompleteDynamics function is different because it takes into account the modifications of state of the actuators and contacts.

• When NextSamplingTime is reached the results of computation are stored in the output variables and a new loop solving the problem for the next time step begins.

The CompleteDynamics function

This function is defined in the Kernel directory and returns xdot = [qdot, qddot]': the velocity qdot is already known, so only the acceleration qddot needs to be computed. This computation is described in section 1.1.5 of the User Documentation and follows these two steps:

• The call to ActuationDynamics

  defined in the ActuationModel/*Name* module/directory. This function is part of the actuation model, it computes the torques (Torques) to apply to the human model in order to follow the desired trajectory.

  The call to LagrangianDynamics

  defined in the LagrangianDynamics/*Name* module/directory. It computes the acceleration qddot by solving the quadratic problem with constraints given by equation 1.7 in the User Documentation:
  with contact constraints.
  The different terms appearing in this equation are provided in the definition of the LagrangianModel/*Name*: the inertia matrix M, the nonlinear effect matrix N and the gravity vector G. Only T(q)u corresponds to the Torques computed at the previous step.

The EventDetection function

This function defined in the Kernel directory returns a float vector Event = EventDetection(t, x). It relies on two functions: first the ActuationEventDetection defined in ActuationModel/*Name* directory that does nothing in our case because there is no actuator dynamics nor events considered. Then ContactEventDetection defined in the LagrangianDynamic/*Name* directory which computes the distance between the contact points and the environment, if the value is negative the point is currently below contact, if positive it is above. A change of sign in this vector helps detecting therefore impacts and lifts-off, allowing to modify the differential equations accordingly.