Metody wyznaczania niezawodności systemów z czasem aktywacji rezerwy
9.6 Eksperyment obliczeniowy
Zatem zmienna D | a określa pozostały od chwili a czas dostawy, przy założeniu, że dostawa nie zaszła do chwili a.
Kolejne awarie pojazdów zdarzają się średnio co E(A). Zatem jeśli j + 1-wsze uszkodzenie zachodzi w chwili t, to wcześniejsze zaszło (statystycznie) w chwili t − E(A). Jeszcze wcześniejsze w chwili t−2E(A), a najstarsze w chwili t−jE(A). Czas oczekiwania na przemieszczenie pojazdu, który początkowo był na pozycji j + 1, na pozycję j wynosi:
Fjj+1= min{D | E(A), D | 2E(A), . . . , D | jE(A)
| {z }
j
}
W oszacowaniu czasu na kolejne przesunięcie, czyli na pozycję j − 1, należy uwzględnić czas, który już upłynął czekając na awans na pozycję j. Średni czas, po którym to następuje to
E(Fjj+1), a więc:
Fj−1j+1 = min{D | (E(A) + E(Fjj+1)), D | (2E(A) + E(Fjj+1)), . . . , D | ((j − 1)E(A) + E(Fjj+1))
| {z }
j−1
}
Uogólniając, jeśli pojazd znajdował się początkowo na j + 1 pozycji, to czas przesunięcia na dowolną pozycję i jest określony przez:
Fij+1 = min{D | (E(A) + i+1 X k=j E(Fkj+1)), . . . , D | (jE(A) − i+1 X k=j E(Fkj+1) | {z } i )}
Uszkodzonemu pojazdowi zostanie przydzielona rezerwa, gdy ten przejdzie na n-tą pozycję. Łączny czas oczekiwania na rezerwę określa zatem suma:
Cj = n X i=j
Fij+1
Zmienna C oraz jej gęstość fC są ostatecznie równe:
C = ∞ X j=n CjPj = ∞ X j=n Pj n X i=j Fij+1 fC = P∞ j=nfD(j−n+1)Pj P∞ j=nPj
Warto zauważyć, że gdy rozkład dostawy będzie jednak wykładniczy, to bazując na własności braku pamięci można zredukować wszystkie zmienne warunkowe, co spowoduje, że Fij+1 stanie się zmienną wykładniczą. Co za tym idzie, całe rozwiązanie zredukuje się do oszacowania Ia.
9.6 Eksperyment obliczeniowy
W kolejnych 3 podrozdziałach zostaną zidentyfikowane grupy systemów z rezerwą czasową, dla których opisane wcześniej algorytmy są efektywnymi metodami wyznaczania niezawodności.
122 ROZDZIAŁ 9. NIEZAWODNOŚĆ SYSTEMÓW Z REZERWĄ ZIMNĄ . . .
Grupy systemów będą oznaczane ciągiem symboli składającym się z trzech liter, w którym pierwsza oznacza typ dystrybuanty czasu pomiędzy awariami, druga czasu wymiany, a ostatnia czasu dostawy. Przykładowo grupa systemów M/W/M to taka, w której czas pomiędzy awariami ma rozkład wykładniczy (litera M), czas wymiany ma rozkład Weibulla (W), a dostawę ponownie wyraża rozkład wykładniczy.
Obliczenia w każdym z podrozdziałów przeprowadzono z wysoką starannością, tak że każda wartość była liczona przynajmniej dwa razy. Wyniki symulacji całego systemu uzyskiwano dwo-ma niezależnymi narzędziami: wolniejszym, będącym produktem transfordwo-macji modelu READ w kod C++ oraz szybszym, utworzonym specjalnie na potrzeby eksperymentu. Każda z pro-ponowanych wcześniej metod estymacji została zaimplementowana dwa razy: numerycznie oraz metodą Monte Carlo samej zmiennej C. Zarówno symulatory, jak i implementacje estymacji zweryfikowano pod kątem zbieżności ze swoimi odpowiednikami dla ok. 400 badanych.
9.6.1 Model M/W/M
Dokładność oszacowań metody Ia oraz Ib wynika z dwu istotnych założeń. Pierwsze z nich po-lega na zastosowaniu kolejki typu M/G/∞ do wyznaczenia P r(N A), czyli prawdopodobieństwa, że rezerwa nie jest dostępna w chwili uszkodzenia regularnego pojazdu. Istotnym ograniczeniem tego modelu jest wymaganie, aby czas pomiędzy kolejnymi żądaniami był opisany rozkładem wykładniczym. Obecnie nie istnieje rozwiązanie modelu G/G/∞, które pozwoliłoby na relaksa-cję tego podejścia. Drugie założenie związane jest z wykładniczym rozkładem dostawy będącym podstawą do opisu czasu do uzyskania dostępu do rezerwy rozkładem hipowykładniczym lub statystyką pozycyjną. Metoda II nie jest ograniczona wykładniczym rozkładem dostawy.
Pierwsze założenie nie będzie ograniczało jakości estymacji, gdy rozkład czasu pomiędzy awariami będzie wykładniczy. Dodatkowo, jeśli czas dostawy również byłby wykładniczy, to metody Ia i Ib dostarczą dokładnych rozwiązań zmiennej C. Warto zauważyć, że żadna z metod nie ogranicza czasu wymiany, a więc z powyższego wywodu wynika, że model M/W/M powinien być rozwiązywany dokładnie. Celem eksperymentu obliczeniowego w niniejszym podrozdziale jest weryfikacja tego spostrzeżenia.
Czas pomiędzy uszkodzeniami oraz dostawa mają w eksperymencie rozkłady wykładnicze, a wymiana ma rozkład Weibulla o kształcie 1,5. Co więcej, przyjęto 2 możliwości rezerwy czaso-wej oraz 5 stopni rezerwowania, co łącznie daje 10 przypadków testowych zebranych w tab. 9.2. Wartości uzyskanego prawdopodobieństwa hazardu oraz średniego jego czasu okazały się zbieżne z estymacjami dla wszystkich przypadków testowych i te wartości umieszczono w ostat-nich dwu kolumnach tabeli. Zatem w przypadku modeli systemów z rezerwą czasową typu M/W/M, symulacja okazuje się niepotrzebna, a parametry hazardu można otrzymać metodami analitycznymi. Wynikający z tego zysk w czasie obliczeń podsumowano w tab. 9.3.
Obydwie metody Ia i Ib okazały się wyraźnie szybsze niż silnie zoptymalizowany symulator, ale to druga z nich jest najwydajniejsza dostarczając wyniki w czasie ok. 0,8s. Metoda oparta o statystyki pozycyjne jest szybsza z następującego powodu.
9.6. EKSPERYMENT OBLICZENIOWY 123
Tabela 9.2: Przypadki testowe dla modelu M/W/M
Nr przypadku A E D RC [min] n Pr. hazardu
P r(H) Średni czas hazardu E(HT ) [min] λA kE λE λD 1 0,02 1,5 0,05 0,01 41 1 0,53134 62,44 2 0,02 1,5 0,05 0,01 41 2 0,32659 41,83 3 0,02 1,5 0,05 0,01 41 3 0,16704 28,89 4 0,02 1,5 0,05 0,01 41 4 0,08304 19,19 5 0,02 1,5 0,05 0,01 41 5 0,04975 12,27 6 0,02 1,5 0,05 0,01 101 1 0,21212 57,30 7 0,02 1,5 0,05 0,01 101 2 0,08022 38,63 8 0,02 1,5 0,05 0,01 101 3 0,02219 28,74 9 0,02 1,5 0,05 0,01 101 4 0,00483 22,68 10 0,02 1,5 0,05 0,01 101 5 0,00088 18,67
Tabela 9.3: Porównanie czasów obliczeń Czas obliczeń [s]
Nr przypadku metoda Ia metoda Ib Symulator
1 3,28 0,81 75,88 2 2,58 0,81 76,95 3 2,17 0,80 76,66 4 1,69 0,80 75,31 5 1,41 0,78 76,84 6 3,08 0,81 76,24 7 2,55 0,80 75,27 8 2,12 0,80 75,47 9 1,74 0,80 75,38 10 1,42 0,78 75,67
Jakkolwiek zmienna C we wzorze 9.21 wynika z ważonej sumy nieskończonego ciągu zmien-nych o rozkładach hipowykładniczych odpowiadającym kolejnym miejscom w kolejce, to imple-mentacja metody Ia kończy obliczenia, gdy iloczyny Cjsą bardzo małe i dalsze obliczenia wnoszą bardzo niewiele. Ponieważ czas wyznaczenia wartości gęstości rozkładu hipowykładniczego jest proporcjonalny do kwadratu liczby stopni tego rozkładu (wzór 9.20), to czas potrzebny na ob-liczenie zmiennej C zależy od trzeciej potęgi maksymalnej głębokości kolejki. Natomiast czas obliczenia statystyki pozycyjnej w metodzie Ib jest proporcjonalny do głębokości, co wynika z operacji potęgowania (wzór 9.22), więc cały algorytm zależy od kwadratu głębokości (9.23).
Z powyższego spostrzeżenia wynika również dlaczego czas obliczeń w metodzie Ia spada, gdy wzrasta liczba pojazdów rezerwowych, co widać osobno w przypadkach 1-5 oraz 6-10. Ponieważ przegląd kolejki rozpoczyna się od miejsca tuż za ostatnim rezerwowym, wymagany nakład obli-czeniowy jest mniejszy. Powyższego efektu nie widać w metodzie Ib ponieważ zysk ze skracania kolejki jest konsumowany przez wyższy stopień operacji potęgowania.
124 ROZDZIAŁ 9. NIEZAWODNOŚĆ SYSTEMÓW Z REZERWĄ ZIMNĄ . . .
9.6.2 Model W/W/M
W ogólności jakość aproksymacji w proponowanym modelu obliczeniowym wynika z popraw-ności szacowania dwu czynników: P r(N A) oraz fC. Model W/W/M to taki, w którym czasy pomiędzy awariami oraz wymiana są opisane rozkładem Weibulla, podczas gdy dostawa ma wciąż rozkład wykładniczy. Zatem jedyną przyczyną błędów w szacowaniu niezawodności syste-mów o naturze W/W/M metodami Ia, Ib oraz II będzie niedokładność P r(N A) wynikająca z zastosowania modelu kolejkowego M/G/∞ zamiast G/G/∞. Gęstość zmiennej C jest w takich przypadkach dokładna, gdyż wymienione metody poprawnie wyznaczają czas oczekiwania na rezerwę, gdy dostawa ma rozkład wykładniczy.
Eksperyment obliczeniowy przeprowadzono przy zmieniającym się parametrze kształtu czasu pomiędzy awariami. Dane wejściowe zebrano w tab. 9.4.
Tabela 9.4: Zmienne w testach modelu W/W/M
Zmienna Wartości A λ = 0, 016, k = 0, 9 . . . 1, 1; krok 0,01 E λ = 0, 026, k = 1, 25 D λ = 0, 07 RC 41 lub 101 n 1 . . . 5
Metoda Ia, Ib oraz II nie różnią się algorytmem wyznaczania P r(N A), więc dostarczają jednakowe błędy względne przedstawione na rys. 9.1-9.4, a zdefiniowane następująco:
δH(kA, RC, n) = estH(kA, RC, n) − symH(kA, RC, n)
symH(kA, RC, n) ∗ 100% δHT(kA, RC, n) = estHT(kA, RC, n) − symHT(kA, RC, n)
symHT(kA, RC, n) ∗ 100%
Gdzie:
• δH, δHT - błąd względny P r(H) oraz E(HT ) w zależności od kształtu zmiennej A, rezerwy czasowej RC oraz liczby elementów rezerwowych n,
• estH, estHT - wyniki estymacji parametrów hazardu,
• simH, simHT - wyniki symulacji parametrów hazardu.
Ponieważ wśród zbioru parametrów kształtu awarii jest wartość 1, 0 w zbiorze testowanych danych istnieją modele M/W/M, które rozwiązywane są dokładnie. Stąd na rysunkach można zaobserwować, że δH = 0 oraz δHT = 0, gdy awaria ma rozkład wykładniczy.
Gdy kształt awarii oddala się od 1, zarówno ku wartościom większym, jak i mniejszym, błędy metody zwiększają się. Warto zauważyć, maksymalny błąd jest mniejszy od 12%, co w większości zastosowań okaże się wystarczające i pozwoli uniknąć czasochłonnej symulacji działania systemu.
9.6. EKSPERYMENT OBLICZENIOWY 125
Rysunek 9.1: Błędy względne w szacowaniu P r(H) w modelu W/W/M, RC=41 min [99]
126 ROZDZIAŁ 9. NIEZAWODNOŚĆ SYSTEMÓW Z REZERWĄ ZIMNĄ . . .
Rysunek 9.3: Błędy względne w szacowaniu E(HT ) w modelu W/W/M, RC=41 min [99]
9.6. EKSPERYMENT OBLICZENIOWY 127
Analizując bardziej precyzyjnie, jeśli czas awarii jest zmienną o rozkładzie Weibulla o tzw. „ciężkim ogonie”, czyli k < 1, to metody dostarczą rozwiązania niedoszacowane. Przeciwnie, gdy „ogon” jest lekki i k > 1, to wartości estymowane będą zbyt duże w stosunku do rezultatów symulacyjnych. Zjawisko wynika z działania modelu M/G/∞, który dostarczy mniejsze P r(N A) niż w rzeczywistości, gdy awaria ma ciężki „ogon” i większe, gdy „ogon” jest lekki. Zważywszy, że P r(H | N A) jest większe niż P r(H | AV ), to, zgodnie ze wzorem 9.3, mniejsze P r(N A) przyczyni się do zmniejszenia P r(H) i odwrotnie. Analogicznie dla średniego czasu hazardu.
Zgodnie z przypuszczeniami prawdopodobieństwo hazardu jest monotoniczne względem licz-by pojazdów rezerwowych i wyraźnie spada, gdy pula rezerwowa jest liczniejsza, ale błąd względ-ny oszacowań wśród badawzględ-nych przypadków testowych ma ekstremum, gdy n = 4.
9.6.3 Model W/W/W i podsumowanie eksperymentu
W najbardziej ogólnym z rozpatrywanych przypadków czas pomiędzy awariami, wymiana oraz dostawa mają rozkłady Weibulla. Taka konfiguracja wejściowa stanowi największe wyzwanie dla proponowanych metod obliczeniowych.
Analogicznie do poprzedniego eksperymentu w obliczeniach przyjęto maksymalne parame-try kształtu: czasu między awariami oraz wymiany. Kształt dostawy natomiast zmieniał się w zakresie 1, 0 . . . 1, 5. Biorąc pod uwagę jedną rezerwę czasową oraz trzy stopnie rezerwowa-nia otrzymuje się łącznie 60 przypadków testowych opisanych w tab. 9.5. Rezultatem obliczeń każdego z nich są błędy względnego prawdopodobieństwa hazardu oraz średniego jego czasu.
Tabela 9.5: Zmienne w testach modelu W/W/W
Zmienna Wartości A λA= 0, 02, kA= 1, 1 E λE = 0, 05, kE = 1, 5 D λD = 0, 014, kD = 1, 0 . . . 1, 5; krok: 0,1 RC 41 lub 101 n 1 . . . 5
Na każdym z 6 wykresów rys. 9.5 metoda II okazała się wyraźnie lepsza nie tylko od schematu dedykowanego systemom o zmiennych wykładniczych, ale również od metody wykorzystującej statystyki porządkowe i dopuszczającej niewykładniczy czas wymiany. Dominacja metody II rośnie wraz ze zwiększaniem współczynnika kształtu dostawy.
Dodatkowo, odstępy na wykresach prawdopodobieństwa oraz średniego czasu hazardu po-między liniami punktowanymi a przerywanymi rosną, gdy zwiększa się stopień rezerwowania. Można to wyjaśnić następująco. Gdy n jest większe, to szansa na otrzymanie elementu rezer-wowego w chwili awarii również rośnie, czyli P r(N A) maleje. Tym samym prawdopodobieństwo wyprzedzenia dostawy przez wymianę rośnie, czyli proces wymiany ma większe znaczenie dla niezawodności. Konsekwentnie prawidłowy opis wymiany w modelu obliczeniowym przyczyni się do poprawy jakości oszacowania, czyli linia reprezentująca wyniki modelu Ib obniża się w
sto-128 ROZDZIAŁ 9. NIEZAWODNOŚĆ SYSTEMÓW Z REZERWĄ ZIMNĄ . . .
sunku do linii metody Ic (np. rys. 9.5e). Przeciwnie, gdy P r(N A) jest bliskie 1, to zysk z metody opartej o statystyki porządkowe jest niewielki, gdyż prawidłowe szacowanie wymiany ma wtedy marginalne znaczenie. W takiej sytuacji metoda Ib jest nieznacznie lepsza od Ic dedykowanej systemom M/M/M (np. rys. 9.5a).
Naturalnie, dodawanie pojazdów rezerwowych zmniejsza prawdopodobieństwo hazardu. Pa-miętając o wniosku wysuniętym w poprzednim paragrafie, należy uznać, że proponowane metody są szczególnie przydatne w systemach o wysokiej lub bardzo wysokiej niezawodności. Prawdopo-dobieństwa hazardu na wykresach 9.5a, 9.5c oraz 9.5e kształtują się na poziomie: 0, 002, 0, 0001 oraz 0, 0001. Stosowanie metody o zmiennych wykładniczych w systemach o wysokim stopniu rezerwowania prowadzi do nadmiernych błędów, które mogą być w dużym stopniu ograniczone poprzez użycie proponowanych metod. Gdy kształt dostawy jest bliski 1, to zarówno metody I, jak i metoda II okażą się odpowiednie. Gdy jednak w badanym systemie kształt dostawy jest daleki od 1, to jedynie metoda II może być rekomendowana, która w takiej sytuacji wyraźnie przewyższy nie tylko podstawową metodę bazującą na zmiennych wykładniczych, ale również bardziej skomplikowaną wykorzystującą statystyki porządkowe.
W ramach podsumowania w tab. 9.6 porównano błędy w szacowaniu niezawodności określo-nej klasy systemu przez każdy z proponowanych algorytmów. Konfrontując jakość oszacowań z czasem potrzebnym na ich uzyskanie (tab. 9.7) można zauważyć, że najlepszymi schematami są: metoda Ib dla wszystkich klas oprócz W/W/W oraz metoda II dla klasy W/W/W.
Tabela 9.6: Jakość estymacji proponowanymi metodami na tle 4 klas systemów z rezerwą czasową
Metoda Model
M/M/M M/W/M W/W/M W/W/W
Ic Wynik dokładny Błąd umiarkowany Błąd duży Błąd bardzo duży Ia Wynik dokładny Wynik dokładny Błąd niski Błąd duży Ib Wynik dokładny Wynik dokładny Błąd niski Błąd duży
II Wynik dokładny Wynik dokładny Błąd niski Błąd umiarkowany
Tabela 9.7: Średnie czasy obliczeń dla jednego przypadku testowego Metoda Przeciętny czas obliczeń [s]
Ib 2,5
Ic 0,8
II 16,1
9.6. EKSPERYMENT OBLICZENIOWY 129
(a) Prawdopodobieństwo hazardu, n = 3 (b) Średni czas hazardu, n = 3
(c) Prawdopodobieństwo hazardu, n = 4 (d) Średni czas hazardu, n = 4
(e) Prawdopodobieństwo hazardu, n = 5 (f) Średni czas hazardu, n = 5
Rysunek 9.5: Błąd względny prawdopodobieństwa hazardu (lewa kolumna) oraz średniego czasu hazardu (prawa kolumna) w zależności od współczynnika kształtu dostawy kD; linia punktowana - metoda Ic, linia przerywana - metoda Ib, linia ciągła - metoda II [99]