В этом примере показан подход на основе уравнений четности модели для обнаружения и диагностики различных типов неисправностей, возникающих в насосной системе. Этот пример расширяет методы, представленные в документе Диагностика неисправностей центробежных насосов с использованием экспериментов в установившемся состоянии, до ситуации, когда работа насоса охватывает несколько рабочих условий.
В качестве примера приведен анализ центробежных насосов, представленный в книге Rolf Isermann «Приложения диагностики неисправностей» [1]. Он использует функциональные возможности от System Identification Toolbox™, Statistics and Machine Learning Toolbox™, Control System Toolbox™ и Simulink™ и не требует Predictive Maintenance Toolbox™.
Установившиеся уравнения напора насоса и крутящего момента не дают точных результатов, если насос работает с быстро изменяющейся или более широким диапазоном скоростей. Трение и другие потери могут стать значительными, и параметры модели будут зависеть от скорости. Широко применимым подходом в таких случаях является создание модели поведения «черного ящика». Параметры таких моделей не обязательно должны быть физически значимыми. Модель используется в качестве устройства для моделирования известных моделей поведения. Выходные сигналы модели вычитаются из соответствующих измеренных сигналов для вычисления остаточных значений. Свойства остатков, такие как их среднее значение, дисперсия и мощность, используются для различения нормальных и неисправных операций.
Используя статическое уравнение напора насоса и динамические уравнения насос-труба, можно вычислить 4 остатка, как показано на рисунке.
Модель состоит из следующих компонентов:
Модель статического насоса:
Динамическая модель трубы: +θ5Qˆ (t-1)
Динамическая модель трубы насоса: (t-1)
Динамическая обратная модель насоса: θ10Mˆmotor (t-1)
Модельные параметры start10 показывают зависимость от скорости насоса. В этом примере вычисляется кусочно-линейное приближение для параметров. Разделить операционный регион на 3 режима :
1500RPM
Нормальный насос работал в диапазоне эталонных скоростей от 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
A. Идентификация модели статического насоса
В уравнении статического насоса, используя измеренные значения скорости (t) насоса и давлений 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. Идентификация динамических моделей труб
Оцените параметры , и в уравнении потока выброса трубы (t-1), используя измеренные значения Q (t) и перепад давления Δp (t) как данные ввода - вывода. См. функцию помощника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» панели инструментов системы управления. См. модель 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. Идентификация модели трубопровода динамического насоса
В модели «Динамический насос-трубопровод» используются те же параметры, что и в предыдущей модели (start5), за исключением того, что для моделирования модели требуется использовать расчетный перепад давления, а не измеренный. Следовательно, новая идентификация не требуется. Проверьте, что расчетные значения, указанные , start5, дают хорошее воспроизведение динамики насоса-трубы.
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. Идентификация модели динамического обратного насоса
Аналогичным образом могут быть идентифицированы параметры start6,..., start10 путем регрессии измеренных значений крутящего момента при предыдущих измерениях крутящего момента и скорости. Однако полное многоскоростное моделирование полученной кусочно-линейной модели не обеспечивает хорошего соответствия данным. Поэтому пробуют другой подход к моделированию черного ящика, который включает в себя идентификацию нелинейной модели 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'}
A. Визуализация разделяемости режимов с использованием графика рассеяния
Начните анализ с визуального осмотра элементов. Для этого рассмотрим Отказ 1: Износ при зазоре. Чтобы просмотреть, какие функции наиболее подходят для обнаружения этой неисправности, создайте график рассеяния характеристик с метками «Healthy» и «HapeGapWear».
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-членный классификатор 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
В предыдущем разделе основное внимание уделялось ручному контролю графиков рассеяния для уменьшения набора функций для определенных типов отказов. Этот подход может быть утомительным и может охватывать не все типы отказов. Можете ли вы разработать классификатор, который сможет обрабатывать все режимы отказов в более автоматизированном режиме? В окне «Статистика и инструментарий машинного обучения» имеется множество классификаторов. Быстрый способ попробовать многие из них и сравнить их показатели - использовать приложение Classification Learner.
Запустите приложение Classification Learner и выберите элемент FeatureTable в рабочей области в качестве рабочих данных для новой сессии. Отложите 20% данных (10 выборок каждого режима) для проверки отсутствия. 
Выберите Все (All) в разделе Тип модели (Model Type) главной вкладки. Затем нажмите кнопку «Поезд». 
За короткое время обучается около 20 классификаторов. Их точность отображается рядом с их именами под панелью истории. Лучше всего работает линейный классификатор SVM, который обеспечивает 86% точность выборок. Этот классификатор имеет определенные трудности в идентификации «HapleGapWear», который он классифицирует как «ImpellerDepose» 40% времени.
Чтобы получить графическое представление о производительности, откройте график «Матрица путаницы» в разделе PLOTS главной вкладки. На графике показана производительность выбранного классификатора (здесь - Линейный классификатор SVM). 
Экспортируйте наиболее эффективный классификатор в рабочую область и используйте его для прогнозирования новых измерений.
Продуманная стратегия диагностики неисправностей позволяет сократить эксплуатационные расходы за счет минимизации простоев и затрат на замену компонентов. Стратегия использует хорошие знания о динамике рабочей машины, которые используются в сочетании с измерениями датчиков для обнаружения и изоляции различных типов неисправностей. В этом примере описан подход на остаточной основе для диагностики неисправностей центробежных насосов. Этот подход является хорошей альтернативой подходам, основанным на оценке и отслеживании параметров, когда задача моделирования является сложной, а параметры модели показывают зависимость от условий эксплуатации.
Подход к диагностике отказов на остаточной основе включает в себя следующие шаги:
Моделирование динамики между измеряемыми входами и выходами системы с использованием физических соображений или методов идентификации системы черного ящика.
Вычисляют остатки как разность между измеренными и модельными полученными сигналами. Остатки могут нуждаться в дополнительной фильтрации для улучшения возможности устранения неисправности.
Извлекайте из каждого остаточного сигнала такие характеристики, как пиковая амплитуда, мощность, куртоз и т.д.
Используйте функции обнаружения и классификации отказов с использованием методов обнаружения и классификации аномалий.
Не все остатки и производные элементы чувствительны к каждому отказу. Представление гистограмм признаков и графиков рассеяния может показать, какие признаки подходят для обнаружения определенного типа неисправности. Этот процесс выбора функций и оценки их эффективности для локализации отказов может быть итеративной процедурой.
Изерманн, Рольф, приложения для диагностики неисправностей. Мониторинг состояния на основе модели: исполнительные механизмы, приводы, машины, установки, датчики и отказоустойчивая система, издание 1, Springer-Verlag Berlin Heidelberg, 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
Динамическое многорежимное моделирование насосно-трубной модели с использованием блока Система ННД.
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