Введение

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

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

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

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

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

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

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

В этом примере мы нашли что: