В этом примере показано, как выполнить многомерное прогнозирование временных рядов данных, измеренных от хищника и популяций добычи в сценарии давки добычи. Движущие силы изменения населения добычи хищника моделируются с помощью линейных и нелинейных моделей временных рядов. Прогнозирование эффективности этих моделей сравнено.
Данные являются двумерными временными рядами, состоящими из популяций с 1 добычей с 1 хищником (в тысячах), собрался 10 раз в год в течение 20 лет. Для получения дополнительной информации о данных, см. Три Экологических Системы Населения: MATLAB и Моделирование Файла MEX на C Timeseries.
Загрузите данные временных рядов.
load PredPreyCrowdingData z = iddata(y,[],0.1,'TimeUnit','years','Tstart',0);
z iddata объект, содержащий два выходных сигнала, y1 и y2, которые относятся к хищнику и популяциям добычи, соответственно. OutputData свойство z содержит данные о населении как 201 2 матрица, такая что z.OutputData(:,1) популяция хищников и z.OutputData(:,2) популяция жертв.
Отобразите данные на графике.
plot(z) title('Predator-Prey Population Data') ylabel('Population (thousands)')
Данные показывают снижение в популяции хищников из-за давки.
Используйте первую половину в качестве данных об оценке для идентификации моделей временных рядов.
ze = z(1:120);
Используйте остающиеся данные, чтобы искать порядки модели и подтвердить результаты прогнозирования.
zv = z(121:end);
Смоделируйте временные ряды как линейный, авторегрессивный процесс. Линейные модели могут быть созданы в полиномиальной форме или форме пространства состояний с помощью команд, таких как ar (только для скалярных временных рядов), arx, armax, n4sid и ssest. Поскольку линейные модели не получают смещения данных (ненулевое условное среднее значение), сначала детрендируют данные.
[zed, Tze] = detrend(ze, 0); [zvd, Tzv] = detrend(zv, 0);
Идентификация требует спецификации порядков модели. Для полиномиальных моделей можно найти подходящие порядки с помощью arxstruc команда. Начиная с arxstruc работает только над моделями одно выхода, выполните поиск порядка модели отдельно каждого выхода.
na_list = (1:10)'; V1 = arxstruc(zed(:,1,:),zvd(:,1,:),na_list); na1 = selstruc(V1,0); V2 = arxstruc(zed(:,2,:),zvd(:,2,:),na_list); na2 = selstruc(V2,0);
arxstruc команда предлагает авторегрессивные модели порядков 7 и 8, соответственно.
Используйте эти порядки модели оценить модель ARMA мультиотклонения, где перекрестные термины были выбраны произвольно.
na = [na1 na1-1; na2-1 na2]; nc = [na1; na2]; sysARMA = armax(zed,[na nc])
sysARMA =
Discrete-time ARMA model:
Model for output "y1": A(z)y_1(t) = - A_i(z)y_i(t) + C(z)e_1(t)
A(z) = 1 - 0.885 z^-1 - 0.1493 z^-2 + 0.8089 z^-3 - 0.2661 z^-4
- 0.9487 z^-5 + 0.8719 z^-6 - 0.2896 z^-7
A_2(z) = 0.3433 z^-1 - 0.2802 z^-2 - 0.04949 z^-3 + 0.1018 z^-4
- 0.02683 z^-5 - 0.2416 z^-6
C(z) = 1 - 0.4534 z^-1 - 0.4127 z^-2 + 0.7874 z^-3 + 0.298 z^-4
- 0.8684 z^-5 + 0.6106 z^-6 + 0.3616 z^-7
Model for output "y2": A(z)y_2(t) = - A_i(z)y_i(t) + C(z)e_2(t)
A(z) = 1 - 0.5826 z^-1 - 0.4688 z^-2 - 0.5949 z^-3 - 0.0547 z^-4
+ 0.5062 z^-5 + 0.4024 z^-6 - 0.01544 z^-7 - 0.1766 z^-8
A_1(z) = 0.2386 z^-1 + 0.1564 z^-2 - 0.2249 z^-3 - 0.2638 z^-4 - 0.1019 z^-5
- 0.07821 z^-6 + 0.2982 z^-7
C(z) = 1 - 0.1717 z^-1 - 0.09877 z^-2 - 0.5289 z^-3 - 0.24 z^-4
+ 0.06555 z^-5 + 0.2217 z^-6 - 0.05765 z^-7 - 0.1824 z^-8
Sample time: 0.1 years
Parameterization:
Polynomial orders: na=[7 6;7 8] nc=[7;8]
Number of free coefficients: 43
Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using ARMAX on time domain data "zed".
Fit to estimation data: [89.85;90.97]% (prediction focus)
FPE: 3.814e-05, MSE: 0.007533
Вычислите 10 шагов вперед (1 год) предсказанный выход, чтобы подтвердить модель по отрезку времени данных об оценке. Поскольку данные были детрендированы для оценки, необходимо задать те смещения для значимых предсказаний.
predOpt = predictOptions('OutputOffset',Tze.OutputOffset');
yhat1 = predict(sysARMA,ze,10, predOpt);predict команда предсказывает ответ по отрезку времени результатов измерений и является инструментом для проверки качества предполагаемой модели. Ответ во время t вычисляется с помощью измеренных значений время от времени t = 0..., t-10.
Постройте предсказанный ответ и результаты измерений.
plot(ze,yhat1)
title('10-step predicted response compared to measured data')
Отметьте, генерация предсказанного ответа и графического вывода его с результатами измерений, может быть автоматизирован с помощью compare команда.
compareOpt = compareOptions('OutputOffset',Tze.OutputOffset');
compare(ze,sysARMA,10,compareOpt)
График, сгенерированный с помощью compare также показывает нормированное среднеквадратичное значение (NRMSE) мера качества подгонки в форме процента.
После проверки данных, прогноз выход модели sysARMA 100 шагов (10 лет) вне данных об оценке, и вычисляют выходные стандартные отклонения.
forecastOpt = forecastOptions('OutputOffset',Tze.OutputOffset');
[yf1,x01,sysf1,ysd1] = forecast(sysARMA, ze, 100, forecastOpt);yf1 предсказанный ответ, возвращенный как iddata объект. yf1.OutputData содержит предсказанные значения.
sysf1 система, похожая на sysARMA но находится в форме пространства состояний. Симуляция sysf1 использование sim команда, с начальными условиями, x01, воспроизводит предсказанный ответ, yf1.
ysd1 матрица стандартных отклонений. Это измеряется, неопределенность предсказывает вследствие эффекта аддитивных воздействий в данных (как измерено sysARMA.NoiseVariance), неопределенность параметра (как сообщается getcov(sysARMA)) и неопределенность сопоставила с процессом отображения прошлых данных к начальным условиям, требуемым для прогнозирования (см. data2state).
Постройте измеренный, предсказанный, и предсказанный выход для модели sysARMA.
t = yf1.SamplingInstants; te = ze.SamplingInstants; t0 = z.SamplingInstants; subplot(1,2,1); plot(t0,z.y(:,1),... te,yhat1.y(:,1),... t,yf1.y(:,1),'m',... t,yf1.y(:,1)+ysd1(:,1),'k--', ... t,yf1.y(:,1)-ysd1(:,1), 'k--') xlabel('Time (year)'); ylabel('Predator population, in thousands'); subplot(1,2,2); plot(t0,z.y(:,2),... te,yhat1.y(:,2),... t,yf1.y(:,2),'m',... t,yf1.y(:,2)+ysd1(:,2),'k--', ... t,yf1.y(:,2)-ysd1(:,2),'k--') % Make the figure larger. fig = gcf; p = fig.Position; fig.Position = [p(1),p(2)-p(4)*0.2,p(3)*1.4,p(4)*1.2]; xlabel('Time (year)'); ylabel('Prey population, in thousands'); legend({'Measured','Predicted','Forecasted','Forecast Uncertainty (1 sd)'},... 'Location','best')
Графики показывают, что прогнозирование использования линейной модели ARMA (с добавленной обработкой смещений) работало несколько, и результаты показали высокую неопределенность по сравнению с фактическими популяциями по отрезку времени этих 12-20 лет. Это указывает, что динамика изменения населения может быть нелинейной.
В этом разделе мы расширяем идентификационный подход черного ящика предыдущего раздела путем добавления нелинейных элементов в модель. Это выполняется с помощью Нелинейной среды моделирования ARX. Нелинейные модели ARX концептуально похожи на линейные единицы (такие как sysARMA выше) в этом они вычисляют предсказанный выход с помощью двухступенчатого процесса:
Сопоставьте измеренные сигналы с конечномерным пробелом регрессора.
Сопоставьте регрессоры с предсказанным выходом с помощью статической функции, которая может быть линейной или нелинейной.
Нелинейность может быть добавлена в двух местах - любой в наборе регрессора (шаг 1), и/или в статической функции отображения (шаг 2). В этом разделе мы исследуем обе из этих опций. Во-первых, мы создаем линейную в регрессоре форму модели Nonlinear ARX, где во всей нелинейности ограничиваются определением регрессора только; статическая функция отображения, которая преобразовывает регрессоры в предсказанный выход, линейна (аффинно, чтобы быть точной). Затем мы создаем модель, которая использует Гауссову функцию Процесса в качестве ее статической функции отображения, но сохраняет ее регрессоры линейными.
Создайте нелинейную модель ARX с 2 выходными параметрами и никакими входными параметрами.
L = linearRegressor(ze.OutputName,{1,1});
sysNLARX = idnlarx(ze.OutputName, ze.InputName, L, []);
sysNLARX.Ts = 0.1;
sysNLARX.TimeUnit = 'years';sysNLARX первый порядок нелинейная модель ARX, которая не использует нелинейной функции; это предсказывает ответ модели с помощью взвешенной суммы его регрессоров первого порядка.
getreg(sysNLARX)
ans = 2x1 cell
{'y1(t-1)'}
{'y2(t-1)'}
Чтобы ввести функцию нелинейности, добавьте полиномиальные регрессоры в модель.
Создайте регрессоры порядка 2 и включайте перекрестные термины (продукты стандартных упомянутых выше регрессоров). Добавьте те регрессоры в модель.
P = polynomialRegressor(ze.OutputName,{1,1},2,false,true);
sysNLARX.Regressors(end+1) = P;
getreg(sysNLARX)ans = 5x1 cell
{'y1(t-1)' }
{'y2(t-1)' }
{'y1(t-1)^2' }
{'y2(t-1)^2' }
{'y1(t-1)*y2(t-1)'}
Оцените коэффициенты (коэффициенты регрессора и смещение) модели с помощью данных об оценке, ze.
sysNLARX = nlarx(ze,sysNLARX)
sysNLARX = Nonlinear multivariate time series model with 2 outputs Outputs: y1, y2 Regressors: 1. Linear regressors in variables y1, y2 2. Order 2 regressors in variables y1, y2 List of all regressors Output functions: Output 1: None Output 2: None Sample time: 0.1 years Status: Estimated using NLARX on time domain data "ze". Fit to estimation data: [88.34;88.91]% (prediction focus) FPE: 3.265e-05, MSE: 0.01048
Вычислите 10 шагов, вперед предсказал выход, чтобы подтвердить модель.
yhat2 = predict(sysNLARX,ze,10);
Предскажите выход шагов модели 100 вне данных об оценке.
yf2 = forecast(sysNLARX,ze,100);
Стандартные отклонения предсказанного ответа не вычисляются для нелинейных моделей ARX. Эти данные недоступны, потому что информация о ковариации параметра не вычисляется во время оценки этих моделей.
Постройте измеренное, предсказанный, и предсказал выходные параметры.
t = yf2.SamplingInstants; subplot(1,2,1); plot(t0,z.y(:,1),... te,yhat2.y(:,1),... t,yf2.y(:,1),'m') xlabel('Time (year)'); ylabel('Predator population (thousands)'); title('Prediction Results for ARX Model with Polynomial Regressors') subplot(1,2,2); plot(t0,z.y(:,2),... te,yhat2.y(:,2),... t,yf2.y(:,2),'m') legend('Measured','Predicted','Forecasted') xlabel('Time (year)'); ylabel('Prey population (thousands)');
Графики показывают, что прогнозирование использования нелинейной модели ARX дало лучшие результаты прогнозирования, чем использование линейной модели. Нелинейное моделирование черного ящика не потребовало предварительных знаний об уравнениях, управляющих данными.
Обратите внимание на то, что, чтобы сократить количество регрессоров, можно выбрать оптимальное подмножество (преобразованных) переменных с помощью анализа главных компонентов (см. pca) или выбор признаков (см. sequentialfs) в Statistics and Machine Learning Toolbox™.
Модель модель sysNLARX использует линейную функцию отображения; нелинейность содержится в регрессорах только. Более гибкий подход включает выбор нетривиальной функции отображения, такой как Гауссова функция Процесса (GP). Создайте нелинейную модель ARX, которая использует только линейные регрессоры, но добавляет Гауссову функцию процесса для отображения регрессоров к выходу.
% Copy sysNLARX so we don't have to create the model structure again sysGP = sysNLARX; % Remove polynomial regressor set sysGP.Regressors(2) = []; % Replace the linear mapping function by a Gaussian Process (GP) function. % The GP uses a squared exponential kernel by default. GP = idGaussianProcess; sysGP.OutputFcn = GP; disp(sysGP)
Nonlinear multivariate time series model with 2 outputs Outputs: y1, y2 Regressors: Linear regressors in variables y1, y2 Output functions: Output 1: Gaussian process function using a SquaredExponential kernel. Output 2: Gaussian process function using a SquaredExponential kernel. Sample time: 0.1 years Status: Estimated using NLARX with prediction focus
У вас теперь есть модель шаблона sysGP, который основан на линейных регрессорах и функции отображения GP. Его параметры являются двумя коэффициентами его экспоненциального ядра в квадрате. Используйте nlarx оценить эти параметры.
sysGP = nlarx(ze, sysGP)
sysGP = Nonlinear multivariate time series model with 2 outputs Outputs: y1, y2 Regressors: Linear regressors in variables y1, y2 List of all regressors Output functions: Output 1: Gaussian process function using a SquaredExponential kernel. Output 2: Gaussian process function using a SquaredExponential kernel. Sample time: 0.1 years Status: Estimated using NLARX on time domain data "ze". Fit to estimation data: [88.07;88.84]% (prediction focus) FPE: 3.369e-05, MSE: 0.01083
getpvec(sysGP)
ans = 10×1
-0.6671
0.0196
0.9558
0.9445
0.0261
-0.0221
-0.5752
0.8479
1.1883
0.0447
Как прежде, вычислите предсказанные и предсказанные ответы и постройте результаты.
Вычислите 10 шагов, вперед предсказал выход, чтобы подтвердить модель.
yhat3 = predict(sysGP,ze,10);
Предскажите выход шагов модели 100 вне данных об оценке.
yf3 = forecast(sysGP,ze,100);
Постройте измеренное, предсказанный, и предсказал выходные параметры.
t = yf3.SamplingInstants; subplot(1,2,1); plot(t0,z.y(:,1),... te,yhat3.y(:,1),... t,yf3.y(:,1),'m') xlabel('Time (year)'); ylabel('Predator population (thousands)'); title('Prediction Results for Gaussian Process Model') subplot(1,2,2); plot(t0,z.y(:,2),... te,yhat3.y(:,2),... t,yf3.y(:,2),'m') legend('Measured','Predicted','Forecasted') xlabel('Time (year)'); ylabel('Prey population (thousands)');
Результаты показывают улучшенную точность по линейной модели. Преимущество использования основанных на GP моделей состоит в том, что они используют очень немного параметров, в противоположность другим формам, таким как Сети Вейвлета (idWaveletNetwork) или сигмоидальные сети (idSigmoidNetwork). Это делает их легкими обучаться с низкой неопределенностью в их оценках параметра.
Если у вас есть предварительные знания системной динамики, можно соответствовать данным об оценке с помощью нелинейной модели серого ящика. Знание о природе динамики может, в целом, помочь улучшить качество модели и таким образом точность прогнозирования. Для динамики добычи хищника, изменений у хищника (y1) и добыча (y2) население может быть представлено как:
Для получения дополнительной информации об уравнениях, см. Три Экологических Системы Населения: MATLAB и Моделирование Файла MEX на C Timeseries.
Создайте нелинейный серый ящик, основанный на модели на этих уравнениях.
Задайте файл, описывающий структуру модели для системы добычи хищника. Файл задает производные состояния и выходные параметры модели в зависимости от времени, состояний, входных параметров и параметров модели. Эти два выходных параметров (хищник и популяции добычи) выбраны в качестве состояний, чтобы вывести нелинейное описание пространства состояний динамики.
FileName = 'predprey2_m';Задайте порядки модели (количество выходных параметров, входных параметров и состояний)
Order = [2 0 2];
Задайте начальные значения для параметров , , , , и , и укажите, что все параметры должны быть оценены. Обратите внимание на то, что требование, чтобы задать исходные предположения для параметров не существовало при оценке моделей sysARMA черного квадрата и sysNLARX.
Parameters = struct('Name',{'Survival factor, predators' 'Death factor, predators' ... 'Survival factor, preys' 'Death factor, preys' ... 'Crowding factor, preys'}, ... 'Unit',{'1/year' '1/year' '1/year' '1/year' '1/year'}, ... 'Value',{-1.1 0.9 1.1 0.9 0.2}, ... 'Minimum',{-Inf -Inf -Inf -Inf -Inf}, ... 'Maximum',{Inf Inf Inf Inf Inf}, ... 'Fixed',{false false false false false});
Точно так же задайте начальные состояния модели и укажите, что должны быть оценены оба начальных состояния.
InitialStates = struct('Name',{'Predator population' 'Prey population'}, ... 'Unit',{'Size (thousands)' 'Size (thousands)'}, ... 'Value',{1.8 1.8}, ... 'Minimum', {0 0}, ... 'Maximum',{Inf Inf}, ... 'Fixed',{false false});
Задайте модель как систему непрерывного времени.
Ts = 0;
Создайте нелинейную модель серого ящика с заданной структурой, параметрами и состояниями.
sysGrey = idnlgrey(FileName,Order,Parameters,InitialStates,Ts,'TimeUnit','years');
Оцените параметры модели.
sysGrey = nlgreyest(ze,sysGrey); present(sysGrey)
sysGrey =
Continuous-time nonlinear grey-box model defined by 'predprey2_m' (MATLAB file):
dx/dt = F(t, x(t), p1, ..., p5)
y(t) = H(t, x(t), p1, ..., p5) + e(t)
with 0 input(s), 2 state(s), 2 output(s), and 5 free parameter(s) (out of 5).
States: Initial value
x(1) Predator population(t) [Size (thou..] xinit@exp1 2.01325 (estimated) in [0, Inf]
x(2) Prey population(t) [Size (thou..] xinit@exp1 1.99687 (estimated) in [0, Inf]
Outputs:
y(1) y1(t)
y(2) y2(t)
Parameters: Value Standard Deviation
p1 Survival factor, predators [1/year] -0.995895 0.0125269 (estimated) in [-Inf, Inf]
p2 Death factor, predators [1/year] 1.00441 0.0129368 (estimated) in [-Inf, Inf]
p3 Survival factor, preys [1/year] 1.01234 0.0135413 (estimated) in [-Inf, Inf]
p4 Death factor, preys [1/year] 1.01909 0.0121026 (estimated) in [-Inf, Inf]
p5 Crowding factor, preys [1/year] 0.103244 0.0039285 (estimated) in [-Inf, Inf]
Status:
Termination condition: Change in cost was less than the specified tolerance..
Number of iterations: 6, Number of function evaluations: 7
Estimated using Solver: ode45; Search: lsqnonlin on time domain data "ze".
Fit to estimation data: [91.21;92.07]%
FPE: 8.613e-06, MSE: 0.005713
More information in model's "Report" property.
Вычислите 10 шагов, вперед предсказал выход, чтобы подтвердить модель.
yhat4 = predict(sysGrey,ze,10);
Предскажите выход шагов модели 100 вне данных об оценке и вычислите выходные стандартные отклонения.
[yf4,x04,sysf4,ysd4] = forecast(sysGrey,ze,100);
Постройте измеренное, предсказанный, и предсказал выходные параметры.
t = yf4.SamplingInstants; subplot(1,2,1); plot(t0,z.y(:,1),... te,yhat4.y(:,1),... t,yf4.y(:,1),'m',... t,yf4.y(:,1)+ysd4(:,1),'k--', ... t,yf4.y(:,1)-ysd4(:,1),'k--') xlabel('Time (year)'); ylabel('Predator population (thousands)'); title('Prediction Results for Grey-Box Model') subplot(1,2,2); plot(t0,z.y(:,2),... te,yhat4.y(:,2),... t,yf4.y(:,2),'m',... t,yf4.y(:,2)+ysd4(:,2),'k--', ... t,yf4.y(:,2)-ysd4(:,2),'k--') legend('Measured','Predicted','Forecasted','Forecast uncertainty (1 sd)') xlabel('Time (years)'); ylabel('Prey population (thousands)');
Графики показывают, что прогнозирование использования нелинейной модели серого ящика дало хорошие результаты прогнозирования и низко прогнозирование выходной неопределенности.
Сравните предсказанный ответ, полученный из идентифицированных моделей, sysARMA, sysNLARX, и sysGrey. Первые два являются моделями дискретного времени и sysGrey модель непрерывного времени.
clf
plot(z,yf1,yf2,yf3,yf4)
legend({'Measured','ARMA','Polynomial AR','GP','Grey-Box'})
title('Forecasted Responses')
xlim([7 23])
Прогнозирование с нелинейной моделью ARX дало лучшие результаты, чем прогнозирование с линейной моделью. Включение знания динамики в нелинейной модели серого ящика далее улучшило надежность модели и поэтому точности прогнозирования.
Обратите внимание на то, что уравнения, используемые в моделировании серого ящика, тесно связаны с полиномиальными регрессорами, используемыми моделью Nonlinear ARX. Если вы аппроксимируете производные в управляющих уравнениях различиями первого порядка, вы получите уравнения, похожие на:
Где, некоторые параметры, связанные с исходными параметрами и к шагу расчета используется для дифференцирования. Эти уравнения предполагают, что только 2 регрессора необходимы для первого выхода (y1 (t-1) и *y1 (t-1) *y2 (t-1)) и 3 для второго выхода при построении модели Nonlinear ARX. Даже в отсутствие таких предварительных знаний, линейные в регрессоре структуры модели, использующие полиномиальные регрессоры, остаются популярный выбор на практике.
Предскажите значения с помощью нелинейной модели серого ящика более чем 200 лет.
[yf5,~,~,ysd5] = forecast(sysGrey, ze, 2000);
Постройте последнюю часть данных (показав 1 sd неопределенность)
t = yf5.SamplingInstants; N = 700:2000; subplot(1,2,1); plot(t(N), yf5.y(N,1), 'm',... t(N), yf5.y(N,1)+ysd5(N,1), 'k--', ... t(N), yf5.y(N,1)-ysd5(N,1), 'k--') xlabel('Time (year)'); ylabel('Predator population (thousands)'); title('Forecasted Population after 200 Years') ax = gca; ax.YLim = [0.8 1]; ax.XLim = [82 212]; subplot(1,2,2); plot(t(N),yf5.y(N,2),'m',... t(N),yf5.y(N,2)+ysd5(N,2),'k--', ... t(N),yf5.y(N,2)-ysd5(N,2),'k--') legend('Forecasted population','Forecast uncertainty (1 sd)') xlabel('Time (years)'); ylabel('Prey population (thousands)'); ax = gca; ax.YLim = [0.9 1.1]; ax.XLim = [82 212];
График показывает, что завышенное население, как предсказывают, достигает установившегося из приблизительно 890, и популяция жертв, как предсказывают, достигает 990.