Анализ надёжности и отказоустойчивости
artpgchАнализ надёжности дроссельных систем
Введение
Срок службы компонента — ключевой показатель его надёжности. Однако прямое измерение времени до отказа часто осложнено цензурированием данных: некоторые образцы не успевают выйти из строя за время испытаний. Как в этих условиях точно оценить долговечность изделия? В данной работе на примере анализа времени наработки на отказ 100 дроссельных заслонок мы исследуем методы обработки цензурированных данных. Сравнивая параметрические и непараметрические подходы, мы покажем, как выявить скрытые закономерности, такие как ранние отказы и эффект старения, для построения адекватных прогнозных моделей.
Присоединим необходимые библиотеки:
using Distributions, Random, Plots, Statistics, StatsBase, StatsPlots, Optim, KernelDensity
Особые свойства данных о наработке на отказ дроссельных заслонок
Анализ надёжности оперирует особыми данными: время жизни всегда положительно, а часть результатов цензурирована. Чтобы изучить это на практике, смоделируем испытания партии из 100 дросселей, создав сценарий с преобладанием долговечных изделий и наличием дефектной группы с ранними отказами.
Random.seed!(2)
lifetime = [rand(Weibull(3, 15000), 90); rand(Weibull(3, 1500), 10)]
В рамках моделирования использованы:
-
Ускоренные испытания (1 тестовый час = 100 часов эксплуатации)
-
Фиксированная длительность теста — 140 часов (эквивалент 14 000 часов)
-
Цензурирование данных для образцов, выдержавших весь испытательный цикл
T = 14000
obstime = sort(min.(T, lifetime))
Из-за ограниченной длительности испытаний для уцелевших дросселей мы можем указать лишь нижнюю границу времени наработки на отказ (>14 000 часов). В статистике надёжности такой тип данных носит название цензурированных наблюдений.
failed = obstime[obstime .< T]
nfailed = length(failed)
survived = obstime[obstime .== T]
nsurvived = length(survived)
censored = obstime .>= T
p1 = plot()
for i in 1:length(obstime)
plot!([0, obstime[i]], [i, i], color=:blue, line=:solid, label="")
end
if nsurvived > 0
for i in 1:nsurvived
plot!([T, 3e4], [nfailed + i, nfailed + i], color=:blue, line=:dot, label="")
end
end
plot!([T, T], [0, nfailed + nsurvived], color=:black, line=:solid, label="")
annotate!(T + 500, 30, text("<-- Неизвестное время работы до отказа", 10, :left))
xlabel!("Время работы до отказа")
ylabel!("Количество наблюдений")
xlims!(0, 3e4)
display(p1)
На данной диаграмме видно, что около 40% наших данных являются цензурированными при значении 14000 часов.
Какие выводы можно сделать по анализу распределения
Для анализа распределения времени жизни используются четыре ключевые функции:
-
Плотность вероятности — показывает вероятность отказа в конкретный момент времени.
-
Отказоустойчивость — вероятность того, что изделие проработает дольше заданного времени.
-
Интенсивность отказов — мгновенная скорость отказа для "выживших" к данному моменту изделий. Её рост указывает на эффект старения.
-
Вероятностный график — функция для визуальной проверки соответствия данных теоретической модели.
x = range(1, 30000, length=100)
p2 = plot(layout=(2,2), size=(800,600))
plot!(p2[1], x, pdf.(Weibull(2, 14000), x), label=nothing, linewidth=2)
plot!(p2[1], x, pdf.(Weibull(2, 18000), x), label=nothing, linewidth=2)
plot!(p2[1], x, pdf.(Weibull(1.1, 14000), x), label=nothing, linewidth=2)
title!(p2[1], "Плотность вероятности")
plot!(p2[2], x, ccdf.(Weibull(2, 14000), x), label=nothing, linewidth=2)
plot!(p2[2], x, ccdf.(Weibull(2, 18000), x), label=nothing, linewidth=2)
plot!(p2[2], x, ccdf.(Weibull(1.1, 14000), x), label=nothing, linewidth=2)
title!(p2[2], "Отказоустойчивочть")
wblhaz(x, a, b) = pdf(Weibull(b, a), x) / ccdf(Weibull(b, a), x)
plot!(p2[3], x, wblhaz.(x, 14000, 2), label=nothing, linewidth=2)
plot!(p2[3], x, wblhaz.(x, 18000, 2), label=nothing, linewidth=2)
plot!(p2[3], x, wblhaz.(x, 14000, 1.1), label=nothing, linewidth=2)
title!(p2[3], "Интенсивность отказов")
weibull_data = rand(Weibull(2, 14000), 40)
sorted_data = sort(weibull_data)
n = length(sorted_data)
empirical_probs = (1:n) ./ (n + 1)
theoretical_quantiles = quantile.(Weibull(2, 14000), empirical_probs)
scatter!(p2[4], theoretical_quantiles, sorted_data, marker=:xcross, markercolor=:blue, markersize=2, label=nothing)
plot!(p2[4], theoretical_quantiles, theoretical_quantiles, linecolor=:black, linestyle=:dash, linewidth=1.5, label=nothing)
title!(p2[4], "Вероятностный график")
#xlabel!(p2[4], "Данные")
#ylabel!(p2[4], "Вероятность")
plot!(p2[4], xticks=[], yticks=[])
display(p2)
Примеры четырёх типов графиков, с иллюстрацией распределения Вейбулла, которое часто используется для модели времени жизни.
Подбор параметров распределения
Распределение Вейбулла, в отличие от экспоненциального с постоянной интенсивностью отказов, позволяет моделировать "старение" оборудования — изменение вероятности отказа во времени. Наряду с ним для анализа надёжности применяются логнормальное, гамма и распределение Бирнбаума-Сандерса.
Для визуальной оценки данных строится эмпирическая функция распределения с 95% доверительными интервалами, показывающая накопленную долю отказов в каждый момент времени.
p3 = plot()
function ecdf_with_ci(times, censored; alpha=0.05)
sorted_indices = sortperm(times)
sorted_times = times[sorted_indices]
sorted_censored = censored[sorted_indices]
n = length(times)
F = zeros(n+1)
F_lo = zeros(n+1)
F_up = zeros(n+1)
x = zeros(n+1)
F[1] = 0.0
F_lo[1] = 0.0
F_up[1] = 0.0
x[1] = 0.0
at_risk = n
for i in 1:n
if !sorted_censored[i]
survival_prob = (at_risk - 1) / at_risk
F[i+1] = 1 - (1 - F[i]) * survival_prob
se = sqrt(F[i+1] * (1 - F[i+1]) / at_risk)
z = quantile(Normal(), 1 - alpha/2)
F_lo[i+1] = max(0, F[i+1] - z * se)
F_up[i+1] = min(1, F[i+1] + z * se)
else
F[i+1] = F[i]
F_lo[i+1] = F_lo[i]
F_up[i+1] = F_up[i]
end
at_risk -= 1
x[i+1] = sorted_times[i]
end
return x, F, F_lo, F_up
end
x_ecdf, empF, empFlo, empFup = ecdf_with_ci(obstime, censored)
function create_stairs(x, y)
x_stairs = Float64[]
y_stairs = Float64[]
for i in 1:length(x)
if i > 1
push!(x_stairs, x[i])
push!(y_stairs, y[i-1])
end
push!(x_stairs, x[i])
push!(y_stairs, y[i])
end
return x_stairs, y_stairs
end
x_stairs, F_stairs = create_stairs(x_ecdf, empF)
x_stairs_lo, F_stairs_lo = create_stairs(x_ecdf, empFlo)
x_stairs_up, F_stairs_up = create_stairs(x_ecdf, empFup)
plot!(x_stairs, F_stairs, color=:blue, linewidth=1.5, label=nothing)
plot!(x_stairs_lo, F_stairs_lo, color=:red, linestyle=:dot, linewidth=1, label=nothing)
plot!(x_stairs_up, F_stairs_up, color=:orange, linestyle=:dot, linewidth=1, label=nothing)
xlabel!("Время")
ylabel!("Соотношение отказов")
title!("Эмпирическая функция распределения")
display(p3)
На графике эмпирического распределения видно, что к 4000 часам отказало около 12% образцов. Из-за цензурирования 50% данных на отметке 14000 часов, функция распределения достигает лишь 0.53, не доходя до 1.0.
Для моделирования используем распределение Вейбулла — гибкий инструмент анализа надёжности. С помощью функции weibull_fit оценим параметры распределения с учетом цензурирования, построим теоретическую функцию распределения и вычислим для неё доверительные границы.
function weibull_fit(data, censored)
function weibull_loglik(θ, data, censored)
scale, shape = θ
if scale <= 0 || shape <= 0
return Inf
end
loglik = 0.0
for i in 1:length(data)
if !censored[i]
loglik += logpdf(Weibull(shape, scale), data[i])
else
loglik += logccdf(Weibull(shape, scale), data[i])
end
end
return -loglik
end
uncensored_data = data[.!censored]
if isempty(uncensored_data)
initial_scale = mean(data)
else
initial_scale = mean(uncensored_data)
end
initial_shape = 2.0
initial_guess = [initial_scale, initial_shape]
lower = [1.0, 0.1]
upper = [1e6, 100.0]
result = optimize(θ -> weibull_loglik(θ, data, censored), lower, upper, initial_guess, Fminbox(LBFGS()))
return result.minimizer
end
function weibull_cdf_ci(x, scale, shape, data, censored; alpha=0.05)
wblF = cdf.(Weibull(shape, scale), x)
n_uncensored = sum(.!censored)
z = quantile(Normal(), 1 - alpha/2)
wblFlo = similar(wblF)
wblFup = similar(wblF)
for i in 1:length(x)
se = z * sqrt(wblF[i] * (1 - wblF[i]) / n_uncensored)
wblFlo[i] = max(0.0, wblF[i] - se)
wblFup[i] = min(1.0, wblF[i] + se)
end
return wblF, wblFlo, wblFup
end
paramEsts = weibull_fit(obstime, censored)
xx = range(1, 2*T, length=500)
wblF, wblFlo, wblFup = weibull_cdf_ci(xx, paramEsts[1], paramEsts[2], obstime, censored)
Мы можем совместить графики эмпирической функции распределения и подобранной функции распределения, чтобы оценить, насколько хорошо распределение Вейбулла описывает надёжность дросселей.
p4 = plot()
plot!(x_stairs, F_stairs, color=:blue, linewidth=1.5, label=nothing)
plot!(xx, wblF, color=:red, linewidth=2, label=nothing)
plot!(xx, wblFlo, color=:red, linewidth=1, linestyle=:dash, label=nothing)
plot!(xx, wblFup, color=:red, linewidth=1, linestyle=:dash, label=nothing)
xlabel!("Время")
ylabel!("Вероятность отказа")
title!("Модели: Эмпирическая и Вейбулла")
xlims!(0, 2*T)
ylims!(0, 1)
display(p4)
Хотя модель Вейбулла позволяет экстраполировать прогнозы за пределы испытаний, её аппроксимация оказалась неудовлетворительной. Наблюдается систематическое расхождение: избыток ранних отказов до 2000 часов и дефицит отказов в интервале 7000-13000 часов относительно теоретической кривой.
Такое расхождение закономерно — оно отражает заложенную в данные специфику профиля отказов, которую классическое распределение Вейбулла не смогло адекватно описать.
Добавление сглаженной непараметрической оценки
Мы можем построить сглаженную непараметрическую кривую через эмпирическую с помощью функции kde. Мы удалим доверительные и добавим две кривые: одну с параметром сглаживания по умолчанию, и одну с параметром сглаживания, равным 1/3 от значения по умолчанию. Меньший параметр сглаживания позволяет кривой точнее следовать данным.
if !@isdefined(empF) || !@isdefined(x_ecdf)
empirical_cdf = ecdf(obstime)
x_ecdf = range(0, maximum(obstime), length=100)
empF = empirical_cdf.(x_ecdf)
end
kde_result = kde(obstime)
u = KernelDensity.default_bandwidth(obstime)
kde_standard = kde(obstime, bandwidth = u)
kde_small = kde(obstime, bandwidth = u/3)
npF_pdf = pdf(kde_standard, xx)
npF3_pdf = pdf(kde_small, xx)
npF_cdf = cumsum(npF_pdf) .* step(xx)
npF_cdf = npF_cdf ./ maximum(npF_cdf)
npF3_cdf = cumsum(npF3_pdf) .* step(xx)
npF3_cdf = npF3_cdf ./ maximum(npF3_cdf)
p5 = plot()
x_stairs, F_stairs = create_stairs(x_ecdf, empF)
plot!(x_stairs, F_stairs, color=:blue, linewidth=1.5, label="Эмпирическая")
plot!(xx, wblF, color=:red, linewidth=2, label="Вейбулла")
plot!(xx, npF_cdf, color=:green, linewidth=2, label="Непараметрическая, стандартная")
plot!(xx, npF3_cdf, color=:magenta, linewidth=2, label="Непараметрическая, 1/3 от стандартной")
xlims!(0, 1.3*T)
title!("Модели")
xlabel!("Время")
ylabel!("Вероятность")
display(p5)
Непараметрическая оценка идеально описывает данные в пределах наблюдений, но непригодна для прогнозирования за границей испытаний — её функция распределения резко обрывается на отметке 14000 часов.
Для дальнейшего анализа построим график интенсивности отказов, вычисленной по этой модели.
hzr = npF3_pdf ./ (1 .- npF3_cdf)
hzr[isinf.(hzr)] .= 0
hzr[isnan.(hzr)] .= 0
p6 = plot(xx, hzr, linewidth = 2, legend = false)
title!("Интенсивность отказов непараметрической модели")
xlims!(0, T)
ylims!(0, 0.00017)
display(p6)
Анализ кривой интенсивности отказов выявил классическую "ваннообразную" форму: пик в районе 1700 часов (ранние отказы), спад и последующий рост (износ). Такая картина характерна для компонентов с двумя уязвимыми периодами — начальной приработкой и конечным старением.
Выводы
Проведённый анализ показал, что классическое распределение Вейбулла недостаточно точно описывает сложный характер отказов дроссельных систем. Непараметрическая модель демонстрирует превосходное соответствие данным, но её главный недостаток — невозможность прогнозирования за пределами периода наблюдений.
Для решения этой проблемы предлагаются два перспективных направления:
-
Использование композитных вероятностных моделей, где отдельные компоненты моделируют ранние отказы и естественное старение
-
Разработка специализированных параметрических моделей, адаптированных к специфическим профилям надежности
Такой подход позволяет сохранить прогнозную силу параметрических методов, одновременно учитывая сложную природу реальных данных о надёжности.