Продолжаем знакомиться с таким методом ускорения проектирования как суррогатное моделирование. Это перевод отличной статьи авторства Шуая Гуо, который мы дополнили своим кодом и некоторыми размышлениями. Если у вас получится повторить пример, или у есть задача, которую можно было бы решить таким образом – зовем высказаться в комментариях, нам будет что обсудить!
В первой части этой серии статей мы обсудили идею использования суррогатных моделей для ускорения такого этапа проектирования изделий как имитационное моделирование. Основной метод здесь – обучение статистической модели, которая будет служить дешевым, но точным заменителем имитационных моделей при выполнении различных задач проектирования, что значительно повышает эффективность анализа.
В части II мы рассмотрим конкретный пример, чтобы продемонстрировать, как использовать суррогатные модели на практике. Дорожная карта для этого примера показана ниже:
Мы начнем с изложения изучаемой проблемы, а затем применим к ней метод суррогатного моделирования. В конце мы проиллюстрируем, как использовать обученную суррогатную модель для выполнения двух типов анализа.
1. Постановка задачи
Турбина – ключевой компонент реактивного двигателя. Сгоревшие газы из камеры сгорания попадают на лопатки турбины и заставляют ее вращаться для создания тяги.
Лопатки турбины работают в чрезвычайно жестких условиях, поскольку они непосредственно сталкиваются с высокотемпературными выхлопными продуктами горения. Экстремальные температуры могут расплавить материал лопаток и поставить под угрозу надежную работу реактивного двигателя.
Как следствие, тепловой анализ является необходимостью для проектирования лопаток турбины. Для этого обычно строится физическая модель для имитации распространения энергии и оценки температурного поля внутри лопатки, учитывая в качестве граничных условий температуры сгоревшего газа и охлаждающего воздуха внутри лопатки.
В данном примере мы используем следующую схему анализа: в качестве входных переменных мы рассматриваем температуры сгоревшего газа и охлаждающего воздуха и соответствующие им коэффициенты теплопроводности. Нас интересует прогнозирование максимальной температуры лопатки (выход) с учетом входных значений. На среднем персональном компьютере один прогон моделирования занимает 1-2 секунды, но помним что это учебный проект.
2. Суррогатное моделирование
Для проектирования лопаток может потребоваться много прогонов модели, чтобы понять, как лопатка реагирует на различные значения температуры газа, температуры охлаждающего воздуха и соответствующие коэффициенты теплопроводности. Одним из способов ускорить этот процесс является обучение суррогатной модели перед началом «медленного» моделирования. Именно это мы и собираемся сделать в данном примере.
Исходные данные для этого примера доступен на сайте MathWorks и входит в стандартную поставку MATLAB. На основе этого примера мы создадим нашу "дорогостоящую и сложную" модель для получения максимальной температуры лопатки при следующих заданных условиях:
заданная температура набегающего потока и охлаждающего потока внутри лопатки
заданная теплопроводность интерфейса воздух-лопатка по "горячей стороне" и по "холодной стороне" лопатки
function max_temp = getMaxTemp(args, thermalmodel)
h_air = args(1); T_air = args(2); h_gas = args(3); T_gas = args(4);
% Создаем модель температурного состояния
thermalmodel = createpde('thermal','steadystate');
% Импортируем геометрию лопатки и создаем сетку
importGeometry(thermalmodel,'Blade.stl');
msh = generateMesh(thermalmodel,'Hmax',0.01);
thermalmodel.Mesh = msh;
% Материал лопатки (в примере это сплав никеля NIMONIC 90)
kapp = 11.5; % зададим теплопроводность в Вт/м/К
thermalProperties(thermalmodel,'ThermalConductivity',kapp);
% Задать граничные условия – теплопроводность на границе воздуха и материала лопатки
thermalBC(thermalmodel,'Face', [15 12 14], 'ConvectionCoefficient',h_air, 'AmbientTemperature',T_air); % Interior cooling
thermalBC(thermalmodel,'Face', [11 10 13 1], 'ConvectionCoefficient',h_gas, 'AmbientTemperature',T_gas); % Pressure side
thermalBC(thermalmodel,'Face',[6 9 8 2 7], 'ConvectionCoefficient',15, 'AmbientTemperature',400); % Root in contact with hot gases
% База лопатки отдаёт тепло другим металлическим частям ротора
thermalBC(thermalmodel,'Face',[3 4 5], 'ConvectionCoefficient',1000, 'AmbientTemperature',300);
% Запустить модель и получить решение (максимальну температуру)
Rt = solve(thermalmodel);
max_temp = max(Rt.Temperature);
endВсе параметры заданы для соответствующих граней модели лопатки (для дальнейших уточнений предлагаем ознакомиться с исходным примером).
2.1 Выборка данных
Для начала мы генерируем несколько обучающих примеров для обучения суррогатной модели.
rng default; % Важно закрепить seed для воспроизводимости результатов
ls_max_points = 100;
n_initial_lhs_points = 8;
var_labels = { 'h_{в.охл}', 'T_{в.охл}', 'h_{г}', 'T_{г}' };
var_limits = [[20, 40]; [120, 180]; [30, 70]; [800, 1200]];
% Распределение по латинскому квадрату
LS = lhsdesign( ls_max_points, 4 );
X_lhs = rescale( LS, ...
repmat( var_limits(:,1)', [ls_max_points,1] ), ...
repmat( var_limits(:,2)', [ls_max_points,1] ));На этом этапе полезно, чтобы точки в распределении были равномерно распределены по всему пространству входных параметров. Для этого мы здесь будем использовать схему выборки латинского гиперкуба. Для примера мы сгенерируем всего 8 точек для начальной выборки и перейдем к следующему шагу.
Теперь мы можем перейти к следующему шагу.
2.2 Обучение модели
На основе первоначально сгенерированных образцов мы можем обучить суррогатную модель. В данном случае в качестве суррогатной модели мы выбрали гауссовский процесс (GP). Используем функцию getMaxTemp для генерации вектора таргетов y.
y = zeros([size(X, 1), 1]);
for i = 1:size(X, 1)
y(i) = getMaxTemp(X(i,:));
endОсновная мотивация в выборе GP-модели заключается в том, что, в отличие от многих других методов машинного обучения с учителем, GP выдает предсказания на неразмеченных образцах, оценивая при этом статистическую неопределенность предсказания в новых точках. Эта особенность является ключом к активному обучению, к которому мы как раз переходим.
2.3 Активное обучение
Итак, мы хотим улучшить первоначальную суррогатную модель с помощью дополнительных обучающих примеров. Мы делаем это итеративно, на каждой итерации генерируем только одну пару «X ; y» (наши признаки – это 4 упомянутых выше проектных параметра, наша целевая переменная: максимальная температура) и добавляем ее к существующему набору обучающих данных.
Чтобы определить, какие точки нужно добавить к датасету, мы ищем в пространстве входных параметров такие области, где текущая суррогатная модель имеет наибольшую оценку ошибки предсказания (confidence interval, CI). Это то место, где текущая суррогатная модель «чувствует себя наименее уверено» относительно своего прогноза.
Есть много способов поиска таких точек, их выбор определяется поведением модели, характером датасета и т.д. Для простоты покажем, как это делать при помощи поиска по сетке (grid search, для комбинаторики был задействован пакет Combinator):
function [max_args, max_yci, max_ypred] = ...
findMaxCIPoint_grid(mdl, limits, grid_steps)
% Все пермутации для 4-х векторов параметров модели (индекс)
idx_grid = combinator(grid_steps, size(limits,1), 'p','r');
% Масштабируем сетку от 0..1 до размеров пространства входных параметров
N = size(idx_grid,1);
X_grid = rescale( idx_grid, ...
repmat(limits(:,1)', [N,1]), ...
repmat(limits(:,2)', [N,1]) );
% Получаем предсказание для всех комбинаций по сетке
[ypred_grid, yci_grid] = predict(mdl, X_grid);
% Находим точку с максимальной неопределенностью
[max_yci,idx] = max(yci_grid);
max_args = X_grid(idx, :);
max_ypred = ypred_grid(idx, :);
endДисклеймер: мы исследовали комбинацию методов поиска точки максимальной неопределенности, добавив рандомизации и несколько шагов градиентого спуска. Это большая тема, о которой стоит поговорить в отдельной статье. Интересный подход использован исходным автором в его статье про суррогатную оптимизацию, но оптимизация и моделирование – слегка разные задачи.
Выбирая новые точки из пространстве параметров в тех областях, где предсказание модели делает наибольшую ошибку, суррогатная модель сможет «обучаться» быстрее всего.
% Для воспроизводимости кода:
n_grid_steps = 31; % шаг разбиения пространства для grid search
n_max_added_points = 20; % Сколько точек мы добавим к датасету
% Установим параметры GP-модели
opts = { 'BasisFunction',basisFunction, ...
'KernelFunction',kernelFunction };
% Датасет для накопления точек будет храниться в переменных X_new и y_new
X_new = zeros(size(X,1)+n_max_added_points, size(X,2));
y_new = zeros(size(y,1)+n_max_added_points, 1);
X_new(1:size(X,1), 1:size(X,2)) = X;
y_new(1:size(y,1),:) = y;
% Обучаем модель на исходных данных
gprMdl = fitrgp(X, y, opts{:});
% Сохраним предсказания и оценку ошибок модели перед добавлением новых данных
[initial_ypred, initial_confidence_intervals] = resubPredict(gprMdl);
for i = 1:n_max_added_points
% Поиск точки с наибольшей потенциалной информацией
[best_args, max_yci, max_ypred] = ...
findMaxCIPoint_grid(gprMdl, var_limits, n_grid_steps);
% Запускаем "тяжелую" модель и дополняем датасет новыми параметрами
X_new(size(X,1)+i, :) = best_args;
y_new(size(y,1)+i, :) = getMaxTemp(best_args, tmodel);
% Заново обучаем модель на дополненном датасете
gprMdl = fitrgp(X_new(1:size(X,1)+i,:), y_new(1:size(y,1)+i,:), opts{:} );
% Сохраняем прогнозы и ошибки новой обученной модели
[ypred, yci] = resubPredict(gprMdl);
endИногда, после добавления всего десятка дополнительных обучающих примеров, мы добивались того, что окончательная ошибка предсказания падала ниже 3% от начального значения. В целом, активное обучение для такого уровня улучшения обычно приводило к добавлению в датасет 10-20 новых точек.
Успех зависит от исходного распределения точек, который зависит от настройки генератора случайных чисел (seed). Это известная проблема, хорошо знакомая все в data science. :)
Для дискуссии, вот примеры графиков поиска точек с наибольшей информативностью в пространстве:
2.4 Испытания
Чтобы оценить точность модели, создаем тестовый набор данных, содержащий 20 тестовых образцов. Построим график чтобы сопоставить прогнозы разных моделей, тяжелой и легкой:
figure
[ypred,yci] = resubPredict(gprMdl);
plot(y_new, ypred, 'bo');
Видно, что данные, полученные из суррогатной модели, прекрасно сопоставляются с данными по исходной высокоточной модели, поскольку все предсказания полностью совпадают с результатами моделирования (грубых ошибок нет, а точность соответствия можно повышать до нужных значений при помощи активного обучения).
3. Развертывание модели
Теперь, когда мы обучили нашу суррогатную модель и оценили ее точность предсказания, пришло время развернуть ее для выполнения нужных нам задач анализа.
3.1 Визуализация
Суррогатная модель может легко быть использована для визуализации «ландшафта» функции «вход-выход». Это помогает аналитикам быстро выявить важные зависимости между выходом и различными входами.
Например, на этом графике показана зависимость между температурой газа и максимальной температурой лопатки при различных значениях коэффициента теплопроводности. Видно, что температура лопатки является почти линейной функцией от температуры газа: при повышении температуры газа температура лопатки пропорционально увеличивается, что имеет интуитивный смысл.
Также мы видим, что связь между температурой лопатки и коэффициентом теплопроводности газа немного нелинейна в режиме высоких температур. Шаг увеличения максимальной температуры лопатки начинает снижаться относительно шага увеличения температуры продуктов горения, если при этом сохраняется коэффициент теплопроводности. То же самое происходит с каждым следующим шагом повышения коэффициента теплопроводности.
3.2 Количественная оценка неопределенности
На практике условия эксплуатации турбины постоянно меняются. То есть мы не уверены в точных значениях температур газа и охлаждающего воздуха и связанных с ними коэффициентов теплопроводности. В результате нам необходимо оценивать характер изменения характеристик нашего продукта, вызванное неопределенными условиями эксплуатации. Именно этого и помогает добиться анализ количественной оценки неопределенности.
Обычно для количественного анализа неопределенности мы используем метод Монте-Карло.
В двух словах, метод Монте-Карло начинается с назначения распределений вероятностей неопределенным входным переменным (распределения выражают наши текущие знания о том, насколько эти переменные «не определены»). Затем мы берем большое количество случайных точек (~o(10⁴)) из входного распределения и вычисляем соответствующие им выходные значения. На основе ансамбля выходных результатов мы можем построить гистограмму выхода, которая статистически полно описывает неопределенность выходных переменных.
Для проведения Монте-Карло анализа последуем вышеупомянутым процедурам. Предположим, что наши четыре входа независимы и имеют нормальное распределение (что на графике). Мы берем в общей сложности 10000 случайных точек, чтобы выполнить количественную оценку неопределенности.
Распределение максимальной температуры лопатки, предсказанное суррогатной моделью, показано на графике. Также мы сравниваем этот результат с распределением, полученным путем применения процедуры Монте-Карло непосредственно к исходной высокоточной модели. Видим, что суррогатная модель точно оценила изменение максимальной температуры лопатки, связанное с неопределенностью входных параметров.
С точки зрения вычислительных затрат, Монте-Карло на основе суррогатной модели занимает 1/100 с на довольно стандартном компьютере, в то время как анализ Монте-Карло на основе симуляции занял около 4 часов. Такого ускорения можно добиться с помощью хорошо обученной суррогатной модели.
4. Выводы
В этой статье мы рассмотрели, как использовать метод суррогатного моделирования для ускорения теплового анализа лопатки турбины реактивного двигателя. Мы продемонстрировали ключевые этапы процесса суррогатного моделирования:
Начальная выборка: метод латинского гиперкуба для создания датасета, хорошо «заполняющего пространство»
Обучение модели: Гауссовский процесс
Активное обучение: последовательное обогащение обучающей выборки для максимизации эффективности обучения модели
Тестирование модели: сравнительный анализ с использованием тестового набора данных
Развертывание модели: визуализация зависимостей и анализ количественной оценки неопределенности
