Models for Event Based Stochastic Processes

If we generalize a finite state automaton by letting the time between two events be stochastic, and possibly also letting the choice of the next state be stochastic, we obtain a model for event based stochastic processes. Such models have many applications, of which one is queuing systems.

1 The Random Package

package Random
import Modelica.Math;
constant Real NV_MAGICCONST=4*exp(-0.5)/sqrt(2.0);
type Seed = Integer[3];

function random "input random number generator with external storage of the seed"
  input Seed si "input random seed";
  output Real x "uniform random variate between 0 and 1";
  output Seed so "output random seed";
algorithm
  so[1] := rem((171 * si[1]),30269);
  so[2] := rem((172 * si[2]),30307);
  so[3] := rem((170 * si[3]),30323);
  // zero is a poor Seed, therfore substitute 1;
  if so[1] == 0 then
    so[1] := 1;
  end if;
  if so[2] == 0 then
    so[2] := 1;
  end if;
  if so[3] == 0 then
    so[3] := 1;
  end if;
  x := rem((so[1]/30269.0 +so[2]/30307.0 + so[3]/30323.0),1.0);
end random;

function normalvariate "normally distributed random variable"
  input Real mu "mean value";
  input Real sigma "standard deviation";
  input Seed si "input random seed";
  output Real x "gaussian random variate";
  output Seed so "output random seed";
protected
  Seed s1, s2;
  Real z, zz, u1, u2;
algorithm
  s1 := si;
  u2 := 1;
  while true loop
    (u1,s2) := random(s1);
    (u2,s1) := random(s2);
    z := NV_MAGICCONST*(u1-0.5)/u2;
    zz := z*z/4.0;
    if zz <= (- Math.log(u2)) then
      break;
    end if;
  end while;
  x := mu + z*sigma;
  so := s1;
end normalvariate;

connector discreteConnector
   discrete Boolean dcon;
end discreteConnector;

end Random;

2 Examples

As a simple example of discrete-time stochastic process simulation we can consider a simple server model with a queue.

2.1 CustomerGeneration Model

The CustomerGenerator model contains an output connector outCustomer that transmits the customer arrival events to the rest of the system. We first start by defining a model which generates customers at random points in time. The time delay until the next customer arrival is assumed to be a normally distributed stochastic variable T, here computed by calling the random number function normalvariate, see above, where the negative part is folded back on the positive side using abs(T) to avoid unrealistic negative service times. Since this function is side-effect free as all Modelica functions, it returns a randomSeed value to be used at the next call to normalvariate, see the help functions above.

model CustomerGeneration
   import Random; // See Math.Random package in Appendix D
   EventPortOutput outCustomer;
   BoolView bv; // Only for plotting to see events
   parameter Real mean = 0;
   parameter Real stDeviation = 1;
   discrete Real Tdiff; // Its absolute value is timediff to next customer
   discrete Real nextCustomerArrivalTime(start=0);
   discrete Random.Seed randomSeed(start={23,87,187});
   Boolean changeflg(start=false);
equation
   when pre(nextCustomerArrivalTime)<time then // Customer arrives
    (Tdiff,randomSeed)=Random.normalvariate(mean,stDeviation,pre(randomSeed));
    nextCustomerArrivalTime = pre(nextCustomerArrivalTime) + abs(Tdiff);
    changeflg = not pre(changeflg); // Changes at each customer arrival
   end when;
   outCustomer = change(changeflg); // Signal true only at event
   bv.inv = outCustomer; // Only for plotting, widening event peaks
// real value comparison to zero not reliable for short durations at events:
// outCustomer = edge(nextCustomerArrivalTime >
// pre(nextCustomerArrivalTime) ); -- wrong, wsignal always false
end CustomerGeneration;

2.1.1 Simulation of CustomerGeneration

We simulate the CustomerGeneration model for a while to observe the customer arrival times.

2.2 ServerWithQueue Model

We can now design a server model ServerWithQueue which through its inCustomer port will accept the arriving customers from the CustomerGeneration model. The nrOfCustomersInQueue represents the current number of customers in the queue. Other state variables of interest are nrOfArrivedCustomers, the number of customers that have so far arrived and nrOfServedCustomers, the number of customers that have so far been served. The first when-equation services the arrival of a new customer, whereas the second when-equation handles the event when a customer is serviced and is ready to leave the queue.

The event when readyCustomer becomes true represents completed service for the current customer. The time needed to serve a customer is custServeTime, which is constant in this model but could be varying according to some formula or be randomly distributed. The variable lastServeStart represents the starting time for servicing the most recent serviced customer.

model ServerWithQueue
   EventPortInput inCustomer;
   parameter Real custServeTime=0.1 "Time needed to service a customer";
   discrete Real nrOfArrivedCustomers;
   discrete Real nrOfServedCustomers;
   discrete Real nrOfCustomersInQueue(start=0);
   discrete Boolean readyCustomer(start=false);
   discrete Boolean serveCustomer(start=false);
   Real lastServeStart "Most recent start servicing";
equation
   readyCustomer = serveCustomer and (custServeTime < time-lastServeStart);
   nrOfCustomersInQueue = nrOfArrivedCustomers - nrOfServedCustomers;
   when inCustomer then // Event: New customer arrives
      nrOfArrivedCustomers = pre(nrOfArrivedCustomers)+1;
   end when;
   when readyCustomer then // Event: Customer finished service
      nrOfServedCustomers = pre(nrOfServedCustomers)+1;
   end when;
   when (pre(nrOfCustomersInQueue)==0 and inCustomer) or
     (pre(nrOfCustomersInQueue)>=1 and pre(readyCustomer)) then
     serveCustomer = true;
     lastServeStart = time;
   end when;
end ServerWithQueue;

2.2.1 Test Model 1

The simulation of the whole system can be done by defining a test model which connects the CustomerGeneration model with the ServerWithQueue model.

2.3 Simulation of testServer1

Using the MathModelica simulation environment we simulate the model for 10 seconds with the customer service time custServeTime = 0.4 seconds.

Number of arrived customers, number of served customers, and number of customers waiting in the queue or being serviced. Customer service time is 0.4 seconds.

2.3.1 Test Model 2

2.4 Simulation of testServer2

We now simulate the model for 10 seconds with the customer service time custServeTime = 1 second.

Number of arrived customers, number of served customers, and number of customers waiting in the queue or being serviced. Customer service time is 1 second.