Exponenta 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 шагов вперед предиктор, i.e., учитывая$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')

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

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

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

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

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

| |

Похожие темы

Документация System Identification Toolbox

Поддержка