Costruzione di macchine

Modulo di:

Progettazione probabilistica e affidabilità

Marco Beghini

Lezione 6:

Combinazioni di variabili aleatorie

Combinazioni di più variabili aleatorie continue

• Distribuzione bivariata

• Variabili Aleatorie indipendenti

• Carico-resistenza (caso statico)

• Metodo di Montecarlo

Combinazioni di più variabili aleatorie continue (1/2)

Distribuzione bivariata o combinata di due V.A.: definizione

y

1x dx+ x

1y dy+

1x

1y

( ) ( ){ } ( )1 1 1 1 1 1, , ,P x x x dx y y y dy x y dxdyϕ∈ + ∩ ∈ + =

Combinazioni di più variabili aleatorie continue (2/2)

Distribuzione bivariata o combinata di due V.A.: proprietà

( ) ( ){ } ( )( )( ) ( )

1 1 1, , ,

,

x

x

f x dx P x x x dx y x y dy dx

f x x y dy

ϕ

ϕ

∞

−∞

∞

−∞

= ∈ + ∀ =

=

∫

∫( ) ( ),y

Distribuzioni marginali

f y x y dxϕ∞

−∞= ∫

Parametri e momenti, , ,x y x yμ μ σ σ

Covarianza e correlazione

( ) ( ) ( ) ( ), ,x yCOV X Y x y x y dxdyμ μ ϕ∞ ∞

−∞ −∞= − ⋅ − ⋅∫ ∫

( ) ( ),1 , 1

x y

COV X YX Yρ

σ σ− ≤ = ≤

Algebra delle variabili aleatorie 1/2

( ),z g x y=Date due V.A. che assumono valori x e y e una funzione:

le proprietà della V.A. che assume valori z sono:

( ) ( ) ( ), , ,z E g x y g x y x y dxdyμ ϕ∞ ∞

−∞ −∞

⎡ ⎤= = ⋅⎣ ⎦ ∫ ∫

( ) ( )22 , ,z zg x y x y dxdyσ μ ϕ∞ ∞

−∞ −∞

⎡ ⎤= − ⋅⎣ ⎦∫ ∫in particolare:

[ ] ( ) ( )E x y E x E y± = ±

[ ] 2 22x x y yVAR x y σ ρσ σ σ± = ± +

[ ] 2 2se scorrelate x yVAR x y σ σ± = +

Algebra delle variabili aleatorie 2/2

Combinazioni lineari

( ) ( )x y x yE a x a y a E x a E y⎡ ⎤± = ±⎣ ⎦

2 2 2 22x y x x x y x y y yVAR a x a y a a a aσ ρ σ σ σ⎡ ⎤± = ± +⎣ ⎦

2 2 2 2se scorrelate x y x x y yVAR a x a y a aσ σ⎡ ⎤± = +⎣ ⎦

Dati due paramenti deterministici non nulli: ,x ya a

Distribuzione combinata di due V.A. gaussiane

Due variabili aleatorie gaussiane

( ) ( )( )( ) 22

21 2

2 1

2

1,2 1

x y yx

x x y y

x y yx

x y

x y e

μ μ μμρ

σ σ σ σρϕ

πσ σ ρ

⎡ ⎤− − ⎛ ⎞−⎛ ⎞−− ⎢ ⎥− +⎜ ⎟⎜ ⎟ ⎜ ⎟⎢ ⎥− ⎝ ⎠ ⎝ ⎠⎣ ⎦=−

Esempio Mathcad sull’effetto dei parametri

xy

50; 13; 33;x y x yμ μ σ σ= = = =

Due variabili aleatorie gaussiane

y

x

0;ρ =

50; 13; 33;x y x yμ μ σ σ= = = =

Due variabili aleatorie gaussiane

x

y

0.9;ρ =

Due variabili aleatorie gaussiane indipendenti

( )

( ) ( )

22

22

12

1122

1,2

1 12 2

yx

x y

yx

yx

yx

x y

yx

x y

x y

x y e

e e

f x f y

μμσ σ

μμσσ

ϕπσ σ

πσ πσ

⎡ ⎤⎛ ⎞−⎛ ⎞−⎢ ⎥− +⎜ ⎟⎜ ⎟ ⎜ ⎟⎢ ⎥⎝ ⎠ ⎝ ⎠⎣ ⎦

⎛ ⎞−⎛ ⎞− − ⎜ ⎟− ⎜ ⎟ ⎜ ⎟⎝ ⎠ ⎝ ⎠

= =

⋅ =

= ⋅

0ρ =

Indipendenti Scorrelate⇒

Calcolo dell’affidabilità per carichi staticiEsempio 6.1Un tirante ha una resistenza definita da una V.A. S gaussiana con parametri: il carico L è anch’esso una V.A.gaussiana con parametri:determinarne l’affidabilità per una singola applicazione del carico.

30kN, 5kNS Sμ σ= =20kN, 6kNL Lμ σ= =

0 20 40 600

0.05

0.1

( )f x

( )kNx

( )Sf x

( )Lf x

Es. 6.1: I soluzione

Si considera la V.A. D differenza resistenza carico:

D S L= −

È una combinazione lineare di due V.A. gaussiane

Proprietà:La differenza di due V.A. gaussiane è una V.A. gaussiana.

Dimostrazione

D1D 1D dD+

Determinare la probabilità che la differenza D sia nell’intervallo specificato:

1S L D− =S

L

1D

( ),S Lϕ

1D dD+

( ) ( )1

11 ,

D dD L

D D Lf D dD S L dS dLϕ

∞ + +

−∞ +

⎡ ⎤= ⎢ ⎥⎣ ⎦∫ ∫

( ) ( )2 2, , 2D S L S L S Lf D N D μ μ σ σ ρσ σ= − + −

In generale:

( ) ( )2 2 2 2, , 2z x x y y x x x y x y y yf z N z a a a a a aμ μ σ ρ σ σ σ= + + +

Dati due paramenti deterministici non nulli: ,x ya a

e due V.A. gaussiane :

con coefficiente di correlazione:

( ) ( )( ) ( )

, ,

, ,x x x

y y y

f x N x

f y N y

μ σ

μ σ

=

=ρ

La V.A.:

Ha la seguente densità di probabilità:

x yz a x a y= +

D S L= −

è la combinazione lineare di due V.A. gaussiane indipendenti

( ) ( ) ( )2 2,10, 5 6 ,10,7.81Df D N D N D= + =

20 0 20 400

0.02

0.04

0.06

0.90R =

D

( )Df D

Es. 6.1: II soluzione

S

L

( ),S Lϕ

L S=

L S>

( ) ( ),L

R P S L L S dS dLϕ∞ ∞

−∞

⎡ ⎤= > = ⎢ ⎥⎣ ⎦∫ ∫

S L>

Es. 6.1: III soluzione (1/2)

( )Lf x ( )Sf x

x1x 1x dx+

per il carico nell’intervallo

Il contributo alla probabilità di non rottura (affidabilità) è

( )1 1,x x dx+

( ) ( ) ( ) ( )

( ) ( )1

1 1 1

1 1

L L Sx

L S

dR f x dx P S x f x dx f x dx

f x dx R x

∞= ⋅ > = ⋅ =

= ⋅

∫

( )1SR x

III soluzione (2/2)

A questo punto è necessario considerare tutte le possibilità per illivello di carico ed effettuarne la somma:

( ) ( )L SR dR f x R x dx∞

−∞= = ⋅∫ ∫

Integrando con inversione di variabili:

( ) ( )S LR f x F x dx∞

−∞= ⋅∫

Definizioni

Safety Margin D

D

SM μσ

= è il reciproco del CV per la differenza S − L

2 2S L

S L

SM μ μ

σ σ

−=

+

D

( )Df D ( )DF D

DDμDSM σ⋅

( )0rottura DP F= ( )1 0DR F= −

Definizioni

Loading Roughness2 2

L L

D L S

LR σ σσ σ σ

= =+

f

x

f

x

Lf

SfLf

Sf

0LR ≅ 1LR ≅

rotturaProtturaP

Probabilità di rottura in funzione del SM

Hyp. S e L variabili aleatorie gaussiane

0 1 2 3 4 5 6 7 8

10-2

0

10-4

10-6

10-8

10-10

10-12

10-14

10-16

SM

.rottP

Coefficiente di sicurezza statico?

Esempio 6.2Dati i seguenti valori di V.A. gaussiane (valori in MPa):

determinare il coefficiente di sicurezza convenzionale, il SM e la probabilità di rottura

( ) ( ),120,42 ; , 270,30L N x S N x= =

0 100 200 300 4000

0.005

0.01

0.015

x

( )Lf x( )Sf x

206MPamaxL = 221MPaamS =

Carico max: 98 percentile alto

Resistenza ammissibile 5 percentile basso

1.07am

max

SL

η = =

2.91; 0.81SM LR= = 31.83 10 0.183%rotturaP −= ⋅ =

Esempio 6.2 bisDati i seguenti valori di V.A. gaussiane (valori in MPa):

determinare il coefficiente di sicurezza convenzionale, il SM e la probabilità di rottura

( ) ( ),500,40 ; ,672,30L N x S N x= =

623 1.07582

am

max

SL

η = = =

3.44; 0.80SM LR= = 40.291 10 0.029%rotturaP −= ⋅ =

Valori tipici di CV per resistenze meccaniche (leghe metalliche)

Resistenza CV

Ultimate stress 0.05

Yielding stress 0.07

Fatigue limit 0.10

Hardness (Brinnel) 0.05

Fracture strength (lower shelf) 0.08

Interpretazione dell’affidabilità statica (Trasformazione di Mellin)

( ) ( )L SR f u R u du∞

−∞= ⋅∫

( )LX F u=

( )SY R u=( ) ( )L Lf u du dF u dX= =

R Y dX= ⋅∫SY R=

LX F=1

1

R

RP

Metodo Montecarlo

• Strumento per effettuare simulazioni numeriche• Algoritmo per calcolare integrali• Permette di simulare in forma numerica situazioni complesse

Esempio Valutare la probabilità che una palla da biliardo lanciata a caso si fermi all’interno di una figura disegnata sul tavolo

• Lanciare un palla a caso significa che la posizione finale è definita da una distribuzione bivariata uniforme nel rettangolo

• Le coordinate della posizione finale della palla sono V.A.indipendenti uniformemente distribuite

• Un generatore di numeri casuali tra 0 e 1 definisce una coordinata

• L’insieme di punti che si ottengono per ripetizione della procedura si chiama campo di Poisson

• Non serve avere l’espressione analitica del contorno della figura, basta una funzione che stabilisca se il punto è dentro o fuori

• Il risultato si ottiene applicando la definizione statistica di probabilità (si prendono tantissimi punti a caso)

• Si comprende l’uso del metodo per il calcolo degli integrali

Generatori di numeri (pseudo)casuali

( )0 1?; mod ; /i i i ix x c a x m r x m−= = + ⋅ =

0m > modulus

0 a m< < multiplier

0 c m≤ < increment

• c e m primi tra loro, • a−1 è divisibile per tutti i fattori primi di m, • a−1 è multiplo di 4 se m è multiplo 4

Esempim a c

• Numerical Recipes 232 1664525 1013904223• Borland C/C++ 232 22695477 1

Campi di Poisson

50n = 200n =

1000n =

Esempio 6.3Stimare π con il metodo Montecarlo

4P π=

4Pπ =

Probabilità che il punto stia dentro lo spicchio di cerchio:

Si genera una distribuzione di Poisson di n punti e si calcola la frazione degli nint punti interni :

4 intnn

π ≅

0 200 400 600 800 10000

1

2

3

4

n

4 intnn

Soluzione

Generazione di numeri casuali con distribuzione voluta (1/2)

( )xdFdy dx f x dxdx

= ⋅ = ⋅

Numero di punti nell’intervallo: dyδ

Densità di corrispondenti valori di x: ( )xdy f xdx

δ =

Se la densità di y è costante e pari a: ( ) 1/1 1yf yδ = = =

dy

dxx

( )F x1

y

Generazione di numeri casuali con distribuzione voluta (2/2)

( )1x F y−=

Applicando a una V.A. distribuita uniformemente in (0,1) la funzione cumulata inversa di una distribuzione nota, si genera una V.A. con tale distribuzione

Esempio 6.3Generare n=20 valori distribuiti secondo una Weibull con parametri:

( ) 0

0

0

0

0 , , ,

1 x x

x xF x x

e x xβ

αα β −⎛ ⎞−⎜ ⎟

⎝ ⎠

≤⎧⎪= ⎨⎪ − >⎩

01, 2, 3xα β= = =

( )

0

0

0

1/0

1/

0

1

ln 1

1ln1

1ln1

x x

F e x x

x x F

x xF

x xF

β

α

β

β

β

α

α

α

−⎛ ⎞−⎜ ⎟⎝ ⎠= − >

−⎛ ⎞− = −⎜ ⎟⎝ ⎠

− ⎡ ⎤⎛ ⎞= ⎜ ⎟⎢ ⎥−⎝ ⎠⎣ ⎦

⎡ ⎤⎛ ⎞= + ⎜ ⎟⎢ ⎥−⎝ ⎠⎣ ⎦

1/

01ln

1x x

y

β

α⎡ ⎤⎛ ⎞

= + ⎢ ⎥⎜ ⎟−⎝ ⎠⎣ ⎦

yk xk0.394 4.001

0.586 4.763

0.480 4.310

0.102 3.215

0.425 4.107

0.692 5.355

0.387 3.977

0.780 6.025

0.420 4.088

0.155 3.338

0.461 4.236

0.547 4.583

0.310 3.743

0.073 3.152

0.040 3.081

0.617 4.921

0.797 6.188

0.057 3.117

0.735 5.655

0.170 3.373

xsk3.081

3.117

3.152

3.215

3.338

3.373

3.743

3.977

4.001

4.088

4.107

4.236

4.31

4.583

4.763

4.921

5.355

5.655

6.025

6.188

k1

2

3

4

5

6

7

8

9

10

11

12

13

14

14

16

17

18

19

20

Fsk0.034

0.083

0.132

0.181

0.230

0.279

0.328

0.377

0.426

0.475

0.525

0.574

0.623

0.672

0.721

0.770

0.819

0.868

0.917

0.966

0.5k

kFsn−

=

3 3.5 4 4.5 5 5.50

0.2

0.4

0.6

0.8

1

Fs

F

F

x

Fs

1ln ln1 F

⎛ ⎞⎛ ⎞⎜ ⎟⎜ ⎟−⎝ ⎠⎝ ⎠

F

( )0ln x x−2 1.5 1 0.5 0 0.5 1

4

2

0

2

Esempio 6.4La sollecitazione di un elemento meccanico è definita da una V.A. di Weibull con parametri:

e la resistenza è una V.A. di Weibull con parametri:

Determinare l’affidabilità per un singolo caricamento statico.

050, 2.0, 100L L Lxα β= = =

0100, 3.0, 120S S Sxα β= = =

( ) ( )0L

L SxR f u R u du

∞= ⋅∫

( )0

min(1, )S

S

S

u x

SR u e

β

α⎛ ⎞−

−⎜ ⎟⎝ ⎠=

( ) ( )01

00

LL L

L

u xLL

L LL L

u xf u e u x

ββαβ

α α

− ⎛ ⎞−−⎜ ⎟⎝ ⎠⎛ ⎞−

= ⋅ ⋅ ≥⎜ ⎟⎝ ⎠

0.948thR =

L

S

0 100 200 3000

0.005

0.01

0.015

0.02

u

( )f u

Distribuzioni

Soluzione con Montecarlo

1 1000 20000.9

0.95

1

R

n

Soluzione numerica in funzione della numerosità del campione di valori del metodo Montecarlo

Precisione e convergenza con Montecarlo

Indichiamo con n il numero di tentativi della simulazione e rispet.:

la probabilità teorica da calcolare e quella ottenuta con Montecarlo

, nP P

È stato dimostrato che:

( )1Prob 2 0.95n

P PP P

n

⎛ ⎞−⎜ ⎟− < =⎜ ⎟⎝ ⎠

Oppure, equivalentemente, che con la confidenza del 95% l’incertezza relativa sulla probabilità (semiampiezza dell’intervallo relativo di previsione) è:

1 12 1nP Pe

P P n−

= = −

Nel caso esaminato: 0.948P =

10 100 1000 10000.

0.01

0.1

0.001

e

n

Stima della semiampiezza dell’intervallo relativo di previsione

La formula mostra che l’incertezza relativa è più alta per la stima di probabilità basse (necessità di n elevati)

2000 4000 6000 8000 1 .1040

0.05

0.1

e

n

Stima della semiampiezza dell’intervallo relativo di previsione:scala lineare

0 1000 20000.9

0.95

1

R

n

Intervalli di confidenza del 95% per il caso esaminato