基于维纳过程模型的剩余使用寿命预测 1.蒙特卡洛方法模拟部件的退化轨迹代码 2.线性维纳模型预测剩余使用寿命的代码及文章参考 3.非线性维纳模型预测剩余使用寿命的代码及文章参考 4.MLE估计算法代码 5.卡尔曼滤波算法更新参数代码。 6.贝叶斯参数更新

蒙特卡洛造点退化数据

先整点仿真数据最实在。维纳过程的核心是漂移+随机波动,咱们用numpy搞个退化轨迹生成器:

import numpy as np

def generate_degradation(drift=0.2, diffusion=0.5, steps=100):
    t = np.arange(steps)
    delta_t = t[1] - t[0]
    increments = drift*delta_t + diffusion*np.sqrt(delta_t)*np.random.randn(steps)
    return np.cumsum(increments)

# 生成10条退化曲线
import matplotlib.pyplot as plt
for _ in range(10):
    plt.plot(generate_degradation())
plt.title("蒙特卡洛退化轨迹")
plt.show()

这段代码的关键在于时间增量的处理——漂移项线性增长,扩散项带着随机噪声。跑出来的曲线那叫一个妖娆,跟实际设备磨损的随机性特匹配。

线性维纳玩转寿命预测

假设退化量超过阈值3.0就算报废,参考Si大神的经典模型[1],咱们可以直接计算首次穿越时间的概率密度:

from scipy.stats import norm

def rul_pdf(t, current_degradation, drift, diffusion):
    threshold = 3.0
    mu = (threshold - current_degradation)/drift
    sigma = diffusion*np.sqrt(t)/drift
    return norm.pdf(t, loc=mu, scale=sigma)

# 示例:当前退化量1.5时的RUL分布
t_range = np.linspace(0, 20, 100)
plt.plot(t_range, rul_pdf(t_range, 1.5, 0.2, 0.5))
plt.title("线性维纳RUL概率密度")
plt.show()

这里用到了首次穿越时间的解析解,计算效率高到飞起。但注意这假设参数已知,现实中咱得接着往下看参数估计。

非线性场景怎么办

设备退化常有加速现象,这时候得用Wang的指数模型[2]。改动生成函数:

def nonlinear_wiener(t_max=100):
    t = np.arange(t_max)
    return 0.001*(t**2) + 0.1*np.cumsum(np.random.randn(t_max))

# 参数估计得用极大似然(后面讲)

参数估计实战

历史数据在手,MLE大法搞起。写个负对数似然函数扔给优化器:

from scipy.optimize import minimize

def mle_estimate(data, dt=1):
    def neg_log_likelihood(params):
        mu, sigma = params
        n = len(data)-1
        diff = np.diff(data)
        return n*np.log(sigma*np.sqrt(dt)) + 0.5*np.sum((diff - mu*dt)**2)/(sigma**2*dt)
    
    result = minimize(neg_log_likelihood, [0.1, 0.5])
    return result.x

# 测试估计效果
data = generate_degradation(steps=500)
print("真实参数(0.2,0.5),估计结果:", mle_estimate(data))

优化结果受初值影响大,多跑几次取稳定值。工业现场数据量大时这个方法挺抗造。

基于维纳过程模型的剩余使用寿命预测 1.蒙特卡洛方法模拟部件的退化轨迹代码 2.线性维纳模型预测剩余使用寿命的代码及文章参考 3.非线性维纳模型预测剩余使用寿命的代码及文章参考 4.MLE估计算法代码 5.卡尔曼滤波算法更新参数代码。 6.贝叶斯参数更新

在线更新用卡尔曼

设备运行时实时数据来了,上卡尔曼滤波动态调整参数。这里把漂移系数当作状态变量:

def kalman_update(prior_mean, prior_var, measurement, dt=1, obs_noise=0.1):
    # 预测步
    pred_mean = prior_mean
    pred_var = prior_var + 0.01  # 过程噪声
    
    # 更新步
    residual = measurement - pred_mean*dt
    K = pred_var*dt / (pred_var*(dt**2) + obs_noise)
    post_mean = pred_mean + K*residual
    post_var = (1 - K*dt)*pred_var
    return post_mean, post_var

# 模拟在线更新过程
true_drift = 0.22
measurements = true_drift + np.random.normal(0, 0.1, 50)
estimate = []
mu, var = 0.2, 0.5
for z in measurements:
    mu, var = kalman_update(mu, var, z)
    estimate.append(mu)
plt.plot(estimate, label='估计值')
plt.axhline(true_drift, color='r', linestyle='--', label='真实值')
plt.legend()
plt.show()

卡尔曼增益自动调节信任程度,数据越多估计越准,这个特性在寿命预测中贼有用。

贝叶斯派怎么玩

最后整个PyMC3实现贝叶斯参数更新,适合小数据场景:

import pymc3 as pm

with pm.Model() as bayesian_model:
    mu = pm.Normal('mu', mu=0.2, sigma=0.5)
    sigma = pm.HalfNormal('sigma', sigma=1)
    likelihood = pm.Normal('obs', mu=mu*1, sigma=sigma, observed=measurements)
    trace = pm.sample(1000, tune=500)
    
pm.plot_posterior(trace)
plt.show()

后验分布直接给出参数的不确定性,做风险决策时比单点估计靠谱多了。

总结

从生成数据到参数动态更新,这套组合拳基本覆盖了维纳过程做寿命预测的核心流程。实际应用时要特别注意模型假设验证——漂移系数是否真恒定?扩散项是不是正态?这些直接关系到预测准头。代码把玩起来,有坑欢迎来唠!

[1] Si X S, et al. IEEE Trans Reliab, 2012

[2] Wang X, et al. Mech Syst Signal Process, 2014

更多推荐