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

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

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

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

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

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

Здоровый насос был запущен в ссылочной области значений скорости 0 - 3 000 об/мин в замкнутом цикле с контроллером с обратной связью. Ссылочный вход является измененным сигналом PRBS. Измерения для моторного крутящего момента, крутящего момента насоса, скорости и давления были собраны в частоте дискретизации на 10 Гц. Загрузите измеренные сигналы и постройте ссылочные и фактические скорости насоса.

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

Образцовая идентификация

A. Статическая идентификация модели насоса

Оцените параметры ${\theta }_1$ и ${\theta }_2$ в статическом уравнении насоса с помощью измеренных значений скорости насоса $\omega \left(\mathit{t}\right)$и перепад давления $\Delta \mathit{p}\left(\mathit{t}\right)$ как данные ввода - вывода. Смотрите, что помощник функционирует 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. Динамическая идентификация модели канала

Оцените параметры ${\theta }_3$, ${\theta }_4$ и ${\theta }_5$ в канале разряжают уравнение потока $$\underset{}{\overset{ˆ}{Q}}(t)={\theta }_3+{\theta }_4\sqrt{\Delta p(t)}+{\theta }_5\underset{}{\overset{ˆ}{Q}}(t-1)$$, использование измеренных значений скорости потока жидкости $\mathit{Q}\left(\mathit{t}\right)$и перепад давления $\Delta \mathit{p}\left(\mathit{t}\right)$ как данные ввода - вывода. Смотрите, что помощник функционирует 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')

C. Динамическая идентификация модели канала насоса

Динамическая модель канала насоса использует те же параметры, идентифицированные выше ($${\theta }_3,{\theta }_4,{\theta }_5$$) за исключением того, что симуляция модели требует использования предполагаемого перепада давлений, а не измеренного. Следовательно никакая новая идентификация не требуется. Проверяйте что ориентировочные стоимости $${\theta }_3,{\theta }_4,{\theta }_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. Динамическая обратная идентификация модели насоса

Параметры $${\theta }_6,...,{\theta }_{\mathrm{}10}$$ может быть идентифицирован подобным образом, путем регрессирования измеренных значений крутящего момента на предыдущем крутящем моменте и измерениях скорости. Однако полная многоскоростная симуляция получившейся кусочной линейной модели не обеспечивает хорошую подгонку к данным. Следовательно различный подход моделирования черного квадрата пробуют, который связал идентификацию модели Nonlinear 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. Набор эксперимента повторяется для каждого из следующих режимов:

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

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
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'

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

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

Начните анализ визуальным осмотром функций. Для этого рассмотрите Отказ 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 (Среднее значение r1) и OneNorm1 (1 норма 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 режима можно отличить на одной или нескольких парах функций. Обучите 20 членских классификаторов Мешконасыпателя Дерева уменьшаемому набору отказов (датчик смещает только), использование уменьшаемого набора функций.

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')

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

B. Классификация мультиклассов с помощью Приложения Classification Learner

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

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

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

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

Основанный на невязке подход диагностики отказа включает следующие шаги:

Ссылки

Поддерживание функций

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

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