В этом примере показано, как выполнить многомерное прогнозирование временных рядов данных, измеренных от хищника и популяций добычи в сценарии давки добычи. Движущие силы изменения населения добычи хищника моделируются с помощью линейных и нелинейных моделей временных рядов. Прогнозирование эффективности этих моделей сравнено.
Данные являются двумерными временными рядами, состоящими из популяций с 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 лет. Это указывает, что динамика изменения населения может быть нелинейной.
Подбирайте нелинейную модель черного ящика к данным об оценке. Вы не требуете предварительных знаний об уравнениях, управляющих данными об оценке. Будет оценена линейная в регрессоре форма модели Nonlinear ARX.
Создайте нелинейную модель ARX с 2 выходными параметрами и никакими входными параметрами.
sysNLARX = idnlarx([1 1;1 1],[],'Ts',0.1,'TimeUnit','years');
sysNLARX
первый порядок нелинейная модель ARX, которая не использует нелинейной функции; это предсказывает ответ модели с помощью взвешенной суммы его регрессоров первого порядка.
getreg(sysNLARX)
Regressors: For output 1: y1(t-1) y2(t-1) For output 2: y1(t-1) y2(t-1)
Чтобы ввести функцию нелинейности, добавьте полиномиальные регрессоры в модель.
Создайте регрессоры до степени 2 и включайте перекрестные условия (продукты стандартных упомянутых выше регрессоров). Добавьте те регрессоры в модель как пользовательские регрессоры.
R = polyreg(sysNLARX,'MaxPower',2,'CrossTerm','on'); sysNLARX.CustomRegressors = R; getreg(sysNLARX)
Regressors: For output 1: y1(t-1) y2(t-1) y1(t-1).^2 y1(t-1).*y2(t-1) y2(t-1).^2 For output 2: y1(t-1) y2(t-1) y1(t-1).^2 y1(t-1).*y2(t-1) y2(t-1).^2
Оцените коэффициенты (коэффициенты регрессора и смещение) модели с помощью данных об оценке, ze
.
sysNLARX = nlarx(ze,sysNLARX)
sysNLARX = Nonlinear ARX model with 2 outputs and 0 inputs Inputs: Outputs: y1, y2 Standard regressors corresponding to the orders: na = [1 1; 1 1] nb = [] nk = [] Custom regressors: For output 1: y1(t-1).^2 y1(t-1).*y2(t-1) y2(t-1).^2 For output 2: y1(t-1).^2 y1(t-1).*y2(t-1) y2(t-1).^2 Nonlinear regressors: For output 1: none For output 2: none Model output is linear in their regressors. 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)'); 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™.
Если у вас есть предварительные знания системной динамики, можно соответствовать данным об оценке с помощью нелинейной модели серого ящика.
Знание о природе динамики может помочь улучшить качество модели и таким образом точность прогнозирования. Для динамики добычи хищника, изменений у хищника (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 шагов, вперед предсказал выход, чтобы подтвердить модель.
yhat3 = predict(sysGrey,ze,10);
Предскажите выход шагов модели 100 вне данных об оценке и вычислите выходные стандартные отклонения.
[yf3,x03,sysf3,ysd3] = forecast(sysGrey,ze,100);
Постройте измеренное, предсказанный, и предсказал выходные параметры.
t = yf3.SamplingInstants; subplot(1,2,1); plot(t0,z.y(:,1),... te,yhat3.y(:,1),... t,yf3.y(:,1),'m',... t,yf3.y(:,1)+ysd3(:,1),'k--', ... t,yf3.y(:,1)-ysd3(:,1),'k--') xlabel('Time (year)'); ylabel('Predator population (thousands)'); subplot(1,2,2); plot(t0,z.y(:,2),... te,yhat3.y(:,2),... t,yf3.y(:,2),'m',... t,yf3.y(:,2)+ysd3(:,2),'k--', ... t,yf3.y(:,2)-ysd3(:,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) legend({'Measured','Linear ARMA','Nonlinear AR','Nonlinear Grey-Box'}) title('Forecasted Responses')
Прогнозирование с нелинейной моделью ARX дало лучшие результаты, чем прогнозирование с линейной моделью. Включение знания динамики в нелинейной модели серого ящика далее улучшило надежность модели и поэтому точности прогнозирования.
Обратите внимание на то, что уравнения, используемые в моделировании серого ящика, тесно связаны с полиномиальными регрессорами, используемыми моделью Nonlinear ARX. Если вы аппроксимируете производные в управляющих уравнениях различиями первого порядка, вы получите уравнения, похожие на:
Где, некоторые параметры, связанные с исходными параметрами и к шагу расчета используется для дифференцирования. Эти уравнения предполагают, что только 2 регрессора необходимы для первого выхода (y1 (t-1) и *y1 (t-1) *y2 (t-1)) и 3 для второго выхода при построении модели Nonlinear ARX. Даже в отсутствие таких предварительных знаний, линейные в регрессоре структуры модели, использующие полиномиальные регрессоры, остаются популярный выбор на практике.
Предскажите значения с помощью нелинейной модели серого ящика более чем 200 лет.
[yf4,~,~,ysd4] = forecast(sysGrey, ze, 2000);
Постройте последнюю часть данных (показав 1 sd неопределенность)
t = yf4.SamplingInstants; N = 700:2000; subplot(1,2,1); plot(t(N), yf4.y(N,1), 'm',... t(N), yf4.y(N,1)+ysd4(N,1), 'k--', ... t(N), yf4.y(N,1)-ysd4(N,1), 'k--') xlabel('Time (year)'); ylabel('Predator population (thousands)'); ax = gca; ax.YLim = [0.8 1]; ax.XLim = [82 212]; subplot(1,2,2); plot(t(N),yf4.y(N,2),'m',... t(N),yf4.y(N,2)+ysd4(N,2),'k--', ... t(N),yf4.y(N,2)-ysd4(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.