Проблемы машиностроения и надежности машин, 2019, № 7, стр. 20-29

1. Современные подходы к исследованию сложных аэрокосмических систем основаны на применении методов математического моделирования [1, 2]. К этому направлению относятся работы сотрудников одной из ведущих отечественных научных школ, которую возглавлял академик РАН К.С. Колесников (1919–2016 гг.) [3]. Оптимальное проектирование, включая междисциплинарную проблематику, предъявляет высокие требования к характеристикам моделей. Прямое моделирование реальных объектов не обеспечивает требуемое качество соответствующих аналитических моделей. Актуальным направлением исследований является разработка эффективных численных методов коррекции математических моделей аэрокосмических систем с использованием модальных данных, полученных при измерениях [4, 5]. Процедуры коррекции существенно связаны с формулировкой и последующим решением обратной задачи на собственные значения [6]. Разработаны достаточно эффективные методики коррекции параметров моделей исследуемых систем. При реализации численных процедур коррекции необходимо учитывать корректность постановки обратной задачи, неполноту и зашумленность измеряемой информации, возможное наличие в спектрах систем кратных частот и др. Во многих приложениях возникают обратные задачи, которые относятся к классу некорректно поставленных, что предполагает применение специальных методов регуляризации [7, 8]. Критериальные функции обратных задач в общем случае являются многомерными, непрерывными, многоэкстремальными и не всюду дифференцируемыми. Поэтому для решения задач коррекции моделей систем и содержательно близких к ним задач вычислительной диагностики (обнаружения повреждений) систем применяют методы глобальной оптимизации [9]. При определении каждого текущего значения критериальной функции в точках допустимой области необходимы значительные вычислительные ресурсы. Этим обусловлена актуальность направления, связанного с созданием эффективных алгоритмов коррекции моделей систем по измеряемым модальным данным. Рассматриваемый подход основан на разработке и применении математических моделей систем, численных методов определения основных динамических характеристик систем, методов теории обратных задач, методов глобальной недифференцируемой оптимизации.

2. Рассмотрим упругие колебания жидкостной ракеты (блока) с топливными баками, частично заполненными жидкостью. Актуальным приближением аналитической модели блока является прямой неоднородный стержень с упруго подвешенными массами [10] (рис. 1).

Продольная ось $Ox$ блока направлена от его вершины к хвостовой части; две другие оси образуют вместе с продольной правую прямоугольную декартову систему координат $Oxyz$. Полная длина блока $l = 25$ м; радиус каждого топливного бака равен $r = 0.75$ м. Исходными данными для вычислений служат функции распределения массы $m(x)$ и изгибной жесткости $EI(x)$ по длине стержня; координата $x$ отсчитывается от левого конца стержня вдоль его осевой линии. Стержень разбит на семь участков. Используются масштабы массы и изгибной жесткости стержня: $m(x) = \bar {m} \cdot {{m}_{0}}$, $EI(x) = \overline {EI} \cdot E{{I}_{0}}$, где ${{m}_{0}} = 147.15$ кг и $E{{I}_{0}} = 103.005 \cdot {{10}^{6}}$ Н м² [10]. Каждая упруго подвешенная масса (осциллятор) является механическим аналогом одного (низшего) тона колебаний жидкости в топливном баке: ${{l}_{j}}$ – расстояние от вершины блока до невозмущенной свободной поверхности жидкости в $j$-м баке; $l_{j}^{*}$ – параметр длины, используемый для определения объема (массы) колеблющейся жидкости; ${{x}_{j}}$ – координата точки подвеса осциллятора; $x_{j}^{*}$ – расстояние от невозмущенной свободной поверхности жидкости до точки подвеса осциллятора; ${{k}_{j}}$ – коэффициент жесткости $j$-го осциллятора; ${{m}_{j}}$ – масса осциллятора; $j = 1,{\text{ }}2$. Значения параметров участков неоднородного стержня, включая длины ${{l}_{i}}$ участков ($i$ – номер участка), приведены в табл. 1 [10].

В верхнем баке ($i = 4$) находится окислитель плотности 11180 кг/м³; в нижнем баке ($i = 6$) – горючее плотности 7848 кг/м³. Сосредоточенные массы ${{m}_{0}}\bar {m}_{j}^{*} = {{m}_{j}}$ осцилляторов расположены в узле 4* между участками 4 и 5 (667.08 кг), а также в узле 7* между участками 7 и 8 (470.88 кг).

Коэффициент жесткости каждого осциллятора ${{k}_{j}}$ выбирают так, чтобы частота собственных колебаний массы ${{m}_{j}}$ была равна частоте собственных колебаний жидкости в неподвижном баке: ${{k}_{1}}$ = 0.16544 ⋅ 10^–6 Н/м; ${{k}_{2}}$ = 0.11678 ⋅ 10^–6 Н/м.

С целью исследования пространственных упругих колебаний блока разработана его конечно-элементная модель. Блок представлен ансамблем одномерных конечных элементов, соединенных между собой в узлах; длина каждого элемента равна длине моделируемого участка разбиения. Инерционные и упругие характеристики конечно-элементной модели для исследования продольных и крутильных колебаний блока получены эквивалентным пересчетом данных, представленных в табл. 1, с учетом геометрических параметров блока. Подвижность жидкости в топливных баках учитывается только при изгибных колебаниях блока в плоскостях $Oxy$, $Oxz$ и моделируется системой осцилляторов.

Обобщенная задача на собственные значения для упругой пространственной системы, дискретизированной с использованием метода конечных элементов, имеет вид [9]

где $M,K$ – соответственно глобальные матрицы инерции (масс) и жесткости системы; $\lambda $ – собственное значение, $\lambda = {{\omega }^{2}}$; $\omega - $ круговая частота собственных колебаний; $\psi $ – собственный вектор, $\psi \in {{{\text{R}}}^{m}}$, ${{{\text{R}}}^{m}}{\text{ }} - {\text{ }}m$-мерное евклидово пространство; собственный вектор нормируется по условию ${{\psi }^{T}}M\psi = 1$, $T$ – символ транспонирования. При исследовании собственных колебаний упругих систем с большим числом степеней свободы достаточно учитывать ограниченное число низших спектральных составляющих, поэтому для численного решения задачи (1) целесообразно использование метода итерации подпространства [11]. Основной целью метода является определение ${{n}_{p}}$ низших собственных значений и соответствующих им собственных векторов (собственных пар), которые удовлетворяют матричному уравнению здесь $\Lambda = {\text{diag}}({{\lambda }_{i}})$, т.е. диагональная матрица размера ${{n}_{p}} \times {{n}_{p}}$, содержащая собственные значения; $\Psi = [{{\psi }^{1}},...,{{\psi }^{p}}]$. Собственные векторы должны также удовлетворять условиям ортогональности ${{\Psi }^{T}}K\Psi = \Lambda $; ${{\Psi }^{T}}M\Psi = I$; где I– единичная матрица порядка ${{n}_{p}}$. Уравнение (2) формулирует необходимые и достаточные условия того, что векторы, содержащиеся в матрице $\Psi $, являются собственными векторами.

Получено численное решение обобщенной задачи на собственные значения (1) (прямой задачи) для конечно-элементной модели блока. Определены – две низшие частоты изгибных колебаний блока в плоскостях $Oxy$ и $Oxz$ с учетом подвижности жидкости, имеющие одинаковое значение 4.0335 Гц. Сравнение со значением аналогичной составляющей спектра частот блока, вычисленным методом начальных параметров и равным 4.036 Гц [10], дает оценку относительной погрешности в 0.062%.

Рассмотрим аналитическую модель ракеты пакетной схемы, составленной из двух одинаковых блоков, характеристики которых приведены выше. Введем узлы межблочных связей, моделируемые пружинными элементами, размещенными в поперечных сечениях, соответствующих стыкам четвертого и пятого, а также седьмого и восьмого участков каждого блока. Верхний и нижний узлы связей имеют номинальные жесткостные характеристики (в системе координат $Oxyz$): $C_{x}^{{\text{в}}}$ = 0; $C_{y}^{{\text{в}}} = C_{z}^{{\text{в}}}$ = 0.8 ⋅ 10^–7 Н/м; $C_{x}^{{\text{н}}} = C_{y}^{{\text{н}}} = C_{z}^{{\text{н}}}$ = 0.8 ⋅ 10^–7 Н/м. Предположим, что значения характеристик могут изменяться, причем их целесообразно задать в виде: $c_{y}^{{\text{в}}} = C_{y}^{{\text{в}}} \cdot {{x}_{1}}{\text{/}}100$; $c_{z}^{{\text{в}}} = C_{z}^{{\text{в}}} \cdot {{x}_{2}}{\text{/}}100$; $c_{x}^{{\text{н}}} = C_{x}^{{\text{н}}} \cdot {{x}_{3}}{\text{/}}100$; $c_{y}^{{\text{н}}} = C_{y}^{{\text{н}}} \cdot {{x}_{4}}{\text{/}}100$; $c_{z}^{{\text{н}}} = C_{z}^{{\text{н}}} \cdot {{x}_{5}}{\text{/}}100$; здесь ${{x}_{i}}$ – переменные модели, определяющие текущие значения характеристик межблочных связей: $0 \leqslant {{x}_{i}} \leqslant 100\% $, $i = \overline {1.{\text{ }}5} $.

Анализ результатов численного решения прямой задачи для конечно-элементной модели ракеты пакетной схемы при различных значениях переменных ${{x}_{i}}$, $i = \overline {1.{\text{ }}5} $ позволяет сделать заключение: определены собственные пары ${{n}_{p}}$ низших тонов собственных колебаний объекта (задано ${{n}_{p}} = 48$, что соответствует учету собственных частот колебаний до 65 Гц). По завершении вычислений решающий алгоритм автоматически упорядочивает собственные значения по возрастанию (имеет также место соответствующее упорядочение собственных пар). При изменении переменных модели отмечена перенумерация (перестановка) собственных пар. Это указывает на принципиальную возможность, при определенных сочетаниях значений переменных модели, появления в собственных спектрах объекта кратных и почти кратных собственных значений, что является характерным для подобных систем и определяет особенности критериальных функций формулируемых обратных задач.

3. Рассмотрим обратную задачу коррекции параметров объекта по модальным данным, которая в рамках выбранной математической модели описывается операторным уравнением

где $X,Y$ – гильбертовы пространства; $A$ – компактный линейный оператор, действующий из $X$ в $Y$. Правая часть возмущенного операторного уравнения представляет приближенные входные данные ${{y}^{\delta }}$, определенные по результатам измерений. Предположим, что погрешность задания входной информации $\delta $ известна и имеет место ${\text{|| }}{{y}^{\delta }} - y{\text{ ||}} \leqslant \delta $. Требуется определить устойчивые приближенные решения ${{x}^{\delta }}$ по заданной приближенной информации ${{y}^{\delta }}$. Во многих приложениях обратные задачи являются некорректно поставленными. Применим метод регуляризации Тихонова [7]. Минимизируемый функционал определен в виде где ${\text{||}}Ax - {{y}^{\delta }}{\text{||}}_{Y}^{2}$ – функционал невязки (квадрат нормы в пространству Y); $\alpha {\text{||}}x{\text{||}}_{X}^{2}$ – стабилизирующий функционал.

Оценка влияния зашумленности входной информации на устойчивость решения обратных задач требует специальных исследований [12].

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

здесь ${{\lambda }_{i}}(x),\lambda _{i}^{*}$ – собственные значения, относящиеся к исходному (текущему) и заданному спектрам; $x$ – вектор переменных модели; $X$ – допустимая область; $n$ – размерность задачи; $J = \{ 1,{\text{ }}...{\text{ }},n\} $; через ${{{\text{R}}}^{n}}$ обозначено $n$-мерное вещественное линейное пространство. Необходимо найти такой вектор переменных модели, который минимизирует отличия между сравниваемыми спектрами (настройка модели на заданный спектр).

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

где где $f(x)$ – критеральная функция; ${{g}_{i}}(x)$ – функции ограничений задачи, $i \in I$; $I = \{ {\text{ }}1,{\text{ }}...{\text{ }},{{n}_{g}}\} $ – конечное множество индексов; ${{n}_{g}}$ – число функций ограничений; $D$ – область поиска; ${{a}_{j}},{{b}_{j}}$ – соответственно нижняя и верхняя границы переменной ${{x}_{j}}$; ${{x}^{*}}$ – глобальное решение. Функции $f(x)$, ${{g}_{i}}(x)$, $i \in I$, задачи (3)–(5) предполагаем непрерывными липшицевыми. Предположим, что действительная функция $f:{{{\text{R}}}^{n}} \to {\text{R}}$ является многоэкстремальной, не всюду дифференцируемой и для нее задана вычислительная процедура, позволяющая определять значения функции в точках допустимой области.

4. К настоящему времени разработаны и находят многочисленные применения эффективные методы глобальной оптимизации. Обзор современных детерминированных методов глобальной оптимизации дан в работе [13]. Эффективность детерминированных алгоритмов существенно ограничена их зависимостью от размерности задачи. В случае большого числа переменных применяют алгоритмы стохастической глобальной оптимизации. Кратный алгоритм столкновения частиц M-PCA (Multi-Particle Collision Algorithm) основан на алгоритме Метрополиса и входит в число наиболее мощных современных стохастических алгоритмов глобальной оптимизации [14]. Применение стохастических алгоритмов глобальной оптимизации требует значительных вычислительных ресурсов. Одним из путей повышения эффективности таких алгоритмов является совершенствование процедуры локального поиска. Этот подход реализуют при разработке гибридных алгоритмов глобальной оптимизации. Например, гибридный алгоритм HJPCA объединяет стохастический алгоритм столкновения частиц PCA, используемый для сканирования пространства переменных, и детерминированный алгоритм локального поиска методом Хука–Дживса [15]. Описание других гибридных алгоритмов глобальной оптимизации приведено в работе [16].

Структуры гибридных алгоритмов глобальной оптимизации, построены на основе стохастических алгоритмов M-PCA и CBM-PCA (кратного алгоритма столкновения частиц с централизованной выборкой) [17]. Работа алгоритма M-PCA основана на использовании аналогии с процессами абсорбции и рассеяния частиц при ядерных реакциях. В его простейшей версии для исследования области поиска используется одна частица: указанная версия алгоритма M-PCA совпадает с алгоритмом PCA, интегрированным в гибридный алгоритм HJPCA. На начальном шаге выбирается пробное решение (Old_Config), которое затем модифицируется посредством стохастического возмущения (Perturbation), что позволяет найти новое решение (New_Config). С помощью функции Fitness дается сравнительная оценка нового и предыдущего решений, на основании которой новое решение можно принять или отвергнуть. Если новое решение отвергнуто, то происходит переход к функции Scattering, реализующей схему Метрополиса. Для сканирования области, перспективной на минимум, применяются функции Perturbation и Small_Perturbation. Эффективность описанного поиска глобального решения алгоритмом значительно повышается за счет одновременного использования большого числа частиц. Такой подход реализует алгоритм M-PCA, который непосредственно ориентирован на применение в среде параллельных вычислений. Наилучшее решение определяется с учетом данных о всех частицах, участвующих в процессе.

В работе [15] приведены результаты сравнения эффективности гибридного алгоритма HJPCA и современных, обладающих высокими характеристиками алгоритмов, реализующих метаэвристики GLOBAL и C-GRASP [18]. По результатам сравнения алгоритма C-GRASP с генетическим алгоритмом Genocop III и алгоритмами, реализующими метаэвристики Scatter Search и DTS_ASP, установлено, что при решении сорока стандартных эталонных тестовых задач глобальной оптимизации алгоритм C-GRASP оказался более эффективным. По результатам тестирования гибридного алгоритма HJPCA, а также алгоритмов GLOBAL и C-GRASP (данные получены для четырнадцати стандартных эталонных тестовых многомерных функций), установлена более высокая эффективность алгоритма HJPCA. При использовании алгоритма HJPCA среднее количество вычислений значений критериальных функций составило 14.5% относительно алгоритма C-GRASP.

Реализованы новые гибридные алгоритмы, интегрирующие алгоритм CBM-PCA и детерминированные методы локальной минимизации. Первый предложенный гибридный алгоритм CBM-PCAHJ объединяет стохастический алгоритм CBM-PCA (сканирование пространства переменных) и метод Хука–Дживса (локальный поиск). Второй гибридный алгоритм CBM-PCALMSI интегрирует стохастический алгоритм CBM-PCA и детерминированный метод линеаризации для сглаживающих аппроксимаций критериальных функций LMSI [19]. Благодаря более эффективному сканированию пространства переменных новый гибридный алгоритм CBM-PCAHJ является более мощным по сравнению с алгоритмом HJPCA. Аналогичная оценка эффективности справедлива и в отношении второго предложенного здесь алгоритма CBM-PCALMSI.

5. Рассмотрим обратную задачу коррекции параметров ракеты пакетной схемы по измеряемым модальным данным. Предположим, что экспериментальное определение собственных форм не обеспечивает требуемой точности, поэтому регистрируемые модальные данные представлены только ограниченным собственным спектром, соответствующим ${{n}_{p}}$ низшим собственным частотам объекта. Входная информация моделируется решением прямой задачи для объекта при значениях коэффициентов жесткости узлов межблочных связей: $c_{x}^{в}$ = 0; $c_{y}^{{\text{в}}}$ = 0.3 ⋅ 10^–7 Н/м; $c_{z}^{{\text{в}}}$ = 0.4 ⋅ 10^–7 Н/м; $c_{x}^{{\text{н}}}$ = = 0.6 ⋅ 10^–7 Н/м; $c_{y}^{{\text{н}}}$ = 0.35 ⋅ 10^–7 Н/м; $c_{z}^{{\text{н}}}$ = 0.45 ⋅ 10^–7 Н/м. Задано ${{n}_{p}} = 48$, при этом учитываются собственные частоты колебаний в интервале до 65 Гц. Предположим, что относительная погрешность входных данных не превышает 0.01%, шумы отсутствуют. Переменные модели представлены вектором $x \in {{{\text{R}}}^{5}}$, определяющим относительные (в %) значения жесткостных характеристик межблочных связей: $0 \leqslant {{x}_{i}} \leqslant 100\% $, $i = \overline {1,{\text{ }}5} $. Критериальная функция обратной задачи определена в виде

где ${{\gamma }_{j}},{\text{ }}{{f}_{j}}(x)$ – весовой коэффициент и частный критерий, соответствующие $j$-му собственному значению, ${{f}_{j}}(x) = (\lambda _{j}^{*} - {{\lambda }_{j}}(x))$; $\lambda _{j}^{*}$ – составляющая заданного входного собственного спектра; ${{\lambda }_{j}}(x)$ – составляющая собственного спектра, определяемого текущим вектором $x \in {{{\text{R}}}^{5}}$; $j = \overline {1,{\text{ }}{{n}_{p}}} $; $\alpha $ – параметр регуляризации; $\alpha {\text{||}}x{\text{|}}{{{\text{|}}}^{2}}$ – стабилизирующий функционал.

Численное решение задачи минимизации критериальной функции $F(x)$ получено с использованием гибридного алгоритма CBM-PCAHJ. После определения области, перспективной на глобальный минимум, локальный поиск реализуется методом Хука–Дживса. Изменение переменных модели ${{x}_{i}}$, $i = \overline {1,{\text{ }}5} $, при возрастании числа итераций в завершающей фазе локального поиска показано на рис. 2; соответствующее изменение значений критериальной функции представлено на рис. 3.

По завершении семнадцати итераций получено $F({{x}^{{17}}}) \approx 9.9263$; целевые $x_{i}^{*}$ и вычисленные приближенные $x_{i}^{{17}}$, $i = \overline {1,{\text{ }}5} $, значения переменных модели приведены в табл. 2.

Значения жесткостных характеристик межблочных связей, восстановленные по приближенным значениям переменных модели, соответственно равны: $c_{y}^{{\text{в}}} \approx 0.300824 \cdot {{10}^{{ - 7}}}$ Н/м; $c_{z}^{{\text{в}}} \approx 0.39774 \cdot {{10}^{{ - 7}}}$ Н/м; $c_{x}^{{\text{н}}} \approx 0.{\text{60019}} \cdot {{10}^{{ - 7}}}$ Н/м; $c_{y}^{{\text{н}}} \approx 0.34970 \cdot {{10}^{{ - 7}}}$ Н/м; $c_{z}^{{\text{н}}} \approx 0.45175 \cdot {{10}^{{ - 7}}}$ Н/м. Наибольшая относительная погрешность решения (определяется восстановленным значением параметра $c_{z}^{в}$) здесь не превышает 0.57%. Можно отметить достаточную согласованность точности полученного приближенного решения модельной задачи и заданной входной информации.

Заключение. Представлена методика коррекции параметров ракеты пакетной схемы по модальным данным. Для математической модели объекта решение прямой задачи получено методом итерации подпространства. Сформулирована обратная задача коррекции параметров объекта, при решении которой используется оптимизационный подход. Предложены два новых гибридных алгоритма глобальной оптимизации. По результатам косвенных оценок установлена достаточно высокая эффективность предложенных алгоритмов. Решение модельной задачи коррекции жесткостных характеристик узлов межблочных связей по модальным данным получено с использованием гибридного алгоритма CBM-PCAHJ. Точность полученного приближенного решения согласована с точностью задания входной информации.

Конфликт интересов. Авторы заявляют об отсутствии конфликта интересов.