Adv Syst Sci Appl 2018; 4:13–38
Published online at http://ijassa.ipu.ru/index.php/ijassa/article/view/624

A Unified Approach to Reliability, Availability, Performability
Analysis Based on Markov Processes with Rewards

Valentina S. Viktorova1∗, Nikolai V. Lubkov1 Armen S. Stepanyants1

1Institute of Control Sciences, Russian Academy of Sciences, Moscow, Russia

Received July 29, 2018; Revised December 17, 2018; Published December 31, 2018

Abstract: This paper discusses a unified approach to reliability, availability and performability
analysis of complex engineering systems. Theoretical basis of this approach is continuous-time
discrete state Markov processes with rewards. From reliability modeling point of view complex
systems are the systems with static and dynamic redundancy, imperfect fault coverage, various
recovery strategies, multilevel operation and varying severity of failure states. We propose a
unified method of calculating the reliability, availability and performability indices based on the
definition of special forms of reward matrix. This method is proved to be effective in calculating
both cumulative and instantaneous measures in steady-state and transient cases. We describe
special analytical software which implements suggested method. We demonstrate the flexibility
of the proposed method and software by analyzing multilevel process unit with protection and
demand-based warm standby system.

Keywords: Availability, Markov process, Markov reward model, multilevel engineering
systems, performability, reliability.

1. INTRODUCTION

At the present state of development of reliability theory the models of analysis can be
subdivided into two categories (classes): the dynamic models and the static models. All
changes of system states in the class of dynamic models are considered as the processes
developing in time. In the class of the static models, states of the system are determined by
the sets of states of system elements at the time moment t.

Markov [1–4], semi-Markov processes, asymptotic methods of the renewal theory and
regenerative processes [5–7], Monte Carlo simulation techniques [8] are used within the
framework of dynamic models. Dynamic models allow to calculate all the main dependability
measures both for repairable and non-repairable systems. These measures are: instantaneous
indices (e.g. availability at the time instant); interval indices (e.g. reliability during the
time interval); time-independent stationary indicators (e.g. mean time between failures).
Known drawbacks inherent in Markov models are the size and stiffness of transient solution
of a system of Kolmogorov-Chapman equations (or ill-conditionality for stationary case).
Possible ways to solve these problems are discussed in [9–11]. Monte Carlo simulation of
modern high reliable systems may require large amount of simulations to obtain results with
desired accuracy. These shortcomings can be eliminated with the help of special acceleration
techniques. There is no limitation to type of distribution of waiting times between the changes
in the class of Monte Carlo models. But the problem of creating a universal Monte Carlo

∗Corresponding author: vsviktorova@gmail.com



14 V.S. VIKTOROVA, N.V. LUBKOV, A.S. STEPANYANTS

model describing the complex reliability behavior of different kind of systems remains
unsolved.

Static models use two main groups of methods – combinatorial-probabilistic [3] and
logical-probabilistic [12]. Combinatorial-probabilistic methods use the basic formulas of
combinatorics and probability theory (probability of the sum and product of events, the
formula of total probability). These formulas are used mainly for serial-parallel and ”m out
of n” redundant schemes. Logical-probabilistic methods, which are based on the construction
of a Boolean function relating the state of the system with the states of its elements. The
resulting Boolean function is transformed to a form that allows to replace logical variables by
corresponding probabilities. Classic failure trees and reliability block diagrams are the main
methodology of logical-probabilistic methods [13, 14]. Classic static models for the case
of repairable systems allow us to calculate only the differential (instantaneous) reliability
indicators determined at the time instant t (e.g. availability, failure frequency, average
efficiency at the time instant t).

In the last decade, combined approaches to reliability models construction are successfully
developing. These approaches are based on the implementation of dynamic properties into
static models. Dynamic fault trees are the most popular practical implementation of these
approaches [15–19].

Modern complex engineering systems with high demands on reliability are characterized
by various features, which include:

• multiple levels of operation efficiency (e.g., performance) and the graceful degradation
of performance in the event of failures [20–23];
• a variety of redundancy implementation techniques (hardware, software, functional,

temporal redundancy) [24–26], types of redundant schemes (parallel, standby, hybrid),
modes of reserve components (cold and warm standby, shared) [3, 4];
• smart recovery strategies, restricted repair resources [27–29];
• imperfect built-in test equipment leading to the presence of undetected failures and false

alarms [30, 31];
• implementation of special multiphase error handling procedures with classification on

permanent and transient faults (mainly for fault tolerant computers) [32–34];
• the possibility of the occurrence of mutually exclusive failures that, with a certain

multiplicity and sequence, can lead to different consequences at the system level [13,14].

Markov modeling allows taking into account the above-mentioned features of the reliability
behavior of systems and calculating most part indicators of reliability, safety, efficiency.
At the design stage, we don’t have, as a rule, objective information about the distribution
functions of random time variables of the model. Therefore, we believe that the assumption
of exponentiality is entirely permissible. The problem of large size of the Markov model can
be solved by decomposition and aggregation of the model parts and by automation of the
model construction.

In this paper, we consider the Markov reward models as common tools for analyzing
complex reliability behavior of technical systems. We propose the extension of Markov
reward modeling for the reliability, abailability, performability (RAP) analysis of the systems
at the design stages. We suggest the method of calculating reliability, availability and
performability indicators based on a special definition of the reward matrix. This method
is quite efficient and well suited for software implementation. It makes calculation both
cumulative and point indices possible only by solving systems of differential or algebraic
equations and does not require the usage of additional numerical integration operations.

The remainder of this paper is organized as follows. A brief summary of the mathematical
foundations of Markov reward model is given in Section II. In Section III, firstly we give
definitions of interval, point, steady-state RAP measures. Next, we present the method for
for calculation of these measures in the framework of Markov reward model. Section IV
presents analytical software implementing the proposed method. The flexibility of described

Copyright c© 2018 ASSA. Adv Syst Sci Appl (2018)



A UNIFIED APPROACH TO RAP ANALYSIS BASED ON MARKOV PROCESSES WITH REWARDS 15

techniques and software is illustrated by investigations of two practical examples in Section
V. In Section VI, we give some concluding remarks.

2. MATHEMATICAL FOUNDATIONS

A system state at the time moment t can be described by stochastic process z(t) of transitions
in system state space Ω. We assume that the state space Ω is discrete, finite and of sizeN . Ω is
defined by states of system components (failed or operational). We define a quality function Φ
on z(t). Changing system states leads to a change in the quality function Φ. z(t) is stochastic
process, therefore, Φ(z(t)) = Φ̃(t) are also stochastic functions. Due to the discrete nature
of Ω, the functions Φ̃(t) are the step functions, whose values are quality levels (for example,
system performability) corresponding to the states. G are indices of system performability. G
can be defined as a measure on the trajectories of the Φ̃(t) . General representation of G takes
the form

G = M{F [Φ(z(t)]}, (2.1)

where M represents mathematical expectation; F is a functional, which is determined by
the type of performability index.

Thus, the model of performability contains three components Ω, z(t),Φ. It is natural and
obvious to present the model of performability by means of transition graph (state graph).
The elements of the state space correspond to the vertices. Possible transitions in the state
space correspond to the arcs of the graph. Φ is defined as a function of the system state si.
This function takes the value of wi corresponding to the performance level of the system in
this state. That is

Φ(si) = wi. (2.2)

In addition, this model can display effects occurring at the transition of the system from
one state to another. This is accomplished by determining reward matrix

W = [wij]
n
n , (2.3)

where wij – is an effect on the system that arises at the transition from state i to state j;
wii = wi.

Random transition process z(t) is most simply set in the case where the value of time of
stay in each state have an exponential distribution, i.e. when the process is a continuous-time
Markov chain (CTMC).

The process z(t) is determined by transition rates, which correspond to arcs of the Markov
graph. Similarly to the case of matrix W , we can create a matrix of transition rates

Λ = [λij]
n
n , (2.4)

where λij – is the rate of transition from state i to state j; λii = −
∑
j,j 6=i

λij .

2.1. Construction of RAP indices
Thus, the Markov model of performability is defined by two matrices Λ and W , i.e.
{R, z(t),Φ} ∼ {Λ,W}. Model {Λ,W} is called Markov reward model (MRM). Markov
reward model describes the behavior of systems with exponential distribution of elements
failure and repair time. States of the model represent the states of the system, which
correspond to different sets of failed and operable items.

If criterion of system failure is specified, then expression (2.1) defines also the indices of
reliability and availability.

Copyright c© 2018 ASSA. Adv Syst Sci Appl (2018)



16 V.S. VIKTOROVA, N.V. LUBKOV, A.S. STEPANYANTS

Let failure criterion is {Φ̃(t) < Φcr} ⇒failure, where Φ̃(t) - is the current value of
the quality function; Φcr - is limit value of the quality, which is allowed by operating
requirements. Then, by setting a particular type of the functional F the following indices
can be determined from Equation (2.1):

• Reliability over the time interval (0, t)R(t):

R(t) = P{Φ̃(τ) > Φcr, ∀τ ∈ (0, t)} under F [·] =

{
0, if ∃τ ∈ (0, t) : Φ̃(τ) < Φcr

1, otherwise
(2.5)

• Mean time to first failure (MTTF) TFF :

TFF =

∞∫
0

R(t)dt underF [·] = min{τ : Φ̃(τ) < Φcr} (2.6)

• Availability A(t) at the time instant t:

A(t) = P{Φ̃(τ) ≥ Φcr} under F [·] =

{
0, Φ̃(τ) < Φcr

1, Φ̃(τ) ≥ Φcr

(2.7)

• Probability of being at the ith level of performance at the time instant t:

P (t) = P{Φ̃(t) = Φi} under F [·] =

{
1, Φ̃(t) = Φi

0, otherwise
(2.8)

• Average accumulated time spent at the ith level of performance:

under F [·] =

∫
t:Φ̃(t)=Φi

dt (2.9)

• Average number of transitions from the ith to the jth performance level during the time
interval (0, t):

under F [·] =
∑

l : Φ̃(tl−) = Φi

Φ̃(tl+) = Φj

1, (2.10)

where tl - the time moment of transition from ith level to jth performance level.
• Time-averaged accumulated reward for the time interval (0, t):

F [·] =
1

t

t∫
0

Φ̃(τ)dτ (2.11)

2.2. Description of MRM Basic Relations
The consept of reward in Markov reward model is a generalized concept. This can be any
effect, loss, costs. According to the provisions of previous section we define the model state
space as a discrete set {si}, i = 1, N . Also we define reward matrix. The elements of the
reward matrix are interpreted as follows: wii is a reward at time unit while the system is

Copyright c© 2018 ASSA. Adv Syst Sci Appl (2018)



A UNIFIED APPROACH TO RAP ANALYSIS BASED ON MARKOV PROCESSES WITH REWARDS 17

in state si (reward rate); wij is impulse reward received in the system under the transition
from state i into state j; N - the number of states of the system. We define the total (or
accumulated [20]) reward obtained by the system for a given time T as the functional
F [Φ(z(t))]:

G̃(T ) =
N∑
i=1

τiwii +
∑
i,j

nijwij, (2.12)

where τi is the residence time in state si during the time interval (0, T ); nij is the number
of system transitions from state i into state j during the time interval (0, T ).

Accumulated reward G̃(T ) is a random variable, since τi and nij are random values.
Average accumulated reward is defined as the mathematical expectation of G̃(T )

Y (T ) = M{G̃(T )}. (2.13)

According to (2.1), Y (T ) is an estimate of the performability index value. If we consider
that members of (2.12) τi and nij depend only on the specific implementation of process z(t),
it is obvious that the type of index is completely determined by the type of matrix W . The
system of Howard differential equations describes the behavior of the average accumulated
reward:

dyi(t)/dt = wii +
∑
j,j 6=i

λijwij +
∑
j

λijyj(t); i, j = 1, N, yi(0) = 0, (2.14)

where yi(t) is the average accumulated reward for the time interval (0, t), given that the initial
state is si(z(0) = si); λij - the rate of transition from state si into state sj .

The matrix form of equation (2.14) is

dY (t)/dt = ΛY (t) +R; Y (0) = 0, (2.15)

where Λ = ‖λij‖ - matrix of transition rates; R is a column vector of constant terms:

ri = wii +
∑
j,j 6=i

λijwij; i, j = 1, N. (2.16)

We assume that for absorbing states sr of the Markov process valuewrr = 0. The absorbing
state is a state that once entered, can not be left. This assumption is quite natural, since the
absorbing states are usually identified with the non-operational states of the system.

3. METHOD OF CALCULATING RAP INDICES BASED ON MARKOV REWARD
MODEL

The expression (2.1) of expected reward will determine the RAP indices in accordance with
(2.5) – (2.11), if the value of the reward for a particular process realization (2.12) coincides
with the functional F[.] in (2.5) – (2.11). It requires to select the appropriate reward matrix
W and optionally to correct the transition rate matrix Λ. Correction of matrix Λ is necessary
in case when the index reflect operation of the system till some event e, for example, till
the system failure. Since all performability indices in the model are interpreted through the
reward, it is obvious that reward should be dismissed after event e has occurred. Formally,
this can be achieved by identification of the event e with a transition to an absorbing state.
The reward associated with this state is assumed to be zero. The correction of matrix Λ is that
we treat the states corresponding to e as absorbing states. To do this correction, it is enough

Copyright c© 2018 ASSA. Adv Syst Sci Appl (2018)



18 V.S. VIKTOROVA, N.V. LUBKOV, A.S. STEPANYANTS

to equate to zero the elements of rows of matrix Λ corresponding to these states. To carry
out further calculations, we partition the set of model states Ω into Ωg, the set of operational
system states, and Ωf , the set of failed states.

All RAP measures can be subdivided into three groups:

1. Interval measures
These metrics are dependent on time and are determined on finite time interval [0, T ],
for example, reliability R(0, T ).

2. Point measures
These metrics are dependent on time and are determined in time point t, for example,
availability A(t).

3. Steady-state measures
These metrics are independent on time and are determined on an unlimited time interval
(0,∞), for example, MTTF.

Calculation of steady-state measures is reduced to finding the stationary solutions of the
system (2.15), i.e., to solving a system of linear algebraic equations:

ΛY +R = 0. (3.1)

We mention that all indices calculated by the expressions (2.15) and (3.1) are vectors. The ith
element of this vector is the estimated value of the index for the initial state si.

3.1. Unreliability
System unreliability over the time interval (0, t)(Q(t)) can be calculated if we define elements
of matrix W as wij = 0 ∀ i, j ∈ Ωg;wij = 1 ∀ i ∈ Ωg and j ∈ Ωf :

wij =

{
1, if i ∈ Ωg, j ∈ Ωf

0, otherwise . (3.2)

Besides we must treat states in Ωf as absorbing states. Namely we should to delete in
Markov graph all arcs leading from states in Ωf to states in Ωg.

The calculation is performed by (2.15).
Complementary index is reliability: R(t) = 1−Q(t).

3.2. Mean Time to First Failure
Mean time to first failure (TFF ) can be calculated if we define diagonal elements of matrix
W corresponding to the operational states as wii = 1 ∀ i ∈ Ωg.

wij =

{
1, if i = j, i, j ∈ Ωg

0, otherwise . (3.3)

Besides we must treat states in Ωf as absorbing states.
Calculation is performed by (3.1). For the given type of reward matrix expression (3.1) is

reduced to

Λ∗TFF =

 −1
...
−1

 , (3.4)

where Λ∗ − k × k matrix derived from matrix Λ by deleting rows and columns corresponding
to nonoperational states; k - total number of operational states.

Copyright c© 2018 ASSA. Adv Syst Sci Appl (2018)



A UNIFIED APPROACH TO RAP ANALYSIS BASED ON MARKOV PROCESSES WITH REWARDS 19

3.3. Standard Deviation of Random Time to First Failure
Standard deviation (σ) and variance (σ2) of mean time to first failure can be calculated in

accordance with the expression σ =
√
M(t2)−M2(t) =

√
T

(2)
FF − T 2

FF .

To calculate vector T (2)
FF , we have to find the solution of algebraic equations

Λ∗T
(2)
FF =

 −2TFF1
...

−2TFFk

 , (3.5)

where TFFi is mean time to first failure given the process starts from the ith state.

3.4. Average accumulated time spent in the operational states
Average accumulated time spent in the operational states during the time interval
(0, t) (T∑(t)) can be calculated if we define diagonal elements of matrix W corresponding
to the operational states as wii = 1 ∀ i ∈ Ωg.

wij =

{
1, if i = j; i, j ∈ Ωg

0, otherwise . (3.6)

Here, we should not artificially create absorbing states.
Calculation is performed by the Equation (2.15).
We can also use this index for analysis of multilevel systems. In this case we should

partition the set of system’s states Ω into m sets Ω1, . . .Ωr, · · · ,Ωm, where Ωr - the set of
operational system states with rth performance level. Average accumulated time spent in Ωr

during the time interval (0, t) (T r∑) can be calculated at wii = 1 ∀ i ∈ Ωr.

3.5. Availability
Availability A(t) is the probability that the system is operational at the time instant t. A(t)
is the point measure. Taking into account that A(t) = dTΣ(t)/dt we can calculate A(t) after
calculation of TΣ(t) as follows:

A(t) = ΛTΣ +R. (3.7)

The formation of reward matrix is carried out in accordance with (3.6).
Unavailability U(t) = 1− A(t).

3.6. Time-averaged availability
General expression of time-averaged availability over the time interval (0, t) (Aav(t)) has the
form

Aav(t) =
1

t

t∫
0

A(τ) dτ. (3.8)

For calculation of this indicator, it is necessary to take the reward matrix, in which the
diagonal elements of the operational system states, are equal to 1/t and all other elements
are zero:

wij =

{
1/t, if i = j; i, j ∈ Ωg

0, otherwise . (3.9)

Artificial creation of absorbing states is not required. Calculations are carried out by (2.15).

Copyright c© 2018 ASSA. Adv Syst Sci Appl (2018)



20 V.S. VIKTOROVA, N.V. LUBKOV, A.S. STEPANYANTS

3.7. Average failures number
Average (expected) failures number during the time interval (0, t)(N(t)) can be calculated if
we define elements of matrixW columns corresponding to the failed states aswij = 1 ∀ i 6= j.
Values of other elements of matrix W are zero:

wij =

{
1, if i ∈ Ωg, j ∈ Ωf

0, otherwise . (3.10)

N(0, t) is failure measure for repairable systems so we should not artificially create absorbing
states. Calculation is performed by (2.15).

3.8. Failure frequency
Failure frequency (ω(t)) is a differential index with respect to N(t). After calculating of
vector N(t) for each initial state, we put it into the right side of (2.15). The resulting vector
dY (t)/dt will be the vector (ω(t)):

ω(t) = ΛN(t) +

 λ1
Σ
...
λNΣ

 , (3.11)

where λiΣ is total transition rate from ith operational state to failure states. For the failure
state, it is zero.

3.9. Average accumulated reward
Average accumulated reward over the time interval (0, t)(ET (t)) integrates reward (loss) of
the system during the time interval (0, t) proportionally to the time of stay in the states and the
quantity of transitions between the states. Reward matrix W has the most general form. The
elements wii, wij are real rewards (losses) measured in terms of units of system performance.
Calculations are carried out by (2.15).

3.10. Average reward
Average reward at the time instant t(Eav(t)) is a point measure. This indicator characterizes
the multi-level systems. Therefore, system states set is subdivided into classes, in accordance
with level of system performance. States of the ith class are characterized by performance
(Ei) per time unit (for example, throughput). The first step of Eav(t) calculation is performed
by the Equation (2.15). As a result of the first calculation step, we get the vector of the average
accumulated reward ET (t). Eav(t) is a differential index with respect to ET (t). So the second
calculation step is to put the vector ET (t) into right side of equations (2.15):

Eav(t) = ΛET (t) +R. (3.12)

3.11. Performability ratio
This index is the ratio of the average accumulated reward to the nominal reward during the
time interval (0, t):

Kp(t) =
ET (t)

En(t)
, (3.13)

where En(t) - nominal (maximum) accumulated reward generated by the system in the
absence of failures. En(t) = Emaxt. Usually Emax = E1.

Copyright c© 2018 ASSA. Adv Syst Sci Appl (2018)



A UNIFIED APPROACH TO RAP ANALYSIS BASED ON MARKOV PROCESSES WITH REWARDS 21

4. ANALYTICAL SOFTWARE BASED ON MARKOV REWARD MODEL

The above uniform approach to calculating performance, availability and reliability measures
was used by us in the development of analytical software. The software is operated under
Windows OS and is written in c] programing language. It has an advanced graphical user
interface (Fig.4.1) and uses OLE automation for generation reports in MS Word formats.

This software allows you to take into account all features of reliability behavior of complex
technical systems, mentioned in sectoion 1:

• multiple levels of operation performance and the graceful degradation of performance
in the event of failures;
• different types of redundancy implementation techniques;
• variety of redundant schemes and backup components load levels;
• recovery strategies and restrictions on repair resources;
• imperfect built-in test;
• possibility of transient faults occurrence;
• common cause failures;
• and so on.

The user of this software can calculate all the RAP indicators described in Section 3:

• Reliability (Unreliability) for a given time interval (0, t) ;
• Mean and Standard Deviation of random time to first failure;
• Average and Accumulated Time spent in the selected subset of states during the time

interval (0, t);
• Availability (Unavailability) at given time instant t;
• Average Failures Number for a given time interval (0, t);
• Failure Frequency at given time instant t;
• Average Reward at given time instant t;
• Average Accumulated Reward for a given time interval (0, t);
• Performability Ratio for a given time interval (0, t).

The software execution process is divided into three steps: creation of the MRM, setting up
of the MRM, numerical solution of the system of Howard equations. The composition and
interrelationships of these processes are shown in Fig.4.2-Fig.4.4.

MRM creation is executed by the user in interactive mode with the help of built-in graphical
editor. The created model can be saved in an external XML file. A saved model can be
reloaded from the XML file (Fig.4.2).

The model’s setting up is performed in the Matrix Configuration Module in accordance
with equations of section 3. Setting up operations are performed automatically after the user
chooses the metric for calculation (Fig.4.3).

Calculation of values of user-selected RAP indices is made in the numerical solution
module (Fig.4.4). In this module, numerical procedures for solving systems of differential and
algebraic Howard equations are realized. Different approaches and methods of computational
mathematics with regard to their software implementation outlined in [35 - 37].

Detailed consideration of numerical evaluation of Markov models transient behavior was
presented in [9], [10]. Three numerical techniques for finding the transient solution of
large Markov models were examined. They are - uniformization, an explicit differential
equation solution method (Runge-Kutta) and special stable implicit method (TR-BDF2). For
the software implementation of numerical solution of the MRM we used efficient method
based on evaluation of matrix exponential at small step. This method had been proposed and
described in [38].

For a stationary case, a numerical solution is obtained by solving the system of algebraic
equations by the Gauss method or the Rotation method [35,36]. The Rotation method works
well for bad conditioned systems, which are often generated by RAP analysis models.

Copyright c© 2018 ASSA. Adv Syst Sci Appl (2018)



22 V.S. VIKTOROVA, N.V. LUBKOV, A.S. STEPANYANTS

Fig. 4.1. Screenshot of main form of the analytical software.

Copyright c© 2018 ASSA. Adv Syst Sci Appl (2018)



A UNIFIED APPROACH TO RAP ANALYSIS BASED ON MARKOV PROCESSES WITH REWARDS 23

Fig. 4.2. Execution process of the analytical software (Step 1).

Copyright c© 2018 ASSA. Adv Syst Sci Appl (2018)



24 V.S. VIKTOROVA, N.V. LUBKOV, A.S. STEPANYANTS

Fig. 4.3. Execution process of the analytical software (Step 2).

Copyright c© 2018 ASSA. Adv Syst Sci Appl (2018)



A UNIFIED APPROACH TO RAP ANALYSIS BASED ON MARKOV PROCESSES WITH REWARDS 25

Fig. 4.4. Execution process of the analytical software (Step 3).Copyright c© 2018 ASSA. Adv Syst Sci Appl (2018)



26 V.S. VIKTOROVA, N.V. LUBKOV, A.S. STEPANYANTS

5. CASE STUDIES

To demonstrate the effectiveness of the software implemented MRM, we present two
examples of reliability, availability, performability evaluation of complex systems. In first
example we investigate technological object with protection system. The second one is
demand-based warm standby system (DB-WSS), described in article [23].

5.1. Markov reward model of multilevel process unit
We consider process unit consisting of technological object (TO) with protection system (PS).
The protection system includes a diagnostic unit (DU) and the actuator (A) (Fig.5.5).

Fig. 5.5. Technological object with protection.

Protection systems are designed to generate control actions to the protected object in
order to prevent transition from process equipment failures to accident. Control actions may
be different, such as change of the operating mode, reduction in productivity, emergency
shutdown of faulty elements and elements belonging to the same processing chain. These
actions prevent the development of accident. Here, we consider the case of technological
process shutdown. However, proposed approach is also suitable for the control actions that
lead to poor performance or mode change.

The main function of protection system is performed sporadically at the time instants of the
occurrence of the object failures, so that the protection system operates in standby temporary
mode. Let us define rewards of operational states and losses of failed states. In this example,
all possible system states can be subdivided into four classes. The first class includes such
technical states in which the object is functioning normally and earns specific rewards per
unit of time of stay in these states (e.g., these states correspond to the nominal performance
of technological object). Let us note that this group includes states with latent failures such
as failure of diagnostic unit. The second class consists of the states of accident-free shutdown
of the technological object. The losses associated with these states are only related to the
downtime of the object. The reward in these states is either zero or negative, if the idle leads to
additional losses per unit time. The third and fourth classes of states are catastrophic failures
of the object. The transition to this class of states from the states of the first class brings a one-
time damage (negative reward) associated with the occurrence of an accident (death of people,
equipment breakdowns, emissions into the atmosphere, etc.). In this model, we consider two
types of catastrophic failures differing in the severity of the consequences – accident I and
accident II. Thus, the normal functioning of the technological object is accompanied by a
linear increase in the accumulated reward in proportion to the time spent in the first class
states. Idle time leads to the preservation of the achieved level of accumulated reward (at
zero values of the reward rate in each state) or to its descent (at negative values of the reward
rates) in proportion to the time spent in states of second, third and fourth classes. When
impulse rewards associated with the transitions between the states are not zero, there is an
abrupt change (more often a decrease) in the accumulated reward. Negative impulse rewards
are due to the costs of recovery from failures or accidents, purchase of equipment, payment
of fines or insurance, etc.

Copyright c© 2018 ASSA. Adv Syst Sci Appl (2018)



A UNIFIED APPROACH TO RAP ANALYSIS BASED ON MARKOV PROCESSES WITH REWARDS 27

We suppose that all accidents occur only due to failures of the technological equipment.
Let protection system in the event of critical situation instantly shut down the technological
process by making the necessary equipment control (for example, switching off). Let the
protection system processes critical situation with coverage β(0 ≤ β ≤ 1). Failures of the
protection system which occur during interval of normal TO operation can lead to different
consequences. Let us single out PS failures of two types: latent failures (no operation) and
explicit failures (false alarm). Latent failures of PS are manifested in the form of absence
of protection actions in the event of critical situation which entails an accident. False alarms
lead to the undesired protection actions in the absence of TO failures.

Parameters of the Markov reward model of the technological object with protection system
are:
wij - impulse losses (negative impulse reward) under the transition from state i into state j;
wii - reward rate (negative or positive) in state i;
β - coverage probability for TO, determined as conditional probability of processing of TO

catastrophic failure by a protection system, provided that the failure has occurred;
αDU - fraction of the protection system latent failures such as non operation of DU;
αA - fraction of the protection system latent failures such as non triggering of actuator;
1− αDU - fraction of the protection system explicit failures such as DU false alarm;
1− αA - fraction of the protection system explicit failures such as actuator false triggering;
ηI - fraction of technological object catastrophic failures of the first kind (accident I);
ηII - fraction of technological object catastrophic failures of the second kind (accident II);
λTO, λDU , λA – failure rates of technological object, diagnostic unit, actuator respectively;
µ - repair rate of technological object after shutdown;
µa - repair rate of technological object after falling into accident II state.
The transition graph of Markov reward model of the technological object with protection

is shown in 5.6.
We partition the set of the MRM states into four subsets:

• normal operation (OP) (states 1,3,4)
• accident-free shutdown (SD) (states 2,6)
• accident I (AI) (state 5)
• accident II (AII) (state 7)

Transition rates matrix Λ of the model is

Λ =


λ11 λ12 λ13 λ14 λ15 0 λ17

λ21 λ22 0 0 0 0 0
0 λ32 λ33 0 λ35 0 λ37

0 0 0 λ44 λ45 λ46 λ47

λ51 0 0 0 λ55 0 0
λ61 0 0 0 0 λ66 0
0 0 0 0 0 0 λ77

 . (5.1)

The diagonal elements of the matrix λii are given by −
∑
j,j 6=i

λij .

Transition rates between states of the graph are determined based on the model parameters
as follows:

λ12 = [1− (1− β)(ηI + ηII)]λTO + (1− αDU)λDU + (1− αA)λA;
λ13 = αDUλDU ;
λ14 = αAλA;
λ15 = (1− β)ηIλTO;
λ17 = (1− β)ηIIλTO;
λ21 = λ61 = µ;
λ32 = (1− ηI − ηII)λTO + (1− αA)λA;

Copyright c© 2018 ASSA. Adv Syst Sci Appl (2018)



28 V.S. VIKTOROVA, N.V. LUBKOV, A.S. STEPANYANTS

λ35 = λ45 = ηIλTO;
λ37 = λ47 = ηIIλTO;
λ46 = (1− ηI − ηII)λTO;
λ51 = µa;

Fig. 5.6. Markov graph of TO with protection.

By specifying the form of reward matrix, it is possible to obtain equations for calculating
various indices. We want to determine the cumulative effect of operation of technological
object taking into account the positive and negative impact of the protection system.
Therefore, it is advisable to calculate the following indices:

1. Q(t) - unreliability over (0, t).
Q(t) can be calculated from equation (2.15).
To calculate Q(t), we must equate to zero λ21, λ51, λ61 and define W as

W =


0 1 0 0 1 0 1
0 0 0 0 0 0 0
0 1 0 0 1 0 1
0 0 0 0 1 1 1
0 0 0 0 0 0 0
0 0 0 0 0 0 0
0 0 0 0 0 0 0

 . (5.2)

2. TΣ(t) -average accumulated time spent in the operational states over the time interval
(0, t).
TΣ(t) is determined by (2.15), setting R as

Copyright c© 2018 ASSA. Adv Syst Sci Appl (2018)



A UNIFIED APPROACH TO RAP ANALYSIS BASED ON MARKOV PROCESSES WITH REWARDS 29

R = [ 1 0 1 1 0 0 0 ]T . (5.3)
3. A(t) - system availability at the time instant t.
A(t) is calculated by (3.7) after calculation of vector TΣ(t).

4. TFF - mean time to first failure.
TFF is determined by (3.4), where matrix Λ∗ is

Λ∗ =

[
λ11 λ13 λ14

0 λ33 0
0 0 λ44

]
. (5.4)

5. σ - standard deviation of random time to first failure.
σ is obtained from Equation (3.5), where vector of free terms is
[ −2TFF1 −2TFF2 −2TFF3 ]T .

6. TA - mean time to first catastrophic failure (accident I or II).
TA is determined by Equation (3.4), where matrix Λ∗ is 5× 5 matrix derived from matrix
(5.1) by deleting rows and columns corresponding to accident I and accident II states
(states 5 and 7).

7. N(t) - expected failures number.
N(t) can be calculated from (2.15) where matrix W takes the form (5.2).

8. ω(t) - failure frequency.
ω(t) can be calculated from (3.11) after N(t) calculation and setting free terms vector
as [ λ12 + λ15 + λ17, 0, λ32 + λ35 + λ37, λ45 + λ46 + λ47, 0, 0, 0 ]T .

9. NA(t) - expected catastrophic failures number.
NA(t) can be calculated from (2.15) where matrix W has unit elements only in the
transitions to catastrophic failures (states 5 and 7). Similarly, we can find separately
NAI(t) and NAII(t).

10. TΣSD(t) –average accumulated time spent in the states of technological object shutdown
over time interval (0, t). TΣSD(t) is determined by Equation (2.15), setting R as
[ 0, 1, 0, 0, 0, 1, 0 ]T .

11. NSD(t) - average number of transitions to states of shutdown.
NSD(t) can be calculated from Equation (2.15) where matrix W has unit elements only
in the transitions to shutdown states (states 2 and 6).

12. ET (t) - average accumulated reward over the time interval (0, t).
Calculation of ET (t) is performed by (2.15).
Matrix of reward W is

W =


w11 w12 w13 w14 w15 0 w17

w21 w22 0 0 0 0 0
0 w32 w33 0 w35 0 w37

0 0 0 w44 w45 w46 w47

w51 0 0 0 w55 0 0
w61 0 0 0 0 w66 0
0 0 0 0 0 0 w77

 . (5.5)

13. Eav(t) - average reward at the time instant t.
In accordance with (3.12) Eav(t) is determined by solving the system of algebraic
equations

 Eav1(t)
...

Eav7(t)

 = Λ

 ET1(t)
...

ET7(t)

+

 w11
...

w77(t)

+


∑
j 6=1

λ1jw1j

...∑
j 6=7

λ7jw7j

 . (5.6)

Copyright c© 2018 ASSA. Adv Syst Sci Appl (2018)



30 V.S. VIKTOROVA, N.V. LUBKOV, A.S. STEPANYANTS

The numerical results of above indices calculation for different values of parameter β are
summarized in table 5.1. All time-dependent metrics were calculated on one year period
(t = 8760hours). We note once again that the solution, obtained by the MRM, is a vector
whose ith component corresponds to the start of the Markov process from state i. The table
and diagrams include the first components of solution vectors, i.e the system starts from
completely operational state S1.

The functions of ET (t), Eav(t) for different values of catastrophic failure coverage
parameter β are shown in Fig.5.7 and Fig.5.8 respectively.

Calculation of indices and plotting are performed with the following system parameters:
λTO = 0.0001(/hour);λDU = 0.00001(/hour);λA = 0.00005(/hour).
µ = 0.1(/hour);µa = 0.004(/hour).
αDU = 0.5;αA = 0.9.
ηI = 0.4; ηII = 0.05.
Numerical values of impulse rewards are
w15 = −105;w17 = −108;w12 = w13 = w14 = 0;
w21 = 0;w32 = 0;w35 = −105;w37 = −108;
w45 = −105;w46 = 0;w47 = −108;
w51 = w61 = 0.
Numerical values of reward rates are
w11 = 15(/hour);w33 = 10(/hour);w44 = 10(/hour);
w22 = w55 = w66 = w77 = 0(/hour).

Fig. 5.7. TO average accumulated reward ET (t) versus time t for different coverage values β.

The first family of curves (Fig.5.7) shows the significant dependence of the average
accumulated reward on the catastrophic failure coverage β. For the model parameters under
consideration β should be greater than 0,9. The second family of curves(Fig.5.8) shows a
sharp decline in the growth of the average accumulated reward in time for all coverage values.
It occurs due to the high rate of PS latent failures. To prevent such a recession, it is advisable

Copyright c© 2018 ASSA. Adv Syst Sci Appl (2018)



A UNIFIED APPROACH TO RAP ANALYSIS BASED ON MARKOV PROCESSES WITH REWARDS 31

Fig. 5.8. TO average reward Eav(t) versus time t for different coverage values β.

Table 5.1. Calculation results over one year period.

Catastrophic failure coverage β
Index 0.8 0.85 0.9 0.95 1
Q(t) 0.612 0.612 0.612 0.612 0.612
TΣ(t) 8670.6 8682.4 8694.3 8706.2 8718.1
A(t) 0.9814 0.9836 0.9858 0.988 0.99
TFF 9360.1 9360.1 9360.1 9360.2 9360.1
σ 9555.6 9555.6 9555.6 9555.6 9555.6
TA 4774.7 5143.4 5573.8 6082.8 6694.1
N(t) 0.942 0.943 0.944 0.946 0.947
NA(t) 0.124 0.108 0.091 0.075 0.058
NAI(t) 0.11 0.096 0.081 0.0663 0.0516
NAII(t) 0.014 0.012 0.01 0.0083 0.0064
TΣSD(t) 8.16 8.34 8.52 8.7 8.88
NSD(t) 0.817 0.835 0.853 0.871 0.889
ET (t) -15489 3146 21834 40575 59368
Eav(t) -6.068 -4.202 -2.33 -0.443 1.45

to carry out preventive maintenance of the process unit (TO+PS), during which latent failures
of the protection system are detected.

5.2. Reliability and Performability Investigation of DB-WSS
We take demand-based warm standby system (DB-WSS), described in Section IV.A of article
[23] as the second object of study. This DB-WSS is composed of two components A1 and A3

(we keep the numbering made by the authors of [23]).
The components A1 and A3 can be in four degradation states, which differ in capacity. The

components capacity levels are 5, 4, 2, 0 and 5,4,1,0, respectively. Component A1 is online.
Component A3 can be online or in warm standby state depends on the system demand. We

Copyright c© 2018 ASSA. Adv Syst Sci Appl (2018)



32 V.S. VIKTOROVA, N.V. LUBKOV, A.S. STEPANYANTS

consider scenario B from [23], when demand determined as 5. So initial state of A3 is warm
standby. The transition time distributions for the degradation processes for each component
are exponential distributions with parameters λ1 = 1/100, λw3 = 1/300, λ3 = λo3 = 1/150
(/day), where (w) indicates warm standby state and (o) – operating online state.

First, we will construct a Markov model completely corresponding to the degradation
process described in [23]. Namely, we will construct a model of a binary-state system
consisting of multistate components. The transition Markov graph of this DB-WSS is shown
in Fig.5.9. The states S1 and S10 are operational states of the system that differ in the capacity
of components. The state (j, k) corresponds to the jth capacity level of element A1(CA1

j ) and
the kth capacity level of element A3(CA3

k ). State S11 is failed system state. CA1
j + CA3

k ≥ 5

for operational states; CA1
j + CA3

k < 5 for failed state.

Fig. 5.9. Markov graph of binary-state DB-WSS.

For the first case, we calculate the following reliability indicators:

1. R(t) reliability during (0, t).
Unreliability Q(t) can be calculated from equation (2.15), where vector R =

[ 0 0 0 0 0 0 λ3 λ1 + λ3 λ1 + λ3 λ1 0 ]T .
Reliability R(t) = 1−Q(t). Table 5.2 shows the results of our calculation of R(t) and
values of R(t) taken from [23]. For the model under study (Fig.5.9), an exact analytical
solution was obtained (see Appendix). The results of R(t) calculations using (A.2),
(A.3) are given in the last column of Table 5.2.

2. TFF -mean time to first failure.

Copyright c© 2018 ASSA. Adv Syst Sci Appl (2018)



A UNIFIED APPROACH TO RAP ANALYSIS BASED ON MARKOV PROCESSES WITH REWARDS 33

Table 5.2. Calculation results of system reliability.

R(t)
Time [day] MRM Method proposed in [23] Monte Carlo from [23] Exact analytical solution
50 0.99317 0.9917 0.9931 0.993170
100 0.93879 0.9333 0.9385 0.938794
200 0.67104 0.6610 0.6730 0.671043
300 0.38090 - - 0.380899
400 0.18958 - - 0.189584
500 0.08858 0.0864 0.0884 0.088584

The TFF is determined in accordance with (3.4) TFF=284.61 [day].
3. σ - standard deviation of random time to first failure.

The σ is obtained from (3.5), where vector of free terms is
[ −2T1 −2T2 −2T3 −2T4 −2T5 −2T6 −2T7 −2T8 −2T9 −2T10 ]T .
σ = 154.05 [day].

Next, we assume the possibility of a multi-level system operation. We assume that the states
of the system with a capacity less than 5 can be subdivided into groups (levels) with capacity
{4, 3, 2, 1, 0}. The model of the multistate system consisting of multistate components is
shown in Fig.(5.10). The properties of the levels are given in Table 5.3.

Table 5.3. The groups properties and measures for multistate model

Levels The group Capacity T i
Σ

1 S1,S2,S3,S4,S5,S6,S7,S8,S9,S10 5 273.10
2 S11,S13 4 63.48
3 S12 3 17.75
4 S15 2 16.38
5 S14 1 52.43
6 S16 0 76.86

For the second case, we calculate the following reliability and performability indicators:

1. T iΣ - average total residence time of the system in the states of the ith group during (0, t).
The T iΣ is determined in accordance with Section 3.4 . For example, T 2

Σ is calculated by
(2.15) with vector R = [ 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 ]T .
Last column of Table 5.3 shows the results of calculation of T iΣ for t = 500[day].

2. Kp(t) - performability ratio.
The Kp(t) is calculated in accordance with (3.13).
To calculate performability ratio, we should determine the reward matrix. For this
system, the reward matrix will consist only of diagonal elements. We will define reward
rates through the capacities, corresponding to the states of operation and degradation
as wii = Ci/(Cmaxt). For Cmax = 5 and t = 500 [day] we have w1 = · · · = w10 =
0.002;w11 = w13 = 0.0016;w12 = 0.0012;w14 = 0.0008;w15 = 0.0004;w16 = 0. As a
result of solving the system of equations (2.15), we obtain Kp(500) = 0.70315.

This performability index is useful in choosing the best project version of a multi-level
system, provided that the reliability requirements are met. In this case, the reliability of
the system during the time interval (0,500) [days] is extremely small (see last raw of Table
5.2). However, the value of the performability index is quite satisfactory, which confirms the
correctness of the choice of the redundancy scheme.

Copyright c© 2018 ASSA. Adv Syst Sci Appl (2018)



34 V.S. VIKTOROVA, N.V. LUBKOV, A.S. STEPANYANTS

Fig. 5.10. Markov graph of multistate DB-WSS.

Copyright c© 2018 ASSA. Adv Syst Sci Appl (2018)



A UNIFIED APPROACH TO RAP ANALYSIS BASED ON MARKOV PROCESSES WITH REWARDS 35

6. CONCLUSION

This paper presents a unified approach for the estimation of reliability, availability,
performability of engineering systems. The approach uses Markov process with rewards
to model transient and steady-state behavior of the systems. By using Markov reward
models, we can derive various instantaneous and cumulative measures and estimate wide
range of indices. We present straightforward calculation technique for evaluation of
reliability, availability and performability indices based on special definition of reward matrix.
The technique implementation requires only solution of system of differential equations
describing the behavior of the average reward. Any additional calculations, such as numerical
integration, are not required. Based on this approach a software for reliability, availability,
performability analysis was created. The software is operated under Windows OS and is
written in c] programing language. It has an advanced graphical user interface, uses effective
numerical method for solving MRM and applies OLE automation for generation reports
in MS Word formats. The efficiency of the proposed approach and accuracy of numerical
solution have been shown by the case studies of the technological system with protection and
demand-based warm standby system.

7. APPENDIX

In this section, we present analytical solution of the model of a binary-state system consisting
of multistate components (Fig.5.9). The column vector p(t) of the state probabilities is
calculated by solving a system of Kolmogorov-Chapman equations:

p′1(t) = −p1(t) · (λ1 + λw3 );

p′2(t) = p1(t) · λ1 − p2(t) · (λ1 + λ3);

p′3(t) = p1(t) · λw3 − p3(t) · (λ1 + λw3 );

p′4(t) = p2(t) · λ1 − p4(t) · (λ1 + λ3);

p′5(t) = p2(t) · λ3 + p3(t) · λ1 − p5(t) · (λ1 + λ3);

p′6(t) = p3(t) · λw3 − p6(t) · (λ1 + λw3 );

p′7(t) = p4(t) · λ1 − p7(t) · λ3;

p′8(t) = p4(t) · λ3 + p5(t) · λ1 − p8(t) · (λ1 + λ3);

p′9(t) = p5(t) · λ3 + p6(t) · λ1 − p9(t) · (λ1 + λ3);

p′10(t) = p6(t) · λw3 − p10(t) · λ1;

(7.1)

The solution is

Copyright c© 2018 ASSA. Adv Syst Sci Appl (2018)



36 V.S. VIKTOROVA, N.V. LUBKOV, A.S. STEPANYANTS

p1(t) = −e−(λ1+λw3 )t;

p2(t) =
λ1

λ3 − λw3
(e−(λ1+λw3 )t − e−(λ1+λ3)t);

p3(t) = λw3 te
−(λ1+λw3 )t;

p4(t) = (
λ1

λ3 − λw3
)
2

(e−(λ1+λw3 )t − e−(λ1+λ3)t)− λ2
1

λ3 − λw3
e−(λ1+λ3)t;

p5(t) =
λ1

λ3 − λw3
(e−(λ1+λw3 )t − e−(λ1+λ3)t)− λ1λ3t

λ3 − λw3
e−(λ1+λ3)t +

λ1λ
w
3 t

λ3 − λw3
e−(λ1+λw3 )t;

p6(t) =
(λw3 t)

2

2
e−(λ1+λw3 )t;

p7(t) =
λ1

λ1 + λw3 − λ3

e−λ3t +
λ2

1 + λ1(λ3 − λw3 )

(λ3 − λw3 )2 e−(λ1+λ3)t − λ3
1

(λ3 − λw3 )2(λ1 + λw3 − λ3)
e−(λ1+λw3 )t

+
λ2

1t

λ3 − λw3
e−(λ1+λ3)t;

p8(t) =
2λ2

1

(λw3 − λ3)2 e
−(λ1+λw3 )t +

λ2
1(λw3 − 2λ3)t

(λw3 − λ3)2 e−(λ1+λ3)t − (λ1t)
2λ3

λ3 − λw3
e−(λ1+λ3)t

+
λ2

1λ
w
3 t

(λw3 − λ3)2 e
−(λ1+λw3 )t − 2λ2

1

(λw3 − λ3)2 e
−(λ1+λ3)t;

p9(t) =
λ1λ

w
3 t

λw3 − λ3

e−(λ1+λ3)t − λ1(λ3t)
2

2(λ3 − λw3 )
e−(λ1+λ3)t − λ1

λw3 − λ3

e−(λ1+λw3 )t − λ1λ
w
3 t

λw3 − λ3

e−(λ1+λw3 )t

− λ1(λw3 t)
2

2(λw3 − λ3)
e−(λ1+λw3 )t +

λ1

λw3 − λ3

e−(λ1+λ3)t;

p10(t) = e−λ1t − (λw3 t)
2 + 2λw3 t+ 2

2
e−(λ1+λw3 )t.

(7.2)
Sought reliability is found as

R(t) =
10∑
i=1

pi(t). (7.3)

REFERENCES

1. Barlow, R. E. & Proschan, F. (1965) Mathematical theory of reliability. NY, Wiley.
2. Gnedenko, B. V., Belyaev, Y. K. & Solovyev, A. D. (1969) Mathematical methods in

reliability theory. New York, Academic Press.
3. Kozlov, B. A. & Ushakov, I. A., (1970) Reliability Handbook. New York, N Y, USA:

Holt, Rinehartand Winston.
4. Viktorova, V. S. & Stepanyants, A. S. (2016) Models and methods of reliability analysis

of technical systems. 2th ed. Moscow: LENAND Publ..
5. Smith, W. L., (1958) Renewal theory and its ramifications. J. Roy. Statist. Soc., ser. B,

vol. 20, no. 2, 243–302.

Copyright c© 2018 ASSA. Adv Syst Sci Appl (2018)



A UNIFIED APPROACH TO RAP ANALYSIS BASED ON MARKOV PROCESSES WITH REWARDS 37

6. Cox, D. R. (1962) Renewal Theory. New York, Wiley.
7. Kalashnikov, V. V. (1994) Topics on Regenerative Processes. Boca Raton: CRC Press.
8. Malinowski, J. (2013) “A Monte Carlo Method for Estimating Reliability Parameters of

a Complex Repairable Technical System with Inter-Component Dependencies,” IEEE
Trans. Reliability, vol. 62, no.1, 256–266.

9. Reibman, A. L. & Trivedi, K. S. (1988) “Numerical transient analysis of Markov
models,” Computers&Operations Research, vol. 15, no. 1, 19–36.

10. Reibman, A. L., Smith, R. & Trivedi, K. S. (1989) “Markov and Markov reward
model transient analysis: An overview of numerical approaches,” European Journal of
Operational Research, vol. 40, 257–267.

11. Lindemann, C., Malhotra, M. & Trivedi, K. S. (1995) “Numerical Methods for
Reliability Evaluation of Markov Closed Fault-Tolerant Systems,” IEEE Trans.
Reliability, vol. 44, no. 4, 694–704.

12. Ryabinin, I. A. (1976) Reliability of engineering systems. M., Mir Pbl..
13. Henley, E. J. & Kumamoto, H. (1996) Fault Tree Construction in Probabilistic Risk

Assessment and Management for Engineers and Scientists, (2nd ed.) IEEEPress, ch. 4,
sec. 4.3, 166–172.

14. Twigg, D. W., Ramesh, D. W., Sandadi, U. R. & Sharma, T. C. (2000) “Modeling
mutually exclusive events in fault trees,” Annual Reliability and Maintainability
Symposium. 2000 Proceedings. International Symposium on Product Quality and
Integrity, Los Angeles, CA, 8–13.

15. Meshkat, L., Dugan, J. B. & Andrews, J. D. (2002) “Dependability analysis of systems
with on-demand and active failure modes, using dynamic fault trees,” IEEE Trans.
Reliability, vol.51, no. 2, 240–251.

16. Tang, Z. & Dugan, J. B. (2006) “BDD-based reliability analysis of phased-mission
systems with multimode failures,” IEEE Trans. Reliability, vol. 55, no. 2, 350–360.

17. Xing, L., Morrissette, B. A. & Dugan, J. B. (2014) “Combinatorial Reliability Analysis
of Imperfect Coverage Systems Subject to Functional Dependence,” IEEE Trans.
Reliability, vol. 63, no. 1, 367–382.

18. Zhu, P., Han, J., Liu, L. & Lombardi, F. (2015) “A Stochastic approach for the analysis of
dynamic fault trees with spare gates under probabilistic common cause failures,” IEEE
Trans. Reliability, vol. 64, no. 3, 878–892.

19. Li, Z., Gu, J., X, J., Fu, L., An, J. & Dong, Q. (2017) “Reliability analysis of complex
system based on dynamic fault tree and dynamic Bayesian network,” 2017 Second
International Conference on Reliability Systems Engineering (ICRSE), Beijing, China,
1–6.

20. Smith, R. M., Trivedi, K. S. & Ramesh, A. V. (1988) “Performability analysis: measures,
an algorithm, and a case study,” IEEE Trans. Computers, vol. 37, no. 4, 406–417.

21. Amari, S. V., Xing, L., Shrestha, A., Akers, J. & Trivedi, K. S. (2010) “Performability
analysis of multistate computing systems using multivalued decision diagrams,” IEEE
Trans. Computers, vol. 59, no. 10, 1419–1433.

22. Viktorova, V. S., Sverdlik, Y. M. & Stepanyants, A. S. (2010) “Analyzing reliability for
systems with complex structure on multilevel models,” Automation and Remote Control,
vol. 71, no. 7, 1410–1414.

23. Jia, H., Ding, Y., Peng, R. & Song, Y. (2017) “Reliability Evaluation for Demand-Based
Warm Standby Systems Considering Degradation Process,” IEEE Trans. Reliability, vol.
66, no. 3, 795–805.

24. Hazra, N. K. & Nanda, A. K., (2014) “Component Redundancy Versus System
Redundancy in Different Stochastic Orderings,” IEEE Trans. Reliability, vol. 63, no.
2, 567–582.

25. Dolega, B., Kopecki, G. & Tomczyk, A. (2016) “Possibilities of using software
redundancy in low cost aeronautical control systems,” 2016 IEEE Metrology for
Aerospace (MetroAeroSpace), Padua, Italy, 33–37.

Copyright c© 2018 ASSA. Adv Syst Sci Appl (2018)



38 V.S. VIKTOROVA, N.V. LUBKOV, A.S. STEPANYANTS

26. Bolvashenkov, I., Kammermann, J. & Herzog, H.-G. (2017) “Methodology for
quantitative assessment of fault tolerance of the multi-state safety-critical systems with
functional redundancy,” 2017 International Conference on Information and Digital
Technologies (IDT), Zilina – Slovakia, 74–83.

27. Pan, X., Yang, Y., Zhang, G. & Zhang, B. (2017) “Resilience-based optimization of
recovery strategies for network systems,” 2017 Second International Conference on
Reliability Systems Engineering (ICRSE), Beijing, China, 1–6.

28. Lim, T.-J. & Lie, C. H. (2000) “Analysis of System Reliability with Dependent Repair
Modes ,” IEEE Trans. Reliability, vol. 49, no. 2, 153–162.

29. Abhilash, Chakka, R. & Challa, R.-K. (2013) “Numerical performance evaluation of
heterogeneous multi-server models with breakdowns and FCFS, LCFS-PR, LCFS-NPR
repair strategies,” 3rd IEEE International Advance Computing Conference (IACC),”
Ghaziabad, India, 566–570.

30. Amari, S. V., Myers, A. F., Rauzy, A. & Trivedi, K. S. (2008) “Imperfect Coverage
Models: Status and Trends,” in Misra K.B. (eds) Handbook of Performability
Engineering, London, Springer.

31. Spiridonov, I. B., Stepanyants, A. S. & Victorova, V. S. (2012), Design testability anal-
ysis of avionic systems. Reliability: theory&applications. [Online]. vol. 7, no. 3 (26).
66–73. Available: http://www.gnedenko-forum.org/Journal/2012/032012/RTA 3 2012-
08.pdf.

32. Ng, Y. W. & Avizienis, A. (1980) “A unified reliability model for fault tolerant
computers,” IEEE Trans. Computers, vol. C-29, no.11, 1002–1011.

33. Makam, S. V. & Avizienis, A. (1982) ”ARIES 81: A reliability and life-cycle evaluation
tool for fault-tolerant systems,” IEEE 12-th Fault-Tolerant Computing Symposium
(FTCS-12), Santa Monica, CA, 267–274.

34. Bavuso, S. J., Dugan, J. B. , Trivedi, K. S., Rothmann, E. M. & Smith, W. E.
(1987) “Analysis of Typical Fault-Tolerant Architectures using HARP,” IEEE Trans.
Reliability, vol. R-36, no. 2, 176–185.

35. Forsythe, G. E., Malcolm, M. A. & Moler, C. B. (1977) Computer methods for
mathematical computations, NJ, Prentice Hall.

36. Ikramov, K.D. (1987) Numerical methods of linear algebra. Moscow: Znanie Publ..
37. Ustinov, S. M. & Zimnitskiy, V. A. (2009) Computational mathematics. SPb., BHV-

Petersburg Publ..
38. Rakitskiy, Y. V., Ustinov, S. M. & Chernorutskiy, I. G. (1979) Numerical methods for

solving stiff systems. Moscow: Nauka Publ., 81–112.

Copyright c© 2018 ASSA. Adv Syst Sci Appl (2018)