Методы и алгоритмы решения задач целочисленного квадратичного программирования на основе линеаризации целевой функции Текст научной статьи по специальности «Математика»

УДК 519.852

МЕТОДЫ И АЛГОРИТМЫ РЕШЕНИЯ ЗАДАЧ ЦЕЛОЧИСЛЕННОГО КВАДРАТИЧНОГО ПРОГРАММИРОВАНИЯ

НА ОСНОВЕ ЛИНЕАРИЗАЦИИ ЦЕЛЕВОЙ ФУНКЦИИ

В. Д. Киселев

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

Ключевые слова: целочисленное квадратичное программирование, теория двойственности, целевая функция.

Задачи целочисленного квадратичного программирования (ЦКП) возникают при проектировании сложных систем, организации вычислений в системах с распределенной обработкой данных, составлении расписаний и др. Основные направления исследований по методам решения задач ЦКП, связаны с изучением следующих возможностей:

- разработка новых алгоритмов решения задач средней и большой размерности;

- создание специализированных алгоритмов решения, учитывающих структурную специфику задач для достижения эффективности вычислений и объемов памяти;

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

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

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

Анализ методов и алгоритмов решения задач ЦКП выявил ряд недостатков:

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

402

- сложность формирования секущих плоскостей в методах отсечения;

- сложность определения коэффициентов линейной функции или существенное увеличение размерности решаемой задачи в методах, основанных на линеаризации целевой функции;

- малая размерность (< 30 переменных для точных алгоритмов) и большое время решения задач.

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

- при методе ветвей и границ - порядок ветвления переменных и оценки границ решения задачи ЦКП на базе двойственной по отношению к линеаризованной задаче;

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

1. Рассмотрим задачу дискретного квадратичного программирования. Необходимо максимизировать [1]

n n n

f(x) = Xcjxj + XX dkjxkxj (1)

j=1 j=1k=1

при ограничениях:

n

Xai jxj <bi, i = 1,2,...,m, (2)

j=1

xj = 0,1,2,..., j = 1,2,...,n, (3)

где D = [dkj ] - неотрицательная матрица порядка n, A = [ai j ] - матрица размерности n x m; C = [c j ] и B = [bi ] - неотрицательные векторы размерности n и m соответственно.

Обозначим через u j ((j = 1, ..., n) максимальное значение целевой

функции следующей задачи при ограничениях (2), (3):

n

uj = max X dkjxk, (4)

k=1

Составим линейную функцию g(x) такую, что

n

n n

f(x) = X CjXj + X j=1 j=1

n

X dkjxk k=1

В(х) = X + uJ)VxJ (5)

j=1

Для любых X = { xj, ] = 1,2,...,п }, удовлетворяющих ограничениям

(2), (3), §(х) > ^х), что подтверждается преобразованием

п п

VxJ £ XcJxJ + X 1 J=1 j=1

Для определения коэффициентов в линейной форме воспользуемся одним из вариантов решения непрерывной задачи линейного программирования (4), (2), в которой условие целочисленности переменных ослабляется и заменяется условием

Xj > 0, j = 1,2,..., п. (6)

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

Двойственной по отношению к (4), (2), (6) является следующая за-

дача:

m

uq = min X biYi, (7)

i=1

m

X ai kyi > dk q 5 k = 1,2,...,n, (8) i=1

yi > 0, i = 1,2,..., m. (9) Алгоритм решения задачи (7)-(9) включает следующие шаги:

10. Определить величину p(r) = —bi—, где k = 1, 2,... - номер ите-

X ai k

kEl(r)

рации; I(r) - множество индексов условий (8), для которых неравенства не выполняются (при r = 1 I(1) ={ k = 1,2,...,n }).

20. Выбрать двойственную переменную y(r), для которой выполняется условие d(r) = min p(r).

i

m+nr-1 (,)

dk q - X Xai kye)

30. Вычислить значение переменной yer) = min-i 1 x 1

k ae k

40. Исключить из множества I(r) индекс уравнения, для которого выполняется условие r = n. Если условие не выполняется, то положить r = r + 1 и перейти к п. 10, в противном случае - к п. 50.

0 Г (r) m

5 . Вычислить yj = Y у. , i = 1,2,...,m и uq = Yь;у;.

t =1 j i=!

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

Таким образом, после решения n задач (7)-(9) определяются величины Uj линейной формы (5). На этом первый этап заканчивается.

На втором этапе при решении задачи ЦКП методом ветвей и границ величина оценки решения определяется как Hsi(xj) = Csi(xj) + Zsi(xj), где

Csl(xl) = Y cjxj + X X dklXkXi - значение целевой функции 1-j e si j k e si

частичного решения задачи (1)-(3); Si ={ (kj), j = 1,2,..., l; 0 < l < n } -

множество индексов переменных задачи (1)-(3), вошедших в l-частичное решение. Для сокращения объема вычислений при определении величины ZSl (xl ) воспользуемся теорией двойственности.

Результаты экспериментальной оценки показали, что наиболее целесообразно использовать способ однократного решения двойственной по отношению к (5), (2), (3) задачи. Для решения двойственной задачи используем рассмотренный выше алгоритм. В этом случае величина Zs1(xi)

при ветвлении всех переменных основной задачи определяется выражени-

m l

ем Zs1(xi) = Xbjyj, где Y = (У1,У2, . ,Ут) - решение двойственной по i=1 j

отношению к (5), (2), (3) задачи, Ь1 = bj - X aijxj, i = 1,2,...,m.

jeSl

Выбор переменной для включения индекса в множество индексов

переменных s1 ведется при условии Hsl(xi) = max Hsi(xk). Для отсечения

k

бесперспективных вариантов используется условие Hs1 < C0, где Со - значение достигнутого рекорда.

Проведена оценка эффективности предложенного алгоритма решения задачи ЦКП, при этом данные для эксперимента генерировались так: коэффициенты cj и pkj (j = 1, ..., m, k = 1, ..., n) - как независимые целые

числа, равномерно распределенные в интервале [0, 100], коэффициенты a^

(i = 1, ... , m, j = 1, ... , n) - как независимые целые числа, равномерно распределенные в интервале [1, 50]; правые части ограничений - как независимые

целые числа в интервале

(i = 1, ..., m). Результаты экспери-

п П

X 2 X ^ >1 ]=1

мента получены по решению 10 задач разной размерности 10(20, 30, 40) х 1(2,3) и показали следующее. Время решения задач каждой размерности невелико и слабо зависит от числа ограничений. Решение задач средней размерности осуществляется за приемлемое время. Отсутствие ярко выраженной зависимости времени решения задач ЦКП от размерности позволяет констатировать, что время во многом определяется характером и структурой исходных данных.

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

П П П

f(x) = XCjXj + XX PkjXkXj , (10)

]=1 j=1k=1

при ограничениях:

n

Xai jxj £ bi, i = 1,2,...,m, (11)

j=1

xj = {0,1}, j = 1,2,...,n, (12)

где P = [pkj ] -неотрицательная матрица порядка n, A = [ai j ] -матрица размерности n x m; C = [cj ] и B = [bi ] - неотрицательные векторы размерности n и m соответственно.

Так как xj = xj для xj = {0,1}, то линейные члены в (10) могут неявно рассматриваться в квадратичной части, поэтому для удобства дальнейших рассуждений линейную часть в (10) опустим. Для j = 1, ..., n обозначим через uj максимальное значение целевой функции следующей задачи:

n

uj = max X pkjxk, j = 1,2,...,n, (13)

k=1

при ограничениях (11)-(12).

Составим линейную функцию g(x) такую, что

n

g(x) = X uj(xj). (14)

j=1

Для любых X = { xj, j = 1,2,...,n }, удовлетворяющих (11)-(12), имеем g(x) > f(x), что подтверждается преобразованием n

.Для определения коэффициентов u j в линейной

n

f (x) = X j=1

X pkjxk k=1

форме (14) воспользуемся решением непрерывной задачи линейного программирования (13), (11), в которой условие целочисленности переменных ослабляется и заменяется условием

0<xj < 1, j = 1,2,...,n. (15)

Применение достаточно простого симплекс-метода для многократного (n раз) решения задач (11), (10), (15) ведет к большим вычислительным трудностям, поэтому предлагается использовать приближенное решение соответствующих двойственных задач (14), для определения которого воспользуемся итерационным алгоритмом, описанным ниже. Двойственной по отношению к (13), (11), (15) при фиксированном значении индекса j = г (г = 1, ..., n) является следующая задача:

, m m+n

uq = min X bjyj + X yi, (16)

i=1 i=m+1

m

X ai kyi + ymk > pk q , k = ^v^n (17)

i=1

yi > 0, i = 1,2,..., m + n. (18)

Таким образом, после решения n задач (16)-(18) определяются величины u j линейной формы (14). На этом первый этап заканчивается.

На втором этапе при решении задачи ЦКП с булевыми переменными методом ветвей и границ, величина оценки решения и порядок ветвления переменных определяется решением двойственной по отношению к (14), (11), (15) задачи, имеющей вид:

Z = min

при ограничениях:

^ m m+n ^

X bi yi + X yi

V i=1 i=m+1 )

(19)

m

X aijУj + Уm+j > ^ , ) = 1,2,...,П, (20)

i=1

yi > 0, i = 1,2,...,m. (21)

Введем следующие обозначения: D = {il,i2,...,ik,...,in} - массив, определяющий порядок ветвления переменных, индекс ^ определяет номер условия (20), выполняющегося на ^ й итерации алгоритма решения задачи (19)-(21); Q = {qei, e = 0,1,...,n, i = 1,2,...,m + n} - матрица, ^я строка которой представляет решение двойственной задачи (19)-(21) при ] = il,i2,...,ik ; S = Sl с So - множество индексов переменных, включенных в s -частичное решение (здесь Sl = {jx j = 1}} So = {jx j = 0}).

Для приближенного решения задачи (19)-(21) и определения порядка ветвления переменных задачи (10)-(12) воспользуемся следующим ите-

рационным алгоритмом, основанным на идеях градиентного метода: 10. Определить величину

_

Z aij

dj _ —, i = 1,2,...,m + n, bm+j _ 1,2,...,n.

где k =1, 2, ... - номер итерации, Ik - множество индексов условий (17), для которых неравенство не выполняется (i1 _ {1,2,...,n}).

20. Выбрать двойственную переменную yk, для которой выполня-

k k ется условие dr _ min di

30. Вычислить значение переменной y и индекс переменной для ветвления на (n - k + 1)-м ярусе дерева ветвлений

m+n k-1 ,

uj- Z Z aijyi

yk _ min-i_1 i_1-,

j arj

где y0 _ 0, i =1, ..., m+k. Индекс iq, определяющий минимум yk, является индексом переменной x-k, для ветвления на (n - k + 1)-м ярусе, запишем как P(k) _ i q .

40. Определить значения элементов k-й строки матрицы Q: qki _ qk-1,i + yk, i = 1,2,...,n + m, где qoj _ 0, i = 1,2,...,m + n.

50. Исключить из множества ik индекс уравнения, для которого выполняется условие (20). Проверить условие k = n, и если оно не выполняется, то положить k = k+1 и перейти к п. 10, в противном случае - к п. 60.

k • k m+n

60. Рассчитать у- _ Zyj, i = 1,2,...,m + n. Z_ Zbjyj + Zyi.

j_1 i_1 i_m+1

Итак, исходя из приведенного алгоритма, при решении задачи (10)-(12) ветвление следует начинать с переменной xin , затем переходить к

xjn 1 и т. д. Тогда верхнюю границу при ветвлении переменной x-k (k = n,

m „ m+n

n -1, ..., 2) можно вычислить как H§(xjk) _ Z Z Pkj + Zb qki + Z qki,

где bW _ bj - Z ajj, i = 1,2,...,m.

k, jeSj i _1 i_m+1

ij 1=10 m

Для отсечения бесперспективных вариантов используется условие Из(х1к) < С0, где С0 - значение достигнутого рекорда. В качестве первого допустимого решения задачи принималось нулевое значение, т. е. С0 =

408

0.

Данные для эксперимента по решению задач ЦКП с булевыми переменными генерировались следующим образом: коэффициенты ^ и р^

( =1, ..., п, к = 1, ..., п) - как независимые целые числа, равномерно распределенные в интервале (0, 100), а коэффициенты а у (1 = 1, ..., т, ] = 1, ..., п) -

как независимые целые числа равномерно распределенные в интервале (1, 50). Правые части ограничений генерировались как независимые целые

числа в интервале

г

n

1, X a

ij

у

(i = 1, ..., m).

j=1

Результаты эксперимента, полученные по решению 10 задач каждой размерности с одним ограничением показали, что время решения для задач при n > 40 предложенным алгоритмом уменьшается в 1,1-3,2 раза по сравнению с существующим. Среднее время решения задачи незначительно зависит от количества ограничений.

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

Для решения задачи (1)-(3) разработан алгоритм на основе метода построения последовательности решений, который применим при выполнении ряда условий:

1) можно найти функцию (миноранту) q(x), такую, что q(x) < f(x) для всех X, удовлетворяющих (2),(3);

2) можно построить алгоритм упорядочения последовательности решений в порядке неубывания q(x). Для i = 1, ..., M обозначим через Uj

минимальное значение целевой функции для следующей задачи:

N

Uj = min X bixü, (22)

l=1

при ограничениях типа (2), (3). Составим такую линейную функцию q(x)

M N

q(x) = XX uibjxij, (23)

i=1j=1

что для любых X = {xij,i = 1,2,..., M, j = 1,2,..., N} удовлетворяющих (2),(3) имеем q(x) < f(x). В качестве коэффициентов Uj функции (23) возьмем Uj = min{bl,l = 1,2,..., N}.

Для построения последовательности решений (23), (2), (3) и упорядочения их по неубыванию q(x) используем метод динамического программирования. Для повышения его эффективности воспользуемся способом встречного решения функциональных уравнений. Возможные назначения задач можно задать сетевым графом с нулевой начальной и N-й конечной вершиной. На основании принципа оптимальности метода динами-

409

ческого программирования составим два функциональных уравнения:

^ Кг-^Ч] )= шш^,] Кь...,^)+ и^ху} 0 £ £ 1 = 1,2,...,Т, 1] е О+,

0 £ £ Б., 1=1;2;...;Т; 1] € О^

(24)

(25)

где О+,О- - множество дуг, входящих и выходящих из >й вершины соответственно.

Функциональные уравнения (24) и (25) отличаются от обычных функциональных уравнений тем, что количество ограничений в них не является постоянной величиной, они могут быть решены при различных значениях 1 = 1, .., Т.

Вычислительный процесс начинается с решения функционального уравнения (24) при 1 = 1, при этом определяются оптимальные последовательности /оп (Б1 п ) и соответствующие зависимости Х1 п (/0 п), п = 1, ...,

N. По значениям зависимостей Х1 п (/о п) последовательно находятся

функции Б}п (/0п), 1 = 2, .., Т. Используя значения последовательностей

/о п (Б п) и Б1 п (/0 п), 1 = 2, .., Т, определяем минимальное значение

/0 п (Б1 п) = Яц, при котором выполняется (2) при 1 = 1 и минимальное

значение /0п^ п) = Я^, при котором выполняется (2) при 1 = 1, ..., Т.

Если Я11 = Я12, то значения переменных х;, п = 1, ..., N соответствующих значению Я11, являются решением задачи (23),(2), (3), в противном случае (Я11 > Я12), переходим ко второй итерации, которая заключается в решении функционального уравнения (25) при 1 = 2.

Функциональные уравнения при 1 = 2, .., Т решим аналогично. На к-й итерации применим дополнительные условия:

шт|/0,п Б0,п,..., ^0,п]+ ^[0 -Б0,п)..., ^-1 -Б0-а1)] }< Я,

п = 1,2,...,N г = 3,5,7,...,К = ш1п{к12,К22,...,Кг-12} при прямом порядке решения (24) и

Ш1п

кл °п^,...,Б & ]+ /0,п [(°1 - °п:ки (б г-1 -Б п:^0)] }< К

п = N,N -1,...,1, г = 2,4,6,...,... при обратном порядке решения (25). Завершив решение задачи (2), (3), (13) методом встречного решения функциональных уравнений динамического программирования, последовательность решений, упорядочим по неубыванию функции (23).

Дальнейшая процедура решения задачи (1)-(3) состоит из следую-

щих шагов:

10. Определяется верхняя /Ц = f(Х0) и нижняя / = g(Xo) границы целевой функции (1). Если ^ = Г, то Х0 - оптимальное решение задачи (1)-(3), иначе Х0 -лучшее текущее решение (1)-(3), где Х0 - оптимальное решение (13), (2), (3); переход к п. 20.

20. Из полученной последовательности решений Х выбирается очередное X;, ; = 1, ..., г. Если g(Xi) > ^ то конец и лучшее текущее решение является оптимальным, в противном случае переход к п 3 .

30. Уточняется значение верхней границы. Если А(Х; ) < ^ ,то ^ = и Х; - лучшее текущее решение (1)-(3). Возврат к п. 20.

Метод встречного решения функциональных уравнений динамического программирования для задач ЦКП оценивался по среднему времени решения задачи и числу задач Р, для которых начальное значение верхней границы совпадает с оптимальным решением квадратичной задачи. Данные задачи формировались как независимые целые случайные числа, равномерно распределенные в интервалах: pj, } = 1, ..., N - [1, 10]; ^к, ) = 1, ...,

К, к = 1, ..., К - [1, 50]; Рдоп. ; и Н;к, ; = 1, ..., М, к = 1, ..., К-

' N Л

1 X ^ . j=1

и

N

1, X^к j=1 .

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

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

Список литературы

1. Денисов Ю.И., Киселев В. Д., Мягков В.Ю., Щербина А.М. Модели и методы решения задач проектирования и испытаний АСУ. М.: Изд-во ВПК, 1997. 249 с.

Киселев Владимир Дмитриевич, д-р техн. наук, проф., kiselevvd@yandex.ru, Россия, Тула, Тульский государственный университет

METHODS AND ALGORITHMS OF THE SOLUTION OF PROBLEMS OF INTEGER SQUARE PROGRAMMIRO- VANIYA ON THE BASIS OF LINEARIZA TION OF CRITERION

FUNCTION

V.D. Kiselev

The effective methods of the decision of the general of tasks and TsKP and problem of TsKP with Boolean variables based by N of linearization of criterion function with use of the theory of a duality are considered.

Key words: integer square programming, theory of a duality, criterion function.

Kiselev Vladimir Dmitrievich, doctor of technical science, professor, kise-levvd@yandex.ru, Russia, Tula, Tula State University

УДК 681.3

ОПРЕДЕЛЕНИЕ КООРДИНАТ ТОЧКИ В СИСТЕМЕ НАБЛЮДЕНИЯ НА БАЗЕ ТУ-МОДУЛЕЙ

А. А. Горшков, А.Ю. Андросов, В.В. Семигорелов, А. А. Аршакян

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

Ключевые слова: TV-модуль, распределенные системы, система координат.

Распределенные информационно-измерительные системы, в которых конечным звеном, потребляющим видеоинформацию, является человек-оператор, выполняющий функции наблюдения сцены и принятия решения о возникновении на ней нештатных ситуаций, в настоящее время широко применяются в различных областях народного хозяйства, в частности в промышленных системах контроля технологических процессов [1, 2, 3].

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