Стеклянный производственный процесс метро

Этот пример показывает линейную образцовую идентификацию стеклянного производственного процесса трубы. Эксперименты и данные обсуждены в:

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)

Существуют ясные трудности иметь дело с высокочастотными компонентами вывода. Это, в сочетании с длительной задержкой, предлагает, чтобы мы десятикратно уменьшили данные четыре (т.е. фильтр нижних частот это, и выберите каждое четвертое значение):

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.