Localizing strategy continued

Last week Walid, Michal and Jan began working on a top-level strategy for the enclosure interpreter. Jan's original idea of reusing the constraint solver developed for the event detection and invariant intersection algorithms turned out not be suitable in its naive form. The problem with the idea is that it would not be able to ignore any events, and therefore not be able to localize the first event effectively. The solution is instead to use a subdivision-based approach that separately localizes the left and right endpoints of a segment of time that may contain the next event following. The segment is analyzed to determine if an event has with certainty been located and if continued sub-division is likely to localize the event further. Currently, an implementation is being finalized and integrated with the Acumen IDE. 

Scala course at Chalmers

During the past week Adam has attended an intensive Scala course at Chalmers, brilliantly presented by David Raymond Christiansen of ITU. The course introduced several aspects of advanced Scala programming, focusing on how its type system can be used with various approaches to embedding DSLs in Scala. Particular weight was given to Scala's implicit conversions mechanism, which e.g. makes it possible to emulate Haskell's type classes (with some additional flexibility) in Scala and can make interoperation with Java libraries smoother. The material of the course built up to the topic covered during the last day, a technique and framework called Lightweight Modular Staging (LMS). Adam found it easier to grasp after having understood the polymorphic DSL embedding technique, covered earlier in the week, of which LMS is an instance. In brief, this technique is based on defining the DSL as operations on an abstract data type. A program in the DSL is a nested call to these operations, and interpreters for the DSL are defined as implementations of the abstract operations which make up the program.

Progression of Zeno systems

During the past week Adam has been implementing a series of models developed by Walid, Jan and Adam, evaluating the enclosure semantics' performance on Zeno systems. An example of such a model is the following, inspired by a formulation of Zeno's paradox:

class Main(simulator)

  private

    d := 1; d' := -1;

    s := 1/2;

    mode := "q"; 

  end

  switch mode assume s <= abs(d)

    case "q"

      if s == d 

        s := s/2; 

        mode := "q" 

      end;

      d' = -1;

  end;

end

Adding a mode invariant (using the assume keyword) enables the interpreter to prune enclosures of impossible solutions, yielding a better overall enclosure. Though, at first glance, s <= d seems to be a suitable invariant for this system, this invariant does not hold after d crosses zero – changing this to s <= abs(d) ensures that the invariant holds for the entire simulation time. 

Adam will spend the coming week participating in a one-week Scala course given at Chalmers, learning about things such as implicit conversions, API design, polymorphic embedding of DSLs, scala-virtualized and LMS.

Towards a localizing strategy

Last week Jan integrated the new atomic step with the enclosure interpreter and started working on the top-level strategy. The idea has been to use the same constraint propagation technique employed to intersect enclosures with mode invariants and event guards during event handling. By composing event guards with the enclosures computed for the free variables in the guard and contracting the domain of the resulting predicate it should be possible to contract time segments around sub-segements where an event may occur. The advantage with this approach is that much of the necessary infrastructure is already in place, but a disadvantage is that it generally cannot show that an event has happened with certainty. 

Simpler atomic step

Last week Walid, Michal and Jan worked on a simplified version of the atomic step, i.e. the computation of an enclosure for a hybrid ODE-IVP over a fixed time segment. Jan implemented the specification and wrapped it in a function that makes it work as a drop-in replacement for the old atomic step, a.k.a. solve-VtE. Currently, the new atomic step produces a single constant enclosure for each variable, i.e. it cannot represent non-contiguous enclosures, such as those arising from discrete transitions in the simulation. This limitation can be lifted by changing the representation of the system state used by the new atomic step, from a single interval per variable and mode, to a set of intervals per variable and mode.

Reviewing literature, implementing interval π, and planning the upcoming year

During the past week Adam has continued reviewing literature on reachability analysis and also worked on developing example models illustrating Zeno behaviour in Acumen. Among these is one of a "bouncing inverted pendulum", where a mass is attached to a stiff massless rod which in turn is attached to the ground, constraining the motion of the mass as it falls and bounces against the ground. Adam also added support for constants to Acumen, used this to add π to all interpreters, and worked on developing a study plan for 2013 which, among other things, is to include the research questions that he will work on during the upcoming year.

Testing, examples, a merged plot and a highlighting editor

During the past week Adam continued reviewing literature on reachability analysis of hybrid systems, implemented a couple of new example models (e.g. one from optimal control theory used by Fuller in 1960 to illustrate Zeno behaviour), added generators and property-based tests for the enclosure semantics, developed a merged plot view for the new plotter and added syntax highlighting (based on RSyntaxTextArea) to the Acumen IDE code editor.

Polynomial intervals, not interval polynomials!

Following Walid's suggestion, Jan implemented splitting of initial condition boxes during IVP solving as a way to obtain narrowing enclosures for stable ODEs, such as x' = -x and x'' = -x'-x. Jan noticed that the current representation of solutions for IVPs as linear interval polynomials does allow for representation of narrowing enclosures. This forces the enclosures to widen between segments and prevents the piece-wise solution from converging. The way to allow narrowing enclosures while retaining a piece-wise affine representation of enclosures is to use pairs of linear interval polynomials to bound the solution.

Debugging using runtime assertions

During the past week Jan presented the ongoing work on Acumen's enclosure semantics at the CERES Open Day. Jan has also been working on debugging the enclosure interpreter. Walid helped to identify and Jan corrected a bug in the way that default field components were handled by the enclosure interpreter.  Adam and Jan identified and corrected a bug in the way the enclosure interpreter was identifying and handling events that lead to inconsistent states. The highly effective method was to inject runtime assertions that encapsulate invariants of the program at boundaries between sub-programs of the interpreter. This way the offering states could easily be reported as assertion failures without the need to sift through what otherwise tend to be long and cluttered logs. 

New plotting, debugging by contract and reachability analysis literature

During the past week Adam has worked on integrating a basic JFreeChart-based plotter into the Acumen IDE with Kevin and added PDF export of such plots. He also reviewed literature on reachability analysis tools and found the recent paper Taylor Model Flowpipe Construction for Non-linear Hybrid Systems by Chen, Ábraahám and Sankaranarayanan especially interesting. While working with Jan on debugging the enclosure interpreter, Adam found that using Scala's require statements (runtime assertions) as contracts for functions simplified the debugging process considerably.

The Acumen-Ptolemy II integration continues

The Acumen-Ptolemy II integration continues 

During the last week, the master student Carlos Fuentes has joined the group, and has set the environment to continue with he integration between Acumen and Ptolemy II project. Currently, he is studying Ptolemy papers to analyse the different models of computation, and how to face the creation of an agent between acumen and those models.

Smart grid models in Acumen

Anita has started collecting and documenting models of several smart grid components. The main goal is to be able to simulate smart micro-grids by connecting several of these component models. This collection of models can also be used to compare Acumen to other tools such as Simulink and Modelica.

Some interesting components for smart grids are illustrated below. The figure shows an electric motor connected to a voltage source which is, in turn, regulated by a PI controller in a current feedback loop. The desired current output is the square signal h(t). The corresponding Acumen code can be found below. This example illustrates how different electromechanical components, sensors and controllers may be modeled individually and connected to form a micro grid.



class RLcircuit(R, L)
  private i = 0; i' = 0; v = 0; e = 0; end
  i' [=] (v - e - R*i) / L;
end

class ElectricMotor(Jeq, kE, kT, TL)
  private i = 0; wm = 0; wm' = 0; e = 0; end
  wm' [=] (kT*i - TL) / Jeq;
  e     [=] wm*kE;
end

class PIcontroller(kP, kI, kPWM)
  private h = 0; int_h = 0; int_h' = 0; i = 0; int_i = 0; int_i' = 0; v = 0; end
  int_i'  [=] i;
  int_h' [=] h;
  v       [=] kPWM*(kP*(h - i) + kI*(int_h - int_i));
end

class Main(simulator)
  private
    motor      = create ElectricMotor(0.0002, 0.1, 0.1, 0);
    circuit     = create RLcircuit(4, 2);
    controller = create PIcontroller(50, 100, 6);
    mode = 1; t = 0; t' = 1;
  end
 
  t' [=] 1;
  controller.i [=] circuit.i;
  circuit.v     [=] controller.v;
  circuit.e    [=] motor.e;
  motor.i      [=] circuit.i;

  simulator.timeStep [=] 0.001;
  simulator.endTime  [=] 1;

  switch mode
    case 1
    controller.h [=] 0;
    if t >= 0.1 mode = 2; end

    case 2
    controller.h [=] 10;
    if t >= 0.2 mode = 3; end

    case 3
    controller.h [=] 0;
  end

end











--
Anita Pinheiro Sant'Anna

Compass Gait in Acumen

Anita just got back from a three-week trip to the US, during which time she visited the Power Systems Lab at Portland State University, the MAHI Lab at Rice University, the Nuclear Robotics Group and the Human Centered Robotics Lab at the University of Texas at Austin, and the AMBER Lab at Texas A&M University.

During her visit to the AMBER Lab, Anita and Ryan Sinnet implemented a passive compass gait model in Acumen, and used the Acumen3D extension to generate a 3D animation. A screen shot of the animation and the corresponding Acumen code can be found below.

class Main(simulator)
  private
    q1=0; q2=0; q1'=-0.5; q2'=2; q1''=0; q2''=0; //theta1 and theta2 angles
    L=1; M=10; m=5; pi=3.1416; g=9.81; //model constants
    s=0.06; //slope of the surface
    pm1x=0; pm1z=0; pM2x=0; pM2z=0; pm3x=0; pm3z=0;
    p2z=0; p2x=0; //position of stance foot
    sfx=0; sfz=0; //position variable for animation
    h=0; h'=0; //guards

    _3D=[["Sphere",[0,0,0],0.01*5,[1,0,0],[0,0,0]],
         ["Sphere",[0,0,0],0.01*10,[1,0,0],[0,0,0]],
         ["Sphere",[0,0,0],0.01*5,[0,1,0],[0,0,0]],
         ["Cylinder",[0,0,-1/2],[0.01,1],[0,0,1],[3.1416/2,0,0]],
         ["Cylinder",[0,0,-1/2-1/2],[0.01,1],[0,0,1],[3.1416/2,0,0]],
         ["Box",[0,0,0],[10,2,2],[0.5,0.5,0],[3.1416/2,0.06,0]]];
  end

  //DINAMICS EQUATIONS
  //theta 1: angle between vertical and stance foot
  q1''[=](2*(2*q1'*q2'*L*m*sin(q2) + q2'^2*L*m*sin(q2) + q1'^2*L*m*(sin(q2) - sin(2*q2)) + g*(-2*(m + M)*sin(q1) + m*sin(q1 + 2*q2))))/(L*(-3*m - 4*M + 2*m*cos(2*q2)));

  //theta 2: angle between legs
  q2''[=](-2*(-2*q1'^2*L*(-3*m - 2*M + 2*m*cos(q2))*sin(q2) +      2*q1'*q2'*L*m*(sin(q2) - sin(2*q2)) + q2'^2*L*m*(sin(q2) - sin(2*q2)) +      g*(-2*(m + M)*sin(q1) + (3*m + 2*M)*sin(q1 - q2) - 2*m*sin(q1 + q2) -        2*M*sin(q1 + q2) + m*sin(q1 + 2*q2))))/   (L*(-3*m - 4*M + 2*m*cos(2*q2)));

  //ADDITIONAL EQUATIONS FOR ANIMATION
  //position of link 1 - middle of stance leg
  pm1x[=] -(L*sin(q1))/2;
  pm1z[=] (L*cos(q1))/2;

  //position of link 2 - hip joint
  pM2x[=]-L*sin(q1);
  pM2z[=]L*cos(q1);

  //position of link 3 - middle of swing leg
  pm3x[=](L*(-2*sin(q1) + sin(q1 + q2)))/2 ;
  pm3z[=]L*cos(q1) - (L*cos(q1 + q2))/2 ;

  //reference position of stance foot - just for animation
  p2x[=]L*(-sin(q1) + sin(q1 + q2));
  p2z[=]L*(cos(q1) - cos(q1 + q2));
 
   _3D[=][["Sphere",[sfx+pm1x-2,0,sfz+pm1z],0.01*m,[1,0,0],[0,0,0]],
         ["Sphere",[sfx+pM2x-2,0,sfz+pM2z],0.01*M,[1,0,0],[0,0,0]],
         ["Sphere",[sfx+pm3x-2,0,sfz+pm3z],0.01*m,[0,1,0],[0,0,0]],
         ["Cylinder",[sfx+pm1x-2,0,sfz+pm1z],[0.01,L],[0,0,1],[3.1416/2,q1,0]],
         ["Cylinder",[sfx+pm3x-2,0,sfz+pm3z],[0.01,L],[0,0,1],[3.1416/2,q1+q2,0]],
         ["Box",[0-2,0,0],[20,0.01,2],[0.5,0.5,0],[3.1416/2,-s,0]]];

  //GUARD
  h[=]2*L*(1/cos(s))*sin(q2/2)*sin(q1 + q2/2 + s);
  h'[=]2*(q1' + q2'/2)*L*cos(q1 + q2/2 + s)*(1/cos(s))*sin(q2/2) +  q2'*L*cos(q2/2)*(1/cos(s))*sin(q1 + q2/2 + s);

  if h <= 0  &&  h' < 0 && q2>0.05

    //RESET MAP
    //reset for dinamics equation
    q1=q1+q2;
    q2=-q2;
    q1'=(q2'*m + q1'*(m - 2*(m + 2*M)*cos(q2)))/(-3*m - 4*M + 2*m*cos(2*q2));
    q2'=(q2'*m*(-1 + 2*cos(q2)) + 8*q1'*(m + M)*(1 + 2*cos(q2))*sin(q2/2)^2)/   (-3*m - 4*M + 2*m*cos(2*q2));
   
    //just for animation
    sfx=sfx+p2x;
    sfz=sfz+p2z;

  end
end