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

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

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

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

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

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

Крутящий момент$$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.

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

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 комплектов измерений от насоса, работающего в безотказном условии, запустив его со скоростью 2900 об/мин для 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 вычисляет расстояние между тестовыми выборками Махаланобиса и распределением набора ссылки выборки (исправный насос набора параметров здесь):

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')

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

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

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

А. Различение разрывов зазоров по тестам коэффициента правдоподобия

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

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}}$ параметры имеют перекрывающиеся функции распределения вероятностей (PDF).

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

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

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

Давайте:

Рассмотрите доступные наборы параметров как тестовые выборки для предсказания режима. Присвойте предсказанный режим как принадлежащий тому, для которого joint 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%), хотя PDF трех режимов значительно перекрывались. Это связано с тем, что на вычисление значения PDF влияет как местоположение (среднее), так и амплитуда ( отклонение).

B. Многоклассовая классификация режимов отказа с использованием древовидной упаковки

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

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

Здесь показано использование классификатора TreeBagger для этой задачи. Tree Bagger - это метод обучения ансамблей, который использует агрегирование (упаковку в мешки) функций для создания деревьев решений, которые классифицируют маркированные данные. 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 может быть вычислена путем изучения ее вероятности неправильной классификации для наблюдений вне мешка как функции от количества деревьев решений.

% 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