CC BY
19
2018 Вычислительные методы в дискретной математике №41
УДК 519.852
МЕТОД ПОСЛЕДОВАТЕЛЬНОЙ АКТИВАЦИИ ОГРАНИЧЕНИЙ В ЛИНЕЙНОМ ПРОГРАММИРОВАНИИ
В. С. Колосов г. Москва, Россия
Представлен алгоритм решения задачи линейного программирования посредством процедуры последовательной активации ограничений (включения в расчёт одного за другим) с удержанием состояния оптимальности на сгенерированной последовательности вложенных многогранников. Вследствие сжатия области допустимых решений при критерии оптимальности «шах» целевая функция на каждом шаге убывает (движение к максимуму сверху), в противоположность росту в других методах (снизу). Компьютерные эксперименты демонстрируют преимущества программной реализации этого алгоритма перед опцией симплекс-метода программы linprog библиотеки MATLAB в скорости и полноте выводимой информации.
Ключевые слова: линейное программирование, активация ограничения, MATLAB, компьютерный эксперимент.
DOI 10.17223/20710410/41/11
METHOD FOR SEQUENTIAL ACTIVATION OF LIMITATIONS IN LINEAR PROGRAMMING
V. S. Kolosov
Moscow, Russia E-mail: [email protected]
In this paper, we present an algorithm for solving a problem in linear programming by means of procedure of sequential activation of limitations (inclusions in calculation are made one after another) with retention of the optimality status in the generated sequence of nested polyhedrons. Due to compression of the admissible solutions area, the objective function decreases on each step at the "max" criterion of an optimality (the movement to a maximum "from above") contrary to its growth in other methods ("from below"). Geometrically, the motion to the maximum starts from the trivially determined starting point and continues along the broken line outside of the polytope of admissible solutions. In the theoretical justification of the algorithm, the signs for the incompatibility of the system of conditions in the problem, for the uniqueness and nonuniqueness of the solution, and for its unboundedness are formulated and proved. Computer experiments demonstrate the advantages of the program implementation of this algorithm in speed and completeness of the output information over the simplex-method option of the MATLAB's library linprog program.
Keywords: linear programming, activation of limitation, MATLAB, computer experiment.
Введение
Число существенно различных методов условной линейной оптимизации невелико. Еще в 1939 г. Л. В. Канторович сформулировал ряд задач линейного программирования и для их решения предложил метод разрешающих множителей [1]. В 1947 г. Дж. Данцигом был создан симплекс-метод, нашедший отражение в отечественной литературе в изложении автора в переводе [2]. Разработанные в то время методы последовательного улучшения плана, сокращения оценок и сокращения невязок классифицируются по типу решаемых ими задач — прямой, двойственной или сразу обеих [3]. В середине 1960-х годов И. И. Дикин в [4] предложил решение задачи линейного программирования методом внутренних точек. Поиск новых подходов к оптимизации дал эффектный, но вычислительно неэффективный метод эллипсоидов Л. Г. Хачияна [5], принципиально важным достоинством которого явилось выражение времени счёта полиномом от размерности задачи. Полиномиальность варианта метода внутренних точек Н. Кармаркара от 1984 г. в совокупности с его вычислительной эффективностью послужила причиной его патентования в США. Следствием изысканий в данном направлении было множество публикаций, в их числе [6-8]. Метод внутренних точек глубоко проработан и как эффективный оптимизатор задействован в ряде программных продуктов, в том числе в МАТЬАБ.
К недавним разработкам относятся метод оператора-проектора [9] и скелетный алгоритм [10]. Суть первого заключается в проектировании произвольной точки на гиперплоскости ограничений вплоть до попадания в область допустимых решений и удержании последующих перемещений в заданных границах согласно критерию оптимальности. Скелетный алгоритм не требует обращения матриц и подменяет исходную задачу рядом вспомогательных задач пониженной размерности. В процессе тестовых расчётов установлено быстрое продвижение к оптимуму в начале вычислительного процесса и его резкое замедление по мере приближения к цели. Метод внутренних точек является приближённым и находит решение со сколь угодно малой заданной погрешностью, но полиномиален; ни один из известных точных алгоритмов таким важным свойством не обладает, за некоторыми исключениями. Так, при условии неотрицательности всех параметров задачи и совместности системы ограничений метод экспоненциальной аппроксимации [11] обещает находить точный оптимум в задаче с п переменными и т ^ п равенствами всего за (п — т) шагов. Алгоритм выявляет и исключает переменные с гарантированно нулевыми оптимальными значениями, и решение исходной задачи завершается решением системы уравнений с квадратной матрицей т-го порядка. Перспективы ослабления жёстких требований этого метода к исходным данным не просматриваются.
Перечисленные выше методы объединяют два общих свойства:
— многогранник допустимых решений подразумевается существующим в своем целостном виде, т. е. предполагается одновременная действенность всех ограничений системы условий;
— на каждом шаге поиск следующей точки производится с умыслом улучшения значения целевой функции в пределах области допустимых решений.
Недавние публикации свидетельствуют, что, несмотря на снижение интереса к проблематике линейного программирования, поиск новых путей в этой сфере продолжается. В [12] автор завершил разработку точного алгоритма решения задачи линейного программирования с использованием параметризации лимитов системы условий, запрограммировал его на языке Еог1гап-4 и успешно апробировал на ЕС-ЭВМ. С учётом
современных возможностей программных и аппаратных средств данный метод им основательно доработан и реализован средствами MATLAB.
Далее изложена доказательная база этого по сути геометрического метода решения общей задачи линейного программирования, отличительная черта которого состоит в следующем. Многогранник допустимых решений не существует изначально, а формируется путём последовательного включения в расчёт (активации с помощью параметризации лимитов) ограничений с удержанием оптимальности целевой функции на генерируемой последовательности вложенных многогранников. Начальным является гиперпараллелограмм, задаваемый ограничениями на переменные с возможно неограниченными лимитами, на котором оптимум определяется элементарно. Затем в его систему условий добавляется первое нетривиальное ограничение, быстрым пересчётом находится оптимум на построенном многограннике (вложенном в предыдущий) и т. д. Активируя неравенства, алгоритм последовательно усекает образующиеся многогранники до окончательной конфигурации, постепенно сужая область допустимых решений, вычисляя и фиксируя на каждой итерации алгебраические и геометрические характеристики текущего оптимума для последующих применений.
Метод оперирует с неравенствами и, исключая переменные через посредство равенств, понижает размерность задачи, но может работать напрямую и с равенствами, легко разбираясь с ситуациями вырожденности и маскировки равенств совокупностью неравенств. В результате сжатия области допустимых решений в задачах с критерием оптимальности «шах» целевая функция на каждом шаге убывает или, по крайней мере, не возрастает (к максимуму сверху), в противоположность её традиционному увеличению в других методах (снизу). Геометрически происходит перемещение в оптимум вдоль некоторой ломаной прямой из тривиально определяемой точки, которая для многогранника допустимых решений может оказаться как граничной, так и внешней. В процессе активации ограничений выявляется несовместность системы условий с указанием виновных в этом ограничений, вычисляются координаты оптимальной вершины и смежных с ней в случае неединственности, а также определяются двойственные переменные как оценки чувствительности.
Разработанный алгоритм реализован автором в виде программы mpao в среде MATLAB и апробирован на обширном экспериментальном материале. Показана простота процедуры ввода исходных данных и запуска программы. Установлено преимущество по быстродействию и широте выводимой информации программы mpao по сравнению с опцией симплекс-метода программы linprog библиотеки MATLAB на примерах задач с размерностями до 1000 х 1000.
1. Теоретическое обоснование алгоритма
В евклидовом пространстве En решается задача линейного программирования в самой общей постановке с разделением ограничений на «многофакторные», содержащие несколько переменных, и «координатные» (однофакторные), устанавливающие диапазоны возможных изменений собственно переменных при допущении отрицательности и неограниченности всех фигурирующих здесь величин включая сами переменные:
L(x) = cx ^ шах, c,x Е En,
— ж ^ b(i) ^ aix ^ ^ а^=(а(1)... a(n)),i=1,... ,m (многофакторные), (1)
— ж ^ ^ x(j) ^ ^ j = 1,..., n (координатные).
Здесь сам вектор нумеруется нижним индексом, а его координаты — верхним в скобках. Многофакторные и координатные ограничения, высекая в En многогранник до-
пустимых решений (МДР) M, геометрически равнозначны. Первым отвечают гиперплоскости с задаваемой нормалями а ориентацией в пространстве, а вторым — ортогональные осям координат. Равенства реализуются уравниванием значений нижней и верхней границ разрешённых изменений и в отдельную группу не выделяются.
Равенствам отвечает размещение МДР в пересечении соответствующих гиперплоскостей, что открывает принципиальную возможность понижения размерности задачи за счёт фактического исключения (в отличие от симплекс-метода) части переменных. Это обстоятельство порождает две алгоритмические ветви. Следуя первой ветви, из равенства a(1)x(1) + ... + a(n)x(n) = b(i) выражается переменная с ненулевым коэффициентом; пусть это x(n) = ai^x(1) + ... + ain 1)x(n-1) + ain). При её исключении из остальных многофакторных ограничений и целевой функции координатное ограничение d(n) ^ x(n) ^ d(n) преобразуется в многофакторное неравенство ^ a(1)x(1) + ... + a(n 1)x(n-1) ^ с меньшим числом переменных. В итоге
размерность пространства решений понижается, но только за счет переменных при сохранении количества многофакторных ограничений. Исключением равенств из системы условий заканчивается подготовительный шаг алгоритма. С этого момента МДР считается полноразмерным (dimM = n), т.е. задаётся только неравенствами. Особенности другой ветви алгоритма без исключения равенств из системы условий будут разъяснены позже.
Совокупности неравенств в En отвечает выпуклый многогранник M. Её любая подсистема ранга n представляет собой неограниченный многогранник M(x) D M с единственной вершиной в точке x, являющейся решением соответствующей системы линейных равенств; пусть это {а^ж = b(i) : i = 1,... ,n} .В параметрическом представлении
М(ж) = {ж = x + £ в.?V(x), ^ 0, j = 1,..., n), (2)
I j=1 J
где Vj (x) — направляющие векторы рёбер Rj (x) = {x = X + tVj (x), t ^ 0} многогранника M, выходящих из вершины xx и образованных пересечением n — 1 гиперплоско-
n
стей GGi — ^ x : а^ж — Тогда Rj (xx) = P| Gi и {xx} = Rj (xx) П Gj. Алгебраически
i=1,i=j
Rj (xx) = {x : aix = b(i), i = 1,... , n, i = j; ajx ^ b(j)}, где оно индексируется номером того единственного ограничения, которое в его образовании не участвует. Возьмём произвольную точку xj G Rj (xx) и построим вектор Vj (xx) = xj — xx, который, будучи направляющим вектором, удовлетворяет однородной системе ранга n — 1:
|aiVj(xx) = aixj — aixx = b(i) — b(i) = 0, i = 1,..., j — 1, j + 1,..., n.
Этот вектор ортогонален нормалям n— 1 гиперплоскостей и может быть найден с помощью алгоритма Грама — Шмидта, который линейно независимую совокупность векторов u1,..., un-1 G En рекуррентными соотношениями ортогонального проектирования
u°=Mi, uk=uk-1— (uk-1uJk-1) uk-1, ui-1=ui-1/ I ui-1 ll, k < i < n (3)
(k — «порядок» проекции) преобразует в ортонормированную систему {Ui-1, i = 1,... , n — 1 : Ui-1 uj-1 = }, — символ Кронекера. Очевидно, что при линейной независимости {uo, u1,... , un-1} решением системы {uiv = 0: i = 1,... , n — 1} будет вектор
n- 1 ( )
v (u0,u1,... ,un-1) = —u0 + (u0ui-1) ui-1, где векторы ui-1 вычислены проектором
i=1
Грама — Шмидта (3). Легко проверяются такие свойства построенного вектора, как
ио'У < 0 и = 0, г € {1,... ,п — 1}. Поскольку вершина ж образована пересечением п гиперплоскостей, то, принимая поочередно в качестве и0 их нормали, можно найти все направляющие векторы. В такой схеме для каждого вектора строится своя матрица ортогонализации. Поэтому более эффективно работает приём с одноразовым обращением матрицы координат нормалей и принятием в качестве направляющих векторов её столбцов, взятых с обратным знаком. С ростом размерности задачи метод ортогона-лизации быстро уступает обращению матрицы, которое представляет собой механизм вычисления направляющих векторов, выходящих из вершины многогранника.
Рассматриваемый метод опирается на общеизвестный факт достижения линейной функцией экстремума на выпуклом замкнутом многограннике по крайней мере в одной из его вершин. Метод реализует идею формирования из имеющихся ограничений начального многогранника М0, с гарантией содержащего оптимум задачи х*, нахождения на нём оптимальной вершины ж0 и построения последовательности оптимальных вершин ж0 ^ х1 ^ ... ^ х*т = х * на промежуточных вложенных многогранниках М0 1Э М1 1Э ... Мт = М. В этой цепочке текущий Мк получен усечением предыдущего Мк-1 следующим к-м ограничением.
В качестве М0 разумно принять ортогональный параллелепипед, образованный координатными ограничениями. В случае неопределённости какой-либо границы её значение сообразно конкретной ситуации полагается равным Если М = 0, то,
безусловно, х * € М0. В силу структурной примитивности М0 максимум на нём линейной функции Ь достигается в вершине х0 = ,... , х0п) ^, координаты которой
ОГ / ' если с(О) ^ 0, • определяются простым правилом: х0 = < ') ') ] = 1,... , п.
I с/ ' , если с << 0,
Определение 1. Полупространство (неравенство) Р = {х : ах ^ Ь} назовём значимым для точки х € Еп, если х € Р (ах > Ь), и незначимым в противном случае (ах ^ Ь). Значимость и незначимость гиперплоскости (равенства) О = {х : ах = Ь} определяются через неравенство ах = Ь и равенство ах = Ь.
В результате усечения М0 первым многофакторным двусторонним неравенством образуется многогранник М1 = М0 П Р1 П Р1, где Р1, Р1 —встречные полупространства с непустым пересечением: Р1 = {х : Ь(1) ^ а1 х}, Р1 = {х : а1х ^ Ь }. Поскольку М1 С М0, то тахЬ(х) ^ тахЬ(х). Характер последующих действий определяется
х€М 1 хЕЫо
значимостью данного ограничения для х0, т. е. расположением величины а1х0 относительно отрезка , Ь ]. Неравенство Ь(1) ^ а1х0 ^ Ь означает, что х0 € М1, и потому максимум Ь на М1 достигается в той же самой вершине, т.е. хЦ = х0. Таким образом, данное ограничение для х0 незначимо и можно перейти к построению многогранника М2, усекая М1 следующим ограничением. В противном случае проверяется М1 = 0, вычисляется х1 и т. д.
Рассмотрим к-й шаг алгоритма, к началу которого уже построен многогранник М*-1 = М0 П Р1 П Р1 П ... П Рк-1 П Рк-1; пусть Мл-1 = 0. Если х^-1 € Мк, то хк хк-
_1 и текущий шаг на этом завершается. Остаётся разобраться с ситуацией
(к) _(к)
х*к_ 1 € Мк. Здесь двустороннее неравенство Ь( ) ^ акх ^ Ь для вершины хк_ 1 может
(к) _(к)
быть не выполнено лишь с одной стороны и либо акх*к_ 1 < Ь( ), либо акх к_ 1 > Ь . Принципиального различия между типами этих двух неравенств нет, так как один сводится к другому умножением неравенства на —1.
производятся основные рас-
В случае значимости ограничения ак1 ф Ь(к),Ь( )
чёты. Пусть для определённости ф Рк, граница которого (к = {ж : акж = Ь(к)}. Параметризуем статическое полупространство Рк = {ж : акх ^ Ь( ^, образовав динамическую структуру Рк(Л) = {ж : акж ^ Л}, Рк(Ь(к)) = Рк, (5к(Л) = {ж : акж = Л}, ¿к (б(к)) = (к.
—(к)
Определение 2. Активацией значимого ограничения ак ж ^ Ь называется изменение параметра Л от до Ь ), в результате чего движение гиперплоскости ((Л) влечёт трансформацию жк-1 ^ жк.
Активация ограничения — динамическая процедура его приведения в активное состояние в случае значимости. Под активностью ограничения Р = {ж : аж ^ Ь} в точке ж традиционно понимается её принадлежность граничной гиперплоскости, т. е. ж Ф 5 = {ж : аж = Ь} или аж = Ь. Вычисление ж* составляет содержание «нулевого» шага алгоритма, на котором активированы только все координатные ограничения.
Обозначим: {ж^- : ] = 1,...,г} — концевые точки рёбер Я (жк_ 1) с направляющими векторами г^ (жк_ 1), являющиеся вершинами Мк-1, смежными с ж*_ 1; Рк(Л) = = {ж : акж ^ Л} — активируемое динамическое полупространство; Л* = акжк_ 1; Лj = = ак:Xj-; Мк(Л) = Мк-1 П /Рк(Л). _
Сигналом к активации многофакторного ограничения Рк является его значимость,
и на к-м шаге этому отвечает решение задачи оптимизации Ь(жк(Л)) = тах Ь(ж),
хемк(л)
~ ~ —(к) ~ при том что Мк(Л) = Мк-1 ПРк(Л), = Л ^ Ь . По построению Мк(Л) при всех Л из
указанного диапазона относительно Мк-1 не расширяется и Мк (Л) С Мк-1, Мк (+^) = = Мк-1, Мк (Ь(к)) = Мк, ж*(Ь(к)) = жк. Полупространство Рк(Л) при Л ^ Л* = акжк-1 для жк_ 1 является незначимым. Его граница (5к (Л) при Л = Л* проходит через ж*_ 1, а при Л < Л* отвечает за трансформацию жк_ 1 ^ жк (Л).
Теорема 1. Если акгj• (жк-1) ^ 0 для всех ] ф {1,...,г}, то для всех Л < Л* система условий задачи (1) несовместна (М = 0).
Доказательство. Для любого Л ^ Л* по построению Мк(Л) С Мк-1, и в соответствии с представлением (2) для любого Р Ф Мк (Л) с учётом в ^ 0 для всех
г _
Ф {1,..., г} находим ак ж = ак жк-1 + ^ в ак гj (жк-1) ^ ак ж*-1 = Л*. Тогда ж ф Мк (Л) при Л < Л*, т. е. Мк(Л) = 0, а вследствие М С Мк(Л) и М = 0. ■
Замечание 1. Здесь сформирован механизм выявления противоречивости системы условий задачи в процессе активации ограничений.
Теорема 2. Пусть ж — вершина выпуклого многогранника М = 0 с направляющими векторами всех выходящих из неё рёбер г1(;г),... , гг(ж). Максимум Ь = еж на М достигается в ж, если егj (ж) ^ 0 для всех ]; эта вершина единственна в случае cгj (ж) < 0 для всех ]. Верно и обратное.
Доказательство. Достаточность следует из (2). Из включения М С М(ж) по-
г
лучаем, что еж = еже + ^ в егj (ж) для всех ж Ф М. Учитывая существование хотя
j=l
бы одного вк > 0 при ж = ж в случае егк (ж) ^ 0, получаем еж ^ еже, а при строгом неравенстве егк (ж) < 0 таким же строгим будет неравенство еж < еже.
Необходимость докажем от противного. Пусть сх > сх для всех х € М, но существует гк(х), для которого сгк(х) > 0. Так как х(£0) = х + ¿0гк(х) € М при некотором малом ¿0 > 0, то сх(£0) = сх + ¿0сгк(х) > сх, что противоречит оптимальности х. Далее, предположив единственность оптимума и существование гк(х), для которого сгк(х) = 0, получаем противоречие со сделанным предположением, поскольку сх(£0) = сх при х(£0) = х. ■
Замечание 2. Так фиксируется неединственность решения.
Теорема 3. Если Мк = 0, то для любого Л € {Л* — е, Л*}, где Л* = акхк-1, при достаточно малом е > 0 имеет место акхк(Л) = Л, чему геометрически соответствует хк(Л) € Ок(Л).
Доказательство. Поскольку х к (Л) € Мк (Л), то ак х к (Л) ^ Л. При допущении противного (х к (Л) € О к (Л)) это неравенство может быть только строгим (ак х * (Л) < Л). Из Мк(Л) С Мк-1 следует хк(Л) € Мк-1 и в силу выпуклости Мк-1 Э [хк(Л),хк_ 1]. По условию Л < Л * и, следовательно, концевые точки отрезка х* (Л) и хк_ 1 лежат строго по разные стороны гиперплоскости Ок(Л), а значит, существует хл € (хк(Л),хк-1), что хл € О к (Л). В параметрическом представлении (хк(Л),х к_ 1) = {х(£) = х * (Л) + ¿(х * _1 —
— хк(Л)) : £ € (0,1)}; тогда существует ¿л € (0,1), при котором хл = х(£л) = х к(Л) +
-1—х к(Л)).Отсюда находим хк_ 1 =х к(Л)+¿0 (хл—х * (Л)) ¿0 = 1/^л >1.Сделанное допущение х к(Л) € О к (Л) влечёт строгое неравенство сх * (Л) > сх для всех х € О к (Л), в том числе схк(Л) > схл или с(хл — х к(Л)) < 0. Посчитав схк-1 = с(хк(Л) + ¿0(хл —
— х к (Л))) = сх к (Л) + ¿0с(хл — х * (Л)) < сх к (Л), приходим к противоречию с оптимальностью хк_ 1 на Мк-1 э х к(Л). ■
Замечание 3. Таким образом, если Мк = 0, то либо сохраняется оптимальность предыдущей вершины и х к = х*_ 1, либо оптимальная вершина удовлетворяет активируемому ограничению как равенству с одной из его границ, что означает её принадлежность гиперплоскости (к = {х : акх = Ь(к)} или Ок = {х : акх = Ь )}.
Теорема 4. Если многогранник Мк = 0 и активируемое ограничение Рк = {х :
—(к) г —(к)
акх ^ Ь } значимо для хк-1, т. е. 7 = : акгО- (хк-1) < 0} = 0 и акхк-1 > Ь , то в малой окрестности [/£(Л *) = {Л :0 <Л * — Л<е} вершина х к (Л), доставляющая целевой функции Ь(х) максимум на Мк(Л), может быть представлена в виде
ЛЛ
х*(Л) = хк-1 + (хЛ )г'о(хк-1), Л € и£(Л*), (4)
где индекс определяется решающим правилом
сг'о(х к-1) • сго(х к-1)
./\ К— 1/ /г\
— min ----. (5)
«k.(xk-i) .'eJ ßkv.(4-i)
Доказательство. В малой окрестности xk_ 1 при Л ^ Л * = ßkxk_ 1 вершинами Mk(Л) являются точки пересечения рёбер R. (xk_ 1) многогранника Mk-1 с гиперплоскостью Gk(Л) = {x : ßkx = Л}, и только они. Поэтому в согласии с теоремой 3 остаётся выбрать ту из них, в которой L достигает максимума. По построению ЛЛ
xk.(Л) = xk-i +--—гv.(xk-i) является вершиной Mk(Л), поскольку она лежит на
ßk(xk-i)
ребре с направляющим вектором v. (xk_ 1) и в то же время оказывается на гиперплос-
кости (5 к (Л), что нетрудно проверить прямым расчётом:
/Л\ ( * /А А^Ч Ч ? (жк_* Л Л>|! Л
акжк,? (Л) = аП жк_ 1 + (Л - Л )-——- = акж^ 1 + Л - Л = Л.
V ак г? (жк-1)) к 1
Целевая функция в этих вершинах принимает следующие значения
/л\ * , /л л*\ ег?'(жк-1) • _ Т
еxk,j(Л) = еж*-1 + (Л - Л ) , * ч, 3 ф ^
ак Ч? (жк-1)
Здесь второе слагаемое неположительно, так как, согласно теореме 2, ег?(жк_ 1) ^ 0 и по условию Л — Л* ^ 0, акг?(жк_ 1) < 0. Поэтому для всех Л Ф ие(Л*) максимум Ь достигается на ребре с направляющим вектором г?0 (жк_ 1), для которого величина отношения (ег?(жк_ 1))/(акг?(жк_ 1)) минимальна, т.е. жк,?0(Л) = жк(Л). ■
Замечание 4. Полученный результат позволяет придать двойственным оценкам следующий смысл. Найдём двойственную оценку активируемого ограничения в вершине жк(Л) прямым расчётом с учётом её аналитического представления (4):
= ±с + Л — Л* г (ж* Л = ег?0 (жк-1) (6)
<*Л ^Л4ж*-1 + акг?о (жк-1)?(жк-1^ акгл(жк-1). (6)
Таким образом, двойственная оценка активируемого ограничения по сути является решающим правилом (4), и наоборот, в чём и проявляется связь прямой и двойственной задач. Теперь формулы (4), (5) можно интерпретировать так: на текущей стадии активации значимого ограничения оптимальными будут точки ребра, вычисляющего двойственную оценку этого ограничения [13]. Правило (5) позволяет после нахождения оптимума вычислять двойственные переменные посредством последовательной псевдоактивации (назначения в качестве активируемого) каждого базового ограничения.
Далее для удобства принята сквозная нумерация ограничений: сначала координатные (1,..., п), а затем многофакторные (п + 1,..., п + т).
Определение 3. Назовём базой вершины ж совокупность номеров всех активных ограничений В (ж) = {¿к : а^ж = Ь^, к = 1,... , г ^ п}, т. е. тех, которым она удовлетворяет как равенствам (принадлежит гиперплоскостям), а сами ограничения (гиперплоскости) — базовыми. Величина д = г — п +1 называется кратностью вершины. Однократная вершина — простая.
Активация ограничения в случае его значимости наряду с изменением оптимума влечёт соответствующую трансформацию базы. Локальный результат теоремы 4 допускает следующее расширение.
Теорема 5. База вершины жк(Л) = ж?(Л) сохраняется в пределах интервала (Л?0, Л*), левая граница которого отвечает концевой точке жж?0 ребра Я?0(жк-1), т.е. Л?0 = акжж?0. Вершина ж?0, как смежная с жк_ 1, вычисляется из расчёта максимального сдвига вдоль вектора г?0 (жк_ 1), допускаемого остальными ограничениями многогранника Мк-1, помимо базовых, а также противоположной гранью 3о-го базового условия.
Доказательство. Очевидно, что кандидатами на пересечение с ребром Я?0 (жк_ 1) являются все активированные к настоящему моменту небазовые ограничения. В дополнение к этому необходимо учитывать также гарантированное пересечение данного ребра с противоположной гранью базового ограничения с номером 30, поскольку из всех базовых гиперплоскостей направляющий вектор г?0 (жк_ 1) не ортогонален только
ей. Сообразно этим обстоятельствам параметры принятого для кандидата обозначения а«х = 6(5) конкретизируются по следующей схеме:
. (хк-1)>0 ^ а=е«, 6(в) .(хк-1)<0^ав=-бв, &(в)=-^(в),
5 = 1) \Ь } ,
для многофакторных ограничений агV.(хк_ 1)>0 ^ а«=аг, 6(8)=6(г),
аг. (4- 1)<0:
=-0г, 6(5)=-Ь(г),
г = 1,..., к, 5 = п + г € В(хк_ 1) \ {п+^0}
Здесь е« — единичный базисный вектор: е() = . Первая группа включает все небазовые координатные ограничения, а вторая — активированные к настоящему моменту небазовые многофакторные ограничения. Базовое ограничение с номером учитывается в той группе, к которой оно принадлежит.
Величина сдвига вдоль ребра Я. (ж^-1) до пересечения с гиперплоскостью С« = = {х : а«х = 6(а)} в точке х.0,« определяется соотношениями
х
.70,«
X
к-1
+ (хк-1), (хк-1 + . (хк-1^ = Ь ),
Ь(«) -
а«хк- 1
70,«
а«г>70(хк_ 1)
Концевой точке ребра хх70, сопряжённой с х£_ 1 по ребру Я. (хк_ 1), отвечает максимально возможный сдвиг в пределах многогранника Мк-1:
¿70,«0 = ^.0,«, ^ ={1,...,п +к —1} \ в(хк-1),х70 = + .,«0.(хк-1)
и А. = акхх70 —левый конец интервала стабильности базы В(х^(А)). ■
Таким образом, х£(А) —оптимум X на Мк(А) для всех А € [А., А*] и характер последующих действий определяется одной из двух ситуаций:
Если Ь(к) ^ А70, то хк = хк (Ь(к))
х
к-1
+
Ь(к) - А*
( * ч.(хк-1). Преобразование ак. (хк-1)
хк-1 ^ хк влечёт трансформацию базы В(хк-1) ^ В(хк) = (В(хк-1) \ Оо}) и и {к}, состоящую в замещении элемента на номер к, чем и завершается текущий шаг алгоритма.
2) В случае Ь( ) < А. активация ограничения продолжается относительно хх70 с базой В(х70) = (В(хк-1) \ Оо}) и {з0} вплоть до занятия гиперплоскостью 6к(А) её геометрически окончательного положения.
На этом заканчивается построение вычислительных конструкций метода последовательной активации ограничений (МПАО), название которого происходит от характера производимых действий. В процессе активации ограничений теорема 1 выявляет несовместность системы условий задачи. Формулы теорем 4 и 5 вычисляют координаты оптимума и двойственные переменные. Неограниченность некоторых координат х* означает неограниченность решения задачи по соответствующим переменным. В случае неединственности оптимума в интересах послеоптимизационного анализа теорема 2 позволяет перечислить смежные вершины в виде множества
X * (М) = {х.
т х т,.
х
+ г.V. € М : сь.(хт) = 0, ] = 1,... , п*} .
Каждая итерация заканчивается замещением в базе одного элемента другим, чему соответствует замена только одной строки матрицы нормалей активных ограничений.
Поэтому при нахождении направляющих векторов посредством инвертирования матриц использование факторизации существенно повышает быстродействие алгоритма.
Следуя первой ветви алгоритма с исключением переменных через равенства, по окончании решения задачи оптимальные значения исключенных переменных восстанавливаются по найденным на предварительном шаге их выражениям через действующие переменные.
Отличие ветви алгоритма без исключения равенств из системы условий состоит в том, что активация многофакторных ограничений производится в следующем порядке: сначала равенства как безусловно активные, затем неравенства. Порядок обращаемых матриц равен размерности пространства решений и, естественно, больше по сравнению с первой ветвью алгоритма. Зато объём других вычислений несколько уменьшается, поскольку по мере поступления равенств в базу промежуточной оптимальной вершины при активации ограничения в расчёт принимаются только те рёбра, которые удерживают последующие перемещения в гиперплоскостях равенств.
Подозрение на возможность зацикливания вычислений в МПАО возникает лишь относительно ситуации прохождения кратных вершин (определение 3), которая в результате накопления ошибок округления может как разрешиться сама собой, так и ложно возникнуть. В [3] в борьбе с зацикливанием предлагается делать ставку на регулярность метода. Следуя этой рекомендации, в случае неединственности оптимального ребра неоднозначность выбора в (5) разрешается по принципу наименьшего номера.
Процедура активации ограничения обладает регулярностью в том смысле, что в соответствии с теоремой 4 движение допускается вдоль направляющего вектора, образующего с нормальным вектором активируемого ограничения исключительно тупой угол. Это обстоятельство обеспечивает безаварийное прохождение вершин кратности больше 1, хотя МПАО настроен на простые вершины. Убедиться в этом помогает умозрительный пример 2-кратной вершины с применением к нему приёма «расклеивания» [14], суть которого состоит в представлении такой вершины в виде двух близко расположенных простых вершин за счёт малого изменения свободного члена со смещением «лишней» гиперплоскости вдоль её нормали во избежание сужения множества допустимых решений. В результате алгоритм если все-таки и пройдёт по фиктивному ребру, то всё равно выйдет на ненулевое оптимальное ребро. Однако тот же результат достижим при виртуальном расклеивании кратной вершины. Величина смещения вдоль фиктивного ребра равна 0, но, несмотря на сохранение координат промежуточной оптимальной вершины, произойдёт преобразование её базы, что означает фактический переход в структурно другую вершину. Повторение элементного состава базы означало бы возврат в ранее пройденную вершину, что возможно лишь при движении под запрещённым острым углом в обратном направлении. Это свидетельствует о безразличии МПАО к кратности вершин и невозможности зацикливания по этой причине. МПАО нечувствителен к вырожденности вершин в традиционной терминологии, легко распознаёт и разрешает случаи дублирования равенств и их маскировки совокупностью неравенств, что подтверждается экспериментальными расчётами.
В отличие от метода внутренней точки (МВТ), являющегося приближённым, МПАО — точный, подобно симплекс-методу (СМ). Однако теоретическая точность метода ещё не признак его безусловного превосходства, так как злую шутку может сыграть накопленная ошибка округления. В то же время приближённый метод может оказаться и точнее, и быстрее, что подтверждается практикой применения МВТ. Вычислительная эффективность МПАО определяется количеством шагов при активации
многофакторных ограничений. Её примитивная оценка сверху на базе подсчёта числа проходимых алгоритмом вершин промежуточных многогранников является ожидаемо экспоненциальной, и возможность её снижения вызывает законное сомнение.
2. Анализ результатов вычислительных экспериментов
МПАО реализован в MATLAB 7 R2011b программой mpao, которая решает задачу линейного программирования с двусторонними ограничениями, допуская неопределённость границ и отрицательность переменных. Условия-равенства оформляются равенством верхней и нижней границ допустимых изменений, в отдельную группу не выделяются и могут перемежаться неравенствами. Исходные данные передаются в программу массивом определённой структуры или текстовым файлом trafaret с вводом с клавиатуры только его номера при запуске программы. Данные из трафаретного файла считываются программой как строчные, что обеспечивает отсутствие жёстких требований к формату их расположения в поле ввода и визуализацию исходных данных сообразно вкусам пользователя. В случаях несовместности, неограниченности, неединственности решения выдаются соответствующие сообщения. Всплывающие меню с активными кнопками обеспечивают управление функциональными возможностями программы в части её запуска, визуализации исходных данных и полученных результатов, вычисления альтернативных решений в случае неединственности.
Рассмотрим результаты применения программ mpao и linprog библиотеки MATLAB в опциях крупномасштабного алгоритма LSA (МВТ) и среднемасштабно-го алгоритма MSA (СМ) на простом примере
L = —ж(1) + ж(2) ^ min, —ж(1) + ж(2) ^ 1, — 2ж(1) + ж(2) ^ 2, 3ж(1) + ж(2) ^ 3
с допустимыми решениями в треугольнике с вершинами А(0,5, 1,5), В(0,2, 2,4), С(—1, 0) и оптимумом на ребре АС.
В программе linprog не предусмотрено выявление неединственности решения, хотя уже сам факт возникновения такой ситуации в практических задачах может быть весьма весомым. Опция LSA вышла на внутреннюю точку ж* = (—0,2743, 0,7257) ребра АС, а MSA — на вершину ж* = (—1, 0) и там остановилась.
Программа mpao нашла обе оптимальные вершины ж1 = (0,5,1,5) и ж2 = (—1, 0), т. е. всё ребро АС. Вычисленные двойственные переменные в обоих решениях имеют одинаковые соответствующие значения (1, 0,0), что объясняется следующим обстоятельством. В первом решении активными являются первое и третье ограничения, а во втором — первое и второе. Превратим активные ограничения в свободные и объединим их правые части A(k) в вектор A = (A(1),A(2),A(3)). При условии сохранения базы оптимальной вершины (теорема 5) прямыми выкладками получаем L(X1(A)) = L(X2(A)) = A(1). Обозначив оптимум двойственных переменных y* = (y(1)* , y(2)* , y(3)*), следуя формуле (6) и [14], находим
(1), = dL(x1(A)) = dL(X2(A)) = _ (fcr = dL(X1(A)) = dL(X2(A)) , = 2 3 У dA(1) dA(1) 1 У dA(k) dA(k) 0 k 2 3
т.е. у двойственной задачи оптимум у* = (1, 0, 0) единственен, что и подтверждается её непосредственным решением.
Замена первого ограничения на —ж(1)+ж(2) ^ 3 превращает систему условий в несовместную. Реакция linprog: решение не найдено. Диагностика трао: совместность нарушает ограничение 3 (правостороннее). Другой порядок активации ограничений может изменить номер нарушителя.
Изъятие первого неравенства влечёт неограниченность решения. Сообщение linprog: решение не ограничено, ж* = (2,46 ■ 1031, —7,39 ■ 1031), Ь = —9,85 ■ 1031. Ответ трао: решение неограничено по переменным ж(1) и ж(2), Ь = —то.
Масштабная апробация программы трао производилась на базе задач, созданных с помощью специальной модификации трао, для которой исходными данными являются количество переменных и ограничений, а также наличие или отсутствие требования совместности системы условий. Используя генератор случайных чисел, эта программа выдаёт коэффициенты целевой функции и ограничений, а также лимиты ограничений и переменных. В таблице приведено время счёта для задач с размерностями в диапазоне 10 х 10 ^ 1000 х 1000 для программ трао и linprog на одноядерном ПК с частотой процессора 1,73 ГГц и ОЗУ 1 Гб. Столь маломощный компьютер был выбран умышленно для обеспечения ощутимости времени счёта на задачах малой размерности.
Размер трао ЬБЛ МБЛ
т х п мин. сред. макс. мин. сред. макс. мин. сред. макс.
10 х 10 0,00 0,00 0,08 0,03 0,04 0,67 0,01 0,025 0,562
20 х 20 0,00 0,01 0,05 0,03 0,05 0,09 0,06 0,074 0,219
30 х 30 0,03 0,04 0,12 0,11 0,13 0,78 0,16 0,185 0,719
40 х 40 0,09 0,11 0,12 0,12 0,13 0,17 0,34 0,369 0,454
50 х 50 0,13 0,14 0,2 0,12 0,14 0,20 0,73 0,758 0,875
80 х 80 0,73 0,76 0,89 0,40 0,44 0,56 2,89 3,067 3,485
100 х 100 1,44 1,47 1,52 0,76 0,79 0,86 7,03 7,333 7,593
200 х 200 16,4 16,6 16,9 6,09 6,21 6,31 116 116,8 118,3
300 х 300 60,2 60,6 60,9 6,75 6,98 7,61 310 310,5 311,1
300 х 500 327 328 329 21,8 21,9 22,4 617 618,3 622,4
500 х 300 98,9 101 101 17,1 17,3 17,8 874 874,7 875,3
500 х 500 437 442 446 36,4 36,7 37,0 2346 2351 2356
500 х 1000 3393 3410 3433 447 471 516 5902 5904 5906
1000 х 500 805 851 942 45,5 46,5 47,7 10895 11113 11505
1000 х 1000 5854 5862 5869 328 333 345 46091 48961 51223
Примечание. Указаны минимальное (мин.), максимальное (макс.) и усреднённое (сред.) время счёта в секундах по 1000 прогонов для каждой из малоразмерных задач с постепенным снижением числа прогонов до трёх по мере роста размерности. Усреднённое время рассчитано как среднеарифметическое по числу прогонов.
В экспериментальных расчётах последовательность активации ограничений заметно влияла на количество итераций, и это обстоятельство логично использовать с максимальной выгодой. Пока что для повышения эффективности алгоритма порядок активации устанавливается исходя из просто реализуемых и по возможности малообременительных правил. Одним из них является выбор кандидата на активацию с максимальной значимостью. Под мерой значимости может пониматься, например, возможность наибольшего снижения максимума целевой функции при активации ограничения, что подразумевает выявление кандидата на активацию на каждом шаге алгоритма. Активацию ограничений можно производить также в порядке убывания угла между нормалями целевой функции и ограничений. Детали такой стратегии устанавливаются один раз изначально, основанием для этого служит следующее соображение. В Е2 одна из нормалей базовой гиперплоскости оптимальной вершины составляет с
нормалью целевой функции минимальный угол. Хотя уже в E3 это вовсе не обязательно, результаты экспериментальных расчётов показывают, что такая последовательность активации, не будучи наилучшей, сокращает количество итераций на десятки процентов и потому вполне оправдана. В случае нахождения оптимальной стратегии активации ограничений традиционным неопровержимым аргументом против принципиальной возможности понижения экспоненциальной оценки числа итераций МПАО будет контрпример.
Анализ полученных результатов свидетельствует, что, в противоположность MSA, быстродействие mpao и LSA в большей мере зависит от числа переменных, нежели от количества ограничений. Это обстоятельство обеспечивает первой ветви алгоритма (с исключением части переменных через равенства) вычислительное преимущество по сравнению со второй ветвью (решение задачи непосредственно с равенствами). На малых размерностях до 50 х 50 уверенно лидирует mpao, а далее при стабильном удержании ею второго места быстро растущее превосходство переходит к опции LSA, которая демонстрирует слабую зависимость времени счёта от размерности задач. На размерности 1000 х 1000 MSA многократно уступает LSA и mpao. На ПК с двухъядерным процессором с частотой 2 ГГц и ОЗУ 2 Гб вычисления ускорились примерно на 10 % с сохранением соотношения величин этого показателя для всех трёх программ.
В условиях наличия высокопроизводительных компьютеров быстродействие далеко не всегда является основным критерием эффективности алгоритма и соответствующего программного продукта, на первое место могут выдвигаться их функциональные и сервисные возможности.
Заключение
Сопоставление СМ и МПАО приводит к следующим выводам. Оба метода — точные и конечные. Решая задачу на максимум, СМ на каждом шаге пытается увеличить текущее значение целевой функции, продвигаясь вдоль одной из координатных осей в границах МДР, т. е. максимум достигается снизу. Для превращения неравенств в равенства и нахождения начальной точки вводятся искусственные переменные, что формально повышает размерность задачи.
МПАО, напротив, предпочитает сокращать количество переменных путём их исключения через посредство равенств, поскольку его вычислительная эффективность обратно пропорциональна числу неизвестных. Последовательно активируя ограничения, алгоритм усекает начальный гиперпараллелепипед до заданного МДР, сужая область допустимых решений и сохраняя состояние оптимальности. В результате этого максимальное значение целевой функции на каждом шаге уменьшается или, по крайней мере, не увеличивается и приближение к максимуму выглядит как спуск сверху. Начальная точка маршрута определяется тривиально, и последующее движение к оптимальной вершине происходит по траектории, представляющей собой ломаную прямую, вершинами которой являются оптимальные вершины промежуточных многогранников. Эта траектория пролегает вне МДР; если в ходе активации ограничений не выявлена их несовместность, то его первой же захваченной алгоритмом точкой как раз и будет искомая оптимальная вершина. В связи с этим напрашиваются такие параллели: симплекс метод — метод внутренней точки, метод последовательной активации ограничений — метод внешней точки. В отличие от геометрического метода, одолевающего задачи максимум с тремя переменными с визуализацией МДР, MnAO, базируясь на геометрических представлениях, решает задачи любой размерности. МПАО активно использует органичную связь прямой и двойственной задач. Это проявляется в том,
что решающее правило выбора оптимального ребра (5) вычисляет двойственную оценку ограничения на текущей стадии его активации. После нахождения оптимума задачи посредством псевдоактивации каждого из базовых ограничений решающее правило находит оптимальное решение двойственной задачи. Последнее действие вычислительно малообременительно, поскольку к его началу известны направляющие векторы рёбер многогранника допустимых решений, выходящих из оптимальной вершины, и активные в ней ограничения. Таким образом, в МПАО двойственность играет решающую роль на всех стадиях вычислительного процесса.
Анализ результатов экспериментальных расчётов свидетельствует о теоретической и практической состоятельности МПАО и превосходстве в быстродействии программы трао над опцией симплекс-метода программы linprog, растущим с размерностью задач. Потенциал программы позволяет: констатировать неограниченность целевой функции и указывать ответственные за это переменные; выявлять ограничения, вызывающие несовместность системы условий; в отличие от linprog, в случае неединственности описывать структуру множества решений в виде смежных оптимальных вершин; определять чувствительность оптимума относительно активных ограничений посредством двойственных оценок; подобно прожектору, высвечивать структуру МДР в окрестности оптимума, что в случае неполной достоверности исходных данных имеет практическое значение. Резервы повышения быстродействия МПАО ещё не исчерпаны, но уже сейчас всего лишь с экспериментальной программой трао он уверенно занимает среднее положение между СМ и МВТ и по ряду параметров более информативен. Стремительное продвижение к оптимуму МВТ вне конкуренции и точным методом пока недостижимо.
ЛИТЕРАТУРА
1. Канторович Л. В. Математические методы в организации и планировании производства. Л.: Изв. Ленингр. гос. ун-та, 1939. 67с.
2. Данциг Дж. Линейное программирование, его применение и обобщения. М.: Прогресс, 1966. 600 с.
3. Юдин Д. Б., Гольштейн Е. Г. Линейное программирование: теория, методы и приложения. М.: Наука, 1969. 424 с.
4. Дикин И. И. Итеративное решение задач линейного и квадратичного программирования // Докл. АН СССР. 1967. Т. 174. С. 747-748.
5. Хачиян Л. Г. Полиномиальный алгоритм в линейном программировании // Докл. АН СССР. 1979. Т. 244. С. 1093-1096.
6. Зоркальцев В. И. Двойственные алгоритмы внутренних точек // Изв. вузов. Математика. 2011. №4. С. 33-53.
7. Зоркальцев В. И., Медвежонков Д. С. Численные эксперименты с вариантами алгоритмов внутренних точек на нелинейных задачах потокораспределения // Управление большими системами. Ин-т систем энергетики им. Л. А. Мелентьева СО РАН, 2013. Вып. 46. С.68-87.
8. Медвежонков Д. С. Экспериментальные исследования алгоритмов внутренних точек на нелинейных задачах потокораспределения // Вестник Бурятского гос. ун-та. 2013. №9. С. 12-16.
9. Вылегжанин О. Н., Шкатова Г. И. Решение задачи линейного программирования с использованием оператора-проектора // Изв. Томского политехн. ун-та. 2009. Т. 314. №5. С.37-40.
10. Бахшиян Б. Ц., Гориянов А. В. Скелетный алгоритм решения задачи линейного программирования и его применение для решения задач оценивания // Вестник МАИ. 2008. Т. 15. №2. С. 5-16.
11. Гордуновский В. М. Метод экспоненциальной аппроксимации для линейного программирования. М.: Эдитус, 2013. 32 с.
12. Колосов В. С. Конечный параметрический метод решения задачи линейного программирования // Науч. труды МЛТИ. 1978. Вып. 103. С. 197-198.
13. Колосов В. С. Роль двойственных оценок в параметрическом методе решения задачи линейного программирования // Науч. труды МЛТИ. 1986. Вып. 183. С. 51-54.
14. Карманов В. Г. Математическое программирование. М.: ФИЗМАТЛИТ, 2004. 264с.
REFERENCES
1. Kantorovich L. V. Matematicheskie metody v organizacii i planirovanii proizvodstva [Mathematical Methods in the Organization and Planning of Production]. Leningrad, News of Leningrad State University, 1939. 67 p. (in Russian)
2. Dantzig J. Linejnoe programmirovanie, ego primenenie i obobshcheniya [Linear Programming, its Application and Generalizations]. Moscow, Progress Publ., 1966. 600p. (in Russian)
3. Yudin D. B. and Golstein E. G. Linejnoe programmirovanie: teoriya, metody i prilozheniya [Linear Programming: Theory, Methods and Applications]. Moscow, Nauka Publ., 1969. 424 p. (in Russian)
4. Dikin 1.1. Iterativnoe reshenie zadach lineynogo i kvadratichnogo programmirovaniya [Iterative solution of problems of linear and quadratic programming]. Report AS USSR, 1967, vol. 174, pp. 747-748. (in Russian)
5. Khachiyan L. G. Polinomial'nyy algoritm v lineynom programmirovanii [Polynomial algorithm in linear programming]. Report AS USSR, 1979, vol.244, pp. 1093-1096. (in Russian)
6. Zorkaltsev V. I. Dvoystvennye algoritmy vnutrennikh tochek [Dual algorithms of interior-points]. Iz. VUZ, Mathematics, 2011, no. 4 pp. 33-53. (in Russian)
7. Zorkaltsev V. I. and Medvejonkov D. S. Chislennye eksperimenty s variantami algoritmov vnutrennikh tochek na nelineynykh zadachakh potokoraspredeleniya [Computational experiments with variants of interior-point algorithms for nonlinear flow distribution problems]. Institute of Energy Systems of SB RAS, 2013, no. 46, pp. 68-87. (in Russian)
8. Medvejonkov D. S. Eksperimental'nyye issledovaniya algoritmov vnutrennikh tochek na nelineynykh zadachakh potokoraspredeleniya [Experimental researches of interior- point algorithms for solution of flow distribution nonlinear problems]. Bulletin of the Buryat State University, 2013, no. 9, pp. 12-16. (in Russian)
9. Vylegzhanin O N. and Skatova G. I. Reshenie zadachi lineynogo programmirovaniya s ispol'zovaniem operatora-proektora [The solution of the linear programming problem using the operator- projector]. News Tomsk Polytechnic University, 2009, vol. 314, no. 5, pp. 37-40. (in Russian)
10. Bakhshiyan B. Ts. and Gurianov A. V. Skeletnyy algoritm resheniya zadachi lineynogo programmirovaniya i ego primenenie dlya resheniya zadach otsenivaniya [The skeletal algorithm for solving the linear programming problem and its application for solving estimation problems]. Bulletin MAI, 2008, vol.15, no. 2, pp. 5-16. (in Russian)
11. Gordunovski V. M. Metod ehksponencial'noj approksimacii dlya linejnogo programmirovaniya [The Method of Exponential Approximation for Linear Programming]. Moscow, Editus Publ., 2013. 32 p. (in Russian)
12. Kolosov V. S. Konechnyy parametricheskiy metod resheniya zadachi lineynogo programmirovaniya [The finite parametric method for solving linear programming problem]. Scientific Works MLTI, 1978, vol. 103, pp. 197-198. (in Russian)
13. Kolosov V. S. Rol' dvoystvennykh otsenok v parametricheskom metode resheniya zadachi lineynogo programmirovaniya [The role of dual estimates in the parametric method for solving the linear programming problem]. Scientific Works MLTI, 1986, vol. 183, pp. 51-54. (in Russian)
14. Karmanov V. G. Matematicheskoe programmirovanie [Mathematical Programming]. Moscow, FIZMATLIT Publ., 2004. 264 p.
