IV. ročník celostátní konference SPOLEHLIVOST KONSTRUKCÍ Téma: Posudek - poruchy - havárie 23.až 24.4.2003 Dům techniky Ostrava
151 ISBN 80-02-01551-7
SIMULAČNÍ POSUZOVÁNÍ SPOLEHLIVOSTI PŘI KORELOVANÝCH VELIČINÁCH Jaroslav Menčík Abstract The paper explains the work with correlated variables in reliability assessment using Monte Carlo simulation technique. After explanation of basic concepts, the relation between correlation and linear regression is shown. Then, expressions are given for generation of two and several linearly correlated quantities. Similar procedure can also be used for nonlinear correlation. The paper is complemented by practical recommendations.
1. Úvod Při vyšetřování chování konstrukcí simulační technikou Monte Carlo počítáme sledovanou veličinu y = f(x1, x2, x3,…) opakovaně pro náhodně generované hodnoty vstupních veličin x1, x2, x3… Z výsledků velkého množství těchto „pokusů“ můžeme sestrojit histogram, který informuje o rozmezí výskytu hodnot y, pravděpodobnosti překročení určitých mezí či vzniku poruchy apod. Výpočty jsou nejjednodušší, jestliže jednotlivé vstupní veličiny jsou navzájem nezávislé (např. modul pružnosti materiálu a plocha průřezu nosníku nebo tíha konstrukce a zatížení větrem). Takovéto případy jsou nejčastější. Mezi některými veličinami však může existovat určitá závislost, a to buď těsná (funkční), jako mezi rozměry průřezu válcované tyče a momentem setrvačnosti, anebo volná, kdy spolu se změnami jedné veličiny se více či méně mění i hodnoty druhé veličiny (např. objemová hmotnost a pevnost betonu, modul pružnosti a pevnost betonu, složky zatížení větrem v různých směrech apod.). V takovém případě hovoříme o korelační nebo statistické závislosti. Zvláštním případem je tzv. autokorelace, kdy náhodná veličina v určitém místě závisí do jisté míry na hodnotách téže veličiny v sousedních bodech nebo v předcházejících časových okamžicích. Příkladem jsou vlastnosti betonu nebo základové půdy, anebo teploty konstrukce, které se den ode dne (i během dne) mění, ale do značné míry závisí také na tom, o jaké roční období se jedná. Pokud ke korelacím při výpočtech nepřihlížíme, můžeme se dopustit chyby. Například k náhodně vygenerovaným velmi nízkým hodnotám modulu pružnosti betonu by počítač mohl generovat velmi vysoké hodnoty pevnosti, což neodpovídá skutečnosti. Při zohlednění korelací odrážejí výpočty skutečnost věrněji a předpovědi jsou přesnější, s menším rozptylem. Kromě toho někdy nemáme k dispozici hodnoty vstupní veličiny, kterou bychom potřebovali, ale pouze hodnoty jiné veličiny, která je s „naší“ veličinou korelována. Například přímé měření meze kluzu nosných částí existující ocelové konstrukce není dost dobře možné; lze ji však přibližně zjistit z vtiskové tvrdosti. O práci s korelovanými nebo autokorelovanými veličinami při posuzování spolehlivosti stavebních konstrukcí existuje řada publikací. Použité postupy jsou však zpravidla zmíněny jen rámcově nebo užívají speciální software. V řadě případů ale vystačíme s velmi jednoduchými postupy, které lze snadno aplikovat i při použití univerzálního programu, jako je např. AntHill [1]. V tomto příspěvku vysvětlíme základní pojmy a ukážeme některé jednoduché postupy pro zahrnutí korelací do výpočtů Monte Carlo. Jaroslav Menčík, Prof. Ing., CSc., Univerzita Pardubice, Dopravní fakulta Jana Pernera, Katedra provozní spolehlivosti, diagnostiky a mechaniky v dopravě, Studentská 95, 532 10 Pardubice, tel.: (+420) 466 036 529, e-mail:
[email protected] .
152
2. Základní pojmy Vzájemnou proměnlivost dvou veličin xi, xj charakterizuje tzv. kovariance, která je pro výběr o rozsahu N definovaná jako [2, 3] 1 N cov ( xi x j ) = (1) ∑ (xi ,K − mi )(x j ,K − m j ) N − 1 K =1 kde K znamená K-tou dvojici veličin xi a xj, a mi, mj jsou jejich aritmetické průměry. Kovariance závisí na jednotkách příslušných proměnných. Lepší představu o těsnosti vazby dvou veličin dává tzv. koeficient korelace, definovaný jako
rij = cov(xixj) /(sisj) (2) kde si, resp. sj je (výběrová) směrodatná odchylka veličiny xi, resp. xj. Korelační koeficient r je bezrozměrný a může nabývat hodnot od –1 do +1. Při r = 0 vzájemná souvislost neexistuje, zatímco r = +1 nebo –1 znamená deterministickou (funkční) závislost; pro r > 0 hodnoty xj s růstem xi rostou, při r < 0 klesají. Za nevýznamnou se obvykle považuje korelace s r < 0,2, za dosti těsnou korelace s r > 0,8. Tři ukázky jsou na obr. 1; některé příklady korelovaných veličin u stavebních konstrukcí, včetně hodnot koeficientů korelace, lze najít v [4, 5].
Obr.1: Příklady souborů se stejnou střední hodnotou a směrodatnou odchylkou a s rozdílnou mírou korelace. Podrobnosti viz kapitola 4. Koeficient korelace (1) je vhodný při normálním rozdělení základních souborů. V některých případech se používá tzv. Spearmanův korelační koeficient, který je neparametrický a navíc méně citlivý na náhodné vybočení některého bodu. Každé hodnotě veličiny xi (a podobně i xj) se přiřadí pořadové číslo podle její velikosti (nejmenší hodnotě odpovídá 1), a Spearmanův korelační koeficient se stanoví jako rS = 1 – [(6ΣdK2)/(N3 – N)] , (3) kde dK je rozdíl mezi pořadovými čísly dvojic hodnot xi,K, xj,K (K = 1, 2,…N), a N je počet těchto dvojic [2]. Obě definice korelačních koeficientů dávají podobné hodnoty r. Kovariance a korelace může existovat i mezi více veličinami, x1x2, x1x3, x2x4 apod. Jednotlivé hodnoty můžeme sestavit do tabulek, resp. matic: c11 c12 ... c1n r11 r12 ... r1n c c ... c r r ... r 21 22 2n 2n [c ] = , [ r ] = 21 22 , ( 4a,b ) ............ ............ cn1 cn 2 ... cnn rn1 rn 2 ... rnn kde cij = cov(xixj), i, j, = 1, 2, …n, kde n je počet proměnných. Z definice (1) je zřejmé, že cii odpovídá rozptylu s2 dané veličiny. Složky rij korelační matice se počítají podle vztahu (2); platí
153
přitom rii = 1. Protože korelační i kovarianční matice jsou symetrické kolem hlavní diagonály, tzn. rij = rji, cij = cji, zapisuje se často jen polovina matice pod diagonálou. Výpočet součinitelů či matic kovariance a korelace z naměřených hodnot nečiní problém, neboť je součástí mnoha programů (viz např. „ „Nástroje“ – „Analýza dat“ v Excelu). Doporučuje se však vždy otestovat, zda vypočtené hodnoty skutečně vyjadřují zákonitou souvislost mezi oběma veličinami a nejsou jen důsledkem náhodnosti výběru o omezeném rozsahu (viz dále).
3. Korelace a lineární regrese Kovariance i korelace, definované vztahy (1) a (2), souvisí s lineární regresní závislostí. Například hodnotu veličiny xj můžeme vyjádřit z veličiny xi jako (5) xj = aji + bji xi , přičemž platí (6) bji = cji /si2 = rji × (sj /si) Při korelaci můžeme také vyjádřit xi jako funkci xj, tzn. xi = aij + bijxj; nyní pro směrnici regresní přímky platí bij = cij/sj2 = rij × (sj/si); rovněž platí bijbji = rij2. Zatímco v prvém případě (vztah 5) při výpočtu regresních konstant aji, bji považujeme veličinu xi za nezávisle proměnnou a xj za závisle proměnnou, ve druhém případě je tomu naopak. Korelační koeficient vyjadřuje těsnost rozložení jednotlivých bodů kolem regresní přímky. Podobně jako můžeme vzdálenost bodu xj od průměru mj vyjádřit jako součet vzdálenosti xj od regresní přímky a vzdálenosti odpovídajícího bodu na regresní přímce od průměru, můžeme celkový rozptyl veličiny xj napsat jako součet rozptylu hodnot regresní funkce kolem průměru mj a rozptylu jednotlivých hodnot xj kolem regresní čáry (tzv. reziduální rozptyl), (7) sj2 = sj,reg2 + sj,res2. 2 Čtverec korelačního koeficientu rji (tzv. koeficient determinace) vyjadřuje, jaká část celkové proměnlivosti xj je způsobená proměnlivostí veličiny xi: rji2 = sj,reg2/sj2 = 1 – sj,res2/sj2 (8) Reziduální rozptyl můžeme tedy vyjádřit z celkového rozptylu jako sj,res2 = sj2 (1 – rji2) . (9) To znamená, že čím je korelační koeficient blíže k jedné, tím je rozptyl kolem regresní přímky menší ve srovnání s celkovým rozptylem, a tím větší chyby bychom se dopouštěli, kdybychom při simulačních výpočtech vzájemnou korelaci nezohlednili – zejména při větších sklonech b regresní přímky.
4. Generování dvou náhodných korelovaných veličin Při generování dvou korelovaných veličin x1, x2 využijeme skutečnosti, že x2 můžeme vyjádřit jako součet hodnoty regresní funkce x2,reg = f(x1) vypočtené pro náhodnou hodnotu prvé veličiny, a náhodné složky charakterizující rozptyl x2 kolem této funkce. Postup sestává ze dvou kroků. V prvém vygenerujeme náhodnou hodnotu veličiny x1, a ve druhém vypočítáme hodnotu x2 jako x2 = f(x1) + ∆x2 . (10) Je-li rozdělení jednotlivých hodnot x2 kolem regresní čáry normální, můžeme použít vztah ∆x2 = u2 s2,res , (11) kde u2 je náhodná hodnota veličiny s normovaným normálním rozdělením. S přihlédnutím k (9) můžeme psát .(12) x2 = f(x1) + u2 s2 (1 – r2)1/2 Jestliže mezi x1, x2 existuje lineární regresní závislost (5) a veličina x1 má normální rozdělení, můžeme při generování užít vztahy
154
x1 = m1 + u1s1 , (13a) 2 1/2 x2 = m2 + u1s1 b + u2 s2 (1 – r ) , (13b) resp. x2 = m2 + u1s2 r + u2 s2 (1 – r2)1/2, (13c) kde u1 a u2 jsou dvě náhodná čísla s normovaným normálním rozdělením. [Ve vztahu (13b) bylo využito výrazů f(x1) = m2 + b(x1 – m1) a x1 – m1 = u1s1; úprava (13c) pak vznikla užitím vztahu (6).] Při N opakováních těchto kroků dostaneme soubory dvou normálně rozdělených náhodných veličin s průměry m1, m2, směrodatnými odchylkami s1, s2 a korelačním koeficientem r. Jednotlivé generované dvojice hodnot x1, x2 potom dosazujeme do simulačních výpočtů vyšetřované veličiny y = y(x1, x2). Podobně by šlo postupovat i v případě, kdy rozdělení hodnot x2 kolem regresní funkce je jiné než normální. Na obr.1 jsou tři soubory korelovaných veličin x1, x2 generovaných podle vztahů (13a, c) pro normální rozdělení s parametry µ1 = 100, σ1 = 30, µ2 = 700, σ2 = 130. Zvolené hodnoty koeficientu korelace byly: ρ = 0 (obr. 1a), ρ = 0,3 (obr. 1b) a ρ = 0,9 (obr. 1c). Pro generované soubory po 200 bodech byly výběrové charakteristiky: a) m1 = 101,44; s1 = 31,78; m2 = 699,86; s2 = 131,08; r = – 0,004; b) m1 = 101,44; s1 = 31,78; m2 = 701,48; s2 = 131,76; r = 0,309; c) m1 = 101,44; s1 = 31,78; m2 = 704,80; s2 = 136,93; r = 0,904; tedy blízké parametrům základních souborů. (Při vyšších počtech generovaných hodnot byla shoda ještě lepší.)
5. Generování většího počtu navzájem korelovaných veličin Náhodné hodnoty pro n normálně rozdělených a navzájem korelovaných veličin x1, x2, … xn dostaneme prostřednictvím následujícího předpisu [6, 7]: x1 = k11u1 + m1 x2 = k12u1 + k22u2 + m2 x3 = k13u1 + k23u2 + k33u3 + m3 (14) ……………… xn = k1nu1 + k2nu2 + k3nu3 + … + knnun + mn , kde u1, u2,… un jsou náhodná čísla s normovaným normálním rozdělením, a kij jsou konstanty, které určíme postupně pomocí kovarianční matice (4a) a vztahů: c11 = k112 c12 = k11k12 c22 = k122 + k222 (15) ……………… cij = k1ik1j + k2ik2j + k3ik3j + … + kiikij ; i≤j?n. Snadno se přesvědčíme, že při n = 2 dostaneme ze (14) a (15) vztahy (13a,c). Poznámka: Konstanty kij jsou prvky trojúhelníkové matice [k] vzniklé tzv. Choleského rozkladem matice kovariancí podle vztahu [c] = [k][k]T. Poznamenejme, že výrazy (14) a (15) platí jen pro určité kombinace korelačních koeficientů. Ani v reálných úlohách se nevyskytují libovolné hodnoty. Není například možné, aby existovala vysoká pozitivní korelace mezi veličinami x1 a x2 a také mezi x1 a x3, a současně vysoká negativní korelace mezi x2 a x3. Může také existovat případ, kdy jedna vstupní veličina (např. x3) závisí na dvou jiných veličinách (x1, x2), které mohou být ještě korelovány navzájem. V takovém případě postupujeme podobně jako v kapitole 4, vztahy (10), (11). Nalezneme regresní závislost x3 = f(x1, x2) a reziduální směrodatnou odchylku sres,3 jednotlivých hodnot x3 od regresní funkce. Potom k náhodně generovaným hodnotám x1 a x2 počítáme hodnoty regresní funkce a ke každé přidáme náhodnou hodnotu u3 sres,3.
6. Nelineární závislost dvou veličin Koeficient korelace (2), vycházející z definice kovariance (1), vyjadřuje míru těsnosti lineárního vztahu. I pro nelineární závislosti se však používá koeficient korelace. V tomto případě je definován jako +√(r2), kde koeficient determinace r2 se počítá podle výrazu (8),
155
přičemž reziduální směrodatná odchylka se stanoví ze součtu čtverců odchylek jednotlivých hodnot od regresní křivky. (Tuto hodnotu r2 uvádí u nelineární regrese i Excel, který ji nazývá Hodnota spolehlivosti; viz „Graf“ – „Přidat spojnici trendu“ – „Možnosti“.) V případě nelineární korelace mezi dvěma (vstupními) veličinami xi, xj můžeme simulační výpočty Monte Carlo zpřesnit podobně jako při lineární korelaci. Nalezneme vzájemný regresní vztah xj = f(xi), koeficient determinace r2 a reziduální směrodatnou odchylku sj,res podle (9), kde sj je celková směrodatná odchylka veličiny xj. Při simulaci generujeme v prvém kroku náhodnou hodnotu té z veličin, kterou jsme při regresní analýze zvolili jako nezávisle proměnnou (xi), a ve druhém kroku vypočítáme odpovídající hodnotu druhé veličiny ležící na regresní křivce, xj = f(xi), a k ní přičteme složku usj,res, charakterizující náhodné kolísání veličiny xj kolem regresní křivky. (Uvedený postup předpokládá, že rozdělení jednotlivých hodnot kolem regresní křivky je normální a reziduální rozptyl v celém oboru xi konstantní. Podobným způsobem, který byl navržen v [8], by šlo postupovat i v jiných případech.)
7. Autokorelace Autokorelace znamená, že určitá hodnota náhodné veličiny závisí do jisté míry i na hodnotách téže veličiny v sousedních bodech nebo v předcházejících časových okamžicích. Příkladem mohou být vlastnosti betonu v konstrukci, vlastnosti základové půdy, anebo teploty konstrukce, které se v průběhu dne mění a každý den jsou jiné, ale do značné míry závisí na ročním období. Při lineární závislosti počítáme autokovarianci nebo součinitel autokorelace opět podle vztahů (1) a (2), do nichž za hodnoty xi dosadíme řadu naměřených hodnot vyšetřované veličiny, a za xj dosadíme hodnoty téže veličiny, xi, ale posunuté o určitý úsek (časový nebo geometrický). Při posunu o jednu hodnotu, tj. při xj,K = xi,K-1, dostaneme autokovarianci nebo autokorelační koeficient prvého řádu, rI. (Při generování příslušné náhodné hodnoty bychom pak užili vztahu xK = xK-1 + zK, kde zK je náhodná složka). Pro xj,K = xi,K-2 dostaneme autokorelační koeficient druhého řádu, rII, atd. Tímto způsobem můžeme získat tzv. autokorelační funkci, r = r(∆x), která vyjadřuje, do jaké míry je určitá hodnota ovlivněna hodnotami vzdálenějšími. U materiálových vlastností (např. pevnost betonu) je pravděpodobné, že v malé oblasti se budou od sebe navzájem lišit méně než od vlastností v místech vzdálených. To znamená že koeficient autokorelace bude se vzdáleností klesat. Při výpočtech se často předpokládá exponenciální charakter autokorelační funkce, např. [9, 10]: ri,j = exp [– (ξ/L)2] , (16) kde ξ je vzdálenost dvou bodů xi a xj, a L je tzv. korelační délka, která v uvedeném případě odpovídá vzdálenosti, v níž koeficient korelace r klesne na hodnotu 1/e, tj. přibližně na 37%. Tento přístup umožňuje modelovat náhodné rozložení vlastností v určité oblasti (tj. náhodná pole), a užívá se u stochastické varianty metody konečných prvků. Blíže viz např. [5, 9].
8. Diskuse Existují i další postupy pro generování korelovaných náhodných veličin, např. s použitím tzv. simulovaného žíhání [5, 11]. Některé byly implementovány do speciálních počítačových programů pro pravděpodobnostní analýzu konstrukcí [12]. Postupy uvedené v oddílech 4 6 tohoto příspěvku jsou však mimořádně jednoduché a nevyžadují zvláštní matematické znalosti ani vybavení; vztahy (13) – (15) lze snadno implementovat i do univerzálních simulačních programů jako je AntHill. Navíc je nutno mít na zřeteli, že i při větším počtu vstupních veličin vykazují výraznou korelaci často jen dvě nebo tři veličiny. Je pravda, že téměř všechny prvky korelační matice naměřených hodnot bývají nenulové. Někdy to ale může být způsobeno čistě náhodným setkáním výběrových hodnot. Je proto vždy vhodné ověřit, zda zjištěné korelační koeficienty r jsou statisticky významné (odlišné od nuly). Pro N dvojic hodnot lze předpokládat, že korelace existuje, jestliže platí | r | (N – 2)1/2 / (1 – r2)1/2 > tα,N-2 ,
(17)
156
kde α je hladina významnosti a tα,N-2 je α-kritická hodnota rozdělení t pro N-2 stupně volnosti. [Důležitou roli hraje, z kolika naměřených hodnot byl korelační koeficient počítán. Vztah (17) také poskytuje určité vodítko pro stanovení potřebného rozsahu měření k prokázání korelace.] Ale i když je nějaký nízký koeficient statisticky významný, stojí za to uvážit, zda je nezbytně nutné při simulačních výpočtech tuto korelaci uvažovat. Všimněme si příkladu z obr. 1b. Testová charakteristika (17) pro r = 0,309 a N = 200 vychází 4,58, což je o dost více než jednoprocentní kritická hodnota rozdělení t (= 2,60). Korelace je tedy statisticky významná. Přesto je z obr. 1b zřejmé, že je z praktického hlediska nevýrazná. Připomeňme, že hlavním důvodem práce s korelacemi je zpřesnění předpovědí tím, že z celkového rozptylu dané veličiny se oddělí část, která je vysvětlitelná změnami jiné veličiny. Jak vyplývá z (9), náhodná složka kolísání, charakterizovaná směrodatnou odchylkou, se při uvažování korelace zmenší z celkové hodnoty s na hodnotu sres = s(1 – r2)1/2. Například pro r = 0,3 klesne sres pouze na 95% „původní“ hodnoty s. Takovéto „zpřesnění“ je často zanedbatelné, někdy i s ohledem na přesnost, s jakou byly získány vstupní údaje. Při r = 0,9 však sres klesá na 44% celkové směrodatné odchylky, a zpřesnění může být významné. Roli hraje i samotná velikost směrodatné odchylky s vzhledem k průměru m. Podobné úvahy jsou velmi důležité také u autokorelovaných veličin, kde počet hodnot potřebný pro spolehlivé stanovení autokorelační funkce může být několikanásobně vyšší než pro pouhé určení rozptylu. Jsou-li potřebná měření nákladná (např. modul pružnosti zeminy u pozemní stavby), je vhodné nejprve matematickým modelováním pro různé korelační délky ověřit, do jaké míry mohou být výsledky konkrétního problému autokorelací ovlivněny [5]. Teprve zjištění výrazného vlivu je východiskem pro uskutečnění doplňujících měření.
Oznámení Příspěvek vznikl v rámci grantového projektu GAČR 103/01/0243.
Literatura [1] Marek, P., Guštar, M.: Ant-Hill, M-Star a další počítačové programy pro simulační analýzu konstrukcí. ARTech, Nad Vinicí 7, Praha 4, 2000. [2] Škrášek, J., Tichý, Z.: Základy aplikované matematiky III. SNTL, Praha, 1990. [3] Pelikán, J.: Modelování a simulace náhodných jevů. Sborník III. celostátní konference Spolehlivost konstrukcí, 10. 4. 2002. Dům techniky Ostrava, 2002, s. 93 – 98. [4] Florian, A., Navrátil, J., Stráský, J.: Moderní metody analýzy mostních konstrukcí. Projekt FRVŠ č. 685/94. VUT FAST Brno, zpráva 1995. [5] Rusina, R.: Metoda stochastických konečných prvků. Disertační práce, VUT Brno, 2002. [6] Virius, M.: Aplikace matematické statistiky (Metoda Monte-Carlo). ČVUT, Praha, 1991. [7] Čačko, J., Bílý, M., Bukoveczky, J.: Meranie, vyhodnocovanie a simulácia náhodných procesov. Veda, Bratislava, 1984. [8] Menčík, J.: Reliability assessment by the Monte Carlo method with small amount of data. Building Research Journal, 47 (1999), No. 1, pp. 27 – 39. [9] Teplý, B., Novák, D.: Spolehlivost stavebních konstrukcí. VUT - CERM, Brno, 1999. [10] Kala, Z.: Respecting the influence of geometrical and material imperfections of a steel frames when calculating their load-carrying capacity. Proc. Int. Conf. Quality and Reliability in Building Industry, Levoča, 24.-26. 10. 2001, TU Košice, 2002, s. 248 – 255. [11] Vořechovský, M., Novák, D., Rusina, R.: A new efficient technique for samples correlation in Latin Hypercube Sampling. Proc. 7th Int. Sci. Conf., Košice, May 22 – 24, s. 102 – 108.
157
[12] Novák, D., Rusina, R., Vořechovský, M.: FREET – software pro pravděpodobnostní posudky výpočtově náročných problémů mechaniky kontinua. Sborník III. celostátní konference Spolehlivost konstrukcí, 10. 4. 2002. DT Ostrava, 2002, s. 71 – 74.
