Inverted pendulum & experience in Acumen

Following standard procedure of classic mechanics, Mingkun set up the Lagrange equations, converted them to ODEs and give some arbitrary driving force to see how the system reacts.  Mingkun resorted to Mathematica to find the optimal force to balance the system under the constraint that minimal strength, minimal tilting, minimal horizontal shifting are considered.  In order to simulate one system with Acumen, Mingkun has to set up the equation of system, which is mostly differential equations. After that, converting them to ODE manually, or symbolic simplification is not very mentally challenging. In the end, Acumen is used to calculate a series of values for each system variable.  The most challenging part, Mingkun thinks, is the first phase, setting up system equations, and in this part, Mingkun could't get much help from Acumen'10 directly.

Parametric Acumen

Last week we managed to make Acumen parametric in the type used for real numbers. We have tested this version of Acumen with both Scala's Double and the library for constructive reals in Java by Boehm. We believe the current implementation provides us with a tool to evaluate accuracy and efficiency of different implementations of real numbers when building models and simulations. Immediate improvements to the current version include 1) parsing constants as written by the Acumen programmer directly into the type of numbers in use, avoiding the current conversion to Double first; and 2) finding a suitable interface for comparing numbers, for example with an extra argument for the tolerance.

On the sensitivity of hybrid models to sampling rate

The way Adam had chosen to evaluate the breaking ball signal function turned out to be the source of the strange behavior described in his paragraph from a week ago. The culprit is the Yampa function reactimate which allows interleaving I/O actions – in this case, rendering – with the sampling of a pure signal function. Such a construct is necessary when implementing an interactive application, but in this case (a deterministic simulation) it was doing more harm than good. Interactive applications (e.g. games) need to adjust the sampling rate with regard to the performance of the system they are executed on (affected by hardware and other applications running on the system) so that simulation time always progresses at the same rate compared to real time. The effects of this dynamic adjustment of the sampling rate can be seen in the figure below.

Above, the breaking ball model was modified not to add a new ball on each impact (reducing it to a bouncing ball simulation), to increase the number of bounces per unit of time the gravity constant was also increased, and the COR of the ball was set to 1. To emulate the effect of handling an exponentially increasing collection of signal functions, which was causing the sampling rate to decrease, Adam added a computationally intensive calculation to the list of actions performed in reactimate's input/sense argument, and toggled this on integer boundaries. The effects of this load on the system are apparent, and they manifest themselves in two ways. Firstly, the number of samples is clearly reduced as can be seen by the sparse distribution of data points along the plot in these intervals. Secondly, though the ball has a COR of 1, the ball also gains in elevation on each bounce in these intervals. The second effect is a consequence of the first, and is caused by the fact that when the sampling rate goes down, the accumulated effect of gravity on the position of the ball also decreases.

To evaluate a signal function with a pre-specified sampling rate, Yampa's embed function can be used. In this case the simulation ground to a halt (in contrast with last week's result) at around t = 4.5, the Zeno point of the sub-system of balls with lower COR.

Termination of interval Picard

During the past week Jan has been continuing to study the Picard method and its generalisation to the case of uncertain initial values. An interesting question that arises is how various representations of interval functions interact with time step and initial value precision in their effect on the termination of the method.

Yingfu's weekly paragraph

Performance issues
This week Yingfu started looking at the performance issues that came up while he was trying to run the competition last week. The problem is that the new version hangs when running the tournament code, and he found out the reason is because the program used up the heap. To solve that problem, Yingfu suggest to increase the heap, and reduce the sampling rate of 3D data (from 100 frame per second to 30).

Breaking ball in Yampa

After last week's initial attempt at implementing a breaking ball simulation in Yampa, Adam set out to familiarize himself with Yampa's switching combinators, starting with the basics. As mentioned in Henrik Nilsson's presentation accompanying the paper Functional Reactive Programming, continued, the simplest way to form a collection of signal functions (SFs) in Yampa is to use parB. This combinator uses the same input for all SFs in the collection, and also does not provide any way to modify the collection during execution. The pSwitchB combinator provides a way to update the collection (through its third parameter), but still uses the same input for each SF in the collection – the suffix B of the names of both of these combinators stands for "broadcast" and indicates the latter. The pSwitch combinator extends pSwitchB by providing a way to match (potentially unique) inputs with each signal function, through its first parameter – the routing function.

pSwitch :: Functor col =>

  -- Routing function, matches inputs with SFs in collection

  (forall sf . (a -> col sf -> col (b, sf)))

  -- Collection of SFs, used until next update

  -> SF (a, col c) (Event d)

  -- Produces an event if SF collection needs updating

  -> col (SF b c)

  -- Function to create updated collection of SFs  

  -> (col (SF b c) -> d -> SF a (col c))

  -- Resulting SF corresponding to pSwitch

  -> SF a (col c)

The system Adam set out to model was one where a ball is dropped from a height of 10 m (with a gravity constant of -10 m/s/s). Upon hitting the ground, the ball breaks into two new balls (bouncing with COR 0.75 and 0.5, respectively), and the process repeats.

Using pSwitch, it was possible to implement the simulation but, as can be seen in the plot below, a puzzling behavior appears in the output beginning at the Zeno point (~ 4.3) of the subsystem consisting of the "bottom" set of balls (corresponding to COR 0.5). Analytically, this is where the number of bounces approaches infinity. Adam will continue to investigate this behavior, he currently suspects that it has to do with the way Yampa samples the collection of SFs.

Interval arithmetic for the uninitiated

During the past week Jan has been working on identifying a self-contained set of definitions that is sufficient for a presentation of the interval Picard algorithm to an audience outside the interval methods community. Next, Jan will write a minimal Scala implementation of the identified operations to serve as reference for future work.

Interval inequalities and Picard continued

Jan has continued to implement the components of the Picard operator for interval polynomials. The additions are evaluation of a polynomial, either in one or all of its variables, as well as the primitive function operator. Once evaluation was there, the interval versions of order relations were easy to add.

One interesting aspect of the interval version of an order relation is that it becomes three-valued. This is because intervals represent any number they contain, thus only disjoint intervals can safely be said to be 'less than' or 'greater than' each other. To see the problem with intersecting intervals consider [1,5] < [3,7]. The interval [1,5] represents any value between 1 and 5, so also 4. If the number represented by [3,7] is 3, then the result should be 'true', if the number is 4 or above the result should be 'false'. Therefore it makes sense to return a third truth value for the interval relation, representing the lack of information. Scala provides a neat way of adding an extra value to a type through the Option type, and Option[Boolean] is therefore the right return type of interval relations, with None representing the lack of information, at least when they are used to represent real numbers.

The next steps towards an interval polynomial Picard operator are to add composition of polynomials. In the case of multivariate polynomials this requires a slightly more general type of composition as the variable to compose into needs to be specified. Luckily, the operation itself is fairly easy to implement as it is built up from the ring operations already present.

Built in Support scala.xml._

For the past few weeks Anil is exploring  XMT and the book  Programming Scala. The chapter 10 "Herding XML in scala" introduces the Scala's built in support for Reading,Exploring,Looping,Matching,Writing XML. He hope and believe that XML built in support from Scala could be used for framework XMT(eXtensible MPEG-4 Textual Format), for representing MPEG-4 scene description using textual syntax. In future he is planning to work with the XML built in support from Scala. 

Floating in Acumen

We have made progress towards making the current implementation of Acumen parametric over the type of real numbers. We have made the necessary changes in most of the functions that needed to be modified (functions dealing with GroundValues) and we have identified the classes that need to be made parametric. The key type is the type of GroundValues used in both the abstract syntax trees generated by the parser and in the values used by the interpreter. From the start it is the class AppModel that initiates parsing and interpretation. Results are then collected in a TraceModel and visualized with TraceView. Up to the TraceModel we will keep parametric types and functions, while in TraceView we have to fall back to Doubles in order to use Java's libraries for 2D-graphics.

A first encounter with pSwitch

A simulation of a breaking ball (a bouncing ball that splits in two on each bounce, see by Globally parallel, locally sequential Paul and Walid for details) requires that the simulator maintains a dynamic collection of objects representing the balls. AFRP supports this scenario through its parallel switching combinators, but initial attempts to express the breaking ball example using pSwitch proved unsuccessful. Adam will continue to familiarize himself with this API aiming for a minimal working example using one of the parallel switching combinators.

Interval polynomials in Scala

Interval arithmetic (IA) provides a way to perform numeric computations that yield safe approximations of the resulting values, i.e. ranges within which the result is guaranteed to lie. One application is for solving of initial value problems (IVPs), where the so-called Interval Picard method may be applied to compute bounds on the solutions of an IVP, even when the initial values are not exactly known, but lie in some interval.

One deficiency of IA is the so-called dependency problem, which arises when variables in expressions are replaced by intervals, effectively treating multiple occurrences of the same variable as distinct. To reduce such effects one can retain some relational information by using interval functions in place of plain intervals.

A rich but computationally tractable class of interval functions are polynomials with interval coefficients and as a step towards an interval IVP solver Jan has implemented the basic ring operations for interval polynomials in Scala.

acumen.util.Random

The acumen.util package contains the Random object, which ports part of the standard Haskell random generators to Scala. The generators are implemented using lazy values and may be used to produce streams of random Int or Double values that can be passed between objects, providing random values without the need for IO calls. The implementation currently suffers from a division by zero defect in the Int generator triggered when called from within the Double generator. The likely cause is that the Haskell version uses unbounded Integer values for most of the arithmetic while Random uses the Scala 32bit Int throughout. In the coming week Jan will make sure BigInt is used whenever Haskell uses Integer.

Grip on Scala

Anil is reading the book Programming Scala. He finished first 3 chapters and worked with the examples of the book  in eclipse IDE with Scala plugin. He is planning to work with all the chapters of the book to gain more grip on Scala.

Hybrid system modeling in AFRP

Continuing the exploration of Arrowized FRP, Adam implemented a model of a furnace (adapted from Example 15 in Michael S. Branicky's Introduction to Hybrid Systems) in Yampa. Together with the bouncing ball and breaking ball examples, these will serve as benchmarks a comparison of FRP languages as hosts for hybrid systems modeling.  Below is a comparison of Acumen and Yampa implementations of the furnace benchmark (note that rendering-specific parts of the FRP code have been omitted). Suggestions for improvements of the Yampa implementation are most welcome!

Acumen

Yampa

class Main(simulator)  private    state = "Off"; x  = 0; x' = 1  end  if state == "On"    if x >= 25      state = "Off"    else      x' [=] 40 - x    end;  end;  if state == "Off"    if x <= 21      state = "On"    else      x' [=] 10 - x    end;  end end

type State = (Bool, Double, Double) thermostat :: State -> SF () State thermostat i =  switch (aux i)  (\(mOut, xOut, vOut) ->           thermostat (not mOut, xOut, vOut)) aux :: State -> SF () (State, Event State) aux (isOn,xIn,x'In) = proc () -> do      rec x   <- (xIn +) ^<< integral -< x'      x'  <- returnA -< if isOn                        then 40 - x                        else 10 - x  tooHot  <- edge    -<     isOn && x >= 25  tooCold <- edge    -< not isOn && x <= 21  let s = (isOn, x, x')  returnA -< (s, (rMerge tooHot tooCold) `tag` s)