Determistic simplifications of stochastic systems in SD

[tahonkan@cc.hut.fi]

Dear List readers, I am a post-graduate student at Helsinkin Universty of Technology. My doctoral dissertation will be reviewed and accepted in February. In the disseratation I present two SD models: one for industrial maintenance and repair, and one for spare parts supply chain. From the initial comments I know that there will be questions about the validity of the models. The maintenance model includes component wear and degradiation model (known as hazard curve) as a function of time. The aim is to model the behavior of the maintenance system and component reliability. Usually, the maintenance simulation and degradiation has been handled with stochastic discrete event monte-carlo simulation. Also, the maintenance and repair times (delays) are usually modelled with stochastic distributions (e.g. gamma, normal). The benefit of this approach is that the distributions have distributions-specific deviations and the simulation really generates these values randomly. Now the question is, can I handle these events (failures & maintenance) with their average values (e.g. failure_rate(t)= hazard_rate(t) and maintenance_delay=4 weeks)? And, are there any relation between 3rd order delay and probabilistic distributions? If anyone has seen this issue been analyzed I would appreciate any links, references and comments. Br, Tuomo Honkanen Helsinki University of Technology tuomo.honkanen@hut.fi

[=?iso-8859-1?Q?Jean-Jacques_Laub]

Hi Tuomo Your problem has too something to do with the general discussion about continuous, discrete time or event. Suppose that every period of time, say a day, there is an average of 100 people coming to make a trip. The capacity of the mean of communication is 100. A continuous deterministic modeler will draw a simple model taking the average of 100 as a rate , to stay simple. The result of the model, is no queue, and the mean of communication is always full. If the model calculates the average use of the capacity of the mean of communication it will be 100%. The reality will not be that depending on the distribution law of the people coming each day. This is obviously wrong. It is necessary to change to a stochastic continuous model. Another example: The same example. There is an average of 1 people coming every week, and the capacity is 1. The period is still one day, as the trip takes one day. If you are a continuous deterministic modeler, you will take 1/7 person every day, and it does not even make sense. Going to a continuous stochastic model does not make sense either, because the data need to be integers. You then need to go to a continuous time discrete event stochastic model. The second problem is evident, because you can see immediately that the problem needs discrete events. With the first problem because it is so simple, you can quickly find that you need some stochastic way of generating samples. But if your problem is more complicated, choosing to work with averages like is done with rates in most continuous SD model, can lead to wrong solutions. Do I have illustrated your problem correctly, if not do not mind. About the way to handle the events and representing them with their averages, it all depends on the problem. The average of the result is not in general equal to the result of the average. F(average(x)) <> average (F(X)) if F represents the model and x your event but it is often nearly equal but it can be very different. I have the same kind of problem. It occurs with stochastic simulation. Is it possible to replace some samples around an average by the average? And when rates for example can be replaced by their average and when not. It is often a question of common sense. When in doubt the solution is to calculate the average of the results it is much more safer, even if it costs a lot of calculation. Another problem of stochastic continuous time simulation with discrete or not events, is when you want to make some optimization. You generally get absurd results corresponding to exotic samples. The solution is to find values for the parameters that optimize the average of the results of a sufficient number of simulations. Using Vensim, I subscript all the variables of the model, or add one more subscript to those already subscripted, choosing for the size of the subscripts the number of simulation I choose. It is important to generate samples of events that are independent for each value of the subscript or each simulation. When you make a simulation the model calculates the values of the variables for every value of the subscript and you can then calculate the average of the results you want. The average will of course need no subscript. It is also easy to use the optimization possibilities of the software. You just choose as pay off the average of the results chosen. A drawback is the number of calculations. If the model is big that can be long. The model has to calculate all the parallel simulations and calculate the average and that is only one simulation. It all depends on the size of the model the number of simulations and on the way the average converges toward the optimum. Another problem when you use SD with stochastic simulation, is that discrete events necessitate generally integer values. Vensim handles integer functions, but does not recognize integer values like would any programming language where you can define a variable as integer. If you have a parameter which is an integer, you are obliged to define it as real, and to work with the integer part of it. You then get multiple solutions of your optimization because different values of the parameter give the same result. Integer(1.2) = integer(1.5) This is not a real problem, but it can eventually increase the number of calculations. On the other side, integer optimization takes too a lot of calculations, and optimizing is easier with continuous values. I think that in Powersim you can define integer variables. But I think it has no optimizer, perhaps because of that. I do not know how works I Think. In a recent message, Bob Eberlein, said that the next version of Vensim would do stochastic optimization. Maybe this will address these kinds of problem. Hope that helps a bit. I have used SD only for 18 months so I am no expert of the field, I only try to do what I can with what I know. So dont take what I say for granted. You could regret it. Experts know much better then me. Regards. J.J. Laublé From: =?iso-8859-1?Q?Jean-Jacques_Laubl=E9?= <JEAN-JACQUES.LAUBLE@WANADOO.FR>

["Nathaniel Osgood"]

Hi Tuomo, A few comments: The issue you have raised is a common one, but doesnt receive as much discussion as it deserves. In general, the key question when substituting mean values for a random variable (or mean realizations for variable driven by a stochastic process) in an SD model is whether the system (and the model that represents it) exhibits nonlinearity with respect to the value. If it does, using the expected value (E[.]) rather than drawing Monte Carlo samples from the distribution will give misleading results. The reason for this is that for a nonlinear function f and random variable X, E[f(X)] can be quite far from f(E[X]). In your case, we are dealing with the issue of a system that exhibits (presumably -- see below) stochastic transitions. As such we are dealing with the difference in model output between a stochastic process and a deterministic process. System evolution will in general be nonlinear with respect to transition parameters. How big a difference this is will depend on several things, notably the presence of positive feedbacks in the broader system that are sensitive directly or indirectly to fluctuations in the transition rates and details of your hazard function. I have done experiments some systems (comparing SD and agent-based models of infectious disease spread) in which the system (or at least the models that represent it) is HIGHLY sensitive to deviations. Without knowing further details about your system, it is hard to assess the level of difference to be expected between the mean of the stochastic process and the deterministic process you have described. Nonetheless, a few basic calculations suggest that the variations of the stochastic process from the deterministic trajectory may be significant. Lacking details on hazard(t), consider the simple case of a stochastic process associated with a hazard function that is uniform (value = lambda) over time. Within this process, failures observe a Poisson process, with the number of failures in a given timestep of length dt distributed according to a Bernoulli distribution (risk of per-component failure lambda * dt, # of trials equal to # of operating components) and with the likelihood of survival of a given component over time exponentially distributed (likelihood of survival of a given component after an interval of time t = Exp[-lambda * t]). If a group of components of size n start operating at time 0, the number of components remaining in operation after time t would be a random variable with mean Exp[-lamda*dt]*n, and the VARIANCE in this quantity would be n*Exp[-lambda*t]*(1-Exp[-lamda*t]). As t rises, this variance eventually reaches a value as large as n*.25. This suggests that the point deviations of the stochastic process from the deterministic process could be rather high. >From your description, it was not clear whether nonlinearity that would boost the effect of deviations in transitions would be present elsewhere in the system, but it seems plausible that it would be. For example, in terms of the dynamics, the length of time required before a given item awaiting repair is fixed would seem likely to depend on the number of items requiring repair -- but how significant this effect is will depend on the specifics of the system. Nonlinearities in the dynamics of a system such as this can lead to magnifications of stochastic fluctuations. You may also have nonlinearity in relationships in other areas of the system -- for example, a nonlinear rise in the cost over the span of a components "downtime". Using average transition times rather than representing the distribution may lead to misleading results here as well. Thinking more broadly, you consider whether the the seemingly stochastic variations failures and repair times simulated in the traditional discrete event models you mrefernece might reflect in part a underlying heterogeneity in the component population. If so, this could significantly challenge the results of those models. In particular, if due to heterogeneity a given component is likely to be associated with a persistent repair profile over a period of time, both dealing with the average trajectory AND treating the failures and repair times as stochastic are likely to miss an important component of the system, and may give misleading results in simulations. Finally, from your description, I wasnt sure I was interpreting the form of the hazard function properly. In your description, you mention that the hazard rate is a "function of time" -- which would imply a non-stationary process (with the hazard uniform across all working components). Perhaps what you meant was that the hazard rate specific to a component changes as a function of the time that particular component has been in circulation? If the latter is correct, this raises the question of how you are keeping track of the time a given component has been in circulation -- through stocks disaggregated (or arrayed) by time, or another mechanism. In the presence of a nonlinear hazard, you will want to take some care in representing this within the model, but given the timing this probably isnt a first-order consideration. Thanks for bringing up this interesting topic. I hope these comments are helpful. Best, n ----------------------------------------------- Nathaniel Osgood, PhD Research Associate (TDP), Senior Lecturer (CEE) Massachusetts Institute of Technology 77 Massachusetts Avenue Rm 1-175 Cambridge, MA 02139 (617) 253-9725 nosgood@mit.edu