We investigate the performance of a production system with correlated demand through diffusion approximation. The key performance metric under consideration is the extreme points that this system can reach. This problem is mapped to a problem of characterizing the joint probability density of a two-dimensional Brownian motion and its coordinate running maximum. To achieve this goal, we obtain the stationary distribution of a reflected Brownian motion within the positive quarter-plane, which is of independent interest, through investigating a solution of an extended Helmhotz equation.

1. Introduction

There are extensive studies on some classic one-dimensional models in the field of operations research, most notably, the work on the behavior and scheduling of single-server queuing system, as well as the understanding of performance and management of single-item inventory system. Probabilistic tools and techniques such as random walks and integral transformations are traditionally used in analyzing them. More recently, concepts and methodologies in dynamical systems and diffusion processes are brought in through fluid and diffusion approximations, and they extend our understanding and capability of analyzing and control of these systems substantially.

While highly desirable, extending these studies to their multidimensional counterparts is a rather difficult problem. Extensions of classic probabilistic methods and techniques, such as random walk and integral transform, introduce a new level of complication that requires deeper understanding in algebraic and complex geometry for their general treatment. For the approximations methods, their multidimensional counterparts, multidimensional dynamical systems, and multidimensional diffusion processes, also pose significant barriers for either qualitative understanding and quantitative computation. While a general and sophisticated methodology is lacking, some very interesting results have been obtained in establishing convergence results that lay the foundation to both fluid and diffusion approximation schemes. Meanwhile, the quantitative aspect of the problem seriously lagged behind, and results in approximation and control of multidimensional systems are very limited. This calls for more focus be put on computational efforts on such system, so that new tools and techniques can be added to the arsenal of attacking these problem.

In this paper, we aim at extend some theoretical and computational understanding to a simple but representative two-dimensional operations research model. Not only can it serve as a building block for analyzing much more complicated systems, but also provides insights to the understanding of basic structure of general systems with correlated demand. This model is motivated by several typical applications in production systems and scheduling. In such a system, it is very common that one resource is demanded by multiple demand processes, similarly, one class of customer demand could contain a combination of different products. The correlation induced by this commonality poses difficulty even for the simplest version of the model. In this paper, after briefly introducing the mathematical model and stochastic processes that characterize the basic underlying relationship, we provide an diffusion approximation to the stochastic processes in the form of a reflected two-dimensional Brownian motion, which is just an abstraction of previous known results in various forms. Then our main effort will be concentrated on the analysis of some key quantities of this diffusion process. More specifically, we are interested at the joint probability density of the Brownian motion and its coordinate running maximum. To produce the computational result, we relate its computation to finding a proper solution to a classic elliptic partial differential equation, the extended Helmhotz equation.

The rest of the paper will be organized as follows, in Section 2, we will introduce the mathematical model and some preliminary results; in Section 3 we will present the a key Brownian bridge argument; then, in Section 4, we solve an extended Helmhotz equation related to the invariant measure; in Section 5, we present the final calculation results for the probability density, and we conclude the paper with come comments on related future works in Section 6.

2. Model

Suppose that there are two types of products, type and type . Each of them is supplied by an independent renewal process. Let us denote them as and , respectively. Statistically, we only require that the interarrival time follows a distribution with finite second moments. Meanwhile, there are two classes of demand, a class- demand requires one unit of type product, a class- demand requires one unit of type product, and a class- demand requires one unit of each product and . Note that class- demand can only be fulfilled when both two types of products are provided. If any of the product requirement can not be satisfied, the demand will be backlogged. Although the model is quite simplistic on the surface, it embodies the basic structure and trade-off for many widely considered systems, for example, assemble-to-order systems in inventory management and queuing systems with flexibility servers. Results on this problem can have many indications to the understanding of those more complicated systems.

Diffusion approximation can be obtained for this system, see, for example [1, 2], the limiting process is a two-dimensional Brownian motion with correlation. One of the most important problems for characterizing the process is the probabilistic distribution of the coordinate running maximum. For many applications, the coordinate running maximums are closely related to the capacity provisioning since they represent how far the process will reach at different directions. For probability theory, it is also a fundamental problem. And this will be the focus of the paper.

From now on, we will focus on the calculation of the probability density function calculation. More specifically, let be a two-dimensional Brownian motion with constant correlation and initial state . Therefore, we have , , for , and . For detailed properties of such process, see, for example, [3]. A quantity that of general interest is, for and and any fixed time where are the running maximum for for respectively, that is . It is easy to see that we can just focus on the case of , general case can be easily obtained through scaling. To be more specific, we have, for any , because and have the same law due to the scaling invariant property of Brownian motion.

In the case of one-dimensional Brownian motion, the density of the Brownian motion and its running maximum is a classical result obtained through reflection principle. In two dimension, the dependence between the two coordinates creates a barrier for generalizing the result easily. In this paper, first, we convert the calculation of the quantity (2.1) to the calculation of density function of the invariant measure for a reflected two-dimensional Brownian motion in the positive orthant with negative drifts through a Brownian bridge argument; then we obtain a Laplace transform of the invariant measure of a reflected two-dimensional Brownian motion in the positive orthant with negative drifts through solving a partial differential equation.

3. A Brownian Bridge Argument

There are many different approaches to the classic problem of the joint distribution of a standard Brownian motion and its running maximum, see, for example, [4]. Here we introduce an approach connecting probability laws of a conditional Brownian motion to a Brownian bridge. To the best of our knowledge, no similar approaches appeared before in the literature. A standard Brownian motion, , which starts at , conditioning on the event that it ends at at time , has the same law as where is a standard Brownian bridge, that is, a Brownian motion satisfies . Therefore, from the representation of Brownian bridge, we have, for , where is a standard Brownian motion, means equal in distribution. Thus, for any , the event of will have the same probability as, Apply a variable transform , we get, or equivalently, This probability of the above event is equivalent to the stationary distribution of reflected at zero, which, in turn, equals Take derivative with respect to , we have, So, This is consistent with the classic results obtained from other approaches such as the reflection principle.

For the two-dimensional Brownian motion, applying the same argument, calculating for and can be reduced to calculating Apply the construction of reflected Brownian motion through maximum of Brownian motion with negative drift, see, for example, [5], as well as the two-dimensional version of P. Lévy theorem, see [3], we can conclude that (3.11) is the same as where is a two-dimensional Brownian motion with negative drifts that reflected at the two positive axes. To be more precise, where is a two-dimensional driftless Brownian motion with covariance matrix, and are the local time defined as where , . Therefore, following the same logic, to obtain (3.10), we need to compute the stationary distribution of .

4. Stationary Distribution

To calculate the stationary distribution of the reflected Brownian motion , we make use the following result from Harrison and Williams [6].

Theorem 4.1. Functions and are the density of the stationary distribution of on the orthant and the boundary, if and only if the following basic adjoint relationship (BAR) is satisfied: where is the infinitesimal generator and is the surface measures defined on the boundary such that the normal vector is pointed inward.

Remark 4.2. Notice that the boundary in the problem we study consists of the two axis. So the surface measure is basically the Lebesgue measure on with the orientation that guarantees that the normal vector points towards the inside of the orthant.

The main idea in this paper is to find a function such that for any pair , therefore, the first term of (4.1) becomes the Laplace transformation of the stationary distribution, and then to investigate the basic properties of to compute the second term. To find , we need to solve the following equation: Let us first consider the homogeneous version of the equation, Apply transform, then, satisfies where and . Let , then is the solution of the following problem, This is so called extended Helmholtz equation. Separation of variables can be applied to obtain the general forms of the solution. Hence, we know that the general solution is for any constant , and . Therefore, For the inhomogeneous version, it is easy to see that will be the solution.

Now plug into the BAR (4.1), we have where and denote the Laplace transform of and . Given that , and can be arbitrary, it is necessary that In fact, they are all in the form of and and satisfy the quadratic equation which is equivalent to For any , denote as the positive conjugate point of with respect to the curve defined by (4.14). Therefore,

Assume that . It is easy to see that the equation has positive solution only when . Hence, we have Therefore, from (4.10), we obtain the following.

Theorem 4.3. The Laplace transform of the stationary bears the following form:

5. Final Calculations

Recall that our goal is to get an expression for which is, of course, The results in the previous section enable us to obtain the Laplace transform of Apply the Laplace transform with respect to , we get where is the Laplace transform of . Hence, we obtain the following.

Theorem 5.1. The laplace transform of with respect has the form,

6. Conclusions

We studied a simple but representative operations research model for performance analysis and capacity planning. Through diffusion approximation, we identify the problem of characterizing the joint distribution of a two-dimensional Brownian motion and its coordinate running maximum, which is of independent interests in the theory of probability, as well as other applications. Using some probabilistic techniques, we are able to reduce the problem to the solution of an extended Helmholtz equation, and a careful exam of the solution to this well-known partial differential equation leads to the calculation of the Laplace transform of the desired joint probability distribution.

Standard inversion methods exist for estimating the probabilities from the Laplace transform, in the ongoing research, we are exploring the special structure of our solution, and tailor the inversion method to our special needs, thus provide a very efficient and accurate way of estimating these fundamentally important quantities for both probability theory and operations research.