Этот пример показывает линейную идентификацию модели стеклянного производственного процесса трубы. Эксперименты и данные обсуждены в:
V. Wertz, Г. Бэстин и М. Хит: Идентификация стеклянного места размещения рисования трубы. Proc. 10-го Конгресса IFAC, Vol 10, стр Бумага 334-339 номер 14.5-5-2. Мюнхенский август 1987.
Выход процесса является толщиной и диаметром произведенной трубы. Входные параметры являются давлением воздуха в трубе и скорости рисования.
Проблема моделирования процесса от входной скорости до выходной толщины описана ниже. Различные варианты для анализа данных и определения порядка модели обсуждены.
Мы начинаем путем загрузки входных и выходных данных, сохраненных как объект iddata:
load thispe25.mat
Данные содержатся в переменной glass
:
glass
glass = Time domain data set with 2700 samples. Sample time: 1 seconds Outputs Unit (if specified) Thickn Inputs Unit (if specified) Speed
Данные имеют 2 700 выборок одного входа (Скорость) и один выход (Thickn). Шаг расчета составляет 1 секунду.
Для цели оценки и перекрестной проверки, разделение это в две половины:
ze = glass(1001:1500); %Estimation data zv = glass(1501:2000,:); %Validation data
Представление крупным планом данных об оценке:
plot(ze(101:200)) %Plot the estimation data range from samples 101 to 200.
Давайте удалим средние значения как первый шаг предварительной обработки:
ze = detrend(ze); zv = detrend(zv);
Шаг расчета данных составляет 1 секунду, в то время как постоянные времени процесса могут быть намного медленнее. Мы можем обнаружить некоторые довольно высокие частоты в выходе. Для того, чтобы подтвердить это, давайте сначала вычислим спектры ввода и вывода:
sy = spa(ze(:,1,[])); su = spa(ze(:,[],1)); clf spectrum(sy,su) axis([0.024 10 -5 20]) legend({'Output','Input'}) grid on
Обратите внимание на то, что вход имеет очень мало относительной энергии выше 1 рад/секунда, в то время как выход содержит относительно большие значения выше той частоты. Существуют таким образом некоторые высокочастотные воздействия, которые могут вызвать некоторую проблему для построения моделей.
Мы вычисляем импульсную характеристику, с помощью части данных, чтобы получить некоторое понимание потенциальной обратной связи и задержки от входа, чтобы вывести:
Imp = impulseest(ze,[],'negative',impulseestOptions('RegularizationKernel','SE')); showConfidence(impulseplot(Imp,-10:30),3) grid on
Мы видим задержку приблизительно 12 выборок в импульсной характеристике (сначала значительное значение ответа вне доверительного интервала), который является довольно существенным. Кроме того, импульсная характеристика является весьма значительной для отрицательных задержек. Это указывает, что существует хорошая вероятность обратной связи в данных, так, чтобы будущие значения выходного влияния (добавляются к), текущие входные параметры. Входная задержка может быть вычислена явным образом с помощью delayest
:
delayest(ze)
ans = 12
Вероятность обратной связи может быть получена с помощью feedback
:
feedback(ze) %compute probability of feedback in data
ans = 100
Таким образом почти бесспорно, что существует обратная связь, существующая в данных.
Мы также, как предварительный тест, вычисляем оценку спектрального анализа:
g = spa(ze);
showConfidence(bodeplot(g))
grid on
Мы отмечаем, среди прочего, что высокочастотное поведение довольно сомнительно. Может быть желательно ограничить область значений модели частотами ниже, чем 1 рад/с.
Давайте сделаем быструю проверку, если мы можем взять хорошую динамику, только вычислив четвертый порядок модель ARX с помощью данных об оценке и симулировать ту модель с помощью данных о валидации. Мы знаем, что задержка составляет приблизительно 12 секунд.
m1 = arx(ze,[4 4 12]); compare(zv,m1);
Близкое представление результатов симуляции:
compare(zv,m1,inf,'Samples',101:200)
Существуют ясные трудности иметь дело с высокочастотными компонентами выхода. Это, в сочетании с длительной задержкой, предлагает, чтобы мы десятикратно уменьшили данные четыре (т.е. lowpass фильтруют его и выбирают каждое четвертое значение):
if exist('resample','file')==2 % Use "resample" command for decimation if Signal Processing Toolbox(TM) % is available. zd = resample(detrend(glass),1,4,20); else % Otherwise, use the slower alternative - "idresamp" zd = idresamp(detrend(glass),4); end zde = zd(1:500); zdv = zd(501:size(zd,'N'));
Давайте найдем хорошую структуру для подкошенных данных. Сначала вычислите импульсную характеристику:
Imp = impulseest(zde);
showConfidence(impulseplot(Imp,200),3)
axis([0 100 -0.05 0.01])
grid on
Мы снова видим, что задержка является приблизительно 3 выборками (который сопоставим с тем, что мы видели выше; 12-секундная задержка с шагом расчета 4 секунд в zde). Давайте теперь попытаемся оценить модель по умолчанию, где порядок автоматически выбран средством оценки.
Mdefault = n4sid(zde); compare(zdv,Mdefault)
Средство оценки выбрало 4-ю модель порядка. Это, кажется, обеспечивает лучшую подгонку, чем это для неподкошенных данных. Давайте теперь систематически оценим, какую структуру модели и порядки мы можем использовать. Сначала мы ищем задержку:
V = arxstruc(zde,zdv,struc(2,2,1:30)); nn = selstruc(V,0)
nn = 1×3
2 2 3
ARXSTRUC также предлагает задержку 3 выборок, которая сопоставима с наблюдениями от импульсной характеристики. Поэтому мы фиксируем задержку с близостью 3 и тестируем несколько различных порядков с и вокруг этой задержки:
V = arxstruc(zde,zdv,struc(1:5,1:5,nn(3)-1:nn(3)+1));
Теперь мы вызываем selstruc
на возвращенной матрице для того, чтобы выбрать самый предпочтительный порядок модели (функция с минимальными потерями, которую показывают в первой строке V).
nn = selstruc(V,0); %choose the "best" model order
SELSTRUC мог быть вызван всего одним входом, чтобы вызвать интерактивный режим выбора порядка (nn = selstruc(V)
).
Давайте вычислим и давайте проверять модель на "лучший" порядок, возвращенный в переменной nn
:
m2 = arx(zde,nn);
compare(zdv,m2,inf,compareOptions('Samples',21:150));
Модель m2
о том же самом как Mdefault
подгонка данных, но использует низшие порядки.
Давайте протестируем остаточные значения:
resid(zdv,m2);
Остаточные значения в области доверительного интервала, указывая, что существенные движущие силы были получены моделью. Что нулевая полюсом схема говорит нам?
clf showConfidence(iopzplot(m2),3) axis([ -1.1898,1.3778,-1.5112,1.5688])
Из нулевого полюсом графика, существует, индикация относительно удалений нулей-полюсов для нескольких пар. Это вызвано тем, что их перекрытие местоположений, в областях уверенности. Это показывает, что мы должны смочь преуспеть с моделями более низкоуровневыми. Попробуйте [1 1 3] модель ARX:
m3 = arx(zde,[1 1 3])
m3 = Discrete-time ARX model: A(z)y(t) = B(z)u(t) + e(t) A(z) = 1 - 0.115 z^-1 B(z) = -0.1788 z^-3 Sample time: 4 seconds Parameterization: Polynomial orders: na=1 nb=1 nk=3 Number of free coefficients: 2 Use "polydata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using ARX on time domain data "zde". Fit to estimation data: 35.07% (prediction focus) FPE: 0.4437, MSE: 0.4384
Симуляция модели m3
сравненный с данными о валидации показывает:
compare(zdv,Mdefault,m2,m3)
Эти три модели поставляют сопоставимые результаты. Точно так же мы можем сравнить с 5 шагами вперед поддержка прогноза моделей:
compare(zdv,Mdefault,m2,m3,5)
Как эти графики показывают, сокращение порядка модели не значительно уменьшает свою эффективность, предсказывает будущие значения.
Для получения дополнительной информации об идентификации динамических систем с System Identification Toolbox� посещают страницу информации о продукте System Identification Toolbox.