Dynamic analysis of the M/G/1 queueing system with multiple phases of operation

1.   Introduction

Queueing systems with multiple phases of operation are very useful in manufacturing systems, transportation systems, financial systems, etc.; see ^[1,2,3,4]. For example, in automobile manufacturing, raw materials first enter the stamping workshop for stamping and forming. Continuing with the welding assembly in the welding workshop is the second stage, then carry out surface treatment such as painting, is the third stage. Finally, the final assembly and quality inspection are carried out. In the airport boarding process, passengers first queue up at the check-in counter to complete check-in procedures, including checking in luggage. Then, queue up through the security checkpoint for security checks, is the second stage. Finally, they queue up at the boarding gate to wait for boarding. In hospital medical systems, patients first queue up at the registration counter to register. Then, the second stage is to go to the waiting area of the corresponding department to queue up for the doctor's diagnosis. If further examinations are needed, such as blood tests, X-rays, etc., one needs to queue outside the examination department to wait for the examination, which is the third stage. Finally, they take the examination results and return to the doctor for diagnosis or treatment, such as prescribing medication, intravenous infusion, etc.

In this paper, we consider the M/G/1 queueing system with multiple phases of operation, where M means that customer arrival follows a Poisson process, G represents the service rate of the server follows a general distribution, 1 indicates the number of servers in the system. Therefore, the M/G/1 queue is a single-server queue where customers arrive according to a Poisson process, and service times are independent and identically distributed general random variables. The process of establishing a mathematical model for the queuing system is as follows:

There are $n+1$ phases in this system, $0$ is the idle period, and $s\; (s = 1, 2, \cdots, n)$ are the operational phases. We assume that $N(t)$ denotes the number of customers in the system at time $t$, $J(t)$ denotes the phase in which the system operates at time $t$, $\; R_{s}(t)\; (s = 1, 2, \cdots, n)$ represents the elapsed service time of the customer currently receiving service during phase $s$, $\; F_s(x) = Prob\{R_s(t)\le x\}\; (s = 1, 2, \cdots, n)$ denotes the probability distribution function corresponding to $R_s(t)$, and $\mu_s(x)dx$ is the service completion rate of the server in the interval $(x, x+dx]$ if the system is in phase $s$ and satisfies $\mu_s(x)\ge 0$ and $\int^{\infty}_0\mu_s(x)dx = \infty$. The service time of the any two phases of service are mutually independent. Based on the properties of conditional probability and differentiation, we have

Then, from this and $F_s(0) = 0$, we obtain the probability distribution function of the service time of the server in the phase $s$,

According to the definition of probability distribution function, we know that $F_s(x)\ge 0$ and $\lim_{x\to\infty}F_s(x) = 1$. Therefore,

We assume that, in phase $s$, the arrivals occur according to a Poisson process of rate $\lambda_s > 0$. In the idle period, the Poisson arrival rate is $\lambda_0 > 0$; upon arrival, the system moves to some operative phase $s$ with probability $q_s$, and the customer service upon arrival begins immediately, where $q_s > 0$ and $\sum^n_{s = 1}q_s = 1$. Then, according to the definition of Poisson process and exponential distribution, we have

where $o(\Delta t)$ denotes the infinitesimal quantity of $\Delta t$. Clearly, the process $\{(N(t), J(t), R_s(t)):t\ge0\}$ is a continuous-time Markov process with state space

We define

Then, consider the changes in the system during $\Delta t$ time. Based on the formula of total probability, the properties of Markov processes, and the above Eqs (1.1a)–(1.2b), we have (for convenience, assuming $\Delta x$ is the same as $\Delta t$)

We consider the boundary conditions as follows:

Assume that the initial values

satisfy

Based on Eqs (1.2a)–(1.5) and the definition of partial differential derivatives, we obtain the following integro-partial differential equations ^[3]:

In ^[3], several static indices for the system (1.6) such as the static queue length and static sojourn time distribution of an arbitrary customer were developed under the following static hypothesis:

● $\lim_{t\to\infty}p_{0, 0}(t) = p_{0, 0}$,

● $\lim_{t\to\infty}p_{k, s}(\cdot, t) = p_{k, s}(\cdot), \; k\ge 1, \; 1\le s\le n$.

From the perspective of partial differential equations, the above hypothesis implies the following two hypotheses:

● ($\mathbf{H1}$): The system (1.6) admits a time-evolving solution.

● ($\mathbf{H2}$): The aforementioned time-evolving solution converges to its static.

In this article, we investigate the aforementioned static hypothesis, that is, ($\mathbf{H1}$) and ($\mathbf{H2}$), and consider the asymptotic behavior of the dynamic indices of the system (1.6). It is worth noting that when $n = 1$, system (1.6) becomes the classical M/G/1 queuing system ^[5], and a detailed dynamic analysis of this classical queuing system was conducted in ^[6,7,8,9,10]. Gupur et al. ^[6] studied the well-posedness of the queueing system ^[5] and obtained the strong convergence of the solution of the system when the service rate is constant. When the service rate is a bounded function, similar results to ^[6] are obtained by ^[7]. In ^[8,9,10], the exponential convergence of the solution of the queueing system ^[5] was studied. Therefore, our results include the above findings.

In the study of partial differential equations, many solutions of partial differential equations can function as semigroups in Banach space. By studying the properties of semigroups, important properties such as the existence, uniqueness, and stability of solutions to partial differential equations can be obtained; see ^{[11,12,13,14,15]}. In this article, first, we convert the system (1.6) into an abstract Cauchy problem in a natural state Banach space. Then, using the semigroup theory on Banach spaces, we show that the system operator generates a positive $C_0-$semigroup of contractions on the state space. Consequently, we verify that the system (1.6) admits a unique positive time-evolving solution, which shows that ($\mathbf{H1}$) holds under certain conditions.

Second, we investigate the asymptotic behavior of the solution. To this aim, we need to know the spectrum of the system operator; see ^[16,17,18]. For example, in order to study the asymptotic stability of multilayer thermal wave systems, Avalos et al. ^[16] conducted spectral analysis on the system and found that the system operator had no spectral points on the imaginary axis. Drogoul and Veltz ^[17] proved that 0 is the unique spectral point of the spike neural network operator on the imaginary axis, thus obtaining the exponential stability of the system. In ^[18], it was proved that 0 is the point spectrum of the system operator of the queuing system and is the unique spectral point, thus obtaining the solution corresponding to this queuing system that strongly converges to its static solution.

Moreover, it is a challenge to find spectrum on the imaginary axis. Here, we apply Greiner's ^[19] boundary perturbation method to fully describe the spectral distribution of the system operator on the imaginary axis. If the service rates are constants, then we obtain that the system operator of the system (1.6) has uncountable eigenvalues on the left-half complex plane. Consequently, these spectral results imply that the time-evolving solution of system (1.6) at most strongly converges to its static solution. In other words, ($\mathbf{H2}$) holds only in the context of strong convergence.

Finally, we discuss the asymptotic behavior of the dynamic indices of the system (1.6). By using cone theory and positive operator theory, we prove that the time-evolving queue length of the system (1.6) converges to its static queue length under some conditions. This asymptotic result includes the result of ^[3].

The remaining part of this article is organized as follows. In the next section, we rewrite the system (1.6) as an abstract Cauchy problem in a Banach space and provide the well-posedness. In Section 3, we provide a complete asymptotic behavior of the solution of system (1.6). In Section 4, we discuss the dynamic indices of the system (1.6). To demonstrate the exponential convergence of the solution, we conduct some numerical analysis in Section 5. We conclude this article in the last section.

2.   Well-posedness

We choose the state Banach space of system (1.6) as follows:

Define the maximal operator of the system (1.6) by

with domain

where $p_{k, s}$ are absolutely continuous functions, $k\ge1, \; 1\le s\le n$, and

We choose the boundary space as $\partial X = \underbrace{l^1\times l^1\times\cdots\times l^1}_{n}$ and define boundary operators $\Psi, \Phi:D(A_m)\rightarrow \partial X$ of the system (1.1a) by

Now, we introduce the system operator $(A_{\Phi}, D(A_{\Phi}))$ of the system (1.6) by

Then, the system (1.6) can be written as an abstract Cauchy problem in the Banach space $X$:

where $g_1(\cdot) = (g_{0, 0}, g_{1, 1}(\cdot), g_{2, 1}(\cdot), \cdots), g_s(\cdot) = (g_{1, s}(\cdot), g_{2, s}(\cdot), \cdots), \; s = 2, 3, \cdots, n$.

Theorem 2.1. Let $A_\Phi$ be defined by Eq (2.1). If $\mu_s(x)\le\overline{\mu}_s = \sup_{x\in[0, \infty)}\mu_s(x) < \infty\; (1\le s\le n)$, then $A_\Phi$ generates a positive $C_0-$semigroup $e^{A_\Phi t}$ of contractions on $X$.

Proof. For self contained and conciseness, we only sketch the proof of Theorem 2.1. To start, we divided operator $A_{\Phi}$ into three parts $A_{\Phi} = A+U+E$, where

where

We can verify that $\|(\gamma I-A)^{-1}\| < \frac{1}{\gamma-\mathbb{M}_0}$ for all $\gamma > \mathbb{M}_0 = \max\{\overline{\mu}_1, \overline{\mu}_2, \cdots, \overline{\mu}_n\}$, where $I$ represents the identity operator. Moreover, it is easy to prove that $D(A)$ is dense in $X$. Then, using the Helle-Yosida theorem (see [20,Theorem 4.20]), we obtain that $A$ generates a $C_0-$semigroup. Due to the operators $U$ and $E$ being linear bounded, with perturbation theory of the semigroup, we know that $A_{\Phi}$ generates a $C_0-$semigroup $e^{A_\Phi t}$. Thus, we complete the proof of this theorem. □

Next, we investigate the isometry of $e^{A_\Phi t}$. It is easy to obtain that $X^*$, the dual space of $X$, is as follows:

If we take a set $X_{+} = \{(p_1, p_2, \cdots, p_n)\in X\mid p_{0, 0}\geq0, \; p_{k, s}(\cdot)\geq 0, \; k\geq 1, \; 1\le s\le n\}$ in $X$, then, Theorem 2.1 ensures $e^{A_\Phi t}X_{+}\subset X_{+}$. Now, for any $(p_1, p_2, \cdots, p_n)\in D(A_\Phi)\cap X_{+}$, we choose

It is not difficult to calculate that for $(p^*_1, p^*_2, \cdots, p^*_n)\in X^*$ and $(p_1, p_2, \cdots, p_n)\in X_+$, we have

This shows that $(p^*_1, p^*_2, \cdots, p^*_n)\in Q(p_1, p_2, \cdots, p_n)$, where

In addition, we obtain for $(p_1, p_2, \cdots, p_n)\in D(A_\Phi)$ and $(p^*_1, p^*_2, \cdots, p^*_n)\in Q(p_1, p_2, \cdots, p_n)$ that

Then, Eq (2.3) implies that $A_\Phi$ is conservative with respect to $Q(\cdot)$. Theorem 3.6.1 of ^[21] stated that: Assume that $A_\Phi$ is densely defined, conservative with respect to $Q(\cdot):D(A_{\Phi})\to X^*$ and $(\gamma I-A_{\Phi})D(A_{\Phi}) = X$ for some $\gamma > 0$. Then, for some $g\in D((A_{\Phi})^2)$, the corresponding semigroup to Cauchy problem (2.2) is isometric. Hence, we obtain the following result.

Theorem 2.2. Let $\mu_s(\cdot)$ satisfy $\mu_s(x)\le\sup_{x\in[0, \infty)}\mu_s(x) < \infty, \; 1\le s\le n$. If the initial value $(p_1, p_2, \cdots, p_n)(\cdot, 0) = (g_{1}(\cdot), g_{2}(\cdot), \cdots, g_n(\cdot))$ of the system (2.2) belongs to $D(A_\Phi^2)$, then semigroup $e^{A_\Phi t}$ is isometric for $(g_{1}(\cdot), g_{2}(\cdot), \cdots, g_n(\cdot))$. That is,

By combining Theorems 2.1 and 2.2, we obtain the main result in this section.

Theorem 2.3. Let $\mu_s(\cdot)$ satisfy $\mu_s(x)\le\sup_{x\in[0, \infty)}\mu_s(x) < \infty, 1\le s\le n$. If the initial value $(g_{1}(\cdot), g_{2}(\cdot), \cdots, g_n(\cdot))$ of system (2.2) belongs to $D(A_\Phi^2)$, then system (2.2) admits a unique positive time-evolving solution $(p_1(\cdot, t), p_2(\cdot, t), \cdots, p_n(\cdot, t))$ which satisfies

Proof. Due to $(g_{1}(\cdot), g_{2}(\cdot), \cdots, g_n(\cdot))\in D(A_\Phi^2)$ and $(g_{1}(\cdot), g_{2}(\cdot), \cdots, g_n(\cdot))\in X_+$, it is easy to see that $(g_{1}(\cdot), g_{2}(\cdot), \cdots, g_n(\cdot))\in D(A_\Phi^2)\cap X_+$. By Theorem 1.81 of ^[12], we see that the system (1.6) has a unique positive time-evolving solution $(p_1(\cdot, t), p_2(\cdot, t), \cdots, p_n(\cdot, t))$, which can be expressed as

From this together with Eq (2.7), we obtain

This illustrates the physical meaning of $(p_1(\cdot, \cdot), p_2(\cdot, \cdot), \cdots, p_n(\cdot, \cdot))$. □

3.   Asymptotic behavior of the time-evolving solution

In this section, our main objective is to address the issue of asymptotic behavior of the time-evolving solution that we stated in Eq (2.8). In this regard, we prove that the time-evolving solution of the system (1.6) strongly converges but not exponentially converges to its static solution. In other words, hypothesis ($\mathbf{H2}$) holds only in the context of strong convergence.

3.1.   Strong convergence of the time-evolving solution

The main result of this subsection is given by the following Theorem 3.1.

Theorem 3.1. Let $\mu_s(x):[0, \infty)\to[0, \infty)$ be a measurable function that satisfies

Then, the time-evolving solution of the system (2.2) strongly converges to its static solution. In other words,

where $(p^*_1, p^*_2, \cdots, p^*_n)$ and $(p_1, p_2, \cdots, p_n)$ are the eigenvectors associated to zero, respectively.

To prove the above Theorem 3.1, we need to find the spectra of $A_{\Phi}$ along the imaginary axis. For this, we first provide the following seven lemmas.

Lemma 3.1. Let $A_{\Phi}$ be defined by Eq (2.1). If $\mu_s(x):[0, \infty)\to[0, \infty)$ is measurable function that satisfies

then, zero is an eigenvalue of $A_{\Phi}$ with geometric multiplicity one.

Proof. We need to solve $A_{\Phi}(p_1, p_2, \cdots, p_n) = 0$ for unknown $(p_1, p_2, \cdots, p_n)$. This equation is equivalent to

Solve Eqs (3.1b) and (3.1c) to obtain

If we take $p_{k, s}(0) = 2^{-(k+1)}p_{1, s}(0), \; p_{1, s}(0) = q_s \lambda_{0}p_{0, 0} > 0$ and define

then $p_{k, s}(0) = 2^{-(k+1)}p_{1, s}(0)$ satisfies the boundary conditions (3.1d) and (3.1e). Therefore, since the Cauchy product of series, the formula $\int_0^{\infty}\mu_s(x)e^{-\int_0^x \mu_s(\xi)d\xi}dx = 1$, and

we have

Eq (3.3) means that zero is an eigenvalue of $A_{\Phi}$. Moreover, by Eqs (3.1a) and (3.1d)–(3.2), we see that the geometric multiplicity of zero is one. □

Now, we use the idea of ^[19] to describe the other spectrum of $A_\Phi$ along the imaginary axis. For this objective, we define the operator $(A_0, D(A_0))$ by

and discuss the inverse of $A_0$.

For given $(y_1, y_2, \cdots, y_n)\in X$, consider $(\gamma I-A_0)(p_1, p_2, \cdots, p_n) = (y_1, y_2, \cdots, y_n)$ of unknown $(p_1, p_2, \cdots, p_n)\in D(A_0)$. This equation can be equivalently written as the following system of equations

By solving Eqs (3.4a)–(3.4c) and using Eq (3.4d), we obtain

Denoting by

then the Eqs (3.5a)–(3.5c) and $\varphi_s f(x) = \int_0^{\infty}f(x)\mu_s(x)dx, f\in L^1[0, \infty)$ give, if the resolvent of $A_0$ exists,

The following Lemma 3.2 indicates the resolvent set $\rho(A_0)$ of $A_0$.

Lemma 3.2. Let $\mu_s(x):[0, \infty)\to[0, \infty)$ be a measurable function that satisfies

Then, $\{\gamma\in \mathbb{C}\mid {\rm Re}(\gamma)+\lambda_0 > 0, {\rm Re}(\gamma)+\underline{\mu}_s > 0\}\subset\rho(A_0)$.

Proof. For all $f\in L^1[0, \infty)$, by performing integration by parts to Eq (3.6), it is easy to obtain that $E_s$ satisfies the following inequality:

Then, using the inequality $\|\varphi_s\|\leq\sup_{x\in[0, \infty)}\mu_s(x)$, we calculate for any $(y_1, y_2, \cdots, y_n)\in X$ that

That is, inequality (3.8) means that the result of this lemma is correct. □

Next, we use the following Lemma 3.3 to provide a specific expression for the Dirichlet operator.

Lemma 3.3. Let $\gamma\in\{\gamma\in \mathbb{C}\mid {\rm Re}(\gamma)+\lambda_0 > 0, {\rm Re}(\gamma)+\underline{\mu}_s > 0\}$. Then, we have $(p_1, p_2, \cdots, p_n)\in \ker(\gamma I-A_m)$ if, and only if,

Proof. If $(p_1, p_2, \cdots, p_n)\in \ker(\gamma I-A_m),$ then $(\gamma I-A_m)(p_1, p_2, \cdots, p_n) = 0$, that is,

By solving Eqs (3.10a)–(3.10c), we obtain

Since $(p_1, p_2, \cdots, p_n)\in\ker(\gamma I-A_m)$, according to the Sobolev embedding theorem ^[22], we can easily obtain

Hence, Eqs (3.11a)–(3.12) show that Eqs (3.9a)–(3.9c) are true.

On the other hand, if Eqs (3.9a)–(3.9c) hold, due to the formula

it holds true for any $c > 0$ and positive integer $k\ge 1$, and performing integration by parts, we deduce

Then, by the Cauchy product of series, we calculate that

for any $\gamma > -\underline{\mu}_s$. Inequalities (3.12) and (3.13) show that $(p_1, p_2, \cdots, p_n)\in X$. In addition, by Eq (3.9b), we have

Combining the above Eqs (3.14a) and (3.14b) with inequality (3.13), we obtain

This inequality implies that

Hence, Eqs (3.13)–(3.15) indicate that $(p_1, p_2, \cdots, p_n)\in D(A_m)$ and

□

Clearly, by the definition it is not difficult to show that the boundary operator $\Psi$ is surjective. In addition, for all $\gamma\in\rho(A_0)$, the operator

is invertible. Now, for any $\gamma\in\rho(A_0)$, we introduce the Dirichlet operator by

Then, using Lemma 3.3, for any $\gamma\in\rho(A_0)$, we can obtain the following specific expression for $D_{\gamma}$:

where

Finally, using the expression of Dirichlet operator (3.16) and the boundary operator $\Phi$, we can calculate the specific expression of $\Phi D_{\gamma}$ as follows:

where

The following Lemma 3.4 was found in ^[23], and we use this lemma along with the above results in this subsection to provide spectrum $\sigma(A_{\Phi})$ of $A_{\Phi}$ on the imaginary axis.

Lemma 3.4. Assume $\gamma\in\rho(A_0)$. If there exists $\gamma_0$ that satisfies $1\notin\sigma(\Phi D_{\gamma_0})$, then $\gamma\in\sigma(A_{\Phi})$ if, and only if, $1\in\sigma(\Phi D_{\gamma})$.

Lemma 3.5. Let $A_{\Phi}$ be defined by Eq (2.1). If $\mu_s(x):[0, \infty)\to[0, \infty)$ is a measurable function that satisfies

then, we have $i\mathbb{R}\cap \sigma(A_{\Phi}) = \{0\}, \; i^2 = -1$.

Proof. If we take $\gamma = ib, \; \; i^2 = -1, \; b\in \mathbb{R}\setminus\{0\}$, then applying the Riemann-Lebesgue lemma, we obtain that there exists $\mathbb{M} > 0$ that satisfies

for all $|b| > \mathbb{M}$. Hence, using inequality (3.18) and the formulas $\sum^n_{s = 1}q_s = 1$ and

we calculate for $p_s(0) = (p_{1, s}(0), p_{2, s}(0), p_{3, s}(0), \cdots)\in l^1\setminus\{0\}$ that

That is, $\|\Phi D_{ib}\| < 1$ for all $|b| > \mathbb{M}$. Since $\lambda_0 > 0$ and $\underline{\mu}_s > 0$, there exists $\gamma_1 = \min\{\lambda_0, \underline{\mu}_s\} > 0$ such that $\{\gamma\in\mathbb{C}\mid{\rm Re}(\gamma) > -\gamma_1\}\subset \rho(A_0)$. This means that $\gamma = ib\in\rho(A_0)$. Then, by the above inequality (3.31), we know that the spectral radius $r(\Phi D_{ib})$ of operator $\Phi D_{ib}$ satisfies $r(\Phi D_{ib})\leq\|\Phi D_{ib}\| < 1$ if $|b| > \mathbb{M}$. In other words, $1\notin\sigma(\Phi D_{ib})$ for $|b| > \mathbb{M}$. This indicates that there must be $\gamma_0 = 2|b|$ satisfying $1\notin\sigma(\Phi D_{\gamma_0})$. Consequently, using Lemma 3.4, we obtain $\gamma = ib\notin\sigma(A_{\Phi})$ for $|b| > \mathbb{M}$, i.e.,

On the other hand, since $e^{A_\Phi t}$ is a positive uniformly bounded semigroup (Theorem 2.1), using Corollary 2.3 of ^[24], we obtain that $\sigma(A_{\Phi})\cap i\mathbb{R}$ is imaginary additively cyclic, which states that $ib\in\sigma(A_{\Phi})\cap i\mathbb{R}$, and we deduce that $ibk\in\sigma(A_{\Phi})\cap i\mathbb{R}$ for every integer $k$. Therefore, combining the above discussion with the inclusion relationship $\{ib\big||b|\leq \mathbb{M}\}\subset\sigma(A_{\Phi})\cap i\mathbb{R}$ and Lemma 3.1, we have $\sigma(A_{\Phi})\cap i\mathbb{R} = \{0\}$. □

Lemma 3.6. The specific expression of the adjoint operator ${A_{\Phi}}^*$ of $A_{\Phi}$ is as follows:

with domain $D(A_{\Phi}^*) = \{(p^*_1, p^*_2, \cdots, p^*_n)\in X^*\mid \frac{dp^*_{k, s}(x)}{dx}\; \mathit{\text{existing and}}\; p^*_{k, s}(\infty) = h\}$, where $\phi_s = \frac{d}{dx}-[\lambda_s+\mu_s(x)]$ and $h$ is a positive constant that is independent of $k$ and $s$.

Proof. For $(p_1, p_2, \cdots, p_n)\in D(A_{\Phi})$ and $(p^*_1, p^*_2, \cdots, p^*_n)\in D(A_{\Phi}^*)$, using integration by parts, we calculate that

Then, from the last equation of the above Eq (3.20), we can obtain ${A_{\Phi}}^*$. □

Lemma 3.7. The zero is an eigenvalue of $A^*_{\Phi}$ with geometric multiplicity one.

Proof. We consider $A^*_{\Phi}(p^*_1, p^*_2, \cdots, p^*_n) = 0$. This equation is equivalent to

By the above equations, it is easy to investigate that

is a positive solution of Eqs (3.21a)–(3.21d). In addition, Eqs (3.21a)–(3.21d) are equivalent to

Clearly, Eqs (3.22a)–(3.22c) imply that the geometric multiplicity of zero is one. □

Proof of Theorem 3.1: Theorem 2.1 shows that semigroup $e^{A_\Phi t}$ is a uniformly bounded $C_0-$semigroup on Banach space $X$. In addition, using Lemmas 3.1, 3.5, and 3.7, we obtain that $\sigma_p(A_\Phi)\cap i \mathbb{R} = \sigma_p(A_{\Phi}^*)\cap i\mathbb{R} = \{0\}$ and $\{\gamma \in \mathbb{C} \mid \gamma = i b, b\neq0, b \in \mathbb{R}\}\subset\rho(A_\Phi)$, and zero is an eigenvalue of $A_\Phi^*$ with algebraic multiplicity one. Hence, due to Theorem 1.96 of ^[12], we obtain that the time-evolving solution of system (2.2) converges strongly to its static solution. In other words,

where $(p^*_1, p^*_2, \cdots, p^*_n)$ and $(p_1, p_2, \cdots, p_n)$ are the eigenvectors associated to zero in Lemmas 3.7 and 3.1, respectively.

3.2.   Exponential convergence of the time-evolving solution

To prove the exponential convergence of the time-evolving solution, we need to find the spectral distribution of $A_\Phi$ on the left-half complex plane. For this objective, we first provide the following Lemma 3.8.

Lemma 3.8. If $\lambda_s < \mu_s, \; 1\le s\le n$, then each point in

is an eigenvalue of $A_\Phi$ with geometric multiplicity one, in particular

Proof. For each $\gamma\in\Lambda$, we consider the equation $A_\Phi(p_1, p_2, \cdots, p_n) = \gamma (p_1, p_2, \cdots, p_n)$ of unknown $(p_1, p_2, \cdots, p_n)\in D(A_\Phi)$. This is equivalent to the following system:

Solving Eqs (3.23b) and (3.23c), we have

From this together with the formula

we obtain

We can thus combine Eqs (3.23e) and (3.25) to obtain

This yields

Clearly, the above equation is equivalent to

For any complex number $\xi_s$ and $\eta_s, \; 1\le s\le n$, if we set

then it is easy to see that $\xi_s$ and $\eta_s$ satisfy the following two equations:

From Eq (3.29), it is easy to determine that

Note that from Eq (3.28), we observe that

Then, by reusing the above equations and organizing it, we obtain

If $\xi_s = \eta_s$, then Eq (3.32) is simplified as

The inequality $|p_{k+2, s}(0)|\le k|\xi_s|^{k-1}|p_{3, s}(0)|+(k-1)|\xi_s|^{k}|p_{2, s}(0)|$ can be obtained by taking the absolute value of the above equation. Then, taking the sum of $k = 2$ to $\infty$ for this inequality, we obtain

If $\xi_s\neq\eta_s$, then Eq (3.32) can be written as

This means that

Moreover, take the $L^1[0, \infty)-$norm on both sides of Eq (3.34), and using the formula

for all ${\rm Re}(\gamma)+\mu_s > 0$, we have

Therefore, for all ${\rm Re}(\gamma)+\mu_s > 0$, by the above inequalities and Cauchy product of series, we obtain

In addition, from Eqs (3.23a), (3.25), (3.23d), and (3.26), it is easy to calculate that

Finally, by Eqs (3.30a) and (3.30b), it is easy to see that $\gamma\in\Lambda$ if, and only if, ${\rm Re}(\gamma)+\mu_s > 0$ and $|\xi_s| < 1, \; |\eta_s| < 1, \; 1\le s\le n$. Therefore, if $\gamma\in\Lambda$, then from Eqs (3.36), (3.33), (3.34), (3.37a), and (3.37d), we obtain

The Eq (3.38) shows that for any $\gamma$ in $\Lambda$ is an eigenvalue of $A_{\Phi}$. Moreover, Eqs (3.24), (3.26), (3.32), (3.37a), and (3.37d) mean that the geometric multiplicity of every $\gamma\in\Lambda$ is one.

Next, we observe the case $\gamma\in\mathbb{R}$. Since Theorem 2.3 implies that $(0, \infty)\subset\rho(A_\Phi)$, the real spectrum of $A_\Phi$ in the interval $(-\infty, 0]$ exists. We discuss the real spectrum of $A_\Phi$ in the following three cases.

Case 1: $\; (\gamma+\lambda_s+\mu_s)^2 > 4\lambda_s\mu_s$ if, and only if, $|\gamma+\lambda_s+\mu_s| > 2\sqrt{\lambda_s\mu_s}$. Since $\gamma+\mu_s > 0$ and $\gamma+\lambda_s+\mu_s > 2\sqrt{\lambda_s\mu_s}$, we have $\gamma > 2\sqrt{\lambda_s\mu_s}-\lambda_s-\mu_s$. From this together with $\lambda_s < \mu_s$ and $\gamma+\mu_s > 0$, it is easy to calculate that

This implies $\left(\max_{1\le s\le n}\{2\sqrt{\lambda_s\mu_s}-\lambda_s-\mu_s\}, 0\right) \subset\sigma(A_\Phi)$. Then, from this together with Lemma 3.1, we obtain

Case 2: $(\gamma+\lambda_s+\mu_s)^2 = 4\lambda_s\mu_s$ if, and only if, $|\gamma+\lambda_s+\mu_s| = 2\sqrt{\lambda_s\mu_s}$. Since $\gamma+\mu_s > 0$ and $\gamma+\lambda_s+\mu_s = 2\sqrt{\lambda_s\mu_s}$, we deduce $\gamma = 2\sqrt{\lambda_s\mu_s}-\lambda_s-\mu_s$. Then, using $\lambda_s < \mu_s$ and $\gamma+\mu_s > 0$, we have

This shows that $\max\{2\sqrt{\lambda_1\mu_1}-\lambda_1-\mu_1, \cdots, 2\sqrt{\lambda_n\mu_n}-\lambda_n-\mu_n\}$ is an eigenvalue of $A_\Phi$.

Case 3: $(\gamma+\lambda_s+\mu_s)^2 < 4\lambda_s\mu_s$ if, and only if, $-2\sqrt{\lambda_s\mu_s} < \gamma+\lambda_s+\mu_s < 2\sqrt{\lambda_s\mu_s}$. Since $\gamma+\mu_s > 0$ and $0 < \gamma+\lambda_s+\mu_s < 2\sqrt{\lambda_s\mu_s}$, we obtain $\gamma < 2\sqrt{\lambda_s\mu_s}-\lambda_s-\mu_s$. Then, from this together with $\lambda_s < \mu_s$, $\gamma+\mu_s > 0$, and $i^2 = -1$, we have

Therefore,

Hence, this implies that

Consequently, by summing up the above three cases, we obtain

□

Let $\omega_0(A_\Phi)$, $\omega_{ess}(A_\Phi)$, $s(A_\Phi)$ represent the growth bound, the essential growth bound, and spectral bound of $A_\Phi$, respectively. The spectral mapping theorem ^[11] means that

Hence, from this property and Lemma 3.8, we obtain that $e^{A_\Phi t}$ has uncountable eigenvalues. Therefore, $e^{A_\Phi t}$ is not compact and it is not eventually compact by Corollary V.3.2 of ^[11].

Additionally, due to $e^{A_\Phi t}$ being a $C_0-$semigroup on $X$ with generator $A_\Phi$, using Corollary IV.2.11 of ^[11], we know that $\omega_0 = \max\{\omega_{ess}, s(A_\Phi)\}$ and $\sigma(A_\Phi)\cap\{\gamma\in\mathbb{C}\; |\; {\rm Re}(\gamma)\ge w\}$ is finite for every $w > \omega_{ess}$. Using Lemma 3.8, we can obtain that the spectrum determined condition $\omega_0 = s(A_\Phi)$ holds and $\omega_0 = s(A_\Phi) = 0$ (we suggest that readers refer to the proof of Theorem 4.1 in ^[25] for similar proofs in this part). Hence, using the aforementioned discussions, we have $\omega_{ess} = 0$. Then, by Proposition 3.5 of ^[11], we derive that $e^{A_\Phi t}$ is not quasi-compact.

The main result of this subsection is given by the following Theorem 3.2.

Theorem 3.2. Let $\mu_s(\cdot): = \mu_s$ be a constant and $\lambda_s < \mu_s, \; 1\le s\le n$. Then, the time-evolving solution of the system (2.2) cannot exponentially converge to its static solution. That is to say, there are no constants $\mathcal M > 0$ and $\varepsilon > 0$ such that

for any $t\ge 0$ and $(p_1, p_2, \cdots, p_n)\in D(A_{\Phi})$, where $(p_1, p_2, \cdots, p_n)_{(0)}$ is the eigenvector associated to zero.

Proof. Assume that $(p_1, p_2, \cdots, p_n)_{(0)}$ and $(p_1, p_2, \cdots, p_n)_{(r)}$ are the eigenvectors of 0 and

in Lemma 3.8, respectively, for any $r\in (0, 1)$. Hence, using $A_{\Phi}(p_1, p_2, \cdots, p_n)_{(0)} = 0$ and $A_{\Phi}(p_1, p_2, \cdots, p_n)_{(r)} = r\beta_s (p_1, p_2, \cdots, p_n)_{(r)}$, we have

Therefore,

That is, there are no constants $\mathcal M > 0$ and $\varepsilon > 0$ such that

for any $t\ge 0$ and $(p_1, p_2, \cdots, p_n)\in D(A_{\Phi})$. □

In neural network ^[17] and reliability model ^[12], it has been proven that the semigroup corresponding to these systems is a quasi-compact strongly continuous semigroup, thus they obtain the dynamic solution of the corresponding system that strongly converges to its steady-state solution. Therefore, the result of Theorem 3.2 is significantly different from those in ^[12,17].

4.   Asymptotic behavior of the time-evolving queue length

Define the time-evolving queue length of system (1.6) by

Then, by combining Theorems 2.1, 2.3 and 3.1 and Lemma 3.1, we can obtain the asymptotic behavior of $L(t)$.

Theorem 4.1. Let $\mu_s(x):[0, \infty)\to[0, \infty)$ be a measurable function that satisfies

If the initial value $(g_1(\cdot), g_2(\cdot), \cdots, g_n(\cdot))$ of the system (1.6) and the eigenvector $(\tilde{p}_1(\cdot), \tilde{p}_2(\cdot), \cdots, \tilde{p}_n(\cdot))$ corresponding to zero satisfying $\tilde{p}_s(\cdot)\ge u_s(\cdot)$, then time-evolving queue length $L(\cdot)$ of system (1.6) converges to its static queue length. That is to say,

Proof. For any $(p_1, p_2, \cdots, p_n)$ and $(y_1, y_2, \cdots, y_n)$ in $X$, we introduce an order relation "$\ge$" by

Then, it is not difficult to show that "$\ge$" is a partial order relation in $X$. Therefore, $(X, \ge)$ is a poset. Let $(\tilde{p}_1, \tilde{p}_2, \cdots, \tilde{p}_n)$ be the positive eigenvector associated to zero of $A_{\Phi}$ (Lemma 3.1). Let $(\tilde{p}_1, \tilde{p}_2, \cdots, \tilde{p}_n)$ and the initial value $(g_1, g_2, \cdots, g_n)$ of the system (2.2) satisfy the aforementioned partial order relation

That is, $(\tilde{p}_1, \tilde{p}_2, \cdots, \tilde{p}_n)\ge (g_1, g_2, \cdots, g_n)$. Due to $e^{A_{\Phi}t}$ being a positive linear operator (Theorem 2.1), it is a monotone increasing operator. In addition, by Theorem 2.3 and Lemma 3.1, we know that

Therefore, from this together with the partial order relation $(\tilde{p}_1, \tilde{p}_2, \cdots, \tilde{p}_n)\ge (g_1, g_2, \cdots, g_n)$, we have

Theorem 3.1 includes the following result:

Hence, by Lemma 3.1 and the Lebesgue theorem, we obtain

This inequality shows that

□

Remark 4.1. In Theorem 4.1, the static queue length is obtained by using the eigenvector that related to zero (Lemma 3.1). This is the same as the result obtained by introducing probability generating functions in ^[3].

Similarly, we can obtain the other time-evolving indicators of system (1.6) such as the convergence of the time-evolving average number of customers to its static average number of customers.

5.   Numerical analysis

To prove the correctness of the exponential convergence results in this article, we perform numerical analysis on the spectral results in Lemma 3.8. The numerical analysis results are shown in Figures 1 and 2. We obtain these numerical results using Matlab.

In Figure 1(a), we consider that system (1.6) only has the idle period and one operational phase, that is, $n = 1$ in system (1.6). In other words, we consider the classical queuing model ^[5], and the point spectrum results of the system operator $A_{\Phi}$ of this queueing model ^[5] are obtained in detail in ^[8,9,10]. If we take $\lambda_1 = 0.1, \; \mu_1 = 0.9$, then it is easy to see that $\frac{\lambda_1}{\mu_1} = \frac{0.1}{0.9} < 1$ and $2\sqrt{\lambda_1\mu_1}-\lambda_1-\mu_1 = -0.4$. In addition, Figure 1(a) means that $0 < \xi_1 < 1$ and $0 < \eta_1 < 1$ for all $\gamma\in(-0.4, 0)$. Hence, by Eqs (3.34) and (3.36), we see that every $\gamma\in(-0.4, 0)$ is the point spectrum of $A_{\Phi}$.

In Figure 1(b), we consider that system (1.6) only has the idle period and two operational phases, that is, $n = 2$ in system (1.6). In this case, we take $\lambda_1 = 0.1, \; \mu_1 = 0.9$ and $\lambda_2 = 0.2, \; \mu_2 = 0.8$. Then, these values satisfy $\frac{\lambda_1}{\mu_1} = \frac{1}{9}, \; \frac{\lambda_2}{\mu_2} = \frac{1}{4}$, and $\max_{s = 1, 2}\{2\sqrt{\lambda_s\mu_s}-\lambda_s-\mu_s\} = -0.2$. Moreover, when $\gamma\in(-0.2, 0)$, from Figure 1(b) we see that $\xi_1, \xi_2, \eta_1$, and $\eta_2$ satisfy $0 < \xi_n < 1$ and $0 < \eta_n < 1$. Therefore, by Eqs (3.34) and (3.36), we know that all $\gamma\in(-0.2, 0)$ are the point spectrum of $A_{\Phi}$. Figure 1 means that the point spectrum results in Lemma 3.8 are correct if $\lambda_s < \mu_s, s = 1, 2$. Therefore, the exponential convergence result of Theorem 3.2 is valid.

In the following, we check whether the condition $\lambda_s < \mu_s, \; 1\le s\le n$ in Lemma 3.8 is necessary. In Figure 2, we consider that system (1.6) only has the idle period and two operational phases, that is, $n = 2$ in system (1.6). If we take $\lambda_1 = 0.1, \; \mu_1 = 0.9$ and $\lambda_2 = 0.8, \; \mu_2 = 0.2$, then we have $\lambda_1 < \mu_1$ and $\lambda_2 > \mu_2$. In this case, Figure 2(a) means that $0 < \xi_1 < 1, \; 0 < \eta_1 < 1$, and $\xi_2 > 1, \; \eta_2 > 1$ for all $\gamma\in(-0.2, 0)$.

If we take $\lambda_1 = 0.9, \; \mu_1 = 0.1$ and $\lambda_2 = 0.8, \; \mu_2 = 0.2$, then we have $\lambda_n > \mu_n$. In this case, Figure 2(b) implies that $\xi_n > 1$ and $\eta_n > 1$ for all $\gamma\in(-0.2, 0)$. Of course, we can choose some different $\lambda_n$ and $\mu_n$, at least one of which satisfies $\lambda_{s_0} > \mu_{s_0}$ for some $s_0 = 1, 2, \cdots, n$, to obtain the same conclusion. As a result, Figure 2 means that we cannot obtain whether it is $\gamma\in\sigma_p(A_{\Phi})$, or even whether it is $\gamma\in\sigma(A_{\Phi})$ under the above circumstances, where $\gamma\in(\max_{1\le s\le n}\{2\sqrt{\lambda_s\mu_s}-\lambda_s-\mu_s\}, 0)$. Therefore, in Lemma 3.8 (or in Theorem 3.2), we must consider the condition $\lambda_s < \mu_s$. This condition is also the stability condition obtained for system (1.6) in reference ^[3].

6.   Conclusions

In this article, we conduct a dynamic analysis of the M/G/1 queueing system with multiple phases of operation. Using operator semigroup theory, we prove that there exists a unique time-evolving solution for this system. We obtain the spectral distribution of the system operator on the imaginary axis and prove that the system operator has an infinite number of eigenvalues on the left-half of the complex plane. As a result, the above solution converges at most strongly to its static solution. We also discuss the compactness of the system's corresponding semigroup by using these spectral results. However, we have not obtained the complete spectrum of the system operator on the left-half of the complex plane. This is the work we will continue to do in the future. Additionally, we obtain that the dynamic queue length of the model strongly converges to its static queue length.

The method described in this article can only be applied to queuing systems established using the supplementary variable method and described by partial differential equations. For example, we cannot use the method proposed in this paper for the queuing systems in ^[1,26,27].

Use of AI tools declaration

The authors declare that they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgements

We are grateful to the anonymous referees, who read carefully the manuscript and made valuable comments and suggestions. This work was supported by the National Natural Science Foundation of China (No: 12301150) and Natural Science Foundation of Xinjiang Uygur Autonomous Region (No: 2024D01C229).

Conflicts of Interest

The authors declare there is no conflict of interest.