Symulacyjna metoda wyznaczania parametrów niezawodnościowych systemu przesyłowego o zależnych elementach

Instytut Bada Systemowych PAN

Wysza szkoła Informatyki Stosowanej i Zarzdzania

Streszczenie

Rozwaany jest ogólny model systemu przesyłowego (danych, sygnałów, energii elektrycznej, surowców, itp.) zbudowanego z wzłów i łczy, w którym wzły ko
cowe pełni rol ródeł lub uj. Komponenty systemu ulegaj losowym uszkodzeniom, których usuwaniem zajmuje si pewna liczba ekip naprawczych. Zakłada si, e czasy do uszkodzenia poszczególnych komponentów maj rozkłady wykładnicze, przy czym dla komponentu połczonego z wzłem ródłowym czas ten ma inny rozkład, ni dla tego samego komponentu odłczonego od wzła ródłowego (std zaleno). Zakłada si te, e czasy do naprawy poszczególnych komponentów maj rozkłady inne od wykładniczych. Złoono modelu powoduje, e parametry niezawodnociowe procesu uszkodze
i napraw rozwaanego systemu s praktycznie niemoliwe do wyznaczenia w sposób analityczny, zatem do ich oszacowania uyto symulacji Monte Carlo oraz estymacji statystycznej. Okazało si przy tym, e nietrywialnym zadaniem jest wyznaczenie przedziału ufnoci dla szacowanych parametrów. W referacie przedstawiono lematy formułujce warunki powracalnoci procesu uszkodze
i napraw, lematy okrelajce granice szukanych przedziałów ufnoci, oraz zaprezentowano przykłady obliczeniowe.

Słowa kluczowe: system przesyłowy, zaleno stochastyczna, parametry niezawodnociowe, polityka napraw, symulacja typu Monte Carlo, estymacja statystyczna.

1. Wprowadzenie

W niniejszym opracowaniu rozwaany jest ogólny model systemu przesyłowego (danych, sygnałów, energii elektrycznej, surowców, itp.) zbudowanego z wzłów i łczy tworzcych struktur drzewiast. Niech {e0,...,em} bdzie zbiorem wszystkich komponentów systemu, tj. wzłów i łczy. S one ponumerowane w taki sposób, e jeli ei jest rodzicem ej (ei ej), to i < j, przy czym e0 jest korzeniem struktury. Zadaniem systemu jest przesyłanie danego towaru od wzła ródłowego (korzenia struktury) do wzłów docelowych (lici struktury). Schemat przykładowego systemu składajcego si z 13 wzłów i 12 łczy jest przedstawiony na rys. 1.

Kady z komponentów moe znajdowa si w jednym z dwóch stanów niezawodnociowych: sprawnoci – 1, lub uszkodzenia – 0. Zakładamy, e e0 jest zawsze w stanie 1. Naprawa uszkodzonego komponentu rozpoczyna si tak szybko jak jest to moliwe, tzn., gdy jest wolna jedna z ekip naprawczych. Ze wzgldu na ograniczon liczb tych ekip od uszkodzenia komponentu do rozpoczcia jego naprawy moe upłyn pewien czas.

Rys. 1. Schemat przykładowego systemu przesyłowego

Kolejno, w której komponenty s poddawane naprawie zaley od stosowanej strategii napraw – bd rozpatrywane dwie takie strategie. Funkcjonowanie komponentu ei, 1  i  m, jest okrelone przez trzy funkcje rozkładu: Fi – dystrybuant czasu do uszkodzenia sprawnego ei połczonego z e0, Gi – dystrybuant czasu do uszkodzenia sprawnego ei odłczonego od e0, oraz Hi – dystrybuant czasu do naprawy uszkodzonego ei. Zakłada si , e Fi i Gi s wykładnicze,

w odrónieniu od Hi, która moe by dowoln funkcj rozkładu na przedziale [0,). Zakłada

si te, e Fi Gi, co w sposób formalny wyraa zasad, e komponenty pracujce "pod obcieniem" s bardziej podatne na awari, ni nieobcione. Wynika std, e czas do uszkodzenia ei zaley od zachowania wszystkich komponentów połoonych midzy ei a e0, nie zaley natomiast od zachowania pozostałych komponentów. Z jednym zastrzeeniem – ei funkcjonuje niezalenie od komponentów połoonych midzy ei a e0 do chwili uszkodzenia jednego z nich. Z kolei czas do naprawy uszkodzonego komponentu nie zaley od zachowania wszystkich pozostałych komponentów. Zauwamy te, e Gi 0, 1  i  m, jeli załoymy, e komponenty odłczone od e0 nie uszkadzaj si.

Przesłanie towaru od e0 do ei jest moliwe tylko wtedy, gdy ei jest sprawny i połczony z e0, tzn., gdy wszystkie komponenty połoone midzy e0 a ei (w tym ei) s sprawne. Poniewa komponenty ulegaj awariom, wic okresy, w których e0 jest połczony ze sprawnym ei przeplataj si z okresami, w których ei jest uszkodzony albo odłczony od e0. Celem niniejszego opracowania jest wyznaczenie dla kadego z komponentów rednich długoci obu tych przedziałów czasowych, czyli odpowiedników parametrów MTBF (redni czas midzy uszkodzeniami) i MTTR (redni czas do naprawy). Jako e złoono modelu wyklucza zastosowanie metod analitycznych, bdziemy szacowa powysze parametry posługujc si estymacj statystyczn, przy czym szczególny nacisk zostanie połoony na okrelenie dokładnoci estymacji.

Jest oczywiste, e długoci okresów istnienia połczenia midzy sprawnym ei a e0, oraz okresów braku tego połczenia zale od dwóch dodatkowych czynników: liczby ekip naprawczych wyznaczonych do obsługi systemu, oraz stosowanej strategii naprawczej. Zwikszenie liczby ekip prowadzi do skrócenia rednich okresów braku połczenia, poniewa rednie okresy oczekiwania na napraw ulegaj wtedy skróceniu. Jeli chodzi o strategie naprawcze, to niewtpliwie ma znaczenie, w jakiej kolejnoci s naprawiane uszkodzone komponenty. Bdziemy tu rozwaali dwie strategie. Zgodnie z pierwsz z nich komponenty s naprawiane w takiej kolejnoci, w jakiej si uszkadzaj, czyli kolejka do naprawy jest typu FIFO. Jeli w jednej chwili uszkodzi si wicej ni jeden komponent (takie zdarzenie zachodzi z zerowym prawdopodobie
stwem, chyba, e jest to uszkodzenie o wspólnej przyczynie), to jako pierwszy jest naprawiany ten z najwikszym indeksem. Strategi t nazwiemy "FIFO albo najwikszy indeks pierwszy". Jeli Gi  0, 1  i  m, to dla kadego liniowo uporzdkowanego zbioru komponentów (tj. wszystkich komponentów połoonych midzy e0 a wzłem-liciem) strategia ta nadaje najwyszy priorytet komponentom najbardziej oddalonym od e0. Istotnie, poniewa tylko komponenty połczone z e0 mog si uszkadza, wic jeli ey znajduje si poniej ex (co implikuje y>x), to ey moe si uszkodzi tylko wtedy gdy ex jest sprawny, czyli uszkodzenie ey moe nastpi tylko przed uszkodzeniem ex albo równoczenie z nim, skd wynika, e ey poprzedzi ex w kolejce do naprawy.

Zgodnie z drug strategi komponenty s naprawiane w kolejnoci odwrotnej ni ta, w jakiej si uszkadzaj, tj. tworz kolejk typu LIFO. Jeli w jednej chwili uszkodzi si wicej ni jeden komponent, to jako pierwszy jest naprawiany ten z najmniejszym indeksem. Strategi t nazwiemy "LIFO albo najmniejszy indeks pierwszy". Jeli Gi  0, 1  i  m, to dla kadego liniowo uporzdkowanego zbioru komponentów strategia ta nadaje najwyszy priorytet komponentom najmniej oddalonym od e0. Istotnie, jeli ey znajduje si powyej ex (co implikuje y<x), to uszkodzenie ey moe nastpi tylko po uszkodzeniu ex albo równoczenie z nim, skd wynika, e ey poprzedzi ex w kolejce do naprawy.

2. Podstawowe definicje i notacja

Zdefiniujmy najpierw podstawowe charakterystyki systemu: Li(1) – czas do uszkodzenia sprawnego ei połczonego z e0 Li

(2)

– czas do uszkodzenia sprawnego ei odłczonego od e0 Ri – czs naprawy uszkodzonego ei

Fi, Gi, Hi – funkcje rozkładu zmiennych losowych Li (1)

, Li (2)

, Ri i – intensywno uszkodze sprawnego ei połczonego z e0 1 – strategia “FIFO albo z najwikszym indeksem pierwszy”

2 – strategia “FIFO albo z najmniejszym indeksem pierwszy”

Aby nie komplikowa oznacze, w pozostałych definicjach zakłada si, e liczba ekip naprawczych wynosi r, a stosowana jest strategia s, s∈{1,2}.

Aj(i) – długo j-ego okresu, podczas którego sprawny ei pozostaje połczony z e0 Bj

(i)

– długo j-ego okresu, podczas którego ei pozostaje uszkodzony lub odłczony od e0 a(i), b(i) – przecitne długoci zdefiniowanych powyej okresów

X – proces uszkodze i napraw rozwaanego systemu, tj. X = {[X1(t),...,Xm(t)], t0}, gdzie zmienna losowa Xi(t) jest zdefiniowana nastpujco:

Xi(t) = –q: w chwili t komponent ei jest na miejscu q w kolejce do naprawy Xi(t) = 0: ei jest w trakcie naprawy

Xi(t) = 1: ei jest sprawny i połczony z e0 Xi(t) = 2: ei jest sprawny i odłczony od e0

!k – chwila k-ego powrotu X do stanu pocztkowego, tj. stanu [1,…,1], k  0, przy czym !0 = 0

Qk (i)

– liczba przywróce połczenia midzy sprawnym ei a e0 w interwale (!k–1, !k], k  1 Uk(i) – całkowity czas, przez który w interwale (!k–1, !k] komponent ei pozostaje uszkodzony lub odłczony od e0, k  1

Formalnie, parametry a(i) i b(i) s definiowane nastpujco:

(1)

gdzie e _j e_i oznacza, e ej znajduje si midzy e0 a ei, albo ej=ei. Wzór (2) jest bezporedni konsekwencj faktu, e ei funkcjonuje niezalenie od komponentów połoonych midzy e0 a ei, ale tylko do chwili uszkodzenia jednego z nich. Z kolei b(i) _{jest w ogólnym przypadku zmienn losow} (do sprecyzowania pozostaje rodzaj zbienoci), lecz przy pewnych załoeniach (o których mowa w dalszym cigu) b(i) jest wartoci stał – granic, do której [B1(i)+...+Bn(i)]/n zbiega według prawdopodobiestwa przy n dcym do nieskoczonoci.

Kady interwał [τk–1, τk), k1, bdzie nazywany cyklem operacyjnym procesu X. Poniewa L1(1),…, Lm(1) s niezalene i maj rozkłady wykładnicze, wic X mona podzieli na niezalene i stochastycznie identyczne podprocesy {X(t), t∈[!k–1, !k)}, k1. W konsekwencji, (Qk(i), k  1) i (Uk(i), k  1) s cigami niezalenych zmiennych losowych o jednakowych rozkładach, co ma istotne znaczenie dla dalszych rozwaa. W zwizku z tym przyjmiemy dla uproszczenia, e Q(i) = Q1

(i)

oraz U(i) = U1 (i)

.

3. Estymacja parametru b(i)

W [5] przedstawiono algorytm symulujcy przebieg procesu X. W trakcie działania tego algorytmu dla kadego i∈{1,...,m} symulowane s zmienne losowe Qk(i), k  1 oraz Uk(i), k  1. W tym samym opracowaniu udowodnione zostały dwa lematy tworzce baz teoretyczn dla estymacji b(i)

. Pierwszy z nich okrela warunki dostateczne do tego, aby proces X był powracajcy. Teza drugiego z tych lematów mówi, e jeli X jest powracajcy, oraz 0 < rmin Ri rmax < , 1  i  m, to [B1(i)+…+Bn(i)]/n zbiega według prawdopodobiestwa do E[U(i)]/E[Q(i)] przy n dcym do nieskoczonoci, czyli

(3) b(i)=E(U(i))/E(Q(i))

Losowy proces " = {"(t), t0} jest nazywany powracajcym, jeli ma nastpujce własnoci: 1. Stan " w chwili t=0 jest ustalony, tj. "(0) ma rozkład jednopunktowy.

2. Z prawdopodobiestwem jeden " powraca do stanu "(0) w skoczonym czasie, czyli Pr(!1 < ∞) = 1, gdzie !1 jest (losow) chwil pierwszego (od chwili 0) powrotu " do stanu "(0).

3. Proces "1 = {"(!1+t), t0} jest stochastycznie identyczny z procesem " (" “zaczyna si na nowo” w chwili !1).

Czasami druga z powyszych własnoci jest zastpowana silniejsz, tzn. E(!1)<. Wicej szczegółów na temat procesów powracajcych mona znale  w [4].

Zgodnie z (3), jako estymatora b(i) uyjemy ilorazu estymatorów E[U(i)] i E[Q(i)]. Powstaje tu naturalne pytanie – jaka jest dokładno takiej estymacji? W kategoriach estymacji przedziałowej problem polega na wyznaczeniu przedziału ufnoci dla ilorazu wartoci przecitnych dwóch zmiennych losowych, które – co naley podkreli – nie s niezalene. Szukany przedział ufnoci okrela nastpujce twierdzenie, którego dowód znajduje si w [5].

Twierdzenie 1. Niech X ≥ 0, Y ≥ ymin > 0 bd zmiennymi losowymi o skoczonych wartociach oczekiwanych µX, µY, oraz skoczonych odchyleniach standardowych σX, σY (ymin jest pewn stał). Niech MX,K i MY,L bd rednimi z prób losowych ze zmiennych X i Y, o rozmiarach K i L. Niech (4) ] L   , K   [ max y q 2 # ₂ X Y Y X min /4 1 = −

gdzie q1–/4 jest kwantylem rzdu 1 – /4 standaryzowanego rozkładu normalnego, czyli (5) 4  1 ) q Pr(Z≤ _1−/4 = −

gdzie Z ma rozkład normalny o wartoci przecitnej równej 0 i wariancji równej 1. Wtedy, dla odpowiednio duych K i L, zachodzi nastpujcy wzór:

(6) # )    M M Pr( _ Y X L Y, K X, ≤ > −

tzn. [MX,K/MY,L – #α, MX,K/MY,L + #α] jest przedziałem ufnoci dla X/Y o poziomie ufnoci 1 – .

W celu dalszego uproszczenia notacji, zamiast U(i) i Q(i) bdziemy uywali oznacze U i Q, przyjmujc, e indeksem domylnego komponentu jest i. Niech U* and Q* bd zmiennymi losowymi o nastpujco zdefiniowanych dystrybuantach:

(7) FU*(u) = Pr(U*<u) = Pr(U<u|Q  1), FQ*(q) = Pr(Q*<q) = Pr(Q<q|Q  1)

Zatem U* i Q* s "warunkowymi" zmiennymi losowymi odpowiadajcymi U i Q, gdzie warunkiem jest, aby ei uszkadzał si albo tracił połczenie z e0 przynajmniej jeden raz w interwale (τ0, τ1). Istotne jest, e Q*  1. Poniewa E(U|Q=0)=0, wic

(8) E(Q*) E(U*) 1) Q | E(Q 1) Q | E(U 0) Q E(Q, 1) Q E(Q, 0) Q E(U, 1) Q E(U, EQ EU = ≥ ≥ = = + ≥ = + ≥ =

W konsekwencji, zamiast E(U)/E(Q) moemy estymowa E(U*)/E(Q*) otrzymujc ten sam wynik. Zauwamy raz jeszcze, e w twierdzeniu 1 nie zakłada si niezalenoci zmiennych losowych X i Y, a wic mona je zastosowa do silnie zalenych zmiennych losowych, takich jak U* i Q*. W rezultacie stwierdzamy, e przedział [MX,K/MY,L – #α, MX,K/MY,L + #α ], gdzie

(9) ,  ) K   ( max 2qz #_{ 1/4} U* Q* Q* U* L − =

jest przedziałem ufnoci dla E(U*)/E(Q*) o poziomie ufnoci 1 – . Zauwamy te, e twierdzenia 1 nie mona zastosowa do U i Q, poniewa Pr(Q(i)=0) > 0 dla niektórych i ze zbioru {1,...,m}, wic nie jest prawd, e Q  qmin > 0.

Podsumowujc, stwierdzamy co nastpuje: jeli liczby całkowite K i L spełniaj ponisze warunki: (10) 2 * U * Q /4 1 2 * Q * U /4 1 2q   L ,   2q K _      ≥       ≥ − − ε ε

to generujc próby losowe z U* i Q* o rozmiarach K i L, a nastpnie wyliczajc iloraz ich rednich, estymujemy E(U*)/E(Q*) z zadan dokładnoci ε (połowa długoci przedziału ufnoci) na zadanym poziomie ufnoci 1 – . W praktyce, wartoci przecitne i odchylenia standardowe zmiennych losowych U* i Q* s zastpowane we wzorze (10) odpowiednimi rednimi lub wariancjami z próby.

Wynikiem powyszych rozwaa jest nastpujcy algorytm estymujcy E(U*)/E(Q*): • Wykonaj Np “symulacji pilotowych” (tj. zasymuluj Np cykli operacyjnych X) w celu

wyznaczenia przyblionych wartoci parametrów $U*, $U*, $Q*, oraz $Q* • Znajd minimalne liczby całkowite K i L spełniajce warunek (10)

• Wykonaj ponisz procedur estymujc E(U*)/E(Q*). Procedura 1

j1 = 0 ; j2 = 0; do {

symuluj pełny cykl operacyjny X; if (Q≥1) {

if (j1 < K) {

j1 = j1 + 1 ; µU* = µU* + (U – µU*) / j1 ; } if (j2 < L) { j2 = j2 + 1 ; µQ* = µQ* + (Q – µQ*) / j2; } } } while (j1 < K OR j2 < L); return (µU* / µQ*);

Uwaga: µU* i µQ* s aktualizowane na podstawie poniszego wzoru: (11)

µ

=

µ

+

(x

−

µ

)/j

gdzie

(12)

_µ

Na bazie Procedury 1 został napisany program komputerowy wyznaczajcy dla kadego i∈{1,...,m} podstawowy parametr niezawodnociowy elementu ei – b(i). Dane do programu s pobierane z pliku zawierajcego parametry dystrybuant Fi, Gi i Hi dla poszczególnych elementów, oraz pliku zawierajcego odwzorowanie struktury systemu. Przykładowe wyniki tego programu s przedstawione w tabelach 1 i 2. Zakłada si, e Li(2) = 0 z prawdopodobiestwem 1 (tzn. Gi

0), co oznacza, e sprawny komponent odłczony od e0 nie moe si uszkodzi. Z kolei Li(1) i Ri maj rozkłady wykładnicze z parametrami i = 0.01 i i = 0.1, 1  i  m. Jednostk czasu jest godzina. K oznacza liczb wszystkich symulowanych cykli operacyjnych, K* – liczb cykli, w których Q  1, natomiast T – całkowity czas oblicze. Obliczenia były prowadzone dla i = 8, r = 2, s = 1, 2. Uywano komputera klasy PC z procesorem Intel Core 2 o czstotliwoci 2,14 GHz.

Tabela 1. Wyniki estymacji dla s = 1

\# # = 0.1 # = 0.2  = 0.99 EU*/EQ* = 15.54 EQ* = 1.37 VU* = 425.98 VQ* = 0.56 K = 4,555,000 K* = 2,512,000 T $ 2' EU*/EQ* = 15.56 EQ* = 1.37 VU* = 426.71 VQ* = 0.56 K = 1,141,000 K* = 629,000 T $ 30''  = 0.95 EU*/EQ* = 15.55, EQ* = 1.37 VU* = 426.37 VQ* = 0.56 K = 2.896.000 K* = 1.598.000 T $ 1’15'' EU*/EQ*= 15.55 EQ* = 1.37 VU* = 427.18 VQ* = 0.56 K = 726,000 K* = 401.000 T $ 20''

Tabela 2. Wyniki estymacji dla s=2

\# # = 0.1 # = 0.2  = 0.99 EU*/EQ* = 14.34, EQ* = 1.49 VU* = 492.78 VQ* = 0.92 K = 6,280,000 K* = 3.467.000 T $ 3' EU*/EQ* = 14.33 EQ* = 1.49 VU* = 494.62 VQ* = 0.92 K = 1,577,000 K* = 871.000 T $ 45''  = 0.95 EU*/EQ* = 14.35 EQ* = 1.49 VU* = 494.88 VQ* = 0.92 K = 4,009,000 K* = 2.214.000 T $ 2' EU*/EQ* = 14.33 EQ* = 1.49 VU* = 492.22 VQ* = 0.92 K = 998,000 K* = 550.000 T $ 30''

%LEOLRJUDILD

[1] Billinton R., Allan R.N.: Reliability Evaluation of Power Systems, Plenum Press 1996. [2] Billinton R., Li W.: Reliability Assessment of Electrical Power Systems Using Monte Carlo

Methods, Springer 1994.

[3] Brown R.E.: Electric Power Distribution Reliability, CRC Press 2002.

[4] Feller W.: An Introduction to Probability Theory and its Applications, Wiley 1968.

[5] Malinowski J.: Simulating failure-repair process and evaluating reliability parameters for single-source multiple sink commodity transportation network witch stochastically dependent components, IBS PAN Warszawa 2009.

[6] Malinowski J.: Determining reliability parameters for a tree-structured commodity transfer system with non-independent components, Reliability, Risk and Safety – Theory and Applications, CRC Press 2010, pp. 1564–1568.

[7] Ridgen S.E., Basu A.P.: Statistical Methods for the Reliability of Repairable Systems, Wiley 2000.

SIMULATION METHOD OF DETERMINING RELIABILITY PARAMETERS OF A TRANSMISSION SYSTEM WITH INTERDEPENDENT ELEMENTS

Summary

A general model of a transmission system is considered (of data, signals, electric energy, raw materials etc.), composed of nodes and links, where end nodes are sources or sinks. System components are subject to random damage, removed by a definite number of repair teams. It is assumed that times to damage of individual components have exponential distribution, but that for a component linked to a source node this time has a different distribution than for the same component not linked to the source node (therefrom the dependence). It is further assumed that times to repair of the damaged components have distributions different from exponential. Complexity of the model makes it impossible to determine the reliability parameters of the damage and repair processes of the system in an analytic manner, and so Monte Carlo simulation and statistical estimation are used for this purpose. It turned out that establishing the confidence interval for the parameters estimated is a non-trivial problem. The paper presents the lemmas containing the conditions of renewal of the damage and repair processes, lemmas defining the bounds on the confidence intervals sought, and the computational examples.

Keywords: transmission system, stochastic dependence, reliability parameters, repair policy, Monte Carlo type simulation, statistical estimation.

Jacek Malinowski

Instytut Bada Systemowych PAN Newelska 6, 01-447 Warszawa

Wysza Szkoła Informatyki Stosowanej i Zarzdzania Newelska 6, 01-447 Warszawa

e-mail: jacek.malinowski@ibspan.waw.pl jacmalin@wsisiz.edu.pl