We develop a high-fidelity simulation model of the patient arrival process to an endocrinology clinic by carefully examining appointment and arrival data from that clinic. The data include the time that the appointment was originally made as well as the time that the patient actually arrived, as well as if the patient did not arrive at all, in addition to the scheduled appointment time. We take a data-based approach, specifying the schedule for each day by its value at the end of the previous day. This data-based approach shows that the schedule for a given day evolves randomly over time. Indeed, in addition to three recognized sources of variability—(i) no-shows, (ii) extra unscheduled arrivals, and (iii) deviations in the actual arrival times from the scheduled times—we find that the primary source of variability in the arrival process is variability in the daily schedule itself. Even though service systems with arrivals by appointment can differ in many ways, we think that our data-based approach to modeling the clinic arrival process can be a guideline or template for constructing high-fidelity simulation models for other arrival processes generated by appointments. The online supplement is available at https://doi.org/10.1287/ijoc.2017.0773 .

An efficient algorithm is developed to calculate the periodic steady-state distribution and moments of the remaining workload Wy at time yc within a cycle of length c, 0 ≤ y < 1, in a single-server queue with a periodic arrival-rate function. The algorithm applies exactly to the GIt/GI/1 model, where the arrival process is a time-transformation of a renewal process. A new representation of Wy makes it possible to apply a modification of the classic rare-event simulation for the stationary GI/GI/1 model exploiting importance sampling using an exponential change of measure. We establish bounds between the periodic workload and the stationary workload with the average arrival rate that enable us to prove that the relative error in estimates of P(Wy > b) is uniformly bounded in b. With the aid of a recent heavy-traffic limit theorem, the algorithm also applies to compute the periodic steady-state distribution of (i) reflected periodic Brownian motion (RPBM) by considering appropriately scaled GIt/GI/1 models and (ii) a large class of general Gt/G/1 queues by approximating by GIt/GI/1 models with the same heavy-traffic limit. Simulation examples demonstrate the accuracy and efficiency of the algorithm for both GIt/GI/1 queues and RPBM. The online supplement is available at https://doi.org/10.1287/ijoc.2017.0766 .

We consider an automatic overload control for two large service systems modeled as multiserver queues such as call centers. We assume that the two systems are designed to operate independently, but want to help each other respond to unexpected overloads. The proposed overload control automatically activates sharing (sending some customers from one system to the other) once a ratio of the queue lengths in the two systems crosses an activation threshold (with ratio and activation threshold parameters for each direction). In this paper, we are primarily concerned with ensuring that the system recovers rapidly after the overload is over, either because (i) the two systems return to normal loading or (ii) the direction of the overload suddenly shifts in the opposite direction. To achieve rapid recovery, we introduce lower thresholds for the queue ratios, below which one-way sharing is released. As a basis for studying the complex dynamics, we develop a new six-dimensional fluid approximation for a system with time-varying arrival rates, extending a previous fluid approximation involving a stochastic averaging principle. We conduct simulations to confirm that the new algorithm is effective for predicting the system performance and choosing effective control parameters. The simulation and the algorithm show that the system can experience an inefficient nearly periodic behavior, corresponding to an oscillating equilibrium (congestion collapse) if the sharing is strongly inefficient and the control parameters are set inappropriately.

Motivated by large-scale service systems with network structure, we introduced in a previous paper a time-varying open network of many-server fluid queues with customer abandonment from each queue and time-varying proportional routing among the queues, and showed how performance functions can be determined. The deterministic fluid model serves as an approximation for the corresponding non-Markovian stochastic network of many-server queues with Markovian routing, experiencing periods of overloading at the queues. In this paper we develop a new algorithm for the previous model and generalize the model to include non-exponential service-time distributions. In this paper we report results of implementing the algorithms and studying their computational complexity. We also conduct simulation experiments to confirm that the algorithms are effective in computing the performance functions and that these performance functions provide useful approximations for the corresponding stochastic models.

This paper investigates ways to create algorithms to invert Laplace transforms numerically within a unified framework proposed by Abate and Whitt (2006). That framework approximates the desired function value by a finite linear combination of transform values, depending on parameters called weights and nodes, which are initially left unspecified. Alternative parameter sets, and thus algorithms, are generated and evaluated here by considering power test functions. Real weights for a real-variable power algorithm are found for specified real powers and positive real nodes by solving a system of linear equations involving a generalized Vandermonde matrix, using Mathematica. The resulting power algorithms are shown to be effective, with the parameter choice being tunable to the transform being inverted. The powers can be advantageously chosen from series expansions of the transform. Experiments show that the power algorithms are robust in the nodes; it suffices to use the first n positive integers. The power test functions also provide a useful way to evaluate the performance of other algorithms.

We introduce and investigate a framework for constructing algorithms to invert Laplace transforms numerically. Given a Laplace transform \hat{f} of a complex-valued function of a nonnegative real-variable, f, the function f is approximated by a finite linear combination of the transform values; i.e., we use the inversion formula f(t) \approx f_n (t) \equiv \frac{1}{t} \sum_{k = 0}^{n}\omega_{k}\hat{f}\biggl(\frac{\alpha_{k}}{t}\biggr),\quad 0 < t < \infty, where the weights ω k and nodes α k are complex numbers, which depend on n, but do not depend on the transform \hat{f} or the time argument t. Many different algorithms can be put into this framework, because it remains to specify the weights and nodes. We examine three one-dimensional inversion routines in this framework: the Gaver-Stehfest algorithm, a version of the Fourier-series method with Euler summation, and a version of the Talbot algorithm, which is based on deforming the contour in the Bromwich inversion integral. We show that these three building blocks can be combined to produce different algorithms for numerically inverting two-dimensional Laplace transforms, again all depending on the single parameter n. We show that it can be advantageous to use different one-dimensional algorithms in the inner and outer loops.

It is often possible to effectively calculate probability density functions (pdf's) and cumulative distribution functions (cdf's) by numerically inverting Laplace transforms. However, to do so it is necessary to compute the Laplace transform values. Unfortunately, convenient explicit expressions for required transforms are often unavailable for component pdf's in a probability model. In that event, we show that it is sometimes possible to find continued-fraction representations for required Laplace transforms that can serve as a basis for computing the transform values needed in the inversion algorithm. This property is very likely to prevail for completely monotone pdf's, because their Laplace transforms have special continued fractions called S fractions, which have desirable convergence properties. We illustrate the approach by considering applications to compute first-passage-time cdf's in birth-and-death processes and various cdf's with non-exponential tails, which can be used to model service-time cdf's in queueing models. Included among these cdf's is the Pareto cdf.

In order to provide an adequate quality of service to large-bandwidth calls, such as video conference calls, service providers of integrated services networks may want to allow some customers to book their calls ahead, i.e., make advance reservations. We propose a scheme for sharing resources among book-ahead (BA) calls (that announce their call holding times as well as their call initiation times upon arrival) and non-BA calls (that do not announce their holding times). It is possible to share resources without allowing any calls in progress to be interrupted, but in order to achieve a more efficient use of resources, we think that it may be desirable to occasionally allow a call in progress to be interrupted. (In practice, it may be possible to substitute service degradation, such as bit dropping or coarser encoding of video, for interruption.) Thus, we propose an admission control algorithm in which a call is admitted if an approximate interrupt probability (computed in real time) is below a threshold. Simulation experiments show that the proposed admission control algorithm can be better (i.e., yield higher total utilization or higher revenue) than alternative schemes that do not allow interruption, such as a strict partitioning of resources.

The notion of effective bandwidths has provided a useful practical framework for connection admission control and capacity planning in high-speed communication networks. The associated admissible set with a single linear boundary makes it possible to apply stochastic-loss-network (generalized-Erlang) models for capacity planning. We consider the case of network nodes that use a priority-service discipline to support multiple classes of service, and we wish to determine an appropriate notion of effective bandwidths. Just as was done previously for the first-in first-out (FIFO) discipline, we use large-buffer asymptotics (large deviations principles) for workload tail probabilities as a theoretical basis. We let each priority class have its own buffer and its own constraint on the probability of buffer overflow. Unfortunately, however, this leads to a constraint for each priority class. Moreover, the large-buffer asymptotic theory with priority classes does not produce an admissible set with linear boundaries, but we show that it nearly does and that a natural bound on the admissible set does have this property. We propose it as an approximation for priority classes; then there is one linear constraint for each priority class. This linear-admissible-set structure implies a new notion of effective bandwidths, where a given connection is associated with multiple effective bandwidths: one for the priority level of the given connection and one for each lower priority level. This structure can be used regardless of whether the individual effective bandwidths are determined by large-buffer asymptotics or by some other method.

Traffic measurements from communication networks have shown that many quantities characterizing network performance have long-tail probability distributions, i.e., with tails that decay more slowly than exponentially. Long-tail distributions can have a dramatic effect upon performance, but it is often difficult to describe this effect in detail, because performance models with component long-tail distributions tend to be difficult to analyze. We address this problem by developing an algorithm for approximating a long-tail distribution by a finite mixture of exponentials. The fitting algorithm is recursive over time scales. At each stage, an exponential component is fit in the largest remaining time scale and then the fitted exponential component is subtracted from the distribution. Even though a mixture of exponentials has an exponential tail, it can match a long-tail distribution in the regions of primary interest when there are enough exponential components.

It is known that probability density functions and probability mass functions usually can be calculated quite easily by numerically inverting their transforms (Laplace transforms and generating functions, respectively) with the Fourier-series method. Other more general functions can be substantially more difficult to invert, because the aliasing and roundoff errors tend to be more difficult to control. In this article we propose a simple new scaling procedure for nonprobability functions that is based on transforming the given function into a probability density function or a probability mass function and transforming the point of inversion to the mean. This new scaling is even useful for probability functions, because it enables us to compute very small values at large arguments with controlled relative error.

We identify conditions under which relatively large buffers will be required in broadband communication networks. For this purpose, we analyze an infinite-capacity stochastic fluid model with a general stationary environment process (without the usual independence or Markov assumptions). With that level of generality, we are unable to establish asymptotic results, but by a very simple argument we are able to obtain a revealing lower bound on the steady-state buffer-content tail probability. The bounding argument shows that the steady-state buffer content will have a long-tail distribution when the sojourn time in a set of states with positive net input rate itself has a long-tail distribution. If a set of independent sources, each with a general stationary environment process, produces a positive net flow when all are in high-activity states, and if each of these sources has a high-activity sojourn-time distribution with a long tail, then the steady-state buffer-content distribution will have a long tail, but possibly one that decays faster than the tail for any single component source. The full buffer-content distribution can be derived in the special case of a two-state fluid model with general high- and low-activity-time distributions, assuming that successive high- and low-activity times come from independent sequences of i.i.d. random variables. In that case the buffer-content distribution will have a long tail when the high-activity-time distribution has a long tail. We illustrate by giving numerical examples of the two-state model based on numerical transform inversion.

The Laguerre method for numerically inverting Laplace transforms is an old established method based on the 1935 Tricomi–Widder theorem, which shows (under suitable regularity conditions) that the desired function can be represented as a weighted sum of Laguerre functions, where the weights are coefficients of a generating function constructed from the Laplace transform using a bilinear transformation. We present a new variant of the Laguerre method based on: (i) using our previously developed variant of the Fourier-series method to calculate the coefficients of the Laguerre generating function, (ii) developing systematic methods for scaling, and (iii) using Wynn's ϵ-algorithm to accelerate convergence of the Laguerre series when the Laguerre coefficients do not converge to zero geometrically fast. These contributions significantly expand the class of transforms that can be effectively inverted by the Laguerre method. We provide insight into the slow convergence of the Laguerre coefficients as well as propose a remedy. Before acceleration, the rate of convergence can often be determined from the Laplace transform by applying Darboux's theorem. Even when the Laguerre coefficients converge to zero geometrically fast, it can be difficult to calculate the desired functions for large arguments because of roundoff errors. We solve this problem by calculating very small Laguerre coefficients with low relative error through appropriate scaling. We also develop another acceleration technique for the case in which the Laguerre coefficients converge to zero geometrically fast. We illustrate the effectiveness of our algorithm through numerical examples.

Although ATM seems to be the wave of the future, one analysis requires that the utilization of the network be quite low. That analysis is based on asymptotic decay rates of steady-state distributions used to develop a concept of effective bandwidths for connection admission control. The present authors have developed an exact numerical algorithm that shows that the effective-bandwidth approximation can overestimate the target small blocking probabilities by several orders of magnitude when there are many sources that are more bursty than Poisson. The bad news is that the appealing simple connection admission control algorithm using effective bandwidths based solely on tail-probability asymptotic decay rates may actually not be as effective as many have hoped. The good news is that the statistical multiplexing gain on ATM networks may actually be higher than some have feared. For one example, thought to be realistic, the analysis indicates that the network actually can support twice as many sources as predicted by the effective-bandwidth approximation. The authors also show that the effective bandwidth approximation is not always conservative. Specifically, for sources less bursty than Poisson, the asymptotic constant grows exponentially in the number of sources (when they are scaled as above) and the effective-bandwidth approximation can greatly underestimate the target blocking probabilities. Finally, they develop new approximations that work much better than the pure effective-bandwidth approximation.

We present a simple algorithm for numerically inverting Laplace transforms. The algorithm is designed especially for probability cumulative distribution functions, but it applies to other functions as well. Since it does not seem possible to provide effective methods with simple general error bounds, we simultaneously use two different methods to confirm the accuracy. Both methods are variants of the Fourier-series method. The first, building on Dubner and Abate (Dubner, H., J. Abate. 1968. Numerical inversion of Laplace transforms by relating them to the finite Fourier cosine transform. JACM 15 115–123.) and Simon, Stroot, and Weiss (Simon, R. M., M. T. Stroot, G. H. Weiss. 1972. Numerical inversion of Laplace transforms with application to percentage labeled experiments. Comput. Biomed. Res. 6 596–607.), uses the Bromwich integral, the Poisson summation formula and Euler summation; the second, building on Jagerman (Jagerman, D. L. 1978. An inversion technique for the Laplace transform with applications. Bell System Tech. J. 57 669–710 and Jagerman, D. L. 1982. An inversion technique for the Laplace transform. Bell System Tech. J. 61 1995–2002.), uses the Post-Widder formula, the Poisson summation formula, and the Stehfest (Stehfest, H. 1970. Algorithm 368. Numerical inversion of Laplace transforms. Comm. ACM 13 479–490 (erratum: 13 624).) enhancement. The resulting program is short and the computational experience is encouraging. INFORMS Journal on Computing, ISSN 1091-9856, was published as ORSA Journal on Computing from 1989 to 1995 under ISSN 0899-1499.

A new algorithm is developed for calculating normalization constants (partition functions) and moments of product-form steady-state distributions of closed queuing networks and related models. The essential idea is to numerically invert the generating function of the normalization constant and related generating functions appearing in expressions for the moments. It is known that the generating function of the normalization constant often has a remarkably simple form, but numerical inversion evidently has not been considered before. For p -dimensional transforms, as occur with queuing networks having p closed chains, the algorithm recursively performs p one-dimensional inversions. The required computation grows exponentially in the dimension, but the dimension can often be reduced by exploiting conditional decomposition based on special structure. For large populations, the inversion algorithm is made more efficient by computing large sums using Euler summation. The inversion algorithm also has a very low storage requirement. A key ingredient in the inversion algorithm is scaling. An effective static scaling is developed for multichain closed queuing networks with only single-server and (optionally) infinite-server queues. An important feature of the inversion algorithm is a self-contained accuracy check, which allows the results to be verified in the absence of alternative algorithms.

The algorithm developed in Choudhury et al. (1994) for computing (exact) steady-state blocking probabilities for each class in product-form loss networks is extended to cover general state-dependent arrival and service rates. This generalization allows to consider, for the first time, a wide variety of buffered and unbuffered resource-sharing models with non-Poisson traffic, as may arise with overflows in the context of alternative routing. As before, the authors consider noncomplete-sharing policies involving upper-limit and guaranteed-minimum bounds for the different classes, but in the present paper both bounds are discussed simultaneously. These bounds are important for providing different grades of service with protection against overloads by other classes. The algorithm is based on numerically inverting the generating function of the normalization constant, which is derived in the present paper. Major features of the algorithm are: dimension reduction by elimination of nonbinding resources and by conditional decomposition based on special structure, an effective scaling algorithm to control errors in the inversion, efficient treatment of multiple classes with identical parameters and truncation of large sums. The authors show that the computational complexity of the inversion approach is usually significantly lower than the alternative recursive approach.>

The authors introduce a deterministic fluid model and two stochastic traffic models for wireless networks. The setting is a highway with multiple entrances and exits. Vehicles are classified as calling or non-calling, depending on whether they have calls in progress. The deterministic model ignores the behavior of individual vehicles and treats them as a continuous fluid, whereas the stochastic traffic models consider the random behavior of each vehicle. However, all three model's use the same two coupled partial (or ordinary) differential equations to describe the system evolution. The call density and call handoff rate (or their expected values in the stochastic models) are readily computable by solving these equations. Numerical examples are presented to illustrate how the models can be used to investigate various aspects of time and space dynamics in wireless networks. These examples also show that the models can serve as useful tools for system engineering and planning.>

Introduces a deterministic fluid model and two stochastic traffic models for wireless networks. The setting is a highway with multiple entrances and exits. Vehicles are classified as calling or noncalling, depending upon whether or not they have calls in progress. The main interest is in the calling vehicles; but noncalling vehicles are important because they can become calling vehicles if they initiate (place or receive) a call. The deterministic model ignores the behavior of individual vehicles and treats them as a continuous fluid, whereas the stochastic traffic models consider the random behavior of each vehicle. However, all three models use the same two coupled partial differential equations (PDEs) or ordinary differential equations (ODEs) to describe the evolution of the system. The call density and call handoff rate (or their expected values in the stochastic models) are readily computable by solving these equations. Since no capacity constraints are imposed in the models, these computed quantities can be regarded as offered traffic loads. The models complement each other, because the fluid model can be extended to include additional features such as capacity constraints and the interdependence between velocity and vehicular density, while the stochastic traffic model can provide probability distributions. Numerical examples are presented to illustrate how the models can be used to investigate various aspects of time and space dynamics in wireless networks.>

Rate-control throttles with token banks or leaky buckets have been used for overload control in telecommunication systems and have been recommended for traffic policing in broadband integrated services digital networks (B-ISDN). Enhancing the token-bank throttle with a buffer to shape the admitted traffic has been suggested. Researchers have shown that the presence of the buffer can dramatically reduce the squared coefficient of variation of the interadmission time. However, the authors show that the impact of the buffer on longer-time-scale characteristics of the admitted traffic is much less dramatic. In particular, they show (primarily through simulations) that the job buffer has much less impact on higher values of the index of dispersion for intervals and on small tail probabilities for the steady-state number in system at a downstream queue (with only this one arrival stream). Indeed, the smoothing benefit of the job buffer decreases as longer-time-scale characteristics become more important. However, if the downstream queue is fed by many sources with throttles, as would be the case in most applications, then the relevant time scale at the downstream queue indeed becomes relatively short. The simulation results show that the benefit of traffic shaping can be much greater. The benefit gained in reduced buffer requirements at the downstream queue, though, is typically significantly less than the sum of all job buffers added to the throttles. A full cost/benefit analysis depends on the relative cost of buffer space in the two places and on details of the relevant application.>

The authors continue an investigation of the way diverse traffic from different data applications affects the performance of packet queues. This traffic often exhibits significant dependence among successive interarrival times, among successive service times, and between interarrival times and service times, which can cause a significant degradation of performance under heavy loads (and often even under moderate loads). This dependence and its effects on performance (specifically, the mean steady-state workload) are partially characterized here by the cumulative correlations in the total input process of work, which is referred to as the index of dispersion for work (IDW). The authors evaluate approximations for the mean steady-state workload based on the IDW by making comparisons with computer simulations.>

Measurements and approximations are proposed to describe the variability of offered traffic to a queue and predict the average workload in the queue. The principal traffic measurement considered is a normalized version of the variance of the total input of work as a function of time, which is called the index of dispersion for work (IDW). Given ample traffic data, the IDW can easily be estimated using sample averages. Given a mathematical model, such as a multiclass queue in which each class has GI/G/1 offered traffic, the IDW can often be calculated analytically, or approximated by using the limits as t approaches 0 and t approaches infinity . The basic premise is that the average workload is primarily determined by the offered traffic, beyond the offered load (the deterministic rate work arrives), through the IDW. Support is provided for this premise, and it is shown how the average workload can be predicted from the IDW or basic model parameters.>

Methods are discussed for computing transient performance measures for the M/M/1 queue. These performance measures are often expressed in terms of modified Bessel functions without any discussion about computation. In fact, a common expression for the probability transition function of the M/M/1 queue length process has an infinite sum of modified Bessel functions. For actually generating numbers, however, it is convenient to use numerical integration with associated integral representations, as was first pointed out by P.M. Morse (Oper. Res., vol.3, p.255-61, 1955).>

The burstiness of the total arrival process has been previously characterized in packet network performance models by the dependence among successive interarrival times. It is shown that associated dependence among successive service times and between service times and interarrival times also can be important for packet queues involving variable packet lengths. These dependence effects are demonstrated analytically by considering a multiclass single-server queue with batch-Poisson arrival processes. For this model and more realistic models of packet queues, insight is gained from heavy-traffic limit theorems. This study indicates that all three kinds of dependence should be considered in the analysis and measurement of packet queues involving variables packet lengths. Specific measurements are proposed for real systems and simulations. This study also indicates how to predict expected packet delays under heavy loads. Finally, this study is important for understanding the limitations of procedures such as the queuing network analyzer (QNA) for approximately describing the performance of queuing networks using the techniques of aggregation and decomposition.>

This paper analyzes a model of a multiplexer for packetized voice and data. A major part of the analysis is devoted to characterizing the aggregate packet arrival process resulting from the superposition of separate voice streams. This is done via the index of dispersion for intervals (IDI), which describes the cumulative covariance among successive interarrival times. The IDI seems very promising as a measurement tool to characterize complex arrival processes. This paper also describes the delays experienced by voice and data packets in the multiplexer using relatively simple two-parameter approximations.

Abstract Customer B overtakes customer A in a queueing system if A arrives before B but B departs first. To understand the phenomenon of overtaking in a network of queues and its impact on sojourn times, it is important to develop concepts of the amount of overtaking. In this paper, two stochastic measures of the amount overtaking are defined: the random number of customers overtaken by an arbitrary customer and the random permutation of the original order of n customers upon departure. Describing these measures, however, often presents rather difficult combinatorial and probabilistic problems. In this paper, we make several conjectures about the amount of overtaking in open Jackson networks and its impact on sojourn times. Then we describe the amount of overtaking in several smaller networks, for example, a single node with the first‐come first‐served discipline and instantaneous Bernoulli feedback.