Kołowrocki Krzysztof, Kuligowska Ewa, Soszyńska-Budny Joanna: Monte carlo simulation for optimization of object operation process and reliability. (Symulacja monte carlo w optymalizacji procesu eksploatacji oraz niezawodności obiektu technicznego. )

(1)

MONTE CARLO SIMULATION FOR OPTIMIZATION

OF OBJECT OPERATION PROCESS

AND RELIABILITY

SYMULACJA MONTE CARLO W OPTYMALIZACJI

PROCESU EKSPLOATACJI ORAZ NIEZAWODNOŚCI

OBIEKTU TECHNICZNEGO

Krzysztof Kołowrocki, Ewa Kuligowska, Joanna Soszyńska-Budny

Maritime University, Faculty of Navigation e-mail: e.kuligowska@wn.am.gdynia.pl

Abstract: This paper presents the computer simulation technique related to the

reliability of an object under variable operation conditions. The Monte Carlo simulation is performed and, consequently, the unknown parameters for the considered object operation process and its main reliability characteristics are identified. The optimal transient probabilities of the object operation process at the particular operation states and its optimal unconditional reliability function are determined.

Keywords: reliability, operation process, optimization, Monte Carlo simulation Streszczenie: Niniejszy artykuł prezentuje komputerową technikę symulacyjną,

związane z niezawodnością obiektu technicznego w zmiennych warunkach eksploatacyjnych. Przeprowadzona jest symulacja Monte Carlo, a w konsekwencji zidentyfikowane są nieznane parametry badanego procesu eksploatacji obiektu

i jego główne cechy niezawodności. Wyznaczone są optymalne

prawdopodobieństwa przebywania procesu eksploatacji obiektu w poszczególnych stanach eksploatacyjnych oraz jego optymalna bezwarunkowa funkcja niezawodności.

Słowa kluczowe: niezawodność, proces eksploatacji, optymalizacja, symulacja

(2)

1. Introduction

The object reliability function analytical determination very often leads to complicated formulae and therefore it is sometimes difficult to implement modeling, prediction, optimization and cost analysis using this way [2]-[7]. The Monte Carlo simulation method is a tool that sometimes allows to simplify solving this problem [8]. The analytical approach to systems reliability is shortly presented and next the computer simulation modelling methods related to the reliability of an object under variable operation conditions. The Monte Carlo method is used for semi-Markov object operation process simulation to obtain its realizations [1]. It allows examining the reliability of objects at variable operation conditions.

2. An object operation process

We assume that an object during its operation at the fixed moment t, t 0,, may be at one of , N, different operations states zb, b1,2,...,. Consequently, we mark by Z(t), t 0,, the object operation process, that is a function of a continuous variable t, taking discrete values at the set {z1,z2,...,zv }

of the object operation states. We assume a semi-Markov model [2], [4] of the object operation process Z(t) and we mark by _bl its random conditional sojourn times at the operation states

z

_b, when its next operation state is ,z l b,l1,2,...,v, bl.

Consequently, the operation process may be described by the following parameters:

 the vector of the initial probabilities [p_b(0)]1__v of the object operation process Z(t) staying at the particular operations states at the moment

t



0 ;

 the matrix of the probabilities [p_bl]_v__v of the object operation process Z(t)

transitions between the operation states z and _b z , l b,l1,2,...,, bl;

 the matrix of the conditional distribution functions [Hbl(t)]___ of the object operation process Z(t) conditional sojourn times θ at the operation states_bl . Having identified the probabilities p of transitions between the operation states bl

and the distributions of conditional sojourn times θ , the mean values bl M of the b

object operation process Z(t) unconditional sojourn times b, b1,2,...,, at the particular operation states can be determined [4]. Further, the limit values of the object operation process Z(t) transient probabilities pb(t)P(Z(t)zb),

, ,..., 2 , 1  

b at the particular operation states

b

p = limpb(t)

t

(3)

can be determined [4]. Other practically interesting characteristics of the object operation process Z(t) when the operation time  is sufficiently large, are its total sojourn times ˆ at the particular operation states ,_b z b b1,2,...,v, during the fixed

object opetation time that have approximately normal distribution with the expected value given by

, ] ˆ [ ˆ _ _ b b b E p M   b1,2,...,, (2) where p are defined by (1). _b

3. Reliability of an object subjected to varying operation conditions

We assume that every operation state of the object operation process Z(t), ),

, 0  

t described in section 2, have an influence on the object reliability [4]. Therefore, the object reliability at the particular operation state zb, b1,2,...,v, can be described using the conditional reliability function R(b)(t), t 0,



,

, ,..., 2 , 1  

b that is the conditional probability that the object conditional lifetime )

(b

T is greater than t , while the object operation process Z(t) is at the operation state zb, b1,2,..., [4]. Further, we denote the object unconditional lifetime by T and the unconditional reliability function of the object by R(t)P(Tt),

). , 0  

t In the case when the object operation time θ is large enough, the unconditional reliability function of the object is approximated by [4]



  v b b b t p t 1 ) ( ), ( ) ( R R t0,) (3)

where p_b,b1,2,...,, are the object operation process limit transient probabilities given by (1). Hence the mean value of the object unconditional lifetime T is given by



     1 b b b p (4)

where _b are the mean values of the object conditional lifetimes T(b) at the

operation state z_b, b1,2,...,.

4. Monte Carlo simulation approach to an object operation process modelling

We denote by zb(g), b1,2,...,, the realization of the object operation process initial operation state at the moment t0 generated from the distribution defined by the vector [p_b(0)]1_v. This realization is generated according to the formula

(4)

                   (0) 1, ) 0 ( ) 0 ( , ), 0 ( ) 0 ( ) 0 ( , ), 0 ( 0 , ) ( 1 2 1 2 1 1 2 1 1 g p p p z p p g p z p g z g zb      (5)

where g is a randomly generated number from the uniform distribution on the interval 0,1 .

We denote by z_bl(g), l1,2,...,, bl, the sequence of the realizations of the object operation process consecutive operation states generated from the distribution defined by the matrix [p_bl]_v__v. Those realizations are generated according to the formula

                   1, , , , 0 , ) ( 1 1 13 12 13 12 12 3 12 2 1 g p p p z p p g p z p g z g zl      (6)                                          , 1 , , , , , , 0 , ) ( 1 2 1 1 2 1 1 2 1 1 1 2 1 2 2 1 1 1 1 g p p p z p p p g p p p z p p p g p p p z p g z g z b b b bb b b bb b b b bb b b bb b b b b l b            (7) for b2,3,...,,                    , 1, , , 0 , ) ( 2 2 1 1 2 1 1 2 1 1 g p p p z p p g p z p g z g z _l              (8)

where g is a randomly generated number from the uniform distribution on the interval 0,1 .

We denote by bl( j), b,l1,2,...,,bl, j1,2,...,nbl, the realizations of the

(5)

distribution H , defined by the matrix _bl [Hbl(t)]___, where n is the number of bl

those sojourn time realizations during the experiment time ~. Those realizations are generated according to the formulae

) ( 1 h Hbl bl    , b,l1,2,...,,bl (9) where Hbl1(h) 

is the inverse function of the distribution function Hbl(t) and h is a randomly generated number from the uniform distribution on the interval 0,1 . Having those realisations, it is possible to determine approximately the total sojourn time at the operation state z during the time of the experiment _b ~

applying the formula



       b l l n j j bl b b l 1 1 ~ , ₍₁₀₎

and the object operation process unconditional mean sojourn times are given by

b b b n M  1 ~ ,



   1 l bl b n n . (11)

Further, the limit transient probabilities defined by (1) can be approximately obtained using the formula

  ~ ~ b b p  , b1,2,...,,



     1 ~ ~ b b . (12)

5. Optimization of an object operation process and reliability at

variable operation conditions

The object operation process has significant influence on its reliability [4]. According to (4), the mean value , of the object unconditional lifetimes is determined by the limit values of transient probabilities pb,b1,2,...,v, of the object operation process at the operation states given by (

12)

12). The corresponding optimal values p_b,b1,2,...,v, of the transient probabilities may be found to maximize the mean value  of the unconditional object lifetimes.

Therefore, the optimization problem can be formulated as a linear programming model with the objective function of the following form

(6)



     1 b b b p , (13)

with the following bound constraints for the unknown transient probabilities p_b

, b b b p p p    b1,2,...,v,



   1 . 1 b b p (14)

The optimal value of the transient probabilities p can be found after completing _b

a few steps described in [4].

6. Applications

An exemplary object operation process

We consider an exemplary object operating at

v

= 4 operation states z1, z2, z and 3 4

z . The probabilities of the initial operation states of this object operation process

are fixed arbitrarily in the following way

]. 40 . 0 , 29 . 0 , 10 . 0 , 21 . 0 [ )] 0 ( [pb  (15)

The probabilities of the object operation process Z(t) transitions between the operation states z and _b z , l b,l1,2,3,4, bl are also fixed arbitrarily [3]

and given in the matrix below

. 0 30 . 0 22 . 0 48 . 0 72 . 0 0 16 . 0 12 . 0 50 . 0 30 . 0 0 20 . 0 46 . 0 32 . 0 22 . 0 0 ] [              bl p (16)

Moreover, we assume that the distribution functions of the exemplary object operation process conditional sojourn times measured in days are as follows:

], 0052 . 0 exp[ 1 ) ( 12 t t H    H31(t)1exp[0.0011t], ], 0021 . 0 exp[ 1 ) ( 13 t t H    H32(t)1exp[0.0021t], ], 005 . 0 exp[ 1 ) ( 14 t t H    H34(t)1exp[0.0033t], ], 0104 . 0 exp[ 1 ) ( 21 t t H    H41(t)1exp[0.0031t], ], 0123 . 0 exp[ 1 ) ( 23 t t H    H42(t)1exp[0.0020t], ], 0182 . 0 exp[ 1 ) ( 24 t t H    H43(t)1exp[0.0023t], t0,). (17)

(7)

From the above, applying (17), the conditional mean values Mbl E[bl], , 4 , 3 , 2 , 1 ,l

b of the exemplary object sojourn times at the particular operation states measured in days are fixed as follows:

,

192

12



M

₁₃



480 ,

M

₁₄



200 ,

,

96

21



M

23



81 ,

M

24



55 ,

, 870 31  M

M

₃₂



480 ,

M

₃₄



300 ,

,

325

41



M

₄₂



510 ,

M

₄₃



438 .

(18) Exemplary object operation process transient probabilities Monte Carlo evaluation

The simulation is performed according to data given in section 6.1. The first step is to select the initial operation state z_b(g), b1,2,3,4, at the moment t0, using formula (5), which is given by

                , 1 60 . 0 , 60 . 0 31 . 0 , 31 . 0 21 . 0 , 21 . 0 0 , ) ( 4 3 2 1 g z g z g z g z g z_b

where g is a randomly generated number from the uniform distribution on the interval

0 ,

1

. The next operation state z , _l l1,2,3,4, is generated according to (6) – (8), from z_bl(g), b1,2,3,4, defined as             , 1 54 . 0 , 54 . 0 22 . 0 , 22 . 0 0 , ) ( 4 3 2 1 g z g z g z g zl             , 1 5 . 0 , 5 . 0 2 . 0 , 2 . 0 0 , ) ( 4 3 1 2 g z g z g z g z l             , 1 28 . 0 , 28 . 0 12 . 0 , 12 . 0 0 , ) ( 4 2 1 3 g z g z g z g zl             . 1 7 . 0 , 7 . 0 48 . 0 , 48 . 0 0 , ) ( 3 2 1 4 g z g z g z g z l

For instance, if zb(g)z1, then the next operation state would be z , 2 z or 3 z 4 generated from z1l(g). Applying (9), the realizations of the empirical conditional sojourn times are generated according to the formulae

(8)

], 1 ln[ 0219 . 0 ) ( 12 t  H  31(t)0.0993ln[1H], ], 1 ln[ 0548 . 0 ) ( 13 t  H  32(t)0.0548ln[1H], ], 1 ln[ 0228 . 0 ) ( 14 t  H  34(t)0.0342ln[1H], ], 1 ln[ 0110 . 0 ) ( 21 t  H  41(t)0.0371ln[1H], ], 1 ln[ 0092 . 0 ) ( 23 t  H  42(t)0.0582ln[1H], ], 1 ln[ 0063 . 0 ) ( 24 t  H  43(t)0.0500ln[1H],

where H is a randomly generated number from the uniform distribution on the interval 0,1 .

The object operation process characteristics are calculated using the Monte Carlo method with time of the experiment fixed as ~18250 days.

Applying (12) the limit values of the object operation process transient probabilities at the operations states z_b are as follows:

, 233 . 0 1 p p20.037, p30.278, p40.452. (19) Based on the formula (11) and applying (19), the object operation process unconditional mean sojourn times b, b1,2,...,v, measured in days at the

particular operation states are given by 45

. 289 1

M , M₂68.71, M3 349.62, M4391.48. (20) Hence, applying (2) and according to (19), the object operation process expected values E[ˆb] of the total sojourn times ˆ , b b1,2,3,4, at the particular operation

states z , _b b1,2,3,4, and during the fixed operation time  1 year = 365 days are given by , 0 . 85 365 233 . 0 ] ˆ [₁    E E[ˆ₂]0.03736513.5, , 5 . 101 365 278 . 0 ] ˆ [3    E E[ˆ₄]0.452365165.0. (21) Exemplary object reliability

For the considered exemplary object, we assume that the conditional reliability functions are different at various operation states and have the exponential forms

] 002067 . 0 exp[ ) ( ) 1 ( t t   R at z , 1 () exp[ 0.00144] ) 2 ( t t   R at z , 2 ] 002611 . 0 exp[ ) ( ) 3 ( t t   R at z , 3 () exp[ 0.003939] ) 4 ( t t   R at z , (22) 4 for t0,).

(9)

Hence, the mean values of the object conditional lifetimes at the particular operation states measured in days are given by

87 . 483 1  days, 2694.44 days, 04 . 383 3  days, 4253.88 days.

The object operation time is large enough to apply the formula (4) and get 88 . 253 452 . 0 04 . 383 278 . 0 44 . 694 037 . 0 87 . 483 233 . 0          359.675days. (23)

Exemplary object reliability optimization

The object characteristics can be improved by changing the parameters of its operations process. According to (23), the objective function defined by (13) takes the form

  p1483.87 p2694.44 p3383.04p4253.88. (24)

where pb, b1,2,3,4, are the transient probabilities we want to optimize. Arbitrarily assumed, the bound constraints of the transient probabilities p b

respectively are: , 351 . 0 201 . 0 p1 0.03p20.105, , 395 . 0 245 . 0 p3 0.309 p40.459,



  4 1 . 1 b b p

The object conditional lifetime mean values _b, b1,2,3,4, arranged in non-increasing order are as follows ₁ ₃ ₄. Further, we substitute x1 p2,

, 1 2 p

x  x3p3, x4  p4,and we maximize with respect to x i, i1,2,3,4, the

linear form (24) that takes the form

 x₁694.44x₂483.87x₃383.04 x₄253.88,

with the following bound constraints , 105 . 0 03 . 0 1 1 1 x x x     x2 0.201 x2 0.351 x2,   _ _ _ _ , 395 . 0 245 . 0 3 3 3 x x x     x4 0.309 x4 0.459 x4,   _ _ _ _



  4 1 . 1 i i x

(10)

Therefore, we calculate [4]



   4 1 , 785 . 0 i i x x   x yˆ1 = 1 – 0.785 = 0.215 , 0 0 0_ _ x x  x1x10.075, x2x2 0.225, , 375 . 0 3 3_ _ x x  x4x4 0.525.

From the above, the inequality from [4] takes the form , 215 . 0   I I x x 

it follows that the largest value I{0,1,2,3,4} such that this inequality holds is .

1



I Therefore, we fix the optimal solution that maximizes linear function (24) and we get 105 . 0 1 1x  x  , 2 1 1 2 yˆ x x x x     0.2150.1050.030.2010.341, , 245 . 0 3 3x  x  x₄x₄0.309. (25) Finally, after making the inverse substitution, we get the optimal transient probabilities , 105 . 0 1 2x  p  p₁x₂0.341, p₃x₃0.245, p₄x₄0.309, (26) Finally, after substituting in (24) the above optimal transient probabilities that maximize the exemplary object mean lifetime  we get its optimal value

  p1483.87 p2694.44 p3383.04p4253.88

0.341483.870.105694.440.245383.040.309253.88



410.21 days, (27)

which is greater than that before optimization given by (23)(.

Substituting the optimal solution (26) into the formula (3) the corresponding optimal unconditional reliability function of the object is of the form

) (t

R 0.341R(1)(t)0.105R(2)(t)0.245R(3)(t)0.309R(4)(t),

(11)

Exemplary object operation process optimal characteristics

Having the values of the optimal transient probabilities determined by (26), it is possible to find the optimal unconditional mean values of the sojourn times of the object operation process at the operation states. Substituting the optimal transient probabilities at operation states

, 341 . 0 1 p p₂0.105, p30.245, p40.309,

determined by (26) and considering the steady probabilities determined in [4] , 236 . 0 1  2 0.169,



3 0.234, 40.361,

we get the following system of equations

0 123101 . 0 079794 . 0 057629 . 0 155524 . 0 ₁ ₂ ₃ ₄  M M M M 0 037905 . 0 02457 . 0 15126 . 0 02478 . 0 M1 M2 M3 M4 0 088265 . 0 17679 . 0 041321 . 0 057702 . 0 M₁ M₂ M₃ M₄ 0 24945 . 0 072306 . 0 052221 . 0 072924 . 0 M1 M2 M3 M4 (28)

with the unknown optimal mean values M . _b

Since the determinant of the main matrix of the system of equations (28) is equal to 0, then its rank is less than 4 and there are non-zero solutions of this system of equations that are ambiguous and dependent on one or more parameters. Thus, we may fix some of them and determine the remaining ones. In our case, according to (20), we conclude that it is sensible to assume M₄ 400. Consequently, from [4], the obtained optimal mean values of the object unconditional sojourn times at the operation states, are as follows

1

M  675, M ₂  290, M ₃  490, M₄ 400. (29) It can be seen that these solutions differ much from the values M , 1 M2, M and 3,

4

M given by (20).

Having these solutions, it is also possible to look for the optimal values M of the bl

mean values M of the conditional sojourn times at operation states. Namely, bl

substituting the probabilities of the object operation process transitions between the operation states, determined by (16) and the optimal mean values M given by b

(12)

12 22 . 0 M 0.32M₁₃0.46M₁₄675 21 20 . 0 M 0.30M₂₃ 0.50M₂₄290 31 12 . 0 M 0.16M₃₂ 0.72M₃₄ 490 41 48 . 0 M 0.22M₄₂0.30M₁₄400 with the unknown optimal values M we want to find. bl

As the solutions of the above system of equations are ambiguous, then we arbitrarily fix some of them because, for instance because of practically important reasons, and we find the remaining ones. In this case we proceed as follows:

- we fix M₁₂200, M₁₃500 and we find M₁₄1024; - we fix M₂₁100, M₂₃100 and we find M₂₄480; - we fix M₃₁900, M₃₂500 and we find M₃₄419;

- we fix M₄₁300, M₄₂500 and we find M₄₃487. (30) It can be seen that these solutions differ much from the mean values of the object conditional sojourn times at the particular operation states before its operation process optimization given by (18).

Another very useful and much easier to be applied in practice tool that can help in planning the operation process of an object are the system operation process optimal mean values of the total sojourn times at the particular operation states during the system operation time that by the assumpion is equal to  1year = 365 days. Under this assumption, after aplying (2), we get their optimal values

, 5 . 124 365 341 . 0 ] ˆ [1  p1    E  E[ˆ2]p20.10536538.3, , 4 . 89 365 245 . 0 ] ˆ [3 p3    E  E[ˆ4]p4 0.309365112.8, (31) that differ much from the values of E[ˆ1], E[ˆ2], E[ˆ3], E[ˆ4] determined by (21). In practice, the knowledge of the optimal values of M _b M and _bl E[ˆb] given respectively by (29)-(31), can be very important and helpful for an object operation process planning and improving in order to make the object operation more reliable and safer.

7. Comments on the object operation and reliability new strategy

The comparison of the selected object characteristics before the object operation process optimization given by (18), (20) and (21) with their values after the object operation process optimization respectively given by (29)-(31) justifies the sensibility of the performed object operation process optimization.

(13)

From the above it can be suggested to organize the object operation process in the way that causes the replacing approximately the conditional mean sojourn times

bl

M of the object at the particular operation states before the optimization given by

(18) by their optimal values M_bl after the optimization given by (30). However, the fulfilling this suggestion of the operation process parameters changing is not easy in practice.

It seems to be practically a bit easier way to change the object operation process characteristics by the reorganizing the operation process that results in replacing approximately the unconditional mean sojourn times _{M of the object at the}_b particular operation states before the optimization given by (20) by their optimal values M_b after the optimization given by (29).

The easiest way of the object operation process reorganizing is that leading to the replacing approximately the total sojourn times E[ˆ_b] before the optimization given by (21) by their optimal values E[ˆ_b] after the optimization given by (31).

8. Conclusions

The achieved results may be considered as an illustration of the possibilities of applications of the proposed methods and procedures to the system operation and reliability optimization. The obtained evaluation may be a very useful example in complex systems reliability characteristics and prediction, especially during the design and when planning and improving their operation processes safety and effectiveness.

9. References

[1] Grabski F., Jaźwiński J.: Funkcje o losowych argumentach w zagadnieniach niezawodności, bezpieczeństwa i logistyki. Wydawnictwa Komunikacji i Łączności, Warszawa, 2009.

[2] Grabski F.: Semi-Markowskie modele niezawodności i eksploatacji. Instytut Badań Systemowych PAN, Warszawa, 2002.

[3] Kołowrocki K. and Soszyńska J.: Reliability modeling of a port oil transportation system’s operation processes. International Journal of Performability Engineering, Vol. 6, No 1, 2010, s. 77-87.

[4] Kołowrocki K. and Soszyńska-Budny J.: Reliability and Safety of Complex Systems and Processes: Modeling – Identification – Prediction – Optimization. Springer, 2011.

[5] Kuligowska E.: Preliminary Monte Carlo approach to complex system reliability analysis. Summer Safety and Reliability Seminars – SSARS 2012, Journal of Polish Safety and Reliability Association, Vol.3,No1,2012,s.59-71.

(14)

[6] Soszyńska J.: Systems reliability analysis in variable operation conditions. International Journal of Reliability, Quality and Safety Engineering. Special Issue: System Reliability and Safety, Vol. 14, No 6, 2007, s. 617-634.

[7] Soszyńska J.: Reliability and risk evaluation of a port oil pipeline transportation system in variable operation conditions. International Journal of Pressure Vessels and Piping Vol. 87, No 2-3, 2010, s. 81-87.

[8] Zio E. and Marseguerra M.: Basics of the Monte Carlo Method with Application to System Reliability. LiLoLe, 2002.

Krzysztof Kołowrocki is a Full Professor and the Head of Mathematics

Department at the Faculty of Navigation of Gdynia Maritime University. His field of interest is mathematical modeling of safety and reliability of complex systems and processes. He has published several books and over 300 scientific articles and papers. He is the President of Polish Safety and Reliability Association. His home site can be found at: http://www.am.gdynia.pl/~katmatkk/

Ewa Kuligowska MSc Eng., is an Assistant at Department of

Mathematics of the Faculty of Navigation of Gdynia Maritime University. Her field of interest is Monte Carlo simulation analysis of complex systems reliability. She has published several publications in scientific journals and conference proceedings.

Joanna Soszyńska is an Assistant Professor at Department of

Mathematics of the Faculty of Navigation of Gdynia Maritime University. Her field of interest is mathematical modelling of safety and reliability of complex systems at variable operation conditions. She has published two books and over 100 papers in scientific journals and conference proceedings.