Метод Стронгина

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

Требуется найти точку ${\displaystyle x^{*}\in [a;\;b]}$ такую, что ${\displaystyle f(x^{*})=\min \left\{f(x)\colon x\in [a;\;b],\;g_{j}(x)\leqslant 0,\;1\leqslant j\leqslant m\right\}}$. Предполагается, что функции ${\displaystyle f(x)}$ и ${\displaystyle g_{j}(x),\;j={\overline {1,\;m}}}$ удовлетворяют условию Липшица на отрезке ${\displaystyle [a;\;b]}$.

Обозначим ${\displaystyle g_{m+1}(x)=f(x)}$, тогда для ${\displaystyle j={\overline {1,\;m+1}}}$ выполняются следующие неравенства:

${\displaystyle |g_{j}(x+\Delta x)-g_{j}(x)|\leqslant L_{j}\Delta x,\;a\leqslant x+\Delta x\leqslant b,}$

где ${\displaystyle L_{j}\geqslant 0}$ — константы Липшица.

Пусть ${\displaystyle Q_{0}=[a;\;b]}$. Ограничение, имеющее номер ${\displaystyle j}$, выполняется во всех точках области ${\displaystyle Q_{j}=\left\{x\in [a;\;b]\colon g_{j}(x)\leqslant 0\right\}}$, которая называется допустимой для этого ограничения. При этом допустимая область ${\displaystyle Q}$ исходной задачи определяется равенством:

${\displaystyle Q=\bigcap _{j=0}^{m}Q_{j}.}$

Испытание в точке ${\displaystyle x\in [a;\;b]}$ состоит в последовательном вычислении значений величин ${\displaystyle g_{1}(x),\;\ldots ,\;g_{\nu }(x)}$, где значение индекса ${\displaystyle \nu }$ определяется условиями:

${\displaystyle x\in Q_{j},\;0\leqslant j<\nu ,\;x\notin Q_{\nu }.}$

Выявление первого нарушенного ограничения прерывает испытание в точке ${\displaystyle x}$. В случае, когда точка ${\displaystyle x}$ допустима, то есть ${\displaystyle x\in Q}$ испытание включает в себя вычисление всех функций задачи. При этом значение индекса принимается равным величине ${\displaystyle \nu =m+1}$.

Пара ${\displaystyle \nu =\nu (x),\;z=g_{\nu }(x)}$, где индекс ${\displaystyle \nu }$ лежит в границах ${\displaystyle 1\leqslant \nu \leqslant m+1}$, называется результатом испытания в точке ${\displaystyle x}$.

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

${\displaystyle \psi (x^{*})=\min _{x\in [a;\;b]}\psi (x),}$

${\displaystyle \psi (x)={\begin{cases}g_{\nu }(x)/L_{\nu }&\nu <M,\\(g_{M}(x)-g_{M}^{*})/L_{M}&\nu =M.\end{cases}}}$

Здесь ${\displaystyle M=\max \left\{\nu (x)\colon x\in [a;\;b]\right\}}$, а ${\displaystyle g_{M}^{*}=\min \left\{g_{M}(x)\colon x\in \bigcap _{i=0}^{M-1}Q_{i}\right\}}$.

В силу определения числа ${\displaystyle M}$, задача отыскания ${\displaystyle g_{M}^{*}}$ всегда имеет решение, а если ${\displaystyle M=m+1}$, то ${\displaystyle g_{M}^{*}=f(x^{*})}$.

Дуги функции ${\displaystyle \psi (x)}$ липшицевы на множествах ${\displaystyle \bigcap _{i=0}^{j}Q_{i},\;0\leqslant j\leqslant M-1}$ с константой 1, а сама ${\displaystyle \psi (x)}$ может иметь разрывы первого рода на границах этих множеств.

Несмотря на то, что значения констант Липшица и величина ${\displaystyle g_{M}^{*}}$ заранее неизвестны, они могут быть оценены в процессе решения задачи.

Пусть ${\displaystyle x^{0}=a,\;x^{1}=b}$. Индексы концевых точек считаются нулевыми, а значения ${\displaystyle z}$ в них не определены. Первое испытание осуществляется в точке ${\displaystyle x^{3}=(a+b)/2}$. Выбор точки ${\displaystyle x^{k+1},\;k\geqslant 3}$ любого последующего испытания определяется следующими правилами:

1. Перенумеровать точки ${\displaystyle x^{0},\;\ldots ,\;x^{k}}$ ${\displaystyle k}$ предшествующих испытаний нижними индексами в порядке увеличения значений координаты: ${\displaystyle a=x_{0}<\ldots <x_{i}<\ldots <x_{k}=b}$ и сопоставить им значения ${\displaystyle z_{i}=g_{\nu }(x_{i}),\;\nu =\nu (x_{i}),\;i={\overline {1,\;k}}}$.
2. Для каждого целого числа ${\displaystyle \nu ,\;1\leqslant \nu \leqslant m+1}$ определить соответствующее ему множество ${\displaystyle I_{\nu }}$ нижних индексов точек, в которых вычислялись значения функций ${\displaystyle g_{\nu }(x)}$:

   ${\displaystyle I_{\nu }=\{i\colon \nu (x_{i})=\nu ,\;1\leqslant i\leqslant k\},\;1\leqslant \nu \leqslant m+1.}$

   Также определить максимальное значение индекса ${\displaystyle M=\max\{\nu (x_{i}),\;1\leqslant i\leqslant k\}.}$

3. Вычислить текущие оценки для неизвестных констант Липшица:

   ${\displaystyle \mu _{\nu }=\max\{|g_{\nu }(x_{i})-g_{\nu }(x_{j})|/(x_{i}-x_{j})\colon i,\;j\in I_{\nu },\;i>j\}.}$

   Если множество ${\displaystyle I_{\nu }}$ содержит менее двух элементов или если значение ${\displaystyle \mu _{\nu }}$ оказывается равным нулю, то принять ${\displaystyle \mu _{\nu }=1}$.

4. Для всех непустых множеств ${\displaystyle I_{\nu },\;\nu ={\overline {1,\;M}}}$ вычислить оценки

   ${\displaystyle z_{\nu }^{*}={\begin{cases}\min\{g_{\nu }(x_{i})\colon x_{i}\in I_{\nu }\}&\nu =M,\\-\varepsilon _{\nu }&\nu <M,\end{cases}}}$

   где вектор с неотрицательными координатами ${\displaystyle \varepsilon _{R}=(\varepsilon _{1},\;\ldots ,\;\varepsilon _{m})}$ называется вектором резервов.

5. Для каждого интервала ${\displaystyle (x_{i-1};\;x_{i}),\;1\leqslant i\leqslant k}$ вычислить характеристику

   ${\displaystyle R(i)={\begin{cases}\Delta _{i}+{\frac {(z_{i}-z_{i-1})^{2}}{(r_{\nu }\mu _{\nu })^{2}\Delta _{i}}}-2{\frac {z_{i}+z_{i-1}-2z_{\nu }^{*}}{r_{\nu }\mu _{\nu }}}&\nu =\nu (x_{i})=\nu (x_{i-1}),\\2\Delta _{i}-4{\frac {z_{i-1}-z_{\nu }^{*}}{r_{\nu }\mu _{\nu }}}&\nu =\nu (x_{i-1})>\nu (x_{i}),\\2\Delta _{i}-4{\frac {z_{i}-z_{\nu }^{*}}{r_{\nu }\mu _{\nu }}}&\nu =\nu (x_{i})>\nu (x_{i-1}),\end{cases}}}$

   где ${\displaystyle \Delta _{i}=x_{i}-x_{i-1}}$.

   Величины ${\displaystyle r_{\nu }>1,\;\nu ={\overline {1,\;m}}}$ являются параметрами алгоритма. От них зависят произведения ${\displaystyle r_{\nu }\mu _{\nu }}$, используемые при вычислении характеристик в качестве оценок неизвестных констант Липшица.

6. Определить интервал ${\displaystyle (x_{t-1};\;x_{t})}$, которому соответствует максимальная характеристика ${\displaystyle R(t)=\max\{R(i),\;1\leqslant i\leqslant k\}}$.
7. Провести очередное испытание в середине интервала ${\displaystyle (x_{t-1};\;x_{t})}$, если индексы его концевых точек не совпадают:

   ${\displaystyle x^{k+1}={\frac {1}{2}}(x_{t}+x_{t-1}).}$

   В противном случае провести испытание в точке

   ${\displaystyle x^{k+1}={\frac {1}{2}}(x_{t}+x_{t-1})-{\frac {z_{t}-z_{t-1}}{2r_{\nu }\mu _{\nu }}},\;\nu =\nu (x_{t})=\nu (x_{t-1}),}$

   увеличить ${\displaystyle k}$ на 1.

8. Если ${\displaystyle x_{t}-x_{t-1}<\varepsilon }$ (${\displaystyle \varepsilon >0}$ — заданная точность метода), то прекратить выполнение алгоритма, иначе перейти на шаг 1.

Пусть исходная задача оптимизации имеет решение ${\displaystyle x^{*}}$ и выполняются следующие условия:

1. каждая область ${\displaystyle Q_{j},\;j={\overline {1,\;m}}}$ представляет собой объединение конечного числа отрезков, имеющих положительную длину;
2. каждая функция ${\displaystyle g_{j}(x),\;j={\overline {1,\;m+1}}}$ удовлетворяет условию Липшица с соответствующей константой ${\displaystyle L_{j}}$;
3. компоненты вектора резервов удовлетворяют неравенствам ${\displaystyle 0\leqslant 2\varepsilon _{\nu }<L_{\nu }(\beta -\alpha )}$, где ${\displaystyle \beta -\alpha }$ — длина отрезка ${\displaystyle [\alpha ;\;\beta ]}$, лежащего в допустимой области ${\displaystyle Q}$ и содержащего точку ${\displaystyle x^{*}}$;
4. начиная с некоторого значения ${\displaystyle k}$ величины ${\displaystyle \mu _{\nu }}$, соответствующие непустым множествам ${\displaystyle I_{\nu }}$, удовлетворяют неравенствам ${\displaystyle r_{\nu }\mu _{\nu }>2L_{\nu }}$.

Тогда верно следующее:

1. точка ${\displaystyle x^{*}}$ является предельной точкой последовательности ${\displaystyle \{x^{k}\}}$, порождаемой методом при ${\displaystyle \varepsilon =0}$ в условии остановки;
2. любая предельная точка ${\displaystyle x^{0}}$ последовательности ${\displaystyle \{x^{k}\}}$ является решением исходной задачи оптимизации;
3. сходимость к предельной точке ${\displaystyle x^{0}}$ является двухсторонней, если ${\displaystyle x^{0}\neq a,\;x^{0}\neq b}$.

Общая схема последовательного метода выглядит следующим образом:

1. Упорядочить точки предшествующих испытаний в порядке возрастания их координат: ${\displaystyle a=x_{0}<\ldots <x_{i}<\ldots <x_{k}=b}$.
2. Вычислить для каждого интервала ${\displaystyle (x_{i-1};\;x_{i}),\;1\leqslant i\leqslant k}$ характеристику ${\displaystyle R(i)}$.
3. Определить интервал ${\displaystyle (x_{t-1};\;x_{t})}$, которому соответствует максимальная характеристика ${\displaystyle R(t)=\max\{R(i),\;1\leqslant i\leqslant k\}}$.
4. Провести следующее испытание в точке ${\displaystyle x^{k+1}=d(t)\in (x_{t-1};\;x_{t})}$, где ${\displaystyle d(t)}$ — правило размещения точки следующего испытания в интервале с номером ${\displaystyle t}$.
5. Проверить выполнение критерия остановки ${\displaystyle x_{t}-x_{t-1}<\varepsilon }$.

Параллельная модификация заключается в том, что на шаге 3 вместо одного интервала с наилучшей характеристикой выбирать ${\displaystyle p>1}$ интервалов в порядке убывания характеристик и параллельно проводить в каждом из них испытания.

Схема параллельного алгоритма:

1. Упорядочить точки предшествующих испытаний в порядке возрастания их координат: ${\displaystyle a=x_{0}<\ldots <x_{i}<\ldots <x_{k}=b}$.
2. Вычислить для каждого интервала ${\displaystyle (x_{i-1};\;x_{i}),\;1\leqslant i\leqslant k}$ характеристику ${\displaystyle R(i)}$.
3. Характеристики интервалов упорядочить по убыванию: ${\displaystyle R(i_{1})>\ldots >R(i_{k})}$.
4. Для всех интервалов с номерами ${\displaystyle i_{1},\;\ldots ,\;i_{p}}$ провести испытания в точках ${\displaystyle x^{k+j}=d(i_{j})\in (x_{i_{j}-1};\;x_{i_{j}}),\;j={\overline {1,\;p}}}$.
5. Проверить выполнение критерия остановки: ${\displaystyle \exists j,\;1\leqslant j\leqslant p\colon x_{i_{j}}-x_{i_{j}-1}<\varepsilon }$.

Такая схема распараллеливания целесообразна, если проведение испытания (то есть вычисление функций задачи) — трудоемкий процесс.

Метод достаточно просто обобщается на случай, когда функции ${\displaystyle g_{j}(x),\;j={\overline {1,\;m+1}}}$ удовлетворяют условию Гёльдера с показателем ${\displaystyle 1/n}$, где ${\displaystyle n\in \mathbb {N} }$, то есть

${\displaystyle |g_{j}(x+\Delta x)-g_{j}(x)|\leqslant H_{j}(\Delta x)^{1/n},\;a\leqslant x+\Delta x\leqslant b}$.

На шаге 3 значения ${\displaystyle \mu _{\nu }}$ вычисляются по формуле:

${\displaystyle \mu _{\nu }=\max\{|g_{\nu }(x_{i})-g_{\nu }(x_{j})|/(x_{i}-x_{j})^{1/n}\colon i,\;j\in I_{\nu },\;i>j\}.}$

На шаге 5 ${\displaystyle \Delta _{i}=(x_{i}-x_{i-1})^{1/n}}$.

На шаге 7 в случае совпадения индексов концевых точек

${\displaystyle x^{k+1}={\frac {1}{2}}(x_{t}+x_{t-1})-\operatorname {sgn}(z_{t}-z_{t-1}){\frac {|z_{t}-z_{t-1}|^{n}}{2r_{\nu }\mu _{\nu }^{n}}},\;\nu =\nu (x_{t})=\nu (x_{t-1}).}$

На шаге 8 критерий остановки принимает вид ${\displaystyle (x_{t}-x_{t-1})^{1/n}<\varepsilon }$.

• Параметры ${\displaystyle r_{\nu }}$ отвечают за надежность метода. Чем больше их значения, тем больше итераций метода требуется для достижения заданной точности и тем вероятнее выполнение условия сходимости 4. Если устремить все ${\displaystyle r_{\nu }}$ к бесконечности, то ${\displaystyle R(i)=\Delta _{i}}$, то есть метод превращается в перебор по равномерной сетке.
• Использование ненулевого вектора резервов позволяет ускорить сходимость метода, однако при этом необходимо оценить возможность выполнения условия сходимости 3.
• Одномерный метод может быть применен для решения многомерных задач без ограничений. Многомерная задача на множестве ${\displaystyle S=\{(x_{1},\;\ldots ,\;x_{n})\in \mathbb {R} ^{n}\colon a_{i}\leqslant x_{i}\leqslant b_{i},\;i={\overline {1,\;n}}\}}$ представляется в виде

${\displaystyle \min _{(x_{1},\;\ldots ,\;x_{n})\in S}f(x_{1},\;\ldots ,\;x_{n})=\min _{a_{1}\leqslant x_{1}\leqslant b_{1}}\min _{a_{2}\leqslant x_{2}\leqslant b_{2}}\ldots \min _{a_{n}\leqslant x_{n}\leqslant b_{n}}f(x_{1},\;\ldots ,\;x_{n}).}$

Для решения задачи ${\displaystyle \min _{a_{1}\leqslant x_{1}\leqslant b_{1}}\phi (x_{1})}$, где ${\displaystyle \phi (x_{1})=\min _{a_{2}\leqslant x_{2}\leqslant b_{2}}\ldots \min _{a_{n}\leqslant x_{n}\leqslant a_{n}}f(x_{1},\;\ldots ,\;x_{n})}$ можно использовать одномерный алгоритм, но, чтобы вычислить значение функции ${\displaystyle \phi (x_{1})}$, необходимо решить задачу оптимизации размерности ${\displaystyle n-1}$.

Если ${\displaystyle n=2}$, то задача ${\displaystyle \min _{a_{2}\leqslant x_{2}\leqslant b_{2}}f(x_{1},\;x_{2})}$ решается одномерным методом (значение переменной ${\displaystyle x_{1}}$ при этом зафиксировано), иначе к ней также применяется процедура снижения размерности. Такой способ решения многомерных задач довольно трудоемкий, поэтому на практике применим при ${\displaystyle n\leqslant 5}$. Наличие нелинейных функциональных ограничений может привести к потере липшицевости во вспомогательных одномерных задачах.

1. Баркалов К. А., СтронгинР. Г. Метод глобальной оптимизации с адаптивным порядком проверки ограничений // Ж. вычисл. матем. и матем. физ. — 2002. — Т. 42. — № 9. — стр. 1338—1350.
2. Городецкий С. Ю., Гришагин В. А. Нелинейное программирование и многоэкстремальная оптимизация. — Нижний Новгород: Издательство Нижегородского Университета, 2007.
3. Стронгин Р. Г. Численные методы в многоэкстремальных задачах (информационно-статистические алгоритмы). — М.: Наука, 1978. — 240 с.
4. Sergeyev Ya. D., Grishagin V. A. Sequential and parallel algorithms for global optimization // Optimization Methods and Software, 3:1-3, 1994, pp. 111—124.
5. Маркин Д. Л., Стронгин Р. Г. Метод решения многоэкстремальных задач с невыпуклыми ограничениями, использующий априорную информацию об оценках оптимума // Ж. вычисл. матем. и матем. физ., 27:1 (1987), стр. 56—62.

• [1] - реализация метода на языке C++.
• [2] - реализация на языке C++ модификации метода метода для решения многокритериальных многомерных задач.