The Modeling Monad Continuation

Earlier I wrote that I discovered a modeling monad for simulating dynamic systems in F#. I denoted the corresponded type by Dynamics<’a> and the corresponded workflow builder by dynamics. Recently I discovered that a monad transformer based on this modeling monad and the continuation monad is suitable for the discrete event-driven and process-driven simulation. The corresponded monad has type: DynamicsCont<’a> = Dynamics<’a -> unit> -> Dynamics<unit>. I denoted the corresponded workflow builder by dynamicscont.

The Dynamics Modeling Monad

I invented a monad which I called the Dynamics modeling monad. This monad allows the modeler to build a model based on differential equations and then simulate it as accepted in System Dynamics. Moreover, my method allows creating a hybrid model, where some parts can be defined as iterative processes, while the rest part will be described in terms of the equations. This is a very interesting mix of declarative and imperative approaches.

Earlier I wrote about the method in my blog in Russian. This message is mainly a translation of that information. Here I'll try to give a complete definition of the Dynamics modeling monad. Also based on this monad I created an open source simulation library written in F#. The library has name Aivika and it can be downloaded by the following link.

Let's consider the following simple dynamic system:

da/dt = -f;     a(t0) = 100;
db/dt = f - g;  b(t0) = 0;
dc/dt = g;      c(t0) = 0;
f = ka * a;
g = kb * b;
ka = 1;
kb = 1.

This system can be integrated with help of Euler's method and the Runge-Kutta method. To simulate, we have to define the initial time starttime, the final time stoptime and the main integration time step dt. The time scale is divided into nodes with step dt: starttime, starttime + dt, starttime + 2*dt, ..., stoptime. In case of the Runge-Kutta method we also create internal sub-nodes. As a result, each node is defined uniquely by a pair of two integers: the iteration number and a phase. Euler's method has only one phase. The 2nd order Runge-Kutta method has two phases. The 4th order method has four phases.

In F# we can introduce the corresponded types:

type Iteration = int
type Phase = int
type Method = Euler | RungeKutta2 | RungeKutta4

The model can be stochastic, i.e. the differential equations can contain random functions. In Aivika we can manage the random number generator with help of parameter that has the following type:

type Randomness = SimpleRnd | StrongRnd

Then the integration parameters, they are also called Simulation Specs, can be described by a value of the next type:

type Specs = {
   spcStartTime: float; spcStopTime: float; spcDT: float;
   spcMethod: Method; spcRandomness: Randomness
}

For our example we can take the following specs: the initial time is 0.0, the final time is 10.0, the integration time step equals 0.1, we apply the 4th order Runge-Kutta method and we use the strong pseudo-random number generator, although the latter is not applied in this particular model, but we have to define this parameter nevertheless.

let specs = {
   spcStartTime=0.0; spcStopTime=10.0; spcDT=0.1;
   spcMethod=RungeKutta4; spcRandomness=StrongRnd
}

We begin to describe the system with constants ka and kb.

let ka = 1.0
let kb = 1.0

Then we declare integrals and define their initial values.

let a = Integ 100.0
let b = Integ 0.0
let c = Integ 0.0

Each of the integrals has three properties: Value, Inflow and Outflow. These properties are constrained by the differential equations, one of which I will show on example of integral a:

d/dt (a.Value) = (a.Inflow - a.Outflow)

The same equations take place for integrals b and c as well. The Value property is read-only. It returns a dynamic process that represents the current value of the integral. The Inflow and Outflow properties are settable. They both define a derivative of the integral. They are also dynamic processes. By default, these two properties Inflow and Outflow represent a zero value. If we modify none of them, then the derivative will be zero too.

After the integrals are declared we define the auxiliary variables that may depend on these integrals.

let f = ka * a.Value
let g = kb * b.Value

Finally, we complete the model definition writing the differential equations themselves.

a.Inflow <- -f
b.Inflow <- f - g
c.Inflow <- g

This is all we need to define the model. Now we can simulate it and receive the results for integral a:

> runDynamics a.Value specs;;
val it : float [] =
 [|100.0; 90.48375; 81.87309014; 74.0818422; 67.03202889; 60.65309344;
   54.88119344; 49.65856187; 44.93292897; 40.65699912; 36.78797744;
   33.28714154; 30.11945393; 27.2532114; 24.65972767; 22.31304633;
   20.18968106; 18.26838054; 16.52991577; 14.95688766; 13.53355284;
   12.24566612; 11.08033792; 10.02590526; 9.071815051; 8.208518451;
   7.427375314; 6.720567711; 6.081021686; 5.50233646; 4.978720367;...|]

We can also ask Aivika to return the results for more number of variables including the current simulation time.

> runDynamics Dynamics.ofArray [| time; a.Value; b.Value; c.Value |] specs;;
[|[|0.0; 100.0; 0.0; 0.0|];
 [|0.1; 90.48375; 9.048333333; 0.4679166667|];
 [|0.2; 81.87309014; 16.37454263; 1.752367234|];
 [|0.3; 74.0818422; 22.22445032; 3.693707481|];
 [|0.4; 67.03202889; 26.81268809; 6.155283021|];
 [|0.5; 60.65309344; 30.32640707; 9.020499487|]; ...|]

As regards to me, for the first time it created a great impression upon me by its conciseness. It looked fantastic that I could define and simulate the dynamic system in a general purpose programming language using only a few lines of code!

The secret of this method is in a new monad that I discovered. This monad is a variation of the idiomatic Reader monad. If you remember, the latter represents a computation that calculates something dependent on the external environment variable that is passed when we want to run the computation. I extended this idea to represent the dynamic processes.

Thus, we simulate some dynamic process that returns a value of generic type 'a at each integration time. Let's denote the corresponded polymorphic type as Dynamics<'a>:

type Dynamics<'a> (iter: Iterator<'a>) =
   member d.Iterator = iter

let iterator (p: Dynamics<'a>) = p.Iterator

This type is just a wrapper for the function that calculates a value at each integration node by the specified simulation specs. If you remember, the integration nodes are defined uniquely by the iteration number and phase.

type Iterator<'a> = Run -> Specs -> Iteration -> Phase -> 'a

Here we can see a new type Run, which is responsible for the non-determinism and that opens door for such techniques as caching. This is not pure functional thing, which makes this monad unusual. After removing some small optimizations, the type looks like this

type Run () =
   let mutable disposed = false
   let disposedEvent = new Event<_>()
   member x.Disposed = disposedEvent.Publish
   member x.IsDisposed = disposed
   interface IDisposable with
       member x.Dispose() =
           if not disposed then
               disposedEvent.Trigger()
               disposed <- true

Each time we start a simulation, we create a new instance of the Run class. During the simulation we pass this object to the dynamic process. Using the Run object, any function can create internal intermediate structures that will contain actual data while this object is alive. For example, the integral function can create a cache-table of the previous values. After the simulation is finished, the Run object is destroyed. But before it the Dispose method is called and the Disposed event is raised. The integral function or another function on more low level catches this event and destroys its cache table. It allows us to achieve the following goal.

Different simulation runs of the dynamic process can return different results. But inside one run we can memorize and use the previous values. At the end of the run we can remove all intermediate data structures. It allows us to use caching and other similar techniques.

The run function for the dynamic process is very simple. We applied it already in the examples above.

let runDynamics d specs = [|
  
   use r = new Run ()
  
   let f1 = iterator d r specs
   let f2 = fun n -> f1 n 0
  
   for i = 0 to (Specs.iterations specs - 1) do yield f2 i
|]

Here the iterations function returns the total number of iterations by the specified simulation specs. We are interested only in those values of the dynamic process, which are calculated in the main integration nodes. Therefore the zero phases are requested. Below is stated the definition of the integ function that returns an integral by the specified derivative and initial value. You can look at this function to gain some insight about the meaning of phases. Also the next value that returns the current simulation time can give an idea:

let time = Dynamics<float> (fun r s n ph ->
   let delta m ph =
       match m, ph with
           | Euler, 0 -> 0.0
           | RungeKutta2, 0 -> 0.0
           | RungeKutta2, 1 -> s.spcDT
           | RungeKutta4, 0 -> 0.0
           | RungeKutta4, 1 -> s.spcDT / 2.0
           | RungeKutta4, 2 -> s.spcDT / 2.0
           | RungeKutta4, 3 -> s.spcDT
           | _  -> failwith "Incorrect method and/or phase"
   in s.spcStartTime + float (n) * s.spcDT + delta s.spcMethod ph)

The rest part of the Dynamics modeling monad can be now defined easily. The return and bind functions are almost the same as in case of the Reader monad.

type DynamicsBuilder() =

   member d.Return (a) =
       new Dynamics<'a> (fun r s n ph -> a)
      
   member d.Bind (m, k) =
       new Dynamics<'b> (fun r s n ph ->
           let a = iterator m r s n ph
           let m' = k a
           iterator m' r s n ph)

Now we can define the dynamics workflow. It will allow us to use a syntactic sugar similar to the do-notation from Haskell.

let dynamics = new DynamicsBuilder()

Following the approach described in excellent book "The Haskell School of Expression" by Paul Hudak, we can define arithmetic operations and functions like sin and cos in accordance with the next pattern:

type Dynamics<'a> with
   static member (+) (a: Dynamics<float>, b: Dynamics<float>) =
       Dynamics.lift2 (+) a b

Here the lift2 function is rather idiomatic. It takes some binary function and two values wrapped in the monad. We extract these values, apply the function to them and then wrap the result in the monad.

module Dynamics =
   let lift2 f m1 m2 = dynamics {
       let! a1 = m1
       let! a2 = m2
       return f a1 a2
   }

But I actually applied in Aivika a more optimal method that uses a knowledge about the internal structure of the monad:

module Dynamics =
   let lift2 f m1 m2 = Dynamics<_> (fun r s n ph ->
       let a1 = iterator m1 r s n ph
       let a2 = iterator m2 r s n ph
       in f a1 a2)

Now we can define the integral and random functions. Also we can add the memo function that caches the previous values of the specified dynamic process, calculating these values strictly sequentially by the integration nodes, which is important for calculation of the integral, for example. The random functions also use a similar caching technique.

module Dynamics =
   val memo: Dynamics<'a> -> Dynamics<'a>

For brevity I will give only a definition of the integ function that returns an integral by the specified derivative and initial value. This definition contains only an implementation of the 4th order Runge-Kutta method. Other methods are simpler and they are implemented in a similar way.

let integ (f: Dynamics<float>) (i: Dynamics<float>) =
  
   let integRK4 (y: Dynamics<_>) (f: Dynamics<_>) (i: Dynamics<_>) r s =
       let vy n = iterator y r s n 0
       let vi n = iterator i r s n 0
       let k1 n = iterator f r s n 0
       let k2 n = iterator f r s n 1
       let k3 n = iterator f r s n 2
       let k4 n = iterator f r s n 3
       fun n ph ->
           match n, ph with
               | 0, 0 -> vi 0
               | n, 0 -> vy (n-1) + s.spcDT/6.0 *
                           (k1 (n-1) + 2.0 * k2 (n-1) + 2.0 * k3 (n-1) + k4 (n-1))
               | n, 1 -> vy n + s.spcDT/2.0 * k1 n
               | n, 2 -> vy n + s.spcDT/2.0 * k2 n
               | n, 3 -> vy n + s.spcDT * k3 n
               | _ -> failwithf "integRK4: incorrect phase = %i" ph

   let rec y = Dynamics.memo z
   and z = Dynamics<_> (fun r s ->
       match s.spcMethod with
           | Euler -> integEuler y f i r s
           | RungeKutta2 -> integRK2 y f i r s
           | RungeKutta4 -> integRK4 y f i r s)
   in y

Unfortunately, the integ function is not sufficient for all cases. The derivative can recursively depend on the integral itself. Therefore I applied a trick and added the Integ type that has properties Value, Inflow and Outflow, although we still have to define the auxiliary variables before differential equations as it was shown in the example above. This is a price of the programming language with a strict evaluation model. For example, in a lazy language the integ function would be sufficient and there would no reasons to define the Integ type. But laziness implies purity, while the Dynamics modeling monad is not pure at all.

We saw that the differential equations are defined declaratively enough. We say what to simulate rather than how to simulate. Only in the strict language we define the order of dependencies between the variables and integrals, which is necessary if the model contains loopbacks.

Nevertheless, during a simulation of such differential equations we use the memo function under the hood. This function is exceptionally imperative. It takes some dynamic process, i.e. a value in the Dynamics modeling monad, and returns a duplicate of this process, where all calculations are performed strictly sequentially by the integration nodes. It allows us to create hybrid models as we know that the first cached process will be calculated in a strict order, if there are no other references to it. Such a process can be implemented as iterative.

val iterative: ('a -> 'b) -> Dynamics<'a> -> Dynamics<'b>
let iterative f m = map f m |> memo

Here the iterative function takes function f with some side-effect, a source of the arguments for that second function and returns already an iterative process generated by function f. The map function is an idiomatic functor map function, which can be optimized for this monad.

let map f m = Dynamics<_> (fun r s n ph ->
   let a = iterator m r s n ph
   in f a)

With help of the Dynamics modeling monad we can simulate the difference equations too. It is even simpler than the differential equations. Also we can simulate the conveyor stock. I think that we should use for the last case such an immutable data structure as heap. We have to keep all the history of the conveyor and there is nothing better than a structure that would use shared data among the different integration nodes. The heap allows us to receive relatively fast the next value of the conveyor, namely a full value, not just a modified state.

As regards to the speed of simulation, Aivika is slow now. I compared it with simulation tool MapSim, which I developed too. MapSim compiles a simulation into a fast byte-code which is executed with more high speed. But I think that the main strength of my new method is that the Dynamics modeling monad allows the modeler to build easily hybrid models of any complexity, which would be almost impossible or hard to implement in more fast but specialized software tools like MapSim.