Алгоритмы квадратичного программирования

Квадратичное определение программирования

Квадратичное программирование - задача нахождения вектора x, который минимизирует квадратичную функцию, возможно подверженную линейным ограничениям:

\underset{}{} \frac{}{}^{} minx12xTHx +^{}

cTx

(1)

так, что A· x ≤ b, Aeq· x = beq, l ≤ x ≤ u.

`interior-point-convex` `quadprog` Алгоритм

interior-point-convex алгоритм выполняет следующие шаги:

Преколве/Постсольве
Создать начальную точку
Предиктор-корректор
Условия остановки
Обнаружение несходимости

Примечание

Алгоритм имеет два кодовых пути. Он принимает один, когда матрица Гессена H является обычной (полной) матрицей двойников, и он принимает другой, когда H является разреженной матрицей. Дополнительные сведения о разреженном типе данных см. в разделе Разреженные матрицы. Как правило, алгоритм быстрее для больших задач, которые имеют относительно мало ненулевых членов, когда вы указываете H как sparse. Аналогично, алгоритм быстрее для небольших или относительно плотных проблем, когда вы указываете H как full.

Преколве/Постсольве

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

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

Если алгоритм обнаруживает неосуществимую или неограниченную проблему, он останавливается и выдает соответствующее сообщение о выходе.

Алгоритм может прийти к одной осуществимой точке, которая представляет решение.

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

Подробнее см. Gould and Toint [63].

Создать начальную точку

Начальная точка x0 для алгоритма:

Инициализировать x0 кому ones(n,1), где n - количество строк в Н.
Для компонентов, имеющих обе верхние границы ub и нижняя граница lb, если компонент x0 не находится строго внутри границ, для компонента установлено значение (ub + lb)/2.
Для компонентов, имеющих только одну границу, при необходимости измените компонент, чтобы он лежал строго внутри границы.
Сделайте шаг предиктора (см. Предиктор-корректор) с незначительными коррекциями для осуществимости, а не полный шаг предиктор-корректор. Это помещает начальную точку ближе к центральному пути, не влекя за собой накладные расходы на полный шаг предиктора-корректора. Для получения подробной информации о центральном пути см. Nocedal and Wright [7], стр. 397.

Предиктор-корректор

Разреженные и полные внутренние выпуклые алгоритмы отличаются главным образом фазой предиктор-корректор. Алгоритмы схожи, но отличаются в некоторых деталях. Описание основного алгоритма см. в разделе Mehrotra [47].

Алгоритмы начинаются с превращения линейных неравенств Ax < = b в неравенства вида Ax > = b умножением A и b на -1. Это не имеет никакого отношения к решению, но делает проблему такой же формы найденной в некоторых литературах.

Редкий предиктор-корректор
Полный предиктор-корректор

Разреженный предиктор - корректор. Аналогично fmincon алгоритм внутренней точки, разреженный interior-point-convex алгоритм пытается найти точку, где сохраняются условия Каруша-Куна-Такера (KKT). Для задачи квадратичного программирования, описанной в Quadratic Programming Definition, существуют следующие условия:

$\begin{matrix} _{}^{} {\overset{}{}}^{} \\ \overset{}{} \overset{}{} \\ _{}_{} \\ _{}_{} Hx+c−AeqTy−A¯Tz=0A¯x−b¯−s=0Aeqx−beq=0sizi=0, i=1,2 ..., \\ ms≥0 \\ z≥0. \end{matrix}$

Здесь

$\overset{}{}$ A - это расширенная линейная матрица неравенства, которая включает границы, записанные как линейные неравенства. $\overset{b}{}$ - соответствующий линейный вектор неравенства, включая границы.
s - вектор слабых элементов, преобразующих ограничения неравенства в равенства. s имеет длину m, число линейных неравенств и границ.
z - вектор множителей Лагранжа, соответствующий s.
y - вектор множителей Лагранжа, связанных с ограничениями равенства.

Алгоритм сначала предсказывает шаг из формулы Ньютона - Рафсона, затем вычисляет шаг корректора. Корректор пытается лучше применить нелинейное ограничение _sizi = 0.

Определения для шага предиктора:

_rd, двойной остаток:
$_{rd} = Hx +_{c}^{} - {\overset{}{}}^{} AeqTy$ − A = Tz.
_req, остаточное ограничение основного равенства:
$_{req} =_{} {Aeqx}_{-}$ beq.
_rineq, основной остаток ограничения неравенства, который включает границы и слабости:
$_{rineq} \overset{}{=} A\nx \overset{}{} -$ b\-s.
_rsz, остаточная комплементарность:
_rsz = Sz.
S - диагональная матрица слабых членов, z - столбчатая матрица множителей Лагранжа.
_rc, средняя комплементарность:
$_{rc} \frac{=^{}}{}$ sTzm.

На шаге Ньютона изменения в x, s, y и z задаются следующим образом:

(\begin{matrix} H0 & -_{}^{} & {\overset{}{AeqT}}^{-} \\ _{A} & = \\ \overset{}{} \\ TAeq000A & = \end{matrix} \begin{matrix} I000ZOS) \end{matrix} (\begin{matrix} _{} \\ _{} \\ {ΔxΔsΔyΔz}_{) =} \\ -_{(} \end{matrix}

rdreqrineqrsz).

(2)

Однако полный шаг Ньютона может оказаться неосуществимым из-за ограничений позитивности на s и z. Поэтому quadprog сокращает шаг, при необходимости, для поддержания позитива.

Кроме того, для поддержания «центрированного» положения во внутреннем пространстве, вместо попытки решить _sizi = 0, алгоритм берёт положительный параметр λ, и пытается решить

_sizi = _startrc.

quadprog заменяет _rsz в уравнении шага Ньютона _{на rsz} + ΔsΔz - _startrc1, где 1 - вектор единиц. Также, quadprog переупорядочивает уравнения Ньютона для получения симметричной, более численно стабильной системы для вычисления шага предиктора.

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

Полный предиктор-корректор. Полный алгоритм предиктора-корректора не объединяет границы в линейные ограничения, поэтому он имеет другой набор переменных слабости, соответствующих границам. Алгоритм сдвигает нижние границы на ноль. И если на переменной имеется только одна граница, алгоритм превращает её в нижнюю границу нуля, отрицая неравенство верхней границы.

$\overset{}{}$ A - это расширенная линейная матрица, которая включает в себя как линейные неравенства, так и линейные равенства. $\overset{b}{}$ - соответствующий линейный вектор равенства $\overset{}{.}$ Кроме того, элемент «s» включает в себя элементы для расширения вектора «x» с переменными «s», которые превращают ограничения неравенства в ограничения равенства:

$\overset{}{A}_x \begin{matrix} _{=} & ( \end{matrix} \begin{matrix} _{Aeq0AI} \\ ) \end{matrix} ($ x0s),

где x0 означает исходный x-вектор.

Условия KKT:

\begin{matrix} {\overset{}{}}^{} \\ \overset{}{} \overset{}{} \\ _{}_{} Hx+c−A¯Ty−v+w=0A¯x=b¯x+t=uvixi=0, i=1,2 ..., \\ _{}_{} mwiti=0, i=1,2 ..., \\ nx, v, w, t≥0. \end{matrix}

(3)

Чтобы найти решение x, переменные слабости и двойственные переменные для уравнения 3, алгоритм в основном рассматривает шаг Ньютона-Рафсона:

(\begin{matrix} H & {\overset{}{-}}^{A} \\ \overset{}{T0} & − \\ IIA \\ 0000 & − \\ I0 & - \end{matrix} \begin{matrix} \begin{array}{l} I00V00X000W0T \end{array} \end{matrix}) (\begin{matrix} {\overset{}{}}^{} ΔxΔyΔtΔvΔw) = \\ \overset{}{-} (\overset{}{Hx} \\ + c - \\ \begin{array}{l} A \\ Ty \end{array} \end{matrix} \begin{matrix} -_{} \\ v_{} \\ +_{} \\ \begin{array}{l} {wA}_{} \\ _{x} \end{array} \end{matrix} -

b u − x − tVXWT) = − (rdrprubrvxrwt),

(4)

где X, V, W и T - диагональные матрицы, соответствующие векторам x, v, w и t соответственно. Остаточные векторы в крайней правой части уравнения:

_rd, двойной остаток
_rp, основной остаток
_руб., остаток верхней границы
_rvx, остаток комплементарности нижней границы
_rwt, остаток комплементарности верхней границы

Алгоритм решает уравнение 4, сначала преобразуя его в симметричную матричную форму

(\begin{matrix} - & {\overset{DA}{}}^{} \\ \overset{}{DTA} \end{matrix} \begin{matrix} = \\ 0) \end{matrix} (\begin{matrix} {ΔxΔy}_{)} \end{matrix} =

− (Rrp),

(5)

где

$\begin{array}{l} D = H^{+} X^{-} \\ 1V +_{T} -^{}_{1WR} =^{-}_{rd} -^{X} -_{} \end{array}$ 1rvx + T − 1rwt + T − 1Wrub.

Все обратные матрицы в определениях D и R просты в вычислении, поскольку матрицы диагональны.

Чтобы вывести уравнение 5 из уравнения 4, обратите внимание, что вторая строка уравнения 5 такая же, как и вторая строка матрицы уравнения 4. Первая строка уравнения 5 является результатом решения двух последних строк уравнения 4 для Δv и Δw, а затем решения для Δt.

Для решения уравнения 5 алгоритм следует существенным элементам Альтмана и Гондцио [1]. Алгоритм решает симметричную систему разложением ЛПНП. Как указывают такие авторы, как Вандербей и Карпентер [2], это разложение является численно устойчивым без какого-либо поворота, поэтому может быть быстрым.

Полное quadprog алгоритм предиктор-корректор в значительной степени такой же, как в linprog 'interior-point' алгоритм, но включает и квадратичные члены. См. раздел Предиктор-корректор.

Ссылки

[1] Альтман, Анна и Дж. Гондзио. Регуляризованные симметричные неопределенные системы в методах внутренних точек для линейной и квадратичной оптимизации. Методы оптимизации и программное обеспечение, 1999. Доступно для загрузки здесь.

[2] Вандербей, Р. Дж. и Т. Дж. Карпентер. Симметричные неопределенные системы для методов внутренних точек. Математическое программирование 58, 1993. стр 1–32. Доступно для загрузки здесь.

Условия остановки

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

$start= max (1, ‖ \overset{H}{} ‖,_{‖} A ‖, \overset{}{‖} Aeq_{‖,} ‖$ c ‖, ‖ b ‖, ‖ beq ‖)

Алгоритм останавливается, когда выполняются все эти условия:

$\begin{matrix} {‖_{}}_{} {_{}}_{} \\ {_{}}_{} \\ _{} rp‖1+‖rub‖1≤ρTolCon‖rd‖∞≤ρTolFunrc≤TolFun \end{matrix}$ ,

где

$_{rc} \underset{}{=} макси ({мин}_{} (_{|}_{xivi} |,_{|} xi |, |_{} {vi}_{} |),_{} мин (_{|}$ tiwi |, | ti |, | wi |)).

_rc по существу измеряет размер остатков комплементарности xv и tw, которые являются векторами нулей в решении.

Обнаружение несходимости

quadprog вычисляет функцию «За заслуги» в каждой итерации. Функция оценки заслуг является показателем осуществимости. quadprog останавливается, если функция заслуг становится слишком большой. В этом случае quadprog объявляет проблему неосуществимой.

Функция качества связана с условиями KKT для задачи - см. Предиктор-корректор. Используйте следующие определения.

$\begin{matrix} start= max (1, ‖ \overset{H}{} ‖,_{‖} A ‖, \overset{}{‖} Aeq_{‖,} ‖ \\ _{c} ‖,_{‖} b_{} \\ ‖_{, ‖} \overset{beq}{} ‖ \overset{}{)} \\ {req}_{} = Aeqx_{−}^{}_{} {\overset{}{}}^{beqrineq} {\overset{=}{}}_{A wetx} \\ \frac{}{−}^{b} -^{} {\overset{}{swsrd}}^{=} {\overset{}{}}_{Hx +} c_{+}^{}_{} \end{matrix}$ AeqTλ eq + A wetλ veneqg = 12xTHx + ct

Обозначение $\overset{}{A}$ = и $\overset{b}{}$ _ n означает коэффициенты линейного неравенства, дополненные слагаемыми для представления границ для разреженного алгоритма. ${\overset{}{Обозначение}}_{λ veneq}$ аналогично представляет множители Лагранжа для линейных ограничений неравенства, включая связанные ограничения. Это называлось z в Predictor-Corrector, $_{а}$ λ eq называлось y.

Функция «За заслуги»

$\frac{}{1group(} max {(_{‖}}_{} {req‖∞}_{, ‖}_{} {_{} rineq‖\infty}_{,} ‖$ rd‖∞) + g).

Если эта функция заслуг становится слишком большой, quadprog объявляет проблему неосуществимой и останавливается с флагом выхода -2.

`trust-region-reflective` `quadprog` Алгоритм

Многие методы, используемые в решателях Optimization Toolbox™, основаны на областях доверия, простой, но мощной концепции оптимизации.

Чтобы понять подход доверительной области к оптимизации, рассмотрим задачу неограниченной минимизации, минимизируйте f (x), где функция принимает векторные аргументы и возвращает скаляры. Предположим, что вы находитесь в точке x в n-пространстве и хотите улучшить, т.е. переместитесь в точку с меньшим значением функции. Основная идея состоит в том, чтобы аппроксимировать f более простой функцией q, которая разумно отражает поведение функции f в окрестности N вокруг точки x. Эта окрестность является областью доверия. Пробный этап s вычисляется путем минимизации (или приблизительно минимизации) над N. Это подпроблема области доверия,

\underset{мин}{} {q (s), s\inN

(6)

Текущая точка обновляется до x + s, если f ( x + s ) < f (x); в противном случае текущая точка остается неизменной, и N, область доверия, уменьшается, и вычисление пробного этапа повторяется.

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

В стандартном способе доверительной области ([48]) квадратичное приближение q определяется первыми двумя членами приближения Тейлора к F при x; окрестность N обычно имеет сферическую или эллипсоидальную форму. Математически обычно указывается подпроблема области доверия

\min \frac{{}{}^{} 12 sTHs^{+} sTg так  , что ‖ Ds‖\leqΔ

(7)

где g - градиент f в текущей точке x, H - гессеновская матрица (симметричная матрица вторых производных), D - диагональная масштабная матрица, Δ - положительный скаляр, а ∥. ∥ - 2-норма. Для решения уравнения 7 существуют хорошие алгоритмы (см. [48]); такие алгоритмы обычно включают вычисление всех собственных значений H и процесс Ньютона, применяемый к светскому уравнению

$\frac{}{} \frac{}{} 1Δ−1‖s‖=0.$

Такие алгоритмы обеспечивают точное решение уравнения 7. Однако они требуют времени, пропорционального нескольким факторизациям Х. Поэтому для масштабных проблем необходим другой подход. Несколько аппроксимационных и эвристических стратегий, основанных на уравнении 7, были предложены в литературе ([42] и [50]). Подход аппроксимации, используемый в решателях Optimization Toolbox, заключается в ограничении подпроблемы области доверия двумерным подпространством S ([39] и [42]). После вычисления подпространства S работа по решению уравнения 7 является тривиальной, даже если необходима полная информация о собственном значении/собственном векторе (поскольку в подпространстве задача является только двумерной). Доминирующая работа теперь перешла к определению подпространства.

Двумерное подпространство S определяют с помощью процесса предварительно кондиционированного сопряженного градиента, описанного ниже. Решатель определяет S как линейное пространство, охватываемое s1 и s2, где s1 находится в направлении градиента g, а s2 является либо приближенным направлением Ньютона, т.е. решением

_{} H⋅s2=−g,

(8)

или направление отрицательной кривизны,

_{}^{}_{} s2T⋅H⋅s2<0.

(9)

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

Эскиз неограниченной минимизации с использованием идей области доверия теперь легко дать:

Сформулируйте двухмерную подпроблему области доверия.
Решите уравнение 7, чтобы определить пробный этап s.
Если f (x + s) < f (x), то x = x + s.
Отрегулируйте Δ.

Эти четыре шага повторяются до сходимости. Измерение Δ области доверия корректируется в соответствии со стандартными правилами. В частности, он уменьшается, если стадия испытания не принимается, то есть f (x + s) ≥ f (x). Для обсуждения этого аспекта см. [46] и [49].

Решатели Optimization Toolbox рассматривают несколько важных особых случаев f со специализированными функциями: нелинейные наименьшие квадраты, квадратичные функции и линейные наименьшие квадраты. Однако лежащие в основе алгоритмические идеи те же, что и для общего случая. Эти особые случаи рассматриваются в последующих разделах.

Метод области доверия подпространства используется для определения направления поиска. Однако вместо ограничения шага (возможно) одним шагом отражения, как в случае нелинейной минимизации, при каждой итерации проводится кусочно-отражающий поиск линии. Для получения подробной информации о поиске строк см. [45].

Метод предварительно кондиционированного сопряженного градиента

Популярным способом решения больших, симметричных, положительных определённых систем линейных уравнений Hp = -g является метод Предварительно кондиционированных сопряжённых градиентов (PCG). Этот итеративный подход требует возможности вычисления матрично-векторных произведений формы H· v, где v - произвольный вектор. Симметричная положительная определенная матрица М является предварительным условием для Н. То есть М = ^C2, где ^C-1HC-1 является хорошо кондиционированной матрицей или матрицей с кластеризованными собственными значениями .

В контексте минимизации можно предположить, что матрица Гессена H симметрична. Однако H гарантированно является положительным определенным только в районе сильного минимизатора. Алгоритм PCG выходит, когда он сталкивается с направлением отрицательной (или нулевой) кривизны, то есть ^dTHd ≤ 0. Направление р выхода PCG является либо направлением отрицательной кривизны, либо приближенным решением для системы Ньютона Hp = -g. В любом случае p помогает определить двумерное подпространство, используемое в подходе «доверительная область», обсуждаемом в документе «Методы доверительной области для нелинейной минимизации».

Ограничения линейного равенства

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

Общая задача минимизации с ограничением линейного равенства может быть записана

\min {f (x ) так , что Ax = b

(10)

где A - матрица m-на-n (m ≤ n). Некоторые решатели Optimization Toolbox предварительно обрабатывают А для удаления строгих линейных зависимостей методом, основанным на факторизации LU ^AT [46]. Здесь A принимается рангом m.

Метод, используемый для решения уравнения 10, отличается от неограниченного подхода двумя существенными способами. Сначала вычисляется начальная осуществимая точка x0 с использованием разреженного шага наименьших квадратов, так что Ax0 = b. Во-вторых, алгоритм PCG заменяется на уменьшенные предварительно кондиционированные сопряженные градиенты (RPCG), см. [46], чтобы вычислить приблизительный уменьшенный шаг Ньютона (или направление отрицательной кривизны в нулевом пространстве A). Ключевой шаг линейной алгебры включает в себя решение систем вида

[\begin{matrix} {\overset{}{}}^{} \\ \overset{}{} & CA˜TA˜0 \end{matrix}] \begin{matrix} [ \end{matrix} st] \begin{matrix} = \end{matrix} [

r0],

(11)

где $\overset{A˜}{}$ аппроксимирует A (малые ненулевые значения A устанавливаются в ноль при условии, что ранг не теряется), а C является разреженным симметричным положительно-определенным приближением к H, т.е. C = H. Подробнее см. [46].

Ограничения коробки

Проблема, связанная с рамкой, имеет вид

\min {f (x ) так , что l≤x≤u},

(12)

где l - вектор нижних границ, а u - вектор верхних границ. Некоторые (или все) компоненты l могут быть равны - ∞, а некоторые (или все) компоненты u могут быть равны ∞. Метод генерирует последовательность строго возможных точек. Два метода используются для поддержания осуществимости при достижении надежного поведения сходимости. Сначала масштабированный измененный шаг Ньютона заменяет неограниченный шаг Ньютона (для определения двумерного подпространства S). Во-вторых, отражения используются для увеличения размера шага.

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

{(D (x)}^{)} - 2

g = 0,

(13)

где

$D (x) = {{diag}_{} (}^{| vk} |$ − 1/2),

и вектор v (x) определяется ниже для каждого 1 ≤ i ≤ n:

Если _gi < 0 и _ui < ∞ то _vi = _xi - _ui
Если _gi ≥ 0 и _li > - ∞ то _vi = _xi - _li
Если _gi < 0 и _ui = ∞ то _vi = -1
Если _gi ≥ 0 и _li = - ∞ то _vi = 1

Уравнение 13 нелинейной системы не дифференцируется везде. Недифференцируемость возникает, когда _vi = 0. Избежать таких моментов можно, сохранив строгую выполнимость, т.е. ограничив l < x < u.

Масштабированный модифицированный шаг Ньютона _sk для нелинейной системы уравнений, заданный уравнением 13, определяется как решение линейной системы.

\overset{}{M}^^{} DsN \overset{=}{}

− g ^

(14)

на k-ой итерации, где

\overset{}{g}^^{=} D − 1g =^{diag} (|

v | 1/2) g,

(15)

\overset{}{M}^^{=} D^{-} 1HD − 1 +^{}

diag (g) Jv.

(16)

Здесь ^Сп играет роль якобианца | v |. Каждая диагональная составляющая диагональной матрицы Jv равна 0, -1 или 1. Если все компоненты l и u являются конечными, Jv = diag (знак (g)). В точке, _где gi = ₀, vi может быть не $_{дифференцируемым .}^{}$ В такой точке определяется Jiiv = 0. Недифференцируемость этого типа не вызывает беспокойства, поскольку для такого компонента не является существенным, какое _{значение} принимает vi. _{Кроме того}, | vi | все еще будет прерывистым в этой точке, но _{функция} | vi |· gi непрерывна.

Во-вторых, отражения используются для увеличения размера шага. (Один) этап отражения определяется следующим образом. Учитывая шаг p, который пересекает ограничивающее ограничение, рассмотрим первое ограничивающее ограничение, пересекаемое p; предположим, что это ограничение i-ой границы (i-я верхняя или i-я нижняя граница). Затем шаг отражения ^pR = p, за исключением i-ой составляющей, где ^pRi = -_pi.

`active-set` `quadprog` Алгоритм

После завершения этапа предварительного решения, active-set алгоритм работает в две фазы.

Фаза 1 - получение возможной точки относительно всех ограничений.
Фаза 2 - итеративное снижение целевой функции при сохранении списка активных ограничений и сохранении выполнимости в каждой итерации.

active-set стратегия (также известная как метод проекции) аналогична стратегии Gill et al., описанной в [18] и [17].

Шаг предварительного решения

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

Алгоритм может прийти к одной осуществимой точке, которая представляет решение.

Подробнее см. Gould and Toint [63].

Алгоритм фазы 1

В фазе 1 алгоритм пытается найти точку x это удовлетворяет всем ограничениям без учета целевой функции. quadprog запускает алгоритм Phase 1, только если предоставленная начальная точка x0 неосуществим.

Для начала алгоритм пытается найти точку, которая осуществима в отношении всех ограничений равенства, таких как X = Aeq\beq. При отсутствии решения x к уравнениям Aeq*x = beqзатем алгоритм останавливается. При наличии решения X, следующим шагом является удовлетворение границ и линейных неравенств. В случае отсутствия ограничений равенства X = x0, начальная точка.

Начиная с X, алгоритм вычисляет M = max(A*X – b, X - ub, lb – X). Если M <= options.ConstraintTolerance, то точка X выполнимо, и алгоритм фазы 1 останавливается.

Если M > options.ConstraintToleranceалгоритм вводит неотрицательную переменную ослабления γ для задачи вспомогательного линейного программирования

$\underset{minx,}{}$ γ γ

такой, что

$\begin{matrix} Ax−γ≤bAeq \\ x=beqlb−γ≤x≤ub+γγ≥−ρ . \end{matrix}$

В данном случае start- это ConstraintTolerance опция, умноженная на абсолютное значение наибольшего элемента в A и Aeq. Если алгоритм находит γ = 0 и точку x, которая удовлетворяет уравнениям и неравенствам, то x является возможной точкой фазы 1. Если нет решения вспомогательной задачи линейного программирования x с γ = 0, то задача фазы 1 неосуществима.

Для решения задачи вспомогательного линейного программирования алгоритм устанавливает γ 0 = M + 1, устанавливает x0 = Xи инициализирует активный набор как фиксированные переменные (если таковые имеются) и все ограничения равенства. Алгоритм переформулирует линейные программные переменные p, чтобы быть смещением x от текущей точки x0, а именно x = x0 + p. Алгоритм решает задачу линейного программирования теми же итерациями, что и в фазе 2 для решения задачи квадратичного программирования, с соответствующим образом модифицированным Гессеном.

Алгоритм фазы 2

В терминах переменной d проблема заключается в

\begin{matrix} \underset{^{}}{} mind∈ℜnq (d) \frac{}{=}^{}^{12dTHd} + \\ _{} cTd,_{} Aid = bi, i_{=} \\ _{1}, . ._{.},_{} meAid\leqbi, i = \end{matrix}

me + 1,..., m.

(17)

Здесь _Ай ссылается на iтретья строка матрицы А «m-на-n».

Во время фазы 2 активный набор ${\overset{}{A _}}_{k,}$ который является оценкой активных ограничений (на границах ограничений) в точке решения.

Алгоритм обновляет ${\overset{}{A _}}_{k}$ на каждой итерации k, формируя основу для направления поиска dk. Ограничения равенства всегда остаются в активном ${\overset{}{наборе}}_{A _}$ k. _{Направление} поиска dk вычисляется и минимизирует целевую функцию, оставаясь на любых активных границах ограничения. Алгоритм формирует возможное подпространство _для dk из _базиса Zk, столбцы которого ортогональны оценке активного ${\overset{}{}}_{множества}$ A _ k ( ${\overset{}{то}}_{} {есть}_{} A$ _ k _ Zk = 0). Поэтому направление поиска, которое формируется из линейного суммирования любой комбинации _{столбцов} Zk, гарантированно остается на границах активных ограничений.

Алгоритм формирует матрицу _Zk из последних m - l столбцов QR-разложения матрицы ${\overset{}{A}}_{wetkT,}^{}$ где l - число активных ограничений и l < m То есть Zk задаётся

_{Zk} = Q [:, l +

1: m],

(18)

где

$^{} {\overset{}{QTA}}_{ПДТ}^{} \begin{matrix} = \end{matrix} [$ R0].

После нахождения _Zk алгоритм ищет новое направление поиска _dk, которое минимизирует q (d), где _dk находится в нулевом пространстве активных ограничений. То есть _dk является линейной комбинацией столбцов _Zk: ${\overset{d}{}}_{}^_{k}$ = Zkp для некоторого вектора p.

Просмотр квадратичного как функции p путем замены на _dk дает

q (p) \frac{}{=}^{}_{}^{}_{}^{12pTZkTHZkp}_{+}

cTZkp.

(19)

Дифференцирование этого уравнения относительно выходов p

∇q (р)_{=}^{}_{}_{ZkTHZkp}^{} +

ZkTc.

(20)

∇q (p) называется спроецированным градиентом квадратичной функции, поскольку он является градиентом, спроецированным в подпространстве, определенном _Zk. Термин $_{}^{}_{}$ ZkTHZk называется проектируемым гессенским. Предполагая, что спроецированная гессеновская матрица $_{}^{}_{}$ ZkTHZk является положительной полуопределённой, минимум функции q (p) в подпространстве, определяемом Zk, возникает, когда ∇q (p ) = 0, что является решением системы линейных уравнений

_{}^{}_{} ZkTHZkp =_{-}^{}

ZkTc.

(21)

Алгоритм затем делает шаг формы

$_{xk +} 1_{} =_{xk}$ + αdk,

где

$_{dk} =_{}$ Zkp.

Из-за квадратичной природы целевой функции на каждой итерации существует только два выбора длины шага α. Шаг единицы вдоль _dk является точным шагом к минимуму функции, ограниченной нулевым пространством ${\overset{}{A _}}_{k}$ . Если алгоритм может сделать такой шаг, не нарушая ограничений, то этот шаг является решением квадратичной программы ( уравнение 18). В противном случае шаг вдоль dk до ближайшего ограничения меньше единицы, и алгоритм включает новое ограничение в активный набор при следующей итерации. Расстояние до границ ограничения в любом направлении dk задается как

$\underset{α=mini∈{1, ..., м}{}} \frac{{-_{} (_{}_{Aixk}}{-_{}_{bi}})$ Айдк},

которая определена для ограничений, не входящих в активный набор, и где направление _dk равно границе ограничения, то есть $_{}_{Aidk} > 0, i = 1, .$ .., m.

Когда активное множество включает в себя n независимых ограничений, без расположения минимума, алгоритм вычисляет множители Лагранжа _{λ k}, которые удовлетворяют неингулярному множеству линейных уравнений

{\overset{}{A}}_{=}^{}_{} kTλ k =_{c}

+ Hxk.

(22)

Если все элементы _{λ k} являются положительными, _xk является оптимальным решением задачи квадратичного программирования Уравнение 1. Однако если какая-либо составляющая _{λ k} отрицательна, и составляющая не соответствует ограничению равенства, то минимизация не является полной. Алгоритм удаляет элемент, соответствующий самому отрицательному множителю, и запускает новую итерацию.

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

Примечание

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

Иногда, active-set алгоритм может иметь трудности с обнаружением, когда проблема является неограниченной. Эта проблема может возникнуть, если направление неограниченности v является направлением, в котором квадратичный член v 'Hv = 0. Для численной стабильности и обеспечения факторизации Холеского, active-set алгоритм добавляет небольшой строго выпуклый член к квадратичной цели. Этот небольшой член заставляет объективную функцию быть ограниченной от -Inf. В этом случае active-set алгоритм достигает предела итерации вместо сообщения о том, что решение является неограниченным. Другими словами, алгоритм останавливается с флагом выхода 0 вместо -3.

Ссылки

[1] Гилл, П. Э., В. Мюррей, М. А. Сондерс и М. Х. Райт. Практическая процедура антициклирования для линейно ограниченной оптимизации. Математическое программирование 45 (1), август 1989, стр. 437-474.

Теплый старт

При запуске quadprog или lsqlin 'active-set' алгоритм с теплым стартовым объектом в качестве начальной точки, решатель пытается пропустить многие шаги фазы 1 и фазы 2. Объект теплого начала содержит активный набор ограничений, и этот набор может быть правильным или близким к исправлению новой проблемы. Поэтому решатель может избежать итераций для добавления ограничений в активный набор. Кроме того, начальная точка может быть близка к решению новой проблемы. Дополнительные сведения см. в разделе optimwarmstart.

Документация