[R-sig-dyn-mod] Discrete event models
Thomas Petzoldt
Thomas.Petzoldt at TU-Dresden.de
Tue Sep 14 06:36:18 CEST 2010
Dear Jose,
many thanks for the interesting examples and sorry for the delay because
of vacation and conferences. Here are some ideas:
1) Both, the activity-based and the event-based method can be
implemented in simecol using class indbasedModel (or a derived class).
For the activity-based approach you may proceed like follows:
- define init as a data structure, e.g. a list that contains the state
of the model, the event queue and the clock.
- times can be a vector with fixed time steps as usual that indicate the
iteration counter, while the real simulation time (the events) are
stored in the init list. The event queue can be a data frame, e.g. like
follows:
data.frame(clock= ..., eventType= ..., furtherArgs = ...)
- vectorizing is possible. So multiple events at the same time (or
within a discrete time window) can efficiently be implemented by
applying subset() on the data frame with the events.
2) A "process-based" approach in your sense would require the definition
of a kind of meta-language. While this is of course possible in R, I
currently don't know someone who has it on his/her task list.
If there is someone else who is interested to do this, it would be fine,
but there may be the risk that this would be rather specific for a
particular model paradigm, because discrete event models form a very
diverse field.
3) There is also the possibility of a hybrid approach so that the
simulation is controlled in R, but execution takes place in another
language. There are several possibilities to communicate between R and
Python and there are also possibilities to communicate with other
event-based environments, e.g.:
Thiele, J.C. and Grimm, V. (2010) NetLogo meets R: Linking agent-based
models with a toolbox for their analysis. Environmental Modelling &
Software 25(8), 972-974. http://dx.doi.org/10.1016/j.envsoft.2010.02.008
4) Depending on your particular question you may also have a look on
other Markov-type models, e.g.:
* Package "pomp" for simulating and analysing partially observed Markov
processes.
* GillespieSSA: Gillespie's Stochastic Simulation Algorithm (SSA), see
also http://www.jstatsoft.org/v25/i12/paper
* and some more on CRAN.
Thomas Petzoldt
On 9/5/2010 5:34 PM, jose romero wrote:
> Hello list:
>
> On the previous topic I asked about using simecol for discrete event,
> process-based simulation. Included bellow is the code for the
> simulation of a simple M:M:1 - GD:oo:oo queue using three approaches
> (activity based, event based and process based). The first two
> scripts are R code, the last is Python using the SimPy library.
>
> Please note the following: 1) The activity-based version is the
> slowest of all. Furthermore, as the ocurrence of arrivals and
> services (both Poisson process) is determined by Bernoulli trials,
> this approach is very inexact. No doubt this simulation model could
> be implemented by the individual based simecol objects (more
> time-efficiently?)
>
> 2) The event based version is more efficient, as the clock jumps from
> event to event. However, the programmer has to take care of all the
> book-keeping details (event scheduling). Consider for example how
> one would have to re-write the code completely if more servers are
> used, there is queue reneging behavior, etc.
>
> 3) the process based approach is the simplest for the programmer. It
> would be nice if something like SimPy could be implemented in R. It
> occurs to me that the fork R-package, which handles processes, could
> be used for this end.
>
> My questions are:
>
> Is this process based approach way too far out of the simecol scope
> (simecol objects being essentially based on time iterators, solvers
> and such)? Is the process-based approach useful for the
> implementation of other ecological models? Or considering the
> event-based approach, could such programs be vectorized? Is there a
> way to implement solvers in simecol that advance the clock to the
> nearest future event?
>
> Thanks and Regards,
>
> jlrp
>
> ##=====================================================================
> ## Simulation of a simple M:M:1 - GD:oo:oo queue
> ## Activity based approach
> ## José Romero - sep 4, 2010
> ##=====================================================================
>
> #---------------------
> # initial parameters
> #---------------------
>
> tmax<- 10000 v #maximum time of simulation
> alpha<- 6 #arrival rate
> lambda<- 8 #service rate
> queue<- NULL #client queue, initially empty
> service_times<- NULL #record of service times
> times<- numeric(0) #record of total time for each state
> step<- 0.001 #time-steps
> p_alpha<- alpha*step #probability of arrival in one time step
> p_lambda<- lambda*step #probability of departure in one time step
> clock<- 0
>
>
> #---------------------
> # function definitions
> #---------------------
>
> not_empty<- function(queue) return(!(length(queue)==0))
>
> empty<- function(queue) return (length(queue)==0)
>
> Update_clock<- function() {
> state<- length(queue) + 1
> if (is.na(times[state])) times[state]<<- 0
> (sample(size=1,c(0,1),prob=c(1-p_alpha,p_alpha))==1)
> times[state]<<- times[state] + step
> clock<<- clock + step
> }
>
> #----------------------------------------------------------------------
> # Main simulation loop:
> #----------------------------------------------------------------------
>
> for (i in 1:(tmax/step)) {
> arrival<- (sample(size=1,c(0,1),prob=c(1-p_alpha,p_alpha))==1)
> if (arrival) queue<- c(queue,clock)
> dispatch<- (sample(size=1,c(0,1),prob=c(1-p_lambda,p_lambda))==1)&
> (not_empty(queue))
> if (dispatch) {
> arrival_time<- queue[1]
> queue<- queue[-1]
> service_times<- c(service_times,clock - arrival_time)
> }
> Update_clock()
> }
>
>
> #----------------------------------------------------------------------
> # Simulation results:
> #----------------------------------------------------------------------
>
> #first convert all NA's in times to 0
> times[is.na(times)]<- 0
>
> #then divide times by clock to obtain probability distribution
> times<- times / clock
>
> #last compare to theoretical state distribution
> #and plot results
> max_x<- length(times)-1
> theoretical<- dgeom(x=0:max_x,prob=1- alpha/lambda)
>
plot(x=rep(0:max_x,2),y=c(times,theoretical),col=c(rep("blue",max_x+1),
> rep("red",max_x+1)))
> mean(service_times)
>
>
> ##=====================================================================
> ## Simulation of a simple M:M:1 - GD:oo:oo queue
> ## Event based approach
> ## José Romero - sep 4, 2010
> ##=====================================================================
>
> #---------------------
> # initial parameters
> #---------------------
>
> tmax<- 10000 #maximum time of simulation
> alpha<- 6 #arrival rate
> lambda<- 8 #service rate
> queue<- NULL #client queue, initially empty
> service_times<- NULL #record of service times
> times<- numeric(0) #record of total time for each state
> clock<- 0
>
>
> #---------------------
> # function definitions
> #---------------------
>
> not_empty<- function(queue) return(!(length(queue)==0))
>
> empty<- function(queue) return (length(queue)==0)
>
> Update_clock<- function(next_event_time) {
> time_dif<- next_event_time - clock
> state<- length(queue) + 1
> if (is.na(times[state])) times[state]<<- 0
> times[state]<<- times[state] + time_dif
> clock<<- next_event_time
> }
>
> Arrival<- function(queue) {
> queue<<- c(queue, clock)
> arrival_time<- clock + rexp(n=1,rate=alpha)
> return(arrival_time)
> }
>
> Dispatch<- function(queue) {
> arrival_time<- queue[1]
> queue<<- queue[-1]
> return(arrival_time)
> }
>
> #----------------------------------------------------------------------
> # Main simulation loop:
> #----------------------------------------------------------------------
>
> next_arrival<- rexp(1,rate=alpha)
>
> while (clock<tmax) {
>
> if (empty(queue)) next_dispatch<- next_arrival+rexp(n=1,rate=lambda)
>
> while (next_arrival<=next_dispatch) {
> Update_clock(next_arrival)
> next_arrival<- Arrival(queue)
> }
> while (not_empty(queue)& (next_dispatch<next_arrival)) {
> Update_clock(next_dispatch)
> service_time<- clock - Dispatch(queue)
> service_times<- c(service_times,service_time)
> if (not_empty(queue)) next_dispatch<- clock +
rexp(n=1,rate=lambda)
> }
> }
>
> #----------------------------------------------------------------------
> # Simulation results:
> #----------------------------------------------------------------------
>
> #first convert all NA's in times to 0
> times[is.na(times)]<- 0
>
> #then divide times by clock to obtain probability distribution
> times<- times / clock
>
> #last compare to theoretical state distribution
> #and plot results
> max_x<- length(times)-1
> theoretical<- dgeom(x=0:max_x,prob=1- alpha/lambda)
>
plot(x=rep(0:max_x,2),y=c(times,theoretical),col=c(rep("blue",max_x+1),
> rep("red",max_x+1)))
> mean(service_times)
>
>
> ##=====================================================================
> ## Simulation of a simple M:M:1 - GD:oo:oo queue
> ## Process based approach
> ## Jose Romero - sep 4, 2010
> ##=====================================================================
>
> #---------------------
> # initial parameters
> #---------------------
>
> from SimPy.Simulation import *
> from random import *
>
>
> class Client(Process):
> served_clients = []
> def go(self,arrival_time,myServer):
> yield request,self,myServer
> yield hold,self,expovariate(lambd=8)
> yield release,self,myServer
>
> class Client_source(Process):
> def go(self,finish,myServer):
> while(now()<finish):
> c=Client()
> activate(c,c.go(arrival_time=now(),myServer=myServer))
> yield hold,self,expovariate(lambd=6)
>
> initialize()
> endtime=1000
> server=Resource(capacity=1,monitored=True,monitorType=Monitor)
> client_source=Client_source()
> activate(client_source,client_source.go(finish=endtime,myServer=server))
> simulate(until=endtime)
> print server.waitMon
>
>
>
>
> [[alternative HTML version deleted]]
>
>
>
>
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
More information about the R-sig-dynamic-models
mailing list