Введение в прогнозирование Динамической системы отклика

Прогнозирование отклика динамической системы является предсказанием будущих выходов системы с помощью прошлых выходных измерений. Другими словами, заданные наблюдения y (<reservedrangesplaceholder5>) = {y(1)..., y(N)} выхода системы, прогнозирование - предсказание выходов y(N+1)..., y(N+H) до будущего периода времени H.

Когда вы выполняете прогнозирование в программном обеспечении System Identification Toolbox™, вы сначала идентифицируете модель, которая подходит для прошлых измеренных данных из системы. Модель может быть линейной моделью временных рядов, такой как AR, ARMA и модели пространства состояний, или нелинейной моделью ARX. Если экзогенные входы влияют на выходы системы, можно выполнить прогнозирование с помощью моделей вход-выход, таких как ARX и ARMAX. После идентификации модели вы затем используете forecast команда для вычисления y(N+1),..., y(N+H). Команда вычисляет прогнозируемые значения следующим образом:

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

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

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

Тулбокс позволяет вам прогнозировать данные временных рядов (только для вывода) с помощью линейных моделей, таких как AR, ARMA и модели пространства состояний. Вот рисунок прогнозирования отклика авторегрессивной модели, за которым следуют шаги прогнозирования для более сложных моделей, таких как модели скользящего среднего и модели пространства состояний.

Авторегрессионные модели

Предположим, что вы собрали данные временных рядов y (t) = {y(1),..., y(N)} стационарного случайного процесса. Если предположить, что данные являются авторегрессионным процессом второго порядка (AR), можно описать динамику следующей моделью AR:

$y (t) + a_{1} y (t - 1) + a_{2} y (t - 2) = e (t)$

Где _a1 и _a2 - коэффициенты аппроксимации, а e (t) - термин шума.

Можно идентифицировать модель, используя ar команда. Программа вычисляет коэффициенты соответствия и отклонение e (t) путем минимизации 1-шаговых ошибок предсказания между наблюдениями {y(1),..., y(N)} и ответом модели.

Принимая, что e инноваций (t) являются нулевой средней белой последовательностью, можно вычислить предсказанный выход $\hat{y} (t)$ использование формулы:

$\hat{y} (t) = - a_{1} y (t - 1) - a_{2} y (t - 2)$

Где y(t-1) и y(t-2) являются либо измеренными данными, если они доступны, либо ранее предсказанными значениями. Для примера прогнозируемые выходы пяти шагов в будущем:

$\begin{array}{l} \hat{y} (N + 1) = - a_{1} y (N) - a_{2} y (N - 1) \\ \hat{y} (N + 2) = - a_{1} \hat{y} (N + 1) - a_{2} y (N) \\ \hat{y} (N + 3) = - a_{1} \hat{y} (N + 2) - a_{2} \hat{y} (N + 1) \\ \hat{y} (N + 4) = - a_{1} \hat{y} (N + 3) - a_{2} \hat{y} (N + 2) \\ \hat{y} (N + 5) = - a_{1} \hat{y} (N + 4) - a_{2} \hat{y} (N + 3) \end{array}$

Обратите внимание, что расчет $\hat{y} (N + 2)$ использует ранее предсказанное значение $\hat{y} (N + 1)$ поскольку измеренные данные недоступны после временного шага N. Таким образом, прямой вклад измеренных данных уменьшается, когда вы прогнозируете дальше в будущее.

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

Модели скользящего среднего

В моделях скользящего среднего (MA) выход зависит от текущих и прошлых инноваций (e (t), e (t-1), e (t-2), e (t-3)....). Таким образом, прогнозирование отклика моделей MA требует знания начальных условий измеренных данных.

Предположим, что данные временных рядов y (t) из вашей системы могут соответствовать модели скользящего среднего второго порядка:

$y (t) = e (t) + c_{1} e (t - 1) + c_{2} e (t - 2)$

Предположим, что y(1) и y(2) являются единственными доступными наблюдениями и их значения равны 5 и 10, соответственно. Можно оценить коэффициенты модели _c1 и _c2, используя armax команда. Предположим, что оценочные _c1 и _c2 значения равны 0,1 и 0,2 соответственно. Затем, принимая, как и прежде, что e (t) является случайной переменной с нулевым средним значением, можно предсказать выходное значение в то время t используя следующую формулу:

$\hat{y} (t) = c_{1} e (t - 1) + c_{2} e (t - 2)$

Где e(t-1) и e(t-2) являются ли различия между измеренной и предсказанной характеристиками в разы t-1 и t-2, соответственно. Если измеренных данных для этих раз не существует, используется нулевое значение, потому что процесс инноваций e (t) принимается как нулевой белый Гауссов шум.

Поэтому прогнозируемый выход в момент времени t = 3 равен:

$\hat{y} (3) = 0.1 e (2) + 0.2 e (1)$

Где, нововведения e(1) и e(2) являются ли различие между наблюдаемым и прогнозируемым значениями выхода в момент времени t равной 1 и 2, соответственно:

$\begin{array}{l} e (2) = y (2) - \hat{y} (2) = y (2) - [0.1 e (1) + 0.2 e (0)] \\ e (1) = y (1) - \hat{y} (1) = y (1) - [0.1 e (0) + 0.2 e (- 1)] \end{array}$

Потому что данные измерялись со времени t равными 1, значения e(0) и e(-1) неизвестны. Таким образом, чтобы вычислить прогнозируемые выходы, значение этих начальных условий e(0) и e(-1) обязательно. Можно либо принять нулевые начальные условия, либо оценить их.

Нулевые начальные условия: Если вы задаете, что e(0) и e(-1) равны 0, значения ошибок и прогнозируемые выходы:
$\begin{array}{l} e (1) = 5 - (0.1 * 0 + 0.2 * 0) = 5 \\ e (2) = 10 - (0.1 * 5 + 0.2 * 0) = 9.5 \\ \hat{y} (3) = 0.1 * 9.5 + 0.2 * 5 = 1.95 \end{array}$
Прогнозные значения в моменты времени t = 4 и 5:
$\begin{array}{l} \hat{y} (4) = 0.1 e (3) + 0.2 e (2) \\ \hat{y} (5) = 0.1 e (4) + 0.2 e (3) \end{array}$
Вот e(3) и e(4) приняты равными нулю, так как нет измерений после времени t = 2. Это предположение дает, $\hat{y} (4) = 0.2 * e (2) = 0.2 * 9.5 = 1.9$ , и $\hat{y} (5) = 0$ .
Таким образом, для этой модели MA второго порядка, предсказанные выходы, которые находятся более чем на два временных шагов за последней измеренной точкой данных (t = 2), все равны нулю. В целом, когда приняты нулевые начальные условия, прогнозируемые значения сверх порядка чистой модели MA без авторегрессивных членов все равны нулю.
Предполагаемые начальные условия: Можно оценить начальные условия, минимизировав квадратную сумму 1-шаговых ошибок предсказания всех измеренных данных.
Для модели MA, описанной ранее, оценка начальных условий e(0) и e(-1) требует минимизации следующей функции затрат методом наименьших квадратов:
$V = e {(1)}^{2} + e {(2)}^{2} = (y (1 ) - [0 .1 e (0 ) + 0 {2 е (-1)])}^{2} + (y (2 ) - [0 .1 e (1 ) + 0 {2 е (0)])}^{2}$
Подстановка a = e(0) и b = e(-1), функцией стоимости является:
$V (a, b) = {(5 - [0.1 a + 0.2 b])}^{2} + {(10 - [0.1 (5 - [0.1 a + 0.2 b]) + 0.2 a])}^{2}$
Минимизация выражений V e(0) = 50 и e(-1) = 0, что дает:
$\begin{array}{l} e (1) = 5 - (0.1 * 50 + 0.2 * 0) = 0 \\ e (2) = 10 - (0.1 * 0 + 0.2 * 50) = 0 \\ \hat{y} (3) = 0 \\ \hat{y}(4) = 0 \end{array}$
Таким образом, для этой системы, если ошибки предсказания минимизированы за доступные две выборки, все будущие предсказания равны нулю, что является средним значением процесса. Если бы было доступно более двух наблюдений, вы бы оценили e(-1) и e(0) использование подхода методом наименьших квадратов для минимизации 1-шаговых ошибок предсказания по всем доступным данным.
В этом примере показано, как воспроизвести эти прогнозные результаты с помощью forecast команда.
Загрузите измеренные данные.
```
PastData = [5;10];
```
Создайте модель MA с полиномиальными коэффициентами A и C, равными 1 и [1 0,1 0,2], соответственно.
```
model = idpoly(1,[],[1 0.1 0.2]);
```
Задайте нулевые начальные условия и спрогнозируйте выход на пять шагов в будущее.
```
opt = forecastOptions('InitialCondition','z');
yf_zeroIC = forecast(model,PastData,5,opt)
```
```
yf_zeroIC = 5×1

    1.9500
    1.9000
         0
         0
         0
```
Укажите, что программное обеспечение оценивает начальные условия и прогнозирует выход.
```
opt = forecastOptions('InitialCondition','e');
yf_estimatedIC = forecast(model,PastData,5,opt)
```
```
yf_estimatedIC = 5×1
10^-15 ×

   -0.3553
   -0.3553
         0
         0
         0
```

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

Модели в пространстве состояний

Модель дискретного времени пространства состояний данных временных рядов имеет вид:

$\begin{array}{l} x (t + 1) = A x (t) + K e (t) \\ y (t) = C x (t) + e (t) \end{array}$

Где, x (t) является вектором состояния, y (t) являются выходами, e (t) является шумом. A, C и K являются матрицами пространства состояний с фиксированным коэффициентом.

Можно представлять любую произвольную структуру линейной модели в форме пространства состояний. Для примера может быть показано, что модель ARMA, описанная ранее, выражена в форме пространства состояний, используя A = [0 0;1 0], K = [0.5;0] и C = [0.2 0.4]. Можно оценить модель пространства состояний из наблюдаемых данных с помощью таких команд, как ssest и n4sid. Можно также преобразовать существующую полиномиальную модель, такую как AR, MA, ARMA, ARX и ARMAX, в форму пространства состояний, используя idss команда.

Преимущество формы пространства состояний заключается в том, что любая авторегрессивная или скользящая модель со многими терминами временной задержки (t-1, t-2, t-3,...) имеет только одну временную задержку (t-1) в переменных состояния, когда модель преобразуется в форму пространства состояний. В результате необходимые начальные условия для прогнозирования переводятся в одно значение для вектора начального состояния X(0). forecast команда преобразует всю линейную модель в форму пространства состояний, а затем выполняет прогнозирование.

Чтобы предсказать ответ модели пространства состояний:

Сгенерируйте модель предиктора на 1 шаг вперед для идентифицированной модели. Модель предиктора имеет вид:
$\begin{array}{l} \hat{x} (t + 1) = (A - K * C) \hat{x} (t) + K y (t) \\ \hat{y} (t) = C * \hat{x} (t) \end{array}$
Где y (t) является измеренным выходом и $\hat{y} (t)$ - предсказанное значение. Измеренный выход доступен до временного шага N и используется как вход в модели предиктора. Начальный вектор состояния $\hat{x} (0) = x_{0}$ .
Присвойте значение вектору начального состояния _x0.
Начальные состояния либо заданы как нули, либо оценены путем минимизации ошибки предсказания за измеренный временной интервал данных.
Задайте нулевое начальное условие, если система находилась в состоянии покоя до сбора наблюдений. Можно также задать нулевые начальные условия, если модель предиктора достаточно стабильна, потому что стабильность подразумевает, что эффект начальных условий быстро уменьшается, когда наблюдения собираются. Модель предиктора стабильна, если собственные значения A-K*C находятся внутри модуля круга.
Вычислить $\hat{x} (N + 1)$ , значение состояний в момент времени t = N+1, момент времени после последней доступной выборки данных.
Для этого моделируйте модель предиктора, используя измеренные наблюдения в качестве входов:
$\begin{array}{l} \hat{x} (1) = (A - K * C) x_{0} + K y (0) \\ \hat{x} (2) = (A - K * C) \hat{x} (1) + K y (1) \\ ⋮ \\ \hat{x} (N + 1) = (A - K * C) \hat{x} (N) + K y (N) \end{array}$
Симулируйте ответ идентифицированной модели для H шаги с использованием $\hat{x} (N + 1)$ в качестве начальных условий, где H - горизонт предсказания. Этот ответ является прогнозируемым ответом модели.

Воспроизведение выхода `forecast` Команда

Открыть Live Script

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

Загрузка данных временных рядов.

load iddata9 z9

z9 является iddata объект, который хранит данные временных рядов (без входов).

Задайте данные, которые будут использоваться для оценки модели.

observed_data = z9(1:128); 
Ts = observed_data.Ts;
t = observed_data.SamplingInstants;
y = observed_data.y;

Ts - шаг расчета измеренных данных, t является временной вектор и y - вектор измеренных данных.

Оцените модель пространства состояний в дискретном времени 4-го порядка.

sys = ssest(observed_data,4,'Ts',Ts);

Прогнозируйте выход модели пространства состояний на 100 шагов в будущее с помощью forecast команда.

H = 100;
yh1 = forecast(sys,observed_data,H);

yh1 - прогнозируемый выход, полученный с помощью forecast команда. Теперь воспроизведите выход путем ручной реализации алгоритма, используемого forecast команда.

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

A = sys.A;
K = sys.K;
C = sys.C;

Сгенерируйте предиктор на 1 шаг вперед, где A матрица Predictor модель A-K*C и B матрица K.

Predictor = idss((A-K*C),K,C,0,'Ts',Ts);

Оцените начальные состояния, которые минимизируют различие между наблюдаемым выходным y и 1-ступенчатый предсказанный ответ идентифицированной модели sys.

x0 = findstates(sys,observed_data,1);

Распространите вектор состояния до конца наблюдаемых данных. Для этого моделируйте предиктор с помощью y как вход и x0 как начальные состояния.

Input = iddata([],y,Ts);
opt = simOptions('InitialCondition',x0);
[~,~,x] = sim(Predictor,Input,opt);
xfinal = x(end,:)';

xfinal - значение вектора состояния в момент времени t(end), последний раз, когда наблюдаемые данные доступны. Прогнозирование 100 временных шагов в будущее начинается на следующем временном шаге, t1 = t(end)+Ts.

Чтобы реализовать алгоритм прогнозирования, значение вектора состояния в момент времени t1 обязательно. Вычислите вектор состояния путем применения уравнения обновления состояния Predictor модель в xfinal.

x0_for_forecasting = Predictor.A*xfinal + Predictor.B*y(end);

Симулируйте идентифицированную модель для H шаги с использованием x0_for_forecasting в качестве начальных условий.

opt = simOptions('InitialCondition',x0_for_forecasting);

Потому что sys является моделью временных рядов, задайте входы для симуляции как H-by-0 сигнал, где H - желаемое количество выходных выборок симуляции.

Input = iddata([],zeros(H,0),Ts,'Tstart',t(end)+Ts);
yh2 = sim(sys,Input,opt);

Сравните результаты forecast командная yh1 с полученными вручную результатами yh2.

plot(yh1,yh2,'r.')

Figure contains an axes. The axes with title y1 contains 2 objects of type line. These objects represent yh1, yh2.

График показывает, что результаты совпадают.

Прогнозирование отклика линейных моделей с экзогенными входами

Открыть Live Script

Когда существуют экзогенные раздражители, влияющие на систему, система не может считаться стационарной. Однако, если эти стимулы измеримы, вы можете рассматривать их как входы для системы и учитывать их эффекты при прогнозировании выхода системы. Рабочий процесс для прогнозирования данных с экзогенными входными параметрами аналогичен для прогнозирования данных временных рядов. Вы сначала идентифицируете модель, чтобы соответствовать измеренным входно-выходным данным. Затем вы задаете ожидаемые входные значения для временного интервала прогнозирования и прогнозируете выход идентифицированной модели, используя forecast команда. Если вы не задаете ожидаемые входные значения, они приняты равными нулю.

В этом примере показано, как спрогнозировать модель ARMAX с экзогенными входами в тулбоксе:

Загрузите входно-выходные данные.

load iddata1 z1

z1 является iddata объект с данными ввода-вывода в 300 временных точках.

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

past_data = z1(1:150);

Идентифицируйте модель ARMAX Ay (t) = Bu (t-1) + Ce (t) порядка [2 2 2 1].

na = 2; % A polynomial order
nb = 2; % B polynomial order
nc = 2; % C polynomial order
nk = 1; % input delay
sys = armax(past_data,[na nb nc nk]);

Прогнозируйте ответ на 100 временных шагов в будущем, после последней выборки наблюдаемых данных past_data. Укажите ожидаемые входы в 100 будущих временных точках.

H = 100;
FutureInputs = z1.u(151:250);
forecast(sys,past_data,H,FutureInputs)
legend('Past Outputs','Future Outputs')

Figure contains an axes. The axes contains 2 objects of type line. These objects represent Past Outputs (y1), Future Outputs.

Прогнозирование отклика нелинейных моделей

Тулбокс также позволяет вам прогнозировать данные с помощью нелинейных моделей ARX, Hammerstein-Wiener и нелинейных серых ящиков.

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

Прогнозирующий ответ нелинейных моделей ARX

Нелинейная модель ARX временных рядов имеет следующую структуру:

$y (t) = f (y (t - 1), y (t - 2), ..., y (t - N)) + e (t)$

Где f - нелинейная функция с входами R (t), регрессоры модели. Регрессорами могут быть задержанные по времени переменные y(t-1), y(t-2),..., y(t-N) и их нелинейные выражения, такие как y (t-1)², y(t-1)y(t-2), abs(y(t-1)). Когда вы оцениваете нелинейную модель ARX из измеренных данных, вы задаете регрессоры модели. Можно также задать структуру f, используя отличные структуры, такие как вейвлеты и древовидные разделы. Для получения дополнительной информации смотрите страницу с описанием для nlarx команда оценки.

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

$R (t) = {[y (t - 1), y (t - 2), y {(t - 1)}^{2}, y {(t - 2)}^{2}, y (t - 1) y (t - 2)]}^{T}$

Затем f(R)=W'*R+c, где _{W=[w1 w2 w3 w4 w5]} является вектором взвешивания, и c - смещение выхода.

Нелинейная модель ARX имеет вид:

$y (t) = w_{1} y (t - 1) + w_{2} y (t - 2) + w_{3} y {(t - 1)}^{2} + w_{4} y {(t - 2)}^{2} + w_{5} y (t - 1) y (t - 2) + c + e (t)$

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

Когда вы используете forecast команда, программное обеспечение вычисляет прогнозируемые выходы модели путем симуляции модели H временных шагов в будущее, используя последние N измеренные выходные выборки в качестве начальных условий. Где N - самая большая задержка в регрессорах, и H - заданный горизонт прогноза.

Для модели линейного регрессора предположим, что вы измерили 100 выборок выхода y, и вы хотите спрогнозировать четыре шага в будущее (H = 4). Самая большая задержка в регрессорах модели - N = 2. Поэтому программное обеспечение берёт последние две выборки данных y(99) и y(100) в качестве начальных условий и прогнозирует выходы как:

$\begin{array}{l} \hat{y} (101) = w_{1} y (100) + w_{2} y (99) + w_{3} y (^{100) 2} + w_{4} y (^{99) 2} + w_{5} y (100) y (99) \\ \hat{y} (102) = w_{1} \hat{y} (101) + w_{2} y (100) + w_{3} \hat{y} (^{101) 2} + w_{4} y (^{100) 2} + w_{5} \hat{y} (101) y (100) \\ \hat{y} (103) = w_{1} \hat{y} (102) + w_{2} \hat{y} (101) + w_{3} \hat{y} (^{102) 2} + w_{4} \hat{y} (^{101) 2} + w_{5} \hat{y} (102) \hat{y} (101) \\ \hat{y} (104) = w_{1} \hat{y} (103) + w_{2} \hat{y} (102) + w_{3} \hat{y} (^{103) 2} + w_{4} \hat{y} (^{102) 2} + w_{5} \hat{y} (103) \hat{y} (102) \end{array}$

Если ваша система имеет экзогенные входы, нелинейная модель ARX также включает регрессоры, которые зависят от входных переменных. Процесс прогнозирования аналогичен процессу для данных временных рядов. Вы сначала идентифицируете модель, sys, используя входно-выходные данные, past_data. Когда вы прогнозируете данные, программное обеспечение моделирует идентифицированную модель H временные шаги в будущем, используя последние N измеренные выходные выборки в качестве начальных условий. Вы также задаете ожидаемые входные значения для временного интервала прогнозирования, FutureInputs. Синтаксис для прогнозирования отклика нелинейных моделей с экзогенными входами такой же, как для линейных моделей, forecast(sys,past_data,H,FutureInputs).

См. также

forecast | predict | sim

Подробнее о

Моделируйте и прогнозируйте идентифицированный выход модели

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