В этом примере показано, как создавать нелинейные модели временных рядов серых полей. Модели временных рядов представляют собой модели, не использующие никаких измеренных входных данных. Изучаются три идеализированные экологические системы, в которых два вида:
конкурировать за одну и ту же еду, или:
находятся в ситуации хищника-добычи
В примере показано моделирование на основе файлов MATLAB и C MEX.
Во всех трех исследованных популяционных системах мы интересны тем, как популяция двух видов меняется со временем. Чтобы смоделировать это, пусть 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 погибают каждый раз. Чистый эффект этих операторов можно суммировать в типе пространства состояний структуры модели: (временной ряд):
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).
В случае, если два вида конкурируют за одну и ту же пищу, то именно общая популяция вида контролирует доступность пищи, и, в свою очередь, их смертность. Простой, но распространенный подход заключается в предположении, что уровень смертности может быть записан как:
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.
Затем мы загружаем (смоделированные, хотя и поврежденные шумом) данные и создаем объект 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');
Моделирование с исходной моделью ясно показывает, что она не может справиться с истинной динамикой населения. См. рисунок графика. Для типа временного ряда модели IDNLGREY обратите внимание, что вывод модели определяется начальным состоянием.
compare(z, nlgr, 1);
Рис. 1: Сравнение истинных результатов и смоделированных результатов начальной модели двух видов.
Чтобы преодолеть довольно плохую производительность начальной модели, мы приступаем к оценке 4 неизвестных параметров и 2 начальных состояний с использованием NLGREYEST. Укажите параметры оценки с помощью NLGREYESTOPTIONS; в этом случае «Display» имеет значение «on», что означает, что информация о ходе оценки отображается в окне хода выполнения. С помощью команды NLGREYESTOPTIONS можно задать основные свойства алгоритма, такие как «GradityOptions», «SeireMethod», «MaxIterations», «Tolerance», «Display».
opt = nlgreyestOptions;
opt.Display = 'on';
opt.SearchOptions.MaxIterations = 50;
nlgr = nlgreyest(z, nlgr, opt);Оценочные значения параметров и начальных состояний хорошо соответствуют значениям, используемым для формирования истинных выходных данных:
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: Сравнение истинных результатов и смоделированных результатов расчетной модели двух видов.
Настоящее изобретение предоставляет дополнительную информацию об оценочной модели, например, о неопределенностях параметров и других связанных с оценкой величинах, таких как функция потерь и мера 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.
Предположим, что первый вид живет на втором. Доступность пищи для вида 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)
который лучше подходит с точки зрения оценки.
На этот раз мы введем эту информацию в C MEX-файл с именем predprey1_c.c. Файл модели структурирован как стандартный файл IDNLGREY C MEX (см. пример «Создание файлов модели 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];
}
Поскольку модель имеет тип временного ряда, ни compute_dx, ни compute_y не включают u во входной список аргументов. Фактически, основная интерфейсная функция predprey1_c.c, даже не объявляет u, даже если пустой u ([]) всегда передается predprey1_c методами IDNLGREY.
Скомпилированный C MEX-файл вместе с исходным вектором параметров, адекватным начальным состоянием и некоторой административной информацией затем подаются в качестве входных аргументов конструктору объекта 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.
Наш следующий шаг - загрузить (смоделированные, хотя и поврежденные шумом) данные и создать объект 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');
Моделирование с исходной моделью указывает на то, что она не может точно справиться с истинной динамикой населения. См. окно печати.
compare(z, nlgr, 1);
Рис. 3: Сравнение между истинными выходами и моделируемыми выходами исходной классической модели хищник-добыча.
Чтобы улучшить производительность начальной модели, мы продолжаем оценивать 4 неизвестных параметра и 2 начальных состояния с использованием NLGREYEST.
nlgr = nlgreyest(z, nlgr, nlgreyestOptions('Display', 'on'));
Оценочные значения параметров и начальных состояний очень близки к тем, которые использовались для генерации истинных выходных данных:
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.
Последнее исследование популяции также посвящено 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)
Интерпретация этих уравнений по сути та же, что и выше, за исключением того, что в отсутствие хищников рост популяции добычи будет удерживаться в страхе.
Новая ситуация моделирования отражена C MEX-файлом под названием 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).
Скомпилированный C MEX-файл вместе с исходным вектором параметров, адекватным начальным состоянием и некоторой административной информацией затем подаются в качестве входных аргументов конструктору объекта 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.
Затем мы загружаем (смоделированные, хотя и поврежденные шумом) данные и создаем объект 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');
Моделирование с исходной моделью ясно показывает, что она не может справиться с истинной динамикой населения. Смотрите рисунок.
clf compare(z, nlgr, 1);
Рис. 6: Сравнение между истинными выходами и моделируемыми выходами исходной модели хищник-добыча с переполненностью добычи.
Для улучшения производительности начальной модели мы приступаем к оценке 5 неизвестных параметров и 2 начальных состояний.
nlgr = nlgreyest(z, nlgr, nlgreyestOptions('Display', 'on'));
Оценочные значения параметров и начальных состояний снова довольно близки к тем, которые использовались для генерации истинных выходных данных:
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(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.
В этом примере показано, как выполнять моделирование временных рядов IDNLGREY на основе файлов модели MATLAB и MEX.