...

Affidabilità dei sistemi - Università degli Studi della Basilicata

by user

on
Category: Documents
33

views

Report

Comments

Transcript

Affidabilità dei sistemi - Università degli Studi della Basilicata
AFFIDABILITA’ DEI SISTEMI STOCASTICI
(semplici)
Un sistema (o uno qualsiasi dei
suoi componenti) può essere soggetto a stress “casuali”.
Es: un fusibile in un circuito;
una trave di acciaio sotto carico;
l’ala di un aereo sotto l’influenza di forze
Collasso del
componente
Lo stato in cui il componente
(o il sistema) cessa di funzionare
Definizione
Se un componente è messo in condizioni di “stress” in un istante temporale
t=0 e viene osservato fino a quando “collassa”, la durata della vita (o il tempo
necessario al collasso) è una variabile aleatoria che indicheremo con T
Definizione di guasto
Per guasto si intende
(a) la cessazione di un dispositivo ad adempiere
la funzione richiesta (guasti totali);
(b) variazione prestazione del dispositivo che lo
renda inservibile (guasti parziali);
Guasti intermittenti
Classificazione dei guasti in base al tempo di vita del sistema
Guasti infantili
Guasti di usura
Guasti casuali
errori montaggio
errori progettazione
Indipendenti
dal tempo
E’ anche nota come v.a. esponenziale negativa. Perché?
Esercizio: calcolare i momenti della v.a. esponenziale
Esercizio: Supponiamo di aver messo in funzione un pezzo alle ore
8 di un certo giorno. Supponiamo che resti costantemente in funzione. Supponiamo che la durata della sua vita sia espressa da una
legge esponenziale (in ore). Calcolare la probabilità che il pezzo cessi
di funzionare durante un turno di lavoro 8-16.
Nelle applicazioni, la v.a. esponenziale viene usata per descrivere il tempo di vita
(aleatorio) di un sistema complesso fino al suo collasso.
FX (t ) ⇒ funzione di guasto (failure function)
RX (t ) = 1 − FX (t ) ⇒ funzione di affidabilità (reliability function)
Prop: Se X ∼ Exp(λ ), allora P ( X > s + t | X > s ) = P ( X > t )
ASSENZA DI MEMORIA
Si dimostra che
(a) Se P ( X > s + t | X > s ) = P ( X > t ) ∀s, t ≥ 0 ⇒ X ∼ Exp (λ )
(b) Se P ( X > s + t | X > s ) = P ( X > t ) ∀s, t ∈ {0,1, 2,…} 0 ⇒ X ∼ Geom( p )
Altra formulazione della proprietà di assenza di memoria
La durata di vita residua ha la stessa distribuzione della durata
di vita del componente considerato.
Prop: Sia X ∼ Exp (λ ) e sia a >0. Posto X a = X − a allora
P ( Xa ≤ t | X > a) = P ( X ≤ t )
Se T è un tempo di vita con legge qualsiasi, qual è la distribuzione della sua vita
residua?
FT ( a + t ) − FT ( a )
FTa ( t ) = P (Ta < t ) = P (T − a ≤ t | T ≥ a ) =
1 − FT ( a )
Nel caso esponenziale, questa proprietà lascia qualche perplessità,
perché implica che il sistema non è soggetto ad usura.
Def : Sia T il tempo di vita di un sistema, e sia T
v.a. assolutamente contina. Si definisce tasso di
gausto (o funzione di rischio ) - failure rate or
hazard function la funzione
f (t )
Z (t ) =
R(t )
dove f (t ) rappresenta la funzione densità di T
e R(t ) la sua funzione di affidabilità.
EX: Nel caso esponenziale, il tasso di guasto è….
Tasso di guasto in funzione dell’età
Significato probabilistico :
t + ∆t
f (s )ds
f (t )∆t
∫
t
Z (t )∆t =
=
= P(t ≤ T < t + ∆t | T ≥ t )
P(T ≥ t )
P(T ≥ t )
probabilità che il componente che è durato fino all'istante t collassi nell'intervallo di tempo
successivo [t , t + ∆t ]
Teorema Se T è una variabile aleatoria continua, con funzione di
distribuzione F ( ⋅) tale che F ( 0 ) = 0, allora
 t

R ( t ) = exp  − ∫ Z ( s)ds  , t ≥ 0
 0

f (t )
R ' (t )
⇒ Z (t ) = −
⇒
DIM: Z (t ) =
R(t )
R(t )
t
∫ Z(s)ds = − log R ( t ) + cos t.
⇒
0
(
t
)
R ( t ) = K exp − ∫ Z ( s )ds dove K dipende dalle condizioni iniziali
0
Esercizio: Stabilire se la funzione R(t)=exp[-(t/b)^k] è una funzione di affidabilità.
Determinare il tasso di guasto.
Proprietà del tasso di guasto
i Z (t ) ≥ 0
i Z (t ) ≥ f (t )
i
∫
t
0
i
Z ( s ) ds < ∞
∫ Z (t ) d t =∞
Dunque non è una probabilità!
ℜ
Contiene informazioni istantanee!
t
H (t ) = − ln[ R(t )] = ∫ Z ( s ) ds ⇒ tasso di guasto cumulativo
0
Proprietà del tasso di guasto cumulativo
i H (t ) ≥ 0
i H ( t ) è funzione crescente nel tempo.
Cosa accade nel caso esponenziale?
Perché si introduce ?
Un problema che spesso va affrontato nelle applicazioni è stabilire quale
sistema è più affidabile.
Def : Un sistema A si dice più affidabile di un sistema B
se RA ( t ) ≥ RB ( t ) per ogni t.
 > 1 se H A ( t ) > H B (t )
RA ( t )

= exp  H A ( t ) − H B (t )  =  1 se H A ( t ) = H B (t )
RB ( t )
< 1 se H ( t ) < H (t )
A
B

Teorema : Se Z A ( t ) > Z B ( t ) ∀t ⇒ RA ( t ) <RB ( t ) ∀t
(vale anche sul segno <)
Non vale il viceversa.
100
Z(t)=10t
Z(t)=(t-3)2
90
80
70
60
50
40
30
20
10
0
0
1
2
3
4
5
6
7
8
9
10
1
R(t) associata a 10t
R(t) associata a (t-3)2
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
Esercizio: Si assuma che un dispositivo abbia un tasso di guasto lineare. Determinare
la funzione di affidabilità e la funzione di guasto.
Esercizio: Si assuma che un dispositivo abbia una funzione di guasto descritta da una
v.a. uniforme sull’intervallo [0,100]. Determinare il tasso di guasto.
Come si descrivono gli altri rami della curva a vasca da bagno?
costante
β = 1

Z (t ) = αβ t β −1 =  β > 1 crescente
 β < 1 decrescente

Esercizio: Il modello di vita di un apparato propulsivo di un sottomarino atomico ha
Un tasso di guasto pari a 2 t^3. Determinare l’affidabilità del sistema a due mesi.
DISTRIBUZIONE DI WEIBULL
Calcolare la funzione densità.
Una v.a. di Weibull ha funzione densità
(αβ ) t β −1 exp ( −α t β ) t > 0
f (t ) = 
0
altrove

con α e β parametri reali nonnegativi.
NB: quella sul MATLAB ha la stessa forma. DISTTOOL:
Esercizio: Calcolare la funzione di ripartizione
Esercizio: Calcolare i momenti centrali
Per β = 2, si ha Z (t ) = 2α t ⇒ Distribuzione di Raylegh
La distribuzione di Rayleigh (Lord Rayleigh - fisico britannico premio nobel),
di parametro α =σ 2 >0, descrive la v.a. Z = X 2 + Y 2 dove X ∼ N (0, σ 2 ) e
Y ∼ N (0, σ 2 ) (indipendenti).
Visualizzare funzione in Matlab
Per β = 1/ 2, si ha Z (t ) = 0.5α t −0.5 ⇒ Curva mortalità infantile
NOTA: I parametri della legge di Weibull
• Il parametro β è detto parametro di forma. Questo perchè variando
questo parametro, varia la forma della pdf
Weibull pdf
5
beta=0.5
beta=1
beta=1.5
beta=5
4.5
4
3.5
3
2.5
2
1.5
1
0.5
0
0
0.5
1
1.5
2
2.5
• Il parametro α è detto parametro di scala. Questo perchè variando
questo parametro, varia il range in cui la pdf è diversa da zero. Però
per vederlo in questa forma è necessario riscrivere la pdf.
Anzicchè
(αβ ) t β −1 exp ( −α t β ) t > 0
f (t ) = 
0
altrove

 1 
rimpiazziamo α con  β  e scriviamo
α 
 β  t  β −1
  t β 
   exp  −    t > 0
f (t ) =  α  α 
 α  

0
altrove

Qual è il vantaggio?
Esercizio: Il tempo di guasto di un radar per la navigazione aerea è
rappresentabile mediante una distribuzione di Weibull con parametro di
scala pari a 2000 ore e parametro di forma pari a 2.1.
Questo significa che il
tempo di vita del sistema evolve su un intervallo temporale piuttosto ampio. Si ha
2.1
α =1/2000 ≈ 10 ^ (−7)
Weibull, scala=2000, alfa=10(-7)
-4
4.5
x 10
4
3.5
3
2.5
2
1.5
1
0.5
Attenzione all’uso del Matlab
0
0
500
1000
1500
2000
2500
3000
3500
4000
4500
5000
Anzicchè
(αβ ) t β −1 exp ( −α t β ) t > 0
f (t ) = 
0
altrove

 1 
rimpiazziamo α con  β  e scriviamo
α 
β

f (t ) =  α


t 
 
α 
β −1
  t β 
exp  −    t > 0
 α  


0
altrove
MATLAB
CON PARAMETRO
DI SCALA
Esercizio: Sia T il tempo di vita di un componente elettronico in un aereo descritto da una v.a. di Weibull di parametro di scala 1100 ore di volo e di parametro di forma 3. Determinare la probabilità di guasto dopo 100 ore.
a ) Trovare il parametro α
b) Calcolare weibcf
Esercizio: Determinare la massima lunghezza di volo tale la probabilità di
guasto sia inferiore a 0.05.
Prob. di guasto
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0
500
1000
1500
2000
2500
3000
Altra legge impiegata per descrivere il tempo di vita di un sistema è la legge
gaussiana per via del teorema del limite centrale.
Prestare attenzione alla legge dei 3 sigma
Talvolta viene
usato il
valore assoluto
della
gaussiana!
Confronti
0.8
Gauss. standard e val. assoluto della gaussiana a confronto
Gauss. Stand.
|Gauss. stand.|
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
-3
-2
-1
0
1
2
3
Tasso di guasto di un modello gaussiano
T ∼ N (100,10)
tasso di guasto di una gaussiana
0.45
0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05
0
60
70
80
90
100
110
120
130
140
Esercizio: Due compagnie producono pneumatici per auto.
Entrambe le compagnie dichiarano che il tempo di vita medio dei
loro prodotti è 2000 miglia. Dopo una indagine statistica risulta che
il tempo di vita dei pneumatici prodotti dalla compagnia A è esponenziale di parametro 0.0005 e invece il tempo di vita dei pneumatici prodotti dalla compagnia B è gaussiano, di media 2000 miglia
e di deviazione standard 200 miglia. Se la politica di un negozio
che affitta le macchine è di cambiare i pneumatici non appena
essi raggiungono le 2000 miglia, quali pneumatici converrebbe
comprare?
RA (2000) = 0.3678, RB (2000) = 0.5
Esercizio: Sia T il tempo di vita di un componente elettronico, una v.a. gaussiana
di deviazione standard 10 ore. Assumendo che dopo un periodo pari a 100 ore,
l’affidabilità stimata per il componente sia 99%, determinare la media di T.
Esercizio: il tempo di vita di un sistema ha un tasso di guasto costante pari a
0.01. Determinare per quante ore deve funzionare affinché abbia un’affidabilità
pari almeno a 90%.
Altri tipi di curve bathtub
Componenti meccanici
Curve reali
Affidabilità condizionata
Il dispositivo in esame ha funzionato fino al tempo s. Qual è la probabilità che
funzioni per un altro intervallo di tempo t?
P (T > s + t | T > s ) = R( s + t | s )
Conoscendo la funzione di affidabilità, come è possibile calcolare la funzione di
affidabilità condizionata?
R( s + t )
R( s + t | s) =
R( s)
►Cosa accade nel caso esponenziale?
► Quanto vale R(t+0|0)?
► Quali sono le proprietà di questa funzione per s fissato e t variabile?
► Quali sono le proprietà di questa funzione per t fissato e s variabile?
Tasso di guasto Weibull di par. 2 e 4
2000
Affidabilità condizionata
per d. Weibull
1800
1600
1400
1200
1000
800
600
1
0.14
400
0.9
0.12
0.8
200
0
0
1
2
3
4
5
6
7
8
9
10
0.7
0.1
0.6
0.08
0.5
Aff. Weibull di parametri 2 e 4
1
0.4
0.06
0.9
0.3
0.04
0.2
0.8
0.7
0.1
0.02
0.6
0
00
0
0.5
0.4
0.3
0.2
0.1
0
Aff.
condizionata
t+1|1
Aff.
condizionata
t+1|t
0
1
2
3
4
5
6
7
8
9
10
1
1
2
2
3
4
3
5
4
6
5
7
6
8
7
9
8
10
9
10
Qual è la probabilità condizionata di guasto di un componente di età t durante l’
intervallo seguente di durata ∆ t ?
F (t + ∆t | t ) =
R (t ) − R(t + ∆t )
R (t )
1 − R ( t + ∆t | t ) =
F (t + ∆t ) − F (t ) P (t < T ≤ t + ∆t )
=
=
= P (t < T ≤ t + ∆t | T > t )
R (t )
P(T > t )
Tempo medio di guasto (MEAN TIME TO FAILURE)
∞
Teorema: Se E [T ] < ∞ ⇒ E[T ] = ∫ R (t ) dt
0
Teorema molto generale


E [ X ] = ∫ xf ( x) dx = ∫  ∫ dt 
0
00

∞
∞
x
∞

f ( x) dx = ∫ dt  ∫ f ( x) dx 
0
t

∞
Attenzione al significato
MTTF2 = MTTF1
Esercizio: Il tempo di guasto di un RWR (radar warning receiver) in un aereo
da combattimento segue una legge di Weibull con parametro di scala 1200 ore
di volo e parametro di forma 3. Lo stesso radar montato su di un elicottero ha
un tempo di vita esponenziale con parametro 0.001. Confrontare le affidabilità
dei due sistemi.
Confrontare il tempo
medio di vita.
1
elicottero
areo guerra
0.9
Se il fornitore dà una garanzia
di 750 ore di volo, calcolare il
rischio per i due sistemi.
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0
500
1000
1500
2000
2500
3000
Tempo medio di vita residua (MEAN RESIDUAL LIFE)
Il tempo medio di guasto per un sistema di età t0 prende il nome di
tempo medio di vita residua
∞
MTTF ( t0 ) = ∫ ( t − t0 ) f (t | t0 )dt
t0
∞
Una formula diversa
per il calcolo del tempo
medio di vita residua:
∫ R(t ) dt
MTTF ( t0 ) =
t0
R ( t0 )
Esercizio: Il tempo di guasto di un radar per la navigazione aerea è rappresentabile
mediante una distribuzione di Weibull con parametro scala pari a 2000 ore e parametro
forma pari a 2.1. Uno di questi radar viene montato su di un aereo. Viene dichiarato che
l’età del radar è di 800 ore. Trovare il tempo medio di vita residua per questo radar.
MBTF= MEAN OPERATING TIME BETWEEN FAILURES
e non MEAN TIME BETWEEN FAILURES
E’ il tempo medio tra due consecutivi guasti
• Se dopo ogni riparazione, il sistema torna "come nuovo" allora MTBF=MTTF,
ignorando il tempo impiegato per la riparazione
• Se il tasso di guasto è esponenziale, allora MBTF=T / n, ignorando il tempo impiegato per la riparazione.
• Più in generale: MBTF=MTTF+MTTR
Tempo medio per la
riparazione
Se si ha un campione casuale di tempi, come si stimano tutte queste funzioni?
Esempio: si misura il tempo di vita di 262 lampadine. Il numero di guasti è riportato
in tabella. Stimare la funzione guasto, la funzione di affidabilità, la densità di guasto,
Il tempo medio tra guasti.
Costruzione dell’istogramma
a) Mettere in un vettore i centri delle classi.
b) Riportare in un altro vettore le frequenze assolute delle classi
c) Costruire l’istogramma (file prova.m)
Per stimare la densità di guasto…
bar(c,fa/taglia)
oppure costruire il poligono di frequenza plot(c,fa/taglia,’--rc’)
Per stimare la funzione di ripartizione di guasto…
cdf=cumsum(fa/taglia)
inf=[0,100,200,300,400,500,600,700,800,900];
stairs(inf,cdf) (file prova1.m)
Per stimare la funzione di affidabilità …
aff=1-cumsum(fa/taglia);
stairs(inf,aff)
Per stimare il tasso di guasto…
Z (t )∆t ≈
P ( t < T ≤ t + ∆t )
P(T > t )
>> for i=1:9
tasso(i)=(cdf(i+1)-cdf(i))/aff(i);
end
F (t + ∆t ) − F (t )
=
R(t )
tasso di guasto
1
0.9
0.8
E per stimare MTTF?
Procedure per dati raggruppati
datigrouped(c,fa)
0.7
0.6
0.5
0.4
0.3
E’ possibile stabilire il modello
di vita delle lampadine?
0.2
0.1
0
TEST
50
150
250
350
450
550
650
750
850
TEST PER LA BONTA’ DI ADATTAMENTO
TEST CHI QUADRATO
TEST DI KOLMOGOROV
Siano Oi le frequenze osservate e siano Ei le frequenze attese (cioè quelle
che ci si attende se il modello probabilistico segue quello ipotizzato).
H 0 : non esiste una differenza tra le frequenze osservate e quelle attese
H1 : esiste una differenza tra le frequenze osservate e quelle attese
2
STASTICA TEST
Oi − Ei )
(
2
χ =∑
≈ χ k2− p −1
Ei
Non esiste una procedura in MATLAB
Bisogna calcolare e!!
Chitest(o,e,alfa,grlibertà)
Ad esempio se ipotizziamo una legge gaussiana, si ha
Ei = NP ( X ∈ Ci )
>>classi=c+50;
>> e(1)=taglia*normcdf(classi(1),mean,dev);
>> for i=1:7
e(i+1)=taglia*(normcdf(classi(i+1),mean,dev)-normcdf(classi(i),mean,dev));
end
>> e(10)=taglia*(1-normcdf(classi(9),mean,dev));
osservate vs attese
Chitest(fa,e,7)
70
60
Si rigetta l’ipotesi di gaussianità
50
40
Altro tipo di test è il test di
Kolmogorv-Smirnov.
30
20
In quel caso i dati servono
singolarmente.
10
0
1
2
Si potrebbe pensare di usare una distribuzione di Weibull.
Come si stimano i parametri di scala e di forma?
Servono i dati singolarmente!
Esistono vari metodi. Quello MATEMATICAMENTE più “corretto” è il metodo
della massima verosimiglianza, che ora andremo a studiare.
Definizione: Sia ( X 1 , X 2 ,… , X n ) un campione casuale di v.a. i.i.d con legge
fϑ ( x). Indichiamo con ( x1 , x2 ,… , xn ) il valore dell'osservazione. Chiamiamo
funzione di verosimiglianza la funzione
n
Lϑ ( x1 , x2 ,… , xn ) = fϑ ( x1 , x2 ,… , xn ) = ∏ fϑ ( xi )
i =1
Definizione: Chiamiamo stimatore di massima verosimiglianza per ϑ il valore ϑˆ che
massimizza la funzione di verosimiglianza Lϑ ( x1 , x2 ,… , xn )
Calcolare il massimo della funzione di verosimiglianza può essere difficile, ma
ricordando che il log è una funzione monotona, il valore che massimizzerà
n
Lϑ ( x1 , x2 ,… , xn ) renderà massimo log  Lϑ ( x1 , x2 ,… , xn )  = ∑ log  fϑ ( xi ) .
i =1
Gli stimatori di massima verosimiglianza non sono necessariamente
corretti, ma per ipotesi abbastanza blande sulla continuità delle leggi pdf sono consistenti ed asintoticamente normali.
Molti stimatori che conosciamo sono stimatori di massima verosimiglianza.
Esempio: Supponiamo di avere i seguenti
tempi di vita
17,21,33,37,39,42,56,98,129,132,140
In MATLAB weibfit(x)
>> weibfit(x)
[alfa,beta] = 0.0011
1.5651
Per determinare quale distribuzione “fitta” meglio i dati, si possono utilizzare metodi grafici che consistono nel costruire i cosiddetti “probability plotting papers”.
Si basano sulla funzione di ripartizione cumulativa che si vuole riconoscere nei dati e
usano assi tali che, quando questi dati effettivamente provengono da una popolazione
con legge simile alla funzione di ripartizione, si distribuiscono lungo una linea retta.
(quindi la tecnica usata è quella della regressione)
Probability plotting papers necessitano di due cose:
a) I dati (ordinati in ordine crescente)
b) La funzione di ripartizione cumulativa
Metodi di stima
Costruzione F(x)
Mean Rank
i/(n+1)
Median Rank
(i-0.3)/(n+0.4) (Bernard)
Symmetric CDF
(i-0.5)/n
Particolarmente
utile per leggi
asimmetriche
Molto usata!
• Studi recenti hanno dimostrato che questo metodo grafico è più accurato
• Non è ancora chiaro quali di queste stime sia la “migliore”
Weibull probability paper


1 

1 
 y = log log

Si ha la seguente retta log log
=
β
log
x
−
β
log
α
.
Poni
1
−
F
(
x
)




1
−
F
(
x
)



z = log x

Un grafico sugli assi ( z , y ) è un Weibull probability paper
Esempio: Supponiamo di avere i seguenti tempi di vita
17,21,33,37,39,42,56,98,129,132,140
Ordinati
Dati
Mean
Rank
Median
Rank
1
17
0.0833
0.0614
2
21
0.1667
0.1491
3
33
0.2500
0.2368
4
37
0.3333
0.3246
5
39
0.4167
0.4123
6
42
0.5000
0.5000
7
56
0.5833
0.5877
mean rank
>> i=[1:1:11];
>> cdfe=(i-0.3)./11.4;
>> y=log(log(1./(1cdfe)));
1.5
1
0.5
0
-0.5
-1
-1.5
median rank
-2
1.5
-2.5
1
-3
0.5
-3.5
2.5
3
3.5
4
4.5 0
5
-0.5
-1
-1.5
RETTA DI REGRESSIONE
LINEARE!!
-2
-2.5
-3
2.5
3
3.5
4
4.5
5
β =b
scale = exp(−a / β )
Vanno trasformati!
>> beta =
1.4354 -6.2267
Confrontiamo questi dati con quelli ottenuti prima usando weibfit
β = 1.4354 (regressione)
β = 1.5651( verosimiglianza )
log scale = -6.2267 (regressione) ⇒ scale
= exp ( -beta(2) / beta(1) ) ⇒ α = 1/ scale ^ β
= 0.0020
α = 0.0011 ( verosimiglianza )
Verificare se i residui hanno legge gaussiana: sia con metodi grafici che con metodi
Inferenziali!!
>> residuals
Normal Probability Plot
0.95
residuals =
0.75
Probability
-0.5989
0.0332
-0.1005
0.1080
0.3359
0.4950
0.3276
-0.2366
-0.3843
-0.1388
0.1595
0.90
0.50
0.25
0.10
0.05
-0.6
-0.4
-0.2
0
0.2
Data
>>normplot(residuals)
TEST
0.4
Costruione norm-plot dei residui
>> ressort=sort(residuals)
>> plot(ressort,z,’*’)
costruzione normplot
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
Ci sono dei software che costruiscono i norm-plot con i quantili
della gaussiana. La forma del diagramma di dispersione non
cambia. Cambia ovviamente l’intervallo di rappresentazione
sull’asse delle y.
>> z=norminv(cdfe,media,dev);
costruzione normplot con i quantili gaussiani
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
weibplot(x)
C’è un modo per capire prima se ha senso cercare una legge
di Weibull che fitta i dati?
Weibull Probability Plot
0.96
0.90
0.75
Probability
0.50
0.25
0.10
0.05
2
10
Data
Effettuiamo un test per verificare se i dati hanno legge di Weibull
con i parametri trovati
TEST DI KOLMOGOROV-SMIRNOV
Se l’ipotesi da testare è che il campione casuale proviene da
una popolazione gaussiana con media 0 e varianza 1, allora la
sintassi è semplicemente
Se l’ipotesi da testare è che il campione casuale proviene da una
popolazione gaussiana non standard, bisogna prima standardizzare i dati e poi usare il kstest:
E’ il caso del vettore residuals su cui abbiamo lavorato prima!!
Standardizziamo:
>> resstd=(residuals-media)/dev;
>> [h,p,kstat,cv]=kstest(resstd)
h= 0
p = 0.9983
kstat = 0.1116
cv = 0.3912
CODA
Null Hypothesis: F(x) equal to CDF for all x.
For TAIL = 0 (2-sided test), alternative: F(x) not equal to CDF.
For TAIL = 1 (1-sided test), alternative: F(x) greater than CDF.
For TAIL = -1 (1-sided test), alternative: F(x) less than CDF.
For TAIL = 0, 1, and -1, the K-S test statistics are T = max|S(x) - CDF|,
T = max[S(x) - CDF], and T = max[CDF - S(x)], respectively.
TORNIAMO AL NOSTRO ESEMPIO.
17,21,33,37,39,42,56,98,129,132,140
Mettiamo i dati in un vettore x
Costruiamo il vettore della funzione di ripartizione teorica
>>alfa=1/(exp(-beta(2)/beta(1)))^beta(1)
>>ccf(:,1)=x;
>>ccf(:,2)=weibcdf(x,alfa,beta(1))
>> [h,p,ksstat,cv]=kstest(x',ccf,0.05,0)
h = 0; p = 0. 7163; ksstat = 0. 2009; cv = 0.3912;
>> weibstat(alfa,beta(1))
ans = 69.50
ATTENZIONE ALL’
USO DEGLI INTERVALLI
DI CONFIDENZA!!
Intervalli di confidenza per i parametri di una distribuzione di Weibull
Usate weibfit!!
>> [par,ci]=weibfit(x)
par = 0.0011
1.5651
ci = -0.0067 -0.0055
0.0090 3.1357
A cosa serve???
Intervalli di confidenza per la affidabilità di una legge di Weibull
1
Stimato
inf
sup
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0
20
40
60
80
100
120
140
160
180
200
Esercizio: anche con i parametri calcolati con il metodo della massima
verosimiglianza il test non rigetta l’ipotesi?.
Exponential probability paper

 1 
 1 
 y = log 

Si ha la seguente retta log 
=
λ
x
.
Poni
1
−
F
(
x
)




1
−
F
(
x
)



z=x

Un grafico sugli assi ( z , y ) è un Exponential probability paper
Esempio: Supponiamo di avere
i seguenti tempi di vita
14,27,32,34,54,57,61,66,67,1
02,134, 152, 209, 230
1
0.9
0.8
0.7
0.6
Stimiamo la prob. di guasto
0.5
0.4
0.3
>> cdfe=(i-0.3)/14.4;
>> stairs(x,cdfe)
0.2
0.1
0
0
50
100
150
200
250
1
Stimiamo la affidabilità
>> stairs(x,1-cdfe)
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0
50
100
150
200
250
7
6
Anzicché lavorare puntualmente
stimiamo la densità di guasto con un
istogramma
5
4
3
>>c=[40,90,140,190]
>>fa=hist(x,c)
>>fa=[7,3,2,2]
2
1
0
40
90
140
190
Stimiamo la densità di guasto
>> classi=[0,65,115,165,230];
>> for i=1:4
dens(i)=fa(i)/14*1/(classi(i+1)-classi(i));
end
-3
8
x 10
7
6
5
4
3
2
40
60
80
100
120
140
160
180
200
Stimiamo il tasso di guasto
>> tasso =
>> cdf=cumsum(fa/14)
>> r(1)=1;
>>for i=2:4 r(i)=1-cdf(i-1); end
0.0077
0.0086
0.0100
0.0154
>> tasso=dens./r;
Regressione
>> i=[1:1:14];
>> cdfe=(i-0.3)/14.4;
>> y=log(1./(1-cdfe));
>> polytool(x,y)
3.5
3
2.5
2
1.5
1
0.5
>> beta
beta =
0.0124 -0.1517
0
-0.5
20
40
60
80
100
120
140
160
180
200
220
>> 1/expfit(x)
ans =
0.0113
>> ccf(:,1)=x';
>> lambda=beta(1);
>> ccf(:,2)=expcdf(x',1/lambda);
>> [h,p,kstat,cv]=kstest(x',ccf,0.05,0)
h = 0, p = 0.5004, kstat = 0.2126, cv = 0.3489
Nel caso gaussiano, conviene stimare la media e la deviazione standard
con i metodi usuali (media campionaria e varianza campionaria).
Normplot(x)
Normal Probability Plot
0.95
0.90
Probability
0.75
0.50
0.25
0.10
0.05
20
40
60
80
Data
100
120
140
Se sui dati dell’esempio sulla distribuzione di Weibull facessimo un normplot…
Test di ipotesi sulla media
Test di ipotesi sulla media
[h,pv,ci]=ttest(residuals,0)
h= 0
pv = 1
ci = -0.2223
0.2223
In entrambi i casi il campione casuale deve provenire da una popolazione
gaussiana. Cosa accade se la popolazione non è gaussiana?
• Si conta il numero di segni positivi (r+)
Esempio: Supponiamo di avere i seguenti tempi di vita
14,27,32,34,54,57,61,66,67,102,134, 152, 209, 230
Abbiamo dimostrato che la popolazione da cui provengono i dati segue una
legge esponenziale. Verifichiamo se la mediana è 100.
Abbiamo 5 segni positivi, che sono in numero inferiore a 7=14/2. Pertanto il
P-value risulta essere
p = 2P ( R ≤ r
+
+
)
>> p=2*binocdf(5,14,0.5)
>> p = signtest(x,100,0.05)
p=
p=
0.4240
0.4240
Si può sostituire
anche un vettore
Questo stesso test può essere utilizzato per verificare se due campioni casuali
hanno la stessa media.
Esempio: Supponiamo di avere i seguenti tempi di vita
14,27,32,34,54,57,61,66,67,102,134, 152, 209, 230
Abbiamo dimostrato che la popolazione da cui provengono i dati segue una
legge esponenziale di mediana 100. Consideriamo un secondo campione
7,29,30,31,44,47,69,70,75,105,134, 152, 210, 240.
>> x-y
ans =
7
-2
2
3
10
10
-8
-4
>> [P,H,STATS] = SIGNTEST(x,y,0.05)
P = 0.7744
H= 0
STATS = sign: 5
-8
-3
0
0
-1 -10
>> 2*binocdf(5,12,0.5)
ans =
0.7744
E se i campioni hanno taglia
diversa?
the Mann-Whitney U test
Test di Mann-Whitney (o della somma dei ranghi): due campioni indipendenti;
è uno dei test non parametrici più potenti e serve a verificare se due gruppi
indipendenti appartengono alla stessa popolazione. È un'alternativa molto valida
al test parametrico T-Student, quando non possono considerarsi le ipotesi del T
test, oppure la scala di misura è più debole di una scala ad intervalli.
Es. tempi di rottura di uno stesso componente in due ambienti diversi
Ambiente A
Rottura
Rango
1000
1
1300
5
1200
3
Ambiente B
Rottura
Rango
1400
6
1600
7
1180
2
1220
4
Nel caso di piccoli campioni, una volta calcolato la somma dei ranghi si confronta
questa con il valore teorico di una tavola in corrispondenza delle dimensioni dei
due campioni
n1n2 (n1 + n2 + 1)
n1 (n1 + n2 + 1)
σ
=
µ
=
u
u


12
2


U = min ∑ R1 , ∑ R2 
j
 i

U − µ u ± 0.5
σu
(+0.5 coda a sx; 0.5 coda a dx)
Calcolo di Skewness (asimmetria) e kurtosis (curtosi)
>> skewness(x)
ans =
0.9702
>> kurtosis(x)-3
ans =
-0.2990
Skewness (asimmetria) e
legge di Weibull
forma=0.5
forma=1.5
forma=2
1
0.8
Per β ∈ (0,2.6) ⇒ sk<0
0.6
0.4
forma=2.6
forma=3
forma=3.3
1.2
0.2
0
1
0
1
2
3
4
5
6
7
8
9
10
0.8
Per β ∈ (2.6,3.4) ⇒ sk ≈ 0
0.6
0.4
4
forma=4
forma=8
forma=12
0.2
3.5
0
0
0.5
1
1.5
2
2.5
3
3.5
4
3
2.5
2
Per β > 3.4 ⇒ sk>0
1.5
1
0.5
0
0
0.5
1
1.5
2
2.5
3
3.5
4
Si ottiene lo stesso risultato del test di Fisher?
INTERVALLI DI CONFIDENZA PER LA MEDIA DI UNA DISTRIBUZIONE
ESPONENZIALE
E' possibile dimostrare che lo stimatore di massima verosimiglianza
n
per λ è costituito dalla media campionaria ∑ X i / n.
i =1
n
n
∑ X i ∼ Gamma(n, λ )
2λ ∑ X i ∼ χ 22n
i =1
i =1
L'intervallo di confidenza è
 2
2
χ
χ
α
 α ,2 n
1− ,2 n
1 − α = P  n2
<λ < n2
2 X
2∑ X i
 ∑ i
i =1
 i =1






>> den=sum(x);
>> a1=chi2inv(0.025,28);
>> a2=chi2inv(0.975,28);
>> sup=a1/(2*den), inf=a2/(2*den)
sup = 0.062 inf = 0.0179
>>1/mean(x)
ans=0.0113
ESERCIZIO:
Considerando un albero rotante con tasso di guasto costante pari a 20 guasti
su 10^5 ore di funzionamento, calcolare il numero (medio) di alberi che si romperanno tra 10.000 e 20.000 ore, considerando un numero di alberi in prova pari a 50.
GUASTI MULTIPLI
Si assuma che il guasto del componente possa avvenire per effetto di due diversi
meccanismi di guasto. La densità di guasto in tal caso si calcola come
f (t ) = f1 (t ) R2 (t ) + f 2 (t ) R1 (t )
Esercizio: Il tempo di vita di un componente per effetto del meccanismo di guasto A è
descritto da una v.a. esponenziale di parametro 0.002 (per ore). Il tempo di vita dello
stesso componente per effetto del meccanismo di guasto B è descritto da una v.a. esponenziale di parametro 0.005 (per ore). Determinare densità di guasto complessiva e
funzione di affidabilità. (effettuare i grafici)
Esercizio: Supponiamo il tasso di guasto di un componente durante il periodo di
rodaggio sia di Weibull con parametro di forma 0.60 e parametro di scala 11000 ore,
durante il periodo di fase adulta sia di Weibull con parametro di forma 1,00 e parametro di scala 4000 ore e durante la fase di usura sia di Weibull con parametro di
forma 2.50 e parametro di scala 4500 ore. Scrivere la funzione densità di guasto
Scrivere la funzione densità di guasto
f (t ) = f1 (t ) R2 (t ) R3 (t ) + f 2 (t ) R1 (t ) R3 (t ) + f 3 (t ) R1 (t ) R2 (t )
Effettuare un grafico
Diversa è la situazione di un componente per il quale le informazioni circa il tasso di
guasto è ripartito lungo l’arco temporale. La affidabilità totale è una somma di affidabilità, ciascuna si riferisce a un periodo della vita, secondo dei pesi che restituiscono
il peso di quel tasso di guasto lungo l’arco temporale
R(t ) = a1 R1 (t ) + a2 R2 (t ) + a3 R3 (t )
con a1 + a2 + a3 = 1
f (t ) = a1 f1 (t ) + a2 f 2 (t ) + a3 f3 (t )
Iperesponenziale
F (t ) = a1 F1 (t ) + a2 F2 (t ) + a3 F3 (t )
MISTURE
Rank
Tempo guasto
# guasti
Guasti cum.
Median rank %
1
600
23
23
25.1
2
1200
9
32
35.1
3
1800
2
34
37.3
4
2400
3
37
40.6
5
3000
11
48
52.8
6
3600
6
54
59.4
7
4200
3
57
62.7
8
4800
6
63
69.4
9
5400
4
67
73.8
10
6600
3
70
77.1
11
6900
4
74
81.5
12
7200
8
82
90.4
13
7800
3
85
93.7
14
8400
0
85
93.7
15
9000
1
86
94.8
16
9600
2
88
97.0
17
10200
2
90
99.2
>> x=[23,9,2,3,11,6,3,6,4,3,4,8,3,0,1,2,2];
>> y=cumsum(x);
>> perc=(y-0.3)/(sum(x)+0.4);
Weibull Probability Plot
0.96
0.90
0.75
0.50
P robability
>> tguasto=[600,
1200,1800, 2400,
3000,3600, 4200,
4800,5400,6600,
6900,7200,7800,
8400,9000,9600,
10200]
>> weibplot(tguasto)
0.25
0.10
Attenzione!! Questo
grafico non è corretto.
Perché?
0.05
0.02
3
4
10
10
Data
>> beta
polytool(tguasto,y)
beta =
0.0071 24.4036
110
100
E se dividessimo
i 3 gruppi?
90
80
70
60
50
40
30
20
10
1000
2000
3000 4000
5000
6000 7000
8000
9000 10000
Rank
Tempo guasto
# guasti
Guasti cum.
Median rank %
1
600
23
23
25.1
2
1200
9
32
35.1
3
1800
2
34
37.3
4
2400
3
37
40.6
5
3000
11
48
52.8
6
3600
6
54
59.4
7
4200
3
57
62.7
8
4800
6
63
69.4
9
5400
4
67
73.8
10
6600
3
70
77.1
11
6900
4
74
81.5
12
7200
8
82
90.4
13
7800
3
85
93.7
14
8400
0
85
93.7
15
9000
1
86
94.8
16
9600
2
88
97.0
17
10200
2
90
99.2
Num.prova
Tempo guasto
# guasti
Guasti cum.
Median rank %
1
600
23
23
60.7
2
1200
9
32
84.8
3
1800
2
34
90.1
4
2400
3
37
98.1
5
3000
11
48
28.6
6
3600
6
54
44.7
7
4200
3
57
52.7
8
4800
6
63
68.7
9
5400
4
67
79.4
10
6600
3
70
87.4
11
6900
4
74
98.1
12
7200
8
82
47.0
13
7800
3
85
65.3
14
8400
0
85
65.3
15
9000
1
86
71.3
16
9600
2
88
83.6
17
10200
2
90
95.8
Num.prova
Tempo guasto
# guasti
Guasti cum.
Median rank %
1
600
23
23
60.7
2
1200
9
32
84.8
3
1800
2
34
90.1
4
2400
3
37
98.1
5
3000
11
11
28.6
6
3600
6
17
44.7
7
4200
3
20
52.7
8
4800
6
26
68.7
9
5400
4
30
79.4
10
6600
3
33
87.4
11
6900
4
37
98.1
12
7200
8
8
47.0
13
7800
3
11
65.3
14
8400
0
11
65.3
15
9000
1
12
71.3
16
9600
2
14
83.6
17
10200
2
16
95.8
>> y1=cumsum(x(1:4))
>>cdfe1=(y1-0.3)/ (sum(x(1:4))+0.4)
>> y2=cumsum(x(5:11))
>>cdfe2=(y2-0.3)/ (sum(x(5:11))+0.4)
>> y3=cumsum(x(12:17))
>>cdfe3=(y3-0.3)/ (sum(x(12:17))+0.4)
>> cdfe1 =
>> cdfe2 =
>> cdfe3 =
0.6070
0.2861
0.4695
0.8476
0.4465
0.6524
0.9011
0.5267
0.6524
0.9813
0.6872
0.7134
0.7941
0.8354
0.8743
0.9573
0.9813
Non possiamo usare weibplot. Facciamo il grafico a mano!!
>> ord1=log(log(1./(1-cdfe1)));
>> ord2=log(log(1./(1-cdfe2)));
>> ord3=log(log(1./(1-cdfe3)));
plot(log(tguasto(1:4)),ord1,'*',log(tguasto(5:11)),ord2,'*',log(tguasto(12:17)),ord3,'*')
Weibull probability papers
1.5
infanzia
adulto
usura
1
0.5
0
-0.5
-1
-1.5
6
6.5
7
7.5
8
8.5
9
Facciamo 3 regressioni diverse!!
9.5
>> polytool(log(tguasto(1:4)),ord1)
>> scale1=exp(-beta(2)/beta(1)), form=beta(1)
scale1 = 652.9613
form = 0.9800
>> polytool(log(tguasto(5:11)),ord2)
>> scale1=exp(-beta(2)/beta(1)), form=beta(1)
scale1 =
4.5330e+003
form = 2.6426
>> polytool(log(tguasto(12:17)),ord3)
>> scale1=exp(-beta(2)/beta(1)), form=beta(1)
scale1 = 1.3675e+004
form = 2.1832
β1
β2
 t 
 t 
 t 
37
37
16
R(t ) = exp  −  + exp  −  + exp  − 
90
90
90
 α1 
 α2 
 α3 
Effettuare un grafico
β3
La fase di raccolta dei dati sperimentali richiede che il prodotto – spesso un prototipo – sia
testato per tempi sufficientemente lunghi
Prove di affidabilità
Time-terminated
Failure-terminated
Le prove vengono eseguite in successione, ovvero sostituendo le unità che si sono guastate
con altre nuove fino allo scadere del tempo di prova.
Campionamento nel caso esponenziale
Failure-terminated
La prova viene fatta finire al primo guasto se questo risulta maggiore di un determinato
tempo di prova. In genere, dopo l’80% delle unità inizialmente inserite sulle macchine
di prova.
ta m ⋅ td
MTTF = =
r
r
r = unità che si sono guastate durante il test
N = m + ( r − 1)
ta = tempo complessivo di sperimentazione
l'ultima unità rotta non
td = tempo complessivo di prova di ogni macchina è stata rimpiazzata
m = macchine di prova
1
2
3
4
5
7
6
8
10
9
11
12
m = 6, r = 7, N = 6 + 6,td = 600, ta = 6*600
MTTF = (6*600) / 7 = 514
⇒ su 10^6 ore si hanno 1945 guasti
Caso failure-terminated: senza rimpiazzo delle unità
guastate
 m 
ta =  ∑ ti  + (m − r )td
 i =1 
Tempi di guasto delle r unità guastate
Tempi di funzionamento delle m-r unità non guastate
m = 8, r = 7,
td = 700, ta = 4180
MTTF=597ore
Caso time-terminated: con rimpiazzo delle unità guastate
m = 6, r = 8, N = 6 + 8, td = 700, ta = 6*700
MTTF = (6*700) / 8 = 625
Caso time-terminated: senza rimpiazzo delle unità
guastate
 m 
ta =  ∑ ti  + (m − r )td
 i =1 
Calcolare MTTF per esercizio
Caso gaussiano: failure-terminated senza rimpiazzo delle unità
guastate
Non è possibile usare dati provenienti da campioni interrotti. Le unità testate vengono
portate tutte al guasto. Vengono sostituite solo quelle che si guastano precocemente.
Campioni incompleti (Censored
sample)
Si definisce campione incompleto, un insieme di dati riguardanti sia componenti che non
hanno funzionato bene durante la prova, che componenti la cui prova è stata interrotta
senza che si sia raggiunto il guasto.
In questo caso non è possibile stimare la funzione di ripartizione perché i guasti sono
inframmezzati da unità testate non arrivate al guasto.
Assegnazione dell’INCREMENTO fittizio del numero ordinale, dipendente dalle
unità non arrivate al guasto.
NO = numero d'ordine
IO = incremento d'ordine
n +1− A
IO =
B +1
IO(1) = 1
A = NO del guasto precedente
B = numero di pezzi ancora da considerare incluso quello rotto
n = dimensione del campione di prova
10 + 1 − 1
IO(2) =
= 1.111
1+ 8
NO(m) = NO(m − 1) + IO(m)
NO (2) = NO(1) + IO(2) = 1 + 1.111 = 2.111
Non ci sono sospensioni
tra il secondo ed il terzo
guasto
IO(2) = IO(3)
NO (3) = NO(2) + IO(3) = 2.111 + 1.111 = 3.222
IO(4) =
10 + 1 − 3.222
= 1.296
1+ 5
E così via…
Num guasti
Tempo guasto
NO(i)
MEDIAN
RANK
1
21
1.0
0.0670
2
40
2.111
0.1739
3
66
3.222
0.2808
4
84
4.518
0.4055
5
150
6.678
0.6133
6
200
8.839
0.8212
NO(i ) − 0.3
n + 0.4
Fly UP