Журнал вычислительной математики и математической физики, 2019, T. 59, № 10, стр. 1666-1680

ВВЕДЕНИЕ

Подходы к решению задач многокритериальной оптимизации (МКО) [1], основанные на аппроксимации и изучении границы Парето (см. [2]–[5]), занимают важное место среди других подходов к решению этой задачи (см., например, [5], [6]). При использовании этого подхода часто возникают затруднения, вызванные тем, что задача аппроксимации границы Парето поставлена некорректно, поскольку граница Парето зачастую неустойчива по отношению к возмущениям параметров задачи (см. [7]). В связи с этим был предложен альтернативный подход к анализу границы Парето, подробно описанный в [8], [9]. Этот подход основан на аппроксимации оболочки Эджворта–Парето (ОЭП) множества достижимых критериальных векторов (которую иногда, особенно в экономических проблемах, называют Free Disposal Hull), и которая, как правило, устойчива по отношению к возмущениям параметров задачи. Аппроксимация ОЭП используется в [8, 9] для визуализации многомерной границы Парето с помощью двумерных сечений ОЭП.

Рассмотрим общую формулировку задачи МКО [1]. Пусть допустимыми решениями задачи МКО являются элементы x компактного множества $X \subset {{{\text{R}}}^{n}}$. Пусть значения компонентов вектора критериев оптимизации $y \in {{{\text{R}}}^{m}}$ задаются вектор-функцией y = f(x), причем в ${\text{R}}_{{}}^{n}$ и ${\text{R}}_{{}}^{m}$ введены метрики, которые будут уточнены в дальнейшем. Для определенности рассмотрим задачу многокритериальной минимизации, т.е. будем считать, что желательно уменьшение значения каждого из отдельных критериев при неизменных остальных. Точнее говоря, будем считать, что критериальная точка y' более предпочтительна, чем точка y (другими словами, y' доминирует y по Парето), если y' ≤ y и y' ≠ y. Решением задачи многокритериальной минимизации является паретова (недоминируемая) граница P(Y) множества достижимых критериальных векторов Y = f(X), определяемая как

и совокупность точек $P(X) \subset X$, порождающих P(Y). Совокупность $P(X)$ называется множеством решений, оптимальных (эффективных) по Парето.

Оболочка Эджворта–Парето (ОЭП) множества Y = f(X) в задачах многокритериальной минимизации определяется как Y* = Y + ${\text{R}}_{ + }^{m}$, где ${\text{R}}_{ + }^{m}$ – неотрицательный конус пространства ${\text{R}}_{{}}^{m}$. ОЭП может быть представлена в виде

Как видно, кроме точек, принадлежащих Y, множество Y* содержит недостижимые точки, доминируемые точками Y.

Аппроксимация ОЭП используется в визуальном методе поиска предпочтительных решений в задачах многокритериальной оптимизации, имеющим название метод достижимых целей/диалоговых карт решений (МДЦ/ДКР). Благодаря тому, что P(Y*) = P(Y) (см. [1], [6]), визуализация Y* позволяет получить информацию о P(Y). Важно, что доминируемые границы множества Y, затрудняющие визуализацию при m > 2, исчезают при переходе к множеству Y*, что в значительной степени облегчает визуализацию P(Y). Основными шагами МДЦ/ДКР являются (см. [8], [9]):

1) аппроксимация ОЭП множества достижимых критериальных векторов;

2) диалоговый анализ границы Парето множества достижимых критериальных векторов экспертом или лицом, принимающим решение (ЛПР), на основе визуализации карт решений – наборов двухмерных сечений ОЭП;

3) выбор наиболее предпочтительной достижимой цели на границе Парето (также осуществляется экспертом или ЛПР);

4) информирование ЛПР о решении, приводящем к выбранной достижимой цели.

В рамках МДЦ/ДКР наиболее сложной вычислительной проблемой является аппроксимация ОЭП. В выпуклом случае (в том числе линейном) разработаны эффективные методы полиэдральной аппроксимации ОЭП (см. [8], [9]), которые применимы для решения задач многокритериальной оптимизации при числе критериев от двух до десяти и большой размерности пространства решений ${\text{R}}_{{}}^{n}$. В случае нелинейных невыпуклых задач многокритериальной оптимизации, которым посвящена настоящая работа, ситуация является более сложной.

Заметим, что внутреннюю аппроксимацию ОЭП в невыпуклом случае разумно строить в виде объединения конусов

задаваемую конечным числом точек $T \subset Y$ (так называемой базой аппроксимации). Обратим внимание на то, что T* – это ОЭП множества T. В идеале точки базы T должны принадлежать P(Y), поэтому в случае использования аппроксимации (0.2) задачи аппроксимации границы Парето и ОЭП близки между собой.

Как известно, если задача многокритериальной оптимизации описана с использованием заданных в явном виде математических функций, для которых можно найти постоянную Липшица, то при достаточно малых величинах постоянной Липшица и малых размерностях пространств решений и критериев для аппроксимации ОЭП можно использовать методы покрытия (см. [3]). Если же размерность пространства решений велика, но число критериев по-прежнему не слишком велико, причем функции (свертки) критериев имеют не слишком большое число локальных экстремумов, то для аппроксимации ОЭП можно использовать глобальную оптимизацию сверток [2]. В более сложных случаях, когда связь между переменными и критериями задается вычислительным модулем (“черным ящиком”), а число локальных экстремумов сверток достигает нескольких десятков, можно использовать так называемые многофазные методы, предложенные в работах [10], [11] и подробно изученные в [12]–[14]. Эти методы представляют собой различные способы распространения на задачи МКО идей мультистарта – известного подхода к глобальной оптимизации в невыпуклых задачах, который основан на градиентной локальной оптимизации, начинающейся из случайных точек множества решений [15].

Как показывает опыт практического использования, методы мультистарта не подходят для аппроксимации ОЭП в случае задач МКО, в которых число локальных экстремумов сверток критериев велико. В таких задачах, в принципе, можно использовать генетические методы аппроксимации границы Парето [16], [17], в рамках которых наличие большого числа локальных экстремумов сверток критериев не вызывает затруднений. Недостаток генетических методов состоит в том, что они могут сходиться очень медленно.

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

Надо отметить, что идея использовать решения задач глобальной скалярной оптимизации отдельных критериев в популяциях генетических методов аппроксимации границы Парето не нова – она была предложена в еще 1985 г. в первом многокритериальном генетическом методе VEGA (vector evaluated genetic algorithm) (см. [18]), где задачи глобальной скалярной оптимизации решались также с помощью генетических методов. Экспериментальные исследования, однако, показали, что такой подход “обычно не находит хорошие компромиссы между значениями критериев” [19, гл. 23]. Судя по всему, нахождению хороших компромиссов в этих экспериментах помешали недостаточно точное решение задач глобальной оптимизации и использование популяций малых объемов. В данной статье мы показываем, что достаточно точное решение задач глобальной оптимизации отдельных критериев с помощью классических методов мультистарта и использование полученных оптимальных решений в популяциях больших объемов позволяют принципиально ускорить сходимость генетического метода. Таким образом, нами предлагается новый метод, который позволяет достаточно точно аппроксимировать ОЭП, по крайней мере, для одного класса сложных нелинейных прикладных задач МКО, а именно задач построения правил управления каскадом водохранилищ с критериями типа надежности выполнения предъявляемых требований к каскаду. Такая задача характеризуется заданием модели в виде вычислительного модуля, наличием большого числа кусочно-постоянных критериев выбора решения и очень большого числа локальных экстремумов сверток критериев.

В первом разделе статьи описывается математическая модель каскада водохранилищ, к правилам управления которыми предъявляются противоречивые требования, что делает задачу выбора правил задачей МКО. Второй раздел посвящен экспериментальному анализу сходимости известного генетического алгоритма NSGA-II (см. [16], [19]) при аппроксимации ОЭП для описанной задачи. В третьем разделе дается описание метода инжекции оптимумов. В четвертом разделе статьи приводятся результаты сравнения скорости сходимости и точности аппроксимации алгоритма NSGA-II и метода инжекции оптимумов на основе экспериментов по аппроксимации ОЭП для задачи управления каскадом водохранилищ.

1. ЗАДАЧА ВЫБОРА ПРАВИЛ УПРАВЛЕНИЯ КАСКАДОМ ВОДОХРАНИЛИЩ

В данной статье предлагаемый метод аппроксимации ОЭП изучается на основе прикладной задачи управления каскадом из $I$ водохранилищ, расположенных в основном русле реки. Водохранилища используются для регулирования уровня реки, производства электроэнергии, предотвращения наводнений, обеспечения речного судоходства и водоснабжения в населенных пунктах и т.д. Пусть i – номер водохранилища, $i = 1,\;...,\;I$, причем водохранилище с номером $i + 1$ находится ниже i-го водохранилища. Математическая модель каскада представляет собой динамическую систему с дискретным временем (подробнее см. в [20], [21]). Считается, что интервал времени t начинается в момент t–1 и заканчивается в момент t. Модель каскада рассматривается в течение периода времени продолжительностью в T лет, причем каждый год в модели разбивается на ${{n}_{0}}$ календарных интервалов, не обязательно равных по продолжительности. Тогда общее число интервалов за T лет составляет ${{t}_{0}} = {{n}_{0}}T$.

Состояние динамической системы в момент t задается величинами $w_{i}^{t}$, $i = 1,\;...,\;I$, где $w_{i}^{t}$ – объем воды в i-м водохранилище. Динамика величин $w_{i}^{t}$ описывается балансовым уравнением

где $r_{i}^{t}$ – попуск из i-го водохранилища за интервал времени t; $q_{i}^{t}$ – боковая приточность к i-му водохранилищу за интервал времени t. Термин “боковая” используется для того, чтобы отличать приточность $q_{i}^{t}$ от приточности $r_{{i - 1}}^{t}$, обеспечиваемой попуском из водохранилища, лежащего выше по течению.

Кроме уравнения состояния, в модели задаются некоторые нелинейные функции (см. [20]), на основе которых по состоянию i-го водохранилища в начале и конце интервала времени t и попуска за этот интервал времени можно рассчитать средние за интервал уровни воды выше и ниже плотины. В свою очередь, на основе уровней можно рассчитать среднюю выработку электроэнергии в единицу времени на i-м гидроузле в этот интервал времени и другие важные показатели функционирования водохранилища.

Управляющими воздействиями в задаче управления каскадом являются попуски воды через плотины $r_{i}^{t}$, $i = 1,\;...,\;I$, $t = 1,\;...,\;{{t}_{0}}$. Поскольку предсказать приточность водохранилищ на весь рассматриваемый период в T лет заранее невозможно, специалисты управляют попуском на основе правил, в которых величина попуска определяется в начале каждого интервала времени как функция имеющейся информации. В используемом нами правиле попуска, которое обсуждается в [22], попуск i-го водохранилища на интервале $t$ определяется в начале интервала, т.е. в момент $t - 1$, согласно соотношению

где ${{\psi }_{i}}( \cdot , \cdot , \cdot , \cdot )$ – некоторый заданный алгоритм расчета попуска i-го водохранилища, а ${{\alpha }_{i}}$ – набор его параметров, значения которых (для всех водохранилищ) определяются в процессе построения правила попуска. Величина $w_{i}^{{t - 1}}$ в момент $t - 1$ известна, а величина боковой приточности $q_{i}^{t}$ прогнозируется в момент $t - 1$ достаточно точно, так что ее можно считать известной. Отметим, что использовать правило (1.2) можно только тогда, когда уже определен попуск $r_{{i - 1}}^{t}$ водохранилища i – 1, лежащего выше по течению.

Предполагается, что совокупность параметров $\alpha = \left( {{{\alpha }_{i}},\;i = 1,...,I} \right)$, содержащихся в алгоритмах расчета попуска всех водохранилищ, должна принадлежать некоторому выпуклому ограниченному множеству $\Xi $, которое определяется особенностями правил управления водохранилищами и здесь уточняться не будет. Подчеркнем, что параметры правил управления для всех водохранилищ должны быть согласованы и определены в рамках одной задачи принятия решений (так называемые синхронизированные правила). Это увеличивает размерность пространства решений и число критериев, а также предъявляет дополнительные требования к типу правил управления (1.2).

Как уже говорилось, управление водохранилищами является многоцелевым – оно должно удовлетворять экономические, экологические, коммунальные и прочие требования к различным показателям водопользования. В том случае, когда требования всех водопользователей не могут быть удовлетворены одновременно, приходится рассматривать задачу многокритериального выбора правил управления, в рамках которой должно быть найдено компромиссное решение о значениях критериев выбора решения, характеризующих степень нарушения различных требований. В качестве критериев берется доля перебоев, т.е. интервалов времени, на которых происходило нарушение требований. Выбор именно таких критериев обусловлен указаниями Минприроды РФ. Это предопределяет сложность решаемой задачи МКО, состоящей в следующем.

Для математической формулировки критериев удобно использовать функцию Хэвисайда $\Theta (z)$, равную нулю при ${\text{ }}z \leqslant 0$ и единице при ${\text{ }}z > 0$. Величина ${\text{ }}z$ – это значение некоторого показателя функционирования водохранилища, причем при ${\text{ }}z \leqslant 0$ требование к этому показателю выполнено, а при ${\text{ }}z > 0$ – имеет место перебой. Рассмотрим примеры.

В качестве примера критерия, связанного с односторонним требованием, приведем критерий, отражающий долю перебоев в производстве электроэнергии на i-й гидроэлектростанции за весь период времени

где $z_{e}^{{i,t}} = e_{i}^{*} - e_{i}^{t}$, где $e_{i}^{*}$ – требование к выработке электроэнергии на интервале t, $e_{i}^{t}$ – фактическая выработка, $t = 1,\;...,\;{{t}_{0}}$. Таким образом, при $z_{e}^{{i,t}} > 0$ величина $z_{e}^{{i,t}}$ является нарушением – отклонением от выполнения требования к выработке электроэнергии на интервале t. Значение критерия (1.3) – это доля интервалов времени, в течение которых имело место нарушение, т.е. недопроизводство электроэнергии.

Двустороннее требование $h_{i}^{{\min }} \leqslant h_{i}^{{}} \leqslant h_{i}^{{\max }}$ к уровню воды $h_{i}^{t}$ выше плотины i-го водохранилища учитывается в виде критерия:

где $z_{h}^{{i,t}} = \max (h_{i}^{t} - h_{i}^{{\max }},h_{i}^{{\min }} - h_{i}^{t})$, а $[h_{i}^{{\min }},h_{i}^{{\max }}]$ – желаемый диапазон значений уровня величин. При $z_{h}^{{i,t}} > 0$ эта величина является отклонением уровня $h_{i}^{t}$ от желаемого диапазона. Аналогичное двустороннее ограничение накладывается на попуски $r_{i}^{t}$, диапазоны значений которых выбираются на основе требований нормальной работы гидроузла, отсутствия наводнений и др. Ясно, что критерии типа (1.3), (1.4) принимают значения от 0 до 1, причем их значения желательно уменьшать.

Для того чтобы сравнить последствия различных наборов параметров ${{\alpha }_{i}},i = 1,\;...,\;I$, необходимо уметь по заданным величинам параметров находить значения критериев типа (1.3), (1.4). Описанная модель пока не представляет такой возможности, поскольку она является открытой – в ней не заданы внешние воздействия на систему, а именно боковая приточность $q_{i}^{t}$ $i = 1,\;...,\;I$, $t = 1,\;...,\;{{t}_{0}}$, в балансовом уравнении (1.1). Для преодоления этой проблемы в водохозяйственной практике принято моделировать будущую приточность на основе приточности, наблюдавшейся в прошлом. В частности, часто рассматривают достаточно продолжительный промежуток времени (от нескольких десятков до сотни лет), для которого собраны исторические значения боковой приточности. Эта информация определяет число интервалов времени ${{t}_{0}}$ в соотношениях (1.1)–(1.4), для которых рассчитываются значения критериев. По начальным значениям объемов воды в водохранилищах, историческим величинам приточности $q_{i}^{t}$, $i = 1,\;...,\;I$, $t = 1,\;...,\;{{t}_{0}}$, и выбранным параметрам ${{\alpha }_{i}},i = 1,\;...,\;I$, на основе соотношений модели можно для всех водохранилищ (начиная с верхних) последовательно рассчитать объемы воды и попуски на весь период времени. По результатам этих расчетов находятся значения тех показателей, для которых сформулированы требования. По реализовавшимся значениям показателей и заданным заранее требованиям рассчитываются соответствующие величины ${\text{ }}z$, что позволяет по формулам типа (1.3), (1.4) найти значения критериев.

Описанная модель использовалась для анализа бассейна реки Ангара. При этом каждый год был разбит на 22 интервала (5 месяцев разбито на три декады каждый, а для 7 месяцев интервалом является месяц целиком). Был подготовлен исторический ряд приточности продолжительностью 103 года, содержащий 2266 отрезков времени. Число параметров правила управления для одного интервала времени было равно шести, так что правило попуска воды для Иркутского водохранилища содержало $6 \times 22 = 132$ параметра, а весь каскад, в который были включены три работающих водохранилища, описывался 284 параметрами правил управления. При выборе параметров было принято во внимание 24 критерия, значения которых можно было рассчитать по заданным наборам параметров.

Для этой модели была поставлена задача аппроксимации ОЭП с точностью в 1% в максимальной метрике (т.е. отклонение по каждому из критериев не должно превышать 0.01, поскольку критерии принимают значения от 0 до 1). Было заранее предъявлено требование к числу расчетов критериальной вектор-функции (далее будем говорить просто о критериальной функции). Дело в том, что в будущих прикладных исследованиях задачу аппроксимации ОЭП потребуется решать многократно с различными вариантами исходной гидрологической информации, требований и ограничений. Кроме того, вместо рассматриваемого в данной работе исторического ряда длиной в 103 года, планируется использовать искусственные ряды приточности для периода в 10 тысяч лет. В связи с этим число расчетов критериальной функции не должно превосходить 10–15 миллионов, что соответствует приблизительно 9–12 часам работы ноутбука среднего класса для периода в 100 лет и такой же продолжительности работы многопроцессорного кластера со 100 процессорами для периода в 10 тысяч лет.

Была сделана попытка решить эту задачу с помощью двухфазного метода, основанного на решении задач локальной оптимизации с помощью градиентных методов. В связи с тем, что каждый из критериев (1.3), (1.4) является суммой большого числа (двух с лишним тысяч) функций Хэвисайда от переменных $z_{{}}^{{i,t}} = 0,\;i = 1,\;...,\;I,\quad t = 1,\;...,\;{{t}_{0}}$, описывающих отклонения от требований во все моменты времени, критерии как функции параметров правил управления являются кусочно-постоянными, причем имеют огромное число поверхностей разрыва. Для того чтобы использовать градиентные методы поиска локальных минимумов критериев, был предложен метод эрзац-функций [23], в рамках которого градиентные методы применяются к вспомогательным функциям, получаемым из исходного критерия типа (1.3), (1.4) заменой функций Хэвисайда на некоторые непрерывные функции (эрзац-функции). В методе эрзац-функций в точках локальных оптимумов вспомогательной функции рассчитываются значения критерия, после чего в качестве приближенного решения задачи глобальной оптимизации критерия выбирается та точка, в которой достигается его наименьшее значение. Исследования, проведенные в [23], показали, что использование степенных эрзац-функций степени от 0.5 до единицы позволяет найти достаточно точное решение задачи минимизации критериев ${{y}_{j}},\;j = 1,\;...,\;m$.

В методах аппроксимации границы Парето или ОЭП (см., например, [2], [10]–[12]) часто используются функции (свертки) критериев. Для свертки критериев можно провести аналогичную замену функций Хэвисайда на степенные функции. При этом, однако, оказалось, что из-за наличия в задаче (1.1)–(1.4) огромного числа локальных экстремумов соответствующей вспомогательной функции найти решение задачи, относительно близкое к решению задачи глобальной оптимизации свертки, обычно не удается. Это стало ясно после проведения большой серии экспериментов по аппроксимации ОЭП в задаче МКО для описанной выше модели каскада водохранилищ реки Ангара с 24 критериями типа (1.3), (1.4) на основе итеративного двухфазного метода [11] с различными свертками критериев.

Все эксперименты дали качественно одинаковый результат. Опишем результаты одного из них. В этом эксперименте для задачи (1.1)–(1.4) была получена аппроксимация ОЭП, для чего понадобилось значительное число расчетов критериальной функции (около 110 миллионов расчетов). Для анализа построенной аппроксимации ОЭП была использована контрольная точка, в качестве которой была взята критериальная точка $y({{\alpha }_{{11}}})$, заведомо принадлежащая ОЭП и, более того, решение ${{\alpha }_{{11}}}$ было выбрано экспертом в качестве предпочтительного решения задачи (1.1)–(1.4) в диалоговой процедуре поиска единственного решения задачи МКО (см. [21]). Уточним, что процедура [21] не использовала аппроксимации ОЭП. После 35 миллионов расчетов критериальной функции отклонение δ точки $y({{\alpha }_{{11}}})$ от аппроксимации ОЭП оказалось равно 0.128 и оставалось неизменным на следующих 75 миллионах итераций двухфазного метода. Поскольку эта величина δ на порядок превосходит желаемую величину в 0.01, а перспективы уменьшения отклонения после завершения 110 миллионов расчетов представлялись весьма сомнительными, то был сделан вывод о том, что двухфазный метод не в состоянии достаточно точно аппроксимировать ОЭП в изучаемой задаче.

В связи со сказанным была сделана попытка использовать в рассматриваемой задаче аппроксимации ОЭП генетические методы аппроксимации границы Парето [16], [17], в рамках которых наличие большого числа локальных экстремумов сверток критериев не вызывает затруднений. Многокритериальные генетические алгоритмы построения точек границы Парето [16], [17] основаны на эволюционных концепциях развития популяции допустимых решений, в частности, на их скрещивании (кроссовер) и мутации, а также на отборе наиболее перспективных из них. Главным достоинством генетических алгоритмов являются простота реализации, отсутствие каких-либо требований к модели исследуемого объекта, в том числе непрерывности и, тем более, дифференцируемости критериальных функций, а также возможность применения к задачам типа «черного ящика», в которых критериальные функции заданы вычислительным модулем. Недостатком генетических методов является их возможная медленная сходимость.

Опишем эксперимент, целью которого является поиск ответа на вопрос о том, можно ли использовать генетический алгоритм для решения задачи построения аппроксимации ОЭП в задаче разработки правил управления каскадом. Поскольку медленная сходимость является основной проблемой при использовании генетических методов, изучалась сходимость генетического алгоритма при аппроксимации ОЭП в задаче (1.1)–(1.4).

2. ЭКСПЕРИМЕНТАЛЬНОЕ ИЗУЧЕНИЕ СХОДИМОСТИ ГЕНЕТИЧЕСКОГО АЛГОРИТМА

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

В данной работе для аппроксимации ОЭП в задаче (1.1)–(1.4) используется итерационный алгоритм многокритериальной оптимизации NSGA-II (см. [16], [24]), считающийся одним из наиболее эффективных современных генетических методов аппроксимации границы Парето [21, гл. 23]. При построении новой популяции, кроме обычных генетических операторов (кроссовер и мутация), применяется удаление из популяции точек, близких по значениям критериев к другим точкам. Эти принципы позволяют алгоритму NSGA-II, с одной стороны, удерживать число точек в популяции на заданном уровне, а с другой – обеспечивать их относительную равномерность в пространстве критериев.

Приведем краткое формальное описание алгоритма NSGA-II для задачи построения ОЭП для общей модели (0.1). Заранее должны быть заданы ${{N}_{{{\text{pop}}}}}$ – достаточно большое число точек в выходном множестве генетического алгоритма и ${{N}_{{{\text{it}}}}}$ – число итераций алгоритма.

Краткое описание алгоритма NSGA-II (подробности см. в [16], [24])

В качестве исходной популяции точек ${{Q}_{0}}$ берется случайная выборка некоторого объема из $X$.

3. МЕТОД ИНЖЕКЦИИ ОПТИМУМОВ

Напомним, что нам требуется метод, предназначенный для аппроксимации ОЭП в нелинейных задачах МКО большой размерности (сотни переменных) с существенным числом кусочно-постоянных критериев (порядка 20), заданных с помощью вычислительного модуля и обладающих очень большим числом локальных экстремумов критериев и их сверток. Предлагаемый далее метод инжекции оптимумов (ИО) состоит из двух шагов:

1) решения задачи глобальной оптимизации отдельных критериев с помощью классических градиентных методов оптимизации и

2) построения аппроксимации ОЭП на основе модифицированного генетического алгоритма аппроксимации границы Парето, в популяцию которого периодически включаются оптимальные решения задач глобальной оптимизации отдельных критериев.

Схема метода инжекции оптимумов

Заранее должны быть заданы параметры: ${{N}_{{{\text{pop}}}}}$ – желаемое число точек в популяции (в том числе и в выходном множестве), K – периодичность инжекции оптимумов и ${{\varepsilon }_{{{\text{gen}}}}} > 0$ – параметр правила остановки.

Шаг 1. Достаточно точно решаются задачи глобальной минимизации всех частных критериев. Эти задачи решаются для каждого критерия независимо от других. Обозначим совокупность точек множества допустимых решений X, являющихся решениями задач глобальной минимизации отдельных критериев, через $R$.

Шаг 2. Исходная популяция ${{Q}_{0}}$ модифицированного генетического алгоритма строится следующим образом: совокупность точек $R \subset X$ пополняется случайными точками из множества $X$ до заранее заданного числа ${{N}_{{{\text{pop}}}}}$. Далее выполняются итерации модифицированного генетического алгоритма NSGA-II, дополненные периодической инжекцией оптимумов.

4. ЭКСПЕРИМЕНТАЛЬНОЕ СРАВНЕНИЕ МЕТОДА ИНЖЕКЦИИ ОПТИМУМОВ С ГЕНЕТИЧЕСКИМ АЛГОРИТМОМ

Сравним метод инжекции оптимумов (ИО) с генетическим алгоритмом аппроксимации границы Парето NSGA-II на основе эксперимента по аппроксимации ОЭП для задачи (1.1)–(1.4). Методы сравнивались между собой при одном и том же числе расчетов критериальной функции.

На шаге 1 метода ИО должно быть найдено множество $R$. Обсудим некоторые особенности используемой нами реализации шага 1 метода ИО для модели (1.1)–(1.4). Рассмотрим способ решения этой задачи для некоторого критерия $j \in [1,\;...,\;m]$, имеющего вид (1.3). Считается, что задача уже решена для критериев с номерами $r = 1,\;...,\;j - 1$. Используется метод мультистарта, т.е. на $\Xi $ генерируется случайная равномерно распределенная выборка ${{H}_{N}}$ стартовых точек процедуры локальной оптимизации, которая реализуется программными средствами, разработанными сотрудниками ВЦ РАН им. А.А. Дородницына А.И. Голиковым и Н.И. Грачевым в рамках системы ДИСО [25] и переработанными в 2000-х годах для персональных компьютеров.

Заметим, что значения критериев неотрицательны, поэтому факт нахождения решения со значением критерия, достаточно близким к нулю, означает, что задача глобальной оптимизации этого критерия решена. Этот факт используется при расчете оптимумов $R$ по модели (1.1)–(1.4).

Итак, перед началом шага 1 должны быть заданы N₁ – число элементов выборки ${{H}_{N}}$ и величина ${{\varepsilon }_{0}}$ – допустимое отклонение от нуля критериев оптимизации. Шаг 1 для некоторого критерия $j \in [1,\;...,\;m]$ состоит из следующих пунктов.

1. Проверяем, не выполняется ли условие ${{y}_{j}}(\alpha _{r}^{{}}) < {{\varepsilon }_{0}}$ на решениях $\alpha _{r}^{{}} \in \Xi $, найденных при минимизации уже рассмотренных критериев $r = 1,\;...,\;j - 1$. Если условие выполняется, то решение для критерия j найдено и переходим к следующему критерию.

2. На множестве $\Xi $ генерируется выборка ${{H}_{N}} = \{ \alpha _{j}^{{(1)}},\;...,\;\alpha _{j}^{{({{N}_{1}})}}\} $ объемом N₁.

3. Выполняется N₁ итераций, на каждой из которых решается задача локальной минимизации на множестве $\Xi $ вспомогательной функции

стартующая на q-й итерации из точки $\alpha _{j}^{{(q)}} \in {{H}_{N}}$. Обратим внимание на то, что вместо критерия (1.3) используется непрерывная вспомогательная функция (4.1), полученная путем замены функций Хэвисайда $\Theta (z_{e}^{{j,t}})$ на непрерывные эрзац-функции $z_{e}^{{j,t}}\Theta (z_{e}^{{j,t}})$, равные $z_{e}^{{j,t}}$ при $z_{e}^{{j,t}} > 0$, и равные нулю при $z_{e}^{{j,t}} \leqslant 0$. Для стартовой точки $\alpha _{j}^{{(q)}} \in {{H}_{N}}$ найденную точку локального минимума обозначим через $\hat {\alpha }_{j}^{{(q)}}$. Если выполняется условие то задача минимизации критерия ${{y}_{j}}$ решена, причем локальный оптимум – это текущая точка $\alpha _{j}^{*} = \hat {\alpha }_{j}^{{(q)}}$. Переходим к следующему критерию. Если условие (4.2) не выполняется, переходим к следующей итерации.

4. Если для всех локальных оптимумов функции (4.1), т.е. $\hat {\alpha }_{j}^{{(1)}},\;...,\;\hat {\alpha }_{j}^{{({{N}_{1}})}}$, условие (4.2) не выполняется, то в качестве приближенного решения задачи глобальной минимизации критерия на множестве $\Xi $ берем точку, на которой достигается минимум критерия (1.3) на множестве $\hat {\alpha }_{j}^{{(1)}},\;...,\;\hat {\alpha }_{j}^{{({{N}_{1}})}}$.

Шаг 1 метода ИО завершается после проведения описанных расчетов для всех критериев.

Вычислительные затраты шага 1 метода ИО при N₁ = 100 составили около 10 миллионов расчетов. На шаге 2 метода ИО выполнялись итерации генетического алгоритма с параметрами ${{N}_{{{\text{pop}}}}} = 10\,000$ и остановкой по отклонению с порогом ${{\varepsilon }_{{{\text{gen}}}}}$ = 0.01 без инжекции в популяцию множества $R$ на промежуточных шагах (т.е. величина K на шаге 2 была взята достаточно большой). Процесс аппроксимации на шаге 2 был завершен за 545 итераций, что потребовало около 5.5 миллионов расчетов критериальной функции. На фиг. 3, на которой приведена зависимость отклонения точки $y({{\alpha }_{{11}}})$ от аппроксимации в зависимости от номера итерации, видно, что достигнутое отклонение δ превосходило 0.075, что не соответствует требованиям к точности аппроксимации. Поэтому величина ${{\varepsilon }_{{{\text{gen}}}}}$ была уменьшена в два раза, до 0.005. Требование (3.1) при таком ${{\varepsilon }_{{{\text{gen}}}}}$ было выполнено на шаге 2 только после 3269-й итераций, для чего потребовалось провести около 32.7 миллионов расчетов критериальной функции. При этом оказалось (см. фиг. 3), что δ = 0.005, т.е. точка практически принадлежит аппроксимации. Всего же метод ИО потребовал 42.7 миллионов расчетов критериальной функции. Такого числа расчетов достаточно для 4270 итераций алгоритма NSGA-II. После проведения такого числа итераций алгоритм NSGA-II достиг отклонения δ = 0.05 (см. фиг. 1). Таким образом, по этому показателю точности аппроксимации метод ИО превзошел алгоритм NSGA-II примерно в 10 раз!

Отметим, что попутно в эксперименте методом ИО была достигнута требуемая величина δ = 0.01, для чего потребовалось 2300 итераций шага 2 (см. фиг. 3). Число расчетов критериальной функции, требуемых для получения такой точности методом ИО, близко к 33 миллионам (учитываются затраты на оба шага). Это позволяет осуществить 3300 итераций генетического алгоритма. Как видно на фиг. 1, при таком числе расчетов критериальной функции величина отклонения δ, достигнутая алгоритмом NSGA-II, равна 0.06. Таким образом, при таком числе расчетов критериальной функции метод ИО превосходит генетический алгоритм по рассматриваемому показателю в шесть раз. Как видно, с ростом точности аппроксимации превосходство метода ИО увеличивается. Этот эффект связан с тем, что вычисления шага 1 являются постоянными затратами метода ИО, доля которых в общих затратах падает с ростом затрат шага 2.

Проанализируем асимптотическую скорость сходимости метода ИО. Преобразуем результаты, приведенные на фиг. 3. На фиг. 4, аналогично фиг. 2, по горизонтальной оси откладывается десятичный логарифм числа итераций k, а по вертикальной оси – десятичный логарифм отклонения δ. Так же как и в случае генетического алгоритма, асимптотическое поведение графика экспериментальной зависимости при $\lg k > 3.25$ может быть приближенно описано линейным соотношением $\lg \delta = a + b\lg k$. На основании экспериментальных данных получаем a = 5.2, b = –2.2, поэтому приближенная асимптотическая зависимость отклонения от числа итераций имеет вид

Сравнивая зависимости (4.3) и (2.1), видим, что, в отличие от первой степени числа расчетов, имеющей место в случае генетического алгоритма, имеет место значительно более быстрая, квадратичная скорость сходимости метода ИО. Этот факт подтверждает вывод о предпочтительности метода ИО при достаточно большом числе итераций.

Проведем теперь прямое сравнение двух полученных аппроксимаций ОЭП с использованием метода функций включения, предложенного в [26]. Обозначим базу аппроксимации, построенную методом ИО за 3270 итераций шага 2 метода через ${{T}_{{{\text{OIM}}}}}$. Базу аппроксимации, построенную генетическим алгоритмом за то же число расчетов критериальной функции, обозначим через ${{T}_{{{\text{NSGA}}}}}$. Пусть (A)_ε есть ε-окрестность множества A в максимальной метрике.

Обозначим через ${{\nu }_{{{{T}_{{{\text{NSGA}}}}}}}}(\varepsilon ,T_{{{\text{OIM}}}}^{*})$ долю точек ${{T}_{{{\text{NSGA}}}}}$, попавших в ε-окрестность множества $T_{{{\text{OIM}}}}^{*}$, т.е.

Зависимость ${{\nu }_{{{{T}_{{{\text{NSGA}}}}}}}}(\varepsilon ,T_{{{\text{OIM}}}}^{*})$ от ε при $\varepsilon \geqslant 0$  будем называть функцией включения  базы ${{T}_{{{\text{NSGA}}}}}$ в ε-окрестность множества $T_{{{\text{OIM}}}}^{*}$. Аналогичным образом определяется функция включения базы ${{T}_{{{\text{OIM}}}}}$ в ε-окрестность множества $T_{{{\text{NSGA}}}}^{*}$, которая обозначается ${{\nu }_{{{{T}_{{{\text{OIM}}}}}}}}(\varepsilon ,T_{{{\text{NSGA}}}}^{*})$ и имеет смысл доли точек ${{T}_{{{\text{OIM}}}}}$, попавших в ε-окрестность множества $T_{{{\text{NSGA}}}}^{*}$.

Фиг. 5 содержит графики функций включения ${{\nu }_{{{{T}_{{{\text{NSGA}}}}}}}}(\varepsilon ,T_{{{\text{OIM}}}}^{*})$ (кривая 1) и ${{\nu }_{{{{T}_{{{\text{OIM}}}}}}}}(\varepsilon ,T_{{{\text{NSGA}}}}^{*})$ (кривая 2). На рисунке видно, что ${{\nu }_{{{{T}_{{{\text{NSGA}}}}}}}}(\varepsilon ,T_{{{\text{OIM}}}}^{*})$ $ \geqslant $ ${{\nu }_{{{{T}_{{{\text{OIM}}}}}}}}(\varepsilon ,T_{{{\text{NSGA}}}}^{*})$ при всех значениях $\varepsilon \geqslant 0$. Это означает, что совокупность точек ${{T}_{{{\text{NSGA}}}}}$ лежит ближе к множеству $T_{{{\text{OIM}}}}^{*}$, чем совокупность точек ${{T}_{{{\text{OIM}}}}}$ к множеству $T_{{{\text{NSGA}}}}^{*}$. В частности, множество $T_{{{\text{NSGA}}}}^{*}$ содержит 70% точек ${{T}_{{{\text{NSGA}}}}}$, а ε-окрестность уже при ε = 0.003 содержит более 95% точек ${{T}_{{{\text{NSGA}}}}}$. При этом радиус полного покрытия ${{T}_{{{\text{NSGA}}}}}$ множеством $T_{{{\text{OIM}}}}^{*}$ составляет менее 0.01. Таким образом, точки ${{T}_{{{\text{NSGA}}}}}$ сконцентрированы в окрестности множества $T_{{{\text{OIM}}}}^{*}$.

Наоборот, $T_{{{\text{NSGA}}}}^{*}$ не содержит точек ${{T}_{{{\text{OIM}}}}}$, а при ε = 0.03 в ${{(T_{{{\text{NSGA}}}}^{*})}_{\varepsilon }}$ содержится не более 20% точек ${{T}_{{{\text{OIM}}}}}$, т.е. большинство точек ${{T}_{{{\text{OIM}}}}}$ не близки к $T_{{{\text{NSGA}}}}^{*}$. Учитывая также то, что радиус полного покрытия ${{T}_{{{\text{OIM}}}}}$ множеством $T_{{{\text{NSGA}}}}^{*}$ составляет 0.077, можно утверждать, что аппроксимация ОЭП, задаваемая множеством точек ${{T}_{{{\text{OIM}}}}}$, существенно лучше аппроксимации, задаваемой множеством ${{T}_{{{\text{NSGA}}}}}$.

Таким образом, метод функций включения также подтверждает, что аппроксимация ОЭП, задаваемая базой ${{T}_{{{\text{OIM}}}}}$, значительно превосходит аппроксимацию, задаваемую базой ${{T}_{{{\text{NSGA}}}}}$.

Заметим также, что графический анализ границы Парето, основанный на визуализации двумерных сечений аппроксимации ОЭП, полученной с помощью метода ИО, позволил ЛПР улучшить решение ${{\alpha }_{{11}}}$ по многим важным критериям за счет уступок по менее важным. Этот результат будет описан в отдельной публикации.

ЗАКЛЮЧЕНИЕ. ВЫВОДЫ И ДАЛЬНЕЙШИЕ ИССЛЕДОВАНИЯ

Результаты проведенного эксперимента показывают, что при достаточном числе итераций в рассматриваемой задаче (1.1)–(1.4) метод ИО имеет больший порядок скорости сходимости и дает лучшую аппроксимацию ОЭП, чем алгоритм NSGA-II при том же числе расчетов критериальной функции. Таким образом, решение задач глобальной оптимизации отдельных критериев с помощью методов классического мультистарта и использование их результатов в исходной популяции генетического алгоритма NSGA-II приводит в случае задачи (1.1)–(1.4) к резкому повышению эффективности работы генетического алгоритма. В дальнейшем предполагается изучить вопрос о том, как влияют на эффективность метода его параметры, а именно число стартовых точек в методе мультистарта и объем популяции генетического алгоритма, а также выяснить, для каких еще классов моделей метод ИО является более эффективным, чем алгоритм NSGA-II.