В этом примере показан основанный на модели подход к обнаружению и диагностике различных типов неисправностей, возникающих в насосной системе. В качестве примера приведен анализ центробежных насосов, представленный в книге Rolf Isermann «Приложения диагностики неисправностей» [1].
Насосы являются основным оборудованием во многих отраслях промышленности, включая энергетику и химическую промышленность, добычу полезных ископаемых и полезных ископаемых, производство, отопление, кондиционирование воздуха и охлаждение. Центробежные насосы используются для транспортировки текучих сред путем преобразования кинетической энергии вращения в гидродинамическую энергию потока текучей среды. Энергия вращения обычно поступает от двигателя внутреннего сгорания или электродвигателя. Текучая среда поступает в рабочее колесо насоса вдоль или вблизи оси вращения и ускоряется рабочим колесом, проходя радиально наружу в диффузор.

Насосы имеют повреждения гидравлических или механических компонентов. Наиболее часто неисправными компонентами являются кольцевые уплотнения скольжения и шарикоподшипники, хотя отказы других компонентов, включая приводной двигатель, лопасти рабочего колеса и подшипники скольжения, также не являются редкостью. В таблице ниже перечислены наиболее распространенные типы отказов.
Кавитация: Образование пузырьков пара внутри жидкости, если статическое давление падает ниже давления пара. Пузыри резко разрушаются, что приводит к повреждению лопастных колес.
Газ в жидкости: падение давления приводит к растворению газа в жидкости. В результате происходит разделение газа и жидкости и нижней части головки.
Сухой прогон: Отсутствие жидкости приводит к отсутствию охлаждения и перегреву подшипника. Важно для начальной фазы.
Эрозия: механическое повреждение стенок из-за твердых частиц или кавитации
Коррозия: повреждение агрессивными жидкостями
Износ подшипника: Механическое повреждение в результате усталости и трения металла, образование пятен и разрывов
Закупорка отверстий рельефного отверстия: приводит к перегрузке/повреждению осевых подшипников
Закупорка уплотнений скользящего кольца: приводит к более высокому трению и меньшей эффективности
Увеличение количества разделенных уплотнений: приводит к потере эффективности
Отложения: отложения органического материала или химические реакции на входе или выходе ротора снижают эффективность и повышают температуру.
Колебания: Дисбаланс ротора через повреждение или отложения на роторе. Может привести к повреждению подшипника.
Доступные датчики
Обычно измеряют следующие сигналы:
Перепад давления между входом и выходом
Частота вращения
Крутящий момент двигателя и крутящий момент насоса
Расход (расход) жидкости на выходе из насоса
Ток, напряжение, температура приводного двигателя (здесь не учитываются)
Температура жидкости, осадки (здесь не рассматриваются)
Сигнал, приложенный к ротору радиального центробежного насоса, приводит к скорости λ вращения и передает увеличение импульса насосной жидкости от входа ротора меньшего радиуса к выходу ротора большего радиуса. Уравнения турбины Эйлера вырабатывают зависимость между перепадом давления , скоростью и скоростью выпуска жидкости (расход) :
h1start2-h2startQ
где Δpαg - теоретический (идеал; без потерь) напор насоса, измеренный в метрах h1h2 - постоянные пропорциональности. При учете конечного числа лопаток рабочего колеса, потерь на трение и ударных потерь из-за нетангенциального потока реальный напор насоса задается:
hnnω2-hnvωQ-hvvQ2
где , и - константы пропорциональности, которые должны рассматриваться как параметры модели. Соответствующий крутящий момент насоса:
hnnωQ-hnvQ2-hvvQ3ω)
Механические части двигателя и насоса вызывают увеличение скорости при приложении крутящего момента в соответствии с:
) -Mf (t)
где - отношение инерции двигателя и насоса, а - фрикционный крутящий момент, состоящий из фрикционной Кулона и вязкой фрикционной ) согласно:
Mf1λ (t)
Насос соединен с трубопроводной системой, которая транспортирует жидкость из нижнего резервуара хранения в верхний. Уравнение баланса импульса дает:
+ Hstatic
где - коэффициент сопротивления трубы, 1gA с трубы l и поперечного сечения А, Hstatic - высота хранилища над насосом. модели
Дефекты могут быть обнаружены путем изучения определенных признаков, извлеченных из измерений, и сравнения их с известными порогами приемлемого поведения. Обнаруживаемость и изолируемость различных неисправностей зависит от характера эксперимента и наличия измерений. Например, анализ постоянной скорости с измерениями давления позволяет обнаружить только неисправности, вызывающие большие изменения давления. Кроме того, он не может достоверно оценить причину сбоя. Однако многоскоростной эксперимент с измерениями разности давлений, крутящего момента двигателя и расхода может обнаружить и изолировать многие источники неисправностей, такие как те, которые происходят из газовых оболочек, сухого хода, больших отложений, дефектов двигателя и т.д.
Подход, основанный на модели, использует следующие методы:
Оценка параметров: С помощью измерений от здоровой (номинальной) работы машины оценивают параметры модели и количественно оценивают их неопределенность. Затем измерения тестовой системы используются для повторной оценки значений параметров, и полученные оценки сравниваются с их номинальными значениями. Эта техника является основной темой этого примера.
Образование остатков: Модель обучена, как и прежде, здоровой машине. Выходные данные модели сравнивают с измеренными наблюдениями из тестовой системы и вычисляют остаточный сигнал. Этот сигнал анализируется на его величину, дисперсию и другие свойства для обнаружения неисправностей. Большое количество остатков может быть разработано и использовано для различения различных источников разломов. Этот метод рассматривается в примере диагностики отказов центробежных насосов с использованием остаточного анализа.
Обычной практикой калибровки и контроля насоса является запуск его с постоянной скоростью и регистрация статического напора насоса и скорости нагнетания жидкости. Изменяя положение клапана в трубопроводной системе, регулируется объем выпуска жидкости (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

Аналогичные изменения могут наблюдаться в характеристиках крутящего момента-потока и для других типов неисправностей.
Для автоматизации диагностики неисправностей следует обратить наблюдаемые изменения в количественную информацию. Надежный способ сделать это - подогнать параметризованную кривую к показанным выше характеристическим данным напора. С помощью уравнений, регулирующих динамику насоса и трубы, и упрощенной зависимости крутящего момента получают следующие уравнения:
k1, k2 - параметры, подлежащие оценке. Если Q, H и Mp, параметр может быть оценен по линейным наименьшим квадратам. Эти параметры являются функциями, которые можно использовать для разработки алгоритма обнаружения и диагностики неисправностей.
Вычислите и постройте график значений параметров, оцененных для вышеуказанных 3 кривых. В качестве данных используйте измеренные значения Q, H и MP, а в качестве номинальной частоты вращения насоса - λ = 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
Таблицы показывают, что значения hnn и k0 уменьшаются, когда зазор велик, в то время как они больше номинальных значений для малого зазора. С другой стороны, значения hvv и k2 увеличиваются для большого зазора и уменьшаются для малого зазора. Зависимость и от зазора менее ясна.
Предварительный анализ показал, как изменения параметров могут указывать на неисправность. Однако даже для здоровых насосов существуют различия в измерениях из-за шума измерений, загрязнения текучей среды и изменения вязкости, а также характеристик скользящего крутящего момента двигателя, работающего на насосе. Эти вариации измерений вносят неопределенность в оценки параметров.
Соберите 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 стандартным отклонениям (≈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
Test Ensemble содержит набор измерений скорости насоса, крутящего момента, напора и расхода в различных положениях клапана. Все измерения содержат провал зазора разной величины.
% 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 с использованием данных параметров исправного насоса. Поскольку метки неисправностей не используются, следует рассматривать все образцы как поступающие из одного и того же (здорового) класса. Поскольку изменения параметров и наиболее характерны для потенциальных неисправностей, используйте только эти параметры для обучения классификатора 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')

Подобный анализ может быть выполнен для обнаружения других видов неисправностей, таких как износ или отложения на выходе из рабочего колеса, как обсуждалось далее в контексте локализации неисправностей.
Если информация о типе (типах) неисправности в тестовой системе доступна, она может использоваться для создания алгоритмов, которые не только обнаруживают неисправности, но и указывают их тип .
А. Разграничение дефектов зазоров по тестам отношения правдоподобия
Изменения в зазоре могут быть разделены на два типа - меньший, чем ожидалось, зазор (желтая линия на графике характеристик напора насоса) и больший, чем ожидаемый зазор (красная линия). Наборы данных испытаний нагрузочного насоса, содержащие дефекты зазоров, в которых заранее известна природа неисправности (большая или малая). Эти метки неисправностей используются для выполнения трехсторонней классификации в следующих режимах:
Режим 1: Нормальный зазор (здоровое поведение)
Режим 2: большой зазор
Режим 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'})
Гистограмма показывает, что обеспечивает хорошую разделяемость между тремя режимами, но параметры hnv, hvv имеют перекрывающиеся функции распределения вероятности (PDF ).
Параметры крутящего момента насоса:
helperPlotHistogram(HealthyTheta2, LargeTheta2, SmallTheta2, {'k0','k1','k2'})
Для параметров крутящего момента индивидуальная разделяемость не очень хороша. Есть еще некоторые различия в среднем и дисперсии, которые могут быть использованы обученным 3-модовым классификатором. Если PDF-файлы показывают хорошее разделение в среднем или расхождении, можно разработать тесты отношения правдоподобия, чтобы быстро назначить набор тестовых данных наиболее вероятному режиму. Это показано ниже для параметров напора насоса.
Давайте:
: гипотеза о том, что параметры напора относятся к режиму здорового насоса
: гипотеза о том, что параметры напора относятся к насосу с большим зазором
: гипотеза о том, что параметры напора относятся к насосу с небольшим зазором
Рассмотрим доступные наборы параметров в качестве тестовых выборок для прогнозирования режима. Назначьте прогнозируемый режим как относящийся к режиму, для которого совместное PDF имеет наибольшее значение (гипотеза выбирается если (H1)). Затем постройте график результатов сравнения истинного и предсказанного режимов в матрице путаницы. Функцияmvnpdf используется для вычисления значения PDF и функций confusionmatrix и heatmap используются для визуализации матрицы путаницы. См. вспомогательную функцию pumpModeLikelihoodTest.
% Pump head confusion matrix
figure
pumpModeLikelihoodTest(HealthyTheta1, LargeTheta1, SmallTheta1)
График путаницы показывает идеальное разделение между тремя режимами, что неудивительно, учитывая четкое разделение между гистограммами для параметров hnn.
% 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» и «BreadBlade» смешиваются только между собой. Это говорит о том, что симптомы для этих категорий отказов недостаточно различимы для данного классификатора.
Продуманная стратегия диагностики неисправностей позволяет сократить эксплуатационные расходы за счет минимизации простоев и затрат на замену компонентов. Стратегия использует хорошие знания о динамике рабочей машины, которые используются в сочетании с измерениями датчиков для обнаружения и изоляции различных типов неисправностей.
В этом примере рассматривается параметрический подход к обнаружению и изоляции неисправностей на основе стационарных экспериментов. Этот подход требует тщательного моделирования динамики системы и использования параметров (или их преобразований) в качестве признаков для разработки алгоритмов диагностики неисправностей. Параметры использовались для обучения детекторам аномалий, выполнения тестов отношения правдоподобия и для обучения мультиклассных классификаторов.
Как использовать методы классификации в реальных испытаниях насосов
Ниже приводится краткое описание процесса диагностики неисправностей.
Запустить тестовый насос на номинальной скорости. Поверните выпускной клапан к различным настройкам для управления расходом. Для каждого положения клапана запишите скорость насоса, расход, перепад давления и крутящий момент.
Оценка параметров для уравнений характеристик напора насоса и крутящего момента насоса (установившийся режим).
Если неопределенность/шум являются низкими и оценки параметров являются надежными, оцененные параметры могут быть непосредственно сопоставлены с их номинальными значениями. Их относительные величины указывали бы на характер разлома.
В общей шумной ситуации используйте методы обнаружения аномалий, чтобы сначала проверить наличие неисправности в системе. Это может быть сделано очень быстро путем сравнения оцененных значений параметров со средними значениями и значениями ковариации, полученными из исторической базы данных здоровых насосов.
Если неисправность указана, используйте методы классификации неисправностей (такие как тесты отношения правдоподобия или выходные данные классификатора) для выявления наиболее вероятной причины (причин). Выбор метода классификации будет зависеть от имеющихся данных датчиков, их надежности, серьезности неисправности и наличия исторической информации о режимах неисправности.
Для получения информации о подходе к диагностике неисправностей, основанном на остаточном анализе, см. пример диагностики неисправностей центробежных насосов с использованием остаточного анализа.
Изерманн, Рольф, приложения для диагностики неисправностей. Мониторинг состояния на основе модели: исполнительные механизмы, приводы, машины, установки, датчики и отказоустойчивая система, издание 1, Springer-Verlag Berlin Heidelberg, 2011.
Линейная посадка для параметров насоса.
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