Три экологических системы населения: 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

Здесь естественно выбрать два состояния в качестве выходных параметров, i.e., мы позволяем 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);

Figure contains 2 axes objects. Axes object 1 contains 2 objects of type line. These objects represent Two species competing for the same food (Population, species 1), Two species competing for the same food: 67.36%. Axes object 2 contains 2 objects of type line. These objects represent Two species competing for the same food (Population, species 2), Two species competing for the same food: 67.83%.

Рисунок 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);

Figure contains 2 axes objects. Axes object 1 contains 2 objects of type line. These objects represent Two species competing for the same food (Population, species 1), Two species competing for the same food: 98.42%. Axes object 2 contains 2 objects of type line. These objects represent Two species competing for the same food (Population, species 2), Two species competing for the same food: 97.92%.

Рисунок 2: Сравнение между истинными выходными параметрами и симулированными выходными параметрами приблизительно двух моделей разновидностей.

PRESENT предоставляет дополнительную информацию о предполагаемой модели, e.g., о неопределенности параметра и других связанных с оценкой количествах, как функция потерь и 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 увеличится когда количество индивидуумов первых увеличений разновидностей, e.g., согласно:

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

Как в предыдущем примере населения, также здесь невозможно однозначно определить шесть отдельных параметров. С тем же видом репараметризации как в вышеупомянутом случае, i.e., 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);

Figure contains 2 axes objects. Axes object 1 contains 2 objects of type line. These objects represent Classical 1 predator - 1 prey system (Population, predators), Classical 1 predator - 1 prey system: 4.842%. Axes object 2 contains 2 objects of type line. These objects represent Classical 1 predator - 1 prey system (Population, preys), Classical 1 predator - 1 prey system: 4.523%.

Рисунок 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);

Figure contains 2 axes objects. Axes object 1 contains 2 objects of type line. These objects represent Classical 1 predator - 1 prey system (Population, predators), Classical 1 predator - 1 prey system: 98.08%. Axes object 2 contains 2 objects of type line. These objects represent Classical 1 predator - 1 prey system (Population, preys), Classical 1 predator - 1 prey system: 97.97%.

Рисунок 4: Сравнение между истинными выходными параметрами и симулированными выходными параметрами предполагаемой классической модели добычи хищника.

Как ожидалось ошибки предсказания, возвращенные PE, малы и случайной природы.

pe(z, nlgr);

Figure contains 2 axes objects. Axes object 1 contains an object of type line. These objects represent Classical 1 predator - 1 prey system (Population, predators), Classical 1 predator - 1 prey system. Axes object 2 contains an object of type line. These objects represent Classical 1 predator - 1 prey system (Population, preys), Classical 1 predator - 1 prey system.

Рисунок 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);

Figure contains 2 axes objects. Axes object 1 contains 2 objects of type line. These objects represent 1 predator - 1 prey system exhibiting crowding (Population, predators), 1 predator - 1 prey system exhibiting crowding: 62.34%. Axes object 2 contains 2 objects of type line. These objects represent 1 predator - 1 prey system exhibiting crowding (Population, preys), 1 predator - 1 prey system exhibiting crowding: 46.38%.

Рисунок 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);

Figure contains 2 axes objects. Axes object 1 contains 2 objects of type line. These objects represent 1 predator - 1 prey system exhibiting crowding (Population, predators), 1 predator - 1 prey system exhibiting crowding: 97.53%. Axes object 2 contains 2 objects of type line. These objects represent 1 predator - 1 prey system exhibiting crowding (Population, preys), 1 predator - 1 prey system exhibiting crowding: 97.36%.

Рисунок 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.

Для просмотра документации необходимо авторизоваться на сайте