Приближенные алгоритмы поиска граничных точек для задачи условной псевдобулевой оптимизации

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

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

Универсальные поисковые алгоритмы, такие как генетический алгоритм (ГА), эволюционный алгоритм, алгоритм имитации отжига, могут применяться для решения любых задач оптимизации. Но для некоторых практических задач они не находят приемлемого решения. Это происходит при наличии сложной системы ограничений на переменные и, как следствие, малой допустимой области. Для таких случаев необходима разработка алгоритмов, учитывающих специфику решаемой задачи.

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

Приведем определения, необходимые для описания работы алгоритмов оптимизации [1]:

• -    псевдобулевой функцией будем называть вещественную функцию на множестве булевых переменных: f : B ₂ ⁿ ^ R * , где B ₂ = { 0,1 } ; B 2 = B 2 х B 2 х ... х B ₂ ;

• -    точки X ¹, X ² е B 2 назовем k-соседними , если они отличаются значением к координат, где к = 1, n ; I- соседние точки будем называть просто соседними ;

• -    множество O k ( X ) , к = 0, n , всех точек B 2 , к -соседних к точке А , назовем к-м уровнем точки А ;

• -    множество точек W ( X °, X l ) = { X °, X 1, . . ., X l } с B 2 - это путь между точками X ⁰ и X¹ , если V i = 1,..., l точка X ' является соседней к X ' ^{- 1} ;

• - множество A с B 2 назовем связным множеством , если V X °, X l е A существует путь W ( X ⁰, X l ) с A ;

точка X * е B 2 , для которой f ( X ) <  f ( X ), V X е O 1 ( X ) , является локальным минимумом псевдобулевой функции / ,

• -    псевдобулевую функцию, имеющую только один локальный минимум, будем называть унимодальной на множестве B 2 функцией;

• -    унимодальную функцию. / "назовем монотонной на B 2 , если V X^k е O k ( X *), к = 1, n :

f ( X k ^{- 1}) <  f ( X k ), V X k ^{- 1} е O k -1 ( X * ) n O 1 ( X^k ) .

Пример . Полином от булевых переменных

* X1,-..,x.) = 2>Д X, j=1    'ев j где вj с {1,..., n}, является унимодальнойи_монотонной псевдобулевой функцией при aj > 0, j = 1, m.

Постановка задачи. Рассмотрим задачу следующего вида:

С ⁽ X ) ^ шах . ,                (1)

X е о с B 2

где С ( X ) - монотонно возрастающая от X ⁰ псевдобу-левая функция; О с B . - некоторая подобласть пространства булевых переменных, определяемая заданной системой ограничений, например:

A j (X ) <  H j , j = 1 m .             (2)

В общем случае множество допустимых решений 5 является несвязным множеством.

Множество допустимых решений задачи (1) обладает следующими свойствами:

• - точка Y _е A является граничной точкой множества А , если существует X е O 1 ( Y ) , для которой X й A ;

• -    точку Y е O ' ( X ⁰) n A будем называть крайней точкой множества А , с базовой точкой X ⁰ е A , если V X е O 1 ( Y ) П O ' _{+ 1} ( X ⁰) выполняется X й A ;

• -    ограничение, определяющее подобласть пространства булевых переменных, будем называть активным , если оптимальное решение задачи условной оптимизации не совпадает с оптимальным решением соответствующей задачи оптимизации без учета ограничения.

Рассмотрим другие свойства множества допустимых решений [2]:

• -    если целевая функция является монотонной унимодальной функцией, а ограничение активно, то оптимальным решением задачи (1) будет точка, принадлежащая подмножеству крайних точек множества допустимых решений 5 с базовой точкой X ⁰ , в которой целевая функция принимает наименьшее значение:

С ( X ⁰) = min С ( X ) ;

• X е B 2.

• -    если функция ограничения (2) является унимодальной функцией, то множество допустимых решений 5 задачи (1) представляет собой связное множество;

• -    количество крайних точек связного множества допустимых решений задачи (1) 5 <  5 = С^{n /2]}, где [ n / 2 ] - целая часть числа n / 2 .В случае 5 = 5 _max все крайние точки принадлежат O _{n / 2}( X °) при четном п и O _{( n} _ _1)/2( X °) (или O _{( n + 1)/2}( X °) ) при нечетном п .

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

Рассмотрим обобщенную функцию

m

F ( X ) = C ( X ) - r £ max{0, A j ( X ) - H j }     (3)

^j =1

для задачи (1) с ограничениями (2).

В этом случае имеет место следующее утверждение [3]:

Утверждение. Крайние точки множества допустимых решений 5 исходной задачи (1) с ограничениями (2) с базовой точкой X0 являются точками локальных максимумов для обобщенной штрафной функции (3) при параметре г, удовлетворяющем условию r >    C (Z) - C (X')

£ max{0, Aj( Z) - H} j=1

для любой крайней точки X' е O _k ( X ⁰) и любой точки Z е O _{k +} i ( X °) n O i ( X ') .

Таким образом, задача (1) идентична задаче

F ( X ) ^ max ,               (4)

X е ВП в которой количество локальных максимумов (в случае связности множества допустимых решений) 5 < 5max = Cn/2]. В общем случае, когда допустимое множество не является связным, число локальных максимумов теоретически может достигать 2”-1.

Задача (4) решается известными поисковыми методами: локального поиска, генетическим алгоритмом, алго ритмом имитации отжига и т. д.

Существенным недостатком такого подхода является потеря свойства монотонности целевой функции C ( X ) . При введении даже простых (например, линейных) ограничений обобщенная функция будет полимодальной немонотонной функцией с экспоненциальным числом точек локальных максимумов.

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

Общая схема прямого алгоритма поиск а следующая:

• -    шаг 1. Полагаем X 1 = X ⁰ , i = 1 ;

• -    шаг 2. В соответствии с правилом выбираем X i _{+ 1} е O i ( X ⁰) n O 1 ( X i ) n S . Если таких точек нет, то переходим на шаг 3; иначе i = i + 1 и повторяем шаг 2;

• - шаг3. X _opt = X +1 .

Общая схема двойственного алгоритма поиска такова:

• -    шаг 1. Полагаем X 1 е O n ( X ⁰) , i = 1 ;

• -    шаг 2. В соответствии с правилом выбираем X i _{+ 1} е O n - j ( X ⁰) n O 1 ( X i ) . Если X, _{+ 1} е S , то переходим на шаг 3; иначе i = i + 1 и повторяем шаг 2;

• - шаг3. X opt = X +1 .

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

Рассмотренные далее алгоритмы поиска граничных точек отличаются друг от друга лишь правилом выбора следующей точки на шаге 2 общих схем, прямого и двой ственного алгоритмов поиска.

Будем использовать следующие правила поиска:

• -    правило 1 - случайный поиск граничных точек (СПГ). Точка X i _{+ 1} выбирается случайным образом. Каждая точка на следующем шаге может быть выбрана с равной вероятностью. При решении практических задач эти вероятности могут быть не равными, а вычисляться на основании специфики задачи до начала поиска;

• -    правило 2 - Гриди-алгоритм. Точка X i _{+ 1} выбирается по условию

^ ⁽ X i +D = ^max ^ ⁽ X^{j )} ,

j где Xj е Oi (X0) n O1 (Xi) n S для прямого алгоритма и Xj е On - i (X0) n O1 (Xi) для двойственного алгоритма.

Функция % ( X ) выбирается исходя из специфики задачи, например: целевая функция % ( X ) = C ( X ) , удельная ценность % ( X ) = C ( X )/ A ( X ) (для одного ограничения) и т. д.;

• -    правило 3 - адаптивный случайный поиск граничных точек (АСПГ). Точка X i _{+ 1} выбирается случайным образом в соответствии с вектором вероятностей

P = (p1, p2,’”, pj ) , где J - количество точек, из которых производится выбор;

i   Ц Xj )     —

P j = j , j = 1, J ,

^ X ( X l )

l =1

здесь Xj е Oi (X0) n O1 (Xi) n S для прямого алгоритма и Xj е On - i (X0) n O1( Xi) для двойственного алгоритма. Алгоритм АСПГ является дополнением к Гриди-алгорит му;

• -    правило 4 - модифицированный случайный поиск граничных точек (МСПГ). Точка X i _{+ 1} выбирается по условию

^ ^{( X} i +1 ) = ^max ^ ^{( X} r ) ,

r где Xr -точки, выбранные по правилу 1, r = 1,R , здесь R - задаваемый параметр алгоритма.

Гриди-алгоритм является регулярным алгоритмом, т. е. при повторном запуске из одной и той же точки он выдает одинаковое решение. Остальные алгоритмы можно запускать несколько раз и из найденных решений выб- рать наилучшее. Трудоемкость каждого запуска алгоритма ограничена:

T 5 n (n +1) " 2 ’ но средняя трудоемкость Гриди-алгоритма и АСПГ значительно выше остальных, так как на каждом шаге они просматривают все соседние точки следующего уровня, в отличие от СПГ и МСПГ, которые в двойственной схеме просматривают на каждом шаге лишь одну и R точек со ответственно.

Рассмотрим применение описанных алгоритмов на одной из практических задач большой размерности с несвязным множеством допустимых решений.

Задача оптимизации загрузки производственных мощностей литейных отделений. В литейных отделениях

(ЛО) производится продукция различного вида. В каждом ЛО имеется специализация по видам продукции, которую могут производить на его литейных машинах (ЛМ). Имеется некоторое количество заказов на производство продукции. Каждому заказу соответствует объем, вид продукции и срок выполнения. Вид продукции характеризуется производительностью за смену (всего три смены в сутки). Замена производимой на ЛМ продукции требует ее перенастройки, занимающей одну смену Необходимо загрузить производственные мощности ЛО таким образом, чтобы выполнялись все заказы, продукция производилась более равномерно во времени и число перенастроек оборудования ЛМ было минимальным.

Входные данные:

- 1 - количество дней планирования;

• -    J - количество ЛО;

  - K j - количество ЛМ в] -мЛО, j = 1, J ;

  - L - количество заказов на продукцию, производи

мую на ЛМ (и соответствующее число видов продукции);

- V l - производительность / -го вида продукции за смену на ЛМ, l = 1, L ;

• -    T l - срок выполнения / -го заказа (на производство I -го вида продукции) на ЛМ, l = 1, L ;

• -    W _l - объем / -го заказа (на производство / -го вида продукции) на ЛМ, l = 1, L ;

z_jt - характеризует специализацию ЛО:

• 1,    если ЛМ j -го ЛО способны

  z jl =1

  
  

производить продукцию l -го вида,

0 в противном случае;

- а - коэффициент жесткости ограничения по требо ванию равномерности производства продукции по дням, 0 < а < 1.

Для построения модели введем следующие булевы переменные:

Y = ^{ y jkl ^{} е} B 2 ” ; где В 2 = { °,1 } , B 2 = B 2 х B 2 хЛх B 2 - множество булевых переменных;

X = { x jkl } е B П ;

1, если в i -й день на k -й ЛМ j -го ЛО

У уи = производится продукция l -го вида,

0 в противном случае;

1, если в г-й день на k-й ЛМ j -го ЛО xjkl = начинается производство продукции l-го вида,

0 в противном случае;

Общая размерность булева вектора У(иХ)

J n = 1 • L X K.

j =1

Дополнительные замечания:

• 1)    j ⁵ Ууи V i , j , k , l ;

• 2)    X jkl = У у_И • ⁽¹ — y , —1, jkl ) V i , j , k , l

(y0jkl = 0 Vj, k, l).~ (5)

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

С(X) ^ min ,(6)

A1(Y) > w„i = й,

Ai2(Y) >a^ W1, i = 1,1, ае (0,1), W = X W I, (8) l=1/

IJKjL       IJKjL где С (X) = XXXX xjkl = XXXX yju (1 — У.-1,jkl);

i =1 j =1 k =1 l =1              i =1 j =1 k =1 l =1

T l J K j

A ? ( ^Y ) = XXX V l • Ууи ( ² + y i-1, jkl ) ;

i =1 j =1 k =1

J ^K j L

A - ²( Y ) = XXX V • y jkl ( 2 + У . -1, jkl ) ;

j =1 k =1 l =1

y 0jkl = 0 V j , k , l ;

• 2)    дополнительные ограничения:

LL

X y k ^{5 1} ( X X jkl ^{5 1} ), i = 1, 1 , j = 1 J , k = L K j , (9) l =1                      l =1

yjkl 5 zjl (xijkl 5 zjl), i = U, j = 1J, k = 1Kj , l = 1L.         (10)

Модель обладает следующими свойствами:

• -    имеется два пространства булевых переменных (обозначим их B X и B Y ), соответствующих векторамКи У . Каждой точке Y е B Y соответствует единственная точка X е B^X , компоненты которой определяются соотношением (5). При этом точке X е B^X может соответствовать несколько точек Y е B Y (с разным значением функций ограничений);

• -    целевая функция (6) является линейной и унимодальной монотонной в пространстве B X с точкой минимума X ⁰ = (0,..., 0) . В пространстве B Y целевая функция является квадратичной и унимодальной немонотонной с точкой минимума Y ⁰ = (0,..., 0) ;

• -    функции ограничений (7) и (8) в пространстве B Y являются квадратичными и унимодальными монотонными псевдобулевыми функциями с точкой минимума Y ⁰ = (0,..., 0) . В пространстве B функции ограничений однозначно не определены;

• -    множество допустимых решений в пространствах

J

B^x и B Y ограничено сверху 1 X K j -м уровнем точки j =1

минимума ( X ⁰ и Y ⁰ ) согласно ограничению (9). В пространстве B Y этот уровень соответствует тому, что в любой день на каждой литейной машине производится продукция;

• -    множество допустимых решений в общем случае является несвязным множеством (в пространстве B Y ).

Таким образом, решение задачи полностью определяется переменными У, чего нельзя сказать про перемен-ныеК. Но целевая функция от А'обладает хорошими конструктивными свойствами, в связи с чем поиск оптиму- ма по А' является более эффективным, чем по Y. Но так как функции ограничений (7), (8) от А'не определены, то следует находить значения этих функций от соответствующей точки Y. Таких точек может быть несколько:

г у у у

• ⁷ 1 , 2 , ..., н ,

причем в некоторых из них решение может быть допустимым, а в некоторых нет. А так как функции ограничений здесь монотонны, то для конкретной А следует выбирать Y , принадлежащую наиболее возможному уровню (с большими значениями функций):

(„    "I

• ^Y = arg y m ax ^ y „ , h         ^ i , j , k ,          J

^Где Y h = ⁽ y llll , ..., yU KJjL )

Ниже приведен один из алгоритмов такого преобразования ( алгоритм 1 ):

• -    шаг 1. Полагаем, что N _jk = 0 , j = 1, J , k = 1, K j ,

i = 1;                      ____ ______ ____

• -    шаг 2. Для j = 1, J , k = 1, K j , l = 1, L выполняем условие: если x _ijkl = 1 , то N _jk = l ;__

• -    шаг 3. Для j = 1, J , k = 1, K j , l = 1, L выполняем условие: если N_jk = l , то y^ y = ¹ , иначе y_ukl = 0 ;

• -    шаг 4. Если i < I , то i = i + 1 и переходим на шаг 2.

В то же время решение Y, полученное из найденного наилучшего вектора Xopt таким способом, может соответствовать ситуации, когда количество выпускаемой продукции намного превосходит требуемое значение (что, в общем-то, не противоречит построенной модели, но может влиять на равномерность загрузки производственных мощностей, которая оптимизируется на следующем этапе). Поэтому после завершения первого этапа поиска Yopt следует определять по правилу

У у у у opt ⁷ ^1 ,^ 2 , ■■■,^ H , f^   "

Y opt = arg , ^min - Z У ум .

Y h : A l ( Y h )> W l , l =1, L ^ i , j^ , _l j

В этом случае алгоритм преобразования несколько отличается от предыдущего ( алгоритм 2 ):

• -    шаг 1. Полагаем, что N _jk = 0 , j = 1, J , k = 1, K j ;

• -    ша г 1а. Полагаем, что y _jjH = 0 , i = 1, I , j = 1, J ,

k = 1, K j , i = 1;       ___ _____ ___

• -    шаг 2. Для j = 1, J , k = 1, K j , l = 1, L выполняем условие: если x_ijkl = 1 , то N _jk = l ;

• -    шаг 3. Для j = 1, J , k = 1, K j , l = 1, L выполняем условие: если N _jk = l и A _l ( Y ) <  W _l , то y _ljkl = 1 ;

• -    шаг 4. Если i < I , то i = i + 1 и переходим на шаг 2.

Здесь в шаг 3 добавлено условие A1 (Y) < Wt, а также введен шаг 1а для возможности такого вычисления. На добности в таком преобразовании во время поиска нет, оно необходимо только для определения конечного Yopt.

Для решения задачи оптимизации загрузки производственных мощностей использовались двойственные алгоритмы: СНГ, Гриди-алгоритм и МСПГ. Алгоритм АСПГ не рассматривался из-за его чрезмерной трудоемкости при многократном запуске. Единичный же запуск АСПГ очень редко может выдать решение лучше, чем решение Гриди-алгоритмом. Стартовой точкой поиска являлась точка безусловного минимума целевой функции X ⁰ = (0,..., 0) . Найденное решение улучшалось соответствующим прямым алгоритмом.

Кроме того, задача решалась генетическим алгоритмом. При реализации ГА была выбрана схема, которая эффективно работала при неоднократном решении других задач комбинаторной оптимизации [4; 5].

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

Осредненные результаты решения 10 задач месячного планирования загрузки производственных мощностей. приведены в таблице. Средние значения входных данных: I = 31 ; J = 3 ; K 1 = 12 ; K 2 = 9 ; K ₃ = 7 ; L = 36 ; а = 0,5 ; V е [40, 50] ; W l е [20, 25 000] . При этом общая размерность булева вектора n = 31 248 .

Количество запусков L алгоритмов выбиралось таким образом, чтобы трудоемкость приблизительно равнялась трудоемкости единичного запуска Гриди-алгоритма. В данном случае трудоемкость T = 8 • 10⁵ .

Модифицированный случайный поиск граничных точек по сравнению с Гриди-алгоритмом является более гибкой процедурой, так как позволяет выбирать параметры L и R , влияющие на трудоемкость алгоритма и точность решения. Гриди алгоритм такой возможности не предоставляет, и при больших размерностях время поиска может быть чрезмерно велико. Что касается эффективности, то при примерно одинаковой трудоемкости точность найденных решений различается несущественно.

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

Алгоритм	Количество запусков L	Найденное решение С opt
спг	2 000	Не найдено
Гриди	1	49
МСПГ, R =1000	12	47
МСПГ, R =100	60	52
МСПГ, R =10	200	Не найдено
ГА*	^^^^^в	Не найдено

* Параметры ГА: турнирная селекция с величиной турнира 5, размер популяции 100, наибольшее число поколений 8 000, вероятность мутации 0,000 1.

смотренной задачи являются двойственные алгоритмы МСПГ и Гриди.