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