Affidabilità dei sistemi - Università degli Studi della Basilicata
by user
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 00 ∞ ∞ 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