Русский     English

Aivika Tutorial for Beginners

Here we will take a simple discrete event simulation model. We will run it to see the results in final time point. Then we will add a simulation experiment. We will instruct Aivika to save data in the CSV files. All the results will be automatically collected together in a single HTML file, which we can open with help of the Internet browser. We will define a Monte-Carlo simulation. We will launch 1000 simulation runs. The Aivika libraries will plot histograms and a deviation chart with confidence intervals by rule 3-sigma. As before, the results will be collected in the HTML file for analysis and observing.

Simulation Model

The prerequisite is that we have to install the aivika package. It will work on Windows, Linux and OS X.

$ cabal install aivika

We will take model MachRep2 described in document Introduction to Discrete-Event Simulation and the SimPy Language.

The model description is as follows.

Two machines, but sometimes break down. Up time is exponentially distributed with mean 1.0, and repair time is exponentially distributed with mean 0.5. In this example, there is only one repairperson, so the two machines cannot be repaired simultaneously if they are down at the same time.

Let us find the long-run proportion of up time. Also let us find the long-run proportion of the time that a given machine has immediate access to the repairperson when the machine breaks down. Output values should be about 0.6 and 0.67.

In Aivika it can be written in the following way.

import Control.Monad
import Control.Monad.Trans

import Simulation.Aivika

meanUpTime = 1.0
meanRepairTime = 0.5

specs = Specs { spcStartTime = 0.0,
                spcStopTime = 1000.0,
                spcDT = 1.0,
                spcMethod = RungeKutta4,
                spcGeneratorType = SimpleGenerator }
     
model :: Simulation Results
model =
  do -- number of times the machines have broken down
     nRep <- newRef 0 
     
     -- number of breakdowns in which the machine 
     -- started repair service right away
     nImmedRep <- newRef 0
     
     -- total up time for all machines
     totalUpTime <- newRef 0.0
     
     repairPerson <- newFCFSResource 1
     
     let machine :: Process ()
         machine =
           do upTime <-
                randomExponentialProcess meanUpTime
              liftEvent $
                modifyRef totalUpTime (+ upTime) 
              
              -- check the resource availability
              liftEvent $
                do modifyRef nRep (+ 1)
                   n <- resourceCount repairPerson
                   when (n == 1) $
                     modifyRef nImmedRep (+ 1)
                
              requestResource repairPerson
              repairTime <-
                randomExponentialProcess meanRepairTime
              releaseResource repairPerson
              
              machine

     runProcessInStartTime machine
     runProcessInStartTime machine

     let upTimeProp =
           do x <- readRef totalUpTime
              y <- liftDynamics time
              return $ x / (2 * y)

         immedProp :: Event Double
         immedProp =
           do n <- readRef nRep
              nImmed <- readRef nImmedRep
              return $
                fromIntegral nImmed /
                fromIntegral n

     return $
       results
       [resultSource
        "upTimeProp"
        "The long-run proportion of up time (~ 0.6)"
        upTimeProp,
        --
        resultSource
        "immedProp"
        "The proption of time of immediate access (~0.67)"
        immedProp]
  
main =
  printSimulationResultsInStopTime
  printResultSourceInEnglish
  model specs

Here the model value returns the simulation results within computation, but the main entry runs the simulation and prints the results in final time point.

Let it be file MachRep2.hs. Let us run it.

The output depends on the random number generator. We will see something similar to

$ runghc MachRep2.hs 
----------

-- simulation time
t = 1000.0

-- The long-run proportion of up time (~ 0.6)
upTimeProp = 0.6030071376098631

-- The proption of time of immediate access (~0.67)
immedProp = 0.6839032527105922

We received the expected results in the final time point. After the model seems to be valid, we can proceed to creating a simulation experiment for more detailed analysis.

Simulation Experiment

To create and run simulation experiments, we have to install the following packages. It will work on Windows, Linux and OS X.

$ cabal install aivika aivika-experiment

Let us begin with a simple experiment. We will save our results in time points that are called integration time points in Aivika. They are defined with help of simulation specs. This is value specs in the code. We start from time 0.0 to 1000.0 with step 1.0. We still have to define the integration time points even if none integration method is actually used in our model. Here we can use these points for saving the results in the CSV files.

We already return the simulation results with help of value model. Therefore, let us copy file MachRep2.hs to a new file Model.hs. Now we have to remove the main entry from it. We will have another main in another file. Also let us remove the specs value from the new file. We will provide with the simulation specs in another file too. It should be useful to separate what is related to the model only. Finally, we have to put the code in a new module Model.

So, we will receive the following file Model.hs. This is a file with the simulation model defined.

module Model (model) where

import Control.Monad
import Control.Monad.Trans

import Simulation.Aivika

meanUpTime = 1.0
meanRepairTime = 0.5

model :: Simulation Results
model =
  do -- number of times the machines have broken down
     nRep <- newRef 0 
     
     -- number of breakdowns in which the machine 
     -- started repair service right away
     nImmedRep <- newRef 0
     
     -- total up time for all machines
     totalUpTime <- newRef 0.0
     
     repairPerson <- newFCFSResource 1
     
     let machine :: Process ()
         machine =
           do upTime <-
                randomExponentialProcess meanUpTime
              liftEvent $
                modifyRef totalUpTime (+ upTime) 
              
              -- check the resource availability
              liftEvent $
                do modifyRef nRep (+ 1)
                   n <- resourceCount repairPerson
                   when (n == 1) $
                     modifyRef nImmedRep (+ 1)
                
              requestResource repairPerson
              repairTime <-
                randomExponentialProcess meanRepairTime
              releaseResource repairPerson
              
              machine

     runProcessInStartTime machine
     runProcessInStartTime machine

     let upTimeProp =
           do x <- readRef totalUpTime
              y <- liftDynamics time
              return $ x / (2 * y)

         immedProp :: Event Double
         immedProp =
           do n <- readRef nRep
              nImmed <- readRef nImmedRep
              return $
                fromIntegral nImmed /
                fromIntegral n

     return $
       results
       [resultSource
        "upTimeProp"
        "The long-run proportion of up time (~ 0.6)"
        upTimeProp,
        --
        resultSource
        "immedProp"
        "The proption of time of immediate access (~0.67)"
        immedProp]

Now we should define the experiment itself that says what to simulate, what should be done with the results of simulation, how we should process the results. Here an idea is that the model and experiment are separated. The model provides with data. The experiment gets data and does something with them.

module Experiment (experiment, generators) where

import Data.Monoid

import Simulation.Aivika
import Simulation.Aivika.Experiment

specs = Specs { spcStartTime = 0.0,
                spcStopTime = 1000.0,
                spcDT = 1.0,
                spcMethod = RungeKutta4,
                spcGeneratorType = SimpleGenerator }

experiment :: Experiment
experiment =
  defaultExperiment {
    experimentSpecs = specs,
    experimentRunCount = 3,
    experimentDescription = "Experiment Description" }

upTime = resultByName "upTimeProp"
immedProp = resultByName "immedProp"

generators :: [WebPageGenerator r]
generators =
  [outputView defaultExperimentSpecsView,
   outputView defaultInfoView,
   outputView $ defaultTableView {
     tableSeries = upTime <> immedProp }]

Let it be file Experiment.hs.

Here we define the specs, experiment and generators. The generators say that we want to show the specs, get the detailed information about simulation data and save the data by specified variables in the CSV files. For demonstrating purposes, we will launch three simulation runs only, but their number can be much greater. Aivika quite easily and quickly processes thousands of runs for this simple model.

Also please note how we can combine different result sources, when defining the table series. This is a general rule. Almost everywhere we can combine different sources to receive a compound source.

Having the model and experiment provided, now we have to define a launcher. Let it be file Main.hs.

import Simulation.Aivika.Experiment

import Model
import Experiment

main = runExperiment experiment generators (WebPageRenderer ()) model

You can find all files for this simulation experiment in the following Zip archive.

Now we will start the launcher. We will receive something that looks like this. It took about 1 second on my MacBook Pro to finish.

$ runghc Main.hs 
Updating directory experiment
Generated file experiment/Table(1).csv
Generated file experiment/Table(2).csv
Generated file experiment/Table(3).csv
Generated file experiment/index.html

It means that Aivika created a new directory experiment with three CSV files for each simulation run and it generated an HTML file.

In my case I received the following HTML file that you can open in your Internet browser.

The CSV files are fine if you are going to analyze them in the R statistics tool, for example. But Aivika can plots charts and histograms itself, which can be useful for quick analysis too. Let us see how we can create them.

Charting on Windows

To plot the charts, we have to use a charting back-end. At present, there is only one option available for Windows: we have to use the Diagrams-based charting back-end.

$ cabal install Chart Chart-diagrams
$ cabal install aivika aivika-experiment aivika-experiment-chart aivika-experiment-diagrams

It will work for happy users of Linux and OS X too, but they can also use the Cairo-based charting back-end that usually works significantly faster and creates files of less size. The Cairo-based charting back-end is used below.

The Model.hs file remains the same. It already provides with all data we need for charting.

The Experiment.hs file has changed.

module Experiment (experiment, generators) where

import Data.Monoid

import Simulation.Aivika
import Simulation.Aivika.Experiment
import Simulation.Aivika.Experiment.Chart

specs = Specs { spcStartTime = 0.0,
                spcStopTime = 1000.0,
                spcDT = 1.0,
                spcMethod = RungeKutta4,
                spcGeneratorType = SimpleGenerator }

experiment :: Experiment
experiment =
  defaultExperiment {
    experimentSpecs = specs,
    experimentRunCount = 1000,
    experimentDescription = "Experiment Description" }

upTime = resultByName "upTimeProp"
immedProp = resultByName "immedProp"

generators :: ChartRendering r => [WebPageGenerator r]
generators =
  [outputView defaultExperimentSpecsView,
   outputView defaultInfoView,
   outputView $ defaultFinalStatsView {
     finalStatsTitle = "Statistics",
     finalStatsSeries = upTime <> immedProp },
    outputView $ defaultDeviationChartView {
     deviationChartTitle = "Deviation Chart",
     deviationChartLeftYSeries = upTime,
     deviationChartRightYSeries = immedProp },
   outputView $ defaultFinalHistogramView {
     finalHistogramPlotTitle  = "Histogram",
     finalHistogramSeries = upTime <> immedProp }]

We launch 1000 simulation runs. We want to see the statistics for variable values in the final time point, plot the deviation chart and plot histograms for the values in the final time point.

The launching file MainUsingDiagrams.hs looks different.

import Simulation.Aivika.Experiment
import Simulation.Aivika.Experiment.Chart.Backend.Diagrams

import Graphics.Rendering.Chart.Backend.Diagrams

import Model
import Experiment

main = runExperiment experiment generators (WebPageRenderer $ DiagramsRenderer SVG loadCommonFonts) model

You can find all source files in the following Zip archive.

Now the simulation experiment took 55 seconds on my MacBook Pro, when interpreting. A binary executable compiled with option -O2 spent 40 seconds on this task.

It printed the following output:

$ runghc MainUsingDiagrams.hs 
Updating directory experiment
Generated file experiment/DeviationChart.svg
Generated file experiment/FinalHistogram.svg
Generated file experiment/index.html

Here you can open the corresponding HTML file with the results of the simulation experiment. It contains the deviation chart and histograms.

Charting on Linux and OS X

As it was mentioned above, the users of Linux and OS X can also use the Cairo-based charting back-end. Probably, this is a preferred option for them. It works faster. It creates the PNG files of less size than the SVG files usually have. PNG is a raster image format, while SVG is a vector graphics format. Therefore, we have such a difference in performance, but the both formats are useful. Moreover, the creation of SVG files is made in Haskell more portable and cross-platform.

The prerequisites for using the Cairo-based charting back-end are as follows.

$ cabal install gtk2hs-buildtools
$ cabal install cairo
$ cabal install Chart Chart-cairo
$ cabal install aivika aivika-experiment aivika-experiment-chart aivika-experiment-cairo

Now we use another launcher MainUsingCairo.hs that looks different too.

import Simulation.Aivika.Experiment
import Simulation.Aivika.Experiment.Chart.Backend.Cairo

import Graphics.Rendering.Chart.Backend.Cairo

import Model
import Experiment

main = runExperiment experiment generators (WebPageRenderer $ CairoRenderer PNG) model

All necessary sources are contained in the same Zip archive.

Now the simulation experiment took 22 seconds on my MacBook Pro, when interpreting. A binary executable compiled with option -O2 spent 10 seconds on this task. This is 1000 simulation runs in a single-threaded mode.

$ runghc MainUsingCairo.hs 
Updating directory experiment(2)
Generated file experiment(2)/DeviationChart.png
Generated file experiment(2)/FinalHistogram.png
Generated file experiment(2)/index.html

As before, it saved the charts and created a new HTML file that you can open by the following link.

Here is the deviation chart that shows the trends and confidence intervals in the integration time points.

Here is the histogram that shows the distribution of values in the final time point.