Method of searching for global extremum of a continuous function on a simplex

Full Text

В широкой области технических приложений возникает необходимость решения нелинейной и, в общем случае, невыпуклой задачи: (1) f (x) → max,   m xj = 1; (2) j=1  xj 0, j = 1, 2, . . . , m, где f (x) - некоторая m-мерная непрерывная функция x = (x1 , x2 , . . . , xm ). Для решения невыпуклых оптимизационных задач наиболее часто используют стохастический подход. Он включает широкий спектр методов - от простого случайного поиска до различных рандомизированных эвристик [1-3]. Все эти методы так или иначе используют глобальный и локальный поиск в различных реализациях и сочетаниях. В методе имитации отжига переход от глобального к локальному поиску осуществляется постепенно [4, 5]. На начальном этапе осуществляется случайный поиск по всей допустимой области. Но затем, оставаясь случайным, он все более и более локализуется, становится избирательным. В методе кроссэнтропии [6-8] идея постепенной локализации поиска реализуется через механизм изменения распределения поисковой выборки с тем, чтобы перспективные варианты выпадали с большей вероятностью. Существуют методы, использующие так называемый «мультистарт» локальных поисков (детерминированных или рандомизированных) сразу из множества точек области определения [9, 10]. Для решения некоторых трудно формализуемых задач эффективны генетические алгоритмы [11-13]. В целом же эти алгоритмы обладают всеми достоинствами и недостатками рандомизированных эвристик [14]. Авторы предлагают для решения задачи (1), (2) модификацию известного метода Ψ-преобразования [15-17], сочетающего детерминированный и стохастический подходы. По сравнению с методами [1-13] предложенный метод менее чувствителен к размерности задачи, поскольку не предполагает прямого поиска. Для задачи max{f (x) | x ∈ E}, где E - некоторое подмножество евклидова пространства, в рассмотрение вводится функция Ψ(ζ), которая определяется как мера объема множества ˜ E(ζ) = {x | x ∈ E, f (x) ζ}. Поскольку Ψ(ζ) является монотонно убывающей, ноль этой функции соответствует решению задачи. Определение функции Ψ(ζ) в аналитическом виде затруднительно, поэтому некоторое ее приближение строится по результатам статистических испытаний. О значении Ψ(ζ) можно судить по относительно˜ му количеству случайных точек, попадающих в E(ζ), при условии, что точки равномерно распределены во множестве E. Поэтому в поставленной задаче (1), (2), прежде чем использовать алгоритм Ψ-преобразования, необходимо при статистических испытаниях решить вопрос генерирования случайных точек, равномерно распределенных в области (2). Правильный (m-1)-мерный симплекс (2) можно поместить в соответствующий (m - 1)-мерный единичный гиперкуб в Rm-1 , получение равномерного 756 Об одном методе поиска глобального экстремума . . . распределения случайных чисел в котором - задача тривиальная. Поскольку симплекс является подмножеством гиперкуба, точки, попадающие в симплекс, также равномерно распределены. Однако для поставленной задачи (1), (2) этот путь непродуктивен. Как известно, объем правильного (m - 1)мерного симплекса со стороной единичной длины равен 1 (m - 1)! m 2m-1 , и при росте m быстро стремится к нулю. Соответственно, стремится к нулю и вероятность того, что случайная точка, сгенерированная в гиперкубе, окажется в симплексе. Поэтому для получения случайных допустимых вариантов используем алгоритм [18, 19], который обеспечивает перенос равномерно распределенных случайных точек в симплекс (2) путем аффинных и линейных преобразований, не влияющих на закон распределения. Алгоритм 1. 1. Генерируются (m - 1) случайных чисел, равномерно распределенных на единичном отрезке. Эти числа нумеруются в порядке возрастания: ξ1 ξ2 ... ξm-1 . 2. Формируется вектор ξ = (ξ1 , ξ2 , . . . , ξm-1 ), где ξ1 = ξ1 , ξ2 = ξ2 - ξ1 , . . . , ξm-1 = ξm-1 - ξm-2 . Этот вектор, как показано в работе [18], равномерно распределен на множестве m-1 M= x∈R m-1 | x1 0, x2 0, . . . , xm-1 0, xi 1 . i=1 3. Все угловые точки многогранника M , кроме нулевой, лежат на рассто√ янии 2 друг от друга. Многогранник M растягивается до правильного ˜ . Для этого перенесем вершину, лежащую в начале координат, так, M чтобы √ расстояние от нее до всех остальных вершин составляло также 2. Это будет точка с координатами (x01 , x02 , . . . , x0m-1 ), где x0i = √ = r = (1 - m)/(m - 1), i = 1, 2, . . . , m - 1. Получился многогранник, конгруэнтный (m - 1)-мерному симплексу в m-мерном пространстве. ˜ , отображающую точку ξ мноОбозначим через ξ˜ точку множества M ˜ жества M . Точку ξ определяем в результате отражения ξ относительно гиперплоскости m-1 x ∈ Rm-1 | xi = 1 . i=1 Для этого сначала находим проекцию ξ на эту гиперплоскость: ξ˘i = ξi + ∆, i = 1, 2, . . . , m - 1, где ∆ определяем из условия m-1 ξ˘i = 1 i=1 757 Л и в ш и ц М. Ю., С и з и к о в А. П. по формуле 1 1- ∆= m-1 m-1 ξi . i=1 Отражение находим с учетом того, что коэффициент растяжения ра√ вен m: √ ˘ ξ˜ = ξ˘ + m(ξ - ξ). Это аффинное преобразование не оказывает влияния на закон распре˜. деления, поэтому точка ξ˜ имеет равномерное распределение в M ˜ 4. Точка ξ переносится в (m - 1)-мерный симплекс m-мерного пространства переменных задачи (1), (2). Для этого представляем ее как выпук˜ m-мерного пространства. лую линейную комбинацию угловых точек M Коэффициенты этой комбинации α1 , α2 , . . . , αm есть не что иное, как новые координаты точки α. Они определяются следующим образом:      ˜  ξ1 α1 1 0 ... 0 r -1  α2   0 1 ... 0 r  ξ˜2        ...   . . .  = . . . . . . . . . . . . . . .   . α  0  ˜ 0 . . . 1 r m-1 ξm-1  αm 1 1 ... 1 1 1 Теперь вернемся к методу Ψ-преобразования. Алгоритм 2. 1. Выполняется серия случайных проб, число которых N должно быть достаточным для получения надежных статистических оценок. На каждой пробе k = 1, 2, . . . , N с помощью алгоритма 1 генерируется допустимое значение аргумента x(k) и находится f (x(k) ). По результатам всех проб определяются среднее арифметическое fc и максимальное fm наблюдаемых значений целевой функции: 1 fc = N N f (x(k) ), fm = max{f (x(1) ), f (x(2) ), . . . , f (x(N ) )}. k=1 Отрезок [fc , fm ] разбивается на L (порядка десяти) уровней ζ = (ζ1 , ζ2 , . . . , ζL ), где ζl = fc + (l - 1)(fm - fc )L-1 . 2. Выполняется вторая серия случайных проб. Для каждого уровня l = = 1, 2, . . . , L подсчитывается nl - число проб, в которых выполнилось условие f (x(k) ) ζl , и средняя координата: 1 x ¯il = nl 758 N (k) δ l xi , k=1 i = 1,2, . . . , m, δl = 1, f (x(k) ) ζl , 0, f (x(k) ) < ζl . Об одном методе поиска глобального экстремума . . . По результатам формируются последовательности: -1 ψ = {ψl }L , l=1 , где ψl = nl N x ˜i = {¯ xil }L l=1 , i = 1, 2, . . . , m. (3) (4) 3. Определяются квадратичные тренды. Пусть   1 1 1 2 4 1 Λ= , . . . . . . . . . 1 L L2 тогда параметры a = (a0 , a1 , a2 ) тываются по формуле тренда ψ(l) = a0 + a1 l + a2 l2 рассчиa = (Λ Λ)-1 Λ ψ . Аналогично рассчитываются параметры bi = (bi0 , bi1 , bi2 ) тренда x ˜i (l) = = bi0 + bi1 l + bi2 l2 : bi = (Λ Λ)-1 Λ x ˜i . 4. Определяется наименьший положительный корень l∗ уравнения ψ(l) = 0. Его существование обусловлено тем, что ψ(l) на отрезке [1, L] представляет собой тренд последовательности (3), монотонно убывающей до значения, близкого к нулю. Подстановкой l∗ в тренды x ˜(l) находится приближенное решение задачи: x∗i = x ˜i (l∗ ), i = 1, 2, . . . , m. Назначение описанного алгоритма 2 - войти в область глобального экстремума. Предположим, что это удалось. Перейдем к локальному поиску самого экстремума. Методы, предполагающие использование производных первого или второго порядков (метод сопряженных градиентов, метод Ньютона), здесь неприменимы по причине возможной негладкости целевой функции f (x). В этом случае необходимо использовать метод, для реализации которого достаточно значений самой целевой функции f (x). Например, последовательное отражение симплексов [20]. Классический метод последовательного отражения симплексов - метод поиска экстремума в m-мерном евклидовом пространстве с помощью m-мерных симплексов [20]. В задаче (1), (2) речь идет о поиске экстремума в (m - 1)мерном симплексе с помощью (m - 1)-мерных же симплексов. Построим в m-мерном евклидовом пространстве правильный m-мерный симплекс, одна из вершин которого (назовем ее нулевой) лежит в начале координат. Пусть длина ребра равна ∆x, тогда координаты вершин соответствуют столбцам матрицы:     (0) (1) (2) (m) x1 x1 x1 . . . x1 0 q q ... q  (0) (1) (2) (m)  0 x2 q q ... q  x2 x2 . . . x2  ,  =  . . . . . . . . . . . . . . .  . . . . . . . . . . . . . . . (0) (1) (2) (m) 0 q q ... q xm xm xm . . . xm 759 Л и в ш и ц М. Ю., С и з и к о в А. П. где ∆x √ ∆x √ √ ( m + 1 + m - 1), q = √ ( m + 1 - 1). m 2 m 2 Перенесем этот симплекс из начала координат так, чтобы все его вершины, кроме x(0) , расположились на гиперплоскости q = x1 + x2 + · · · + xm = 1 и чтобы (m - 1)-мерный симплекс с вершинами x(1) , x(2) , . . . , x(m) оказался в области начала поиска. При параллельном переносе нулевой вершины в точку x ˜(0) координаты остальных вершин становятся следующими: x ˜(k) = x ˜(0) + x(k) , k = 1, 2, . . . , m. Найдем x ˜(0) из условия (k) (k) x ˜1 + x ˜2 + · · · + x ˜(k) m = 1. Достаточно взять какую-нибудь одну вершину, например, первую: m m (1) x ˜i (0) = i=1 x ˜i + i=1 m (0) = x ˜i + i=1 ∆x √ √ ( m + 1 + m - 1) + m 2 ∆x √ √ m+1+m-1+ m 2 m i=2 m-1 ∆x √ √ ( m + 1 - 1) = m 2 √ ( m + 1 - 1) i=1 m (0) = x ˜i + i=1 Отсюда получаем m (0) x ˜i i=1 = ∆x √ √ m m + 1 = 1. m 2 ∆x √ = 1 - √ m + 1. 2 Определим начальное расположение m-мерного симплекса, привязав его к точке x∗ , полученной методом Ψ-преобразования. В этом случае координаты нулевой вершины симплекса определятся следующим образом: (0) x ˜i ∆x √ = x∗i 1 - √ m + 1 , 2 i = 1, 2, . . . , m. Координаты остальных вершин, которые и образуют начальный «рабочий» симплекс, представлены столбцами матрицы:   (1) (2) (m) x ˜1 x ˜1 ... x ˜1  (1) (2) (m)  x x ˜2 ... x ˜2   ˜2 ,  ... ... ... ...  (1) (2) (m) x ˜m x ˜m . . . x ˜m 760 Об одном методе поиска глобального экстремума . . . где (k) x ˜i = x∗i √ ∆x √ ( m m 2 √ ∆x √ ( m m 2 ∆x √ 1- √ m+1 + 2 + 1 + m - 1), i = k, + 1 - 1), i = k. Процесс поиска состоит в реализации последовательности шагов, на каждом из которых происходит построение нового «рабочего» симплекса путем отражения вершины с минимальным значением целевой функции f (x) от противоположной грани исходного симплекса. Пусть v ∈ [1, m] - номер такой вершины, тогда отраженная вершина определяется по формуле x ˘(v) = m 2 m-1 x(k) - x(v) - x(v) . k=1 Покажем, что x ˘(v) лежит на гиперплоскости x1 + x2 + · · · + xm = 1. Действительно, m (v) x ˘i i=1 = = 2 m-1 2 m-1 m m m (k) xi i=1 k=1 m m (k) xi k=1 i=1 (v) - xi (v) - xi=1 = i=1 m m (v) - xi i=1 (v) - xi i=1 = 2 (m - 1) - 1 = 1. m-1 Поскольку поиск проводится в ограниченном пространстве, может оказаться, что x ˘(v) выходит за допустимую область. Предположим для простоты, что на текущем шаге алгоритма выполняется f (x(1) ) < f (x(2) ) < · · · < f (x(m) ). Тогда нужно последовательно перебирать вершины до тех пор, пока не получим допустимого отражения. Отражение лучшей точки будет признаком окончания текущей серии. Следующая серия начинается с уменьшения длины ребра ∆x и переноса нового «рабочего» симплекса в область лучшей точки предыдущей серии. Описанная методика может быть использована для поиска оптимального компонентного состава различных смесей, в частности, для расчета рецептур смешения товарных нефтепродуктов. Некоторые показатели качества смесей рассчитываются по специальным, иногда довольно сложным, эмпирическим формулам [21, 22]. Пусть Q - множество индексов контролируемых показателей некоторой смеси, и пусть φq (x) - зависимость (линейная или нелинейная) q-того показателя смеси от компонентов смешения, взятых в пропорциях x = (x1 , x2 , . . . , xm ). Тогда задачу смешения можно сформулировать так: m C= cj xj → min, j=1 761 Л и в ш и ц М. Ю., С и з и к о в А. П.  + φq (x) ∈ [p- q ∈ Q,  q , pq ],   m xj = 1,  j=1   xj 0, j = 1, 2, . . . , m, + где cj - цена j-того компонента; p- q , pq - граничные значения для q-того показателя. Требования по качеству смеси удобно учесть в виде штрафной составляющей целевой функции: m f (x) =   j=1 m (5) ϕ2q (x) → min, cj xj + µ q∈Q xj = 1, j=1  xj 0, j = 1, 2, . . . , m, где µ - весовой коэффициент штрафной составляющей; ϕq (x) - процентное отклонение расчетного значения q-того показателя от границы заданного интервала, определяемое по формуле  2 - -  10 (pq - φq (x))/p- q , pq - φq (x) > 0, + 2 + 10 (φq (x) - pq )/pq , p+ ϕq (x) = q - φq (x) < 0,  + 0, φq (x) ∈ [p- q , pq ]. Для примера решена задача составления трехкомпонентной смеси со следующими условиями. Стоимость C = 1.5x1 + 1.3x2 + 1.0x3 . Требования к показателям (октановому числу, плотности, содержанию серы): φ1 (x) = 90.2x1 +88.5x2 +73.6x3 -14.3x1 x2 -9.8x1 x3 +28.4x1 x2 x3 φ2 (x) = 0.82x1 + 0.59x2 + 0.67x3 0.68, φ3 (x) = 0.001x1 + 0.003x2 + 0.002x3 0.002. 85, (6) При составлении функции f (x) в (5) используется весовой коэффициент µ = 10. По результатам предварительного статистического исследования области изменения f (x) сформировано 10 уровней с шагом ∆ζ = -70, начиная с ζ1 = 700. Для каждого из них выполнено 1000 статистических испытаний. Получены последовательности (3), (4) и построены тренды (см. рисунок): a) ψ(l) = 0.6916 - 0.02906l - 0.00374l2 , b) x ˜1 (l) = 0.33041 - 0.00145l + 0.000098l2 , c) x ˜2 (l) = 0.327667 + 0.013837l - 0.00072l2 , d) x ˜3 (l) = 0.341923 - 0.012392l + 0.000626l2 . Из уравнения ψ(l) = a0 + a1 l + a2 l2 = 0 получено l∗ = 10.2554. Подстановкой l∗ в тренды x ˜1 (l), x ˜2 (l), x ˜3 (l) найдено приближенное решение: x∗1 = 0.3259, 762 x∗2 = 0.3934, x∗3 = 0.2807. Об одном методе поиска глобального экстремума . . . Тренды результатов статистических испытаний [Trends of results of statistical tests] Значение целевой функции (5): f (x∗ ) = 115.3. Уточнение решения осуществляется по описанному выше модифицированному методу отражения правильных симплексов. Проведено 20 серий с уменьшением длины ребра симплекса по формуле ∆xn = 0.1n-1 , где n - номер серии. Результаты первого приближения x∗ и уточненного x ˘ решения представлены в таблице ниже. ∗ x x ˘ x1 x2 x3 φ1 φ2 φ3 f (x) 0.3259 0.3702 0.3934 0.3990 0.2807 0.2308 83.16 83.71 0.687 0.693 0.00209 0.00206 115.3 87.6 Как видно из таблицы, уже на первом этапе получено достаточно хорошее приближение. В результате реализации алгоритма выявлена несовместность системы ограничений (6) и получено решение по критерию (5), основная составляющая которого - взвешенная сумма квадратов процентных отклонений показателей от требуемых значений. 763 Л и в ш и ц М. Ю., С и з и к о в А. П. Заключение. В статье предложен метод Ψ-преобразования, модифицированный для решения задачи поиска глобального экстремума непрерывной функции на симплексе, имеющей важное прикладное значение. Метод является альтернативой прямому случайному поиску, эффективность которого существенно падает с увеличением числа переменных. Объем вычислений, реализуемых в рамках метода Ψ-преобразования, напротив, мало зависит от размерности задачи. Он определяется в основном числом статистических испытаний, которое, по данным [16], не превышает в общей сложности нескольких тысяч и не требует больших затрат машинного времени. Декларация о финансовых и других взаимоотношениях. Исследование не имело спонсорской поддержки. Все авторы принимали участие в разработке концепции статьи и в написании рукописи. Авторы несут полную ответственность за предоставление окончательной рукописи в печать. Окончательная версия рукописи была одобрена всеми авторами. Авторы не получали гонорар за статью.

About the authors

Aleksandr P Livshits

Samara State Technical University

Email: mikhaillivshits@gmail.com
(Dr. Techn. Sci.; mikhaillivshits@gmail.com), Head of Dept., Dept. Of Management and System Analysis of Thermal Power and Socio-Technical Systems 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

Aleksandr Pavlovich Sizikov

Samara State Technical University

Email: apsizikov@mail.ru
(Cand. Econ. Sci.; apsizikov@mail.ru; Corresponding Author), Doctoral Student, Dept. of Management and System Analysis of Thermal Power and Socio-Technical Systems 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

References

Supplementary files