Ctrl+K
  • REPEAT
en
  • REPEAT
en

Решение двумерного нестационарного уравнения теплопроводности с помощью PINN#

Внешние аппроксимационные модели машинного обучения позволяют учитывать 2D и 3D зависимости и уравнения при моделировании изделий и систем на платформе REPEAT VISION.

В данном практическом руководстве приведено подробное математическое описание моделей физически-информированных нейронных сетей (PINN), а также приведены идеи из актуальных научных статей по достижению лучшей сходимости решений уравнений.

В процессе выполнения данного практического руководства вы получите следующие навыки:

  • Навыки создания и обучения ML модели для решения задачи теплопереноса в двумерном объекте;

  • Навыки выполнения экспорта модели машинного обучения в формат ONNX для последующей интеграции на платформу REPEAT VISION;

  • Навыки интеграции ML модели на платформу REPEAT VISION для использования модели объекта в проекте системы или изделия.

1.  Описание задачи#

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

Параметры системы

  • Область моделирования задана как \( x \in [0, L_x] \) и \( y \in [0, L_y] \).

  • В качестве материала используется медь с температурапроводностью \( \alpha = 1.12 \cdot 10^{-4} \) м\( ^2 /\)с.

Временной интервал моделирования

  • Моделирование происходит во временном интервале \( t \in [0, 12.0] \).

Начальные и граничные условия

  • Использованы граничные условия вида Дирихле.

  • Температурное распределение задаётся аналитически как \( u(x,y,0) = \sin(\pi \frac{x}{L_x}) \sin(\pi \frac{y}{L_y}) \).

Поскольку максимальная температура в такой постановке достигается в центре области и равна 1 °C (в то время как на границе 0 °C) имеет смысл рассматривать малую область. Например, пусть \( L_x = L_y = 0.1\) м.

Постановку задачи в начальный момент времени можно визуализировать следующим образом:

(png, hires.png)

../../../_images/problem_plotter_plot_problem_setup.png

Распределение температуры в начальный момент времени.#

2.  Математическая постановка задачи#

Модель теплопроводности в двумерном пространстве описывается уравнением (1).

(1)#\[\frac{\partial u}{\partial t} = \alpha \Delta u = \alpha \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right),\]

где:

  • \( u(x,y,t) \) – температура в точке с координатами \( (x,y) \) в момент времени \( t \);

  • \( \frac{\partial u}{\partial t} \) – частная производная температуры по времени, отражающая скорость изменения температуры;

  • \( \frac{\partial^2 u}{\partial x^2} \) и \( \frac{\partial^2 u}{\partial y^2} \) – вторые частные производные по координатам \( x \) и \( y \) соответственно, которые описывают диффузионное распределение температуры по пространству;

  • \( \alpha \) – коэффициент теплопроводности, характеризующий, насколько быстро тепло распространяется в материале.

Часто при записи уравнений в частных производных, для краткости, опускают символы \( \partial \). Тогда переменная, по которой берется производная, становится нижним индексом. Например, уравнение теплопроводности в такой нотации будет иметь вид:

(2)#\[u_t = \alpha \left(u_{xx} + u_{yy} \right),\]

Для решения подобного рода уравнений требуются дополнительные условия, позволяющие найти частное решение уравнения. В начальный момент времени \( t = 0 \) температура задана аналитически функцией:

(3)#\[u(x,y,0) = \sin(\pi x) \sin(\pi y),\]

На границе области \( \partial \Omega \), где \( \Omega = [0,1] \times [0,1] \), температура фиксирована и равна нулю:

(4)#\[u(x,y,t) = 0 \quad \text{при } (x,y) \in \partial \Omega, \quad t > 0.\]

Для данной постановки задачи аналитическое решение можно получить методом разделения переменных (методом Фурье). Предположим, что решение имеет вид:

\[ u(x,y,t) = X(x) \, Y(y) \, T(t). \]

Подставляя это представление в исходное уравнение, после разделения переменных и применения граничных условий (где \( X(0)=X(1)=0 \) и \( Y(0)=Y(1)=0 \)) получаем собственные функции для пространственной части:

\[ X(x) = \sin(\pi x), \quad Y(y) = \sin(\pi y). \]

Временная функция \( T(t) \) определяется из обыкновенного дифференциального уравнения:

\[ \frac{dT}{dt} = -2\pi^2\alpha \, T(t), \]

откуда

\[ T(t) = e^{-2\pi^2\alpha\, t}. \]

Таким образом, полное аналитическое решение записывается следующим образом:

(5)#\[u(x,y,t) = e^{-2\pi^2\alpha\, t} \sin(\pi x) \sin(\pi y).\]

Это решение удовлетворяет начальному условию и нулевым граничным условиям. Эволюция решения во времени имеет следующий вид:

Evolution

Эволюция решения.#

3.  Описание моделей PINN#

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

В случае рассматриваемой задачи дано дифференциальное уравнение в области \( \Omega \) с граничными и начальными условиями. Аналитическое решение \( u(x,y,t) \) аппроксимируется нейронной сетью \( \hat{u}(x,y,t; \theta) \), где \( \theta \) – набор параметров сети. При этом модель принимает на вход пространственные и временную координаты \( (x, y, t) \), а возвращает приближенное значение температуры \( u \).

В отличии от большинства других задач где используются нейронные сети, модель PINN не имеет обобщающих способностей в смысле решения разнообразных задач. Она лишь описывает решение конкретной задачи в заданной области в любой точке пространства и времени. Таким образом, для нахождения величины в момент времени \( t_i \) не нужно находить значения в моменты \( t_{i-1}, t_{i-2}, ... \), как это нужно было бы, к примеру, при использовании конечных разностей или объемов – метод сразу возвращает значение для поданного набора данных. Более того, метод является бессеточным, а значит может быть получено значение параметра в абсолютно любой точке расчетной области. Можно считать модель PINN своеобразным методом шифрования данных о процессе в конкретной области в виде весов обученной модели.

Функция потерь в PINN обычно состоит из нескольких слагаемых:

  1. Потеря на уравнении

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

    \[ f(x,y,t;\theta) = \hat{u}_t - \alpha \left( \hat{u}_{xx} + \hat{u}_{yy} \right). \]

    В случае отсутствия источника тепла в области, как в поставленной задаче, \( \lim\limits_{\theta \to \hat{\theta}} f(x,y,t;\theta) \to 0\), где \(\hat{\theta} \) набор параметров, при котором ошибка на уравнении равна нулю (абсолютное соответствие). Функция потерь для остатка может быть записана как среднеквадратичная ошибка:

    \[ \mathcal{L}_{\text{PDE}} = \frac{1}{N_f} \sum_{i=1}^{N_f} \left| f(x_i, y_i, t_i;\theta) \right|^2. \]

  2. Потеря на граничных условиях

    Для точек, принадлежащих границе \( \partial \Omega \), модель должна удовлетворять граничным условиям, например, \( \hat{u}(x,y,t) = g(x,y,t) \). Потеря граничных условий:

    \[ \mathcal{L}_{\text{BC}} = \frac{1}{N_b} \sum_{i=1}^{N_b} \left| \hat{u}(x_i, y_i, t_i;\theta) - g(x_i, y_i, t_i) \right|^2. \]

  3. Потеря на начальных условиях

    Если заданы начальные условия \( \hat{u}(x,y,0) = h(x,y) \), то:

    \[ \mathcal{L}_{\text{IC}} = \frac{1}{N_i} \sum_{i=1}^{N_i} \left| \hat{u}(x_i, y_i, 0;\theta) - h(x_i, y_i) \right|^2. \]

Общая функция потерь имеет вид:

\[ \mathcal{L}(\theta) = w_{\text{PDE}}\,\mathcal{L}_{\text{PDE}} + w_{\text{BC}}\,\mathcal{L}_{\text{BC}} + w_{\text{IC}}\,\mathcal{L}_{\text{IC}}, \]

где \( w_{\text{PDE}} \), \( w_{\text{BC}} \) и \( w_{\text{IC}} \) – весовые коэффициенты, регулирующие вклад каждого слагаемого.

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

Далее проводится обучение модели, при котором параметры \( \theta \) оптимизируются с использованием выбранного алгоритма оптимизации таким образом, чтобы минимизировать общую функцию потерь \( \mathcal{L}(\theta) \).

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

3.1.  Модификации метода PINN#

Для повышения скорости сходимости и общей точности решения для моделей PINN применяются следующие доработки:

Использование нестандартных функций активации

Стандартные функции активации (например, Tanh или ReLU) могут быть недостаточно эффективными для аппроксимации сложных функций с высокочастотными компонентами. Исходя из того, что множество решений уравнений описывается в виде гармонических функций (или рядов по ним), было предложено использовать функцию sin в качестве активационной. Пример такого слоя:

\[ \text{SinLayer}(x) = \sin\Big(\omega_0 \, (W x + b)\Big), \]

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

Гибридные оптимизаторы

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

  • Adam — градиентный метод первого порядка, он быстро сходится в начале обучения за счёт адаптивного шага.

  • LBFGS — квазиньютоновский метод, он часто используется на финальных этапах для тонкой настройки параметров, что позволяет достичь более высокой точности за счёт второго порядка информации о функции потерь.

Функция потерь минимизируется сначала оптимизатором Adam[1], а затем модель дообучается оптимизатором LBFGS[2] для повышения точности[3].

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

4.  Построение модели PINN#

Примечание

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

Также вы можете скачать файл PINN модели и сразу использовать её в проекте на платформе REPEAT VISION с помощью блока Аппроксимационная модель.

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

Для написания модели будет использован язык Python и фреймворк PyTorch, поэтому необходимы базовые знания этого языка и навыки работы с PyTorch.

Важно

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

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

Импорт библиотек#
1import torch
2import torch.nn as nn
3import torch.optim as optim
4
5import numpy as np
6import matplotlib.pyplot as plt
7
8import onnx
9import json

Для проверки обученной модели получем функцию аналитического решения. Она позволит как посчитать ошибку по \( L^2 \) норме в конкретный момент времени, так и провести сравнительный анализ эволюции решения полученного с помощью PINN и аналитическим решением.

Аналитическое решение#
 1def analytical_solution(x: np.ndarray, 
 2                        y: np.ndarray, 
 3                        t: np.ndarray,
 4                        bounds: dict,
 5                        alpha: float) -> np.ndarray:
 6    """
 7    Аналитическое решение
 8
 9    Параметры:
10        x: np.ndarray
11            Вектор координат x.
12        y: np.ndarray
13            Вектор координат y.
14        t: np.ndarray
15            Вектор координат t.
16        bounds: dict
17            Границы расчетной области.
18        alpha: float
19            Коэффициент температурапроводности материала.
20
21    Возвращает:
22        np.ndarray
23            Вектор значений температур, являющихся аналитическим решением.
24    """
25    u = (np.sin(np.pi * x/bounds['x']) * np.sin(np.pi * y/bounds['y'])
26         * np.exp(-2 * np.pi**2 * alpha * t/(bounds['x']*bounds['y'])))
27    
28    return u

  • Строки 1-5. Определение функции analytical_solution, которая вычисляет аналитическое решение уравнения теплопроводности в двумерной области.

  • Строка 25-26. Вычисление аналитического решения согласно формуле (5) с учетом размерности:

    \[ u(x, y, t) = \sin\left(\frac{\pi x}{L_x}\right) \sin\left(\frac{\pi y}{L_y}\right) \exp\left(-2 \pi^2 \alpha \frac{t}{L_x L_y}\right). \]

    где \( L_x \) и \( L_y \) — размеры области, берущиеся из bounds['x'] и bounds['y'] соответственно.

Функция возвращает массив температурных значений в зависимости от координат, времени и материала.

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

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

Слой с функцией активации sin#
 1class SinLayer(torch.nn.Module):
 2    """
 3    Слой нейронной сети с синусоидальной функцией активации.
 4
 5    Параметры:
 6        in_features: int
 7            Количество входных признаков.
 8        out_features: int
 9            Количество выходных признаков.
10        bias: bool
11            Использовать ли смещение в линейном слое.
12        is_first: bool
13            Является ли слой первым в сети.
14        omega_0: float
15            Частота для функции синуса.
16    """
17    def __init__(self, 
18                 input: float, 
19                 output: float, 
20                 bias: bool=True,
21                 is_first: bool=False, 
22                 omega_0: float=30):
23        super().__init__()
24        self.omega_0 = omega_0
25        self.is_first = is_first
26        self.in_features = input
27        self.linear = torch.nn.Linear(input, output, bias=bias)
28        self.init_weights()
29    
30    def init_weights(self):
31        """
32        Инициализация весов линейного слоя.
33        """
34        with torch.no_grad():
35            if self.is_first:
36                self.linear.weight.uniform_(-1 / self.in_features, 
37                                            1 / self.in_features)
38            else:
39                self.linear.weight.uniform_(-np.sqrt(6 / self.in_features) 
40                                            / self.omega_0, 
41                                            np.sqrt(6 / self.in_features) 
42                                            / self.omega_0)
43        
44    def forward(self, input):
45        """
46        Прямой проход через слой.
47
48        Параметры:
49            input: torch.Tensor
50                Входной тензор.
51
52        Возвращает:
53            torch.Tensor
54                Результат применения синусоидальной функции активации.
55        """
56        return torch.sin(self.omega_0 * self.linear(input))

  • Строка 1. Определение класса SinLayer, который представляет собой слой нейронной сети с синусоидальной функцией активации.

  • Строки 17-22. Определение конструктора класса __init__, который принимает параметры слоя.

  • Строки 23-27. Инициализация параметров omega_0, is_first, in_features и создание линейного слоя (torch.nn.Linear).

  • Строка 28. Вызов метода init_weights(), который задает начальные значения весов.

  • Строка 30. Определение метода init_weights(), который устанавливает начальные значения весов линейного слоя.

  • Строки 35-37. Если слой является первым (is_first = True), веса инициализируются равномерно в диапазоне \( [\frac{-1}{\text{in}\_\text{features}}, \frac{1}{\text{in}\_\text{features}}] \).

  • Строка 38-42. В противном случае веса инициализируются в диапазоне \( \left[-\frac{\sqrt{6 / \text{in}\_\text{features}}}{\omega_0}, \frac{\sqrt{6 / \text{in}\_\text{features}}}{\omega_0}\right] \), что улучшает сходимость при обучении сети.

  • Строка 44. Определение метода forward(), который выполняет прямой проход через слой.

  • Строка 56. Применение линейного слоя к входным данным, умножение результата на частоту omega_0 и передача в функцию torch.sin(), реализующую синусоидальную активацию.

Подобные слой используется в архитектурах нейросетей, таких как SIREN (Sinusoidal Representation Networks), и позволяют эффективно аппроксимировать широкий спектр частот зависимостей в данных различных типов.

Теперь создадим архитектуру нейронной сети. Класс SinPINN представляет собой полносвязную сеть, состоящую из нескольких линейных слоёв с синусоидальной функцией активации(экземпляры класса SinLayer). Модель принимает на вход координаты \( (x, y) \) и время \( t \) и выдаёт предсказанную температуры \( u \) точки, с координатами \( (x, y) \) в момент времени \( t \).

Нейронная сеть PINN#
 1class SinPINN(nn.Module):
 2    """
 3    Физически-информированная нейронная сеть (PINN) с функцией активации синус.
 4
 5    Параметры:
 6        layers: nn.ModuleList
 7            Список слоёв нейронной сети.
 8    """
 9    def __init__(self, 
10                 layers: list):
11        super().__init__()
12        self.layers = nn.ModuleList()        
13        for i in range(len(layers) - 1):
14            if i == 0:
15                self.layers.append(SinLayer(layers[i], 
16                                            layers[i+1], is_first=True))
17            elif i == len(layers) - 2:
18                self.layers.append(nn.Linear(layers[i], 
19                                             layers[i+1],))
20            else:
21                self.layers.append(SinLayer(layers[i], 
22                                            layers[i+1], is_first=False))
23    
24    def forward(self, 
25                x: torch.Tensor, 
26                y: torch.Tensor, 
27                t: torch.Tensor) -> torch.Tensor:
28        """
29        Прямой проход через сеть.
30
31        Параметры:
32            x: torch.Tensor
33                Тензор координат x.
34            y: torch.Tensor
35                Тензор координат y.
36            t: torch.Tensor
37                Тензор время t.
38
39        Возвращает:
40            torch.Tensor
41                Выход сети.
42        """
43        inputs = torch.cat([x.reshape(-1, 1), 
44                            y.reshape(-1, 1), 
45                            t.reshape(-1, 1)], dim=1)
46        for layer in self.layers:
47            inputs = layer(inputs)
48        return inputs

  • Строка 1: Объявление класса SinPINN, который наследует функциональность от базового класса nn.Module библиотеки PyTorch. Это обеспечивает возможность использовать все встроенные механизмы для построения и обучения нейронных сетей.

  • Строка 9: Объявление конструктора класса, который принимает layers — список определяющий количество нейронов для каждого слоя сети.

  • Строка 11: Вызов конструктора базового класса nn.Module для корректной инициализации всех внутренних механизмов PyTorch, необходимых для работы с нейронными сетями.

  • Строка 12: Инициализация объекта ModuleList для хранения слоёв нейронной сети. Использование ModuleList позволяет автоматически регистрировать параметры всех добавляемых слоёв.

  • Строки 13–22: Цикл по индексам списка layers создающий слои нейросети. Отдельно проверяется, является ли слой первым в сети. Если это условие выполнено, создаётся SinLayer с флагом is_first=True. Если же встречается последний слой — создаётся обычный Linear слой без синусоидальной функции активации. В остальных случаях создаётся SinLayer с is_first=False.

  • Строка 24: Объявление метода forward, который задаёт порядок прохождения входных данных через сеть. Этот метод принимает на вход тензоры x, y и t и возвращает предсказание модели (температуру).

  • Строки 43-45: Объединение входных данных (координаты x, y и время t) в один тензор. Данные преобразуются из векторов в 2D-тензоры размерности (-1, 1) и конкатенируются вдоль оси dim=1, то есть вдоль оси столбцов. В итоге получится матрица, в строках которой записаны определенные состояния системы \( x_i, y_i, t_i \).

Примечание

x.reshape(-1, 1):

  • -1 означает, что PyTorch автоматически подставит количество строк, равное количеству элементов в x.

  • 1 задает количество столбцов, превращая одномерный тензор (N,) в двумерный (N,1).

  • Строки 46–47: Цикл для прохождения данных через все слои.

Внимание

Обратите внимание что функция активации в данном случае интегрирована в слои.

  • Строка 48: Возврат выходного тензора inputs, содержащего предсказанные значения, полученные в результате прямого прохода через всю сеть.

4.1.  Определение функции потерь#

Функция потерь включает три компонента:

  • Потери на начальных условиях — сравнение предсказанного значения температуры в момент времени \( t=0 \).

  • Потери на граничных условиях — сравнение предсказанного значения температуры в точках \( (x, y) \in \partial \Omega \), где \( \Omega = [0,L] \times [0,L] \).

  • Потери для уравнения — невязка, посчитанная по дифференциальному уравнению \( u_t - \alpha (u_{xx} + u_{yy}) = 0 \).

Для начала рассмотрим компонент потерь для уравнения и напишем функцию для его расчета.

Функция потерь#
 1def pde_residual(model: SinPINN, 
 2                 x: torch.Tensor, 
 3                 y: torch.Tensor, 
 4                 t: torch.Tensor, 
 5                 alpha: float):
 6    """
 7    Вычисляет остаток по уравнению теплопроводности для модели PINN.
 8    Вид уравнения:
 9        u_t - alpha * (u_xx + u_yy) = 0
10
11    Параметры:
12        model: torch.nn.Module
13            Модель, вычисляющая u(x, y, t).
14        x: torch.Tensor
15            Тензор координаты x.
16        y: torch.Tensor
17            Тензор координаты y.
18        t: torch.Tensor
19            Тензор времени t.
20        alpha: float
21            Коэффициент теплопроводности.
22
23    Возвращает:
24        torch.Tensor
25            Остаток уравнения f = u_t - alpha*(u_xx + u_yy).
26    """
27    u = model(x, y, t)
28
29    u_x, u_y, u_t = torch.autograd.grad(u, [x, y, t], 
30                                        grad_outputs=torch.ones_like(u), 
31                                        create_graph=True)
32    
33    u_xx = torch.autograd.grad(u_x, x, grad_outputs=torch.ones_like(u_x), 
34                               create_graph=True)[0]
35    u_yy = torch.autograd.grad(u_y, y, grad_outputs=torch.ones_like(u_y), 
36                               create_graph=True)[0]
37    
38    f = u_t - alpha * (u_xx + u_yy)
39    return f

  • Строка 1. Объявление функции pde_residual, которая принимает в качестве аргументов модель, тензоры точек по координатам и времени для вычисления остатков (невязки) уравнения. Функция возвращает итоговое значение функции потерь на уравнении.

  • Строка 27. Вычисление значения выхода модели \( u(x, y, t) \).

  • Строка 29-31. Вычисление первых производных \( u_x, u_y, u_t \) с помощью функции автодифференцирования torch.autograd.grad. Параметр grad_outputs=torch.ones_like(u) указывает, что производная рассчитывается для каждого элемента тензора u. Флаг create_graph=True позволяет строить вычислительный граф для дальнейшего дифференцирования, что необходимо для вычисления второй производной.

  • Строка 33-35. Вычисление вторых производных \( u_{xx} \) и \( u_{yy} \) по x и y соответственно, используя torch.autograd.grad аналогичным образом.

  • Строка 38. Вычисление остатка уравнения теплопроводности. Если предсказания модели удовлетворяет физическому уравнению (2), \( f = u_t - \alpha (u_{xx} + u_{yy}) \) должно стремиться к нулю.

  • Строка 39. Возвращение результата f, который затем будет использован для обучения PINN.

4.2.  Подготовка данных и обучение модели#

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

Точки коллокации для уравнения#
 1def generate_pde_points(n: int, 
 2                        bounds: dict, 
 3                        device: torch.device):
 4    """
 5    Генерирует точки коллокации для уравнения PDE.
 6
 7    Параметры:
 8        n: int
 9            Число точек коллокации.
10        bounds: dict
11            Границы области.
12        device: torch.device
13            Устройство (CPU или GPU).
14
15    Возвращает:
16        Tuple[torch.Tensor, torch.Tensor, torch.Tensor]
17            Тензоры x, y, t с requires_grad=True.
18    """
19    x = (bounds['x'] * torch.rand(n, 1, device=device)).requires_grad_()
20    y = (bounds['y'] * torch.rand(n, 1, device=device)).requires_grad_()
21    t = (bounds['t'] * torch.rand(n, 1, device=device)).requires_grad_()
22    return x, y, t

  • Строка 1. Объявление функции generate_pde_points, которая генерирует точки коллокации для подсчета остатков решения дифференциального уравнения.

  • Строка 19-21. Генерация случайных значений для координат x, y, t в пределах указанных границ bounds в два этапа:

    • Создает случайные значения от 0 до 1 с помощью torch.rand(n, 1, device=device).

    • Умножение на bounds['x'], bounds['y'], bounds['t'] масштабирует точки в заданные границы расчетной области.
      вызов метода requires_grad_() позволяет вычислять градиенты по этим переменным.

  • Строка 22. Возвращение сгенерированных тензоров x, y, t.

Точки коллокации для начальных условий#
 1def generate_initial_condition(n: int, 
 2                               bounds: dict, 
 3                               device: torch.device):
 4    """
 5    Генерирует точки для начального условия.
 6
 7    Параметры:
 8        n: int
 9            Число точек начального условия.
10        bounds: dict
11            Границы расчетной области.
12        device: torch.device
13            Устройство (CPU или GPU).
14
15    Возвращает:
16        Tuple[torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor]:
17            Тензоры x, y, t, u, 
18            где t — нулевой тензор, а u = sin(pi * x/L) * sin(pi * y/L).
19    """
20    x = bounds['x'] * torch.rand(n, 1, device=device)
21    y = bounds['y'] * torch.rand(n, 1, device=device)
22    t = torch.zeros_like(x, device=device)
23    u = torch.sin(np.pi * x/bounds['x']) * torch.sin(np.pi * y/bounds['y'])
24    return x, y, t, u

  • Строка 1. Объявление функции generate_initial_condition, создающей точки для проверки начального условия.

  • Строки 20-21. Генерация случайных значений для x и y в пределах границ, аналогично способу в функции generate_pde_points.

  • Строка 22. Установка t = 0, так как начальные условия задаются в момент времени \( t=0 \).

  • Строка 23. Вычисление u, основываясь на поставленном начальном условии (3) \( u = \sin(\pi x/L) \sin(\pi y/L) \).

  • Строка 24. Возвращение сгенерированных тензоров x, y, t, u.

Точки коллокации для граничных условий#
 1def generate_boundary_conditions(n: int, 
 2                                 bounds: dict, 
 3                                 device: torch.device):
 4    """
 5    Генерирует точки граничных условий.
 6
 7    Параметры:
 8        n_b: int 
 9            Число точек для граничных условий.
10        bounds: dict
11            Границы области, например, {'x': 1.0, 'y': 1.0, 't': 1.0}.
12        device: torch.device
13            Устройство (CPU или GPU).
14
15    Возвращает:
16        Tuple[torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor]:
17            x_b, y_b, t_b, u_b, где u_b = 0.
18    """
19    x = bounds['x'] * torch.rand(n, 1, device=device)
20    y = bounds['y'] * torch.rand(n, 1, device=device)
21    t = bounds['t'] * torch.rand(n, 1, device=device)
22
23    border_selector = torch.randint(0, 4, (n, 1), device=device)
24    x[border_selector == 0] = 0.0 # Left bound
25    x[border_selector == 1] = bounds['x'] # Right bound
26    y[border_selector == 2] = 0.0 # Bottom bound
27    y[border_selector == 3] = bounds['y'] # Upper bound
28
29    u_b = torch.zeros_like(t, device=device)
30    return x, y, t, u_b

  • Строка 1. Объявление функции generate_boundary_conditions, создающей точки для проверки граничных условий.

  • Строки 19-21. Генерация случайных точек x, y, t в пределах границ, аналогично способу в функции generate_pde_points.

  • Строка 23. Создание border_selector, случайно выбирающего граничные точки.

  • Строки 24-27. Присвоение граничных значений координат x=0, x=L, y=0, y=L в зависимости от border_selector.

  • Строка 29. Установка нулевого граничного условия u_b = 0 в соответствии с условием (4).

  • Строка 30. Возвращение сгенерированных тензоров x, y, t, u_b.

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

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

Параметры задачи#
 1torch.manual_seed(42)
 2np.random.seed(42)
 3
 4device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
 5
 6N_F, N_B, N_I = 10000, 5000, 2000
 7BOUNDS = {'x': 0.1, 'y': 0.1, 't': 12.0}
 8ALPHA = 1.12 * 1e-4
 9
10model = SinPINN([3, 64, 1]).to(device)
11mse_loss = nn.MSELoss()
12optimizer = optim.Adam(model.parameters(), lr=1e-3)
13
14x_f, y_f, t_f = generate_pde_points(N_F, BOUNDS, device)
15x_i, y_i, t_i, u_i = generate_initial_condition(N_I, BOUNDS, device)
16x_b, y_b, t_b, u_b = generate_boundary_conditions(N_B, BOUNDS, device)
17
18w_pde, w_bc, w_ic = 10.0, 1.0, 1.0

  • Строка 1. Установка torch.manual_seed(42), чтобы обеспечить воспроизводимость эксперимента при использовании PyTorch.

  • Строка 2. Установка np.random.seed(42), чтобы обеспечить воспроизводимость для библиотеки NumPy.

  • Строка 4. Определение вычислительного устройства device (GPU, если доступно, или CPU).

  • Строка 6. Определение количества точек коллокации для уравнения (N_F), для граничных условий (N_B), для начального условия (N_I).

  • Строка 7. Определение границ расчетной области по оси x, по оси y и максимального значения времени моделирования в виде словаря BOUNDS.

  • Строка 8. Определение коэффициента теплопроводности, соответствующего меди ALPHA = 1.12 * 1e-4.

  • Строка 10. Создание модели SinPINN с одним скрытым слоем в 64 нейрона и перемещение модели на устройство device.

  • Строка 11. Определение функции потерь mse_loss, использующей среднеквадратичную ошибку (MSE).

  • Строка 12. Определение оптимизатора Adam для обновления параметров модели.

  • Строки 14-16. Генерация обучающих данных (точек коллокации) для соответствующий частей постановки задачи.

  • Строка 18. Формальное определение весов для различных частей функции потерь. Они могут быть скорректированы для лучшей сходимости или учета определенных значений.

Перейдем к обучению модели. Основной цикл обучения проходит в течение заданного количества эпох, где на каждой эпохе:

  • Обнуляются градиенты.

  • Вычисляются компоненты функции потерь.

  • Производится обратное распространение ошибки.

  • Обновляются параметры модели.

  • Сохраняется история значений компонент функции потерь.

  • Периодически выводится значение потерь для контроля процесса обучения.

Цикл обучения#
 1epochs = 20000
 2
 3loss_history = {'pde': [],
 4                'ic': [],
 5                'bc': []}
 6
 7for epoch in range(epochs+1):
 8    optimizer.zero_grad()
 9    
10    u_pred_i = model(x_i, y_i, t_i)
11    loss_ic = mse_loss(u_pred_i, u_i)
12    
13    u_pred_b = model(x_b, y_b, t_b)
14    loss_bc = mse_loss(u_pred_b, u_b)
15    
16    loss_pde = mse_loss(pde_residual(model, x_f, y_f, t_f, ALPHA), 
17                        torch.zeros_like(x_f))
18    
19    loss = w_ic * loss_ic + w_bc * loss_bc + w_pde * loss_pde
20    
21    loss.backward()
22    optimizer.step()
23
24    loss_history['ic'].append(loss_ic.item())
25    loss_history['bc'].append(loss_bc.item())
26    loss_history['pde'].append(loss_pde.item())
27    
28    if epoch % 1000 == 0:
29        print(f"Epoch {epoch}, Loss: {loss.item():.3e}, "
30              f"IC: {loss_ic.item():.3e}, "
31              f"BC: {loss_bc.item():.3e}, "
32              f"PDE: {loss_pde.item():.3e}")

  • Строка 1. Устанавливается количество эпох, определяющее число итераций обучения модели.

  • Строки 3-5. Создается словарь loss_history для хранения значений компонент функции потерь.

  • Строки 7-30. Основной цикл обучения модели. В нем на каждой из эпох происходит:

    • Обнуление градиентов оптимизатора.

    • Вычисление ошибки для начального условия. Модель предсказывает значения u для точек начального условия, затем по ним вычисляется среднеквадратичная ошибка относительно истинных значений u.

    • Аналогичным образом вычисляется ошибка для граничных условий.

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

    • Все компоненты функции потерь складываются с соответствующими весам.

    • Вычисляются градиенты параметров модели с учетом всех составляющих потерь.

    • Выполняется шаг оптимизации, обновляющий параметры модели.

    • Далее происходит сохранение значений компонент функции потерь и периодический вывод.

Для визуальной оценки качества процесса обучения можно провести визуализацию значений компонент функции потерь от эпох в логарифмической шкале.

График функции потерь#
 1def plot_loss_history(loss_history: dict, 
 2                      step: int=200):
 3    """
 4    Функция для построения графика функции потерь.
 5
 6    Параметры:
 7        loss_history : dict
 8            Словарь, содержащий данные для потерь.
 9        step : int, optional
10            Шаг выборки данных для отображения, по умолчанию 100.
11    """
12    epochs = range(0, len(loss_history['pde']), step)
13
14    plt.plot(epochs, loss_history['pde'][::step], label='PDE Loss')
15    plt.plot(epochs, loss_history['bc'][::step], label='BC Loss')
16    plt.plot(epochs, loss_history['ic'][::step], label='IC Loss')
17    
18    plt.yscale('log')
19    plt.xlabel('Epoch')
20    plt.ylabel('Loss')
21    plt.grid()
22    plt.legend()
23    plt.show()

  • Строки 1-2. Объявление функции plot_loss_history, предназначенная для построения графика функции потерь.

  • Строка 12. Определяет список epochs, который содержит номера эпох, с шагом step. Нужно для корректного отображения номера эпохи по оси \( x \).

  • Строки 14-16. Строятся графики для каждой из компонент потерь, при этом значения берутся с шагом step для большей гладкости.

Вызовем эту функцию, передав ей историю обучения loss_history.

Построение графика функции потерь#
1plot_loss_history(loss_history)

Будет построен график следующего вид:

../../../_images/loss_function-05c2d8c4.png

График функции потерь.#

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

Из-за несколько отличающегося интерфейса взаимодействия в реализации метода LBFGS в библиотеке PyTorch от остальных оптимизаторов, необходимо объявить вспомогательную функцию closure, которая будет вызываться при работе оптимизатора.

Дообучение модели#
 1def closure():
 2    """
 3    Вспомогательная функция для оптимизатора LBFGS.
 4    """
 5    optimizer_lbfgs.zero_grad()
 6
 7    loss_ic = mse_loss(model(x_i, y_i, t_i), u_i)
 8
 9    loss_bc = mse_loss(model(x_b, y_b, t_b), u_b)
10
11    loss_pde = mse_loss(pde_residual(model, x_f, y_f, t_f, ALPHA), 
12                           torch.zeros_like(x_f))
13    
14    loss = w_ic * loss_ic + w_bc * loss_bc + w_pde * loss_pde
15
16    loss.backward()
17    
18    return loss
19
20optimizer_lbfgs = optim.LBFGS(model.parameters(), lr=1.0, 
21                              max_iter=500, history_size=50)
22optimizer_lbfgs.step(closure)

  • Строка 1-18. Создается функция closure, которая будет использоваться для оптимизации с использованием метода LBFGS. В ней происходят следующие действия:

    • Обнуление градиентов оптимизатора.

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

    • Полученные компоненты суммируются с весами, получая итоговую функцию потерь loss.

    • Вычисляются градиенты модели с помощью вызова loss.backward().

  • Строки 20-21. Создается объект оптимизатора LBFGS, в который предаются параметры модели, шаг обучения, максимальное число итераций и размером истории оптимизации.

  • Строка 22. Выполняется шаг оптимизации с использованием функции closure, обновляя параметры модели для минимизации функции потерь.

Таким образом мы получили обученную модель для решения поставленной задачи.

5.  Экспорт модели в формат ONNX#

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

Экспорт модели в ONNX#
 1model.double()
 2
 3# Создание тестового входного тензора (x, y, t)
 4x_temp = torch.randn(1, 1, device=device, dtype=torch.float64) * BOUNDS['x']
 5y_temp = torch.randn(1, 1, device=device, dtype=torch.float64) * BOUNDS['y']
 6t_temp = torch.randn(1, 1, device=device, dtype=torch.float64) * BOUNDS['t']
 7
 8input_names = ['x', 'y', 't']
 9output_names = ['u']
10
11torch.onnx.export(model,
12                  (x_temp, y_temp, t_temp),
13                  "pinn_heat2d.onnx",
14                  opset_version=12,
15                  input_names=input_names,
16                  output_names=output_names)

  • Строка 1. Перевод модели в режим двойной точности float64. Это важно для точных вычислений при работе с уравнениями в частных производных.

  • Строки 4-6. Создание тестовых входных тензоров x_temp, y_temp, t_temp со случайными значениями, масштабированными по границам BOUNDS.

  • Строки 8-9. Определение имен входных \( (x, y, t) \) и выходного \( (u) \) параметров модели.

  • Строки 11-16. Экспорт модели в формат ONNX с параметрами:

    • model — объект класса SinPINN.

    • (x_temp, y_temp, t_temp) — кортеж входных тензоров.

    • «pinn_heat2d.onnx» — имя выходного файла.

    • opset_version=12 — версия ONNX API.

    • input_names и output_names — имена входных и выходного параметров.

Дополнительно, для интеграции в REPEAT VISION нужно добавить метаданные о границах расчетной области.

:caption: Добавление метаданных в файл
:lineno-start: 1
# Загрузка экспортированной модели ONNX
onnx_model = onnx.load("pinn_heat2d.onnx")

# Создание словаря метаданных с границами параметров x, y, t
metadata_dict = {
    "x": f"{BOUNDS['x']}",
    "y": f"{BOUNDS['y']}",
    "t": f"{BOUNDS['t']}"
}
metadata_json = json.dumps(metadata_dict)

metadata_prop = onnx.StringStringEntryProto()
metadata_prop.key = "bounds"
metadata_prop.value = metadata_json

# Добавление метаданных и сохранение обновленной модели
onnx_model.metadata_props.append(metadata_prop)
onnx.save(onnx_model, "pinn_heat2d.onnx")

  • Строка 2. Загрузка ранее экспортированной модели ONNX.

  • Строки 5-9. Создание словаря metadata_dict, содержащего границы входных параметров.

  • Строка 10. Конвертация словаря в JSON-формат.

  • Строки 12-14. Создание объекта StringStringEntryProto для хранения пар «ключ-значение».

  • Строка 17. Добавление метаданных в модель ONNX.

  • Строка 18. Сохранение обновленной модели в файл «pinn_heat2d.onnx».

6.  Создание проекта на платформе REPEAT VISION#

Для начала необходимо создать новый проект в REPEAT VISION и присвоить ему название, например, 2D теплопроводность с PINN. После создания проекта установите параметр Время моделирования равным 12 [с] в соответствии с использованным ранее при обучении модели.

Time_of_modeling

Рис. 1 Время моделирования#

6.1.  Импорт модели в формте ONNX#

В первую очередь, импортируем обученную модель, сохранённую в формате ONNX, в проект REPEAT VISION. Для этого откройте библиотеку Внешние модели и перенесите на рабочую область блок Аппроксимационная модель.

Block_init

Рис. 2 Инициализация блока Аппроксимационная модель#

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

Init_done

Рис. 3 Инициализация модели выполнена#

Важно

Обратите внимание что после импорта модели, количество входных и выходных портов было определено автоматически. Оно должно соответствовать заданным при экспорте модели в ONNX. Если это не так — в результате экспорта была допущена ошибка.

Inputs_and_outputs

Рис. 4 Входные и выходные порты#

В качестве координат могут быть поданы как любые значения из области моделирования, как в виде констант, так и каких-либо зависимостей. Создадим блоки координат, для простоты обозначим их постоянными значениями и поставим в центр расчетной области — точку \( (0.05, 0.05) \). Для этого добавим два блока Константа из библиотеки Автоматика и соединим их с соответствующими входами блока аппроксимации.

SimulationTime

Рис. 5 Добавление блоков Константа#

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

SimulationTime

Рис. 6 Добавление блока Время моделирования#

Для получения выходного значения модели воспользуемся блоком Выходной порт из библиотеки Внешние модели.

OutPort

Рис. 7 Добавление блока Выходной порт#

Добавим визуализацию выходного порта, для построения графика зависимости координаты от времени.

Plot_creation

Рис. 8 Создание графика#

7.  Запуск расчета#

Запустите расчёт кнопкой «Запустить моделирование»play. Расчёт должен успешно завершиться следующим сообщением:

Расчет остановлен со стороны решателя. Достигнут лимит модельного времени.

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

Results

Рис. 9 График результатов расчета#

8.  Выводы#

В рамках данного практического руководства вы:

  • Получили более подробное представление об использовании нейронных сетей для решения физических задач.

  • Применили актуальные техники по улучшению сходимости моделей PINN в виде синусоидальной функции активации и использования нескольких оптимизаторов.

  • Улучшили навыки создания и обучения моделей нейронных сетей с использованием PyTorch.

  • Освоили экспорт в формат ONNX и последующую интеграцию в проект на платформе REPEAT VISION.

  • Научились добавлять учет двумерных зависимостей в решатель REPEAT VISION.

9.  Список литературы#

На этой странице