Bartłomiej Piekarski

171160

Data utworzenia: 01.06.2010r.

Łukasz Tkacz

171032

Łukasz Przywarty

171018

Zadanie projektowe:

Niezawodność i diagnostyka układów cyfrowych

Temat: „Ocena niezawodności systemu pomiarowego typu '2z3' z zawod-

nym układem diagnostyki”

Prowadzący:

dr inż. Kazimierz Kapłon

Strona 1

Spis treści

I. Szczegółowy opis zadania projektowego . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .3

II. Koncepcja symulacji . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

a) Wykresy czasowe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

b) Tabela stanów . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

c) Graf stanów . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .6

d) Schemat blokowy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

III. Implementacja programowa modelu symulacyjnego w środowisku Matlab . . . . . . . . . . . . 9

a) Listing kodu. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

b) Wyniki symulacji. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

IV. Szacowanie charakterystyk niezawodnościowych . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

V. Podsumowanie i wnioski. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

VI. Bibliografia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

Strona 2

I. Szczegółowy opis zadania projektowego

W projekcie należy określić charakterystyki niezawodnościowe systemu pomiarowego,

który składa się z trzech elementów pomiarowych (czujników) i układu diagnostyki i przełączania.

Czujniki są funkcjonalne i niezawodnościowo takie same. Do poprawnej pracy systemu

pomiarowego konieczna jest sprawność co najmniej dwóch czujników. System wykonuje swoją

pracę w sposób ciągły. W ujęciu niezawodnościowym wszystkie elementy są zawodne i napra-

wialne. Znane są: rozkład prawdopodobieństwa czasu poprawnej pracy czujnika do uszkodzenia

F

c

(t) i układu diagnostyki F

d

(t) oraz rozkłady G

c

(t) i odpowiednio G

d

(t) czasu do zakończenia

naprawy. Elementy systemu są naprawiane przez K ekip remontowych. Ponadto:

•

układ diagnostyki i przełączania czujników w stanie uszkodzenia nie zakłóca pracy

czujników (tzn. następuje utrata zdolności do wykrywania uszkodzeń i przełączania

elementów, ale jeśli co najmniej 2 z 3 czujników są w stanie sprawności, to system nadal

poprawnie wykonuje pomiary),

•

czasy do uszkodzenia elementu τ

c

, τ

d

mają rozkłady wykładnicze o parametrach λ

c

, λ

d

•

czasy naprawy Θ

c

, Θ

d

mają rozkłady normalne obcięte (Θ > 0) o parametrach μ

c

, σ

c

, μ

d

, σ

d

(wartości dobrać indywidualnie)

W ramach projektu należy:

•

opracować model symulacyjny działania systemu uwzględniający występowanie uszkodzeń

i napraw elementów (przedstawić ten model w postaci schematu blokowego)

•

opracować implementację programową modelu symulacyjnego w środowisku Matlab,

•

dla wybranych wartości parametrów wykonać na modelu niezbędną liczbę eksperymentów

•

na podstawie otrzymanych wyników oszacować charakterystyki niezawodności badanego

systemu

II. Koncepcja symulacji

Zadanie polega na utworzeniu modelu symulacyjnego systemu '2z3' z zawodnym układem

diagnostyki. Szczegółowe założenia projektowe zawarte są w punkcie I. niniejszej pracy.

Na samym początku należy przyjrzeć się zachowaniu systemu w poszczególnych chwilach czasu.

Zakładając, że dostępny jest jeden konserwator, w jednej chwili czasu może on naprawiać tylko

jeden element systemu. W momencie gdy uszkodzi się inny element, a konserwator będzie zajęty,

element pozostaje uszkodzony i czeka na naprawę. Szczególnie traktujemy układ diagnostyki:

naprawiamy go w pierwszej kolejności (gdy jest dostępny konserwator, który taką naprawę może

Strona 3

wykonać. W przeciwnym razie układ naprawiany jest zaraz po ukończeniu naprawy innego

elementu).

a) Wykresy czasowe

Przedstawimy teraz przykładowe wykresy czasowe, zakładając, że korzystamy z usług jednego

konserwatora a stany elementów systemu oznaczone są następująco:

•

0 – sprawny

•

1 - uszkodzony w naprawie

•

2 – uszkodzony bez naprawy

Na rysunku poniżej c

1,

c

2,

c

3,

oznaczają kolejne czujniki, natomiast D – układ diagnostyki

Jak można zauważyć na Rys. 1 uzyskaliśmy wstępne wyniki symulacji dla konkretnych zdarzeń

Strona 4

Rysunek 1: Wykresy czasowe dla rozpatrywanego systemu

(uszkodzeń i napraw). Oznaczyliśmy literą Q czasy naprawy, natomiast literą τ – czasy między

uszkodzeniami.

b) Tabela stanów

Pomocne będzie rozpatrzenie wszystkich kombinacji stanów poszczególnych elementów wraz

z uwzględnieniem stanu całego systemu. Niestety, wiele z kombinacji jest niemożliwych, tzn., że

nie spełniają warunków określonych w treści zadania projektowego. Kombinacje stanów

przedstawimy w postaci tabeli. Kolorem szarym zaznaczono niemożliwe kombinacje stanów

poszczególnych elementów., wyróżniono natomiast stany możliwe. Zakładamy, że K = 1.

Stan S

sys

Stan elementu

Stan

systemu

Stan S

sys

Stan elementu

Stan

systemu

D

S1 S2 S3

D

S1 S2 S3

S

0000

0

0

0

0

0

S

0010

0

0

1

0

0

S

0001

0

0

0

1

0

S

0011

0

0

1

1

S

0002

0

0

0

2

S

0012

0

0

1

2

1

S

0020

0

0

2

0

S

1121

1

1

2

1

S

0021

0

0

2

1

1

S

1122

1

1

2

2

S

0022

0

0

2

2

S

1200

1

2

0

0

0

S

0100

0

1

0

0

0

S

1201

1

2

0

1

S

0101

0

1

0

1

S

1202

1

2

0

2

1

S

0102

0

1

0

2

1

S

1210

1

2

1

0

S

0110

0

1

1

0

S

1211

1

2

1

1

S

0111

0

1

1

1

S

1212

1

2

1

2

S

0112

0

1

1

2

S

1220

1

2

2

0

1

S

0120

0

1

2

0

1

S

1221

1

2

2

1

S

0121

0

1

2

1

S

1222

1

2

2

2

1

S

0122

0

1

2

2

1

S

2000

2

0

0

0

S

0200

0

2

0

0

S

2001

2

0

0

1

0

S

0201

0

2

0

1

1

S

2002

2

0

0

2

S

0202

0

2

0

2

S

2010

2

0

1

0

0

S

0210

0

2

1

0

1

S

2011

2

0

1

1

S

0211

0

2

1

1

S

2012

2

0

1

2

1

S

0212

0

2

1

2

1

S

2020

2

0

2

0

S

0220

0

2

2

0

S

2021

2

0

2

1

1

Strona 5

Stan S

sys

Stan elementu

Stan

systemu

Stan S

sys

Stan elementu

Stan

systemu

S

0221

0

2

2

1

1

S

2022

2

0

2

2

S

0222

0

2

2

2

S

2100

2

1

0

0

0

S

1000

1

0

0

0

0

S

2101

2

1

0

1

S

1001

1

0

0

1

S

2102

2

1

0

2

1

S

1002

1

0

0

2

0

S

2110

2

1

1

0

S

1010

1

0

1

0

S

2111

2

1

1

1

S

1011

1

0

1

1

S

2112

2

1

1

2

S

1012

1

0

1

2

S

2120

2

1

2

0

1

S

1020

1

0

2

0

0

S

2121

2

1

2

1

S

1021

1

0

2

1

S

2122

2

1

2

2

1

S

1022

1

0

2

2

1

S

2200

2

2

0

0

S

1100

1

1

0

0

S

2201

2

2

0

1

1

S

1101

1

1

0

1

S

2202

2

2

0

2

S

1102

1

1

0

2

S

2210

2

2

1

0

1

S

1110

1

1

1

0

S

2211

2

2

1

1

S

1111

1

1

1

1

S

2212

2

2

1

2

1

S

1112

1

1

1

2

S

2220

2

2

2

0

S

1120

1

1

2

0

S

2221

2

2

2

1

1

S

2222

2

2

2

2

Tabela 1: Tabela stanów systemu

Okazuje się, że większość kombinacji jest niepoprawna. Możliwe zestawienia stanów elementów są

mniej liczne. Nasz system, przy założeniu, że korzystamy z usług jednego konserwatora przez

większa część czasu pozostaje niesprawny. Stan sprawności systemu zachodzi tylko w 11

przypadkach.

c) Graf stanów

Aby lepiej zobrazować zachowanie systemu w momencie naprawy lub uszkodzenia elementu

przedstawimy graf stanów. Graf stanów przedstawia tylko możliwe stany systemu, oznaczone

kółkiem, w którym widnieje nazwa stanu. Poszczególne przejścia między stanami oznaczone są

strzałkami. Adnotacje przy strzałkach (n

i

lub u

i

) oznaczają naprawę lub uszkodzenie

poszczególnych elementów. Stan sprawności całego systemu wyróżniono pogrubieniem stanów.

Stany niesprawności zaznaczono przerywaną linią.

Strona 6

Zależności między stanami są dość skomplikowane, ale jednoznaczne. W momencie uszkodzenia

lub naprawy elementu przechodzimy do kolejnego stanu. Procedura ta powtarza się przez cały czas

trwania symulacji.

c) Schemat blokowy algorytmu symulacji

W tym miejscu możemy skupić się na algorytmie symulacji w postaci schematu blokowego.

Schemat blokowy zawiera wyczerpujące dane na temat działania systemu.

Strona 7

Rysunek 2: Graf stanów

Strona 8

Rysunek 3: Model symulacyjny systemu - schemat blokowy

Legenda:

•

F

c

(t) oraz F

d

(t) - rozkład prawdopodobieństwa czasu poprawnej pracy do uszkodzenia

czujnika oraz układu diagnostyki

•

G

c

(t) oraz G

d

(t) - rozkład prawdopodobieństwa czasu naprawy czujnika i układu

diagnostyki

•

K – liczba konserwatorów

•

T

sym

– całkowity czas symulacji

•

L

sym

– liczba symulacji

•

T

1

, T

2

, T

3

oraz T

d

– moment uszkodzenia poszczególnych czujników oraz układu

diagnostyki

•

T

x

– czas najbliższego uszkodzenia (wydarzenia)

•

S

x

oraz S

d

– stany czujników i układu diagnostyki

•

L

0 –

liczba oczekujących na naprawę

Opiszemy pokrótce algorytm działania systemu. Na samym początku określamy niezbędne

parametry, takie jak rozkłady prawdopodobieństw, liczba symulacji itd. Losujemy czasy poprawnej

pracy i wybieramy minimum z tych czasów. W momencie gdy uszkodzeniu uległ konkretny

element systemu naprawiamy go, zakładając, ze są dostępni konserwatorzy. Priorytet (jeśli chodzi

o kolejność wykonywanych napraw) ma układ diagnostyki. Czas naprawy jest również parametrem

losowym. Gdy nie ma dostępnych konserwatorów zwiększamy liczbę oczekujących na naprawę

i staramy się o rozpoczęcie naprawy w innym momencie czasu. Jeśli wszystkie elementy zostały

naprawione zwiększamy czas od startu symulacji i zmniejszamy czasy uszkodzeń. Wyniki

poszczególnych prób zapisujemy w postaci plików wyjściowych o konkretnym numerze.

Wszystkie kroki powtarzamy do momentu aż przekroczymy maksymalny czas symulacji

i maksymalną ilość symulacji.

Dokładny opis działania symulatora znajduje się w punkcie IV niniejszej pracy.

IV. Implementacja programowa modelu symulacyjnego w środowisku Matlab

Program symulacyjny oprócz tego, że generuje losowe zdarzenia (uszkodzenia, naprawy)

będzie także obliczał następujące miary niezawodnościowe:

•

TTF – czas do uszkodzenia,

•

TTR – czas do naprawy,

•

MTBF – średni czas między uszkodzeniami,

Strona 9

•

MTTF – średni czas do uszkodzenia,

•

MTTR – średni czas do naprawy.

a) Listing kodu

Opiszemy działanie symulatora na podstawie dołączonych fragmentów kodu. Na samym początku

ustalamy maksymalny czas symulacji oraz liczbę tych symulacji. Dobieramy również rozkłady

prawdopodobieństwa dla czasów poprawnej pracy i czasów naprawy.

clear

Tsym =

500

;

Lsym =

1000

;

F =

13

;

G =

10

;

Listing 1: Podstawowe parametry

Ustalamy ścieżki zapisu plików zawierających czasy do uszkodzenia – TTF oraz czasy do naprawy

TTR. W tych plikach przechowywać będziemy obliczone w dalszej części programu wartości.

sciezka_ttf =

sprintf

(

'wyniki\\ttf.txt'

);

[file_ttf, message] =

fopen

(sciezka_ttf,

'w'

);

%Sprawdzenie poprawności otwarcia

if

file_ttf == -1

disp

(wiadomosc)

return

;

end

;

sciezka_ttr =

sprintf

(

'wyniki\\ttr.txt'

);

[file_ttr, wiadomosc] =

fopen

(sciezka_ttr,

'w'

);

if

file_ttr == -1

disp

(wiadomosc)

return

;

end

;

Listing 2: Ścieżki zapisu plików ttf i ttr

W tym miejscu możemy zaimplementować schemat blokowy działania systemu – w postaci kilku

zagnieżdżonych pętli 'for'. Definiujemy zmienne i losujemy momenty uszkodzenia

for

L = 1:Lsym

%Inicjalizacja zmiennych.

K = 1;

%liczba konserwatorów

Strona 10

T = 0;

%czas symulacji

tx = 0;

%bieżąca chwila (minimum z wylosowanych momentów)

LO = 0;

%liczba oczekujących na naprawę

S = [0,0,0,0];

%stany poszczególnych elementów systemu - pierwszy układ

%diagnostyki, następnie czujniki

Tx = [1,2,3,10];

%czasy uszkodzenia (najbliższego wydarzenia)

%Losowanie momentów uszkodzenia

for

i = 1:4

Tx(i) = exprnd(F);

end

Listing 3: Zmienne i momenty uszkodzenia

Sprawdzimy również czy system jest sprawny. Wykorzystamy do tego celu dodatkową zmienną

sprawne, która będzie przechowywała liczbę elementów sprawnych.

%Zmienna przechowująca ilość elementów sprawnych

sprawne = 0;

for

i = 2:4

%Jeśli stan elementu jest równy 0 (jest sprawny) zwiększamy

%zmienną sprawne o 1

if

S(i) == 0

sprawne = sprawne + 1;

end

end

%W przypadku gdy liczba czujników sprawnych jest równa przynajmniej 2

%cały system jest sprawny (zasada '2z3')

if

sprawne > 1

stan = 0;

%w przeciwnym wypadku układ jest niesprawny

else

stan = 1;

end

Listing 4: Sprawdzenie stanu systemu

Mając stany poszczególnych elementów jak i całego systemu możemy zapisać wyniki do pliku.

Definiujemy ścieżkę zapisu pliku i wpisujemy do niego wartości:

Strona 11

%Ustalenie ścieżki zapisu dla plików wynikowych symulacji

sciezka_proba =

sprintf

(

'wyniki\\sym%d.txt'

,L);

[file_proba, wiadomosc] =

fopen

(sciezka_proba,

'w'

);

if

file_proba == -1

disp

(wiadomosc)

return

;

end

;

%Zapisujemy wyniki do pliku w postaci: stan systemu, stan układu

%diagnostyki, stany czujników, czas od momentu startu symulacji.

fprintf

(file_proba,

'%d\t%d\t%d\t%d\t%d\t%f\r\n'

, stan, S(1), S(2), S(3),

S(4),T);

Listing 5: Zapis stanów elementów i systemu do pliku

Dochodzimy do najważniejszej części programu – symulacji. Dopóki czas jest mniejszy od

maksymalnego czasu symulacji wykonujemy procedury zobrazowane w schemacie blokowym.

Rozpoczynamy od wyszukania minimum z czasów uszkodzeń.

while

T < Tsym

%zwiększamy bieżącą chwilę

tx = 2 * Tsym;

%Szukamy minimum z czasów uszkodzeń

for

i = 1:4

if

Tx(i) < tx

x = i;

tx = Tx(i);

end

end

Listing 6: Szukanie minimum z czasów uszkodzeń

Następnie sprawdzamy czy urządzenie w chwili tx_ było sprawne. Jeśli tak – sprawdzamy czy

układ diagnostyki był sprawny i czy liczba konserwatorów jest większa od 0. W przypadku

uzyskania twierdzącej odpowiedzi wiemy, że urządzenie uległo uszkodzeniu, musimy je więc

naprawić. Zmniejszamy liczbę dostępnych konserwatorów i losujemy czas naprawy. Zaznaczamy

stan elementu jako 1 (w naprawie). W przeciwnym razie zwiększamy liczbę oczekujących na

naprawę i oznaczamy stan elementu jako 2 (czekający na naprawę).

if

S(x) == 0

if

S(1) == 0 && K>0

Strona 12

K = K-1;

Tx(x) = exprnd(G) + tx;

S(x) = 1;

else

Tx(x) = 5

* Tsym;

LO = LO + 1;

S(x) = 2;

end

Listing 7: Sprawdzenie sprawności w chwili tx_

Jeśli element w chwili tx_ był uszkodzony losujemy czas poprawnej pracy, zwiększamy liczbę

konserwatorów i oznaczamy stan elementu jako 0 (sprawny)

else

K = K

+1

;

Tx(x) = exprnd(F) + tx;

S(x) =

0

;

end

Listing 8: Sprawdzenie sprawności w chwili tx_ c.d.

Szczególnie traktujemy układ diagnostyki. Jeśli układ ten oczekuje na naprawę i są dostępni

konserwatorzy naprawiamy go w pierwszej kolejności.

if

S(

1

) ==

2

&& K >

0

K = K -

1

;

Tx(

1

) = exprnd(G) + tx;

S(

1

) =

1

;

LO = LO -

1

;

end

Listing 8: Próba naprawy układu diagnostyki

Nie tylko układ diagnostyki ulega uszkodzeniu. Podobnie dzieje się z czujnikami, dopóki liczba

oczekujących na naprawę czujników jest większa od 0, układ diagnostyki jest sprawny i są dostępni

konserwatorzy, naprawiamy poszczególne elementy w kolejności 1-wszy, 2-gi, 3-ci czujnik.

i

=

2

;

while

LO >

0

&& S(

1

) ==

0

&& K >

0

if

S(

i

) ==

2

K = K -

1

;

Tx(

i

) = exprnd(G) + tx;

S(

i

) =

1

;

Strona 13

LO = LO -

1

;

end

i

=

i

+

1

;

end

Listing 9: Naprawa poszczególnych czujników systemu

Możemy teraz zwiększyć czas od startu systemu a także zmniejszyć czasy najbliższego wydarzenia

i kolejny raz policzyć ile elementów jest sprawnych. Określenie liczby sprawnych elementów

systemu wygląda identycznie jak w Listingu 4, więc pominiemy ten fragment kodu. Wyniki

zapisujemy do pliku w sposób zaprezentowany w końcowej części Listingu 5.

%Zwiększamy czas od startu systemu

T = T + tx;

%Zmniejszamy czasy najbliższego wydarzenia

for

i

=

1

:

4

Tx(

i

) = Tx(

i

) - tx;

end

Listing 10: Korekta czasów

Nasz symulator wygenerował do tej pory uszkodzenia i naprawy a także czasy tych zdarzeń. Mając

do dyspozycji te dane jesteśmy w stanie obliczyć parametry TTF oraz TTR. Rezultaty zapisujemy

do plików tekstowych.

%Otwieramy plik z wynikami konkretnej próby

[file_proba, wiadomosc] =

fopen

(sciezka_proba,

'r'

);

if

file_proba ==

-1

disp

(wiadomosc)

return

;

end

%Z pliku zczytujemy wartości do tablicy

tab =

fscanf

(file_proba,'

%d%d%d%d%d%f'

,[ 6 inf ]);

fclose

(file_proba);

tab = tab';

rozmiar =

size

(tab,

1

);

ostatniStan = tab(

1

,

1

);

%Pierwszy stan systemu

ostatniRazSprawny =

0

;

%Moment ostatniej sprawności

ostatniRazUszkodzony =

0

;

%Moment ostatniego uszkodzenia

%W tym miejscu obliczymy poszczególne czasy do uszkodzenia i czasy do

Strona 14

%naprawy - TTF oraz TTR

for

i

=

2

:rozmiar

%Zczytujemy czas od momentu startu symulacji

czas = tab(

i

,

6

);

%oraz stan systemu

stan = tab(

i

,

1

);

%Jeśli ostatni odczytany stan różni się od bieżącego sprawdzamy czy

%jest równy 1 jeśli tak zapisujemy do pliku obliczoną wartość TTF.

%W przeciwnym razie - wartość TTR

if

ostatniStan ~= stan

if

stan ==

1

czasSprawnego = (czas - ostatniRazSprawny);

ostatniRazUszkodzony = czas;

fprintf

(file_ttf,

'%f\r\n'

,

czasSprawnego);

else

czasUszkodzonego = (czas - ostatniRazUszkodzony);

ostatniRazSprawny = czas;

fprintf

(file_ttr,

'%f\r\n'

,

czasUszkodzonego);

end

end

ostatniStan = stan;

end

end

fclose

(file_ttf);

fclose

(file_ttr);

Listing 11: Obliczenie parametrów TTF i TTR

Aby umożliwić późniejszą analizę wyników umieścimy wartości TTF oraz TTR w tablicach:

%Otwieramy pliki z wartościami TTF oraz TTR

sciezka_ttf =

sprintf

(

'wyniki\\ttf.txt'

);

[file_ttf, wiadomosc] =

fopen

(sciezka_ttf,

'r'

);

if

file_ttf ==

-1

disp

(wiadomosc)

return

;

end

;

sciezka_ttr =

sprintf

(

'wyniki\\ttr.txt'

);

[file_ttr, wiadomosc] =

fopen

(sciezka_ttr,

'r'

);

if

file_ttr ==

-1

disp

(wiadomosc)

Strona 15

return

;

end

;

%Zapisujemy wartości TTF oraz TTR w tablicy

ttf =

fscanf

(file_ttf,'

%f'

, [ 1 inf ]);

ttf = ttf';

ttr =

fscanf

(file_ttr,'

%f'

,[ 1 inf ]);

ttr = ttr';

fclose

(file_ttf);

fclose

(file_ttr);

Listing 12: Zapis wartości TTF i TTR w tablicach

Na samym końcu obliczamy pozostałe parametry MTBF, MTTF, MTTR

fprintf

(

'MTBF: %f\n'

,mean(ttf) + mean(ttr));

fprintf

(

'MTTF: %f\n'

,mean(ttf));

fprintf

(

'MTTR: %f\n'

,mean(ttr));

Listing 13: Obliczenie pozostałych parametrów niezawodnościowych.

b) Wyniki symulacji

Uruchamiając program uzyskaliśmy pliki wynikowe, których fragmenty przedstawiono w Tabeli 2.

Stan

układu

Stan

diagno-

styki

Stan 1.
czujni-

ka

Stan 2.
czujni-

ka

Stan 3.
czujni-

ka

Czas od

początku

symulacji

TBF

MTBF

MTTF

MTTR

0
0
0
0
0
1
0
0
1
1
1
0
1
1
1
1

0
1
0
0
2
2
1
0
0
0
0
0
0
2
2
1

0
0
0
0
0
2
2
1
1
1
0
0
2
2
2
2

0
0
0
1
1
1
0
0
0
2
1
0
0
0
2
2

0
0
0
0
0
0
0
0
2
2
2
1
1
1
1
0

0.000000
7.687255
8.781624
9.834445
9.864509

17.490966
29.341057
32.566749
37.286037
37.961525
44.649317
55.324184
56.910808
58.083068

66.521111

70.361735

-
-
-
-
-

17.490966

-
-

19.795071

0.675488
6.687792

-

12.261491

1.172260
8.438043
3.840624

36.870919 8.776564 28.094354

Tabela 2: Fragment plików wynikowych

Strona 16

IV. Szacowanie charakterystyk niezawodnościowych

Podsumujmy nasze dotychczasowe działania: analizując treść zadania projektowego

utworzyliśmy schemat blokowy działania systemu. Na tej podstawie napisaliśmy program

w środowisku Matlab. Program symuluje działanie systemu, zapisuje wyniki do poszczególnych

plików i oblicza wartości TTF oraz TTR. Wykorzystamy je do szacowania charakterystyk

niezawodnościowych.

Aby ukazać z jaką częstością pojawiają się kolejne wartości TTF i TTR wykorzystamy histogramy.

Histogram w programie Matlab można utworzyć wydając polecenie: hist(dane). Jeśli

w miejsce danych podstawimy ttf uzyskamy wykres:

Postępując podobnie w przypadku ttr wyświetlimy drugi histogram:

Strona 17

Rysunek 4: Histogram częstości wystąpień wartości TTF

Rysunek 5: Histogram częstości wystąpień wartości TTR

Z Rys. 4 oraz Rys. 5 wynika, że zdecydowanie częściej występują małe wartości TTF oraz TTR. Im

większa wartość tym rzadziej ona występuje. Wiedząc jak wyglądają histogramy wartości TTF

i TTR możemy przystąpić do właściwej części zadania – dopasowywania rozkładów.

W trakcie poszukiwań informacji na temat dopasowywania rozkładów do zmiennych losowych

natrafiliśmy na bardzo pomocne narzędzie zaimplementowane w programie Matlab – dfittool.

Na podstawie danych (w naszym przypadku tablic z wartościami TTF oraz TTR) generuje

histogramy, do których można dopasować funkcje rozkładu np. Weibulla, gamma, normalnego.

Narzędzie korzysta z odrębnego graficznego interfejsu użytkownika, wszystkie operacje

wykonujemy za pomocą dostępnych przycisków. Dfittool automatycznie generuje wykresy

oraz potrzebne parametry np. dla rozkładu Weibulla – parametr skali oraz kształtu (odpowiednio

a i b w programie Matlab). Użycie narzędzia jest bardzo proste wystarczy wpisać (np. w oknie

poleceń): dfittool(dane) w naszym przypadku: dfittool(ttf) lub dfittool(ttr)

Po uruchomieniu narzędzia otwiera się osobne okno przedstawione na Rys. 6. Okno to zawiera

kilka ważnych elementów.

Najistotniejsze dla nas są

przyciski znajdujące się nad

wykresem. Aby zaimportować

zbiór danych należy skorzystać

z przycisku Data. Kolejny

raz wyświetli się nowe okno,

które pozwala wybrać, do któr-

ych zmiennych dopasowywać

będziemy rozkłady. Operacje

dodawania danych akceptu-

jemy przyciskiem Create
data set.

Przycisk New Fit... po-

zwala dodać nowe dopaso-

wanie. Podobnie jak w przypadku tworzenia zbioru zmiennych korzystamy z odrębnego okna

wyboru, które umożliwia nazwanie konkretnego dopasowania (Fit name) oraz zdefiniowanie po-

żądanego rozkładu (Distribution). Warto zaznaczyć, że do jednego wykresu możemy dodać

wiele dopasowań, pozwala to łatwiej określić, który rozkład jest najbardziej odpowiedni.

Przyciski Manage fits... służą do zarządzania dopasowaniami, natomiast Evaluate...

Strona 18

Rysunek 6: Główne okno narzędzia dfittool

oraz Exclude... kolejno do badania konkretnych zmiennych oraz wykluczania wartości,

których nie chcemy obejmować dopasowywaniem.

Ważnym elementem okna narzędzia dfittol jest również pole Display type. Pozwala ono

wybrać sposób wyświetlania wykresów np. Density oznacza wyświetlanie w postaci funkcji

gęstości prawdopodobieństwa natomiast Cumulative probability – funkcji rozkładu

prawdopodobieństwa.

Spróbujemy teraz dopasować histogram wartości TTF oraz TTR do kilku rozkładów. Uznaliśmy,

że najodpowiedniejsze będą rozkłady:

•

normalny,

•

eksponencjalny,

•

Weibulla.

Uzyskaliśmy następujące rezultaty (w postaci funkcji gęstości prawdopodobieństwa):

Z analizy przeprowadzonego doświadczenia (Rys. 7 oraz Rys. 8) wynika , że najlepiej dopasowany

do histogramu jest rozkład eksponencjalny. Rozkład normalny kompletnie nie pasuje, natomiast

rozkład Weibulla osiąga zbyt duże gęstości.

Strona 19

Rysunek 7: Dopasowania rozkładów do wartości TTF

Narzędzie dfittol wygenerowało dla rozkładów następujące parametry:

•

dla rozkładu normalnego:

–

mu (średnia μ) – 8.77656

–

sigma (odchylenie standardowe σ) – 9.48085

•

dla rozkładu Weibulla

–

a (parametr skali λ)– 8.45563

–

b (parametr kształtu k) - 0.92366

•

dla rozkładu eksponencjalnego

–

mu (parametr λ) – 8.77656

Niestety – tylko na podstawie dopasowań do histogramów nie możemy jednoznacznie określić,

który rozkład najlepiej pasuje do naszych wartości. Przedstawimy więc te same dane i rozkłady na

wykresie funkcji prawdopodobieństwa. W tym celu w polu Display type narzędzia

dfittol wybieramy opcję Probability plot.

Strona 20

Rysunek 8: Dopasowania rozkładów do wartości TTR

Widzimy, że o ile dla wartości TTF rozkład eksponencjalny jest dość dobrze dopasowany to dla

wartości TTR – nie (linia czerwona oznacza dopasowanie idealne). Porównamy więc wyniki

z pozostałymi rozkładami.

Strona 21

Rysunek 9: Wykres probabilistyczny dla wartości TTF i rozkładu eksponencjalnego

Rysunek 10: Wykres probabilistyczny dla wartości TTR i rozkładu eksponencjalnego

Rozkład normalny:

Podobnie jak w przypadku dopasowania do histogramu rozkład normalny odrzucamy, ze względu

na to, że zdecydowanie odbiega od ideału (Rys. 11 oraz Rys 12)

Strona 22

Rysunek 11: Wykres probabilistyczny dla wartości TTF i rozkładu normalnego

Rysunek 12: Wykres probabilistyczny dla wartości TTR i rozkładu normalnego

Rozkład Weibulla:

W tym momencie pojawia się problem. Z analizy dopasowań do histogramu wynika, że nasz

rozkład najlepiej przybliża rozkład eksponencjalny. Zaprzeczeniem tego są wykresy

Strona 23

Rysunek 13: Wykres probabilistyczny dla wartości TTF i rozkładu Weibulla

Rysunek 14: Wykres probabilistyczny dla wartości TTR i rozkładu Weibulla

probabilistyczne, na których jednoznacznie widać, że najbardziej dopasowany jest rozkład Weibulla

(Rys. 12 oraz Rys. 13). Uznajemy więc, że rozkład Weibulla będzie najlepszym wyborem.

V. Podsumowanie i wnioski

W ramach projektu opracowaliśmy model symulacyjny systemu uwzględniający

występowanie uszkodzeń i napraw elementów. Model ten został przedstawiony w postaci

schematu blokowego. Na tej podstawie utworzyliśmy implementację programową modelu

w środowisku Matlab. Poszczególne zdarzenia – uszkodzenia i naprawy wygenerowaliśmy

funkcjami programu Matlab. Symulator tworzy pliki wynikowe z rezultatami, które są

następnie opracowywane w postaci histogramów i wykresów a także miar

niezawodnościowych – TTF, TTR, MTBF, MTTF, MTTR. Szacowanie charakterystyk

niezawodnościowych pozwoliło określić, że najlepiej dopasowanym rozkładem jest rozkład

Weibulla. Funkcje rozkładów prawdopodobieństwa liczone są według następujących

wzorów:

•

dla rozkładu eksponencjalnego:

f  x = e

−

x

, dla x 0

f  x =0,

dla x0

gdzie: λ - parametr

•

dla rozkładu normalnego

f  x =

1





2

e

−

x−

2

2 

2

,

dla x0

f  x =0,

dla x0

gdzie μ – średnia, σ – odchylenie standardowe

•

dla rozkładu Weibulla

f  x = k  x

p−1

e

−

x

k

, dla x0

f  x =0,

dla x0

gdzie λ – parametr skali, k – parametr kształtu

Uwaga 1: W pracy często pojawia się pojęcie 'uszkodzenia' systemu. Przypominamy:

system jest niesprawny gdy nie działa więcej niż jeden czujnik lub układ diagnostyki jest

niesprawny. Niesprawność systemu oznaczaliśmy zawsze cyfrą 1, sprawność – 0.

Uwaga 2: Kod w postaci pliku Matlab`owego oraz sprawozdanie w wersji elektronicznej są

zawarte na dołączonej do projektu płycie CD.

Strona 24

VI. Bibliografia

Podczas tworzenia niniejszego projektu korzystaliśmy z następujących źródeł:

1. Jarnicki J. „Wykład o estymacji rozkładów zmiennych losowych”

2. Praca zbiorowa, „Niezawodność i eksploatacja systemów”, pod redakcją W. Zamojskie-

go, Wrocław 1981

3. Pratap R., „Matlab 7 dla naukowców i inżynierów”, Warszawa 2006

Miejsce na płytę CD

Strona 25