Экспериментальные данные

Мы начинаем путем загрузки входных и выходных данных, сохраненных как объект 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)

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

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

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