Performance Analysis of Production Systems with Correlated Demand via Diffusion Approximations

Abstract

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 1 and type 2. 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-(1) demand requires one unit of type 1 product, a class-(2) demand requires one unit of type 1 product, and a class-(1,2) demand requires one unit of each product 1 and 2. Note that class-(1,2) 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 𝜌>0 and initial state (𝑋₁(0),𝑋₂(0))=(0,0). Therefore, we have 𝖤[𝑋_𝑖(𝑡)]=0, 𝖤[𝑋²_𝑖(𝑡)]=𝑡, 𝖤[𝑋_𝑖(𝑠)𝑋_𝑖(𝑡)]=min{𝑠,𝑡} for 𝑖=1,2, and 𝖤[𝑋₁(𝑠)𝑋₂(𝑡)]=𝜌min{𝑠,𝑡}. For detailed properties of such process, see, for example, [3]. A quantity that of general interest is, for 𝑧≥𝑥 and 𝑤≥𝑦 and any fixed time 𝑡₀>0𝖯𝑋₁𝑡₀∈𝑑𝑥,𝑀₁𝑡₀∈𝑑𝑧;𝑋₂𝑡₀∈𝑑𝑦,𝑀₂𝑡₀∈𝑑𝑤, (2.1) where 𝑀_𝑖(𝑡) are the running maximum for 𝑋_𝑖(𝑡) for 𝑖=1,2 respectively, that is 𝑀_𝑖(𝑡)=sup_{0≤𝑠≤𝑡}{𝑋_𝑖(𝑠)}. It is easy to see that we can just focus on the case of 𝑡₀=1, general case can be easily obtained through scaling. To be more specific, we have, for any 𝑡₀, 𝖯𝑋₁𝑡₀∈𝑑𝑥,𝑀₁𝑡₀∈𝑑𝑧;𝑋₂𝑡₀∈𝑑𝑦,𝑀₂𝑡₀𝑋∈𝑑𝑤=𝖯₁(1)∈𝑑𝑥𝑡₀,𝑀₁(1)∈𝑑𝑧𝑡₀;𝑋₂(1)∈𝑑𝑦𝑡₀,𝑀₂(1)∈𝑑𝑤𝑡₀, (2.2) 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 0, conditioning on the event that it ends at 𝑥 at time 1, has the same law as 𝑊(𝑡)+𝑥𝑡, (3.1) where 𝑊(𝑡) is a standard Brownian bridge, that is, a Brownian motion satisfies 𝑊(0)=𝑊(1)=0. Therefore, from the representation of Brownian bridge, we have, for 𝑡∈[0,1], 𝑊(𝑡)+𝑥𝑡^d1=𝑡𝐵𝑡−1+𝑥𝑡, (3.2) where 𝐵(𝑡) is a standard Brownian motion, =^d means equal in distribution. Thus, for any 𝑧≥𝑥, the event of 1𝑡𝐵𝑡−1+𝑥𝑡≤𝑧,∀𝑡∈[0,1], (3.3) will have the same probability as, 𝐵1𝑡≤𝑧−1𝑡−𝑥,∀𝑡∈[0,1]. (3.4) Apply a variable transform 𝑠=1/𝑡−1, we get, {𝐵(𝑠)≤𝑧𝑠+(𝑧−𝑥),∀𝑠≥0}, (3.5) or equivalently, {𝐵(𝑠)−𝑧𝑠≤𝑧−𝑥,∀𝑠≥0}. (3.6) This probability of the above event is equivalent to the stationary distribution of 𝐵(𝑠)−𝑧𝑠 reflected at zero, which, in turn, equals 1−exp[−2𝑧(𝑧−𝑥)]. (3.7) Take derivative with respect to 𝑧, we have, 2(2𝑧−𝑥)exp[−2𝑧(𝑧−𝑥)]. (3.8) So, 𝖯[𝑋(1)(1)][𝑀(1)(1)]=2()∈𝑑𝑧,𝑀∈𝑑𝑎=𝖯∈𝑑𝑎∣𝑋∈𝑑𝑧2𝑧−𝑥√−()2𝜋exp2𝑧−𝑥²2. (3.9) 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 𝖯𝑋₁(1)∈𝑑𝑥,𝑀₁(1)≤𝑧,𝑋₂(1)∈𝑑𝑦,𝑀₂(1)≤𝑤, (3.10) for 𝑧≥𝑥 and 𝑤≥𝑦 can be reduced to calculating 𝖯𝑌₁(𝑡)−𝑧𝑡≤𝑧−𝑥,𝑌₂(𝑡)−𝑤𝑡≤𝑤−𝑦,∀𝑡≥0. (3.11) 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 𝖯𝑌₁(𝑡)≤𝑧−𝑥,𝑌₂(𝑡)≤𝑤−𝑦,∀𝑡≥0, (3.12) where (𝑌₁(𝑡),𝑌₂(𝑡)) is a two-dimensional Brownian motion with negative drifts that reflected at the two positive axes. To be more precise, 𝑌₁(𝑡)=−𝑧𝑡+𝑊₁(𝑡)+𝐿₁(𝑡),𝑌₂(𝑡)=−𝑤𝑡+𝑊₂(𝑡)+𝐿₂(𝑡), (3.13) where (𝑊₁,𝑊₂) is a two-dimensional driftless Brownian motion with covariance matrix, Γ=1𝜌𝜌1, (3.14) and 𝐿_𝑖(𝑡),𝑖=1,2 are the local time defined as 𝐿_𝑖(𝑡)≜sup_{0≤𝑠≤𝑡}−𝜃_𝑖𝑡+𝑊_𝑖(𝑡)⁻, (3.15) 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:ℒ𝑓(_ℝ2+𝑧)𝜋𝑥,𝑦𝑑𝑥𝑑𝑦+²_{𝑗=1𝐹𝑗}𝐷_𝑗𝑓(𝑧)𝜎_𝑖(𝑧)𝑑𝑣_𝑗(𝑧)=0, (4.1) whereℒ=𝐷^𝑡Γ𝐷+𝜃₁,𝜃₂^𝑡𝐷 (4.2) 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: 𝜕²𝑢𝜕𝑥²𝜕+2𝜌²𝑢+𝜕𝜕𝑥𝜕𝑦²𝑢𝜕𝑦²+𝜃₁𝜕𝑢𝜕𝑥+𝜃₂𝜕𝑢𝜕𝑦=𝑒^{−𝑧𝑥−𝑤𝑦}. (4.3) Let us first consider the homogeneous version of the equation, 𝜕²𝑢𝜕𝑥²𝜕+2𝜌²𝑢+𝜕𝜕𝑥𝜕𝑦²𝑢𝜕𝑦²+𝜃₁𝜕𝑢𝜕𝑥+𝜃₂𝜕𝑢𝜕𝑦=0. (4.4) Apply transform, 𝜌𝜉=𝑥,𝜂=−1−𝜌²1𝑥+1−𝜌²𝑦, (4.5) then, 𝑢(𝜉,𝜂) satisfies Δ_𝜉𝜂𝑢+𝜃₁𝑢_𝜉+𝜃₂𝑢_𝜂=0, (4.6) where 𝜃₁=𝜃₁ and 𝜃₂=(𝜃₁−𝜌𝜃₂)/1−𝜌². Let 𝑣(𝜉,𝜂)=𝑒𝜃₁^𝜉+𝜃₂^𝜂𝑢(𝜉,𝜂), then 𝑣(𝜉,𝜂) is the solution of the following problem, Δ_𝜉𝜂𝑣−𝜃₁𝜃₂𝑣=0. (4.7) 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 𝑣𝜉,𝜂=𝑐₁𝑒^{−√𝑘𝜉}+𝑐₂𝑒^√𝑘𝜉𝑐₃𝑒⁻𝜃₁𝜃_2−𝑘𝜂+𝑐₄𝑒𝜃₁𝜃_2−𝑘𝜂, (4.8) for any constant 𝑐₁,𝑐₂,𝑐₃,𝑐₄, and 𝑘≤𝜃²₁/4+𝜃²₂/4. Therefore, 𝑢𝑥,𝑦=𝑐₁𝑐₃exp−√𝑘+𝜃₁+𝜌𝜃₁𝜃₂−𝑘+𝜃₂1−𝜌²𝑥−𝜃₁𝜃₂−𝑘+𝜃₂1−𝜌²𝑦+𝑐₁𝑐₄exp−√𝑘+𝜃₁−𝜌𝜃₁𝜃₂−𝑘+𝜃₂1−𝜌²𝑥+𝜃₁𝜃₂−𝑘−𝜃₂1−𝜌²𝑦+𝑐₂𝑐₃exp√𝑘−𝜃₁+𝜌𝜃₁𝜃₂−𝑘+𝜃₂1−𝜌²𝑥−𝜃₁𝜃₂−𝑘+𝜃₂1−𝜌²𝑦+𝑐₂𝑐₄exp√𝑘−𝜃₁−𝜌𝜃₁𝜃₂−𝑘−𝜃₂1−𝜌²𝑥+𝜃₁𝜃₂−𝑘−𝜃₂1−𝜌²𝑦. (4.9) For the inhomogeneous version, it is easy to see that 𝑢−(1/(𝑧²+2𝜌𝑧𝑤+𝑤²−𝜃𝑧−𝜃𝑤))𝑒^{−𝑧𝑥−𝑤𝑦} will be the solution.

Now plug 𝑢(𝑥,𝑦) into the BAR (4.1), we have𝜋(𝑧,𝑤)+𝑐₁𝑐₃−√𝑘+𝜃₁+𝜌𝜃₁𝜃₂−𝑘+𝜃₂1−𝜌²𝜎𝜃₁𝜃₂−𝑘+𝜃₂/21−𝜌²−𝜃₁𝜃₂−𝑘+𝜃₂1−𝜌²𝜎−√𝑘+𝜃₁+𝜌𝜃₁𝜃₂−𝑘+𝜃₂1−𝜌²+𝑐₁𝑐₄−√𝑘+𝜃₁−𝜌𝜃₁𝜃₂−𝑘−𝜃₂1−𝜌²𝜎𝜃₁𝜃₂−𝑘−𝜃₂1−𝜌²−𝜃₁𝜃₂−𝑘−𝜃₂1−𝜌²𝜎−√𝑘+𝜃₁−𝜌𝜃₁𝜃₂−𝑘−𝜃₂1−𝜌²+𝑐₂𝑐₃√𝑘−𝜃₁+𝜌𝜃₁𝜃₂−𝑘+𝜃₂1−𝜌²𝜎𝜃₁𝜃₂−𝑘+𝜃₂1−𝜌²−𝜃₁𝜃₂−𝑘+𝜃₂1−𝜌²𝜎√𝑘−𝜃₁+𝜌𝜃₁𝜃₂−𝑘+𝜃₂1−𝜌²+𝑐₂𝑐₄√𝑘−𝜃₁−𝜌𝜃₁𝜃₂−𝑘−𝜃₂1−𝜌²𝜎𝜃₁𝜃₂−𝑘−𝜃₂1−𝜌²−𝜃₁𝜃₂−𝑘−𝜃₂1−𝜌²𝜎√𝑘−𝜃₁−𝜌𝜃₁𝜃₂−𝑘−𝜃₂1−𝜌²−𝜎(𝑧)−𝜎(𝑤)𝑧²+2𝜌𝑧𝑤+𝑤²−𝜃₁𝑧−𝜃₂𝑤=0, (4.10) where 𝜋 and 𝜎 denote the Laplace transform of 𝜋(𝑥,𝑦) and 𝜎(𝑥). Given that 𝑐_𝑖, 𝑖=1,…,4 and 𝑘 can be arbitrary, it is necessary that −√𝑘+𝜃₁+𝜌𝜃₁𝜃₂−𝑘+𝜃₂1−𝜌²𝜎𝜃₁𝜃₂−𝑘+𝜃₂1−𝜌²−𝜃₁𝜃₂−𝑘+𝜃₂1−𝜌²𝜎−√𝑘+𝜃₁+𝜌𝜃₁𝜃₂−𝑘+𝜃₂1−𝜌²−√=0,𝑘+𝜃₁−𝜌𝜃₁𝜃₂−𝑘−𝜃₂1−𝜌²𝜎𝜃₁𝜃₂−𝑘−𝜃₂1−𝜌²−𝜃₁𝜃₂−𝑘−𝜃₂1−𝜌²𝜎−√𝑘+𝜃₁−𝜌𝜃₁𝜃₂−𝑘−𝜃₂1−𝜌²√=0,𝑘−𝜃₁+𝜌𝜃₁𝜃₂−𝑘+𝜃₂1−𝜌²𝜎𝜃₁𝜃₂−𝑘+𝜃₂1−𝜌²−𝜃₁𝜃₂−𝑘+𝜃₂1−𝜌²𝜎√𝑘−𝜃₁+𝜌𝜃₁𝜃₂−𝑘+𝜃₂1−𝜌²√=0,𝑘−𝜃₁−𝜌𝜃₁𝜃₂−𝑘−𝜃₂1−𝜌²𝜎𝜃₁𝜃₂−𝑘−𝜃₂1−𝜌²−𝜃₁𝜃₂−𝑘−𝜃₂1−𝜌²𝜎√𝑘−𝜃₁−𝜌𝜃₁𝜃₂−𝑘−𝜃₂1−𝜌²=0. (4.11) In fact, they are all in the form of 𝐴𝜎(𝐵)−𝐵𝜎(𝐴)=0, (4.12) and 𝐴 and 𝐵 satisfy the quadratic equation 𝐴+𝜌𝐵+𝜃₁²+1−𝜌²𝐵+𝜃₂²=𝜃²₁+𝜃²₂4, (4.13) which is equivalent to 𝐴²+2𝜌𝐴𝐵+𝐵²+2𝜃₁𝐴+2𝜃₂3𝜃𝐵+²¹+𝜃²₂−2𝜌𝜃₁𝜃₂41−𝜌²=0. (4.14) For any 𝑥∈ℝ, denote 𝜂(𝑥) as the positive conjugate point of 𝑥 with respect to the curve defined by (4.14). Therefore, 𝜂(𝑥)=4𝜌²−1𝑥²+8𝜌𝜃₂−8𝜃₁𝑥+4𝜃²₂+3𝜃²₁+𝜃²₂−2𝜌𝜃₁𝜃₂/41−𝜌²−2𝜃₂+2𝜌𝑥2. (4.15)

Assume that 𝜎(𝑥)=𝜂(𝑥)⁻¹. It is easy to see that the equation has positive solution only when 𝜃_𝑖<0. Hence, we have 𝜎(𝑧)=24𝜌²−1𝑧²+8𝜌𝜃₂−8𝜃₁𝑧+4𝜃²₂+3𝜃²₁+𝜃²₂−2𝜌𝜃₁𝜃₂/41−𝜌²−2𝜃₂.+2𝜌𝑧 (4.16) Therefore, from (4.10), we obtain the following.

Theorem 4.3. The Laplace transform of the stationary 𝜋(𝑥,𝑦) bears the following form:𝜋()=2𝑧,𝑤𝑧²+2𝜌𝑧𝑤+𝑤²−𝜃₁𝑧−𝜃₂𝑤×𝑤4𝜌²−1𝑤²+8𝜌𝜃₂−8𝜃₁𝑤+4𝜃²₂+3𝜃²₁+𝜃²₂−2𝜌𝜃₁𝜃₂/41−𝜌²−2𝜃₂−𝑧+2𝜌𝑤4𝜌²−1𝑧²+8𝜌𝜃₂−8𝜃₁𝑧+4𝜃²₂+3𝜃²₁+𝜃²₂−2𝜌𝜃₁𝜃₂/41−𝜌²−2𝜃₂.+2𝜌𝑧 (4.17)

5. Final Calculations

Recall that our goal is to get an expression for 𝖯𝑋₁(𝑡)∈𝑑𝑥,𝑀₁(𝑡)∈𝑑𝑧;𝑋₂(𝑡)∈𝑑𝑦,𝑀₂(𝑡)∈𝑑𝑤, (5.1) which is, of course, 𝜕²𝖯𝑋𝜕𝑧𝜕𝑤¹(𝑡)∈𝑑𝑥,𝑀₁(𝑡)≤𝑑𝑧;𝑋₂(𝑡)∈𝑑𝑦,𝑀₂(𝑡)≤𝑑𝑤. (5.2) The results in the previous section enable us to obtain the Laplace transform of 𝖯𝑋₁(𝑡)∈𝑑𝑥,𝑀₁(𝑡)≤𝑑𝑧;𝑋₂(𝑡)∈𝑑𝑦,𝑀₂(𝑡)≤𝑑𝑤. (5.3) Apply the Laplace transform with respect to 𝑥,𝑦, we get ^∞_−∞𝑒^{𝑠𝑥+𝑡𝑦}𝖯𝑋₁(1)∈𝑑𝑥,𝑀₁(1)≤𝑧;𝑋₂(1)∈𝑑𝑦,𝑀₂(1)=≤𝑤_∞−∞𝑒^{𝑠𝑥+𝑡𝑦}𝖯𝑀₁(1)≤𝑧,𝑀₂(1)≤𝑤∣𝑋₁(1)∈𝑑𝑥,𝑋₂(1)𝑋∈𝑑𝑦×𝖯₁(1)∈𝑑𝑥,𝑋₂(1)=∈𝑑𝑦_∞−∞𝑒^{𝑠𝑥+𝑡𝑦}𝖯𝑌₁≤𝑧−𝑥,𝑌₂(1)=≤𝑤−𝑦_∞−∞𝑒^{𝑠𝑥+𝑡𝑦𝑧−𝑥}_−∞^𝑤−𝑦_−∞𝜋_𝑧,𝑤(𝑢,𝑣)𝑑𝑢𝑑𝑣𝖯𝑋₁(1)∈𝑑𝑥,𝑋₂(1)=1∈𝑑𝑦𝑑𝑥𝑑𝑦𝑒𝑠𝑡^{−𝑠𝑧+𝑡𝑤}𝜋(−𝑠,−𝑡;𝑧,𝑤)∗𝜙(𝑠,𝑡), (5.4) where 𝜙(𝑠,𝑡) is the Laplace transform of 𝖯[𝑋₁(1)∈𝑑𝑥,𝑋₂(1)∈𝑑𝑦]. Hence, we obtain the following.

Theorem 5.1. The laplace transform of𝖯𝑋₁(1)∈𝑑𝑥,𝑀₁(1)≤𝑑𝑧;𝑋₂(1)∈𝑑𝑦,𝑀₂(1)≤𝑑𝑤 (5.5) with respect (𝑥,𝑦) has the form, 1𝑒𝑠𝑡^{−𝑠𝑧+𝑡𝑤}𝜋(−𝑠,−𝑡;𝑧,𝑤)∗𝜙(𝑠,𝑡). (5.6)

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.