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

Данные об измеренном отношении ток-мощность хранятся в 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');

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

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

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

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

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

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

Создайте предиктор для произвольного состояния. Идентифицированные ковариации модели должны быть переведены в модель предиктора с помощью translatecov() команда. The 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')

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

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

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

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