В этом примере показано, как использовать подход моделирования на основе данных для обнаружения неисправностей.
Раннее обнаружение и изоляция аномалий в работе машины может помочь снизить аварийность, сократить время простоя и, таким образом, снизить эксплуатационные расходы. Этот подход включает в себя обработку измерений в режиме реального времени от работы системы, чтобы отметить любое неожиданное поведение, которое указывает на вновь возникшую неисправность.
В этом примере рассматриваются следующие аспекты диагностики неисправностей:
Определение аномального поведения системы с помощью остаточного анализа
Обнаружение изношенности по строительным моделям поврежденной системы
Изменения системы слежения с помощью оперативной адаптации параметров модели
В модельном подходе к обнаружению динамическая модель соответствующей системы сначала строится с использованием измеренных входных и выходных данных. Хорошая модель способна точно предсказать реакцию системы для определенного будущего временного горизонта. Когда предсказание не является хорошим, остатки могут быть большими и могут содержать корреляции. Эти аспекты используются для обнаружения случаев отказа.
Рассмотрим здание, подверженное ударам и вибрациям. Источником вибраций могут быть различные типы стимулов в зависимости от системы, такие как порывы ветра, контакт с работающими двигателями и турбинами или наземные колебания. Воздействия являются результатом испытаний системы на импульсный удар, которые добавляются для достаточного возбуждения системы. Модель Simulink pdmMechanicalSystem.slx является простым примером такой структуры. Возбуждение происходит от периодических ударов, а также от колебаний земли, моделируемых фильтрованным белым шумом. Выходной сигнал системы собирается датчиком, который подвергается измерительному шуму. Модель способна моделировать различные сценарии, включающие структуру в здоровом или поврежденном состоянии.
sysA = 'pdmMechanicalSystem'; open_system(sysA) % Set the model in the healthy mode of operation set_param([sysA,'/Mechanical System'],'LabelModeActiveChoice','Normal') % Simulate the system and log the response data sim(sysA) ynormal = logsout.getElement('y').Values;
Входной сигнал не измерялся; все, что мы записали, это ответ ynormal. Поэтому мы строим динамическую модель системы с использованием методов «слепой идентификации». В частности, мы строим ARMA модель записанного сигнала как представление системы. Этот подход работает, когда предполагается, что входной сигнал является (отфильтрованным) белым шумом. Так как данные подвержены периодическим ударам, мы разделили данные на несколько частей каждый, начиная с возникновения удара. Таким образом, каждый сегмент данных содержит отклик на один удар плюс случайные возбуждения - ситуацию, которая может быть зафиксирована с помощью модели временных рядов, где эффект удара приписывается подходящим начальным условиям.
Ts = 1/256; % data sample time nr = 10; % number of bumps in the signal N = 512; % length of data between bumps znormal = cell(nr,1); for ct = 1:nr ysegment = ynormal.Data((ct-1)*N+(1:500)); z = iddata(ysegment,[],Ts); znormal{ct} = z; % each segment has only one bump end plot(znormal{:}) % plot a sampling of the recorded segments title('Measured Response Segments')
Разбейте данные на части оценки и проверки.
ze = merge(znormal{1:5});
zv = merge(znormal{6:10});
Оценка модели временного ряда 7-го порядка в форме state-space с помощью ssest() команда. Порядок модели был выбран перекрестной проверкой (проверка соответствия данным проверки) и остаточным анализом (проверка того, что остатки не коррелированы).
nx = 7; model = ssest(ze, nx, 'form', 'canonical', 'Ts', Ts); present(model) % view model equations with parameter uncertainty
model =
Discrete-time identified state-space model:
x(t+Ts) = A x(t) + K e(t)
y(t) = C x(t) + e(t)
A =
x1 x2 x3
x1 0 1 0
x2 0 0 1
x3 0 0 0
x4 0 0 0
x5 0 0 0
x6 0 0 0
x7 0.5548 +/- 0.04606 -2.713 +/- 0.2198 5.885 +/- 0.4495
x4 x5 x6
x1 0 0 0
x2 0 0 0
x3 1 0 0
x4 0 1 0
x5 0 0 1
x6 0 0 0
x7 -8.27 +/- 0.5121 9.234 +/- 0.3513 -7.956 +/- 0.1408
x7
x1 0
x2 0
x3 0
x4 0
x5 0
x6 1
x7 4.263 +/- 0.02599
C =
x1 x2 x3 x4 x5 x6 x7
y1 1 0 0 0 0 0 0
K =
y1
x1 1.025 +/- 0.01401
x2 1.444 +/- 0.0131
x3 1.907 +/- 0.01271
x4 2.385 +/- 0.01203
x5 2.857 +/- 0.01456
x6 3.26 +/- 0.0222
x7 3.552 +/- 0.0336
Sample time: 0.0039062 seconds
Parameterization:
CANONICAL form with indices: 7.
Disturbance component: estimate
Number of free coefficients: 14
Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Termination condition: Near (local) minimum, (norm(g) < tol)..
Number of iterations: 7, Number of function evaluations: 15
Estimated using SSEST on time domain data "ze".
Fit to estimation data: [99.07 99.04 99.15 99.05 99.04]% (prediction focus)
FPE: 0.6242, MSE: [0.5974 0.6531 0.5991 0.5871 0.6496]
More information in model's "Report" property.
Отображение модели показывает относительно небольшую неопределенность в оценках параметров. Мы можем подтвердить надежность, вычисляя 1-sd (99,73%) достоверность, связанную с оценочным спектром измеряемого сигнала.
h = spectrumplot(model); showConfidence(h, 3)
Доверительная область невелика, хотя существует около 30% неопределенности в ответе на более низких частотах. Следующим шагом в проверке является просмотр того, насколько хорошо модель предсказывает ответы в наборе данных проверки. zv. Мы используем горизонт прогнозирования на 25 шагов вперед.
compare(zv, model, 25) % Validation against one dataset
График показывает, что модель способна предсказать отклик в первом эксперименте набора данных проверки 25 временных шагов (= 0,1 сек) в будущем с точностью > 85%. Для просмотра соответствия другим экспериментам в наборе данных используйте контекстное меню осей печати, открывающееся щелчком правой кнопкой мыши.
Последним шагом при проверке модели является анализ генерируемых ею остатков. Для хорошей модели эти остатки должны быть белыми, то есть показывать статистически незначительные корреляции для ненулевых лагов:
resid(model, zv)
Остатки в основном некоррелированы при ненулевых лагах. После получения модели нормального поведения мы переходим к изучению того, как модель может быть использована для обнаружения неисправностей.
Обнаружение неисправности - это маркировка нежелательных или неожиданных изменений в наблюдениях системы. Неисправность вызывает изменения в динамике системы вследствие постепенного износа или внезапных изменений, вызванных отказом датчика или поломкой деталей. При появлении неисправности модель, полученная при нормальных рабочих условиях, не может предсказать наблюдаемые реакции. Это приводит к увеличению разницы между измеренным и прогнозируемым откликом (остатками). Такие отклонения обычно отмечаются большой квадратной суммой остатков или наличием корреляций.
Поместите модель Simulink в вариант поврежденной системы и смоделируйте. Мы используем один удар в качестве входного, так как остаточный тест требует белого ввода с, возможно, переходным из-за начальных условий.
set_param([sysA,'/Mechanical System'],'LabelModeActiveChoice','DamagedSystem'); set_param([sysA,'/Pulse'],'Period','5120') % to force only one bump sim(sysA) y = logsout.getElement('y').Values;
resid(model, y.Data) set_param([sysA,'/Pulse'],'Period','512') % restore original
Остатки теперь больше и показывают корреляции при ненулевых лагах. Это основная идея обнаружения неисправностей - создание остаточной метрики и наблюдение за тем, как она изменяется с каждым новым набором измерений. Здесь используется простой остаток, основанный на одношаговой ошибке прогнозирования. На практике генерируются более совершенные остатки, адаптированные к потребностям приложений.
Более детальный подход к обнаружению неисправностей заключается также в определении модели неисправного (поврежденного) состояния системы. Затем мы можем проанализировать, какая модель с большей вероятностью объяснит измерения в реальном времени из системы. Эта компоновка может быть обобщена на модели для различных типов неисправностей и, таким образом, используется не только для обнаружения неисправности, но и для определения того, какая из них («изоляция»). В этом примере используется следующий подход:
Мы собираем данные с системой, работающей в нормальном (здоровом) и известном состоянии, вызванном износом.
Мы определяем динамическую модель, представляющую поведение в каждом состоянии.
Мы используем подход кластеризации данных, чтобы провести четкое различие между этими состояниями.
Для обнаружения неисправностей мы собираем данные с запущенной машины и определяем модель ее поведения. Затем мы прогнозируем, в каком состоянии (нормальном или поврежденном) наиболее вероятно объяснить наблюдаемое поведение.
Мы уже смоделировали систему в режиме нормальной эксплуатации. Теперь мы моделируем модель pdmMechanicalSystem в режиме «конец жизни». Это сценарий, при котором система уже ухудшилась до конечного состояния допустимой эксплуатации.
set_param([sysA,'/Mechanical System'],'LabelModeActiveChoice','DamagedSystem'); sim(sysA) y = logsout.getElement('y').Values; zfault = cell(nr,1); for ct = 1:nr z = iddata(y.Data((ct-1)*N+(1:500)),[],Ts); zfault{ct} = z; end
Теперь мы создаем набор моделей, по одной для каждого сегмента данных. Как и прежде, мы строим модели временных рядов 7-го порядка в форме state-space. Выключите вычисление ковариации для скорости.
mNormal = cell(nr,1); mFault = cell(nr, 1); nx = order(model); opt = ssestOptions('EstimateCovariance',0); for ct = 1:nr mNormal{ct} = ssest(znormal{ct}, nx, 'form', 'canonical', 'Ts', Ts, opt); mFault{ct} = ssest(zfault{ct}, nx, 'form', 'canonical', 'Ts', Ts, opt); end
Убедитесь, что модели mFault являются хорошим представлением неисправного режима работы:
compare(merge(zfault{:}), mFault{:}, 25)
Нормальные и неисправные оценочные спектры показаны ниже.
Color1 = 'k'; Color2 = 'r'; ModelSet1 = cat(2,mNormal,repmat({Color1},[nr, 1]))'; ModelSet2 = cat(2,mFault,repmat({Color2},[nr, 1]))'; spectrum(ModelSet1{:},ModelSet2{:}) axis([1 1000 -45 40]) title('Output Spectra (black: normal, red: faulty)')
График спектра показывает разницу: поврежденный режим имеет свои первичные резонансы усилены, но спектры в противном случае перекрываются. Далее создадим способ количественного различения нормального и неисправного состояния. Мы можем использовать такие подходы к кластеризации и классификации данных, как:
Нечеткая кластеризация C-средств. Посмотрите fcm() в панели инструментов нечеткой логики.
Поддержка векторного машинного классификатора. Посмотрите fitcsvm () в окне «Статистика и инструментарий машинного обучения».
Самоорганизующиеся карты. Посмотрите selforgmap() в инструментарии глубокого обучения.
В этом примере мы используем метод классификации Support Vector Machine. Кластеризация информации из двух типов моделей (mNormal и mFault) может основываться на различных видах информации, которую могут предоставить эти модели, таких как местоположения их полюсов и нулей, их местоположения пиковых резонансов или их список параметров. Здесь мы классифицируем моды по местоположениям полюсов, соответствующим двум резонансам. Для кластеризации мы помечаем полюса моделей здорового состояния «хорошим», а полюса моделей неисправного состояния «неисправным».
ModelTags = cell(nr*2,1); % nr is number of data segments ModelTags(1:nr) = {'good'}; ModelTags(nr+1:end) = {'faulty'}; ParData = zeros(nr*2,4); plist = @(p)[real(p(1)),imag(p(1)),real(p(3)),imag(p(3))]; % poles of dominant resonances for ct = 1:nr ParData(ct,:) = plist(esort(pole(mNormal{ct}))); ParData(nr+ct,:) = plist(esort(pole(mFault{ct}))); end cl = fitcsvm(ParData,ModelTags,'KernelFunction','rbf', ... 'BoxConstraint',Inf,'ClassNames',{'good', 'faulty'}); cl.ConvergenceInfo.Converged
ans = logical 1
cl - классификатор SVM, разделяющий обучающие данные. ParData в хорошие и неисправные регионы. Использование predict способ этого классификатора может назначать входной вектор nx-by-1 одной из двух областей.
Теперь мы можем протестировать классификатор для его прогнозирования (нормальный или поврежденный) собирать пакеты данных из системы, параметры которой изменяются таким образом, что они переходят от работоспособности (mode = 'Normal') к полному повреждению (mode = 'TaxingSystem') непрерывным образом. Чтобы смоделировать этот сценарий, переведем модель в режим «Ухудшенная система».
set_param([sysA,'/Mechanical System'],'LabelModeActiveChoice','DeterioratingSystem'); sim(sysA) ytv = logsout.getElement('y').Values; ytv = squeeze(ytv.Data); PredictedMode = cell(nr,1); for ct = 1:nr zSegment = iddata(ytv((ct-1)*512+(1:500)),[],Ts); mSegment = ssest(zSegment, nx, 'form', 'canonical', 'Ts', Ts); PredictedMode(ct) = predict(cl, plist(esort(pole(mSegment)))); end I = strcmp(PredictedMode,'good'); Tags = ones(nr,1); Tags(~I) = -1; t = (0:5120)'*Ts; % simulation time Time = t(1:512:end-1); plot(Time(I),Tags(I),'g*',Time(~I),Tags(~I),'r*','MarkerSize',12) grid on axis([0 20 -2 2]) title('Green: Normal, Red: Faulty state') xlabel('Data evaluation time') ylabel('Prediction')
График показывает, что классификатор предсказывает поведение примерно до средней точки как нормальное и после этого находится в состоянии отказа.
В предыдущем анализе использовались пакеты данных, собранные в разное время во время работы системы. Альтернативным, часто более удобным, способом мониторинга работоспособности системы является создание адаптивной модели её поведения. Новые измерения обрабатываются непрерывно и используются для рекурсивного обновления параметров модели. На влияние износа или неисправности указывает изменение значений параметров модели.
Рассмотрим сценарий износа еще раз. По мере старения системы происходит большее «бряцание», которое проявляется как возбуждение нескольких резонансных режимов, а также увеличение пиковой реакции системы. Этот сценарий описан в модели. pdmDeterioratingSystemEstimation который совпадает с режимом «Ухудшенная система» pdmMechanicalSystem за исключением того, что импульсные удары, добавленные для автономной идентификации, отсутствуют. Ответ системы передается в блок «Recursive Polynomial Model Estimator», который сконфигурирован для оценки параметров структуры модели ARMA. Фактическая система запускается в исправном состоянии, но ухудшается до состояния окончания срока службы в течение 200 секунд.
initial_model = translatecov(@(x)idpoly(x),model);
sysB = 'pdmDeterioratingSystemEstimation';
open_system(sysB);
Блок «ARMA-модель» инициализирован с использованием параметров и ковариационных данных из оценочной модели нормального поведения, полученной в предыдущем разделе после преобразования в формат полинома (ARMA). translatecov() функция используется так, что данные ковариации параметра также преобразуются. Блок использует алгоритм «Коэффициент забывания» с коэффициентом забывания, установленным чуть меньше 1, для обновления параметров в каждый момент выборки. Выбор фактора забывания влияет на скорость обновления системы. Небольшое значение означает, что обновления будут иметь высокую дисперсию, в то время как большое значение затруднит адаптацию оценщика к быстрым изменениям.
Оценка параметров модели используется для обновления выходного спектра и его доверительной области 3-sd. Система будет явно изменена, когда доверительная область спектра не будет перекрывать доверительную область здоровой системы на интересующих частотах. Порог обнаружения неисправности показан с использованием черной линии на графике, обозначающей максимально допустимые коэффициенты усиления на определенных частотах. По мере накопления изменений в системе спектр смещается через эту линию. Это служит визуальным индикатором неисправности, которая может использоваться для выполнения ремонтных работ в реальных системах.
Запустите моделирование и просмотрите график спектра по мере его обновления.
sim(sysB)
Текущие оценки параметров модели также используются для вычисления местоположений полюсов системы, которые затем подаются в классификатор SVM для прогнозирования того, находится ли система в «хорошем» или «неисправном» состоянии. Это решение также отображается на графике. Когда нормированный балл предсказания меньше 0,3, решение считается ориентировочным (близким к границе различения). См. сценарий pdmARMASpectrumPlot.m для получения подробной информации о том, как вычисляется текущая оценка предсказания спектра и классификатора.
Можно реализовать процедуру адаптивной оценки и печати вне Simulink с помощью recursiveARMA() функция. Как блок «Recursive Polynomial Model Estimator», так и блок recursiveARMA() создание кода поддержки функций для целей развертывания.
Схема классификации может быть обобщена на случай, когда существует несколько известных способов отказа. Для этого нам понадобятся многогрупповые классификаторы, где режим относится к определенному типу отказа. В этом примере эти аспекты не рассматриваются.
Этот пример показал, как системные схемы идентификации в сочетании с подходами к кластеризации и классификации данных могут помочь в обнаружении и изоляции неисправностей. Были обсуждены как последовательный пакетный анализ, так и схемы адаптации в режиме онлайн. Была идентифицирована модель структуры ARMA измеряемого выходного сигнала. Подобный подход может быть принят в ситуациях, когда имеется доступ как к входным, так и к выходным сигналам, и хотел бы использовать другие типы модельных структур, такие как модели State-space или Box-Jenkins.
В этом примере мы обнаружили, что:
Корреляции в остатках, основанные на модели нормальной работы, могут указывать на начало отказа.
Постепенно ухудшающиеся неисправности могут быть обнаружены путем использования непрерывно адаптирующейся модели поведения системы. Заданные пороговые значения характеристик модели, такие как границы ее выходного спектра, могут помочь визуализировать начало и прогрессирование отказов.
Когда источник неисправности необходимо изолировать, жизнеспособный подход заключается в том, чтобы заранее создать отдельные модели соответствующих режимов отказа. Затем можно использовать классификационный подход для присвоения прогнозируемого состояния системы одному из этих режимов.