Обнаружение отказов центробежных насосов с использованием анализа невязок

Этот пример использует:

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

Пример следует анализу центробежного насоса, представленному в книге «Приложения диагностики отказа» Рольфа Изермана [1]. Он использует функциональность от System Identification Toolbox™, Statistics and Machine Learning Toolbox™, Control System Toolbox™ и Simulink™ и не требует Predictive Maintenance Toolbox™.

Многоскоростные запуски - диагностика путем остаточного анализа

Установившаяся головка насоса и уравнения крутящего момента не дают точных результатов, если насос запускается с быстро изменяющейся или более широкой областью значений скоростей. Трение и другие потери могут стать значительными, и параметры модели показывают зависимость от скорости. Широко применимым подходом в таких случаях является создание модели поведения черного ящика. Параметры таких моделей не должны быть физически значимыми. Модель используется как устройство для симуляции известного поведения. Выходы модели вычитаются из соответствующих измеренных сигналов, чтобы вычислить невязки. Свойства невязок, такие как их среднее значение, отклонение и степень, используются, чтобы различать нормальные и дефектные операции.

Используя статическое уравнение головки насоса и динамические уравнения насоса-трубопровода, можно вычислить 4 невязки, как показано на рисунке.

Модель имеет следующие компоненты:

Статическая модель насоса: $Δ p_{}^{ˆ} (t) = θ_{1} ω^{2} (t) + θ_{2} ω (t)$
Динамическая модель трубопровода: $Q_{}^{ˆ} (t) = θ_{3} + θ_{4} \sqrt{Δ p (t)} + θ_{5} Q_{}^{ˆ} (t - 1)$
Динамическая модель насоса-трубопровода: ${Q_{}^{ˆ}}_{}^{ˆ} (t) = θ_{3} + θ_{4} \sqrt{{Δ p (t)}_{}^{ˆ}} + θ_{5} {Q_{}^{ˆ}}_{}^{ˆ} (t - 1)$
Динамическая модель обратного насоса: ${M_{}^{ˆ}}_{m o t o r} (t) = θ_{6} + θ_{7} ω (t) + θ_{8} ω (t - 1) + θ_{9} ω^{2} (t) + θ_{10} {M_{}^{ˆ}}_{m o t o r} (t - 1)$

Параметры модели $θ_{1}, . . ., θ_{10}$ показать зависимость от скорости насоса. В этом примере вычисляется кусочно-линейное приближение для параметров. Разделите операционную область на 3 режима:

$ω \leq 900 R P M$
$900 < ω \leq 1500 R P M$
$ω > 1500 R P M$

Исправный насос вращался над задающей скоростью области значений 0-3000 об/мин в закрытом цикле с помощью контроллера с обратной связью. Входной параметр ссылки является измененным сигналом PRBS. Измерения крутящего момента двигателя, крутящего момента насоса, скорости и давления были собраны с частотой дискретизации 10 Гц. Загрузите измеренные сигналы и постройте график ссылки и фактических скоростей насоса (это большие наборы данных, ~ 20 МБ, которые загружаются с сайта файлов поддержки MathWorks).

url = 'https://www.mathworks.com/supportfiles/predmaint/fault-diagnosis-of-centrifugal-pumps-using-residual-analysis/DynamicOperationData.mat';
websave('DynamicOperationData.mat',url);
load DynamicOperationData
figure
plot(t, RefSpeed, t, w)
xlabel('Time (s)')
ylabel('Pump Speed (RPM)')
legend('Reference','Actual')

Задайте рабочие режимы на основе областей значений скорости насоса.

I1 = w<=900;            % first operating regime
I2 = w>900 & w<=1500;   % second operating regime
I3 = w>1500;            % third operating regime

Идентификация модели

А. Идентификация статической модели насоса

Оцените параметры $θ_{1}$ и $θ_{2}$ в уравнении статического насоса с использованием измеренных значений скорости насоса $ω (t)$ и перепада давления $Δ p (t)$ в качестве входно-выходных данных. Смотрите функцию helper staticPumpEst который выполняет эту оценку.

th1 = zeros(3,1);  
th2 = zeros(3,1); 
dpest = nan(size(dp));  % estimated pressure difference
[th1(1), th2(1), dpest(I1)] = staticPumpEst(w, dp, I1);  % Theta1, Theta2 estimates for regime 1
[th1(2), th2(2), dpest(I2)] = staticPumpEst(w, dp, I2);  % Theta1, Theta2 estimates for regime 2
[th1(3), th2(3), dpest(I3)] = staticPumpEst(w, dp, I3);  % Theta1, Theta2 estimates for regime 3
plot(t, dp, t, dpest) % compare measured and predicted pressure differential
xlabel('Time (s)')
ylabel('\Delta P')
legend('Measured','Estimated','Location','best')
title('Static Pump Model Validation')

B. Идентификация динамической модели трубопровода

Оцените параметры $θ_{3}$ , $θ_{4}$ и $θ_{5}$ в уравнении расхода разряда трубопровода $Q_{}^{ˆ} (t) = θ_{3} + θ_{4} \sqrt{Δ p (t)} + θ_{5} Q_{}^{ˆ} (t - 1)$ , с использованием измеренных значений скорости потока жидкости $Q (t)$ и перепада давления $Δ p (t)$ в качестве входно-выходных данных. Смотрите функцию helper dynamicPipeEst который выполняет эту оценку.

th3 = zeros(3,1);  
th4 = zeros(3,1); 
th5 = zeros(3,1);
[th3(1), th4(1), th5(1)] = dynamicPipeEst(dp, Q, I1); % Theta3, Theta4, Theta5 estimates for regime 1
[th3(2), th4(2), th5(2)] = dynamicPipeEst(dp, Q, I2); % Theta3, Theta4, Theta5 estimates for regime 2
[th3(3), th4(3), th5(3)] = dynamicPipeEst(dp, Q, I3); % Theta3, Theta4, Theta5 estimates for regime 3

В отличие от случая со статической моделью насоса, динамическая модель трубопровода показывает динамическую зависимость от значений скорости потока жидкости. Чтобы симулировать модель в меняющихся режимах скорости, в Simulink создается кусочно-линейная модель с использованием блока «LPV System» Control System Toolbox. Смотрите модель Simulink LPV_pump_pipe и вспомогательную функцию simulatePumpPipeModel который выполняет симуляцию.

% Check Control System Toolbox availability
ControlsToolboxAvailable = ~isempty(ver('control')) && license('test', 'Control_Toolbox');
if ControlsToolboxAvailable
    % Simulate the dynamic pipe model. Use measured value of pressure as input
    Ts = t(2)-t(1);
    Switch = ones(size(w));
    Switch(I2) = 2;
    Switch(I3) = 3;
    UseEstimatedP = 0;
    Qest_pipe = simulatePumpPipeModel(Ts,th3,th4,th5);
    plot(t,Q,t,Qest_pipe) % compare measured and predicted flow rates
else
    % Load pre-saved simulation results from the piecewise linear Simulink model
    load DynamicOperationData Qest_pipe
    Ts = t(2)-t(1);
    plot(t,Q,t,Qest_pipe)
end
xlabel('Time (s)')
ylabel('Flow rate (Q), m^3/s')
legend('Measured','Estimated','Location','best')
title('Dynamic Pipe Model Validation')

С. Идентификация модели трубопровода динамического насоса

В модели Dynamic pump-pipe используются те же параметры, идентифицированные выше ( $θ_{3}, θ_{4}, θ_{5}$ ) за исключением того, что симуляция модели требует использования расчетного давления различия а не измеренного. Следовательно, новая идентификация не требуется. Проверяйте, что оценочные значения $θ_{3}, θ_{4}, θ_{5}$ обеспечивают хорошее воспроизведение динамики насоса-трубопровода.

if ControlsToolboxAvailable
    UseEstimatedP = 1;
    Qest_pump_pipe = simulatePumpPipeModel(Ts,th3,th4,th5);
    plot(t,Q,t,Qest_pump_pipe) % compare measured and predicted flow rates
else
    load DynamicOperationData Qest_pump_pipe 
    plot(t,Q,t,Qest_pump_pipe)
end

xlabel('Time (s)')
ylabel('Flow rate Q (m^3/s)')
legend('Measured','Estimated','location','best')
title('Dynamic Pump-Pipe Model Validation')

Подгонка фактически идентична той, которая получена с использованием измеренных значений давления.

D. Идентификация модели динамического обратного насоса

Параметры $θ_{6}, . . ., θ_{10}$ можно идентифицировать аналогичным образом, путем регрессии измеренных значений крутящего момента на предыдущих измерениях крутящего момента и скорости. Однако полная многоскоростная симуляция полученной кусочно-линейной модели не обеспечивает хорошей подгонки данных. Следовательно, используется другой подход моделирования черного ящика, который включает идентификацию нелинейной модели ARX с рациональными регрессорами, чтобы соответствовать данным.

% Use first 300 samples out of 550 for identification
N = 350;
sys3 = identifyNonlinearARXModel(Mmot,w,Q,Ts,N)

sys3 =
Nonlinear ARX model with 1 output and 2 inputs
 Inputs: u1, u2
 Outputs: y1
 Standard regressors corresponding to the orders:
   na = [2]
   nb = [2 1]
   nk = [0 1]
 Custom regressors:
   u1(t-2)^2
   u1(t)*u2(t-2)
   u2(t)^2
 Nonlinear regressors:
  none
 Model output is linear in regressors.
Sample time: 0.1 seconds

Status:                                          
Estimated using NLARX on time domain data.       
Fit to estimation data: 50.55% (simulation focus)
FPE: 1759, MSE: 3214

Mmot_est = sim(sys3,[w Q]);
plot(t,Mmot,t,Mmot_est) % compare measured and predicted motor torque
xlabel('Time (s)')
ylabel('Motor Torque (Nm)')
legend('Measured','Estimated','location','best')
title('Inverse pump model validation')

Генерация остатков

Задайте остаток модели как различие между измеренным сигналом и соответствующим выходным выходом модели. Вычислите четыре невязки, соответствующие четырем компонентам модели.

r1 = dp - dpest;
r2 = Q - Qest_pipe;
r3 = Q - Qest_pump_pipe;

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

r4 = Mmot - movmean(Mmot_est,[1 5]);

Взгляд на остатки обучения:

figure
subplot(221)
plot(t,r1)
ylabel('Static pump - r1')
subplot(222)
plot(t,r2)
ylabel('Dynamic pipe - r2')
subplot(223)
plot(t,r3)
ylabel('Dynamic pump-pipe - r3')
xlabel('Time (s)')
subplot(224)
plot(t,r4)
ylabel('Dynamic inverse pump - r4')
xlabel('Time (s)')

Редукция данных остатка

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

Рассмотрим набор из 20 экспериментов по системе насосов с использованием входных реализаций PRBS. Набор экспериментов повторяется для каждого из следующих режимов:

Исправный насос
Отказ 1: Износ зазора в зазоре
Отказ 2: Небольшие отложения на выходе рабочего колеса
Отказ 3: Отложения на входном отверстии рабочего колеса
Отказ 4: Абразивный износ на выходе рабочего колеса
Отказ 5: Сломанный блейд
Отказ 6: Кавитация
Отказ 7: Смещение датчика скорости
Отказ 8: Смещение расходомера
Отказ 9: Смещение датчика давления

Загрузите экспериментальные данные

url = 'https://www.mathworks.com/supportfiles/predmaint/fault-diagnosis-of-centrifugal-pumps-using-residual-analysis/MultiSpeedOperationData.mat';
websave('MultiSpeedOperationData.mat',url);
load MultiSpeedOperationData
% Generate operation mode labels
Labels = {'Healthy','ClearanceGapWear','ImpellerOutletDeposit',...
    'ImpellerInletDeposit','AbrasiveWear','BrokenBlade','Cavitation','SpeedSensorBias',...
    'FlowmeterBias','PressureSensorBias'};

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

% HealthyR = helperComputeEnsembleResidues(HealthyEnsemble,Ts,sys3,th1,th2,th3,th4,th5); % Healthy data residuals

% Load pre-saved data from "helperComputeEnsembleResidues" run
url = 'https://www.mathworks.com/supportfiles/predmaint/fault-diagnosis-of-centrifugal-pumps-using-residual-analysis/Residuals.mat';
websave('Residuals.mat',url);
load Residuals

Функция остатков, которые имели бы наибольшую степень дискриминации видов, не известна априори. Поэтому сгенерируйте несколько функции кандидата: среднюю, максимальную амплитуду, отклонение, куртоз и 1-норму для каждого остаточного. Все функции масштабируются с помощью области значений значений в здоровом ансамбле.

CandidateFeatures = {@mean, @(x)max(abs(x)), @kurtosis, @var, @(x)sum(abs(x))};
FeatureNames = {'Mean','Max','Kurtosis','Variance','OneNorm'};
% generate feature table from gathered residuals of each fault mode
[HealthyFeature, MinMax] = helperGenerateFeatureTable(HealthyR, CandidateFeatures, FeatureNames);
Fault1Feature  = helperGenerateFeatureTable(Fault1R,  CandidateFeatures, FeatureNames, MinMax);
Fault2Feature  = helperGenerateFeatureTable(Fault2R,  CandidateFeatures, FeatureNames, MinMax);
Fault3Feature  = helperGenerateFeatureTable(Fault3R,  CandidateFeatures, FeatureNames, MinMax);
Fault4Feature  = helperGenerateFeatureTable(Fault4R,  CandidateFeatures, FeatureNames, MinMax);
Fault5Feature  = helperGenerateFeatureTable(Fault5R,  CandidateFeatures, FeatureNames, MinMax);
Fault6Feature  = helperGenerateFeatureTable(Fault6R,  CandidateFeatures, FeatureNames, MinMax);
Fault7Feature  = helperGenerateFeatureTable(Fault7R,  CandidateFeatures, FeatureNames, MinMax);
Fault8Feature  = helperGenerateFeatureTable(Fault8R,  CandidateFeatures, FeatureNames, MinMax);
Fault9Feature  = helperGenerateFeatureTable(Fault9R,  CandidateFeatures, FeatureNames, MinMax);

В каждой таблице функций 20 функций (5 функции для каждого сигнала остатка). Каждая таблица содержит 50 наблюдений (строк), по одному из каждого эксперимента.

N = 50; % number of experiments in each mode
FeatureTable = [...
   [HealthyFeature(1:N,:), repmat(Labels(1),[N,1])];...
   [Fault1Feature(1:N,:),  repmat(Labels(2),[N,1])];...
   [Fault2Feature(1:N,:),  repmat(Labels(3),[N,1])];...
   [Fault3Feature(1:N,:),  repmat(Labels(4),[N,1])];...
   [Fault4Feature(1:N,:),  repmat(Labels(5),[N,1])];...
   [Fault5Feature(1:N,:),  repmat(Labels(6),[N,1])];...
   [Fault6Feature(1:N,:),  repmat(Labels(7),[N,1])];...
   [Fault7Feature(1:N,:),  repmat(Labels(8),[N,1])];...
   [Fault8Feature(1:N,:),  repmat(Labels(9),[N,1])];...
   [Fault9Feature(1:N,:),  repmat(Labels(10),[N,1])]];
FeatureTable.Properties.VariableNames{end} = 'Condition';

% Preview some samples of training data
disp(FeatureTable([2 13 37 49 61 62 73 85 102 120],:))

     Mean1       Mean2      Mean3      Mean4      Max1        Max2        Max3        Max4       Kurtosis1    Kurtosis2    Kurtosis3    Kurtosis4     Variance1    Variance2    Variance3    Variance4    OneNorm1    OneNorm2    OneNorm3    OneNorm4            Condition        
    ________    _______    _______    _______    _______    ________    ________    _________    _________    _________    _________    __________    _________    _________    _________    _________    ________    ________    ________    ________    _________________________

     0.32049     0.6358     0.6358    0.74287    0.39187     0.53219     0.53219     0.041795      0.18957    0.089492     0.089489              0     0.45647       0.6263       0.6263      0.15642      0.56001    0.63541     0.63541      0.24258    {'Healthy'              }
     0.21975    0.19422    0.19422    0.83704          0     0.12761     0.12761      0.27644      0.18243     0.19644      0.19643        0.20856    0.033063      0.16772      0.16773     0.025119     0.057182    0.19344     0.19344     0.063167    {'Healthy'              }
      0.1847    0.25136    0.25136    0.87975    0.32545     0.13929      0.1393     0.084162      0.72346    0.091488     0.091485       0.052552    0.063757      0.18371      0.18371     0.023845      0.10671    0.25065     0.25065      0.04848    {'Healthy'              }
           1          1          1          0    0.78284     0.94535     0.94535       0.9287      0.41371     0.45927      0.45926        0.75934     0.92332            1            1      0.80689      0.99783          1           1       0.9574    {'Healthy'              }
     -2.6545    0.90555    0.90555    0.86324     1.3037      0.7492      0.7492      0.97823    -0.055035     0.57019      0.57018         1.4412      1.4411      0.73128      0.73128      0.80573       3.2084    0.90542     0.90542      0.49257    {'ClearanceGapWear'     }
    -0.89435    0.43384    0.43385     1.0744    0.82197     0.30254     0.30254    -0.022325      0.21411    0.098504       0.0985      -0.010587     0.55959      0.29554      0.29554     0.066693       1.0432    0.43326     0.43326      0.13785    {'ClearanceGapWear'     }
     -1.2149    0.53579    0.53579     1.0899    0.87439     0.47339     0.47339      0.29547     0.058946    0.048138     0.048137        0.17864     0.79648      0.48658      0.48658      0.25686       1.4066     0.5353      0.5353      0.26663    {'ClearanceGapWear'     }
     -1.0949    0.44616    0.44616      1.143    0.56693     0.19719     0.19719    -0.012055     -0.10819     0.32603      0.32603     -0.0033592     0.43341      0.12857      0.12857     0.038392       1.2238    0.44574     0.44574     0.093243    {'ClearanceGapWear'     }
     -0.1844    0.16651    0.16651    0.87747    0.25597    0.080265    0.080265      0.49715       0.5019     0.29939      0.29938         0.5338    0.072397     0.058024     0.058025       0.1453      0.20263    0.16566     0.16566      0.14216    {'ImpellerOutletDeposit'}
     -3.0312     0.9786     0.9786    0.75241      1.473     0.97839     0.97839       1.0343    0.0050647      0.4917       0.4917        0.89759      1.7529       0.9568       0.9568       0.8869       3.6536    0.97859     0.97859      0.62706    {'ImpellerOutletDeposit'}

Проект классификатора

А. Визуализация разделяемости режима с помощью графика поля точек

Начните анализ путем визуального осмотра функций. Для этого рассмотрите Отказ 1: Износ зазора. Чтобы просмотреть, какие функции наиболее подходят для обнаружения этого отказа, сгенерируйте график поля точек признаков с метками 'Healthy' и 'ClearanceGapWear'.

T = FeatureTable(:,1:20);
P = T.Variables;
R = FeatureTable.Condition;
I = strcmp(R,'Healthy') | strcmp(R,'ClearanceGapWear');
f = figure;
gplotmatrix(P(I,:),[],R(I))
f.Position(3:4) = f.Position(3:4)*1.5;

Хотя функции в столбцах 1 и 17 четко не видны, они обеспечивают наибольшее разделение. Анализируйте эти функции более тщательно.

f = figure;
Names = FeatureTable.Properties.VariableNames;
J = [1 17];
fprintf('Selected features for clearance gap fault: %s\n',strjoin(Names(J),', '))

Selected features for clearance gap fault: Mean1, OneNorm1

gplotmatrix(P(I,[1 17]),[],R(I))

График теперь четко показывает, что функции Mean1 и OneNorm1 могут использоваться, чтобы отделить исправный режим от режима отказа зазора. Аналогичный анализ может быть выполнен для каждого режима отказа. Во всех случаях можно найти набор функций, которые различают режимы отказа. Следовательно, обнаружение дефектного поведения всегда возможно. Однако изоляция отказов сложнее, поскольку на те же функции влияет несколько типов отказов. Для примера функций Mean1 (Mean of r1) и OneNorm1 (1-norm of r1) показывают изменение для многих типов отказов. Тем не менее, некоторые отказы, такие как смещения датчика, легче изолировать, когда отказ разделим во многих функциях.

Для трех отказов смещения датчика выберите функции из ручного осмотра графика поля точек.

figure;
I = strcmp(R,'Healthy') | strcmp(R,'PressureSensorBias') | strcmp(R,'SpeedSensorBias') | strcmp(R,'FlowmeterBias');
J = [1 4 6 16 20]; % selected feature indices
fprintf('Selected features for sensors'' bias: %s\n',strjoin(Names(J),', '))

Selected features for sensors' bias: Mean1, Mean4, Max2, Variance4, OneNorm4

gplotmatrix(P(I,J),[],R(I))

График поля точек выбранных функций показывает, что 4 режима можно различить на одной или нескольких парах функций. Обучите классификатор Tree Bagger представителей для уменьшенного набора отказов (только смещения датчика) с помощью уменьшенного набора функций.

rng default % for reproducibility
Mdl = TreeBagger(20, FeatureTable(I,[J 21]), 'Condition',...
   'OOBPrediction','on',...
   'OOBPredictorImportance','on');
figure
plot(oobError(Mdl))
xlabel('Number of trees')
ylabel('Misclassification probability')

Ошибка неправильной классификации менее 3%. Таким образом, можно выбрать и работать с меньшим набором функций для классификации определенного подмножества отказов.

B. Многоклассовая классификация с использованием приложения Classification Learner

Предыдущий раздел был посвящен ручному осмотру графиков поля точек, чтобы уменьшить набор функций для конкретных типов отказов. Этот подход может стать утомительным и может не охватывать все типы отказов. Можете ли вы спроектировать классификатор, который может обрабатывать все режимы отказа более автоматизированным образом? В Statistics and Machine Learning Toolbox доступно множество классификаторов. Быстрый способ попробовать многих из них и сравнить их выступления - использовать приложение Classification Learner.

Запустите приложение Classification Learner и выберите FeatureTable из рабочей области в качестве рабочих данных для нового сеанса. Выделите 20% данных (по 10 выборок каждого режима) для валидации отключения.
Выберите Все в разделе Тип модели (Model Type) на основной вкладке. Затем нажмите кнопку Train.
За короткое время обучено около 20 классификаторов. Их точность отображается рядом с именами на панели истории. Линейный классификатор SVM выполняет самое лучшее, получая 86% точности на выборки задержки. Этот классификатор испытывает некоторые трудности с определением «ClearanceGapWear», который он классифицирует как «FlellerOutletDeposit» 40% времени.
Чтобы получить графическое представление о эффективности, откройте график Confusion Matrix из раздела PLOTS на основной вкладке. График показывает эффективность выбранного классификатора (здесь классификатор линейного SVM).

Экспортируйте классификатор лучших результатов в рабочую область и используйте его для предсказания новых измерений.

Сводные данные

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

Подход к диагностике остаточного отказа включает следующие шаги:

Моделируйте динамику между измеряемыми входами и выходами системы с помощью физических факторов или методов системы идентификации черного ящика.
Вычисление остатков как различия между измеренными и моделированными сгенерированными сигналами. Остатки могут потребоваться для дальнейшей фильтрации, чтобы улучшить изоляционность отказа.
Извлеките функции, такие как пиковая амплитуда, степень, куртоз и т. Д. Из каждого остаточного сигнала.
Используйте функции для обнаружения и классификации отказов с помощью методов обнаружения аномалий и классификации.
Не все остатки и выведенные функции чувствительны к каждому отказу. Представление гистограмм функций и графиков поля точек может выявить, какие функции подходят для обнаружения определенного типа отказа. Этот процесс выбора функций и оценки их эффективности для изоляции от отказа может быть итерационной процедурой.

Ссылки

Изерманн, Рольф, Приложения диагностики отказа. Основанный на модели контроль условия: приводы, двигатели, оборудование, объекты, датчики, и отказоустойчивая система, издание 1, Springer-Verlag Берлин Гейдельберг, 2011.

Вспомогательные функции

Оценка параметра статического уравнения насоса

function [x1, x2, dpest] = staticPumpEst(w, dp, I)
%staticPumpEst Static pump parameter estimation in a varying speed setting
% I: sample indices for the selected operating region.

w1 = [0; w(I)];
dp1 = [0; dp(I)];
R1 = [w1.^2 w1];
x = pinv(R1)*dp1;
x1 = x(1);  
x2 = x(2);  

dpest = R1(2:end,:)*x;
end

Оценка динамического параметра трубопровода

function [x3, x4, x5, Qest] = dynamicPipeEst(dp, Q, I)
%dynamicPipeEst Dynamic pipe parameter estimation in a varying speed setting
% I: sample indices for the selected operating region.

Q = Q(I);
dp = dp(I);
R1 = [0; Q(1:end-1)];
R2 = dp; R2(R2<0) = 0; R2 = sqrt(R2);
R = [ones(size(R2)), R2, R1];

% Remove out-of-regime samples
ii = find(I);
j = find(diff(ii)~=1);
R = R(2:end,:); R(j,:) = [];
y = Q(2:end); y(j) = [];
x = R\y;

x3 = x(1);
x4 = x(2);
x5 = x(3);

Qest = R*x;
end

Динамическая многофункциональная симуляция модели насоса-трубопровода с использованием блока LPV System.

function Qest = simulatePumpPipeModel(Ts,th3,th4,th5)
%simulatePumpPipeModel Piecewise linear modeling of dynamic pipe system.
% Ts: sample time
% w: Pump rotational speed
% th1, th2, th3 are estimated model parameters for the 3 regimes.
% This function requires Control System Toolbox.

ss1 = ss(th5(1),th4(1),th5(1),th4(1),Ts);
ss2 = ss(th5(2),th4(2),th5(2),th4(2),Ts);
ss3 = ss(th5(3),th4(3),th5(3),th4(3),Ts);
offset = permute([th3(1),th3(2),th3(3)]',[3 2 1]);
OP = struct('Region',[1 2 3]');
sys = cat(3,ss1,ss2,ss3);
sys.SamplingGrid = OP;

assignin('base','sys',sys)
assignin('base','offset',offset)
mdl = 'LPV_pump_pipe';
sim(mdl);
Qest = logsout.get('Qest');
Qest = Qest.Values;
Qest = Qest.Data;
end

Определите динамическую модель для динамики обратного насоса.

function syse = identifyNonlinearARXModel(Mmot,w,Q,Ts,N)
%identifyNonlinearARXModel Identify a nonlinear ARX model for 2-input (w, Q), 1-output (Mmot) data.
% Inputs:
%  w: rotational speed
%  Q: Flow rate
%  Mmot: motor torque
%  N: number of data samples to use
% Outputs:
%  syse: Identified model
%
% This function uses NLARX estimator from System Identification Toolbox.

sys = idnlarx([2 2 1 0 1],'','CustomRegressors',{'u1(t-2)^2','u1(t)*u2(t-2)','u2(t)^2'});
data = iddata(Mmot,[w Q],Ts);
opt = nlarxOptions;
opt.Focus = 'simulation';
opt.SearchOptions.MaxIterations = 500;
syse = nlarx(data(1:N),sys,opt);
end

Документация