Журнал вычислительной математики и математической физики, 2019, T. 59, № 10, стр. 1649-1665

ВВЕДЕНИЕ

В данной статье излагается история формирования семейства эффективных алгоритмов метода внутренних точек. Основой этих алгоритмов послужила идея 1965 г. Л. Канторовича оценки множителей Лагранжа при неоптимальном плане методом наименьших квадратов [1]–[4]. В развитие этой идеи И. Дикиным был создан и исследован алгоритм [5]–[7], получивший впоследствии название “affine scaling method”. До середины 80-х годов алгоритмы этого типа развивались в России, в Сибирском энергетическом институте (И. Дикин, В. Зоркальцев) и в Вычислительном центре (Ю. Евтушенко, В. Жадан) Академии наук СССР [4], [8]–[10]. Накоплен большой опыт использования этих алгоритмов в моделях энергетики [4], [11]–[13]. Широкую популярность во многих странах эти алгоритмы получили после многообещающей статьи 1984 г. Н. Кармаркара [14]. Частными случаями рассматриваемых в данной статье алгоритмов является affine scaling method, dual affine scaling method [15]–[17] и другие родственные алгоритмы, имеющие свои преимущества и недостатки.

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

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

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

Заданы матрица A размера $m \times n$, векторы $c \in {{R}^{n}}$, $b \in {{R}^{m}}$. Векторы $x \in {{R}^{n}}$, $u \in {{R}^{m}}$ составляют переменные этих задач.

1. СЕМЕЙСТВО АЛГОРИТМОВ

Вычисления можно начинать с любого вектора ${{x}^{0}} \in {{R}^{n}}$, у которого все компоненты положительные, что будем представлять в виде неравенства ${{x}^{0}} > 0$. Вырабатывается последовательность векторов ${{x}^{k}} > 0$, где $k = 0,1,2,\; \ldots $, – номер итерации.

На итерации k вычисляется вектор невязок ограничений-равенств задачи (1)

Затем решается система линейных уравнений относительно вектора переменных $u \in {{R}^{m}}$: Решение этой системы обозначим через ${{u}^{k}}$.

В (4) матрица ${{D}_{k}} = \operatorname{diag} {{d}^{k}}$ – диагональная матрица, образованная вектором весовых коэффициентов ${{d}^{k}} \in {{R}^{n}}$. Компоненты вектора ${{d}^{k}}$должны удовлетворять неравенствам

Здесь $\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\sigma } $, $\bar {\sigma }$ – любые непрерывные функции от положительного аргумента, удовлетворяющие двум условиям: для некоторых $M > 0$, $\delta > 0$. Условиям (5)–(7) удовлетворяет много конкретных правил задания весовых коэффициентов. Поэтому описываемый вычислительный процесс представляет семейство алгоритмов.

Будем считать, что система $Ax = b$ совместна и $\operatorname{rank} A = m$. Тогда существует и единственное решение системы (4). Вычислив ${{u}^{k}}$, находим направление и шаг корректировки решения:

В (9), (10) имеем при заданном $\gamma \in (0,1)$. Если ${{s}^{k}} \geqslant 0$ при ${{r}^{k}} \ne 0$, то полагаем ${{\lambda }_{k}} = 1$.

При ${{r}^{k}} = 0$ ситуация ${{s}^{k}} = 0$ означает, что полученное решение ${{x}^{k}}$ из X является оптимальным, $\bar {X} = X$. Ситуация ${{s}^{k}} \geqslant 0$, ${{s}^{k}} \ne 0$ при ${{r}^{k}} = 0$ означает, что $\bar {X} = \phi $, целевая функция задачи (1) неограниченна снизу на множестве X. Далее считаем, что $\bar {X} \ne \phi $, $\bar {X} \ne X$.

Итеративный переход осуществляется по правилу

Из (3)–(12) вытекают следующие свойства вычислительного процесса:

Поскольку ${{x}^{k}} > 0$, $0 < \gamma < 1$, то из (13) получаем ${{x}^{{k + 1}}} > 0$.

В изложенном вычислительном процессе можно выделить два этапа.

Этап 1. Ввод в область допустимых решений. Если ${{r}^{k}} \ne 0$, то, согласно (14), абсолютные величины невязок ограничений-равенств задачи (1) сокращаются в $(1 - {{\lambda }_{k}})$ раз. Этим объясняется, почему при вводе в область допустимых решений величина шага ${{\lambda }_{k}}$ в правиле вычисления (9) ограничивается сверху единицей. Если ${{\lambda }_{k}} = 1$, то вектор ${{x}^{{k + 1}}}$ будет удовлетворять ограничениям-равенствам задачи (1), ${{r}^{{k + 1}}} = 0$.

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

Значения дополнительной переменной образуют последовательность $x_{{n + 1}}^{0} = 1$, $x_{{n + 1}}^{k} = (1 - {{\lambda }_{{k - 1}}})x_{{n + 1}}^{{k - 1}}$, $k = 1,2,\; \ldots $. Вектор ${{x}^{k}}$ и величина $x_{{n + 1}}^{k}$, $k = 1,2,\; \ldots $, составляют допустимое, монотонно улучшающееся по итерациям решение задачи (15).

Отметим, что вектор ${{s}^{k}}$ можно рассматривать как результат решения следующей вспомогательной задачи с вектором переменных $s \in {{R}^{n}}$:

В условии этой задачи отражена приоритетная цель – минимизация абсолютных значений невязок ограничений-равенств задачи (1). Вектор ${{u}^{k}}$ составляют множители Лагранжа ограничений данной вспомогательной задачи. В целевой функции вспомогательной задачи (16) учитывается целевая функция исходной задачи (1). Это позволяет получать первое допустимое решение, близкое к оптимальному решению.

Этап 2. Оптимизация в области допустимых решений. Если ${{r}^{k}} = 0$, то

на данной и последующих итерациях осуществляется процесс улучшения решения в области допустимых решений.

Подмножества алгоритмов. Из рассматриваемого семейства выделим подмножество алгоритмов, обладающих более сильным, чем (7), свойством: при некотором $M > 0$ для любых $\beta \in (0,1]$, $\alpha \in (0,\beta ]$ имеем

Выделим еще более узкий класс алгоритмов: для изменяющихся $\alpha > 0$, $\beta \in (0,1]$ имеем

Примеры правил задания весовых коэффициентов. Условиям (5), (6), (17) удовлетворяет правило

при заданном $p \geqslant 1$. При $p > 1$ выполняется условие (18).

Условиям (5), (6), (17) удовлетворяет также следующее правило: при заданном $\varepsilon > 0$ для $k \geqslant 1$

Условие невырожденности. Если для любого $x \in X$ либо не существует, либо существует и только единственный вектор $u \in {{R}^{m}}$, при котором

то задачу (1) будем называть невырожденной. Аналогично определяется невырожденность задач (2), (15).

В условии невырожденности (21) не требуется, чтобы вектор u был допустимым решением задачи (2). Соотношение (21) для допустимых решений обоих задач (1), (2), т.е. для $x \in X$, $u \in U$ является условием дополняющей нежесткости, гарантирующим оптимальность векторов x и u. Известно, что, если задачи (1), (2) имеют оптимальные решения, то условие (21) для некоторых $x \in \bar {X}$, $u \in \bar {U}$ выполняется в строгой форме, когда

Это означает, что $x \in \operatorname{ri} \bar {X}$, $u \in \operatorname{ri} \bar {U}$, где $\operatorname{ri} \bar {X}$, $\operatorname{ri} \bar {U}$ – подмножества относительно внутренних точек оптимальных решений задач (1), (2).

Приводимые ниже теоремы содержат утверждения, анонсированные в [18] и доказанные в [19].

Теорема 1. Если задача (15) невырожденная и имеет оптимальное решение с ${{x}_{{n + 1}}} = 0$, $x > 0$, то через конечное число итераций будет получено решение ${{x}^{k}} > 0,\,\,\,{{x}^{k}} \in X.$

Если задача (1) невырожденная, $\bar {X} \ne \phi $, $\bar {X} \ne X$, ${{x}^{k}} \in X$ на некоторой итерации k, то верно следующее:

1) при выполнении условий (5)–(7), для некоторых $\bar {x} \in \bar {X}$, $\bar {u} \in \operatorname{ri} \bar {U}$, ${{x}^{k}} \to \bar {x}$, ${{u}^{k}} \to \bar {u}$ при $k \to \infty $;

2) если выполняется условие (17), то $\bar {x} \in \operatorname{ri} \bar {X}$, величины $\left\| {{{x}^{k}} - \bar {x}} \right\|$, $\left\| {{{u}^{k}} - \bar {u}} \right\|$ сходятся к нулю при $k \to \infty $ линейно;

3) если выполняется условие (18), то

И, начиная с некоторой итерации, на всех последующих

Причем в алгоритме можно использовать итеративно изменяющийся параметр $\gamma $, сходящийся к единице, что в силу (23) позволяет получать сверхлинейную скорость сходимости.

Более полное теоретическое обоснование [20]–[24], не использующее предположение о невырожденности, имеется для алгоритмов с весовыми коэффициентами (19) для $p \in [1,3]$. При этом требуется выполнение дополнительного условия

Верно следующее утверждение.

Теорема 2. Пусть в вычислительном процессе (3), (4), (8)–(12) используются весовые коэффициенты (19) при $p \in [1,3]$ и при параметре $\gamma $, удовлетворяющем условию (24).

Если у задачи (1) имеется решение $x \in X$, $x > 0$, то через конечное число итераций будет получено такое решение ${{x}^{k}} \in X$, ${{x}^{k}} > 0$.

Если $\bar {X} \ne \phi $, $\bar {X} \ne X$ и, начиная с некоторой итерации, ${{x}^{k}} \in X$, осуществляется процесс оптимизации в области допустимых решений задачи (1), то существуют $\bar {x} \in {\text{ri}}\bar {X}$, $\bar {u} \in {\text{ri}}\bar {U}$ такие, что при некоторых $M > 0$, $\alpha \in (0,1)$ верно

При $p > 1$ справедливо (22) и следующее свойство:

В [20] было дано обоснование процесса оптимизации в области допустимых решений алгоритмом с $p = 1$. В [21] было дано обоснование для случая $p = 2$. В [23] дано обоснование для случая $p \in (1,3]$. Обоснование процесса ввода в область допустимых решений (первая часть утверждения в теореме 2) было дано в препринтах [23], [24].

2. МЕТОДИКА КАНТОРОВИЧА ОЦЕНКИ МНОЖИТЕЛЕЙ ЛАГРАНЖА ПРИ НЕОПТИМАЛЬНОМ ПЛАНЕ И АЛГОРИТМ ДИКИНА

При оптимизации в области допустимых решений, когда ${{r}^{k}} = 0$, вспомогательная задача (4) равносильна проблеме минимизации квадратичной выпуклой функции от вектора $u \in {{R}^{m}}$:

где Действительно, приравняв градиент функции ${{\Phi }_{k}}$ к нулю, получим систему (4).

В 1965 г. Л. Канторович предложил использовать задачу (25) для приближенной оценки множителей Лагранжа при неоптимальном допустимом решении ${{x}^{k}} \in X$, ${{x}^{k}} > 0$. При этом в качестве весовых коэффициентов Канторовичем предлагалось использовать $d_{j}^{k} = x_{j}^{k}$, $j = 1,\; \ldots ,\;n$. Это равносильно (19) при $p = 1$. Эта методика Канторовича направлена на минимизацию методом наименьших квадратов невязок в соотношении (21), в целях достижения наилучшего приближения в выполнении условия дополняющей нежесткости. Данная методика Канторовича предлагалась применительно к проблеме формирования рентных оценок сельскохозяйственных земель [1]–[4].

Из свойств метода наименьших квадратов следует, что вычисляемый по правилу (8) вектор ${{s}^{k}}$, на основе получаемого в результате решения задачи (25) вектора ${{u}^{k}}$, не выводит из области допустимых решений, $A{{s}^{k}} = 0$.

Причем, если $\bar {X} \ne X$, то ${{s}^{k}}$ будет направлением уменьшения значения целевой функции задачи (1), ${{c}^{{\text{т}}}}{{s}^{k}} < 0$.

На основе этих свойств И. Дикиным и С. Анцизом были разработаны и экспериментально апробированы первые варианты алгоритмов оптимизации в области допустимых решений рассматриваемого здесь типа [7].

Экспериментальные расчеты И. Дикина и С. Анцыза показали, что вместо весовых коэффициентов (19) при $p = 1$ лучше использовать весовые коэффициенты (19) при $p = 2$, т.е. лучше положить

В качестве шага корректировки решения первоначально использовалась величина

Для такого алгоритма И. Дикиным была доказана линейная сходимость векторов ${{x}^{k}}$ к точке $\bar {x} \in \operatorname{ri} \bar {X}$ при предположении о невырожденности задачи (1).

Исследования по развитию и теоретическому обоснованию непрерывных аналогов этих алгоритмов с 70-х годов осуществлялись также Ю. Евтушенко и В. Жаданом в Вычислительном центре АН СССР [8], [9].

Интерес к алгоритмам рассматриваемого здесь типа проявился за пределами России только с середины 80-х годов после многообещающей статьи Н. Кармаркара [14]. Анонсированная Н. Кармаркаром высокая вычислительная эффективность его алгоритма не подтвердилась. Вместе с тем его статья привела к огромному количеству публикаций по рассматриваемым в данной статье вычислительным процессам с весовыми коэффициентами (19) с $p = 2$ (например, [15]–[17]). Такой вариант алгоритмов внутренних точек получил название “affine scaling method” (ASM). Фактически он почти полностью соответствует исходному алгоритму И. Дикина.

В первоначальном алгоритме И. Дикина для ввода в область допустимых решений использовался прием, аналогичный применявшемуся в симплекс-методе – введение m дополнительных переменных по одной для каждого уравнения. Приведенный выше способ ввода в область допустимых решений и способ вычисления шага ${{\lambda }_{k}}$ были предложены и реализованы в алгоритме [10]. Алгоритмы рассматриваемого здесь типа уже с 70-х годов прошлого века применялись в ряде моделей энергетики Сибирского энергетического института СО АН СССР [4], [11]–[13]. В том числе они использовались для реализации моделей управления режимами функционирования электроэнергетических систем [12], в вычислительных комплексах анализа надежности электроснабжения [11], в расчетах термодинамического равновесия [13].

3. СПОСОБЫ ПРОТИВОДЕЙСТВИЯ ПОГРЕШНОСТЯМ ВЫЧИСЛЕНИЙ

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

Погрешности в решении системы линейных уравнений (4) проявятся, в частности, в том, что для вычисленного направления улучшения решения ${{s}^{k}}$ требуемое равенство $A{{s}^{k}} = {{r}^{k}}$ будет выполняться в приближенном виде

где ${{\rho }^{k}}$ – вектор, отражающий погрешности вычислений. При использовании эффективных алгоритмов решения системы (4), например, метод квадратного корня, абсолютные значения компонент вектора ${{\rho }^{k}}$ будут, как правило, относительно небольшими. Могут ли небольшие погрешности вычислений усиливаться по итерациям и приводить к большим погрешностям в вычислительном процессе?

При вводе в область допустимых решений не будет происходить накопление по итерациям погрешностей вычислений, поскольку рассчитываемый на очередной итерации по формуле (3) вектор невязок ограничений-равенств задачи (1) будет учитывать и погрешность вычислений предыдущей итерации. Естественно, что вследствие погрешностей вычислений невязки будут итеративно изменяться не в точном соответствии с формулой (14), а по несколько иному правилу. А именно, это будет происходить по следующей формуле:

Так как при вводе в область допустимых решений величина шага ${{\lambda }_{k}}$ находится в ограниченном интервале $(0,1]$, влияние погрешностей ${{\rho }^{k}}$ на отклонение от формулы (3) будет незначительным, не превышающим абсолютные значения компонент вектора ${{\rho }^{k}}$:

Оптимизация в области допустимых решений при реализации вычислительного процесса на ЭВМ должна осуществляться, конечно, не при точной реализации условия ${{r}^{k}} = 0$. Необходимо использовать более реалистичное условие близости всех компонент вектора ${{r}^{k}}$ нулю, например из [4], [10], условие

где $\varepsilon $ – заданная положительная малая величина. При выполнении такого условия следует на данной итерации переопределять значение вектора ${{r}^{k}}$, положить его равным нулевому вектору.

В процессе такой оптимизации в области допустимых решений невязки балансовых ограничений будут меняться из-за погрешностей вычислений. Через t итераций вектор невязок будет иметь значение:

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

Корректировка в области допустимых решений может осуществляться [4] по правилам, используемым при вводе в область допустимых решений. Поскольку при выборе направления улучшения решения в процессе ввода в область допустимых решений учитывается целевая функция исходной задачи (1), то выбираемое направление корректировки обычно близко к направлениям оптимизации в области допустимых решений.

Выбор способа задания весовых коэффициентов. Величина и последствия погрешностей вычисления направления улучшения решения существенно зависят от выбора правила задания весовых коэффициентов $d_{j}^{k}$. Для уменьшения погрешностей вычислений в решении вспомогательной задачи (4) целесообразно по возможности уменьшать разброс значений коэффициентов $d_{j}^{k}$. В этой связи вместо правила (19) более эффективна следующая его модификация:

при заданных $N > 0$, $p \geqslant 1$. Такая модификация была рекомендована в [4], [10]. В этом случае сохраняются все теоретические результаты для алгоритмов с весовыми коэффициентами (19). Параметр N рекомендуется задавать примерно равным среднеарифметическому значению переменных задачи (1). Можно использовать изменяющийся по итерациям параметр N. Правила его задания нуждаются в дополнительных теоретических и экспериментальных исследованиях.

Рассмотрим правила (19) или (27) при $p > 1$. Это включает и известный affine scaling method, соответствующий случаю $p = 2$. Для таких алгоритмов теоретически обоснована высокая, асимптотически независимая от исходных данных задачи скорость сходимости процесса оптимизации в области допустимых решений, при $k \to \infty $ имеем

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

Большим недостатком таких вариантов алгоритмов в вычислительном отношении является бесконечное возрастание по итерациям величины шага ${{\lambda }_{k}}$ при оптимизации в области допустимых решений. Это теоретически доказанный и наблюдаемый в расчетах факт.

С ростом ${{\lambda }_{k}}$ малые погрешности в решении вспомогательной задачи поиска направления улучшения решения будут усиливаться в соответствии с (26) и будут быстро приводить к большим нарушениям ограничений-равенств задачи (1). Из-за этого возникает необходимость частой, порой через каждую итерацию оптимизации, корректировки в область допустимых решений.

Для алгоритмов с весовыми коэффициентами (19) или (7) с $p = 1$ величина шага ${{\lambda }_{k}}$ при оптимизации в области допустимых решений ограничена сверху некоторым положительным числом на всех итерациях. Это является также теоретически установленным и реально наблюдаемым фактом. Этот факт можно рассматривать как достоинство алгоритмов c $p = 1$. Недостатком таких алгоритмов является медленная скорость сходимости к оптимальному решению. Это подтверждают экспериментальные расчеты и имеющиеся теоретические результаты.

Для алгоритмов с $p = 1$ свойство (28) не выполняется. Асимптотическая при $k \to \infty $ скорость сходимости в этом случае зависит от исходных данных задачи. Пусть задача (1) имеет оптимальные решения и невырожденна. Тогда существует и единственное решение двойственной задачи (2) $\bar {u} \in \bar {U}$, к которому сходится последовательность векторов ${{u}^{k}}$, $k = 0,1,2,\; \ldots $. Введем множество номеров компонент вектора $g(\bar {u})$ с положительными значениями

Согласно теореме 1 алгоритм с весовыми коэффициентами (19) или (27) с $p = 1$ будут приводить последовательность векторов ${{x}^{k}}$, $k = 0,1,2,\; \ldots ,$ к некоторому решению $\bar {x} \in \operatorname{ri} \bar {X}$. Для этого вектора нулевые значения будут иметь компоненты с номерами из I: Пусть Можно доказать, что для любого $j \in I$ при $k \to \infty $ имеем Это означает, что алгоритмы с $p = 1$ имеют более медленную асимптотическую скорость сходимости, чем алгоритмы с $p > 1$, поскольку величины ${{g}_{j}}(\bar {u}){\text{/}}{v}$ для некоторых $j \in I$ могут быть и, обычно бывают, меньше единицы. Соотношение (30) означает также, что асимптотическая скорость сходимости для алгоритмов с $p = 1$, в отличие от алгоритмов с $p > 1$, зависит от исходных данных задач (1), (2).

В целях объединения преимуществ, рассмотренных выше двух вариантов весовых коэффициентов, было введено правило (20). Это правило целесообразно дополнить алгоритмом итеративного изменения коэффициента $\varepsilon $, т.е. модифицировать правило (20) к виду

где ${{\varepsilon }_{k}}$ – положительные величины. Они должны итеративно изменяться таким образом, чтобы не быть слишком малыми и при этом в пределе не превосходить ${{g}_{j}}(\bar {u})$ для j, при которых эти значения положительные. Здесь $\bar {u}$ – предельное значение последовательности векторов ${{u}^{k}}$, $k = 0,1,2,\; \ldots $. Например, можно воспользоваться правилом Это правило для невырожденной задачи (1) будет приводить к следующему предельному значению последовательности величин ${{\varepsilon }_{k}}$: где $I$ – множество номеров (29).

Для алгоритмов с весовыми коэффициентами (20) или (31) при оптимизации в области допустимых решений шаг ${{\lambda }_{k}}$ ограничен сверху. При этом, если выполняется (33), то имеет место асимптотическая скорость сходимости (28). Первоначально введенное в [18] правило (20) рекомендовалось использовать после нескольких итераций алгоритмов с весовыми коэффициентами (27) с $p = 2$. Как показали результаты экспериментальных расчетов, правилом (20) можно пользоваться и с самого начала вычислительного процесса.

В дополнение к правилу (32) возможно введение ограничения сверху на значения весовых коэффициентов, как в формуле (27). Правила задания и итеративного управления параметрами $N$ и ${{\varepsilon }_{k}}$ нуждаются в дополнительных теоретических и экспериментальных исследованиях.

4. ДВОЙСТВЕННЫЕ АЛГОРИТМЫ

В изложенном в разд. 2 вычислительном процессе происходит монотонное улучшение по итерациям решения задачи (1). Сначала уменьшаются абсолютные величины невязок ограничений-равенств. Затем происходит монотонное улучшение целевой функции в области допустимых решений задачи (1). При этом итеративно вырабатываемые приближения ${{u}^{k}}$ для задачи (2), хотя и сходятся к оптимальному решению этой задачи, но не гарантированно монотонно в каком-либо смысле.

Рассматриваемый в данном разделе вычислительный процесс является симметричным аналогом вычислительного процесса раздела (2). Монотонное улучшение решения осуществляется для двойственной задачи (2). Начальным приближением могут служить любые векторы ${{u}^{0}} \in {{R}^{m}}$, ${{g}^{0}} \in {{R}^{n}}$, ${{g}^{0}} > 0$.

Исходными на итерации $k = 0,1,2,\; \ldots $ являются векторы ${{u}^{k}} \in {{R}^{m}}$, ${{g}^{k}} \in {{R}^{n}}$, ${{g}^{k}} > 0$. Пусть

где вектор ${{d}^{k}}$ из ${{R}^{n}}$ должен задаваться таким, чтобы выполнялись неравенства Здесь $\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\sigma } $, $\bar {\sigma }$ – некоторые непрерывные функции от положительного аргумента, удовлетворяющие условиям (6), (7).

Вычисляем векторы

Находим при заданном $\gamma \in (0,1)$_,значения Осуществляем итеративный переход

Из приведенных расчетных формул следует, что на всех итерациях

При ${{\delta }^{k}} \ne 0$ имеет место этап ввода в область допустимых решений задачи (2). При ${{\delta }^{k}} = 0$ реализуется этап оптимизации в области допустимых решений. И, если $\bar {U} \ne U$, $\bar {U} \ne \phi $, то

Приведенные здесь двойственные алгоритмы в более общем виде впервые были изложены в [18], где были анонсированы утверждения по их обоснованию аналогичные утверждениям теоремы 1. Доказательства этих утверждений опубликованы в [25].

Аналогом (19) является следующее правило вычисления весовых коэффициентов при $p \geqslant 1$

Алгоритм оптимизации в области допустимых решений с весовыми коэффициентами (38) при $p = 2$ получил название “dual affine scaling method” [26].

Аналогом правила (20) является следующий способ задания весовых коэффициентов: при $k \geqslant 1$ для некоторого $\varepsilon > 0$ имеем

Для правил (38), (39) можно высказать рекомендации по их уточнению и развитию, аналогичные высказанным в разд. 3 для правил (19), (20).

Альтернативный путь вычислений. Основной вычислительной проблемой при реализации изложенного здесь двойственного вычислительного процесса является решение системы линейных уравнений с симметричной неотрицательно определенной матрицей размера m, к чему сводится правило (35) вычисления вектора ${{r}^{k}}$. Вместо (35)–(37) можно использовать иной путь вычислений, основной проблемой в котором является поиск решения системы линейных уравнений с симметрической положительно-определенной матрицей размера n. Этот путь представляется в виде следующей, равносильной (35)–(37), последовательности вычислений

Вспомогательная задача поиска направления улучшения решения. Представленные выше два варианта вычислений при реализации двойственных алгоритмов внутренних точек представляют два способа решения следующей вспомогательной задачи поиска векторов $r \in {{R}^{m}}$, $z \in {{R}^{n}}$ в результате минимизации квадратичной выпуклой функции при линейных ограничениях в форме равенств:

где Решением этой задачи являются векторы ${{r}^{k}}$ и ${{z}^{k}}$. Вектор ${{x}^{k}}$ состоит из множителей Лагранжа ограничений (44).

Исходя из (44), можно выразить вектор z через вектор r и подставить это выражение в целевую функцию (45). Придем к задаче безусловной минимизации квадратичной выпуклой функции от вектора переменных r. Приравняв градиент этой функции нулю, приходим к расчетным формулам (35)–(37).

Другой путь вычислений основывается на использовании для решения задачи (43), (44) метода неопределенных множителей Лагранжа. Обозначим через x вектор искомых (неопределенных) множителей Лагранжа ограничений (44). Из условий оптимальности следует, что

Выразим векторы z и r через x:

Подставив эти выражения в условие (44), получим систему линейных уравнений относительно вектора переменных xРешая эту систему, согласно (40), находим вектор ${{x}^{k}}$ и затем, подставляя ${{x}^{k}}$ в выражение (46), (47), по правилам (41), (42) находим ${{r}^{k}}$ и ${{z}^{k}}$.

Возможные способы повышения точности решения вспомогательной задачи поиска направления улучшения решения. Одним из способов “противодействия” погрешностям вычислений при реализации как прямых, так и двойственных алгоритмов внутренних точек является повышение точности решения системы линейных уравнений (4) или, для двойственных алгоритмов, системы (35) либо системы (40). Эти системы имеют симметрические неотрицательно-определенные матрицы. Причем у системы (40) матрица положительно-определенная. Для систем (4), (35) матрицы будут положительно-определенными, если $\operatorname{rank} A = m$. Для решения таких систем линейных уравнений эффективен, как известно, метод квадратного корня. Обычно на базе этого метода и реализуются рассмотренные здесь алгоритмы внутренних точек. Можно предложить три направления повышения точности решения вспомогательной задачи и “противодействия” погрешностям вычислений.

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

2. Метод квадратного корня не использует имеющееся приближение к решению. Вместе с тем у систем (4), (35) и (40) имеются, особенно в конце вычислительного процесса, хорошее исходное приближение – решение этих систем на предыдущей итерации. Есть смысл использовать только итеративный метод решения указанных систем с имеющихся приближений. Для того, чтобы решения, полученные на предыдущих итерациях, были более близкими к требуемым на данной итерации, есть смысл уменьшать величину шага ${{\lambda }_{k}}$ и за счет этого уменьшать расхождения между величинами $d_{j}^{k}$ и $d_{j}^{{k + 1}}$.

3. Условие (34) означает, что для рассматриваемого вычислительного процесса в качестве значения $d_{j}^{k}$ подходит любая величина из интервала $[\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\sigma } (g_{j}^{k}),\bar {\sigma }(g_{j}^{k})]$. Поэтому весовые коэффициенты $d_{j}^{k}$ можно рассматривать как переменные, искомые величины из указанных интервалов. Этот факт можно эффективно использовать для повышения точности решения системы (40). Приведем один из возможных вариантов правила вычисления, использующего эту возможность.

Пусть ${{\tilde {D}}_{k}}$ – диагональная матрица размера n, состоящая из коэффициентов

где $N$ – заданная положительная величина.

Эта правило является развитием (39) и может рассматриваться как уточненный симметричный аналог правила (32). Введем нижнюю и верхнюю границу диапазонов, в которых могут находиться весовые коэффициенты. Пусть

Рассмотрим систему (48), у которой в качестве матрицы ${{D}_{k}}$ используется определенная здесь матрица ${{\tilde {D}}_{k}}$. Пусть вектор ${{x}^{k}}$ является решением этой системы с некоторой погрешностью. Если бы решение было точным, то выполнялись бы при ${{d}_{j}} = \tilde {d}_{j}^{k}$ равенства

где Равенство (49) при ${{d}_{j}} = \tilde {d}_{j}^{k}$ может точно не выполняться, но, за счет корректировок величины ${{d}_{j}}$, можно сократить, даже порой до нуля, невязки требуемого равенства (49). Пусть w – заданная положительная малая величина. Если $\left| {x_{j}^{k}} \right| \leqslant w$, то полагаем $d_{j}^{k} = \tilde {d}_{j}^{k}$. В этом случае весовые коэффициенты и невязки условия (49) не изменяются.

Пусть $\left| {x_{j}^{k}} \right| > w$. Вычисляем

Возможно три случая:

1) если $\underline d _{j}^{k} \leqslant \hat {d}_{j}^{k} \leqslant \bar {d}_{j}^{k}$, то полагаем $d_{j}^{k} = \hat {d}_{j}^{k}$;

2) если $\hat {d}_{j}^{k} < \underline d _{j}^{k}$, то полагаем $d_{j}^{k} = \underline d _{j}^{k}$;

3) если $\hat {d}_{j}^{k} > \bar {d}_{j}^{k}$, то полагаем $d_{j}^{k} = \bar {d}_{j}^{k}$.

Во всех трех случаях абсолютные значения невязок условия (49) будут сокращены. Причем в первом случае до нуля.

Критерий отсутствия решения. Для выявления ситуаций отсутствия решений у задач (1), (2) можно воспользоваться теоремами об альтернативных системах линейных неравенств [27]. Известно, что задачи (1), (2) не имеют решения в том и только том случае, если пусто хотя бы одно из множеств X или U. Множество X пусто в том и только том случае, если не пусто множество

Множество U пусто в том и только том случае, если не пусто множество

Для выявления ситуации отсутствия решения у задач (1), (2) алгоритмами первого раздела следует на каждой итерации проверять, находится ли вектор ${{u}^{k}}$ в множестве $\hat {U}$ и вектор ${{s}^{k}}$ в множестве $\hat {X}$. Для двойственных алгоритмов критерием отсутствия решений у задач (1), (2) может служить принадлежность на очередной итерации вектора ${{r}^{k}}$ множеству $\hat {U}$ и принадлежность вектора ${{z}^{k}}$ множеству $\hat {X}$. Как показывают результаты расчетов, по этим критериям обычно уже на начальных итерациях выявляется ситуация отсутствия решения.

5. ПРИМЕР ИСПОЛЬЗОВАНИЯ АЛГОРИТМОВ

С 70-х годов прошлого века рассмотренные здесь алгоритмы активно используются в Институте систем энергетики СО РАН (ранее называвшемся Сибирским энергетическим институтом) при реализации ряда моделей энергетики [4], [11]–[13], [28]. Экспериментальные исследования вариантов алгоритмов внутренних точек на конкретных моделях являются важным дополнением к теоретическим исследованиям.

В табл. 1 представлены результаты сравнительного исследования, выполненного А.Ю. Филатовым, четырех вариантов алгоритмов внутренних точек на линеаризованных задачах расчета допустимых режимов электроэнергетической системы. Исходная модель допустимых режимов нелинейная. Для ее реализации используется процедура итеративной линеаризации. На каждой итерации требуется найти решение или выявить несовместность системы линейных уравнений и неравенств следующего вида:

Здесь заданными являются: матрица A размера $m \times n$, векторы $b \in {{R}^{m}}$, $\underline x $, $\bar {x}$ из ${{R}^{n}}$. Переменные составляют вектор x из ${{R}^{n}}$. Исходные модели расчета допустимых режимов и их линеаризации были подготовлены О.Н. Войтовым.

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

Рассматривается два варианта прямых и два варианта двойственных алгоритмов внутренних точек. В алгоритме ${{F}_{1}}$ используются весовые коэффициенты (19) при $p = 2$. В алгоритме ${{F}_{2}}$ используются весовые коэффициенты (20). В алгоритмах ${{G}_{1}}$ и ${{G}_{2}}$ реализуются двойственные вычислительные процессы с весовыми коэффициентами (38) при $p = 2$ и (39). Алгоритм ${{F}_{1}}$ можно идентифицировать как ASM. Алгоритм ${{G}_{1}}$ можно считать соответствующим dual affine scaling method.

Рассматривалась система (50) с m, изменяющимся в диапазоне от 30 до 41, и n, равным во всех случаях 80. Эксперименты проводились с 14 совместными и 14 несовместными системами. В табл. 1 указано среднее значение числа итераций по всем 14 задачам и в скобках – диапазоны числа итераций по 14 задачам.

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

Данные табл. 1 наглядно демонстрируют практическую эффективность обсуждавшихся выше критериев несовместности. Несовместность ограничений нередко выявляется на первых итерациях. Особенно хорошо проявили себя в этом двойственные алгоритмы внутренних точек. Эти алгоритмы также лучше себя проявили при решении совместных систем. Наилучшим из четырех алгоритмов по результатам этих экспериментов можно признать алгоритм ${{G}_{2}}$.

6. ОТНОСИТЕЛЬНО ВНУТРЕННИЕ ТОЧКИ ОПТИМАЛЬНЫХ РЕШЕНИЙ

Обозначим ${{J}_{0}}(x)$, $J(x)$, ${{J}_{ + }}(x)$, ${{J}_{ - }}(x)$ множества номеров компонент вектора $x \in {{R}^{n}}$ с нулевыми, ненулевыми, положительными и отрицательными значениями. Множество $J(x)$ называется носителем вектора x.

В теории линейных неравенств большая роль отводится [29], [30] решениям системы неравенств с максимальными (нерасширяемыми) наборами активных ограничений. Для множества решений задачи линейного программирования (1) это будут оптимальные решения с минимальным (не сужаемым) носителем. К таким решениям обычно приводят алгоритмы симплекс-метода.

Алгоритмы внутренних точек приводят к “антиподам” таких решений, а именно к оптимальным решениям задачи (1) с максимальным (не расширяемым) носителем, к относительно внутренним точкам множества оптимальных решений.

Если для векторов $x \in X$, $u \in U$ условие дополняющей нежесткости, гарантирующее их оптимальность, состоит в выполнении соотношений

то для решений $\bar {x} \in \operatorname{ri} \bar {X}$, $\bar {u} \in \operatorname{ri} \bar {U}$ выполняются равенства

Это означает, что для произвольных оптимальных решений $x \in \bar {X}$, $u \in \bar {U}$ выполняются соотношения

Относительно внутренние точки оптимальных решений могут быть полезны, в частности, для анализа устойчивости множества оптимальных решений к варьированию исходных данных задачи (1). Располагая относительно внутренними точками оптимальных решений, можно эффективно описывать множество оптимальных решений, что полезно для лексикографической оптимизации. Если мы имеем произвольное оптимальное решение $\tilde {x} \in \bar {X}$, то множество оптимальных решений обычно задается путем введения дополнительного ограничения на значение целевой функции:

Если же известно, что мы располагаем относительно внутренней точкой оптимальных решений $\bar {x} \in \operatorname{ri} \bar {X}$, то множество оптимальных решений можно определять просто исключением из рассмотрения компонент вектора переменных, получивших нулевые значения для вектора $\bar {x}$:

При этом не вводится дополнительное ограничение-равенство на значение целевой функции задачи (1).

Описание (51) удобно, в частности, для решения многоэкстремальных задач лексикографической оптимизации [30]. В частности, такой была модель оценки дефицита мощности ЭЭС [4], [11]. В этой модели после решения задачи минимизации суммарного дефицита по системе возникала неоднозначность распределения этого суммарного дефицита по узлам системы. Поэтому в качестве задачи второго этапа вычислений осуществляется равномерное (по специальному линейному критерию) распределение дефицита между узлами ЭЭС. При этом использовалось представление (51), при котором в задаче второго этапа учитывались только потенциально дефицитные узлы ЭЭС с неоднозначным распределением этого дефицита.

Чебышёвская проекция на линейное многообразие. Пусть L – линейное многообразие в ${{R}^{n}}$, h – заданный вектор положительных весовых коэффициентов, $\alpha (x) = \mathop {\max }\limits_{j \in J} {{h}_{j}}\left| {{{x}_{j}}} \right|$ – взвешенная чебышёвская норма вектора $x \in {{R}^{n}}$, $J = \{ j = 1,\; \ldots ,\;n\} $. Рассмотрим задачу поиска чебышёвской проекции начала координат на линейное многообразие

В таком виде представляются многие прикладные проблемы, в том числе линейная чебышёвская аппроксимация.

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

Рассмотрим пример. Пусть $n = 2$, линейное многообразие задается условием ${{x}_{2}} = 1$, т.е. это прямая, параллельная оси абсцисс. Множество чебышёвских проекций начала координат на такое многообразие составляет интервал значений ${{x}_{1}}$ от $ - {{h}_{2}}{\text{/}}{{h}_{1}}$ до ${{h}_{2}}{\text{/}}{{h}_{1}}$. Причем по содержательным соображениям удовлетворительным является только значение ${{x}_{1}} = 0$.

Можно отметить, что именно в точке ${{x}_{1}} = 0$, ${{x}_{2}} = 1$ будет достигаться минимум на прямой ${{x}_{2}} = 1$ функции

при любом $p \geqslant 1$ и любых ${{h}_{1}} > 0$, ${{h}_{2}} > 0$. При $p = 1$ это будет октаэдральная проекция начала координат на линейное многообразие. При $p = 2$ – евклидова проекция, при $p \in (1,\infty )$ – гёльдеровская проекция. И только в предельном случае $p = \infty $, что соответствует чебышёвской проекции, будем иметь указанный выше интервал решений.

Прежде всего следует определиться с общим понятием “ближайшей к началу координат точки линейного многообразия”. В качестве самого общего подхода предлагается использовать парето-оптимальные решения многокритериальной задачи минимизации абсолютных значений всех компонент. Множество векторов линейного многообразия с парето-минимальными абсолютными значениями компонент обозначим через

В рассматриваемом выше примере множество Q будет состоять только из одного вектора с компонентами ${{x}_{1}} = 0$, ${{x}_{2}} = 1$.

Пусть S – линейное подпространство, параллельное L, т.е. множество, состоящее из векторов $s = x - y$, где x и y – векторы из L. Критерием принадлежности вектора множеству Q может служить следующее, доказываемое стандартно, утверждение.

Лемма 1. Вектор $x \in L$ находится в Q в том и только том случае, если не существует вектора $s \in S$, $s \ne 0$, при котором

Для преодоления проблемы возможной неоднозначности чебышёвских проекций обычно используется специальное ограничение на линейное многообразие L, на правила которыми оно задается. Это так называемое условие Хаара, фактически сводящееся к требованию единственности чебышёвской проекции. Условие Хаара ограничивает применение чебышёвской проекции. Оно порой трудно проверяемое и иногда нарушаемое. Предлагается пойти по другому пути – использовать в дополнение к задаче (52) лексикографическую оптимизацию с алгоритмами, приводящими к относительно внутренним точкам оптимальных решений.

Процесс вычисления чебышёвской проекции представим в виде конечной последовательности поиска относительно внутренних точек оптимальных решений задач линейного программирования. Пусть ${{J}_{1}} = J$, ${{S}^{1}} = S$. При $t = 1$ задача (52) представима в следующем виде:

Пусть ${{\alpha }^{t}}$ – оптимальное значение целевой функции этой задачи, ${{x}^{t}}$ – относительно внутренняя точка оптимальных решений, т.е. это оптимальное решение с минимальным набором ограничений (54), выполняющихся в виде равенства.

Пусть $I_{ + }^{t}$ – множество номеров $j \in {{J}_{t}}$, для которых

Пусть $I_{ - }^{t}$ – множество номеров $j \in {{J}_{t}}$, для которых Поскольку ${{x}^{t}}$ – решение задачи (53), (54) с минимальным набором активных ограничений, то не существует вектора $s \in {{S}^{t}}$, $s \ne 0$, при котором

Зафиксируем значения компонент вектора переменных с номерами из $I_{ + }^{t}$ и $I_{ - }^{t}$:

Положим

Если ${{J}_{{t + 1}}} \ne \phi $, то переходим к решению задачи (53), (54) при $t: = t + 1$ и при фиксированных по условию (56) значениях компонент вектора переменных на последнем и всех предыдущих этапах вычислений.

Поскольку множество ${{J}_{{t + 1}}}$ строго включено в ${{J}_{t}}$, то через конечное число этапов таких вычислений процесс завершится и множество ${{J}_{{t + 1}}}$ окажется пустым. Пусть это будет этап с номером $T$. Вектор ${{x}^{{\text{т}}}}$ определяется описанным вычислительным процессом однозначно. При этом справедлива

Теорема 3. Вектор ${{x}^{{\text{т}}}}$ находится в Q.

Доказательство. Это вытекает из (55) для $t = 1,\; \ldots ,\;T$ и леммы 1.

Поиск относительно внутренней точки оптимальных решений задачи (53), (54) может осуществляться рассмотренными в данной статье алгоритмами внутренних точек.