TECHNICKÁ UNIVERZITA V LIBERCI Fakulta mechatroniky, informatiky a mezioborových studií
DISERTAČNÍ PRÁCE MODELOVÁNÍ DYNAMICKÉ SPOLEHLIVOSTI UŽITÍM MARKOVSKÉ ANALÝZY
Vypracoval:
Ing. Josef Chudoba
Školitel: Školitel specialista:
Ing. Pavel Fuchs, CSc. Ing. David Vališ, Ph.D.
Technická univerzita v Liberci, Studentská 2, 461 17 Liberec
Tel.: 485 353 763
http://www.fm.tul.cz
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Anotace Disertační práce na téma Modelování dynamické spolehlivosti užitím markovské analýzy má za hlavní cíl rozšířit markovskou analýzu ve spolehlivosti o prostředek, který umožňuje zjistit a popsat časovou a výkonovou dynamiku složitých systémů především se síťovou strukturou. Základním předpokladem markovské analýzy jsou konstantní intenzity poruch a oprav mezi stavy. V práci je uveden prostředek na časovou dynamiku, který umožňuje řešit úlohy s nekonstantní intenzitou poruch a oprav mezi stavy. Dále je schopen řešit modely systémů s preventivní údržbou. Matematickým základem řešení je sestavení obyčejných diferenciálních rovnic s nekonstantními parametry a následné řešení metodou Runge-Kutta pro komponenty, které jsou stejně spolehlivé jako staré. Nebo pro komponenty systému, které jsou stejně spolehlivé jako nové, řešit úlohu pomocí metody Monte Carlo. Výkonovou dynamiku lze popsat u systémů s násobným počtem shodných jednotek, u kterých je výsledná spolehlivost závislá na objemu výroby, kombinací jednotek v provozu a typu zálohování (aktivní nebo pasivní). Poznatky z této disertační práce lze použít například při modelování spolehlivosti elektrické, silniční a železniční sítě. Pro aplikační řešení projektu bylo použito modelování kompresorové stanice a přilehlých linií tranzitního plynovodu RWE Transgas. Vypočítáním okamžité nepohotovosti systému se zjišťuje pravděpodobnost, že daná síť nebude schopna uspokojit požadavky zákazníka. To znamená u plynárenské společnosti nedodání plynu v odpovídající kvalitě nebo u ostatních uvedených příkladů např. přetížení sítě. Tuto problematiku dosud neřeší dostupné softwarové nástroje používané ve spolehlivosti. Tato práce by měla výrazně přispět k efektivnímu stanovení pravděpodobnosti vzniku kritických a nežádoucích výpadků při modelování složitých systémů. Vyhodnocením a porovnáním výsledků několika scénářů modelů spolehlivosti se umožní snížit pravděpodobnost vzniku havarijních situací, například účinnější údržbou systému nebo početnějším zálohováním komponent. Účinným snížením pravděpodobnosti vzniku havarijní situace se následně ušetří finanční náklady plynoucí z frekvence výpadku a celkové doby při nedodávání plynu.
Klíčová slova: Dynamická spolehlivost, markovská analýza, výkonová analýza, údržba, Weibullovo rozdělení, kompresorová stanice, RWE Transgas
Strana: 2 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Abstract This dissertation thesis on the topic of „Modelling of dynamical dependability by using stochastic processes“ has the main aim of expanding the Markov analysis (in dependability) by adding an instrument, which allows learning and describing time and performance dynamics of complicated systems, especially network structure. The base hypothesis of the Markov analysis is that failure and repair rates between two postures are constant. The instrument for time dynamics can solve tasks whose failure and repair rates are not constant and also solve repairs in predetermined maintenance. Mathematical solving of this problem is based on a construction of differential equations with non-constant parameters and their solution. The Runge-Kutta method is used for components, which are as „good as old“. The Monte Carlo method is usually used for components of a system, which are as “good as new”. Performance dynamics can be described as systems with multiple counts of the same items. Resultant dependability of these items depends on the production volume and the combination of items, which are in use, or in active and passive redundancy. An example can be mentioned as modelling of dependability of an electrical network and systems for railways and roads. Practical solving of this project was shown on the modelling of a compressor station and adjacent lines of the gas pipeline RWE Transgas. After instantaneous unavailability, a solution is possible to determine the probability that this gas pipeline system is not able to provide gas distribution to customers in required amounts. This could mean that the provider can’t provide gas distribution in the required qualities or another example can be that gas pipeline becomes overloaded etc. Conventional software engines in dependability didn’t solve these problems at that time. This dissertation thesis may help on a large scale by an efficient evaluation of probability of the creation on catastrophic failures by the modelling of complicated systems. Thanks to the revaluation of modelled causes of dependability it is possible to decrease the probability of cataleptic failure. This can be achieved, for example, by more efficient maintenance, or multiple component redundancy. If the probability of cataleptic failure is effectively decreased, it would bring a reduction of cost resulting from the frequency of gas distribution failure and also the total time of failure.
Key words: Dynamic reliability, Markov analysis, performance analysis, maintenance, Weibull distribution, compressor station, RWE Transgas
Strana: 3 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Poděkování: Tato disertační práce byla zrealizována díky finanční podpoře státního rozpočtu České republiky prostřednictvím projektů:
AV ČR číslo 1ET401940412 s názvem Modelování a kvantifikace spolehlivosti dynamických systémů, MŠMT Výzkumná centra (PP2-DP01) číslo 1M0554 s názvem Výzkumné centrum Pokročilé sanační technologie a procesy
Prohlašuji, že jsem tuto předkládanou disertační práci vypracoval samostatně a uvedl jsem veškeré prameny, kterých jsem použil. Datum:
Podpis:
Chtěl bych na tomto místě upřímně poděkovat všem, kteří mi radou, pomocí, konzultací či podporou pomohli při tvorbě disertační práce a bez jejichž přispění by tato práce nikdy nevznikla.
Strana: 4 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Obsah: Anotace........................................................................................................................2 1.
Úvod .................................................................................................................12
2.
Cíle disertační práce .......................................................................................14
3.
Teoretické poznatky........................................................................................16
3.1 3.2 3.3 3.4 3.5
4. 4.1 4.2 4.3 4.4 4.5
5. 5.1 5.2 5.3 5.4
Stochastické procesy ............................................................................................. 17 Základní metody spolehlivosti vhodné k zjištění pohotovosti ........................... 23 Markovská analýza ve spolehlivosti...................................................................... 25 Statistická rozdělení popisující dobu do poruchy komponenty ......................... 30 Softwarové prostředky spolehlivosti využitelné v markovské analýze ............. 39
Markovská analýza a neexponenciální rozdělení doby do přechodu .........40 Modelování opravovaných a neopravovaných objektů....................................... 41 Objekt po obnově stejně spolehlivý jako starý .................................................... 44 Objekt po obnově stejně spolehlivý jako nový .................................................... 50 Příklad ...................................................................................................................... 55 Zaměnitelnost modelů - systém je stejně spolehlivý jako nový/starý ............... 57
Modelování údržbových zásahů ....................................................................58 Modelování údržby pomocí markovské analýzy .................................................. 60 Preventivní údržba pro systém stejně spolehlivý jako starý .............................. 63 Preventivní údržba pro systém stejně spolehlivý jako nový .............................. 66 Příklad ...................................................................................................................... 68
6.
Výkonové konfigurace systému.....................................................................70
6.1 6.2 6.3
Přechod z RBD modelu k přechodovému diagramu............................................ 71 Sestavení matice pravděpodobností přechodů pro systém s více výkonovými režimy.......................................................................................................................... 78 Příklad ...................................................................................................................... 81
7.
Další možné aplikace ......................................................................................83
7.1 7.2
Modelování způsobilosti systému ......................................................................... 83 Modelování zisků a ztrát......................................................................................... 84
8.
Analýza spolehlivosti kompresorové stanice tranzitního plynovodu.........86
9.
Otevřené problémy disertační práce .............................................................87
10.
Závěr.................................................................................................................88
11.
Literatura..........................................................................................................92
Strana: 5 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Seznam obrázků Obr. 1: Vanová křivka ..................................................................................................................... 31 Obr. 2: Závislost intenzity/proudu poruch na čase pro objekty stejně spolehlivé jako starý a stejně spolehlivé jako nový .................................................................................................. 42 Obr. 3: Přechodový diagram řešeného příkladu ............................................................................. 49 Obr. 4: Porovnání nepohotovostí pomocí softwaru Isograph (vlevo) a pomocí MATLABu® (vpravo)................................................................................................................................ 49 Obr. 5: Přechodový diagram řešeného systému............................................................................. 55 Obr. 6: Graf nepohotovosti systému, doba simulace 1000 h, systém stejně spolehlivý jako nový (vlevo), systém stejně spolehlivý jako starý (vpravo) .................................................. 56 Obr. 7: Graf nepohotovosti systému, doba simulace 10 000 h, systém stejně spolehlivý jako nový (vlevo), systém stejně spolehlivý jako starý (vpravo) .................................................. 56 Obr. 8: Rozdělení typů údržby ........................................................................................................ 58 Obr. 9: Graf nepohotovosti systému s preventivní údržbou, systém stejně spolehlivý jako starý... 69 Obr. 10: Graf nepohotovosti systému s preventivní údržbou, systém stejně spolehlivý jako nový ..................................................................................................................................... 69 Obr. 11: RBD diagram a přechodový diagram sériového modelu pro dvě sériové komponenty s rozdílnou intenzitou oprav ................................................................................................. 72 Obr. 12: Přechodový diagram pro sériové komponenty s obdobnou intenzitou oprav ................... 72 Obr. 13: RBD a přechodový diagram pro paralelně aktivně zálohované komponenty.................... 73 Obr. 14: RBD a přechodový diagram pro dvě paralelně zálohované komponenty a poruchou se společnou příčinou .......................................................................................................... 74 Obr. 15: RBD diagram pro paralelně pohotovostně zálohované komponenty ................................ 74 Obr. 16: Přechodový diagram pro paralelně pohotovostně zálohované komponenty .................... 75 Obr. 17: Markovský přechodový diagram pro dvě komponenty v pohotovostní záloze .................. 75 Obr. 18: Přechodový diagram pro m z n aktivně zálohované komponenty .................................... 76 Obr. 19: RBD diagram pro m z n pohotovostně zálohované komponenty...................................... 77 Obr. 20: Přechodový diagram pro komponenty zálohované m z n s pohotovostní zálohou........... 77 Obr. 21: Struktura matice pravděpodobností přechodů pro dvě výkonové konfigurace ................. 79 Obr. 22: RBD diagram řešeného příkladu....................................................................................... 81 Obr. 23: Přechodové diagramy řešeného příkladu ......................................................................... 81 Obr. 24: Funkce okamžité nepohotovosti řešeného příkladu provozovaného ve dvou odlišných výkonových režimech........................................................................................... 82
Strana: 6 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Seznam tabulek Tabulka 1: Chyba při výpočtu střední intenzity přechodu ......................................................47 Tabulka 2: Vzorce pro generování doby do přechodu...........................................................53
Seznam příloh Příloha 1: Softwarové prostředky spolehlivosti - markovská analýza Příloha 2: Software pro modelování systému stejně spolehlivého jako starý Příloha 3: Software pro modelování systému stejně spolehlivého jako nový Příloha 4: Analýza spolehlivosti kompresorové stanice tranzitního plynovodu
Seznam zkratek AV ČR
Akademie věd České republiky
ETA
analýza stromu událostí
FMEA
analýza způsobu a důsledků poruch
FTA
analýza stromu poruchových stavů
LBD
logický blokový diagram
LCC
analýza nákladů životního cyklu
MA
markovská analýza
PN
analýza Petriho sítí
RBD
blokový diagram bezporuchovosti
RCM
údržba zaměřená na bezporuchovost
TUL
Technická univerzita v Liberci
Strana: 7 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Použité značení a zjištění parametrů A
součinitel asymptotické pohotovosti
A(t )
funkce okamžité pohotovosti
A(t ) = ∑ pi (t ) ΩU
t
1 2 A(τ )dτ t2 − t1 ∫t1
součinitel střední pohotovosti
A (t1 , t2 ) =
Amin t1 ,t2
minimální pohotovost v intervalu t ∈ t1 ,t2
Amin t1 , t2 = min A(t )
α ,β
parametry Weibullova rozdělení
C (t )
funkce okamžité způsobilosti
1
A (t1 , t 2 )
N
C (t ) = ∑ Ci ⋅ pi (t ) i =1
způsobilost i -tého stavu
Ci
t
C (t1 , t2 )
2
součinitel střední způsobilosti
1 2 C (t1 , t2 ) = C (τ )dτ t2 − t1 t∫1 t2
CT (t1 , t2 )
3
celková způsobilost
CT (t1 , t2 ) = ∫ C (τ )dτ t1
Costi
zisk za jednotku času ve stavu i
Costij
zisk za každý přechod ze stavu i do stavu j
Costuij
zisk za každý preventivní údržbový zásah ze stavu i do stavu j t
Cost (0, t )
součinitel středního zisku
CostT (t1 , t 2 )
celkový zisk v intervalu < t1 ,t2 >
1 Cost (0, t ) = ∫ Cost (τ )dτ tt
t
1
1 Často A (0, t ) = ⋅ ∫ A(τ )dτ t 0 t
2
1 C (0, t ) = ⋅ ∫ C (τ )dτ t 0 t
3
CT (0, t ) = ∫ C (τ )dτ 0
Strana: 8 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
N
N
N
i =1
i , j =1 i≠ j
i , j =1 i≠ j
CostT (t1 , t2 ) = ∑ pi (t1 ) ⋅ Costi ⋅ (t2 − t1 ) + ∑ pi (t1 ) ⋅ pij (t1 , t2 ) ⋅ Costij + ∑ pi (t1 ) ⋅ uij (t1 , t2 ) ⋅ Costuij
D(T )
rozptyl
E (T )
střední hodnota
F (t )
distribuční funkce
... h1N (t ) ⎤ ... h2 N (t ) ⎥⎥ ... ... ⎥ ⎥ ... hNN (t )⎦
⎡ h11 (t ) h12 (t ) ⎢ h (t ) h (t ) 22 h(t ) = ⎢ 21 ⎢ ... ... ⎢ ⎣hN 1 (t ) hN 2 (t )
h(t )
matice intenzit přechodů
h(t )
intenzita přechodu
hij (t )
intenzita přechodu ze stavu i do stavu j , prvek matice h(t ) t
hij (t1 , t2 )
střední intenzita přechodu
i, j , k
označení stavů
λ (t )
okamžitá intenzita poruch
1 2 h(t1 , t2 ) = h(τ )dτ t2 − t1 ∫t1
t
λ (t1 , t 2 )
4
střední intenzita poruch
l
délka kroku
MTBF
střední doba provozu mezi poruchami
MTTF
střední doba do poruchy
MTTR
střední doba do obnovy
µ (t )
okamžitá intenzita oprav
λ (t1 , t2 ) =
1 2 λ (τ )dτ t2 − t1 t∫1
t
µ (t1 , t 2 )
střední intenzita oprav
n
počet sledovaných komponent
N
počet stavů popisující systém
ν
počet realizací
pi (t )
4
1 2 µ (t1 , t 2 ) = µ (τ )dτ t2 − t1 ∫t1
pravděpodobnost stavu i , prvek pravděpodobnostního vektoru p(t )
λ (0, t ) Strana: 9 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
pravděpodobnost přechodu ze stavu i do stavu j v intervalu < t1 , t2 ) ,
pij (t1 , t 2 )
prvek matice pravděpodobnosti přechodů P(t )
⎡ p1 (t ) ⎤ ⎢ p (t ) ⎥ p(t ) = ⎢ 2 ⎥ ⎢ M ⎥ ⎢ ⎥ ⎣ pN (t )⎦
p(t )
pravděpodobnostní vektor
P(t1 , t2 )
matice pravděpodobnosti přechodů
⎡ p11 (t1 , t2 ) ⎢ p (t , t ) P(t1 , t2 ) = ⎢ 21 1 2 ⎢ ... ⎢ ⎣ pN 1 (t1 , t2 )
p12 (t1 , t2 )
...
p22 (t1 , t2 )
...
... ... p N 2 (t1 , t2 ) ...
R (t )
pravděpodobnost bezporuchového provozu
r
počet poruch
t
čas
ti
časový okamžik poruchy
tr
časový okamžik r -té poruchy
T
celková doba zkoušky
T*
celková kumulovaná doba provozu
TDT (t )
kumulovaná doba nepoužitelného stavu (total down time)
p1N (t1 , t2 ) ⎤ p2 N (t1 , t2 ) ⎥⎥ ⎥ ... ⎥ pNN (t1 , t2 )⎦
t
TDT (t ) = ∫ U (τ )dτ 0
kumulovaná doba použitelného stavu (total up time)
TUT (t )
t
TUT (t ) = ∫ A(τ )dτ 0
U
součinitel asymptotické nepohotovosti
U (t )
funkce okamžité nepohotovosti
U (t ) = ∑ pi (t ) ΩD
t
U (t1 , t 2 )
5
5
součinitel střední nepohotovosti
U (t1 , t2 ) =
1 2 U (τ )dτ t2 − t1 t∫1
U (0, t ) Strana: 10 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
U max t1 ,t 2
maximální nepohotovost v intervalu t ∈ t1 ,t2
U max t1 , t2 = max U (t ) , t ∈ t1 ,t2
U(t1 , t2 )
matice údržby
⎡ u11 (t1 , t2 ) u12 (t1 , t2 ) ⎢ u (t , t ) u (t , t ) 22 1 2 U(t1 , t2 ) = ⎢ 21 1 2 ⎢ ... ... ⎢ ⎣u N 1 (t1 , t2 ) u N 2 (t1 , t2 )
uij (t1 , t2 )
... u1N (t1 , t2 ) ⎤ ... u2 N (t1 , t2 ) ⎥⎥ ⎥ ... ... ⎥ ... u NN (t1 , t2 )⎦
pravděpodobnost preventivní údržby ze stavu i do stavu j v intervalu
< t1 , t2 ) , prvek matice údržby U (t1 , t2 )
X (t )
náhodná veličina
z (t )
parametr proudu poruch
z (t ) = lim
z (t1 , t2 )
střední parametr proudu poruch
1 2 z (t1 , t2 ) = z (τ )dτ t2 − t1 t∫1
ΩU
množina provozních a částečně poruchových stavů
ΩD , Ω A
množina poruchových stavů a množina absorpčních stavů
E [r (t + ∆t ) − r (t )] ∆t → 0 ∆t t
Strana: 11 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
1. Úvod Firmy, které vyrábějí nebo chtějí vyrábět nové technologické celky, musí podle platné legislativy zajistit, aby jejich výrobky byly bezpečné. V rámci udržení si dobrého jména na trhu a také z důvodů ekonomických musí zajistit, aby jejich výrobky splňovaly určité předem dané parametry spolehlivosti a tím pádem i pohotovosti, bezporuchovosti, udržovatelnosti a zajištěnosti údržby. V poslední době lze sledovat trend, při kterém je zřejmý stále se zvyšující zájem o vytváření složitějších modelů spolehlivosti. Zatímco dříve se používal pro popis funkce systému obvykle dvoustavový model, dnes se již často pracuje s modely vícestavovými. Dříve se pravděpodobnost vzniku poruchy komponenty do času t popisovala pomocí exponenciálního rozdělení - F (t ) = 1 − e − λt , při kterém se předpokládá, že komponenta s dobou používání nedegraduje. Dnes se již využívá při popisu pravděpodobnosti vzniku poruchy komponenty do času t i jiných rozdělení. Velmi často například Weibullovo t −( ) β
rozdělení - F (t ) = 1 − e α , pomocí kterého lze popsat degradaci komponenty v čase. Nadto je vhodné uvést, že technologické systémy mají stále sofistikovanější politiku údržby. Při sjednávání obchodních kontraktů je již zcela běžné, že jsou ve smlouvách uváděny i parametry týkající se spolehlivosti systémů. Jako příklad stále složitějších systémů a řešených analýz spolehlivosti jsou uvedeny dva typy zakázek od dvou různých firem pro TUL, které autor vyhodnocoval. V prvním případě bylo potřeba zjistit, na základě dat o poruchách, zda je vhodné komponenty po poruše stále opravovat nebo bude doporučena jejich výměna [101] a [104]. Ve druhém případě se zkoumalo, zda výrobní zařízení (každé v řádové hodnotě sto tisíc euro) budou schopny vytvořit dvakrát víc výrobků než se předpokládalo v návrhu, přičemž žádné z obdobných provozovaných zařízení nevyrobilo ani polovinu předpokládaného počtu výrobků [105]. V obou předkládaných případech je zjištění výsledků dáno především teoretickými základy spolehlivosti, vytvořeným modelem spolehlivosti, znalostí statistiky a samozřejmě nutných interdisciplinárních konzultací odborníků. Jak již bylo výše naznačeno, je nutné stále prohlubovat základní metody spolehlivosti tak, aby byly schopny popsat (z pohledu spolehlivosti, bezpečnosti, ekonomiky a také spokojenosti odběratele a dodavatele) stále složitější systémy. V názvu disertačního tématu je uveden termín „dynamická spolehlivost“. Tento termín je v poslední době hojně využíván a je často chápán ve více významech. Zatím neexistuje ustálená definice tohoto termínu. V [1] na straně 147 je uvedena následující definice: „Dynamická spolehlivost zahrnuje studium stochastických procesů, které popisují dynamiku systému, kdy důraz je kladen zejména na poruchy charakterizované vystoupením nějaké procesní proměnné z bezpečné domény“. Lze popsat více typů dynamických úloh, které jsou ve světě řešeny. Termín dynamická spolehlivost se objevuje často například v následujících úlohách: • •
Parametry spolehlivosti jsou závislé nejenom na době provozu, ale navíc na dalších fyzikálních veličinách, které se sledují a nabývají spojitých hodnot v určitém intervalu a navíc je lze pouze obtížně zdiskrétnit. Doba do poruchy komponenty je popsána pomocí speciálních rozdělení. Dynamická spolehlivost se v tomto případě zabývá použitelností nových modelů a zjištěním optimálních parametrů speciálních rozdělení. Strana: 12 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
• • •
Řešení spolehlivosti systémů pomocí metody sítí. Systém se popisuje pomocí více než dvou stavů, kde doba do poruchy není popsána pomocí exponenciálního rozdělení. Odlišuje se od druhého bodu tím, že se předpokládá nutnost popsat systém více než dvěma stavy. Zjišťování spolehlivosti systému v závislosti na jeho výkonnosti.
V rámci disertační práce se autor bude zabývat posledními dvěma body z výše uvedeného seznamu. Disertační práce si klade za cíl připojit do již užívané markovské analýzy ve spolehlivosti prostředek, který bude popisovat časovou degradaci komponent, a dále pomocí kterého bude možné popsat preventivní údržbu systému. Nadto bude řešit jeho výkonovou dynamiku, kdy spolehlivost systému je závislá na jeho produkovaném výkonu. Tento model je složitější, ale zároveň poskytuje i přesnější výsledky proti současnému stavu tím, že obvykle musí být popsán pomocí nehomogenních markovských procesů. Přínosy disertační práce lze využít pro stanovení spolehlivosti přepravních linií například elektrického vedení, silniční a železniční sítě nebo transitních plynovodů. Dynamická spolehlivost se u přepravních linií projevuje tím, že výsledná pohotovost systému je závislá i na objemu přepravy. V rámci disertační práce bude řešen aplikační příklad, kterým bude popsán model spolehlivosti kompresorové stanice a přilehlých linií tranzitního plynovodu RWE Transgas, a.s. Cílem autora je v aplikačním příkladě použít všech modelovacích prostředků spolehlivosti uvedených v disertační práci a tím přispět ke snížení maximální nepohotovosti systému v čase U max (t1 , t2 ) a pravděpodobnosti vzniku závažných nežádoucích událostí.
Strana: 13 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
2. Cíle disertační práce Hlavním cílem disertační práce je vytvořit model spolehlivosti, který umožní rozšířit modelování pohotovosti vícestavových systémů, u kterých se nepředpokládá asymptotická (ustálená) pohotovost A . Matematickým základem pro výpočet funkce okamžité pohotovosti A(t ) (nepohotovosti U (t ) ) vícestavových systémů je v této práci teorie stochastických procesů. Spolehlivostním základem je metoda markovské analýzy. Vytvořený model umožní (na základě znalosti topologie přechodového diagramu) rozšířit stávající vícestavový markovský model o popis doby do poruchy komponent pomocí neexponenciálního rozdělení (se zaměřením především na Weibullovo) - kapitola 4 a dále o popis preventivní údržby - kapitola 5. Tím bude možné stanovit, na základě vstupních dat, hodnotu funkce okamžité pohotovosti A(t ) systému v čase. Obě kapitoly jsou časově dynamickou složkou pohotovosti systému. V současné době je Weibullovo rozdělení obvykle využíváno pro popis doby do poruchy dvoustavových systémů a zřídkakdy vícestavových systémů. Tím, že intenzity přechodů mezi stavy nejsou u Weibullova rozdělení konstantní, nelze úlohu řešit pomocí markovských homogenních procesů. Za předpokladu, že jsou komponenty po obnově stejně spolehlivé jako staré (as good as old), se úloha řeší pomocí soustavy diferenciálních rovnic s nekonstantními parametry. V druhém případě, kdy komponenty mají po obnově vlastnost, že jsou stejně spolehlivé jako nové (as good as new), je nutné využít metod řešení pomocí stochastických procesů - například metody Monte Carlo. Oba dva přístupy jsou natolik odlišné, že první řešící úlohu pomocí soustavy diferenciálních rovnic je uveden v kapitole 4.2, druhý pomocí metody Monte Carlo je uveden v kapitole 4.3. Kapitola 5 se zabývá modelováním systémů s preventivní údržbou. Obnova komponenty po poruše se v přechodovém diagramu modeluje pomocí přechodu ze stavu i do stavu j . Přechod je ohodnocen intenzitou opravy hij (t ) . Na systémech se však velice často provádí preventivní údržba, kdy se komponenty v definovaném čase vymění, nebo se zkontroluje jejich provozuschopnost. Tímto typem údržby přejde v předem daném čase t model systému s určitou stanovenou pravděpodobností ze stavu i do stavu j . Tento typ údržby způsobí nespojitost funkce okamžité pohotovosti A(t ) systému. Na systém po preventivní údržbě lze pohlížet jako na nový systém s jinými pravděpodobnostmi jednotlivých stavů, ale se shodnou topologií modelu přechodového diagramu. Pro popis preventivní údržby je zavedena matice údržby U (t1 , t2 ) ( t1 < t2 ), která stanovuje pravděpodobnosti přechodů mezi jednotlivými stavy v časovém K matematickému řešení úloh bude využito markovských řetězců.
intervalu
< t1 , t2 ) .
Další řešenou problematikou je modelování výkonových konfigurací systému - kapitola 6. Výkonovou dynamiku lze sledovat u systémů, kde se připojují resp. odpojují subsystémy nebo pracují na nižší výkon na základě okamžitého množství objemu výroby. Pro tyto systémy je možné sestavit větší množství modelů spolehlivosti, které jsou závislé na požadovaném objemu výkonu. Pro každý výkonový scénář se vytvoří přechodový diagram, který se nemusí odlišovat v topologii od ostatních, ale který se odlišuje především v intenzitách přechodů mezi stavy a v definici, zda je daný stav provozní, částečně poruchový, nebo poruchový. Celková pohotovost systému je potom funkcí objemu výroby a způsobem zálohování jednotlivých subsystémů. Dalším rozšířením modelu bude další aplikace, která umožní zjistit způsobilost systému - kapitola 7.1. Dále model umožní ohodnotit stavy a přechody mezi stavy (porucha, Strana: 14 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
obnova nebo údržba komponenty) ziskem nebo ztrátou - kapitola 7.2. Tím lze, krom funkce okamžité pohotovosti A(t ) , určit i funkci okamžité způsobilosti C (t ) . Dále pomocí ohodnocení přechodů a stavů lze stanovit funkci celkového zisku v zadaném časovém intervalu CostT (t1 , t 2 ) . Tato veličina může být jedním ze vstupů pro metodu nákladů životního cyklu LCC (analýza nákladů životního cyklu). Definováním určitých předpokládaných scénářů popisujících chování systému v čase lze zjistit celkový vliv vložených finančních prostředků na pohotovost a způsobilost systému. Tento prostředek je vhodný především pro průmyslovou aplikační činnost. Zatím se autor nesetkal v literatuře se spojením neexponenciálního rozdělení doby do poruchy vícestavových systémů, preventivní údržbou a popisem výkonových konfigurací, které jsou řešeny pomocí markovské analýzy, kromě prací napsané autorem ([89], [91], [92], [93], [100] a [112]). V aplikačním příkladě, který je součástí disertační práce, je vytvořen a vyhodnocen model spolehlivosti reálné kompresorové stanice a přilehlých trasových linií tranzitního plynovodu RWE Transgas, a.s. - příloha 4. Přepravní síť plynu, obdobně jako jiné přepravní sítě, v případě poruchy, výpadku, napadení útočníkem, představuje ohrožení pro velkou část populace. Následky v případě poruchy mohou být významné z pohledu zdraví, ekonomiky a životního prostředí. Kompresorová stanice tranzitního plynovodu byla vybrána, protože se jedná o zařízení, které je zálohované, komponenty jsou složitě udržovány a s časem z pohledu spolehlivosti degradují. Celý systém je navíc schopen přepravovat rozdílné objemy plynu. Počet turbokompresorů a vstupních/výstupních potrubních linií v provozu je funkcí právě přepravovaného objemu plynu. V dnešní době je v České republice, ale i v celé Evropské unii financováno mnoho projektů, které se zabývají problematikou spolehlivosti, bezpečností tranzitních sítí, ochranou systému před napadením, problematikou životního prostředí a krizového řízení. Z tohoto pohledu je řešené téma aktuální. Tato disertační práce vznikla díky finanční podpoře projektu AV ČR číslo 1ET401940412 s názvem Modelování a kvantifikace spolehlivosti dynamických systémů a projektu MŠMT Výzkumná centra (PP2-DP01) číslo 1M0554 s názvem Výzkumné centrum Pokročilé sanační technologie a procesy.
Strana: 15 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
3. Teoretické poznatky V kapitole 3 jsou uvedeny obecné poznatky, které jsou využity při tvorbě disertační práce. V kapitole 3.1 jsou představeny základy stochastických procesů. V kapitole je definován rozdíl mezi stochastickými a markovskými řetězci/procesy s uvedením základní markovské vlastnosti. Markovské řetězce/procesy jsou následně rozděleny na homogenní a nehomogenní. Této teorie je využito pro stanovení ukazatelů spolehlivosti systému s využitím metody markovské analýzy. Mostem mezi teorií stochastických procesů a ukazateli spolehlivosti je výpočet pravděpodobnosti stavů v čase pi (t ) a následné stanovení funkce okamžité nepohotovosti systému U (t ) =
∑ p (t ) . ΩD
i
V kapitole 3.2 jsou rozvedeny základní charakteristiky metod spolehlivosti, které literatura uvádí jako vhodné pro zjištění funkce okamžité nepohotovosti systému U (t ) . Metoda markovské analýzy je vzhledem k rozsahu využití v disertační práci popsána v kapitole 3.3. Tam jsou uvedeny i základní způsoby numerického řešení diferenciálních rovnic metodou Runge-Kutta a metodou Monte Carlo. V kapitole 3.4 jsou uvedeny způsoby zjištění bodových odhadů parametrů rozdělení pro dvoustavový systém, kde doba do poruchy je popsána buď exponenciálním, nebo Weibullovým rozdělením. Způsoby zjištění bodových odhadů parametrů pro jiná rozdělení jsou pouze naznačena. Weibullovo rozdělení se v praxi využívá pro popis období časných poruch a období degradace. V disertační práci je spojena markovská analýza spolehlivosti (popisující vícestavové systémy) a využití Weibullova rozdělení doby do přechodu. V kapitole 3.5 jsou uvedeny základní charakteristiky softwarových prostředků, které se zabývají markovskou analýzou spolehlivosti.
Strana: 16 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
3.1 Stochastické procesy [2], [3], [4], [5], [6] Stochastický proces - se nazývá každá funkce X (t ) , jejíž hodnota při dané hodnotě argumentu t (nabývajícího hodnot z určitého definičního oboru J ) je náhodnou veličinou. Obsahuje-li definiční obor J jen konečně nebo spočetně mnoho hodnot, mluví se o stochastickém procesu s diskrétním argumentem. V opačném případě o stochastickém procesu se spojitým argumentem t . Parametrem t se pro potřebu disertační práce, pokud není uvedeno jinak, rozumí čas. Stochastický proces lze chápat jako funkci dvou proměnných a to jednak argumentu (času t ), jednak elementárního jevu E , patřícího do úplného prostoru Ω elementárních jevů. Tedy X (t ) = Φ ( E , t ) , kde E ∈ Ω a t ∈ J . Jako příklad náhodné veličiny ve spolehlivosti lze uvést poruchu/obnovu komponenty. Průsek stochastického procesu - Pro každou pevnou hodnotu t = t0 ∈ J závisí Φ ( E0 , t ) jen na elementárním jevu E , takže jde o náhodnou veličinu. Ta se nazývá průsekem stochastického procesu X (t ) v bodě t0 a označuje se X (t0 ) . Realizace stochastického procesu - Při pevném elementárním jevu E = E0 závisí Φ ( E0 , t ) jen na proměnné t , proto představuje nenáhodnou funkci argumentu t . To jest Φ ( E0 , t ) = x(t ) . Každá taková funkce x(t ) se nazývá realizací daného stochastického procesu X (t ) . Realizaci x(t ) je možno považovat za výsledek provedeného pozorování jednoho a téhož stochastického procesu. Rozložení stochastického procesu - Pro určení rozložení stochastického procesu se vezme v úvahu ν realizací x1 (t ), x2 (t ),..., xν (t ) stochastického procesu X (t ) . Zvolí se t = t1 , pro které přechází stochastický proces v náhodnou veličinu X (t1 ) , zvanou průsek stochastického
procesu X (t ) v bodě t1 . Z uvažovaného počtu ν realizací se vybere právě ν 1 , jejichž
hodnoty pro t = t1 jsou nejvýše rovny zvolenému číslu x1 , které se nazvou hladinou x1 . Podle
ν1 . Tedy P ( X (t1 ) ≤ x1 ) značí pravděpodobnost, ν →∞ ν
zákona velkých čísel platí P ( X (t1 ) ≤ x1 ) = lim
že hodnota náhodné veličiny X (t1 ) nebude větší než x1 . Tato pravděpodobnost je funkcí dvou argumentů t1 a x1 . Je tedy F ( x1 , t1 ) = P ( X (t1 ) ≤ x1 ) . Funkce jednorozměrná distribuční funkce stochastického procesu X (t ) .
F
se nazývá
Střední hodnota stochastického procesu X (t0 ) - se nazývá náhodná funkce x(t ) , která se pro každou hodnotu času t z definičního oboru J rovná střední hodnotě odpovídajícího průseku stochastického procesu X (t ) v bodě t . Tedy x(t ) = E ( X (t )) . Stochastický proces může mít následující základní podoby: •
s diskrétními stavy a spojitým argumentem, Strana: 17 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
• • •
s diskrétními stavy a diskrétním argumentem, se spojitými stavy a spojitým argumentem, se spojitými stavy a diskrétním argumentem.
Disertační práce využívá pouze stochastické procesy s diskrétními stavy, kterých může být konečný, případně nejvýše spočetný počet stavů. 3.1.1
Markovské řetězce
Markovské řetězce - Řekneme, že posloupnost realizací pokusu nebo příslušná posloupnost náhodných veličin X n (pro n = 1,2,.. ) tvoří markovský řetězec, platí-li pro každý index
m = 2,3,.. a pro všechny možné hodnoty j , i,..., r těchto náhodných veličin vztah: P[ X m = j X m−1 = i ] = P[ X m = j X m−1 = i,..., X 1 = r ] = pij
( m)
(3.1)
Hodnoty j = 1,2,.. (a tedy také jevy E1 , E2 ,... ) se nazývají stavy uvažovaného systému. Podmíněná pravděpodobnost pij
(m)
se nazývá pravděpodobností přechodu daného
systému po m -té realizaci pokusu ze stavu i do stavu j . Příslušná matice se nazývá matice přechodu po m krocích ( m = 1,2,3,... ) a značí se Pm . Pro ∀i, j ∈ N platí:
P( m )
⎡ p11( m ) ⎢ (m) p = ⎢ 21 ⎢ ... ⎢ (m) ⎢⎣ pN 1
p12
(m)
...
(m)
p22 ... ... ... (m) pN 2 ...
⎤ (m) ⎥ p2 N ⎥ ... ⎥ (m) ⎥ pNN ⎥⎦ p1N
(m)
Homogenní markovský řetězec - Řekne se, že posloupnost
(3.2)
{X n }
náhodných veličin X n
tvoří homogenní markovský řetězec, jestliže podmíněná pravděpodobnost pij
(m)
nezávisí
na hodnotě m , tj. jestliže pro všechny hodnoty i, j , m platí:
pij
(m)
= pij
()
(3.3a)
Pro homogenní markovský řetězec platí, že matice pravděpodobností přechodů P je tvořena prvky pravděpodobností přechodů mezi stavy p ij pro ∀i, j ∈ N .
⎡ p11 ⎢p P = ⎢ 21 ⎢ ... ⎢ ⎣ p N1 Nechť
p12
...
p 22
...
...
...
pN 2
...
{X m }
p1N ⎤ p 2 N ⎥⎥ ... ⎥ ⎥ p NN ⎦
je homogenní markovský řetězec a nechť pij
(3.3b)
(m)
je pravděpodobnost
přechodu za m kroků ze stavu i do stavu j . Pak platí že
Pm = P1
m
(3.4). Strana: 18 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Stochastická matice - matice pravděpodobností přechodů P , jejíž prvky p ij ( ∀i, j ∈ N ) jsou nezáporné a součet prvků v každém řádku se rovná jedné, se nazývá stochastická. Nebo-li
∀i, j ∈ N platí, že 0 ≤ pij ≤ 1 a
N
∑p j =1
ij
=1
Matice pravděpodobností přechodů P , která je využívána v disertační práci, je stochastická matice. Krom pravděpodobností přechodu pij je třeba znát také počáteční stav, v němž se systém nachází na začátku procesu. Tento stav může být náhodný, a potom je třeba znát příslušnou pravděpodobnost náhodné veličiny X 0 . Nebo může být pevně určen, potom příslušný stav lze vyjádřit pravděpodobností P ( X 0 = k ) = 1 . Tyto počáteční pravděpodobnosti jsou prvky pravděpodobnostního vektoru p(0) . Chapman-Kolmogorova rovnice - Nechť X n je homogenní markovský řetězec a nechť pijm je pravděpodobnost přechodu za m kroků ze stavu i do stavu j . Nechť je dále m1 + m2 = m , kde m1 , m2 , m jsou celá čísla. Pak platí vztah:
pijm1 + m2 = ∑ pikm1 p kjm2
(3.5)
k
Pravděpodobnosti stavů homogenních markovských řetězců po m krocích pi (m) lze zjistit na základě znalosti prvků pravděpodobnostního vektoru p(m) z počátečního pravděpodobnostního vektoru p(0) a matice pravděpodobnosti přechodů P pomocí vztahu (3.6):
p(m) = p(0)P m
(3.6)
Nehomogenní markovský řetězec - Řekne se, že posloupnost {X n } náhodných veličin X n tvoří nehomogenní markovský řetězec, jestliže podmíněná pravděpodobnost pij
(m )
závisí
na hodnotě m . Pravděpodobnosti stavů nehomogenních markovských řetězců po m krocích lze zjistit p(m) z počátečního na základě znalosti prvků pravděpodobnostního vektoru pravděpodobnostního vektoru p(0) a matice pravděpodobnosti přechodů P pomocí vztahu (3.7):
p(m) = p(0) ⋅ P(0) ⋅ P(1) ⋅ ... ⋅ P(m) 3.1.2
(3.7)
Markovské procesy
Markovský proces - Stochastický proces X (t ) se nazve markovským procesem, právě když pro n = 1,2,.. a pro libovolné hodnoty t m ( m = 0,1,2,..., n ), pro něž platí nerovnosti
t0 < t1 < ... < tn , platí pro všechny reálná čísla x a y vztah: P[ X (t n ) < y X (tn−1 ) = x] = P[ X (tn ) < y X (t n−1 ) = x, X (t n−2 ) = xn−2 ,..., X (t0 ) = x0 ] , (3.8) kde x0 , x1 ,..., xn − 2 jsou libovolná reálná čísla. Strana: 19 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Zjednodušeně jsou markovské procesy takové stochastické procesy, u nichž v každém okamžiku τ podmíněné rozložení pravděpodobnosti v budoucím okamžiku t ( t > τ ) závisí jen na hodnotě veličiny X (τ ) a nezávisí na hodnotách nabytých v okamžicích t < τ . U obecného stochastického procesu závisí příslušná funkce podmíněného rozložení stochastického procesu X (τ ) i na hodnotách nabytých v okamžicích t < τ .
Podmíněná distribuční funkce - Nechť X (t ) je markovský proces. Podmíněnou distribuční funkci ve tvaru F (t 2 , y t1 , x) = P( X (t 2 ) < y X (t1 ) = x) (3.9)
X (t 2 ) < y za předpokladu, že ( t1 < t 2 ) náhodná veličina X (t1 ) nabyla hodnoty x . Tato
vyjadřující pravděpodobnost, že v okamžiku t2
bude
v předchozím okamžiku t1 podmíněná distribuční funkce se nazývá funkce přechodu.
Homogenní markovský proces - Závisí-li funkce přechodu F (t 2 , y t1 , x) jen na délce intervalu t 2 − t1 , nikoli na okamžicích t1 a t2 , pak se příslušný markovský proces nazývá homogenní vzhledem k času t . Distribuční funkce se pak značí F (t , x, y ) . Diskrétní markovský proces - se nazývá diskrétní proces X (t ) , jestliže pro n = 1,2,... , pro libovolné okamžiky t m ( m = 0,1,2,..., n ), přičemž t0 < t1 < ... < t n , a pro libovolná celá čísla
i, j , i0 , i1 ,..., in−2 platí vztah: P ( X (t n ) = j X (t n−1 ) =i, X (t n−2 ) = in−2 ,..., X (t0 ) = i0 ) = P( X (t n ) = j X (t n−1 ) =i ) . (3.10a) Funkce ve tvaru:
pij (t1 , t 2 ) = P( X (t 2 ) = j X (t1 ) =i )
(3.10b),
definovanou pro libovolné okamžiky t1 , t 2 ( t1 < t 2 ) a pro libovolná celá čísla i , j , se nazývá pravděpodobnost přechodu ze stavu i do stavu j . Je-li diskrétní markovský proces homogenní, závisí pravděpodobnost přechodu ze stavu i do stavu j , tj. pij (t1 , t 2 ) , pouze na rozdílu t = t 2 − t1 a lze napsat:
pij (t1 , t 2 ) = pij (t ) Chapman-Kolmogorovy rovnice lze zobecnit pro diskrétní markovský proces, za předpokladu, že t1 ,t 2 jsou libovolné okamžiky a dále τ ∈ (t1 , t 2 ) , na rovnici:
pij (t1 , t 2 ) = ∑ pik (t1 ,τ ) pkj (τ , t 2 ) .
(3.11a)
k
Pro diskrétní markovský proces homogenní jsou Chapman-Kolmogorovy rovnice ve tvaru:
Strana: 20 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
pij (t + τ ) = ∑ pik (t ) pkj (τ )
(3.11b)
k
Intenzita přechodu markovského procesu - Nechť X (t ) je diskrétní markovský proces, který má nejvýše spočetně mnoho stavů i = 0,1,2,... , a nechť pro 0 < τ < t značí pij (τ , t ) funkci pro pravděpodobnost přechodu ze stavu i do stavu j . Funkce hii (t ) a hij (t ) jsou definované vztahem:
lim
pij (t , t + ∆t )
∆t
∆t →0 i≠ j
lim
= hij (t ) ≥ 0
p jj (t , t + ∆t ) − 1
∆t
∆t →0 i= j
= hii (t ) ≤ 0
(3.12a)
(3.12b)
a jsou nazvány intenzitou přechodu mezi stavy. Intenzity přechodu mezi stavy hii (t ) a hij (t ) ( i ≠ j ) jsou prvky matice intenzit přechodů h(t ) . Platí, že
∑ h (t ) = 0 ij
j
Chapman-Kolmogorovy diferenciální rovnice - Nechť pij (τ , t ) jsou pravděpodobnosti přechodů v markovském procesu X (t ) s nejvýše konečným počtem stavů i = 0,1,2,... . Nechť funkce hii (t ) a hij (t ) ( i ≠ j ) tohoto procesu existují a jsou spojité. Potom funkce pij (τ , t ) splňuje následující soustavu diferenciálních rovnic:
∂pij (τ , t ) ∂t
= ∑ pik (τ , t )hkj (t ) ,
(3.13a)
k
při počátečních podmínkách pii (τ ,τ ) = 1 , pij (τ ,τ ) = 0 ( i ≠ j ). Za předpokladu, že uvažovaný markovský proces X (t ) je homogenní, jsou funkce intenzity přechodu nezávislé na času t a lze je zapsat ve tvaru hii a hij ( i ≠ j ). Parciální diferenciální rovnice přejdou v obyčejné diferenciální rovnice ve tvaru:
dpij (t ) dt
= ∑ pik (τ , t )hkj
(3.13b).
k
Při počátečních podmínkách pii (0) = 1 , pij (0) = 0 ( i ≠ j ). Z intenzit přechodu markovského procesu lze zjistit matici pravděpodobnosti přechodů
P (t1 , t2 ) za předpokladu t2 − t1 = ∆t → 0 podle vzorce:
Strana: 21 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
⎡1 ⎢0 P(t1 , t2 ) = I + h(t1 ) ⋅ ∆t = ⎢ ⎢... ⎢ ⎣0
0 1 ... 0
... ... ... ...
⎡ N h12 (t1 ) ⎢− ∑ h1i (t1 ) i =1, i ≠1 0⎤ ⎢ N ⎢ h (t ) 0 ⎥⎥ h2i (t1 ) − ∑ 21 1 + ∆t ⎢ i = 1 , i ≠ 2 ...⎥ ⎢ ... ... ⎥ ⎢ 1⎦ ⎢ hN 1 (t1 ) hN 2 (t1 ) ⎢⎣
⎤ ⎥ ⎥ ... h2 N (t1 ) ⎥ ⎥ ⎥ ... ... ⎥ N ... − ∑ hNN (t1 )⎥ ⎥⎦ i =1, i ≠ N ...
h1N (t1 )
(3.14). Pravděpodobnosti stavů markovských procesů v čase t lze stanovit na základě znalosti počátečních podmínek p(0) a matice intenzit přechodů h(t ) popisujících systém pomocí vztahu:
dp(t ) = p(t )h dt
pro homogenní markovský proces,
(3.15a),
dp(t ) = p(t )h(t ) dt
pro nehomogenní markovský proces
(3.15b).
Strana: 22 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
3.2 Základní metody spolehlivosti vhodné k zjištění pohotovosti 3.2.1
Základní definice V disertační práci bude použito termínů, které jsou definovány v normách ČSN zabývajících se spolehlivostí. Především se jedná o texty [67], [68], [74] a [75]. Jedním z dílčích cílů autora disertační práce je aplikovat teoretické poznatky do praktického příkladu modelování kompresorové stanice tranzitního plynovodu. Tento příklad je uveden v příloze 4. Normy ČSN jsou často používány jako dorozumívací jazyk v praxi, proto budou v textu uvedeny odkazy, mimo jinou literaturu, i na ně. V názvu disertační práce je uveden termín „spolehlivost“, kterou lze definovat: Spolehlivost (cit. [67] def. 2-3) - je souhrnný termín používaný pro popis pohotovosti a činitelů, které ji ovlivňují: bezporuchovost, udržovatelnost a zajištěnost údržby. Pohotovost (cit. [67] def. 2-5) - je schopnost objektu být ve stavu schopném plnit požadovanou funkci v daných podmínkách, v daném časovém okamžiku nebo v daném časovém intervalu, za předpokladu, že jsou zajištěny požadované vnější prostředky. Bezporuchovost (cit. [67] def. 2-6) - je schopnost objektu plnit požadovanou funkci v daných podmínkách a v daném časovém intervalu. Udržovatelnost (cit. [67] def. 2-7) - je schopnost objektu v daných podmínkách používání setrvat ve stavu nebo vrátit se do stavu, v němž může plnit požadovanou funkci, jestliže se údržba provádí v daných podmínkách a používají se stanovené postupy a prostředky. Zajištěnost údržby (cit. [67] def. 2-8) - schopnost organizace poskytující údržbářské služby zajišťovat podle požadavků v daných podmínkách prostředky potřebné pro údržbu podle dané koncepce údržby. 3.2.2
Základní kvantitativní ukazatele spolehlivosti Metody analýzy spolehlivosti se používají k stanovení ukazatelů bezporuchovosti, pohotovosti, udržovatelnosti a bezpečnosti systému. Analýzy spolehlivosti se provádějí ve všech etapách života objektu. Aplikují se na různých úrovních a stupních rozčlenění systému pro vyhodnocení a stanovení ukazatelů spolehlivosti systému nebo investičního celku. Používají se též pro porovnání výsledků analýzy se specifikovanými požadavky. Analytické metody umožňují vyhodnotit kvalitativní a kvantitativní ukazatele spolehlivosti a činitelů, které ji ovlivňují. Lze odhadnout hodnoty ukazatelů, které jsou definovány v [67]: • • •
ukazatele pohotovosti v kapitole 11, ukazatele bezporuchovosti v kapitole 12, ukazatele udržovatelnosti a zajištěnosti údržby v kapitole 13.
Způsoby výpočtu a metody zjištění jednotlivých ukazatelů lze nalézt například v [7], [8], [9], [10], [11], [12], [13], [101] a [107] a jsou také uvedeny na straně 8 disertační práce. 3.2.3
Základní kvantitativní analýzy spolehlivosti Pro zjištění a ověření funkce okamžité pohotovosti A(t ) systému se dnes v praxi nejčastěji využívá následujících základních metod spolehlivosti. Jejich seznam je uveden například v [12], [13], [73]. Odkazy na typy softwaru užívaného ve spolehlivosti lze najít v [14]. Strana: 23 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Mezi základní kvantitativní analýzy spolehlivosti lze zařadit například: • • • • • • •
předpověď intenzity poruch, [15], [16], [17], [18], [19], [79], [114] pravdivostní tabulka - analýza funkční struktury, [70] RBD - analýza blokového diagramu bezporuchovosti, [20], [69], [70] FTA - analýza stromu poruchových stavů, [21], [71], [72] ETA - analýza stromu událostí, [22], [23] MA - markovská analýza, 6 PN - analýza Petriho sítí,
Žádná jednotlivá metoda analýzy spolehlivosti není dostatečně vyčerpávající a pružná, aby se vypořádala s možnými složitostmi modelu požadovanými k vyhodnocení význačných rysů praktických systémů. Aby bylo zajištěno řádné zpracování složitých nebo multifunkčních systémů, může být nezbytné uvážit použití několika vzájemně se doplňujících metod analýzy.
6
Základy metody jsou vysvětleny v kapitole 3.3, kde jsou uvedeny i odkazy na literaturu
Strana: 24 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
3.3 Markovská analýza ve spolehlivosti [24], [25], [26], [27], [28], [29], [74] a [75] Tato kapitola popisuje metodu markovské analýzy, což je jedna z metod analýz spolehlivosti, jejichž seznam je uveden v kapitole 3.2.3. Protože tato metoda je využívána v celé disertační práci, je jí věnována celá kapitola. Metoda využívá matematických poznatků z kapitoly 3.1. 3.3.1
Podstata metody Markovské modely ve spolehlivosti představují metodu, která umožňuje řešit dynamickou závislost charakteristik poruchy či obnovy jednotlivých součástek a přizpůsobit je stavům přechodového diagramu systému. Markovskými modely lze zachytit vlivy jak poruch komponent závislých na pořadí, tak změny intenzit přechodů vyplývající z namáhání či jiných faktorů. Z tohoto důvodu je markovská analýza vhodná metoda pro hodnocení spolehlivosti funkčně složitých konstrukcí systému a složitých strategií oprav a údržby. Markovská metoda je založena na teorii markovských procesů (kapitola 3.1), kde argument je obvykle čas. Pro výpočet ukazatelů spolehlivosti jsou obvykle využívány homogenní markovské procesy, které vyžadují, aby byly intenzity přechodů mezi stavy konstantní. K reprezentaci chování systému pomocí markovské analýzy je nutné stanovit všechny možné stavy systému, znázorněné graficky v přechodovém diagramu. Každému stavu je určeno, zda se jedná o provoz, částečně poruchový stav nebo poruchový stav. Přechody mezi stavy v přechodovém diagramu jsou ohodnoceny intenzitami poruch/oprav. Následně se pro zjištění výsledné hodnoty parametrů pohotovosti nebo bezporuchovosti převede přechodový diagram pomocí vzorců uvedených v kapitole 3.1 (vzorec 3.14) na matematický model. Vyhodnocují se pravděpodobnosti, kdy systém je v čase t v jednotlivých stavech přechodového diagramu. Typickými výstupy markovského modelu jsou pravděpodobnosti, s jakými se systém nachází v daném čase a v daném stavu. Příkladem zpracovaného výsledku pravděpodobností stavů je ukazatel funkce okamžité pohotovosti A(t ) , kdy se sčítají pravděpodobnosti stavů reprezentující provoz a částečnou poruchovost v čase:
A(t ) =
N
∑ p (t ) .
i =1,ΩU
i
(3.16)
Z použití této metodiky plynou následující výhody. Metoda poskytuje pružný pravděpodobnostní model pro analýzu chování systému. Je možné ji přizpůsobit pro složité redundantní konfigurace, pro složitou koncepci údržby, složité modely ošetření poruchových stavů, degradované režimy provozu a poruchy se společnou příčinou. Metoda poskytuje pravděpodobnostní řešení pro moduly, které se mají vložit do jiných modelů, jako jsou blokové diagramy bezporuchovosti - RBD a stromy poruchových stavů - FTA. Metoda umožňuje přesné modelování posloupností událostí se specifickým typem nebo pořadím výskytu. Metoda však má i jistá omezení. Zvyšujícím se počtem komponent systému exponenciálně roste počet stavů, což vede ke zvýšení pracnosti analýzy. Pro uživatele může být obtížné model sestavit a ověřit jeho správnost.
Strana: 25 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
3.3.2
Sestavení markovského přechodového diagramu a matice intenzit Markovská analýza využívá přechodový diagram mezi stavy, který je grafickou reprezentací vlastností spolehlivosti systému. Tento diagram modeluje spolehlivostní hlediska chování systému v čase. Na systém se pohlíží jako na určitý počet jednotek, z nichž každá může existovat pouze v jednom ze dvou stavů: provozním nebo poruchovém. Systém jako celek však může existovat v mnoha různých stavech, z nichž každý je určen danou kombinací nefunkčních a funkčních jednotek. Definují se tyto základní stavy systému: provozní, částečně poruchový stav, poruchový stav. V některé literatuře jsou rozdělovány stavy na funkční, degradované a nefunkční.
Počet stavů reprezentujících provoz (funkční stav), částečně poruchových (degradovaných) a poruchových (nefunkčních), může být v systému větší než jeden. Hlavní důvod spočívá v tom, že počet sledů poruch, které vedou k poruše celého systému je obecně také více. Pokud u některé jednotky systému dojde k poruše, nebo na některé jednotce dojde k obnově, „přechází“ systém v přechodovém diagramu z jednoho stavu do následujícího stavu. Tento druh modelu se obecně nazývá model s diskrétními stavy a spojitým argumentem. Na počátku analýzy musí být definovány stavy, které jsou předmětem zájmu analýzy spolu s intenzitami přechodu z jednoho stavu do jiného (intenzity poruch/oprav). Při modelování je nutné splnit podmínku, že namodelovaný systém je bez paměti - vzorec (3.8). Tato podmínka znamená, že budoucí chování systému závisí pouze na přítomném stavu, a ne na minulosti. Tato podmínka je splněna, jestliže: • •
intenzity přechodů mezi stavy jsou časově konstantní - homogenní markovský proces, intenzity přechodů mezi stavy nejsou časově konstantní. Platí však, že intenzity přechodů v čase tn jsou spjaty k pevnému časovému okamžiku tn −1 ( tn > tn −1 ) a nikoliv k časovým okamžikům předcházejícím tn− 2 ,...,t0 ( t0 < t1 < ... < tn − 3 < tn − 2 < tn −1 ) nehomogenní markovský proces (podmínka vzorce 3.8).
Kvalitativní analýza vyžaduje určení všech možných stavů systému, znázorněných v podobě přechodového diagramu. Kvantitativní analýza využívá intenzity poruch/oprav a vztahy mezi stavy znázorněnými přechodovým diagramem, které umožňují sestavit žádanou matici intenzit přechodů h . Matice je základem pro matematický model, popsaný v rovnici (3.15a), k zjištění ukazatelů pohotovosti (rovnice 3.16) nebo bezporuchovosti systému. Obecně, pro markovské procesy se spojitým časem, prvek hij čtvercové matice h představuje intenzitu přechodu (intenzitu poruch/oprav) ze stavu i do stavu j . Přechody mezi stavy lze interpretovat jako poruchu resp. obnovu komponent systému. Jestliže číselné označení stavů se zvětšuje se zvětšující se degradací systému, potom prvky nad diagonálou představují intenzity poruch, pod diagonálou jsou intenzity oprav. Hodnoty nediagonálních prvků matice intenzit přechodů h jsou známy z přípravných analýz (přechody mezi stavy) a jsou vždy nezáporné. Jednotlivé intenzity přechodů v matici h lze interpretovat: • •
nenulový nediagonální prvek matice h označuje intenzitu přechodu ze stavu i do stavu j, nulový nediagonální prvek matice h označuje nemožnost okamžitého přechodu ze stavu i do stavu j, Strana: 26 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
•
nulový diagonální prvek označuje absorpční stav, tzn. pravděpodobnost setrvání v daném stavu je jednotková.
Diagonální prvky matice intenzit přechodů h se dopočítávají tak (rovnice 3.12), že součet prvků v každém řádku je 0, tzn. že diagonální prvky jsou vždy nekladné. Obdobně jako pro markovské procesy lze definovat pro markovské řetězce pravděpodobnosti přechodů pij , které určují pravděpodobnost přechodu mezi stavy i a j . Pravděpodobnosti přechodů pij mezi stavy i a j se zapisují do matice pravděpodobnosti přechodů P . Viz vzorec (3.3b). 3.3.3
Vyhodnocení markovského přechodového diagramu Pomocí matice intenzit přechodů h lze řešit soustavu diferenciálních rovnic
dp(t ) = p(t )h 7 , která je použitelná pro markovské homogenní procesy. Rovnice dt dp(t ) (3.15b) = p(t )h(t ) - je určena pro markovské nehomogenní procesy. dt Rovnice (3.15) představuje soustavu lineárních diferenciálních rovnic prvního řádu. V počátečních podmínkách (vektor p(0) ) se většinou předpokládá, že pravděpodobnost plně funkčního stavu je 1, všechny ostatní stavy jsou rovny 0. Pro pravděpodobnostní vektor p(t ) platí, že ∀i ∈ {1, K , N }, pi ≥ 0 ,
∑p
i
= 1.
i
Pro každý prvek ∀i, j ∈ {1,K, N } matice intenzit přechodů h(t ) je dáno, že hij (t ) ≥ 0 ( i ≠ j ) a h jj (t ) ≤ 0 . Dále pro ∀i platí součet
∑ h (t ) = 0 . ij
j
Řešením soustavy diferenciálních rovnic se získá vektor průběhu pravděpodobností stavů v čase. Výsledná funkce okamžité pohotovosti systému A(t ) se získá součtem pravděpodobností stavů přes všechny provozní a částečně poruchové stavy A(t ) =
N
∑ p (t ) .
i =1,ΩU
i
Funkce okamžité pohotovosti A(t ) se ustálí na asymptotické hodnotě pohotovosti A tehdy, když matice pravděpodobností/intenzit přechodů má konstantní prvky po celou dobu simulace a celá matice je regulární. Diferenciální rovnice se obvykle řeší pomocí numerických metod, mezi kterými lze vyzdvihnout metody Runge-Kutta a metodu Monte Carlo. Softwarové prostředky obvykle využívají metody Runge-Kutta. Obtížnost řešení úlohy spočívá v řešení soustavy diferenciálních rovnic. Splněna také
7
d [ p1 (t ) L dt
pi (t )] = [ p1 (t ) L
⎡h11 L h1i ⎤ pi (t )] ⋅ ⎢⎢ M O M ⎥⎥ ⎢⎣ hi1 L hii ⎥⎦ Strana: 27 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
musí být podmínka (základní markovská vlastnost rovnice 3.8), že je systém bez paměti. To například znamená, že intenzity přechodů mezi jednotlivými stavy jsou v čase konstantní. Konstantní intenzita přechodu mezi 2 stavy je splněna, když je doba do přechodu popsána exponenciálním rozdělením. 3.3.4
Řešení obyčejných diferenciálních rovnic metodou Runge-Kutta [5], [30], [31], [32]
Obyčejné diferenciální rovnice s počátečními podmínkami lze řešit několika metodami, které lze rozlišit na jednokrokové a k-krokové metody. U jednokrokových metod výpočet řešení diferenciální rovnice probíhá tak, že po provedení jednoho kroku se vlastně řeší nový problém s počáteční podmínkou zadanou v tom bodě, jehož řešení je známo z předešlého kroku. Jednokroková metoda má své výhody, protože se jednoduše programuje a dále není složitý přechod k jiné délce integračního kroku během výpočtu. Na druhou stranu má metoda i své nevýhody. Nejzávažnější je pravděpodobně ta, že lze velice obtížně odhadnout lokální chybu. U k-krokové metody se pro řešení diferenciální rovnice v bodě ym + k využívají funkční hodnoty k předcházejících bodů ym , ym +1 ,..., ym + k −1 . V disertační práci je pro výpočet funkce okamžité nepohotovosti U (t ) systému použita jednokroková metoda. Nejjednodušší jednokroková metoda je Eulerova metoda. V této metodě se rozdělí
t2 − t1 . Dělícími body jsou ti = t1 + i ⋅ l , i = 1,..., m . m Přibližná hodnota řešení daná pravděpodobnostním vektorem p(ti ) je vypočtena z rekurence interval t ∈ t1 ,t2
na m dílků délky l =
dané vzorcem (3.17)
p(ti +1 ) = p(ti ) + l ⋅
dp(ti ) , i = 1,..., m − 1 . dt
(3.17)
Počáteční podmínka je předpokládána ve tvaru p(0) . V disertační práci jsou výpočty provedeny pomocí Eulerovy metody, protože při modelování preventivní údržby jsou pravděpodobnosti jednotlivých stavů nespojité funkce. Dále se v disertační práci často využívají přechody mezi markovskými procesy a diskrétními markovskými procesy například při modelování preventivní údržby či přechody mezi výkonovými konfiguracemi systému - viz kapitoly 5 a 6. Výpočet pomocí Eulerovy metody je snadněji interpretovatelný než při použití metody Runge-Kutta např. čtvrtého řádu. 3.3.5
Metoda Monte Carlo [28], [29], [33], [34] a [35]
Metoda Monte Carlo je založena na mnohonásobném opakování simulačních realizací. Každá simulační realizace se sestává z generování náhodného pokusu, který popisuje chování systému v čase. Ve výsledku se sleduje rozložení stochastického procesu . Každá realizace startuje z předem daného stavu, který je dán počátečními podmínkami simulace. Posléze se pomocí vygenerování náhodných čísel z příslušné předepsané distribuční funkce určí doba do přechodu ze stavu i do jednoho z dalších stavů, kam umožňuje přejít systému přechodový diagram. Toto se opakuje pro všechny možné přechody ze stavu i . Časový okamžik přechodu a stav do kterého systém přejde ze stavu i , Strana: 28 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
je dán minimalizací časových okamžiků přes všechny možné přechody ze stavu i . Prvním krokem analýzy je tvorba vstupních souborů a řídícího programu na základě generování události podle požadované distribuční funkce. Tato funkce může být stanovena postupem nezávislým na výpočetním programu nebo je někdy součástí komplexnějších modelových nástrojů. Dalším krokem je aplikování takto získaných parametrů na celkový systémový model včetně uvažovaných podprogramů. Tento cyklus se opakuje tak dlouho, až je vyplněno kritérium maximálního počtu událostí, eventuelně kritérium přesnosti statistického výběru, které je průběžně kontrolováno řídícím programem. Posledním krokem pravděpodobnostní analýzy je statistické vyhodnocení výsledků. K řešení analýz pomocí metody Monte Carlo lze využít matematického softwaru například MATLAB® - http://www.mathworks.com/, nebo Goldsim - http://www.goldsim.com/.
.
Strana: 29 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
3.4 Statistická rozdělení popisující dobu do poruchy komponenty Doba do poruchy komponent se nejčastěji popisuje pomocí exponenciálního rozdělení. Jeho výhodou oproti jiným rozdělením je, že doba do poruchy se nemění v závislosti na době používání komponenty, neboli intenzita poruch komponent je konstantní. Tento předpoklad platí u komponent, které nejsou namáhány například silou, momentem nebo elektricky [15], [16], [19]. Další příklady jsou uvedeny v kapitole 3.4.1. V případě, že je shromážděno dostatek dat o době do poruchy a existují statisticky zjištěné předpoklady, že data o době do poruchy nejsou z exponenciálního rozdělení, používá se v technické praxi nejčastěji Weibullovo rozdělení. Pomocí Weibullova rozdělení se popisuje doba do poruchy namáhaných komponent. Mezi příklady z jednotlivých odvětví lze zařadit: • • •
strojní komponenty - ložiska, automobily, stavební prvky - mosty, elektromagnetické komponenty - relé.
Tato kapitola se zabývá způsoby zjištění bodových odhadů parametrů Weibullova rozdělení, které se používá pro popis doby do poruchy komponenty. Vzorce pro výpočet parametrů jsou uváděny v příslušných kapitolách pro: • •
censorovaná i necensorovaná data, komponenty, které jsou po poruše obnoveny výměnou/opravou.
Poznatků uvedených v této kapitole bude využito v disertační práci při modelování doby do přechodu (v přechodovém diagramu) s využitím neexponenciálního rozdělení. Stanovení relevantních bodových odhadů intenzit přechodů je důležitým základem pro modely spolehlivosti a následné technické využití. 3.4.1
Vanová křivka [1], [8], [12], [13], [39]
„Vanová křivka“ je zidealizovaná křivka, u které se na vodorovné ose vynáší doba provozu komponenty a na svislou osu její okamžitá intenzita poruch λ (t ) . Podle charakteru změny okamžité intenzity poruch komponenty lze křivku rozdělit do tří období. V prvním období se okamžitá intenzita poruch λ (t ) komponenty snižuje. Toto období se nazývá období časných poruch. Příčiny zvýšené intenzity poruch na počátku provozu jsou poruchy v důsledku výrobních poruch, nesprávné montáže, skrytých poruch, chyb při návrhu a při výrobě apod. Dobu do poruchy komponenty lze popsat v tomto období pomocí Weibullova rozdělení. V druhém období je okamžitá intenzita poruch λ (t ) komponenty přibližně konstantní. Toto období se nazývá období normálního života. Dochází k běžnému využití komponenty, k poruchám dochází náhodně. Nedochází k opotřebení, které by změnilo funkční vlastnosti komponenty. Dobu do poruchy komponenty lze popsat v tomto období pomocí exponenciálního rozdělení. V závěrečném třetím období se okamžitá intenzita poruch λ (t ) komponenty výrazně zvyšuje. Toto období se nazývá období degradace. V komponentě dochází k opotřebení. Procesy stárnutí a opotřebení mění funkční vlastnosti komponenty. Dobu do poruchy komponenty lze popsat v tomto období pomocí Weibullova rozdělení. Strana: 30 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
U mnoha komponent se nevyskytuje (často nelze prokázat) první období, kdy se okamžitá intenzita poruch snižuje (období časných poruch). Je to například u komponent, které jsou: • • •
zaběhnuty u výrobce, vnitřně zálohovány, jsou vysoce spolehlivé a není tak dostatek dat o poruchách, ze kterých lze snižování okamžité intenzity poruch prokázat.
Pro statistické prokázání období časných poruch komponenty je obvykle třeba získat velké množství dat o poruchách. U mnoha komponent se nevyskytuje třetí období, kdy se okamžitá intenzita poruch zvyšuje (tj. období degradace). Je to především u těch, které jsou vyřazeny dříve než se začnou projevovat degradační vlivy. Dále je to u komponent, které nejsou namáhány. Jako příklad lze uvést většinu elektronických komponent. Graficky je zobrazen zidealizovaný tvar vanové křivky na obr. 1.
intenzita poruch
λ(t) [hod]
Typická vanová křivka objektu
doba provozu
t [hod]
Obr. 1: Vanová křivka 3.4.2
Obnovované komponenty opravou / výměnou Pro zjišťování parametrů rozdělení při sledování doby do poruchy je nutné rozlišovat následující případy: • •
Komponenta se po poruše skutečně opravuje (opravovaný objekt). Doba do poruchy je popsána pomocí parametru proudu poruch z (t ) . Komponenta se po poruše neopravuje (neopravovaný objekt). Doba do poruchy je popsána pomocí intenzity poruch λ (t ) . Strana: 31 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Protože rozdíl v termínech „okamžitá intenzita poruch“ a „parametr proudu poruch“ se ne vždy důsledně rozlišuje, jsou zde uvedeny jejich definice. V disertační práci bude termín „okamžitá intenzita poruch“ využit pro neopravované komponenty, nebo pro opravované komponenty které jsou použity pouze do první poruchy. Termín „parametr proudu poruch“ bude použit pro opravované komponenty. Termín „intenzita přechodu mezi stavy i a j “ (zkráceně i „intenzita přechodu“) s označením hij (t ) představuje prvek matice intenzit přechodů h(t ) . Tento termín slučuje, pro potřebu disertační práce, termíny okamžitá intenzita poruch λij (t ) a okamžité intenzity oprav
µij (t ) pro neopravované objekty (nebo opravované komponenty používané do první poruchy). A dále slučuje termíny parametr proud poruch z (t ) a okamžité intenzity oprav µij (t ) pro opravované objekty. Okamžitá intenzita poruch (cit. [67] def. 12-02) - λ (t ) - limita poměru podmíněné pravděpodobnosti, existuje-li, že časový okamžik T vzniku poruchy objektu leží v daném časovém intervalu (t , t + ∆t ) k délce časového intervalu ∆t , jestliže ∆t se blíží nule, za podmínky, že na začátku časového intervalu je objekt v použitelném stavu Parametr proudu poruch (cit. [67] def. 12-04) - z (t ) - limita poměru, existuje-li, středního počtu poruch opravovaného objektu v časovém intervalu (t , t + ∆t ) , k délce tohoto intervalu ∆t , jestliže se délka časového intervalu blíží nule
E [r (t + ∆t ) − r (t )] ∆t → 0 ∆t
z (t ) = lim
V kapitolách 3.4.4, 3.4.5, 3.4.6 jsou uvedeny způsoby zjištění parametrů exponenciálního a Weibullova rozdělení, jestliže je/není po poruše komponenta opravována. O komponentě, která se po poruše skutečně opravuje, lze usuzovat, že má velmi podobné parametry rozdělení doby do poruchy jako shodná komponenta, u které nedošlo k poruše. V anglicky psané literatuře se tato vlastnost nazývá termínem „as good as old“. V česky psané literatuře neexistuje jednotný překlad. V této práci bude použito termínu „stejně spolehlivý jako starý“. O komponentě, která se po poruše neopravuje, lze usuzovat, že má velmi podobné parametry rozdělení doby do poruchy jako shodná komponenta, která je zcela nová. V anglicky psané literatuře se tato vlastnost nazývá termínem „as good as new“. V této práci bude použito termínu „stejně spolehlivý jako nový“. 3.4.3
Způsob získání neznámých parametrů předpokládaných rozdělení [1], [30]
Bodové odhady neznámých parametrů rozdělení z dat popisujících dobu do poruchy lze získat: • •
početní metodou - maximální věrohodností, momentů, graficky - pravděpodobnostní papír.
Strana: 32 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Metoda maximální věrohodnosti - Metoda maximální věrohodnosti je založena na maximalizaci věrohodnostní funkce. Věrohodnostní funkce je sdružená hustota pravděpodobnosti daného náhodného výběru a je brána jako funkce neznámých parametrů. Nechť X 1 , X 2 ,..., X n je náhodný výběr popisující dobu do poruchy z rozdělení s distribuční funkcí F ( x; Θ) . Tvar distribuční funkce je znám a Θ je/jsou neznámé parametr/parametry. Problém úlohy spočívá v nalezení takové statistiky, která bude funkcí X 1 , X 2 ,..., X n a bude použita jako odhad Θ .
f ( x; Θ) je hustota pravděpodobnosti vyšetřované funkce Θ , nazveme ji věrohodnostní funkcí založenou na x( x1 , x2 ,..., xn ) a označenou jako L(x; Θ) . Věrohodnostní Nechť
funkci obdržíme pomocí vzorce 3.19.
L( x1 , x2 ,..., x n ; Θ) = f1 ( x1; Θ) ⋅ ... ⋅ f n ( xn ; Θ)
(3.19)
Pro jednodušší získání maxima věrohodnostní funkce lze využít toho, že se logaritmováním maximum funkce nezmění.
ln L( x1 , x2 ,..., x n ; Θ) = ln[ f1 ( x1; Θ) ⋅ ... ⋅ f n ( xn ; Θ)] Maximum této funkce se získá derivováním podle neznámých parametrů Θ . První derivace se položí rovna 0 a druhá derivace musí být záporná. Aplikační příklad je uveden v poznámce 8. Vzorec 3.19 je uveden pro necensorovaná data. Pro data censorovaná časem je věrohodnostní funkce ve tvaru (3.20a):
L( x1 , x2 ,..., x n ; Θ) =
r n! ⋅ ∏ f i ( xi ; Θ) ⋅ R n − r (T ) (n − r )! i =1
(3.20a)
Pro data censorovaná poruchou je věrohodnostní funkce ve tvaru (3.20b):
L( x1 , x2 ,..., x n ; Θ) =
8
Náhodný výběr
f (t ) = λ ⋅ e
− λt
r n! ⋅ ∏ f i ( xi ; Θ) ⋅ R n − r ( xr ) (n − r )! i =1
X( X 1 , X 2 ,..., X n ) je z exponenciálního rozdělení. Hustota pravděpodobnosti je
. Potom věrohodnostní funkce je
funkce se obdrží
(3.20b)
L(t ; Θ) = λ ⋅ e n
−λ
n
∑ ti i =1
. Logaritmováním věrohodnostní
n
ln L(t ; Θ) = n ⋅ ln λ − λ ∑ ti . Lokální extrém se určí derivováním podle neznámého i =1
d ln L(t ; λ ) n d 2 ln L(t ; λ ) n parametru λ . = − ∑ ti = 0 . Druhá derivace je rovna = − 2 < 0. 2 dλ dλ λ i =1 λ n Z rovnice pro první derivaci je bodový odhad parametru λˆ = n . ∑ ti n
i =1
Strana: 33 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Metoda momentů Dalším způsobem, jak lze získat neznámé parametry určitého rozdělení, je metoda momentů. Tato metoda spočívá v tom, že se porovnávají výběrové momenty získaných dat s odpovídajícími výběrovými momenty teoretického rozdělení. Pro získání Θ neznámých parametrů je nutné řešit Θ rovnic. Nutný je předpoklad, že až do Θ -ho momentu jsou řešení rovnice konečné a dále, že soustava rovnic je regulární. Metoda pravděpodobnostního papíru Další metodou - grafickou - jak získat neznámé parametry předpokládaného rozdělení, je metoda pravděpodobnostního papíru. Jestliže vynesené okamžiky poruch na vodorovnou osu pravděpodobnostního papíru a příslušné pravděpodobnosti poruchy vynesené na svislou osu tvoří přímku, potom lze učinit závěr, že data jsou z příslušného rozdělení. Nejčastější typy pravděpodobnostních papírů jsou: • • • •
exponenciální, Weibullův, normální, log-normální.
Pravděpodobnostní papíry jsou uvedeny např. v [40], [41], kde lze i nalézt návod jak vyčíst parametry příslušných rozdělení. 3.4.4
Exponenciální rozdělení Pomocí exponenciálního rozdělení se popisuje doba do poruchy komponenty v druhém období - normálního života. Základní charakteristiky exponenciálního rozdělení jsou uvedeny ve vzorcích (3.21)
R(t ) = e − λt
F (t ) = 1 − e − λt
h(t ) = λ
(3.21)
Parametr λ se získá pomocí metody maximální věrohodnosti pomocí vzorců (3.22)
λ=
r , T*
(3.22a),
kde T * je kumulovaná doba provozu. Kumulovaná doba se stanoví: • • •
všechny doby do poruchy komponent jsou známé T * =
n
∑t i =1 r
censorování I. typu - časem
i
(3.22b),
T * = ∑ ti + (n − r ) ⋅ T (3.22c), i =1 r
censorování II. typu - poruchami
T * = ∑ ti + (n − r ) ⋅ tr (3.22d), i =1
kde
ti - časový okamžik, při kterém došlo k poruše komponenty, T - celková doba zkoušky, tr - čas, při kterém nastala r -tá porucha a zkouška byla zároveň ukončena.
Konfidenční meze bodového odhadu parametru λ lze určit pomocí [42], [77].
Strana: 34 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
V praxi je nutné předpoklad použitelnosti popisu pravděpodobnosti bezporuchového provozu pomocí exponenciálního rozdělení ověřit. Statistické testy jsou podrobně uvedeny například v [42] a [76]. Vzorec (3.22a) lze upravit pro opravované komponenty za předpokladu, že střední doba provozu mezi poruchami je konstantní - kapitola 3.4.6, na:
z (t ) =
r T*
(3.22e)
3.4.5
Weibullovo rozdělení pro neopravované objekty Pro popis doby do poruchy komponenty v období časných poruch a degradace je vhodné použít Weibullovo rozdělení. Dále se Weibullovým rozdělením obvykle popisují komponenty s dobou do poruchy, které jsou namáhány. Weibullovo rozdělení se především používá v případech, kdy bylo statisticky prokázáno, že pro popis doby do poruchy nelze využít exponenciálního rozdělení. Základní charakteristiky Weibullova rozdělení jsou uvedeny ve vzorcích (3.23).
R (t ) = e
t −( ) β
F (t ) = 1 − e
α
t −( ) β
α
h(t ) =
β ⋅ t β −1 αβ
α,β > 0
(3.23)
Ze vzorců (3.23) je zřejmé: • • •
parametr β = 1 , potom Weibullovo rozdělení přechází na exponenciální rozdělení, parametr β > 1 , potom komponenta degraduje a má zvýšené množství poruch na konci fyzického života, parametr 0 < β < 1 , potom má komponenta zvýšené množství poruch v počátečních fázích provozu.
Před zjišťováním bodových odhadů parametrů je nutné zjistit, zda data o poruchách pocházejí z Weibullova rozdělení. Vhodné testy jsou uvedeny například v [45], [80]. Parametry Weibullova rozdělení lze získat pomocí metod, které jsou popsány v odstavci 3.4.2. Pomocí metody momentů (pouze pro necensorovaná data) se parametry rozdělení získají řešením rovnic (3.24).
E (T ) = αΓ(1 +
1
β
)
⎡ 2 1 ⎤ D(T ) = α 2 ⎢Γ(1 + ) − Γ 2 (1 + )⎥ , β β ⎦ ⎣
(3.24)
kde E (T ) je střední hodnota, D(T ) je rozptyl a Γ(x) je gamma funkce. Zjištěná střední hodnota a rozptyl z dat o poruchách komponent se zadají do uvedených rovnic a vypočtou se bodové odhady neznámých parametrů αˆ a βˆ Weibullova rozdělení. Metodou největší věrohodnosti se získají bodové odhady parametrů Weibullova rozdělení pomocí rovnic (3.25a, 3.25b, 3.25c). Literární odkazy jsou uvedeny v [1], [43] a [80].
Strana: 35 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
•
všechny doby do poruchy komponent jsou známé
⎛ n β ⎜ ∑ ti ⎜ i =1 ⎜ n ⎜ ⎝ •
1
⎞β ⎟ ⎟ =α ⎟ ⎟ ⎠
∑ (t β ln(t ) i
i
i =1
n
∑t
β
−
1
β
−
1 n ∑ lnti = 0 n i =1
(3.25a),
i
i =1
censorování I. typu - časem
⎛ r β ⎜ ∑ ti + (n − r )T β ⎜ i =1 ⎜ r ⎜ ⎝ •
n
1
⎞β ⎟ ⎟ =α ⎟ ⎟ ⎠
r
∑ (t β ln(t ) + (n − r )T β ln(T ) i
i
i =1
r
∑ t β + (n − r )T β
−
1
β
−
1 r ∑ lnti = 0 (3.25b), r i =1
−
1 r ∑ lnti = 0 (3.25c), r i =1
i
i =1
censorování II. typu - poruchami
⎛ r β β ⎜ ∑ ti + (n − r )tr ⎜ i =1 ⎜ r ⎜ ⎝
1
⎞β ⎟ ⎟ =α ⎟ ⎟ ⎠
r
∑ (t β ln(t ) + (n − r )t i
i
i =1
β r
r
∑ t β + (n − r )t β i
i =1
ln(tr )
−
1
β
r
kde ti jsou okamžiky, kdy došlo k poruše; n počet komponent; T celková doba zkoušky;
α , β bodové odhady Weibullova rozdělení.
Stanovení, zda vstupní data pro neopravované objekty splňují exponenciální nebo Weibullovo rozdělení, a následné zjištění bodových odhadů parametrů bylo autorem použito například v [90], [101] a [104]. 3.4.6
Weibullovo rozdělení pro opravované objekty Při testování konstantnosti proudu poruch lze využít testy, které jsou uvedeny v [76] 9. V těchto testech se zjišťuje, zda časovým prodlužováním zkoušky, doba mezi poruchami roste, je konstantní, nebo klesá Softwarové prostředky většiny firem, jako je např. Itemsoft či Reliasoft, umožňují zjistit parametry Weibullova rozdělení pouze pro neopravované objekty. Pro testování vhodnosti rozdělení je nutný předpoklad, že je nutné mít alespoň data o 6 poruchách. Nechť X je testovací veličina, dále doba do i-té poruchy je ti a kumulovaná doba zkoušky je T * . Potom, jestliže T * > tr , se testovací veličina zjistí: r
X=
∑t i =1
T
9
i
*
−r
T* 2
r 12
(3.26a).
Vzorce jsou značně obsáhlé, proto je uveden pouze odkaz.
Strana: 36 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Jestliže T * = tr , potom r −1
X=
∑t i =1
i
− (r − 1)
tr 2
(3.26b).
r −1 12
tr
Testovací veličina X má normované normální rozdělení. Zjištění, zda vstupní data pro opravované objekty splňují exponenciální nebo Weibullovo rozdělení, je použito například v práci [105]. 3.4.7
Nevýhody Weibullova rozdělení Nevýhoda Weibullova rozdělení s parametrem β > 1 spočívá v tom, že intenzita poruch v čase t = 0 je nulová. Degradace mnoha systémů se projevuje až při delším čase provozování. Z tohoto důvodu je vhodné popsat intenzitu poruch komponent pomocí superpozice Weibullova a exponenciálního rozdělení ve tvaru:
h(t ) = λ +
β ⋅ t β −1 . αβ
(3.27)
Tím je zajištěno, že v čase t = 0 je intenzita poruch rovna h(t ) = λ . Nevýhoda takto zvolené intenzity poruch spočívá v obtížném zjištění parametrů λ , α a β , nicméně intenzita poruch v tomto tvaru lépe charakterizuje chování vanové křivky. [46], [47] Druhý problém Weibullova rozdělení vyplývá z hodnot jeho parametrů, které jsou získávány ze všech dat o poruchách - tedy i z těch, kde se degradace zcela neprojevovala. Tím se získají hodnoty parametrů, které nejlépe odpovídají datům o poruchách od začátku sledování až do okamžiku výpočtu parametrů. Při opětovném sledování dat o poruchách a jejich následné analýze je stále stejně dlouhý čas, při kterém k degradaci nedocházelo. Ale zároveň je delší čas, při kterém k degradaci již docházelo. To má za následek, že se zkrátí střední doba do poruchy komponent a zvyšuje se hodnota parametru β . [48] Ve [49] a [50] je popsáno modifikované Weibullovo rozdělení, pomocí kterého lze popsat všechny tři části vanové křivky. Jeho charakteristiky jsou uvedeny ve vzorci (3.28). V literatuře [97] bylo ověřeno, že toto rozdělení je použitelné pro popis doby do poruchy, jestliže existuje dostatečný počet vstupních dat. V práci [97] je analyzováno 25 000 dat o poruchách se zjištěním, že rozdělení velice dobře popisuje všechny tři části vanové křivky. Nevýhoda tohoto rozdělení je v tom, že jednotlivé bodové odhady parametrů nelze (na rozdíl od vzorce 3.27) jednoduše interpretovat.
R(t ) = e − at
b λt
e
F (t ) = 1 − e − at
b λt
e
h(t ) = a(b + λt ) ⋅ t b −1 ⋅ eλt
(3.28)
a > 0; b, λ ≥ 0
Strana: 37 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
3.4.8
Použití dalších neexponenciálních rozdělení Nejčastěji používaná neexponenciální rozdělení doby do poruchy jsou Weibullovo rozdělení, Gamma rozdělení a jeho speciální případ Erlangova rozdělení, normální rozdělení a log-normální rozdělení. Použití Weibullova rozdělení je popsáno v předcházejících kapitolách. Gamma rozdělení se využívá pro popis doby do poruchy systémů, které jsou tvořeny shodnými zálohovanými komponentami. Dále se předpokládá, že doba do poruchy každé komponenty je popsána exponenciálním rozdělením. Pravděpodobnost bezporuchového provozu je: τ
∞
R (t ) =
1 τ α −1e β dτ α ∫ β Γ(α ) t
(3.29).
Speciálním případem gamma rozdělení je Erlangovo rozdělení, které vznikne pro β celé. Erlangovo rozdělení popisuje dobu do β -té poruchy komponenty, která má exponenciální rozdělení doby do poruchy. Dále doba do obnovy je nulová nebo popsaná exponenciálním rozdělením. Pravděpodobnost bezporuchového provozu je: β −1
(λ t ) j j! j =0
R(t ) = e − λt ⋅ ∑
(3.30).
Normální rozdělení se často používá pro popis doby do obnovy komponenty a dále pro komponenty, které silně degradují. Pravděpodobnost bezporuchového provozu je: ∞
1 R(t ) = e σ 2π ∫t
− (τ − µ ) 2
2σ 2
dτ
(3.31).
Strana: 38 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
3.5 Softwarové prostředky spolehlivosti využitelné v markovské analýze V [14] je uveden příspěvek, který se zabývá počítačovou podporou ve spolehlivosti. V tomto článku je uveden soupis jednotlivých softwarových produktů různých výrobců. V příloze 1 je uvedena stručná charakteristika spolehlivostního softwaru popisující markovskou analýzu. Jedná se o produkty následujících firem: • • • • •
Item Software, [53] Relex Software, [54] IsographDirect, [55] SoHaR, [56] Logan, [57]
Ke každému programu, který je uveden v seznamu, byla prostudována trial verze.
Strana: 39 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
4. Markovská analýza a neexponenciální rozdělení doby do přechodu V současné době se modeluje markovskou analýzou ve spolehlivosti doba do poruchy/obnovy především pomocí exponenciálního rozdělení. Hlavním důvodem je to, že intenzita přechodu mezi stavy je konstantní. V této kapitole bude vysvětlen způsob jak modelovat pomocí markovských technik neexponenciální rozdělení doby do přechodu. Tato kapitola shrnuje poznatky, které jsou uvedeny v kapitole 3. Na jedné straně používá matematický aparát uvedený v kapitole 3.1. Na druhé straně, z pohledu spolehlivosti, využívá markovskou analýzu viz kapitola 3.3 a dále pro popis neexponenciálních rozdělení a stanovení jejich parametrů kapitolu 3.4. Při analýze spolehlivosti užitím markovských technik se využívá základního předpokladu, který je uveden ve vzorci (3.8). Tento předpoklad je ve většině odborných knih popsán tak, že navržený model musí být bez paměti. Předpoklad je splněn, jestliže doba do přechodu je popsána pomocí exponenciálního rozdělení, potom jsou intenzity přechodů hij (t ) konstantní a je možné využít pro řešení analýzy postupu uvedeného v kapitole 3.3. V disertační práci se autor zabývá případy, kdy toto není splněno. Jestliže nelze využít pro popis doby do přechodu exponenciálního rozdělení, hledají se složitější modely, které popisují dynamiku intenzity přechodu - např. Weibullovo rozdělení. Způsoby zjištění parametrů Weibullova rozdělení pro dvoustavový opravovaný a neopravovaný systém jsou uvedeny v kapitole 3.4.5 a 3.4.6. Často se využívá při popisu doby do obnovy komponent normální rozdělení. Způsob zjištění parametrů normálního rozdělení je uveden například v [6], [7], [8] a [9]. Předpoklad uvedený ve vzorci (3.8) a popis doby do přechodu pomocí neexponenciálního rozdělení jsou dva navzájem si odporující předpoklady. Pro analýzu spolehlivosti systému je nutné buď: • •
použít zjednodušené předpoklady, které zaručují markovskou vlastnost, nebo použít pro výpočet nehomogenní markovské procesy.
V následujících kapitolách bude uveden způsob, jak lze modelovat úlohy spolehlivosti pomocí markovské analýzy spolehlivosti, když intenzity přechodů nejsou konstantní. V kapitole 4.1 je uveden rozdíl v modelování a následném řešení opravovaných a neopravovaných objektů. V kapitole 4.2 je vysvětlen způsob modelování, který neporušuje markovskou vlastnost a je založen na diskrétních markovských procesech. Tato metoda je vhodná pro modelování opravovaných objektů. V kapitole 4.3 je uveden způsob modelování založený na metodě Monte Carlo. Tato metoda je vhodná i pro modelování neopravovaných objektů. Při modelování pomocí metody markovské analýzy se v kapitole 4 využívají markovské nehomogenní procesy, nebo stochastické procesy. Přestože není u těchto procesů zajištěna základní markovská vlastnost, že systém je bez paměti - vzorec (3.8), bude tato metoda spolehlivosti i nadále nazývána markovskou analýzou. Hlavní důvod je ten, že pojmenování této analýzy spolehlivosti je vžité. Další důvod je ten, že autor disertační práce se necítí povolán navrhnout nový termín.
Strana: 40 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
4.1 Modelování opravovaných a neopravovaných objektů Markovská analýza (jejíž princip je uveden v kapitole 3.3) předpokládá, že všechny bodové odhady intenzit přechodů mezi stavy jsou konstantní. Tento předpoklad, zajišťující, že systém je tzv. bez paměti, je uveden ve vzorci (3.8). U řady technických systémů však vyvstává otázka jak modelovat systém, který obsahuje více než 2 stavy a zároveň jednotlivé komponenty nesplňují podmínku, že intenzity přechodů mezi stavy jsou konstantní. Tato problematika byla pro dvoustavové systémy nastíněna v kapitole 3.4.5 a 3.4.6. V této kapitole je uveden prostředek, který umožní propojit na jedné straně markovskou analýzu a na druhé straně nekonstantní intenzity přechodů mezi stavy. V kapitole 3.4.4 je uveden způsob zjištění bodového odhadu parametru exponenciálního rozdělení. Tímto rozdělením lze popsat dobu do poruchy komponenty, která
r (3.22a) se T* r vypočte bodový odhad intenzity poruch pro neopravované komponenty. Dle vzorce z (t ) = * T
nedegraduje (komponenta je v období normálního života). Dle vzorce λ =
(3.22e) se získá bodový odhad proudu poruch pro opravované objekty. Způsob výpočtu bodového odhadu parametru exponenciálního rozdělení - λ není tedy závislý na tom, zda je objekt opravovaný nebo neopravovaný. V kapitolách 3.4.5 a 3.4.6 jsou uvedeny způsoby zjištění parametrů Weibullova rozdělení pro neopravovaný a opravovaný objekt. Pro každý typ objektů jsou způsoby výpočtu bodových odhadů parametrů α , β odlišné. V kapitole 4.1.1 a 4.1.2 jsou definovány a popsány dva základní termíny: • •
systém stejně spolehlivý jako starý (as good as old), systém stejně spolehlivý jako nový (as good as new).
V prvním případě, kdy objekt je po poruše opravován, lze předpokládat, že parametr proudu poruch po obnově komponenty je shodný jako před poruchou. Definuje se termín, že komponenta je „as good as old“. Ustálený překlad do češtiny pro tento výraz neexistuje, proto je využit v této práci termín „stejně spolehlivý jako starý“. Tento termín se týká především opravovaných objektů. Ve druhém případě, kdy objekt je po poruše nahrazen za nový, lze předpokládat, že intenzita poruch po obnově je shodná jako u nové komponenty. Definuje se termín, že komponenta je „as good as new“. Ustálený překlad do češtiny pro tento výraz neexistuje, proto je využit v této práci termín „stejně spolehlivý jako nový“. 4.1.1
Systém stejně spolehlivý jako starý U komponenty, která je stejně spolehlivá jako stará, je po obnově parametr proudu poruch z (t ) (intenzita přechodu hij (t ) ) shodný, jako kdyby k poruše nedošlo. Existuje tím vztažný bod - počátek t = 0 , ke kterému je vztažen veškerý čas intenzit přechodů hij (t ) mezi stavy. Tento model je vhodný pro systémy, kde je komponenta po poruše opravena. Úloha stanovení pohotovosti systému se řeší pomocí obyčejných diferenciálních rovnic s nekonstantními koeficienty. U úlohy je splněna základní markovská podmínka (3.8), že systém je bez paměti, proto lze úlohy řešit numericky, například metodou Runge-Kutta. Tato úloha je řešena v kapitole 4.2. Strana: 41 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Okamžité hodnoty proudu poruch
z (t ) (intenzity přechodu mezi stavy hij (t ) )
opravované komponenty jsou zobrazeny na obr. 2. Autor vyzkoušel v řadě prací [86], [89], [91], [92], [93], [96] a [100], že je pro řešení markovské analýzy ve spolehlivosti pro systém stejně spolehlivý jako starý dostačující Eulerova metoda. Její nevýhoda spočívá v tom, že z důvodu numerické stability musí být menší časový krok než u složitějších a z pohledu doby výpočtu výhodnějších metod RungeKutta. Tím je prodloužena doba výpočtu. Eulerova metoda má však tu výhodu, že lze snadněji a intuitivněji modelovat preventivní údržbu (viz kapitola 5), než je tomu při výpočtu metodou Runge-Kutta. 4.1.2
Systém stejně spolehlivý jako nový U objektu, který je stejně spolehlivý jako nový, je po poruše komponenta nahrazena novou. Intenzita poruch λ (t ) nové komponenty v okamžiku výměny γ je shodná jako intenzita poruch komponenty staré v čase t = 0 . Časový posun γ je vztažen k časovému okamžiku obnovy. Tento okamžik γ není přesně znám. Tím vznikají velké potíže při řešení úlohy, protože systém je s pamětí a nelze použít pro řešení popis obyčejnými diferenciálními rovnicemi 1. řádu. Úloha se řeší pomocí metody Monte Carlo. Tato úloha je řešena v kapitole 4.3. Na obr. 2 je uvedena závislost okamžité intenzity poruch λ (t ) pro objekt po poruše stejně spolehlivý jako nový a proudu poruch z (t ) pro objekt po poruše stejně spolehlivý jako starý. Způsoby řešení systémů, které jsou buď „stejně spolehlivý jako starý“, nebo „stejně spolehlivý jako nový“ se navzájem zcela vylučují. Z tohoto důvodu jsou v disertační práci v jednotlivých kapitolách uvedeny odděleně oba způsoby. Stejně tak i zdrojové texty v softwaru MATLAB®, které jsou uvedeny v příloze 2 a příloze 3, jsou odlišné pro obě metody.
Obr. 2: Závislost intenzity/proudu poruch na čase pro objekty stejně spolehlivé jako starý a stejně spolehlivé jako nový Existují i jiné modely, kde výsledná intenzita přechodu komponenty hij (t ) je omezená okamžitou intenzitou poruch λ (t ) a parametrem proudu poruch
z (t ) . Tento model Strana: 42 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
intenzit/proudu poruch není popsán v disertační práci. Řešení takového modelu by bylo možné metodou Monte Carlo. Nová intenzita přechodu komponenty po obnově by byla vypočtena pomocí vzorců, které by zadal autor analýzy. Při řešení takovéto úlohy je pro účinnou analýzu spolehlivosti nutné velké množství vstupních dat o poruchách.
Strana: 43 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
4.2 Objekt po obnově stejně spolehlivý jako starý Jak již bylo popsáno v úvodní kapitole 4.1 bude v tomto odstavci popsán model, kdy objekt je po obnově stejně spolehlivý jako starý. Úloha se řeší pomocí soustavy obyčejných diferenciálních rovnic - kapitola 4.2.1. Tento model je vhodný pro opravované systémy. V prostředí MATLAB® byl vytvořen software, který řeší tuto úlohu. Zdrojový text je uveden a popsán v příloze 2. 10 4.2.1
Popis systému pomocí markovských řetězců Ideou při řešení modelu je nahrazení parametru proudu poruch z (t )
11
nebo obecně
okamžitých intenzit přechodů mezi stavy hij (t ) střední intenzitou přechodů mezi stavy
hij (tm , tm +1 ) podle vzorce (4.1). Následně ze znalostí středních intenzit přechodů mezi stavy hij (tm , tm +1 ) lze zjistit matici pravděpodobností přechodu P(m) dle vzorce (4.2). Prvky této matice představují diskrétní markovský nehomogenní proces.
1 ∀i, j ∈ (1,..., N ) hij (tm , tm +1 ) = (tm +1 − tm )
t m +1
∫ h (t )dt
(4.1)
ij
tm
P(tm , tm +1 ) = I + h ⋅ (tm , tm +1 ) = ⎡1 ⎢0 =⎢ ⎢... ⎢ ⎣0
0 1 ... 0
... ... ... ...
⎡ N h12 (tm , tm +1 ) ⎢− ∑ h1i (tm , tm +1 ) 0⎤ ⎢ i =1,i ≠1 N ⎢ h (t , t ) 0 ⎥⎥ − ∑ h2i (tm , tm +1 ) + 21 1 m m + (tm +1 − tm ) ⎢ i =1, i ≠ 2 ...⎥ ⎢ ... ... ⎥ ⎢ 1⎦ ⎢ hN 1 (tm , tm +1 ) hN 2 (tm , tm +1 ) ⎢⎣
⎤ ⎥ ⎥ ... h2 N (tm , tm +1 ) ⎥ ⎥ ⎥ ... ... ⎥ N ... − ∑ hNN (tm , tm +1 )⎥ ⎥⎦ i =1, i ≠ N
...
h1N (tm , tm +1 )
(4.2) Z důvodu numerické stability musí být matice P(tm , tm +1 ) stochastická, tj. prvky matice
⎡ p11 (tm , tm +1 ) ⎢ p (t , t ) P (tm , tm +1 ) = ⎢ 21 m m +1 ⎢ ... ⎢ ⎣ p N 1 (tm , tm +1 )
p12 (tm , tm +1 )
...
p22 (tm , tm +1 ) ... ... ... pN 2 (tm , tm +1 ) ...
p1N (tm , tm +1 ) ⎤ p2 N (tm , tm +1 ) ⎥⎥ ⎥ ... ⎥ pNN (tm , tm +1 )⎦
hodnotou 0 ≤ pij (tm , tm +1 ) ≤ 1 . Dále musí být splněno
∀i, j ∈1,..., N
pro
N
∑ p (t j =1
ij
m
omezeny
, tm +1 ) = 1 . Doporučuje se
z důvodu přesnosti řešení, není však nutné, aby hodnoty prvků pij (tm , tm +1 ) byly buď 10
Softwarový prostředek zároveň řeší i přínosy uvedené v dalších kapitolách. Proto je uveden až na konci teoretické části. 11
Je možno uvést i intenzitu poruch, za přijmutí předpokladu, že neopravovaný objekt má po nahražení shodnou intenzitu poruch jako před poruchou.
Strana: 44 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
pij (t m , t m+1 ) → 0 + , nebo pij (tm , tm +1 ) → 1− . pij (tm , tm +1 ) představuje pravděpodobnost přechodu ze stavu i do stavu j v časovém intervalu (t m , t m +1 ) . Ze
znalosti
matice
pravděpodobnosti
přechodů
P(tm , tm +1 )
se
určí
hodnota
pravděpodobnostního vektoru p(tm +1 ) pomocí vzorce 4.3.
p(tm +1 ) = p(tm )P(tm , tm +1 )
(4.3)
Nechť je doba do přechodu popsána Weibullovým rozdělením s parametry α a β ,
β ⋅ t β −1 kde α , β > 0 . Potom je intenzita přechodu hij dána podle vzorce (3.21) hij (t ) = . Dle αβ vzorce (4.1) je střední intenzita přechodu na intervalu tm , tm +1 :
hij (tm , tm +1 ) =
tmβ +1 − tmβ α β (tm +1 − tm )
(4.4a)
V odstavci 3.4.7 je uvedeno rozdělení s intenzitou přechodu h(t ) = λ +
β ⋅ t β −1 . Toto αβ
rozdělení odstraňuje nevýhodu, kdy při β > 1 je h(0) = 0 . Střední intenzita poruch je na intervalu tm , tm +1 :
tmβ +1 − tmβ hij (tm , tm +1 ) = λ + β α (tm +1 − tm )
(4.4b)
Protože je zjišťování střední hodnoty na daném intervalu tm , tm +1
pro celou matici a
v každém časovém intervalu složité, je snaha výpočet zjednodušit. Pro výpočet je jednodušší a rychlejší použít hodnotu okamžité intenzity poruch ve středu intervalu tm , tm +1 = tm + 0.5 ⋅ l . Relativní chybu η , vzniknuvší změnou střední hodnoty intenzity přechodu v intervalu
tm , tm +1 na hodnotu intenzity přechodu v bodě h(
η=
t m +1 + t m ) , lze vyjádřit vzorcem (4.5): 2
t m+1 + tm ) 2 t +t h( m+1 m ) 2
h(t m , t m+1 ) − h(
(4.5a)
Pro intenzitu přechodu danou Weibullovým rozdělením s intenzitou přechodu uvedenou ve vzorci (4.4a) a (4.4b) je relativní chyba závislá na velikosti parametrů, které vstupují do vzorce 4.5 - α , β , λ , t a h . Dosazením rovnice (4.4b) do rovnice (4.5a) a následnou úpravou, se zjistí, že není závislá na intenzitě poruch λ . To znamená, že místo
t mβ +1 − t mβ , tj. vzorcem vzorce (4.4b) lze dosadit střední hodnotu ve tvaru hij (t m , t m+1 ) = β α (t m+1 − tm ) (4.4a). Pro relativní chybu η se obdrží následující vzorec:
Strana: 45 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
β ⋅ (tm + 0.5 ⋅ l ) β −1 tmβ +1 − tmβ − α β (tm +1 − tm ) αβ η= β ⋅ (tm + 0.5 ⋅ l ) β −1
(4.5b)
αβ
Jestliže je:
t m +1 + t m ) 2 β > 2 , potom je výraz (vzorec 4.5) kladný (konvexní funkce), t m +1 + t m h( ) 2 t +t h(t m , t m +1 ) − h( m +1 m ) 2 β < 2 , potom je výraz záporný (konkávní funkce). t m +1 + t m h( ) 2 h(t m , t m +1 ) − h(
•
•
Relativní chybu lze vyjádřit pomocí totálního diferenciálu za použití parciálních derivací. Ty jsou uvedeny v rovnicích (4.6). V případě, že β < 2 , je nutné vzorce (4.6) násobit koeficientem -1.
∂η =0 ∂α
(4.6a)
β
β
∂η tm +1 ⋅ ln(tm +1 ) − tm ⋅ ln tm − l ⋅ (tm + 0.5l ) β −1 − β ⋅ l ⋅ (tm + 0.5l ) β −1 ⋅ ln(tm + 0.5l ) = − β ∂β t m +1 − t β β
β
β
β
− tm − β ⋅ l ⋅ (tm + 0.5l ) β −1 ) ⋅ (tm +1 ⋅ ln(tm +1 ) − tm ⋅ ln t ) (t − m +1 β β (tm +1 − tm ) 2 ∂η β ⋅ tm +1 = ∂h β
β −1
− β ⋅ (tm + 0.5l ) β −1 − 0.5β ⋅ ( β − 1) ⋅ l ⋅ (tm + 0.5l ) β − 2 − β β t m +1 − t m
β
(tm +1 − tm − β ⋅ l ⋅ (tm + 0.5l ) β −1 ) ⋅ ( β ⋅ tm +1 − β β (tm +1 − tm ) 2 ∂η β ⋅ tm +1 = ∂t β
β −1
β
− β ⋅ tm
β −1
(4.6b)
(4.6c)
)
β −1
− β ⋅ ( β − 1) ⋅ l ⋅ (tm + 0.5l ) β − 2 − β β t m +1 − t m
− tm − β ⋅ l ⋅ (tm + 0.5l ) β −1 ) ⋅ ( β ⋅ tm +1 (t − m +1 β β (tm +1 − tm ) 2
β −1
− β ⋅ tm
β −1
(4.6d)
)
Strana: 46 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Z důvodu obtížnosti odhalení vlivu jednotlivých parciálních derivací bude uveden příklad. V něm se sleduje rozdíl ve výsledcích střední intenzity přechodů pro Weibullovo rozdělení, které se nejčastěji používá pro popis nekonstantních intenzit poruch hij (t ) . Cílem je
stanovit relativní chybu η .
Jako referenční model pro porovnání výsledných středních intenzit jsou použity následující parametry α = 10000 , β = 1,5 , tm = 100 h a tm +1 = 101 h. Je sledována chyba
výpočtu, jestliže se parametry jednotlivě mění na: α = 5000 , α = 20000 , β = 0,7 , β = 1,2 ,
β = 1,8 , tm = 10 h, tm = 1000 h a délka kroku tm +1 − tm = 5 h, tm +1 − tm = 20 h. Výsledky stability jsou uvedeny v tabulce 1. V prvních čtyřech sloupcích tabulky 1 jsou uvedeny parametry Weibullova rozdělení. V pátém sloupci je uvedena střední intenzita poruch vypočtená pomocí vzorce 4.1. V dalším sloupci je uvedena okamžitá intenzita poruch zjištěná ve středu intervalu. V posledním sedmém sloupci je zjištěna relativní chyba intenzit při výpočtu. Relativní chyba zjištěná z předem zadaných hodnot je 1,03.10-6. Tabulka 1: Chyba při výpočtu střední intenzity přechodu
α
β
tm
t m +1
h(tm , tm +1 )
[h]
[h]
[h-1]
h(
t m +1 + t m ) 2
t m +1 + t m ) 2 t +t h ( m +1 m ) 2
h(tm , tm +1 ) − h(
[h-1]
10000
1,50
100
101
1,504.10-5
1,504.10-5
1,03.10-6
5000
1,50
100
101
4,253.10-5
4,253.10-5
1,03.10-6
20000
1,50
100
101
5,317.10-6
5,317.10-6
1,03.10-6
10000
0,70
100
101
2,783.10-4
2,783.10-4
1,61.10-6
10000
1,20
100
101
4,782.10-5
4,782.10-5
6,60.10-7
10000
1,80
100
101
4,539.10-6
4,539.10-6
6,60.10-7
10000
1,50
10
11
4,860.10-6
4,861.10-6
9,45.10-5
10000
1,50
1000
1001
4,745.10-5
4,745.10-5
1,04.10-8
10000
1,50
100
105
1,519.10-5
1,519.10-5
2,48.10-5
10000
1,50
100
120
1,573.10-5
1,573.10-5
3,45.10-4
Z tabulky 1 lze zjistit, že rozdíl v intenzitách je prakticky na jeden krok zanedbatelný. Relativní chyba je významně závislá pouze na délce intervalu, což je způsobeno použitou numerickou metodou. Snížení lze dosáhnout buď zmenšením kroku, nebo použitím metody Runge-Kutta.
Strana: 47 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Výpočet v MATLABu® Pomocí softwaru MATLAB® je naprogramován prostředek, který umožňuje mimo jiné modelovat neexponenciální dobu do přechodu. Předpokladem je, že při analýze je pevně zakotven čas k počátku simulace. Při výpočtu je využito změny z matice intenzit přechodu h(t ) na matici pravděpodobnosti přechodu P(tm , tm +1 ) . Tato změna je popsána v kapitole 4.2.1. 4.2.2
Zdrojový text a popis je uveden pro systém, který je stejně spolehlivý jako starý v příloze 2. Ukazatele pohotovosti, které software stanovuje a týkají se popisu komponent s neexponenciální dobou do přechodu, jsou následující: •
funkce okamžité pohotovosti
• • • •
součinitel střední pohotovosti kumulovaná doba nepoužitelného stavu kumulovaná doba použitelného stavu funkce okamžité nepohotovosti
A(t ) , 12 A (0, t ) , TDT (t ) , TUT (t ) , U (t ) ,
•
součinitel střední nepohotovosti
U (0, t ) .
4.2.3
Řešený příklad Příklad, který je zde uveden, bude analyzován pomocí prostředku uvedeného v příloze 2. Pro porovnání výsledků bude shodný příklad vytvořen a analyzován v softwarovém prostředku IsographDirect, který umožňuje modelovat systémy s neexponenciální dobou do přechodu. Zadání příkladu: Nechť existuje systém, který je popsán pomocí sedmi stavů, které jsou označeny 1 až 7. Stav 1 je stav provozu, stavy 2 až 6 jsou částečně poruchové a stav 7 je poruchový. Nechť intenzita přechodu mezi dvěma po sobě označenými stavy vychází
1,5 ⋅ t 1,5 −1 . Po poruše je systém opravován, tedy je 10001,5 stejně spolehlivý jako starý a intenzita přechodu z každého stavu do stavu 1 je h(t ) = 0,01 h-1.
z Weibullova rozdělení a je h(t ) = 0,001 +
Přechodový diagram je zobrazen na obr. 3. Úkolem je zjistit ukazatele pohotovosti. Vyjádření výsledků analýzy ze softwaru Isograph (vlevo) a MATLAB® (vpravo) jsou následující: • • • •
A(10000) = 0,9977 , součinitel střední pohotovosti A (0,10000) = 0,9991 , kumulovaná doba nepoužitelného stavu TDT (10000) = 8,54 TDT (10000) = 8,53 h, kumulovaná doba použitelného stavu TUT (t ) = 9990 h TUT (t ) = 9991,5 h, funkce okamžité nepohotovosti U (10000) = 0,00227 U (10000) = 0,0023 ,
•
součinitel střední nepohotovosti
•
funkce okamžité pohotovosti
A(10000) = 0,998
U (0,10000) = 0,000853
Graf nepohotovosti systému získaného z Isographu a MATLABu® jsou uvedeny na obr. 4. Oba grafy i výsledné hodnoty se shodují.
12
Způsob výpočtu jednotlivých parametrů jsou uvedeny v kapitole použité značení.
Strana: 48 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Obr. 3: Přechodový diagram řešeného příkladu
Obr. 4: Porovnání nepohotovostí pomocí softwaru Isograph (vlevo) a pomocí MATLABu® (vpravo)
Strana: 49 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
4.3 Objekt po obnově stejně spolehlivý jako nový Komponenta, která je po poruše nahrazena (je stejně spolehlivá jako nová), má intenzitu poruch vztaženu k počátku provozu komponenty. Tím, že čas obnovy není znám, ztrácí se základní markovská vlastnost (vzorec 3.8), že systém je bez paměti. Analyzovaný systém nelze řešit pomocí diferenciálních rovnic s nekonstantními koeficienty, které jsou uvedeny ve vzorci 3.15. Úloha se řeší například pomocí metody Monte Carlo (kapitola 3.3.5). Její princip je založen na mnohonásobném opakování stejného náhodného pokusu a na sledování jeho náhodného procházení v čase. Na analyzovaný systém se vytvoří přechodový diagram. Tento přechodový diagram je svojí topologií shodný s přechodovým diagramem systému, který splňuje markovskou vlastnost (komponenta po obnově je stejně spolehlivá jako stará). Odlišná jsou však pravidla intenzit přechodů mezi jednotlivými stavy. Každému přechodu se definují podmínky pro změny do dalších stavů. Tyto podmínky mohou být vztaženy buď k počátečnímu času simulace (obdobně jako v kapitole 4.2), nebo k jinému časovému okamžiku (nejčastěji k obnově komponenty). Určí se, kdy dojde k jednotlivým přechodům v systému a jaká je nová konfigurace dosažená tímto přechodem. Každý přechod, který je popsán v přechodovém diagramu má ve zdrojovém textu programu vlastní proceduru. Jednotlivé typy vkládaných procedur jsou odlišné podle typu přechodu a budou uvedeny v příslušných kapitolách 4.3.1, 5.3 a 6.2. Jestliže se do přechodového diagramu přidá/odebere některý přechod ze stavu i do stavu j , vloží/odstraní se i v zdrojovém textu řešícího programu příslušná procedura. V každé realizaci metody Monte Carlo dochází k náhodnému procházení jednotlivými stavy přechodového diagramu systému. Nechť je systém v čase t ve stavu i . Pro všechny možné přechody vycházející ze stavu i se vypočte doba do přechodu podle pravidel daných přechodem. Výsledný stav j (z množiny stavů vycházejících ze stavu i ) je dán přechodem, reprezentujícím minimální okamžik přechodu, z množiny přechodů vycházejících ze stavu i . Fyzický čas, kdy dojde k přechodu ze stavu i , je dán součtem minimální doby do přechodu ze stavu i a časem, kdy došlo k přechodu do stavu i . 4.3.1
Generování doby přechodu mezi stavy Čas do přechodu ze stavu i do dalšího stavu se zjistí pomocí vygenerovaného náhodného čísla, které reprezentuje pravděpodobnost v distribuční funkci. Distribuční funkce je funkcí parametrů statistického rozdělení a času F (t ) = f ( parametry, t ) . Řeší se inverzní úloha, kdy jsou známy parametry rozdělení, a náhodným číslem je vygenerována F (t ) . Není znám čas do přechodu t . Tedy t = f ( F (t ), parametry ) . V tabulce 2 jsou uvedeny vzorce pro některá nejčastěji používaná statistická rozdělení pro generování časů do přechodu. Do zdrojového textu se přidá procedura pro generování časů do přechodu. Každá procedura představuje jeden přechod ze stavu i do stavu j . Procedura má následující základní strukturu t2=t+inten1; if (t2<=t1) t1=t2; state=j; end
Strana: 50 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Tato procedura volá funkci, která je v tomto případě nazvána inten1. Pokud je funkce nazvána jiným názvem, bude použito i jiného názvu v příkazu. Pro generování doby do přechodu se využije jedné z následujících možností: náhodný generátor z daného rozdělení implementovaný v softwaru MATLAB®, matematické vzorce, pokud není pro dané statistické rozdělení v MATLABu® generátor implementován a nelze použít matematické vzorce pro obtížné vyjádření t = f ( F (t ), parametry) , lze využít numerického výpočtu. Použití jednotlivých rozdělení je ve stručnosti popsáno v kapitole 3.4.8. • • •
V prvním případě, kdy se využívá pro zjištění potenciálního času přechodu ze stavu i do stavu j (pomocí funkce se určí pouze okamžik možného přechodu ze stavu i do stavu j ) náhodného generátoru z MATLABu®, se využije následující funkce: function t3=inten1 definování jména funkce inten1 a výstupní proměnné t3 t3=exprnd(10000); zjištění, za jakou dobu by došlo k přechodu
(generátor z exponenciálního rozdělení) end
V druhém případě, kdy se využívá pro zjištění potenciálního času přechodu ze stavu i do stavu j matematického vzorce se využije v MATLABu® naprogramovaná následující funkce: function t3=inten1 definování jména funkce inten1 a výstupní proměnné t3 t3=-log(rand())/0.001; zjištění, za jakou dobu by došlo k přechodu 13 end
Ve třetím případě, kdy se využívá pro zjištění potenciálního času přechodu ze stavu i do stavu j numerického výpočtu, se využije v MATLABu® naprogramovaná následující funkce: function t3=inten1 definování jména funkce inten1 a výstupní proměnné t3 nahcislo=rand(); definování náhodného čísla funkce=0; hodnota distribuční funkce t3=0; čas, při kterém by došlo k přechodu cyklus while - zjištění prvního odhadu možného přechodu cyklus je ukončen při prvním nesplnění podmínky funkce=1-exp(-0.000001*t3); zadání distribuční funkce t3=t3+100; přidání času k okamžiku možného přechodu end konec cyklu while while (nahcislo>funkce)
13
t3=t3-200; funkce=0;
odpočet času po prvním nesplnění podmínky vynulování distribuční funkce pro další cyklus
while (nahcislo>funkce)
shodný cyklus while jako nahoře přesné zjištění okamžiku možného přechodu
V MATLABu® funkce
log( x) vypočítá přirozený logaritmus čísla x Strana: 51 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
funkce=1-exp(-0.000001*t3); t3=t3+1; end t3=t3-1;
odpočet času po prvním nesplnění podmínky
end
Při modelování pomocí výše uvedených příkladů je nutné zajistit, aby doba do možného přechodu t3 nebyla záporná. Tento předpoklad se týká především normálního rozdělení. Funkce s podmínkou, že doba do přechodu nebude záporná je uvedena níže. function t3=inten1 definování jména funkce inten1 a výstupní proměnné t3 t3=normrnd(10000,1000); zjištění za jakou dobu by došlo k přechodu if (t3<0) ošetření zápornosti doby do přechodu t3=0; end end
V tabulce 2 je uveden seznam nejčastějších rozdělení popisující dobu do poruchy a vzorce, pomocí kterých lze daný přechod mezi stavy řešit.
Strana: 52 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Tabulka 2: Vzorce pro generování doby do přechodu Funkce
F (t )
t = f (Pr, par ) v MATLABu®
h(t )
Rozdělení
2. způsob řešení
3. způsob řešení
h(t )
Numericky
1. způsob řešení
Numericky
-----
− ln(Pr) − λt exponenciální 1 − e
Weibullovo
1− e
λ β t β −1 αβ
t −( ) β
α
λ
exprnd(λ)
α β − ln(Pr)
wblrnd(α,β)
Numericky
normrnd(µ,σ)
Numericky
lognrnd(µ,σ)
− (t − µ ) 2 2σ 2
e t
Normální
1 e σ 2π ∫0
− (τ − µ ) 2 2σ
dτ
2
∫e
Log-normální
1 σ 2π
∫
e
2σ 2
1 e t dτ
τ
0
∞
∫
− (ln t − µ ) 2 2σ 2
− (ln τ − µ ) 2
e
2σ 2
dτ
τ
t
λ β −1
Erlangovo
dτ
2σ 2
t
− (ln τ − µ ) 2 t
− (τ − µ ) 2
∞
(λ t ) j! j =0
1 − e − λt ⋅ ∑
j
β −1
1 j j = 0 ( β − 1 − j )!(λt )
( β − 1)!∑
β
∑ j =1
− ln(Pr j )
λ
Gamrnd(α,β)
β ∈Z
t
t α −1e β τ
t
Gamma
1 τ α −1e β dτ α β Γ(α ) ∫0 1− e 1− e
t ⎡t ⎤ − ⎢ f (τ ) dτ + g (τ ) dτ ⎥ ⎢⎣ 0 ⎥⎦ 0
∫
∫
t − λt − ( ) β
α
∞
∫τ
τ α −1 β
e dτ Numericky
t
f (t ) + g (t )
λ+
Gamrnd(α,β)
Numericky
-----
Numericky
-----
β −1
βt αβ
Výpočet v MATLABu® Pomocí softwaru MATLAB® je naprogramován prostředek, který umožňuje modelovat systém, kde objekty jsou po poruše stejně spolehlivé jako nové. Zdrojový text a popis je uveden v příloze 3. Ukazatele spolehlivosti, které software zjišťuje a týkají se popisu komponent s neexponenciální dobou do přechodu, jsou následující: 4.3.2
•
funkce okamžité pohotovosti
•
součinitel střední pohotovosti
A(t ) , A (0, t ) , Strana: 53 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
• • •
kumulovaná doba nepoužitelného stavu kumulovaná doba použitelného stavu funkce okamžité nepohotovosti
TDT (t ) , TUT (t ) , U (t ) ,
•
součinitel střední nepohotovosti
U (0, t )
Strana: 54 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
4.4 Příklad Příklad, který je zde uveden, bude analyzován pro objekty stejně spolehlivé jako nové (určen pro neopravované objekty). Tento prostředek je uveden v kapitole 4.3.2, software popsán zdrojový text popsán v příloze 3. Pro porovnání výsledků bude shodný příklad analyzován pomocí softwaru, který analyzuje objekty stejně spolehlivé jako staré (opravované objekty) a byl uveden v kapitole 4.2.2. Mějme systém, který je tvořen čtyřmi identickými komponentami. Komponenty jsou v pohotovostní záloze. Intenzita poruch každé komponenty je ve tvaru h(t ) = λ +
β t β −1 . αβ
Parametry rozdělení jsou: λ = 0,0001 , α = 10000 , β = 1,5 . Komponenty mají intenzitu oprav µ = 0,01 . Komponenty jsou provozovány shodnou dobu. Cílem příkladu je porovnat nepohotovost systému U (t ) , jestliže komponenty jsou po poruše buď stejně spolehlivé jako staré nebo stejně spolehlivé jako nové. Přechodový diagram je zobrazen na obr. 5.
1
µ
λ
2
µ
λ
3
µ
4
λ
Obr. 5: Přechodový diagram řešeného systému Doba simulace byla zvolena 1000 a 10 000 h. Výsledkem jsou grafy střední nepohotovosti systému U (t ) . Při analýze metodou Monte Carlo, řešení pro objekty stejně spolehlivé jako nové, bylo použito 106 simulačních realizací. Při analýze metodou RungeKutta, pro objekty stejně spolehlivé jako staré, byl použit časový krok 0,1 h. Výsledky jsou uvedeny na: • • • •
obr. 6 - doba simulace 1000 h, systém stejně spolehlivý jako nový, obr. 6 - doba simulace 1000 h, systém stejně spolehlivý jako starý, obr. 7 - doba simulace 10 000 h, systém stejně spolehlivý jako nový, obr. 7 - doba simulace 10 000 h, systém stejně spolehlivý jako starý. Z grafů jsou zřejmé následující výsledky pro oba typy systémů:
•
systém stejně spolehlivý jako starý, doba simulace 1000 h U (1000) = 1,55 ⋅ 10−6 ,
•
systém stejně spolehlivý jako nový, doba simulace 1000 h U (1000) = 0,85 ⋅ 10−6
•
systém stejně spolehlivý jako starý, doba simulace 10 000 h
U (10000) = 8,1 ⋅ 10−6
•
systém stejně spolehlivý jako nový, doba simulace 10 000 h
U (10000) = 2,1 ⋅ 10−6
Strana: 55 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Obr. 6: Graf nepohotovosti systému, doba simulace 1000 h, systém stejně spolehlivý jako nový (vlevo), systém stejně spolehlivý jako starý (vpravo)
Obr. 7: Graf nepohotovosti systému, doba simulace 10 000 h, systém stejně spolehlivý jako nový (vlevo), systém stejně spolehlivý jako starý (vpravo)
Strana: 56 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
4.5 Zaměnitelnost modelů - systém je stejně spolehlivý jako nový/starý V kapitole 4.2 byl popsán model, který popisuje komponenty, které jsou po obnově stejně spolehlivé jako staré. Tento model je vhodný pro opravované objekty. V kapitole 4.3 byl popsán model popisující komponenty, kterou jsou stejně spolehlivé jako nové. Tento model je vhodný pro neopravované objekty. Oba tyto modely, které jsou svoji povahou zcela odlišné, lze za splnění určitých podmínek zaměňovat, protože mohou mít téměř shodné výsledky. Zaměnitelnost obou metod je možná, jestliže: • •
střední doba do první poruchy systému je výrazně delší než doba analýzy systému (vysoce spolehlivé systémy), intenzity přechodů mezi stavy hij (t ) jsou konstantní funkce.
Jestliže je intenzita přechodu klesající funkce (například Weibullovo rozdělení s parametrem β < 1 ), potom hodnota funkce okamžité nepohotovosti U (t ) je větší pro systém stejně spolehlivý jako nový než pro systém stejně spolehlivý jako starý. Jestliže je intenzita přechodu rostoucí funkce (například Weibullovo rozdělení s parametrem β > 1 ), potom hodnota funkce okamžité nepohotovosti U (t ) je větší pro systém stejně spolehlivý jako starý než pro systém stejně spolehlivý jako nový.
Strana: 57 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
5. Modelování údržbových zásahů V této kapitole bude do markovské analýzy ve spolehlivosti vložen prostředek, který umožňuje modelovat systémy s preventivní údržbou. Markovská analýza popsaná v kapitole 3.3 14 předpokládá, že jednotlivé přechody se odehrávají ve spojitém čase - porucha nebo obnova komponenty. Pro řešení se sestavila soustava diferenciálních rovnic, která se následně řešila. Výsledkem analýzy byla spojitá funkce okamžité pohotovosti A(t ) . Při preventivní údržbě dochází k přechodům mezi stavy podle předem daných pravidel daných koncepcí údržby. Dochází k nespojitostem funkce okamžité pohotovosti A(t ) . Terminologie a rozdělení typů údržby jsou normativně popsány v [67] a [68]. Základní rozdělení údržby je uvedeno na obr. 8 , který je převzat z [68] - str. 28.
Porucha
Doba mezi poruchami
Porucha
Údržba Před zjištěním poruchového stavu
Po zjištění poruchového stavu
Preventivní údržba
Údržba podle stavu
Plánovaná, nepřetržitá nebo na požádání
Údržba po poruše
Údržba s předem stanovenými intervaly
Plánovaná
Odložená
Okamžitá
Obr. 8: Rozdělení typů údržby Údržba (cit. [67] def. 7-1) - kombinace všech technických a administrativních činností, včetně činností dozoru, zaměřených na udržení ve stavu nebo navrácení objektu do stavu, v němž může plnit požadovanou funkci. Preventivní údržba (cit. [68] def. 7-1) - údržba prováděná v předem stanovených intervalech nebo podle předepsaných kritérií a zaměřená na snížení pravděpodobnosti poruchy nebo degradace fungování objektu. Údržba po poruše (cit. [68] def. 7-6) - údržba prováděná po zjištění poruchového stavu a zaměřená na uvedení objektu do stavu, ve kterém může vykonávat požadovanou funkci. Údržba podle stavu (cit. [68] def. 7-4) - preventivní údržba, která se skládá z monitorování výkonnosti a/nebo parametrů a z následných opatření. Údržba s předem stanovenými intervaly (cit. [68] def. 7-3) - preventivní údržba prováděná 14
Popis analýzy uvedený v kapitole 3.3 byl rozšířen v kapitole 4
Strana: 58 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
v souladu se stanovenými časovými intervaly nebo stanoveným počtem jednotek používání, avšak bez předchozího zjišťování stavu objektu.
Strana: 59 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
5.1 Modelování údržby pomocí markovské analýzy Složité systémy mají obvykle i komplexní politiku údržby, která se skládá z mnoha dílčích úkolů. Základní typy údržby v daných časech jsou dány plány technické údržby. Jednotlivé kroky technické údržby lze charakterizovat podle různých cílů a přístupů k zajištění udržovatelnosti objektu. V současné době jsou využívány především následující typy [12]: 1. 2. 3. 4. 5. 6. 7. 8. 9.
údržba po poruše, údržba po prohlídce, systém standardní preventivní údržby, preventivní periodické opravy, údržba postupnou výměnou skupin, údržba diferencovanou proporcionální péčí, údržba podle skutečného stavu, údržba za použití optimalizace celkových nákladů, údržba podle posouzení důsledků poruch.
V následujícím textu bude k jednotlivým typům údržeb (viz seznam výše) popsán způsob modelování pomocí markovské analýzy. Zvláštní důraz je kladen na to, zda použít model pro systém stejně spolehlivý jako starý (kapitola 5.2), nebo pro systém stejně spolehlivý jako nový (kapitola5.3). 1. Údržba po poruše Při této údržbě se připouští riziko poruchy tím, že objekt se obnoví až po jejím vzniku. Nepředpokládá se preventivní údržba. Modelování údržby po poruše je popsáno pomocí markovských procesů v kapitole 3.3. Z údajů o poruchovosti a z přechodového diagramu se zjistí intenzity přechodů reprezentující poruchy a obnovy zařízení. Vytvoří se matice intenzit přechodů h(t ) , určí se počáteční podmínky p(0) a řeší se soustava diferenciálních rovnic (3.15a). Základní markovská vlastnost předpokládá, že matice intenzit přechodů h(t ) není závislá na čase. To znamená, že doba do poruchy/obnovy splňuje exponenciální rozdělení. V kapitole 4 je uveden prostředek, který umožňuje modelovat i neexponenciální rozdělení doby do přechodu, tzn. intenzita přechodu je závislá na čase. Následně se řeší soustava diferenciálních rovnic (3.15b). Pro řešení diferenciálních rovnic (3.15) je typické, že pravděpodobnosti všech stavů jsou při údržbě po poruše spojité funkce. Změny pravděpodobností jednotlivých stavů mohou nastat v každém časovém okamžiku. Ostatní typy údržbových zásahů, které jsou uvedeny v seznamu, nelze v současné době pomocí markovské analýzy dostatečně podrobně řešit. V přechodovém diagramu dojde vlivem preventivní údržby v předem daném čase k přechodu ze stavu i do stavu j s danou pravděpodobností podle předem daných pravidel. K přechodu mezi stavy dojde v diskrétním čase a není tím splněn předpoklad, že změna pravděpodobností stavů je spojitá. Z tohoto důvodu nelze použít pro stanovení funkce okamžité nepohotovosti systému U (t ) rovnici (3.15). Proto jsou v disertační práci vytvořeny prostředky - viz kapitola 5.2 a 5.3, které umožňují modelovat preventivní údržbu. Oba prostředky se liší v tom, zda je systém po obnově stejně spolehlivý jako starý - kapitola 5.2, nebo stejně spolehlivý jako Strana: 60 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
nový - kapitola 5.3. V současné době je v teorii spolehlivosti modelována preventivní údržba především pro dvoustavové systémy. V disertační práci bude tento typ údržby pro vícestavové systémy vysvětlen pomocí metody markovské analýzy. V příloze 2 a 3 jsou uvedeny zdrojové texty naprogramované v MATLABu®, které umožňují analyzovat systémy a který je stejně spolehlivý jako starý a stejně spolehlivý jako nový. Tyto programy umožňují modelovat také systémy s preventivní údržbou a dále systémy, které mají nekonstantní intenzity přechodů mezi stavy. 2. Údržba po prohlídce - jsou předepsány intervaly prohlídek objektu, při kterých se zjistí skutečný technický stav a podle výsledku se provede oprava systému. Mohou nastat oba dva modelové případy, tzn. kdy objekt je po obnově stejně spolehlivý jako starý i jako nový. 3. Systém standardní preventivní údržby - je založen na znalosti normované životnosti jednotlivých dílů, které jsou opravovány v předepsaném rozsahu a lhůtách zásadně preventivním způsobem bez ohledu na jejich aktuální skutečný technický stav. K modelování standardní preventivní údržby lze použít model, kdy komponenta je po obnově stejně spolehlivá jako nová, protože komponenty jsou nahrazovány v předem daných časech bez ohledu na jejich aktuální skutečný technický stav. 4. Preventivní periodická oprava - je založena na organizačně technických opatřeních obsluhy a dozoru nad objektem, který poskytuje prostředky k preventivní a periodické údržbě podle daného plánu. Při modelování tohoto typu údržby lze použít podle druhu úlohy oba základní modely, tj. kdy systém je stejně spolehlivý jako nový i kdy systém je stejně spolehlivý jako starý. 5. Postupná výměna skupin Údržba, při které jsou postupně vyměňovány skupiny komponent, je používána u objektů, jejichž charakter, konstrukční a technologické provedení dovoluje snadnou montáž a demontáž celých skupin a podskupin. Tento typ údržby je možné řešit především pomocí modelu, kdy je systém stejně spolehlivý jako nový. 6. Diferencovaná proporcionální péče - rozlišuje péči o objekt podle důležitosti jeho funkce jako základního výrobního prostředku. Připouští určitý počet poruch a oprav s tím, že jsou posouzeny důsledky těchto poruch. Tomu je podřízena i diferencovaná péče o objekt. Tento údržbový zásah lze modelovat pomocí obou typů prostředků, tj. kdy systém je stejně spolehlivý jako starý i stejně spolehlivý jako nový. Při modelování pomocí markovské analýzy je nutné brát v potaz, že stavy výsledného přechodového diagramu popisují nejenom kombinaci poruch jednotlivých komponent, ale zároveň i jejich důležitost k udržení provozuschopnosti, bezpečnosti a ekonomičnosti. Z tohoto důvodu je vytvoření přechodového diagramu u tohoto typu údržby složitější a je nutné sledovat více vazeb než u jiných typů údržeb. 7. Údržba podle skutečného technického stavu - údržbové úkony se provádí na základě objektivně zjištěného skutečného stavu v okamžiku, kdy hodnota kontrolované veličiny charakterizující technický stav dosáhne mezní hodnoty. Tento údržbový zásah lze modelovat pomocí obou typů prostředků, tj. kdy systém je stejně spolehlivý jako starý i stejně spolehlivý jako nový. Strana: 61 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Při modelování přechodového diagramu je nutné rozlišit dva typy poruch - bezpečné a nebezpečné. Mohou nastat následující čtyři případy: • • • •
hodnota kontrolované veličiny charakterizující technický stav je vyšší než mezní hodnoty a nedojde následně k poruše - bezpečná porucha (kontrolované veličiny), hodnota kontrolované veličiny charakterizující technický stav je vyšší než mezní hodnoty a dojde následně k poruše - nebezpečná porucha, hodnota kontrolované veličiny charakterizující technický stav je nižší než mezní hodnoty a nedojde k poruše - provoz, hodnota kontrolované veličiny charakterizující technický stav je nižší než mezní hodnoty a dojde k poruše - nebezpečná porucha.
8. Optimalizací celkových nákladů na pořízení a údržbu objektu - údržba je založena na sledování jednotkových nákladů na pořízení a udržení objektu v provozu. Náklady jsou objektivním ukazatelem technického stavu, protože integrují vlastnosti spolehlivosti jako jsou bezporuchovost, pohotovost, opravitelnost, diagnostikovatelnost a složku nákladů na údržbu a opravu. Při tvorbě přechodového diagramu se bere ohled na ocenění ziskem nebo ztrátou jednotlivých stavů. Dále se oceňují přechody mezi stavy, které jsou reprezentovány poruchou, obnovou nebo údržbou systému. Tento typ údržby je již blízký spolehlivostně-ekonomickým metodám jako jsou například metoda nákladů životního cyklu (LCC) či údržba zaměřená na bezporuchovost (RCM) Při modelování tohoto typu se využívá obou typů modelů, tj. kdy systém je stejně spolehlivý jako starý, nebo stejně spolehlivý jako nový. 9. Podle posouzení důsledků poruch - tento typ údržby se zaměřuje na účelnost a hospodárnost se zdůrazněním požadavků na spolehlivost a bezpečnou funkci systému. Opět i tento typ údržbového zásahu lze modelovat pomocí obou typů modelů, tj. kdy systém je stejně spolehlivý jako nový, nebo kdy je systém stejně spolehlivý jako starý. Protože v seznamu u řady údržbových zásahů je uvedeno, že je možné úlohu řešit oběma základními typy modelů, je nutné popsat pravidla, kdy je vhodnější použít model systému stejně spolehlivého jako starý, nebo stejně spolehlivého jako nový. V prvním případě, kdy je použit model systému, který je po obnově stejně spolehlivý jako starý, jsou porouchané komponenty opraveny. Tento model lze využít, jestliže počet opravovaných komponent je výrazně větší než počet vyměňovaných komponent. Zdrojový text tohoto modelu je uveden a popsán v příloze 2. V druhém případě, kdy je použit model systému, který je po obnově stejně spolehlivý jako nový, jsou porouchané komponenty vyměněny za nové. Tento model lze využít, jestliže počet vyměňovaných komponent je výrazně větší než počet opravovaných komponent. Zdrojový text tohoto modelu je uveden a popsán v příloze 3. Disertační práce nebere v úvahu případ, kdy systém je mezi těmito dvěma případy. To znamená v případech kdy je významný počet komponent, které jsou opravovány i nahrazovány za nové. Potom je výsledná intenzita přechodu po obnově komponenty omezena hodnotami intenzit přechodu systému stejně spolehlivého jako nový a stejně spolehlivého jako starý.
Strana: 62 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
5.2 Preventivní údržba pro systém stejně spolehlivý jako starý 5.2.1
Modelování Systém, který je po poruše opraven (je stejně spolehlivý jako starý) má intenzitu poruch vztaženu k počátku analýzy. Tím je zachována markovská vlastnost, že systém je bez paměti. Analyzovaný systém se nahradí soustavou diferenciálních rovnic, které se řeší například metodou Runge-Kutta. Pro popis pravděpodobnosti stavů pi , při údržbě po poruše, je využito matice pravděpodobností přechodů P (t1 , t2 )
15
. Přechod mezi maticí intenzit poruch h(t ) a maticí
pravděpodobností přechodů P (t1 , t2 ) je dán vzorcem 3.8 - P (t1 , t2 ) = I + h(t )(t2 − t1 ) , t1 < t2 ,
t2 − t1 → 0 . Hodnoty pravděpodobnostního vektoru v čase se zjistí pomocí rovnic 3.7. Jestliže se v systému provádí preventivní údržba, nelze přechodový diagram systému popsat pomocí intenzity poruch/oprav. Je nutné připojit prostředek, který zajistí časově nespojitý popis pravděpodobnosti stavů - matici údržby U (t1 , t2 ) :
⎡ u11 (t1 , t2 ) u12 (t1 , t2 ) ⎢ u (t , t ) u (t , t ) 22 1 2 U (t1 , t2 ) = ⎢ 21 1 2 ⎢ ... ... ⎢ ⎣u N 1 (t1 , t2 ) u N 2 (t1 , t2 )
... u1N (t1 , t2 ) ⎤ ... u2 N (t1 , t2 ) ⎥⎥ . ⎥ ... ... ⎥ ... u NN (t1 , t2 )⎦
Prvek uij (t1 , t2 ) matice údržby U (t1 , t2 ) představuje pravděpodobnost preventivní údržby v intervale < t1 , t2 ) , kdy konfigurace systému odpovídající stavu i přejde údržbou do stavu j . V přechodovém diagramu se prvky uij (t ) zadávají společně s intenzitami přechodu
λij (t ) a µij (t ) . V zásadě lze vymezit dva základní případy preventivní údržby [60], [61]: • •
údržbovým zásahem jsou v daném časovém intervalu pokryty všechny komponenty a pravděpodobnost, že údržba je provedena správně, je 1, údržbovým zásahem nejsou v daném časovém intervalu pokryty všechny komponenty nebo pravděpodobnost, že údržba je provedena správně, je menší než 1.
V přechodovém diagramu reprezentuje každý stav konfiguraci provozních a poruchových komponent. Každá komponenta má určeny politikou údržby časy, kdy dojde k preventivní údržbě. Tím dojde v čase t údržbovým zásahem k přechodu ze stavu i do stavu j . Stav j je reprezentován stejnou konfigurací provozních a poruchových komponent až na ty, které byly právě obnoveny.
15
⎡ p11 (t1 , t2 ) L p1 j (t1 , t 2 )⎤ ⎢ ⎥ P(t1 , t2 ) = ⎢ M O M ⎥ , prvek pij (t1 , t2 ) představuje pravděpodobnost přechodu ze ⎢ pi1 (t1 , t 2 ) L pij (t1 , t 2 ) ⎥ ⎣ ⎦
stavu i do stavu j v časovém intervalu < t1 , t2 ) .
Strana: 63 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Mohou nastat následující varianty podob matice U (t1 , t2 ) : • •
•
•
V každém řádku matice U (t1 , t2 ) je prvek na diagonále roven 1. V systému není provedena v daném časovém intervalu žádná preventivní údržba. V některém řádku matice U (t1 , t2 ) je prvek na diagonále nulový. Údržbovým zásahem je v daném časovém intervalu pokryta příslušná konfigurace komponent a pravděpodobnost, že údržba je provedena správně, je 1. V některém řádku matice U (t1 , t2 ) je prvek na diagonále v intervalu 0 a 1. Údržbovým zásahem je v daném časovém intervalu pokryta příslušná konfigurace komponent a pravděpodobnost, že údržba je provedena správně, je menší než 1. V některém řádku matice U (t1 , t2 ) je prvek na diagonále roven jedné. Údržbovým zásahem není v daném intervalu pokryta příslušná konfigurace komponent. Pro matici U (t1 , t2 ) platí ∀i
∑u
ij
(t1 , t2 ) = 1 , ∀i, j uij (t1 , t2 ) ≥ 0 .
j
Maticově
lze popsat provedenou preventivní údržbu systému v čase t: p(t + u ) = p(t ) ⋅ U(t ) . 16 Matici přechodů mezi stavy kombinovanou s maticí údržby lze zapsat do jediné rovnice 5.1
p(t + ∆t ) = p(t )(P(t , t + ∆t ) ⋅ U(t , t + ∆t ))
∆t → 0
(5.1).
V rovnici 5.1 se předpokládá, že údržbový zákrok je proveden na konci časového intervalu (kroku). V případě, že údržbový zákrok je proveden na začátku časového intervalu o délce ∆t , přejde rovnice 5.1 na tvar 5.2:
p(t + ∆t ) = p(t )(U(t , t + ∆t ) ⋅ P(t , t + ∆t ))
∆t → 0
(5.2).
Je tedy nutný předpoklad, z důvodů přesnosti, že matice P (t1 , t2 ) má na diagonále prvky, které jsou blízké hodnotě 1. Tento předpoklad je většinou splněn, protože časový krok, kterým je násobena matice intenzit přechodů h , je z důvodu stability malý. 5.2.2
Určení hodnot prvků matice údržby Zdrojový text je naprogramován v MATLABu®. Při výpočtu se předpokládá, že časové okamžiky preventivní údržby jsou periodické. Zadává se perioda a pravděpodobnosti reprezentující míru pokrytí údržbovým zásahem. Při běhu programu se sleduje, zda v daném časovém kroku dojde k preventivní údržbě. Nechť systém je složen z komponent k1 , k2 ,...kn různého typu. Každý stav i ,
i = 1,..., N , přechodového diagramu je dán kombinací počtu komponent k1 , k2 ,...kn v provozu a komponent v poruše. Nechť je na těchto komponentách prováděna preventivní údržba s různou periodou. Potom lze popsat pro každou komponentu/skupinu komponent matici údržby U1 (t ), U 2 (t ),..., U n (t ) . Každá z matice U1 (t ), U 2 (t ),..., U n (t ) má shodné vlastnosti jako matice údržby U (t ) , tj. stochastické matice. Pravděpodobnostní vektor v čase t po provedené údržbě se vypočte:
16
p(t + u ) je pravděpodobnostní vektor po provedené údržbě v čase t . Strana: 64 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
p(t + u ) = p(t ) ⋅ U1 (t ) ⋅ U 2 (t ) ⋅ ... ⋅ U n (t )
(5.3).
Nechť stav i , i ∈ 1,..., N , přechodového diagramu je dán kombinací komponent
k1 , k 2 ,...,k l ,..., k n ve stavu provozu a v poruše. Nechť se sestavuje matice údržby l -té komponenty - tj. matice U l (t1 , t2 ) . Nechť stav i obsahuje r porouchaných komponent s označením l . Nechť jsou preventivní údržbou pokryty všechny porouchané komponenty s pravděpodobností p . Nechť systém přejde v čase t opravou jedné komponenty ze stavu i do stavu j a následně opravou dalších komponent do dalších stavů. Potom i -tý řádek matice U l (t ) je tvořen nenulovými prvky, které jsou na pozicích i , j ,…, za předpokladu, že 0 < p < 1 . Jestliže p = 1 , potom jsou nenulové prvky v i -tém řádku na pozicích j ,… . Hodnota těchto prvků je dána binomickým rozdělením s parametry Bi (r , p ) . Z důvodu stability úlohy nesmí být žádný z prvků matice P záporný. Jestliže je některý prvek matice P záporný, doporučuje se zmenšit časový krok analýzy. Při analýze je nutné sledovat, zda některý prvek vektoru p(t ) není větší než 1 nebo menší než 0. Dále se sleduje, zda 1 − ε <
∑ p (t ) < 1 + ε , kde ε i∈Ω
i
je předem definovaná chyba výpočtu.
Chybu výpočtu při opakovaném násobení matice P (t1 , t2 ) lze odhadnout provedením shodné analýzy při delším časovém kroku než (t2 − t1 ) . Jestliže oba výsledky funkce okamžité pohotovosti A(t ) jsou, až na povolenou chybu řešení, obdobné, lze učinit závěr, že výsledek analýzy je přesný. Pokud je rozdíl řešení funkce okamžité pohotovosti významný, doporučuje se zmenšit časový krok analýzy.
Strana: 65 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
5.3 Preventivní údržba pro systém stejně spolehlivý jako nový Systém, který je po poruše nahrazen (je stejně spolehlivý jako nový) má intenzitu poruch vztaženu k počátku provozu komponenty. Tím se ztrácí základní markovská vlastnost, že systém je bez paměti. Analyzovaný systém se řeší například metodou Monte Carlo, jejíž princip je uveden v kapitole 3.3.5. V kapitole 4.3.1 je uveden způsob výpočtu funkce okamžité nepohotovosti U (t ) pomocí metody Monte Carlo pro systémy, které mají dobu do přechodu popsanou pomocí Weibullova rozdělení. Přechody mezi stavy se vkládají pomocí procedur a funkcí do zdrojového textu programu. Základní typy procedur a funkcí pro popis doby dob přechodu jsou uvedeny ve stejné kapitole. Analyzovaný systém má přechodový diagram se shodnou topologií jako systém, který je stejně spolehlivý jako starý. Jestliže se v systému provádí preventivní údržba, nelze přechodový diagram systému popsat pomocí intenzity přechodu. Je nutné připojit další typ procedur, která zajistí časově nespojitý popis pravděpodobnosti stavů. U systému, který je stejně spolehlivý jako starý se přechody způsobené preventivní údržbou uskutečňují pomocí matice údržby U (t1 , t2 ) . V kapitole 5.2 jsou vymezeny dva základní případy. Buď jsou údržbovým zásahem v daném časovém intervalu < t1 , t2 ) pokryty všechny komponenty a pravděpodobnost, že údržba je proveden správně, je 1. Potom je prvek matice údržby uij (t1 , t2 ) = 1 . Nebo nejsou údržbovým zásahem v daném časovém intervalu < t1 , t2 )
pokryty všechny komponenty nebo pravděpodobnost, že údržba je
provedena správně je menší než 1. Potom je prvek matice údržby 0 ≤ uij (t1 , t2 ) < 1 . U systému, který je stejně spolehlivý jako nový se přechody způsobené preventivní údržbou uskutečňují pomocí vložených procedur a funkcí. Do hlavního programu se k příslušnému stavu i přidá následující procedura: x=rand(); if(x<0.18)
prvek matice údržby uij , obvykle je odlišná hodnota
t2=t+periodaudrzby1(t); volání procedury „periodaudrzby1“ if (t2<=t1) zjištění, zda dojde k přechodu t1=t2; state=2; stav j end end if((x>0.18)&(x<0.99))
prvek matice údržby uij
t2=t+periodaudrzby1(t); volání procedury „periodaudrzby1“ if (t2<=t1) zjištění zda dojde k přechodu t1=t2; state=1; stav k end end
Protože při modelování systému, který je stejně spolehlivý jako nový se vkládají procedury do hlavního programu, je nutné i prvky matice údržby uij (t1 , t2 ) vkládat pomocí Strana: 66 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
procedur. V ukázce procedury, která je uvedena výše, se nejdříve generuje náhodné číslo. Dále jsou do procedury vloženy podmínky (například if((x>0.18)&(x<0.99))), které určují do jakého nového stavu systém přejde. Nechť stav i , i ∈ 1,..., N , přechodového diagramu je dán kombinací komponent
k1 , k 2 ,...,k l ,..., k n ve stavu provozu a v poruše. Nechť se sestavuje procedura l -té komponenty, která popisuje její preventivní údržbu. Nechť stav i obsahuje r porouchaných objektů s označením l . Nechť jsou preventivní údržbou pokryty všechny porouchané komponenty s pravděpodobností p . Nechť systém přejde v čase t obnovou jedné komponenty ze stavu i do stavu j a následně obnovou dalších komponent do dalších stavů. Potom procedura, popisující preventivní údržbu, přechází s určenou pravděpodobností do nových stavů i , j ,…, za předpokladu, že 0 < p < 1 . Jestliže p = 1 , potom procedura přejde s jistotou do stavu j . Hodnota pravděpodobností těchto prvků je dána binomickým rozdělením s parametry Bi(r , p) . Procedura, která byla uvedena výše, volá funkci „periodaudrzby1“, která má následující strukturu: function [t3] = periodaudrzby1(t) t3=(fix(t/10000)+1)*10000-t;
zadání periody údržby 10 000 h
Strana: 67 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
5.4 Příklad V této kapitole je uveden vysvětlující příklad dvoustavového systému. Způsob řešení je obdobný jako u aplikačního příkladu, který je v příloze 4. Aplikační příklad se zabývá modelováním spolehlivosti kompresorové stanice tranzitního plynovodu RWE Transgas. Systém je tvořen dvěma stavy: • •
stav 0 - komponenta je v provozním stavu, stav 1 - komponenta je v poruchovém stavu. Intenzita poruch mezi stavem 0 a 1 je ve tvaru h = λ +
β ⋅ t β −1 , kde jednotlivé αβ
parametry jsou λ = 10−6 h-1, α = 2 ⋅ 106 a β = 1.5 . Intenzita oprav mezi stavem 1 a 0 je
⎡ β ⋅ t β −1 λ − − µ = 0,001 h . Matice intenzit přechodů je h(t ) = ⎢ αβ ⎢ µ ⎣
β ⋅ t β −1 ⎤ α β ⎥⎥ . Řešením −µ ⎦ soustavy diferenciálních rovnic (3.15b) lze stanovit funkci okamžité pohotovosti A(t ) . -1
λ+
Každých 10 000 h se všechny komponenty otestují. V případě, že jsou v provozním stavu jsou ponechány na stejné pozici. V druhém případě, kdy je zjištěno, že komponenta je v poruchovém stavu, je vyměněna. Pravděpodobnost, že porouchaná komponenta je testem rozpoznána a vyměněna, je 0,95. Dále pravděpodobnost, že neporouchaná komponenta je testem rozpoznána jako porouchaná, je nulová. Délka časového kroku je stanovena na 1 h. Matice pravděpodobností přechodů mezi
⎡ 1,5 ⋅ t 0,5 1,5 ⋅ t 0,5 ⎤ −6 −6 1 10 10 − − + stavy je pomocí vzorce P (t1 , t2 ) = I + h(t ) ⋅ (t2 − t1 ) = ⎢ (2 ⋅ 106 )1,5 (2 ⋅ 106 )1,5 ⎥⎥ , ⎢ 10− 3 0,999 ⎥⎦ ⎢⎣ t −t kde 2 1 = t . Hodnoty parametrů matice údržby U (t1 , t2 ) lze rozdělit pro dva případy. Za 2 prvé pro intervaly, které neobsahují časový okamžik k ⋅ 10000 h, kde k je celé kladné číslo. ⎡1 0⎤ V tomto případě je matice údržby U (t1 , t2 ) = ⎢ ⎥ . V druhém případě, pokud obsahují ⎣0 1 ⎦ časový okamžik k ⋅ 10000 , kde k je celé číslo. Příslušná matice údržby je 1 0 ⎡ ⎤ U(t1 , t2 ) = ⎢ ⎥ . Přechodový diagram systému je: ⎣0,95 0,05⎦ 0
µ10 u10
λ01
1 .
Pro stanovení funkce okamžité nepohotovosti systému U (t ) je využit model v softwarovém prostředku MATLAB®. Na obr. 9 (systém stejně spolehlivý jako starý) a obr. 10 (systém stejně spolehlivý jako nový) je uveden graf U (t ) a U (t1 , t2 ) systému, který je analyzován v této kapitole. Strana: 68 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Obr. 9: Graf nepohotovosti systému s preventivní údržbou, systém stejně spolehlivý jako starý
Obr. 10: Graf nepohotovosti systému s preventivní údržbou, systém stejně spolehlivý jako nový Z grafu funkce okamžité nepohotovosti U (t ) lze zjistit: •
hodnoty funkce okamžité nepohotovosti U (t ) ,
•
hodnoty součinitele střední nepohotovosti U (t1 , t2 ) ,
•
hodnotu maximální nepohotovosti systému U max (t1 , t2 ) .
Často je vhodné při analýzách systému sledovat součinitel střední nepohotovosti. Tento parametr lze snížit například prodloužením střední doby do poruchy komponent, zkrácením střední doby do obnovy, úpravou politiky údržby. Grafy nepohotovostí lze využít i pro analyzování vlivu periody preventivní údržby na celkovou nepohotovost systému. Graf funkce okamžité nepohotovosti může být také vstupem pro další analýzy, například pro ekonomickou analýzu nákladů životního cyklu. Ekonomické modelování je uvedeno v kapitole 7. Strana: 69 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
6. Výkonové konfigurace systému Mnoho zálohovaných zařízení může, podle požadavků, pracovat nejenom na maximální provozní výkon, ale i na odlišný výkon. Ke každé komponentě je možné přiřadit informaci, jaký podíl z maximálního možného výkonu systému je schopna sama realizovat. Změnou výkonu systému se následně může změnit také funkce okamžité pohotovosti A(t ) sledovaného systému. Z toho důvodu je nutné popsat systém odděleně pro plný výkon a snížený výkon, přičemž možných scénářů sníženého výkonu (výkonové konfigurace systému) je obvykle více. Lze předpokládat, že systém, který vyrábí nižší výkon, bude mít i vyšší pohotovost. Jako příklad systému s více výkonovými konfiguracemi lze uvést výrobní systém tvořený linkami. U tohoto systému se podle požadavků na objem výroby mohou jednotlivé linky postupně uvádět z nevyužitého stavu do provozu a naopak. Dalším příkladem mohou být tranzitní soustavy, např. silniční a železniční doprava, transport plynu nebo ropy, kdy je spolehlivost systému závislá na objemu přepravy. Úloha na zjištění parametrů spolehlivosti systému, při přepravě plynu tranzitním plynovodem je uvedena v kapitole 8 a příloze 4. Pro popis výkonových konfigurací systému je nutné, aby spolehlivost systému byla závislá na jeho výkonu. Protože se jedná typicky o zálohovaný systém, jsou při maximálním provozním výkonu obvykle všechny komponenty v provozu. Při nižším výkonu systému jsou v provozu buď všechny komponenty, které však pracují na nižší výkon, potom jsou zálohovány v aktivní záloze [67] - definice 15-02.17 Nebo jsou v provozu pouze některé komponenty a ostatní jsou v nevyužitém stavu připraveny na provoz v případě poruchy, potom je systém v pohotovostní záloze [67] - definice 15-03. 18 Pro každou výkonovou konfiguraci analyzovaného systému lze vytvořit blokový diagram bezporuchovosti systému (RBD diagram). RBD diagramy jsou v této disertační práci využívány z toho důvodu, že lze u nich relativně snadným způsobem zobrazit různé typy zálohování a tím pádem i různé výkonové konfigurace systému. Jestliže systém pracuje na menší výkon, potom může být jiný RBD model systému. Změny RBD diagramů závisí především na typu zálohování komponent, zda jsou v aktivní nebo v pohotovostní záloze. Výkonové konfigurace popisující systém mohou mít vlivem různého zálohování komponent odlišné hodnoty intenzit poruch λ (t ) / proudu poruch z (t ) . Následně pro každou výkonovou konfiguraci a získaného RBD diagramu se sestaví přechodový diagram systému. Přechody mezi stavy se ohodnotí intenzitou poruch λ (t ) (proudem poruch z (t ) ) a intenzitou oprav µ (t ) (vztaženou k počátku simulace nebo k okamžiku obnovy). Případně je možné ohodnotit přechody pravděpodobnostmi přechodů mezi stavy pij (t1 , t2 ) . Model se následně řeší v závislosti na tom, zda je systém stejně spolehlivý jako starý soustavou diferenciálních rovnic, nebo stejně spolehlivý jako nový procedurami a funkcemi popisující přechody mezi stavy. Daný model se následně řeší numericky.
17
Jako příklad lze uvést silniční sítě a elektronické komponenty v systému.
18
Jako příklad lze uvést výrobní linky, kdy v případě nižšího objemu výroby lze jednu linku vypnout.
Strana: 70 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
6.1 Přechod z RBD modelu k přechodovému diagramu V této kapitole jsou zobrazeny základní, nejjednodušší markovské diagramy pro paralelní a m z n zálohování komponent. Každý typ zálohování lze dále rozdělit na zálohy aktivní a pohotovostní ([67] def. 191-15-02, 191-15-03) . K řešení je využito článků [58] a [59]. Články se zabývají modelováním spolehlivosti elektrické sítě pomocí nehomogenních markovských procesů. Základem je vytvoření přechodového diagramu pro každou výkonovou konfiguraci. Každá výkonová konfigurace je propojena s ostatními, tím lze modelovat přechody z jedné do ostatních výkonových konfigurací. RBD diagramy lze převést na přechodové diagramy. Obráceně to neplatí. Všechny RBD diagramy hodnoceného systému, kterými jsou popsány různé výkonové konfigurace, se převedou na přechodové diagramy. Intenzity přechodů nejsou mezi ekvivalentními stavy různých výkonových konfigurací shodné, proto jsou v této kapitole uvedeny základní způsoby zálohování i s uvedením násobících koeficientů intenzit přechodů. Obdobně stupeň degradace (provozní, částečně poruchový, poruchový stav) příslušných stavů různých výkonových konfigurací není shodný. Mezi základní typy zálohování patří: • • • • •
sériové, paralelní - aktivní záloha, paralelní - pohotovostní záloha, m z n - aktivní záloha, m z n - pohotovostní záloha.
V přechodových diagramech na obr. 11, 12, 13, 14, 16, 18 a 20 se předpokládá, že všechny komponenty mají shodnou intenzitu poruch λ a intenzitu oprav µ . Provoz či částečně poruchové stavy jsou označeny bílou barvou, poruchové stavy červeně. Pokud není uvedeno jinak, stav 0 na následujících přechodových diagramech představuje 0 porouchaných z celkového množství komponent, stav 1 představuje jednu porouchanou komponentu atd. Dále budou popsány základní přechodové diagramy pro výše uvedené základní typy zálohování. Sériové zálohování Blokový diagram bezporuchovosti RBD sériového modelu tvořeného dvěma komponentami lze zobrazit podle obr. 11 vlevo. Porucha komponenty 1 převede model systému z provozního stavu 0 do poruchového stavu 1 s intenzitou poruch λ1 . V poruchovém stavu 1 může dojít pouze k obnově a přechodu do provozního stavu 0 s intenzitou oprav µ1 . Obdobně pro popis poruchy komponenty 2. Při poruše převede model systému z provozního stavu 0 do poruchového stavu 2 - intenzita λ 2 . Zpět do stavu 0 se systém dostane obnovou -
intenzita µ 2 . Obdobně lze rozšířit i pro n komponent. Provozní stav sériového systému předpokládá, že všechny komponenty jsou v provozu. Jestliže na některé komponentě vznikne porucha, je v poruchovém stavu celý systém. Přechodový diagram je zobrazen na obr. 11 vpravo.
Strana: 71 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
1
2
Obr. 11: RBD diagram a přechodový diagram sériového modelu pro dvě sériové komponenty s rozdílnou intenzitou oprav Jestliže n komponent má shodnou nebo podobnou intenzitu oprav, za předpokladu λi << µ i ∀i ∈ 1,.., n , lze stavy 1 až n spojit do jediného stavu 1. Získá se přechodový diagram systému, který je uveden na obr. 12. Výsledná intenzita poruch se získá λ = Intenzita oprav se stanoví ze vzorce µ =
0
∑λ i
∑λ . i =1
i
i
λ ∑i µi i µ
n
.
λ
1
Obr. 12: Přechodový diagram pro sériové komponenty s obdobnou intenzitou oprav Paralelní zálohování aktivní V následujících odstavcích je uveden způsob modelování komponent, které jsou paralelně a aktivně zálohované. Aktivní záloha znamená, že jsou všechny provozuschopné komponenty v provozu. Systém je tvořen n shodnými komponentami. Komponenty mají intenzitu poruch λ a intenzitu oprav µ . Jednotlivé typy stavů nastanou: • • •
provozní stav - všechny komponenty jsou v provozním stavu, částečně poruchový stav - některé, nikoliv však všechny, komponenty jsou v provozuneschopném stavu, poruchový stav - všechny komponenty jsou v provozuneschopném stavu.
Blokový diagram bezporuchovosti je zobrazen na obr. 13 . Porucha jedné komponenty převede systém v přechodovém diagramu z provozního stavu 0 do částečně poruchového stavu 1. Ve stavu 1 může dojít k poruše druhé komponenty a tím k přechodu do částečně poruchového stavu 2 nebo k obnově a přechodu do provozního stavu 0. Poruchou všech komponent přejde přechodový model do poruchového stavu n . Přechodový diagram je zobrazen na obr. 13.
Strana: 72 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
0
µ
nλ
1
µ
(n-1)λ
2
µ
λ
n
Obr. 13: RBD a přechodový diagram pro paralelně aktivně zálohované komponenty V následujících několika řádcích je uveden způsob jak modelovat poruchu se společnou příčinou. Na obr. 14 je zobrazen přechodový diagram systému, který je tvořen dvěma shodnými komponentami a je paralelně aktivně zálohován. Na komponentě se vyskytují dva módy poruch. První je porucha komponenty, která má intenzitu poruch λ1 . Druhá je porucha, která ohrozí chod obou komponent (například napájení), v tomto případě je intenzita poruch λ 2 .
Porucha, která ohrozí provoz jedné komponenty, způsobí změnu z provozního stavu 0 do částečně poruchového stavu 1 s intenzitou poruch 2λ1 , protože jsou v provozu 2 shodné
komponenty. Intenzita poruch ze stavu 1 do stavu 2 je λ1 , protože se může porouchat pouze zbývající jedna komponenta.
Při poruše se společnou příčinou vyřadí porucha více zálohovaných komponent z provozu. Intenzita vzniklé poruchy se společnou příčinou je λ2 . Intenzita poruch se společnou příčinou není závislá na počtu provozovaných komponent, ale pouze na intenzitě jevu, který vyřadí z provozu všechny komponenty. V uvedeném příkladě, tvořeném dvěma komponentami, je nutné k intenzitě přechodu ze stavu 1 do stavu 2 přičíst i intenzitu poruch se společnou příčinou λ 2 . Celková intenzita přechodu ze stavu 1 do stavu 2 je λ1 + λ2 . Intenzita opravy jakékoliv poruchy je µ .
Strana: 73 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Obr. 14: RBD a přechodový diagram pro dvě paralelně zálohované komponenty a poruchou se společnou příčinou Paralelní zálohování pohotovostní V této podkapitole je popsána tvorba přechodového diagramu pro paralelní pohotovostní zálohu. Při aktivní záloze jsou v provozu všechny komponenty. Při pohotovostní záloze jsou v provozu pouze ty komponenty, které jsou nutné pro výrobu. Náhradní komponenty se použijí až po poruše provozované komponenty. Systém je tvořen n shodnými komponentami. Komponenty mají intenzitu poruch λ a intenzitu oprav µ . Jednotlivé typy stavů nastanou: • • •
provozní stav - jedna komponenta je v provozu, ostatní jsou v provozuschopném stavu, částečně poruchový stav - některé, nikoliv však všechny, komponenty jsou v poruchovém stavu, poruchový stav - všechny komponenty jsou v poruchovém stavu.
Blokový diagram bezporuchovosti je zobrazen na obr. 15.
Obr. 15: RBD diagram pro paralelně pohotovostně zálohované komponenty Porucha jedné komponenty převede systém v přechodovém diagramu z provozního stavu 0 do částečně poruchového stavu 1. Ve stavu 1 může dojít k poruše druhé komponenty a tím k přechodu do částečně poruchového stavu 2 nebo k obnově a přechodu do provozního Strana: 74 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
stavu 0. Poruchou všech komponent přejde přechodový model do poruchového stavu n . Přechod ze stavu 0 do stavu 1 má okamžitou intenzitu poruch λ (t ) , protože v provozu je pouze jedna komponenta, ostatní jsou v pohotovostní záloze. Přechod ze stavu 1 do stavu 2 má okamžitou intenzitu poruch také λ (t ) , protože v provozu je jedna komponenta, jedna komponenta je v poruchovém stavu a ostatní jsou v pohotovostní záloze. Předpokládá se, že jednotlivé komponenty jsou přibližně stejně využívány pro provoz. Přechodový diagram systému je uveden na obr. 16. 0
µ
λ
1
µ
λ
2
µ
λ
n
Obr. 16: Přechodový diagram pro paralelně pohotovostně zálohované komponenty Při modelování paralelního pohotovostního zálohování je nutné předpokládat, že všechny komponenty jsou přibližně stejně využívány. Pokud tomu tak není, je nutné modelovat systém pomocí odlišného přechodového diagramu. Nechť paralelní systém je tvořen dvěma komponentami A a B. Komponenta A je v provozu, komponenta B je v nevyužitém stavu - pohotovostní záloze - stav 0. Poruchou komponenty A přejde systém v diagramu do stavu 1. Po poruše komponenty A mohou nastat dva případy. Buď bude komponenta B převedena z nevyužitého stavu do provozu (systém je v provozu) - stav 2. Nebo komponenta B nebude v provozu a dojde k poruše celého systému - stav 4. Po poruše komponenty B model systému přejde do stavu 3. Ze stavu 3 může model systému přejít buď do stavu 0, nebo do stavu 4. Blíže viz obr. 17.
Obr. 17: Markovský přechodový diagram pro dvě komponenty v pohotovostní záloze Sloučením stavů 0 a 2 do jednoho provozního stavu, a obdobně stavů 1 a 3 do částečně provozuschopného stavu, se získá přechodový diagram, který je uveden na obr. 18. Přechodové intenzity poruch mezi provozním a částečně poruchovým stavem zůstanou shodné, protože v provozu je v daném čase pouze jediná komponenta, druhá je mimo provoz a její intenzita poruch je v daném okamžiku nulová.
Strana: 75 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
m z n zálohování aktivní Při m z n aktivním zálohování je systém tvořen n shodnými komponentami, které mají intenzitu poruch λ a intenzitu oprav µ . Jednotlivé typy stavů nastanou: • •
provozní stav - všechny komponenty jsou v provozním stavu, částečně poruchový stav - některé komponenty jsou v poruchovém stavu a zároveň je více než m komponent v provozním stavu, poruchový stav - méně než m komponent je v provozním stavu.
•
Blokový a přechodový diagram je zobrazen na obr. 18. Porucha jedné komponenty převede systém v přechodovém diagramu z provozního stavu 0 do částečně poruchového stavu 1. Ve stavu 1 může dojít k poruše druhé komponenty a tím k přechodu do opět částečně poruchového stavu 2 nebo k obnově a přechodu do provozního stavu 0. Jestliže je v provozu méně než m komponent, systém přechází do poruchového stavu.
0
µ
nλ
1
µ
(n-1)λ
2
mλ
nm+1
µ (m-1)λ
λ
n
Obr. 18: Přechodový diagram pro m z n aktivně zálohované komponenty m z n zálohování pohotovostní Při m z n pohotovostním zálohování je systém tvořen n shodnými komponentami, kde je v provozu právě m komponent. Komponenty mají intenzitu poruch λ a intenzitu oprav µ . Jednotlivé typy stavů nastanou: • •
provozní stav - m komponent je v provozu, ostatní jsou v nevyužitém stavu, částečně poruchový stav - některé komponenty jsou v poruchovém stavu a zároveň je m komponent v provozním stavu, • poruchový stav - méně než m komponent je v provozním stavu. Blokový diagram bezporuchovosti je zobrazen na obr. 19.
Strana: 76 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Obr. 19: RBD diagram pro m z n pohotovostně zálohované komponenty Intenzita přechodu ze stavu 0 do stavu 1 je mλ , protože v každém čase je právě m komponent v provozu, ostatní jsou v nevyužitém stavu. Obdobně intenzita poruch ze stavu 1 do stavu 2 je mλ , protože je opět m komponent v provozu. Intenzita poruch z provozuneschopného stavu n − m + 1 do stavu n − m je pouze (m − 1)λ , protože počet komponent v provozu je m − 1 . Přechodový diagram je zobrazen na obr. 20. 0
µ
mλ
1
µ
mλ
2
mλ
nm+1
µ (m-1)λ
λ
n
Obr. 20: Přechodový diagram pro komponenty zálohované m z n s pohotovostní zálohou
Strana: 77 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
6.2 Sestavení matice pravděpodobností přechodů pro systém s více výkonovými režimy Funkce okamžité pohotovosti výkonové konfigurace systému, která je popsaná přechodovým diagramem majícím N stavů, lze vyřešit pomocí: • •
soustavy N obyčejných diferenciálních rovnic prvního řádu, pro systém stejně spolehlivý jako starý, metody Monte Carlo, pro systém stejně spolehlivý jako nový.
Obecně analyzovaný systém je charakterizován několika výkonovými konfiguracemi, které mají shodnou topologii přechodového diagramu, ale které se odlišují hodnotami okamžitých intenzit přechodů a definicí zda jednotlivé stavy jsou provozní, částečně poruchové nebo poruchové. Z okamžitých intenzit přechodů se vytvoří pro každou výkonovou konfiguraci: • matice pravděpodobnosti přechodů P (t1 , t2 ) - systém stejně spolehlivý jako starý, • procedury vložené do zdrojového textu výpočtového programu - systém stejně spolehlivý jako nový. Systém stejně spolehlivý jako starý Model systému může přecházet poruchou, obnovou či údržbou komponent v rámci jedné výkonové konfigurace mezi stavy podle pravidel přechodového diagramu. Přechody jsou dány maticí pravděpodobnosti přechodů P (t1 , t2 ) . Dále model systému může přecházet mezi odpovídajícími stavy různých výkonových konfigurací podle předem daných pravidel. Z tohoto důvodu je nutné, aby jednotlivé výkonové konfigurace měly shodnou topologii přechodového diagramu. Přechody mezi konfiguracemi mohou být buď časově spojité, tzn. že mohou nastat kdykoliv v čase. Tento systém lze popsat pomocí matice intenzit přechodů h(t ) a následně pomocí vzorce (4.1) převést soustavu na matici pravděpodobností přechodů P (t1 , t2 ) . Nebo mohou nastat pouze v určitých diskrétních časech a pak je nutné popsat systém pomocí matice pravděpodobnosti přechodů P (t1 , t2 ) . Systém s v různými výkonovými scénáři se řeší maticí pravděpodobnosti přechodů, která má v bloků kolem diagonály. Každý scénář s příslušnou výkonovou konfigurací reprezentuje jeden blok. Každý výkonový scénář reprezentovaný blokem je čtvercový a má tolik řádků a sloupců, jako je počet stavů přechodového diagramu každé výkonové konfigurace. Dále jsou v matici přechody mezi různými výkonovými konfiguracemi odpovídajících stavů. Tyto pravděpodobnosti přechodů jsou topologicky rovnoběžné s diagonálou matice pravděpodobnosti přechodů P (t1 , t2 ) . Na obr. 21 je ukázána struktura pravděpodobností matice přechodu s nenulovými prvky. Systém má pět stavů a vyskytuje se ve dvou výkonových konfiguracích. Červeně jsou vyznačeny prvky ij matice pravděpodobnosti přechodů, které mohou nabývat nenulových hodnot. Ostatní bílé prvky jsou vždy nulové.
Strana: 78 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Přechod mezi první a druhou výkonovou konfigurací pro stav 3
Intenzita přechodu mezi stavem tři a čtyři druhé výkonové konfigurace
Obr. 21: Struktura matice pravděpodobností přechodů pro dvě výkonové konfigurace Systém stejně spolehlivý jako nový Model systému, který je stejně spolehlivý jako nový, může přecházet v rámci jedné výkonové konfigurace mezi stavy podle přechodového diagramu. Přechody jsou dány pravidly, které jsou vkládány ve formě procedur do zdrojového textu programu (příloha 3). Dále může přecházet mezi různými výkonovými konfiguracemi. V tomto případě model systému přechází mezi odpovídajícími stavy různých výkonových konfigurací. Z tohoto důvodu je nutné, aby jednotlivé výkonové konfigurace měli shodnou topologii přechodového diagramu. Přechody mezi konfiguracemi mohou být buď časově spojité, nebo mohou nastat pouze v určitých diskrétních časech. V prvním případě, kdy dochází k časově spojitým přechodům, se přechod mezi odpovídajícími stavy různých výkonových konfigurací modeluje obdobně jako intenzity přechodů. Do hlavního programu je vložena procedura: t2=t+vykonkonfig; if (t2<=t1) t1=t2; state=j; end
která volá funkci: function t3=vykonkonfig definování
jména
funkce
vykonkonfig1
a
výstupní
proměnné t3 t3=exprnd(10000); zjištění, za jakou dobu by došlo k přechodu mezi
konfiguracemi end
V druhém případě, kdy dochází k přechodu pouze v určitých diskrétních časech, se přechody mezi odpovídajícími stavy různých výkonových konfigurací modelují obdobně jako preventivní údržba. Do hlavního programu je vložena procedura: x=rand(); if(x<0.18) t2=t+vykonkonfig1(t); if (t2<=t1) t1=t2;
generování náhodného čísla, určující do jaké výkonové konfigurace systém přejde změna do jedné výkonové konfigurace volání procedury vykonkonfig1 zjištění, zda dojde k přechodu Strana: 79 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
state=2; end end if((x>0.18)&(x<0.99)) t2=t+vykonkonfig1(t); if (t2<=t1) t1=t2; state=1; end end
stav j změna do jiné výkonové konfigurace volání procedury vykonkonfig1 zjištění, zda dojde k přechodu stav k
Procedura volá funkci periodaudrzby, která má následující strukturu: function [t3] = vykonkonfig1(t) t3=(fix(t/10000)+1)*10000-t; zjištění doby přechodu, perioda změny 10 000 h
Strana: 80 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
6.3 Příklad Podnik má 3 výrobní linky. Na ranní a odpolední směně jsou v provozu dvě linky a na noční směně je v provozu pouze jedna. Při poruše jakékoliv provozované linky je k dispozici v pohotovostní záloze okamžitě další linka. Pravděpodobnost, že při najetí linky vznikne porucha, je nulová. Obdobně pravděpodobnost výpadku více než jedné linky je nulová. Intenzita poruch výrobní linky je stanovena na λ = 10 −6 h-1 a intenzita oprav je µ = 0,001 h-1. Na obr. 22 jsou zobrazeny dva RBD diagramy. Vlevo je diagram ranní a odpolední směny, kdy jsou v provozu 2 výrobní linky. Vpravo je diagram noční směny, kdy je v provozu pouze 1 výrobní linka. 1
1
2 ze 3
2
1 ze 3
2
3
3 Obr. 22: RBD diagram řešeného příkladu
Pomocí postupu uvedeného v této kapitole se získá přechodový diagram, který je na obr. 23. Pro ranní a odpolední směnu dojde k poruše, když jsou dvě nebo tři linky v poruše. Při noční směně dojde k poruše při výpadku všech tří linek současně. 0
µ
2λ
1
µ
2λ
2
µ
λ
3
0
µ
λ
1
µ
λ
2
µ
λ
3
Obr. 23: Přechodové diagramy řešeného příkladu Na obr. 24 je zobrazen graf nepohotovosti systému, který je popsán výše. Nespojitost funkce je způsobena přechody mezi nepohotovostí jednotlivých konfigurací. Z výsledků lze vypočítat součinitel střední pohotovosti, který reálně popisuje skutečnou nepohotovost systému.
Strana: 81 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Detail dolní části grafu
Obr. 24: Funkce okamžité nepohotovosti řešeného příkladu provozovaného ve dvou odlišných výkonových režimech
Strana: 82 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
7. Další možné aplikace 7.1 Modelování způsobilosti systému Podle normy [67] je způsobilost definována: Způsobilost (cit. [67] def. 02-04) - je schopnost objektu plnit požadavky na služby s danými kvantitativními charakteristikami při daných vnitřních podmínkách. Hlavní nevýhodou ukazatelů pohotovosti tkví v tom, že stavy v přechodovém diagramu jsou buď provozní, částečně poruchové, nebo poruchové 19. Výsledná funkce okamžité pohotovosti A(t ) systému je dána součtem pravděpodobností (popisujících systém) přes všechny provozní a částečně poruchové stavy. V praktickém případě při zjišťování ukazatelů pohotovosti se předpokládá, že provozní a částečně poruchový stav má plnou funkčnost, tzn. že případná porucha se neprojeví na plánované výrobě. Poruchový stav představuje poškození takového rozsahu, že případná porucha se projeví na plánované výrobě. Způsobilost jednotlivých stavů představuje podíl objemu skutečné výroby (dané příslušným stavem v přechodovém diagramu) k objemu požadované výroby. Pro účely disertační práce jsou definovány následující ukazatele způsobilosti: Funkce okamžité způsobilosti C (t ) - ekvivalent pro pohotovost funkce okamžité pohotovosti - Pravděpodobnost, že objekt je ve stavu schopném plnit požadavky na služby s danými kvantitativními charakteristikami při daných vnitřních podmínkách. N
C (t ) = ∑ Ci ⋅ pi (t )
(7.1)
i =1
Součinitel střední způsobilosti C (t1 , t2 ) - ekvivalent pro pohotovost součinitel střední pohotovosti - Střední hodnota funkce okamžité způsobilosti v daném časovém intervalu (t1 , t 2 ) . t
1 2 C (t1 , t2 ) = C (τ )dτ t2 − t1 t∫1
(7.2)
Celková způsobilost CT (0, t ) - ekvivalent pro kumulovanou dobu použitelného stavu TUT (t ) . 19
Stav provozu předpokládá, že všechny komponenty jsou v provozním stavu a jejich provozuschopnost není snížena. Částečně poruchový stav předpokládá, že existuje porucha některé komponenty (některých komponent), která nemá vliv na provozuschopnost celého systému. Poruchový stav předpokládá, že existuje porucha některé komponenty (některých komponent), která má vliv na provozuschopnost celého systému. Strana: 83 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
t2
CT (t1 , t 2 ) = ∫ C (τ )dτ
(7.3)
t1
Stavy provozu a částečně poruchové mají způsobilost stavů Ci = 1 . Stavy, které jsou poruchové mají způsobilost stavů omezenou 0 ≤ Ci ≤ 1 . Vstupními parametry pro výpočet funkce
okamžité
způsobilosti
C (t )
jsou
způsobilosti
jednotlivých
stavů
Ci , dále
pravděpodobnosti stavů pi (t ) , které se získají při výpočtu úloh pohotovosti.
7.2 Modelování zisků a ztrát [63], [64], [66], [81], [82] Zdrojové texty, které jsou uvedené v příloze 2 a 3, umožňují ocenit finanční částkou (ziskem nebo ztrátou) každý stav k a každý přechod (daný poruchou, obnovou nebo údržbou) ze stavu i do stavu j . Každý provozní stav, poruchový stav, porucha, obnova a údržba v sobě nese finanční zisk nebo ztrátu, kterou lze uvést v měnové jednotce. Výsledkem finančního modelování je určení celkového zisku v daném časovém intervalu pro daný zvolený scénář. Pomocí vytvořených scénářů popisujících různé chování systému v čase lze zjistit výsledné ekonomické a spolehlivostní parametry jako funkci času. Měnící se parametry mohou být například následující: • • • •
změna počtu zálohovaných komponent, jiná perioda při preventivní údržbě, změna koncepce údržby nebo změna střední doby do obnovy, dřívější nebo pozdější výměna komponent.
Pro každý vytvořený scénář, odlišující se od ostatních změnou výše uvedených vstupních parametrů, lze stanovit výslednou pohotovost systému a tím i celkové náklady (zisky) CostT (t1 , t 2 ) za sledovaný čas analýzy. Pomocí vytvořených scénářů lze optimalizovat, z pohledu ukazatelů spolehlivosti, náklady životního cyklu (LCC) komplexních systémů. Dalším uplatněním metody je její použití při optimalizování údržby systémů metodou údržby zaměřenou na bezporuchovost (RCM). Při zjišťování spolehlivostně ekonomických ukazatelů jednotlivých scénářů je nutné mít na paměti, že: • •
Zvyšující se spolehlivostí systému obvykle rostou náklady na provoz, údržbu a bezpečnostní prvky systému. Pravděpodobnost, že nastanou finančně náročné poruchy nebo poruchy s následky na zdraví, je však nižší. Snižující se spolehlivostí systému obvykle rostou náklady na údržbu po poruše, stejně tak za nedodání výrobků v požadovaných termínech. Zároveň pravděpodobnost, že nastanou finančně náročné poruchy nebo poruchy s následky na zdraví, je vyšší. Nižší jsou však náklady na provoz, preventivní údržbu a bezpečnosti prvky systému. Celkový zisk CostT (t1 , t2 ) v časovém intervalu t1 ,t2 , za předpokladu, že t2 − t1 → 0 ,
se stanoví: Strana: 84 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
N
CostT (t1 , t2 ) = ∑ pi (t1 ) ⋅ Costi ⋅ (t2 − t1 ) + i =1
N
∑ pi (t1 ) ⋅ hij (t1 ) ⋅ Costij ⋅ (t2 − t1 ) +
i , j =1 i≠ j
N
∑ p (t ) ⋅ u
i , j =1 i≠ j
i
1
ij
(t1 , t2 ) ⋅ Costuij
nebo N
N
N
i =1
i , j =1 i≠ j
i , j =1 i≠ j
CostT (t1 , t2 ) = ∑ pi (t1 ) ⋅ Costi ⋅ (t2 − t1 ) + ∑ pi (t1 ) ⋅ pij (t1 , t2 ) ⋅ Costij + ∑ pi (t1 ) ⋅ uij (t1 , t2 ) ⋅ Costuij (8.1), kde
CostT (t1 , t2 )
celkový zisk v intervalu t1 ,t2
Costi
zisk za jednotku času ve stavu i
Costij
zisk za každý přechod ze stavu i do stavu j
Costuij
zisk za každý preventivní údržbový zásah ze stavu i do stavu j
pi (t ) hij (t )
pravděpodobnost stavu i , prvek pravděpodobnostního vektoru p(t )
uij (t1 , t2 )
pravděpodobnost preventivní údržby ze stavu i do stavu j v intervalu
intenzita přechodu ze stavu i do stavu j , prvek matice h(t )
(t1 , t 2 ) , prvek matice údržby U (t1 , t2 ) pravděpodobnost přechodu ze stavu i do stavu j v intervalu (t1 , t 2 ) ,
pij (t1 , t2 )
prvek matice pravděpodobnosti přechodů P (t ) Pomocí vzorce (8.2) lze stanovit součinitel středního zisku Cost (0, t ) . t
Cost (0, t ) =
1 Cost (τ )dτ t ∫t
(8.2)
Jako důkaz nutnosti provádět ekonomické úvahy při modelování údržby po poruše a preventivní údržby je uvedena následující citace z [12] str. 178: „Vzhledem k tomu, že dosud nebyly provedeny žádné seriozní kvantitativní studie, objasňující vzájemné souvislosti (ekonomické důsledky) mezi preventivní a nápravnou údržbou, je obtížné „optimalizovat“ rozsah pouze preventivní údržby bez uvážení důsledků pro provozní poruchy a jejich opravy. Proto každá seriózní optimalizace musí být postavena na pevných principech a to: 1. Buď na technicko - ekonomické optimalizaci nákladů, spojených s pořízením objektu a jeho udržením v provozuschopném stavu za celý jeho technický život; 2. Nebo na posouzení provozních, technických, bezpečnostních nebo ekonomických důsledků, které mohou mít pro uživatele jednotlivé provozní poruchy, důvodně očekávané u objektu během celé doby jeho provozu.“ Obdobně důležité je modelování vícestavových systémů a zjišťování parametrů spolehlivosti pro nekonstantní intenzity přechodů.
Strana: 85 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
8. Analýza spolehlivosti kompresorové stanice tranzitního plynovodu Při řešení projektu Akademie věd číslo 1ET401940412 s názvem Modelování a kvantifikace spolehlivosti dynamických systémů bylo jedním z cílů vytvořit a následně řešit dynamický model spolehlivosti kompresorové stanice tranzitního plynovodu RWE Transgas. Tento model, který je popsán dále, využívá přínosů disertační práce, které jsou uvedeny v kapitolách 4, 5, 6, 7. Pro řešení je použito softwarového prostředku uvedeného v přílohách 2 a 3. V rámci projektu AV byla vytvořena analýza spolehlivosti kompresorové stanice a přilehlých potrubních linií tranzitního plynovodu [112], která využívá všech výše uvedených přínosů disertační práce. Její zestručněná verze je uvedena v příloze 4. Projekt AV č. 1ET401940412 byl úspěšně obhájen.
Strana: 86 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
9. Otevřené problémy disertační práce •
Kapitola 3.4 - Statistická rozdělení popisující dobu do poruchy komponenty V této kapitole jsou uvedeny základní způsoby zjištění bodových odhadů exponenciálního a Weibullova rozdělení. Dále je naznačeno, že v případě kdy nelze použít pro popis doby do poruchy ani exponenciální ani Weibullovo rozdělení je nutné hledat složitější statistické modely. V programech, které jsou uvedeny a popsány v příloze 2 a 3 řešící funkci okamžité nepohotovosti U (t ) systému se předpokládá, že intenzity přechodu mezi stavy jsou ve tvaru hij (t ) = λ +
βt β −1 . αβ
Možné rozšíření softwarového prostředku je především v tom, že bude řešit dobu do přechodu komponent, které nejsou popsány pomocí superpozice exponenciálního a Weibullova rozdělení. V softwaru by se tato úloha řešila buď vložením upravených volaných funkcí, pro systém stejně spolehlivý jako nový - viz diskuze v kapitole 4.3. Nebo pro systém stejně spolehlivý jako starý upravením zdrojového textu programu. • •
Kapitola 4.1 - Modelování opravovaných a neopravovaných objektů Kapitola 5.1 - Modelování údržby pomocí markovské analýzy V této kapitole jsou popsány rozdíly v intenzitách přechodu pro systémy, které jsou stejně spolehlivé jako starý a stejně spolehlivé jako nový. Pro oba modely byly vytvořeny zdrojové texty v prostředí Matlab, které řeší funkci okamžité nepohotovosti systému. V disertační práci není řešen případ, kdy výsledná intenzita přechodu po obnově komponenty je odlišná od intenzit přechodů systémů modelů stejně spolehlivého jako nový a stejně spolehlivého jako starý. •
Kapitola 5.1 - Modelování údržby pomocí markovské analýzy Zdrojový text programu, který je uveden v příloze 2 a 3, předpokládá, že preventivní údržba je s periodickými intervaly. V případě, že nedochází k údržbovým zásahům v periodických intervalech, potom není možné pomocí přiloženého softwaru úlohu řešit. Obdobně není řešen případ, kdy k údržbě dochází v náhodných časech. •
Kapitola 7.2 - Modelování zisků a ztrát U řady analyzovaných systémů nemusí být finanční ohodnocení přechodů mezi stavy, způsobené poruchou, obnovou či údržbou konstantní. Software, který je uveden v příloze 2 a 3 však neumožňuje zadávat ohodnocení stochasticky podle předem daných pravidel. Stejně tak nelze ohodnotit přechod, kdy finanční částka je funkcí doby výpadku. To je potřebné, jestliže markovská analýza má být vstupní či pomocnou analýzou pro metodu RCM. • •
Příloha 2 - Popis softwaru pro řešení systému stejně spolehlivého jako starý Příloha 3 - Popis softwaru pro řešení systému stejně spolehlivého jako nový Software, který vypočítává funkci okamžité nepohotovosti systému byl naprogramován pro účely disertační práce. V případě jeho dalšího využívání je vhodné vytvořit více komfortní verzi. •
Příloha 4 - Analýza spolehlivosti kompresorové stanice tranzitního plynovodu Z důvodu rozsahu disertační práce nejsou v příloze uvedeny grafy, které se zabývají způsobilostí systému a finančním ohodnocením stavů a přechodů mezi stavy. Tyto grafy jsou uvedeny v práci [113]. Strana: 87 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
10.
Závěr
V disertační práci byly v rámci markovské analýzy ve spolehlivosti dosaženy následující hlavní výsledky: 1. sestavení modelu neexponenciálního rozdělení doby do přechodu, 2. sestavení modelu preventivní údržby, 3. sestavení modelu výkonových konfigurací systému. Ad 1: Použití markovské analýzy předpokládá, že doba do přechodu mezi stavy splňuje exponenciální rozdělení. Tento předpoklad je velmi omezující pro tvorbu analýz. V kapitole 4 jsou řešeny dva základní modelové případy. • V prvním případě je intenzita poruch po obnově komponenty shodná jako u komponenty, u které k poruše nedošlo - systém je stejně spolehlivý jako starý (as good as old). Při řešení byla sestavena soustava obyčejných diferenciálních rovnic s nekonstantními parametry a následně vypočtena numericky. • V druhém případě je intenzita poruch po obnově je shodná jako u nové komponenty systém je stejně spolehlivý jako nový (as good as new). Úloha se řeší pomocí metody Monte Carlo, kdy jednotlivé přechody z přechodového diagramu jsou do hlavního programu vkládány pomocí funkcí a procedur. Případy, kdy intenzita poruch po obnově komponenty je mezi těmito dvěma limitními případy, není v disertační práci řešena. Ad 2: Při řešení obvyklé úlohy markovskou analýzou se sestaví soustava diferenciálních rovnic, která se následně řeší. Výsledkem jsou spojité funkce okamžité pohotovosti A(t ) systému. V kapitole 5 je uveden prostředek, který připojuje k markovské analýze preventivní údržbu. Tento prostředek má za následek, že funkce okamžité pohotovosti není spojitá. Řešení tohoto problému je pak závislé na tom, zda je systém stejně spolehlivý jako starý či jako nový. V prvním případě se sestaví diferenciální rovnice. V čase, kdy dochází k preventivní údržbě je pravděpodobnostní vektor p(t ) vynásoben maticí údržby U(t , t ) , čímž se získá nová hodnota vektoru p(t + u ) . Nová hodnota vektoru p(t + u ) charakterizuje hodnotu vektoru po údržbě. V druhém případě se sestaví procedury, které jsou součástí řešícího softwaru. Tyto procedury zajistí přechod z jednoho stavu do dalšího. Software neřeší úlohu, kdy dochází k preventivní údržbě, která však není periodická. Ad 3: V kapitole 6 je vytvořen prostředek, který umožňuje modelovat výkonové konfigurace systému. Tento prostředek je využitelný, jestliže je systém zálohován a jeho spolehlivost je závislá na objemu výroby - například při modelování spolehlivosti tranzitních sítí. Princip je založen na vytvoření několika přechodových diagramů popisujících různý objem výroby. Tyto přechodové diagramy mají shodnou topologii, ale jsou odlišné v intenzitách přechodu a definici provozuschopnosti stavů. Všechny uvedené výsledky jsou implementovány do zdrojového textu naprogramovaného v programu MATLAB®. Zdrojové texty jsou přílohou disertační práce. Strana: 88 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
V kapitole 8 je řešen aplikační příklad - modelování kompresorové stanice tranzitního plynovodu. V této kapitole je využito vše, co bylo teoreticky popsáno v kapitolách 4 až 7. Na systému je definováno 10 různých výkonových konfigurací. Dle charakteru analyzovaného systému se předpokládá, že komponenty mají neexponenciální dobu do poruchy a mají navíc preventivní údržbu. Bylo vytvořeno a řešeno mnoho scénářů, které popisovaly chování systému v závislosti na změně: • • • •
výkonové konfigurace, střední doby do poruchy, typu rozdělení doby do poruchy, preventivní údržby.
Modelování kompresorové stanice tranzitního plynovodu RWE Transgas bylo jedním z dílčích cílů úspěšně obhájeného projektu AV 1ET401940412 s názvem Modelování a kvantifikace spolehlivosti dynamických systémů.
Strana: 89 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Conclusion The following results were concluded using the Markov analysis in dependability: 1. construction of model with non-exponential distribution of time to transition 2. construction of maintenance model with forward intervals 3. construction of model of system performance configuration 1) The use of the Markov analysis supposes that the time to transition between two postures fulfil the conditions of exponential distribution. This hypothesis is very limited in analysis construction. There are two basic examples solved in chapter 4 • In first case the failure rate after the component recovery is equally as limiting as the rate by component, which was fault-free. The system is “as good as old”. The mathematical solving of this case was realized after the construction of an ordinary differential equation system with non-constant parameters. • In the second case the failure rate after component recovery is the same as the rate by component, which was new. System is “as good as new”. The Monte Carlo method can be used to solve this case, where each transition is inserted from a transition diagram to the main program with the help of functions and procedures. In such cases where the failure rate between two causes is limited, this thesis does not attempt to make a solution. 2.) The system of differential equation is constructed through the ordinary tasks of the Markov analysis. The results of it are continuous function of actual system emergency A(t ) . Forward intervals as an instrument, which extends Markov analysis is mentioned in chapter 5. This instrument states that the function of an actual system emergency is not continuous. The solving of this problem depends on the fact that the system is „good as old“ or „good as new“. In first cause differential equations are constructed. In the time of forward intervals maintenance is dependability vector p(t ) multiplied by maintenance array U (t , t ) , which causes a new value of vector p(t + u ) . This new value characterises the vector value after maintenance. In second cause the procedures, which are components of solving software are constructed. These procedures ensure transition from one to another posture. The software doesn’t solve the task, when the forward intervals maintenance is not periodic. 3.) In chapter 6 an instrument is mentioned, which allows modelling of systems performance configuration. This instrument is useful, when the system is backed up and its dependability depends on the production volume - for example modelling of transit networks. The principle of modelling is based on the creation of some transition diagrams, which describe different production volume. These transition diagrams have the same topology, but different intensity rates and definition of postures serviceability. All the mentioned results are implemented into source text through the software, MATLAB®. Source texts are mentioned in appendix of this thesis. In chapter 8 the application task is solved - modelling of transit networks compressor station. In this task all benefits were applied, which were theoretically described in chapters 4 - 7. It was defined 10 performance configuration of system. On the base of the analyzed system Strana: 90 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
character is supposed, that failure rates of components are non - exponential and also has maintenance with forward intervals. Many cases were created and solved, which described the behaviour in context of the following changes. • • • •
performance configuration, average failure rate, distribution type of failure rate, maintenance with forward intervals.
Modelling of a compressor station of a gas pipeline RWE Transgas was one of the fractional aims of the project AV 1ET401940412 (modelling and quantification of dynamic systems dependability), which was successfully upheld.
Strana: 91 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
11.
Literatura
Knihy, články a časopisy [1]
Briš R., Inovační metody pro ocenění spolehlivosti prvků a systémů, Monografie 1. vydání, Ostrava, VŠB-Technická univerzita Ostrava, 2007, ISBN 978-80-248-1596-1
[2]
Kořenář V., Stochastické procesy, Praha, Vysoká škola ekonomická v Praze, 1998, ISBN 80-7079-813-0
[3]
Walter J., Stochastické modely v ekonomii, Praha, SNTL, 1970, 196 s.
[4]
Walter J., Stochastické procesy, Praha, SPN, 1983
[5]
Škrášek J., Tichý Z., Základy aplikované matematiky III, Praha, SNTL, 1990
[6]
Tichý Z. a kol., Aplikovaná matematika I a II, Praha, SNTL, 1978, 2386 s.
[7]
Šor J. B., Statističeskije metody analiza i kontrolja kačestva i nadějnosti, Radio Moskva 1962
[8]
Šor J. B., Praktické problémy teorie spolehlivosti, ČVUT Praha, 1967
[9]
http://www.weibull.com/systemrelwebcontents.htm
[10]
MIL-HDBK 338B, Electronic reliability design handbook, 1998
[11]
Schüller J. C. H., Brinkman J. L., Van Gestel P. J., van Otterloo R. W, Methods for determinining and processing probabilities „Red Book“, Publication Series on Dangerous Substances 4,The Commitee fort he Prevention of Disasters by Hazardous materials, Nederland, 2005
[12]
Fuchs P., Vališ D., Metody analýzy a řízení rizika, TUL, Liberec, 2004 http://flow.kmo.tul.cz/~www/czech/predmet.php?id=13&t=soubory [22. 2. 2008]
[13]
Fuchs P., Využití spolehlivosti v provozní praxi, TUL, Liberec, http://flow.kmo.tul.cz/~www/czech/predmet.php?id=13&t=soubory [22. 2. 2008]
[14]
Vališ D., Software pro podporu spolehlivosti, materiály Odborné skupiny pro spolehlivost 27. setkání, ČSJ Praha, 2007
[15]
MIL-HDBK-217F, Reliability prediction of electronic equipment, 1991
[16]
http://src.alionscience.com/spidr/
[17]
Vintr M., Současné přístupy k predikci bezporuchovosti prvků, ČSJ Praha, Perspektivy jakosti, Vol. 2 (2005), No.3, pp. 27-31, ISSN 1214-8865
[18]
Vintr M., Přehled metod a nástrojů pro odhad bezporuchovosti prvků, ČSJ Praha, 27. setkání Odborné skupiny pro spolehlivost, pp. 26-57
[19]
Matějček, Katalog Tesla, k. p., 1985
[20]
System analysis reference: Reliability, Availability and Optimization
spolehlivosti
elektrotechnických
a
elektronických
2004
výrobků,
http://www.weibull.com/systemrelwebcontents.htm [21]
Vesely W. E., Goldberg F. F., Fault tree Handbook, U.S. Nuclear Regulatory Commission, NUREG 492
[22]
http://www.event-tree.com/ Strana: 92 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
[23]
http://www.eventtreeanalysis.com/
[24]
Barlow R.E., Proschan F., Statistical Theory of Reliability and Life Testing. Probabilistic Models, New York, Holt, Rinehart and Winston, 1975
[25]
Billinton R., Allan R.N., Reliability Evaluation of Engineering Systems. Concepts and Techniques. Second Edition, New York, Plenum Press, 1992
[26]
Birolini A., Quality and Reliability of Technical Systems. Theory - Practise Management, Berlin, Springer Verlag, 1994
[27]
Lisnianski A., Levitin G., Multi-state system reliability, World Scientic Publishing Co., 2003, ISBN 981 238-306-9
[28]
Levitin G., The Universal Generating Function Optimalization, Springer, ISBN 185 233-927-6
[29]
Tůma J., Rusek S., Martínek Z., Chemišinec I., v elektroenergetice, ČVUT Praha, ISBN 80-239-6483-6
[30]
Rektorys K., Přehled užité matematiky, Prométheus, Praha 1995, ISBN 80-7196-179-5
[31]
Vitásek E., Numerické metody, SNTL 1987
[32]
Ralston A., Rabinowitz P., A first course in numerical analysis, Dover Publications 2001, ISBN 0-486-41454-X
[33]
Fishman G. S., Monte Carlo: Concepts, Algorithms and Applications, Volume 1 od Springer Series in Operations Research, Springer-Verlag, New York, 1996
[34]
Virius M., Aplikace matematické statistiky: metoda Monte Carlo, ČVUT Praha 1998
[35]
Fabian F., Kluiber Z.: Metoda Monte Carlo a možnosti jejího uplatnění, Prospektrum s.r.o., Praha 1998
[36]
System analysis reference: Reliability, Availability and Optimization
in
Reliability Goňo
Analysis
R.:
and
Spolehlivost
http://www.weibull.com/systemrelwebcontents.htm [37]
Hammersley J.M., Handscomb D.C., Monte Carlo methods, Methuen &Co Ltd., 1965
[38]
Using Weibull++ 7 and Monte Carlo simulation for probabilistic design, Reliability HotWire č. 56, 2005 http://www.weibull.com/hotwire/issue56/hottopics56.htm
[39]
Wilkins J. D.: The bathtub curve and produkt failure behavior, Reliability HotWire, č. 21 a 22, 2002 http://www.weibull.com/hotwire/issue21/hottopics21.htm http://www.weibull.com/hotwire/issue22/hottopics22.htm
,
[40]
Jiang R., Murthy D.N.P.: The exponentiated Weibull family: A graphical approach, IEEE Transactions on Reliability, 1999, p. 68-72
[41]
http://www.weibull.com/GPaper/index.htm
[42]
http://www.weibull.com/lifedatawebcontents.htm
[43]
Farnum N.R., Booth P., Uniqueness of maximum likelihood estimators of the 2parameter Weibull distribution, IEEE Transaction on Reliability, 1997, p. 523-525
[44]
Likelihood ratio konfidence bounds, Reliability http://www.weibull.com/hotwire/issue42/relbasics42.htm
HotWire,
č.
42,
2004
Strana: 93 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
[45]
Shimokawa T., Liao M.: Goodness-of-fit tests for type-I extreme-value and 2-parameter Weibull distributions, IEEE Transactions on Reliability, 1999, p. 79-86
[46]
Praks P., Bacarizo H. F., Labeau P. E., On the modeling of ageing using Weibull models: Case studies, Safety, Reliability and Risk Analysis - ESREL 2008, p. 559-565, ISBN 13-978-0-415-48513-5
[47]
Hsieh H. K.: Prediction intervals for Weibull observations, based on early failure data, IEEE Transactions on Reliability, 1996, p. 666-670
[48]
Klutke G. A., Kiessler P.C., Wortman M.A., A Critical look at the bathtub curve, IEEE Transactions on Reliability, 2003 p. 125-129
[49]
Lai C. D., Xie M., Murthy D. N. P., A Modified Weibull Distribution, IEEE Transtions on Reliability, Vol. 52 No. 1, 2003, p. 33-37
[50]
Ng H. K. T., Parameter Estimation for a Modified Weibull Distribution for Progressively Type-II Censores Samples, IEEE Transactions on Reliability, Vol. 54 No. 3, 2005, p. 374-380
[51]
Vintr Z., Co lze najít na internetu o spolehlivosti, Univerzita obrany Brno, 2005
[52]
Using Weibull++ 7 and Monte Carlo simulation for probabilistic design, Reliability HotWire č. 56, 2005 http://www.weibull.com/hotwire/issue56/hottopics56.htm
[53]
www.itemsoft.com
[54]
www.relex.com
[55]
www.isograph-software.com
[56]
www.sohar.com
[57]
www.rmclogan.co.uk
[58]
Platis A, Limnios N., Du M. L., Performability of electric-power systems modeled by non-homogeneous Markov chains, IEEE Transactions on reliability, 1996, p. 605-610
[59]
Kyandoghere K., Reliability modeling & performance of variable link-capacity networks, IEEE Transactions on reliability, 1998, p. 44-45
[60]
Cassady C.R., A Generic Model of Equipment Availability Under Imperfekt Maintenance, IEEE Transactions on reliability, 2005, p. 564-571
[61]
Sheu H.S., Lin Y.B., Liao G.L., Optimum policies for a system with general imperfect maintenance, Reliability engineering & system safety, 2006, p. 362-369
[62]
Zequeira R.I., Bérenguer C., Periodic imperfect preventive maintenance with two categories of competing failure modes, Reliability engineering & system safety, 2006, p. 460-468
[63]
Financial applications for Weibull analysis, Reliability HotWire, č. 14, 2002 http://www.weibull.com/hotwire/issue14/hottopics14.htm
[64]
Moubray J. M.: Reliability-centered Maintenance, New-York,1999
[65]
Svoboda A. a kolektiv, Plynárenská příručka - 150 let plynárenství v Čechách a na Moravě, GAS s.r.o, Praha 1997, ISBN 80-902339-6-1, 1308 stran
[66]
Fuchs P., Údržba zaměřená na bezporuchovost, 21. setkání OSS, Praha 2005 Strana: 94 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Technické normy [67]
ČSN IEC 50(191) (01 0102), Mezinárodní elektrotechnický slovník, kapitola 191: Spolehlivost a jakost služeb, Federální úřad pro normalizaci a měření, 1993
[68]
ČSN EN 13306 (01 0660) Terminologie údržby, ČNI Praha 2002
[69]
ČSN IEC 1078 (01 0677):1993 Metody analýzy spolehlivosti metoda blokového diagramu bezporuchovosti. ČSN Praha, 1993
[70]
ČSN EN 61078 (01 0677):2007, Techniky analýzy spolehlivosti - Blokový diagram bezporuchovosti a booleovské metody, ČSN Praha 2007
[71]
ČSN IEC 1025 (01 0676):1993, Analýza stromu poruchových stavov, ČNI Praha, 1993
[72]
ČSN EN 61025(01 0676):2007, Fault tree analysis, 2006
[73]
ČSN IEC 60300-3-1 (01 0690), Management spolehlivosti - Část 3-1: Pokyn k použití Techniky analýzy spolehlivosti - Metodický pokyn, 2003
[74]
ČSN IEC 1165 (01 0691):1996, Použití Markovových metod, ČNI Praha, 1996
[75]
ČSN EN 61165 (01 0691):2007 Použití Markovových metod. ČNI Praha, 2007
[76]
ČSN IEC 60605-6 (01 0644-6) Zkoušení bezporuchovosti zařízení - Část 6: Testy platnosti předpokladu konstantní intenzity poruch nebo konstantního parametru proudu poruch, ČNI 1998
[77]
ČSN IEC 60605-4 (01 0644-4) Zkoušení bezporuchovosti zařízení - Část 4: Statistické postupy pro exponenciální rozdělení - Bodové odhady, konfidenční intervaly, předpovědní intervaly a toleranční intervaly, ČNI 2002
[78]
ČSN EN 61164 (01 0647) Růst bezporuchovosti - Metody statistických testů a odhadů, ČNI Praha 2005
[79]
ČSN IEC 61609 (01 0649), Elektronické součástky - Bezporuchovost - Referenční podmínky pro intenzity poruch a modely namáhání pro přepočty
[80]
ČSN IEC 61649 (01 0653) Testy dobré shody, konfidenční intervaly a dolní konfidenční meze pro data s Weibullovým rozdělením, ČNI Praha 1998
[81]
ČSN IEC 60300-3-11 (01 0690-3-11) Management spolehlivosti - Část 3-11: Návod k použití - Údržba zaměřená na bezporuchovost, ČNI Praha, 1999
[82]
ČSN EN 60300-3-3 (01 0690-3-3) Management spolehlivosti - Část 3-3: Pokyn k použití - Analýza nákladů životního cyklu, ČNI Praha 2005
Strana: 95 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
Publikační činnost autora Konference [83]
Chudoba J., Evaluation of dependability by using Markov analysis, International Congress ECMS 2005 Electronics, Controle, Modeling, Measurement and Signále, Toulouse 2005, Evaluation of dependability by using Markov analysis
[84]
Chudoba J.,What is possible to find out from railroad akcident statistics, konference Věda a krizové situace, TUO Ostrava 2005, ISBN 80-248-0944-3
[85]
Chudoba J., Zpracování dat o poruchách výrobků z provozu, konference CQR, Hejnice 2006
[86]
Praks P., Bris R., Chudoba J.: Motivace k modelování spolehlivosti kompresorové stanice tranzitního plynovodu s přihlédnutím k dynamickému chování dispečera, Sborník abstraktů semináře aplikované matematiky SAMO’06, TUO Ostrava 2006, ISBN 80-248-1151-0
[87]
Chudoba J., Analysis of the probability of car accidents, Reliability, safety and diagnostics of transport structures and means 2005, Pardubice 2005, ISBN 80-7194-769-5
[88]
Chudoba J.: Hodnocení přesnosti výsledků z metody FMECA, In: Request ’06, Ed. Centrum pro jakost a spolehlivost výroby, ČVUT Praha 2007, ISBN 978-80-01-03709-6
[89]
Praks. P., Chudoba J., Briš R.: Modelování spolehlivosti kompresorové stanice tranzitního plynovodu užitím spojitých markovských procesů, Ed. Centrum pro jakost a spolehlivost výroby, ČVUT Praha 2007, ISBN 978-80-01-03709-6
[90]
Chudoba J.: Determination of cause for increased failure rate of electrical components, 9th international scientific konference Applied Mechanics 2007, TUO Ostrava 2007, ISBN 978-80-248-1389-9
[91]
Praks P., Chudoba J., Briš R., Koucký M.: Reliability analysis of a natural gas compression station and surrounding gas pipeline network with assuming of performance changes by a dispatcher, In: Proceedings of the European Safety and Reliability Conference 2007 (ESREL 2007). Ed. Terje Aven&Jan Erik Vanen, London: Tailor&Francis Group, 2007, ISBN 978-0-415-44786-7 Chudoba J., Model spolehlivosti kompresorové stanice tranzitního plynovodu s neexponenciální dobou do poruchy komponent, Věda a krizové situace 2007, Ostrava, 23.10.2007, ISBN 978-80-7385-011-1
[92]
[93]
Chudoba J., Model spolehlivosti kompresorové stanice tranzitního plynovodu s užitím markovských procesů, Opotřebení spolehlivost diagnostika 2007, Brno, 2007, ISBN 978-80-7231-294-8
[94]
Chudoba J., Praks P.: Reliability analysis of a natural gas compression station and probability estimation of non-delivery of agreed amount of the gas, Současnost a budoucnost krizového řízení 2007, Praha 2007
[95]
Chudoba J., Metody spolehlivosti a jejich využití v ekonomii, In. VII. Mezinárodní konference studentů doktorských studijních programů IMEA, Ed. Kocourek Aleš, Technická univerzita v Liberci, 2008, ISBN 978-80-7372-331-6
[96]
Chudoba J., Modelování spolehlivosti kompresorové stanice tranzitního plynovodu metodou markovských procesů za využití metody Monte Carlo, In.13. vedeckej Strana: 96 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
konference Riešenie krízových situácií v specifickém prostredí, Ed. Fakulta špeciálneho inžinierstava Žilinskej univerzity v Žiline, ISBN 978-80-8070-846-7, 2008 [97]
Chudoba J., Případová studie využitelnosti modifikovaného Weibullova rozdělení při popisu doby do poruchy, konference Request 2008, Brno 2008
[98]
Chudoba J., Reliability analysis of a natural gas compression station, prediction of a profit and loss arise from estimation of non-delivery of agreed amount of the gas, In. 26th International Conference Mathematical Methods in Economics 2008, Ed. Kocourek Aleš, Technická univerzita v Liberci, 2008 ISBN 978-80-7372-387-3
[99]
Chudoba J., Reliability Modelling of Gas Line Compression Station by Markov Process Method and by further Apllication of Monte-Carlo Method, Fakulta bezpečnostního inženýrství, TUO Ostrava, 2008 Publikované články
[100]
Chudoba J., Modelování spolehlivosti průmyslové soustavy, Odborná skupina pro spolehlivost, 28. setkání, Česká společnost pro jakost, Praha 2007, ISBN 978-80-0201965-7 Technické zprávy
[101]
Ságl P., Chudoba J., Rozbor degradace typů komponent v rámci M3-5, ČEZ, a.s., JE Dukovany, celek 1, zpráva Technické univerzity v Liberci, 2007
[102]
Chudoba J., Posouzení spolehlivosti a důsledků změny intervalu zkoušení pojistných ventilů rafinérie Kralupy, zpráva Technické univerzity v Liberci, 2005
[103]
Chudoba J., Stanovení pravděpodobnosti nehody při transportu osobním automobilem, zpráva Technické univerzity v Liberci, 2004
[104]
Chudoba J., Ságl P., Určení příčin zvýšené poruchovosti signálek řady IS, KIS v jaderné elektrárně Dukovany, zpráva Technické univerzity v Liberci, 2007
[105]
Fuchs P., Chudoba J., Řešení problematiky spolehlivosti vstřikovacích forem Cadence Innovation při zvýšeném počtu zdvihů, zpráva Technické univerzity v Liberci, 2007
[106]
Fuchs P., Vališ D., Chudoba J., Kamenický J., Zajíček J., Bezporuchovost a životnost Techniky analýzy bezporuchovosti, skripta Technické univerzity v Liberci, 2006
[107]
Fuchs P., Vališ D., Chudoba J., Kamenický J., Zajíček J., Řízení spolehlivosti, elektronický učební text k předmětu Řízení jakosti a spolehlivosti, Liberec, 2006
[108]
Chudoba J., Dynamické modelování kritičnosti komponent a zařízení složitých průmyslových soustav - Komplexní rešerše přístupů k vytváření a využívání pravděpodobnostně orientovaných spolehlivostních modelů složitých technologických soustav s dynamickými prvky, zpráva Technické univerzity v Liberci, 2006
[109]
Chudoba J., Alternativní výpočet optimálního počtu náhradních dílů systému kontroly a řízení Jaderné elektrárny Dukovany, zpráva Technické univerzity v Liberci, 2005
[110]
Chudoba J., Dynamické modelování kritičnosti komponent a zařízení složitých průmyslových soustav - Příprava a zpracování komplexní informace o aktuálním stavu a možnostech využití metodologie markovských procesů při řešení otázek dynamického modelování spolehlivosti, zpráva Technické univerzity v Liberci, 2006
[111]
Chudoba J., Metodický pokyn k softwaru Hodnocení spolehlivosti - DPD, zpráva Technické univerzity v Liberci, 2007 Strana: 97 z 98
Technická univerzita v Liberci Disertační práce Ing. Josefa Chudoby
[112]
Chudoba J., Dynamický model spolehlivosti kompresorové stanice tranzitního plynovodu, zpráva Technické univerzity v Liberci, 2007
[113]
Chudoba J., Modelování dynamické spolehlivosti užitím stochastických procesů, disertační teze, Technická univerzita v Liberci, 2008 Software
[114]
Balatka M., Chudoba J., Intepor-intenzita poruch, verze 2.0 (poslední aktualizace 3. 10. 2005), Technická univerzita v Liberci Vedení ročníkových projektů a bakalářských prací Obhájené
[115]
Bezr J., Určení pohotovosti vstřikovacího stroje v Cadence Innovation k.s., bakalářská práce, Technická univerzita v Liberci, 2007
[116]
Dolina J., Statistické a spolehlivostní vyhodnocení reálné komponenty, ročníkový projekt, Technická univerzita v Liberci, 2007
[117]
Polák J., Určení pohotovosti vstřikovací formy v Cadence Innovation k.s., bakalářská práce, Technická univerzita v Liberci, 2007 Zvolené ročníkové projekty a bakalářské práce pro školní rok 2008 až 2009
[118]
Nguyen B. H., Vytvoření slovníku spolehlivosti, ročníkový projekt
[119]
Šarman A., Tvorba formuláře na sběr dat po poruchách, bakalářská práce
[120]
Taneček M., Zpracování dat o poruchách z provozu RWE Transgas, bakalářská práce
[121]
Teubelová H., Statistické zjištění vztahů mezi kurzy akcií a indexem PX na Burze cenných papírů, ročníkový projekt
[122]
Trégl P., Vytvoření softwaru zajišťující zjištění parametrů neexponenciálních rozdělení, bakalářská práce
[123]
Váňa M., Zjištění použitelnosti modifikovaného Weibullova rozdělení, bakalářská práce
Strana: 98 z 98