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

Abstract


A non-convex problem of mathematical programming is considered, which permissible region is a simplex. A two-stage algorithm is proposed for approximate solution of the problem. The region of global optimum is determined using the Ψ-transform method at the first stage; local “fine-tuning” of the solution is performed at the second stage. The Ψ-transform was modified taking into account the special features of the problem under consideration. Ψ-function is determined according to the results of statistical tests implemented using the generator of random points uniformly distributed over the simplex. The proposed method of reflection of regular simplexes is used for fine-tuning of the solution. An example of application of the developed algorithm for solving the problem of optimization of component composition of the hydrocarbon mixture is presented.

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
244, Molodogvardeyskaya st., Samara, 443100, Russian Federation
(Dr. Techn. Sci.; mikhaillivshits@gmail.com), Head of Dept., Dept. Of Management and System Analysis of Thermal Power and Socio-Technical Systems

Aleksandr Pavlovich Sizikov

Samara State Technical University

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

References

  1. Zhigljavsky A., Žilinskas A. Stochastic Global Optimization / Springer Optimization and Its Applications. vol. 9. Berlin: Springer, 2008. xiii+262 pp. doi: 10.1007/978-0-387-74740-8.
  2. Marti K. Stochastic Optimization Methods / Applications in Engineering and Operations Research. Berlin: Springer, 2015. xxiv+368 pp. doi: 10.1007/978-3-662-46214-0.
  3. Пантелеев А. В. Метаэвристические алгоритмы поиска глобального экстремума. М.: МАИ Принт, 2009. 159 с.
  4. M. Tim Jones AI Application Programming / Programming Series. Boston: Charles River Media, 2003. 496 pp.
  5. Савин А. Н., Тимофеева Н. Е. Применение алгоритма оптимизации методом имитации отжига на системах параллельных и распределённых вычислений // Изв. Сарат. унта. Нов. сер. Сер. Математика. Механика. Информатика, 2012. Т. 12, № 1. С. 110-116.
  6. Botev Z. I., Kroese D. P. Global likelihood optimization via the cross-entropy method, with an application to mixture models / Proceedings of the 2004 Winter Simulation Conference. Washington: IEEE, 2004. pp. 529-535. doi: 10.1109/wsc.2004.1371358.
  7. Ernst D., Glavic M., Stan G.-B., Mannor S., Wehenkel L. The cross-entropy method for power system combinatorial optimization problems / 2007 IEEE Lausanne Power Tech. Washington: IEEE, 2007. pp. 1290-1295. doi: 10.1109/pct.2007.4538502.
  8. Evans G. E., Keith J. M., Kroese D. P. Parallel cross-entropy optimization / 2007 Winter Simulation Conference. Washington: IEEE, 2007. pp. 2196-2202. doi: 10.1109/wsc.2007.4419854.
  9. Zhigljavsky A. Theory of Global Random Search / Mathematics and Its Applications (Soviet Series). vol. 65. Berlin: Springer, 1991. xviii+341 pp. doi: 10.1007/978-94-011-3436-1.
  10. Феоктистов А. Г., Горский С. А. Реализация метода мультистарта в пакете Градиент // Вестник НГУ Серия: Информационные технологии, 2007. Т. 5, № 2. С. 78-82.
  11. Cohoon J., Karro J., Lienig J. Evolutionary Algorithms for the Physical Design of VLSI Circuits / Advances in Evolutionary Computing / Natural Computing Series; eds. A. Ghosh, S. Tsutsui. Berlin: Springer, 2003. pp. 683-712. doi: 10.1007/978-3-642-18965-4_27.
  12. Гладков Л. А., Курейчик В. В., Курейчик В. М. Генетические алгоритмы. М.: Физматлит, 2006. 320 с.
  13. Гладков Л. А., Гладкова Н. В. Особенности использования нечетких генетических алгоритмов для решения задач оптимизации и управления // Известия ЮФУ. Технические науки, 2009. № 4(93). С. 130-136.
  14. Skiena S. S. The Algorithm Design Manual. London: Springer, 2008. xvi+730 pp. doi: 10.1007/978-1-84800-070-4.
  15. Chichinadze V. K. Solution of nonlinear nonconvex optimization problems by Ψ-transformation method // Computers & Mathematics with Applications, 1991. vol. 21, no. 6-7. pp. 7-15.
  16. Чичинадзе В. К. Решение невыпуклых нелинейных задач оптимизации. Метод Ψпреобразования. М.: Наука, 1983. 256 с.
  17. Ахмадиев Ф. Г., Гильфанов Р. М. Математическое моделирование и оптимизация “состав-свойство” многокомпонентных смесей // Известия КГАСУ, 2012. № 2(20). С. 289-297.
  18. Новоселов А. А. Равномерное распределение на стандартном симплексе в Rn , http: //risktheory.novosyolov.com/lectures/unifs.pdf (дата обращения: 23.10.2016).
  19. Русин Ю. В. Алгоритмы статистического моделирования вероятностных распределений. Ярославль: ЯрГУ, 2006. 58 с.
  20. Зедгинидзе И. Г. Планирование эксперимента для исследования многокомпонентных систем. М.: Наука, 1976. 390 с.
  21. Чинакал В. О. Оптимизация рецептуры светлых нефтепродуктов / Оптимизация, исследование операций, бионика. М.: Наука, 1973. С. 198-205.
  22. Никитин В. А., Мусаев А. А. Оптимизация компаундирования углеводородных смесей // Тр. СПИИРАН, 2007. Т. 4. С. 327-336.

Statistics

Views

Abstract - 11

PDF (Russian) - 9

Cited-By


PlumX

Dimensions

Refbacks

  • There are currently no refbacks.

Copyright (c) 2016 Samara State Technical University

Creative Commons License
This work is licensed under a Creative Commons Attribution 4.0 International License.

This website uses cookies

You consent to our cookies if you continue to use our website.

About Cookies