Этот пример показывает, как оценить передаточную функцию от измеренных данных сигнала.
В этом примере мы оцениваем передаточную функцию для теплообменника. Теплообменник состоит из температуры хладагента, температуры продукта и температуры окружающей среды воздействия. Мы оценим хладагент к передаточной функции температуры продукта.
Результаты измерений хранятся в файле 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', который мы оценили, хранится как числитель и значения знаменателя в свойстве "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.