Три экологических системы генеральной совокупности: MATLAB и моделирование файла MEX на C timeseries
Этот пример показывает, как создать нелинейные серые модели временных рядов поля. Модели временных рядов являются моделями, не используя измеренных входных параметров. Три идеализированных экологических системы изучены где две разновидности также:
конкурируйте за ту же еду, или:
находятся в ситуации добычи хищника
Пример показывает моделирование и на основе MATLAB и на основе файлов MEX на C.
Экологические системы генеральной совокупности
Во всех трех системах генеральной совокупности, исследованных, мы интересны в том, как генеральная совокупность двух разновидностей меняется в зависимости от времени. Чтобы смоделировать это, позвольте x1 (t), и x2 (t) обозначают количество людей соответствующих разновидностей во время t. Позвольте l1, и l2 обозначают уровень рождаемости, сопоставленный с x1 (t) и x2 (t), соответственно, оба принятые, чтобы быть константами в зависимости от времени. Показатели смертности разновидностей зависят и от доступности еды и, если хищники присутствуют на риске того, чтобы быть съеденным. Довольно часто (и в целом), показатель смертности для разновидностей i (i = 1 или 2) может быть записан ui (x1 (t), x2 (t)), где ui(.) является некоторой соответствующей функцией. На практике это означает, что ui (x1 (t), x2 (t)) *xi (t) животные разновидностей i умирает каждая единица измерения времени. Результирующий эффект этих операторов может быть получен в итоге в типе пространства состояний образцовой структуры: (timeseries):
d
-- x1(t) = l1*x1(t) - u1(x1(t), x2(t))*x1(t)
dt
d
-- x2(t) = l2*x2(t) - u2(x1(t), x2(t))*x2(t)
dt
Здесь естественно выбрать два состояния в качестве выходных параметров, т.е. мы позволяем y1 (t) = x1 (t) и y2 (t) = x2 (t).
1. Две Разновидности, которые Конкурируют за Ту же Еду
В случае, если две разновидности конкурируют за ту же еду, затем это - полная генеральная совокупность разновидности, которая управляет доступностью еды, и в свою очередь их показателей смертности. Простой все же общий подход должен принять, что показатель смертности может быть записан как:
ui(x1(t), x2(t)) = gi + di*(x1(t) + x2(t))
для обеих разновидностей (i = 1 или 2), где gi и di являются неизвестными параметрами. В целом это дает структуру пространства состояний:
d
-- x1(t) = (l1-g1)*x1(t) - d1*(x1(t)+x2(t))*x1(t)
dt
d
-- x2(t) = (l2-g2)*x2(t) - d2*(x1(t)+x2(t))*x2(t)
dt
Мгновенная проблема с этой структурой состоит в том, что l1, g1, l2, и g2 не может быть идентифицирован отдельно. Мы можем только надеяться идентифицировать p1 = l1-g1 и p3 = l2-g2. Также позволяя p2 = d1 и p4 = d2, каждый получает повторно параметризованную образцовую структуру:
d
-- x1(t) = p1*x1(t) - p2*(x1(t)+x2(t))*x1(t)
dt
d
-- x2(t) = p3*x2(t) - p4*(x1(t)+x2(t))*x2(t)
dt
y1(t) = x1(t)
y2(t) = x2(t)
В этом первом примере генеральной совокупности мы обращаемся к моделированию файла MATLAB. Уравнения выше затем вводятся в файл MATLAB, preys_m.m, со следующим содержимым.
function [dx, y] = preys_m(t, x, u, p1, p2, p3, p4, varargin)
%PREYS_M Two species that compete for the same food.
% Output equations.
y = [x(1); ... % Prey species 1.
x(2) ... % Prey species 2.
];
% State equations.
dx = [p1*x(1)-p2*(x(1)+x(2))*x(1); ... % Prey species 1.
p3*x(2)-p4*(x(1)+x(2))*x(2) ... % Prey species 2.
];
Файл MATLAB, наряду с начальным вектором параметра, соответствующим начальным состоянием и некоторой административной информацией затем питается как входные параметры конструктора Object IDNLGREY. Заметьте, что начальные значения и параметров и начальных состояний заданы как массивы структур с Npo (количество объектов параметра = количество параметров, если все параметры являются скалярами), и Nx (количество состояний) элементы, соответственно. Через эти массивы структур возможно полностью присвоить значения свойств не по умолчанию ('Name', 'Модуля', 'Значение, 'Минимум', 'Максимум', и 'Фиксированный') к каждому параметру и начальному состоянию. Здесь мы присвоили 'Минимальное' значение каждого начального состояния, чтобы обнулить (население положительно!), и также заданный, что оба начальных состояния должны быть оценены по умолчанию.
FileName = 'preys_m'; % File describing the model structure. Order = [2 0 2]; % Model orders [ny nu nx]. Parameters = struct('Name', {'Survival factor, species 1' 'Death factor, species 1' ... 'Survival factor, species 2' 'Death factor, species 2'}, ... 'Unit', {'1/year' '1/year' '1/year' '1/year'}, ... 'Value', {1.8 0.8 1.2 0.8}, ... 'Minimum', {-Inf -Inf -Inf -Inf}, ... 'Maximum', {Inf Inf Inf Inf}, ... 'Fixed', {false false false false}); % Estimate all 4 parameters. InitialStates = struct('Name', {'Population, species 1' 'Population, species 2'}, ... 'Unit', {'Size (in thousands)' 'Size (in thousands)'}, ... 'Value', {0.2 1.8}, ... 'Minimum', {0 0}, ... 'Maximum', {Inf Inf}, ... 'Fixed', {false false}); % Estimate both initial states. Ts = 0; % Time-continuous system. nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts, ... 'Name', 'Two species competing for the same food', ... 'OutputName', {'Population, species 1' 'Population, species 2'}, ... 'OutputUnit', {'Size (in thousands)' 'Size (in thousands)'}, ... 'TimeUnit', 'year');
Команда PRESENT может использоваться, чтобы просмотреть информацию о первоначальной модели:
present(nlgr);
nlgr =
Continuous-time nonlinear grey-box model defined by 'preys_m' (MATLAB file):
dx/dt = F(t, x(t), p1, ..., p4)
y(t) = H(t, x(t), p1, ..., p4) + e(t)
with 0 input(s), 2 state(s), 2 output(s), and 4 free parameter(s) (out of 4).
States: Initial value
x(1) Population, species 1(t) [Size (in t..] xinit@exp1 0.2 (estimated) in [0, Inf]
x(2) Population, species 2(t) [Size (in t..] xinit@exp1 1.8 (estimated) in [0, Inf]
Outputs:
y(1) Population, species 1(t) [Size (in thousands)]
y(2) Population, species 2(t) [Size (in thousands)]
Parameters: Value
p1 Survival factor, species 1 [1/year] 1.8 (estimated) in [-Inf, Inf]
p2 Death factor, species 1 [1/year] 0.8 (estimated) in [-Inf, Inf]
p3 Survival factor, species 2 [1/year] 1.2 (estimated) in [-Inf, Inf]
p4 Death factor, species 2 [1/year] 0.8 (estimated) in [-Inf, Inf]
Name: Two species competing for the same food
Status:
Created by direct construction or transformation. Not estimated.
More information in model's "Report" property.
2. Данные ввода - вывода
Мы затем загружаем (моделируемый, хотя поврежденный шум) данные, и создайте объект IDDATA, описывающий одну конкретную ситуацию, где две разновидности конкурируют за ту же еду. Этот набор данных содержит 201 выборку данных, покрывающую 20 лет эволюции.
load(fullfile(matlabroot, 'toolbox', 'ident', 'iddemos', 'data', 'preydata')); z = iddata(y, [], 0.1, 'Name', 'Two species competing for the same food'); set(z, 'OutputName', {'Population, species 1', 'Population, species 2'}, ... 'Tstart', 0, 'TimeUnit', 'Year');
3. Производительность первоначальных двух моделей разновидностей
Симуляция с первоначальной моделью ясно показывает, что это не может справиться с истинной демографической динамикой. См., что график фигурирует. Поскольку тип timeseries модели IDNLGREY замечает, что образцовый вывод определяется начальным состоянием.
compare(z, nlgr, 1);
Рисунок 1: Сравнение между истинными выходными параметрами и моделируемыми выходными параметрами первоначальных двух моделей разновидностей.
4. Оценка параметра
В порядке преодолеть довольно низкую производительность первоначальной модели мы продолжаем оценивать 4 неизвестных параметра и эти 2 начальных состояния с помощью NLGREYEST. Задайте опции оценки с помощью NLGREYESTOPTIONS; в этом случае 'Отображение' установлен в 'on', что означает, что информация о прогрессе оценки отображена в окне прогресса. Можно использовать NLGREYESTOPTIONS, чтобы задать свойства основного алгоритма, такие как 'GradientOptions', 'SearchMethod', 'MaxIterations', 'Допуск', 'Отображение'.
opt = nlgreyestOptions;
opt.Display = 'on';
opt.SearchOptions.MaxIterations = 50;
nlgr = nlgreyest(z, nlgr, opt);5. Производительность приблизительно двух моделей разновидностей
Ориентировочные стоимости параметров и начальных состояний хорошо соответствуют используемым, чтобы сгенерировать истинные выходные данные:
disp(' True Estimated parameter vector');True Estimated parameter vector
ptrue = [2; 1; 1; 1];
fprintf(' %6.3f %6.3f\n', [ptrue'; getpvec(nlgr)']); 2.000 2.004
1.000 1.002
1.000 1.018
1.000 1.010
disp(' ');disp(' True Estimated initial states');True Estimated initial states
x0true = [0.1; 2]; fprintf(' %6.3f %6.3f\n', [x0true'; cell2mat(getinit(nlgr, 'Value'))']);
0.100 0.101
2.000 1.989
Чтобы далее оценить качество модели (и проиллюстрировать улучшение по сравнению с первоначальной моделью), мы также моделируем предполагаемую модель. Моделируемые выходные параметры сравниваются с истинными выходными параметрами в окне графика. Как видно, предполагаемая модель довольно хороша.
compare(z, nlgr, 1);
Рисунок 2: Сравнение между истинными выходными параметрами и моделируемыми выходными параметрами приблизительно двух моделей разновидностей.
PRESENT предоставляет дополнительную информацию о предполагаемой модели, например, о неуверенности параметра и других связанных с оценкой количествах, как функция потерь и FPE Акэйка (Итоговая Ошибка Прогноза) мера.
present(nlgr);
nlgr =
Continuous-time nonlinear grey-box model defined by 'preys_m' (MATLAB file):
dx/dt = F(t, x(t), p1, ..., p4)
y(t) = H(t, x(t), p1, ..., p4) + e(t)
with 0 input(s), 2 state(s), 2 output(s), and 4 free parameter(s) (out of 4).
States: Initial value
x(1) Population, species 1(t) [Size (in t..] xinit@exp1 0.100729 (estimated) in [0, Inf]
x(2) Population, species 2(t) [Size (in t..] xinit@exp1 1.98855 (estimated) in [0, Inf]
Outputs:
y(1) Population, species 1(t) [Size (in thousands)]
y(2) Population, species 2(t) [Size (in thousands)]
Parameters: ValueStandard Deviation
p1 Survival factor, species 1 [1/year] 2.00429 0.00971109 (estimated) in [-Inf, Inf]
p2 Death factor, species 1 [1/year] 1.00235 0.00501783 (estimated) in [-Inf, Inf]
p3 Survival factor, species 2 [1/year] 1.01779 0.0229598 (estimated) in [-Inf, Inf]
p4 Death factor, species 2 [1/year] 1.0102 0.0163506 (estimated) in [-Inf, Inf]
Name: Two species competing for the same food
Status:
Termination condition: Near (local) minimum, (norm(g) < tol)..
Number of iterations: 5, Number of function evaluations: 6
Estimated using Solver: ode45; Search: lsqnonlin on time domain data "Two species competing for the same food".
Fit to estimation data: [98.42;97.92]%
FPE: 7.747e-09, MSE: 0.0001743
More information in model's "Report" property.
B.1. Классическая система добычи хищника
Примите теперь, когда первая разновидность живет на второй. Доступность еды для разновидностей 1 затем пропорциональна x2 (t) (количество людей разновидностей 2), что означает, что показатель смертности разновидностей 1 уменьшается, когда x2 (t) увеличивается. Этот факт получен простым выражением:
u1(x1(t), x2(t)) = g1 - a1*x2(t)
где g1 и a1 являются неизвестными параметрами. Точно так же показатель смертности разновидностей 2 увеличится когда количество людей первых увеличений разновидностей, например, согласно:
u2(x1(t), x2(t)) = g2 + a2*x1(t)
где g2 и a2 являются еще двумя неизвестными параметрами. Используя линейный уровень рождаемости, принятый выше, каждый получает структуру пространства состояний:
d
-- x1(t) = (l1-g1)*x1(t) + a1*x2(t)*x1(t)
dt
d
-- x2(t) = (l2-g2)*x2(t) - a2*x1(t)*x2(t)
dt
Как в предыдущем примере генеральной совокупности, также здесь невозможно однозначно определить шесть отдельных параметров. С тем же видом репараметризации как в вышеупомянутом случае, т.е. p1 = l1-g1, p2 = a1, p3 = l2-g2, и p4 = a2, получена следующая образцовая структура:
d
-- x1(t) = p1*x1(t) + p2*x2(t)*x1(t)
dt
d
-- x2(t) = p3*x2(t) - p4*x1(t)*x2(t)
dt
y1(t) = x1(t)
y2(t) = x2(t)
который лучше подходит с точки зрения оценки.
На этот раз мы вводим эту информацию в файл MEX на C, названный predprey1_c.c. Образцовый файл структурирован как стандартный IDNLGREY файл MEX на C (см. пример, названный, "Создавая Файлы Модели IDNLGREY" или idnlgreydemo2.m), с состоянием и выводят функции обновления, compute_dx и compute_y, можно следующим образом.
void compute_dx(double *dx, double t, double *x, double **p,
const mxArray *auxvar)
{
/* Retrieve model parameters. */
double *p1, *p2, *p3, *p4;
p1 = p[0]; /* Survival factor, predators. */
p2 = p[1]; /* Death factor, predators. */
p3 = p[2]; /* Survival factor, preys. */
p4 = p[3]; /* Death factor, preys. */
/* x[0]: Predator species. */
/* x[1]: Prey species. */
dx[0] = p1[0]*x[0]+p2[0]*x[1]*x[0];
dx[1] = p3[0]*x[1]-p4[0]*x[0]*x[1];
}
/* Output equations. */
void compute_y(double *y, double t, double *x, double **p,
const mxArray *auxvar)
{
/* y[0]: Predator species. */
/* y[1]: Prey species. */
y[0] = x[0];
y[1] = x[1];
}
Поскольку модель имеет тип timeseries, ни compute_dx, ни compute_y не включают u в список входных параметров. На самом деле, основная функция интерфейса predprey1_c.c, даже не объявляет u даже при том, что пустой u ([]) всегда передается predprey1_c методами IDNLGREY.
Скомпилированный файл MEX на C, наряду с начальным вектором параметра, соответствующим начальным состоянием и некоторой административной информацией затем питается как входные параметры конструктора Object IDNLGREY:
FileName = 'predprey1_c'; % File describing the model structure. Order = [2 0 2]; % Model orders [ny nu nx]. Parameters = struct('Name', {'Survival factor, predators' 'Death factor, predators' ... 'Survival factor, preys' 'Death factor, preys'}, ... 'Unit', {'1/year' '1/year' '1/year' '1/year'}, ... 'Value', {-1.1 0.9 1.1 0.9}, ... 'Minimum', {-Inf -Inf -Inf -Inf}, ... 'Maximum', {Inf Inf Inf Inf}, ... 'Fixed', {false false false false}); % Estimate all 4 parameters. InitialStates = struct('Name', {'Population, predators' 'Population, preys'}, ... 'Unit', {'Size (in thousands)' 'Size (in thousands)'}, ... 'Value', {1.8 1.8}, ... 'Minimum', {0 0}, ... 'Maximum', {Inf Inf}, ... 'Fixed', {false false}); % Estimate both initial states. Ts = 0; % Time-continuous system. nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts, ... 'Name', 'Classical 1 predator - 1 prey system', ... 'OutputName', {'Population, predators', 'Population, preys'}, ... 'OutputUnit', {'Size (in thousands)' 'Size (in thousands)'}, ... 'TimeUnit', 'year');
Модель добычи хищника затем дословно просматривается через команду PRESENT.
present(nlgr);
nlgr =
Continuous-time nonlinear grey-box model defined by 'predprey1_c' (MEX-file):
dx/dt = F(t, x(t), p1, ..., p4)
y(t) = H(t, x(t), p1, ..., p4) + e(t)
with 0 input(s), 2 state(s), 2 output(s), and 4 free parameter(s) (out of 4).
States: Initial value
x(1) Population, predators(t) [Size (in t..] xinit@exp1 1.8 (estimated) in [0, Inf]
x(2) Population, preys(t) [Size (in t..] xinit@exp1 1.8 (estimated) in [0, Inf]
Outputs:
y(1) Population, predators(t) [Size (in thousands)]
y(2) Population, preys(t) [Size (in thousands)]
Parameters: Value
p1 Survival factor, predators [1/year] -1.1 (estimated) in [-Inf, Inf]
p2 Death factor, predators [1/year] 0.9 (estimated) in [-Inf, Inf]
p3 Survival factor, preys [1/year] 1.1 (estimated) in [-Inf, Inf]
p4 Death factor, preys [1/year] 0.9 (estimated) in [-Inf, Inf]
Name: Classical 1 predator - 1 prey system
Status:
Created by direct construction or transformation. Not estimated.
More information in model's "Report" property.
B.2. Данные ввода - вывода
Наш следующий шаг должен загрузить (моделируемый, хотя поврежденный шум) данные, и создайте объект IDDATA, описывающий эту конкретную ситуацию добычи хищника. Этот набор данных также содержит 201 выборку данных, покрывающую 20 лет эволюции.
load(fullfile(matlabroot, 'toolbox', 'ident', 'iddemos', 'data', 'predprey1data')); z = iddata(y, [], 0.1, 'Name', 'Classical 1 predator - 1 prey system'); set(z, 'OutputName', {'Population, predators', 'Population, preys'}, ... 'Tstart', 0, 'TimeUnit', 'Year');
B.3. Производительность первоначальной классической модели добычи хищника
Симуляция с первоначальной моделью указывает, что это не может точно справиться с истинной демографической динамикой. Смотрите окно графика.
compare(z, nlgr, 1);
Рисунок 3: Сравнение между истинными выходными параметрами и моделируемыми выходными параметрами первоначальной классической модели добычи хищника.
B.4. Оценка параметра
Чтобы улучшать производительность первоначальной модели, мы продолжаем оценивать 4 неизвестных параметра и эти 2 начальных состояния с помощью NLGREYEST.
nlgr = nlgreyest(z, nlgr, nlgreyestOptions('Display', 'on'));
B.5. Производительность предполагаемой классической модели добычи хищника
Ориентировочные стоимости параметров и начальных состояний очень близко к тем, которые использовались, чтобы сгенерировать истинные выходные данные:
disp(' True Estimated parameter vector');True Estimated parameter vector
ptrue = [-1; 1; 1; 1];
fprintf(' %6.3f %6.3f\n', [ptrue'; getpvec(nlgr)']); -1.000 -1.000
1.000 1.000
1.000 1.000
1.000 0.999
disp(' ');disp(' True Estimated initial states');True Estimated initial states
x0true = [2; 2]; fprintf(' %6.3f %6.3f\n', [x0true'; cell2mat(getinit(nlgr, 'Value'))']);
2.000 2.002
2.000 2.000
Чтобы далее оценить качество модели (и проиллюстрировать улучшение по сравнению с первоначальной моделью), мы также моделируем предполагаемую модель. Моделируемые выходные параметры сравниваются с истинными выходными параметрами в окне графика. Как видно снова, предполагаемая модель довольно хороша.
compare(z, nlgr, 1);
Рисунок 4: Сравнение между истинными выходными параметрами и моделируемыми выходными параметрами предполагаемой классической модели добычи хищника.
Как ожидалось ошибки прогноза, возвращенные PE, являются небольшими и случайной природы.
pe(z, nlgr);
Рисунок 5: ошибки Прогноза получены с предполагаемой классической моделью добычи хищника IDNLGREY.
C1 . Система добычи хищника с давкой добычи
Последнее исследование генеральной совокупности также посвящено 1 хищнику и 1 системе добычи, с различием, являющимся этим, мы здесь вводим термин-p5*x2 (t) ^2 представление промедления роста добычи из-за давки. Повторно параметризованная образцовая структура от предыдущего примера затем только дополнена с этим сроком давки:
d
-- x1(t) = p1*x1(t) + p2*x2(t)*x1(t)
dt
d
-- x2(t) = p3*x2(t) - p4*x1(t)*x2(t) - p5*x2(t)^2
dt
y1(t) = x1(t)
y2(t) = x2(t)
Интерпретация этих уравнений является по существу тем же самым как выше, за исключением того, что в отсутствие хищников рост генеральной совокупности добычи будет держаться в страхе.
Новая ситуация с моделированием отражается файлом MEX на C, названным predprey2_c.c, который является почти тем же самым как predprey1_c.c. Кроме изменений, связанных с количеством параметров модели (5 вместо 4), основное различие заключается в функции обновления состояния compute_dx:
/* State equations. */
void compute_dx(double *dx, double t, double *x, double **p,
const mxArray *auxvar)
{
/* Retrieve model parameters. */
double *p1, *p2, *p3, *p4, *p5;
p1 = p[0]; /* Survival factor, predators. */
p2 = p[1]; /* Death factor, predators. */
p3 = p[2]; /* Survival factor, preys. */
p4 = p[3]; /* Death factor, preys. */
p5 = p[4]; /* Crowding factor, preys. */
/* x[0]: Predator species. */
/* x[1]: Prey species. */
dx[0] = p1[0]*x[0]+p2[0]*x[1]*x[0];
dx[1] = p3[0]*x[1]-p4[0]*x[0]*x[1]-p5[0]*pow(x[1],2);
}
Заметьте, что добавленный срок промедления вычисляется как p5[0]*pow (x [1], 2). Степень функциональной головы (..) задан в библиотеке C math.h, который должен быть включен наверху образцового файла через оператор #include "math.h" (это не необходимо в predprey1_c.c, когда это только содержит стандарт C математика).
Скомпилированный файл MEX на C, наряду с начальным вектором параметра, соответствующим начальным состоянием и некоторой административной информацией затем питается как входные параметры конструктора Object IDNLGREY:
FileName = 'predprey2_c'; % File describing the model structure. Order = [2 0 2]; % Model orders [ny nu nx]. 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}); % Estimate all 5 parameters. InitialStates = struct('Name', {'Population, predators' 'Population, preys'}, ... 'Unit', {'Size (in thousands)' 'Size (in thousands)'}, ... 'Value', {1.8 1.8}, ... 'Minimum', {0 0}, ... 'Maximum', {Inf Inf}, ... 'Fixed', {false false}); % Estimate both initial states. Ts = 0; % Time-continuous system. nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts, ... 'Name', '1 predator - 1 prey system exhibiting crowding', ... 'OutputName', {'Population, predators', 'Population, preys'}, ... 'OutputUnit', {'Size (in thousands)' 'Size (in thousands)'}, ... 'TimeUnit', 'year');
Путем введения имени объекта модели (nlgr) основная информация о модели отображен в командном окне. Обратите внимание на то, что, как прежде, подарок (nlgr) предоставляет более всесторонние сводные данные модели.
nlgr
nlgr =
Continuous-time nonlinear grey-box model defined by 'predprey2_c' (MEX-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).
Name: 1 predator - 1 prey system exhibiting crowding
Status:
Created by direct construction or transformation. Not estimated.
C2 . Данные ввода - вывода
Затем мы загружаем (моделируемый, хотя поврежденный шум) данные, и создайте объект IDDATA, описывающий этот тип давки ситуации добычи хищника. Этот набор данных содержит 201 выборку данных, покрывающую 20 лет эволюции.
load(fullfile(matlabroot, 'toolbox', 'ident', 'iddemos', 'data', 'predprey2data')); z = iddata(y, [], 0.1, 'Name', '1 predator - 1 prey system exhibiting crowding'); set(z, 'OutputName', {'Population, predators', 'Population, preys'}, ... 'Tstart', 0, 'TimeUnit', 'Year');
C.3. Производительность первоначальной модели добычи хищника с давкой добычи
Симуляция с первоначальной моделью ясно показывает, что это не может справиться с истинной демографической динамикой. Смотрите фигуру.
clf compare(z, nlgr, 1);
Рисунок 6: Сравнение между истинными выходными параметрами и моделируемыми выходными параметрами первоначальной модели добычи хищника с давкой добычи.
C.4. Оценка параметра
Чтобы улучшать производительность первоначальной модели, мы продолжаем оценивать 5 неизвестных параметров и эти 2 начальных состояния.
nlgr = nlgreyest(z, nlgr, nlgreyestOptions('Display', 'on'));
C.5. Производительность предполагаемой модели добычи хищника с давкой добычи
Ориентировочные стоимости параметров и начальных состояний снова вполне близко к тем, которые использовались, чтобы сгенерировать истинные выходные данные:
disp(' True Estimated parameter vector');True Estimated parameter vector
ptrue = [-1; 1; 1; 1; 0.1];
fprintf(' %6.3f %6.3f\n', [ptrue'; getpvec(nlgr)']); -1.000 -1.000
1.000 1.001
1.000 1.002
1.000 1.002
0.100 0.101
disp(' ');disp(' True Estimated initial states');True Estimated initial states
x0true = [2; 2]; fprintf(' %6.3f %6.3f\n', [x0true'; cell2mat(getinit(nlgr, 'Value'))']);
2.000 2.003
2.000 2.002
Чтобы далее оценить качество модели (и проиллюстрировать улучшение по сравнению с первоначальной моделью), мы также моделируем предполагаемую модель. Моделируемые выходные параметры сравниваются с истинными выходными параметрами в окне графика. Как видно снова, предполагаемая модель довольно хороша.
compare(z, nlgr, 1);
Рисунок 7: Сравнение между истинными выходными параметрами и моделируемыми выходными параметрами предполагаемой модели добычи хищника с давкой добычи.
Мы завершаем третий пример генеральной совокупности путем представления информации модели, возвращенной PRESENT.
present(nlgr);
nlgr =
Continuous-time nonlinear grey-box model defined by 'predprey2_c' (MEX-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) Population, predators(t) [Size (in t..] xinit@exp1 2.00281 (estimated) in [0, Inf]
x(2) Population, preys(t) [Size (in t..] xinit@exp1 2.00224 (estimated) in [0, Inf]
Outputs:
y(1) Population, predators(t) [Size (in thousands)]
y(2) Population, preys(t) [Size (in thousands)]
Parameters: Value Standard Deviation
p1 Survival factor, predators [1/year] -0.999914 0.00280581 (estimated) in [-Inf, Inf]
p2 Death factor, predators [1/year] 1.00058 0.00276684 (estimated) in [-Inf, Inf]
p3 Survival factor, preys [1/year] 1.0019 0.00272154 (estimated) in [-Inf, Inf]
p4 Death factor, preys [1/year] 1.00224 0.00268423 (estimated) in [-Inf, Inf]
p5 Crowding factor, preys [1/year] 0.101331 0.0005023 (estimated) in [-Inf, Inf]
Name: 1 predator - 1 prey system exhibiting crowding
Status:
Termination condition: Change in cost was less than the specified tolerance..
Number of iterations: 8, Number of function evaluations: 9
Estimated using Solver: ode45; Search: lsqnonlin on time domain data "1 predator - 1 prey system exhibiting crowding".
Fit to estimation data: [97.53;97.36]%
FPE: 4.327e-08, MSE: 0.0004023
More information in model's "Report" property.
Заключения
Этот пример показал, как выполнить моделирование timeseries IDNLGREY на основе MATLAB и файлов модели MEX.