Наблюдение за насосом и выявление его поломок

Насосы являются существенным оборудованием во многих отраслях промышленности включая степень и химический, минерал и горная промышленность, производство, нагревание, кондиционирование воздуха и охлаждение. Центробежные насосы используются, чтобы транспортировать жидкости преобразованием вращательной кинетической энергии к гидродинамической энергии потока жидкости. Вращательная энергия обычно прибывает из двигателя внутреннего сгорания или электродвигателя. Жидкость вводит рабочее колесо насоса вперед или близко к вращающейся оси и ускоряется рабочим колесом, текущий радиально исходящий в диффузор.

Насосы испытывают повреждение своих гидравлических или механических компонентов. Наиболее часто дающие сбой компоненты являются скользящими кольцевыми изоляциями и шарикоподшипниками, несмотря на то, что отказы к другим компонентам включая ведущий двигатель, лопатки рабочего колеса и скользящие подшипники также весьма распространены. Таблица ниже приводит наиболее распространенные типы отказов.

Доступные датчики

Следующие сигналы обычно измеряются:

Математические модели системы насоса и трубопровода

Крутящий момент$$M$$примененный ротор радиального центробежного насоса приводит к скорости вращения $$\omega$$ и передает увеличение импульса жидкости насоса от входа ротора меньшего радиуса к выходу ротора большего радиуса. Турбинные уравнения Эйлера дают к отношению между перепадом давления $$\Delta p$$, скорость $$\omega$$ и жидкость разряжает уровень (скорость потока жидкости) $$Q$$:

где $$H_{th}=\frac{\Delta p}{\rho g}$$ теоретическое (идеал; без потерь) крышка насоса, измеренная в метрах и ${\mathit{h}}_1$, ${\mathit{h}}_2$ коэффициенты пропорциональности. При составлении конечного числа лопаток рабочего колеса, потерь на трение и потерь удара из-за нетангенциального потока, действительной крышкой насоса дают:

где ${\mathit{h}}_{\mathrm{nn}}$, ${\mathit{h}}_{\mathrm{n}\mathit{v}}$ и ${\mathit{h}}_{\mathrm{vv}}$ коэффициенты пропорциональности должны быть обработаны как параметры модели. Соответствующий крутящий момент насоса:

Механические детали двигателя и насоса заставляют скорость увеличиваться, когда крутящий момент применяется согласно:

где $$J_P$$ отношение инерции двигателя и насоса, и $$M_f$$ фрикционный крутящий момент, состоящий из трения Кулона $$M_{f0}$$ и вязкое трение $$M_{f1}\omega (t)$$ согласно:

Насос соединяется с трубопроводной системой, которая транспортирует жидкость от более низкого резервуара для хранения до верхнего. Выражения уравнения баланса импульса:

где $$h_{rr}$$ коэффициент сопротивления трубопровода, $$a_F=\frac{l}{gA}$$ с длиной трубопровода $$l$$ и площадь поперечного сечения $$A$$, и $$H_{static}$$ высота устройства хранения данных по насосу. Параметры модели $$h_{nn},h_{nv},h_{vv}$$ или известны от физики или может быть оценен путем подбора кривой измеренным сигналам датчика к входным параметрам/выходным параметрам модели. Тип используемой модели может зависеть от условий работы, под которыми запущен насос. Например, полная нелинейная модель системы трубопровода насоса не может требоваться, если насос всегда запускается на постоянной угловой скорости.

Методы обнаружения отказа

Отказы могут быть обнаружены путем исследования определенных функций, извлеченных из измерений и сравнения их с известными порогами приемлемого поведения. Обнаружительная способность и изонеустойчивость различных отказов зависят от природы эксперимента и доступности измерений. Например, анализ постоянной скорости с измерениями давления только может обнаружить отказы, вызывающие большие скачки давления. Кроме того, это не может надежно оценить причину отказа. Однако многоскоростной эксперимент с измерениями перепада давления, крутящего момента двигателя и скорости потока жидкости может обнаружить и изолировать много источников отказов, таких как те, которые происходят из газовых корпусов, пробного прогона, больших депозитов, моторные дефекты и т.д.

Основанный на модели подход использует следующие методы:

Экспериментирование постоянной скорости: анализ отказа по оценке параметра

Установившаяся практика для калибровки насоса и контроля должна запустить его на постоянной скорости и записать статический главный и жидкий уровень выброса насоса. Меняя положение клапана в трубопроводной системе, жидкий объем выброса (GPM) отрегулирован. Увеличение уровня выброса заставляет крышку насоса уменьшаться. Измеренные главные характеристики насоса могут быть сравнены с поставляемыми изготовителем значениями. Любые различия указали бы на возможность отказов. Измерения для головы доставки и уровня выброса потока были получены симуляциями системной модели трубопровода насоса в Simulink.

На номинальной скорости 2 900 об/мин идеальные характеристики крышки насоса для здорового насоса, предоставленного производителем, как показано.

url = 'https://www.mathworks.com/supportfiles/predmaint/fault-diagnosis-of-centrifugal-pumps-using-steady-state-experiments/PumpCharacteristicsData.mat';
websave('PumpCharacteristicsData.mat',url);
load PumpCharacteristicsData Q0 H0 M0 % manufacturer supplied data for pump's delivery head
figure
plot(Q0, H0, '--');  
xlabel('Discharge Rate Q (m^3/h)')
ylabel('Pump Head (m)')
title('Pump Delivery Head Characteristics at 2900 RPM')
grid on
legend('Healthy pump')

Отказы, которые вызывают заметное изменение в характеристиках насоса:

Для анализа неисправных насосов скорость, крутящий момент и измерения скорости потока жидкости были собраны для насосов, затронутых различными отказами. Например, когда введенный отказ находится в звонке разрешения, измеренные главные характеристики для насосов показывают, что ясное переключает характеристическую кривую на нижний регистр.

load PumpCharacteristicsData Q1 H1 M1  % signals measured for a pump with a large clearance gap
hold on
plot(Q1, H1);  

load PumpCharacteristicsData Q2 H2 M2  % signals measured for a pump with a small clearance gap
plot(Q2, H2);  
legend('Healthy pump','Large clearance','Small clearance')
hold off

Подобные изменения видны в характеристиках потока крутящего момента и для других типов отказа.

Для автоматизации диагностики отказа вы поворачиваете наблюдаемые изменения в количественных данных. Надежный способ сделать так состоит в том, чтобы соответствовать параметрированной кривой к данным о характеристике главного потока, отображенным на графике выше. Используя управляющие уравнения динамики трубопровода насоса и использование упрощенного отношения крутящего момента, получены следующие уравнения:

${\mathit{h}}_{\mathrm{nn}},{\mathit{h}}_{\mathrm{nv}},{\mathit{h}}_{\mathrm{vv}},{\mathit{k}}_0,{\mathit{k}}_1,{\mathit{k}}_2$ параметры должны быть оценены. Если вы измеряетесь $\omega ,\mathit{Q},\mathit{H}$и ${\mathit{M}}_{\mathit{p}}$, параметр может быть оценен линейным методом наименьших квадратов. Эти параметры являются функциями, которые могут быть использованы, чтобы разработать обнаружение отказа и алгоритм диагноза.

Предварительный Анализ: Сравнение значений параметров

Вычислите и постройте значения параметров, оцененные для вышеупомянутых 3 кривых. Используйте измеренные значения $\mathit{Q},\mathit{H}$ и ${\mathit{M}}_{\mathit{p}}$ как данные и $$\omega =2900RPM$$ как номинальная скорость насоса.

w = 2900; % RPM
% Healthy pump
[hnn_0, hnv_0, hvv_0, k0_0, k1_0, k2_0] = linearFit(0, {w, Q0, H0, M0});
% Pump with large clearance
[hnn_1, hnv_1, hvv_1, k0_1, k1_1, k2_1] = linearFit(0, {w, Q1, H1, M1});
% Pump with small clearance
[hnn_2, hnv_2, hvv_2, k0_2, k1_2, k2_2] = linearFit(0, {w, Q2, H2, M2});
X = [hnn_0 hnn_1 hnn_2; hnv_0  hnv_1  hnv_2; hvv_0  hvv_1  hvv_2]';
disp(array2table(X,'VariableNames',{'hnn','hnv','hvv'},...
    'RowNames',{'Healthy','Large Clearance', 'Small Clearance'}))

                          hnn           hnv           hvv   
                       __________    __________    _________

    Healthy            5.1164e-06    8.6148e-05     0.010421
    Large Clearance     4.849e-06     8.362e-05     0.011082
    Small Clearance    5.3677e-06    8.4764e-05    0.0094656

Y = [k0_0 k0_1 k0_2; k1_0  k1_1  k1_2; k2_0  k2_1  k2_2]';
disp(array2table(Y,'VariableNames',{'k0','k1','k2'},...
    'RowNames',{'Healthy','Large Clearance', 'Small Clearance'}))

                           k0           k1           k2    
                       __________    ________    __________

    Healthy            0.00033347    0.016535    2.8212e-07
    Large Clearance    0.00031571    0.016471    3.0285e-07
    Small Clearance    0.00034604    0.015886    2.6669e-07

Таблицы показывают это ${\mathit{h}}_{\mathrm{nn}}$ и ${\mathit{k}}_0$ значения уменьшают, когда разрыв разрешения является большим, в то время как они больше, чем номинальная стоимость для маленького разрешения. С другой стороны, ${\mathit{h}}_{\mathrm{vv}}$ и ${\mathit{k}}_2$ значения увеличиваются для большого разрыва разрешения и уменьшения для маленького разрыва. Зависимость ${\mathit{h}}_{\mathrm{nv}}$ и ${\mathit{k}}_1$ на разрешении разрыв менее ясен.

Слияние неопределенности

Предварительный анализ показал, как изменения параметра могут указать на отказ. Однако даже для здоровых насосов существуют изменения измерений вследствие шума измерения, жидкого загрязнения и изменений вязкости и характеристик крутящего момента промаха двигателя, запускающего насос. Эти изменения измерения вводят неопределенность в оценках параметра.

Соберите 5 наборов измерений от насоса, действующего при условии без отказов путем выполнения его на уровне 2 900 об/мин для 10 положений клапана дросселя выброса.

url = 'https://www.mathworks.com/supportfiles/predmaint/fault-diagnosis-of-centrifugal-pumps-using-steady-state-experiments/FaultDiagnosisData.mat';
websave('FaultDiagnosisData.mat',url);
load FaultDiagnosisData HealthyEnsemble
H = cellfun(@(x)x.Head,HealthyEnsemble,'uni',0);
Q = cellfun(@(x)x.Discharge,HealthyEnsemble,'uni',0);
plot(cat(2,Q{:}),cat(2,H{:}),'k.')
title('Pump Head Characteristics Ensemble for 5 Test Runs')
xlabel('Flow rate (m^3/hed)')
ylabel('Pump head (m)')

График показывает изменение характеристик даже для здорового насоса при реалистических условиях. Эти изменения должны быть учтены для того, чтобы сделать диагностику отказа надежной. Следующие разделы обсуждают обнаружение отказа и методы изоляции для зашумленных данных.

Обнаружение аномалии

Во многих ситуациях измерения только здоровых машин доступны. В этом случае статистическое описание здорового состояния, инкапсулировавшего средним значением и ковариацией вектора параметра, может быть создано с помощью доступных измерений. Измерения тестового насоса могут быть сравнены с номинальной статистикой, чтобы протестировать, вероятно ли, что тестовый насос является здоровым насосом. Неисправный насос ожидается к обнаруженному как аномалия в представлении функций обнаружения.

Оценочное среднее значение и ковариация крышки насоса и параметров крутящего момента.

load FaultDiagnosisData HealthyEnsemble
[HealthyTheta1, HealthyTheta2] = linearFit(1, HealthyEnsemble);
meanTheta1 = mean(HealthyTheta1,1);
meanTheta2 = mean(HealthyTheta2,1);
covTheta1  = cov(HealthyTheta1);
covTheta2  = cov(HealthyTheta2);

Визуализируйте неопределенность параметра как 74% областей доверия, который соответствует 2 стандартным отклонениям ($$\sqrt{chi2inv(0.74,3)}\approx 2$$). Смотрите, что помощник функционирует helperPlotConfidenceEllipsoid для деталей.

% Confidence ellipsoid for pump head parameters
f = figure;
f.Position(3) = f.Position(3)*2;
subplot(121)
helperPlotConfidenceEllipsoid(meanTheta1,covTheta1,2,0.6);   
xlabel('hnn')
ylabel('hnv')
zlabel('hvv')
title('2-sd Confidence Ellipsoid for Pump Head Parameters')
hold on

% Confidence ellipsoid for pump torque parameters
subplot(122)
helperPlotConfidenceEllipsoid(meanTheta2,covTheta2,2,0.6); 
xlabel('k0')
ylabel('k1')
zlabel('k2')
title('2-sd Confidence Ellipsoid for Pump Torque Parameters')
hold on

Серые эллипсоиды показывают области доверия здоровых параметров насоса. Загрузите непомеченные тестовые данные для сравнения со здоровой областью.

load FaultDiagnosisData TestEnsemble

TestEnsemble содержит набор скорости насоса, крутящего момента, головы и измерений скорости потока жидкости в различных положениях клапана. Все измерения содержат отказ разрыва разрешения различных величин.

% Test data preview
disp(TestEnsemble{1}(1:5,:)) % first 5 measurement rows from the first ensemble member

      Time       Run    ValvePosition    Speed      Head     Discharge    Torque 
    _________    ___    _____________    ______    ______    _________    _______

    180 sec       1          10          3034.6    12.367     35.339      0.35288
    180.1 sec     1          10          2922.1    9.6762     36.556       4.6953
    180.2 sec     1          10          2636.1    11.168     36.835       9.8898
    180.3 sec     1          10          2717.4    10.562      40.22      -12.598
    180.4 sec     1          10          3183.7     10.55     40.553       14.672

Вычислите тестовые параметры. Смотрите, что помощник функционирует linearFit.

% TestTheta1: pump head parameters
% TestTheta2: pump torque parameters
[TestTheta1,TestTheta2] = linearFit(1, TestEnsemble);
subplot(121)
plot3(TestTheta1(:,1),TestTheta1(:,2),TestTheta1(:,3),'g*')
view([-42.7 10])
subplot(122)
plot3(TestTheta2(:,1),TestTheta2(:,2),TestTheta2(:,3),'g*')
view([-28.3 18])

Каждый зеленый маркер-звездочка внесен одним тестовым насосом. Маркеры, которые находятся вне доверительных границ, могут быть обработаны как выбросы, в то время как те внутри или от здорового насоса или вышли из обнаружения. Обратите внимание на то, что маркер от конкретного насоса может быть отмечен как аномальный в представлении крышки насоса, но не в представлении крутящего момента насоса. Это могло быть вследствие других источников отказов, обнаруживаемых этими представлениями или базовой надежностью давления, и закрутить измерения.

Определение количества обнаружения аномалии Используя области доверия

В этом разделе обсужден способ использовать информацию об области доверия для обнаружения и оценить серьезность отказов. Метод должен вычислить "расстояние" тестовой выборки от среднего значения или медианы здорового распределения области. Расстояние должно быть относительно нормального "распространения" здоровых данных о параметре, представленных его ковариацией. Функциональный MAHAL вычисляет расстояние Mahalanobis тестовых выборок от распределения ссылочного демонстрационного набора (здоровый набор параметров насоса здесь):

ParDist1 = mahal(TestTheta1, HealthyTheta1);  % for pump head parameters

Если вы принимаете 74%-ю доверительную границу (2 стандартных отклонения) как приемлемое изменение для здоровых данных, любые значения в ParDist1, которые больше 2^2 = 4, должны быть помечены как аномальные и следовательно показательные из дефектного поведения.

Добавьте значения расстояния в график. Красные линии отмечают аномальные тестовые выборки. Смотрите, что помощник функционирует helperAddDistanceLines.

Threshold = 2;
disp(table((1:length(ParDist1))',ParDist1, ParDist1>Threshold^2,...
    'VariableNames',{'PumpNumber','SampleDistance','Anomalous'}))

    PumpNumber    SampleDistance    Anomalous
    __________    ______________    _________

         1            58.874          true   
         2            24.051          true   
         3             6.281          true   
         4            3.7179          false  
         5             13.58          true   
         6            3.0723          false  
         7            2.0958          false  
         8            4.7127          true   
         9            26.829          true   
        10           0.74682          false  

helperAddDistanceLines(1, ParDist1, meanTheta1, TestTheta1, Threshold);

Так же для крутящего момента насоса:

ParDist2 = mahal(TestTheta2, HealthyTheta2);  % for pump torque parameters
disp(table((1:length(ParDist2))',ParDist2, ParDist2>Threshold^2,...
    'VariableNames',{'PumpNumber','SampleDistance','Anomalous'}))

    PumpNumber    SampleDistance    Anomalous
    __________    ______________    _________

         1            9.1381          true   
         2            5.4249          true   
         3            3.0565          false  
         4             3.775          false  
         5           0.77961          false  
         6            7.5508          true   
         7            3.3368          false  
         8           0.74834          false  
         9            3.6478          false  
        10            1.0241          false  

helperAddDistanceLines(2, ParDist2, meanTheta2, TestTheta2, Threshold);
view([8.1 17.2])

Графики теперь не только показывают обнаружение аномальных выборок, но также и определяют количество их серьезности.

Определение количества обнаружения аномалии Используя классификатор одного класса

Другой эффективный метод, чтобы отметить аномалии должен создать классификатор одного класса для здорового набора данных параметра. Обучите классификатор SVM с помощью здоровых данных о параметре насоса. С тех пор нет никаких используемых меток отказа, не обрабатывают все выборки как прибывающий из того же (здорового) класса. Начиная с изменений в параметрах ${\mathit{h}}_{\mathrm{nn}}$ и ${\mathit{h}}_{\mathrm{vv}}$ являются самыми показательными из потенциальных отказов, используют только эти параметры для обучения классификатор SVM.

nc = size(HealthyTheta1,1);
rng(2)  % for reproducibility
SVMOneClass1 = fitcsvm(HealthyTheta1(:,[1 3]),ones(nc,1),...
    'KernelScale','auto',...
    'Standardize',true,...
    'OutlierFraction',0.0455);

Постройте тестовые наблюдения и контур решения. Отметьте векторы поддержки и потенциальные выбросы. Смотрите, что помощник функционирует helperPlotSVM.

figure
helperPlotSVM(SVMOneClass1,TestTheta1(:,[1 3]))
title('SVM Anomaly Detection for Pump Head Parameters')
xlabel('hnn')
ylabel('hvv')

Контур, разделяющий выбросы от остальной части данных, происходит, где значение контура 0; это - кривая уровня, отмеченная "0" в графике. Выбросы отмечены красными кругами. Подобный анализ может быть выполнен для параметров крутящего момента.

SVMOneClass2 = fitcsvm(HealthyTheta2(:,[1 3]),ones(nc,1),...
    'KernelScale','auto',...
    'Standardize',true,...
    'OutlierFraction',0.0455);
figure
helperPlotSVM(SVMOneClass2,TestTheta2(:,[1 3]))
title('SVM Anomaly Detection for Torque Parameters')
xlabel('k0')
ylabel('k2')

Подобный анализ может быть выполнен для обнаружения других видов отказов, таких как износ или депозиты при выходе рабочего колеса, как обсуждено затем в контексте изоляции отказа.

Изоляция отказа Используя установившиеся параметры как функции

Если информация о типе отказа (отказов) в системе тестирования доступна, это может использоваться, чтобы создать алгоритмы, которые не только обнаруживают отказы, но также и указывают на их тип.

A. Различение отказов разрыва разрешения тестами отношения правдоподобия

Изменения в разрыве разрешения могут быть разделены в два типа - меньший, чем ожидаемый разрыв (желтая линия в графике характеристики крышки насоса) и больше, чем ожидаемый разрыв (красная линия). Загрузите тестовые наборы данных насоса, содержащие отказы разрыва разрешения, где природа (большого или маленького) отказа известна заранее. Используйте эти метки отказа, чтобы выполнить классификацию с 3 путями среди следующих режимов:

url = 'https://www.mathworks.com/supportfiles/predmaint/fault-diagnosis-of-centrifugal-pumps-using-steady-state-experiments/LabeledGapClearanceData.mat';
websave('LabeledGapClearanceData.mat',url);
load LabeledGapClearanceData HealthyEnsemble LargeGapEnsemble SmallGapEnsemble

Ансамбли содержат данные из 50 независимых экспериментов. Подбирайте установившиеся линейные модели как прежде, чтобы параметрировать данные о крутящем моменте и крышка насоса.

[HealthyTheta1, HealthyTheta2] = linearFit(1,HealthyEnsemble);
[LargeTheta1, LargeTheta2]     = linearFit(1,LargeGapEnsemble);
[SmallTheta1, SmallTheta2]     = linearFit(1,SmallGapEnsemble);

Постройте гистограммы параметра, чтобы проверять, существует ли отделимость среди этих 3 режимов. Функциональный histfit используется, чтобы построить гистограммы и соответствующие подходящие кривые нормального распределения. Смотрите, что помощник функционирует helperPlotHistogram.

Параметры крышки насоса:

helperPlotHistogram(HealthyTheta1, LargeTheta1, SmallTheta1, {'hnn','hnv','hvv'})

Гистограмма показывает это ${\mathit{h}}_{\mathrm{nn}}$ предложения хорошая отделимость среди этих трех режимов, но ${\mathit{h}}_{\mathrm{nv}},{\mathit{h}}_{\mathrm{vv}}$ параметры имеют перекрывающиеся функции распределения вероятностей (PDFs).

Параметры крутящего момента насоса:

helperPlotHistogram(HealthyTheta2, LargeTheta2, SmallTheta2, {'k0','k1','k2'})

Для параметров Крутящего момента отдельная отделимость не очень хороша. Существует все еще некоторое изменение среднего значения и отклонений, которые могут быть использованы обученным классификатором с 3 режимами. Если PDFs показывают хорошее разделение в среднем значении или отклонении, можно спроектировать тесты отношения правдоподобия, чтобы быстро присвоить тестовый набор данных наиболее вероятному режиму. Это показывают затем для параметров крышки насоса.

Пусть:

Рассмотрите доступные наборы параметров как тестовые выборки для предсказания режима. Присвойте предсказанный режим как принадлежащий одному, для которого объединенная PDF имеет самое высокое значение (Гипотеза ${{\rm H}}_0$ выбран ${{\rm H}}_1$ если $\mathit{p}\left({\mathit{H}}_0\right)>\mathit{p}\left({\mathit{H}}_1\right)$). Затем постройте результаты, сравнивающие истинные и предсказанные режимы в матрице беспорядка. Функциональный mvnpdf используемо в вычислениях значение PDF и функции confusionmatrix и heatmap используются для матричной визуализации беспорядка. Смотрите, что помощник функционирует pumpModeLikelihoodTest.

% Pump head confusion matrix
figure
pumpModeLikelihoodTest(HealthyTheta1, LargeTheta1, SmallTheta1)

График беспорядка показывает совершенное разделение между этими тремя режимами, которое не удивляет, учитывая ясное разделение среди гистограмм для ${\mathit{h}}_{\mathrm{nn}}$ параметры.

% Pump torque confusion matrix
pumpModeLikelihoodTest(HealthyTheta2, LargeTheta2, SmallTheta2)

Результаты немного хуже для параметров крутящего момента. Однако, показатель успешности довольно высок (97%) даже при том, что PDFs этих трех режимов, перекрытых значительно. Это вызвано тем, что вычисление значения PDF затронуто оба местоположением (среднее значение), а также амплитуда (отклонение).

B. Классификация мультиклассов режимов отказа Используя древовидное укладывание в мешки

В этом разделе другой метод классификации обсужден, который более подходит, когда классификация среди большего числа режимов требуется. Рассмотрите следующие режимы работы отказа насоса:

Проблема классификации более трудна, потому что существует только три вычисленные параметра, и необходимо различать 7 режимов работы. Таким образом не только необходимо сравнить предполагаемые параметры для каждого режима отказа к здоровому режиму, но также и друг другу - оба, направление (увеличение или сокращение значения) и величина (10%-е изменение по сравнению с 70%-м изменением) изменения параметра должно быть учтено.

Здесь использование классификатора TreeBagger для этой проблемы показывают. Древовидный Мешконасыпатель является ансамблем, изучающим метод, который использует агрегацию начальной загрузки (укладывание в мешки) функций, чтобы создать деревья решений, которые классифицируют маркированные данные. 50 помеченных наборов данных собраны для этих 7 режимов работы. Оцените параметры крышки насоса для каждого набора данных и обучите классификатор с помощью подмножества оценок параметра от каждого режима.

url = 'https://www.mathworks.com/supportfiles/predmaint/fault-diagnosis-of-centrifugal-pumps-using-steady-state-experiments/MultipleFaultsData.mat';
websave('MultipleFaultsData.mat',url);
load MultipleFaultsData
% Compute pump head parameters
HealthyTheta = linearFit(1, HealthyEnsemble);
Fault1Theta  = linearFit(1, Fault1Ensemble);
Fault2Theta  = linearFit(1, Fault2Ensemble);
Fault3Theta  = linearFit(1, Fault3Ensemble);
Fault4Theta  = linearFit(1, Fault4Ensemble);
Fault5Theta  = linearFit(1, Fault5Ensemble);
Fault6Theta  = linearFit(1, Fault6Ensemble);

% Generate labels for each mode of operation
Label = {'Healthy','ClearanceGapWear','ImpellerOutletDeposit',...
    'ImpellerInletDeposit','AbrasiveWear','BrokenBlade','Cavitation'};
VarNames = {'hnn','hnv','hvv','Condition'};
% Assemble results in a table with parameters and corresponding labels 
N = 50; 
T0 = [array2table(HealthyTheta),repmat(Label(1),[N,1])];
T0.Properties.VariableNames = VarNames;
T1 = [array2table(Fault1Theta), repmat(Label(2),[N,1])];
T1.Properties.VariableNames = VarNames;
T2 = [array2table(Fault2Theta), repmat(Label(3),[N,1])];
T2.Properties.VariableNames = VarNames;
T3 = [array2table(Fault3Theta), repmat(Label(4),[N,1])];
T3.Properties.VariableNames = VarNames;
T4 = [array2table(Fault4Theta), repmat(Label(5),[N,1])];
T4.Properties.VariableNames = VarNames;
T5 = [array2table(Fault5Theta), repmat(Label(6),[N,1])];
T5.Properties.VariableNames = VarNames;
T6 = [array2table(Fault6Theta), repmat(Label(7),[N,1])];
T6.Properties.VariableNames = VarNames;

% Stack all data
% Use 30 out of 50 datasets for model creation
TrainingData = [T0(1:30,:);T1(1:30,:);T2(1:30,:);T3(1:30,:);T4(1:30,:);T5(1:30,:);T6(1:30,:)];

% Create an ensemble Mdl of 20 decision trees for predicting the 
% labels using the parameter values
rng(3) % for reproducibility
Mdl = TreeBagger(20, TrainingData, 'Condition',...
   'OOBPrediction','on',...
   'OOBPredictorImportance','on')

Mdl = 
  TreeBagger
Ensemble with 20 bagged decision trees:
                    Training X:              [210x3]
                    Training Y:              [210x1]
                        Method:       classification
                 NumPredictors:                    3
         NumPredictorsToSample:                    2
                   MinLeafSize:                    1
                 InBagFraction:                    1
         SampleWithReplacement:                    1
          ComputeOOBPrediction:                    1
 ComputeOOBPredictorImportance:                    1
                     Proximity:                   []
                    ClassNames:  'AbrasiveWear'   'BrokenBlade'    'Cavitation' 'ClearanceGapWear'       'Healthy' 'ImpellerInletDeposit' 'ImpellerOutletDeposit'

  Properties, Methods

Эффективность модели TreeBagger может быть вычислена путем изучения ее misclassification вероятности для наблюдений из сумки как функция количества деревьев решений.

% Compute out of bag error
figure
plot(oobError(Mdl))
xlabel('Number of trees')
ylabel('Misclassification probability')

Наконец, вычислите эффективность предсказания модели на тестовых выборках, которые никогда не использовались для роста деревьев решений.

ValidationData = [T0(31:50,:);T1(31:50,:);T2(31:50,:);T3(31:50,:);T4(31:50,:);T5(31:50,:);T6(31:50,:)];
PredictedClass = predict(Mdl,ValidationData);
E = zeros(1,7);
% Healthy data misclassification
E(1) = sum(~strcmp(PredictedClass(1:20), Label{1}));
% Clearance gap fault misclassification
E(2) = sum(~strcmp(PredictedClass(21:40), Label{2}));
% Impeller outlet deposit fault misclassification
E(3) = sum(~strcmp(PredictedClass(41:60), Label{3}));
% Impeller inlet deposit fault misclassification
E(4) = sum(~strcmp(PredictedClass(61:80), Label{4}));
% Abrasive wear fault misclassification
E(5) = sum(~strcmp(PredictedClass(81:100), Label{5}));
% Broken blade fault misclassification
E(6) = sum(~strcmp(PredictedClass(101:120), Label{6}));
% Cavitation fault misclassification
E(7) = sum(~strcmp(PredictedClass(121:140), Label{7}));
figure
bar(E/20*100)
xticklabels(Label)
set(gca,'XTickLabelRotation',45)
ylabel('Misclassification (%)')

График показывает, что абразивный износ и поврежденные блейд-отказы неправильно классифицируются для 30% выборок валидации. Более внимательное рассмотрение в предсказанных метках показывает, что в неправильно классифицированных случаях, марки 'AbrasiveWear' и 'BrokenBlade' смешаны друг между другом только. Это предполагает, что признаки для этих категорий отказа не достаточно различимы для этого классификатора.

Сводные данные

Хорошо спроектированная стратегия диагностики отказа может сократить эксплуатационные затраты путем минимизации сервисного времени простоя и заменяющих затрат компонента. Стратегия извлекает выгоду из хорошего знания о динамике операционной машины, которая используется в сочетании с измерениями датчика, чтобы обнаружить и изолировать различные виды отказов.

Этот пример обсудил параметрический подход для обнаружения отказа и изоляции на основе установившихся экспериментов. Этот подход требует тщательного моделирования системной динамики и использования параметров (или преобразования этого) как функции разработки алгоритмов диагностики отказа. Параметры использовались для учебных детекторов аномалии, выполняя тесты отношения правдоподобия и для обучения классификаторов мультикласса.

Как использовать методы классификации в реальном тестировании насосов

Сводные данные рабочего процесса диагностики отказа следуют.

Для подхода диагностики отказа на основе остаточного анализа смотрите пример Обнаружения отказов центробежных насосов с использованием анализа невязок.

Ссылки

Вспомогательные Функции

Линейная подгонка, чтобы накачать параметры.

function varargout = linearFit(Form, Data)
%linearFit Linear least squares solution for Pump Head and Torque parameters.
%
% If Form==0, accept separate inputs and return separate outputs. For one experiment only.
% If Form==1, accept an ensemble and return compact parameter vectors. For several experiments (ensemble).
if Form==0
   w = Data{1};
   Q = Data{2};
   H = Data{3};
   M = Data{4};
   n = length(Q);
   if isscalar(w), w = w*ones(n,1); end
   Q = Q(:); H = H(:); M = M(:);
   Predictor = [w.^2, w.*Q, Q.^2];
   Theta1 = Predictor\H;
   hnn =  Theta1(1);
   hnv = -Theta1(2);
   hvv = -Theta1(3);
   Theta2 = Predictor\M;
   k0 =  Theta2(2);
   k1 = -Theta2(3);
   k2 =  Theta2(1);
   varargout = {hnn, hnv, hvv, k0, k1, k2};
else
   H = cellfun(@(x)x.Head,Data,'uni',0);
   Q = cellfun(@(x)x.Discharge,Data,'uni',0);
   M = cellfun(@(x)x.Torque,Data,'uni',0);
   W = cellfun(@(x)x.Speed,Data,'uni',0);
   N = numel(H);

   Theta1 = zeros(3,N);
   Theta2 = zeros(3,N);
   
   for kexp = 1:N
      Predictor = [W{kexp}.^2, W{kexp}.*Q{kexp}, Q{kexp}.^2];
      X1 = Predictor\H{kexp};
      hnn =  X1(1);
      hnv = -X1(2);
      hvv = -X1(3);
      X2 = Predictor\M{kexp};
      k0 =  X2(2);
      k1 = -X2(3);
      k2 =  X2(1);
      
      Theta1(:,kexp) = [hnn; hnv; hvv];
      Theta2(:,kexp) = [k0; k1; k2];
   end
   varargout = {Theta1', Theta2'};
end
end

Тест вероятности членства.

function pumpModeLikelihoodTest(HealthyTheta, LargeTheta, SmallTheta)
%pumpModeLikelihoodTest Generate predictions based on PDF values and plot confusion matrix.

m1 = mean(HealthyTheta);
c1 = cov(HealthyTheta);
m2 = mean(LargeTheta);
c2 = cov(LargeTheta);
m3 = mean(SmallTheta);
c3 = cov(SmallTheta);

N = size(HealthyTheta,1);

% True classes
% 1: Healthy: group label is 1.
X1t = ones(N,1);
% 2: Large gap: group label is 2.
X2t = 2*ones(N,1);
% 3: Small gap: group label is 3.
X3t = 3*ones(N,1);

% Compute predicted classes as those for which the joint PDF has the maximum value.
X1 = zeros(N,3); 
X2 = zeros(N,3); 
X3 = zeros(N,3); 
for ct = 1:N
   % Membership probability density for healthy parameter sample
   HealthySample  = HealthyTheta(ct,:);
   x1 = mvnpdf(HealthySample, m1, c1);
   x2 = mvnpdf(HealthySample, m2, c2);
   x3 = mvnpdf(HealthySample, m3, c3);
   X1(ct,:) = [x1 x2 x3];
   
   % Membership probability density for large gap pump parameter
   LargeSample  = LargeTheta(ct,:);
   x1 = mvnpdf(LargeSample, m1, c1);
   x2 = mvnpdf(LargeSample, m2, c2);
   x3 = mvnpdf(LargeSample, m3, c3);
   X2(ct,:) = [x1 x2 x3];
   
   % Membership probability density for small gap pump parameter
   SmallSample  = SmallTheta(ct,:);
   x1 = mvnpdf(SmallSample, m1, c1);
   x2 = mvnpdf(SmallSample, m2, c2);
   x3 = mvnpdf(SmallSample, m3, c3);
   X3(ct,:) = [x1 x2 x3];
end

[~,PredictedGroup] = max([X1;X2;X3],[],2);
TrueGroup = [X1t; X2t; X3t];
C = confusionmat(TrueGroup,PredictedGroup);
heatmap(C, ...
    'YLabel', 'Actual condition', ...
    'YDisplayLabels', {'Healthy','Large Gap','Small Gap'}, ...
    'XLabel', 'Predicted condition', ...
    'XDisplayLabels', {'Healthy','Large Gap','Small Gap'}, ...
    'ColorbarVisible','off');
end