Метод множителей Лагранжа

Лекция 7 При выборе оптимального варианта действий (оптимального управления) в задачах, связанных с проектированием объектов всегда задаются ограничения, накладываемые на искомые оптимальные значения параметров или характеристик объектов. Эти ограничения задаются, как правило, функциями. Рассмотрим общий подход к решению подобных задач, а именно, метод множителей Лагранжа. Суть метода состоит в следующем. Пусть требуется найти экстремум функции при ограничениях . Составляется функция: (6.3) Эту функцию называют функцией Лагранжа, а числа  — множителями Лагранжа. Далее можно считать, что все переменные в поставленной задаче независимы, т. е. можно применить теорему Ферма. Следовательно, имеем систему уравнений: (6.4) где  — число переменных (размерность вектора );  — число функций, определяющих ограничения. Число уравнений может быть сокращено на одно уравнение за счет того, что значение может быть принято равным единице. Тогда метод множителей Лагранжа сводится к алгоритму: Шаг 1. Формализация задачи. Составление уравнения (6.3). Шаг 2. Составление системы уравнений (6.4). Шаг 3. Нахождение всех стационарных точек. При этом следует заметить, что поскольку все могут быть умножены на любое число, то значение  можно принять равным единице. Шаг 4. Перебор всех стационарных точек и выбор тех, в которых функция (6.3) достигает максимума или минимума в зависимости от существа решаемой задачи. Пример 6.1 Пусть требуется сконструировать резервуар, состоящий из днища, крыши и боковых стенок, объемом из листовой стали. Толщина листов и плотность стали фиксированы. Необходимо выбрать такую конструкцию, для которой расход стали будет минимальным. Масса расхода стали на конструкцию: (6.5) где  — площадь поверхности резервуара. Таким образом, достаточно выбрать конструкцию с минимальным значением площади поверхности , что будет адекватно выбору конструкции с минимальным расходом металла. Пусть рассматриваются два варианта резервуаров: прямоугольной и цилиндрической формы. Для резервуара в форме прямоугольного параллелепипеда площадь поверхности определяется по формуле: (6.6) где  — высота параллелепипеда;  — стороны параллелепипеда. Требуется найти такие значения , которые доставят минимум функции (6.6). Составим функцию Лагранжа: (6.7) Составляем систему уравнений: (6.8) Решение составленной системы уравнений, произведенное с помощью пакета Mathcad, приведено ниже в протоколе 6.2. Опуская мнимые и отрицательные корни, получим, что минимальным затратам материала в данном случае отвечает куб со стороной: Зададимся значением  м3, тогда получим площадь поверхности куба: (6.9) Затраты металла при толщине стенки 1 см и плотности стали 7800 кг/м3 составят: Рассмотрим второй вариант конструкции, а именно, цилиндр, высотой  и радиусом . Необходимо выбрать такие и , которые доставят минимум площади поверхности цилиндра: При ограничениях: (6.10) Составим функцию Лагранжа: (6.11) Решение поставленной задачи в Mathcad приведено в протоколе 6.3. Для значения  м3, получим значение площади поверхности цилиндрического резервуара: высота  м, радиус  м, суммарная площадь поверхности резервуара: (6.12) Сравнив значения (6.9) и (6.12), приходим к выводу, что резервуар в форме цилиндра по критерию расхода стали выгоднее. При толщине стенки 0,01 м и плотности 7800 кг/м3, получим: масса металла для резервуара в форме прямоугольного параллелепипеда составит 136812 кг, а для резервуара в форме цилиндра — 126282 кг. Таким образом, разность в затратах стали составит 12420 кг. Протокол 6.2. Решение системы уравнений 6.8, произведенное с помощью пакета Mathcad Окончание протокола 6.3 Протокол 6.3. Решение задачи об оптимальных размерах резервуара цилиндрической формы в Mathcad Сформулированная выше задача с физической точки зрения страдает существенным недостатком. Если рассматриваемый резервуар предназначен для хранения жидкости, например, нефти, то не учтено основное ограничение. Кроме заданного объема необходимо обеспечить достаточную прочность резервуара. Для цилиндрических резервуаров равной толщины стенок на основании уравнения Лапласа толщина стенки должна быть не менее : (6.13) где  — плотность хранимой жидкости, кг/м3;  — ускорение силы тяжести, м/с2;  — радиус цилиндра, м;  — высота цилиндра, м;  — допустимое напряжение материала стенок резервуара, Па. С учетом ограничения (6.13) задача проектирования цилиндрического резервуара с равной толщиной металла формулируется следующей постановкой. Найти такие значения радиуса цилиндра, его высоты и толщины металла , для которых масса затраченного металла будет минимальна при ограничениях, накладываемых на заданный объем резервуара и на прочность резервуара, определяемую (6.13). Модель цилиндрического резервуара, его площади поверхности при толщине стенок, крышки, и днища равных , по условию удовлетворения прочности для хранения жидкости плотностью имеет вид: Формализованная постановка задачи имеет вид: (6.14) Решение задачи (6.14) для стального резервуара с плотностью стали 7800 кг/м3, динамической прочностью 200 МПа, предназначенного для хранения 5000 м3 нефти плотностью 850 кг/м3, выполненное в Mathcad (протокол 6.4) дает следующие результаты: – радиус резервуара 14,71 м, высота 7,35 м, – толщина стали 4,5 мм, – объем затраченного металла составляет 11,19 м3, – масса стали 71,57 т, – расчетный объем резервуара при этом равен 4996 м3. Различие результатов решения задачи в постановке (6.11) и (6.14) объясняется изменением ограничений на неравенства и введением третьего аргумента — минимальной толщины стенки металла, обеспечивающего прочность резервуара. Кроме того, если в задаче (6.11) толщина стали выбрана наудачу (1 см), то в задаче (6.14) ее значение (0,45 см) получено в результате решения оптимизационной задачи. Расхождение значений параметров цилиндрических резервуаров при выбранной наудачу толщине металла и полученной в результате решения, объясняют различием моделей поставленных задач по их физическому смыслу. Рассмотрим еще одну нелинейную модель определения параметров. Пример 6.2 Рассматривается система с нагрузочным резервированием и простейшим потоком отказов. Протокол 6.4. Решение задачи определения оптимальной конструкции резервуара цилиндрической формы с учетом прочности в Mathcad Техническая система состоит из одного основного и трех одинаковых с основным резервных устройств (рис. 6.2). Основное устройство подвергается простейшему потоку отказов с интенсивностью . Каждое из резервных устройств до своего включения в работу на полную нагрузку подвергается простейшему потоку отказов с интенсивностью и после включения в работу эта интенсивность для соответствующего устройства мгновенно «подскакивает» до значения . Система функционирует по алгоритму: после отказа основного устройства включается в работу устройство , затем и т. д. Отказ системы адекватен отказу последнего из всех устройств. Переключающее устройство принимается идеальным. Требуется определить надежность технической системы. Рис. 6.2. Вид технической системы с нагрузочным резервированием, состоящей из одного основного и трех резервных устройств Состояние системы будем нумеровать двумя индексами. Первый индекс будем считать равным единице, если основное устройство находится в работоспособном состоянии, и равным нулю, если основное устройство отказало. Второй индекс будем принимать равным числу исправных резервных устройств. Тогда состояния системы (см. рис. 6.3) принимают следующие обозначения:  — в исправном состоянии находятся все четыре устройства технической системы (основное и три резервных);  — исправно основное устройство, из трех резервных устройств одно отказало, два находятся в исправном состоянии;  — основное устройство исправно и исправно одно резервное;  — исправно только основное устройство, все три резервных отказали;  — основное устройство отказало, все три резервных исправны;  — основное устройство отказало, исправны два резервных устройства;  — основное устройство отказало, работает только одно из резервных;  — отказали все четыре устройства (произошел отказ системы). Рис. 6.3 Размеченный граф состояний технической системы с «нагруженным» резервированием Система уравнений Колмогорова для вероятностей состояний рассматриваемой системы имеет вид: (6.15) К этим уравнениям необходимо добавить условие: (6.16) которое позволяет исключить одно (любое) из уравнений Колмогорова. Интегрирование системы уравнений (6.15) может быть осуществлено в следующей последовательности: из первого уравнения находится функция , равная: (6.17) Это выражение подставляется во второе уравнение, которое теперь содержит только одну неизвестную функцию , которую достаточно легко можно определить и подставить в третье уравнение, и т. д. На каждом шаге такого процесса новые выражения функций мы находим через уже известные, пока не доходим до функции , которую выражаем через все остальные: (6.18) Теперь можно определить надежность рассматриваемой технической системы: (6.19) Пример изменения вероятностей нахождения системы в перечисленных выше состояниях, рассчитанных с помощью пакета Mathсad, приведен на рис. 6.4. На рис. 6.4 введены обозначения (решение системы уравнений): , , , , , , , Расчеты проводились при , , . Решение задачи в среде Mathcad приведено в протоколе 6.5, где интенсивности отказов имеют иные обозначения для удобства программирования. Пояснение к решению. По оси абсцисс рис. 6.4 отложено значение , определяющее число шагов интегрирования системы и входящее аргументом встроенной функции . Весь рассматриваемый интервал времени (1 год), определенный вторым и третьим аргументом функции разделен на 100 интервалов (при отсчете от 0 имеем значение 99). Для получения значения времени , соответствующего значению , необходимо в данном случае выполнить действия: ( 6.20) Если бы рассматривался интервал времени продолжительностью 5 лет: (6.21) где  — значение, выбранное на оси абсцисс рис. 6.4, а 100 — заданное предельное значение изменения . Рис. 6.4. Изменение вероятностей состояния системы Кривые на рис. 6.4 — функции изменения вероятности нахождения системы в соответствующем состоянии. К задачам нелинейного программирования часто относят задачи вариационного исчисления. Отличие задач вариационного исчисления состоит в том, что в роли управления выступает не число, а функция, т. е. значения целевой функции , которую теперь называют функционалом, зависят от выбранного управления , где множество представляет собой множество функций. Одной из старинных задач, послужившей началом в развитии вариационного исчисления, является задача о брахистохроне. Задача о брахистохроне В июньском номере журнала «Акта Эрудиторум» за 1696 год была помещена заметка знаменитого швейцарского ученого Иоганна Бернулли. Вот как формулирует постановку задачи сам И. Бернулли. В вертикальной плоскости даны точки А и В (рис. 6.5). Определить путь АМВ, спускаясь по которому под действием собственной силы тяжести, тело М, начав двигаться из точки А, достигнет точки В в кратчайшее время. Рис. 6.5. Типовой пример вариационной задачи. Задача о брахистохроне Решение поставленной задачи И. Бернулли произвел следующим образом. Пусть уравнение кривой, по которой должно двигаться тело , задано функцией . Вычисляется время, за которое тело массой (без учета трения) скатится из точки в точку по желобу . Согласно закону Галилея скорость тела в точке с координатами при движении только под действием силы тяжести без учета трения зависит лишь от ординаты и не зависит от формы кривой между точками и . Кинетическая энергия тела в точке равна и эта энергия равна разности потенциальных энергий в рассматриваемых точках. В итоге, скорость в точке оказывается равной , где  — ускорение силы тяжести. Далее рассматривается участок пути между точками: где  — сколь угодно малое приращение абсциссы. Длина этого участка пути приближенно равна: Воспользовавшись приближенным равенством: получается: Полагая скорость движения на малом участке постоянной и равной в итоге, время, необходимое для прохождения малого участка окажется приближенно равным: Окончательно, общее время прохождения телом пути от точки до точки по кривой запишется в виде интеграла: (6.22) где верхний предел интегрирования «» представляет собой координату по оси ординат конечного движения тела , а именно, координату точки . Таким образом, получается следующая аналитическая форма задачи о брахистохроне. Требуется найти такую функцию , для которой значение  будет минимальным. При этом множество функций, претендующих на участие в решении задачи, ограничено, эти функции должны удовлетворять следующим условиям: , . Решение подобных задач относится к классу вариационного исчисления. Кратко остановимся на сущности вариационного исчисления. Пусть задана некоторая непрерывная функция трех переменных . Рассмотрим следующий функционал: (6.23) Будем рассматривать этот функционал в пространстве . Тогда можно дать следующее определение. Вариационным исчислением называется раздел теории экстремальных задач, где изучаются максимумы и минимумы таких функционалов при различного рода ограничениях. Чтобы по функции из получить число необходимо продифференцировать . Затем подставить во второй аргумент, а  — в третий аргумент функции . Тогда получим функцию одного переменного, которая значению ставит в соответствие число . Задача о брахистохроне относится к числу простейших задач классического вариационного исчисления. В общем виде подобная задача принимает следующий вид. (6.24) где «» означает экстремум (в частности, максимум или минимум) функционала, а множество ограничений: (6.25) Решение поставленной классической задачи вариационного исчисления основано на известной теореме Эйлера. Пусть  — непрерывно дифференцируемая функция трех переменных. Тогда, если функция доставляет локальный экстремум функционалу (6.24), то удовлетворяется следующее уравнение: (6.26) где Если в задаче (6.24) появляются дополнительные ограничения, типа , то для ее решения необходимо составить функцию Лагранжа и затем считать, что ее переменные независимы, т. е. можно применить теорему Ферма. Вновь вернемся к задаче о брахистохроне. Во времена И. Бернулли уже было известно, что уравнение: представляет собой уравнение циклоиды. Осуществим подстановку: После преобразования, получим окончательно: Заменяя на и на получим уравнение циклоиды. Из полученных формул вытекает способ построения искомой кривой И. Бернулли (брахистохроны). График решения задачи И. Бернулли показан на рис. 6.6. Следует заметить, что все циклоиды гомотетичны и выпуклы, следовательно, можно взять любую из них, имеющую точку своей левой вершиной и соединить эту вершину с точкой прямой линией, которая пересечет построенную циклоиду в точке . Затем следует сделать гомотетию с коэффициентом , таким образом, будет получена искомая брахистохрона. Построение брахистохроны от точки (0; 0) до точки (-8; 12,5) приведено в протоколе 6.6. К сожалению, авторам монографии не известны компьютерные апробированные программы решения широкого класса задач вариационного исчисления. Возможно, они существуют, но степень их использования мало известна. Рис. 6.6. Пояснение к решению задачи И. Бернулли Протокол 6.6. Вид брахистохроны, построенной в среде Mathcad К классу задач вариационного исчисления можно отнести так называемые изопериметрические задачи, т. е. задачи выбора функции, доставляющей максимум или минимум площади (объему, стоимости) при заданном периметре (площади поверхности, уравнении состояния объекта или процесса). Например, задачу Дидоны. Рассмотрим популярную при изучении динамики средних и проблем экологической безопасности задачу, теоретическая часть которой подробно изложена Е. С. Вентцель [7]. Здесь эта задача будет доведена до решения в среде Mathcad. Требуется ее решение в Matlab. Пусть в некоторой местности обитают животные двух видов: (А) — хищники, (В) — травоядные (жертвы). Животные (А) питаются животными вида (В), которые довольствуются растительной пищей. Состояние каждого животного будем характеризовать только по двум признакам: животное живо или погибло. Граф состояний системы (рис. 6.7) разделим на два графа, соответствующие каждому виду в отдельности. Здесь стрелки, ведущие из состояния в и из в , учитывают смертность животных, причем для травоядных (В) смертность может наступать по двум причинам: естественная смерть и гибель от хищника. Двойные стрелки, входящие в соответствующие состояния символизируют пополнение численности животных за счет рождаемости. Рис. 6.7. Размеченные графы изменения численности животных Пусть ,  — численности хищников в состоянии (живы) и в состоянии (мертвы) соответственно. Аналогично и  — соответствующие состояния для травоядных. Математические ожидания в данный момент времени для численности хищников будем обозначать через и , а для травоядных — , . Введем обозначения: ,  — интенсивности гибели хищников и травоядных; ,  — интенсивности прироста (рождаемости) хищников и травоядных. Для составления дифференциальных уравнений изменения численности животных необходимо задаться видом зависимостей , , , . Рассмотрим динамику численности травоядных. Предположим, что запасы пищи для травоядных не зависят от их численности и от численности хищников. Тогда можно предположить, что средняя рождаемость в единицу времени животных вида В (в пересчете на одно животное) остается постоянной. Обозначим эту постоянную через и получим прирост численности травоядных (состояние ) в виде: (6.27) Естественно предположить, что число встреч животных А и В, заканчивающихся тем, что животное А поедает животное В, пропорционально числу () живых животных вида А и числу () имеющихся животных вида В, т. е. число таких встреч (число животных вида В, которых съели животные вида А) равно: где  — константа (эффективность охоты хищника). Обозначим через постоянную смертности травоядных вследствие естественной смерти, тогда число погибших травоядных в единицу времени будет состоять из двух слагаемых: (6.28) при этом второе слагаемое в составе  — это число животных, поедаемых в единицу времени хищниками. Разделим (6.28) на и получим относительную численность гибели травоядных: В пересчете на одно животное вида В, т. е. при  = 1 интенсивность гибели травоядных будет равна: (6.29) Рассмотрим динамику численности животных вида А. Поскольку согласно принятой модели единственной пищей хищников являются травоядные, то и рождаемость и смертность животных вида А будут зависеть от числа поедаемых в единицу времени животных вида В. Число поедаемых травоядных определено числом эффективных встреч хищников с травоядными, а в пересчете на одного животного вида А составит величину (напомним, что — число живых травоядных). Таким образом, интенсивность смертности хищников будет представлять собой некоторую функцию аргумента , но поскольку  — константа, а изменение пока остается неизвестным, можно записать: (6.30) Очевидно, что функция (6.30) будет убывающей функцией своего аргумента, т. е. чем больше у хищников будет пищи, тем меньше их будет умирать. Аналогично, рождаемость хищников будет также зависеть от наличия травоядных. Обозначим через функцию рождаемости хищников приходящейся на одно животное в единицу времени. Легко заметить, что  — должна быть возрастающей функцией своего аргумента. Рисунок 6.8 отражает изменение функций смертности и рождаемости хищников в зависимости от значений аргумента. Средний прирост популяции хищников за счет рождаемости (в единицу времени) определится выражением: (6.31) Теперь представляется возможным составить дифференциальные уравнения для средних численностей состояний: (6.32) где  — изменение во времени математического ожидания численности живых хищников;  — изменение во времени математического ожидания численности живых травоядных. Пользуясь принципом квазирегулярности и заменяя численности , , от которых зависят интенсивности (6.27), (6.29), (6.30), (6.31), их средними значениями , , получим: (6.33) где  — постоянная смертности травоядных вследствие естественной смерти;  — эффективность охоты хищников;  — число трагичных для травоядных встреч;  — средняя рождаемость травоядных в единицу времени в пересчете на одно животное;  — интенсивность смертности хищников;  — средний прирост хищников за счет рождаемости. Рис. 6.8. Изменение функций смертности и рождаемости хищников в зависимости от значений аргумента В первое уравнение системы (6.33) входит только разность двух функций, которую можно обозначить через : (6.34) где  — текущее значение численности травоядных (). Поскольку функция (6.34) представляет собой разность возрастающей функции и убывающей функций , то сама функция  — возрастающая и в отличие от положительных функций , функция может менять знак (рис. 6.9). Содержательное значение функции (6.34) сводится к тому, что она определяет относительный прирост численности хищников в единицу времени в зависимости от числа травоядных. Предположим, что функция (6.34) меняет знак в точке , которая имеет смысл той численности травоядных, при которой рождаемость и смертность хищников в среднем уравновешиваются. При увеличении численности травоядных сверх численность хищников будет в среднем возрастать, при уменьшении численности животных вида В — численность хищников в среднем будет убывать. Назовем критической численностью травоядных (животных вида В). По аналогии можно ввести понятие критической численности хищников, при которой число травоядных животных в среднем остается неизменным. Рис. 6.9. Пример изменения функции Из второго уравнения системы (6.33) при и , получим: Откуда, учитывая, что , имеем после деления на : (6.35) Если число хищников начнет превышать , то численность травоядных будет убывать, а при числе хищников меньшем  — численность травоядных будет возрастать. Пользуясь введенными обозначениями, уравнение динамики средних (6.33) можно переписать в следующем виде: (6.36) Анализ средних численностей популяций Решение системы (6.36) в заданный момент времени — это точка на плоскости , . Введем обозначение: и рассмотрим движение во времени точки , изображающей решение системы уравнений (6.36) на плоскости в координатах , , обращаясь к рис. 6.10, где нанесены критические значения численности хищников и травоядных . Рис. 6.10. Траектория движения точки решения системы (6.36) на плоскости Выберем в качестве начальной точку A на оси . Если при этом численность хищников убывает , что возможно только в случае , т. е. если численность травоядных меньше критической, то с течением времени движение точки будет направлено вниз (численность травоядных уменьшается). Если при этом и (численность травоядных убывает), что возможно в случае , то движение точки будет происходить влево. Т. е. общее направление движения точки  — влево и вниз. При (точка пересекла прямую , прошла точку В) производная согласно второму уравнению системы (6.36) меняет знак и становится положительной, т. е. точка , продолжая движение влево (производная остается отрицательной), начинает одновременное движение вверх к точке С. После прохождения точкой прямой (прохождение точки С) функция меняет знак и становится положительной, т. е. точка  начинает движение вправо, поскольку, как следует из первого уравнения системы (6.36), производная становится положительной. Одновременно точка продолжает движение вверх к точке D вследствие положительного значения производной . После того, как траектория движения вновь пересечет прямую (проходит точку D), левая часть второго уравнения в (6.36) становится отрицательной и движение точки происходит вниз к точке Е. Следовательно, точка , изображающая решение системы дифференциальных уравнений (6.36), вращается вокруг точки с координатами (,) по часовой стрелке, причем это вращение в общем случае происходит не по окружности, приближаясь к этой точке или удаляясь от нее. Точка (,) представляет положение равновесия системы. Если в качестве начальных условий для решения системы (6.36) выбрать (,), то численности , не будут изменяться и решение системы (5.36) не будет зависеть от времени. Действительно, при функция и правая часть первого уравнения в системе (6.36) принимает нулевое значение. Аналогично, при обращается в ноль правая часть второго уравнения: Следовательно, независимо от времени, имеем: (6.37) Если при иных начальных условиях с течением времени точка неограниченно будет приближаться к этой точке, то рассматриваемая система имеет устойчивое статическое решение. Решение системы не будет устойчивым, если с течением времени точка  будет удаляться от точки (,). В общем случае возможны три варианта: a. после одного «оборота» точка , соответствующая решению системы (6.36), окажется ближе к точке (,); b. после одного «оборота» точка решения системы окажется дальше от точки (,); c. после одного «оборота» точка решения системы в точности вернется к своему начальному положению. В последнем случае все последующие обороты точки будут совершаться по одной и той же замкнутой траектории. Такой путь в теории дифференциальных уравнений называется циклом. Рассмотрим варианты дальнейшего развития процесса, изображенного на рис. 6.10. Следующий оборот точки вокруг положения равновесия еще больше приблизит ее к точке (,), и с каждым следующим оборотом точка будет все больше приближаться к положению равновесия, т. е. к точке (,). Здесь возможны два варианта: – точка , соответствующая решению системы (6.36), будет неограниченно приближаться к положению равновесия (к точке (,)), в этом случае положение равновесия будет устойчивым; – точка будет приближаться к точке (,), но не неограниченно, т. е. точка будет неограниченно приближаться снаружи к какому-то циклу, который, естественно, будет устойчивым. Аналогичные рассуждения можно провести для случая, когда после первого оборота точка удалится от положения равновесия. В этом случае: – точка может неограниченно удаляться от положения равновесия (,); – точка может неограниченно приближаться изнутри к какому-то устойчивому циклу. Рассмотрим пример решения системы (6.36), изображенный траекторией на рис. 6.11. В том случае, когда строго реализуется вариант «с» замкнутому циклу соответствует периодическое решение , т. е. если число хищников возрастает, то число травоядных начинает сокращаться, достигая критической численности . После этого начинает убывать число хищников. После того, как число хищников минует критическое значение , численность травоядных начнет возрастать. В конце концов, численности хищников и травоядных достигнут своих первоначальных значений, и процесс возобновится. По аналогии с положением равновесия, цикл также может быть устойчивым, если решения, начинающиеся вблизи его внешней границы, неограниченно приближаются к циклу, и неустойчивым, если они удаляются от цикла. Такие решения будут соответствовать устойчивому и неустойчивому динамическому равновесию системы. Рис. 6.11. Неординарный пример траектории точки решения системы уравнений (6.36) с одним устойчивым циклом «С» Рассмотрим еще один вариант возможного предельного поведения решений системы (6.36). Положение равновесия может быть неустойчиво, хотя и имеется один устойчивый цикл (рис. 6.11). В этом случае предельный режим выражается в периодическом возрастании и убывании численностей, и независимо от начальных условий точка  будет удаляться от точки (,), соответствующей устойчивому решению, приближаясь к внутренней границе цикла. Рассмотренные варианты решений системы дифференциальных уравнений (6.36) могут составлять различные комбинации, образуя иногда довольно сложную картину. В частности, положение равновесия может быть неустойчиво, и циклов нет. При таком варианте никакого предельного режима не существует. Численность животных каждого вида то убывает почти до нуля, то быстро возрастает, причем с каждым «оборотом» эти колебания становятся все больше, возрастая неограниченно. Естественно, что практическая реализация такого случая не встречается, поскольку начнут сказываться не учтенные в постановке задачи и в рассматриваемой модели факторы: ограниченность запасов пищи или пространства, ошибки принципа квазирегулярности (возрастающие при малых значениях численностей), миграция животных, способность разнообразного питания и т. д. Возможны и более сложные случаи, один из примеров которых изображен на рис. 6.12. Здесь имеем устойчивое положение равновесия, вокруг которого существует неустойчивый цикл . За внешней границей цикла существует еще один, но уже устойчивый цикл . Рис. 6.12. Пример стремления решений системы (6.36) к различным циклам В этом варианте два предельных режима: – положение равновесия (,); – периодический режим, соответствующий циклу . К первому предельному режиму численности животных , будут приближаться, если начальные условия , находились внутри цикла (заштрихованная область). Ко второму режиму численности , будут стремиться, если исходная точка находилась вне цикла . Для нахождения предельных режимов системы (6.36) необходимо исследование этой системы дифференциальных уравнений при конкретных значениях входящих в нее числовых параметров и конкретном виде функции . Пример 6.3 Пусть функция , определяющая относительный прирост численности хищников в единицу времени в зависимости от числа травоядных, имеет вид: (6.38) Откуда значение критической численности травоядных, полученное при будет равно: . Пусть , тогда при , , получим критическое значение численности хищников . В этом случае система уравнений динамики средних (6.36) будет иметь вид: (6.39) Численное решение системы (6.39) при различных начальных условиях показывает, что существует неустойчивое положение равновесия (, ) и один устойчивый цикл. Результаты решения представлены графиками рис. 6.13, рис. 6.14. Рис. 6.13. Изменение во времени численности особей хищников и травоядных , полученное в результате решения системы (6.39) Для анализа полученных результатов рассмотрим графики рис. 6.13 и выберем в качестве начальной точку условных единиц времени. В этот момент численность хищников находится в минимуме и составляет примерно 563 единицы, а численность травоядных равна своему критическому значению . С течением времени численность травоядных возрастает, так как численность хищников все еще меньше своего критического значения ( животных). В момент времени примерно равный 14 единицам численность хищников достигает критического значения, а численность травоядных достигает максимума, (около 2089 голов). Далее с течением времени численность хищников продолжает возрастать, так как численность травоядных превышает критическое значение и пищи для хищников вполне достаточно. В тот момент (примерно 20 условных единиц времени), когда численность травоядных снижается до критического значения, численность хищников, достигнув максимума (около 1500 единиц), перестает возрастать и далее начинает убывать, поскольку число травоядных становится меньше критического значения и продолжает убывать. В момент времени, равный примерно 27 условным единицам численность хищников становится равной критическому значению и продолжает убывать вследствие нехватки пищи (численность травоядных к этому моменту достигает своего минимума и начинает возрастать). В момент времени, соответствующий 36 единицам, численность травоядных достигает критического значения, и весь цикл повторяется заново. Таким образом, через 36 единиц времени весь цикл изменения численности травоядных и хищников повторится. Вид устойчивого цикла решения системы уравнений (6.39) представлен на рис. 6.14. Рис. 6.14. Устойчивый цикл развития популяций хищников и травоядных Решение системы (6.39) производилось в системе Mathcad при начальных данных: численность хищников составляла 800 единиц при численности травоядных — 500 единиц. Устойчивый цикл представляет собой аналог решения проблемы экологической безопасности.