Analýza spolehlivosti tlakové nádoby metodou Monte Carlo Jakub Nedbálek
Abstrakt: Cílem práce je ukázat možnost využití Monte Carlo simulace pro studium úloh z oblasti spolehlivosti. V našem případě máme k dispozici systém, ve kterém se deterministicky vyvíjí procesní proměnná tlak. Časový průběh proměnné vychází z diferenciálních rovnic, jejichž řešení je nám známé. Všechny změny, resp. poruchy komponent považujeme za stochastické události. Vzhledem k dynamice procesu se snažíme vyčíslit pravděpodobnost poruchy vztaženou ke vstupním parametrům a vytvořit funkci pravděpodobnosti poruchy v čase. Pro komplikovanější typy úloh se numerické řešení pomocí metody Monte Carlo zdá být efektivní, nicméně jistou nevýhodou může být časová náročnost řešení díky výpočetní době a konstruování příslušného algoritmu. V naší konkrétní úloze byl čas, potřebný k modelování 2 .106 simulací, menší než 1000 s, což je zanedbatelná hodnota.
1.
Úvod
V této úloze je konkrétně nasimulováno a vyšetřeno chování uzavřeného systému, ve kterém se mění tlak v závislosti na čase. Cílem modelování je určit spolehlivost této soustavy vzhledem ke komponentám, jenž mají za úkol udržet tlak v patřičných mezích. Prostředkem k řešení je metoda Monte Carlo (MC), jako alternativa k analytickému postupu v [1].
2. Popis benchmark procesu Mějme k dispozici systém, v němž je hlavní procesní proměnou x tlak. Předpokládejme, že její počáteční hodnota v čase t0 je x0 a dále že se x vyvíjí exponenciálně. Dosáhle-li tlak úrovně x = l, začne působit komponenta C1, která může s pravděpodobností p0 selhat, s pravd. 1- p0 –p1 pracovat korektně, nebo se spustit částečně chybně s pravd. p1. Poslední dva případy vedou ke snížení x v systému k počáteční hodnotě x0. V případě, že x > L > l, systém považujeme za porouchaný. Situace odpovídá případu, kdy tlak v nádobě je regulován ventilem, který se při hodnotě tlaku l má otevřít, jinak dojde k roztržení nádoby - tlak vzroste nad bezpečnou mez L. Systém je doplněn komponentou C2, která akceleruje růst x, díky čemuž se systém dostane rychleji na úroveň x = l s následným zapojením komponenty C1 podle scénáře podobného předchozímu. Stavovou veličinu je možno popsat následujícími dif. rovnicemi:
a x, i = 1,4,6 dx ={ i ai > 0 ∀i , − ai x, i = 2,3,5 dt
(1)
kdy platí a2 > a3, a4 = a1 + b, a5 = a2 - b, a6 = - a3 + b, a b = 0.15 Intenzitu poruch λ považujme konstantní v celém zledovaném rozmezí < x0,L>. Pro koeficienty dosazujeme následující numerické hodnoty: x0 l
1 3
p0 p1
0.02 0.04
a1 a2
0.2 0.25
a4 a5
0.35 0.1
L
4
λ
4.10-2
a3
0.1
a6
0.05
Pro řešení úlohy předpokládáme, že pro danou možnost se stavová proměnná vyvíjí deterministicky a k případným změnám v průběhu vývoje dochází stochasticky. Protože naším úkolem není rozebírat postup výpočtu (1), podívejme se rovnou, jak vypadá řešení diferenciálních rovnic s příslušnými konstantními koeficienty a1 – a6 :
Obr. 1. Řešení (1) – vývoj stavové proměnné pro jednotlivé koeficienty
3. Řešení Monte Carlo metodou Úkolem je nalézt funkci pravděpodobnosti poruchy systému. Základem řešení naší testovací studie bude sestavení optimálního algoritmu MC metody. Za prostředí, v němž budeme případ simulovat, zvolme Matlab. Při sestavování algoritmického řešení vycházíme z následujících informací: - vzrůst x nad L - porucha - pokles x na x0 - úspěch - vzrůst x nad l - změna koeficientu a, změna hraničních časů a úrovní - hraniční úrovně x=xf, resp. časy t=tf zjistíme z rovnice x(t)=xe.exp(ae.(t-te)) – řešení vztahu (1) – kde xe a te jsou počáteční (inicializační) hodnoty v simulačním cyklu - uvažujeme, že při jedné simulaci můžou nastat obě poruchy u C1 i C2 – vybíráme dřívější čas
-
-
-
čas T1 je náhodně generován z exponenciálního rozložení s parametrem λ - odpovídá přechodu C2 do poruchy a změně ai o b. pravděpodobnosti p jsou generovány náhodně a platí: s pravděpodobností p0 přejde vývoj do stavu a1, 1- p0 –p1 do stavu a2 a p1 do a3. počáteční hodnoty pro každý cyklus simulace jsou dány te = 0 , xe = x0 = 1 a ai = a1 = 0.2
Algoritmus v Matlabu vypadá následovně: function simstav_x global X por usp m
% globalni promenne
lambda = 0.04; beta = 0.15; %vstupni koeficienty %ze zadani X1 = 1; X2 = 3; X3 = 4; A = [ 0.2 -0.25 -0.1 ]; p0 = 0.02; p1 = 0.04; % -----------------------------------------------------ms = input('pocet milionu simulaci = '); if isempty(ms) ,error('neudano'), end if ms~=round(ms) | ms<1 ,error('neprirozene'), end doba = zeros(ms*1000000,1) ;% doby trvani % odpredu porucha % odzadu uspech MS = length(doba) ; % pocet simulaci usp = 0 ;% pocet uspechu por = 0 ;% pocet poruch while
usp+por < MS
%cyklus pres pocet simulaci
% provedeni jedne simulace T1 = -log(rand)/lambda ; %doba fungovani prvku C2, %generovano nahodne
if elseif else end
nah = rand ; nah < p0 , nah < p1 , , %
an = A(1) ; an = A(3) ; an = A(2) ; funkce prvku C1
% %provadeni dilcich prechodu pocinaje pocatecnim stavem % tf = 0; xf = X1; af = A(1); %pocatecni inicializace while 1 %dokud neskonci jedna simulace te = tf ; xe = xf ; ae = af ; %dosazeni "hranicniho" casu if T1>te , tf_T = T1 ; else , tf_T = Inf ; end %kdyby pokles|vzrust pres "hranicni" uroven L if ae>0 %nalezeni nejblizsi vyssi urovne pro akci "vzrust" if xe < X2 , xf_X = X2 ; CO = 2 ; else , xf_X = X3 ; CO = 3 ; end else %nalezeni nejblizsi nizsi urovne pro akci "pokles" xf_X = X1 ; CO = 1 ; end %vypocet odpovidajiciho casu pro prislusny stav x tf_X = te + log(xf_X/xe)/ae ; %rozliseni co nastane drive (porucha od C2/zasah C1) if tf_T < tf_X %realizace "prechodove" akce == porucha C2 %dojde k akceleraci rustu tf = tf_T ; an = an + beta ; xf = xe*exp(ae*(tf-te)) ; af = ae + beta ; else % % realizace "vzrustove|poklesove" akce % switch CO case 1 % konec- uspech jak = +1 ; cas = tf_X ; break case 2 % akce zasah C1 tf = tf_X ; xf = xf_X ; af = an ;
case 3 %konec- porucha jak = -1 ; cas = tf_X ; end end %rozliseni co nastane
break
end %provadeni dilcich prechodu (while 1) %dosazeni vysledku jedne simulace %zapis cas uspechu/poruchy if jak>0 , usp = usp + 1 ; doba(MS+1-usp) = cas ; else , por = por + 1 ; doba( por) = cas ; end end %provadeni vsech simulaci % %prezentace simulacniho reseni X = sort(doba(1:por)) ; m = length(X) ; Y = (por/(por+usp))*(0:m)/m ;
%delka vektoru %pro kazde m vypocti %pravdepodobnost
plot([0 X'],Y,'r-' ) title([num2str(por) '+' num2str(usp) ' histories']) %popisy os a grafu ylabel('failure probability') xlabel('time')
Výsledkem je graf pravděpodobnosti poruchy systému. 59260+1940740 histories 0.03
0.025
failure probability
0.02
0.015
0.01
0.005
0
0
5
10
15
20
25
30
35
40
45
time
Obr. 2. Grafický výstup pro 2 .106 simulací MC. (nahoře – poč.poruch + poč.úspěchů)
4. Závěr Testovací studie byla provedena pro 2 .106 simulací MC metodou. Pro jiný počet simulací v daném řádu se tvar výsledné funkce příliš nemění. Výstup byl konfrontován s výsledkem [1] , kde bylo uvažováno řešení analyticky. Matematický průběh pravděpodobnostní funkce zjištěné MC metodou je v dobré shodě s analytickým výsledkem - MC simulaci můžeme tedy považovat za dostatečně spolehlivou pro naši benchmarkovou studii. Pro 2 .106 cyklů zabralo modelování na počítači PII 0,5GHz, 256 MB RAM čas t<1000s, což je zanedbatelná hodnota.
Literatura [1]
Labeau, E.P.: Pressurisation test-cases for dynamic reliability, Université Libre de Bruxelles, Brussels, 2002
[2]
Virius, M.:Základy výpočetní techniky (Metoda Monte Carlo), ČVUT, Praha, 1985
[3]
Briš R., Praks P.: Special Case of Dynamic Reliability Analysis Based on Time Dependent Acyclic Graph, The International Symposium on STOCHASTIC MODELS in RELIABILITY, SAFETY, SECURITY and LOGISTICS, February 15-17, 2005 Beer Sheva, Israel, ISBN 9984-668-79-7, str. 69-70
Adresa autora: Ing. Jakub Nedbálek, Vysoká škola báňská – Technická univerzita Ostrava, fakulta elektrotechniky a informatiky, katedra aplikované matematiky, Tř.17.listopadu 15, 708 33 Ostrava-Poruba. e-mail:
[email protected] Tato práce byla vytvořena v rámci projektu MŠMT 1M06047 - CQR.
