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

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

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

Используя статическое уравнение крышки насоса и динамические уравнения трубопровода насоса, могут быть вычислены эти 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