В этом примере показано, как использовать основанный на данных подход моделирования для обнаружения отказа.
Раннее обнаружение и изоляция аномалий в операции машины может помочь уменьшить аварии, уменьшить время простоя и, таким образом, сэкономить эксплуатационные расходы. Подход включает обработку живых измерений от операции системы, чтобы отметить любое неожиданное поведение, которое будет указывать на недавно разработанный отказ.
Этот пример исследует следующие аспекты диагностики отказа:
Обнаружение ненормального поведения системы путем остаточного анализа
Обнаружение ухудшения при построении моделей поврежденной системы
Отслеживая изменения системы с помощью онлайн-адаптации параметров модели
В основанном на модели подходе к обнаружению сначала строят динамическую модель соответствующей системы с использованием измеренных входных и выходных данных. Хорошая модель способна точно предсказать ответ системы для определенного будущего временного горизонта. Когда предсказание не является хорошим, невязки могут быть большими и могут содержать корреляции. Эти аспекты используются для выявления случаев отказа.
Рассмотрите создание, подверженное влияниям и вибрации. Источником вибраций могут быть различные типы раздражителей в зависимости от системы, такие как порывы ветра, контакт с работающими двигателями и турбинами или колебания земли. Эти влияния являются результатом импульсных тестов на отбойник в системе, которые добавляются, чтобы возбуждать систему в достаточной степени. Модели 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- модели зарегистрированного сигнала как представления системы. Этот подход работает, когда входной сигнал принимается как (фильтрованный) белый шум. Поскольку данные подвержены периодическим ударам, мы разделяем данные на несколько частей каждый, начиная с частоты отбоя. Таким образом, каждый сегмент данных содержит ответ на один bump плюс случайные возбуждения - ситуацию, которая может быть получена с помощью модели временных рядов, где эффект bump приписывается подходящим начальным условиям.
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});
Оцените модель timeseries 7-го порядка в форме пространства состояний с помощью 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
Невязки теперь больше и показывают корреляции при ненулевых лагах. Это основная идея обнаружения отказов - создание невязки метрики и наблюдение за тем, как она изменяется с каждым новым набором измерений. Здесь используется простая невязка, основанный на 1-ступенчатой ошибке предсказания. На практике генерируются более совершенные невязки, адаптированные к потребностям приложения.
Более подробный подход к обнаружению отказа состоит также в том, чтобы идентифицировать модель неисправного (поврежденного) состояния системы. Затем мы можем проанализировать, какая модель с большей вероятностью объяснит живые измерения из системы. Эта схема может быть обобщена для моделей для различных типов отказов и, таким образом, используется не только для обнаружения отказа, но и для определения того, какой из них («изоляция»). В этом примере мы используем следующий подход:
Мы собираем данные с системой, работающей в нормальном (исправном) и известном износостойком состоянии, вызванном окончанием срока службы.
Мы идентифицируем динамическую модель, представляющую поведение в каждом состоянии.
Мы используем подход кластеризации данных, чтобы провести четкое различие между этими состояниями.
Для обнаружения отказа мы собираем данные с работающей машины и идентифицируем модель ее поведения. Затем мы прогнозируем, какое состояние (нормальное или поврежденное) с наибольшей вероятностью объяснит наблюдаемое поведение.
Мы уже моделировали систему в ее режиме нормальной операции. Теперь мы симулируем модель 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-го порядка в форме пространство состояний. Отключите ковариационные расчеты скорости.
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-Means. См. fcm() в Fuzzy Logic Toolbox.
Машина опорных векторов. См. fitcsvm () в Statistics and Machine Learning Toolbox.
Самоорганизующиеся карты. См. selforgmap() в Deep Learning Toolbox.
В этом примере мы используем Машину опорных векторов метод классификации. Кластеризация информации из двух типов моделей (mNormal и mFault) может основываться на различных видах информации, которую эти модели могут предоставить, таких как местоположения их полюсов и нулей, их местоположения пиковых резонансов или их список параметров. Здесь мы классифицируем режимы по положениям полюса, соответствующим двум резонансам. Для кластеризации мы пометим полюсы моделей исправного состояния с 'good' и полюсы моделей дефектного состояния с 'faulty'.
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 на 1 одной из двух областей.
Теперь мы можем протестировать классификатор на его предсказание (нормальный по сравнению с поврежденным), собирать пакеты данных из системы, параметры которой изменяются таким образом, чтобы он непрерывно переходил от исправного (режим = 'Normal') к полному повреждению (режим = 'DamagedSystem'). Чтобы симулировать этот сценарий, мы поместили модель в режим 'DeterioratingSystem'.
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 который аналогичен режиму 'DeterioratingSystem' pdmMechanicalSystem за исключением того, что импульсных ударов, которые были добавлены для оффлайн-идентификации, нет. Ответ системы передается в блок «Recursive Polynomial Model Estimator», который был сконфигурирован, чтобы оценить параметры структуры модели ARMA. Фактическая система начинается в исправном состоянии, но ухудшается до условий окончания срока службы в течение периода времени 200 секунд.
initial_model = translatecov(@(x)idpoly(x),model);
sysB = 'pdmDeterioratingSystemEstimation';
open_system(sysB);
Блок «Модель ARMA» был инициализирован с помощью параметров и ковариационных данных предполагаемой модели нормального поведения, выведенной в предыдущем разделе после преобразования в полиномиальный (ARMA) формат. The translatecov() функция используется так, что параметрические ковариационные данные также преобразуются. Блок использует алгоритм «Забывающего фактора» с коэффициентом забывания, установленным на чуть меньше 1, чтобы обновить параметры в каждый момент дискретизации. Выбор забывающего фактора влияет на то, как быстро система обновляется. Небольшое значение означает, что обновления будут иметь высокое отклонение, в то время как большое значение будет труднее для оценщика адаптироваться к быстрым изменениям.
Оценка параметров модели используется для обновления спектра выхода и его 3-sd доверием области. Система явно изменится, когда доверительная область спектра не совпадает с доверительной областью здоровой системы на интересующих частотах. Порог обнаружения отказа показан с использованием черной линии на графике, помечающем максимально допустимые усиления на определенных частотах. Когда изменения в системе накапливаются, спектр дрейфует через эту линию. Это служит визуальным индикатором отказа, который может использоваться для вызова ремонта в реальных системах.
Запустите симуляцию и смотрите график спектра по мере обновления.
sim(sysB)
Текущие оценки параметров модели также используются, чтобы вычислить местоположения системного полюса, которые затем подаются в классификатор SVM, чтобы предсказать, находится ли система в «хорошем» или «отказе» состоянии. Это решение также отображается на графике. Когда нормированный счет предсказания меньше 3, решение рассматривается ориентировочным (близким к контуру различения). Смотрите скрипт pdmARMASpectrumPlot.m для получения дополнительной информации о том, как вычисляется текущая оценка предсказания спектра и классификатора.
Возможно реализовать процедуру адаптивной оценки и графического изображения вне Simulink с помощью recursiveARMA() функция. Блок «Recursive Polynomial Model Estimator», а также recursiveARMA() генерация кода поддержки функций в целях развертывания.
Схема классификации может быть обобщена в случае, когда существует несколько известных режимов отказа. Для этого нам понадобятся многогрупповые классификаторы, где режим относится к определенному типу отказа. Эти аспекты не исследуются в этом примере.
Этот пример показал, как схемы системы идентификации в сочетании с подходами к кластеризации и классификации данных могут помочь в обнаружении и изоляции отказов. Обсуждались как последовательный пакетный анализ, так и онлайновые схемы адаптации. Была идентифицирована модель структуры ARMA измеренного выходного сигнала. Аналогичный подход может быть принят в ситуациях, когда один имеет доступ как к входу, так и к выходу сигналам, и хотел бы использовать другие типы структур модели, таких как модели полинома State-space или Box-Jenkins.
В этом примере мы обнаружили, что:
Корреляции в невязках, основанные на модели нормальной операции, могут указывать на начало отказа.
Постепенно ухудшающиеся отказы могут быть обнаружены с помощью постоянно адаптирующейся модели поведения системы. Предустановленные пороги характеристик модели, такие как ограничения на ее выходном спектре, могут помочь визуализировать начало и прогрессирование отказов.
Когда источник отказа должен быть изолирован, жизнеспособный подход состоит в том, чтобы заранее создать отдельные модели соответствующих типов отказа. Затем классификационный подход может использоваться, чтобы присвоить предсказанное состояние системы одному из этих режимов.