Analysis of single-cell RNA-sequencing data identifies a hypoxic tumor subpopulation associated with poor prognosis in triple-negative breast cancer

1.   Introduction

Semi-open queueing networks (SOQNs), also referred to as open queues with restricted capacity, have garnered the attention of the stochastic modeling research community in recent years; see the survey ^[1] and the paper ^[2]. SOQNs are characterized by users, who would like to receive service, arriving from an infinite population of sources and permanently departing from the network after service, as in open queueing networks. However, only a restricted number of users can receive service in the network at the same time, as in closed networks.

A SOQN can be organized into a complex consisting of a set of orders (permissions, windows, threads, tokens, operators, robots, automated guided vehicles, etc.) required to begin service and a core (inner) network, whose nodes are dedicated to providing service for an admitted user. If an arriving user finds an available order in the set, the user enters the core network and sequentially receives service in its nodes. After the end of service in the required sequence of nodes, the user leaves the core network and returns the order, that was received at the entrance to the set of orders. If the set is empty at a user arrival instant, then the user joins a buffer and waits for an order release (this kind of service organization is called back-ordering SOQN, see ^[3]) or is lost (called user loss SOQN).

As stated in ^[1], the utilization of SOQNs originates from the study ^[4] conducted by Avi-Itzhak and Heyman, who examined the efficacy of a multi-programming computer system. In that system, multiple tasks can be executed concurrently, with the central processor managing one task while peripheral devices attend to other tasks.

Brief overviews of the research on SOQN analysis may be found, for example, in ^[1,5,6,7] and more recent publications ^{[3,8,9,10,11,12,13]} and references therein. Examples of possible applications of SOQNs are listed in ^[1], as follows: manufacturing material control rules commonly implemented in manufacturing shops to manage the inventory of both finished goods and work-in-process; modeling the operation of vehicle rental providers, where users borrow vehicles from rental depots and return them after the rental time; modeling systems for storage and retrieval based on autonomous vehicles; modeling communication networks based on window flow control; and modeling health-care systems to accurately depict the dynamics of patient-resource waiting, see also ^[14]. Unmanned manufacturing factories (UMFs) and robotic mobile fulfillment systems (RMFSs) are also very common applications of SOQNs; see, for example, ^{[3,15,16,17,18,19,20,21,22,23,24]}.

The vast bulk of analytical models of SOQNs currently available in the literature, as well as open and closed queueing networks, assume that the arrival flow of users is determined by the stationary Poisson process. However, it is widely acknowledged that modern real-world arrival processes are poorly fitted by such process. It may be too optimistic to predict real-world systems' performance metrics assuming that the stationary Poisson process describes the arrival flow. This is due to its constant arrival rate, relatively low variance, and zero coefficient of correlation for the sequential inter-arrival periods. Therefore, any irregularity in the arrival process, which is the inherent feature of the majority of real flows, is ignored. Namely, the existence of periods of peak rates leads to congestion and worsens the performance characteristics of the system.

As a much more adequate model of real-world arrival processes, M. Neuts ^[25] introduced the versatile arrival process, later called the Markov arrival process (MAP), see ^[26]. More information about the MAP and its extension, such as the batch Markov arrival process (BMAP), can be found, e.g., in ^{[27,28,29,30,31]}. To our knowledge, a SOQN with the MAP is only taken into consideration in a small number of articles; see the corresponding references in ^[7,8,9].

The MAP assumes that arriving users are of the same type. A more general model of the arrival process is the MMAP (marked Markov arrival process), in which arriving users are heterogeneous; see, e.g., ^[32]. A SOQN with multi-server nodes and the MMAP, where the type of user predefines the node in the network where the first service of the user is performed, was recently considered in ^[10].

In this paper, we analyze a SOQN with the MMAP and single-server nodes. The novelty of the model, compared to all existing literature that examines SOQN, consists of the possibility of changing the rate of service in the nodes of the network depending on the number of users admitted to the network.

The literature devoted to the analysis of various controlled queues, including control by the arrival rate, service rate, or number of active servers, is huge. This is easily explained by the fast technological development and consequent possibilities to monitor the state of a real-world system or network and dynamically change the service regime correspondingly. Such control can provide a significant economic effect due to the possibility of regime management in such a way that fast but expensive service regimes are temporarily used only in the case of high system congestion.

Control by the service rate in the nodes can be implemented via the use of equipment with distinct performance (e.g., leasing of channels or routers with different bandwidths in telecommunication networks) or the use of various numbers of workers providing service in a node. High-performance equipment and (or) more staff can be used when the quantity of users receiving service in the network is large, which may imply a long delay of the user in the network, possible dissatisfaction with the quality of service, and possible choice of another service provider in the future. When the network is less congested, it makes sense to use lower-performance equipment and (or) less staff to spend less money on service provision.

As early papers that examine queues with controlled service rates, we can mention ^[33,34]. When the system has several available service regimes listed in ascending order of service rate (and, correspondingly, increased cost), it is proved that the optimal control policy is monotone, i.e., of the multi-threshold type. This policy is defined by the set of integer numbers (thresholds). Selection of the service rate at decision moments is performed based on the relationship between the system's current user count and the predetermined thresholds. An application of such a type of policy to the BMAP/G/1 queue with a controlled service time, which can be varied at moments of service completion, was implemented in ^[35]. The BMAP/SM/1-type queueing system with a multi-threshold control by the service rate, semi-Markov service process, and user retrials in the case of a busy server was analyzed in ^[36].

One drawback of the threshold policy is that when the current number of users in the system approaches a certain threshold, there is a chance that the service rate will fluctuate often, a phenomenon known as oscillation. To overcome this negative phenomenon (due to the possible time loss or charge for the switching), the so-called hysteresis strategy has been offered. This type of strategy assumes that each policy threshold can be split into two thresholds. The rate of service increases when the larger of these two thresholds is exceeded. When the length of a queue drops to a smaller threshold, the service slows. The system stays in the prior regime when the queue length falls between the thresholds. This delay in switching (called hysteresis in Greek) leads to a decrease in the frequency of regime switching.

The M^x/G/1-type queueing system with two available service regimes, a batch stationary Poisson arrival process, a general distribution of service time, an account of switch-over times, and a hysteresis policy of control was considered in ^[37,38]. The result was generalized to the BMAP/G/1-type queueing system in ^[39]. An analysis of the queue with hysteresis control is more tricky than that of the same type of queue with threshold control. This is because, given the fixed values of the thresholds of the threshold policy, the choice of a service rate at a decision instant is uniquely determined by the current number of users in the network. In the case of the hysteresis policy of control, selection depends also on the service rate used before the moment of decision-making.

In this paper, we consider a user loss SOQN with arbitrary topology, an arriving flow defined by the MMAP, the availability of a finite number $L, \; L\ge 2,$ of service regimes, stochastic routing, and the hysteresis type of control by the service rates in the nodes. An exact algorithmic study of the steady-state dynamic of the network is implemented and numerically illustrated.

The structure of the paper is as follows. The description of the SOQN under study is presented in Section 2. The topology of the network is explained, a necessary notation is presented, and the hysteresis-type control policy is defined in detail. In Section 3, a multidimensional continuous-time Markov chain defining the system behavior under the fixed set of thresholds of the control strategy is introduced. The methodology tracing back to ^[40] is used here for keeping track of the number of users residing in all nodes of the network. The explicit form of the infinitesimal generator of this Markov chain is derived. Formulas for the calculation of a variety of performance measures of the network via the vectors of the steady-state probabilities of the Markov chain are given in Section 4. A numerical illustration of how the algorithmic results can be used for the analysis of the network consisting of 3 nodes, having 3 different service regimes, and an opportunity to simultaneously accommodate up to 40 users is presented in Section 5. Section 6 summarizes the results of the implemented study and briefly discusses possible directions for further research.

2.   An explanation of the model

Let us consider a SOQN that consists of $K$ nodes. The illustration of the network structure is presented in Figure 1.

Arrivals of users follow a marked Markov arrival process (MMAP). The irreducible underlying Markov chain $\{v_t, \; t \ge 0\}$ with the state space $\{1, \dots, V\}$ governs this process. The generator $H$ of this chain is represented by the form

where the diagonal entries of the matrix $H_0$ are negative. The rates at which the chain $\{v_t, \; t \ge 0\}$ exits the appropriate states are determined by their moduli. The non-diagonal entries of the matrix $H_0$ are nonnegative and specify the rates at which the chain $\{v_t, \; t \ge 0\}$ makes transitions without users arriving. The nonnegative entries in the matrix $H_k$ specify the rates at which the chain $\{v_t, \; t \ge 0\}$ transits when a user arrives at the network's $k$-th node, $k = \overline{1, K}.$ The expression $k = \overline{1, K}$ indicates that the integer variable $k$ takes values in the set $\{1, 2, \dots, K\}.$

The average rate of arrivals to the $k$-th node of the network is calculated as $\lambda_k = \boldsymbol \theta H_k {\bf e}, \; k = \overline{1, K},$ where the row vector $\boldsymbol \theta$ defines the stationary probability distribution of the states of the chain $\{v_t, \; t \ge 0\}$. The vector $\boldsymbol \theta$ satisfies the equations $\boldsymbol \theta H = {\bf 0}, \; \boldsymbol \theta {\bf e} = 1.$ Here, ${\bf 0}$ represents the row vector with 0s, and ${\bf e}$ represents the column vector with 1s. The total rate of arrivals is defined as $\lambda = \sum\limits_{k = 1}^{K} \lambda_k.$

Various details on the MMAP, its properties, and attributes, such as the distribution of moments of inter-arrival times, the coefficients of variation, and the correlation of inter-arrival times, can be found, for example, in ^{[27,28,29,30,32]}.

The network can process a maximum of $N$ users simultaneously. When there are $N$ users on the network, the user who arrives is regarded as lost. Every network's node functions as a queueing system having a single-server and a buffer with a capacity that is large enough to ensure that users are never lost.

The network can serve users in $L, \, L\ge 2,$ regimes. If the current service regime is $l$, the distribution of service time in the $k$-th node is exponential with intensity $\mu_k^{(l)}, \; l = \overline{1, L}.$ After service completion at the $k$-th node, with probability $p_{k, k'}$ the user transits for service to the $k'$-th node, $k' = \overline{1, K}, \, k'\neq k,$ or with probability $p_{k, 0}$ ends service in the network. Here, $p_{k, k} = 0, \; \sum\limits_{k' = 0}^K p_{k, k'} = 1$ for all $k, \; k = \overline{1, K}.$

The switching between the available service regimes is implemented by the following control strategy.

Let two sets of integer numbers (thresholds), $\{L_1^+, L_2^+, \dots, L_{L-1}^+\}$ and $\{L_1^-, L_2^-, \dots, L_{L-1}^-\},$ such that

be fixed.

We assume that the network instantly transitions to operating in the $(l+1)$-th regime if the number of users exceeds the threshold $L^+_l$ and the current regime is $l$. For example, if the network operates in the first service regime, the number of users in the network is $L^+_1,$ a new user is admitted to the network, and the service is immediately switched to the second regime.

If some user leaves the network when it operates in the $l$-th regime, $2\le l\le L$, and the number of users in the network drops to the value $L_{l-1}^-,$ the network immediately switches to operation in the $(l-1$)-th regime.

Let us present the sets listed below:

The network functions in the $l$-th regime if the quantity of users in the network equals $n$ and $n\in W_l, \, l = \overline{1, L},$ as stated in the model description. If $n\in F_l, \, l = \overline{1, L-1},$ the network can operate in both the $l$-th and $(l+1)$-th regime, and the used regime has to specified.

The impatience of users (leaving the queue without receiving service because they waited in the buffer for too long) is an essential feature of many real-world systems. For a survey of the literature dealing with queues with impatient users, see, e.g., ^[41]. To take the impatience phenomenon into account, we assume the impatience of the users waiting in the buffers of the nodes of the network. A user in the $k$-th buffer reneges after an amount of time that is exponentially distributed with intensity $\beta_k, \, k = \overline{1, K},$ independently from other users who are waiting. Such a user permanently reneges from the network (is lost).

We aim to analyze the invariant distribution of the network states and derive expressions for the calculation of the main performance characteristics of the network under the arbitrarily fixed values of the thresholds. As a result, a variety of optimization issues can be formulated and resolved.

3.   The process that describes the network states' dynamics

The continuous-time Markov chain

can be used to define the behavior of the considered network.

Here, during the instant $t,$

● the component $n_t$ defines the quantity of users in the network, $n_{t} = \overline{0, N};$

● the component $i_t$ is an indicator of the current regime of the network operation. It is defined only for the values $n_{t}$ such that $n_{t}\in F_l, \, l = \overline{1, L-1}.$ It admits values 0 or 1. Namely, the value $i_{t} = 0$ indicates the operation in the $l$-th regime, and $i_{t} = 1$ indicates the operation in the $(l+1)$-th regime;

● the component $v_{t}$ specifies the status of the underlying process of the MMAP, $v_{t} = \overline{1, V};$

● the component $m_{t}^{(k)}$ defines the quantity of users in the $k$-th node, $m_{t}^{(k)} = \overline{0, n_t},$ $\sum\limits_{k = 1}^{K} m_{t}^{(k)} = n_t,$ $k = \overline{1, K}.$

Let us call the set $\{ n_t, i_t, \, v_{t}, m^{(1)}_{t}, \dots, m^{(K)}_{t}\}$ of the states of Markov chain $\zeta_{t}, \, t\ge 0$ enumerated in the reverse lexicographic order of the processes $m_{t} ^ {(1)}, \dots, m_ {t} ^ {(K)}$ and the direct lexicographic order of the processes $(i_t, v_t)$ as level $n_t.$

The Markov chain $\{\zeta_ {t}, \; t \ge 0\}$ is regular and irreducible. Its state space is finite. Therefore, the invariant probabilities of the states of this chain

exist for all possible values of the queuing network's parameters.

The probability of the states that correspond to the level $n, \, n = \overline {0, N },$ can be formed as row vectors ${ \boldsymbol \pi}_n$. The probability vectors ${{ \boldsymbol \pi}}_n, \, n = \overline {0, N},$ can be computed as the solution of the system of linear algebraic equations (known as Chapman-Kolmogorov, balance, or equilibrium equations):

where the matrix $Q$ represents the generator of the Markov chain $\{\zeta_ {t}, \; t \ge 0\}$.

To calculate the vectors ${{ \boldsymbol \pi}}_n, \, n = \overline {0, N},$ the exact expression for the generator $Q$ must be obtained.

Let us present the following notation for future use in this paper:

● $O$ is a zero matrix of the suitable size, and $I$ is the identity matrix. If the size is not clear from context, it is indicated by a suffix, e.g., $I_K$ is the identity matrix of size $K$;

● the symbols $\otimes$ and $\oplus$ represent the Kronecker product and sum of matrices, respectively; refer to ^[42];

● ${\bf b}^{(k)}$ is a vector of size $K$ defined as ${\bf b}^{(k)} = \{\underbrace{0, 0, \dots, 0, 1}_k, 0, \dots, 0\}, \; k = \overline{1, K};$

● ${\rm diag}\{\dots \}$ denotes the diagonal matrix with the diagonal entries shown in brackets;

● $\Omega^{(l)}, \, l = \overline{1, L},$ is the square matrix of dimension $K$ defined as:

where the matrix $P$ with the entries $p_{k, k'}, \, k, k' = \overline{1, K},$ defines the one-step transition probabilities of a user within the network;

● ${\bf p}^{(l)}, \, l = \overline{1, L},$ is the column vector defined as ${\bf p}^{(l)} = (p_{1}^{(l)}, p_{2}^{(l)}, \dots, p_{K}^{(l)})^T$ where $p_{k}^{(l)} = p_{k, 0}\mu_k^{(l)}, \; {k = } \overline{1, K}$;

● ${ \boldsymbol \beta}$ is the column vector defined as ${ \boldsymbol \beta} = (\beta_{1}, \beta_{2}, \dots, \beta_{K})^T$.

Let us first introduce the set of matrices that describe the behaviour of the $K$-dimensional stochastic process ${\bf m}_t = \{m^{(1)}_{t}, \dots, m^{(K)}_{t}\}, {\; t \ge 0, }$ defining the quantity of users in each node of the network. To define the transition rates of this process, conditioned on the knowledge that $n$ users are present in the network working in the $l$-th regime at instant $t$, we need the following matrices:

$a)$ the matrix ${\cal T}_n (\Omega^{(l)})$ establishes the process ${\bf m} _ {t}, {\; t \ge 0, }$ of transition intensities at the precise instant when a user ends the service in one network's node and moves on to another node, $n = \overline{1, N}, \, l = \overline{1, L}$.

The following six steps make up the procedure for calculating the matrices ${\cal T}_n (\Omega^{(l)}), \; n = \overline{1, N}, \, l = \overline{1, L}:$

$(1)$ Determine the matrices $\Omega_k^{(l)}, \, k = \overline{1, K-2},$ that result from subtracting the $K-2-k$ first rows and columns from the matrix $\Omega^{(l)}.$

$(2)$ Calculate the set of matrices $X^{(r, l)}_{n, k}$ using the recursive formulas:

where $\omega_{i_1, i_2}^{k, l}$ is the $(i_1, i_2)$th entry of the matrix $\Omega_k^{(l)}$, and $r_k^{(l)}$ represents how many rows there are in the matrix $\Omega_k^{(l)}$.

$(3)$ Utilizing the recursive formulas, determine the set of the matrices $Y^{(r, l)}_{n, k}$:

$(4)$ Calculate the matrices $X_{n, k}^{(l)} = X^{(r_k-2, l)}_{n, k}, \, n = \overline{1, N},$ and $Y_{n, k}^{(l)} = Y^{(r_k-2, l)}_{n, k}$, $n = \overline{0, N-1}, \, k = \overline{1, K-2}, \, l = \overline{1, L}.$

$(5)$ Calculate the matrices ${\cal T}^{(k, l)}_{n} = {\cal T}^{(k, l)}_{n}(\Omega^{(l)})$ using the recursive formulas:

$(6)$ Calculate the matrices ${\cal T}_n(\Omega^{(l)})$ as

b) The matrix ${\cal S}_{n}({\bf p}^{(l)})$ defines the intensities of the process ${\bf m} _ {t}, \, t\ge 0,$ of transitions that occur at the moment of a user service completion in some node and its departure from the network, $n = \overline{1, N}, \, l = \overline{1, L}$.

The matrices ${\cal S}_{n}({\bf p}^{(l)}),$ $n = \overline{1, N}, \, l = \overline{1, L},$ can be found as

where the matrices ${\cal S}^{(K-1, l)}_{n}({\bf p}^{(l)})$ are recursively computed as:

c) The matrix ${\cal I}_{n}({ \boldsymbol \beta})$ defines the intensities of the process ${\bf m} _ {t}, \, t\ge 0,$ of transitions that happen at the moment of some user loss due to impatience, $n = \overline{2, N}$. The matrices ${\cal I}_{n}({ \boldsymbol \beta}),$ $n = \overline{2, N},$ can be found as ${\cal I}_n({ \boldsymbol \beta}) = {\cal I}^{(K-1)}_{n}({ \boldsymbol \beta}), \, n = \overline{2, N},$ where the matrices ${\cal I}^{(K-1)}_{n}({ \boldsymbol \beta}), \; \, n = \overline{2, N},$ are recursively obtained as

Hereinafter, we assume that ${\cal I}_{1}({ \boldsymbol \beta})$ is a zero matrix.

d) The probabilities of the process ${\bf m}_t$ transition at the time of a type-$k$ user admittance to the network are defined by the matrix ${\cal P}_n ({\bf b}^ {(k)}),$ $k = \overline{1, K}, \; n = \overline{0, N-1}.$ The paper ^[43] contains the algorithm for computing matrices ${\cal P}_n({\bf b}^{(r)})$. Using the denotations of our paper, we represent that algorithm below for the benefit of the reader. The matrices ${\cal P}_n({\bf b}), \, n = \overline{0, N-1},$ where ${\bf b} = (b_1, \dots, b_K), \; {\bf b}\in\{{\bf b}^{(1)}, {\bf b}^{(2)}, \dots, {\bf b}^{(K)}\},$ can be computed as

where the matrices ${\cal P}^{(k)}_{n}$ of size $(n+1)\times(n+2), \;$$n = \overline{1, N-1},$ are defined as

Here, the vectors $\tilde {\bf b}^{(k)}$ are defined as

e) The overall rate at which the process ${\bf m}_{t}$ exits the respective state is given by the moduli of the diagonal components of the diagonal matrix ${\cal E}_ {n}^{(l)},$ $n = \overline{1, N}, \, l = \overline{1, L}$. The matrices ${\cal E}_n^{(l)}$ are given by the formula

The presented formulas and algorithms used for the description of transition rates and probabilities of the process ${\bf m}_{t}$ are based on the ideas presented in ^[40] for the simultaneous description of the phase-type service processes in several independent servers.

Remark 1. A fact worth mentioning is that here, we compute the superfluous number of matrices ${\cal T}_n^{(l)}(\Omega^{(l)}),$ ${\cal E}_n^{(l)},$ and ${\cal S}_n^{(l)}({\bf p}^{(l)})$ because if the network operates in the $l$-th regime, then the number $n$ of users in the network should admit not an arbitrary value in the range from 0 to $N$ but should belong to the interval $[0, L^+_1]$, if $l = 1,$ $[L^-_{l-1}+1, L^+_{l}]$, if $l = \overline{2, L-1},$ and $[L^-_{L-1}+1, N]$, if $l = L.$ However, because we intend to build up the surfaces illustrating certain dependencies of the network performance indicators on the thresholds and solve the optimization problem, we would like to compute the generator $Q$ of the Markov chain $\zeta_ {t}, \; t \ge 0,$ for all possible sets of the thresholds $L_l^+$ and $L_l^-, \, l = \overline{1, L-1}.$ Therefore, if we compute all the matrices ${\cal T}_n^{(l)}(\Omega^{(l)}),$ ${\cal E}_n^{(l)},$ and ${\cal S}_n^{(l)}({\bf p}^{(l)})$ for every value of $n$ in the range from 0 to $N$ and all values of $l$ in the range from 1 to $L$ from the early beginning, we will avoid redundant repeated computations during the building of the surfaces and the solution of the optimization problem. If it is required to make calculations only for the fixed set of the thresholds, computation of these matrices can be implemented only for the values on $N$ from the corresponding diapason.

After calculating the matrices mentioned above, which completely characterize the service process of users ${\bf m}_t, {\; t \ge 0, }$ in all nodes of the network, we are ready to formulate the following statement.

Theorem 1. The infinitesimal generator $Q$ of the Markov chain $\zeta_ {t}, \; t \ge 0,$ under study has a block-tridiagonal structure.

The diagonal blocks $Q_{n, n}, \, n = \overline{0, N},$ of the generator are given by formulas:

The subdiagonal blocks $Q_{n, n-1}, \, n = \overline{1, N},$ are given by formulas:

The updiagonal blocks $Q_{n, n+1}, \, n = \overline{0, N-1},$ are given by formulas:

Proof. To implement the theorem's proof, every potential transition of the Markov chain $\zeta_ {t}$ over an infinitesimally short interval is analyzed, and the transition intensities are rewritten in block matrix form. The simultaneous transition rates of two independent Markov chains, $v_t$ and ${\bf m}_t,$ are defined via the Kronecker product of matrices.

The meaning of the blocks of the generator is transparent taking into account the described above probabilistic meaning of the matrices $H_0, \; H_1$ and ${\cal T}_n^{(l)}(\Omega^{(l)}),$ ${\cal E}_n^{(l)},$ and ${\cal S}_n^{(l)}({\bf p}^{(l)}).$ The off-diagonal components of the diagonal blocks $Q_{n, n}$ define the rates of transition of the Markov chain $\zeta_ {t}$ inside the level $n.$ The modules of the negative diagonal entries of the blocks $Q_{n, n}$ define the rates of the departure of the Markov chain $\zeta_ {t}$ from the states that belong to the level $n.$ The entries of the blocks $Q_{n, n-1}$ define transition intensities of the Markov chain $\zeta_ {t}$ from the level $n$ to the level $n-1.$ Such transitions can occur during service completion of some user in the network or user loss due to impatience. The entries of the blocks $Q_{n, n+1}$ define transition intensities of the Markov chain $\zeta_ {t}$ from the level $n$ to the level $n+1.$ Such transitions can occur at the moments of a new user arrival and admission. Different sizes of some blocks are explained by the possibility of transitions of the number $n$ of a level between the sets $W_l$ and $F_l$, which imply the change (between one and two) of the cardinality of the state space of the component $i_t.$ □

Remark 2. The number of blocks in the matrix $Q$ depends on the maximum quantity $N$ of users in the network but does not depend on the number $L$ of the available service regimes. But the size of these blocks may depend on $L$ and the relations between the thresholds $L_l^+$ and $L_l^-, \, l = \overline{1, L-1}.$ This size is minimal when the service rate control strategy is of the threshold, but not the hysteresis, type, i.e., $L_l^+ = L_l^-, \, l = \overline{1, L-1}.$ The square matrix $Q$ may have a huge size. Even in the simpler case of the threshold strategy, it is equal to $V\sum\limits_{n = 0}^N J_n$ where $J_m$ is the quantity of variants to distribute $m$ users among $K$ existing nodes, which is defined by the formula

Therefore, certain efficient algorithms that use the generator's sparse structure are required to solve the system (1). Specifically, the approach from ^[44] can be suggested to determine the queueing network's and Markov chain's stationary probability distribution.

4.   Performance measures

Following the computation of the probability vectors ${ \boldsymbol \pi}_ {n}, \; n = \overline{0, N},$ we may compute various performance metrics of the queuing network under study.

The mean number of users in the network is given by the following expression

The intensity of the output flow of successfully serviced users in the network can be found as

The output rate of users successfully serviced in the network from the $k$-th node is equal to

where ${\bf p}^{(l, k)}$ is a column vector of dimension $K$ that has all zero entries except the $k$-th entry $({\bf p}^{(l, k)})_k$, which is equal to ${p}^{(l)}_k$. The matrices ${\cal S}_n({\bf p}^{(l, k)})$ can be found using the same algorithm as for the matrices ${\cal S}_n({\bf p}^{(l)}).$

The average quantity of users in the $k$-th node, $k = \overline{1, K},$ is given by the following expression

where the matrices ${\cal J}_{n}({\bf b}), \, n = \overline{1, N},$ for the vectors ${\bf b} = (b_1, \dots, b_K), \;$ which take values from the set ${\bf b}\in\{{\bf b}^{(1)}, {\bf b}^{(2)}, \dots, {\bf b}^{(K)}\},$ can be found as

with the matrices ${\cal J}^{(K-1)}_{n}({\bf b}), \; \, n = \overline{1, N},$ recursively obtained as

The average quantity of busy servers in the $k$-th node, $k = \overline{1, K},$ is obtained using the formula

The average quantity of busy servers in the network is given by the following expression

The average quantity of users in the $k$-th node's buffer, $k = \overline{1, K},$ is given by the following expression

The average quantity of users waiting in all buffers of the network can be calculated using the formula

The probability that the network operates in the $l$-th regime at an arbitrary epoch is determined by the following expression

The intensity of the regime increase is given by the formula

The intensity of the regime decrease is given by the following expression

where

The average intensity of regime switching is equal to

The probability of an arbitrary user loss upon arrival caused by the residence of $N$ users in the network is calculated by the expression

The probability of an arbitrary type-$k$ user loss upon arrival due to the residence of $N$ users in the network is computed as

The probability of an arbitrary user loss upon arrival at the $k$-th node due to the residence of $N$ users in the network is given by the following expression

The probability of an arbitrary user loss due to impatience is equal to

The probability of an arbitrary user loss due to impatience in the $k$-th node is computed as

where ${ \boldsymbol \beta}^{(k)}$ is a column vector of size $K$ with all zero entries except the $k$-th entry $({ \boldsymbol \beta}^{(k)})_k$, which is equal to ${\beta}_k$.

The probability of an arbitrary user loss is given by the following expression

Controlling the accuracy of the calculation of the stationary distribution of the network states is made easier by the existence of two distinct expressions for computing the probability $P^{loss}$ and the average intensity of regime switching.

An arbitrary user's loss probability in the $k$-th node is determined as

An arbitrary user's probability of receiving successful service within the network is determined by the formula

5.   Numerical example

Let us show the numerical example that verifies the viability of the suggested methods and formulas, and partially highlights the impact of variation of the thresholds on the value of the key performance indicators of the system and the potential to apply the outcomes to managerial objectives.

In this example, we examine a queueing network with $K = 3$ nodes.

The MMAP flow of users arriving at the network is determined by the following matrices:

The average arrival rate for this arrival flow is $\lambda = 4.8606.$ The average arrival rate $\lambda_k$ and the coefficients of variation $c_{var}^{(k)}$ and correlation $c_{cor}^{(k)}$ of successive inter-arrival times to the $k$-th node, $k = 1, 2, 3,$ have the following values:

We assume that there are $L = 3$ possible service regimes of the network operation. Under the first regime, the service times in the nodes are exponentially distributed with parameters

respectively. Under the second regime, the parameters of the exponential distribution of the service time are

and under the third regime, the parameters are

The transition probabilities of users after the service completion in the nodes are defined as

The rates of the users' departure from the buffers of the nodes due to impatience are defined as follows:

In this numerical example, we assume that up to $N = 40$ users can obtain service in the network at the same time. The regime switching is defined by the four parameters $L^-_1, \, L^-_2, L^+_1$ and $L^+_2$.

The purpose of the numerical example is to show the dependence of the main network's performance measures on the switching parameters. However, there is no opportunity to build 5 D figures. Therefore, for better visualization of the results, let us fix the thresholds defining the rule of the switching between the first and the second regimes as $L^-_1 = 5$ and $L^+_1 = 10$ and vary the thresholds $L^-_2$ and $L^+_2$, defining the rule of the switching between the second and third regimes as follows: The threshold $L^+_2$ varies in the interval $[L^+_1+1, N)$, and the threshold $L^-_2$ varies over the interval $[L^+_1+1, L^+_2]$ with the same step 1.

Figures 2 and 3 illustrate the dependence of the average quantity $N^{network}$ and the average total number $N^{buf}$ of users in buffers on the parameters $L^-_2$ and $L^+_2.$

Complementary to Figure 2, the dynamics of $N^{network}$ are illustrated in Table 1 where values of $N^{network}$ are presented for values of $L^-_2$ and $L^+_2$ in some smaller range.

It is seen from these figures that the minimal values of $N^{network}$ and $N^{buf}$ are achieved for small values of the thresholds $L^-_2$ and $L^+_2.$ When these thresholds increase (this means that the third regime is used only for a larger number of users in the network), the values of $N^{network}$ and $N^{buf}$ increase quite sharply.

This is easily understandable, as the service rate during the use of the third regime is three times higher than during the use of the first regime and is 1.5 times higher than during the use of the second regime.

Figures 4–6 illustrate the dependence of the probabilities $P^{regime}_l$ on the fact that, at any arbitrary moment, the network operates in the $l$-th regime, $l = 1, 2, 3,$ on the parameters $L^-_2$ and $L^+_2.$

The maximum value of $P^{regime}_1$ is achieved for small values of $L^-_2$ and $L^+_2$ because such small values imply a more frequent use of the fastest, the third, regime of operation (this is confirmed by Figure 6) and higher chances that the number of users in the network will drop below the value $L_1^-+1$ and, therefore, the first regime will be used. The maximum value of $P^{regime}_2$ is achieved for large values of $L^-_2$ and $L^+_2$ because the second regime is used until the number of users in the network drops below the value $L_2^-+1,$ which implies the longer use of the second regime.

Figure 7 illustrates the dependence of the average intensity $\phi$ of regime switching on the parameters $L^-_2$ and $L^+_2.$

Figures 8–11 illustrate the dependence of the probabilities $P^{ent-loss}$ of a user loss at the entrance to the network and the probabilities $P_l^{ent-loss-arb}$ of a user loss at the entrance to the $l$-th node of the network, $l = 1, 2, 3,$ on the parameters $L^-_2$ and $L^+_2.$

Figures 12–15 illustrate the dependence of the probabilities $P^{imp-loss}$ of a user loss due to impatience in the network and the probabilities $P_l^{imp-loss}$ of a user loss due to impatience from the buffer of the $l$-th node of the network, $l = 1, 2, 3,$ on the parameters $L^-_2$ and $L^+_2.$

Figure 16 illustrates the dependence of the loss probability $P^{loss}$ of an arbitrary user (due to all reasons) on the parameters $L^-_2$ and $L^+_2.$

Complementary to Figure 16, the dynamics of $P^{loss}$ are illustrated in Table 2.

Because we have fixed $L^+_1 = 10$, and the thresholds $L^-_2$ and $L^+_2$ have to satisfy the inequalities $L^+_1 < L^-_2\le L^+_2 < N,$ the possible values of the thresholds $L^-_2$ and $L^+_2$ range from 11 to 39. The loss probability $P^{loss}$ reaches its minimum value of 0.07887 when $L^-_2 = L^+_2 = 11$. When $L^-_2 = L^+_2 = 39$, the maximum value of the loss probability $P^{loss}$ is reached and equals 0.23454. This fact is obvious because when $L^-_2 = L^+_2 = 11$, the network starts operation in the fastest, the third, service regime as soon as possible. When $L^-_2 = L^+_2 = 39,$ the third regime is switched on only when the number of users in the system is equal to the maximum admissible value of 40. Therefore, many users are lost at the entrance to the network and due to impatience. However, when deciding on the selection of the values of the thresholds, it is necessary to take into account that the use of faster service regimes by default is more costly compared to slower service regimes. Also, it is undesirable to frequently change service regimes because such a change in a real-world network can require some expenditures related to switching to another equipment or inviting or dismissing staff.

Therefore, optimization of the operation of the network requires an exact definition of the cost criterion. Let us assume that the following cost criterion is used to define the network's operational quality:

Here, $a$ is a profit obtained by the system via the service of one user, $b$ is a penalty paid by the network for one user loss upon arrival, $c$ is a penalty paid by the network for one user loss due to impatience, $e_l, \, l = \overline{1, L},$ is the cost of maintaining the $l$-th operation regime per unit time, and $d$ is a charge paid by the network for one switch of an operation regime.

Thus, the cost criterion $E$ represents the average network's revenue per unit of time. Our aim is to find the values of the thresholds $L^-_2$ and $L^+_2$ providing the maximum to the function $E(L^-_2, L^+_2).$

The cost coefficients in this numerical example are fixed at the following values:

Figure 17 shows how the thresholds $L^-_2$ and $L^+_2$ affect the cost criteria $E$.

When $L^{+}_2 = 20$ and $L^{-}_2 = 15,$ the cost criterion reaches its optimal value of $E^* = 5.19909.$ Thus, to obtain the maximal revenue under the fixed above values of the network parameters, it is necessary to switch from the second to the third service regime when the number of users in the network becomes equal to $(L^+_2)+1 = 21$ and switch back to the second service regime when the number of users drops to $L^-_2 = 15.$

All figures presented above illustrate the dependence of the cost criterion $E$ on the thresholds $L^-_2$ and $L^+_2$ under the fixed values of the thresholds $L^-_1 = 5$ and $L^+_1 = 10.$ Let us now assume that $L^-_1 = L^+_1 = L_1$ and $L^-_2 = L^+_2 = L_2,$ i.e., the hysteresis strategy turns into the threshold strategy. This means that if the number $n_t$ of users in the network does not exceed $L_1,$ then the network operates in the first regime. If the number $n_t$ belongs to the interval $(L_1+1, L_2]$, then the network operates in the second regime. If the number $n_t$ exceeds $L_2,$ then the network operates in the third regime.

Figures 18 and 19 illustrate the dependence of the loss probabilities $P^{ent-loss}$ of an arbitrary user upon arrival and $P^{imp-loss}$ of an arbitrary user due to impatience on the parameters $L_1$ and $L_2.$

Figure 20 illustrates the dependence of the cost criterion $E(L_1, L_2)$ on the parameters $L_1$ and $L_2.$

The optimal values of the thresholds are as follows: $L_1 = 0$ and $L_2 = 15.$ This means that the network uses the first regime only when it is empty; once a user arrives, it is serviced in the second regime. When the number of users in the network reaches the value of 16, the network starts operation in the third regime and maintains this regime until the number of users in the network drops to the value of 15. The maximal value of the cost criterion is 5.13969.

Let us return to the hysteresis-type control by the network operation. A numerical solution to the problem of determining the optimal value of the cost criterion under the network's fixed parameters can be found using the results obtained for the computation of the stationary distribution of the network states, the primary performance characteristics, and the cost criterion value. In our example, the value of the cost criterion is the function of four thresholds, $L^-_1, \; L^+_1, \; L^-_2$, and $L^+_2.$ It can be computed that the maximum value of the cost criterion is equal to 5.31252 and is achieved for the following values of the thresholds: $L^-_1 = 0, L^+_1 = 2, L^-_2 = 13, L^+_2 = 18.$ The achieved value of 5.31252 of the cost criterion is higher than the optimal value of 5.13969 of the criterion under the use of the optimal threshold strategy due to the more seldom switching of regimes and the presence of a charge for the switch of the regimes in the cost criterion.

In a more general case, when the number $L$ of available regimes is higher than two (and the number of the thresholds is 2$L$), the problem of finding the maximal value of the criterion and the optimal values of the corresponding thresholds can be deeply complicated by the existence of a huge number of possible combinations of the threshold values. Therefore, this problem deserves a separate consideration. The use of some derivative-free methods of optimization, see, e.g., ^[45,46], can be recommended. As mentioned above, the value of our results consists in providing the possibility to exactly compute the value of the cost criterion for any fixed set of control strategies during the implementation of the search for the optimal values of the thresholds.

It is worth noting that we assumed that the values of the service rates $\mu_k^{(l)}, \; k = \overline{1, K}, \; l = \overline{1, L},$ in the nodes of the network under the fixed regimes of the network operation are fixed. In potential real-world applications, these values can also not be fixed but have to be chosen from some set. Our results can be used to optimize the selection of these rates' values and the corresponding threshold values.

Remark 3. Described above computations were implemented using Wolfram Mathematica on a Lenovo notebook with an Intel(R) Core(TM) i7-1165G7 2.80GHz and 16 GB RAM. Running time for computation of the optimal value of the cost criterion $E(L^-_2, L^+_2)$ defined by formula (3) was equal to 8419 seconds, i.e., 18 seconds for one point $(L^-_2, L^+_2)$ (the total number of points is 465). Because this computation time was acceptable for preparation of the presented examples, no optimization of the code was made. Computation time can be significantly reduced via such an optimization and the use of a more powerful notebook or PC.

A quite long computation time is explained by the large size of the generator $Q$. As mentioned above in the simplest case of the threshold strategy, this size is equal to $V\sum\limits_{n = 0}^N J_n,$ where the numbers $J_n, \; n = \overline{0, N},$ are defined by formula (2). In the considered example of $SOQN$, where we assume $K = 3$ nodes and admission of up to $N = 40$ customers to the network simultaneously, the number $J_N$ is equal to 861. If we consider the network consisting of four nodes and decreased $N$ to 15, we will have about the same (816) size of the block $J_N$ and a similar computation time.

6.   Conclusions

We considered a user loss SOQN with the bursty MMAP-type arrival process, single-server nodes with a controlled service regime, a hysteresis-type control policy, and impatient users. Under the fixed parameters of the control policy, the stationary behavior of this SOQN is described by a multidimensional Markov chain, whose components define the total number of users in the network, the used service regime (if it is not uniquely defined by the number of users in the network), the state of the underlying process of the MMAP, and the number of users in each node of the network. The generator of this chain is obtained as a block tridiagonal matrix. Formulas for computation of the key performance measures of the network are derived. Numerical illustrations of the algorithmic and analytical results obtained are provided.

The results can be applied to managerial objectives, such as the selection of possible variants of service organization in the nodes, including the choice of the suitable equipment and the corresponding staff; routing of users in the network; pricing; and the optimal dynamical scheduling of the variants of service organization depending on the current load of the network.

The results can be generalized into several directions, such as the consideration of back-ordering SOQNs having an infinite or finite buffer for storing users who did not succeed in entering the core network upon arrival because it was completely busy; an account of possible user's impatience during waiting in this buffer; the presence of an orbit for user retrials; or the arrival and impatience rate control. The results from ^[47] are planned to be used for implementing these generalizations.

Author contributions

Ciro D'Apice: Conceptualization, Methodology, Validation, Investigation, Writing-original draft, Writing-review & editing, Project administration; Alexander Dudin: Conceptualization, Methodology, Formal analysis, Investigation, Writing-original draft, Writing-review & editing, Project administration; Sergei Dudin: Methodology, Software, Formal analysis, Investigation, Supervision, Writing-original draft; Rosanna Manzo: Conceptualization, Software, Validation, Formal analysis, Investigation, Writing-original draft, Writing-review & editing. All authors have read and approved the final version of the manuscript for publication.

Use of AI tools declaration

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

Conflict of interest

Rosanna Manzo is the Guest Editor of special issue "Managing complex systems by simulation and optimization techniques" for AIMS Mathematics. Rosanna Manzo was not involved in the editorial review and the decision to publish this article.