Spiralhead 20 авг в 10:01

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

14 мин

3.1K

Блог компании ГК ЛАНИТCAD/CAM*Инженерные системы*

Туториал

Даже небольшое преимущество продукта может принести существенную выгоду. Инженеры постоянно ищут новые способы оптимизации конструкций в существующих ограничениях, чтобы добиться лучших результатов. Так, Airbus в 2006 году анонсировал программу, которая позволила добиться прироста на ~4 процента в показателях эффективности для самолета A320. Отчасти благодаря этому с 2009 по 2016 года (до появления A320neo с новыми двигателями) компания нарастила продажи A320 на ~40% по сравнению с основным конкурентом. В мире будут доминировать те, кто смогут проводить оптимизацию быстрее и эффективнее. Так можно ли ускорить сам процесс оптимизации? В этой статье в блоге ЛАНИТ я бы хотел поделиться одним подходом к оптимизации конструкции, который позволит это сделать.

Для оценки показателей эффективности на этапе проектирования часто применяются CAE- и CFD-программы. Причем обычно они используются в режиме так называемого поверочного расчета: у нас есть изделие с геометрией X, зависящее от исходных требований D, надо выяснить, как оно себя поведет в реальном мире при граничных условиях (ГУ) и начальных условиях (НУ) Y. Для этого мы используем современный CAE-продукт: делаем сетку, задаём ГУ, НУ и параметры решателя. Затем запускаем расчет, достигаем сходимости, получаем рассчитанное распределение полей физических величин Q и вычисляем целевую функцию f. Профит! Расчет CAE используют в настоящем варианте как инструмент подтверждения рабочих гипотез.

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

Берём какой-нибудь хороший оптимизатор, прикручиваем к нашему расчету и запускаем! Среди коммерческих продуктов можно отметить HEEDS и отечественный PSeven. Если вы скряга, как я, и умеете что-то в программировании, то можно собраться на свободно распространяемых оптимизаторах (Scipy, GNU Octave, Gradient-Free-Optimizers).

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

Расчеты могут выполняться от нескольких минут до нескольких месяцев. Навскидку (авторские наблюдения) нормальный средний cfd-расчет в промышленности занимает от часа до шести. Пусть час. Если у вас оптимизируемых параметров меньше трех, то такая оптимизация более-менее эффективна. А если больше 10?

Приведу известный пример от scipy.

def rosen(x):
    """The Rosenbrock function"""
    return sum(100.0 * (x[1:] - x[:-1] ** 2.0) ** 2.0 + (1 - x[:-1]) ** 2.0)
x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])
res = optimize.minimize(
    rosen, x0, method="nelder-mead", options={"xatol": 1e-8, "disp": True}
)

Задача оптимизации в безградиентной постановке для пяти исходных параметров требует >500 вызовов функции.

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

Грустно? А то!

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

«Но позвольте!», ― скажете вы. Мы же умеем оптимизировать (обучать) супербольшие нейронные сети. Там миллиарды параметров, и все успешно обучается. Неужели нельзя как-то использовать те же подходы только для оптимизации конструкции?

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

Как мы обучаем нейросети

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

При обучении нейросети используется секретный ингредиент ― градиент.

Градиент ― матрица чувствительности выходных параметров от входных, которую с помощью обратного распространения ошибки (или цепного правила дифференцирования), мы можем рассчитать. Вектор показывает направление локального экстремума. Это достаточно сильный сигнал, который позволяет оптимизировать системы с таким гигантским количеством параметров, как у современных языковых моделей (1,76 трлн параметров у GPT-4).

А нам такой подход подойдет? Что у нас, у инженеров?

Процесс можно условно представить так:

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

Давайте попробуем разобраться. Но прежде чем окунуться с головой в математическое мамбо-джамбо, предлагаю сначала придумать задачку, на которой будем экспериментировать.

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

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

Соответственно, нашу геометрию задают 12 параметров ― значения 12 координат оси Z для точек сплайна. Остальные считаются константами.

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

$f=\left\{\begin{matrix}C_y-5\cdot C_y(C_x-0.1) \; \text{if}\; C_x>0.1\\C_y\ \text{otherwise}\end{matrix}\right.$

Она следует методу множителя Лагранжа для задачи с ограничением. Коэффициент $\lambda=5$ выбран эмпирически. Такой подход не позволит точно обеспечить условие $C_x \leq 0.1$ , но точное решение и не нужно. Хотелось просто получить аэродинамический профиль из пластины так, чтобы , не росли до бесконечности.

Этапы решения

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

Все выглядит не очень сложно, но как найти градиенты у каждого шага?

Как найти градиент на геометрию

Что мы знаем о процессе построения трехмерной геометрии на современных CAD-пакетах? Конструкторы используют инструментарий различных операций, чтобы получать поверхности. Если они герметично ограничивают некий объем, то считается, что это твердое тело (solid). Сами поверхности ― это, по сути, формулы, которые преобразуют формирующие их параметры в координаты облака точек, лежащих на этих поверхностях.

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

В нашей задаче мы будем вычислять градиент для срединной поверхности крыла от координаты Z для 12 точек. Это так называемая Non-Uniform Rational B-Spline (NURBS) поверхность или поверхность, образованная протяжкой NURBS-кривых.

Для расчета градиента мы воспользуемся библиотекой расчета градиента NURBS-кривых и поверхностей, которая написана для всем известного фреймворка обучения нейросетей PyTorch. Она позволяет рассчитывать координаты NURBS-кривых и поверхностей по координатам ключевых точек, а также использовать градиенты этих функций в оптимизаторе Torch. Фактически операция построения NURBS ведет себя как слой нейронной сети при таком подходе.

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

from NURBSDiff.surf_eval import SurfEval
import torch
def fun(self, input_parameters, num_ctrl_pts1=6, num_ctrl_pts2=2):
    # Превращение входных параметров в дифференцируемый тензор
    inp_ctrl_ys = torch.tensor(input_parameters).flatten()
    inp_ctrl_ys.requires_grad = True
    # Создание слоя по расчету поверхности
    surf_eval = SurfEval(
        # Количество точек по направлению u
        num_ctrl_pts1,
        # Количество точек по направлению v
        num_ctrl_pts2,
    )
    # Создание тензора точек полюсов срединной поверхности
    # и их весов
    x_ctrl = torch.linspace(0, 0.5, num_ctrl_pts1)
    z_ctrl = torch.linspace(0, 1, num_ctrl_pts2)
    X_ctrl, Z_ctrl = torch.meshgrid(x_ctrl, z_ctrl)
    X_ctrl = X_ctrl[:, :, None]
    Z_ctrl = Z_ctrl[:, :, None]
    inp_ctrl_pts = torch.cat((X_ctrl, inp_ctrl_ys, Z_ctrl), 2).reshape(
        (1, num_ctrl_pts1, num_ctrl_pts2, 3)
    )
    weights = torch.ones(1, num_ctrl_pts1, num_ctrl_pts2, 1)
    # Возвращение тензора точек срединной поверхностей с графом вычислений
    return surf_eval(torch.cat((inp_ctrl_pts, weights), -1))

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

Функция дифференцируема. Все тип-топ. Но для нашего расчета еще понадобится настоящая 3D-геометрия. Да еще и строить ее надо автоматически. Где ее взять?

Я не стал искать еще одну удобную библиотеку на Python, потому что под рукой был Siemens NX. С помощью него можно построить геометрию, которая зависит от параметров, и перестраивать ее автоматически, используя язык макросов NX Journaling.

В таком случае схема работы с геометрией выглядит так:

*Схема подготовки геометрической модели и дифференцируемого уравнения срединной поверхности*

Большинство операций над геометрией и так параметрические, так что перестраивалась она лихо. Была только загвоздка с сохранением геометрии в формат x_t. NX Journaling позволяет записывать макросы по мере работы пользователя над геометрией, но операцию сохранения в Parasolid записать нельзя. Пришлось искать код на сайте https://www.nxjournaling.com/ и адаптировать под свои нужды. Сайт очень советую. Он несколько раз меня выручал.

Код функции перестроения приведен здесь.

Как найти градиент у 3D-симуляции

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

$\frac{ \partial}{ \partial t} \int\limits_V WdV + \oint [F-G] da = \int\limits_V HdV$

Векторы определены как

$W = \begin{bmatrix} \rho \\ \rho \upsilon \\ \rho E\end{bmatrix}$ , $F = \begin{bmatrix} \rho \upsilon \\ \rho \upsilon \upsilon + \rho I \\\rho \upsilon H + \rho \upsilon\end{bmatrix}$ ,

$G = \begin{bmatrix} 0 \\ T \\ T\cdot \upsilon + \dot{q} \end{bmatrix}$ , $H = \begin{bmatrix} S_u \\ f_r + f_g + f_p + f_u + f_w + f_L \\S_u \end{bmatrix}$ ,

где:

$\rho$ - плотность жидкости,
$\upsilon$ - скорость жидкости,
- давление,
- удельная полная энергия,
- тензор вязких напряжений,
- тепловой поток,
- вектор объемных сил.

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

Adjoint Solver решает задачу в виде: $\frac {dL} {dD} = [\frac{\partial L} {\partial X} +\frac{\partial L} {\partial Q}\frac{\partial Q} {\partial X}] \frac{\partial X} {\partial D}$ , где:

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

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

Здесь я не готов глубоко погружаться в описание присоединенного решателя. Это тема для отдельной статьи. Если интересно, напишите в комментариях.

Среди известных мне решателей Adjoint Solver поддерживают Fluent и Star-CCM+. Среди отечественных я таких кодов не знаю. Если что-то все же есть, то напишите, пожалуйста, в комментариях. Я для решения задачи буду использовать Star-CCM+.

Еще стоит прояснить моменты, связанные с пробрасыванием градиента до точек, созданных NURBSDiff. На конечном этапе Adjoint Solver происходит расчет матрицы чувствительности целевой функции к искомым параметрам в виде:

$\frac{dL^T}{dD} = \frac{dX^T}{dD} \frac{dL^T}{X}$ ,

где:

$\frac{dL^T}{dD}$ - матрица чувствительности целевой функции к сеточной модели,
$\frac{dX^T}{dD}$ - матрица чувствительности сетки к оптимизируемым параметрам.

$\frac{dX^T}{dD}$ рассчитывается из учета того, что изменение сетки происходит (в Star-CCM+) через радиальные базисные функции (Radial Basis Functions):
$x(x^0) = \alpha + \displaystyle\sum_{j=1}^{N} \beta_j\phi_j (r_j(x^0))$ ,

где:
$r_j (x^0) = \| x^0 - x_j^0\|$
$\phi_j (r_j) =\sqrt{r_j^2 + c_j^2}$

- начальная координата узла сетки,
- новая координата узла сетки.

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

Схема получения градиента для CFD-расчета выглядит так:

*Схема вычисления градиента на этапе CFD-решения*

В Star-CCM+ мне достаточно настроить только целевую функцию и пробросить решение до точек срединной поверхности, которые мне генерирует моя функция (1). Конечно, общие вещи по построению сетки и заданию граничных и начальных условий остаются за пользователем.

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

*Задание Report (а), его определение (б) и внесение его как целевой функции (в)*

Для того, чтобы найти градиент целевой функции до наших точек, их достаточно внести в систему в виде Point Set (набор точек). Задается он через таблицу, которую, в свою очередь, можно импортировать из внешнего файла.

AllSensitivity в папке с таблицами появилось не напрасно. Это уже способ экспортировать собранные градиенты и сохранить их в файл. Star-CCM+ позволяет превратить любые поля в таблицы, которые затем можно сохранить в формате csv.

Конечно, все действия надо делать автоматически, и для этого в Star-CCM+ есть макроязык, основанный на Java. Для проведения этапа трехмерного расчета в глобальной схеме графа вычислений необходимо:

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

Для этого написан java-скрипт, который приведен здесь.

Собираю все вместе

Оркестрировать весь процесс буду на Python. Для общения между программами буду использовать самый примитивный интерфейс ― файлы. Общая схема процесса расчета показана на рисунке:

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

Проверяю идею на оптимизаторах первого порядка

Прежде чем переходить к результатам с учётом градиента, я провел пуск с привлечением коммерческого супер-пупер-оптимизатора SHERPA от Siemens HEEDS, который, тем не менее, не использует градиент целевой функции. Он задействует гибридный метод оптимизации (ссылка на статью про SHERPA), при котором одновременно запускаются различные оптимизационные алгоритмы, а потом из них выбирается лучший.

Так вот результаты:

Да, за ~100 итераций не было особого движения в сторону приемлемого аэродинамического профиля. Так как расчет занимал довольно много времени, я решил дальше не считать. Давайте перейдем к результатам с градиентом. Здесь в качестве оптимизатора я сначала использую обычный градиент с шагом обучения .

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

Да, результаты не идеальны. Присутствует постоянное дребезжание профиля около сошедшегося решения. Можно ли улучшить решение, если использовать более продвинутый оптимизатор? Например, известный ADAM.

ADAM справился хуже ожидаемого. Из-за накопления движения он проскакивает минимум и довольно долго прыгает около него. Скорее всего, можно настроить его гиперпараметры получше, но это потребует продолжительных запусков, а этого хотелось бы избежать.

Можно ли как-то ускорить сходимость? Можно, и для этого как раз придумали оптимизаторы второго порядка.

Как ускорить расчет: оптимизатор второго порядка

Оптимизаторы второго порядка редко используются в обучении нейронных сетей из-за больших требований к памяти, а вот для задач до <1000 параметров они подходят замечательно. Я экспериментировал с несколькими такими оптимизаторами, например, с LBFGS. Однако лучше всего сработал SLSQP (Sequential Least Squares Programming ― последовательное программирование по методу наименьших квадратов) из обоймы оптимизаторов Scipy.

Как видно из решения, хороший результат появляется уже на 500 итерации а лучший ― на 2651 итерации (это также лучший результат для всех оптимизаторов). Не исключаю, что SGD со временем мог бы выдать результат и качественней, но по критериям скорость-эффективность побеждает однозначно SLSQP. Кроме того, он позволяет учитывать ограничивающие функции, то есть я мог бы ввести в параметрах ограничение и использовать в качестве целевой функции только коэффициент подъемной силы.

Также весьма любопытны резкие скачки параметров на итерациях 546, 1262 и т.д. Я полагаю, что это работа алгоритма поиска по линии (Line Search), который по градиентам сначала определяет направление поиска, а затем по этому же направлению ищет локальный максимум.

Сравнение результатов

Для сравнения результатов воспользуемся таблицей, в которой собраны лучшие значения целевой функции, номер итерации CFD решателя и итераций оптимизатора для каждого метода:

Метод оптимизации	Кол-во итераций оптимизатора	Кол-во итераций CFD	Лучшее значение			$K=\frac{C_y}{C_x}$
SHERPA	100	30255	0.849	0.919	0.210	4.379
ADAM	5	744	0.927	1.008	0.115	8.765
SGD	27	2651	0.944	1.042	0.119	8.768
SLSQP	12	2306	0.947	1.031	0.116	8.861

Лучшие результаты показали алгоритмы SLSQP и SGD. При этом первому понадобилось еще и меньшее значение итераций оптимизатора. Но еще более очевидна разница между градиентными и безградиентными методами. Скорость сходимости у градиентных методов на порядок превышает таковую у безградиентных методов. Предполагаю, в случае моей задачи нужно >500 итераций, чтобы оптимизатор начал выдавать удовлетворительный результат.

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

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

Забавный факт. Аэродинамическое качество 8,86 превышает, например, показатели F-4E Фантом II. Конечно, сравнивать напрямую показатели самолета и изогнутой доски неправильно. Я хочу скорее указать, что из наложенных ограничений (пластина с равномерной толщиной, перемещение узлов разрешено только по вертикали, большое количество оптимизируемых параметров) удалось получить результат, приближенный к максимально возможному.

Заключение

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

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

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

Что касается проблем градиентных методов, можно разработать гибридные оптимизаторы с использованием подходов безградиентных оптимизаторов, которые принимают сигнал градиента для ускорения сходимости. Пока что гибридные оптимизаторы не учитывают градиент, а работают как безградиентные. Это возможно. В качестве примера приведу использование Population-based training (обучение, основанное на популяции ), с помощью которого можно оптимизировать гиперпараметры и обучать/оптимизировать модель одновременно.

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

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

Я бы хотел, чтобы представленные методы нашли более широкое распространение, чем сейчас. Уверен, от этого бы выиграли все.

*Статья написана в рамках Хабрачелленджа 2.0, который прошел в ЛАНИТ весной 2024 года.

Теги:

Хабы: