Wienerprosess - AssetLife verktøykasse

Der det er tilgang til informasjon om komponentens degradering på en kontinuerlig skala (for eksempel overvåkning med sensor) kan modellering av degraderingsprosessen ved hjelp av en Wienerprosess være et godt alternativ. En fordel med Wienerprosessen er at det finnes enkle metoder for å modellere fremtidig sannsynlighet for svikt.

Dette kapittelet presenterer hvilke forutsetninger som må være på plass for å kunne benytte Wienerprosessen for degraderingsmodellering og estimering av fremtidig sviktsannsynlighet.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import invgauss
from scipy.stats import norm
from numpy.random import default_rng
rg = default_rng()
%config InlineBackend.figure_format = 'retina'

Generering av syntetiske tilstandsdata¶

Vi antar at degraderingen måles med en helseindikator (HI) x ved tidsinkrement t = [1,2,3, ... ,tau]. Videre har komponenten et definert sviktnivå ved x = L.

x0 = 0 # HI ved tid 0

nu_true = 0.1   # drift av degraderingsprosess
sigma_true = 1  # standardavvik for degraderingsprosess

tau = 1000
n = 10 # antall tidsserier

L = 40 # terskel for svikt

# Visualisering av n degraderingsforløp:
[plt.plot(np.cumsum([x0]+[rg.normal(loc = nu_true,scale = sigma_true) for t in range(tau)])) for i in range(n)];
plt.xlabel('Tid [dager]')
plt.ylabel('Helseindikator (HI)')
plt.hlines(y = L, xmin=0, xmax= tau, color = 'r',label = f'Sviktnivå, D ={L}');
plt.legend();

Vi velger ut ett degraderingsforløp som vi analyserer videre:

# Eksempel degraderingsforløp for videre analyse
arrData = np.cumsum([x0]+[rg.normal(loc = nu_true,scale = sigma_true) for t in range(tau)])

plt.plot(arrData, label = 'Eksempel degraderingsforløp');
plt.hlines(y = L, xmin=0, xmax= tau, color = 'r', label=f'Sviktnivå, D ={L}');
plt.xlabel('Tid [dager]')
plt.ylabel('Helseindikator (HI)');plt.legend();

Forutsetning for bruk av Wienerprosessen for degraderingsmodellering¶

Forutsetningene som er relevant i denne sammenhengen er at endringen i en måling av x til den neste

er stabile over tid
er uavhengige av hverandre
er tilnærmet normalfordelte

Sjek om endringen i degraderingsinkrement er stabil over tid¶

# lager en array med endringen i en måling til den neste:
arrDiff = np.diff(arrData)
arrDiff[:5]

plt.plot(arrDiff);
plt.hlines(y=0,xmin=0,xmax=len(arrDiff),color='k')
plt.title('Stabilitet av inkrementer over tid')
plt.xlabel('Tid [dager]')
plt.ylabel(r'$\Delta$(HI)');

Sjekk om påfølgende måling av helseindikator er statistisk uavhengige¶

For at denne typen statistisk modell skal kunne benyttes for å predikere levetidsfordeling må degraderingsinkrementene komme fra en uavhengig og identisk fordelt variabel (se Uavhengige, identisk fordelte variabler).

plt.acorr(arrDiff,usevlines = True,label = 'autocorrelation function (ACF) for degraderingsinkrementene') #,marker = 'o'
plt.xlim(-1,11)
plt.axhline(1.96/np.sqrt(tau), ls = '--', color = 'r',label='95% KI');
plt.axhline(-1.96/np.sqrt(tau), ls = '--', color = 'r');
plt.ylabel('ACF')
plt.xlabel('Lag/antall tidssteg');plt.legend();

Sjekk om inkrementene er normalfordelte:¶

from scipy.stats import norm

plt.ylabel('Antall');plt.xlabel('Inkrement størrelse');
counts, bins, _ = plt.hist(arrDiff, bins=40, alpha=0.6, label='Data')
mu, sigma = norm.fit(arrDiff)
x = np.linspace(min(arrDiff), max(arrDiff), 200)
bin_width = bins[1] - bins[0]

plt.plot(
    x,
    norm.pdf(x, mu, sigma) * len(arrDiff) * bin_width,
    'r-',
    lw=2,
    label='Tilpasset normalfordeling'
)
plt.legend();plt.show()

Parameterestimering for degraderingsprosessen¶

# Gjennomsnittet for degraderingsinkrementene betegnes som degraderingsprosessens drift.
# Driftparameteren kan regnes ut ved å ta gjennomsnittet av degraderingsinkrementene
delta_x_bar = np.mean(arrDiff)
print("delta_x_bar =",np.round(delta_x_bar,3))

delta_x_bar = 0.126

MTTF = L/delta_x_bar # Mean Time To Failure
print(f'Gjennomsnittlig tid-til-svikt: {np.round(MTTF,0)}')

Gjennomsnittlig tid-til-svikt: 317.0

# Diffusjonskoeffisienten for Wienerprosessen kan estimeres ved å regne ut standardavvikte for degraderingsprosessen.
delta_x_std = np.std(arrDiff)

list_norm = np.linspace(-3*sigma_true,3*sigma_true+nu_true,100)

plt.plot(list_norm,[norm.pdf(x,nu_true,sigma_true) for x in list_norm],label = 'true');
plt.plot(list_norm,[norm.pdf(x,delta_x_bar,delta_x_std) for x in list_norm],label = 'estimert');
plt.hist(arrDiff,density = True,bins = 50);

plt.vlines(x = 0, ymin = -1, ymax= 1,color= 'k', ls = '--');
plt.ylim(0,0.6)
plt.title('Endring i HI fra et tidspunkt til neste')
plt.legend()

Fordeling restlevetid¶

En fordel med Wienerprosessen er at tid til første kryssing av en definert terskel (L) følger invers Gauss fordeling (IG). Hvis svikt er definert når HI overstiger terskel L, med andre ord, svikt inntreffer første gang X(t)≥L, så kan fordelingen for svikttidspunktet modelleres som

T \sim IG(t;μ,λ)

(1)

hvor $μ=(L-x_0)/ν$ og $λ=(L-x_0)^2/σ^2.$

Dette er beskrevet nærmere i prosjektnotatet “Alternativer for videreutvikling av modeller for fremtidig sviktsannsynlighet”. For mer informasjon om IG-fordelingen se: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.invgauss.html

mu_true = (L-x0)/nu_true
lmbda_true = (L-x0)**2/sigma_true**2

Estimering av drift og diffusjonskoeffisient basert på degraderingsforløpet¶

mu_bar = (L-x0)/delta_x_bar
lmbda_bar = (L-x0)**2/delta_x_std**2

list_t = np.arange(0,tau+1)

plt.figure()
plt.plot(list_t,invgauss.pdf(list_t,mu = mu_true/lmbda_true, scale=lmbda_true),label = 'PDF IG faktisk');
plt.plot(list_t,invgauss.pdf(list_t,mu = mu_bar/lmbda_bar, scale=lmbda_bar),label = 'PDF IG estimert');

plt.xlabel('tid [dager]')
plt.ylabel('f(t)');plt.legend()

plt.figure()
plt.plot(list_t,invgauss.cdf(list_t,mu = mu_true/lmbda_true, scale=lmbda_true),label = 'CDF IG faktisk');
plt.plot(list_t,invgauss.cdf(list_t,mu = mu_bar/lmbda_bar, scale=lmbda_bar),label = 'CDF IG estimert');

plt.xlabel('tid [dager]')
plt.ylabel('F(t)');plt.legend();