<figure><img src="image_1.png" /></figure>

HAL
open science

# Operational availability prediction for a fleet of large systems; an approximate solution

Raymond A Marie

<figure><img src="image_2.png" /></figure>

## ► To cite this version:

Raymond A Marie. Operational availability prediction for a fleet of large systems; an approximate solution. ICIL 2023 - International Conference on Industrial Logistics, May 2023, Zadar, Croatia. pp.1-8. hal-04389976

HAL Id: hal-04389976

https://inria.hal.science/hal-04389976v1

Submitted on 12 Jan 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# OPERATIONAL AVAILABILITY PREDICTION FOR A FLEET OF LARGE SYSTEMS; AN APPROXIMATE SOLUTION

Raymond A. MARIE

University of Rennes, Inria & IRISA, Campus de Beaulieu, 35042 Rennes Cedex, France

Raymond.Marie@irisa.fr

**Abstract:** We consider a fleet of systems (e.g. a fleet of aircraft of a certain model) which are returned to operational status by exchanging an easily replaceable unit on site or nearby (provided that an identical unit in good condition is available). Such units are called LRUS (on Line Replaceable Units). These LRUS must be exchanged by qualified personnel who are multi-skilled but not experts in the internal functioning of the different types of LRUS, and this in a relatively short time (a few hours at most). The major goal of this study concerns the number of systems blocked by unavailability of spare LRU (in the form of the probability distribution or expectation of this variable). The explosion in the number of states required for some real-world applications limits the use of a classical Markovian approach to model and evaluate the performance of fleet systems. This is the reason why we propose here an approximate recurrent method to calculate the operational unavailability. Its low relative complexity allows it to be used for applications encountered in the maintenance domain. Although approximated, the method provides the result with a small relative error, ranging from $10^{-2}$ to $10^{-8}$. Although generally less significant, the downtime associated with disassembly and reassembly is also taken into account in this study.

**Keywords:** Performance evaluation, Integrated logistic system, Stock shortage, Line replaceable unit, Operational availability, Intrinsic availability, Proactive queues, Non product-form queueing network.

## 1 INTRODUCTION

We consider a fleet of systems (e.g. a fleet of aircrafts of a certain model) which are returned to operational status by exchanging an easily replaceable unit on site or nearby (provided that an identical unit in good condition is available). Such units are called LRUS (on Line Replaceable Units). These LRUS must be exchanged by qualified personnel who are multi-skilled but not experts in the internal functioning of the different types of LRUS, and this in a relatively short time (a few hours at most). These time and skill constraints mean that LRUS are generally large and expensive subsystems (such as an aircraft engine). After the exchange, the defective LRU is directed to a repair center where it is either brought back to operational level or declared non-repairable and redirected to recycling.

Considering the high cost of most LRUS, and possibly their good reliability, they are supplied in (very) small quantities at the operational site of the systems or at a nearby site. This results in system unavailability due to a lack of LRUS. When looking for the reasons for a poor operational availability of a fleet of systems, LRU waiting times take precedence over manpower shortages.

In order to predict the performance of such organizations, two approaches are a priori possible. The first is to build complex models with a large number of variables solved by simulation. But this approach requires significant human resources to build the model and test it; it is then necessary to simulate a large number of events to obtain metrics with small confidence intervals (because some types of events are

in the category of rare events). The second approach is to look for a more analytical model such as a queueing network. Unfortunately, the class of product-form queueing networks, for which known algorithms can efficiently handle high-dimensional models, does not allow for accurate model modeling the phenomena of no available LRUS. Of course, Markovian modeling can theoretically describe any discrete event system, but the systems we wish to evaluate here would require dealing with continuous time Markov chains (CTMCs) with state space cardinality in the millions or billions of states. Faced with what we can currently call dead ends, we must consider the search for an approximate solution with the hope of finding one that provides results close to the exact results (on examples of modest dimensions to be computable). Note that many articles exist on inventory management in the context of performance evaluation of maintenance or production systems, and which do not use simulation (cf. [1]); but, unless the author is mistaken, they do not take into account the interdependence of unavailable LRU stocks on the operational availability of systems.

This paper, which extends previous results (cf. [2] and [3]), is organized as follows. In Section 2, we specify the assumptions and the overall model which consists of a queueing network that does not belong to the class of product-form networks. In Section 3, we introduce a Markovian representation providing the exact stationary distribution of the number of systems blocked by unavailability of spare LRU. The approximate method for dealing with large models is discussed in Section 4. In Section 5, the results provided by the approximate method are used to finally obtain an estimate of the op-
erational availability of the system fleet. Section 6 is devoted to numerical tests to ensure that the proposed approximate method provides results close to those provided by the reference model. Finally, Section 7 constitutes the conclusion of this study.

## 2 ASSUMPTIONS AND REFERENCE MODEL

We assume that each system consists of $N$ LRU in series. We also assume that $M$ systems are deployed at the operational site. Given the serial structure of a system, we assume here that when a system fails due to the failure of one LRU, the other LRUs cannot fail until the system is back in service. The fleet of systems is considered in steady state and the study focuses on the asymptotic behavior of the features.

$R_i$ denotes the number of LRUs of type $i$ initially allocated as spares at the operational site, $i = 1, 2, ..., N$. Let us introduce the random variable (RV) $S_i$, $i = 1, 2, ..., N$, defined as the quantity of LRU of type $i$ available on the site at a given time. The possible values of the RV $S_i$ are thus the integers $\{0, 1, ..., R_i\}$. Given the high price of the LRU, the failure of a LRU leads to the ordering of a new LRU from the storage center in charge of the operational site. At a given time, several systems may be blocked due to a lack of LRU of type $i$. The RV $X_i$, $i = 1, 2, ..., N$, designates here the number of systems blocked due to a lack of LRU of type $i$ at a given time. Each RV $X_i$ takes its values in the set $\{0, 1, 2, ..., M\}$.

To model the behavior of LRUS stocks, we use the proactive queueing model shown in Figure 1. The presented queue is a $M/M/r < C >$ where $C$ denotes the maximal value of service anticipation (cf. [4]) and $r$ is the number of servers. The variable $j$ (resp. $i$) denotes the number of clients present, (resp. the number of anticipated services available).

In our case, for the LRUS of type $i$, the number of clients present corresponds to the RV $X_i$ and the number of anticipated services available corresponds to the RV $S_i$. Note that, by definition, the product $S_i \times X_i$ is always zero. Now consider the RV $W_i$, $i = 1, 2, ..., N$, defined as the random quantity $W_i = R_i - (S_i - X_i)$. At any time, this quantity corresponds to the number of type $i$ LRUS waiting to be delivered to the operational site.

<figure><img src="image_4.png" /><figcaption>Figure 1: Representation of the proactive queue $M/M/r < C >$.</figcaption></figure>

Our global model is the queueing network represented in Figure 2. The lower sub-network models the $N$ stocks relative to the $N$ types of LRU. This lower sub-network is not a product-

form sub-network. But we assume that the upper sub-network is a product-form sub-network. For that we assume that the uptime of the LRUS follow exponential distributions and that all the other stations of the upper sub-network satisfy the necessary conditions. The failure rate for type $i$ LRU is denoted by $\lambda_i$, $i = 1, 2, ..., N$.

Let us name NBS the RV "number of systems blocked by lack of available LRUS"; we have $NBS = \sum_{i=1}^{N} X_i$. The expectation of NBS is the portion of unavailability attributable to the lack of available LRUS.

Let us put $K_i = M + R_i$, $i = 1, 2, ..., N$; it is easy to see that the possible values of the $W_i$ are the integers $\{0, 1, ..., K_i\}$ and that the knowledge of the couple $(S_i, X_i)$ gives us the value of the RV $W_i$, and vice versa. Note that if $R_i = 0$, we simply get $W_i = X_i$.

We assume here that the delivery rate of a LRU of type $i$ on the site depends only on the type $i$ of the LRU and on the quantity waiting for delivery (quantity equal to $W_i$); which will allow us to consider a homogeneous Markovian model. The delivery rate of a LRU of type $i$ to the site, knowing that $W_i = k_i$, is denoted $\beta_i(k_i)$, $k_i = 1, 2, ..., K_i$, $i = 1, 2, ..., N$. The values of these rates are here assumed to be known (although we are aware of the required experimental work).

<figure><img src="image_5.png" /><figcaption>Figure 2: Structure of the global queueing network.</figcaption></figure>

## 3 EXACT STEADY STATE DISTRIBUTION OF THE RANDOM VARIABLE NBS

The first step is now to aggregate the upper sub-network admitting a product-form into a single station admitting an output rate as a function of the number of systems present in this subnetwork. Let $\nu(n)$ denotes this output rate, $n =$
0, 1, 2, ..., M; with of course $\nu(0) = 0$. The determination of these output rates is known for many years (cf. [5] and [6]): the expression $\nu(n)$ equals $\gamma_{out}G(n-1)/G(n)$. $G(\cdot)$ is the discrete time convolution function built on the product-form upper sub-network and $\gamma_{out}$ is, up to an arbitrary constant, the output routing rate obtained thanks to the routing probability matrix $P_r = (r_{ij})$, where $r_{ij}$ is the probability that a customer leaving station i joins station j, and thanks to the vector $\mathbf{x}$, solution of $\mathbf{x} = \mathbf{x}P_r$.

The second step deals with the non product-form sub-network. The continuous time Markov chain (CMTC) whose states are the n-tuples $(W_1, W_2, ..., W_N)$ constitutes our reference model. Let us denote its states by the n-tuples $(k_1, k_2, ..., k_N)$. The state space $\mathbb{E}$ corresponds to the set :

$$ \mathbb{E} = \{(k_1, k_2, ..., k_N) | 0 \le k_j \le K_j, j = 1, ..., N \\ \text{and } \sum_{j=1}^{N} k_j \le (M + \sum_{j=1}^{N} R_j)\}. \quad (1) $$

To a state $(k_1, k_2, ..., k_N)$ will correspond a vector $(X_1, X_2, ..., X_N)$ whose component values are given by the relation $X_i = \max(0, (k_i - R_i))$, $i = 1, 2, ..., N$. In this paper, we will use the simplified notation $(k_i - R_i)^+ \triangleq \max(0, (k_i - R_i))$, $i = 1, 2, ..., N$. When the CMTC is in a state $(k_1, k_2, ..., k_N)$, the number of blocked systems, noted here $b(k_1, k_2, ..., k_N)$, is given by the relation :

$$ b(k_1, k_2, ..., k_N) = \sum_{j=1}^{N} (k_j - R_j)^+. \quad (2) $$

Thus, when the system is in a $(k_1, k_2, ..., k_N)$ state of $\mathbb{E}$, the failure rate of LRU of type i is equal to

$$ \nu(M - b(k_1, k_2, ..., k_N)) \frac{\lambda_i}{\Lambda}, \text{ where } \Lambda = \sum_{i=1}^{N} \lambda_i. \quad (3) $$

Recall again that, by assumption, the LRU of blocked systems do not continue to fail.

Let us introduce two vector notations: $\mathbf{k}$ denotes the (line) vector $(k_1, k_2, ..., k_N)$ and $\mathbf{e}_i$ denotes a vector whose components are all zero except for the i-th which is equal to unity, the dimension of $\mathbf{e}_i$ being defined by the context. Let us denote $\mathbf{A}$ the infinitesimal generator of the CMTC and $a_{\mathbf{k},\mathbf{k}'}$ the elements of $\mathbf{A}$. The values of the rates $a_{\mathbf{k},\mathbf{k}'}$, for $\mathbf{k}' \neq \mathbf{k}$ are given by the relation :

$$ a_{\mathbf{k},\mathbf{k}'} = \begin{cases} \beta_i(k_i) & \text{if } \mathbf{k}' = \mathbf{k} - \mathbf{e}_i \text{ and } k_i > 0, \\ & i = 1, 2, ..., N \\ \nu(M - b(\mathbf{k})) \frac{\lambda_i}{\Lambda} & \text{if } \mathbf{k}' = \mathbf{k} + \mathbf{e}_i, \\ & i = 1, 2, ..., N \\ 0 & \text{else} \end{cases} \quad (4) $$

By letting us put by convention $\beta_i(0) = 0$, we obtain for the

diagonal elements $a_{\mathbf{k},\mathbf{k}}$:

$$ a_{\mathbf{k},\mathbf{k}} = - \sum_{i=1}^{N} \left( \beta_i(k_i) + \nu(M - b(\mathbf{k})) \frac{\lambda_i}{\Lambda} \right). \quad (5) $$

The knowledge of the infinitesimal generator $\mathbf{A}$ allows us to obtain the stationary distribution of the states $(k_1, k_2, ..., k_N)$ of $\mathbb{E}$ and thus the stationary distribution of the random variable NBS. This model seems to us to be a relatively accurate representation of the studied system. Unfortunately, the drastic increase of the number of states with the number of types of LRU constitutes the limit of this reference model. For example, when $N = 6$, $R_i = 1$ and $K_i = 40$, $i = 1, 2, ..., N$, the number of states of the CMTC is equal to 16 154 399. It is this limitation that led us to look for the approximate method exposed in the following Section. But this reference resolution will allow to test the accuracy of the approximated method on more modest examples.

## 4 DETERMINATION OF THE APPROXIMATE METHOD

The proposed model is constructed in a recursive manner. Schematically, in step $n$, we integrate a fictitious system consisting of $n$ types of LRU in series with a $(n+1)$-th LRU. At the beginning of the step, the behavior of the fictitious system consisting of $n$ types of LRU is characterized by a set of $M$ birth and death processes. At the end of the step, the behavior of the new fictitious system consisting of $(n+1)$ types of LRU is characterized by a new set of $M$ birth and death processes. The $(N-1)$-th and last step, which is simpler, allows to determine an approximate stationary distribution of the random variable NBS.

Slightly particular, the first step constructs the fictitious system comprising only the first two types of LRU. In a first phase, this system is modeled by a CMTC whose states correspond to the possible values of the couples $(W_1, W_2)$, i.e. to the couples $(k_1, k_2)$ of the set $\mathbb{E}_2 = \{(k_1, k_2) | 0 \le k_1 \le K_1, 0 \le k_2 \le K_2, \text{ et } k_1 + k_2 \le (M + R_1 + R_2)\}$. For a state $(k_1, k_2)$ of this fictitious system, the number of blocked systems, noted $b(k_1, k_2)$, is equal to :

$$ b(k_1, k_2) = \sum_{i=1}^{2} (k_i - R_i)^+. \quad (6) $$

The cardinality of the state space of this CMTC is equal to $|\mathbb{E}_2| = (K_1 + 1)(K_2 + 1) - \frac{M(M+1)}{2}$. The quantity $-\frac{M(M+1)}{2}$ corresponds to the number of forbidden states $(k_1, k_2)$ due to the constraint $k_1 + k_2 \le (M + R_1 + R_2)$.

Let $A_2$ be the infinitesimal generator of the CMTC. Figure 3 shows the transition rates in the general case.

The Figure 4 gives the structure of the transition graph in the particular case where $M = 4$, $R_1 = 2$, and $R_2 = 1$; the edges in solid line schematize the existence of transitions related to both failures and deliveries of LRU. The vertically
hatched states are the states of the CMTC for which no system is blocked due to lack of LRU and the horizontally hatched states are the states of the CMTC for which no further failures can occur because all systems are already non-operational. The dashed edges schematize the transitions related to LRU deliveries only.

<figure><img src="image_6.png" /><figcaption>Figure 3: Synthetic representation of the transition rates of the CMTC associated with the couple $(k_1, k_2)$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 4: Structure of the transition graph of the CMTC associated to the study of the pairs $(W_1, W_2)$. Special case $M = 4$, $R_1 = 2$, and $R_2 = 1$.</figcaption></figure>

Once the matrix $A_2$ is known, we determine the vector of asymptotic state probabilities of the CMTC by solving the usual linear system.

Let us denote $p(k_1, k_2)$ the asymptotic probability of the state $(k_1, k_2)$ of the CMTC. These marginal probabilities of $X_1$ and $X_2$ are written:

$$p_1(j) = \begin{cases} \sum_{k_1=0}^{R_1} \sum_{k_2=0}^{K_2} p(k_1, k_2) & \text{if } j=0, \\ \sum_{k_2=0}^{K_2-j} p((R_1+j), k_2) & \text{if } 0 < j \le M \end{cases} \quad (7)$$

and

$$p_2(j) = \begin{cases} \sum_{k_2=0}^{R_2} \sum_{k_1=0}^{K_1} p(k_1, k_2) & \text{if } j=0, \\ \sum_{k_1=0}^{K_1-j} p(k_1, (R_2+j)) & \text{if } 0 < j \le M \end{cases} \quad (8)$$

Thanks to Figure 4, it is easy to understand that for any state $(k_1, k_2)$ such that $k_1 \le R_1$, $X_1$ will be null (this explains the expression of $p_1(0)$), while if $k_1 > R_1$, the maximal value of $k_2$ is $(K_2 - j)$.

At this stage, the basic idea consists in constructing from this CMTC associated to the couple $(W_1, W_2)$, a birth and death process (BDP) in order to approximate the behavior of this couple by a unique RV that we will note $\tilde{W}_2$. Noting $R_{min} = \min(R_1, R_2)$, we notice that we can build a birth and death process by considering the partition $\{\mathbb{E}_{0,0}, \mathbb{E}_{0,1}, \dots, \mathbb{E}_{0,R_{min}}, \mathbb{E}_{1}, \dots, \mathbb{E}_{M}\}$ of the set of states of the CMTC such that

$$\begin{aligned} \mathbb{E}_{0,n} &= \{(k_1, k_2) | (k_1, k_2) \in \mathbb{E}, \\ & \quad \text{and } k_1 = R_1 - R_{min} + n, \\ & \quad \text{and } k_2 \le R_2 - R_{min} + n\} \\ & \cup \{(k_1, k_2) | (k_1, k_2) \in \mathbb{E}, \\ & \quad \text{and } k_2 = R_2 - R_{min} + n, \\ & \quad \text{and } k_1 \le R_1 - R_{min} + n\} \\ & n = 0, 1, 2, \dots, R_{min} - 1. \end{aligned} \quad (9)$$

and

$$\mathbb{E}_i = \{(k_1, k_2) | (k_1, k_2) \in \mathbb{E}, \\ \text{and } b(k_1, k_2) = i\} \quad i = 1, 2, \dots, M. \quad (10)$$

On Figure 5 we have specified some examples of classes in order to illustrate the construction of the selected birth and death process.

We can verify that if the CMTC leaves the subset $\mathbb{E}_i$ at time $t$, it ends up at time $t^+$ either in the subset $\mathbb{E}_{i+1}$ or in the subset $\mathbb{E}_{i-1}$. The transition rates of the birth and death process are identified with the conditional transition frequencies from one class $\mathbb{E}_i$ to the class $\mathbb{E}_{i-1}$ or to the class $\mathbb{E}_{i+1}$. In the same way, if the CMTC leaves the subset $\mathbb{E}_{0,n}$ at time $t$, it finds itself at time $t^+$ either in the subset $\mathbb{E}_{0,n+1}$ (or $\mathbb{E}_1$ if $n = R_{min}$), or in the subset $\mathbb{E}_{0,n-1}$ (if $n > 0$). We have thus defined a birth and death process of cardinality $(M + R_{min} + 1)$.
<figure><img src="image_8.png" /><figcaption>Figure 5: Illustration of the classes of the partition on the structure of the transition graph of the CMTC associated with the study of the pairs $(W_1, W_2)$. Special case $M = 4$, $R_1 = 2$, et $R_2 = 1$.</figcaption></figure>

Since the asymptotic state probabilities $p(k_1, k_2)$ of the CMTC are known, we can compute the transition rates of the birth and death process of cardinality $(M + R_{min} + 1)$ associated with the RV $\tilde{W}_2$. Let us denote $\tau^{(M)}$ (resp. $\theta^{(M)}$) the vector whose $(M + R_{min})$ components correspond to the transition rates associated with failures (resp. deliveries of LRU).

Then we have to repeat $(M - 1)$ times this first phase of step 1 to determine the vectors $\tau^{(m)}$ and $\theta^{(m)}$ associated to a number $m$ of (pseudo) systems formed by the two first types of LRU, $m = 1, 2, \dots, (M - 1)$. At the end of the first step, $M$ birth and death processes of respective cardinality $(m + R_{min} + 1)$ and their $2m$ vectors $\tau^{(m)}$ and $\theta^{(m)}$, $m = 1, 2, \dots, M$, have been computed. To simplify the presentation of the second step, we let us put $\tilde{R}_2 = R_{min}$ and $\tilde{K}_2 = R_{min} + M$.

The next step consists in building a CMTC $(\tilde{W}_2, W_3)$ whose states correspond to the couples $(\tilde{k}_2, k_3)$ of the set $\mathbb{E} = \{(\tilde{k}_2, k_3) | 0 \le \tilde{k}_2 \le \tilde{K}_2, 0 \le k_3 \le K_3, \text{ and } \tilde{k}_2 + k_3 \le (M + \tilde{R}_2 + R_3)\}$. The transition graph of this CMTC is similar to the one constructed to study the couple $(W_1, W_2)$ and represented in Figure 4. Note that for a value $k_3$ of the variable $W_3$, the possible values of the variable $\tilde{W}_2$ are the values $\tilde{k}_2 = 0.1, \dots, (\tilde{K}_2 - (k_3 - R_3))^+$.

A subtlety of this CMTC $(\tilde{W}_2, W_3)$ is the following. For a fixed value $k_3$ of the variable $W_3$ such that $k_3 < K_3$, the transition rates associated with the variable $\tilde{W}_2$ are those of the birth and death process obtained with the value $m(k_3) = \tilde{K}_2 - (k_3 - R_3)^+ - \tilde{R}_2$, the one whose cardinality of states is equal to $(\tilde{K}_2 + 1 - (k_3 - R_3)^+)$ and whose transition rates

are given by the vectors $\tau^{(m(k_3))}$ and $\theta^{(m(k_3))}$. Note that the $M$ birth and death processes determined earlier are needed to study the CMTC $(\tilde{W}_2, W_3)$. Since $m(k_3)$ is zero when $k_3 = \tilde{K}_3$, there is no birth and death process to provide the transition (death) rates associated with the variable $\tilde{W}_2$ and necessary if $\tilde{R}_2$ is positive. In this case, we use the $\tilde{R}_2$ first components of the vector $\theta^{(1)}$. Figure 6 shows the transition rates in the general case.

The idea behind the approach is to think that the variable $\tilde{W}_2$, which can be seen describing the behavior of a super-LRU, will behave from the point of view of availability in a way quite close to the first two real LRU.

<figure><img src="image_9.png" /><figcaption>Figure 6: Synthetic representation of the transition rates of the CMTC associated to the couple $(\tilde{k}_2, k_3)$.</figcaption></figure>

After determining the stationary distribution of the CMTC $(\tilde{W}_2, W_3)$, we can consider the birth and death process of cardinality $(M + R_{min} + 1)$, where $R_{min} = \min(\tilde{R}_2, R_3)$, describing the behavior of the $\tilde{W}_3$, and determine the new vectors $\tau^{(M)}$ and $\theta^{(M)}$. As in the previous step, we must study the $M$ CMTC $(\tilde{W}_2, W_3)$ associated with the values $m, m = 1, 2, \dots, M$.

Having thus the vectors $\tau^{(m)}$ and $\theta^{(m)}$, $m = 1, 2, \dots, M$, we can study the stationary distribution of the new CMTC $(\tilde{W}_3, W_4)$, with parameters $\tilde{R}_3 = R_{min}$ and $\tilde{K}_3 = \tilde{R}_3 + m$, $m = 1, 2, \dots, M$ in order to determine the new vectors $\tau^{(m)}$ and $\theta^{(m)}$, $m = 1, 2, \dots, M$.

In general, we study the stationary distribution of new CMTC $(\tilde{W}_i, W_{i+1})$, with the parameters $\tilde{R}_i = \min(\tilde{R}_{i-1}, R_i)$ and $\tilde{K}_i = \tilde{R}_i + m$, $m = 1, 2, \dots, M$ in order to determine the new vectors $\tau^{(m)}$ and $\theta^{(m)}$, $m = 1, 2, \dots, M$. For the last step, $i = N - 1$, it is sufficient to determine the stationary distribution of the CMTC $(\tilde{W}_{N-1}, W_N)$ associated with $m = M$ to obtain the desired metrics, such as the distribution of the RV NBS and its expectation corresponding to an approximation of the unavailability in asymptotic regime. The skeleton of the
general algorithm of the method is represented table 1.

<caption>Table 1: General algorithm skeleton.</caption>

Begin
Initialize the vectors $\tau^{(m)}$ and $\theta^{(m)}$, using the BDP associated with RV $W_1$ when the $m$ pseudo-systems include only the LRU of type 1, $m = 1, 2, ..., M$.
$n = 1$
While $n < (N - 1)$ Do
For $m = 1, 2, ..., M$ Do
Determine the vectors $\tau^{(m)}$ and $\theta^{(m)}$ relative to the BDP deducted from CMTC $(\widetilde{W}_n, W_{n+1})$ for which the number of pseudo-systems is equal to $m$.
EndFor
$n = n + 1$
EndWhile
Study the CMTC $(\widetilde{W}_{N-1}, W_N)$ to deduce the distribution and the expectation of the NBS.
End

This algorithm is much faster and less memory consuming than the algorithm providing the exact solution of the NBS distribution. Moreover, the low relative error of this algorithm (which will be discussed in Section 6) offers the possibility to perform a dynamic allocation of LRUs in order to minimize the NBS expectation under a budget constraint. Concerning the complexity of this new algorithm, it is of the order $O(N.(M + \bar{R})^2)$, where $\bar{R}$ is the average value of the sequence $R_1, R_2, ..., R_N$.

Let us specify that if $R_i = 0$ for $i = 1, 2, ..., (N - 1)$, this method provides the exact solution (identical to that of the reference model). Moreover, some simplifications are possible since in this case the common components of the vectors $\tau^{(m)}$ (resp. $\theta^{(m)}$) and $\tau^{(m')}$, (resp. $\theta^{(m')}$) take identical values.

# 5 DETERMINING THE OPERATIONAL AVAILABILITY OF THE FLEET OF SYSTEMS

At this stage, we do not yet know the operational availability of the systems because we must also take into account the phases of removing the failed LRUs and reassembling the new ones on the systems. This is the purpose of this last step.

Here, we consider that the mean number of non blocked systems is equals $(M - \text{NBS})$. This is the expected number of systems present in the exponential sub-network. Remembering that $\Lambda = \sum_{i=1}^{N} \lambda_i$, the mean up time of a system is equal to $\Lambda^{-1}$. We now consider an optimistic situation where the Man power is large enough such that there is no waiting time in the exponential stations modeling the phases of removing

and reassembling the LRUs. In that case, the mean disassembling and reassembling time, denoted $\gamma^{-1}$, is equal to:

$$ \gamma^{-1} = \sum_{i=1}^{N} \frac{\lambda_i}{\Lambda} (\phi_i + \delta_i), $$

where $\phi_i$ (resp. $\delta_i$) denotes the inverse of the mean disassembling (resp. reassembling) time of LRU of type $i$. So the quantity $(\Lambda^{-1} + \gamma^{-1})$ is the mean sojourn time in the exponential sub-network (in fact the mean sojourn time of the "customer" modeling the behavior of the system in the queuing model). Thus we obtain the corresponding fleet availability thanks to the formula:

$$ A_{\text{fleet}}^{u} = (M - \mathbb{E}[\text{NBS}]) \frac{\Lambda^{-1}}{[\Lambda^{-1} + \gamma^{-1}]} \quad (11) $$

In a less optimistic situation, in deed a pessimistic one, we overload the disassembling and reassembling stations. In order to obtain a pseudo lower bound on the operational availability, we increase the load of the exponential sub-network by considering the LRU's inventories to be infinite. The stations of the non-exponential sub-network disappear and the average disassembling and reassembling time becomes maximal. This new mean time is obtained thanks to classical algorithms for exponential queueing networks. Let $\Gamma^{-1}$ denote this maximal average disassembling and reassembling response time. With this pessimistic situation, $\Gamma^{-1}$ is greater than $\gamma^{-1}$ and the expression of the new approximation of the fleet availability is given by the formula:

$$ A_{\text{fleet}}^{l} = (M - \mathbb{E}[\text{NBS}]) \frac{\Lambda^{-1}}{[\Lambda^{-1} + \Gamma^{-1}]} \quad (12) $$

Of course, when using only the approximate approach for the determination of the probability distribution of NBS, we don't know the exact value of the expectation $\mathbb{E}[\text{NBS}]$ and so, we cannot say that the exact fleet operational availability value is inside the interval $[A_{\text{fleet}}^{l}, A_{\text{fleet}}^{u}]$. Note also that in general, this interval is not relatively large because $\Lambda^{-1}$ is generally much greater than $\gamma^{-1}$ and even than $\Gamma^{-1}$. Let close this section by considering a little example where $M = 10, N = 4, R_1 = 2, R_2 = 3, R_3 = 2, R_4 = 4$ for which we get:

$$ \Lambda^{-1} = 23.8095; \gamma^{-1} = 0.0486; \Gamma^{-1} = 0.0489; \\ \mathbb{E}[\text{NBS}] = 2.6531; A_{\text{fleet}}^{u} = 7.3319; A_{\text{fleet}}^{l} = 7.3318. $$

The exact markovian model give $\mathbb{E}[\text{NBS}] = 2.6506$. So, for this little example, the relative error on $\mathbb{E}[\text{NBS}]$ is $9.4310^{-4}$ while the relative difference between $A_{\text{fleet}}^{u}$ and $A_{\text{fleet}}^{l}$ is $1.3310^{-5}$.
# 6 NUMERICAL TESTS

We first examine the variation of normalized availability ((M-E[NBS])/M) as a function of M when the product $M\lambda_i$ remains constant. To simplify the description of the experimental parameters, the delivery rates $\beta_i(k_i)$ are taken such that $\beta_i(k_i) = k_i\beta$, $i = 1, 2, ..., N$.

We consider a version where each system consists of 4 types of LRU in series ($N = 4$). Such a reduced number of LRU types allows to compare the answer obtained with the original method to the exact solution for a reasonable number of M values.

On Figure 7, M varies from 1 to 30 with the product $M\lambda_i$ remaining equal to 0.75, $i = 1, 2, 3, 4$ and the parameter $\beta$ being equal to 0.6. The availability calculated by this proposed approximate method is represented for the following three situations: i) $R_i = 0, i = 1, 2, 3, 4$, ii) $R_i = 1, i = 1, 2, 3, 4$, iii) $R_i = 2, i = 1, 2, 3, 4$. As expected, the availability increases with the value of $R_i$.

In the cases $R_i = 1$ and $R_i = 2$, the crosses correspond to the exact solutions calculated for $M \le 18$ when $R_i = 1$ and for $M \le 12$ when $R_i = 2$. When the exact values are available, one can notice that the values provided by the proposed method are almost indistinguishable from the exact values. For $M = 10$, the number of states of the reference model is equal to 2,586 and 5,241 respectively when $R_i = 1, i = 1, 2, 3, 4$ and $R_i = 2, i = 1, 2, 3, 4$. Still for $M = 10$, the approximate method gives the result with a relative error of $7.8844 \times 10^{-5}$ when $R_i = 1, i = 1, 2, 3, 4$ and of $3.9457 \times 10^{-5}$ when $R_i = 2, i = 1, 2, 3, 4$.

<figure><img src="image_10.png" /><figcaption>Figure 7: Availability for $R = 0, 1, 2$. Influence of M when the product $M\lambda_i$ remains equal to 0.75. Case $N = 4$ and $\beta = 0.6$.</figcaption></figure>

The second version considered is one where each system consists of 6 types of LRU in series. With this increase in $N$, we have to reduce the range of the exact results. The product

$M\lambda_i$ is here equal to 0.25, $i = 1, 2, 3, 4$. Keeping the other data the same as in the first version, we obtain the results in Figure 8. The exact solutions are computed for $M \le 8$ when $R_i = 1$ and for $M \le 5$ when $R_i = 2$. Again, the values provided by the proposed method are almost indistinguishable from the exact values provided.

<figure><img src="image_11.png" /><figcaption>Figure 8: Availability for $R = 0, 1, 2$. Influence of M when the product $M\lambda_i$ remains equal to 0.25. Case $N = 6$ and $\beta = 0.6$.</figcaption></figure>

For $M = 5$, the number of states of the reference model is equal respectively to 5, 336 and 26, 262 when $R_i = 1, i = 1, 2, 3, 4$ and $R_i = 2, i = 1, 2, 3, 4$. For this value of M, the approximate method gives the result with a relative error of $8.4051 \times 10^{-5}$ when $R_i = 1, i = 1, 2, 3, 4$ and $3.8500 \times 10^{-6}$ when $R_i = 2, i = 1, 2, 3, 4$. Note that the availability obtained with this version is better than in the previous version whereas here $N = 6$ because the product $M\lambda_i$ is here three times weaker.

Figure 9 illustrates the variation of the relative error of the operational availability as a function of $\beta$ for $R_i = 1, 2$ and $3, i = 1, 2, ..., N$, with here $N = 4$ the parameter $\beta$ varying from 0.01 to 0.7.

Recall that the relative error is zero when $R_i = 0$. We can see that if the precision is acceptable when the value of $\beta$ is lower than 0.1, it is remarkably weak when $\beta$ takes values higher than 0.1. Let us add that for $\beta = 0.01$, the operational availability is equal to 0.0703 if $R_i = 1$ and to 0.1142 if $R_i = 3$. For $\beta = 0.1$, the respective values of the operational availability are equal to 0.4779 and 0.7174. While for $\beta = 0.5$, the respective values of the operational availability are equal to 0.8993 and 0.9950. Thus, it is in the useful area of the $\beta$ parameter that the accuracy of the method is the best. In fact, it is possible that an imprecision on the value of an input parameter creates a greater variation of the operational availability than the one induced by the use of the proposed method. For example, if $\beta = 0.101$ instead of $\beta = 0.1$, the
relative variation of the result is 6 10⁻³.

<figure><img src="image_12.png" /><figcaption>Figure 9: Relative error of the approximate availability as a function of $\beta$, for $R = 1, 2, 3$. Case $N = 4, M = 8$ and $\lambda_i = 0.05$.</figcaption></figure>

# 7 CONCLUSION

With this study, we addressed the potential risk of stock shortages of LRU on operational sites of fleet of large systems. To do this, we proposed an original and accurate approximated solution base on non product-form queueing networks giving the expectation of the number of blocked systems due to lake of LRUS. In order to estimate the accuracy of the proposed solution, a Markovian model has been elaborated; this allows us to obtain an evaluation for numerical examples of reasonable size. An additional step allows us to transform with a good accuracy the proactive queues into product-form queues. This finally allows us to determine a new product-form network and to obtain also the operational availability

of the systems with a good accuracy. This study shows in addition that for the studied systems, it is the LRU stock-outs that are mainly responsible for low operational availability. Nevertheless, more investigation needs to be dedicated to deeper analyse the effects of large differences between the respective values of the different LRU allocations on the accuracy of the method.

# References

[1] Qiwei Hu, John E Boylan, Huijing Chen, and Ashraf Labib. Or in spare parts management: A review. *European Journal of Operational Research*, 266(2):395-414, 2018.

[2] Raymond A. Marie. The impact of the unavailability of line replaceable units on the operational availability of system. In *ICIL 2018 - Fourteenth International Conference on Industrial Logistics*, pages 1-6, Beer Sheva, Israel, May 2018.

[3] Raymond A. Marie. Méthode récurrente approchée pour calculer l'indisponibilité opérationnelle induite par les ruptures de stock de LRU. In *21e Colloque National de Maîtrise des risques et sûreté de fonctionnement*, pages 1-9, Reims, France, October 2018.

[4] Raymond A. Marie. Performance Evaluation of proactive Queues. Technical Report RR-9496, Inria Rennes - Bretagne Atlantique & IRISA, January 2023.

[5] Jeffrey P Buzen. Computational algorithms for closed queueing networks with exponential servers. *Communications of the ACM*, 16(9):527-531, 1973.

[6] Raymond A. Marie. "Modélisation par réseaux de files d'attente". Theses, Université de Rennes, November 1978.