Rta 2 2012

ELECTRONIC JOURNAL
OF INTERNATIONAL GROUP
ON RELIABILITY
JOURNAL IS REGISTERED
IN THE LIBRARY
OF THE U.S. CONGRESS
Gnedenko Forum Publications
RELIABILITY:
ISSN 1932-2321
THEORY&APPLICATIONS
VOL.1 NO.2 (25)

JUNE, 2012
Vol.1
No.2 (25)
June
2012
San Diego
ISSN 1932-2321
"Reliability: Theory & Applications", 2006, 2010, 2011
" Reliability & Risk Analysis: Theory & Applications", 2008
I.A.Ushakov, 2009
A.V.Bochkov, 2009
http://www.gnedenko-forum.org/Journal/index.htm
All rights are reserved

The reference to the magazine "Reliability: Theory & Applications"
at partial use of materials is obligatory.
Journal of International Group on Reliability
RELIABILITY:
THEORY & APPLICATIONS
Vol.1 No.2 (25),
June, 2012
Special Issue 3
100th anniversary of Boris Vladimirovich
Gnedenko's birthday
San Diego
2011
RT&A # 02 (25)
(Vol.1) 2012, June
Journal Council
Editor-in-Chief :
Send your paper
Ushakov, Igor (USA)

e-mail: igusha22@gmail.com
e-Journal Reliability: Theory & Applications

publishes papers, reviews, memoirs, and bibliographical
materials on Reliability, Quality Control, Safety,
Survivability and Maintenance.
Scientific Secretary:
Bochkov, Alexander (Russia)
e-mail: a.bochkov@gmail.com
Theoretical papers have to contain new problems,

finger practical applications and should not be
overloaded with clumsy formal solutions.
Deputy Editors:
Gertsbakh, Eliahu (Israel)
e-mail: elyager@bezeqint.net
Priority is given to descriptions of case studies.
Koowrocki, Krzysztof (Poland)

e-mail: katmatkk@am.gdynia.pl
Krishnamoorthy, Achyutha (India)
General requirements for presented papers
e-mail: krishna.ak@gmail.com
Shybinsky Igor (Russia)
1. Papers have to be presented in English in

MSWord format. (Times New Roman, 12 pt ,
1.5
intervals).
2. The total volume of the paper (with illustrations)
can be up to 15 pages.
3. presented paper has to be spell-checked.
4. For those whose language is not English, we
kindly recommend to use professional linguistic proofs
before sending a paper to the journal.
e-mail: christian.paroissin@univ-pau.fr
Singpurwalla, Nozer (USA)

e-mail: nozer@gwu.edu
Editorial Board:
Belyaev, Yuri (Sweden)
e-mail: Yuri.Belyaev@math.umu.se
Chakravarthy, Srinivas (USA)

e-mail: schakrav@kettering.edu
Dimitrov, Boyan (USA)
***
The Editor has the right to change the paper title and
make editorial corrections.
The authors keep all rights and after the publication
can use their materials (re-publish it or present at
conferences).
Publication in this e-Journal is equal to publication
in other International scientific journals.
Papers directed by Members of the Editorial Boards
are accepted without referring.
The Editor has the right to change the paper title and
make editorial corrections.
e-mail: BDIMITRO@KETTERING.EDU
Genis, Yakov (USA)

e-mail: yashag5@yahoo.com
Kaminsky, Mark (USA)

e-mail: katmatkk@am.gdynia.pl
Kovalenko, Igor (Ukraine)

e-mail: kovigo@yandex.ru
Levitin, Gregory (Israel)

e-mail: levitin@iec.co.il
Limnios, Nikolaos (France)

e-mail: Nikolaos.Limnios@utc.fr
Nikulin, Mikhail
e-mail: M.S.Nikouline@sm.u-bordeaux2.fr
Nelson, Wayne (USA)
The authors keep all rights and after the publication

can use their materials (re-publish it or present at
conferences).
Send your papers to
e-mail: WNconsult@aol.com
Popentiu, Florin (UK)

e-mail: Fl.Popentiu@city.ac.uk
Rykov, Vladimir (Russia)

e-mail: rykov@rykov1.ins.ru
the Editor-in-Chief ,
Igor Ushakov
igusha22@gmail.com
Wilson, Alyson (USA)

e-mail: agw@lanl.gov
Wilson, Simon (Ireland)

e-mail: swilson@tcd.ie
Yastrebenetsky, Mikhail (Ukraine)
or
e-mail: ma_yastreb@mail.ru
Zio, Enrico (Italy)
the Deputy Editor,

Alexander Bochkov
a.bochkov@gmail.com
e-mail: zio@ipmce7.cesnef.polimi.it
Technical assistant
Ushakov, Kristina
e-mail: kudesigns@yahoo.com
RT&A # 02 (25)
(Vol.1) 2012, June
Table of Contents
Table of Contents
Igor Ushakov
SENSITIVITY ANALYSIS OF OPTIMAL REDUNDANCY SOLUTIONS .............................................................. 7
Solving practical optimal redundancy problems, one can ponder: what is the sense of optimizing if input
data are taken from the ceiling? Indeed, statistical data are so unreliable (especially in reliability
problems that such doubts have a very good ground. Not found any sources after searching the answer for
this question, the author decided to make some investigation of optimal solutions sensitivityunder influence
of data scattering. Thus, one can notice that input parameters variation may influence enough significantly
enough on probability of failure-free operation and the total system cost from run to run of Monte Carlo
simulation though optimal solution remains more or less stable.
Shubinsky Igor B., Zamyshlyaev Alexey M.
TOPOLOGICAL SEMI-MARKOV METHOD FOR CALCULATION OF STATIONARY
PARAMETERS OF RELIABILITY AND FUNCTIONAL SAFETY OF TECHNICAL SYSTEMS ............................ 12
The paper offers the method for calculation of reliability parameters and functional safety of technical
systems, differing from known methods by an opportunity of obtaining strict formula expressions of
stationary parameters directly from a system state graph. The method is suitable for solution of both
Markov, and semi-Markov models of reliability and safety. In addition the paper presents some examples
of determining safety and availability factors, as well as time parameters of safety and reliability of the twochannel safety related device.
G. Tsitsiashvili, M. Osipova
RECURRENT SEQUENCE OF PARALLEL-SERIAL CONNECTIONS.................................................................... 23
In this paper a sequence of parallel-serial connections is considered. In this sequence next connection is
obtained by parallel or serial linking of new arc to obtained connection. Distributions of random numbers
of connectivity components are analyzed. These distributions are considered intensively now. Central limit
theorem is proved for these distributions and parameters (mean and variance) of normal limit distribution
are calculated.
Igor Ushakov
UNIVERSAL GENERATING FUNCTION & OPTIMAL REDUNDANCY PROBLEMS .......................................... 28
The Method of U-functions, or the Method of the Universal Generating Function (UGF), was introduced in
[(1986) Ushakov, (1987) Ushakov] and later developed in [(1988) Ushakov; (1995) Gnedenkon&
Ushakov]. Actually this is a generalization and algebraic formalization of the well-known Kettelles
Algorithm [(1962) Kettelle]. In turn, Kettelles Algorithm, is a form of presentation of convolution of
discrete random variables. The method of U-functions is very convenient for computerized calculations.
Last years, this method was significantly developed by G. Levitin and A. Lisnianski
Alexander V. Bochkov, Valery V. Lesnykh, Nikolay N. Zhigirev
DYNAMIC MULTI-CRITERIA DECISION MAKING METHOD FOR SUSTAINABILITY RISK
ANALYSIS OF STRUCTURALLY COMPLEX TECHNO-ECONOMIC SYSTEMS ................................................. 36
The paper considers the characteristics of the functioning and sustainability of the structurally complex
techno-economic systems (SCTES) in terms of different types of risk. The validity of the application to
describe the behavior of this class of systems of semi-empirical mathematical models, which are based on a
vector description of the system states, using the criteria approach for assessing the quality of its
functioning, is demonstrated. Under discussion is conceptual model for the interaction of the object and its
environment, allowing to estimate the "optimal" allocation ratio between the productive system and its
development potential. The concept of non-formalizable threats for the sustainable functioning of this
systems class was introduced. Expert procedure to account non-formalizable threats in case of risk
assessment was proposed. For the construction of indicators for assessing the status the methods of
quantitative analysis based on the theory and multi-criteria utility was used. Multi-criteria utility as an
indicator of sustainability in the form of dimensionless complex hierarchy of indicators was proposed.
Computing for through the convolution of the primary indicators.
RT&A # 02 (25)
(Vol.1) 2012, June
Table of Contents
Fedol Ghazi, Pr: Hamouine Abdelmadjid, Dr: Zouaoui Chikr el Mezouar

BEAM BEHAVIOUR UNDER MONOTONIC LOADS ............................................................................................. 43
In this work we interest to study the beam behavior under monotonic loads in four point, to improve the
mechanical properties of a concrete beam fiber and establish an identification card of the new concrete
beams were comparing these beams witnesses.
G. Tsitsiashvili, M. Osipova
BASIS OF CONTINUUM APPROXIMATION IN MODELS OF GROWING RANDOM
NETWORKS .............................................................................................................................................................. 52
Accuracy asymptotics for differences between prelimit and limit distributions of nodes powers in models of
random growing networks are constructed. A rate of a convergence in these relations is power. Obtained
formulas allow to ground the continuum approximations for considered models.
A. Blokus-Roszkowska, L. Smolarek
COLLISION RISK ESTIMATION FOR MOTORWAYS OF THE SEA ..................................................................... 58
The Motorways of the Sea is rather a new concept and, thus, is still in the process of development by the
European Commission and in the Baltic Sea Region. The Baltic Sea Region is one of the most dynamic
growth areas. Due to this fact the new ideas and technologies are needed to optimize the sea transport
system. In the paper the simulation model for the system safety state evaluation is presented. The
simulation program can constitute a base for decision-support tool, on the level of safety management,
especially to optimally plan the safety transport system.
Kumar Pardeep, Chaturvedi D.K., Pahuja G.L.
A HEURISTIC METHOD FOR RELIABILITY REDUNDANCY OPTIMIZATION OF FLOW
NETWORKS .............................................................................................................................................................. 69
In flow networks, from the quality and service management point of view, measurement of the transmission
ability of a network to meet the customers demand is very important. To meet the ultra high reliable
requirements of such networks, a heuristic method for reliability redundancy optimization of flow networks
using composite performance measure (CPM) integrating reliability and capacity has been proposed. The
method is based upon the selection of main flow paths and backup paths and then optimizing main paths on
priority basis. Thus, the reduced computation work makes the proposed algorithm suitable for designing of
large, reliable telecommunications networks.
A.N. Patowary, J. Hazarika and G. L. Sriwastav
ESTIMATION OF RELIABILITY IN INTERFERENCE MODELS USING MONTE CARLO
SIMULATION............................................................................................................................................................ 78
This paper presents estimation of reliability R P(X Y) of a system, for the cases when its strength (X) and
stress (Y) follow exponential, normal or gamma distributions, using Monte Carlo simulation (MCS). First
the parameters of strength and / stress are estimated and substituting them in the reliability expressions, in
different cases, the estimates of reliability are obtained. Normal distribution is fitted to various sets of
estimated reliability R , generated by MCS. The goodness of fit is tested using Kolmogorov-Smirnov one
sample test.
V. M. Chacko, M. Manoharan
MULTISTATE COHERENT SYSTEMS WITH MULTIPLE STATE TRANSITIONS AT A TIME ........................... 85
A key requirement in defining a multistate coherent system (MCS) is the relevance condition of its
components. A new class of MCSs is introduced with a new component relevance condition. Also we
introduce a more general relevance condition. They are compared with some existing component relevance
conditions. Based on the two new relevance conditions, two component importance measures for MCSs are
defined. They are most appropriate for comparing components when certain type of system improvement is
sought. We introduce new joint importance measures for two or more components with respect to the
proposed relevance conditions. The new MCS classes include several existing MCSs as special case. An
illustrative example of the proposed MCSs is also provided.
RT&A # 02 (25)
(Vol.1) 2012, June
I.Ushakov - SENSITIVITY ANALYSIS OF OPTIMAL REDUNDANCY SOLUTIONS
SENSITIVITY ANALYSIS OF OPTIMAL REDUNDANCY SOLUTIONS

Igor Ushakov
San Diego, California, USA

Solving practical optimal redundancy problems, one can ponder: what is the sense of
optimizing if input data are taken from the ceiling? Indeed, statistical data are so unreliable
(especially in reliability problems) that such doubts have a very good ground.
Not found any sources after searching the answer for this question, the author decided to
make some investigation of optimal solutions sensitivityunder influence of data scattering.
A simple series system of six units has been considered (see Figure 1). For reliability
increasing, one uses a loaded redundancy, i.e. if a main unit k has xk redundant units, its reliability
is found from
Pk ( x k ) 1 (1 p k ) xk 1
where pk is a probability of failure free operation (PFFO) of a single unit k. And the total cost of xk
redundant units is equal to ckxk, where ck is the cost of a single unit of type k.
Figure 1. Series system underwent analysis

Units parameters are presented in Table 1.
Table 1. Input data
Unit 1 Unit 2 Unit 3 Unit 4 Unit 5 Unit 6
0.8
0.8
0.9
0.9
0.9
pk 0.8
5
5
1
1
1
ck 5
Assumed that units are mutually independent, i.e. systems reliability is defined as
P (x
PSystem ( x k , 1 k 6 )
1 k 6
And the total systems redundant units cost is:

C System ( x k , 1 k 6 )
xk
1 k 6
Below presented solutions of both problems of optimal redundancy: direct:

min C ( xk , 1 k 6 ) P ( xk , 1 k 6 ) P *
1 xk
RT&A # 02 (25)
(Vol.1) 2012, June
Ushakov I. - SENSITIVITY ANALYSIS OF OPTIMAL REDUNDANCY SOLUTIONS
and inverse:
max P( x
C ( x k , 1 k 6 ) C *
, 1 k 6 )
1 xk
For finding the optimal solutions, the Steepest Descent Method was applied. For this base
system the solutions for several sets of parameters are presented for Direct Problem in Table 2 and
for Inverse Problem in Table 3. (Numbers are given with high accuracy only for demonstration
purposes; in practice, one has to use only significant positions after a row of nines.)
Table 2. Solution for Direct problem
P*
0.95
0.99
0.995
0.999
x1
3
4
5
6
x2
3
4
4
5
x3
3
3
4
5
x4
3
3
4
4
x5
2
3
3
4
x6
2
3
3
4
Achieved P
0.9559520
0.991187
0.995229
0.999218
System C
52
69
75
93
Table 3. Solution for Inverse problem

C*
50
75
100
x1
3
4
5
x2
3
4
4
x3
2
3
4
x4
2
3
4
x5
2
3
3
x6
2
3
3
Achieved C
46
75
99.5
System P
0.931676
0.995229
0.999602
The questions of interest are: how optimal solution will change if input data are changed?
Two types of experiments have been performed: in the first series of experiments, different units
costs with fixed probabilities were considered (see Figure 2) and in another one different units
probabilities with fixed costs were considered (see Figure 3).
Figure 2. Input data for the first series of experiments
Figure 3. Input data for the second series of experiments

The results of calculations are as follows:
RT&A # 02 (25)
(Vol.1) 2012, June
Table 4. Values of Probabilities of Failure-free operations
Initial
Various C
Various P
0.999
0.999218
0.998996
0.999218
0.995
0.99566
0.99566
0.99566
0.99
0.9922
0.9922
0.9922
0.95
0.955952
0.955952
0.955952
In conclusion , there was performed a Monte Carlo simulation where parameter of the PFFO
and cost were changed simultaneously. In this case, parameters of probabilities of each unit were
calculated (in Excel) as:
pk=0.8pk+0.4pk*RAND()
and
ck= 0.8ck+0.4*RAND(),
i.e. considered a random variation of the values within 20% limits..
The final results for this case are presented in Tables 5 8.
Table 5. Results of Monte Carlo simulations for P*=0.999
No.
1
2
3
4
5
6
7
8
9
10
P* = 0.999
P
C
0.999352 100
0.999218 102
0.999313 102
0.999212 97
0.999182 102
0.999171 97
0.999171 103
0.999596 100
0.999526 100
0.999399 100
x1
6
6
6
5
6
6
6
6
6
6
x2
6
6
6
6
6
6
6
6
6
6
x3
6
6
6
6
6
5
6
6
6
6
x4
4
5
4
4
4
4
4
4
4
4
x5
4
4
4
4
4
4
5
4
4
4
x6
4
4
4
4
4
4
4
4
4
4

P* = 0.995
No. P
C
1
0.995478 84
2
0.996755 85
3
0.995026 85
4
0.996777 79
5
0.996777 84
6
0.995525 79
7
0.996732 85
8
0.996732 85
9
0.995645 84
10 0.99567 84
x1
5
5
5
4
5
5
5
5
5
5
x2
5
4
4
5
5
5
5
5
5
5
x3
5
4
5
5
5
4
5
5
5
5
x4
3
4
4
3
3
3
3
3
3
3
x5
3
3
3
3
3
3
4
4
3
3
x6
3
3
3
3
3
3
3
3
3
3
RT&A # 02 (25)
(Vol.1) 2012, June
No.
1
2
3
4
5
6
7
8
9
10
P* = 0.99
P
0.990147
0.990965
0.990229
0.99185
0.990389
0.99107
0.992185
0.990422
0.990893
0.990466
C
69
70
70
69
71
69
74
71
71
69
x1
4
4
4
4
4
4
5
4
5
4
x2
4
4
4
4
4
4
4
4
4
4
x3
4
4
4
4
4
4
4
4
4
4
x4
3
4
4
3
4
3
3
3
3
3
x5
3
3
3
3
4
3
3
4
3
3
x6
3
3
3
3
3
3
3
3
3
3
No.
1
2
3
4
5
6
7
8
9
10
P* = 0.95
P
0.950045
0.955842
0.951936
0.951711
0.957883
0.951908
0.962227
0.962227
0.95261
0.950393
C
52
52
52
54
50
51
51
51
50
52
x1
3
3
3
3
3
3
3
3
3
3
x2
3
3
3
3
3
3
3
3
3
3
x3
3
3
3
3
3
3
3
3
3
3
x4
3
3
3
2
3
2
2
3
3
3
x5
2
2
2
2
3
2
2
3
2
2
x6
2
2
2
2
2
2
2
2
3
2
Analysis of data presented in Tables 5-8 shows relatively significant difference in numerical
results (see Figure 4).
Figure 4. Deviation of maximum and minimum values of probability

of failure-free operation in results of Monte Carlo simulation
10
RT&A # 02 (25)
(Vol.1) 2012, June
However, the problem is not in coincidence of final values of PFFO or cost. The problem is:
how change of parameters influences on the optimal values of x1, x2, .
However, one can observe that even with a system of six units (redundant groups) a visual
analysis of sets (x1, x2, , x6) is extremely difficult and, at the same time, deductions based on
some averages or deviations of various xk are almost useless.
The author was forced to invent some kind of a special presentation of sets of xks. Since
there is no official name for such kind of graphical presentation, it is called multiple polygons (in
Russian -). On such multiple polygon there are numbers of rays
corresponding to the number of redundant of units (groups). Each ray has several levels
corecponded to the number of calculated redundant units for considered case (see Figure 5).
Figure 5. Multiple polygon axes with numbered levels

These multiple polygons give a perfect visualization of close-to-optimal solutions and
characterize observed deviation of particular solutions. Such multiple polygons for considered
example are given in Figure 6. (Here bold lines re used for connection xk obtained as optimal
solution for units with parameters given in Table 1.)
P* = 0.999
P* = 0.999
P* = 0.999
P* = 0.999
Figure 6. Deviations of optimal solutions for randomly variation of parameters

from the optimal solution obtained for parameters given in Table 1
Thus, one can notice that input parameters variation may influence enough significantly
enough on probability of failure-free operation and the total system cost from run to run of Monte
Carlo simulation though optimal solution remains more or less stable.
11
I. Shubinsky, A. Zamyshlyaev - TOPOLOGICAL SEMI-MARKOV METHOD FOR CALCULATION OF STATIONARY PARAMETERS OF RELIABILITY AND
FUNCTIONAL SAFETY OF TECHNICAL SYSTEMS
RT&A # 02 (25)
(Vol.1) 2012, June
TOPOLOGICAL SEMI-MARKOV METHOD FOR CALCULATION OF

STATIONARY PARAMETERS OF RELIABILITY AND FUNCTIONAL
SAFETY OF TECHNICAL SYSTEMS
Shubinsky Igor B.
Doctor of Technical Sciences, professor, director

of the company " IBtrans", The Public Corporation RZD, Moscow,
Russian Federation
e-mail: igor-shubinsky@yandex.ru
Zamyshlyaev Alexey M.
Ph.D., The Head of Scientific-Technical Complex of JSC "NIIAS",

The Public Corporation RZD, Moscow, Russian Federation
e-mail: A.Zamyshlyaev@gismps.ru
ABSTRACT
The paper offers the method for calculation of reliability parameters and functional safety of
technical systems, differing from known methods by an opportunity of obtaining strict formula
expressions of stationary parameters directly from a system state graph. The method is suitable for
solution of both Markov, and semi-Markov models of reliability and safety. In addition the paper
presents some examples of determining safety and availability factors, as well as time parameters of
safety and reliability of the two-channel safety related device.
Keywords: reliability, functional safety, parameters of reliability and safety, Markov and
semi-Markov models of reliability, loop weight, graph breakdown weight.
1. Introduction
When solving problems of reliability and functional safety of technical systems
mathematical tools of Markov and semi-Markov random processes are widely applied.
Development and solution of Markov and semi-Markov reliability models by traditional methods in
general terms is brought into making up a system of the homogeneous differential equations
describing behavior of the investigated system, their operational transformation, solving system of
equations in the operational form, inverse transformation and finding the required reliability
parameters. Such a way is always fraught with mathematical difficulties, especially when the
number of equations exceeds ten and is problematic to execute correctly inverse transformations of
solutions of system equations obtained in the operational form. Therefore in the majority research
people and, especially practical workers, are compelled to introduce a lot of assumptions which
radically simplify solution of reliability models and allow obtaining reliability parameters of
considered systems in the analytical or numerical form. However these results are already far from
true and there is a natural question: whether it is necessary to aspire to realization of the traditional
plan of construction and solving reliability models of systems.
In many problems of reliability calculation it is enough to be limited by stationary reliability
and availability factors (parameters). In these cases it is necessary to switch over from the model of
differential equations to the model of algebraic equations describing system behavior in the steadystate mode, to solve them, find stationary probabilities of staying system in each of possible states.
12
RT&A # 02 (25)
(Vol.1) 2012, June
Then based on system failure criteria with the help of the specified probabilities probabilistic
system availability and unavailability should be found. Thereafter stationary parameters of system
non-failure operation and maintainability should be defined. Such problems are not connected with
necessity of operational calculus application for development and solution of reliability models. The
required stationary reliability parameters are calculated sufficiently strictly. However, alongside
with the fact that the given plan does not provide definition of a full list of reliability parameters,
there is also the unsolved problem of the big dimension of algebraic equations model. Therefore
even at rather small number of states it is not possible in many cases to analytically describe
required reliability parameters of the system. This circumstance does not depend on a degree of
system graph model connectivity. Dimensions of a system matrixes for algebraic equations
representing reliability of investigated technical system model do not vary both with weak
connectivity and with strong connectivity.
At the same time, graph models of complex systems reliability, as a rule, are poorly
connected. This circumstance has stimulated us to switch from the traditional plan of solving linear
algebraic equations by Kramer's rule to the scheme of breakdown initial graph to the constituent sub
graphs which are not containing single out nodes (for example, disabled states of model or states
which are being on the way from one node to another, or an initial system state). At application of
such a scheme (plan) it has turned out sufficient for solving the system of algebraic equations.
Moreover it has turned out sufficient to be limited by finding ways and loops on the graph, what is
now well formalized.
2. Problem definition
Stationary parameters of technical systems reliability are factors of availability K and
unavailability K , mean time to failure T , dispersion of mean time to failure D , an average
time between failures T , mean idle time average of a T . As functional dependence
K f (T, T ) is known methods of calculation of these three parameters are expedient for
considering simultaneously. Similarly it is necessary to simultaneously consider methods for
calculation of parameters T and D .
Stationary parameters of functional safety of safety related systems is a factor of safety K
, mean time to dangerous failure T , dispersion of time to dangerous failure D , mean time to
protective failure T , dispersion of time to protective failure D , mean time to dangerous
(hazardous) failure T .
Methods of calculation of corresponding groups of reliability and functional safety
parameters are identical. Difference is only in the separation of initial system states on efficient and
disabled (concerning reliability) subsets and nonhazardous, hazardous and protective subsets of
states (relating to functional safety). So in the study [1] the following formula of calculation of
system availability factor which behavior is described by semi-Markov random process is
determined
PiTi
iS
(1)
K
(i, j S ; S S) ,
Pj T j
jS
where S is a subset of efficient system states , S is full set of system states; Pi is a final
probability of staying Markov chain in i-th state; T i is an expectancy of unconditional time of
system staying in i-th state.
By turn, formulas of the calculation of mean time to failure (time between failures) and the
average idle time of the system which behavior is described by semi-Markov random process, and
13
RT&A # 02 (25)
(Vol.1) 2012, June
determined according to the study [2] are the following:

P iTi
iS
T0
;
Pi p ij
iS
jS
P T
j
(2)
jS
P p
j
jS
(3)
ji
iS
where it is implied, that transition from a subset S into a subset S can be carried out not from
any working state capacity, but only from boundary conditions (subset) S . Similarly transition
from S into S can be carried out from a subset belonging to the subset of boundary disabled
states S .
Practical methods for calculation of functional safety parameters of recoverable safety
related systems nowadays are poorly developed.
The purpose of this paper consists in development of practical methods of calculation listed
above stationary parameters of reliability and functional safety of complex technical systems. It is a
question of formalization of calculations on the basis of the graph theory.
3. Calculation of availability and safety factors
3.1. Topological concepts:
Path chain of consistently connected unidirectional arcs starting from a state i and
ending in a state j, path weight l ij pir p rj , where pir - probability of one step transition for i
i ,r , jS
th state in a state r;
Closed loop is a chain of consistently connected unidirectional arcs in which the output of
final vertex in the chain is connected to starting vertex of the chain;
Weight of j th loop C j pij p ji ; self-loop is a special case of the closed loop
i , jS
(entering and leaving arcs in self-loop merge into one arch), weight of a self-loop C j p jj ;
Graph decomposition - a graph part which is not containing assigned vertices and arcs
connected with them; weight of a decomposition G i is calculated taking into account the
exclusion from the graph a vertex i and the arcs connected with it; the weight of a decomposition
G Si is calculated taking into account the additional exclusion from the graph the vertices of a set
of disabled states S and the arcs connected with them; the weight of a decomposition G ki is
calculated taking into account the exclusion from the graph the vertex i, as well as the vertices
located on k-th path from starting vertex into a vertex i and arcs connected with them;
The decomposition weight (determinant) is calculated under Masons formula
G 1 C j C r C j Ci C r C j ...
j
rj
(4)
irj
Application of Masons formula allows to considerably reducing labor input of calculations

of minors on the rarefied matrixes, and matrix G, as a rule, is rarefied.
14
RT&A # 02 (25)
(Vol.1) 2012, June
3.2. Topological formulas of calculation of availability and functional safety factors

The statement 1. If system reliability is modeled by means of the graph states and semiMarkov random process on this set of states, specified by transition probability matrix and a vector
of unconditional expectances of staying time in each graph state the factor of system availability in
the topological form is equal:
G i Ti
iS
G T
(5)
jS
where G i - decomposition weight of the graph without a state i, Ti - expectance of unconditional

staying time of the system in states s i .
Proof. The stationary probability of enclosed homogeneous Markov chain staying in a state
i , is equal
D
Pi n i ,
Dj
j 1
where n - number of states in initial set of states of the system S

2 n
1- p11 - p12 -
p1n
- p 21 1- p 22 - p 2 n
.
pn1 - pn 2 1- pnn
and Di ( D j ) - a minor obtained by deletion of i (j) line and of i (j) column in matrix D . In turn,
both the determinant D and minors Di ( D j ) can be strictly or with acceptable accuracy calculated
under Masons formula (4). Hence, the stationary probability of enclosed homogeneous Markov
chain staying in a state i , is equal to the following
Pi
Di
D
j 1
G i
G i
G j
(6)
jS
j 1
Substituting the formula (6) in the expression (1) we obtain the formula (5), as was to be
shown.
The consequence 1. If S S - a subset of nonhazardous states of safety related system.
The system safety factor is determined as the following
15
RT&A # 02 (25)
(Vol.1) 2012, June
G T
iS1
(7)
G T
jS
The formula (7) is obtained by analogy to the formula (5) concerning set of nonhazardous
states.
4. Topological formulas for calculation of stationary time parameters of reliability

Mean time to system failure
i
G T
T0
iS
(8)
G p
iS
ij
jS
Average idle time of a system

j
G T
T
jS
(9)
ji
jS iS
Formulas (8) and (9) are obtained from formulas (2) and (3) by substitution in them the
formula (6).
Mean time between hazardous failures
G i Ti
iS1
T
,
(10)
G i p ij
iS1
jS1
where S 1 - a subset of nonhazardous states, S 1 - a subset of hazardous states S 1 S 1 S , S 1 - a

subset of boundary nonhazardous states ( S1 S1 ).
Mean time to system failure and dispersion of mean time to failure
With a view of development of the formalized engineering methods for determining these
parameters we shall prove the following statement.
The statement 2. If system reliability is modeled by means of the graph states and semiMarkov random process on this set of states then the confidence curve to system failure in Laplace
transformations at i-th initial state is determined by the following expression
~ ij
~
lk ( z ) G kj ( z )
~
jS k
i ( z)
,
(11)
~
G S ( z )
~
where lk ij ( z ) - k-th path in Laplace transformations, leading from an efficient state of the graph
~
i S into the failure state j S ; G kj ( z ) - graph decomposition weight in Laplace
~
transformations without j-th vertex and the graph vertices located on the k-th path; G S ( z) without vertices, graph decomposition weight of a set of failure states
16
RT&A # 02 (25)
(Vol.1) 2012, June
Proof
In the study [1] it is shown, that the function of time distribution of system staying in a set
of efficient conditions S in Laplace transformations can be obtained from the following equation
~
~
~
~
i ( S ) Qij ( S ) j ( S ) Qil ( S ) .
jS
lS
Let's transform this equation to a matrix form, keeping in mind, that the right part of the
equation is a vector-column of free terms of semi-Markov transitions probabilities for one step
from vertices i, j ,...z S into the vertex l S .
~
~
~
~
(S ) Q ( S ) (S ) Q * ( S ) ,
~
~
~
~
where Q ( S ) (Qij ( S )) - is a matrix of semi-Markov probabilities; Q * (S ) (Qil ( S )) - is a vectorcolumn.
~
In system of the equations the unknown elements are those of the vector-column ( S ) .
After their grouping in the left part we shall obtain
~
~
~
( S ) I Q ( S ) Q * ( S ).
(S )
~
~
Then by Kramer's rule we can find i ( S ) i
, where ( S ) I Q ( S ) , and i (S ) ( S )
~
the determinant, obtained by replacement of i-th column in the matrix I Q ( S ) on a vector~
column of free terms Q * ( S ) provided that i (S ) and (S ) are not equal to zero.
The determinant i (S ) differs from the determinant (S ) G S by the fact that in the
column i the element ~p ( S ) is replaced with the element ~
p ( S ) where i, j S , and l S . As a
il
ij
result we obtain the following
i ( S ) GSi
~ il
( S ) G kl ( S ) .
lS k
Hence,
~ il
~
i ( S )
l S k
( S ) Gkl ( S )
~
GS
and at replacement of an index l on j the required result is obtained. The statement is proved.
From the formula (11) follows
- Mean time to system failure at initial state i=1
~
1 ( z )
z 0
- Dispersion of the mean time to system failure at initial state i=1

2
~
~

2 1 ( z)
1 ( z)
D
z 0 .
z 0
z 2
z
17
(12)
(13)
RT&A # 02 (25)
(Vol.1) 2012, June
Consequence 2 The function of time distribution (confidence curve) to system hazardous

failure in Laplace transformations at i-th initial state is determined by the following expression
~
where lk ij ( z ) - k-th path in Laplace transformations leading from a nonhazardous state of
the graph i S into hazardous failure state j S ;
~ ij
~
lk ( z ) G kj ( z )
~
jS k
i ( z )
,
~
G S ( z )
(14)
From the formula (14) follows

- Mean time to system failure at an initial state i=1
~
1 ( z )
T
z
z 0
(15)
- A dispersion mean time to system failure at an initial state i=1
~
2 1 ( z )
z 2
z 0
2
~
1 ( z )
z 0
(16)
5. Examples
Example 1
The two-channel device is analyzed. It contains two identical and independent channels, as
well as diagnostics tools which check with acceptable frequency for good safety the functioning
state of each channel and compare their output results. Failures of channels are asymmetrical. When
diagnostics tools are sound the fact of failure of any one channel is detected and then the device
transition in a state of protective failure is carried out. In the case of diagnostics tools failure only a
nonhazardous failure of the device can occur. The subsequent behind this event failure of a channel
leads to hazardous failure of the device.
Graph states of reliability and safety of the two-channel device with diagnostics tools
without channels restart is shown on fig. 1.
The description of states:
1 - Serviceable state;
2 - Diagnostics tool failure;
3 - protective failure of the device caused by detected failure of one of the channels;
detection was carried by regular diagnostics tools with probability ;
4 - Not detected failure of one of the channels, owing to failure or insufficient efficiency of
diagnostics tools (hazardous failure of the device).
For presentation of an illustration of opportunities of the offered method we assume, that
failure and recovery flows, as well as a flow of detected failures of one channel are the simple
flows with rates , , .. Restoration is carried out in the state of protective failure 2.
Graph edges on fig. 1 are marked by following parameters: - failure rate of diagnostics
tools; 2 - failure rate of the two working channels; - recovery rate of failures by one repair
team.
Transition from a hazardous state 3 into initial state 0 is shown. The edge 3-0 is marked by
parameter c - recovery rate of hazardous failure of the device, where the factor 0 c 1 . If for the
18
RT&A # 02 (25)
(Vol.1) 2012, June
elimination of hazardous failure there is no need to update the device then c=1 and the rate of
hazardous failure elimination is equal to recovery rate of the device. If it is required to update the
device depending on duration of updating time the given factor will have the value c 1 which
is much less 1. The opportunity of failure of one more channel when the device is in the condition
3, is not considered, as hazardous failure has already taken place and either one channel or two
channels are subjects to recovery.
2(1-v)
2
2v
3
Figure 1. Graph of safety states of the two-channel device without channels restart
The model of reliability and functional safety of the two-channel device on fig.1 provides
the following logic of device operation: an initial state 0 (all elements of the device are serviceable).
In case of diagnostics tools failure there is a transition into a state 1. If at serviceable diagnostics
tools any one channel (a state 2) has failed, and the channel failure is detected in due time out with
probability the device is transferred into a state of protective failure (the device does not function,
the channel is under repair). At the latent failure of the channel probability 1 or at failure of one
channel after the failure of diagnostics tools (the path 0 - 1 - 3) there is a transition into transition a
state 3 of hazardous failure.
Failure criterion: S {0,1} S {2,3} S S S .
Hazardous failure criterion: S {0,1,2} S {3} S S S .
It is required by means of formulas (5) and (7) to determine availability and functional
safety factors of the two-channel device
Solution
- Initial parameters should be defined:
T0 exp[ (2 l ) t ]dt
0
p 01 l exp[ (2 l ) t ]dt
0
1
1
1
1
; T1
; T2 ; T3
;
2 l
2
l
2(1 )
2
; p 02
; p 03
;
2 l
2 l
2 l
p13 p 20 p 30 1
G 0 1; G 1 1 p 02 p 20 p 03 p 30 1 p 02 p 03 ; G 2 1 p 01 p 03 ;
G 3 1 p 02
- Find availability factor
19
G i Ti
iS
Tj
G 0T0 G 1T1
G 0T0 G 1T1 G 2 T2 G 3T3
RT&A # 02 (25)
(Vol.1) 2012, June
jS
2 2 ( 2 )
2(1 )
1
2
2 2 ( 2 ) ( 2 )
c (2 )
(2 )
,
(2 ) 4 2 (1 (1 )) 2
If for elimination of hazardous failure there is no need to update the device then the factor
=1 and expression for availability factor of the device will be transformed in to the following
form:
.
2
- Find safety factor of the device
G i Ti
iS
Tj
G 0T0 G 1T1 G 2T2
G 0T0 G 1T1 G 2T2 G 3T3
jS
1
2
2 2 ( 2 ) ( 2 )
2(1 )
1
2
2 2 ( 2 ) ( 2 )
c (2 )
At =1 the safety factor of the device is determined by means of the following expression:
4 2
.
K
2 ( 2 ) ( 2 )
If the two-channel system is inefficient (in extreme case 0 ) the safety factor of the
device is equal, as one would expect, to its availability factor.
Example 2
In conditions of the example 1 it is required by means of formulas (11), (12), (14), (15) to
determine time parameters of reliability and functional safety of the two-channel device
Solution
From the formula (11) follows, that functions of time distribution to system failure at an
initial state 0 in operational transformations has the following form:
~
p ( z) ~
p13 ( z ) ~
p 03 ( z ) ~
p 02 ( z )
~
0 ( z ) 01
,
1
~
~
~
~
~
where l1 03 ( z ) ~p 01 ( z ) ~p13 ( z ); l2 03 ~
p 03 ( z ); G13 ( z ) 1; G 23 1; G S ( z ) 1 .
20
RT&A # 02 (25)
(Vol.1) 2012, June
In Laplace transformations at exponential distributions of random variables
zt
( 2 ) t
zt
e
dF
(
t
)
e
d
[
1
e
]
0
0
2 0
2 z
~
p01 ( z )
~
p13 ( z ) e zt dF1 (t ) e zt d [1 e 2 t ]
0
2
2 z
2(1 ) zt
2(1 )
~
p 03 ( z )
e dF0 (t ) e zt d [1 e ( 2 )t ]
2 0
2 z
0
2
2
e zt dF0 (t ) e zt d [1 e ( 2 )t ]
2 0
2 z
0
~
p 02 ( z )
~
p20 ( z ) e zt dF2 (t ) e zt d [1 e t ]
Hence,
~
0 ( z)
~
0 ( z )
2
T
2 z
z
z 0
1
.
2
The formula of function of time distribution to system failure in Laplace transformations

under conditions of the given example has the following form:
~
p ( z) ~
p ( z) ~
p 03 ( z )
~
0 ( z ) 01 ~ 13
,
~
1 p 02 ( z ) p 20 ( z )
~
as conditions 0,1 and 2 are non hazardous and G S 1 ~
p 02 ( z ) ~
p 20 ( z )
Hence,
~
0 ( z )
2 [ (1 )(2 z )]( z )
[2 (2 ) (2 ) z z 2 ] (2 z )
From here
~
2 ( 2 )
0 ( z )
.
z 0
z
2 [2 (1 ) ]
If to take into account, that ; , with an margin error less than the first
infinitesimal order then the following expression is true
1
.
T
2 (1 )
At high efficiency of detection of hazardous failures on the basis of two-channel architecture
of the device ( 1 ) its safety depends only on the reliability of the built in diagnostics tools and
the comparator (i.e. on failure rate ).
The conclusion
The offered topological semi-Markov method for calculation of reliability and safety
parameters of technical systems allows determining directly on the states graph the strict or
approximates formula expressions of typical reliability and safety parameters of technical systems
which behavior is described by both Markov, and semi-Markov random processes. Mathematical
21
RT&A # 02 (25)
(Vol.1) 2012, June
positions of the method are illustrated by examples which show simplicity and rigor of finding out
the required reliability and safety parameters.
The literature
1. Korolyuk V.S., Turbins A.F. Semi-Markov processes and their applications. - Kiev:
Naukova dumka, 1976, 179 p..
2. Rainshke K., Ushakov I.A. Estimation of system reliability with use of graphs /under
under the editorship of. I.A.Ushakov.-.: Radio and communication, p.p.1988.-209.
22
RT&A # 02 (25)
(Vol.1) 2012, June
G. Tsitsiashvili, M. Osipova, N. Markova REQURRENT SEQUENCE OF PARALLEL-PARALLEL-SERIAL CONNECTIONS
RECURRENT SEQUENCE OF PARALLEL-SERIAL CONNECTIONS

G. Tsitsiashvili*, M. Osipova**
* IAM FEB RAS, Vladivostok, Russia, ** FEFU, Vladivostok, Russia

* guram@iam.dvo.ru, ** mao1975@list.ru
ABSTRACT
In this paper a sequence of parallel-serial connections is considered. In this sequence next
connection is obtained by parallel or serial linking of new arc to obtained connection.
Distributions of random numbers of connectivity components are analyzed. These
distributions are considered intensively now. Central limit theorem is proved for these
distributions and parameters (mean and variance) of normal limit distribution are calculated.
1. INTRODUCTION
In the reliability theory parallel-serial connections play important role [1] [6] etc. These
connections are widely used in electrotechnics, in computer networks etc. A specific of these
connections is a possibility to calculate their reliability by algorithms with linear complexity by a
number of arcs.
Last years large interest is called to characteristics of networks sparseness. It means that
powers of nodes (a number of incident arcs) is bounded by some positive number (see [7] and large
bibliography in this article). Stochastic modeling and statistical processing of internet type networks
data showed that nodes powers have distribution with heavy tails [8]. Last circumstance makes
actual to consider parallel-serial connections which are free of this lack.
Last time a distribution of numbers of connectivity components in different random networks
are analyzed intensively now [9] [11]. In this paper numbers of connectivity components in
recurrent sequence of connections obtained by parallel or serial linking of new arc is considered.
For this sequence central limit theorem is proved and parameters of limit normal distribution are
calculated.
A problem to calculate a mean and mainly a variance of limit normal distribution in this
model is technically sufficiently complicated. In this paper it is based on central limit theorem for
discrete Markov chains [12] and on a construction of special and sufficiently fast algorithm of such
calculations.
2. MODEL DESCRIPTION
Consider the sequence
, 1, of ports defined recursively by a sequential or parallel
connection of new arc
to the port
. Denote a type of connection by || or accordingly.
Suppose that random variable
characterizes a type of the arc
connection to the port
and
put
= (
=),
||
= (
= ||) = 1
23
0<
< 1.
RT&A # 02 (25)
(Vol.1) 2012, June
Here random variable

(
characterizes a state of the arc
= 1) = (
in working state) = , (
= 0) = 1
= , 0 <
< 1.
The sequences of random variables { , 1}, { , 1} are independent and each of them
consists of independent and identically distributed random variables.
The port
with randomly working arcs is characterized by random vector ( , ) there
is an indicator of a connectivity between initial and final nodes of parallel-sequential
connection
and
is a number of connectivity components in
. Introduce auxiliary
random variables
=
,
=
+1 ,
(1)
=
(2)
then
(
)= (
=)(
)+ (
= ||)(
),
(3)
where ( ) is an indicator of an event .

3. LIMIT THEOREM FOR MARKOV CHAIN CHARACTERIZING CONNECTIVITY
OF PARALLEL-SERIAL CONNECTIONS
Denote
=
, then the sequence
= ( , ),
the states set = {( , ), = 0, 1, = 1, 0, 1} as follows
(
)= (
=)(
,1
)+ (
= ||)(
1, is Markov chain with
).
From the equalities (1) - (3) and the conditions 0 < < 1, 0 < < 1 we see that Markov chain
, 1, states are interconnected. Consequently from the central limit theorem for discrete
Markov chains with finite states set [12, chapters V,VI} there are normally distributed random
vector (0, ) with the dimension six and with zero mean and with covariance matrix and real
numbers ( ), , which do not depend on initial state
so that for any real ( ), ,
( )
( )
> ( ( ),
) ( (0, ) > ( ( ),
)),
(4)
Here ( ) =
( = ) and the inequalities are defined componentwise.
Introduce auxiliary numbers ( ), :
( , 0) = 0,
( , 1) = 1,
( , 1) = 1,
= 0,1.
From the formula (4) it is simple to obtain that there is normally distributed random variable
(0, ) with zero mean and with the covariance > 0 so that for any real
( )(
Using obvious equality

follows
( )
( )
( )) >
( )=
( (0, ) > ),
=
24
(5)
1, rewrite the formula (5) as
RT&A # 02 (25)
(Vol.1) 2012, June
( (0, ) > ),
>
( ) ( ).
(6)
Remark 1. A calculation of the vector ( ( ), )) and especially of covariance matrix in

the formula (4) is sufficiently complicated procedure [12, chapters V, VI}. So to define the mean
and the covariance we use following limit formulas
= lim
, = lim
(7)
which are corollaries of the formula (6) with special initial distribution of
4. CALCULATION OF LIMIT NORMAL DISTRIBUTION PARAMETERS
Choose random vector ( , ) = ( , ) which does not depend on random sequences

1}, { , 1} and satisfies the equalities
((
) = (1,1)) =
||
||
((
) = (0,2)) =
with ( = 1) , ( = 0) . Random sequence

Theorem 1. The equalities
=
,
=
(1
are true.
Proof. To define the constants ,
=

=
,
1, is stationary Markov chain.

(9)
(10)
from (7) we construct recurrent algorithm. Denote

|
(8)
) > 0
+2
=1
= 1),
|
= 1),
= 0),
= 0),
(11)
(12)
where
= 1,
= 2,
Using the formulas (1) - (3), (11) obtain for
=
=
=
||
+(
= 4.
1:
= 1,
||
+(
+ 1)
1)
+(
||
+ 1)
||
||
,
,
=1+ .
Then from (7) we obtain the equality (9).

And

=(
) 2
(13)
||
1, =
||
< 1.
so
+ 2
||
1
1
25
1 ,
RT&A # 02 (25)
(Vol.1) 2012, June
consequently

+ [
Begin now a calculation of
1 ],
1 ],
1.
(14)
. Using the formulas (1) - (3), (12) obtain for 1:

+
+(
||
+ 1)
+(
+2
||
||
,
+(
+2
+2
+ 1)
+2
||
||
+ 1)
)+
||
(1 + ).
So from (14) we obtain

=
+2
=
+2
+2
+2 (
(1 + )
( (1 + ) +
||
||
(1 + ) + 2
( 1)/2) 2
(1 + ) + 2
||
(1
),
(1 + ) =
1
1
(1 + )
||
= 1+3 .
Consequently
=
[1 +
(1 + )] + 2
(1
)+
Then from (7), (13) we have

(1 + ) = (1
= 1+
+ 2 ) > 0.
2
Theorem is proved.
Remark 2. From Remark 1 is possible to replace the condition (8) by more natural suggestion
(
) = (1,1) = ,
) = (0,2) =
so that the equalities (6), (9), (10) are true also.
REFERENCES
1.
2.
3.
Barlow R.E., Proschan F. 1965. Mathematical Theory of Reliability. London and New York:
Wiley.
Ushakov I.A. et al. 1985. Reliability of technical systems. Handbook. Moscow: Radio and
Communication (in Russian).
Riabinin I.A. 2007. Reliability and safety of structural complicated systems. Sankt-Petersberg:
Edition of Sankt-Petersberg university (in Russian).
26
4.
RT&A # 02 (25)
(Vol.1) 2012, June
Solojentsev E.D. 2006. Specific of logic-probability risk theory with groups of incompatible
events. Automatics and remote control, no. 7: 187-203 (in Russian).
5. Satyanarayana A., Wood R.K. 1985. A linear time algorithm for computing k-terminal
reliability in series-parallel networks. SIAM, J. Computing, v. 14: 818-832.
6. Ball M.O., Colbourn C.J., Provan J.S. 1995. Network Reliability. Network Models. Handbook
of Operations Research and Management Science, v. 7: 673-762.
7. Raigorodskiy A. 2010. Models of random graphs and their application. Proceedings of MPhTI,
v. 2, no. 4: 130-140 (in Russian).
8. Kumar R., Raghavan P., Rajagopalan S., Sivakumar D., Tomkins A., Upfal E. 2000. Stochastic
model for the web graph. Proc. 41-st Symposium on Foundation of Computer Science.
9. Sukach E., Ratobilskaia D., Demuskov A. 2011. Computer system of probability-algebraic
modeling of systems with many states. Mathematical mashines and systems, no. 3: 32-39 (in
Russian).
10. Sukach E., 2011. Approaches to decrease dimension of multi component and structurally
complicated systems with many states in estimate of their reliability. Problems of physics,
mathematics and technique. no. 3: 99-103 (in Russian).
11. Timashev A. 2011. Asymptotic expansions for distributions of numbers of connectivity
components in random images and decompositions. Discrete mathematics, v. 23, no. 2: 66-75
(in Russian).
12. Romanovsky V. 1949. Discrete Markov chains. Moscow-Leningrad: State Publishing House of
Technical-Theoretical Literature (in Russian).
27
RT&A # 02 (25)
(Vol.1) 2012, June
I. Ushakov UNIVERSAL GENERATING FUNCTION & OPTIMAL REDUNDANCY PROBLEMS
UNIVERSAL GENERATING FUNCTION & OPTIMAL REDUNDANCY

PROBLEMS
Igor Ushakov
San Diego, California, USA

1.
Introduction
The Method of U-functions, or the Method of the Universal Generating Function (UGF), was
introduced in [(1986) Ushakov, (1987) Ushakov] and later developed in [(1988) Ushakov; (1995)
Gnedenkon& Ushakov]. Actually this is a generalization and algebraic formalization of the wellknown Kettelles Algorithm [(1962) Kettelle]. In turn, Kettelles Algorithm, is a form of presentation
of convolution of discrete random variables. The method of U-functions is very convenient for
computerized calculations.
Last years, this method was significantly developed by G. Levitin and A. Lisnianski
2.
Briefly about Generating Function
Everybody knows that Generating Function (GF) is very convenient mathematical tool for
finding a convolution of discrete random variables.
Consider two non-negative discrete random variables X1 and X2 that are characterized by
discrete distributions
P{ X (1) x1(1) } p1(1) ,

P{X ( 2) x1( 2) } p1( 2) ,
(1)
(1)
(1)
( 2)
( 2)
( 2)
P{X x2 } p 2
P{X x2 } p2 ,
(1)
( 2)
P ( x)
and P
...
...
(
1
)
(
1
)
(
1
)
(
2
)
P{X xn (1) } pn (1) ,
P{X xn( 2( 2) ) } pn( 2( 2) ) ,
(1)
correspondingly where n(1) and n(2) are numbers of discrete realizations of values of each type.
If we are interested in the distribution of r.v. X = X1 + X2, we perform product of generating
functions, perform collecting terms, and get:
(1)
a ( 2)
( z ) (1) ( z ) ( 2) ( z ) pi(1) z ai p (j2) z j p k z Ak
1i n(1)
1 j n ( 2)
1k n
(2)
where Ak is a convolution of two r.v.s. ai(1) and ai(1) .

Thus, this transform suggests multiplication of polynomial coefficients and summation of
polynomial powers. The method of U-functions suggests a transparent and convenient method of
computerized solutions of various enumeration problems where variables are subjects to operations
beyond multiplication and summation, for instance, finding distribution of minimum, maximum,
geometrical summation, etc., depending on physical nature of the analyzed problem.
28
RT&A # 02 (25)
(Vol.1) 2012, June
Table 1. Examples of object interaction depending on physical nature of unit parameters

Name
( 1 , 2 )
( 1 , 2 )
SERIES
PARALLEL
PFFO (for hot redundancy)

Cost
Weight
el. Resistance
1 2
1 2
1 2
1 2
(1 1 ) (1 2 )
1 2
1 2
el. Capacity
(( 1 ) 1 ( 2 ) 1 ) 1
el. Conductivity
(( 1 ) 1 ( 2 ) 1 ) 1
pipeline capacity
min(1 , 2 )
random time to failure
min(1 , 2 )
...
...
number of different redundant units ( 1 , 2 )
(( 1 ) 1 ( 2 ) 1 ) 1
1 2
1 2
1 2
max(1 , 2 )
...
( 1 , 2 )
Here by symbol " " we denote interaction of parameters of various physical nature. In
particular, the method of U-functions can be effectively applied to solving the optimal redundancy
problem.
3.
Method of U-functions
Let us consider GF from another viewpoint. Each k-th discrete distribution one can represented
as a set of triplets:
S ( k ) ( p1( k ) , c1( k ) ,1), ( p2( k ) , c1( k ) , 2),..., ( p1( k ) , cn( k( k) ) , n(k ))
(3)
(k )
where p j
c ( k ) are the probability of failure-free operation (PFFO) of unit k with j-th variant of
and j
redundancy and the cost of this variant, correspondingly. The third component is the number of
redundant units of type k (or, in more general case, the ordering number of variant of unit k).
Indeed, product of two GFs is equivalent to Descartes interaction of two sets S(1) and S(2) ,
i.e. each triplet of set S(1) interacts with all triplets of set sets of S(2). Interaction of two triplets can be
conditionally written as follows:
( pi(1) , ci(1) , i (1) ) ( p (j2) , c (j2) j ( 2) )
(4)
In turn, interaction of triplets consists of interactions of its components that produce a new
triplet
(1)
i
Here interaction
p (j2) ; ci(1) c (j2) ; i (1) j ( 2) p ( 2*) , c ( 2*) , j ( 2 ) )
means product, operator
29
does summation, and operator
(5)
does
RT&A # 02 (25)
(Vol.1) 2012, June
union ( a vector with corresponding components), i.e.

p ( 2*) pi(1) p (j2) pi(1) p (j2) ;
c ( 2*) ci(1) c (j2 ) ci(1) c (j2) ;

( 2*)
(1)
( 2)
(1)
(6)
( 2)
j i j (i , j )
One can easily see that Descartes interaction of duplets that belongs sets S(1) and S(2) is
completely equivalent to product of two generating functions (1) ( z ) and ( 2 ) ( z ).
Analogously with the product of GFs one has to collect terms for getting the final set
S = S1S2.
Naturally, operator possesses commutativity property, i.e.
( a , b )= (b , a )
(7)
and associativity property, i.e.

(a, b, c) = (a ( b , c) )= ((a b ), c ).
4.
(8)
Using U-function for solving of optimal redundancy problems
Let us consider a series system consisting of n units, each of which has PFFO equals pk and
costs ck units. For increasing reliability of each unit, one can use redundancy of individual units.
Each unit k is represented by set of triplets
S k {R0( k ) ; C 0( k ) ; 0}, {R1( k ) ; C1( k ) ;1}, {R2( k ) ; C 2( k ) ; 2}, ..., {R s( k ) ; C s( k ) ; s}...
(9)
where s is the number of redundant units ( any natural number); C s(k ) is the total cost of s redundant
units (usually, a linear function of the number s); and Rs(k ) is the PFFO of unit k with s redundant units.
It is well-known that for loaded redundancy of group including one main and s identical redundant
units:
R s( k ) 1 (1 p k ) s 1 ;
and for an unloaded redundant (spare) units:
( t ) j
Rs( k ) k
exp( k t );
j!
0 j s
Now consider a general procedure of optimal redundancy with the use of U-functions. First
of all, take units 1 and 2 and arrange the Descartes interaction procedure between sets S1 and S2. In
our case
Ri(1) R (j 2 ) Ri(1) R (j 2 ) R K( 2*) ;
C i(1) C (j 2 ) C i(1) C (j 2 ) C K( 2*) ;
i j (i, j ) K ,
i.e. interaction between two numbers produces vector, containing numbers of redundant units of the 1st
and 2nd types..
30
RT&A # 02 (25)
(Vol.1) 2012, June
Here symbol * relating to the number means that this aggregated unit includes all previous
units.
At the next step of sets S1 and S2 interaction one takes aggregated unit 2* and unit 3:
R K( 2*) Rk( 3) R K( 2*) Rk( 3) ( Ri(1) R (j 2 ) ) R k( 3) R L( 3*) ;
C K( 2*) C k( 3) C K( 2*) C k( 3) C i(1) C (j 2 ) C k( 3) C L( 3*) ;
K k ( K , k ) (i, j, k ) L .
Vector L shows that in a series system of 3 units the 1st unit has i redundant ones, the 2nd
unit has j redundant ones and the 3rd units has k redundant ones.
This procedure continues until necessary final triplets will have been generated. Instead of
further abstract presentation of the procedure, let us turn to a simple illustrative numerical example.
The result of interaction is presented in the table below.
Example 1. Consider a series system of four units with parameters given in the table below.
Assume that hot redundancy is used for the system reliability improvement.
Table 2. System unit parameters
Unit-1 Unit-2 Unit-3 Unit-4
PFFO 0.6
0.6
0.7
0.7
Cost 3
5
3
5
In accordance with the description given above, the block diagram of the using U-functions
in this particular case can be presented as follows (see Figure 1).
Figure 1. Block-diagram of the solution procedure for Example 1
31
RT&A # 02 (25)
(Vol.1) 2012, June
Let us solve two problems of optimal redundancy:

(a) Direct problem: Find the optimal allocation of redundant units to reach required PFFO
level of the system equals to 0.97;
(b) Inverse problem: Find the optimal allocation of redundant units to reach maximum
possible PFO level under condition that the total cost of the system does not exceed 70 units of cost.
In this case, the UGF for each unit is defined by a set of triplets (Cost, PFFO, Number of
redundant units). Solution for the first step of the solution (interaction of set 1 and set 2) can be
presented in the form of the table below.
First of all, restrict ourselves with possible solutions for the Direct Problem: since the
system PFFO has to be not less than 0.97, it means that PFFO of each of four redundant groups has
to be not less than 0.97.
Since cost restriction equals 70 cost units, the total cost of redundant units in each redundant
group has to be not larger than, say 20-30.
Keeping this in mind let us construct the table with triplets obtained in the result of
interaction sets for Unit-1 and Unit -2.
Table 3. Result of interaction of UGFs for Unit-1 and Unit-2
S2
15
0.936
3
20
0.9744
4
25
0.9898
5
30
0.9959
6
35
0.9984
7
40
0.9993
8
S1
9
0.936
3
24
0.8761
(3; 3)
29
0.912
(3; 4)
34
0.9264
(3; 5)
39
0.9322
(3; 6)
44
0.9345
(3; 7)
49
0.9354
(3; 8)
12
0.9744
4
27
0.912
(4; 3)
32
0.9495
(4; 4)
37
0.9644
(4; 5)
42
0.9704
(4; 6)
47
0.9728
(4; 7)
52
0.9738
(4; 8)
15
0.9898
5
30
0.9264
(5; 3)
35
0.9644
(5; 4)
40
0.9796
(5; 5)
45
0.9857
(5; 6)
50
0.9881
(5; 7)
55
0.9891
(5; 8)
18
0.9959
6
33
0.9322
(6; 3)
38
0.9704
(6; 4)
43
0.9857
(6; 5)
48
0.9918
(6; 6)
53
0.9943
(6; 7)
58
0.9953
(6; 8)
21
0.9984
7
36
0.9345
(7: 3)
41
0.9728
(7; 4)
46
0.9881
(7; 5)
51
0.9943
(6; 7)
56
0.9967
(7; 7)
61
0.9977
(7; 8)
24
0.9993
8
39
0.9354
(7:3)
44
0.9738
(8; 4)
49
0.9891
(8; 5)
54
0.9953
(6;8)
59
0.9977
(7; 8)
64
0.9987
(8; 8)
In this table triplets that are dominated by others are marked with grey shadowing. One can
observe that dominating sequence occupies an area around diagonal of the table. This property can
be successfully used for minimizing the calculations: as soon as a dominated triplet appears below
this diagonal area, the further calculation in cells located below this cell can be stopped.
Analogously, if a dominated triplet appears upper this diagonal area , the further calculation in
cells located to the right from this cell can be also stopped. We will use this property in further
calculating.
Thus, the dominating sequence characterizing an equivalent Unit-2*is presented in nonshadowed area of table 1. On the basis of data for Unit-2*, we can construct an analogous table for
equivalent Unit-3* (see Table 4).
32
RT&A # 02 (25)
(Vol.1) 2012, June
Table 4. Result of interaction of UGFs for Unit-2* and Unit-3
S2*
S3
9
0.973
3
33
0.8524
(3; 3; 3)
36
0.8874
(4; 3; 3)
39
0.9014
(5; 3; 3)
41
0.9239
(4; 4; 3)
44
0.9384
(5; 4; 3)
12
0.9919
4
36
24
0.869
0.8761
(3; 3; 4)
(3; 3)
39
27
0.9046
0.912
(4; 3; 4)
(4; 3)
42
30
0.9189
0.9264
(5; 3; 4)
(5; 3)
44
32
0.9418
0.9495
(4; 4; 4)
(4; 4)
47
35
0.9566
0.9644
(5; 4; 4)
(5; 4)
50
38
0.9704 0.9442 0.9625
(6; 4) (6; 4; 3) (6; 4; 4)
52
40
0.9532
0.9717
0.9796
(5; 5) (5; 5; 3) (5; 5; 4)
55
43
0.9857 0.9591 0.9777
(6; 5) (6; 5; 3) (6; 5; 4)
58
46
0.9801
0.9881
(7; 5; 4)
(7; 5)
49
0.9891
(8; 5)
15
0.9976
5
18
0.9993
6
21
0.9998
7
24
0.9999
8
42
0.9098
(4; 3; 5)
45
0.9242
(5; 3; 5)
47
0.9472
(4; 4; 5)
50
0.9621
(5; 4; 5)
53
0.9681
(6; 4; 5)
55
0.9772
(5; 5; 5)
58
0.9833
(6; 5; 5)
61
0.9857
(7; 5; 5)
64
0.9867
(8; 5; 5)
61
0.9850
(6; 5; 6)
64
0.9874
(7; 5; 6)
67
0.9884
(8; 5; 6)
67
0.9879
(7; 5; 7)
70
0.9889
(8; 5; 7)
73
0.9885
(8; 5; 8)
Table 5. Final result of calculating
S3*
41
0.9239
(4; 4; 3)
44
0.9418
(4; 4; 4)
47
0.9566
(5; 4; 4)
50
0.9625
S4
15
0.973
3
56
0.899
(4; 4; 3; 3)
59
0.9164
(4; 4; 4; 3)
62
0.9308
(5; 4; 4; 3)
65
0.9365
20
0.9919
4
61
0.9164
(4; 4; 3; 4)
64
0.9342
(4; 4; 4; 4)
67
0.9489
(5; 4; 4; 4)
70
0.9547
25
0.9976
5
30
0.9993
6
35
0.9998
7
40
0.9999
8
69
0.9395
(4; 4; 4; 5)
72
0.9543
(5; 4; 4; 5)
75
0.9602
33
RT&A # 02 (25)
(Vol.1) 2012, June
(6; 4; 4)
52
0.9717
(5; 5; 4)
55
0.9777
(6; 5; 4)
58
0.9833
(6; 5; 5)
61
0.9857
(7; 5; 5)
64
0.9874
(7; 5; 6)
(6; 4; 4; 3)
67
0.9455
(5; 5; 4; 3)
70
0.9513
(6; 5; 4; 3)
67
0.9884
(8; 5; 6)
0.9568
(6; 5; 5; 3)
0.9591
(7; 5; 5; 3)
(6; 4; 4; 4)
72
0.9638
(5; 5; 4; 4)
75
0.9698
(6; 5; 4; 4)
78
0.9753
(6; 5; 5; 4)
81
0.9777
(7; 5; 5; 4)
84
0.9794
(7; 5; 6; 4)
(6; 4; 4; 5)
77
0.9694
(5; 5; 4; 5)
80
0.9754
(6; 5; 4; 5)
83
0.9809
(6; 5; 5; 5)
86
0.9833
(7; 5; 5; 5)
89
0.985
(7; 5; 6; 5)
87
0.9804
(8; 5; 6; 4)
92
0.9860
(8; 5; 6; 5)
85
0.9770
(6; 5; 4; 6)
88
0.9826
(6; 5; 5; 6)
91
0.985
(7; 5; 5; 6)
94
0.9867
(7; 5; 6; 6)
97
0.9877
(8; 5; 6; 6)
99
0.9872
(7; 5; 6; 7)
102
0.9882
(8; 5; 6; 7)
107
0.9883
(8; 5; 6; 8)
All calculations have been done with a simple Excel program.

Solutions of the problems above can be easily found from the last table. First time PFFO
exceed level of 0.97 when X=(6; 5; 5; 4) and the corresponding system cost is 78 cost units. The
inverse problem solution for restriction on the cost equals 70 cost units reaches when X=(4,4,5,3)
and corresponding PFFO is equal to 0.9547.
Notice that due to associativity property of U-functions it is possible to get the same solution
using another order of units interaction.
Figure 2. Second type of units interaction procedure
Remark. By the way, this example shows with transparency that one can consider not only
redundancy as a method of system reliability increase. For instance, one can consider a set of
variants of the units with various reliability and cost. Actually, Unit-2* and Unit-4* can be
considered as black boxes that are characterized by corresponding dominating sequences of
34
RT&A # 02 (25)
(Vol.1) 2012, June
triplets {Qs( k ) ; C s( k ) ; s} where s is just a number of variants of considered Unit-k.
Bibliography in chronological order

(This bibliography on the subject is presented in chronological order to show the steps of
the U-functions method development).
(1962) Kettele, J. D. , Jr. Least-coast allocation of reliability investment. Operations
Research, No. 2.
(1986) Ushakov, I.A. Universal Generating Function (in Russian). Engineering Cybernetics,
No. 5.
(1987) Ushakov, I.A. Universal generating function. Soviet Journal Computer Systems
Science, No.3.
(1987) Ushakov, I.A. Optimal standby problem and a universal generating function. Soviet
Journal Computer and System Science, No.4.
(1987) Ushakov, I.A. Solution of multi-criteria discrete optimization problems using a
universal generating function. Soviet Journal of Computer and System Sciences , No. 5.
(1988) Ushakov, I.A. Solving of optimal redundancy problem by means of a generalized
generating function. Elektronische Informationsverarbeitung und Kybernetik, No.4-5
(1995) Gnedenko, B.V., and Ushakov I.A. Probabilistic Reliability Engineering. John
Wiley & Sons.
(1998) Levitin, G., Lisnianski, A., Ben Haim, H., and Elmakis, D. Redundancy
optimization for series-parallel multi-state systems, IEEE Transactions on Reliability, , No. 2.
(1999) Levitin, G., and Lisnianski, A. Importance and sensitivity analysis of multi-state
systems using universal generating functions method, Reliability Engineering & System Safety, No.
65.
(2000) Levitin, G., and Lisnianski, A. Optimal Replacement Scheduling in Multi-state
Series-parallel Systems. Quality and Reliability Engineering International, , No. 16.
(2000) Ushakov, I.A. The method of generating sequences. European Journal of Operational
Research, No. 2.
(2001) Levitin, G. Redundancy optimization for multi-state system with fixed resourcerequirements and unreliable sources. IEEE Transactions on Reliability, No. 50.
(2001) Levitin, G., and Lisnianski, A. A new approach to solving problems of multi-state
system reliability optimization, Quality and Reliability Engineering International, No. 47.
(2002) Levitin, G. Optimal allocation of multi-state elements in linear consecutivelyconnected systems with delays. International Journal of Reliability Quality and Safety Engineering,
No. 9.
(2003) Levitin, G. Optimal allocation of multi-state elements in linear consecutivelyconnected systems. IEEE Transactions on Reliability, No. 2.
(2005) Levitin, G. The Universal Generating Function in Reliability Analysis and
Optimization. Springer.
35
A Bochkov, V. Lesnykh, N. Zhigirev DYNAMIC MULTI-CRITERIA DECISION MAKING METHOD FOR SUSTAINABILITY RISK ANALYSIS OF
STRUCTURALLY COMPLEX TECHNO-ECONOMIC SYSTEMS
RT&A # 02 (25)
(Vol.1) 2012, June
DYNAMIC MULTI-CRITERIA DECISION MAKING METHOD FOR

SUSTAINABILITY RISK ANALYSIS OF STRUCTURALLY COMPLEX
TECHNO-ECONOMIC SYSTEMS
Alexander V. Bochkov1, Valery V. Lesnykh, Nikolay N. Zhigirev,
Science Research Institute of Economics and Management

in Gas Industry (LLC NIIGAZECONOMIKA)
1
Corresponding author e-mail: a.bochkov@gmail.com
Abstract
The paper considers the characteristics of the functioning and sustainability of the
structurally complex techno-economic systems (SCTES) in terms of different types of risk. The
validity of the application to describe the behavior of this class of systems of semi-empirical
mathematical models, which are based on a vector description of the system states, using the criteria
approach for assessing the quality of its functioning, is demonstrated.
Under discussion is conceptual model for the interaction of the object and its environment,
allowing to estimate the "optimal" allocation ratio between the productive system and its
development potential.
The concept of non-formalizable threats for the sustainable functioning of this systems class
was introduced. Expert procedure to account non-formalizable threats in case of risk assessment
was proposed. For the construction of indicators for assessing the status the methods of quantitative
analysis based on the theory and multi-criteria utility was used. Multi-criteria utility as an indicator
of sustainability in the form of dimensionless complex hierarchy of indicators was proposed.
Computing for through the convolution of the primary indicators.
The hierarchical model proposed to calculate the integral index of multi-criteria preference
of one embodiment of the system over the other. Some results of case study are discussed.
Key words: structurally complex techno-economic systems, risk analysis, sustainability,
multi-criteria method
SCTES are characterized by distribution in space, big variety and interaction of objects
types, non-uniform structure of processing chains, unique conditions of influence of risks of the
various nature on objects of the subsystem and the system as a whole.
In the idea of situation management of SCTES principles of changeable (adaptive) behavior
in terms of possible risks and uncertainties are initially put. Presence of such risks generated by
different circumstances is capable to brake or change this or that way of movement, to force the
system to live under another scenario, different from all variety of plans shaped before.
If as sustainability of SCTES functioning to understand the plan performance of its
development with admissible variation on volumes and terms of problems performance then
situation safety management in this system is reduced to minimization of hazardous losses at
extraordinary situations and to carrying out of actions for their prediction. The success of such
tactics depends substantially on intuition and talent of management of the company, on its ability to
expect the possibility of weakly formalized threats outgrowth into notable risks and losses.
Under weakly formalized threats we understand here the threat for criticality estimation of
36
RT&A # 02 (25)
(Vol.1) 2012, June
which the development of original algorithm of the decision depending on a concrete situation is
required and for which uncertainty and dynamism of the initial data and knowledge can be
characteristic.
However, in the absence in enough of adequately estimated information necessary for
decision-making, tactics of adaptive management quite often turns to a continuous chain of the
"emergency" scenarios leading to disruption of controllability of all system. Hence, company
management should be engaged not only the current work bringing quite notable results which
utility is measured in economic factors, but also to care of creation the company condition
monitoring system and the world surrounding it, watching dynamics of internal and external threats
to its growth and development.
That is optimal control of SCTES aimed at reception of profit on its activity, consists in
ability to find balance of redistribution of the available resources (material, human, information) by
proprietor of the company between productive activity and maintenance of development
potential .
The elementary model illustrating the abovementioned and allowing to estimate "optimal"
proportions of resources distribution between "useful" (production) system and its development
potential is the model of interaction of developing object and its environment (Klykov, 1970;
Zhigirev, etc., 1983).
Let's present the activity of some SCTES, consisting of two subsystems (fig. 1).
Figure 1. Activity of productive system

The first, productive subsystem, brings profit, proportionally to quantity of received
1 (t )
,
resources ( t ) x (t ) with some positive resources increase speed coefficient
( t )
available in the system.
The second subsystem - development potential plays the role of the accelerator (retarder)
of resources reproduction speed in the system.
Herein (t) is proportion of resources distribution between productive subsystem and its
development potential, 1 (t ), 2 (t ) - coefficient of resources exchange intensity between
investigated system and some external in relation to it system in the process of coexistence.
37
RT&A # 02 (25)
(Vol.1) 2012, June
Actually the difference ( t ) 1 ( t ) 2 ( t ) is the share of resources deduced from the cycle of
reproduction in the form of losses of one kind or another, for example, of final consumption, taxes,
etc.
1 ( t )
In the elementary representation "potential" influence is the value of function

(t )
and coefficient (t ) for large-scale systems we consider independent of times in an explicit form as
constants. In this case system development is described by the homogeneous linear equation on a
variable x(t ) at parameters and
1 (t )
dx (t )
x( t )
( t ) x( t )
dt
( t )
(1)
The optimum proportion * (t ) between productive system and its development potential is
defined from the condition
1 *
max
(2)
*
*

At the natural assumption that is monotone function with saturation (fig. 2) there is a
1
.
simple way of its optimum definition, as
1
Figure 2.
According to fig. 2 it is clear that this optimum is reached in some point * , having quite
certain sense. So, if resources for development potential are allocated excessively much * ,
then means * are incorrectly withdrawn from the current reproduction and there is a situation
when efforts to studying and counteraction to numerous risks which the developing system can
never come across are spent.
The point 0 corresponds to the situation when all resources are spent exclusively for
growth of productive system. The potential of similar system is low because of constant losses
which it is possible to avoid if there is a potential for prediction of arising risks and struggle against
them.
38
RT&A # 02 (25)
(Vol.1) 2012, June
The segment 0, 1 shows that if the means allocated for studying and counteraction to
threats and risks are small, then return from similar researches and done actions less than the
resources allocated for them. Information gathering, research of internal and external threats on a
low level doesn't allow to receive an adequate estimation for improvement of decision-making
quality in most cases anyhow developing circumstances.
On the segment 1 , * the contribution to development potential starts to give positive
return, however only in the point 2 the level of "self-support" of expenses for development of
"potential" of system will be reached 2 0 .
Therefore it is expedient to consider this point as the point of "critical" position.
Decrease of potential to the level 2 threatens that "in accordance with the
circumstances" economically expediently there will be strategy of survival strategy of full
refusal of expenses for the decision of problems of prediction and anticipation of threats and risks
and reproduction maintenance only at the expense of escalating inefficient capacities in productive
subsystem 0 .
In spite of the fact that the conceptual model stated above is schematical, it, nevertheless,
provides guidance that threats and risks as a matter of fact are "antipotentials" of development, that
is they are retarders of speed of all system reproduction.
Since SCTES, as a rule, are non-uniform, they are subject to risks various by the nature and
on influence levels. The received expert estimations of optimum proportions, certainly, need
updating if to consider balanced development of the system consisting from many productively and
territorially connected subsystems.
The logic of optimum proportional development in this case also remains. Received
estimations should be considered only as "reference points" for the further researches, otherwise
struggle for escalating of development potential will be carried out only in those territories and only
in those productive-technological chains for which it by theoretical estimations is "economically
expedient" that will lead to destruction of integrity of the system (connectivity loss), withdrawal
from unified state and branch standards.
Let's suppose that the exit of investigated system runway from admissible corridor (a
component of the vector of functioning efficiency indicators) can be caused to four reasons:
a) owing to increase of importance of the problems put before the system to such level that
default of these problems at occurrence of extraordinary situation (and furthermore in a normal
mode of functioning) appears critical for system existence, up to necessity of its re-structuring as a
whole;
b) owing to system simplification or destruction at which locally arising extraordinary
situations are really capable to outgrow in events with large-scale losses under scenarios of cascade
type;
c) owing to dramatic or long deterioration of operating conditions of objects and subsystems
in one or several territorial formations formed, including, as a result of non-formalizable threats
increase;
d) owing to decrease in level of industrial and fire safety and (or) physical protection for
technological blocks and objects of various type.
It is offered for an estimation of extraordinary situation threat level in SCTES to use the
following hierarchical multi-criteria model (Russman, 1991).
Integrated risk of extraordinary situation R(r1 , ri ,, rn ) represents function of risks of
private extraordinary situations occurrence ri (i 1,, n) . The kind of dependence R on the
arguments gets out proceeding from conditions:
0 R(r1,, ri ,, rn ) 1 ;
(3)
R(0,,0,,0) 0 ;
(4)
39
RT&A # 02 (25)
(Vol.1) 2012, June
R(0, , ri , ,0) ri ;
(5)
0 R (r1 , ,1,, rn ) 1 for ri 1 irrespective of values of other arguments.
(6)
Continuous function R(r1 , ri ,, rn ) , meeting (3)-(6), has the following general view
R(r1,ri ,, rn ) 1 1 ri g (r1 ,ri ,, rn ) ,

i 1
(7)
where g (0,, ri ,,0) 1 .

If in special case g (r1, ri ,, rn ) 1 , then formula (7) is of the form
R(r1 ,ri ,, rn ) 1 1 ri
(8)
i 1
states the underestimated estimation of integrated risk from calculation that the stream of
extraordinary situations represents a mix of ordinary events taken from homogeneous, but differing
with values ri (i 1,, n) samples.
But for real systems risks, as a rule, are dependent.
Then we have
n 1
g (r1 , ri ,, rn ) 1
C r r
ij
ij
ij
(9)
i 1 j i 1
n 1
ij
1 Cij 0
i 1 j i 1
ij 0
(10)
ij 0
where Cij - coefficients of risks connection of i and j extraordinary situation; ij and ij - positive
coefficients of elasticity of replacement of corresponding risks, allow to consider the facts of risks
replacement, caused mainly by that simultaneously effective actions for decrease in all risks can't be
done owing to limitation of time and resources.
The current risks values ri (i 1,, n) , entering integral risks factor R are values changed in
time with various speeds (for example depending on the seasonal factor priorities of solved
technological problems essentially change).
Thereof classical calculation of risks equation leads to problems of combinatory complexity
on the initial data having objectively casual, uncertain, often qualitative (semiquantitative) nature.
The decision of problems of risks analysis becomes complicated also because nonformalizable threats can play considerable role.
For the account of these factors it is offered to form values ri as product of four
components:
ri ri( a) ri(b) ri( c) ri( d ) .
(11)
The component ri(a ) in (11) is estimated through categorizing of problems the performance
of which is cancelled or delayed owing to the arisen extraordinary situation (for example, in gas
supply systems categorizing can be defined through percentage distribution of categories of power
users, affected in case of extraordinary situation because of the termination of gas supply).
(b )
The component ri is estimated through maximum permissible losses (MPL) in
extraordinary situations at existing technological levels and materials (subjectively established)
calculations of such losses.
Before achievement the level of MPL ri(b ) can be considered as linear function
40
ri ( b )
Li
,
MPL
RT&A # 02 (25)
(Vol.1) 2012, June
(12)
where L i - current level of losses.

At exceeding of the level of MPL ri(b ) fixed with value 1.
The multiplier ri( ) is estimated as dimensionless value, is calculated under empirically
picked up statistical data about characteristics of objects in a binding to their territorial placing and
has the meaning of indicator of absolute vulnerability of object on which the scenario of i
extraordinary situation is initiated.
In general ri( ) can be considered as the indicator of environment aggression in which the
object functions.
For each territory owing to geographical factors, features of productive structure,
sociocultural, ethnic and other differences the construction of unique models of calculation
significantly found on personal assessment of the experts familiar with this specificity is required.
And a last, valuation of ri(d ) is in terms of ranging of objects types. It reflects quality of
relative "susceptibility" of objects of the set type on a wide range of external changes of the factors
defining ri(c ) . Values ri(d ) are used so that to result estimations of risks of extraordinary situations
initiated by events on objects of various types to a uniform scale.
The offered scheme of calculation of integrated index R is mainly intended for the
preliminary analysis of variants of system development on the basis of hierarchy of the indicators
characterizing all aspects of extraordinary situations including both estimations of consequences
ri(a) , ri(b ) , and estimations of causes ri(c ) , ri(d ) .
Lines of levels of R values allow to build borders of reaction zones to changes of all
spectrum of risks: R R (a "red" zone demanding change of the existing mode of functioning or
additional means for decrease of risks ri(c ) and ri(d ) ), R R R (an orange zone demanding
balancing of contract activity, carrying out of diagnostic and other actions), R R (a "green" zone
in which pertinently to speak about growth and the further development of the system, introduction
of new capacities and new risks connected with their occurrence).
It is obvious that between developing object what is any of SCTES and its development
potential, one of the components of which is the subsystem of safety (risk) management there
should be a balance. High-yielding system with the underestimated risk is doomed to inefficient
functioning owing to losses and on the contrary, the excessive safety concern leads to withdrawal of
resources from a cycle of reproduction.
Thus, for SCTES, having diversified multiphasic production and difficult space-territorial
topology of estimations having only economic character are unacceptable. Expediently complex
development of development potential control system which can be realized, for example, within
the limits of audit formly occurred crisis and precritical situations, estimation and generalization of
experience, development of system of indicators of early detection of threats to steady functioning
of objects (groups of objects).
Let's notice that during development of such indicators system for realization of situational
approach to management of the company, it is necessary to accept a number of "reconciliatory
agreements.
The first main agreement consists that preservation of integrity of system and knowledge
(information) about it is more significant, rather than economic success of separately taken
productive element or productive-territorial formation.
The second agreement: threats and risks are considered as factors, braking development
potential and, accordingly, use of the device of an estimation of the analysis and risks management
41
RT&A # 02 (25)
(Vol.1) 2012, June
without taking charge for risks corresponding to their competence by regulatory bodies is
impossible. Concealment of risks or their revaluationll become a subject of economic auction,
inappropriate in the conditions of approaching threats.
The third agreement: pure one-criteria approach when to every disovered risk (social,
economic, foreign policy) "cost" estimation of consequences of its realizations and (or) prevention
of scenarios of their expansion is offered, isn't universal.
Activity of any person separately, groups of people, labor collective of the company as a
whole many-sided and various, the use of multi-criteria approach with elements of "indistinct"
logic, with use (as far as data permits) detection device of the latent laws in conjunction with and
mutual strengthening of numerous factors therefore is the most pertinent.
The methodological approach taken as a principle offered method has advantage in
comparison with the cost approaches, expressed that multi-criteria utility "absorbs" in itself in
"share" participation all factors, but not just having cost expression ones (that, however, doesn't
allow to remove all uncertainties).
Multi-criteria utility is formed on the indicators having in the basis different dimensions,
units and scales of measurements which are easily arranged in the presence of computing resources
to specific features of investigated objects and risks generated by numerous factors at different
circumstances of place and time.
Classical schemes of multi-criteria analysis, based on linear and multiplicative convolutions
are successfully enough used in design and predesign analysis, at the decision of some problems of
situational management in marketing, at risks estimation of continuation or termination of research
and developmental works but as SCTES is dynamic system, it is offered to use more developed
model of dynamic multi-criteria analysis which will allow to include the situations generated socalled non-formalizable threats and risks-factors into consideration.
References
Klykov, Yu. I. (1974) Big systems situational management. Energiya, Moscow 130 p.
Zhigirev, N.N. (1983), Vorobiev, .I., Kuzmin, I.I., Reznichenko, V.Yu. Model of
interaction of developing object with environment. // Preprint Institute of Atomic Energy 3799/16. 69 p.
Russman, I.B. (1991), Kaplinskiy, .I., Umyvakin, V.. Modeling and algorithmization
non-formalizable problems of the best variants of system choice. Voronezh: VSU - 168 p.
42
Fedol Ghazi, Hamouine Abdelmadjid, Zouaoui Chikr el Mezouar BEAM BEHAVIOUR UNDER MONOTONIC LOADS
RT&A # 02 (25)
(Vol.1) 2012, June
BEAM BEHAVIOUR UNDER MONOTONIC LOADS

Fedol Ghazi
Bechar University B. P. 417, 08000 Bechar, Algeria

E-mail: gfedol@yahoo.fr
Pr: Hamouine Abdelmadjid

E-mail: ahamouine@gmail.com
Dr: Zouaoui Chikr el Mezouar

E-mail: chikrtime@yahoo.fr
Tel: 213 49 81 55 81/91; Fax: 213 49 81 52 44
Abstract
In this work we interest to study the beam behavior under monotonic loads in four point, to
improve the mechanical properties of a concrete beam fiber and establish an identification card of
the new concrete beams were comparing these beams witnesses.
Keywords: matrix materials, reinforcements, fiber characterization, charge arrow
monotonous beam, ductility.
1. Introduction
From earlier the research were based on how to reinforce materials with fiber plantin order
to increase there mechanical resistance and improve there stability. In the past they used the
Torchis was of clay reinforced with straw put in place by compression.
(Kriker et al 2005) used the date palm fibers as building blocks in cement matrix
composites. They showed that the increase in length and percentage of fibers improve the flexural
strength and post-elastic hardness of the composite, but decrease the compressive strength.
The work presented here is an analysis of the behavior of fiber-reinforced beams armed as
bellow:
Study the effect of the incorporation of fibers in a cement matrix

Study the recovery of fibers as reinforcements
Determine the increase in flexural strength under monotonic loading increasing
Observe the mode of cracking
The first set of standardized test used to determine the compressive strength and tensile
strength in bending.
The second set concerns the 4-points bending beams (15x10x60 cm) that will be subjected
to a monotonically increasing load until failure.
43
RT&A # 02 (25)
(Vol.1) 2012, June
2. Presentation of the materials tested

2.1 Materials in base
For any current use of concrete fiber in the building was used aggregate crushers in the
region with a dosage of 350 kg/m3 cement
2.1.1 Cement
The cement used is a type of cement CEMII 42.5 NA 442, physical properties are given in
the table below
Table 1: The results of standard tests carried out on this cement are given in Table 2
Test
start
2h50mn
End
4h06mn
Chatelier expansion
specific surface
(Blaine)mm2/g
consistency normal
hot
2.90
3891
26.91
cold
1.65
Table 2: Mechanical resistance cement (bar)

tests
Compression (b)
Flexion
(b)
Age
7days
266.4
58.3
2days
143.2
35.7
28 days
433.1
77.9
2.1.2 Mixing water

The water used in mixing is the tap (dam Djorf ettorba), the results of physicochemical
analysis are as follows:
Table 3: Results of analysis of water-chemical physic
PH
8,13
Matter in
suspension
Null
Chlorides
Mg/l
234,3
Sulfates
Mg/l
123,02
Residue sec
105 C0
800,00
Conductivity
25 C0 s/cm
0,93
2.1.3 Aggregates
The gravels are Petro graphically micritic limestone partly dolomitize sandstone, sand
consists mainly of these proportions: Silica and limestone following rigorous testing of these
materials are characterized as follows:
Table 4: Particle size analysis has given us the following composition for a dosage of 350Kg/M3
Designation
Class d/D
Mass volumique Absolute
Apparent
Surface properties%
Equivalent sable %
Product
Sand 0/3
2,5 t/m3
1,85 t/m3
Gravel 3/8
1,53 t/m3
1,0 %
67 %
44
Gravel 8/15
2,66 t/m3
1,41 t/m3
0,80 %
RT&A # 02 (25)
(Vol.1) 2012, June
Coefficient LA
Coefficient M.D.E
Nature
Fineness modulus
21 % 40 %
17 % 35 %
calcareous silico
1,95
Sandy dolomitic limestone
Table 5 Composition of 1 m3 of concrete

%
Ciment
Sand (0/4)
Gravel(3/8)
Gravel(8/15)
Concrete sec
Water
13,85
37,15
16
33
100
/
Volume
absolute
112,9
302,7
130,4
268,
815
185
Mass volumique
absolut [t/m3]
3,1
2,5
2,66
2,66
/
1
Mass [Kg] Mass volumique

Apparent [t/m3]
350
1
756,7
1,85
346,9
1,53
715,2
1,41
2168,8
/
185
1
According to this composition were obtained concrete firm for a 1cm subsidence ratio (E / C
= 0.53), the collapse is obtained for a de6cm (E / C = 0.60), for a Mix = 1.5%; Mix in means the
ratio of fiber relative to the aggregates,. Here is the composition of the concrete practice
Table 6: Composition of Concrete Practice Mix for 1 m3 of 1.5%
%
Ciment
Sand (0/4)
Gravel(3/8)
Gravel(8/15)
Concrete sec
Water
14,29
36,71
16
33
100
/
Volume
absolute
112,9
290
126,4
260,7
790
210
Mass volumique
absolut [t/m3]
3,1
2,5
2,66
2,66
/
1
Mass
[Kg]
350
725
336,2
693,46
2104,66
210
Mass volumique
Apparent [t/m3]
1
1,85
1,53
1,41
/
1
The identifications of the various physical and mechanical aggregates showed conformance
to specifications of standard NF P18 301 (Georges 1990), Also, aggregates (G1, G2 and S) shows
no abnormality in their grading curves, and Based on these results in these fractions fall into classes
0 / 3, 3 / 8 and 8/15selon NF P18 560.
2.1.4 Fibre Plant
We use the leaves of palms date palms of the type of Taghit Oisis "Fegousse" to the
saturated state antecedent research have shown a clear difference between concrete and full of dry
fiber Fiber identification:
A)
Property-mechanical: test is used for axial tractions a sample of 20 fiber
dimensions 0.35 mm thick, with a length of 30mm and a width of 6 mm,
45
RT&A # 02 (25)
(Vol.1) 2012, June
Figure 1: Tensile force in function of the elongation

The results are given in the following table:
Table 7: characteristics Mechanical of the fiber
Resistance
traction (Mpa)
114
B)
to The longer %
6,5
Coeff of absorption %
132
Modulus of elasticity
(Mpa)
17,58
Adhesion Fibre Concrete
The Essay adhesion pull-out test To evaluate the bond stress fiber matrix, we use the direct
test method (pull-out), so we anchor the fiber in a cement matrix and then applying a tearing force
on the fiber during the test is fixed fiber dimensions (length 150mm, width 7 mm, thickness
0.55mm), the only variable is the anchorage length we take the length and 3 cm respectively
1.5,2,2.5 are used specimens
Figure 2: Breakout force versus slip

Analysis of the pull-out curve:
Pcr: the end of the elastic behavior, it is a critical shift cr (elastic zone)
Pmax is the maximum force before detachment, which corresponds to a shift max
0: the shift corresponding to the total detachment
After several direct axial pullout tests in the following figure summarizes the adhesion
Force based on the anchorage length
46
RT&A # 02 (25)
(Vol.1) 2012, June
Figure 3-curve of the maximum force based on the anchorage length

3. Experimental Method
Cont de compression (Mpa)
3.1 Compression test: Compression tests are carried out on cubes of 10x10x10 cm 3 after
28 days
Compression Test (Mpa)
2
3
Mix %
Figure 4: test on cubes of 10x10x10 cm 3 after 28 days
3.2 Tensile test 4-point bending

They summers 7X7X28 performed on samples after 28 days
Cont Traction (Mpa)
Tensille Test (Mpa)
2
3
Mix %
Figure 5: 7X7X28 performed on samples after 28 days
47
RT&A # 02 (25)
(Vol.1) 2012, June
After compression tests and tensile Note that there is a considerable drop constraint (over
50%) are increasing the fiber content (Mix) as is consistent with studies Similar [2]
Reinforcement
For the longitudinal reinforcement is used under the terms of 4T10 non fragile (Section
A.4.2 of BAEL)
For transverse reinforcement conditions of shear imposes 6 spacing of 7 cm at the
supports.
We chose the minimum percentage to see the behavior of concrete in tension and
compression.
4. Results of static tests:
4.1 Charge-arrow diagram
Analysis by the Mix
The results obtained can be classified our sample into 2 categories, the first consisting of the
control concrete, and those of Mix 0.5%, 1.0%, the second category that of1, 5%, 2.0% and 3.0%.
In the first category we find that the curves of the three components are combined in a first zone
(zone without degradation) is the elastic zone and an area with a slight shift of the beginning of
cracking (phase elasto-plstique) and finally a plastic phase that ends in failure. The tensile strength
in bending of fiber-reinforced beams is 1.21% times more than the control concrete, for against the
influence of fiber length appears in the arrow registered; it to an arrow of 2.46 mm for a fiber length
of 6 cm and 2.54 mm for the fiber to 4cm Mix 0.5% in the second category there is clearly brought
on beam ductility (an arrow of 3.8 mm for the Mix 3% to 6 cm fiber) against it by a break for a load
of 55 KN of course this is a break due to compression
Module de young (Mpa)
4.2 Module of elasticity

The modulus of elasticity is a constant mechanical stress of materials, is given by the slope
of the first part of the diagram (, %) is an instantaneous modulus calculated first threshold. We
note that for Mix 0.5% fiber is 4 cm in the same behavior materials by the cons there is a drop in
module for the other
4 cm
6 cm
Mix %
Figure 6 : Module of elasticity
48
RT&A # 02 (25)
(Vol.1) 2012, June
Indice de ductulit
4.3 Ductility Index

Among the important advantages of fiber concretes that improve the ductility of materials,
which plays an important role in seismic areas (avoids sudden destruction of a building) is an
important research in the future, so we see that the contribution of fibers to dramatically improve its
ductility than the control concrete (Mix for 0.5%, 1% and 1.5%)
2
3
Mix %
Figure 7 : Ductility Index

5. Conclusion
The incorporation of the fibers of date palms wet brings a significant improvement in the
behavior of beams in four point bending ductility and its bearing capacity, noting that actual
performance summers were obtained without requiring a particular choice of cement (CPJ 45)
commonly used with a dosage of 350 kg/m3 or natural aggregates after the crusher, recalling that
the compressive stress of the control concrete was 21.5 MPa, our goal was to explore the possibility
of using fiber in the current building, after analyzing the different steps it is clear that the Mix 0.5%
of the fibers of 4 cm gives the best results, the tensile strength was 1.21% higher than the control
concrete by cons ductility was lower than the control concrete, the increase in fiber actually
increased its plastic deformation for the Mix 0.5% (4 cm) is the strain at break is greater than that of
control concrete by 50% but in a sharp decrease its compressive stress and tensile
Results of tests:
Witness
Mix 0,5 %
0,011
Constraint Mpa
Constraint Mpa
0,03
0,057
0,078
0,099
0,12
Witness
4cm
6cm
0,144
Relatif deplacement %
Relatif deplacment %
0,178
49
RT&A # 02 (25)
(Vol.1) 2012, June
Mix 1,5 %
Constraint Mpa
Constraint Mpa
Mix 1,0 %
witness
4cm
6cm
witness
4cm
6 cm
relatif deplacement %
Constraint Mpa
Mix 3%
witness
4cm
6cm
References:
[1]
Abbani Said. Caractrisation et durabilit des btons renforcs de fibres vgtales
de palmier dattiers.2003.Mmoire de magister, Universit dOuargla.
[2]
Khenfer.M.M, MorlierP. 'Effet de la longueur des fibres sur les proprits
mcaniques des ciments renforcs de fibres cellulosiques'. Materials and structures (24)(1991)185190.
[3]
A. Kriker, G. Debicki, A. Bali, M. M. Khenfer, M. Chabannet Mechanical
properties of date palm fibres and concrete reinforced with date palm fibres in hot-dry climate
Cement and Concrete Composites, Volume 27 (2005) 554-564
[4]
A.Kriker, A.Bali, G. Debicki, M. Bouziane, M. Chabannet Durability of date palm
fibres and their use as reinforcement in hot dry climates Cement and Concrete Composites Volume
30 (2008) 639-648.
[5]
Nadia Benmansour Etude des performances de produits renouvelables et locaux
adapts aux applications de l'isolation thermique dans le btiment -Thse de Magister Universit
Hadj Lakhdar, Batna (2011)
[6]
G.Dreux et J.Festa Nouveau guide du bton et de ses constituants 8iem dition
Eyrolles 1998
[7]
Cesar Juarez, Gerardo Fajardo Caracterisation microstructurale des fibres
naturelles pour des materiaux composites a base de ciment Canadian Journal of Civil Engineering
March 1, 2009.
50
RT&A # 02 (25)
(Vol.1) 2012, June
[8]
Guide de bonnes pratiques pour lutilisation des fibres dans le bton -Ralis par
lassociation bton Qubec en collaboration avec les manufacturiers de fibres mtalliques et de
fibres synthtiques-janvier 2005.
[9]
Hachani Imane et Khelif Rima Etude de leffet de la gomtrie des fibres
synthtiques sur les proprits physico synthtiques sur les proprits physico-mcaniques des
mortiers sable calcaire mmoire dingnieur dtat (2009) Universit Ammar Telidji Laghouat.
[10] Mourad Chikhi Mtrologie et Modlisation des Transferts dans les Composites
Naturels Faible Cot pour lIsolation Thermique dans le Btiment en Algrie Thse de Magister
Universit Hadj Lakhdar, Batna (2010)
[11]
Pascal Casanova Bton renforces de fibres thse de Doctorat l'Ecole Nationale
des Ponts et Chausses 26 juin 1995.
[12]
Ali Kaflou Etude du comportement des interfaces et des interphases dans les
composites a fibres et a matrices cramiques thse de Doctorat LInstitut National des Sciences
Appliques de Lyon, 20 Mars 2006.
[13] J.M.Alwan , A.E.Naaman et W.Hansen Pull-out Work of Steel Fibers From
Cementitious Composites :Analyse Investigation Department of Civil Engineering, 2340 G.G.
Brown Building, College of Engineering, University of Michigan, Ann Arbor, MI 48109-2125,
USA, 24 may 1991.
[14] H.Houari. Contribution ltude du comportement du bton renforc de fibres
mtalliques soumis
laction des charges maintenues et cycliques. 1993. INSA Lyon.243 p.
[15]
Valade. T. 'Comportement des btons arms de fibres' ; Mmoire de D.E.A,
L.C.P.C, (Dcembre 1987).
[16] Amina Djoudi, Mohammed Mouldi Khenfer, Abderrahim Bali Etude dun nouveau
composite en platre renforce avec les fibres vgtale du palmier -1st International Conf erence on
sustainable Built Environment Infrastructures in Developing Countries ENSET Oran (Algeria)
Octobre 12-14, 2009.
[17] Richane SACI Prise en compte de lorientation des fibres dans une poutre continue
matrice de bton -Thse de Magister Universit Mohamed Boudiaf Msila (juillet2004).
51
G. Tsitsiashvili, M. Osipova BASIS OF CONTINUUM APPROXIMATION IN MODELS OF GROWING RANDOM NETWORKS
RT&A # 02 (25)
(Vol.1) 2012, June
BASIS OF CONTINUUM APPROXIMATION

IN MODELS OF GROWING RANDOM NETWORKS
G. Tsitsiashvili*, M. Osipova**
* IAM FEB RAS, Vladivostok, Russia, e-mail: guram@iam.dvo.ru
** FEFU, Vladivostok, Russia, e-mail: mao1975@list.ru
ABSTRACT
Accuracy asymptotics for differences between prelimit and limit distributions of nodes
powers in models of random growing networks are constructed. A rate of a convergence in these
relations is power. Obtained formulas allow to ground the continuum approximations for considered
models.
1.
INTRODUCTION
Main aim of this paper is to estimate a rate of a convergence to limit distributions in models
of growing random networks [1]. One of the most convenient methods to define limit distribution of
node exponent is the continuum approximation [2], [3]. It is based on asymptotic behavior of
considered distribution when a number of steps tends to the infinity. In [4, p. 124] it is marked that
``a problem of strict formal description of statistical ensemble of random networks with fixed
distribution of nodes exponents is not yet solved``. An absence of correct mathematical base of the
continuum approximation makes susceptible results obtained using this approximation in spite of its
calculation convenience.
In this paper a ground of the continuum approximation in a calculation of limit distributions
for main models of growing random networks is made. Here the model of growing exponential
network, the model of Barabasi and the model of Dorogovtsev are analyzed. Exact asymptotic of a
difference between prelimit and limit distributions when a number of steps tends to the infinity is
obtained and this asymptotic is power. Such results are based on a construction of recurrent
relations for prelimit distributions of nodes exponents distributions and on asymptotic series in
these relations.
2.
EXPONENTIAL NETWORKS
Consider a model of growing exponential network in which new node is connected with
each existing node with equal probabilities. Denote p k,s, t the probability that on the step t 1
the node s,1 s t, is connected with k arcs of non oriented graph of exponential network. Then
k is called the power of the node s . In [2] the following relations are obtained (here ij is the
Kroneker index):
p k 1,s, t t 1 p k, s, t
p k, t, t k1 , p k,s, t 1 =
, k 1.
t
Designate
1 t
P k, t p k, s, t
(1)
t s1
52
RT&A # 02 (25)
(Vol.1) 2012, June
then the recurrent relation

t 1 P k, t 1 P k 1, t t 1 P k, t k1 , k 1 , t 1 ,
is true and P 1,1 1 , P k, t 0 , k t 1 . So for k 1
P 1, t 1
t t t 1 P 1, t
t 1 t
t t 1 t 1 t 2 P 1, t 1
t 1 t
...
1 ... t 1
, t 1,
t 1 t 2
(2)
and for k 2
P k, t 1
t P k 1, t t t 1 P k, t
t 1 t
t P k 1, t t 1 P k 1, t 1 t 1 t 2 P k, t 1
...
t 1 t
t
1
jP k 1, j .
t 1 t j1
(3)
From the last formula
1 ... t 1
, t 1,
2t t 1 4
(4)
2 ... t 1
1

, t 1.
4t t 1 8 4t t 1
(5)
P 2, t 1
P 3, t 1
Denote f k t P k, t 1 2 k , k 1 , then from the formulas (2) (5) we have
f 1 t f 2 t 0 , f3 t ~
1
, t .
4t 2
(6)
Theorem 1. For t the relations

C ln k 3 t
1
fk t ~ k 2
, Ck
, k 3,
(7)
4 k 3 !
t
take place.
Proof. For k 3 the formula (7) is a corollary of the formula (6) then by an induction we
obtain from (3) that
t
1
1
1
f k 1 t P k 1, t 1 k 1
jP k, j k 1
t t 1 j1
2
2
t 1
t
1
1
1
C k ln k 2 t
f
j
jf
j
1
~
, t .

k 1
k
k
t t 1 j1 2 k
t t 1 j1
2
k 2 t2
The formula (7) is proved.

Remark 1. In all sections of this paper the continuum approximation is based on the limit
t P k 1, t P k, t t f k t f k t 1 0 , t .
In this section the formula (8) is a corollary of the formula (7).
53
(8)
3.
RT&A # 02 (25)
(Vol.1) 2012, June
MODEL OF BARABASI-ALBERT
Consider Barabasi-Albert model of growing network [1] in which new node is connected
with each existing node with probability proportional to a power of existing node. Denote p k,s, t
the probability that on the step t 1 the node s, 1 s t, is connected with k arcs of non oriented
graph of Barabasi-Albert network. In [2] the following relations are obtained
k 1
k
p k, t, t k1 , p k,s, t 1 =
p k 1,s, t 1 p k,s, t , k 1 .
2t
2t
From the formula (1) we have
k 1
k
P k 1, t t P k, t k1 , k 1 , t 1 ,
t 1 P k, t 1
2
2
P 1,1 1 , P k, t 0 , k t 1 , P 0, t 0 , t 1 .
Analogously with (2), (3) it is not difficult to obtain
t t
1
1
P 1, t 1
1 1 , t 1 ,
(9)
t 1 j1s j 2s
t
t
k 1 t
k
P k 1, j 1 , k 2 , t 1 , 1 .
2 t 1 j1
2s
t 1
s j1
Lemma 1. For A 0 the equalities
t t
t
1 A t
A tA
, t
,
1
2
s 1 A 1 A 1 A
t 1
j1 s j
P k, t 1
t
t
t
A t 1
1
s 1 A 1 A 1 A
j1 s j1
(10)
(11)
(12)
are true. Here z is the gamma function.

t t
A
Proof. Denote S t 1 then
s
j1s j
1
1 A
S 1 1 A
,
1 A 1 A 2 1 A
consequently the formula (11) for t 1 is proved. Suppose that this formula is true for t and prove
it for t 1
t 1 t 1
A t 1 A
S t 1 1
S t 1
s
t 1
j1 s j
t 1 A
t
t 1
t 1 A t A
.
t 1 1 A 1 A 2 1 A
t 1
1 A 2 1 A
The relation (11) is proved for all natural t . The formula (12) may be proved similar.
Remark 2. In left sides of the formulas (11), (12) we have functions defined for A 0 ,
t 0 and on right sides - the gamma functions which may be non defined for some A 0 , t 0 .
But ratios of these gamma functions in such points may be redefined using limit transition to these
points.
54
RT&A # 02 (25)
(Vol.1) 2012, June
Designate
4
, k 1.
(13)
k k 1 k 2
Theorem 2. For t and k 1 the following relations are true
t 3 / 2
fk t ~
.
(14)
3
Proof. It is known [5, subsection 1.18] that for A 0
1 A t A
(15)
t
~ t ,t .
t 1
Assume that k 1 then the relation (14) is a corollary of the formulas (9), (10), (15) for
A 1/ 2
2
t
1 t 1/ 2
1
~
f1 t P 1, t 1 1
2
3
t 1
3/ 2
3/
2
3/
2
1/ 2
3/ 2
4t
t
~
, t .
9 3 / 2 t 3
f k t P k, t 1 k , k
Suppose that the formula (14) is true for k 1 and prove it for k , k 1 . Represent f k t in
the form f k t a k t b k t where
ak t
t
k 1 t
k
k 1 1 k ,
2 t 1 j1
2s
s j1
t
k 1 t
k
f k 1 j 1 1 .
2 t 1 j1
2s
s j1
Consequently from the formulas (12), (13) for A k / 2 we have
t
t
1
4
k
1
ak t
1
k k 1 2 t 1 j1 s j1 2s k 2
1 t 1
t
4
1

k k 1 2 t 1 1 k / 2 1 k / 2 1 k / 2 k 2
bk t
4 t
k t 1 k / 2
~
, t .
2 t 1 k k 11 k / 2 1 k / 2
1 k / 2
t 3 / 2
, t , and from the gamma function
3
properties and from the formula (15) for A k / 2 we obtain
t
k 1 t
k 1 t
jk / 2
t 3/ 2
bk t
~
d
j
~
, t .
f k 1 j 1
2 t 1 j1
j 2 t 1 1 3 j3/ 2 t k / 2
3
Then asymptotic relation (14) is true for arbitrary natural k .
From the induction assumption f k 1 t ~
55
4.
RT&A # 02 (25)
(Vol.1) 2012, June
MODEL OF DOROGOVTSEV
Consider Dorogovtsev model of growing network in which new node is connected with each
existing node with the probability proportional to a sum of its power (a number of arcs incoming to
existing node) and some constant a 0 . Here a 0 is model parameter. Denote p k,s, t the
probability that on the step t 1 the node s, 1 s t, has the power k . In [2] we obtain the
relations
k 1
ka
p k, s, t 1
p k 1,s, t 1
p k,s, t , p k, t, t k0 , k 0 .
t 1 a
t 1 a
Designate
1 t
P k, t p k, s, t
t s1
then
1
P k, t a 1 t 1 k a P k 1, t k 1 a 1 a k0 ,
P k, t 1
t a 1
P 0,1 1 , P k, t 0 , k t , P 1, t 0 , t 0 .
Analogously with the formulas (9), (10) it is not difficult to obtain
1
a t 1
a t 1 t 1
a
P 0, t 1 1
1
1
,
t a 1 j1 j a 1 s1 js j a 1
P k, t 1
t 1
k 1 a t
ka
P k 1,s 1
, k 0.
t a 1 s 1
a 1 j
js
(16)
(17)
Denote A k k a / a 1 , k 0 ,
k 1 a
1 2a k a
a k 2 2a
Theorem 3. The formulas
f k t ~ Ck t 1 A0 , t , k 0 ,
C0
1
2
1 A0 1 A0
, f k t P k, t 1 k , k 0.
C k 1 a
A0
, Ck k 1
, k 0,
1 A0
a 1 A k A 0
are true.
Proof. From the formulas (11), (15), (16) for A A 0 we have
t 1
1 A t 1 t 1 A 0
f0 t P 0, t 1 0 1 0
2
t 1 A 0
1 A0
A
1
0
0
1 a t 1 C0
~ C0 t 1A0 , t .
1 2a
t
For k 0 we seek f k t P k, t 1 k in the form f k t a k t b k t ,
t 1
A
k 1 a t
k 1 1 k k ,
t a 1 s 1
j
js
t 1
A
k 1 a t
bk t
f k 1 s 1 1 k .
t a 1 s 1
j
js
ak t
56
(18)
RT&A # 02 (25)
(Vol.1) 2012, June
Then from the formula (15) for t and A A k we have
ak t
(k 1 a) k 1 t 1 t 1 A k
1
t a 1
j
s 1 js

1 k

k 1 a k 1 t Ak 1
k 1 a k 1 t 1
~
.
2
t a 1 1 A k 1 A k
1 A k 2 a 1 1 Ak
From the induction assumption f k 1 t ~ C k 1t 1 A and from
the formula A k A 0 and
from the formula (15) for A A k we obtain

bk t
t 1
A
k 1 a t
f k 1 s 1 1 k
t a 1 s 1
j
js
C k 1 k 1 a
1 A 0
t 1
k 1 a t
~
f k 1 s 1
s 1
t a 1 s 1
C k t 1 A0 , t .
t
a 1 A k A0
Consequently asymptotic relation (18) is proved for arbitrary natural k .
Remark 3. A consideration of Dorogovtsev model [2] is connected with its wide application
to modern models of growing random networks. For small a this model gives sufficiently simple
and convenient description of the Internet network with power distribution network nodes
exponents
1 a 1 2a , k
k ~ Dk 2a , D
a
with the parameter close to two [6].
REFERENCES
1. Barabasi L.A., Albert. 1999. Emergence of scaling in random networks. Science, v. 286.
509-512.
2. Dorogovtsev S.N., Mendes J.F.F. 2002. Evolution of Networks. Adv. Phys., v. 51, no 4.
1079-1187.
3. Ginzburg S.L., Nakin A.V., Savitskaya N.E. 2009. Influence of complex network
structure on properties of its dynamical processes. Letters to JETPh, v. 90, no 12. 873-878. (In
Russian).
4. Evin I. 2010. Introduction to theory of complicated networks. Computer researchs and
modeling, v. 2, no 2. 121-141 (In Russian).
5. Beitman G., Erdeii A. 1965. Supreme transcendental functions. Moscow: Science. (In
Russian).
6. Raigorodskiy A. 2010. Models of random graphs and their application. Proceedings of
MPhTI, v. 2, no. 4. 130-140 (in Russian).
57
A. Blokus-Roszkowska, L. Smolarek - COLLISION RISK ESTIMATION FOR MOTORWAYS OF THE SEA
RT&A # 02 (25)
(Vol.1) 2012, June
COLLISION RISK ESTIMATION FOR MOTORWAYS OF THE SEA

A. Blokus-Roszkowska, L. Smolarek
Gdynia Maritime University, Gdynia, Poland

e-mail: ablokus@am.gdynia.pl
ABSTRACT
The Motorways of the Sea is rather a new concept and, thus, is still in the process of
development by the European Commission and in the Baltic Sea Region. The Baltic Sea Region is
one of the most dynamic growth areas. Due to this fact the new ideas and technologies are needed
to optimize the sea transport system. In the paper the simulation model for the system safety state
evaluation is presented. The simulation program can constitute a base for decision-support tool, on
the level of safety management, especially to optimally plan the safety transport system.
INTRODUCTION
The Motorways of the Sea is (MoS) still in the process of development by the European
Commission in the Baltic Sea Region. Cargo transported in containers and trailers have increased
rapidly, also oil tanker traffic has seen a noticeable increase in the Baltic. In the Motorways of the
Sea the focus should be put on the safety of shipping. Sea motorways main elements can be point
out as a part of the European transport corridors network, see Paulauskas & Bentzen (2007).
The safety of MoS routes and quality assurance of the main fairway system is the most
important problem. Traditional risk analysis approach calculates the possibility of ship collision
with historical data, mathematical models and opinions of experts, which evaluates the present risk
level. However, there are not any historical data of MoS crossing that can be used for the analyses
to describe the future of the studied area.
Figure 1. Proposals from the joint Baltic Call for MoS, www.pisil.pl
58
RT&A # 02 (25)
(Vol.1) 2012, June
So it is necessary to do a risk analysis based on a mathematical model, with a combination

of forecasting and simulating system, which should be verified by a concrete example.
2
STOCHASTIC MODEL
The model can be used to protect safety of navigation in the congested areas. This model
presents an approach for modeling both spatial interactions and detailed succession dynamics in the
MoS crossing by placing the semi-Markov processes within a class of stochastic process called
piecewise deterministic Markov (PDM) processes.
A piecewise-deterministic Markov process is a stochastic process that evolves
deterministically until a random time when the process jumps to a new (random) state. PDM
processes were introduced by Davis (1984). These processes are readily amenable to computer
simulations. PDM processes have been used in a variety of settings, including storage processes,
capacity expansion problems, and financial investment models.
To define a PDM process we have to define four basic components such as:
the state space,
a description of the deterministic motion between jumps,
the rate at which jumps occur,
the distribution of the state after a jump has occurred.
Each of these components is described below within the context of ships dynamics. The
deterministic motion of the PDM processes given here is only used to keep track of the times that
ships have been in their current states. In fact, the processes used here are piece-wise linear. A MoS
intersection is divided into several cells (Fig. 2). The area within a cell is assumed to be
environmentally homogeneous belongs to only one MoS or crossing place.
Figure 2. An example of a MoS intersections grid of three types of cells
The state of a MoS intersection is determined by the distribution of ships within it and the
respective functional roles of those ships. This paper focuses on modeling the dynamics of cells.
Assume that a cell may be in any of 5 states and denote the set of possible states by S. The Table 1
specifies five status types or states defined by a ship appearance.
59
RT&A # 02 (25)
(Vol.1) 2012, June
Table 1. Status types and state numbers

State Number
1
2
3
4
5
Status Type
Gap
A ship entering a cell
A ship in a cell
A ship living a cell
Collision alert
The state space of the MoS crossing is defined to be, (Monticino et al. 2002)
E ( S [0, )) N .
(1)
The state of the MoS crossing at time t is

X t (( s1 (t ), 1 (t )), , ( s N (t ), N (t )))
where N = number of cells; si(t) = is the state of cell i at time t; i(t) = is the time that the plot has
been in state si(t) since the last time it changed states.
The deterministic portion of the PDM processes is used here only to keep track of how long
cells have been in their current states. Thus, between jump times, the state of the each cell remains
constant, while the time in that state evolves at unit rate.
3
COLLISION FREQUENCY MODELS
There is one major difference between two collision types (Fig. 3). At X-shaped
intersection the traces of two ships always intersect, whereas they on average only intersect in one
out of two cases at Y-shaped intersection. This means that the geometrical collision probability
needs to be corrected by a factor of 0.5 in case of a Y-shaped intersection.
Figure 3. The two collision type, www.sofartsstyrelsen.dk/SiteCollectionDocuments/CMR/

Sejladssikkerhed,%20GMDSS%20og%20SAR/Bassy/BASSY%20Evaluation.pdf
The Pedersens model considers the crossing of two waterways, (Pedersen 1995). Ships are
grouped by their type and length in order to utilize the different characteristics of vessel groups like
the average speed or manoeuvrability which varies significantly from one ship group to another.
60
RT&A # 02 (25)
(Vol.1) 2012, June
Model of Fowler and Srgrd suggest that the frequency of critical situations is calculated
assuming that traffic movements are uncorrelated, (Jutta 2010). A critical situation denotes that two
ships are crossing within half a nautical mile from each other. Encounter frequency is estimated by
a pair-wise summation across all shipping lanes at the considered location. They do not present a
practical procedure to calculate the number of critical situations.
Macduffs Model is build on molecular collision theory, (MacDuff 1974). Ships on a
shipping lane are regarded as a homogenous group: they are navigating at the same speed and they
have similar dimensions.
Cowi Crossing Collision Model defines crossing collision as a collision that includes ships
sailing along different waterways. Two ships can theoretically collide if their traces intersect. The
possibility of a collision between two ships navigating at intersecting routes can be expressed by
critical time interval.
4
SIMULATION APPROACH
The necessary time to make a decision by navigator usually amounts 3 to 6 min and this
decision time can be also considered in the collision checking range (Xue et al. 2009). However as
we consider only the crossing situation, in which may be involved more than two ships in close
proximity, the average checking range on the waterway should be smaller. For example, nowadays,
for security reasons it is recommended that Automatic Identification System AIS, thereby
controlling system, allows for transmitting minimum 2000 messages per minute.
First we have to find the cell size, as we will determine the simulation time step with the use
of it. We assume, similarly as in the Nagel-Schreckenberg model presented in (Nagel &
Schreckenberg 1992, Wahle et al. 2001), that the waterway is divided into cells with a length of
1
CS cross
. We determine the total crossing density cross according to the following formula:
cross max{12 , 13 }, where 12 1 2 2 2 , 13 1 2 3 2
LOA3
LOA1
LOA2
, 2
, 3
,
and 1
T1 V1 LOA1
T2 V2 LOA2
T3 V3 LOA3
(2)
LOA1 = a length of a ship on the main waterway; LOA2, LOA3 = a length of ships on the
lateral waterways 2 and 3; T1, T2, T3 = mean times between ships starting for waterways no 1, 2
and 3.
The article deals with different levels of risk depending on the mutual distance of vessels
that are on collision courses. Distance is measured in taxicab metric, according to the sequence of
grid cells. We denote a size of a cell by CS. To determine the threshold values of safety distances
defining the states of collision risk we have to consider following results. A ship domain is the area
around the vessel that should be avoided by other vessels and an overlap of two vessels domains is
concerned with a very high risk of collision. A length of a ship domain is assumed to be equal
4LOA (Fujii & Tanaka 1971) and we consider this values as a critical distance of high risk of
collision to define the first risk level. To define the second risk level we assume from references a
distance of passing clear CPA (Closest Point of Approach) equal to 1 nm (1852m). This distance of
passing clear will be considered in the simulation as a critical value between high and low risk of
collision. Further we take a distance of 2 nm as a safety distance corresponding to the low risk of
collision.
Then we define following risk levels:
high risk of collision both ships are entering a cell and number blank cells between two
ships is equal at least d1-2, one ship is entering a cell and second ship is in a cell or leaving a cell
and number blank cells between two ships is equal at least d1-1, both ships are in a cell or leaving a
61
RT&A # 02 (25)
(Vol.1) 2012, June
cell and number blank cells between two ships is equal at least d1; if none of these conditions if
fulfilled we define this situation as a collision alert; the distance d1 is determined from the equation:
4 LOA1 4 LOAmax
d1
CS
(3)
where d denotes the integer part of number d plus 1; LOA1 = a length of a ship on the main
waterway; LOA2, LOA3 = a length of ships on the lateral waterways 2 and 3; LOAmax = max{LOA2,
LOA3}; CS = size of a cell;
low risk of collision both ships are entering a cell and number blank cells between two
ships is equal at least d2-2, one ship is entering a cell and second ship is in a cell or leaving a cell
and number blank cells between two ships is equal at least d2-1, both ships are in a cell or leaving a
cell and number blank cells between two ships is equal at least d2, were distance d2 is determined
from the equation:
1852
d2
;
CS
(4)
negligible risk of collision both ships are entering a cell and number blank cells between
two ships is equal at least d3-2, one ship is entering a cell and second ship is in a cell or leaving a
cell and number blank cells between two ships is equal at least d3-1, both ships are in a cell or
leaving a cell and number blank cells between two ships is equal at least d3, were distance d3 is
determined from the equation:
2 1852
d3
.
CS
(5)
We assume that the simulation time step corresponds to time of vessel moving with the
largest velocity from one cell to the next. Thus the simulation step time, denoted by t , is
determined from the formula:
t
CS
, where V max max{V1 , V2 , V3 }
V max
(6)
and size of a cell CS was described with use of crossing density given in (2). Velocities of vessels
on each waterway V1 , V2 , V3 in knots in the program corresponds to the speed measured in cells per
time step.
4.1
Simulation language and environment
The computer program is written in Java language using SSJ V2.1.3 library with support of
stochastic simulations. The documentation of SSJ can be found in Simard (2012). The Java platform
is the object-oriented programming language that provide several standard packages. To perform
the simulation we use a javaSimulation package that is devoted to process-based discrete event
simulation. The package is a Java implementation of the simulation facilities provided by the
programming language SIMULA. The javaSimulation package provides three different approaches
to discrete event simulation: event-based, activity-based and process-based. The description in
details along with appendices containing Java source code and documentation are given in the
62
RT&A # 02 (25)
(Vol.1) 2012, June
report (McNab 1996).

Java-based simulation tools are very popular because it is the only object-oriented
programming environment that effectively supports standardized components (Kilgore et al. 1998).
For example the book (Garrido 2001) concentrates on object-oriented modeling of simulation using
Java and practical simulation techniques.
4.2
Simulation model
Discrete-event system is a system completely determined by random event times and by the
changes in state taking place at these moments. Basic approaches for constructing a discrete-event
simulation model are event scheduling, activity scanning and process interaction approach. The
event scheduling approach focuses on event, i.e. the moment in time when state changes occur,
while process interaction focuses on processes, i.e. the flow of each entity through the system. In
the activity scanning approach in each cycle of simulation there are independently checked
conditions of all events occurring. For a comprehensive description of basic methods and
techniques related to computer simulation of discrete event systems, a reader is referred to Tyszner
(1990). In the paper we will concentrate on simulation of discrete-event systems by event
scheduling approach.
The event oriented simulation concentrates on handling and sending events. The activity
following each event is implemented as an event routine and the event routine may schedule new
events and re-schedule existing event. In this approach we have to define states, events, rules telling
what will happen when an event occurs and some parameters.
In the simulation model there are considered three sea motorways with four points of
collision marked at the scheme (Fig. 4). Collision points number 1 and 3 are X-shaped intersection
while collision points number 2 and 4 are the type of Y-shaped intersection. These two types of
collision are described in Section 3. The main flow is on the sea motorway 1. The collision problem
includes the characteristics of the ships and their motion before, during and after collision.
Figure 4. The scheme of the collision problem
4.3
Computer program
We consider following major ship types: a tanker, a container carrier, a passenger ship, a
RoPax, a general cargo ship and fast ferry. The mean velocity and length for different types of
63
RT&A # 02 (25)
(Vol.1) 2012, June
vessels at each waterway is given by default from empirical data for maritime traffic in Gulf of
Finland presented in (Montewka et a. 2010). The user can change manually, given by default while
computer starting, values of the vessels velocity and length (Fig. 5).
Figure 5. The starting window of the simulation program

In the simulation the ships starting times can follow different distribution. We can choose
these distribution on the main and lateral waterways from the following list: deterministic, uniform,
exponential, Erlang, normal, log-normal, Beta, gamma and triangular. Thus using the program there
is possible to simulate also non Poisson streams. In Fig. 6 there is presented a program window that
is showing the moment of reading data and choosing the distributions of ships starting times. The
mean waiting time between ships departing for each sea motorway is given by the user. Given by
default the percentage of the starboard ships at each waterway can be modified by the user. In the
simulation the user can also change the distances to collision points. In the case the flow of vessels
into waterway is uniform there is also necessary the minimal accepted distance between vessels,
that is given by default in the computer program depending on the vessels length. It is assumed by
default that the intervals times between successive notifications at the sea motorways 2 and 3 have
exponential distribution and at the sea motorway 1 an uniform distribution.
Figure 6. The program window for reading parameters concerned with vessels move
64
RT&A # 02 (25)
(Vol.1) 2012, June
The scheme of the computer program used for simulation is presented in Fig. 7 and the rules
applied in the computer simulation can be described as follows:
START
DATA
BASE
initial data
run
simulation
SIMULATION MODEL
object ship
object ship
object ship
event list for

waterway 1
starboard
according to the selected

event
scheduler distributions
event list for
waterway 2
headway
headway
starboard
event list for

waterway 3
starboard
headway
control
command
event list to
collision
point 2
main waterway
lateral waterways
event list to
collision
point 1
event list to
collision
point 3
RUN
time
scaning
event list to
collision
point 4
execute and
re-schedule current
events list
check risk of
ships collision
end
simulation
statistical
results
END
Figure 7. The scheme of the simulation program

the system saves data and gives results from the point of view vessels being on the main
waterway 1; we assume that vessels on the main waterway 1 are stand-on vessels, while vessels on
the waterways 2 and 3 are give-way vessels,
at each time step the system examines the ship closest to the crossing on each waterway
and if the safety state of a ship on the waterway 1 has changed the system saves this transition and
the time during which the ship was in the previous safety state,
65
RT&A # 02 (25)
(Vol.1) 2012, June
if the ship on the main waterway 1 is not fully safe the system also examines the next ships
on this waterway,
in the considered crossing situation the starboard ships on the waterway 3 are safe as they
do not cross any other waterway,
in the vessel on the waterway 3 follows headway this waterway then the system checks
situation at the collision point 1 and 2 and the safety state of the vessel, if there is on the main
waterway 1; after passing the collision point 1 if there is a ship on the waterway 1 going starboard
the system checks the situation at the collision point 2, if there is no starboard ship on the waterway
1 the vessel on the waterway 3 is safe and it is not controlled by the system,
the starboard ships on waterway 1 are examined at the collision point 2, while the headway
ships are first checking at the collision point 1, next at the collision point 3 and 4; after passing the
collision point 4, there is assumed the negligible risk of collision,
the starboard ships on the waterway 2 are examined at the collision point 4, while the
headway ships are examined at the collision point 3,
the process is repeated throughout the entire simulation time.
In the computer program, according to the accepted before risk levels and their critical
distances defined by (3)-(5), there are assumed following states:
state 0 collision alert;
state 1 high risk of collision;
state 2 low risk of collision;
state 3 negligible risk of collision.
We denote by p(i), i = 0,1,2,3, probability of system being in the safety state i.
4.4
Results of the program
As a result of the program we obtain following data: the matrix of the system transitions
number between the states and the realizations of the conditional sojourn times at the state until the
transition to the other state. From these results there are also determined the matrix of probabilities
of the systems transitions between the states and the vector of probabilities of the systems being in
the particular states during the simulation time. The obtained from the simulation results can be
used for further identification and safety analysis of the system (Blokus-Roszkowska et al. 2011).
The proposed simulation model is sensitive to its changing parameters, that is depicted at the
graph (Fig. 8). For this reason the input data of the program must be properly selected. The obtain
from the simulation results show the dependency of probabilities of system being in the safety states
during the simulation time on the mean time between vessels departure i.e. on the intensity, at the
sea motorway 1. In the presented example we assume that the departure time of a ship on the main
waterway follows an uniform distribution and the interval times between ships departure on the
crossing waterways are log-normal. We assume that the mean time between ships departure on the
sea motorway 2 equals 1.2 h with standard deviation 0.9 h and on the sea motorway 3 equals 1.3 h
with standard deviation 1 h, equivalently.
66
RT&A # 02 (25)
(Vol.1) 2012, June
0,7
0,6
p(0)
0,4
p(2)
p(3)
0,3
probability
0,5
p(1)
2,8
mean time between vessels' departure

at sea motorway 1
0
2,5
2,2
0,1
1,9
1,6
1,3
0,2
p(0)
p(1)
p(2)
p(3)
Figure 8. Probability of being in the safety states depending on the mean time between vessels
departure at the sea motorway 1
CONCLUSIONS
Java can expose the benefit of computer simulation to a larger audience of problem-solvers,
decision-makers and trainers. Java-based simulation components could be easily distributed,
executed and modified throughout the word over the internet. The main advantage of Java is crossplatform compatibility that can eliminate the need to maintain different versions of the software.
Considering this Java-based programs can provide solutions or assistance in safety transportation
system planning and support cooperation or exchange of information between ships, ports and
terminals. The presented program can serve a base to create a new service that supports the
development of the concept of sea motorways. The simulation programs can optimize the logistical
transportation system with integration of sea motorways.
The idea of the paper was to develop the simulation environment to test the features of
computer-controlled sea motorways. The proposed simulation model for safety state evaluation can
be helpful on the level of safety management. The model allows to predict safety state depending on
changing traffic strength at sea motorways. With the concept of safety management such simulation
models could face the problem of sea transportation development and provide new solutions for
operational optimization and safety transportation system planning.
REFERENCES
Blokus-Roszkowska, A., Koowrocki, K. & Fu X. 2011. Integrated software tools
supporting decision making on operation, identification, prediction and optimization of complex
technical systems reliability and safety, Part 2. Summer Safety & Reliability Seminars Journal of
Polish Safety and Reliability Association (5) 2: 279-288.
Davis, H.H.A. 1984. Piecewise deterministic Markov processes: a general class of nondiffusion stochastic models. Journal of the Royal Statistical Society Series B, 46(3): 353-388.
Fujii, Y., Tanaka, K. 1971. Traffic capacity. Journal of Navigation 24: 543-552.
67
RT&A # 02 (25)
(Vol.1) 2012, June
Garrido, J.M. 2001. Object-Oriented Discrete-Event Simulation with Java. New York:
Kluwer Academic/Plenum Publisher.
Jutta Y. 2010. Modelling Marine Accident Frequency. Masters thesis submitted in partial
fulfillment of the requirements for the degree of Master of Science in Technology in the Degree
Programme in Engineering Physics.
Kilgore, R., Healy, K. & Kleindorfer, G. 1998. The future of Java-based simulation.
Proceedings of the 1998 Winter Simulation Conference, IEEE Computer Society Press: 1707-1712.
MacDuff, T. 1974. The probability of vessel collisions. Ocean Industry: 144-148.
McNab, R. 1996. SimJava: a discrete event simulation library for Java. University of
Newcastle upon Tyne. Available from http://www.javaSim.ncl.ac.uk
Montewka, J., Goerlandt, F., Kujala, P. 2012. Determination of collision criteria and
causation factors appropriate to a model for estimating the probability of maritime accidents. Ocean
Engineering 40: 50-61.
Montewka, J., Hinz T., Kujala, P., Matusiak, J. 2010. Probability modelling of vessel
collisions. Reliability Engineering and System Safety 95: 573-589.
Monticino, M.G., Cogdill, T. & Acevedo, M.F. 2002. Cell Interaction in Semi-Markov
Forest Landscape Models, Integrated Assessment and Decision Support, Proceedings of the First
Biennial Meeting of the International Environmental Modelling and Software Society iEMSs 1:
334- 339.
Nagel, K., Schreckenberg, M. 1992. A cellular automaton model for freeway traffic. Journal
de Physique I France 2(12): 2221-2229.
Paulauskas, V. & Bentzen, K. 2007. Scenario Analysis of New Maritime Corridors, BaSIMBaltic Sea Information Motorways, Aalborg 27.
Pedersen, P. T. 1995. Collision and Grounding Mechanics. Proceedings WEMT 1995, 1:
125-157.
Pietrzykowski, Z. & Uriasz, J. 2009. The Ship Domain A Criterion of Navigational Safety
Assessment in an Open Sea Area. Journal of Navigation, 62(01): 93108.
Simard R. 2012. SSJ: Stochastic Simulation in Java, Document available online at:
http://www.iro.umontreal.ca/~simardr/ssj/indexe.html.
Smalko Z. & Smolarek L. 2010. Modelling a ship safety according to collision threat for
ship routes crossing. Scientific Journals Maritime University of Szczecin, 20(92): 120127.
Tyszner, J. 1990. Digital simulation (in Polish). Wydawnictwo Naukowo-Techniczne,
Warsaw.
Wahle, J., Neubert, L., Esser, J., Schreckenberg, M. 2001. A cellular automaton traffic flow
model for online simulation of traffic. Parallel Computing 27: 719-735.
Xue, Y., Lee, B.S., Han, D. 2009. Automatic collision avoidance of ships. Proceedings of
the Institution of Mechanical Engineers, Part M: Journal of Engineering for the Maritime
Environment 233(1): 33-46.
68
Kumar Pardeep and Chaturvedi D.K., Pahuja G.L. - A HEURISTIC METHOD FOR RELIABILITY REDUNDANCY OPTIMIZATION OF FLOW NETWORKS
RT&A # 02 (25)
(Vol.1) 2012, June
A HEURISTIC METHOD FOR RELIABILITY REDUNDANCY

OPTIMIZATION OF FLOW NETWORKS
Kumar Pardeep * and Chaturvedi D.K. **
Department of Instrumentation,
Kurukshetra University, Kurukshetra, 136 119, India
e-mail: pardeep_kuk@rediffmail.com*, dhiraj_chaturvedi@yahoo.co.in**
*Corresponding author
Pahuja G.L.
Department of Electrical Engineering,

National Institute of Technology, Kurukshetra 136 119, India
e-mail: pahuja_gl@yahoo.co.in
ABSTRACT
In flow networks, from the quality and service management point of view, measurement of
the transmission ability of a network to meet the customers demand is very important. To meet the
ultra high reliable requirements of such networks, a heuristic method for reliability redundancy
optimization of flow networks using composite performance measure (CPM) integrating reliability
and capacity has been proposed. The method is based upon the selection of main flow paths and
backup paths and then optimizing main paths on priority basis. Thus, the reduced computation work
makes the proposed algorithm suitable for designing of large, reliable telecommunications
networks.
Keywords: flow networks;
optimization; heuristic algorithm.
1
capacity
related
reliability;
constrained
redundancy
INTRODUCTION
Constrained reliability redundancy optimization of networks has generally been studied with
reliability as connectivity measure. The practical systems such as computer networks,
telecommunication networks, transportation systems, electrical power transmission networks,
internet etc. can only transport limited amount of flow therefore these are termed as flow limited
networks. To meet the ultra high reliability requirements of such networks, a heuristic method for
reliability redundancy optimization of flow networks has been proposed.
Literature is enriched with reliability redundancy optimization of networks with connectivity
only as a measure of performance (Sharma & Venkateswaran 1971, Aggarwal 1976, Golden &
Magnanti 1977, Gopal et al. 1978, Lee 1980, Bodin et al. 1982, Xue 1985, Dinghua 1987, Fredman
& Tarjan 1987, Kim & Yum 1993, Shen 1995, Martins & Santos 1997, Schrijver 1998, Ahuja 1998,
Kuo & Prasad 2000, Kuo et al. 2001, Park et al. 2004, Pascoal et al. 2005, Kumar et al. 2009,
2010a, b, 2011). However, reliability redundancy optimization of flow networks has rarely been
studied. In present days context existing methods do not fulfil the requirements of management of
quality of service. The reliability redundancy optimization techniques discussed in Gopal et al.
(1978), Dinghua (1987), Kim & Yum (1993), Shen (1995), Park et al. (2004), Kumar et al. (2009,
2010a, b, 2011) are not suited to flow networks. This paper presents a technique for reliability
redundancy optimization of flow networks using combined performance measure named capacity
related reliability (CRR).
69
RT&A # 02 (25)
(Vol.1) 2012, June
A path is a sequence of arcs and nodes connecting a source to a sink. All the arcs and nodes
of network have its own attributes like delay, reliability and capacity etc.. From the quality and
service management point of view, measurement of the transmission ability of a network to meet
the customers demand is very important (Lin 2006). When a given amount of flow is required to be
transmitted through a flow network, it is desirable to optimize the network reliability to carry the
desired flow. In such cases, the system reliability is the measure of quality of the system capability
to transmit desired flow. The capacity of each arc (the maximum flow passing the arc per unit time)
has two levels, 0 and/or a positive integer. The system reliability is the probability that the
maximum flow through the network between the source and the sink is not less than the demand
(Golden & Magnanti 1977, Lee 1980, Bodin et al. 1982, Fredman & Tarjan 1987, Ahuja 1998, Lin
2003, 2004, Pahuja 2004, Lin 2006, 2007a, b). For determining the reliability it is generally
assumed that network is capable of transmitting any required amount of flow between source (s)
and terminal (t) nodes of the network. This presumption is neither valid nor justifiable for real life
systems as links and nodes can carry only limited amount of flow. In 1980 Lee did the pioneer
work of integrated both capacity and reliability & named it combined performance measure as
capacity related reliability (CRR) and also termed such networks as flow networks. Max-FlowMin-Cut theorem has been used to determine the capacity of the network (Sharma &
Venkateswaran 1971, Xue 1985, Shen 1995, Martins & Santos 1997, Schrijver 1998, Lin 2003,
2004, Park et al. 2004, Lin 2006, 2007a, b).
2
COMPOSITE PERFORMANCE MEASURE
Reliability under flow constraint is a more realistic performance measure for flow networks.
The concept of weighted reliability introduced by Aggarwal (1976) requires that all the successful
states qualifying connectivity measure of the network be enumerated. The probability of each
success state is evaluated and is multiplied by the normalized weight (Aggarwal 1985). Rushdi
(1988) evaluated the same performance index as evaluated by Aggarwal 1985 using decomposition
approach. Methods given by Aggarwal (1985), Rushdi (1988), and Shakti (1995) generate both
cancelling (failed) and non-cancelling (success) terms.
The following section presents a heuristic algorithm for reliability redundancy optimization
of flow networks using composite performance measure (CPM) and the method has been utilized to
determine the capacity related reliability performance index.
2.1
Notation
al(X) Sensitivity factor of lth minimal path set

bi(xi) Subsystem selection factor for ith subsystem with xi components
Cj Total amount of resource j available
g ij ( xi ) /
Amount of resources-j consumed in subsystem-i with xi components / Cost of
c ji ( xi )
subsystem i for jth constraint.
h(.) Function yielding system reliability; dependent on number of subsystems (n) and
configuration of subsystems
k Number of constraints, j = 1, 2,, k
L(x) ( L x1 , L x2 ,..., L xn ) , Lower limit of each subsystem i,
m Number of main minimal path sets, l = 1, 2,., m

n Number of subsystems, i = 1, 2,, n
Pl lth minimal path set of the system
1 2
min
1
min
Ps (l ,l ,,l ): priority vector s.t. l and l are the number of minimal path sets arranged
70
RT&A # 02 (25)
(Vol.1) 2012, June
in decreasing order of path selection parameter al(X).

Qi (xi ) Unreliability of subsystem i with xi components.
ri Reliability of a component at subsystem i.
Ri (xi ) Reliability of subsystem i with xi components.
Rr Residual resources [total resource available (Cj) - resources consumed (gi j xi)]
Rs ( X ) System reliability
S ( x ) Set of variables that have been used as key-elements in a given decomposed
expressions
U (x) (U x1 , U x2 ,..., U xn ) , Upper limit of each of subsystem i,
x* Optimal solution
xi Number of components in subsystem i; i = 1,2,.n
X A vector (x1,xn)
Ri Increment in ith stage reliability when a unit is added in parallel to the ith stage
2.2
Assumptions
Following are the assumptions for the rest of the sections:

1.
The system and all its subsystems are coherent.
2.
Subsystem structures (other than coherence) are not restricted.
3.
The networks are modelled with the help of graphs, the paths (ordered pair of arcs
and the members of the ordered pair are reliability and capacity respectively) where in are assigned
as the weight of each link.
4.
Each link can have only two stages up and down.
5.
The network nodes are perfect. If the nodes are not perfect, the method needs to be
modified to deal with nodes failures.
6.
All component states are mutually and statistically independent.
7.
All constraints are separable and additive among components.
8.
Each constraint is an increasing function of xi for each subsystem.
9.
Redundant components cannot cross subsystem boundaries.
2.3
Composite Performance Measure (CPM)
The weighted reliability measure i.e. composite performance measure (CPM), integrating
both capacity and reliability may be stated as by [27, 28]:
CPM
ti Ri
(1)
iS ( x )
Where ti is the normalized weight and is defined as:
ti = Capi / Capmax
i.e. the ratio of capacity in the ith state to the maximum capacity (Capmax) of the system and Ri
probability of the system being in state Si and is computed as:
Ri Pr Si
j / S ij 1
71
pj
qk
k / S ik 0
(2)
2.4
RT&A # 02 (25)
(Vol.1) 2012, June
Capacity Functions of Networks
The capacity function of different arcs connected in parallel is [25]:

C X Par Capi
(3)
i x
and the capacity function of different arcs connected in series is:

C X Ser min Capi
(4)
The rules for connecting series and parallel arcs to integrate capacity and reliability to give
composite performance measure are expressed as:
n
(5)
CR X min Cap r
i
Ser
CR X
Par
ix
n
i
i 1
n
Cap i . ri
i 1
(6)
i 1
CPM for series and parallel networks can be defined as:

CPMPar = CR(X)Par / Capmax
and
CPMSer = CR(X)Ser / Capmax
3
PROBLEM FORMULATION AND HEURISTIC METHOD
3.1
Problem Formulation
(7)
(8)
The general constrained redundancy optimization problem in complex systems can be

reduced to the following integer programming problem (Kuo et al. 2001):
Maximise
(9)
Rs ( X ) h ( R1 ( x1 ), ..., Rn ( x n )),
n
(10)
subject to g i j ( x i ) C j ,
j = 1, 2, , k
i 1
and
3.2
1 xi U x i ,
i = 1, 2, , n.
Proposed Heuristic Method
In real life systems all the arcs are not simultaneously connected to carry flow from source
to sink. Hence a flow path set is the arcs and nodes that actually carry traffic. In these practical
systems all the path sets are not utilized for transfer of information (Hayashi & Abe 2008). The
flow is transmitted through the main path(s) and in case of failure of this path(s), a backup path
completes the task of main path. The backup path(s) come in operation only when the main paths
fail thus, enhancing the reliability of the network, it is presumed that the main path(s) and backup
path(s) for the network are known. The proposed algorithm first optimizes the main path(s) and then
back up path(s) using redundancy optimization technique. This leads to more efficient use of
resources which are generally limited. Unlike existing heuristics, a stopping criterion has been
applied to switch from reliability redundancy optimization of main path(s) to reliability redundancy
optimization of backup paths. The algorithm considers sensitivity factor as the criteria for selecting
the main path from the list of given main paths and then a subsystem for applying redundancy
within the chosen flow path.
3.3
Steps of the Proposed Method
Step1: Firstly select the main path sets m and back up flow paths of the flow Network.
72
RT&A # 02 (25)
(Vol.1) 2012, June
Step2: Let xi = 1 for all i; i = 1, 2,, n.

Step3: Calculate al ( X ) , l = 1, 2,, m for each minimal path set and find l* such that
al* ( X ) = max [ al (X ) ]
and
Ri ( xi )
i Pl
al ( X )
,
j
( g i ( xi ) / k C j )
i Pl
j 1
l 1,2,3.........m .
Step4: For the chosen minimal path set l* find
bi ( xi )
Ri
k
( g i ( xi ) / k C j )
i*
such that bi* ( xi* ) = max [ bi ( xi ) ] and
, for each i P
l
j 1
where Ri 1 Q i ( xi ) Ri ( xi )
*
Step5: Check, by adding one redundant subsystem to unsaturated subsystem i :

i)
if no constraints are violated, add one redundant subsystem to unsaturated subsystem
i*
by replacing xi* with xi* 1 , and go to step 3.

ii) if at least one constraint is exactly satisfied and others are not violated, then add one
*
redundant subsystem to unsaturated subsystem i by replacing xi* with xi* 1 . The x* = X is the
optimal solution. Go to step 6.
iii) if at least one constraint is violated, then remove minimal path set l* from further
consideration and consider the next path having maximum
al* ( X ) value and go to step
4.
iv) if all minimal path sets are now excluded from further consideration, then x = X is the
optimal solution; else go to step 3
Step6: In case any resources are still available optimize the backup paths as discussed in
Step2.
Step4: Evaluate the composite performance measure (CPM) for each subsystem of the
network.
Step5: Evaluate the system reliability using the CPM of the each subsystem.
4
COMPUTATION AND RESULTS
To illustrate the performance of the proposed algorithm a network having six arcs {x1, x2, x3,
x4, x5, x6} and five minimal path sets {y1, y2, y3, y4, y5} as shown in the Figure 1 is considered and
solved for capacity related redundancy reliability optimization using CPM (7 and 8). System
reliability is determined using Bayes method. The network shown in Figure 1 is a bench mark
problem, considered by Hayashi & Abe (2008).
Using Bayes method, the Reliability of the above system can be expressed as:
73
Rs(X) = R3 [1- Q6 {1-(1-Q1Q2)(1-Q4Q5)}]

+ Q3[1-(1-R2R5)(1-R1R4)]*Q6
(11)
y2
y1
y5
x2
y4
RT&A # 02 (25)
(Vol.1) 2012, June
x3
x5
x6
x1
x4
y3
Figure 1 Illustration Network
The simple minimal path sets of the Network are
P1 = {1, 3, 5}, P2 = {2, 3, 4}, P3 = {1, 4}, P4 = {6}, P5 = {2, 5}
The problem is solved for data given in Table 1. Using this initial data the general problem
of constrained reliability redundancy allocation has been solved using the steps discussed in Section
3.3 above. The problem is solved by considering the flow path sets P1, P2, P3 and P4 as main path
sets and P5 as backup path. The proposed algorithm gives the optimal solution (3, 1, 1, 2, 1, 4) with
system reliability Rs = 0.999998458, the optimized subsystem reliability probability Ri and
unreliability probabilities Qi are shown in Table 2.
Table 1 Data for Fig. 1
i
ri
gi1 / c1i
C1
1
2
3
4
5
6
0.70 0.75 0.8 0.85 0.7 0.9
2
3
2
3
1
3
30
Table 2 Optimized subsystem reliability/unreliability for Fig. 1

i
X*
Ri
Qi
x1
3
0.973
0.027
x2
1
0.75
0.25
x3
1
0.8
0.2
x4
2
0.9775
0.0225
x5
1
0.7
0.3
x6
4
0.9999
0.0001
The capacity of each subsystem of the flow path is taken as 100 and the capacity of flow
paths of the network is determined using proposed approach as:
74
P1 =min {3*100, 100, 100} =100

P2 =min {100, 100, 100} =100
P3 =min {3*100, 2*100} =200
P4 =min {4*100} =400
RT&A # 02 (25)
(Vol.1) 2012, June
(12)
Next the composite performance measure CPM expression (13) is derived using (7 and 8)
and the value for CPM for an assumed flow of 200 is suppose to pass through the flow path and it
comes out to be 1.0000 .
CPM
P1
CPM
P2
CPM
P3
CPM
P4
so
min Cap i
[ R1 R 3 R 5 ]
Cap max
= (100/200) x 0.973 x 0.8 x 0.7 = 0.27244
min Cap i
[ R 2 R3 R 4 ]
Cap max
= (100/200) x 0.75 x 0.8 x 0.9775 = 0.29325
min Cap i
[ R1 R 4 ]
Cap max
= (200/200) x 0.973 x 0.9775 = 0.9511075
min Cap i
[ R6 ]
Cap max
= (400/200) x 0.9999
= (2) x 0.9999 as 0 (min Capi / Capmax ) 1
= 1x 0.9999 = 0.9999
(13)
(14)
(15)
(16)
Composite performance measure integrating the reliability with capacity is calculated as:
CPMNetwork
= 1- (1-CPMP1) (1-CPMP2) (1-CPMP3) (1-CPMP4)

= 1- (1-0.27244) x (1-0.29325) x (1- 0.9511) x (1- 0.9999)
= 1.0000
The above result shows that proposed method is capable of optimizing the flow network to
transport the desired capacity through the network with highest reliability. However, the selection
of main paths and backup paths will affect the quality of composite performance measure. Hence
the proper choice of these paths may be done using cardinality criteria (Kumar et al. 2010b) or any
other hierarchical measures of importance.
5
CONCLUSIONS
This paper has presented a new model for designing reliable flow networks capable of
transmitting required flow. The proposed algorithm utilizes the concept of main and backup flow
paths. The choice of backup and flow paths is application specific and paths with minimum
cardinality may be selected as main path and disjoint paths can be the backup paths. The numerical
example demonstrates that the proposed algorithm is fast for designing large, reliable
telecommunications networks because the task of optimization is reduced, as only few paths are
selected as main paths.
75
RT&A # 02 (25)
(Vol.1) 2012, June
REFERENCES
Sharma J. and Venkateswaran K. V. 1971, A direct method for maximizing the system
reliability, IEEE Trans. Reliability, vol. R-20, no. 1, pp. 256259.
Aggarwal K. K. 1976, Redundancy optimization in general system IEEE Trans.
Reliability, vol R-25, pp. 330-332.
Golden B. L. and Magnanti T. L. 1977, Deterministic network optimization: A
bibliography, Networks, vol. 7, no. 2, pp. 149183.
Gopal K., Aggarwal K. K., and Gupta J. S. 1978, An improved algorithm for reliability
optimization, IEEE Trans. Reliability, vol. R-27, no. 5, pp.325328.
Lee S. H. 1980, Reliability Evaluation of Flow Networks, IEEE Tr. Rel. Vol. R-29, pp.
24-26.
Bodin L. D., Golden B. L., Assad A. A., and Ball M. O. 1982, Routing and scheduling of
vehicles and crews: the state of the art, Computers and Operations Research, vol. 10, no. 2, pp.
63211.
Aggarwal K. K. 1985, Integration of Reliability and Capacity in Performance Measure of
Telecommunication Network, IEEE Trans. Reliability, vol. R-34, pp. 184-186.
Xue J. 1985, On multistate system analysis, IEEE Trans. Reliability, vol. 34, no. 4, pp.
329337.
Dinghua S. 1987, A new heuristic algorithm for constrained redundancy-optimization in
complex systems, IEEE Trans. Reliability, vol. R-36, no.5, pp. 621623.
Fredman M. L. and Tarjan R. E. 1987, Fibonacci heaps and their uses in improved
network optimization algorithms, J. ACM, vol. 34, no. 3, pp. 596615.
Rushdi M. 1988, Performance Indexes of Telecommunication Networks, IEEE Trans.
Reliability, vol. R-37, pp. 57-64.
Rai S., and Soh S. 1991, A computer Approach for Reliability Evaluation of
Telecommunication networks with heterogeneous Link capacities, IEEE trans. On Reliability, Vol.
40, pp 441-451.
Kim J. H. and Yum B. J. 1993, A heuristic method for solving redundancy optimization
problems in complex systems, IEEE Trans. Reliability, vol. 42, no. 4, pp. 572578.
Varshney, Joshi, and Chang A. R. 1994, Reliability modeling and performance evaluation
of variable link capacity networks, IEEE Trans. Reliability, vol. 43(3), pp. 378-382.
Shen Y. 1995, A new simple algorithm for enumerating all minimal paths and cuts of a
graph, Microelectronics and Reliability, vol. 35, no. 6, pp. 973976.
Shakti K. 1995, Some Aspects of Reliability Related Performance Evaluations of
Computer Communication Networks, Ph.D. Thesis, Kurukshetra University, Kurukshetra.
Martins E. D. Q. V. and Santos J. L. E. D. 1997, An algorithm for the quickest path
problem, Operations Research Letters, vol. 20, no. 4, pp. 195198.
Schrijver 1998, Theory of Linear and Integer Programming. NY: John Wiley & sons, pp.
155156.
Ahuja R. K. 1998, Minimum cost-reliability ratio problem, Computers and Operations
Research, vol. 16, no. 1, pp. 8389.
Kuo W. and Prasad V. R. 2000, An Annotated Overview of System-Reliability
Optimization, IEEE Trans. Reliability, vol. 49(2), pp.176-187.
Kuo, W. Prasad, R., Tillman, F.A. and Hwang, C.L. 2001 Optimal Reliability Design
Fundamentals and Applications, Cambridge University Press.
Lin Y. K. 2003, Extend the quickest path problem to the system reliability evaluation for a
stochastic-flow network, Computers and Operations Research, vol. 30, no. 4, pp. 567575.
Park C. K., Lee S., and Park S. 2004, A label-setting algorithm for finding a quickest
path, Computers and Operations Research, vol. 31, no. 14, pp. 24052418.
Lin Y. K. 2004, Reliability of a stochastic-flow network with unreliable branches & nodes
76
RT&A # 02 (25)
(Vol.1) 2012, June
under budget constraints, IEEE Trans. Reliability, vol. 53, no. 3, pp. 381387.
Pahuja G. L. 2004, Reliability Evaluation and Optimization Recent and New
Approaches, Ph.D. thesis, Kurukshetra University, Kurukshetra.
Pascoal M. M. B., Captivo M. E. V., and Clmaco J. C. N. 2005, An algorithm for
ranking quickest simple paths, Computers and Operations Research, vol. 32, no. 3, pp. 509520.
Lin Y. K. 2006, Reliability of a computer network in case capacity weight varying with
arcs, nodes and types of commodity, Reliability Engineering and System Safety, vol. 2(5), pp. 1-7.
Lin Y. K. 2007,
a.
Reliability evaluation for an information network with node failure under cost
constraint, IEEE Trans. Systems, Man and Cybernetics-Part A: Systems and Humans, vol. 37, no.
2, pp. 180188.
b.
On a multi-commodity stochastic-flow network with unreliable nodes subject to
budget constraint, European Journal of Operational Research, vol. 176, no. 1, pp. 347360.
Hayashi M. and Abe T. 2008, Evaluating Reliability of Telecommunications Networks
Using Traffic Path Information, IEEE Trans. Reliability, vol. 57(2), pp. 283-294, 2008. pp. 283294.
Kumar P., Chaturvedi D. K. and Pahuja G. L. 2009, Heuristic algorithm for constrained
redundancy reliability optimization and performance evaluation. Proc. IMechE, Part O: J. Risk and
Reliability, vol. 223 (4), 381-386.
Kumar P., Chaturvedi D. K. and Pahuja G. L. 2010,
a. Heuristic Methods for Solving Redundancy Allocation in Complex Systems, Int. J. of
Reliability and Safety, vol. 4, No. 2/3, pp. 285-298.
b. An efficient heuristic algorithm for determining optimal redundancy allocation of
complex networks, Reliability Theory and Applications, vol. 3(18), pp.15-28.
Kumar P., Chaturvedi D. K. and Pahuja G. L. 2011, Constrained Reliability
Redundancy Optimization of Complex Systems using Genetic Algorithm, MIT Int. J of Electrical
and Instrumentation Engineering, vol. 1, No. 1, pp. 41-48.
77
A.N. Patowary, J. Hazarika AND G. L. Sriwastav. - ESTIMATION OF RELIABILITY IN INTERFERENCE MODELS USING MONTE CARLO SIMULATION
RT&A # 02 (25)
(Vol.1) 2012, June
ESTIMATION OF RELIABILITY IN INTERFERENCE MODELS

USING MONTE CARLO SIMULATION
A.N. Patowary, J. Hazarika and G. L. Sriwastav
Department of Statistics, Dibrugarh University, India

e-mail: a_patowary@yahoo.com, jitenhazarika@rediffmail.com,
wastav2@yahoo.com
Abstract
This paper presents estimation of reliability R P(X Y) of a system, for the cases when
its strength (X) and stress (Y) follow exponential, normal or gamma distributions, using Monte
Carlo simulation (MCS). First the parameters of strength and / stress are estimated and substituting
them in the reliability expressions, in different cases, the estimates of reliability are obtained.
Normal distribution is fitted to various sets of estimated reliability R , generated by MCS. The
goodness of fit is tested using Kolmogorov-Smirnov one sample test.
Keywords: Stress-Strength; Monte-Carlo Simulation; Kolmogorov-Smirnov one sample test.
1. Introduction
In interference theory of reliability, reliability and other reliability characteristics of a system
can be expressed as some functions of the parameters of the distributions of the random variables
(r.v.s), strength (X) and stress (Y) associated with the functioning of the system. We estimate these
parameters and substitute these values in the expressions for reliability and other characteristics to
get their estimates. The estimates of parameters used here are maximum likelihood estimators and
as such from the invariance property of MLEs, the corresponding estimators of reliability are also
MLEs. In absence of hard data the numerical values of the estimators can be obtained from
simulation. There exists extensive literature for estimation of reliability analytically for single
component systems e.g. Mazumder [12], Church and Harris [4] etc. But the reliability expressions
for multi-component systems are not simple enough to facilitate analytical estimation of reliability
and its other characteristics. Also due to lack of hard data, one way out is simulation, in particular
Monte Carlo simulation.
With simulation technique it is possible to estimate reliability (or probability of failure) or
other reliability characteristics without going into the analytical techniques. The availability of
personal computer and software makes the process comparatively simple. In fact, to evaluate the
accuracy of the sophisticated analytical techniques or to verify a new technique, simulation is
routinely used to independently evaluate the underlying probability distributions.
1.1 Monte Carlo Simulation:
The Monte Carlo Simulation method is an artificial sampling method which may be used
for solving complicated problems in analytic formulation and for simulating purely statistical
problems. In the simplest form of simulation, each r.v. in a problem is sampled several times to
78
RT&A # 02 (25)
(Vol.1) 2012, June
represent its real distribution. Each realization of r.v.s in the problem produces a set of numbers
that indicates one realization of the problem itself. Solving the problem deterministically for each
realization is known as a simulation cycle, a trial or a run. Using many simulation cycles we get the
overall probabilistic characteristics of the problem, particularly when the number of cycles N is
sufficiently large. Simulation, using a computer, is an inexpensive way (compared to laboratory
testing) to study the uncertainty in a problem.
The primary components of a Monte Carlo Simulation include the followings:
(i) Probability distribution function (or probability density function): The physical (or
mathematical) system must be described by a set of probability distribution functions.
(ii) Random number generator: A source of random numbers uniformly distributed on the unit
interval must be available.
(iii) Sampling rule: A prescription of sampling from specified distribution function, assuming the
availability of random numbers on the unit interval.
(iv) Scoring (or tallying): The outcome must be accumulated into overall tallies or scores for the
quantities of input.
In this paper except exponential distribution we have not used uniformly distributed random
numbers, rather obtained random numbers following particular distribution directly from
MATLAB.
In Section 2, we have estimated reliability of an n-standby system (n=1, 2, 3), through
Monte Carlo Simulation technique. Simulation is performed for exponential stress-strength, normal
stress-strength and gamma stress-strength. In Section-3 we have considered fitting of normal
distribution to estimated reliability in each case, for different true values of the parameters. The
goodness of fit is tested by K-S one sample test (Seigel [18]). Since we have taken a small sample,
20, only, when using 2 test, the number of classes becomes too few, due to pooling. We have
considered the fitting of normal distribution to check whether normal approximation is good enough
for a small of 20.
Some literatures on the topic which we have come across are:
Kamat and Riley [8] presented MCS for a complex system for time to failure (TTF) models.
Some of the others studies of reliability estimation using MCS for TTF models includes
Pulido et.al. [15], Goel [5], Hong and Lind [6], Landis et.al. [10], Tunak et.al. [23], Naess et.al.
[13], Wu et.al. [24] etc.
Stancampiano [21] applied simulation to interference models. Manders et.al. [11], Aldrisi
[1], Stumpf and Schwartz [22], Zhang et.al. [25] have simulated stress-strength. Paul and
Borhanuddin [14], Rezaei et.al. [17] estimated reliability of stress-strength model, using MCS.
Ahmad et.al.[2] obtain Bayes estimates of P (Y< X) using MCS. Borhanuddin et.al. [3] estimated
reliability for multicomponent system using MCS. Rao et al [14] compared reliability estimates for
multicomponent systems evaluated by different methods such as method of moments, modified ML
method and Best Linear Unbiased Estimator through MCS technique.
Kakati and Sriwastav [7] and Sriwastav [20], used simple simulation by taking random
exponential numbers to represent stress-strength. They considered very small samples. From these
samples they first estimate the parameters and substituting these in the expressions of reliability
they get estimated reliability.
2 Reliability Estimation through Monte Carlo Simulation:
Let us consider an n-standby system. Let X1, X2,,Xn be the strengths of the n components
in the system arranged in the order of activation. Let Y1, Y2,,Yn be the stresses faced,
respectively, by 1st, 2nd,,nth component, when they are activated; Xis and Yis are all
independent. For a detailed description of such a system one may refer (Sriwastav and Kakati, [19]).
The reliability Rn of an n-standby system for a single impact of stress is given by,
Rn = R(1) + R(2) + + R(n),
(2.1)
79
RT&A # 02 (25)
(Vol.1) 2012, June
where R(i) is the increment in the system reliability due to the ith component, defined as
R i P X 1 Y 1 , X 2 Y 2 ,..., X i 1 Y i 1 , X i Y i
(2.2)
Here, we have assumed that all the components are having the same strength distributions
and are working under the same environment (stress), i.e. all Xis and Yis are i.i.d. with probability
density functions (pdfs) f(x) and g(y), respectively.
In this section, we use MCS to estimate reliability. The programs are developed in
MATLAB, separately for exponential, normal and gamma. First a set of 5000 values of the
particular r.v. viz. (exponential, normal or gamma) are generated for a particular value of the
parameter(s). Using these values an estimate of the parameters involved is obtained. Substituting
this estimate(s) in the expression of reliability we get an estimate of the reliability. This process is
repeated j times to give j estimates of the parameter(s) and subsequently j estimates of reliability.
The whole process is repeated for different true values of the parameters; for a particular true value
of the parameter(s) j is the sample size.
2.1 (a) Exponential Stress-Strength:
Let us assume that the components strength follows exponential distribution with mean
and the stress on it follows exponential distribution with mean strength unity, without loss of
generality. Then the marginal reliability expression due to the n th component is, (ibid)
n 1
1
R(n) = R(n) =
(2.3)
1 1
So,
R1 =
(2.4)
1
R2 = R1 + (1 R1) R1,
(2.5)
2
R3 = R1 + (1 R1) R1+ (1 R1) R1.
(2.6)
For MCS, let U be the uniform r.v. over (0, 1). Then by following inverse transformation we
can generate exponential random variable with mean as:
Let
U = F(x) = 1 exp (x / )
X = log (1 U)
Now if U is uniform over (0, 1), (1 U) is also uniform over (0, 1). So
X = log (U)
(2.7)
From uniform r.v. U we can generate exponential r.v. with parameter using the above
transformation (2.7). We generate 5000 of U. Then from (2.7), for each U, log U gives a value of
the exponential r.v. X with mean unity. Thus we get 5000 values of X. Multiplying each of these
5000 values of X by ( 0.5, 2,3) we get 5000 values of exponential r.v. (say X1) with mean . The
mean of these 5000 values give an estimate of for a particular true value of . Substituting these
estimates in (2.4), (2.5) and (2.6) we get an estimate of R1, R2 and R3, respectively. For each true
value of the whole process is repeated j times there by giving j estimates of and Rs for a
particular true value of . Here, we have taken j = 20.
(b) Normal Stress-Strength: In case of normal stress-strength let X N , 2 and stress Y ~N (0,
1) by without loss of generality. The reliability expressions are (ibid)
n1
R (n) = [1
.
(2.8)
]
2
2
1

R1=
,
(2.9)
2
1
and R2 and R3 are given by (2.5) and (2.6), respectively.
80
RT&A # 02 (25)
(Vol.1) 2012, June
For generating normal random numbers for particular true values of and 2 , we first
generate standardized normal numbers (Z) by the MATLAB. Next by the following transformation
we generate normal random numbers particular true values of and 2 .
X = + Z where Z ~ N (0, 1).
We estimate and 2 from the sample of X of size 5,000. Substituting these estimates in
R1, R2 and R3, given above, we get estimates of the system reliabilities Rs. This process is repeated
j times for a particular set of set of true values of and 2 . We have taken ( , ) = (-1, .5), (0, .5),
(1, .5), (2, .5), (-1,1), (0, 1), (1, 1), (2, 1), (-1, 2), (0, 2), (1, 2), (2, 2) and j = 20.
(c) Gamma Stress-Strength: We assume that the strength X 1, m and stress Y 1, k .
Then if m and/ or k is an integer (ibid)
m1
(m k i 1)
R(n) = 1
m k i1
i0 k(m i 1)!2
m1
So,
R1 =
n 1
m1
(m k i 1)
k(m i 1)!2mk i1
(2.10)
i0
(m k i 1)
k(m i 1)!2mk i1
(2.11)
i0
Here also R2 and R3 are given by (2.5) and (2.6).

If m and k are not necessarily integers (Kapur and Lamberson, [9])
n 1
R(n) = 1 Beta(m, k)Beta1/ 2 (m, k) Beta(m, k)Beta1/ 2 (m, k)

(2.12)
where Beta (. , .) is Beta function and Beta1/2(. , .) is incomplete beta function, with
parameters m and k.
So,
R1 = Beta(m, k)Beta1/ 2 (m, k)
(2.13)
and then R2 and R3 are given by (2.5) and (2.6).
Here also, by MATLAB, we directly generate gamma random numbers X for strength
population and Y for stress population, of sizes 5000 each with different true values of the
parameters m and k. Substituting these estimates in the above expressions of Rs we get a
reliability estimate of the systems. This process is repeated j times for a particular set of true values
of m and k. Then for each set of true values of m and k the above process is repeated. We have
taken (m, k) = (1,1), (1,2), (2,1), (2,2) and j = 20.
3. Fitting of Normal Distribution to Systems Reliability:
From different expressions of reliability in Section-2 we have obtained the estimates of
reliability substituting the estimated values of the parameters in respective cases. To these estimated
values of reliability for different cases we have fitted normal distribution and tested the goodness of
fit by one sample K-S test. The tabulated value of D for sample size 20 at 5% level of significance
is 0.294 (see Seigel [18]).
Let us first consider the case of exponential stress-strength. For each ( 0.5, 2,3) j =20, the
values of R 1 , R 2 and R 3 are obtained by substituting the corresponding values of obtained in

Section-2, in the expressions (2.4), (2.5) and (2.6). Then for each value of true , mean and s.d. of
R 1 , R 2 and R 3 are calculated. In each case normal distribution is fitted and the goodness of fit is
tested by K-S test. The values are tabulated in Table-3.1. True values of R1, R2, and R3 are also
given in the same table for comparison.
For K-S test, if calculated value is of D < 0.294, the fit is good. From Table- 3.1, columns 5,
9, 13 it is clear that normal distribution gives good fit to the values of R 1 , R 2 and R 3 .
81
RT&A # 02 (25)
(Vol.1) 2012, June
Table -3.1: Exponential Stress-Strength

True
True
R1
Mean
SD
R1
R1
D
for
True
R2
Mean
SD
D
for
R2
R2
.333
.500
.667
.750
.333
.500
.667
.750
.003
.003
.003
.003
.090
.072
.076
.071
Mean
SD
D
for
R3
R1
.5
1
2
3
True
R3
R3
R2
.555
.750
.889
.938
.556
.750
.889
.948
.004
.003
.002
.001
R3
.067
.072
.063
.088
.704
.875
.926
.953
.704
.875
.963
.984
.004
.002
.001
.000
.116
.075
.063
.023
N.B.: The entry .000 in the SD column indicates that the SD is very small. This is the situation for
all the tables.
Next let us consider the case of normal stress-strength. The above procedure is repeated for
different set of ( , 2 ) and their corresponding estimated values from Section-2 are used in
expressions (2.9), (2.5) and (2.6). The results are tabulated in Table- 3.2. From values of D (see
Seigel [18]) in column 6, 10, 14 we see that normal distribution gives good fits to the distributions
of R 1 , R 2 and R 3 .
Table -3.2: Normal Stress-Strength
True
True
True
R1
Mean
SD
R1
R1
D
for
True
R2
Mean
R2
SD
R2
R1
.5
-1
0
1
2
-1
0
1
2
-1
0
1
2
.212
.500
.788
.945
.308
.500
.692
.841
.421
.500
.579
.655
.186
.500
.814
.963
.240
.501
.760
.922
.327
.501
.672
.813
.001
.003
.002
.001
.003
.005
.004
.002
.004
.005
.004
.005
.073
.800
.083
.037
.141
.064
.047
.123
.077
.061
.107
.075
D
for
True
R3
Mean
R3
SD
R3
R2
.379
.750
.955
.997
.522
.750
.905
.975
.665
.750
.823
.899
.337
.750
.966
.999
.424
.751
.942
.994
.547
.751
.892
.965
.002
.003
.000
.000
.004
.005
.002
.000
.006
.005
.003
.002
.117
.080
.073
.021
.115
.064
.043
.105
.044
.039
.078
.082
D
for
R3
.510
.875
.991
.999
.669
.875
.971
.996
.805
.875
.926
.977
.460
.875
.994
.999
.561
.876
.986
.999
.695
.876
.965
.993
.003
.002
.000
.000
.004
.003
.000
.000
.006
.004
.001
.001
.119
.091
.149
.058
.205
.077
.026
.129
.088
.055
.129
.097
Here we would like to point out that for = 0 and = 0.5, the fit is not good.
Similarly for gamma stress-strength, for different sets of true values of stress-strength
parameters (m, k) and taking their corresponding estimates from Section-2 and using this in
expressions (2.11), (2.5) and (2.6) we obtain estimates of R1, R2, and R3 in different situations and
calculate D statistics in each case. All these values are tabulated in Table- 3.3. Comparing the
values of D in columns 6, 10 and 14 with the tabulated values (ibid) we see that normal distribution
gives good fit to reliabilities of systems for gamma stress-strength also.
82
RT&A # 02 (25)
(Vol.1) 2012, June
Table-3.3: Gamma Stress-Strength (m and/ or k are Integer)

True
m
True
k
True
R1
Mean
SD
R1
R1
D
for
True
R2
Mean
R2
SD
R2
R1
1
2
1
2
1
2
.500
.250
.750
.500
.500
.250
.745
.500
.007
.003
.012
.004
.088
.118
.103
.092
D
for
True
R3
Mean
R3
SD
R3
R2
.750
.438
.938
.750
.750
.438
.935
.750
.007
.005
.006
.004
.063
.126
.087
.092
D
for
R3
.875
.579
.984
.875
.875
.579
.983
.875
.005
.006
.002
.003
.085
.123
.119
.102
When neither m nor k is an integer then the corresponding estimates of R1, R2, and R3 are
obtained by substituting the values of m and k from Section-2 in the expressions (2.13) etc. and the
corresponding values are tabulated in Table-3.4.
Table- 3.4: Gamma Stress-Strength (m and K are not necessarily Integer)
True
m
True
k
True
R1
Mean
SD
R1
R1
D
for
True
R2
Mean
R2
SD
R2
R1
1
2
1
2
1
2
.500
.250
.750
.500
.499
.250
.750
.500
.007
.005
.005
.006
.090
.092
.069
.095
D
for
True
R3
Mean
R3
SD
R3
R2
.750
.438
.938
.750
.749
.438
.938
.750
.007
.007
.002
.006
.090
.005
.060
.070
D
for
R3
.875
.579
.984
.875
.874
.579
.984
.875
.006
.008
.001
.004
.076
.065
.071
.096
Conclusion: In this paper, we have estimated the reliability through MCS. We have seen that
normal distribution is fitted to the data sets of estimated reliability obtained by MCS. Once we
know the distribution it is easy for us to obtain the other characteristics of the reliability data sets.
Acknowledgment: This paper was supported by UGC (India) under Major Research Project (Sc.).
References:
[1]
[2]
[3]
[4]
[5]
[6]
Aldrisi, M. M., (1987): A Simulation Approach for Computing System Reliability,

Microelectronics Reliability, Vol. 27, Issue: 3, pp. 463-467.
Ahmad, K.E., Fakhry, M.E. and Jaheen, Z.F., (1997): Empirical Bayes Estimation of
P(X>Y) and characterizations of Burr-type X model, Journal of Statistical Planning and
Inference, Vol. 64, Issue: 2. pp. 297-308.
Borhan Uddin, Md., Pandey, M., Ferdous, J. and Bhuiyan, M.R., (1993): Estimation of
Reliability in a Multicomponent Stress-Strength Model, Microelectronics and Reliability,
Vol.33, Issue: 13, pp.2043-2046.
Church, J.D. and Harris, V., (1970): The Estimation of Reliability from Stress-Strength
Relationship, Technometrics, Vol.12, No.1, pp.49-54.
Goel, L., (2000): Monte Carlo Simulation Based Reliability Studies of a Distribution Test
System, Electric Power System Research, Vol. 54, Issue: 1, pp.55-65.
Hong, H.P. and Lind, N.C., (1996): Approximate Reliability Analysis using Normal
Polynomial and Simulation Results, Vol.18, Issue:4, pp.329-339.
83
[7]
[8]
[9]
[10]
[11]
[12]
[13]
[14]
[15]
[16]
[17]
[18]
[19]
[20]
[21]
[22]
[23]
[24]
[25]
RT&A # 02 (25)
(Vol.1) 2012, June
Kakati, M.C. and Sriwastav, G.L., (1981): Accelerated Life-Testing under Stress-Strength
Model, IAPQR Trans. Vol.11, No.1-2, pp.51-56.
Kamat, S.J. and Riley, M.W., (1975): Determination of Reliability Using Event-Based Monte
Carlo Simulation, IEEE Trans. on Reliab., Vol. R-24, No.1, pp.73-75.
Kapur, K.C. and Lamberson, L.R., (1977): Reliability in Engineering Design, John Wiley
and Sons, New York.
Landis, C.M., Beyerlein, I.J. and McMeeking, R.M., (2000): Micromechanical Simulation of
the Failure of Fiber Reinforced Composites, Journal of the Mechanics and Physics of Solids,
Vol.48, Issue:3, pp.621-648.
Manders, W.P., Bader, M.G. and Chou, T.W., (1982): Monte Carlo Simulation of the
Strength of Composite Fiber Bundles, Fiber Science and Technology, Vol. 17, Issue:3,
pp.183 204.
Mazumder, M.,(1970): Some Estimates of Reliability Using Interference Theory, NRLQ,
Vol.17, pp.159-165.
Naess, A., Leira, B.J. and Batsevych, O., (2009): System Reliability Analysis by Enhanced
Monte Carlo Simulation, Structural Safety, Vol.31, Issue: 5, pp.349 355.
Paul, R.M. and Borhar Uddin, MD., (1997): Estimation of Reliability of Stress-Strength
Model with Non-identical Component Strengths, Microelectronics and Reliability, Vol. 37,
Issue: 6, pp.923-927.
Pulido, J.E., Jacobs, T.L. and Parters de Lima, E.C., (1992): Structural Reliability Using
Monte- Carlo Simulation with Variance Reduction Techniques on Elastic-Plastic Structures,
Computers and Structures, Vol. 43, Issue: 3, pp.419-430.
Rao, G.S. and Kantam, RRL, (2010): Estimation of Reliability in Multicomponent StressStrength Model: Log Logistic Distribution Eletron. J. App. Stat. Annal., Vol.3, Issue: 2 pp.
75-84.
Rezaei, S., Tahmasbi, R. and Mahmoodi, M., (2010): Estimation of P [ Y < X ] for
Generalized Pareto Distribution, Journal of Statistical Planning and Inference, Vol. 140,
Issue: 2, pp. 480-494.
Seigel, S. (1956): Non-Parametric Statistics for Behavioural Sciences, McGraw Hill Book.
Co. Ltd., N.Y.
Sriwastav, G.L. and Kakati, M.C., (1981): A Stress-Strength Model with Redundancy,
IAPQR Trans. Vol.6, No.1, pp.21-27.
Sriwastav, G.L., (1988): Reliability Estimation in Stress-Strength Model with Ancillary
Information, Assam Statistical Review., Vol. 2, No.1-2, pp.16-20.
Stancampiano, P.A., (1978): Simulation Approaches to Probabilistic Structural Design at
Component Level, Clear Engineering, Vol. 50, Issue:1, pp.123-129.
Stumpf, H. and Schwartz, P., (1993): A Monte-Carlo Simulation of the Stress-Rupture of
Seven-Fiber Microcomposites, Vol.49, Issue:3, pp.251-263.
Tunak, M., Linka, A. and Volf, P., (2009): Load-Sharing and Monte Carlo Models of Defects
in a Bundle fibres, Composites Science and Technology, Vol.69, Issue: 9, pp. 1417-1421.
Wu, L., Yi, L., Xue-qin, G. and Deng-hai, W., (2009): Monte-Carlo Simulation for the
Reliability Analysis of Multi-Status Network System Based on Breadth First Search,
Information and Computing Science, pp.280-283.
Zhang, H., Chandrangsu, T. and Ramussen, K.J.R., (2010): Probabilistic Study of the
Strength of Steel Scaffold System, Structural Safety, Vol.32, Issue: 6, pp.393-401.
84
Chacko V.M., Manoharan M. MULTISTATE COHEREHT SYSTEMS WITH MULTIPLE STATE TRANSITION AT A TIME

RT&A # 02 (25)
(Vol.1) 2012, June

MULTISTATE COHERENT SYSTEMS WITH MULTIPLE STATE

TRANSITION AT A TIME
V. M. Chacko
Department of Statistics, University of Calicut,
Kerala -673635, India
chackovm@gmail.com
M. Manoharan
Department of Statistics, University of Calicut,
Kerala-673635, India
Abstract
A key requirement in defining a multistate coherent system (MCS) is the relevance
condition of its components. A new class of MCSs is introduced with a new component relevance
condition. Also we introduce a more general relevance condition. They are compared with some
existing component relevance conditions. Based on the two new relevance conditions, two
component importance measures for MCSs are defined. They are most appropriate for comparing
components when certain type of system improvement is sought. We introduce new joint
importance measures for two or more components with respect to the proposed relevance
conditions. The new MCS classes include several existing MCSs as special case. An illustrative
example of the proposed MCSs is also provided.
Keywords: Reliability, MCS, relevance condition, component importance, joint importance.
1. Introduction
Let us consider a coherent system with n components C={1,2,,n}. Furthermore, suppose

that each component can be in one of M+1 states, {0,1,2,,M}, where 0 is the failed state and
M is the maximal or perfect state. To describe such a multistate system (MSS), a general
theory has been developed in the literature. 5,9,11 A binary state system (BSS) of n components can
be described by a structure function : {0,1}n {0,1}, which presents the state of system as a function
of states of its n components. 4 A binary system is statistically coherent if it satisfies the following
conditions; 4
(i) (x ) is non-decreasing in each argument, where x ( x1 , x2 ,..., xn )
with xi {0,1}, and
(ii) for each i, there exist a vector (.i , x ) , such that (1i , x ) (0i , x ),
where (.i , x ) ( x1 , x2 ,..., xi 1,.i , xi 1,..., xn ).
Note that the condition (i) and (ii) gives, ( j ) j, j 0,1 where j ( j, j ,..., j ) .
In practice, a system and its components often have more than two states of performance. 5
The structure function of the MSS is : S n S , S {0,1,..., M } , which relates the level of performance
of system to level of performance of each of its components. There are various approaches which
extends the structure function from the binary case to the multistate case. 5,9,10,11 The effort resulted
in, extension of the requirement of non-decreasing binary structure function to MSS structure
function. Also note the condition ( j ) j, j 0,1 of the binary coherent system (BCS) is extended to
85

RT&A # 02 (25)
(Vol.1) 2012, June

the MSS requiring ( j ) j, j {0,1,..., M } . 11 The condition (ii) of relevancy in BCS is extended in
various different ways. Some extensions can be seen in Refs. 1, 2 and 12.
In this paper, we extend the relevance condition to the MSS case in a general way, which
includes several existing relevance conditions as special cases. Section 2 introduces the new class of
multistate coherent system(MCS)s and its generalization by introducing a reasonable component
relevance condition. The two new classes are compared with the some existing classes. Section 3
introduces two new component importance and joint importance measures for the proposed MCSs.
Section 4 provides an example of an offshore electrical power generation system. Discussion and
conclusion are given in section 5.
2. Component relevancy and the new classes of MCSs
In this section we discuss the new relevance condition and its generalization on which two
new classes of MCSs are defined. Consider the following component relevance conditions.
NAT 13 : For every component i and level j>0, there exist (.i , x ) such that
( ji , x ) j and (( j 1)i , x ) j.
GRI.1 11 : For every component i and level j>0, there exist
GRI.2 11 : For every component i, there exist
and
(.i , x )
such that ( ji , x ) (( j 1)i , x ).
such that (0i , x ) (M i , x ).

EP 9 : For every component i and level j1, there exist (.i , x ) such that
(.i , x )
( ji , x ) (0i , x ).
NAT and GRI.1 indicate degree of relevance of each component to every level of
performance; while GRI.2 merely states that is not a constant in any of its arguments.
Now consider a situation in which some component is not relevant to every level of
performances, i.e., the system degrades from state j to j-1 or j-2 etc when the component degrades
only from state j to j-2 or j-3 etc. In order to degrade the system, component must degrade more
than one level of performance. For example, 14 let S {0,1,2,3,4}, and the component can take 0, 2,
and 4 when the system can take 0, 1, 2, 3 and 4. Consider the structure function 2 having 5
components in Ref.14. From the minimal path vectors of 2 , we have,
th
2 (41 ,42 ,23 ,44 ,25 ) 4 2 (41 ,42 ,23 ,24 ,25 ) 3 , when the 4 component degrades from state 4 to state 2, the
system degrades from state 4 to state 3. Now consider the structure function 1 with three
components in Ref.14. We have, 1 (41 ,02 ,43 ) 4 1 (41 ,02 ,23 ) 2, when the third component degrades
from state 4 to state 2, the system degrades from state 4 to state 2. Here fourth component must
degrade from state 4 to state 2 for the system to degrade from state 4 to state 3 with respect to 2 .
The third component must degrade from state 4 to state 2 for the system to degrade from state 4
with respect to 1 .
We define a new component relevance condition as, degrading a component from state j to
state j-2 can cause system failure or degradation while degradation of the component from state j
to j-1 cannot cause system failure or degradation.
Now the new class of MCSs can be defined as follows.
Definition.1.: A multistate system of n components with structure function belonging to
class CM.1 if is non-decreasing, ( j ) j, and for each component, there exist (.i , x ) such that
( ji , x ) (( j 2)i , x ).
Now consider the generalization of the new relevance condition, one or more than one level
of degradation of the component can cause the system degradation, i.e., when the component i
degrades from state j to state j {j-1, j-2, j-3,., 1, 0}, the system degrades from state j to any
lower state. Thus we define the generalized class of MCSs with this relevance condition.
86

RT&A # 02 (25)
(Vol.1) 2012, June

Definition 2.: A multistate system of n components with structure function belonging to

class CM.2 if is non-decreasing, ( j ) j, and for each component, there exist (.i , x ) such that
( ji , x ) ( j 'i , x ) where j {j-1, j-2,j-3,.,1,0}.
In the following section we introduce the component importance and joint importance
measures to the new classes of MCSs.
3. Component importance and Joint importance measures
We consider the problem of measuring the reliability importance and structural importance
of individual components, and the joint reliability importance and joint structural importance of two
or more components in the new classes of the MCSs. The main advantage of defining a new
relevance condition is to obtain the importance measures. 6 At the reliability design phase, the joint
importance can improve system designers understanding of the relationship between the
components and the system, and among the components, 3 which is quite desirable. Birnubaum
measure provides the importance of a component in the BSS. 6 It is further extended to the MSS. 7 ,15
Now we consider X ( X 1 , X 2 ,..., X n ) as a random vector with component states X i as random
variables and pij Pr{ X i j} where j S {0,1,..., M }. For the BSS with structure function , the
Birnubaum reliability importance 6 of component i is
I i ( B) P( (1i , x ) (0i , x ) 1) h(1i , p ) h(0i , p ).
where h( p ) , p ( p1 , p2 ,..., pn ) and i pi ( xi 1) , is the reliability function of the BCS,

h( p ) pi h(1i , p ) (1 pi )h(0i , p ) pi I i ( B) h(0i , p ).
Therefore,
h( p )
I i ( B).
pi
We propose the following component importance measures for the two classes of new
MCSs.
1. I i (CM .1) P( ( ji , x ) (( j 2)i , x )).
2. I i (CM .2) P( ( ji , x ) ( j 'i , x )), j ' { j 1, j 2,...,1,0}.
Let the distribution of X i be described by pi ( pi 0 , pi1 ,..., piM ). The reliability function of the
MCS with minimum satisfactory system level j , is P( ( x ) j ) P( ( ji , x ) j ), since
j S
pi 0 pi1 pi 2 ... piM 1 .
Now we prove the following theorems.
Theorem 1. For the CM.1 class,

respect to pij .
I i (CM .1)
is the rate of improvement of P( ( x ) j ) with
Proof. Clearly,
P ( ( x ) j )
p [ P( ( j , x j ) P( (( j 2) , x ) j )] P( (( j 2) , x ) j ),
S \{ j 2}
since
ij
1 pij 2 pi 0 pi1 ... pij 3 pij 1 ... piM .
Differentiating P( ( x ) j ) partially with respect
to pij , we get
I i (CM .1)
P( ( x ) j )
P( ( ji , x ) j ) P( (( j 2) i , x ) j ) P( (( j 2) i , x ) ( ji , x )).
pij
87
Theorem 2. For the CM.2 class,

respect to pij .
I i (CM .2) is

RT&A # 02 (25)
(Vol.1) 2012, June

the rate of improvement of P( ( x ) j ) with
Proof. Clearly,
P( ( x ) j )
p [ P( ( j , x ) j) P( ( j' , x ) j)] P( ( j ' , x ) j ),

ij
S \{ j ' }
where
1 pij ' pi 0 pi1 ... pij '1 pij '1 ... piM and
differentiating the P( ( x ) j ) partially with
respect to pij we get,

I i (CM .2)
P( ( x ) j )
P( ( j 'i , x ) ( ji , x )),
pij
j '{ j 1, j 2,...,1,0}.
Now define the structural definition of the component importance (when reliabilities of
components are not given) with respect to the new relevance conditions.
Consider ( x ) 1 if ( x ) j and 0 otherwise. We define the structural importance of a
component as follows.
Definition.3.: Let : S n S be the MCS structure function in CM.1 class. Then is said to
have the following measures of structural importance for the level j of component i:
I ij (CM .1)
1
Max{0, ( ji , x ) (( j 2)i , x )}.
( M 1) n 1 { x : xi j}
Definition.4.: Let : S n S be the MCS structure function in CM.2 class. Then is said to
have the following measures of structural importance for the level j of component i:
I ij (CM .2)
1
Max{0, ( ji , x ) ( j 'i , x )},
( M 1) n 1 { x : xi j}
j ' { j 1, j 2,...,1,0}.
In order to define the joint importance measures for two or more components in the new
classes of MCSs, we recall the joint structural importance measure(JSIM)s 8 and joint reliability
importance measure(JRIM)s 8 for the MSS with relevance condition in GRI.1. The JSIM (i,j) for
two components i and j with the new relevance conditions can be obtained by replacing m with m-2
or m' {m 1, m 2,...,2,1,0} in,
M M
JSIM ( i ,l ) { SIM ( i ,l ; m,k ) SIM ( i ,l ; m,k )}

m 1 k 1
( (m , k , x
where SIM (i, l ; m, k ) = X
il
q 1
il
) j, (mi , kl , xil ) j q)
( M 1)n 2
Here (true)=1 and (false)=0, and (mi , kl , xil ) j, (mi , kl , xil ) j q where xij ( x1 ,..., mi ,..., kl ,..., xn )
determines the critical path vector to the level j with state m of component i. The JSIM (i,j,k) for
three components can be obtained as, for n n 2 or n {n 1, n 2,...,2,1,0} ,
M
JSIM (i, l , r ) {JSIM (i, l , r ; m, k , n) JSIM (i, l , r ; m, k , n)},

k 1 n 1 m 1
where JSIM ( i ,l , r ; m , k , n ) = JSIM ( i , l , r ; m , k , n ) JSIM ( i , l , r ; m , k, n ).

Thus we can find JSIM of any number of components w. r. t. both relevance conditions in
the MCS classes CM.1 and CM.2. Thus JSIM 8 holds with new relevance conditions by replacing
m with suitable m-2 or m' {m 1, m 2,...,2,1,0} .
Now we consider the JRIM for k components in the MSS. 8 The JRIM of state b1 of
component a1 , state b2 of component a ,..., state bk of the component a k (k n) of the MSS is
2
88

RT&A # 02 (25)
(Vol.1) 2012, June

JRIM ( a1, ,..., ak ; b1 ,..., bk )
k Es
, k 2,..., n,
Ra1 b1Ra 2 b2 ...Ra k bk
where Es P( ( x) j ) is the expected system performance and Ra bi is the reliability function with
i
j 0
respect to level bi of component ai .

Here m =m-1 in the expansion of Es . 8 The results also holds true with the other values of
m =m, m {m-1,m-2,m-3,,0}. Hence the JRIM with the new MCSs can be obtained with
appropriate m values for m , i.e., m=m-2 or m' {m 1, m 2,...,2,1,0} .
Now consider some implications based on the new relevancy definitions. In fact one can
easily prove the following implications.
Theorem 3. CM.1=>EP, CM.2=>NAT=>GRI.1=>GRI.2=>EP,
CM.2=>GRI.1,
CM.2=>GRI.2 and CM.2=>EP.
It is clear that all the existing relevance conditions are special cases of CM.2 (or CM.1)
relevance condition. Hence the existing MCSs are special cases of the proposed MCSs.
4. Example
Ref.14 considered an offshore electrical power generation system, which supply two nearby
oilrings with electrical power. Both oilrings have their own main generation, represented by
equivalent generators A1 and A3 each having capacity of 50MW. In addition the oilrings has a
standby generator A2 that is switched into the network in case of outage of A1 or A3 , or may be used
in extreme load situations in either of the oilrings. The A2 also has capacity 50MW. The control
unit, U , continuously supervises the supply from each of the generators with automatic control of
the switches. If for instance the supply from A3 to oilring 2 is not sufficient, whereas the supply
from A1 to oilring 1 is sufficient, U can activate A2 to supply oilring 2 with electrical power through
the subsea cables L . The components have states {0, 2, 4} and the system have states {0, 1, 2, 3,
4}, where 0, 1, 2, 3 and 4 represents the states of the system at capacities 0MW, 12.5MW,
25MW, 37.5MW, and 50MW respectively.
Table I .Minimal path vectors of
U
2
4
2
4
4
Levels
2
2
4
4
4
A1
2
0
4
0
2
A2
0
2
0
4
2
Table II .Minimal path vectors of 2

Levels
A1
A2
A3
1
1,2
2
2
3
3,4
4
2
4
4
4
2
4
0
4
4
4
0
2
0
2
4
2
0
2
0
4
2
2
0
0
2
0
0
2
4
89

RT&A # 02 (25)
(Vol.1) 2012, June

4
4
4
4
4
4
4
4
4
2
4
4
4
2
4
2
2
0
The minimal path vectors to the levels are given in table 1 and table II, of the structure
functions,
1 (U , A1 , A2 ) I (U 0) min( A1 A2 I (U 4),4), the amount of power that can be supplied to
platform 1, and
the amount of power that can be
2 (U , A1 , L, A2 , A3 ) I (U 0) min( A3 A2 I (U 4) I ( A1 4) L / 4,4),
supplied to platform 2, when I (.) is the indicator function. One may easily verify from the tables
that the new relevance condition of CM.1 and CM.2 are found to be holding good w. r. t. the
structure functions considered in the example. As done for JSIM and JRIM for two components 11
and JSIM for three components, 8 we can compute the concerned joint importance measures for any
number of components in the power generation system.
5. Discussion and Conclusion
The theory of MSS reliability models has been developed to cope with many real-life
situations. The present paper introduced two classes of MCSs with a new relevance condition and
its generalization. It is shown that many MCSs introduced earlier in the literature are included in the
new classes. Structural definitions of importance and joint structural importance measures are
given, and new reliability importance and joint reliability importance measures are introduced. In
system engineering, a practical and difficult problem is the identification of those groups of
components that mostly influence the system behavior with respect to safety and reliability. In this
respect, the main advantage of our importance and joint importance measure with respect to the new
MCS models is the information provided by them for the reliability theoreticians and design
analysts. It gives useful information for safe and efficient operation of the system, where existing
importance measures gives information about individual component importance and joint
importance of components with some limited number of relevance conditions.
Acknowledgments
Authors are thankful to the financial support from CSIR 1 and UGC 2 , India.
References
1. A. M. Abouammoh, and Al-kadi, On measures of importance for components in
multistate coherent system, Microel. Reliab., 31(2) (1991) 109-122.
2. A. M. Abouammoh, and Al-kadi, Multistate coherent system of order k, Microel. Reliab.,
35(11) (1995) 1415-1421.
3. M. J. Armstrong, Joint reliability importance of components, IEEE Trans. Reliab., 44(3)
(1995) 408-412.
4. R. E. Barlow, and F. Proschan, Statistical theory of reliability and Life testing, (Holt,
Rinheart & Winston, NewYork, 1975).
5. R. E. Barlow, and A. Wu, Coherent system with multistate components, Math. Oper.
Res., 3 (1978) 275-281.
6. Z. W. Birnbaum, On the importance of different components in a multi-component
system, in Multivariate analysis II, ed. P. R. Krishnaiah (NewYork, Academic Press, 1969), p. 591592.
7. V. C. Bueno, On the importance of components for multistate monotone systems, Stat.
Prob. Lett., 7 (1989) 51-59.
90

RT&A # 02 (25)
(Vol.1) 2012, June

8. V. M. Chacko and M. Manoharan, Joint importance measures for the multistate system,
Advances in Performance and Safety of Complex systems, eds. A. K., Verma, P. K, Kapur, S. G.,
Ghadge, (Macmillan, New Delhi, 2008), p. 308-314.
9. E. El-Neweihi and F. Proschan, Degradable systems; a survey of multistate system
theory, Commu. Statist. Theory and Methods, 13 (1984) 405-432.
10. E. El-Neweihi, F. Proschan, and J. Sethuraman, Multistate coherent systems, J. Appl.
Prob., 15 (1978) 675688.
11. W. S. Griffith, Multistate reliability models, J. Appl. Prob., 17(1980) 735744.
12. F. C. Meng, Component relevancy and characterization results in multistate systems,
IEEE Trans. Reliab., 42 (1993) 478-483.
13. B. Natvig, Two suggestions of how to define a multistate coherent system, Adv. Appl.
Prob., 14 (1982) 434-455.
14. B. Natvig, S. Sormo, A. T. Holen, and G. Hogasen, Multistate reliability theory -a case
study, Adv. Appl. Prob, 18(1986) 921- 932.
15. S. Wu, Joint importance of multistate systems, Comp. Indust. Eng., 49(2005) 6375.
91

Rta 2 2012

Uploaded by

Copyright:

Available Formats

Rta 2 2012

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Rta 2 2012

Uploaded by

Copyright:

Available Formats

ELECTRONIC JOURNAL

Gnedenko Forum Publications

VOL.1 NO.2 (25)

All rights are reserved

Journal of International Group on Reliability

Send your paper

Ushakov, Igor (USA)

e-Journal Reliability: Theory & Applications

Theoretical papers have to contain new problems,

Priority is given to descriptions of case studies.

Koowrocki, Krzysztof (Poland)

Krishnamoorthy, Achyutha (India)

General requirements for presented papers

Shybinsky Igor (Russia)

1. Papers have to be presented in English in

Singpurwalla, Nozer (USA)

Chakravarthy, Srinivas (USA)

Dimitrov, Boyan (USA)

Genis, Yakov (USA)

Kaminsky, Mark (USA)

Kovalenko, Igor (Ukraine)

Levitin, Gregory (Israel)

Limnios, Nikolaos (France)

Nelson, Wayne (USA)

The authors keep all rights and after the publication

Popentiu, Florin (UK)

Rykov, Vladimir (Russia)

Wilson, Alyson (USA)

Wilson, Simon (Ireland)

Yastrebenetsky, Mikhail (Ukraine)

Zio, Enrico (Italy)

the Deputy Editor,

Fedol Ghazi, Pr: Hamouine Abdelmadjid, Dr: Zouaoui Chikr el Mezouar

I.Ushakov - SENSITIVITY ANALYSIS OF OPTIMAL REDUNDANCY SOLUTIONS

SENSITIVITY ANALYSIS OF OPTIMAL REDUNDANCY SOLUTIONS

San Diego, California, USA

Figure 1. Series system underwent analysis

And the total systems redundant units cost is:

Below presented solutions of both problems of optimal redundancy: direct:

Ushakov I. - SENSITIVITY ANALYSIS OF OPTIMAL REDUNDANCY SOLUTIONS

Table 3. Solution for Inverse problem

Figure 2. Input data for the first series of experiments

Figure 3. Input data for the second series of experiments

Ushakov I. - SENSITIVITY ANALYSIS OF OPTIMAL REDUNDANCY SOLUTIONS

Table 4. Values of Probabilities of Failure-free operations

Table 6. Results of Monte Carlo simulations for P*=0.995

Ushakov I. - SENSITIVITY ANALYSIS OF OPTIMAL REDUNDANCY SOLUTIONS

Table 7. Results of Monte Carlo simulations for P*=0.99

Table 7. Results of Monte Carlo simulations for P*=0.95

Figure 4. Deviation of maximum and minimum values of probability

Ushakov I. - SENSITIVITY ANALYSIS OF OPTIMAL REDUNDANCY SOLUTIONS

Figure 5. Multiple polygon axes with numbered levels

Figure 6. Deviations of optimal solutions for randomly variation of parameters

TOPOLOGICAL SEMI-MARKOV METHOD FOR CALCULATION OF

Doctor of Technical Sciences, professor, director

Ph.D., The Head of Scientific-Technical Complex of JSC "NIIAS",

determined according to the study [2] are the following:

Application of Masons formula allows to considerably reducing labor input of calculations

3.2. Topological formulas of calculation of availability and functional safety factors

p (j2) ; ci(1) c (j2) ; i (1) j ( 2) p ( 2) , c ( 2) , j ( 2 ) )