Введение

Раннее обнаружение и изоляция аномалий в операции машины может помочь уменьшить аварии, уменьшить время простоя и, таким образом, сэкономить эксплуатационные расходы. Подход включает обработку живых измерений от операции системы, чтобы отметить любое неожиданное поведение, которое будет указывать на недавно разработанный отказ.

Этот пример исследует следующие аспекты диагностики отказа:

Идентификация динамической модели поведения системы

В основанном на модели подходе к обнаружению сначала строят динамическую модель соответствующей системы с использованием измеренных входных и выходных данных. Хорошая модель способна точно предсказать ответ системы для определенного будущего временного горизонта. Когда предсказание не является хорошим, невязки могут быть большими и могут содержать корреляции. Эти аспекты используются для выявления случаев отказа.

Рассмотрите создание, подверженное влияниям и вибрации. Источником вибраций могут быть различные типы раздражителей в зависимости от системы, такие как порывы ветра, контакт с работающими двигателями и турбинами или колебания земли. Эти влияния являются результатом импульсных тестов на отбойник в системе, которые добавляются, чтобы возбуждать систему в достаточной степени. Модели 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)')

График спектра показывает различие: поврежденный режим имеет свои первичные резонансы усилены, но спектры в противном случае перекрываются. Затем мы создадим способ количественного различения нормального и дефектного состояний. Мы можем использовать кластеризацию данных и классификационные подходы, такие как:

В этом примере мы используем Машину опорных векторов метод классификации. Кластеризация информации из двух типов моделей (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.

В этом примере мы обнаружили, что: