Три экологические системы населения: MATLAB и C MEX-файловое моделирование файла 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 вместе с вектором начального параметра, адекватным начальным состоянием и некоторой административной информацией подаются в качестве входов в конструктор объекта IDNLGREY. Заметьте, что начальные значения параметров и начальных состояний заданы как массивы структур с элементами Npo (количество объектов параметра = количество параметров, если все параметры скаляры) и Nx (количество состояний), соответственно. Через эти массивы структур можно полностью присвоить значения свойств, отличные от свойств по умолчанию ('Name', 'Unit', 'Value', 'Minimum', 'Maximum' и 'Fixed') каждому параметру и начальному состоянию. Здесь мы присвоили значение 'Minimum' каждого начального состояния нулю (населения положительны!), а также указали, что оба начальных состояния должны быть оценены по умолчанию.
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; в этом случае 'Display' устанавливается на ' on ', что означает, что информация о ходе оценки отображается в окне прогресса. Можно использовать NLGREYESTOPTIONS, чтобы задать основные свойства алгоритма, такие как 'GradientOptions', 'SearchMethod', 'MaxIterations', 'Tolerance', 'Display'.
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 Акайке (Final Prediction Error).
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. Файл модели структурирован как стандартный файл MEX на C IDNLGREY (см. Пример, озаглавленный «Создание файлов модели 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 вместе с вектором начального параметра, адекватным начальным состоянием и некоторой административной информацией далее подаются как входные параметры на конструктор объекта 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.
C.1. Система хищник-добыча с переполненностью добычи
Последнее население исследование также посвящено 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). Степень функции pow (.,.) задана в библиотеке C math.h, которая должна быть включена в верхней части файла модели через оператор # include «math.h» (это не обязательно в predprey1_c.c, так как она содержит только стандартную математику C).
Скомпилированные файлы MEX на C вместе с вектором начального параметра, адекватным начальным состоянием и некоторой административной информацией далее подаются как входные параметры на конструктор объекта 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.
C.2. Входные-выходные данные
Далее мы загружаем (моделируемые, хотя и поврежденные шумом) данные и создаем объект 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.