В этом примере показано, как использовать основанный на данных подход моделирования для обнаружения отказа.
Раннее обнаружение и изоляция аномалий в работе машины могут помочь уменьшать несчастные случаи, уменьшать время простоя и таким образом сократить операционные затраты. Подход включает обработку живые измерения от операции системы, чтобы отметить любое неожиданное поведение, которое указало бы на недавно разработанный отказ.
Этот пример исследует следующие аспекты диагностики отказа:
Обнаружение аварийного поведения системы остаточным анализом
Обнаружение ухудшения путем создавания моделей поврежденной системы
Изменения системы слежения с помощью онлайновой адаптации параметров модели
В основанном на модели подходе к обнаружению динамическая модель заинтересованной системы сначала создана с помощью измеренных входных и выходных данных. Хорошая модель может точно предсказать ответ системы в течение определенного будущего периода времени. Когда предсказание не хорошо, остаточные значения могут быть большими и могли содержать корреляции. Эти аспекты используются, чтобы обнаружить падение отказа.
Считайте создание удовлетворяющим ударам и колебаниям. Источник колебаний может быть различными типами стимулов в зависимости от системы, такими как порывы ветра, контакт с рабочими механизмами и турбинами, или основать колебания. Удары являются результатом импульсивных тестов удара в системе, которые добавляются, чтобы взволновать систему достаточно. Модель Simulink pdmMechanicalSystem.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%-я точность. Чтобы просмотреть подгонку к другим экспериментам в наборе данных, используйте контекстное меню щелчка правой кнопкой осей графика.
Последний шаг в проверке модели должен анализировать остаточные значения, сгенерированные им. Для хорошей модели эти остаточные значения должны быть белыми, i.e., покажите статистически незначительные корреляции для ненулевых задержек:
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 шагом. На практике более усовершенствованные остаточные значения сгенерированы, которые нестандартны к потребностям приложения.
Более подробный подход к обнаружению отказа должен также идентифицировать модель дефектного (поврежденного) состояния системы. Мы можем затем анализировать, какая модель, более вероятно, объяснит живые измерения от системы. Это расположение может обобщаться к моделям для различных типов отказов и таким образом использоваться для того, чтобы не только обнаружить отказ, но также и идентифицировать который один ("изоляция"). В этом примере мы проявляем следующий подход:
Мы собираем данные с системой, действующей в (здоровом) нормальном, и известный износ вызвал состояние конца жизненного цикла.
Мы идентифицируем динамическую модель, представляющую поведение в каждом состоянии.
Мы используем подход кластеризации данных, чтобы чертить ясное различие между этими состояниями.
Для обнаружения отказа мы собираем данные от рабочей машины и идентифицируем модель ее поведения. Мы затем предсказываем, какое состояние (нормальный или поврежденный), скорее всего, объяснит наблюдаемое поведение.
Мы уже симулировали систему в ее режиме нормального функционирования. Мы теперь симулируем модель pdmMechanicalSystem в "конце жизни" режим. Это - сценарий, где система уже ухудшилась к ее конечному состоянию допустимой операции.
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')
График показывает, что классификатор предсказывает поведение приблизительно до средней точки, чтобы быть нормальным и в состоянии отказа после этого.
Предыдущий анализ использовал пакеты данных, собранных в разное время во время операции системы. Альтернатива, часто более удобная, способ контролировать здоровье системы, должна создать адаптивную модель своего поведения. Новые измерения обрабатываются постоянно и используются, чтобы обновить параметры модели рекурсивным способом. Эффект износа или отказа обозначается изменением в значениях параметра модели.
Рассмотрите сценарий износа снова. Как системные возрасты, существует больший "грохот", который проявляется как возбуждение нескольких резонирующих режимов, а также повышение максимальной чувствительности системы. Этот сценарий описан в модели pdmDeterioratingSystemEstimation который является тем же самым как режим 'DeterioratingSystem' pdmMechanicalSystem за исключением того, что импульсивные удары, которые были добавлены для оффлайновой идентификации, не присутствуют. Ответ системы передается блоку "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.
В этом примере мы нашли что:
Корреляции в остаточных значениях на основе модели нормального функционирования могут указать на начало отказа.
Постепенно ухудшающиеся отказы могут быть обнаружены путем использования постоянно адаптирующейся модели поведения системы. Предварительно установленные пороги на характеристиках модели, таких как границы на ее выходном спектре могут помочь визуализировать начало и прогрессию отказов.
Когда источник отказа должен быть изолирован, жизнеспособный подход должен создать отдельные модели заинтересованных типов отказа заранее. Затем подход классификации может использоваться, чтобы присвоить предсказанное состояние системы к одному из этих режимов.