Этот пример показывает основанный на уравнениях четности подход модели для обнаружения и диагноза различных типов отказов, которые происходят в системе накачки. Этот пример расширяет методы, представленные в Диагностике отказа Центробежных Насосов Используя Эксперименты Устойчивого состояния к ситуации, где работа насоса охватывает несколько условий работы.
Пример следует за центробежным анализом насоса, представленным в книге Приложений Диагностики отказа Рольфа Изерманна [1]. Это использует функциональность от System Identification Toolbox™, Statistics and Machine Learning Toolbox™, Control System Toolbox™ и Simulink™ и не требует Predictive Maintenance Toolbox™.
Установившаяся крышка насоса и уравнения крутящего момента не приводят к точным результатам, если насос запущен при быстром варьировании или более широкой области значений скоростей. Трение и другие потери могли стать значительными, и параметры модели показывают зависимость от скорости. Широко применимый подход в таких случаях должен создать модель черного квадрата поведения. Параметры таких моделей не должны быть физически значимыми. Модель используется в качестве устройства для симуляции известных поведений. Выходные параметры модели вычтены из соответствующих измеренных сигналов вычислить остаточные значения. Свойства остаточных значений, такие как их среднее значение, отклонение и степень используются, чтобы различать нормальные и дефектные операции.
Используя статическое уравнение крышки насоса и динамические уравнения трубопровода насоса, могут быть вычислены эти 4 остаточных значения как показано на рисунке.
Модель имеет следующих компонентов:
Статическая модель насоса:
Динамическая модель трубопровода:
Динамическая модель трубопровода насоса:
Динамическая обратная модель насоса:
Параметры модели покажите зависимость от скорости насоса. В этом примере вычисляется кусочная линейная аппроксимация для параметров. Разделите операционную область на 3 режима:
Здоровый насос был запущен в области значений задающей скорости 0 - 3 000 об/мин в замкнутом цикле с контроллером с обратной связью. Ссылочный вход является модифицированным сигналом PRBS. Измерения для крутящего момента двигателя, крутящего момента насоса, скорости и давления были собраны в частоте дискретизации на 10 Гц. Загрузите измеренные сигналы и постройте ссылочные и фактические скорости насоса (это большие наборы данных, ~20MB, которые загружаются с сайта файлов поддержки 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. Статическая идентификация модели насоса
Оцените параметры и в статическом уравнении насоса с помощью измеренных значений скорости насоса и перепад давления как данные ввода - вывода. Смотрите, что помощник функционирует 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. Динамическая идентификация модели трубопровода
Оцените параметры , и в трубопроводе разряжают уравнение потока , использование измеренных значений скорости потока жидкости и перепад давления как данные ввода - вывода. Смотрите, что помощник функционирует 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. Динамическая идентификация модели трубопровода насоса
Динамическая модель трубопровода насоса использует те же параметры, идентифицированные выше () за исключением того, что симуляция модели требует использования предполагаемого перепада давлений, а не измеренного. Следовательно никакая новая идентификация не требуется. Проверяйте что ориентировочные стоимости дайте хорошее воспроизведение динамики трубопровода насоса.
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. Динамическая обратная идентификация модели насоса
Параметры может быть идентифицирован подобным образом, путем регрессирования измеренных значений крутящего момента на предыдущем крутящем моменте и измерениях скорости. Однако полная многоскоростная симуляция получившейся кусочной линейной модели не обеспечивает хорошую подгонку к данным. Следовательно различный подход моделирования черного квадрата пробуют, который связал идентификацию модели 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 Regressors: 1. Linear regressors in variables y1, u1, u2 2. Custom regressor: u1(t-2)^2 3. Custom regressor: u1(t)*u2(t-2) 4. Custom regressor: u2(t)^2 List of all regressors Output function: Linear Function Sample time: 0.1 seconds Status: Estimated using NLARX on time domain data. Fit to estimation data: 49.2% (simulation focus) FPE: 1798, MSE: 3392
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' и '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.
Запустите Приложение Classification Learner и выберите FeatureTable из рабочей области как рабочие данные для нового сеанса. Отложите 20% данных (10 выборок каждого режима) для валидации затяжки. 
Выберите All под разделом Model Type основной вкладки. Затем нажмите кнопку Train. 
В скором времени приблизительно 20 классификаторов обучены. Их точность отображена рядом с их именами под панелью истории. Линейный классификатор SVM выполняет лучшую, производящую 86%-ю точность на хранении выборки. Этот классификатор испытывает некоторые затруднения в идентификации "ClearanceGapWear", который это классифицирует как "ImpellerOutletDeposit" 40% времени.
Чтобы получить графическое представление эффективности, откройте Матричный график Беспорядка от раздела 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