BIPOP

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 (nondynamical) 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.
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.
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 liftsoff, allowing to modify the differential equations
accordingly.
