В этом примере показано, как оценить передаточную функцию по измеренным данным сигнала.
В этом примере мы оцениваем передаточную функцию для теплообменника. Теплообменник состоит из температуры хладагента, температуры продукта и возмущения температуры окружающей среды. Мы оценим функцию передачи температуры хладагента в продукт.
Измеренные данные хранятся в файле 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
compare и resid команды позволяют исследовать, насколько хорошо оцененная модель соответствует измеренным данным.
set(gcf,'DefaultAxesTitleFontSizeMultiplier',1,... 'DefaultAxesTitleFontWeight','normal',... 'Position',[100 100 780 520]); resid(sysTF,data);
clf compare(data,sysTF)
Рисунок показывает, что остатки сильно коррелируют, подразумевая, что в измеренных данных имеется информация, которая не была адекватно захвачена оценочной моделью.
Ранее мы оценивали передаточную функцию из данных, но отдельно для системного заказа не включали много apriori знаний. Из физики проблемы мы знаем, что система стабильна и имеет положительный выигрыш. Проверяя измеренные данные, мы также знаем, что задержка ввода/вывода составляет около 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», который мы оценили, хранится как числитель и знаменатель в свойстве «NoureTF» модели.
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%, что подтверждает, что наши расчетные модели не хуже реальной модели в подгонке измеренных данных - следствие сильного эффекта возмущения температуры окружающей среды.
Дополнительные сведения об идентификации динамических систем с помощью инструментария идентификации систем см. на информационной странице инструментария идентификации систем.