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

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

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

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

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

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

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

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

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

$y (t)_{+} a1y (t_{−} 1) + a2y ($ t − 2) = e (t)

Где _а1 и _а2 - коэффициенты соответствия, а е (t) - член шума.

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

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

$\overset{}{y}^(t)_{=} - a1y (_{t} - 1)$ − a2y (t − 2)

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

$\begin{array}{l} \overset{}{y}^(N + 1_{)} = -_{} a1y (N) \\ \overset{}{-} a2y (N -_{} \overset{)}{1} y^(N_{} + 2) \\ \overset{=}{} - a1y^_{} \overset{N}{(} + 1) -_{} \overset{}{} a2y (N) \\ \overset{y}{}^(N +_{3} \overset{}{)} = - {a1y}_{} \overset{}{^} (N + \\ \overset{)}{2} - a2y^_{} \overset{N}{(} + 1) y_{} \overset{}{^} (N + \end{array}$ 4) = − a1y

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

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

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

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

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

$y (t) = e (_{t}) + c1e_{(} t - 1)$ + c2e (t − 2)

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

$\overset{}{y}^(t)_{} = c1e (t_{} - 1) +$ c2e (t − 2)

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

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

$\overset{}{y}^(3) = 0,1e (2) +$ 0,2e (1)

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

$\begin{array}{l} e (2) = y \overset{2}{(}) - y^(2) = y (2) - [0,1e \\ (1) + 0,2 \overset{}{e} (0)] e (1) = y (1) -y^(1) \end{array}$ = y (1) - [0,1e (0) + 0,2 e (-1)]

Поскольку данные измерялись от времени 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) = 5e (2) = 10 - ( \\ \overset{}{} 0,1 * 5 + 0,2 * 0) = 9 . \end{array}$ 5y ^ (3) = 0,1 * 9,5 + 0,2 * 5 = 1,95
Прогнозируемые значения в моменты времени t = 4 и 5:
$\begin{array}{l} \overset{}{y}^(4) = 0,1e (3) + \\ \overset{}{} 0,2e (2) y^(5) = \end{array}$ 0,1e (4) + 0,2e (3)
Здесь e(3) и e(4) принимают за ноль, так как нет измерений за пределами времени t = 2. Это предположение дает $\overset{}{y}^(4) = 0,2 * e (2) =$ 0,2 * $\overset{}{9,5} = 1,9$ , и y ^ (5) = 0.
Таким образом, для этой модели МА второго порядка все прогнозируемые выходные сигналы, которые более чем на два временных шага превышают последнюю измеренную точку данных (t = 2), равны нулю. В общем случае, когда предполагаются нулевые начальные условия, прогнозируемые значения за пределами порядка чистой модели МА без авторегрессионных членов равны нулю.
Предполагаемые начальные условия: Можно оценить начальные условия, минимизируя квадратную сумму 1-ступенчатых ошибок прогнозирования всех измеренных данных.
Для модели MA, описанной ранее, оценка начальных условий e(0) и e(-1) требует минимизации следующей функции затрат на наименьшие квадраты:
$V = e^{(} 1) 2^{+} e (2) 2 = (y (1) - [{0,1 e (0) +}^{} 0,2 e (-1 )] ) 2 + (y ( 2) - {[ 0,1 e (}^{1}$ ) + 0,2 e (0)]) 2
Замена a = e(0) и b = e(-1), функция затрат:
$V (a, b {) = (5- [0 . 1a}^{+} {0 . 2b]) 2 + (10 - [0,1 (5 - [0 .}^{}$ 1a + 0 .2b]) + 0 .2a]) 2
Минимизация выхода V e(0) = 50 и e(-1) = 0, что дает:
$\begin{array}{l} e (1) = 5- (0,1 * \\ 50 + 0,2 * 0) = 0e (2) = 10 - \\ \overset{}{(} 0,1 * \\ \overset{}{0}+ 0,2 * \end{array}$ 50) = 0y ^ (3) = 0y ^ (4) = 0
Таким образом, для этой системы, если ошибки прогнозирования минимизируются по имеющимся двум выборкам, все будущие прогнозы равны нулю, что является средним значением процесса. Если бы было более двух доступных наблюдений, вы бы оценили e(-1) и e(0) использование метода наименьших квадратов для минимизации одношаговых ошибок прогнозирования по всем доступным данным.
В этом примере показано, как воспроизвести эти прогнозируемые результаты с помощью 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
```

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

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

Дискретная пространственно-временная модель данных временных рядов имеет вид:

$\begin{array}{l} x (t + 1) = Ax (t) \\ + Ke (t) y (t) \end{array}$ = Cx (t) + e (t)

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

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

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

Для прогнозирования реакции модели «состояние-пространство»:

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

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

Открыть сценарий в реальном времени

В этом примере показано, как вручную воспроизводить результаты прогнозирования, полученные с помощью 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);

Прогнозирование выхода модели state-space на 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.

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

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

Открыть сценарий в реальном времени

Когда существуют экзогенные стимулы, влияющие на систему, система не может считаться стационарной. Однако, если эти стимулы измеримы, то вы можете рассматривать их как входные данные для системы и учитывать их влияние при прогнозировании выхода системы. Поток операций для прогнозирования данных с внешними входами аналогичен потоку операций для прогнозирования данных временных рядов. Сначала необходимо определить модель в соответствии с измеренными данными ввода-вывода. Затем необходимо указать ожидаемые входные значения для периода времени прогнозирования и спрогнозировать выходные данные идентифицированной модели с помощью 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) 2, 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=[w₁ w₂ w₃ w₄ w₅] является вектором взвешивания, и c - смещение выходного сигнала.

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

$y (t)_{=} w1y (_{t-1}) + w2y_{(} {t-2) +}^{}_{w3y} {(t-1)}^{} 2_{+} w4y (t-2) 2 + w5y$ (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} \overset{}{y}^(101)_{} = w1y (_{} 100) +_{w2y} (^{99) +}_{} w3y^{(100})_{2} + w4y (99) \\ \overset{2}{} + w5y_{(} \overset{}{} 100) y (_{99}) y ^(_{} \overset{)}{102}^{= w1y}^_{} (^{101) +}_{} \overset{}{w2y} (100) + w3y \\ \hat{} (101)_{2} \overset{+}{} w4y (_{} \overset{)}{100} 2 +_{w5y} \hat{}^{(101)} y_{} \overset{}{(}^{100) y}^_{(} \overset{}{} 103 \overset{}{} \\ \overset{}{}_{} \overset{}{}_{} \overset{}{}_{} \overset{}{}^{}_{} \overset{}{}^{}_{} \overset{}{} \overset{}{} \end{array}$ )

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

См. также

forecast | predict | sim

Связанные примеры

Подробнее

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

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