Queueing Theory in a World where most Queueing Problems are Solved by Simulation (Keynote)
CC Attribution - NonCommercial - NoDerivatives 4.0 International:
You are free to use, copy, distribute and transmit the work or content in unchanged form for any legal and non-commercial purpose as long as the work is attributed to the author in the manner specified by the author or licensor.
Grassmann, Winfried
Monte Carlo simulation is one of the most successful techniques, not only in operations research and performance evaluation, but in science in general. One reason for this extraordinary success is its flexibility. In contrast, most queueing models are rather specialized. In this talk, we suggest methods to make queueing theory more flexible. In particular, we suggest an event-based approach, which provides great flexibility for the modeller. We also show how to convert such event- based models into Markov chains, which can then be solved by classical numerical methods. The suggested method is particularly suited for small models, where its execution times are much lower than Monte-Carlo simulation. For larger problems, the curse of dimensionality takes over, and the execution times based on classical numerical methods increase exponentially. This means that for complex models, simulation finds numerical solutions with less computer time than classical numerical methods.
 And here you see again the title and my name as you see is Winfield Grasmund and I talk about the Lyme Monte Carlo simulation is done so frequently. Monte Carlo simulation is really one of the most successful applications of operations in the field search and even beyond. However the bulk of Turing theory uses deterministic methods that is methods not suggesting to randomness with random variables and outside Turing theory many Turing models are solved by Monte Carlo simulation and the results of Turing theory are often ignored. And how can we make change that? How can we make classical Turing theory more successful? Now the first question is why is simulation so successful? And I see three reasons for this. You see maybe others or more. One for large problems Monte Carlo simulation needs far less computer time. Two simulation is much much easier. Most people have no difficulty doing simulation and three simulation is much more flexible. Now how can Turing theory come to this? Let us look at the different issues, the first issues the execution times. In Monte Carlo simulation execution times start at a high level but they increase linearly with more complexity as measured say by the number of variables you use. In deterministic methods execution times start at a low level but they increase exponentially with more complexity. So Turing theory is an advantage for small models but not for large models. Now to show you this I have here this graph. The red line is the deterministic methods. They start at a low level and go up very steeply and very soon they take more computer time than simulation, the blue line. Our simulation needs a lot of computer time initially and in fact before the advent of computer simulation was not a practical method. However the computer time increases very slowly with a number of values. Now the next thing is mathematics. How can we simplify mathematics? I think we should express numeric methods. Maybe that is a little bit against the stream to say that people want results rather than formulas. Now if you use numeric methods you have a choice. You can use the standard methods or you can use methods tailored to deal with probabilities. Now for standard methods you have lots of software available including software included in simple manipulation languages. However in some cases it is advantageous to exploit the special properties of probabilities mainly that probabilities are also always greater than or equal to zero. And this enables us to avoid subtraction and just reduce rounding errors. Now I should be clear here, subtractions do not cause rounding errors but they increase rounding errors committed earlier, cases quite rapidly often by orders of magnitude. Now if you use special algorithms to deal with probabilities they often have a probabilistic interpretation and that provides additional insights. Now the third thing is can we increase the flexibility? One of the most used simulation paradigms is discrete event simulation. Can we copy this paradigm? Now we can use discrete event systems as they were described by Cassandra and Laportzli in 2008 in a book. And what are discrete event systems? Well the discrete event systems has two main components, the state and events. Now the state is a set of fibers, for instance, q-legs. And all changes occur in discrete steps and discrete event simulation. Now every change is called by an event, because an event can even be defined as something that changes the state. Now q's are clearly discrete event systems, why not use this paradigm during q-legs? Now first let us look at the state. First the state consists of physical state fibers, fibers that you are interested in, it could be q's because we assume that the state fibers are x1, x2, up to xd. And we also assume that the xi are between 0 and ni. In addition to this we need supplementary state fibers, if you have a schedule of future events you need to keep track of the schedule and for this you use a supplementary state fiber or you may record when past events have happened, again that would be supplementary fibers. What you need is enough supplementary values such that the distribution of the future state fibers can be calculated. And you do that and the system becomes a Markov chain consisting of physical and supplementary state fibers. Next we come to the events. Events come in different types such as arrivals, departures, switches of rating lines and so forth. And each event is characterized by three things, an event function, an event condition and the process that generates the events. The first the event function we call it fkx and if the state before event k is x, the state after the event k is fkx. For instance if your state is x1, x2 and you have an event function fkx1, x2 which increases x1 by 1 leaving x2 unchanged and this would be the event function of something like an arrival to line 1. Next we have an event condition and events are prevented from occurring if the condition is not met. And then we have an event process that could be a process, renewal process, a phase start process and hold action. Now here we have an example of a discrete event systems. There are two bins, bin 1 and bin 2. They have finite size, bin 1 the size n1, bin 2 the size n2 and arrivals to bin 1 are lambda 1, arrivals to bin 2 are lambda 2. The bins are used to do repairs and each repair always needs one unit from bin 1 and two units from bin 2 repairs occur at a late new and all these events are possible. Now this system has two state parameters x1, the number in bin 1 and x2, the number in bin 2. And what are the events? We have three events, we finish one, we have two and repair and each event has an event function, an event condition and a rate. So we finish one, we have a path, we see u increase x1 by 1 leaving x2 unchanged and this has the event condition x1 is less than n1 because if the thing is full you don't accept arrivals anymore and the rate is lambda 1. Similarly you have event number 2, it simulates as a rate of lambda 2 and repair has the event function x1 minus 1, x2 minus 2 and this event can only happen if you have enough in the bins, that is x1 must be greater than 0 and x2 must be greater than 1 and if that is 2 then the event happens at a late new. So now you have formulated this event and now you need a transition matrix because the whole thing is a Markov thing. So we take the transition matrix Q and I do it here for any Quasso events because we have to find Qij, the rates of going from state i to state j except for the diagonal which you all know must be determined such that the sum across the row must be 0. Because you take a state out of the state space and you apply all events and for event k you do the following, if the event condition holds, if Ck holds and you have a new state after the event you call it x new which is xk of the present state and the transition again is x, x new which is set to the rate of event k and this you do for all events and for each state. Now this is not really a transition matrix because the x are vectors rather than numbers because we have to associate our vectors with numbers because we have to find for each state a number and this number is used to find both the row and the column of the state and this number must be in a one to one relation with the state. Now if the state space is as simple as this which is true for many events including the for many discrete event systems including the event systems I talked about then you can use the alphanumeric order. This that is the event number is the rank in the alphabetic order the state has and this event number can be found by the following formula. You take the sum of uixi where the ui are the products ntray plus one or tray equals i plus one to d. Okay so now we have events event number and in some cases you can trust enumerate the event numbers and you get a matrix gather quickly. Now things are more complicated if events are generated by renewal process then we must add supplementary variables such as the time since the previous event of the same type or the time to the next event if the event is scheduled. Now the times between events are often continuous random variables and in this case we need to discretize unless we use phase type variables. Now discretization has a disadvantage it forces us to deal with multiple events that occur simultaneously and that increases the number of entries in the transition matrix by considerably. Now if you look at equilibrium solutions we can save computer time by embedding the system at the points where the events do occur. Now once you have a generated your transition matrix you want to find the transcendent equilibrium probabilities that is we want to find the ijt the probabilities to be in state j at country and you probably all know how this is done. We have a recursive formula here that ijt is ijt plus one is a function of iitpij. Now as you all know if t is large enough then the pijt corresponds to an equilibrium probability pij and in order to find the equilibrium probability we have a system of situations which is obtained from this system by just popping the argument t. Now for continuous time of trains we have a similar relationship we have a rate matrix or the increment generator q. Now instead of a recurrence relation we have a dependency equation because pij prime t is the derivative of pij with respect to t and here you have a formula for pij prime t. Again it is large enough that converges to an equilibrium probability which is written down here. And now the next step is to solve these equations both as transient equations and equilibrium equations numerically because this is what I propose and let us do that. Now first however let us look at the properties of the transition matrix. The matrices are huge but sparse. If all events are Poisson processes the matrix size is given by the ranges of the different state drivers and the product of the range is of the different state drivers. Now the density however is low because the entries per row in the transition matrix is equal to the number of events and let the number of events be e. Consequently the density is e divided by n. Now if the events are renewal processes you also have to add the products of the ranges of your supplementary drivers which I call yi and then you get the matrix size of this huge product. Now in each row we actually have two to the e events. It is two to the e because we have to look at all combinations of events that can have an entity. This is an approximation but it is good for m I know too small. Because we have huge matrices. The matrices are also banded. In fact if x1 changes at most by 1 the state number and that represents the band and it changes by u1 which is the product of the ranges of the different state drivers except for x1. It is n1 plus 1. Now here I have an example. I have four state drivers and the ranges of the state drivers are from 0 to 9 and if you have all the elements plus or that gives 10 to the fourth of 10,000 states. However the transition matrix has 10,000 to the square entries which is 100 million. And if you assume five events as it's true for the tandem queue you have only five times 10,000 that is 15,000 non-zero entries. Now the same is true for a lesser to a lesser extent also for men's generated value of u reposes. This suggests that we store only the non-zero entries of the transition matrix. And I apply this idea now in order to find transient solutions of the discrete time mark of time. Thus we use three arrays, one for the row, one for the column and one for the probabilities. Thus we want to apply now this method to calculate pjt plus 1 according to this formula. This is a formula for transient probabilities in this discrete time mark of time. We assume that the total number of non-zero entries we call the dn-top which in some cases in most cases is 10 times the number of events. Thus this is the number of non-zero entries and we want to do the calculation for p going from 1 to t top and here you have the problem to do that. We go from t for 1 to t top for n from 1 to n top. We go to all the non-zero events. And this formula really represents this formula here. The pi nu tray is ijt plus 1, the pi ti is iit and this is the probability ij. Now i is of course what you find in row i, tray is what you find in row n, tray is what you find in column n and the probability is the next. Just to do that for all n. And the remaining statement should be clear because I don't want to go into it. But the important thing is here is that the complexity is proportional to n. This is the important thing. Now we can look at Planck's solution in continuous time. In this case you could solve the differential equations but we suggest another method called randomization or uniformization. It is pretty standard now. In this method you create a matrix p which is q divided by f plus the identity matrix. Where f is chosen such that p is a stochastic matrix. And this we can look at the Markov chain we have described by p. And after n plus 1 steps to the probability to be in state tray is this formula here. But there is now a difference. The steps are not steps of one. They are steps of exponential random bias with rate of f and rate of f. And consequently the pjg the probability to be in the continuous Markov chain in state tray at p is given by this formula which has a pj super scheme n multiplied by the force of distribution. And as I said this is the randomization or uniformization. This method can also be derived by a small modification of the Taylor series. So here we have an effect on mathematics. Now Taylor expansions of matrix exponentials tend to be numerically unstable. And it was argued that no that first hand that this method is potentially unstable. But this is not true and the reason is that we have no subtractions. And if there is no subtractions rounding error think is very slow. So randomization or uniformization as this method is called works for millions of states and it works for very high heat and out problems. And what is even better it performs better than standard algorithm for solving the effect. So far transient solution and now we come to steady state solutions. We discuss only discrete and Markov chains. First we have here the formula for the steady state of a discrete Markov chain. This is a system of equations and we can use Gaussian elimination. Now the mathematicians say Gaussian elimination is unstable if you have something like 10,000 or 100,000 states. But this is not necessarily true because we can reformulate the system such that subtractions are avoided. And we do that when the resulting algorithm has a probabilistic interpretation. Now if you look at linear algebra box then if you solve linear equations and linear equations then you start with the first equation and solve this first equation for x1. We don't do this. We start with the last equation and solve it for pi n. Well why do we do that? Well mathematicians assure us that routing errors are minimized if you start with the smallest n to be first. Just using this mathematic is a mathematical recipe, we eliminate pi n first and then you get the following result. You have Gaussian elimination. Here you have the system of linear equations. We set the tray equal to n, that is that it becomes pi n and we solve for pi n. That's the solution of the last equation for pi n. Then we do the substitutions, that is we first split the sum into the sum from pi equals 1 to n minus 1. Then we have here the last term which only includes pi n. We plug in this value here and we simplify and what you have is a new system of equations which has n minus 1 linear equations covered in n and the coefficients are now given by pi j plus pi n pi nj divided by 1 minus pi nn. Now you can continue. First we can call these new coefficients pi j n minus 1 which is pi j plus pi n pi nj divided by pi 1 minus pi nn. This gives us a system of equations which we could solve for pi n minus 1. Then we can get a new system of equations which is sorted for pi n minus 2 and so forth and so forth. The coefficients after we have solved all the pi's from n up to a capital N from small n up to capital N are given by the following formula here. Pi ij super skipped n minus 1 is pi ij n plus pi i n super skipped n pi nj super skipped n divided by 1 minus pi nn super skipped n. In order to obtain this formula we also need pi n which is this equation. It expresses pi n in terms of the pi i where the i is less than n. This can be used in the black substitution phase for calculating the pi n. Now I told you there is a probabilistic interpretation of the whole thing and to discuss it we first need to use embedded. Now what is embedded? Suppose you have a sequence x1, x2, x3 and so forth which would be a Markov chain and in addition we have a set c. And now we create a subsequence of this sequence here and elements xb of the original sequence are included in the subsequence only if they are elements of c. Because if the set c is 1, 2 and 3 and if xb is 1, 3, 5, 2, 4, 2 then what happens? Well 1 is in c so we keep it, 3 is in c so we keep it, 5 is not in c so we omit it, 2 is in c we keep it, 4 is not in c we omit it and 2 is in c. Now we can actually embed into embedded Markov chains. First we can use this new sequence we call it x superskip c and we can embed it into the set d where d is say 1 or 2 and so we get 1, 2, 2. Now if you look at this you will see that instead of embedding d into this sequence into the sequence x superskip c you could have embedded d also in the original sequence and we would have the same result. First what happens is if you have repeated embedding then you could as well do embedding in one step or vice versa in order to do embedding in for many, many different states you can break it down into one step embedding sets. Because embedding can be done in either in several small steps or equivalently in one large steps. Because we can reduce the state space by one state at a time until only one state is left and this I call state data. Now if you have a Markov chain and if you do embedding by set c then the corresponding embedded chain is also a Markov chain. The problem is now to find the transition probabilities pij superskip c and what is the solution? The solution is to add all the probabilities of all the paths that start in i ending j and avoid any state c in between. Now what happens is that elimination Gaussian elimination is really embedded. First we start with the original sequence states 1, 2, 3, we embed it into a set c which goes from 1, 2 up to n minus 1 and in order to find c we add the probabilities of all the possible sequences to go from i to j avoiding any state less than n. We can go from i to j either indirectly that's pij we can go from i to n and from n to j that's pij, pij, we can go from i to n, we can stay in n, go to j and the probability of this path is pij, pij, pij and so on. Now you go on, for this sequence here you find the probability pij, pij and n square pn, pij. Now we have to add all these probabilities together and we get these probabilities for the embedding set 1 to n minus 1 which is equal to pij, pij, pn, n to the k, you stay k times in state n, pnj and as you see this is a geometric series and the sum of the geometric series of these terms is 1 minus pnn, first what we get is pij plus pij divided by 1 minus pnn and this is nothing else but pij superskipped n minus 1 and the rest follows the same holds now for all pij superskipped n. To state it out you reduce the state space by 1 into 2 step and for the general term you get pij n minus 1 is pijn, pij and n to the n to the n to the n minus pnn and this is the transition matrix for the Markov chain embedded in the set of 1 to n minus 1. Now here is a subtraction and I said we should avoid subtractions and this can be done because all the pnj here are transition probabilities of Markov chains that is the sum across the j is equal to 1 or alternatively if you look at state n, 1 minus pnn superskipped n is the sum j equals 1 to n minus 1 pnj n such that the sum here is 1 if you move that over and it follows that this formula can be written without using subtractions and you get the 3th method the method of TACS-Altas-Mantak-Sahay. Now this method though has a problem it is impossible to exploit specific well it is possible but it is very difficult. Now if the matrix is tense the number of floating point operations goes up with the 3rd power of n. In the example I presented before we had n10 to the 4th state that is 10,000 because this means according to the formula here that we get 2 3rds 10 to the 12th flops and if you do that or if you do these many flops on a laptop as I have here as you have you have then that means around 100,000 seconds compute the time that is one day on the laptop. First this is not a good thing really. Now you may say well we cannot exploit bandedness well what do you get? You get as you can divide by the bandwidth which is denominator here but still you get an n to the 3rd in the numerator. First the algorithm are always n to the 3rd. First in the example you would reduce the time from one day to 100 so but they which is still you don't like to wait too long until you get this one. What is even worse is if you look at transient solutions they increase this n rather than with n to the 3rd. Those transient solutions in a way are easier to obtain than steady state solutions. Now this is kind of strange but for large-stance matrices this is good. In fact in order to find equilibrium solutions you can just iterate the transient solutions until you have reached the equilibrium. Now you can increase convergence a little bit by different tips such as using the trumpets and so forth. You can also go back to classical methods such as Gauss Seidel I use that extensively and the 3rd by 12. Now in Gauss Seidel the order of the states is important and as has been shown by Mitter and Zuckers if you order the states such that the events downstream changing state by the downstream are done first then the convergence is faster. First you look at the flow. Again you have to look at the meaning of the model and by exploiting the properties of the model you can improve your methods. Now state reduction is still great for some theoretical results and I have here two theoretical results that mention the 3rd one. First it allows you to deal with infinite state spaces. Secondly it allows you to deal with matrices that have a discrete cluster and third they can allow you to do good embedding techniques which I don't have the time to discuss. Now first infinite state spaces. Here there is the following result. If J is recurrent and the state space can be cut in such a way that the Pij superskips the N. The only thing that you need to calculate the Pij N does not change by more than epsilon as long as xd plus 1 minus xp the up comes about it. And here is the proof. If J is still recurrent most parts from i to j ending as soon as j is reached have a finite length. Because we now only look at most parts and most parts as we can say these are the parts less than u and if you look only at those parts less than u then the number of states reachable between two visits of J is finite and the states that cannot be reached on such a path can be ignored and we can cut them without changing Pij superskips N by more than epsilon. First in other words you can always cut the state space to a finite size if the Markov chain is recurrent and if the up comes about it. Now the next thing is let's look at transition matrices with repeating columns. Now here I have a very simple proof but you find more complex proofs in the literature. If Pij is equal to Pij j plus k for almost all ij which is true for all non methods then we have Pij N equals Pij j plus k j plus n plus k and the proof is simple. If you have any path from i to j superskips N you can match that with a path from i plus k j plus k n plus k by just shifting the path upwards by k units and when doing this you have you meet the same transition over the next. Consequently this is our next. Now this essentially finishes my talk and here I have some conclusions to go back to the beginning so to say. You make a previous conclusions you may easily what I think these things should be discussed. During this flexibility I propose this big event approaches to reduce computer time. You can try to reduce the number of state values and to make mathematics easier use your method possibly modified to do that for a previous case. Okay that is the all here we have number of references which I found useful and of course my own papers I will use at least for you. Okay any question? Pij Jank please. Yes I am here. Yes I am here. I am hearing you. So we are afraid it is very nice very interesting talk. So when we do numerical and also the Markov for Mendecavo what is your comments on a statistical approach in the future for queuing models? You see I mean what is my suggestion for the future of queuing theory is that the question? Nowadays say the advances of this telecommunications smartphones and all this real-time apps and a lot of times we do not have really have models and we have input we have something and then we have output but in order to do the analysis or to design the system in an optimum way of course sometimes based on the knowledge or other same information we may assume here is the model but even in that case very often a lot of information are still incomplete for example parameters could be missing and for example if you want to say try to show this remaining but you do not really have that information upon arrival maybe information is delayed or information is missing so therefore we need to deal with that one. I see this yes okay that is a talk yesterday I think that is something that has to be done in general you see in order to make the things applicable during the year we should also look at real problems and actual problems. This is not something I did well I did it in the past and there is also they used to be quite long for long time methods which was a practical problem in airline maintenance but I haven't looked at this at this moment but yes this is very important that you look at the statistical method that you do that type of thing that you look at actually during problems and so forth and what I look at is easy methods that once you have formulated such a system how can you solve it without going in a lot of mathematics so this would be the second part kind of a problem and the first part which you alluded to the statistics and so forth that would be the first part. Thank you Winfried oh by the way I have some I have few words but not related to Winfried that today's talk and most of you know Winfried actually when I was doing my PhD Winfried was my PhD supervisor and actually my career greatly impacted by the supervision of Winfried and also the directions I have been doing also greatly influenced impacted by Winfried work for example the state of reduction method that the the GTH method randomization and factorization or other things so I really appreciate I like to use this opportunity to thank Winfried once again thank you Winfried oh thank you for this nice comment. Hello hi yes hi Winfried very nice insight for talk thank you very much my question is about the accuracy of simulation method and related the issue is is it possible to use or when you use simulation how do you identify system singularity and do you have comments on those things. Okay yes well okay well first of all simulation since the error increases with the square root of the number experiments simulation in order to be very accurate in spite a lot of runs but this is with the computer speeds that you have now this is no great problem however there are cases the simulation is doing poorly and that is exactly if you have rare events or I don't remember you mentioned the simulation in the event is poor because the error is kind of well it's an absolute error and if you estimate small probabilities you really need very very many experiments because this is yes this is a problem it's a major problem and there is a group to address exactly this point. Okay thank you. Hello my run Hlemka here you were talking about these huge state spaces sometimes I have an issue that I don't even know how to order the states is that an issue that you come up with. That is an issue yes and if you have ideas how to do that I would like to hear about it. You see I mean there's two orders that I tried and one order is the alphanomeric order which I presented and the alternative is to order as you create the states as you start the state run and see which events bring you to the different states you use these new states again to find further states and so forth. I initially thought that this is a poor method because it takes long to generate all states but this is not true you can find the states in linear time. Thank you.
 And here you see again the title and my name as you see is Winfield Grasmund and I talk  about the Lyme Monte Carlo simulation is done so frequently.  Monte Carlo simulation is really one of the most successful applications of operations  in the field search and even beyond.  However the bulk of Turing theory uses deterministic methods that is methods not suggesting to  randomness with random variables and outside Turing theory many Turing models are solved  by Monte Carlo simulation and the results of Turing theory are often ignored.  And how can we make change that?  How can we make classical Turing theory more successful?  Now the first question is why is simulation so successful?  And I see three reasons for this.  You see maybe others or more.  One for large problems Monte Carlo simulation needs far less computer time.  Two simulation is much much easier.  Most people have no difficulty doing simulation and three simulation is much more flexible.  Now how can Turing theory come to this?  Let us look at the different issues, the first issues the execution times.  In Monte Carlo simulation execution times start at a high level but they increase linearly  with more complexity as measured say by the number of variables you use.  In deterministic methods execution times start at a low level but they increase exponentially  with more complexity.  So Turing theory is an advantage for small models but not for large models.  Now to show you this I have here this graph.  The red line is the deterministic methods.  They start at a low level and go up very steeply and very soon they take more computer time  than simulation, the blue line.  Our simulation needs a lot of computer time initially and in fact before the advent of  computer simulation was not a practical method.  However the computer time increases very slowly with a number of values.  Now the next thing is mathematics.  How can we simplify mathematics?  I think we should express numeric methods.  Maybe that is a little bit against the stream to say that people want results rather than  formulas.  Now if you use numeric methods you have a choice.  You can use the standard methods or you can use methods tailored to deal with probabilities.  Now for standard methods you have lots of software available including software included  in simple manipulation languages.  However in some cases it is advantageous to exploit the special properties of probabilities  mainly that probabilities are also always greater than or equal to zero.  And this enables us to avoid subtraction and just reduce rounding errors.  Now I should be clear here, subtractions do not cause rounding errors but they increase  rounding errors committed earlier, cases quite rapidly often by orders of magnitude.  Now if you use special algorithms to deal with probabilities they often have a probabilistic  interpretation and that provides additional insights.  Now the third thing is can we increase the flexibility?  One of the most used simulation paradigms is discrete event simulation.  Can we copy this paradigm?  Now we can use discrete event systems as they were described by Cassandra and Laportzli  in 2008 in a book.  And what are discrete event systems?  Well the discrete event systems has two main components, the state and events.  Now the state is a set of fibers, for instance, q-legs.  And all changes occur in discrete steps and discrete event simulation.  Now every change is called by an event, because an event can even be defined as something  that changes the state.  Now q's are clearly discrete event systems, why not use this paradigm during q-legs?  Now first let us look at the state.  First the state consists of physical state fibers, fibers that you are interested in,  it could be q's because we assume that the state fibers are x1, x2, up to xd.  And we also assume that the xi are between 0 and ni.  In addition to this we need supplementary state fibers, if you have a schedule of future  events you need to keep track of the schedule and for this you use a supplementary state  fiber or you may record when past events have happened, again that would be supplementary  fibers.  What you need is enough supplementary values such that the distribution of the future state  fibers can be calculated.  And you do that and the system becomes a Markov chain consisting of physical and supplementary  state fibers.  Next we come to the events.  Events come in different types such as arrivals, departures, switches of rating lines and so  forth.  And each event is characterized by three things, an event function, an event condition and  the process that generates the events.  The first the event function we call it fkx and if the state before event k is x, the  state after the event k is fkx.  For instance if your state is x1, x2 and you have an event function fkx1, x2 which increases  x1 by 1 leaving x2 unchanged and this would be the event function of something like an  arrival to line 1.  Next we have an event condition and events are prevented from occurring if the condition  is not met.  And then we have an event process that could be a process, renewal process, a phase start  process and hold action.  Now here we have an example of a discrete event systems.  There are two bins, bin 1 and bin 2.  They have finite size, bin 1 the size n1, bin 2 the size n2 and arrivals to bin 1 are  lambda 1, arrivals to bin 2 are lambda 2.  The bins are used to do repairs and each repair always needs one unit from bin 1 and two units  from bin 2 repairs occur at a late new and all these events are possible.  Now this system has two state parameters x1, the number in bin 1 and x2, the number in  bin 2.  And what are the events?  We have three events, we finish one, we have two and repair and each event has an event  function, an event condition and a rate.  So we finish one, we have a path, we see u increase x1 by 1 leaving x2 unchanged and  this has the event condition x1 is less than n1 because if the thing is full you don't  accept arrivals anymore and the rate is lambda 1.  Similarly you have event number 2, it simulates as a rate of lambda 2 and repair has the event  function x1 minus 1, x2 minus 2 and this event can only happen if you have enough in the  bins, that is x1 must be greater than 0 and x2 must be greater than 1 and if that is 2  then the event happens at a late new.  So now you have formulated this event and now you need a transition matrix because the  whole thing is a Markov thing.  So we take the transition matrix Q and I do it here for any Quasso events because we  have to find Qij, the rates of going from state i to state j except for the diagonal  which you all know must be determined such that the sum across the row must be 0.  Because you take a state out of the state space and you apply all events and for event  k you do the following, if the event condition holds, if Ck holds and you have a new state  after the event you call it x new which is xk of the present state and the transition  again is x, x new which is set to the rate of event k and this you do for all events  and for each state.  Now this is not really a transition matrix because the x are vectors rather than numbers  because we have to associate our vectors with numbers because we have to find for each state  a number and this number is used to find both the row and the column of the state and this  number must be in a one to one relation with the state.  Now if the state space is as simple as this which is true for many events including the  for many discrete event systems including the event systems I talked about then you  can use the alphanumeric order.  This that is the event number is the rank in the alphabetic order the state has and this  event number can be found by the following formula.  You take the sum of uixi where the ui are the products ntray plus one or tray equals i plus  one to d.  Okay so now we have events event number and in some cases you can trust enumerate the  event numbers and you get a matrix gather quickly.  Now things are more complicated if events are generated by renewal process then we must  add supplementary variables such as the time since the previous event of the same type  or the time to the next event if the event is scheduled.  Now the times between events are often continuous random variables and in this case we need  to discretize unless we use phase type variables.  Now discretization has a disadvantage it forces us to deal with multiple events that occur  simultaneously and that increases the number of entries in the transition matrix by considerably.  Now if you look at equilibrium solutions we can save computer time by embedding the system  at the points where the events do occur.  Now once you have a generated your transition matrix you want to find the transcendent equilibrium  probabilities that is we want to find the ijt the probabilities to be in state j at country  and you probably all know how this is done.  We have a recursive formula here that ijt is ijt plus one is a function of iitpij.  Now as you all know if t is large enough then the pijt corresponds to an equilibrium probability  pij and in order to find the equilibrium probability we have a system of situations which is obtained  from this system by just popping the argument t.  Now for continuous time of trains we have a similar relationship we have a rate matrix  or the increment generator q.  Now instead of a recurrence relation we have a dependency equation because pij prime t  is the derivative of pij with respect to t and here you have a formula for pij prime  t.  Again it is large enough that converges to an equilibrium probability which is written  down here.  And now the next step is to solve these equations both as transient equations and equilibrium  equations numerically because this is what I propose and let us do that.  Now first however let us look at the properties of the transition matrix.  The matrices are huge but sparse.  If all events are Poisson processes the matrix size is given by the ranges of the different  state drivers and the product of the range is of the different state drivers.  Now the density however is low because the entries per row in the transition matrix is  equal to the number of events and let the number of events be e.  Consequently the density is e divided by n.  Now if the events are renewal processes you also have to add the products of the ranges  of your supplementary drivers which I call yi and then you get the matrix size of this  huge product.  Now in each row we actually have two to the e events.  It is two to the e because we have to look at all combinations of events that can have  an entity.  This is an approximation but it is good for m I know too small.  Because we have huge matrices.  The matrices are also banded.  In fact if x1 changes at most by 1 the state number and that represents the band and it  changes by u1 which is the product of the ranges of the different state drivers except  for x1.  It is n1 plus 1.  Now here I have an example.  I have four state drivers and the ranges of the state drivers are from 0 to 9 and if you  have all the elements plus or that gives 10 to the fourth of 10,000 states.  However the transition matrix has 10,000 to the square entries which is 100 million.  And if you assume five events as it's true for the tandem queue you have only five times  10,000 that is 15,000 non-zero entries.  Now the same is true for a lesser to a lesser extent also for men's generated value of u  reposes.  This suggests that we store only the non-zero entries of the transition matrix.  And I apply this idea now in order to find transient solutions of the discrete time mark  of time.  Thus we use three arrays, one for the row, one for the column and one for the probabilities.  Thus we want to apply now this method to calculate pjt plus 1 according to this formula.  This is a formula for transient probabilities in this discrete time mark of time.  We assume that the total number of non-zero entries we call the dn-top which in some cases  in most cases is 10 times the number of events.  Thus this is the number of non-zero entries and we want to do the calculation for p going  from 1 to t top and here you have the problem to do that.  We go from t for 1 to t top for n from 1 to n top.  We go to all the non-zero events.  And this formula really represents this formula here.  The pi nu tray is ijt plus 1, the pi ti is iit and this is the probability ij.  Now i is of course what you find in row i, tray is what you find in row n, tray is what  you find in column n and the probability is the next.  Just to do that for all n.  And the remaining statement should be clear because I don't want to go into it.  But the important thing is here is that the complexity is proportional to n.  This is the important thing.  Now we can look at Planck's solution in continuous time.  In this case you could solve the differential equations but we suggest another method called  randomization or uniformization.  It is pretty standard now.  In this method you create a matrix p which is q divided by f plus the identity matrix.  Where f is chosen such that p is a stochastic matrix.  And this we can look at the Markov chain we have described by p.  And after n plus 1 steps to the probability to be in state tray is this formula here.  But there is now a difference.  The steps are not steps of one.  They are steps of exponential random bias with rate of f and rate of f.  And consequently the pjg the probability to be in the continuous Markov chain in state  tray at p is given by this formula which has a pj super scheme n multiplied by the force  of distribution.  And as I said this is the randomization or uniformization.  This method can also be derived by a small modification of the Taylor series.  So here we have an effect on mathematics.  Now Taylor expansions of matrix exponentials tend to be numerically unstable.  And it was argued that no that first hand that this method is potentially unstable.  But this is not true and the reason is that we have no subtractions.  And if there is no subtractions rounding error think is very slow.  So randomization or uniformization as this method is called works for millions of states  and it works for very high heat and out problems.  And what is even better it performs better than standard algorithm for solving the effect.  So far transient solution and now we come to steady state solutions.  We discuss only discrete and Markov chains.  First we have here the formula for the steady state of a discrete Markov chain.  This is a system of equations and we can use Gaussian elimination.  Now the mathematicians say Gaussian elimination is unstable if you have something like 10,000  or 100,000 states.  But this is not necessarily true because we can reformulate the system such that subtractions  are avoided.  And we do that when the resulting algorithm has a probabilistic interpretation.  Now if you look at linear algebra box then if you solve linear equations and linear equations  then you start with the first equation and solve this first equation for x1.  We don't do this.  We start with the last equation and solve it for pi n.  Well why do we do that?  Well mathematicians assure us that routing errors are minimized if you start with the  smallest n to be first.  Just using this mathematic is a mathematical recipe, we eliminate pi n first and then you  get the following result.  You have Gaussian elimination.  Here you have the system of linear equations.  We set the tray equal to n, that is that it becomes pi n and we solve for pi n.  That's the solution of the last equation for pi n.  Then we do the substitutions, that is we first split the sum into the sum from pi equals  1 to n minus 1.  Then we have here the last term which only includes pi n.  We plug in this value here and we simplify and what you have is a new system of equations  which has n minus 1 linear equations covered in n and the coefficients are now given by  pi j plus pi n pi nj divided by 1 minus pi nn.  Now you can continue.  First we can call these new coefficients pi j n minus 1 which is pi j plus pi n pi nj  divided by pi 1 minus pi nn.  This gives us a system of equations which we could solve for pi n minus 1.  Then we can get a new system of equations which is sorted for pi n minus 2 and so forth  and so forth.  The coefficients after we have solved all the pi's from n up to a capital N from small  n up to capital N are given by the following formula here.  Pi ij super skipped n minus 1 is pi ij n plus pi i n super skipped n pi nj super skipped  n divided by 1 minus pi nn super skipped n.  In order to obtain this formula we also need pi n which is this equation.  It expresses pi n in terms of the pi i where the i is less than n.  This can be used in the black substitution phase for calculating the pi n.  Now I told you there is a probabilistic interpretation of the whole thing and to discuss it we first  need to use embedded.  Now what is embedded?  Suppose you have a sequence x1, x2, x3 and so forth which would be a Markov chain and  in addition we have a set c.  And now we create a subsequence of this sequence here and elements xb of the original sequence  are included in the subsequence only if they are elements of c.  Because if the set c is 1, 2 and 3 and if xb is 1, 3, 5, 2, 4, 2 then what happens?  Well 1 is in c so we keep it, 3 is in c so we keep it, 5 is not in c so we omit it, 2  is in c we keep it, 4 is not in c we omit it and 2 is in c.  Now we can actually embed into embedded Markov chains.  First we can use this new sequence we call it x superskip c and we can embed it into  the set d where d is say 1 or 2 and so we get 1, 2, 2.  Now if you look at this you will see that instead of embedding d into this sequence  into the sequence x superskip c you could have embedded d also in the original sequence  and we would have the same result.  First what happens is if you have repeated embedding then you could as well do embedding  in one step or vice versa in order to do embedding in for many, many different states you can  break it down into one step embedding sets.  Because embedding can be done in either in several small steps or equivalently in one  large steps.  Because we can reduce the state space by one state at a time until only one state is left  and this I call state data.  Now if you have a Markov chain and if you do embedding by set c then the corresponding  embedded chain is also a Markov chain.  The problem is now to find the transition probabilities pij superskip c and what is  the solution?  The solution is to add all the probabilities of all the paths that start in i ending j and  avoid any state c in between.  Now what happens is that elimination Gaussian elimination is really embedded.  First we start with the original sequence states 1, 2, 3, we embed it into a set c which goes  from 1, 2 up to n minus 1 and in order to find c we add the probabilities of all the  possible sequences to go from i to j avoiding any state less than n.  We can go from i to j either indirectly that's pij we can go from i to n and from n to j  that's pij, pij, we can go from i to n, we can stay in n, go to j and the probability  of this path is pij, pij, pij and so on.  Now you go on, for this sequence here you find the probability pij, pij and n square  pn, pij.  Now we have to add all these probabilities together and we get these probabilities for  the embedding set 1 to n minus 1 which is equal to pij, pij, pn, n to the k, you stay  k times in state n, pnj and as you see this is a geometric series and the sum of the geometric  series of these terms is 1 minus pnn, first what we get is pij plus pij divided by 1 minus  pnn and this is nothing else but pij superskipped n minus 1 and the rest follows the same holds  now for all pij superskipped n.  To state it out you reduce the state space by 1 into 2 step and for the general term  you get pij n minus 1 is pijn, pij and n to the n to the n to the n minus pnn and this  is the transition matrix for the Markov chain embedded in the set of 1 to n minus 1.  Now here is a subtraction and I said we should avoid subtractions and this can be done because  all the pnj here are transition probabilities of Markov chains that is the sum across the  j is equal to 1 or alternatively if you look at state n, 1 minus pnn superskipped n is  the sum j equals 1 to n minus 1 pnj n such that the sum here is 1 if you move that over  and it follows that this formula can be written without using subtractions and you get the  3th method the method of TACS-Altas-Mantak-Sahay.  Now this method though has a problem it is impossible to exploit specific well it is possible  but it is very difficult. Now if the matrix is tense the number of floating point operations  goes up with the 3rd power of n. In the example I presented before we had n10 to the 4th state  that is 10,000 because this means according to the formula here that we get 2 3rds 10  to the 12th flops and if you do that or if you do these many flops on a laptop as I have  here as you have you have then that means around 100,000 seconds compute the time that  is one day on the laptop. First this is not a good thing really. Now you may say well  we cannot exploit bandedness well what do you get? You get as you can divide by the  bandwidth which is denominator here but still you get an n to the 3rd in the numerator.  First the algorithm are always n to the 3rd. First in the example you would reduce the  time from one day to 100 so but they which is still you don't like to wait too long until  you get this one. What is even worse is if you look at transient solutions they increase  this n rather than with n to the 3rd. Those transient solutions in a way are easier to  obtain than steady state solutions. Now this is kind of strange but for large-stance matrices  this is good. In fact in order to find equilibrium solutions you can just iterate the transient  solutions until you have reached the equilibrium. Now you can increase convergence a little  bit by different tips such as using the trumpets and so forth. You can also go back to classical  methods such as Gauss Seidel I use that extensively and the 3rd by 12. Now in Gauss Seidel the  order of the states is important and as has been shown by Mitter and Zuckers if you order  the states such that the events downstream changing state by the downstream are done  first then the convergence is faster. First you look at the flow. Again you have to look  at the meaning of the model and by exploiting the properties of the model you can improve  your methods. Now state reduction is still great for some theoretical results and I have  here two theoretical results that mention the 3rd one. First it allows you to deal with  infinite state spaces. Secondly it allows you to deal with matrices that have a discrete  cluster and third they can allow you to do good embedding techniques which I don't have  the time to discuss. Now first infinite state spaces. Here there is the following result.  If J is recurrent and the state space can be cut in such a way that the Pij superskips  the N. The only thing that you need to calculate the Pij N does not change by more than epsilon  as long as xd plus 1 minus xp the up comes about it. And here is the proof. If J is still  recurrent most parts from i to j ending as soon as j is reached have a finite length.  Because we now only look at most parts and most parts as we can say these are the parts  less than u and if you look only at those parts less than u then the number of states  reachable between two visits of J is finite and the states that cannot be reached on such  a path can be ignored and we can cut them without changing Pij superskips N by more  than epsilon. First in other words you can always cut the state space to a finite size  if the Markov chain is recurrent and if the up comes about it. Now the next thing is let's  look at transition matrices with repeating columns. Now here I have a very simple proof  but you find more complex proofs in the literature. If Pij is equal to Pij j plus k for almost  all ij which is true for all non methods then we have Pij N equals Pij j plus k j plus n  plus k and the proof is simple. If you have any path from i to j superskips N you can  match that with a path from i plus k j plus k n plus k by just shifting the path upwards  by k units and when doing this you have you meet the same transition over the next. Consequently  this is our next. Now this essentially finishes my talk and here I have some conclusions to  go back to the beginning so to say. You make a previous conclusions you may easily what  I think these things should be discussed. During this flexibility I propose this big event  approaches to reduce computer time. You can try to reduce the number of state values and  to make mathematics easier use your method possibly modified to do that for a previous  case. Okay that is the all here we have number of references which I found useful and of  course my own papers I will use at least for you. Okay any question?  Pij Jank please. Yes I am here. Yes I am here. I am hearing you. So we are afraid it is very  nice very interesting talk. So when we do numerical and also the Markov for Mendecavo  what is your comments on a statistical approach in the future for queuing models? You see  I mean what is my suggestion for the future of queuing theory is that the question?  Nowadays say the advances of this telecommunications smartphones and all this real-time apps and  a lot of times we do not have really have models and we have input we have something  and then we have output but in order to do the analysis or to design the system in an  optimum way of course sometimes based on the knowledge or other same information we may  assume here is the model but even in that case very often a lot of information are still  incomplete for example parameters could be missing and for example if you want to say  try to show this remaining but you do not really have that information upon arrival maybe  information is delayed or information is missing so therefore we need to deal with that one.  I see this yes okay that is a talk yesterday I think that is something that has to be done  in general you see in order to make the things applicable during the year we should also  look at real problems and actual problems. This is not something I did well I did it  in the past and there is also they used to be quite long for long time methods which  was a practical problem in airline maintenance but I haven't looked at this at this moment  but yes this is very important that you look at the statistical method that you do that  type of thing that you look at actually during problems and so forth and what I look at is  easy methods that once you have formulated such a system how can you solve it without  going in a lot of mathematics so this would be the second part kind of a problem and the  first part which you alluded to the statistics and so forth that would be the first part.  Thank you Winfried oh by the way I have some I have few words but not related to Winfried  that today's talk and most of you know Winfried actually when I was doing my PhD Winfried  was my PhD supervisor and actually my career greatly impacted by the supervision of Winfried  and also the directions I have been doing also greatly influenced impacted by Winfried  work for example the state of reduction method that the the GTH method randomization and  factorization or other things so I really appreciate I like to use this opportunity to thank  Winfried once again thank you Winfried oh thank you for this nice comment.  Hello hi yes hi Winfried very nice insight for talk thank you very much my question is  about the accuracy of simulation method and related the issue is is it possible to use  or when you use simulation how do you identify system singularity and do you have comments  on those things.  Okay yes well okay well first of all simulation since the error increases with the square  root of the number experiments simulation in order to be very accurate in spite a lot  of runs but this is with the computer speeds that you have now this is no great problem  however there are cases the simulation is doing poorly and that is exactly if you have  rare events or I don't remember you mentioned the simulation in the event is poor because  the error is kind of well it's an absolute error and if you estimate small probabilities  you really need very very many experiments because this is yes this is a problem it's  a major problem and there is a group to address exactly this point.  Okay thank you.  Hello my run Hlemka here you were talking about these huge state spaces sometimes I have an  issue that I don't even know how to order the states is that an issue that you come  up with.  That is an issue yes and if you have ideas how to do that I would like to hear about  it.  You see I mean there's two orders that I tried and one order is the alphanomeric order which  I presented and the alternative is to order as you create the states as you start the  state run and see which events bring you to the different states you use these new states  again to find further states and so forth.  I initially thought that this is a poor method because it takes long to generate all states  but this is not true you can find the states in linear time.  Thank you.
202008220832-Grassmann_lrv_0_203.png 202008220832-Grassmann_lrv_204_347.png 202008220832-Grassmann_lrv_348_1595.png 202008220832-Grassmann_lrv_1596_2495.png 202008220832-Grassmann_lrv_2496_3443.png 202008220832-Grassmann_lrv_3444_3863.png 202008220832-Grassmann_lrv_3864_4163.png 202008220832-Grassmann_lrv_4164_4331.png 202008220832-Grassmann_lrv_4332_6803.png 202008220832-Grassmann_lrv_6804_8531.png 202008220832-Grassmann_lrv_8532_10331.png 202008220832-Grassmann_lrv_10332_12431.png 202008220832-Grassmann_lrv_12432_14063.png 202008220832-Grassmann_lrv_14064_14375.png 202008220832-Grassmann_lrv_14376_14423.png 202008220832-Grassmann_lrv_14424_14615.png 202008220832-Grassmann_lrv_14616_14759.png 202008220832-Grassmann_lrv_14760_15203.png 202008220832-Grassmann_lrv_15204_15443.png 202008220832-Grassmann_lrv_15444_16115.png 202008220832-Grassmann_lrv_16116_16847.png 202008220832-Grassmann_lrv_16848_16895.png 202008220832-Grassmann_lrv_16896_16955.png 202008220832-Grassmann_lrv_16956_17399.png 202008220832-Grassmann_lrv_17400_18611.png 202008220832-Grassmann_lrv_18612_18671.png 202008220832-Grassmann_lrv_18672_19199.png 202008220832-Grassmann_lrv_19200_19487.png 202008220832-Grassmann_lrv_19488_20579.png 202008220832-Grassmann_lrv_20580_20879.png 202008220832-Grassmann_lrv_20880_21491.png 202008220832-Grassmann_lrv_21492_23963.png 202008220832-Grassmann_lrv_23964_26039.png 202008220832-Grassmann_lrv_26040_26267.png 202008220832-Grassmann_lrv_26268_26771.png 202008220832-Grassmann_lrv_26772_28427.png 202008220832-Grassmann_lrv_28428_31703.png 202008220832-Grassmann_lrv_31704_34271.png 202008220832-Grassmann_lrv_34272_35723.png 202008220832-Grassmann_lrv_35724_38207.png 202008220832-Grassmann_lrv_38208_39935.png 202008220832-Grassmann_lrv_39936_42143.png 202008220832-Grassmann_lrv_42144_46067.png 202008220832-Grassmann_lrv_46068_46319.png 202008220832-Grassmann_lrv_46320_47243.png 202008220832-Grassmann_lrv_47244_47291.png 202008220832-Grassmann_lrv_47292_50867.png 202008220832-Grassmann_lrv_50868_52907.png 202008220832-Grassmann_lrv_52908_56075.png 202008220832-Grassmann_lrv_56076_57491.png 202008220832-Grassmann_lrv_57492_57851.png 202008220832-Grassmann_lrv_57852_58067.png 202008220832-Grassmann_lrv_58068_58403.png 202008220832-Grassmann_lrv_58404_58763.png 202008220832-Grassmann_lrv_58764_61331.png 202008220832-Grassmann_lrv_61332_63095.png 202008220832-Grassmann_lrv_63096_64163.png 202008220832-Grassmann_lrv_64164_64751.png 202008220832-Grassmann_lrv_64752_64847.png 202008220832-Grassmann_lrv_64848_64931.png 202008220832-Grassmann_lrv_64932_65903.png 202008220832-Grassmann_lrv_65904_66251.png 202008220832-Grassmann_lrv_66252_68483.png 202008220832-Grassmann_lrv_68484_70355.png 202008220832-Grassmann_lrv_70356_70475.png 202008220832-Grassmann_lrv_70476_72131.png 202008220832-Grassmann_lrv_72132_73199.png 202008220832-Grassmann_lrv_73200_73283.png 202008220832-Grassmann_lrv_73284_75155.png 202008220832-Grassmann_lrv_75156_75287.png 202008220832-Grassmann_lrv_75288_75659.png 202008220832-Grassmann_lrv_75660_77051.png 202008220832-Grassmann_lrv_77052_77135.png 202008220832-Grassmann_lrv_77136_77172.png