Strojniški vestnik - Journal of Mechanical Engineering 53(2007)10, 621-634 UDK - UDC 621.9.019 Izvirni znanstveni članek -Original scientific paper (1.01) Napoved zbirnega števila okvar popravljivega izdelka na podlagi poteka delovanja Prediction of the Cumulative Number of Failures for a Repairable System Based on Past Performance Boštjan Veber - Marko Nagode - Matija Fajdiga (Fakulteta za strojništvo, Ljubljana) Napoved zbirnega števila okvar popravljivega izdelka je pomembna tema v teoriji zanesljivosti. Popravljiv izdelek je lahko po opravljenem popravilu v treh možnih stanjih: 'dober kot nov', 'slab kot star' in 'boljši kot star toda slabši kot nov'. Običajni verjetnostni modeli za napovedovanje pričakovanega števila okvar temeljijo na prvih dveh stanjih, zato so neprimerni za opis zadnjega stanja, ki je bolj pogosto v praksi. V prispevku je predstavljena robustna rešitev verjetnostnega modela, t.i. splošni obnovitveni proces (SOP), ki je omogoča predstavitev vseh treh stanj izdelka po popravilu. Raziskave kažejo, da je s SOP na osnovi mešanice m-tih Weibullovih porazdelitev omogoča splošen pristop k opisu kompleksnih popravljivih izdelkov in podaja uporabo algoritmom matematičnega pričakovanja (EM) za oceno neznanih parametrov. V prispevku je predstavljen tudi standardni SOP na podlagi dvoparametrične Weibullove porazdelitve. SOP na osnovi mešane porazdelitve z m komponentami in standardni SOP primerjamo z izračunom pričakovanega števila okvar in funkcije napake. Na podlagi rezultatov lahko sklepamo, da predlagan SOP na osnovi mešanice Weibullovih porazdelitev točno opisuje izmerjene vrednosti okvar in je primeren za napovedovanje okvar na podlagi predhodnega poteka delovanja izdelka, kljub omejenemu številu razpoložljivih podatkov o okvarah. © 2007 Strojniški vestnik. Vse pravice pridržane. (Ključne besede: napovedi okvar, popravljivi izdelki, numerično modeliranje, ocenjevanje parametrov) The prediction of the cumulative number of failures for a repairable system is an important topic in reliability theory. A repairable system may end up in one of the three possible states after a repair: 'as good as new', 'as bad as old' and 'better than old but worse than new'. Current probabilistic models used in repairable system analysis account for the first two states, but they do not properly apply to the last one, which is more common in practice. In this paper, a robust solution to a probabilistic model that is applicable to all of the three after repair states, called generalized renewal process (GRP), is presented. This research demonstrates that the GRP based on an m-fold Weibull mixture offers a general approach to modeling complex repairable systems and discusses application of the EM algorithm to estimation of the GRP parameters. This paper also presents a review of the standard GRP based on two-parameter Weibull distribution. The GRP with m mixture components distributions is compared to the standard GRP by calculating the expected cumulative number of failures and the error function. It is shown that the proposed GRP solution with a Weibull mixture accurately describes the failure data and it is suitable for predicting failures based on the past performance of the system, even when a small amount of failure data is available. © 2007 Journal of Mechanical Engineering. All rights reserved. (Keywords: failure prediction, repairable systems, numerical modeling, parameter estimations) 0 UVOD 0 INTRODUCTION Napoved zbirnega števila okvar The prediction of the cumulative number of popravljivega izdelka je pomembna tema v teoriji failures of a repairable system is an important topic 621 Strojniški vestnik - Journal of Mechanical Engineering 53(2007)10, 621-634 zanesljivosti. Omogoča načrtovanje preventivnih vzdrževalnih posegov, stroškov kroga trajanja, ocenjevanje razpoložljivosti izdelka itn. Popravljiv izdelek je lahko v enem izmed treh mogočih stanjih po opravljenem popravilu: ‘dober kakor nov’, ‘slab kakor star’ in ‘boljši kakor star toda slabši kakor nov’. Trenutni verjetnostni modeli za analizo popravljivih izdelkov, kakršna sta obnovitveni proces (OP) in nehomogeni Poissonov proces (NHPP), so primerni za opis le prvih dveh stanj. Za prikaz zadnjega stanja izdelka po popravilu ni učinkovitega in dovolj natančnega postopka. Kijima in Sumita [1] sta predlagala nov verjetnostni model, imenovan splošni obnovitveni proces (SOP), ki omogoča predstavitev vseh mogočih stanj izdelka po popravilu. Postopek s SOP sta izpopolnila Kaminskiy in Krivtsov [2] in predlagala približno rešitev na temelju metode Monte Carlo (MC). Predpostavka, ki omogoča simulacijo SOP z metodo Monte Carlo, je poznavanje porazdelitve verjetnosti časa do prve okvare (ČDPO) in kakovosti popravila q oziroma možnost njune ocene na podlagi razpoložljivih podatkov ([2] in [3]). Čas popravila zanemarimo, tako lahko zaporedje okvar opazovanega izdelka obravnavamo kot naključni točkovni postopek. V glavnem je bil postopek Monte Carlo razvit za uporabo v primerih z veliko množico podatkov. Razpoložljivost večjega števila podatkov omogoča oceno porazdelitve verjetnosti ČDPO in q z visoko stopnjo natančnosti [4]. Seveda je težko zagotoviti zadostno količino podatkov v primerih, ko je na voljo omejeno število enakih izdelkov. Postopek s simulacijo Monte Carlo ima pri SOP nekaj prednosti in pomanjkljivosti. Postopek omogoča rešitev za vse vrste porazdelitev, vključno z empirično, ki je nepristranska in dosledna. Na drugi strani postopek MC potrebuje veliko množico podatkov in je časovno zelo potraten. Ti razlogi omejujejo uporabo postopka Monte Carlo Kaminskega in Krivtsova zunaj avtomobilske industrije. Krivtsov [3] se je zavedal zapletenosti in težavnosti razvoja matematično prilagodljivega verjetnostnega modela za SOP, zato je predlagal drugo pot na podlagi metode največje verjetnosti (MNV), vendar brez ustrezne rešitve za SOP Kasneje je Yanez [4] na podlagi MNV izpeljal rešitev za oceno parametrov SOP Razvoj MNV za in reliability theory. It enables the planning of preventive maintenance actions and costs, an estimation of the system’s availability, etc. A repairable system may end up in one of three possible states after a repair: ‘as good as new’, ‘as bad as old’ and ‘better than old, but worse than new’. Current probabilistic models used in repairable system analysis, such as the renewal process (RP) and the non-homogeneous Poisson process (NHPP), account for the first two states, respectively. However, no practical and accurate approach exists to address the remaining after-repair state. A new probabilistic model to address all the after-repair states called the ‘generalized renewal process’ (GRP) has been proposed by Kijima and Sumita [1]. This GRP approach has been extended by Kaminskiy and Krivtsov [2] and they have offered a Monte Carlo (MC) based approximate solution. The assumption that makes Monte Carlo simulation of the GRP possible is that the time to first failure (TTFF) distribution and the quality of the repair q are known and can be estimated from the available data ([2] and [3]). Furthermore, the repair time is assumed to be negligible so that the failures can be viewed as point processes. The Monte Carlo approach was developed mainly for cases where a large set of data is available. The availability of such data allows for an estimation of the TTFF distribution and q with a high degree of accuracy [4]. However, it would be difficult to obtain the same amount of data for cases where only a limited number of identical systems are present. There are some advantages and disadvantages of using the Monte Carlo approach to the GRP. The approach offers a solution for all kinds of distributions, including empirical ones, which is unbiased and consistent. In contrast, besides the need for large amounts of data, the approach is extremely time consuming. For these reasons the direct application of Kaminskiy and Krivtsov’s Monte Carlo approach outside the automotive industry would be limited. Krivtsov [3] recognized the complexities and the difficulties of developing a mathematically tractable probabilistic model to the GRP, and discussed an alternative maximum-likelihood (ML) estimation approach to solve the GRP without offering any solution. Later, Yanez [4] developed a solution based on ML to estimate the GRP parameters. The 622 Veber B. - Nagode M. - Fajdiga M. Strojniški vestnik - Journal of Mechanical Engineering 53(2007)10, 621-634 oceno parametrov SOP je odpravil potrebo po večji količini podatkov za izvedbo analize okvar (kakor v primeru postopka MC Kaminskega in Krivtsova pri določevanju SOP). Yanezov postopek z MNV je izpeljan s predpostavko, da porazdelitev ČDPO ustreza dvoparametrični Weibullovi porazdelitvi verjetnosti ter da naslednji časi med okvarami (ČMO) sledijo pogojni dvo-parametrični Weibullovi porazdelitvi. Dvo- ali triparametrična Weibullova porazdelitev je najbolj splošno uporabna porazdelitev na področju modeliranja zanesljivosti, ker s parametroma ß (oblika) in 0 (velikost) omogoča zelo različne oblike. Zato je primerna za modeliranje raznolikih podatkov in obratovalnih značilk izdelka, t.i. porazdelitev ČDPO. Oblika porazdelitve funkcije dobe trajanja izdelka je pogosto sestavljena iz več osnovnih oblik, zato je smiselna vpeljava mešane porazdelitve kot osnovne porazdelitev SOP Značilna pomanjkljivost, skupna vsem mešanim porazdelitvam, je težavna ocena neznanih parametrov. Dvoparametrična Weibullova porazdelitev velja za zelo pomembno v teoriji zanesljivosti, zato je primerna njena uporaba kot komponente mešane porazdelitvene ČDPO okvare SOP Veliko število prispevkov je na temo uporabnosti mešanice Weibullovih porazdelitev, posebej pri modeliranju zanesljivosti izdelka, obratovalni trdnosti in analizi preživetja. Predlagane so bile številne tehnike ocenjevanja neznanih parametrov ([5] do [9]). Namen prispevka je ugotavljanje primernosti končne mešanice Weibullovih porazdelitev s pozitivnimi utežmi komponent kot osnovne porazdelitve ČDPO za SOP, kljub temu da ‘dejanska’ porazdelitev ni mešana porazdelitev ali prave uteži niso negativne. Dovolj so le podatki o časih okvar posameznih, če neznane parametre SOP ocenimo z ustrezno metodo. V ta namen je izpeljan postopek EM za SOP z osnovno mešanico m Weibullovih porazdelitev. Prispevek ima naslednjo sestavo. SOP in funkcija največje verjetnosti sta podana v prvem poglavju. Drugo poglavje predstavlja model na temelju m-kratne mešanice Weibullovih porazdelitev in metodo za oceno parametrov na podlagi algoritma EM. Sledi tretje poglavje s prikazom uporabe in primerjave na primerih za predlagani model na podlagi mešane Weibullove porazdelitvene funkcije in osnovni triparametrični SOP na podlagi ocene parametrov z MNV. Četrto development of the ML parameter estimation removes the need to have a large set of failure data in order to perform the analysis (as is the case with Kaminskiy and Krivtsov’s Monte Carlo simulation approach to the GRP). Yanez’s ML approach was derived using the assumption that the TTFF distribution is a two-parameter Weibull distribution and that the subsequent times between failures (TBFs) follow the conditional two-parameter Weibull distribution. The two- and three-parameter Weibull distributions are some of the most commonly used distributions in reliability engineering because of the many shapes they attain for various values of the parameters ß (shape) and ? (scale). It can, therefore, model a wide variety of data and life characteristics, i.e., the TTFF distribution. Since the shape of the life distribution is often composed of more than one basic shape, a reasonable step is to introduce the mixture distribution as the genuine distribution for the GRP. A significant difficulty common to all mixed distributions during their application is the estimation of unknown parameters. Since the Weibull distribution is considered very important in reliability studies, it is natural to consider it as a base distribution in finite mixtures for the distribution of the TTFF in the GRP. There are a number of papers dealing with the usefulness of m-fold Weibull mixture distributions, especially in reliability, fatigue and survival analysis. In addition, several estimation techniques for unknown parameters have been proposed ([5] to [9]). The aim of this article is to prove that a finite Weibull mixture, with positive component weights only, can be used as the underlying distribution of the TTFF of the GRP in spite of the ‘true’ distribution function of the TTFF not being a mixture distribution or the true weights being negative. In this way it is sufficient to have only the observed failure data set of the system, if the unknown parameters of the GRP are estimated by a proper method. For this purpose, the EM procedure is derived for the GRP with an underlying m-fold Weibull mixture. The paper is structured as follows. The GRP and the likelihood function are given in Section 1. Section 2 presents a model based on the m-fold Weibull mixture and an estimation method based on the EM algorithm. Following this, two failure data sets are used in Section 3 to illustrate the application of the proposed GRP based on the Weibull mixture and to compare it with the standard three-parameter GRP based on the ML estimation. Napoved zbirnega števila okvar - Prediction of the Cumulative Number of Failures 623 Strojniški vestnik - Journal of Mechanical Engineering 53(2007)10, 621-634 poglavje sklenemo z razpravo ter s komentarji rezultatov. Predpostavimo, da popravljiv izdelek prične delovati pri času t=0. Zaporedje časov delovanj {X;}™ sestavlja naključni točkovni obnovitveni postopek. Zaradi lažje nadaljnje izpeljave predpostavimo, da je pravilo takojšnje. Kijima [10] je vpeljal pojem navideznega časa delovanja. Če ima izdelek takoj po opravljenem (w-l)-tem popravilu navidezno dobo delovanja V =y, potem ima dolžina «-tega kroga oziroma čas med okvarama (ČMO) X naslednjo zbirno porazdelitveno funkcijo (ZPF), t.i. osnovna porazdelitev: 1.1 Funkcija zanesljivosti vzorca Če izvedemo oceno parametrov SOP za trenutek, ko se pojavi naslednja okvara, potem je We conclude with a discussion and some comments in Section 4. Consider a repairable system that starts functioning at time t=0. A sequence of operating times {Xj}™ forms a renewal-type stochastic point process. For simplicity, the repair is assumed to be instantaneous. Kijima [10] introduced the notion of virtual age. If a system has a virtual age of Vn =y immediately after the («-l)th repair, then the duration of the «th cycle or the time between failures (TBF) X has the following cumulative distribution function (CDF) (the so-called underlying distribution): 1.1 Likelihood function If an estimation of the GRP is made for the moment when the very next failure occurs, then the 1 OSNOVNA DOLOČITEV 1 BASIC DEFINITION l-F(v) R(x + y) (1), F(y) R(y) kjer je F(t) ZPF časa do prve okvare (ČDPO) where F(t) is a CDF of the time to first failure (TTFF) novega izdelka in R(t)=l-F(t) je zanesljivost pri of a new system and R(t)=l-F(t) is the reliability at posameznih časih. Vsoto: the respective time. The summation: r»=X> (2), !=1 s T=Q imenujemo dejanska doba delovanja izdelka. Predpostavimo, da «-to popravilo odpravi poškodbe, ki so se nakopičile pri obratovanju med (w-l)-to in «-to okvaro. S to predpostavko zapišemo navidezno dobo delovanja izdelka po «-tem popravilu z: with T0=0, is called the real age of the system. It is assumed that the nth repair would only compensate for the damage accumulated during the time between the (n-1)th and the nth failure. With this assumption the virtual age of the system after the nth repair is: Vn=Vn_i+q-Xn, « = 1,2,... (3), kjer je q parameter učinkovitosti popravila (ali where q is the repair-effectiveness parameter (or pomladitveni parameter) in V=0. Če opazujemo the rejuvenation parameter) and V=0. If the times čase med zaporednimi okvarami, potem lahko between the successive failures are considered, then navidezno dobo delovanja izrazimo kot [3]: the virtual age can be expressed as [3]: n Vn=q-Y,Xi=fTn Glede na ta model in sliko 1, privzeta vrednost q=0 vodi do OP (‘dober kakor nov’), medtem ko predpostavka q=l vodi do NHPP (‘slab kakor star’). Vrednosti parametra q iz obdobja 0<9<1 predstavljajo stanje izdelka po popravilu, ko je ta ‘boljši kakor star, toda slabši kakor nov’. i=\ (4). According to this model and Fig. 1, the result of assuming a value of q=0 leads to a RP (as good as new), while the assumption of q=1 leads to a NHPP (as bad as old). The values of q that fall in the interval 0----------o------> v3 v2 Vi v3 v2 Vi Navidezni čas Rrtua/ tirne q =\ Dejanski čas Real tirne 0 h h h o h h h o h h h Sl. 1. Primerjava navidezne in dejanske dobe delovanja za različne vrednosti parametra q Fig. 1. Virtual age versus real age for varying q values funkcija zanesljivosti vzorca definirana z [11 ] : likelihood function is defined by [11 ] : n L = f(xl\ß,ß)Ylg(xi\ß,ß,q) (5), ! = 2 kjer sta f(x) in g(x) gostoti porazdelitve verjetnosti (GPV) za ČDPO in pogojna GPV za ČMO. Takoj je pomembno poudariti, da se prva okvara ne pokorava prej omenjeni pogojni verjetnosti. Torej lahko negativni logaritem funkcije zanesljivosti vzorca, imenovan tudi funkcija napake, podamo z izrazom: where f(x) and g(x) are the probability density function (PDF) of the TTFF and the conditional PDF of the TBF. At this point it is important to recall that the first failure does not obey the conditionality mentioned previously. Therefore, the negative logarithm of the likelihood function, which can be regarded as an error function, is given as: L = -log(I) = - log{/(xi |/?,0)} +J>g{g(x, \ß,e,q)} (6). 1.2 Triparametrični SOP 1.2 Three-parameter GRP Običajni SOP je utemeljen na predpostavki, da je osnovna porazdelitev ČMO pogojna dvoparametrična Weibullova porazdelitvena funkcija [4]. Kakor je bilo predhodno omenjeno, se čas do prve okvare ne podreja pogojnosti ČMO. Dvoparametrična Weibullova ZPF za ČDPO ima obliko: The standard GRP is based upon the assumption that the underlying TBF distribution is a conditional two-parameter Weibull probability distribution [4]. As we mentioned before, the time to first failure does not comply with the conditionality of the TBF. The two-parameter Weibull CDF for the TTFF is of the form: F(*,|/?,0) = l-exp ß e (7), kjer sta 0 in ß parametra velikosti in oblike. Za naslednje okvare je enačba (1) podana kot: G(xi|/?,0,?) = l-exp kjer sta x čas med (i-1)-to in i-to okvaro in t je zbirni čas delovanja (dejanska doba delovanja) do (i-1)-te okvare (sl. 2): where ? in ß are the scale and shape parameters, respectively. For the subsequent failures Eq. (1) becomes: q-ti- fi e q-h- e ß (8), where xi is the period of time between the (i-1)th and the ith failure and ti- 1 is the cumulative operating time (i.e., the actual age) up to the (i-1)th failure (see Fig. 2): Napoved zbirnega števila okvar - Prediction of the Cumulative Number of Failures 625 Strojniški vestnik - Journal of Mechanical Engineering 53(2007)10, 621-634 Xi o h i x2 2 A X'-l A X' A -A------AA> 1-2 !-l i t U-l SI. 2. Diagram časov do posameznih okvar Fig. 2. Time-to-failure diagram of the failure events Z odvajanjem obeh ZPF in z nekaj algebre lahko odvod po velikostnem parametru ? logaritma By differentiating the CDF functions and with some algebra the derivatives of the log- funkcije zanesljivosti vzorca (enačba (6)) zapišemo likelihood function (Eq. (6)) with respect to the v samostojni izraz samo z oblikovnim parametrom scale parameter ? can be expressed in terms of the ß in s parametrom učinkovitosti popravila q. S tem shape parameter ß and the repair-effectiveness se bistveno poenostavi reševanje. Sistem dveh parameter q in a closed-form expression, thus enačb za oceno parametrov ß in q je rešujemo sočasno in je prikazan spodaj: simplifying the solution. The system of equations for the parameters ß and q should be solved simultaneously and is shown below: -1 log(x! ) + ^ log(q ¦ ti_x + Xi ) ---------------- ! = 2 xf log(x!) + ^[(q-h_x +xif-log(0 and the parameters 0=(w ,...,w ,0 ,...,0 ) are such za (/-l,...,m) in S™=1w, = 1. Stalnico w.imenujemo utežni količnik in/ je komponenta mešanice gostot porazdelitev verjetiiosti ter je določena s parametri 0. Ker zavzema Weibullova porazdelitev pomembno mesto v teoriji zanesljivosti, je pri našem postopku zvezna gostota porazdelitve verjetnosti (GPV) časa do prve okvare modelirana kot mešanica dvoparametričnih Weibullovih porazdelitev: that wj>0 for (j=1,...,m) and constant wj is called a weighting factor and fj a component density function parameterized by ?j. In our approach a continuous probability density function (PDF) of times to the first failure is modelled as a mixture of two-parameter Weibull distributions, because the Weibull distribution is considered to be very important in reliability studies: f (t1 m |0) = J> 7=1 bj j q kjer sta ß in 6 oblikovni in velikostni parameter komponente mešanice Weibullovih porazdelitev. Vrednosti parametrov so omejene glede na lastnosti Weibullove porazdelitve, tako da velja /?>0 in 0>O za (j=\,...,m). Če imamo mešanico Weibullovih porazdelitev za osnovno GPV za ČDPO, potem je pripadajoča funkcija zanesljivosti določena kot utežna vsota posameznih komponent funkcije zanesljivosti [12]: ßj-i /3,1 h h expj - { J j { J j (13), where the constants ßj in ?j stand for the Weibull shape and scale parameters of the component densities. They are limited according to the characteristics of the Weibull distribution ßj>0 and ?j>0 for (j=1,...,m). If the Weibull mixture is considered as underlying the PDF of the TTFF, the corresponding reliability function can be defined as a weighted sum of component reliability functions [12]: R(t\@) = Y]wjRj(t\Qj) (14). 7=1 Z upoštevanjem zgornjega izraza za funkcijo Considering the expression for the reliability zanesljivosti lahko pogojno ZPF za čas med prvo function, the conditional CDF of the times between in drugo okvaro, za čas med drugo in tretjo okvaro the first and second failure, the second and third itn., zapišemo ob upoštevanju enačbe (1) kot: failure, etc. can be written according to Eq. (1) as: G(Xi \0,q) = Pr {Xt < Xi | Vt_x = qtt_x} = 1 L" H^ta-i+*iM> L>^-iiM> (15). m Napoved zbirnega števila okvar - Prediction of the Cumulative Number of Failures 627 Strojniški vestnik - Journal of Mechanical Engineering 53(2007)10, 621-634 2.2 Algoritem EM 2.2 EM algorithm Algoritem EM je splošna metoda za iskanje parametrov porazdelitvene funkcije, ki dajo največjo vrednost funkcije zanesljivosti vzorca glede na izmerjene vrednosti. Standardni algoritem EM obravnava izmerjene podatke kot nepopolno informacijo o opazovanem naključnem pojavu. Pri algoritmu EM vstavimo izmerjene podatke v večji ‘popoln’ vzorčni prostor. Pripadajoča funkcija zanesljivosti vzorca zagotavlja celoten opis generacije oziroma vzorčenja podatkov glede na iskane parametre. Večji vzorčni prostor seveda ne moremo vzorčiti ali trenutno ni na voljo. Torej računamo pričakovano vrednost logaritma funkcije zanesljivosti vzorca celovitih podatkov, odvisno od izmerjenih vrednosti {x.;i=(1,...,n)} in trenutno oceno parametrov 0(K). Pri SOP na podlagi mešanice Weibullovih porazdelitev je negativni logaritem funkcije zanesljivosti vzorca za analizo do zadnje zaznane okvare podan z: logOL) The EM algorithm is a general method for finding the maximum likelihood estimate of the parameters of an underlying distribution from a given data set [13]. The standard EM algorithm regards the measured data as incomplete information about the underlying stochastic process. For the EM algorithm, one embeds the measured data in a larger ‘complete’ data space. The corresponding likelihood function provides a complete description of the data generation, given the parameters in question. The larger data set, however, cannot be sampled or is not available at the moment. One computes, therefore, the expected value of the log-likelihood function of the complete data conditioned by the observed data {xi;i=(1,...,n)} and the current parameter estimates ?(k). For the GRP based on the Weibull mixture distribution the failure-terminated negative logarithm of the likelihood function is defined by: log {/(X! | 0)} + J] log {^x; | ©,C/1 ** ' ©(i), ék) ) • iog(w, ) 7=1 1=2 7=1 (17), nt rt hi Y, fU I *i,©(i)) • log { fj (*i I e,.) } + Y, Y,sU I xi>0(i),?(i) ) • !°g { gj (xi I Qj,ii) } 7=1 !=2 7=1 kjer smo vpeljali obrobne verjetnosti f(j|x1,®(k)) in g(J|X1, e"),?")), ki jih izpeljemo z uporabo Bayesovega pravila [14]. Vrednotenje zgornjega izraza imenujemo korak E algoritma. Pri drugem koraku (korak M) algoritma EM maksimiziramo pričakovanje, ki smo ga izračunali v prvem koraku, glede na parametre, tako da dobimo novo oceno parametrov 0(i+1). Maksimiziranje enačbe (17) lahko izvedemo z neodvisno maksimizacijo člena, ki vsebuje w., in člena, ki vsebuje 6, saj nista medsebojno povezana. Z nekaj algebre dobimo naslednje izraze: where we introduce the posterior probabilities f(j|x1,?(k)) and g(j|x1,?(k),q(k)), which can be expressed using Bayes’s rule (see Ref. [14]). The evaluation of this expectation is called the E-step of the algorithm. The second step (the M-step) of the EM algorithm is to maximize the expectation we computed in the first step with respect to the parameters to obtain new parameter estimations ?(k+1). To maximize Eq. (17), we can maximize the term containing wj and the term containing ?j independently since they are not related. After some algebra we get the next expressions: ° = EE*t/l*,e(iVi)) ;=2 7=1 iß (*+l) IX- ß (*+i) (k+1) q ti -1 +xi e(k+i) wx(,u)^(?(i+Vi lor1»,?^) (19) k +1) (20), + 1 (k+1) bj kjer (*+!) RM h-\ 18) where expL[fo<*+H-i)/^+1f } (i+l) (*+!) (21). E^-^f^-U)/^} Dobili smo sistem (2m+l) enačb z (2m+l) neznanimi parametri. Na podlagi Newton-Raphsonove metode [15] smo razvili numerično metodo reševanja sistema zahtevnih nelinearnih enačb. Oceno novih parametrov 0(i+1), q{k+r> dobimo s posodobitvijo enačb (18), (19) in (20) z vrednostmi parametrov iz prejšnjega koraka 0(i), g«. Koraka pričakovanja in največje vrednosti se izvajata simultano. Algoritem napreduje z uporabo novih parametrov kot začetnim približkom naslednje iteracije. Algoritem EM ponavlja korake, dokler se ne približa rešitvi. 2.3 Ocena pričakovanega števila okvar Dogodek okvare pri popravljivem izdelku lahko simuliramo z metodo Monte Carlo, pri kateri z uporabo zbirne porazdelitvene funkcije za ČDPO vnesemo prvo okvaro, medtem ko za naslednje okvare uporabimo pogojno zbirno porazdelitveno funkcijo za ČMO. Predpostavimo, da je verjetnost pojava prve okvare F{xx) = P in naslednjih okvar G(x) = P naključno število iz enakomerne porazdelitve na območju [0,1], potem lahko zapišemo enačbi za simuliranje časov okvar: We obtain a system of (2m+1) equations with (2m+1) unknown variables. We have developed a solution based on the Newton-Raphson method to solve these extremely complex equations [15]. Update Eq. (18), (19) and (20) for the estimation of the new parameters ?(k+1), q(k+1), in terms of the old parameters ?(k), q(k) perform both the expectation step and the maximization step simultaneously. The algorithm proceeds by using the newly derived parameters as the guess for the next iteration. The EM steps are iterated until the algorithm converges. 2.3 Estimation of the expected number of failures The occurrence of failure in repairable systems can be simulated with a Monte Carlo method using the TTFF cumulative distribution function to generate the first failure, while for the subsequent failures, a conditional TBF cumulative distribution is used. Assume that the probability of the first failure F(x1) = P and for subsequent failure times G(xi) = P are random numbers from a uniform distribution on the interval [0,1], then we get the equations for the generation of failure times: !=2 !=2 o + =2 =2 1=2 Napoved zbirnega števila okvar - Prediction of the Cumulative Number of Failures 629 Strojniški vestnik - Journal of Mechanical Engineering 53(2007)10, 621-634 in 1- E 7=1 and = p E> ¦RM-h-x +*i) E (22) (23). 7=1 ¦RM-h-i) Predpostavimo obdobje opazovanja [0,7], za katerega želimo oceniti pričakovano število okvar. Na začetku simulacije MC določimo enakomerno porazdelitev za nastajanje naključnih števil v mejah [0,1] za P. Potem uporabimo numerično metodo, kakršni sta Newton-Raphsonova ali sekantna metoda, za rešitev enačb (22) in (23). Naključne vrednosti {x.;Hl,2,...)} določimo, kakor je prikazano na sliki 3. Naključno vrednost x prištejemo vsoti predhodnih naključnih vrednosti ter primerjamo z dobo opazovanja T. Postopek ponavljamo, dokler vsota naključnih vrednosti x. ne preseže T. Vzemimo, da smo ta postopek ponovili k-krat, pri tem je « število zapisanih ali simuliranih okvar pri j-ti ponovitvi. Pričakovano število okvar (ti. zbirna funkcija intenzivnosti CIF) pri času T je podana z: E[N(T)] = Assume a time period [0,T] as the period of interest to estimate the expected number of failures. The MC simulation starts with the definition of a uniform distribution that will generate random numbers in the range [0,1] for P. Then, a numerical method, such as the Newton-Raphson or the secant one, is used to solve Eqs. (22) and (23). Random values of {xi;i=(1,2,...)} are generated, as depicted in Fig. 3. A random value for xi is then added to the sum of the past random xis and compared with the period of interest T. This procedure is repeated until the sum of xis does not exceed T. Assume that this procedure is repeated k times and nj is the number of failures registered or generated at the jth repetition. The expected number of failures (i.e., the cumulative intensity function (CIF)) at time T is given by: E ;=1 (24). 3 NUMERIČNI PRIMERI 3 NUMERICAL EXAMPLES Na primeru iz literature in na simuliranih podatkih sta uporabljena dva SOP, običajni triparametrični model in model na temelju mešanice WeibuUovih porazdelitev. Grafični prikaz rezultatov omogoča bralcu lažje razumevanje predstavljenih pojmov, razlag rezultatov in sklepov. Za oba predstavljena SOP je izračunana vrednost funkcije napake E [14]. Two GRP, a standard three-parameter model and a model based on the Weibull mixture distribution, are applied in two examples, one from the literature and one that is simulated. The use of a graphical display enables the reader to gain a perspective on the various meanings and associated interpretations. In addition, the error function E [14] is calculated for both GRPs. SI. 3. Simulacija MC naključnih časov med okvarami za SOP Fig. 3. MC simulation of the random times between failures for the GRP m m k k 630 Veber B. - Nagode M. - Fajdiga M. Strojniški vestnik - Journal of Mechanical Engineering 53(2007)10, 621-634 Preglednica 1. Ocenjeni parametri SOP za posamezne primere Table 1. Estimated parameters for the GRP in the numerical examples Preglednica 2. Vrednosti E in AIC za SOP uporabljene v primerih Table 2. E and AIC values for the GRP used in the numerical examples Ocenjeni parametri / Estimated parameters Primer 2 SOP model GRP model Primer 1 Primer 2 SOP model Primer 1 Example No. 1 Example No. GRP model Example No. 1 Example No. 2 Standardni SOP E 460,814 31,938 « = 71 « = 33 Standard GRP Standardni SOP A = 3,12, 0o = 3649 A = 2,70, 0„ = 3,920 Standard GRP q =0,409 9 = 0,160 SOP z mešano porazdelitvijo GÄ73 based on mixture L 458,471 31,444 AE -2,343 -0,494 SOP z mešano porazdelitvijo W] = 0,923, 0] = 2664 W] = 0,503,0] = 3,001 (m = 2) GRP based on mixture ßi = 4,26 A = 2,62 (m = 2) w2 = 0,077, 02 = 5006 w2 = 0,497, 02 = 4,999 /i = 4,26 /4 = 3,77 q = 0,428