Этот пример показывает основанный на модели подход для обнаружения и диагноза различных типов отказов, которые происходят в системе накачки. Пример следует за центробежным анализом насоса, представленным в книге Приложений Диагностики отказа Рольфа Изерманна [1].
Насосы являются существенным оборудованием во многих отраслях промышленности включая степень и химический, минерал и горная промышленность, производство, нагревание, кондиционирование воздуха и охлаждение. Центробежные насосы используются, чтобы транспортировать жидкости преобразованием вращательной кинетической энергии к гидродинамической энергии потока жидкости. Вращательная энергия обычно прибывает из двигателя внутреннего сгорания или электродвигателя. Жидкость вводит рабочее колесо насоса вперед или близко к вращающейся оси и ускоряется рабочим колесом, текущий радиально исходящий в диффузор.
Насосы испытывают повреждение своих гидравлических или механических компонентов. Наиболее часто дающие сбой компоненты являются скользящими кольцевыми изоляциями и шарикоподшипниками, несмотря на то, что отказы к другим компонентам включая ведущий двигатель, лопатки рабочего колеса и скользящие подшипники также весьма распространены. Таблица ниже приводит наиболее распространенные типы отказов.
Кавитация: Разработка пузырей пара в жидкости, если статическое давление падает ниже давления пара. Пузыри выходят из строя резко ведущий к повреждению в шкивах ленточной пилы.
Газ в жидкости: перепад давления приводит к растворенному газу в жидкости. Разделение газа и жидкости и более низких главных результатов.
Пробный прогон: Недостающая жидкость приводит к отсутствию охлаждения и перегрева подшипника. Важный для стартовой фазы.
Эрозия: Механическое повреждение стенок из-за трудных частиц или кавитации
Коррозия: Повреждение агрессивными жидкостями
Износ подшипников: Механическое повреждение через усталость и металлическое трение, генерацию точечной коррозии и слез
Включение вспомогательных буровых скважин: Приводит к перегрузке/повреждению упорных подшипников
Включение скользящих кольцевых изоляций: Приводит к более высокому трению и меньшему КПД
Увеличение изоляций разделения: Приводит к снижению эффективности
Депозиты: Депозиты органического материала или посредством химических реакций во входе ротора или выхода уменьшают КПД и увеличивают температуру.
Колебания: неустойчивость Ротора через повреждение или депозиты в роторе. Может нанести ущерб подшипника.
Доступные датчики
Следующие сигналы обычно измеряются:
Перепад давлений между входом и выходом
Скорость вращения
Крутящий момент двигателя и накачайте крутящий момент
Жидкий выброс (поток) уровень при выходе насоса
Ведущий моторный ток, напряжение, температура (не рассмотренный здесь)
Температура жидкости, отложения (не рассмотренный здесь)
Крутящий моментпримененный ротор радиального центробежного насоса приводит к скорости вращения и передает увеличение импульса жидкости насоса от входа ротора меньшего радиуса к выходу ротора большего радиуса. Турбинные уравнения Эйлера дают к отношению между перепадом давления , скорость и жидкость разряжает уровень (скорость потока жидкости) :
где теоретическое (идеал; без потерь) крышка насоса, измеренная в метрах и , коэффициенты пропорциональности. При составлении конечного числа лопаток рабочего колеса, потерь на трение и потерь удара из-за нетангенциального потока, действительной крышкой насоса дают:
где , и коэффициенты пропорциональности должны быть обработаны как параметры модели. Соответствующий крутящий момент насоса:
Механические детали двигателя и насоса заставляют скорость увеличиваться, когда крутящий момент применяется согласно:
где отношение инерции двигателя и насоса, и фрикционный крутящий момент, состоящий из трения Кулона и вязкое трение согласно:
Насос соединяется с трубопроводной системой, которая транспортирует жидкость от более низкого резервуара для хранения до верхнего. Выражения уравнения баланса импульса:
где коэффициент сопротивления трубопровода, с длиной трубопровода и площадь поперечного сечения , и высота устройства хранения данных по насосу. Параметры модели или известны от физики или может быть оценен путем подбора кривой измеренным сигналам датчика к входным параметрам/выходным параметрам модели. Тип используемой модели может зависеть от условий работы, под которыми запущен насос. Например, полная нелинейная модель системы трубопровода насоса не может требоваться, если насос всегда запускается на постоянной угловой скорости.
Отказы могут быть обнаружены путем исследования определенных функций, извлеченных из измерений и сравнения их с известными порогами приемлемого поведения. Обнаружительная способность и изонеустойчивость различных отказов зависят от природы эксперимента и доступности измерений. Например, анализ постоянной скорости с измерениями давления только может обнаружить отказы, вызывающие большие скачки давления. Кроме того, это не может надежно оценить причину отказа. Однако многоскоростной эксперимент с измерениями перепада давления, крутящего момента двигателя и скорости потока жидкости может обнаружить и изолировать много источников отказов, таких как те, которые происходят из газовых корпусов, пробного прогона, больших депозитов, моторные дефекты и т.д.
Основанный на модели подход использует следующие методы:
Оценка параметра: Используя измерения от здоровой (номинальной) работы машины, оцениваются параметры модели, и их неопределенность определяется количественно. Измерения системы тестирования затем используются, чтобы повторно оценить значения параметров, и получившиеся оценки сравнены со своей номинальной стоимостью. Этот метод является основной темой этого примера.
Генерация остатка: модель обучена как прежде для здоровой машины. Выход модели сравнен с измеренными наблюдениями от системы тестирования, и вычисляется остаточный сигнал. Этот сигнал анализируется для его величины, отклонения и других свойств обнаружить отказы. Большое количество остатков может проектироваться и использоваться, чтобы отличить другие источники отказов. Этот метод обсужден в примере Обнаружения отказов центробежных насосов с использованием анализа невязок.
Установившаяся практика для калибровки насоса и контроля должна запустить его на постоянной скорости и записать статический главный и жидкий уровень выброса насоса. Меняя положение клапана в трубопроводной системе, жидкий объем выброса (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
Подобные изменения видны в характеристиках потока крутящего момента и для других типов отказа.
Для автоматизации диагностики отказа вы поворачиваете наблюдаемые изменения в количественных данных. Надежный способ сделать так состоит в том, чтобы соответствовать параметрированной кривой к данным о характеристике главного потока, отображенным на графике выше. Используя управляющие уравнения динамики трубопровода насоса и использование упрощенного отношения крутящего момента, получены следующие уравнения:
параметры должны быть оценены. Если вы измеряетесь и , параметр может быть оценен линейным методом наименьших квадратов. Эти параметры являются функциями, которые могут быть использованы, чтобы разработать обнаружение отказа и алгоритм диагноза.
Вычислите и постройте значения параметров, оцененные для вышеупомянутых 3 кривых. Используйте измеренные значения и как данные и как номинальная скорость насоса.
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
Таблицы показывают это и значения уменьшают, когда разрыв разрешения является большим, в то время как они больше, чем номинальная стоимость для маленького разрешения. С другой стороны, и значения увеличиваются для большого разрыва разрешения и уменьшения для маленького разрыва. Зависимость и на разрешении разрыв менее ясен.
Предварительный анализ показал, как изменения параметра могут указать на отказ. Однако даже для здоровых насосов существуют изменения измерений вследствие шума измерения, жидкого загрязнения и изменений вязкости и характеристик крутящего момента промаха двигателя, запускающего насос. Эти изменения измерения вводят неопределенность в оценках параметра.
Соберите 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 стандартным отклонениям (). Смотрите, что помощник функционирует 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 с помощью здоровых данных о параметре насоса. С тех пор нет никаких используемых меток отказа, не обрабатывают все выборки как прибывающий из того же (здорового) класса. Начиная с изменений в параметрах и являются самыми показательными из потенциальных отказов, используют только эти параметры для обучения классификатор 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 путями среди следующих режимов:
Режим 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'})
Гистограмма показывает это предложения хорошая отделимость среди этих трех режимов, но параметры имеют перекрывающиеся функции распределения вероятностей (PDFs).
Параметры крутящего момента насоса:
helperPlotHistogram(HealthyTheta2, LargeTheta2, SmallTheta2, {'k0','k1','k2'})
Для параметров Крутящего момента отдельная отделимость не очень хороша. Существует все еще некоторое изменение среднего значения и отклонений, которые могут быть использованы обученным классификатором с 3 режимами. Если PDFs показывают хорошее разделение в среднем значении или отклонении, можно спроектировать тесты отношения правдоподобия, чтобы быстро присвоить тестовый набор данных наиболее вероятному режиму. Это показывают затем для параметров крышки насоса.
Пусть:
: гипотеза, что главные параметры принадлежат здоровому режиму насоса
: гипотеза, что главные параметры принадлежат насосу с большим разрывом разрешения
: гипотеза, что главные параметры принадлежат насосу с маленьким разрывом разрешения
Рассмотрите доступные наборы параметров как тестовые выборки для предсказания режима. Присвойте предсказанный режим как принадлежащий одному, для которого объединенная PDF имеет самое высокое значение (Гипотеза выбран если ). Затем постройте результаты, сравнивающие истинные и предсказанные режимы в матрице беспорядка. Функциональный mvnpdf
используемо в вычислениях значение PDF и функции confusionmatrix
и heatmap
используются для матричной визуализации беспорядка. Смотрите, что помощник функционирует pumpModeLikelihoodTest
.
% Pump head confusion matrix
figure
pumpModeLikelihoodTest(HealthyTheta1, LargeTheta1, SmallTheta1)
График беспорядка показывает совершенное разделение между этими тремя режимами, которое не удивляет, учитывая ясное разделение среди гистограмм для параметры.
% 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' смешаны друг между другом только. Это предполагает, что признаки для этих категорий отказа не достаточно различимы для этого классификатора.
Хорошо спроектированная стратегия диагностики отказа может сократить эксплуатационные затраты путем минимизации сервисного времени простоя и заменяющих затрат компонента. Стратегия извлекает выгоду из хорошего знания о динамике операционной машины, которая используется в сочетании с измерениями датчика, чтобы обнаружить и изолировать различные виды отказов.
Этот пример обсудил параметрический подход для обнаружения отказа и изоляции на основе установившихся экспериментов. Этот подход требует тщательного моделирования системной динамики и использования параметров (или преобразования этого) как функции разработки алгоритмов диагностики отказа. Параметры использовались для учебных детекторов аномалии, выполняя тесты отношения правдоподобия и для обучения классификаторов мультикласса.
Как использовать методы классификации в реальном тестировании насосов
Сводные данные рабочего процесса диагностики отказа следуют.
Запустите тестовый насос на его номинальной скорости. Поверните выпускной клапан к различным настройкам, чтобы управлять скоростью потока жидкости. Для каждого положения клапана запишите скорость насоса, скорость потока жидкости, перепады давления и крутящий момент.
Оцените параметры для крышки насоса и накачайте характеристику крутящего момента (устойчивое состояние) уравнения.
Если неопределенность/шум является низкой, и оценки параметра надежны, предполагаемые параметры могут быть непосредственно по сравнению с их номинальной стоимостью. Их относительные величины указали бы на природу отказа.
В общей шумной ситуации используйте методы обнаружения аномалии, чтобы сначала проверять, существует ли отказ, существующий в системе вообще. Это может быть сделано очень быстро путем сравнения предполагаемых значений параметров против среднего значения и значений ковариации, полученных из исторической базы данных здоровых насосов.
Если отказ обозначается, используйте методы классификации отказов (такие как тесты отношения правдоподобия или выход классификатора), чтобы изолировать наиболее вероятную причину (причины). Выбор метода классификации зависел бы от доступных данных датчика, их надежности, серьезности отказа и доступности исторической информации относительно режимов отказа.
Для подхода диагностики отказа на основе остаточного анализа смотрите пример Обнаружения отказов центробежных насосов с использованием анализа невязок.
Изерманн, Рольф, приложения диагностики отказа. Основанный на модели мониторинг состояния: приводы, диски, машинное оборудование, объекты, датчики, и отказоустойчивая система, выпуск 1, Springer-Verlag Берлин Гейдельберг, 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