Этот пример показывает, как использовать основанный на данных подход моделирования для обнаружения отказа. Этот пример использует функциональность от Signal Processing Toolbox™ и System Identification Toolbox™, и не требует Predictive Maintenance Toolbox™.
Раннее обнаружение и изоляция аномалий в работе машины могут помочь уменьшать несчастные случаи, уменьшать время простоя и таким образом сократить операционные затраты. Подход включает обработку живые измерения от операции системы, чтобы отметить любое неожиданное поведение, которое указало бы на недавно разработанный отказ.
Этот пример исследует следующие аспекты диагностики отказа:
Обнаружение аварийного поведения системы остаточным анализом
Обнаружение ухудшения путем создавания моделей поврежденной системы
Изменения системы слежения с помощью онлайновой адаптации параметров модели
В основанном на модели подходе к обнаружению динамическая модель заинтересованной системы сначала создана с помощью измеренных входных и выходных данных. Хорошая модель может точно предсказать ответ системы в течение определенного будущего периода времени. Когда прогноз не хорош, невязки могут быть большими и могли содержать корреляции. Эти аспекты используются, чтобы обнаружить падение отказа.
Считайте создание подвергающимся влиянию и колебаниям. Источник колебаний может быть различными типами стимулов в зависимости от системы, такими как порывы ветра, контакт с рабочими механизмами и турбинами, или основать колебания. Влияние является результатом импульсивных тестов удара в системе, которые добавляются, чтобы взволновать систему достаточно. Модель Simulink idMechanicalSystem.slx является простым примером такой структуры. Возбуждение прибывает из периодических ударов, а также наземных колебаний, смоделированных отфильтрованным белым шумом. Вывод системы собран датчиком, который подвергается шуму измерения. Модель может моделировать различные сценарии, вовлекающие структуру в здоровое или поврежденное состояние.
sysA = 'pdmMechanicalSystem'; open_system(sysA) % Set the model in the healthy mode of operation set_param([sysA,'/Mechanical System'],'OverrideUsingVariant','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-ю модель timeseries порядка в форме пространства состояний с помощью команды 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'],'OverrideUsingVariant','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 шагом. На практике более усовершенствованные невязки сгенерированы, которые нестандартны к потребностям приложения.
Более подробный подход к обнаружению отказа должен также идентифицировать модель дефектного (поврежденного) состояния системы. Мы можем затем анализировать, какая модель, более вероятно, объяснит живые измерения от системы. Это расположение может обобщаться к моделям для различных типов отказов и таким образом использоваться для того, чтобы не только обнаружить отказ, но также и идентифицировать который один ("изоляция"). В этом примере мы проявляем следующий подход:
Мы собираем данные с системой, действующей в (здоровом) нормальном, и известный износ вызвал состояние конца жизненного цикла.
Мы идентифицируем динамическую модель, представляющую поведение в каждом состоянии.
Мы используем подход кластеризации данных, чтобы чертить ясное различие между этими состояниями.
Для обнаружения отказа мы собираем данные от рабочей машины и идентифицируем модель ее поведения. Мы затем предсказываем, какое состояние (нормальный или поврежденный), скорее всего, объяснит наблюдаемое поведение.
Мы уже моделировали систему в ее режиме нормального функционирования. Мы теперь моделируем модель idMechanicalSystem в "конце жизни" режим. Это - сценарий, где система уже ухудшилась к ее конечному состоянию допустимой операции.
set_param([sysA,'/Mechanical System'],'OverrideUsingVariant','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-средств. Смотрите fcm() в Fuzzy Logic Toolbox.
Поддержите векторный классификатор машины. Смотрите fitcsvm () в Statistics and Machine Learning Toolbox.
Самоорганизующиеся карты. Смотрите selforgmap() в Deep Learning Toolbox.
В этом примере мы используем метод классификации Машин Вектора Поддержки. Кластеризация информации от двух типов моделей (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-1 вектор в одну из этих двух областей.
Теперь мы можем протестировать классификатор на его прогноз (нормальный по сравнению с поврежденным), собирают пакеты данных из системы, параметры которой изменяются способом, что он идет от того, чтобы быть здоровым (режим = 'Normal') к тому, чтобы быть полностью поврежденным (режим = 'DamagedSystem') непрерывным способом. Чтобы моделировать этот сценарий, мы помещаем модель в режим 'DeterioratingSystem'.
set_param([sysA,'/Mechanical System'],'OverrideUsingVariant','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')
График показывает, что классификатор предсказывает поведение приблизительно до средней точки, чтобы быть нормальным и в состоянии отказа после этого.
Предыдущий анализ использовал пакеты данных, собранных в разное время во время операции системы. Альтернатива, часто более удобная, способ контролировать здоровье системы, должна создать адаптивную модель своего поведения. Новые измерения обрабатываются постоянно и используются, чтобы обновить параметры модели рекурсивным способом. Эффект износа или отказа обозначается изменением в значениях параметра модели.
Рассмотрите сценарий износа снова. Как системные возрасты, существует больший "грохот", который проявляется как возбуждение нескольких резонирующих режимов, а также повышение максимальной чувствительности системы. Этот сценарий описан в модели idDeterioratingSystemEstimation, которая является тем же самым как режим 'DeterioratingSystem' idMechanicalSystem за исключением того, что импульсивные удары, которые были добавлены для оффлайновой идентификации, не присутствуют. Ответ системы передается блоку "Recursive Polynomial Model Estimator", который был сконфигурирован, чтобы оценить параметры структуры модели ARMA. Фактическая система запускается в здоровом состоянии, но ухудшается к условиям конца жизненного цикла по отрезку времени 200 секунд.
initial_model = translatecov(@(x)idpoly(x),model);
sysB = 'pdmDeterioratingSystemEstimation';
open_system(sysB);
Блок "ARMA model" был инициализирован с помощью параметров и данных о ковариации из предполагаемой модели нормального поведения, выведенного в предыдущем разделе после преобразования в полином (ARMA) формат. Функция translatecov() используется так, чтобы данные о ковариации параметра были также преобразованы. Блок использует "Упущение факторного" алгоритма с набором фактора упущения к немного меньше чем 1, чтобы обновить параметры в каждый момент выборки. Выбор упущения фактора влияет как быстро системные обновления. Маленькое значение означает, что обновления будут иметь высокое отклонение, в то время как большое значение сделает его тяжелее для средства оценки, чтобы адаптироваться к быстрым изменениям.
Оценка параметров модели используется, чтобы обновить выходной спектр и его 3-sd область уверенности. Система ясно изменится, когда область уверенности спектра не перекроет область здоровой системы на частотах интереса. Порог обнаружения отказа показывают с помощью черной линии в графике, отмечающем максимальные позволенные усиления на определенных частотах. Когда изменения в системе накапливаются, дрейфы спектра через эту строку. Это служит визуальным индикатором отказа, который может использоваться, чтобы призвать к восстановлениям в реальных системах.
Запустите симуляцию и см. график спектра, когда это обновляет.
sim(sysB)
Рабочие оценки параметров модели также используются, чтобы вычислить системные местоположения полюса, которые затем поданы в классификатор SVM, чтобы предсказать, находится ли система в "хорошем" состоянии или состоянии "отказа". Это решение также отображено на графике. Когда нормированный счет прогноза - меньше чем.3, решение рассматривается предварительным (близко к контуру различия). См. скрипт pdmARMASpectrumPlot.m для получения дополнительной информации о том, как рабочая оценка спектра и прогноза классификатора вычисляется.
Возможно реализовать адаптивную процедуру оценки и графического вывода вне Simulink с помощью функции recursiveARMA(). Оба блок "Recursive Polynomial Model Estimator", а также recursiveARMA() функционируют генерация кода поддержки в целях развертывания.
Система классификации может быть обобщена к случаю, где существует несколько известных режимов отказа. Для этого нам будут нужны классификаторы мультигруппы, где режим относится к определенному типу отказа. Эти аспекты не исследуются в этом примере.
Этот пример показал, как схемы системы идентификации, объединенные с кластеризацией данных и подходами классификации, могут помочь в обнаружении и изоляции отказов. Оба последовательных пакетных анализа, а также онлайновые схемы адаптации были обсуждены. Модель структуры ARMA измеренного выходного сигнала была идентифицирована. Аналогичный подход может быть принят в ситуациях, где каждый имеет доступ и к сигналам ввода и вывода и хотел бы использовать другие типы образцовых структур, такие как модели полинома Пространства состояний или Поля-Jenkins.
В этом примере мы нашли что:
Корреляции в невязках на основе модели нормального функционирования могут указать на начало отказа.
Постепенно ухудшающиеся отказы могут быть обнаружены путем использования постоянно адаптирующейся модели поведения системы. Предварительно установленные пороги на характеристиках модели, таких как границы на ее выходном спектре могут помочь визуализировать начало и прогрессию отказов.
Когда источник отказа должен быть изолирован, жизнеспособный подход должен создать отдельные модели заинтересованных типов отказа заранее. Затем подход классификации может использоваться, чтобы присвоить предсказанное состояние системы к одному из этих режимов.