Введение При решении задач принятия решений часто возникают ситуации, когда результирующее решение выбирается из большого набора вариантов (альтернатив). Формализованное описание подобных задач часто осуществляется путем введения бинарных переменных по каждой альтернативе, которые могут принимать только два значения, например 0 и 1. Если целевые функции и ограничения в подобной задаче могут быть представлены аналитическими выражениями, то полученные формализованные модели называются задачами математического бинарного (или булевого) программирования. Наиболее важными задачами бинарного программирования являются задачи линейного и квадратичного бинарного программирования, когда целевые функции и ограничения описывают линейными и квадратичными функциями соответственно, поскольку во многих приложениях задач математического программирования использование более сложных функций не позволяет повысить качество результатов ввиду наличия больших погрешностей в исходных данных, сравнимых с погрешностями самой модели. Вышеупомянутые задачи являются частным случаем задач целочисленного математического программирования, и поэтому, казалось бы, могут быть решены методами, используемыми при решении задач целочисленного программирования. Однако одна важная особенность делает все известные методы решения задач целочисленного программирования практически неприемлемыми и бесполезными при решении указанных задач: число альтернатив часто бывает настолько большим, что существующие вычислительные средства оказываются неспособными за разумное время найти решение указанных задач. Так, например, в работах [1, 2] приведены формализованные постановки задачи составления учебного расписания, когда в качестве бинарных переменных берутся следующие: , если занятие по d-й дисциплине у g-й подгруппы проводится в a-й аудитории на p-й паре, и в противном случае. Это означает, что в данной задаче индекс суммирования i представлен набором pgad. Тогда, как показано в [1], в этом случае число всех возможных переменных для типового вуза с общим числом студентов порядка 7 500 достигает 150 миллионов переменных , что делает практически невозможным автоматизированное решение указанной задачи на основе алгоритмов целочисленного программирования. Бинарное программирование как самостоятельное научное направление исследований было выделено еще в 60-е гг. ХХ в.; наиболее важные методы решения задач бинарного программирования, применявшиеся в то время, можно найти в [3]. В дальнейшем сколь-нибудь значимых результатов по решению задач бинарного программирования с большим числом переменных, по-видимому, получить не удалось. Состояние исследований в данном направлении в настоящее время отражено в [4, 5]. Отметим также отсутствие значимых результатов по использованию современных эвристических методов оптимизации (генетических алгоритмов, нейронных сетей, методов муравьиных или пчелиных колоний и др.), которые доказали свою эффективность при решении многих других типов задач оптимизации. Ниже нами предлагаются набор условий для нахождения оптимального решения и построенный на их основе эвристический алгоритм решения. 1. Формализация задачи с использованием метода штрафных функций В анализируемой ниже задаче квадратичного программирования ограничения заданы в виде равенств, в то время как в классической постановке эти ограничения имеют форму неравенств. Однако способы сведения ограничений в виде неравенств к ограничениям в виде равенств и обратно известны, поэтому предлагаемая постановка охватывает все классические постановки задач квадратичного программирования. Приведем формализованную постановку задачи. Пусть задан набор переменных xi (), каждая из которых принимает только значения 0 или 1, т. е. является бинарной переменной; . Тогда задачей квадратичного бинарного программирования является следующая оптимизационная задача: (1) при условии , . (2) Для нас представляет интерес также несколько другая запись целевой функции (1): (3) т. е. квадратичное выражение в целевой функции (1) представляется как сумма из S отдельных квадратичных выражений. Именно в таком виде представлена целевая функция в задаче формирования расписания занятий в [1], где S = 19, каждое квадратичное выражение в сумме по S имеет относительно простой вид и соответствует разным ограничениям задачи, число которых составляет порядка 30. Отметим, что поскольку общее число возможных вариантов решений задачи (1), (2) ограничено (не превосходит 2N), то задача (1), (2) всегда имеет решение. Без ограничения общности можно считать C = 0. Как было указано выше, для решения задачи (1), (2) предлагается использовать метод штрафных функций - т. е. решение задачи (1), (2) заменяется решением следующей задачи безусловной оптимизации (т. е. задачи без ограничений): (4) где () - вспомогательные переменные; Lk > 0 трактуется как величина штрафа за нарушение k-го ограничения: если хотя бы одно из слагаемых в последней сумме не равно нулю (т. е. положительно), то значение целевой функции при больших значениях Lk резко уменьшается. Поэтому оптимальное решение вынуждено стремиться обеспечить выполнение ограничений (2), при которых все слагаемые последней суммы равны нулю. Известно, что при достаточно больших значениях Lk для всех k решения задач (4) и (1), (2) совпадают, поэтому вместо решения задачи (1), (2) можно решать задачу (4). В записи задачи (4) учтены все ограничения задачи (1), (2), кроме самого важного - переменные xi являются бинарными. Для учета данного ограничения предлагается добавить в целевую функцию задачи еще по одному слагаемому на каждую переменную: (5) где L >> max(Li). Действительно, если для некоторого индекса i выражение ) отличается значимо от нуля, в целевой функции появляется большое по модулю отрицательное слагаемое (последняя сумма в правой части), которое сильно уменьшает значение целевой функции по сравнению со случаем ) = 0. В результате исходная задача квадратичного программирования переходит в задачу полиномиального программирования четвертой степени. Однако, как будет следовать из приводимых ниже рассуждений, задача (5) по ряду важных свойств аналогична задаче квадратичного программирования. 2. Базовый случай одной переменной Покажем вначале, что при достаточно больших L решение задачи (5) получается из решения задачи (4) путем целочисленного округления и при этом обеспечивается выполнение требования бинарности переменных xi. Рассмотрим вначале случай функции одной переменной (т. е. N = 1). Тогда левую часть соотношения (5) после деления обеих частей на L можно записать в виде или, после обозначений , , , получаем Проведем замену переменных , т. е. . Тогда, полагая , , из последнего соотношения выводим следующее представление задачи (5): (6) В силу необходимого условия экстремума из (6) следует уравнение Таким образом, точки максимума функции f(z) являются корнями кубического уравнения (7) Уравнение (7) имеет в комплексной плоскости три корня zi(L) (i = 1, 2, 3), которые упорядочим лексикографически, т. е. в порядке возрастания их действительной части; если же действительные части равны, то в порядке возрастания их мнимых частей. Покажем, что корни z1(L) и z3(L) симметричны относительно корня z2(L). Для этого воспользуемся формулами Кардано [6, с. 235]: если дано кубическое уравнение , то (напомним, кубический корень в поле комплексных чисел имеет три возможных значения): (8) причем первый и второй кубические корни в правой части выбираются таким образом, чтобы их произведение равнялось -p/3. В нашем случае , . Кроме того, оба квадратных корня в правой части должны быть одинаковыми. Использование формулы (8) для анализа симметричности корней, ввиду указанных ограничений, неудобно. В частности, необходимо перебрать до девяти пар чисел - все комбинации пар трех значений кубического корня первого слагаемого в правой части (8) и трех - второго слагаемого, поэтому преобразуем выражение (8) в более удобную для анализа форму. Умножим второе слагаемое в правой части (8) на сопряженное выражение - по сути, сопряженным является первое слагаемой в (8). Получим (9) Пусть вначале < 0; отсюда очевидно вытекает p < 0. Тогда . Из (8) следует, что второе слагаемое в правой части равно комплексно сопряженному значению , откуда выводим и, ввиду (9), и Обозначим: - главный аргумент числа w, т. е. и , , где . Для вычисления φ можно использовать следующие соотношения: , если w лежит в первом или втором квадранте, т. е. , и , если w лежит в третьем или в четвертом квадранте, т. е. . Тогда, с учетом (9), корни могут быть записаны в виде (k = 0, -1, 1). Отсюда получаем следующие выражения для корней: 1) k = 0: ; 2) k = -1: ; 3) k = 1: . Так как , то , . Отсюда заключаем: , , и с учетом введенной выше лексиграфической упорядоченности аргументов (т. е. z1 < z2 < z3) выводим: ; ; . (10) Очевидно, что в данном случае все корни - действительные числа. При этом выполняется условие , которое вытекает также из того, что в кубическом уравнении нет квадратичного слагаемого. Случай ≥ 0 рассматривается аналогично. В этом случае w - действительное число, = 0, если w > 0, и = π, если w < 0. Получаем: , , . (11) Нетрудно проверить выполнение равенства . Приведем примеры вычисления корней на основе формул (10) и (11). Пусть вначале p = -1, q = 0, т. е. необходимо найти корни уравнения . Очевидно, z1 = -1, z1 = 0, z1 = 1. Так как , то корни необходимо вычислять по формулам (10). Имеем: , , , откуда выводим: , , и, следовательно,,, , т. е. корни совпали с указанными выше их значениями. Для демонстрации формул (11) рассмотрим случай p = 1, q = 0, т. е. необходимо найти корни уравнения . Очевидно, z1 = -i, z1 = 0, z1 = i. Так как , то корни необходимо вычислять по формулам (11). Имеем: , откуда , , и, следовательно, ввиду (11), , , , т. е. корни совпали с их истинными значениями. Поскольку при из определения a(L) и b(L) следует и , то в силу формулы Кардано для корней кубического уравнения (см. (8)), в которых зависимость корней от коэффициентов кубического уравнения является непрерывной, заключаем: корни уравнения (7) сходятся к корням уравнения z3 - z = 0, корни которого равны -1, 0 и 1, т. е. , , . Отметим, что значения z = -1, z = 0 и z = 1 соответствуют значениям x = 0, x = 0,5 и x = 1 соответственно. Рис. 1. Графики функций Для того чтобы выявить, какой из корней уравнения (7) является точкой максимума, построим графики функций (см. (6)) и параболы - функция f(z) является суммой этих функций. Для определенности предположим, что a(L) > 0. Случаи a(L) < 0 и a(L) = 0 рассматриваются совершенно также. Так как корни функции y2 симметричны относительно оси Оy, то при сложении функций значение суммы в более удаленной от вершины параболы точке z1 будет больше, чем в более близкой к вершине точке z2. Отметим: значение абсциссы вершины не зависит от L. Таким образом, максимум достигается в той из точек z1 и z3, которая находится дальше от вершины параболы. Данный принцип ниже обосновывается на основе оценок. Как видно из рис. 1, при a(L) > 0, точка z3 не может быть точкой максимума. Сравним между собой значения f(z1) и f(z3). Обозначим: и исследуем поведение ψ при . Поскольку при имеем и , , , то ввиду соотношения при выводим: при Отсюда получаем: при . На основе (10) выводим: при , (12) и, аналогично, . На основе (12) имеем: при , и, аналогично, . Отсюда выводим: при . Следовательно, для точки максимума zmax функции f(z) имеем: если > 0, то ; если же < 0, то . При = 0 (т. е. q = 0, p = 1 + a(L)/2 > 0 при достаточно больших L) для функции f(z) справедливо представление , откуда вытекает: , и , т. е. и в качестве точки максимума можно выбирать любую из точек - z1 или z3. Последние соотношения можно получить и на основе соотношений (10) - в данном случае . На основе всего вышесказанного, определения функций a(L) и b(L), полагая , , приходим к следующему правилу выбора решения задачи (6): если A0 + B0 > 0, то ; если A0 + B0 < 0, то ; при A0 + B0 = 0 имеем . Поскольку , , то, полагая (i = 1, 2, 3), приходим к следующему выражению для точек максимума в задаче (6) в терминах переменной x: , (13) где - индикаторная функция, равная 1, когда условие С выполняется, и равная нулю в противном случае. Отметим, что в приведенном правиле полагаем = 0 в случае A0 + B0 = 0, что позволяет упростить исходное выражение целевой функции в случае многих переменных. Отметим, что условие (13) зависит также от параметров метода штрафных функций коэффициентов {Lk, } и введенных дополнительно переменных {tk, }. Общие результаты по методу штрафных функций [7, гл. 5, § 12] позволяют утверждать, что при достаточно больших значениях параметров {Lk, } конечный результат не зависит от конкретных значений этих параметров, хотя, безусловно, целесообразно привести соотношение (13) к виду, который не зависел бы от параметров метода, а зависел бы только от входных параметров задачи (1), (2). 3. Вывод системы уравнений для нахождения решения общей задачи Вернемся к решению общей задачи (1), (2), которая при достаточно больших L равносильна задаче (5). Если точка является точкой максимума задачи (5), то для каждого i, , точка является точкой максимума одномерной задачи по переменной xi, когда значения остальных переменных совпадают со значениями этих же переменных в точке максимума. Последняя задача исследована в разделе 2. Запишем решение полученной одномерной задачи на основе результатов, полученных в разделе 2. Для фиксированного i из (4) получаем: , . Тогда условие оптимальности значения переменной (13) по переменной xi в точке максимума запишется в следующем виде: . Получаем следующую систему уравнений для нахождения оптимального решения задачи (4): , (14) Таким образом, решение исходной задачи (1), (2) квадратичного бинарного программирования свелось к нахождению решения системы (14), состоящей из N уравнений и содержащей N + K переменных {xi. tk }, удовлетворяющих ограничениям (2), tk ≥ 0, и K параметров {Lk}. В частном случае задач бинарного линейного программирования, когда все коэффициенты перед квадратичными слагаемыми равны нулю, система (14) перепишется в виде , (15) Условия (15) при отсутствии ограничений принимают следующий простой вид: , т. е. xi = 1, если Bi > 0, и xi = 0, если Bi < 0, что согласуется с непосредственным решением на основе анализа соответствующей оптимизационной задачи. Метода нахождения решения указанной системы в доступной литературе найти не удалось. Ниже предлагается эвристический алгоритм ее решения, строгое обоснование которого требует дальнейших исследований. 4. Эвристический алгоритм решения задачи Пусть имеется некоторое текущее решение системы (14), удовлетворяющее следующим условиям: для , (16) для , (17) для и для , (18) где I = I(,), , - некоторое подмножество индексов из набора {1, .2, …, N}, для которых условие (14) выполняется, а - множество остальных индексов, для которых условие (14) не выполнено; аналогично определение множества R= R(,). Тогда, если не пусто, то к вектору прибавляем некоторый специальным образом построенный вектор смещения , т. е. получаем . Затем компоненты вектора округляем до ближайшего из чисел 0 или 1; полученный вектор и является новым значением вектора решений . Таким образом, необходимо так выбрать вектор , чтобы в результате описанных построений число элементов множеств и было бы по крайней мере не больше числа элементов множества и . Ниже описывается процедура построения указанного вектора . Попытаемся уменьшить число элементов во множестве на основе следующих эвристических соображений. Обозначим: - выражение, записанное под знаком индикаторной функции в (14). Тогда направление максимального роста функции равно Если (т. е. для индекса i условие (14) выполнено), то желательно так изменять вектор , чтобы для индекса i условие (14) по-прежнему выполнялось. Это означает, что если в (16) xi = 1, т. е. > 0, то желательно, чтобы вектор был бы по направлению как можно ближе к вектору . Если же в (16) xi = 0, т. е. ≤ 0, то вектор должен быть направлен в противоположную сторону, т. е. как можно ближе к вектору -. Выполнение данных условий дает определенные гарантии того, что условие (16) для индекса i сохранится в новом решении . В случае же выполнения (17) для индекса i, т. е. при , требования к вектору в каждом из двух случаев (xi = 0 и xi = 1) заменяются на противоположные по соображениям, аналогичным предыдущему случаю: если в (17) xi = 0, т. е. > 0, то желательно, чтобы вектор был бы направлен в сторону, противоположную направлению роста функции для того, чтобы добиться выполнения противоположного неравенства ≤ 0 - по направлению как можно ближе к вектору -; если же в (17) xi = 1, т. е. ≤ 0, то вектор должен быть направлен в противоположную сторону, т. е. как можно ближе к вектору . Аналогичные рассуждения и обоснования можно использовать и для выбора значений последних K компонентов вектора , относящихся к переменным tk, поскольку необходимо обеспечить выполнение условия tk ≥ 0. Вектор сформируем как линейную комбинацию векторов по всем i. Встает вопрос о величине коэффициента перед каждым слагаемым в этой линейной комбинации. Отметим: можно принять, что вклад слагаемых, для которых выполняется условие (14) (выполнение этих условий надо сохранить), равен вкладу слагаемых, для которых это условие не выполняется, поскольку для решения задачи одинаково важны как сохранение части уже выполняемых условий, так и обеспечение изменения невыполняемых условий. Важно также, чтобы векторы и были равноценны при формировании вектора , поскольку влияние обоих слагаемых одинаково важно при вычислении . Последнее требование мы представим в виде соотношение , где есть длина вектора . Для вычисления коэффициентов указанной линейной комбинации воспользуемся следующим рассуждениями. Если значение достаточно велико по модулю, то значения функции могут быть подвергнуты более резким изменениям с относительно большими шансами, что при этом условие (15) для индекса i сохранится при , условие (17) заменится на противоположное при , и, аналогично, для условий (18). Сказанное означает, что вектор в линейной комбинации должен учитываться в меньшей степени (т. к., как было сказано выше, предполагается вектор максимально «запараллелить» с векторами ), т. е. его коэффициент в линейной комбинации может быть малым, чтобы меньше мешать изменению условий в (17). Если же значение мало, то при формировании вектора в случае (16) (т. е. при ) значение вектора - в составе должно учитываться в большей степени, чтобы обеспечить, по возможности, сохранение условия (16), т. е. коэффициент должен быть большим. В случае же условия (17), т. е. при , при малом абсолютном значении значение коэффициента перед вектором может быть большим, чтобы способствовать изменению условий (17) на противоположные. Отметим: чем больше значение , тем в большей степени должно учитываться значение вектора . Построим векторы и на основе вектора путем замены на нулевые значения последних K координат (относящихся к переменным tk) соответственно первых N координат (относящихся к переменным xi). При этом (N + k)-я координата последнего вектора умножается на tk/2 для всех k от 1 до K. Вышесказанное позволяет предложить следующее выражение для вектора : , (19) где задает знак i-го слагаемого в соответствии с представленными выше обоснованиями для случаев xi = 0 и xi = 1, sgn(tk) =1, если tk ≥ 0, sgn(tk) = -1, если tk < 0, параметры > 0, > 0, > 0 и > 0, > 0, > 0 выбираются так, чтобы и . Можно взять β = 0,1 (на порядок меньше промежутка [0; 1] поиска решений) и δ = 0,01 (большая часть переменных tk равна нулю в оптимальном решении). Тогда линейно выражается через , что после подстановки в выражение для в (19) приводит к сокращению числителя и знаменателя на . Это означает, что вектор от конкретных значений не зависит, и поэтому можно положить = 1. Тогда Аналогичное выражение выводится и для : Необходимо также описать процесс построения начального значения вектора решений . Отметим, что в качестве начального варианта может быть взят любой набор из нулей и единиц длины N, например набор, состоящий из одних нулей. Однако в этом случае в процессе поиска оптимального решения может потребоваться большое число дополнительных итераций, поэтому выбор наиболее приемлемого начального варианта также является важным элементом процедуры решения системы (14). В ряде случаев можно выбрать решение на основе следующей процедуры. Предположим, что целевую функцию исходной задачи (1) удалось представить в виде суммы S отдельных относительно простых слагаемых; например, в задаче формирования учебного расписания [1] целевая функция представляет собой сумму 19 относительно простых слагаемых. Тогда исходная задача (1) запишется в виде (3). По описанной процедуре с нулевыми начальными данными решаем независимо S задач (2), (3), в k-й из которых целевая функция равна k-му квадратичному слагаемому целевой функции. Пусть в результате решения k-й задачи получено оптимальное решение . Тогда в качестве начального решения исходной задачи (1), (2) можно взять следующий вектор: полагаем ; тогда вектор , где , если , и , если . Полученное описанным образом решение во многих случаях будет существенно ближе к искомому оптимальному решению, чем при случайном начальном варианте решения, что значимо уменьшит число требуемых итераций для нахождения оптимального решения . Описанная процедура нахождения решения системы (14) может быть записана в виде следующего общего алгоритма. 1. Задается начальное решение . Можно положить либо взять любой другой набор из N чисел 0 и 1, либо найти начальное значение на основе описанной выше процедуры. 2. Пусть построено решение для некоторого k ≥ 0. Тогда на основе соотношений (17) и (18) находим вектор . 3. Находим вектор : . Затем находим новый вектор решений , где = 1, если , и = 0, если (. 4. Если для i выполнены условия (14), то процесс поиска прекращаем и полагаем . В противном случае переходим к шагу 2. Укажем на следующий важный элемент решения задачи (1), (2). В процессе преобразования указанной задачи в соответствии с методом штрафных функций были введены параметры {Lk} (см. (5)). Выбор указанных параметров метода штрафных функций предполагается выполнить на основе существующих методов, приведенных во многих источниках; например в [7]. Однако специфический тип целевой функции (мы рассматриваем квадратичные функции), по-видимому, позволяет уточнить правила выбора параметров штрафных функций, что предполагается выполнить в дальнейшем. Заключение В ходе исследования рассматривалось обобщение классической задачи квадратичного программирования, когда наряду с линейными ограничениями допускается наличие квадратичных ограничений, в частности случай, когда переменные могут принимать только бинарные значения - 0 или 1. Подобные задачи важны при выборе оптимальных структур систем, состоящих из большого числа вариантов, и выбор или отказ от выбора отдельных вариантов равносилен значениям 1 или 0 соответствующих бинарных переменных. Описана процедура сведения указанной задачи на основе метода штрафных функций к задаче полиномиального программирования четвертой степени, когда все слагаемые, кроме одного, имеющего четвертую степень, имеют степень не более второй. Решение полученной оптимизационной задачи сведено к решению системы уравнений, содержащих только бинарные переменные. Предложен эвристический рекуррентный алгоритм решения полученной системы уравнений, описаны варианты построения начального варианта решения, а также параметров, используемых в предложенном методе. В процессе сведения использованы классические формулы Кардано для корней кубического уравнения, для практического нахождения которых получены соотношения и описана процедура вычисления. Полученный алгоритм может быть использован также при решении многих задач дискретной математики. Возможным продолжением исследований в рамках рассмотренной задачи может стать обоснование приведенного в работе алгоритма решения системы (14), а также поиск и разработка новых, более эффективных методов решения указанной системы, в том числе на основе различных современных эвристических алгоритмов. в частности генетических алгоритмов, алгоритмов муравьиных колоний, пчелиного роя и т. п. Важным направлением является также использование разработанного алгоритма при решении конкретных прикладных задач в различных профессиональных сферах.