Прогноз временных рядов и предсказывающий для прогноза

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

Загрузите и отобразите результаты измерений на графике

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

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');

Образцовая идентификация

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

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

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

Используйте эти 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');

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

Прогнозирование используется, чтобы далее проверить модель. Прогнозирование использует запись результатов измерений, чтобы вычислить образцовое состояние на временном шаге 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. Как мы преобразовываем координаты состояния так, чтобы состояние модели соответствовало (зависящему от времени) размеру слота? Решение состоит в том, чтобы полагаться на фактические, прямые измерения размера слота, взятого периодически. Это весьма распространено на практике, где стоимость проведения прямых измерений высока и только периодически делаться (такой как тогда, когда компонент заменяется).

А именно, преобразуйте состояние предиктора, к, так, чтобы, где измеренное отношение текущей степени, и печной размер слота. В этом примере четыре прямых измерения печного размера слота, sizeMeasured, и печного отношения текущей степени, ySizeMeasured, используются, чтобы оценить. В преобразовании предиктора также должны быть преобразованы ковариации предиктора. Следовательно мы используем команду 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')

Используя прогноз и предсказывающий для прогноза

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

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

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

Смотрите также

| |

Похожие темы