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