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

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

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.

Figure contains 2 axes objects. Axes object 1 with title Thickn contains an object of type line. This object represents untitled1. Axes object 2 with title Speed contains an object of type line. This object represents untitled1.

Предварительный анализ данных

Давайте удалим средние значения как первый шаг предварительной обработки:

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

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent Output, Input.

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

Мы вычисляем импульсную характеристику, с помощью части данных, чтобы получить некоторое понимание потенциальной обратной связи и задержки от входа, чтобы вывести:

Imp = impulseest(ze,[],'negative',impulseestOptions('RegularizationKernel','SE'));
showConfidence(impulseplot(Imp,-10:30),3)
grid on

Figure contains an axes object. The axes object with title From: Speed To: Thickn contains 2 objects of type line. This object represents Imp.

Мы видим задержку приблизительно 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

Figure contains 2 axes objects. Axes object 1 with title From: Speed To: Thickn contains an object of type line. This object represents g. Axes object 2 contains an object of type line. This object represents g.

Мы отмечаем, среди прочего, что высокочастотное поведение довольно сомнительно. Может быть желательно ограничить область значений модели частотами ниже, чем 1 рад/с.

Параметрические модели поведения процесса

Давайте сделаем быструю проверку, если мы можем взять хорошую динамику, только вычислив четвертый порядок модель ARX с помощью данных об оценке и симулировать ту модель с помощью данных о валидации. Мы знаем, что задержка составляет приблизительно 12 секунд.

m1 = arx(ze,[4 4 12]);
compare(zv,m1);

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent zv (Thickn), m1: 11.14%.

Близкое представление результатов симуляции:

compare(zv,m1,inf,'Samples',101:200)

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent zv (Thickn), m1: 12.62%.

Существуют ясные трудности иметь дело с высокочастотными компонентами выхода. Это, в сочетании с длительной задержкой, предлагает, чтобы мы десятикратно уменьшили данные четыре (i.e. 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

Figure contains an axes object. The axes object with title From: Speed To: Thickn contains 2 objects of type line. This object represents Imp.

Мы снова видим, что задержка является приблизительно 3 выборками (который сопоставим с тем, что мы видели выше; 12-секундная задержка с шагом расчета 4 секунд в zde). Давайте теперь попытаемся оценить модель по умолчанию, где порядок автоматически выбран средством оценки.

Mdefault = n4sid(zde);
compare(zdv,Mdefault)

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent zdv (Thickn), Mdefault: 32.73%.

Средство оценки выбрало 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));

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent zdv (Thickn), m2: 33.44%.

Модель m2 о том же самом как Mdefault подгонка данных, но использует низшие порядки.

Давайте протестируем остаточные значения:

resid(zdv,m2);

Figure contains 2 axes objects. Axes object 1 with title AutoCorr contains 2 objects of type line. This object represents m2. Axes object 2 with title XCorr (Speed) contains 2 objects of type line. This object represents m2.

Остаточные значения в области доверительного интервала, указывая, что существенные движущие силы были получены моделью. Что нулевая полюсом схема говорит нам?

clf
showConfidence(iopzplot(m2),3)
axis([ -1.1898,1.3778,-1.5112,1.5688])

Figure contains an axes object. The axes object with title From: Speed To: Thickn contains 4 objects of type line. This object represents m2.

От диаграммы нулей и полюсов, существует, индикация относительно удалений нулей-полюсов для нескольких пар. Это вызвано тем, что их перекрытие местоположений, в областях доверия. Это показывает, что мы должны смочь преуспеть с моделями более низкоуровневыми. Попробуйте [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)

Figure contains an axes object. The axes object contains 4 objects of type line. These objects represent zdv (Thickn), Mdefault: 32.73%, m2: 34.77%, m3: 33.04%.

Эти три модели поставляют сопоставимые результаты. Точно так же мы можем сравнить с 5 шагами вперед поддержка предсказания моделей:

compare(zdv,Mdefault,m2,m3,5)

Figure contains an axes object. The axes object contains 4 objects of type line. These objects represent zdv (Thickn), Mdefault: 31.91%, m2: 34.74%, m3: 33.04%.

Как эти графики показывают, сокращение порядка модели не значительно уменьшает свою эффективность, предсказывает будущие значения.

Для просмотра документации необходимо авторизоваться на сайте