Exponenta Banner
exponenta event banner

Прогнозирование и прогнозирование временных рядов для прогноза

Открыть сценарий

В этом примере показано, как создать модель временных рядов и использовать ее для прогнозирования, прогнозирования и оценки состояния. Измеренные данные поступают из индукционной печи, размер щели которой со временем уменьшается. Размер щели не может быть измерен непосредственно, но измеряется ток печи и потребляемая мощность. Известно, что по мере увеличения размера щели сопротивление щели уменьшается. Таким образом, отношение измеренного тока в квадрате к измеренной мощности пропорционально размеру слота. Измеренное отношение ток-мощность (измерения тока и мощности являются шумными) используется для создания модели временных рядов и модели для оценки текущего размера слота и прогнозирования будущего размера слота. При физическом осмотре размер щели индукционной печи известен в некоторые моменты времени.

Нагрузка и график измеренных данных

Данные измеренного отношения ток-мощность хранятся в iddata_TimeSeriesPrediction Файл MATLAB. Данные измеряются через почасовые интервалы и показывают, что с течением времени отношение увеличивается, указывая на эрозию щели печи. На основе этих данных разрабатывается модель временных рядов. Начните с разделения данных на сегмент идентификации и проверки.

load iddata_TimeSeriesPrediction
n = numel(y);
ns = floor(n/2);
y_id = y(1:ns,:);
y_v = y((ns+1:end),:);
data_id = iddata(y_id, [], Ts, 'TimeUnit', 'hours');
data_v  = iddata(y_v, [], Ts, 'TimeUnit', 'hours', 'Tstart', ns+1);

plot(data_id,data_v)
legend('Identification data','Validation data','location','SouthEast');

Идентификация модели

Эрозия щели может быть смоделирована как система «состояние-пространство» с входным шумом и измеренным отношением ток-мощность в качестве выходного сигнала. Измеренное отношение ток-мощность пропорционально состоянию системы, или

$x_{n+1} = Ax_n + Ke_n$

$y_n = Cx_n + e_n$

Где$x_n$ вектор состояния содержит размер слота;$y_n$ - измеренное отношение ток-мощность;$e_n$ шум и$A, C, K$ должны быть идентифицированы.

Используйте ssest() для идентификации модели дискретного состояния-пространства по измеренным данным.

sys = ssest(data_id,1,'Ts',Ts,'form','canonical')

sys =
  Discrete-time identified state-space model:
    x(t+Ts) = A x(t) + K e(t)
       y(t) = C x(t) + e(t)
 
  A = 
          x1
   x1  1.001
 
  C = 
       x1
   y1   1
 
  K = 
            y1
   x1  0.09465
 
Sample time: 1 hours
  
Parameterization:
   CANONICAL form with indices: 1.
   Disturbance component: estimate
   Number of free coefficients: 2
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                             
Estimated using SSEST on time domain data "data_id".
Fit to estimation data: 67.38% (prediction focus)   
FPE: 0.09575, MSE: 0.09348

Идентифицированная модель минимизирует прогноз на 1 шаг вперед. Проверьте модель с использованием предиктора на 10 шагов вперед, то есть с использованием$y_0,...,y_n$ модели для прогнозирования. $y_{n+10}$Следует отметить, что погрешность между измеренными и прогнозируемыми значениями используется $y_0-\hat{y_0},...,y_n-\hat{y_n}$для прогнозирования$y_{n+10}$.

Используйте предиктор на 10 шагов вперед для идентификационных данных и независимых данных проверки.

nstep = 10;
compare(sys,data_id,nstep)  % comparison of 10-step prediction to estimation data
grid('on');


figure; compare(sys,data_v,nstep)  % comparison to validation data
grid('on');

Вышеприведенное упражнение Оба набора данных показывают, что предиктор соответствует измеренным данным.

Прогнозирование используется для дальнейшей проверки модели. Прогнозирование использует запись измеренных данных$y_0,y_1,...,y_n-\hat{y_n}$ для вычисления состояния модели на временном шаге N. Это значение используется в качестве начального условия для прогнозирования реакции модели на будущий временной промежуток. Мы прогнозируем отклик модели на временной промежуток данных проверки, а затем сравниваем их. Мы также можем вычислить неопределенность в прогнозах и построить график +/- 3 sd их значений.

MeasuredData = iddata(y, [], Ts, 'TimeUnit', 'hours'); % = [data_id;data_v]
t0 = MeasuredData.SamplingInstants;

Horizon = size(data_v,1); % forecasting horizon
[yF, ~, ~, yFSD]  = forecast(sys, data_id, Horizon);
% Note: yF is IDDATA object while yFSD is a double vector
t = yF.SamplingInstants; % extract time samples
yFData = yF.OutputData;      % extract response as double vector
plot(MeasuredData)
hold on
plot(t, yFData, 'r.-', t, yFData+3*yFSD, 'r--', t, yFData-3*yFSD, 'r--')
hold off
title('Forecasted response over the validation data''s time span')
grid on

График показывает, что реакция модели с доверительными интервалами (обозначенными красными пунктирными кривыми) перекрывает измеренное значение для данных проверки. Объединенные результаты прогнозирования и прогнозирования показывают, что модель представляет измеренное отношение ток-мощность.

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

Наконец, мы используем модель для прогнозирования ответа 200 шагов в будущем на период времени в 502-701 часов.

Horizon = 200; % forecasting horizon
[yFuture, ~, ~, yFutureSD] = forecast(sys, MeasuredData, Horizon);
t = yFuture.SamplingInstants; % extract time samples
yFutureData = yFuture.OutputData;      % extract response as double vector
plot(t0, y,...
   t, yFutureData, 'r.-', ...
   t, yFutureData+3*yFutureSD, 'r--', ...
   t, yFutureData-3*yFutureSD, 'r--')
title('Forecasted response (200 steps)')
grid on

Синяя кривая показывает измеренные данные, охватывающие более 1-501 часов. Красная кривая - это прогнозируемый отклик в течение 200 часов за пределами временного диапазона измеряемых данных. Красные пунктирные кривые показывают неопределенность 3 sd в прогнозируемом ответе на основе случайной выборки идентифицированной модели.

Оценка состояния

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

Создайте предиктор для произвольного состояния. Идентифицированные ковариации модели должны быть преобразованы в модель предиктора с использованием translatecov() команда. createPredictor() функция просто извлекает третий выходной аргумент predict() функция для использования с translatecov().

type createPredictor
est = translatecov(@(s) createPredictor(s,data_id),sys)

function pred = createPredictor(mdl,data)
%CREATEPREDICTOR Return 1-step ahead predictor.
%
%   sys = createPredictor(mdl,data)
%
%   Create a 1-step ahead predictor model sys for the specified model mdl
%   and measured data. The function is used by
%   |TimeSeriedPredictionExample| and the |translatecov()| command to
%   translate the identified model covariance to the predictor.

% Copyright 2015 The MathWorks, Inc.
[~,~,pred] = predict(mdl,data,1);

est =
  Discrete-time identified state-space model:
    x(t+Ts) = A x(t) + B u(t)
       y(t) = C x(t) + D u(t)
 
  A = 
           x1
   x1  0.9064
 
  B = 
            y1
   x1  0.09465
 
  C = 
       x1
   y1   1
 
  D = 
       y1
   y1   0
 
Sample time: 1 hours
  
Parameterization:
   CANONICAL form with indices: 1.
   Feedthrough: none
   Disturbance component: none
   Number of free coefficients: 2
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                         
Created by direct construction or transformation. Not estimated.

Модель est - предиктор на 1 шаг вперед, выраженный в тех же координатах состояния, что и исходная модель; sys. Как преобразовать координаты состояния таким образом, чтобы состояние модели соответствовало (зависящему от времени) размеру слота? Решение состоит в том, чтобы полагаться на фактические прямые измерения размера слота, выполняемые с перерывами. Это не редкость на практике, когда затраты на проведение прямых измерений высоки и проводятся только периодически (например, при замене компонента).

В частности, преобразуют предикторное состояние,, $x_n$в, $z_n$где$y_n = Cz_n$$y_n$ измеренное отношение ток-мощность, и является$z_n$ размером щели печи. В этом примере четыре прямых измерения размера щели печи, sizeMeasuredи отношение ток-мощность печи, ySizeMeasured, используются для оценки. $C$При преобразовании предиктора ковариации предиктора также должны быть преобразованы. Следовательно, мы используем translatecov() для выполнения преобразования координат состояния.

Cnew = sizeMeasured\ySizeMeasured;
est  = translatecov(@(s) ss2ss(s,s.C/Cnew),est)

est =
  Discrete-time identified state-space model:
    x(t+Ts) = A x(t) + B u(t)
       y(t) = C x(t) + D u(t)
 
  A = 
           x1
   x1  0.9064
 
  B = 
           y1
   x1  0.9452
 
  C = 
           x1
   y1  0.1001
 
  D = 
       y1
   y1   0
 
Sample time: 1 hours
  
Parameterization:
   CANONICAL form with indices: 1.
   Feedthrough: none
   Disturbance component: none
   Number of free coefficients: 2
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                         
Created by direct construction or transformation. Not estimated.

Теперь предиктор выражается в требуемых координатах состояния. Он имеет один вход, который является измеренным выходом системы (отношение тока к мощности печи), и один выход, который является прогнозируемым выходом системы (размер щели печи). Предиктор моделируется для оценки выходного сигнала системы и состояния системы.

opts = simOptions;
opts.InitialCondition = sizeMeasured(1);
U = iddata([],[data_id.Y; data_v.Y],Ts,'TimeUnit','hours');
[ye,ye_sd,xe] = sim(est,U,opts);

Сравните расчетный выходной сигнал и размер слота с измеренными и известными значениями.

yesdp   = ye;
yesdp.Y = ye.Y+3*ye_sd;
yesdn   = ye;
yesdn.Y = ye.Y-3*ye_sd;
n = numel(xe);
figure, plot([data_id;data_v],ye,yesdp,'g',yesdn,'g')
legend('Measured output','Estimated output','99.7% bound','location','SouthEast')
grid('on')
figure, plot(tSizeMeasured,sizeMeasured,'r*',1:n,xe,1:n,yesdp.Y/est.C,'g',1:n,yesdn.Y/est.C,'g');
legend('Measured state','Estimated state','99.7% bound','location','SouthEast')
xlabel('Time (hours)')
ylabel('Amplitude');
grid('on')

Использование прогнозирования и прогнозирования для прогноза

Комбинация модели предиктора и прогнозирования позволяет выполнить прогноз на индукционной печи.

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

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

См. также

| |

Связанные темы

Документация по инструментам идентификации системы

Поддержка