Algorithms for optimizing systems with multiple extremum functionals

全文:

1. ВВЕДЕНИЕ

Для численного решения задач оптимизации сложных систем в частных производных представляется целесообразным использовать прямой экстремальный подход [1–3]. Суть похода заключается в прямой минимизации (максимизации) экстремальными алгоритмами целевого функционала

$\begin{array}{l}J\left(u\right)={\int }_{\omega }I\left(v\mathrm{,}u\right)d\mathrm{\omega }\to \min \mathrm{,}\mathrm{\omega }\subset \overline{\Omega }\mathrm{,}\\ \mathrm{при} \mathrm{условии} \mathbb{D}\left(\tau \mathrm{,}v\mathrm{,}u\right)v\mathrm{\tau }\in \overline{\Omega }\mathrm{,}\end{array}$ (1.1)

где функция управления $u\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo>}\in U\mathrm{<mo\; stretchy="false">(</mo>}S\mathrm{<mo\; stretchy="false">)</mo>}$, $S\subset \overline{\Omega }$, τ $–$ пространственно-временная переменная, U $–$ допустимое множество управлений, $\overline{\Omega }$ $–$ замкнутая область функционирования системы с состоянием $v\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo>}\in V\mathrm{<mo\; stretchy="false">(</mo>}\overline{\Omega }\mathrm{<mo\; stretchy="false">)</mo>,}V$ $–$ пространство (допустимое множество) состояний системы. Здесь оператор D, действующий на v, включает в себя не только конкретный вид дифференциальных уравнений на Ω, но и краевые условия где-либо на границе ∂Ω. Функция цели $I\mathrm{<mo\; stretchy="false">(</mo>}v\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo>,}u\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}$ определена на множестве ω, а ее значение зависит от параметра v и возможно u.

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

$u_{\mathrm{<mo>*</mo>}}\mathrm{<mo>=</mo>}\arg \min J\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>,}$ (1.2)

где $u_{\mathrm{<mo>*</mo>}}\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo>}$ $–$ оптимальное управление, оптимум, оптималь, оптимальный параметр, доставляющий глобальный минимум функционалу J(u). Задача (1.2) является бесконечномерной задачей оптимизации, поскольку управление u – функция. Для ее решения, обычно, используют методы градиентного спуска или методы сопряженных градиентов, обобщенные на бесконечномерные пространства. Сходимость этих методов за конечное число итераций обоснована только для конечномерных пространств управлений, когда u $–$ вектор, а J(u) при этом $–$ функция. Доказательства сходимости в таких задачах основываются, если не на квадратичности функции J(u), то на ее выпуклости. Для минимизации невыпуклых, многоэкстремальных функций существует большое количество алгоритмов (с той или иной степенью обоснования сходимости), но они не обобщаются на бесконечномерные задачи.

Традиционные градиентные методы, применяемые в бесконечномерных пространствах, имеют вид:

$u^{k+1}\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{u}^k\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo>}-b^k\nabla J^k\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo>;</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo>,}\tau \in S\mathrm{,}k\mathrm{<mo>=</mo>0,1}\dots ,$ (1.3)

где k $–$ номер итерации, $b^k$ $–$ шаговый множитель, регулирующий глубину спуска к минимуму вдоль антиградиента $-\nabla J^k$. Градиент вычисляется с использованием уравнений $\mathbb{D}\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{,}v\mathrm{,}u\mathrm{<mo\; stretchy="false">)</mo>}v\mathrm{<mo>=</mo>0}$. Все разнообразие методов (1.3) заключается в выборе шагового множителя $b^k$ [4–6] для ускорения сходимости к минимуму J(u).

Общепринятая стратегия выбора b^k $–$ это выбор заранее, возможно постоянного, шагового множителя на всех итерациях:

$b^k=b>0\text{\hspace{1em}или\hspace{1em}}{b}^k={\mathrm{\beta }}_1\left(k\right)b\mathrm{,}k=0,1\dots ,$ (1.4)

где ${\beta }_1\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}$ $–$ некоторая положительная, заранее известная функция. Данная стратегия нередко используется при минимизации выпуклых J(u) [5]. Для таких задач оптимизации стратегия является релаксационной, т.е. она ослабляет, уменьшает шаги так, что $‖b^k\nabla J^k‖\stackrel{k\to \infty }{\to }0$ по мере приближения к u_*. Если бы была известна аналитическая зависимость J(u), то можно было бы найти оптимальное b [6]. При этом используются знания о выпуклости J(u), такие как константа липшицевости градиентов, собственные числа матрицы Гессе. Понятно, что при бесконечномерной оптимизации получить что-либо подобное о поведении J(u) вряд ли возможно.

К достоинствам указанной стратегии следует отнести минимальность усилий в ее реализации и отсутствие дополнительных вычислений J или $\nabla J$ на каждой итерации k. Поэтому стратегия может оказаться эффективной в некоторых частных случаях при некоторых исходных данных.

Остается открытым вопрос, как исследовать, решать экстремальные задачи с невыпуклыми, многоэкстремальными целевыми функционалами? В настоящей работе предлагаются градиентные алгоритмы, основанные на специфическом выборе шагового множителя b^k, для исследования и off-line решения задач оптимизации (не непрерывного on-line управления) сложных систем в частных производных.

2. АЛГОРИТМЫ

Рассмотрим следующую (вторую) стратегию выбора b^k, которая использует адаптивную релаксацию шага:

$b^k\mathrm{<mo>=</mo>}{\beta }_2\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}{b}^{k-1}\mathrm{,}k\mathrm{<mo>=</mo>1.2,}\dots ,$

где β₂(k) $–$ положительная, заранее неизвестная функция, значение которой на каждой итерации определяется (адаптируется) на основании полученной ранее информации о поведении J(u). Например, для минимизации выпуклых J(u), это может быть градиентный метод вида [1]:

$\begin{array}{l}\text{Если\hspace{0.33em}}{J}^k<J^{k-1},\text{\hspace{0.33em}тогда:}{b}^k=b_1{b}^{k-1}\mathrm{,}{u}^{k+1}=u^k-b^k\nabla J^k\mathrm{.}\\ \text{Иначе повторять до }{J}^k<J^{k-1}:\\ b^{k-1}\leftarrow b_2{b}^{k-1}\mathrm{,}{u}^{k+1}=u^{k-1}-b^{k-1}\nabla J^{k-1}\mathrm{,}\\ \text{\hspace{1em}если\hspace{0.33em}}{b}^{k-1}\approx \mathrm{0,}\text{\hspace{0.33em}тогда стоп итерации\hspace{0.33em}}k\mathrm{.}\end{array}$ (2.1)

Здесь $k\mathrm{<mo>=</mo>1,2}\dots \mathrm{,}{b}_1\ge \mathrm{1,}{b}_2\mathrm{<mo><</mo>1}$. Начальное значение b⁰ задается из условия J¹ < J⁰.

Алгоритм (2.1) посредством параметра b₁ может увеличивать шаги и тем самым усиливать сходимость к минимуму, а параметром b₂ $–$ контролировать и не допускать расходимость метода из-за чрезмерно больших шагов. Условие завершения «стоп» в цикле повторений предыдущего шага контролирует зацикливание и излишнее количество чрезмерно малых шагов при уменьшении b^k–1. Если b^k–1 ≈ 0, то следует считать, что процесс минимизации J(u) достиг своего предела. Итерации k прекращаются. Очевидно, алгоритм (2.1), в общем случае, для минимизации функционалов, эффективнее алгоритмов описанной ранее в первой стратегии, хотя он требует дополнительного хранения в памяти массивов $u^{k-1}\mathrm{,<mi>\; </mi>}\nabla J^{k-1}$.

Если $J\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}$ не выпуклый, но имеет единственный экстремум $–$ минимум, то в алгоритме минимизации (2.1) следует предусмотреть и возможность ослабления сходимости, если обнаруживается рост нормы градиента (здесь и далее норма вычисляется в $L_2$ ). Это можно сделать следующим градиентным методом:

$\begin{array}{l}\text{Если\hspace{0.33em}}{J}^k<J^{k-1},\text{\hspace{0.33em}тогда:}\\ \cdot \text{\hspace{0.33em}если\hspace{0.33em}}||\nabla J^k||\text{<}||\nabla J^{k-1}||\text{,}\text{\hspace{0.33em}тогда\hspace{0.33em}}{b}^k=b_1{b}^{k-1}\mathrm{,}\\ \text{\hspace{0.33em}иначе\hspace{0.33em}}{b}^k=\frac{||\nabla J^{k-1}||}{b_3||\nabla J^k||}{b}^{k-1}\\ \cdot u^{k+1}=u^k-b^k\nabla J^k\mathrm{.}\\ \text{Иначе повторять до\hspace{0.17em}}{J}^k<J^{k-1}\\ b^{k-1}\leftarrow b_2{b}^{k-1}\mathrm{,}{u}^{k+1}=u^{k-1}-b^{k-1}\nabla J^{k-1}\mathrm{,}\\ \text{\hspace{1em}если\hspace{0.33em}}{b}^{k-1}\approx \mathrm{0,}\text{\hspace{0.33em}тогда стоп итерации\hspace{0.33em}}k\mathrm{.}\end{array}$  (2.2)

Здесь b₃≥1. Если J(u) убывает, а норма градиента не убывает (возрастает), то это означает, что текущий шаг либо пришелся на вогнутую часть J(u), либо перескочил через локальный минимум и попал в область повышенной выпуклости J(u). В любом случае следующий шаг $b^k\nabla J^k$ не должен быть больше предыдущего. Величина такого шага регулируется параметром b₃. Надо быть внимательным при выборе b₃>1, поскольку можно существенно замедлить сходимость алгоритма.

Если J(u) не выпуклый и при этом имеет несколько локальных экстремумов, то алгоритм (2.2) следует модифицировать, посредством адаптации параметра b₁. Для того, чтобы сходимость не завершилась в каком-либо физически неудовлетворительном локальном экстремуме (или хочется проверить существуют ли другие экстремумы), следует приближаться к нему с большим значением b₁ $–$ для перешагивания через экстремум. И наоборот, чтобы не перешагнуть через желаемый экстремум, следует приближаться к нему осторожно, с уменьшенным b₁.

Алгоритм адаптации b₁ существенным образом зависит от конкретной задачи оптимизации. Данная процедура требует кропотливого участия исследователя и вряд ли может быть формализована. Тем не менее, она позволяет исследовать многоэкстремальные задачи бесконечномерной оптимизации.

Обсуждаемый алгоритм выбора шагового множителя имеет вид:

$\begin{array}{l}\text{Если\hspace{0.33em}}{J}^k<J^{k-1},\text{\hspace{0.33em}тогда:}\\ \cdot \text{\hspace{0.33em}если}||\nabla J^k||\text{<}||\nabla J^{k-1}||,\text{\hspace{0.33em}тогда}\\ \left(\text{если\hspace{0.33em}}{u}^k\text{рядом с экстремумом } -\text{ изменить }{b}_1\right)b^k=b_1{b}^{k-1}\mathrm{,}\\ \text{\hspace{0.33em}иначе\hspace{0.33em}}{b}^k=\frac{||\nabla J^{k-1}||}{b_3||\nabla J^k||}{b}^{k-1};\\ \cdot u^{k+1}=u^k-b^k\nabla J^k\mathrm{.}\\ \text{Иначе повторять до\hspace{0.33em}}{J}^k<J^{k-1}:\\ b^{k-1}\leftarrow b_2{b}^{k-1}\mathrm{,}{u}^{k+1}=u^{k-1}-b^{k-1}\nabla J^{k-1}\mathrm{,}\\ \text{\hspace{1em}если\hspace{0.33em}}{b}^{k-1}\approx \mathrm{0,}\text{\hspace{0.33em}тогда стоп итерации\hspace{0.33em}}k\mathrm{.}\end{array}$ (2.3)

Существуют и другие стратегии поиска b^k для минимизации функционалов. Например, стратегии линейного поиска на каждой итерации, которая осуществляется в два этапа. Сначала ищется (задается) отрезок, в направлении $-\nabla J^k$, потом на этом отрезке ищется удовлетворительное b^k. Здесь самым известным представителем является градиентный метод наискорейшего спуска – это стратегия полной релаксации, когда на каждом шаге вдоль направления $-\nabla J^k$, на отрезке заданной длины, выбирается оптимальное значение множителя:

$b^k\mathrm{<mo>=</mo>}\arg \underset{b\mathrm{<mo>></mo>0}}{\min }J\left(u^k-b\nabla J^k\right)\mathrm{.}$

Это замечательная стратегия. Однако, если управление ограничено допустимым множеством или другими специфическими особенностями задачи, то выбор отрезка, содержащего минимум функции J^k(b), может стать невозможным. Любой «конец» отрезка, начавшегося в точке u^k, может выйти за ограничения и сделать абсурдным решение уравнений $\mathbb{D}\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{,}v\mathrm{,}u\mathrm{<mo\; stretchy="false">)</mo>}v\mathrm{<mo>=</mo>0}$.

Для пояснения и иллюстрации актуальности и работоспособности предложенных алгоритмов второй стратегии рассмотрим следующий пример.

3. ПРИМЕР

3.1. Постановка задачи

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

Фиг. 1.Схема поршневой гидропушки

Движение воды в цилиндрическом сопле можно описать следующей квазиодномерной, квазилинейной гиперболической системой уравнений [7]:

$\begin{array}{c}\frac{\partial \rho }{\partial t}+w\frac{\partial \rho }{\partial x}+\rho \frac{\partial w}{\partial x}+\frac{\rho w}{\sigma }\frac{d\sigma }{dx}\mathrm{<mo>=</mo>0,}\\ \frac{\partial w}{\partial t}+\frac{Bn{\rho }^{n-2}}{{\rho }_0^n}\frac{\partial \rho }{\partial x}+w\frac{\partial w}{\partial x}\mathrm{<mo>=</mo>0.}\end{array}$ (3.1)

Состояние системы характеризуется вектор-функцией $v\mathrm{<mo>=</mo><mo>\{</mo>}\rho \mathrm{,}w\mathrm{<mo>\}</mo>}$, где $\rho$ $–$ плотность потока, w $–$ скорость. Состояние определено на пространственновременной области $\overline{\Omega }$ втекания воды в сопло и истечения из него. B и n $–$ постоянные в уравнении состояния Тэта. Управление

$u\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\frac{1}{\sigma \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}}\frac{d\sigma \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}}{dx}\mathrm{,}x\in \left(x_a\mathrm{,}{x}_b\right)\mathrm{,}$ (3.2)

где $\sigma$ $–$ площадь поперечного сечения сопла, $\sigma \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\sigma }_a\text{exp}\left({\displaystyle {\int }_{x_a}^x}u\mathrm{<mo\; stretchy="false">(</mo>}\zeta \mathrm{<mo\; stretchy="false">)</mo>}d\zeta \right)\mathrm{,}x\in \left[x_a\mathrm{,}{x}_b\right]$, ${\sigma }_a$ $–$ площадь ствола гидропушки.

Система (3.1) $–$ это значение оператора $\mathbb{D}v$, которое дополняется краевыми условиями. Слева – это уравнение движения поршня массой $m_p$ : $\frac{dw}{dt}+\frac{{\sigma }_{a}B}{m_p}\left({\left(\frac{\rho }{{\rho }_0}\right)}^n-1\right)\mathrm{<mo>=</mo>0}$. На правой границе, взаимодействующей с атмосферой, плотность воды $\mathrm{\rho }={\mathrm{\rho }}_0={10}^3\text{кг}/{\mathrm{м}}^3$. Начальные условия – это $w_0=\text{\hspace{0.33em}76 м}/\text{с}\mathrm{,}\mathrm{\rho }={\mathrm{\rho }}_0$. Радиус ствола $R_a\mathrm{<mo>=</mo>}$ $33\cdot {10}^{-3}\text{м}$, длина сопла $0.253\text{\hspace{0.33em}м}$, начальная длина столба вода 0.28 м.

Определим целевой функционал. Будем максимизировать среднюю силу действия импульса струи на преграду:

$J\mathrm{<mo>=</mo>}{\left.\frac{1}{t_2-t_1}{\displaystyle {\int }_{t_1}^{t_2}}\rho \sigma w^{2}dt\right|}_{x_b}\mathrm{,}$

где $t_1$ $–$ начало истечения струи из сопла, $t_2$ $–$ конечное время формирования струи. В частности, $t_1\approx 2.7\cdot {10}^{-3}\text{с}$, а время истечения $t_2-t_1=3\cdot {10}^{-4}\text{с}$.

Подобная задача рассматривалась и ранее [8, 9] в рамках подхода классического вариационного исчисления. Были предложены различные необходимые условия оптимальности, но оптимальная форма сопла так и не была получена. Только применение экстремального подхода с анализом управляемости и развитие адаптивных градиентных алгоритмов позволило решить эту сложную задачу оптимизации.

Далее будем обращаться к функционалу $J$ в форме:

$J\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle {\int }_{\omega }}I\mathrm{<mo\; stretchy="false">(</mo>}w\mathrm{,}u\mathrm{<mo\; stretchy="false">)</mo>}dt\to \max \mathrm{,}\omega \mathrm{<mo>=</mo>}{x}_b\times \left(t_1\mathrm{,}{t}_2\right)\mathrm{,}$ (3.3)

где подынтегральная функция цели

${\left.I\mathrm{<mo\; stretchy="false">(</mo>}w\mathrm{,}u\mathrm{<mo\; stretchy="false">)</mo>}\right|}_{\omega }\mathrm{<mo>=</mo>}\frac{{\rho }_0{\sigma }_b{w}_b^2}{t_2-t_1}\mathrm{<mo>=</mo>}\frac{{\rho }_0{\sigma }_a\text{exp}\left({\displaystyle {\int }_{x_a}^{x_b}}u\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}dx\right)w_b^2}{t_2-t_1}\mathrm{.}$

Уточним область определения управления $u\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}$ $–$ это множество $S\mathrm{<mo>=</mo>}$ $\left(x_a\mathrm{,}{x}_b\right)$, при этом форма сопла $\sigma \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}$ определена на $\left[x_a\mathrm{,}{x}_b\right]$. Областью значений управления будем считать полупространство допустимых значений:

$U\mathrm{<mo\; stretchy="false">(</mo>}S\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\left\{u\mathrm{<mo>:</mo>}u\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}\le 0\forall x\in \left(x_a\mathrm{,}{x}_b\right)\right\}$. (3.4)

Такому $U$ соответствует $\sigma \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}\le {\sigma }_a$, что является физически разумным. Необходимо иметь в виду и еще одно ограничение. При сверхзвуковом истечении, когда $w\ge c_0=1475\text{\hspace{0.33em}м}/\text{с}$, струя разрушается. Это ограничение не входит в постановку задачи оптимизации, а является индикатором возможной неудачной оптимизации.

Для прямой максимизации целевого функционала (10) необходимо найти его градиент $\nabla J$. Техника аналитического определения градиента неявно заданного функционала $–$ это самостоятельная сложная задача. Здесь можно использовать подход [1], что было нами сделано в работе [10]. Мы получим следующее выражение градиента:

$\nabla J\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo>;</mo>}x\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle {\int }_{t_1}^{t_2}}\rho w{f}_{1}dt+\kappa J\mathrm{,}x\in \left(x_a\mathrm{,}{x}_b\right)\mathrm{,}$

где f₁ $–$ это компонента сопряженного, по отношению к v, состояния системы, $\kappa \mathrm{<mo>=</mo>0.17}$ $–$ весовой коэффициент выравнивания вычислительных помех решения на $\overline{\Omega }$ исходной нелинейной задачи $\mathbb{D}\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{,}v\mathrm{,}u\mathrm{<mo\; stretchy="false">)</mo>}v\mathrm{<mo>=</mo>0}$ и линейной сопряженной задачи на $\left[x_a\mathrm{,}{x}_b\right]\times \left[t_1\mathrm{,}{t}_2\right]\subset \overline{\Omega }$:

$\begin{array}{c}-\frac{\partial f_1}{\partial t}-w\frac{\partial f_1}{\partial x}-\frac{Bn{\rho }^{n-2}}{{\rho }_0^n}\frac{\partial f_2}{\partial x}+wu{f}_1\mathrm{<mo>=</mo>0,}\\ -\frac{\partial f_2}{\partial t}-\rho \frac{\partial f_1}{\partial x}-w\frac{\partial f_2}{\partial x}+\rho u{f}_1\mathrm{<mo>=</mo>0.}\end{array}$

На левой границе f₁ = 0. На правой, при истечении с t₁ и до t₂ имеем $f_1+\frac{w}{{\rho }_0}{f}_2+\frac{2{\sigma }_b{w}_b}{t_2-t_1}\mathrm{<mo>=</mo>0}$. Начальное (терминальное) условие при t₂ $–$ это f₁ = 0, f₂ = 0. Задача решается в обратном по времени направлении, начиная с нулевого терминального состояния.

3.2. Оценка выпуклости целевого функционала

Для получения представления о выпуклости функционала (3.3), рассмотрим управление в классе функций-конусов. При этом $J\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}J\left(R_b\right)$, т.е. J будет функцией радиуса среза сопла R_b, а $\sigma \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\pi {\left(\frac{R_b-R_a}{x_b-x_a}x+R_a\right)}^2\mathrm{,}{R}_a\mathrm{<mo>=</mo>}\sqrt{{\sigma }_a\mathrm{<mo>/</mo>}\pi }$. Полученная функция J представлена на фиг. 2.

Фиг. 2.Зависимость функционала-функцииJ(R_b)

Правая точка на графике соответствует соплу в виде трубы, R_b = R_a. Левая точка соответствует минимально возможному радиусу среза сопла, при котором еще не наступает сверхзвуковое истечение. Рядом с этой точкой имеется maxJ, с большой вогнутостью и очень малой окрестностью дозвукового течения. Функционал, в процессе изменения радиуса сопла от трубы до минимально допустимого сужения, меняет совою начальную выпуклость на вогнутость. Причем, вогнутость на много сильнее начальной выпуклости.

Если максимизировать функцию J(R_b), начиная с R_b = R_a, классическим конечномерным градиентным методом с первой стратегией выбора b = b⁰, то небольшие начальные шаги $b\nabla J^k$ при движении влево сменятся на большие шаги (из-за возрастания градиента). Поскольку вогнутость имеет место в относительно малой окрестности максимума, то большие шаги могут приводить к перепрыгиванию через экстремум вплоть до появления сверхзвукового истечения. Это недопустимо. Попасть в maxJ, возможно, удастся только с очень малым b и, естественно, при чрезмерно большом количестве итераций. Применение градиентного метода наискорейшего спуска, метода сопряженных градиентов и других методов, использующих линейный поиск наилучшего b^k, здесь невозможно, поскольку такой поиск, в общем случае, будет начинаться со сверхзвукового истечения.

В этой ситуации, для достижения maxJ, целесообразно использовать метод (6) с адаптивной стратегией выбора b^k. Давайте с него и начнем реальную оптимизацию произвольной формы сопла. Решение задачи программировалось автором в среде Delphi 7, вычисления проводились на компьютере с индексом производительности 3.5 Windows на пространственновременной сетке 40 × 500 в течение времени менее минуты.

3.3. Первый локальный максимум

Пусть начальное приближение

$u^0\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>0,}$

что соответствует трубе ${\sigma }^0\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\pi R_a{}^2$. Начальное значение шагового множителя

$b^0\mathrm{<mo>=</mo>}\frac{0.5}{‖\nabla J^0‖}\mathrm{,}$

что соответствует первому шагу $‖u^1-u^0‖\mathrm{<mo>=</mo>0.5}$. Сначала рассмотрим задачу оптимизации без условия (3.4) контроля расширения сопла.

Применим градиентный метод (2.2) для максимизации J(u):

 $\begin{array}{l}\text{Если\hspace{0.33em}}{J}^k>J^{k-1},\text{\hspace{0.33em}тогда:}\\ \cdot \text{\hspace{0.33em}если\hspace{0.33em}}||\nabla J^k||\text{<}||\nabla J^{k-1}||,\text{\hspace{0.33em}тогда\hspace{0.33em}}{b}^k=b_1{b}^{k-1}\mathrm{,}\\ \text{\hspace{0.33em}иначе\hspace{0.33em}}{b}^k=\frac{||\nabla J^{k-1}||}{b_3||\nabla J^k||}{b}^{k-1};\\ \cdot u^{k+1}=u^k+b^k\nabla J^k\mathrm{.}\\ \text{Иначе\hspace{0.33em}повтор}\stackrel{⏜}{\text{я}}\text{ть\hspace{0.33em}до\hspace{0.33em}}{J}^k>J^{k-1}:\\ b^{k-1}\leftarrow b_2{b}^{k-1}\mathrm{,}{u}^{k+1}=u^{k-1}+b^{k-1}\nabla J^{k-1}\mathrm{,}\\ \text{\hspace{1em}если\hspace{0.33em}}{b}^{k-1}\approx \mathrm{0,}\text{\hspace{0.33em}тогда\hspace{0.33em}стоп\hspace{0.33em}итерации\hspace{0.33em}}k\mathrm{.}\end{array}$  (3.5)

Были подобраны значения параметров метода: b₁ = 1.05, b₂ = 0.3, b₃ = 1.05. Подбор осуществлялся из соображений: b₁ должен быть небольшим для «осторожного» подхода к ближайшему сильно вогнутому экстремуму; b₂ во всех расчетах удовлетворительно выполняло свою задачу при 0.3; b₃ $–$ чуть больше единице, иначе, при больших значениях сходимость сильно замедлялась, а при b₃ = 1 происходил слишком большой переход (из-за резкой смены выпуклости) через первый экстремум. Остановка итераций осуществлялась при $‖u^k-u^{k-1}‖\mathrm{<mo>/</mo>}‖u^{k-1}‖\le {10}^{-6}$.

Результаты оптимизации представлены на фиг. 3 в виде радиуса сопла для разных площадей σ. Оптимальной формой сопла оказался практически конус со значением целевого функционала J = 1.26·10⁵.

Фиг. 3.Оптимальное сопло в первом экстремуме

Традиционным бесконечномерным градиентным методом (1.3), с постоянным шаговым множителем, не удалось попасть в рассматриваемый локальный maxJ. При $b\gtrsim 10b^0$ перешагивание экстремума сопровождалось сверхзвуковым истечением. При уменьшении b до 10^–2b⁰ сходимость не заканчивалась даже поле нескольких тысячах итераций, наблюдалось неконтролируемое расширение сопла после перехода локального maxJ.

3.4. Второй локальный максимум

Выясним, является ли первый максимум целевого функционала глобальным. Для этого необходимо посредством метода (3.5) грубо перейти и через максимум, и через близлежащий за ним минимум, чтобы опять выполнилось условие роста функционала J^k > J^k–1. Увеличим в окрестности максимума шаги посредством параметра b₁ от 1.05 до b₁ = 1.16. При этом был реализован переход через локальные максимум и минимум в область роста $J$, но, если не применять ограничение (3.4), то далее происходило неестественное расширение сопла за пределы ствола гидропушки, что приводит к аварийным завершениям расчета состояния потока. Ниже на фиг. 4 точечными линиями показано промежуточное расширение сопла для итерации k = 41.

Ограничение (3.4) легко реализуется проецированием управления на допустимое множество U. После шага любым алгоритмом адаптивной релаксации полученное управление u^k+1 дополнительно корректируется:

$\text{если\hspace{0.33em}}{u}^{k+1}\left(x\right)>0,\text{\hspace{0.33em}тогда\hspace{0.33em}}{u}^{k+1}\left(x\right)=0,x\in \left(x_a\mathrm{,}{x}_b\right)\mathrm{.}$

Дальнейшие расчеты оптимизации ограниченного сопла методом (3.5) и традиционным методом (1.3), при любых b⁰, приводили к сверхзвуковому истечению. Для (3.5) этот «неудачный результат» означает, что параметр увеличения шагов b₁ был недопустимо большим, и возможный второй локальный maxJ был грубо пройден. Здесь вместо метода (3.5) необходимо применять метод типа (2.3) и в окрестности ожидаемого maxJ уменьшать b₁.

Окрестность предполагаемого максимума можно определить по скорости истечения $w_b\mathrm{<mo>=</mo>}w\left(x_b\mathrm{,}t\right)$, $t\in \left(t_1\mathrm{,}{t}_2\right)$. Она должна быть близкой к скорости звука c₀. В частности, такой скоростью была принята $w_b>800\mathrm{м}/\mathrm{с}$. При ее превышении усиление шагов было отменено, т.е. b₁ = 1.

Метод (2.3), для максимизации целевого функционала J, имеет вид:

$\begin{array}{l}\text{Если\hspace{0.33em}}{J}^k>J^{k-1},\text{\hspace{0.33em}тогда:}\\ \cdot \text{\hspace{0.33em}если\hspace{0.33em}}||\nabla J^k||\text{<}||\nabla J^{k-1}||,\text{\hspace{0.33em}тогда}\\ \left(\text{если\hspace{0.33em}}{w}_b>800,\text{\hspace{0.33em}тогда\hspace{0.33em}}{b}_1=1.0\right)b^k=b_1{b}^{k-1}\mathrm{,}\\ \text{\hspace{0.17em}иначе\hspace{0.33em}}{b}^k=\frac{b^{k-1}||\nabla J^{k-1}||}{b_3||\nabla J^k||};\\ \cdot u^{k+1}=u^k+b^k\nabla J^k\mathrm{.}\\ \text{Иначе\hspace{0.33em}повторять\hspace{0.33em}до\hspace{0.33em}}{J}^k>J^{k-1}:\\ b^{k-1}\leftarrow b_2{b}^{k-1}\mathrm{,}{u}^{k+1}=u^{k-1}+b^{k-1}\nabla J^{k-1}\mathrm{,}\\ \text{\hspace{1em}если\hspace{0.33em}}{b}^{k-1}\approx \mathrm{0,}\text{\hspace{0.33em}тогда\hspace{0.33em}стоп\hspace{0.33em}итерации\hspace{0.33em}}k\mathrm{.}\end{array}$ (3.6)

Параметры метода были теми же: начальное $b_1\mathrm{<mo>=</mo>1.16,}{b}_2\mathrm{<mo>=</mo>0.3,}{b}_3\mathrm{<mo>=</mo>1.05}$.

На фиг. 4 сплошными линиями показан радиус сопла полученной оптимальной площади σ⁶⁸ на последней итерации k=68. Блок «Иначе повторять» в (3.6) срабатывал только в конце итераций при k≥66, что потребовало дополнительно 12 вычислений J. Значение целевого функционала составило J = 3.75 · 10⁵, т.е. в 3 раза больше (лучше), чем в первом локальном максимуме.

Фиг. 4.Оптимальное сопло во втором экстремуме

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

4. ВЫВОДЫ

Предложенные градиентные методы минимизации (максимизации) целевых функционалов, с шаговыми множителями адаптивной релаксации шагов при приближении к экстремуму, позволили решить сложную многоэкстремальную задачу оптимизации формы сопла гидропушки. Были выявлены два локальных максима у целевого функционала, характеризующего среднюю силу действия импульса струи на преграду. Второй максимум имел значение целевого функционала в 3 раза больше (лучше), чем первый. Полученное сопло, соответствующее глобальному максимуму, следует считать оптимальным. Применение традиционного градиентного метода не позволило найти ни одного локального экстремума целевого функционала. То есть, решить поставленную задачу градиентным методом с постоянным шаговым множителем не удалось, и по всей вероятности $–$ невозможно.

作者简介

V. Tolstykh

编辑信件的主要联系方式.
Email: mail@tolstykh.com
俄罗斯联邦, ul. Universitetskaya, 24, Donetsk, 283001

参考

补充文件

附件文件

动作

1. JATS XML

下载

2. Fig. 1. Piston hydraulic gun diagram

下载 (50KB)

索引源数据

3. Fig. 2. Dependence of the functional-function J(Rb)

下载 (47KB)

索引源数据

4. Fig. 3. Optimal nozzle at the first extremum

下载 (51KB)

索引源数据

5. Fig. 4. Optimal nozzle at the second extremum

下载 (94KB)

索引源数据