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