Регулирование многовидового рыболовства на основе математического моделирования (стр. 20 )

Для оценки начальных приближений параметров модели систему уравнений (4-1) рассматриваем как регрессионные соотношения, коэффициенты которых должны быть определены, а исходными данными для этого служат результаты последнего этапа подготовки рядов данных предыдущего раздела. Всего следует определить 8 параметров.

Эту процедуру оценки параметров можно выполнить без каких-либо отклонений от приведенного алгоритма только в случае независимости троек "гладких" переменных {u1, u2, E1} и {u1, u2, E2}. Однако известно, что во многих случаях, когда улов на усилие является индексом запаса, существует достаточно сильная линейная связь между усилием с одной стороны и уловом на усилие – с другой. Для нашего примера были построены уравнения линейной регрессии вида

E1=g1+g2u1+g3u2 и E2=h1+h2u1+h3u2

по всем 17 точкам (с 1974 по 1982 г. с полугодовым шагом), коэффициенты множественной корреляции для них составили 0,85 и 0,9 соответственно.

Чтобы повысить определенность задачи, следует оценить коэффициенты улавливаемости q1 и q2 каким-либо независимым методом. Для рассматриваемых запасов ставриды и капского хека оценки запасов этих популяций в районах 1.3+1.4. ИКСЕАФ (Babayan and Bulgakova, 1983a; 1983b; Bulgakova et al., I984; Babayan et al, 1985; ICSEAF, 1983) были сопоставлены с уловами на усилие, в результате получены приближенные оценки: q1= 0,9*10-6 и q2 = 2,5*10-6. Эти значения и использованы в дальнейшем. В тексте опускаем порядковые множители 10-6, поскольку в уравнения переменная qi входит в произведение qiEi, а переменная Ei имеет порядковый множитель 106.

Поскольку в реальных запасах с многовозрастной структурой не может быть очень резких локальные изменений запаса во времени, решено ограничиться расчетами коэффициентов регрессионных уравнений по 13 точкам (полугодовые интервалы) – от середины 1976 г. до середины 1982 г. Последнее полугодие использовалось для тестирования прогностических возможностей модели. На этапе тестирования значения параметров уточнялись.

Рис. 4.3. Фазовая диаграмма системы запасов ставрида - хек: I − сглаженные фактические данные; 2 – результат тестирования модели (4-1); тонкие пунктирные линии – изолинии равновесных уловов (в I06 т в год); крестики и кружки, как на рис. 4.1.

Начальные приближения параметров a1 , b1, c1, a2, b2, c2 при заданных q1 и q2 определяются с помощью программы линейной множественной регрессии из следующих уравнений:

(4-2)

Расчеты, однако, показали, что коэффициент множественной корреляции для первого из соотношений (4-2) (около 0,75) меньше коэффициента корреляции между левой частью этого соотношения и u2 (равного 0,85). Другими словами, в данном случае следует заранее положить b1=0 и строить линейную регрессию величин (du1/(u1dt)+g1E1) от одного фактора и2. Это соображение подтверждается еще и тем, что при учете обоих факторов величина b1 мала и плохо обеспечена.

Для второго соотношения (4-2) коэффициент множественной корреляции достаточно высок (около 0,8), причем исключение любого из факторов снижает величину коэффициента корреляции.

Итак, описанные процедуры дают следующие начальные приближения параметров модели:

a1= I,76; b1=0,00; c1=1,00; a2= 0,56; b2=0,34; c2= 0,18,

причем t-критерий для c1, b2 и c2 достаточно высок и превышает 5, обеспеченность же параметров a1 и a2 всегда высокая. Отметим, что коэффициенты уравнений (4-2) (кроме b1) получились положительными, так что все слагаемые входят в уравнения именно с теми знаками, которые соответствуют начальным представлениям о внутри - и межвидовых отношениях. Кроме того, фактор внутривидовой плотностной регуляции запаса оказался отчетливо выраженным именно для хека, поскольку b2>b1. Значимость этого параметра подтверждает, что параметры модели имеют биологический смысл – ведь каннибализм хека в этом районе – широко известное явление (Lleonart et al., 1985; Konchina, 1986; Crawford et al., 1987). Параметры c1 и c2 получены положительными, что означает, что хищничество мерлузы по отношению к ставриде является более сильным фактором, чем конкуренция. Дело в том, что некоторые ученые (Krzeptowski, 1982) указывали на совпадение состава пищи капского хека и ставриды, основой которой являются эвфаузииды и миктофиды, и на этом основании делали вывод о конкуренции между ставридой и хеком. По данным же (Konchina, 1986), в 1985 г. ставрида составляла до трети пищи хека длиной 24-44 см и до 100% пищи хека длиной 32-40 см, что говорит в пользу более сильного влияния хищничества.

Теперь проводим тестирование модели и уточняем ее параметры. Если задаться зависимостью E1 и E2 от времени, то по уравнению (4-1) можно рассчитывать изменения u1 и u2 (а, значит, и запасов) от некоторого их начального состояния на определенный срок. Проверку качества описания моделью динамики запасов (сглаженных значений u1 и u2) проводим, используя уже найденные представления E1 и E2 в виде полиномов четвертой степени от времени. Решение системы дифференциальных уравнений (4-1) проводим методом Рунге-Кутта, приняв в качестве начальных значений значения, соответствующие середине 1976г., расчеты велись до конца рассматриваемого периода.

Для уточнения значений коэффициентов системы уравнений (4-1) их начальные значения округляли до пяти сотых (a1= 1,75, b1= 0,00, c1= 1,00, a2= 0.55, b2= 0.35, c2= 0.20), а затем значения коэффициентов меняли в некоторых диапазонах с шагом 0,025. При этом для каждой комбинации значений коэффициентов по уравнениям (4-1) численно (с шагом 0,05 года) определяли величины u1 и u2. Различные варианты сравнивали по среднеквадратическому отклонению расчетных величин du1/du2 в исходных точках от сглаженных фактических значений; можно и визуально сравнивать их соответствие на фазовой диаграмме (рис. 4.3), где приведены сглаженные фактические данные и рассчитанные траектории. В конкретном примере наилучшей оказалась комбинация: a1=I,70, b1= 0,05, c1= 0,9, a2=0,5, b2=0,35, c2=0,2.

Таким образом, получена окончательная система уравнений для описания динамики запасов ставриды и хека:

(4-3)

Результаты тестирования модели (4-3) (на рис.4.3 кривые 1 и 2) показали, что модель демонстрирует способность описывать главные закономерности в динамике взаимодействующих запасов. Отметим, что в окончательном варианте параметр b1 больше нуля, хотя и значительно меньше b2.

Использование стандартных программ избавляет от необходимости разработки специальной программы целевого назначения и имеет ряд преимуществ, в том числе позволяет работать в интерактивном режиме, легко меняя вариант расчетов – тип сглаживания данных, выбор целевой функции, способ определения параметров и т. д.

4.4. Анализ чувствительности модели к вариациям параметров

Настройка модели и исследование модели на чувствительность к изменениям параметров – задачи, тесно связанные между собой. В обоих случаях необходимо провести серию численных экспериментов, варьируя параметры модели. При изучении чувствительности модели следует проанализировать различные прогоны, чтобы выяснить, варьирование каких параметров более всего влияет на поведение модели, т. е. какие параметры следует оценивать с большой точностью.

Поведение данной динамической модели характеризуется двумя выходными величинами, которые зависят от времени – u1(t) и u2(t). Поэтому в качестве показателя чувствительности модели к изменению какого-либо из ее параметров (обозначим его через а ) выбираем следующую величину

. (4-4)

Имеется в виду, что при изменении значения выбранного параметра от a (1) до a (2) величина u1(t) для некоторого фиксированного значения t изменится от u1(1) до u1(2), a u2(t) - от u2 (1) до u2(2). На фазовой плоскости {u1,u2} значению параметра a (1) соответствует точка (u1(1),u2(2)), а значению a(2) – точка (u1(2),u2(2)). Тогда – расстояние между этими точками, а радикал в знаменателе равен длине вектора {u1,u2}, соответствующего средней точке между двумя вышеуказанными; а aср=(a(1)+a(2))/2 – среднее значение параметра. Таким образом, выражение (4-4) представляет собой отношение относительного изменения выходной величины к относительному изменению параметра a. При таком подходе остальные параметры следует оставлять неизменными.

Очевидно, что чувствительность модели в смысле (4-4) зависит от момента времени t. Целесообразно рассчитывать чувствительность в тех точках, где наблюдается наибольший разброс решений уравнений (4-1) при варьировании параметра или где решение значительно отклоняется от "реальных" данных.

При настройке модели было сделано около 30 прогонов – решений дифференциальных уравнений (4-1) с различными наборами параметров, но при одних и тех же начальных условиях. Среди этих экспериментов каждому параметру модели соответствуют 2-3 варианта расчетов, в которых меняется только этот параметр. Именно они выбраны для расчетов чувствительности к изменению параметров ai; bi , сi (i= 1,2). Эти эксперименты соответствовали набору исходных данных с 1975 по 1984 г. Для расчетов чувствительности к варьированию q1 и q2 были использованы исходные данные с 1974 по 1983 г., и при каждом сочетании q1 и q2 методом множественной регрессии однозначно оценивались все остальные параметры.