Этот пример показывает, как оценить передаточную функцию из данных измеренного сигнала.
В этом примере мы оцениваем передаточную функцию для теплообменника. Теплообменник состоит из температуры хладагента, температуры продукта и температуры окружающего воздуха нарушения порядка. Мы оценим теплоноситель к передаточной функции температуры продукта.
Измеренные данные хранятся в файле MATLAB и включают измерения изменений температуры хладагента вокруг номинальной и изменения температуры продукта вокруг номинальной.
load iddemo_heatexchanger_data
Соберите измеренные данные с помощью iddata команда и постройте график.
data = iddata(pt,ct,Ts); data.InputName = '\Delta CTemp'; data.InputUnit = 'C'; data.OutputName = '\Delta PTemp'; data.OutputUnit = 'C'; data.TimeUnit = 'minutes'; plot(data)
Из физики задачи мы знаем, что теплообменник может быть описан системой первого порядка с задержкой. Используйте tfest команда, задающая один полюс, отсутствие нулей и неизвестную задержку ввода/вывода для оценки передаточной функции.
sysTF = tfest(data,1,0,nan)
sysTF =
From input "\Delta CTemp" to output "\Delta PTemp":
1.467
exp(-0.0483*s) * --------
s + 1.56
Continuous-time identified transfer function.
Parameterization:
Number of poles: 1 Number of zeros: 0
Number of free coefficients: 2
Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using TFEST on time domain data "data".
Fit to estimation data: 36%
FPE: 0.02625, MSE: 0.02622
The compare и resid команды позволяют нам исследовать, насколько хорошо оцененная модель соответствует измеренным данным.
set(gcf,'DefaultAxesTitleFontSizeMultiplier',1,... 'DefaultAxesTitleFontWeight','normal',... 'Position',[100 100 780 520]); resid(sysTF,data);
clf compare(data,sysTF)
Рисунок показывает, что невязки сильно коррелируют, что означает, что в измеренных данных есть информация, которая не была адекватно получена оценочной моделью.
Ранее мы оценивали передаточную функцию из данных, но отдельно для системного порядка не включали много знаний априори. Из физики задачи мы знаем, что система стабильна и имеет положительный коэффициент усиления. Проверяя измеренные данные, мы также знаем, что задержка ввода/вывода составляет около 1/5 минуты. Мы используем эту информацию, чтобы создать начальную систему и оценить передаточную функцию, используя эту систему как начальное предположение.
sysInit = idtf(NaN,[1 NaN],'ioDelay',NaN); sysInit.TimeUnit = 'minutes';
Ограничьте числитель передаточной функции и условия знаменателя так, чтобы система была стабильной с положительным усилением.
sysInit.Structure.num.Value = 1; sysInit.Structure.num.Minimum = 0; sysInit.Structure.den.Value = [1 1]; sysInit.Structure.den.Minimum = [0 0];
Ограничьте задержку ввода/вывода области значений [0 1] минутой и используйте 1/5 минуту в качестве начального значения.
sysInit.Structure.ioDelay.Value = 0.2; sysInit.Structure.ioDelay.Minimum = 0; sysInit.Structure.ioDelay.Maximum = 1;
Используйте систему как начальное предположение для задачи оценки
sysTF_initialized = tfest(data,sysInit)
sysTF_initialized =
From input "\Delta CTemp" to output "\Delta PTemp":
1.964
exp(-0.147*s) * ---------
s + 2.115
Continuous-time identified transfer function.
Parameterization:
Number of poles: 1 Number of zeros: 0
Number of free coefficients: 2
Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using TFEST on time domain data "data".
Fit to estimation data: 54.09%
FPE: 0.01351, MSE: 0.01349
resid(sysTF_initialized,data);
clf compare(data,sysTF,sysTF_initialized)
В вышесказанном мы рассматривали задачу оценки как задачу оценки передаточной функции, но знаем, что существует некоторая дополнительная структура, которую мы можем навязать. В частности, известно, что система теплообменника является процессом 1-го порядка с задержкой или моделью 'P1D' формы:
Используйте procest команда для дальнейшего навязывания этой структуры задаче.
sysP1D = procest(data,'P1D')
sysP1D =
Process model with transfer function:
Kp
G(s) = ---------- * exp(-Td*s)
1+Tp1*s
Kp = 0.90548
Tp1 = 0.32153
Td = 0.25435
Parameterization:
{'P1D'}
Number of free coefficients: 3
Use "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using PROCEST on time domain data "data".
Fit to estimation data: 70.4%
FPE: 0.005614, MSE: 0.005607
resid(sysP1D,data);
clf compare(data,sysTF,sysTF_initialized,sysP1D)
Остаточные графики всех оценок, выполненных до сих пор, показывают, что остаточная корреляция все еще высока, что означает, что модель недостаточно богата, чтобы объяснить всю информацию в измеренных данных. Ключом к отсутствующей части является температура окружающего воздуха нарушения порядка, которую мы еще не включили в модель.
Создайте модель процесса 'P1D' с ограничением по значениям задержки и постоянной времени и используйте ее как начальное предположение для задачи оценки.
sysInit = idproc('P1D','TimeUnit','minutes');
Ограничьте модель положительным усилением и задержкой в область значений [0 1] минуту.
sysInit.Structure.Kp.Value = 1; sysInit.Structure.Kp.Minimum = 0; sysInit.Structure.Tp1.Value = 1; sysInit.Structure.Tp1.Maximum = 10; sysInit.Structure.Td.Value = 0.2; sysInit.Structure.Td.Minimum = 0; sysInit.Structure.Td.Maximum = 1;
Задайте опцию для использования модели первого порядка ('ARMA1') для компонента нарушения порядка. Используйте модель шаблона sysInit наряду с набором опций для оценки модели.
opt = procestOptions('DisturbanceModel','ARMA1'); sysP1D_noise = procest(data,sysInit,opt)
sysP1D_noise =
Process model with transfer function:
Kp
G(s) = ---------- * exp(-Td*s)
1+Tp1*s
Kp = 0.91001
Tp1 = 0.3356
Td = 0.24833
An additive ARMA disturbance model has been estimated
y = G u + (C/D)e
C(s) = s + 591.6
D(s) = s + 3.217
Parameterization:
{'P1D'}
Number of free coefficients: 5
Use "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using PROCEST on time domain data "data".
Fit to estimation data: 96.86% (prediction focus)
FPE: 6.307e-05, MSE: 6.294e-05
resid(sysP1D_noise,data);
clf compare(data,sysTF,sysTF_initialized,sysP1D,sysP1D_noise)
Остаточный график ясно показывает, что невязки некоррелированы, что означает, что у нас есть модель, которая объясняет измеренные данные. Компонент нарушения порядка 'ARMA1', который мы оценили, сохранен как значения числителя и знаменателя в свойстве «NoiseTF» модели.
sysP1D_noise.NoiseTF
ans =
struct with fields:
num: {[1 591.6038]}
den: {[1 3.2172]}
Хотя мы определили модель, которая объясняет измеренные данные, мы отмечаем, что модель соответствует измеренным данным составляет около 70%. Потеря значения подгонки является следствием сильного эффекта нарушения порядка окружающей температуры, которое показано следующим образом.
Измеренные данные были получены из модели Simulink со следующими точными значениями (d является нарушением порядка Гауссова шума)
y = (1-pi/100)*exp(-15s)/(21.3s+1)*u + 1/(25s+1)*d
Создайте модель 'P1D' с этими значениями и посмотрите, насколько хорошо эта модель соответствует измеренным данным.
sysReal = idproc('P1D','TimeUnit','minutes'); sysReal.Kp = 1-pi/100; sysReal.Td = 15/60; sysReal.Tp1 = 21.3/60; sysReal.NoiseTF = struct('num',{[1 10000]},'den',{[1 0.04]});
compare(data,sysReal,sysTF,sysTF_initialized,sysP1D,sysP1D_noise);
График сравнения показывает, что истинная система, подобранная к измеренным данным, также составляет около 70%, что подтверждает, что наши оценочные модели не хуже реальной модели в подборе кривой измеренных данных - следствие сильного эффекта нарушения порядка окружающей температуры.
Для получения дополнительной информации об идентификации динамических систем с помощью System Identification Toolbox посетите информационную страницу System Identification Toolbox.