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

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

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

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

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

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