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

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

Данные имеют 2700 выборки одного входа (Speed) и одного выхода (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

ARXSTRING также предлагает задержку в 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

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

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