Доклады Российской академии наук. Математика, информатика, процессы управления, 2023, T. 514, № 2, стр. 28-38

2.1. Нелинейное уравнение Шредингера

Распространение лазерного импульса в оптоволоконном кабеле в системе координат, движущейся с групповой скоростью импульса, описывается нелинейным уравнением Шредингера [3]:

где $A(z,T)$ – комплексная временная огибающая электрического поля, $z$ – пространственная координата в оптоволокне, ${{\beta }_{2}}$ – коэффициент дисперсии второго порядка, $T$ – временная развертка импульса, $\alpha $ – коэффициент поглощения, а $\gamma $ – коэффициент Керра. Заметим, что $A = u + iv$ является комплексной функцией, и поэтому уравнение (1) является системой из двух уравнений. Учитываются три физических эффекта, а именно затухание, дисперсия и самофазная модуляция через члены при $\alpha $, ${{\beta }_{2}}$ и $\gamma $ соответственно. Каждый из этих эффектов проявляет себя на расстоянии, большем определенного, называемого характеристической длиной. Можно ввести характеристическую длину затухания ${{L}_{a}}$, дисперсионную длину ${{L}_{d}}$и нелинейную длину ${{L}_{n}}$где ${{T}_{0}}$ и ${{P}_{0}}$ – начальные ширина и мощность в пике пульса соответственно. По соотношению этих длин можно оценить, насколько сильно будет проявляет себя соответствующий эффект. Например, если ${{L}_{n}} \gg {{L}_{d}}$ и $z \sim {{L}_{d}}$, то нелинейный эффект будет пренебрежимо мал по сравнению с дисперсионным эффектом. В данной работе без потери общности ограничимся анализом гауссовского импульса в качестве граничных условий: где импульс на входе оптоволокна, ${{T}_{0}}$ – полуширина импульса, а $[ - {{T}_{{{\text{max}}}}}$, ${{T}_{{{\text{max}}}}}]$ – временное окно. Начальные условия имеют вид

2.2. Замена переменных

При решении уравнения (1) часто прибегают к преобразованию системы координат [14]:

где и

Также $A$ нормируется на пиковую мощность ${{P}_{0}}$:

Уравнение для нормированного сигнала $U$ в преобразованной системе координат $(\xi ,t)$ имеет вид:

где

В данной работе расчеты проводятся для длины волны $\lambda = 1030$ нм и коэффициента ${{\beta }_{2}}$ = = 0.019164. Параметр поглощения $\alpha $ принят равным нулю. Временной интервал T_max составляет 100 пс. Параметр Керра $\gamma $ равен 0.0011465. Энергия импульса и длина волокна варьируются в диапазонах 1e-6 – 10 нДж и 5 – 10 000 м соответственно.

2.3. Метод PINN

Сформулируем метод PINN к нашей задаче. PINN является нейронной сетью, на вход которой подаются значения z и t. На выходе нейронной сети имеем действительную $u$ и мнимую ${v}$ части $A(z,t)$:

В процессе обучения минимизируется функция потерь, имеющая вид

где ${{E}_{{pde}}}$ – невязка уравнения, рассчитанная в ${{N}_{{pde}}}$ точках расчетной области, а ${{E}_{b}}$ и ${{E}_{{in}}}$ – невязки граничных и начальных условий, соответственно, т.е. где ${{N}_{b}}$ и ${{N}_{{in}}}$ – количества точек, взятых на z = 0 и $T = \pm {{T}_{{{\text{max}}}}}$ соответственно. Частные производные $\frac{{\partial {\text{PINN}}({{z}_{j}},{{T}_{j}})}}{{\partial z}}$ и $\frac{{{{\partial }^{2}}{\text{PINN}}({{z}_{j}},{{T}_{j}})}}{{\partial {{T}^{2}}}}$ рассчитываются с помощью автоматического дифференцирования. Минимизация функции E производится по N = = ${{N}_{{pde}}} + {{N}_{b}} + {{N}_{{in}}}$ точкам $\{ {{z}_{j}},{{t}_{j}}\} $, $j = 1,...,N$ из расчетной области. Структура PINN показана на рис. 1.

Эффективность PINN зависит от количества слоев и нейронов в сети, функций активации, выбора точек, в которых производится расчет функции потерь. Кроме того, как отмечено в работе [15], точность метода можно повысить с помощью весов при слагаемых функции потерь:

где ${{\lambda }_{{pde}}}$, ${{\lambda }_{b}}$ и ${{\lambda }_{{in}}}$ – соответствующие веса.

3.1. Конфигурация PINN

Мы разработали модель PINN на основе PyTorch. Нейронная сеть может включать линейные и ResNet-блоки с разными количеством нейронов и функциями активации. В работе [16] предлагается использовать оптимизатор Adam совместно с алгоритмом L-BFGS. Однако в наших расчетах мы обнаружили, что алгоритм Adamax имеет схожую эффективность по сравнению с комбинацией двух вышеупомянутых минимизаторов. Во всех рассмотренных сценариях функция активации Tanhshrink работает лучше, чем функции активации ReLU и Tanh. Поэтому, чтобы упростить анализ, во всех расчетах мы используем по умолчанию Adamax и Tanhshrink.

Фактическое значение скорости обучения минимизатора регулируется планировщиком. Он умножает текущее значение скорости обучения на коэффициент (по умолчанию 0.9), если функция потерь не уменьшается после определенного количества итераций (по умолчанию 30 итераций). Минимальное значение скорости обучения установлено на уровне 1e-7. Частные производные вычисляются с помощью встроенных функций автоматического дифференцирования PyTorch. Веса сети инициализируются с помощью метода Ксавье. Обучение PINN производится на GPU NVidia Tesla V100 16GB. Если не указано иное, время обучения установлено на уровне двух часов для одного сценария. В модели могут быть использованы преобразования (8)–(11) и нормировка (12), а также балансировка слагаемых в функции потерь (19). Для оценки точности PINN используются референсные решения, полученные с помощью метода разделения шагов. Методом сеточного поиска было найдено, что для большинства сценариев наилучшие результаты PINN получаются при использовании конфигурации с 8 ResNet блоками по 150 нейронов в каждом блоке.

Количество точек сэмплирования существенно влияет на время обучения и точность модели. Время обучения одной эпохи линейно зависит от количества точек выборки. На рис. 2, иллюстрирующем функцию потерь, можно выделить области плавного снижения функции потерь, резкого снижения и области сильных колебаний с последующим плато. Увеличивая количество точек с 240 000 до 720 000, мы получаем улучшение результатов, а именно: потери уменьшаются, спад происходит раньше, обучение имеет тенденцию заканчиваться в стабильной точке до начала плато. Однако слишком большое количество точек (в нашем случае 1 040 000) приводит к негативным последствиям, а именно: спад наступает позже и обучение прекращается в области сильных колебаний функции потерь. Таким образом, при оптимальном количестве точек выборки процедура обучения останавливается, когда функция потерь достигает плато. В представленных расчетах используются 720 000–960 000 точек.

3.2. Применение PINN к нелинейному уравнению Шредингера

Применим описанную конфигурацию PINN к уравнению (1). Результаты показаны в табл. 1. Заметим, что с увеличением длины оптоволокна, а также с ростом энергии импульса теряется сходимость. С физической точки зрения это можно объяснить тем, что с ростом этих двух параметров решение для действительной и мнимой частей становится все более осциллирующим. А согласно частотному принципу глубокого обучения (F-principle) [17], наличие высоких частот в решении увеличивает время обучения. Недавно в работе [18] был обнаружен аналогичный эффект. Это наблюдение позволяет предположить, что существует максимальный уровень колебаний (приводящий к верхним ограничениям на длину волокна), выше которого PINN не сходится. В нашем случае увеличение числа слоев и количества нейронов не приводило к сходимости в проблемных случаях.

3.3. Влияние замены переменных на сходимость

Далее, применим PINN к уравнению (13), полученному с помощью преобразований (8)–(11) и нормализации (12). В этом случае сходимость получена для мощности 1 мВт при длине волокна 1 км (см. табл. 2 и рис. 3). Заметим, что точность решения PINN несколько ниже в тех областях, где наблюдаются сильные колебания, что опять же объясняется частотным принципом.

В табл. 3 сравнивается эффективность подхода PINN с эффективностью традиционных алгоритмов, а именно метода разделения шагов (split-step Foutier method – SSFM) и его распараллелиной модификации (massive parallel algorithm – MPA). Детали реализации описаны в работе [19]. Расчеты выполнены с шагом интегрирования 5 × 10^–4 м. Несмотря на то что цифры, приведенные в табл. 3, являются ориентировочными, можно сделать вывод, что время обучения PINN сопоставимо с временем вычислений традиционных методов, хотя и не превосходит результаты SSFM и MPA.

В линейном случае (в котором нелинейным слагаемым при коэффициенте $\gamma $ пренебрегается) PINN сходится для всех энергий при длине оптического волокна до 2 км, как показано в табл. 4. При этом, несмотря на то, что функция потерь имеет относительно малые значения, ошибки в решении могут превышать 5%. В этом случае косвенным фактором, позволяющим оценить точность решения, является выполнение закона сохранения энергии импульса.

Сравнивая табл. 2 и 4, мы видим, что максимальные длины оптического волокна, при которых PINN еще сходится, различны при энергиях 10 и 0.1 нДж. Для понимания этого факта приведем численные значения ${{L}_{n}}$ и ${{L}_{d}}$ в табл. 5. Заметим, что если ${{L}_{n}}$ меньше ${{L}_{d}}$, то именно ${{L}_{n}}$ определяет максимальную длину волокна, которая зависит от энергии. Если ${{L}_{d}}$ < ${{L}_{n}}$, то ${{L}_{d}}$ определяет длину волокна, которая в этом случае не зависит от энергии. Эмпирическим путем мы получили, что область сходимости PINN по длине волокна определяется как $ \sim {\kern 1pt} \min \left( {50{{L}_{d}},50{{L}_{n}}} \right)$. Отметим, что во всех моделях параметр затухания α = 0, что соответствует ${{L}_{a}} \to \infty $. При значениях $\alpha > 0$ сигнал затухает после прохождения расстояния $ \sim {\kern 1pt} {{L}_{a}}$ и тогда вопрос о сходимости PINN снимается.

Большая разница в коэффициентах при членах уравнения может вызвать численные трудности для PINN. Рассмотрим коэффициенты ${{k}_{1}}$, ${{k}_{2}}$, $b$ и $c$, которые получаются после замены переменных (8)–(12). Численные значения приведены в табл. 6, из которой следует, что решение сходится, если коэффициенты при членах NLSE отличаются менее чем на 9 порядков. В противном случае вклад слагаемого при малом коэффициенте в функцию потерь при расчете ${{E}_{{pde}}}$ становится пренебрежимо мал, и PINN не может корректно решить уравнение.

3.4. Синусоидальное отображение

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

применяется на выходе первого слоя. В наших расчетах мы видели, что для решения нелинейного уравнения Шредингера более точные результаты получаются, если разделить первый слой на узлы, относящиеся к $x$ и t, и применить (20) только к узлам t. Соответствующая структура PINN показана на рис. 4. При этом значения MSE уменьшаются в среднем в два раза. Здесь интересно провести аналогию с методом SSFM, в котором применяется преобразование Фурье к координате t. Ранее в работе [21] было предложено использовать блоки преобразования Фурье для повышения точности PINN. Таким образом, использование синусоидальной функции активации для t может быть частично оправдано, а при проектировании архитектуры нейронной сети можно руководствоваться математическими свойствами задачи, как предложено в работе [22].

3.5. Подбор весов к слагаемым функции потерь

Рассмотрим влияние коэффициентов при слагаемых функции потерь. Сразу заметим, что адаптивный алгоритм подборки весов, описанный в работе [15], не привел к существенному улучшению результатов. Поэтому мы рассмотрим подбор веса $\lambda $ перед невязкой по начальным условиям:

В процессе обучения мы меняем $\lambda $ так, чтобы ${{E}_{{pde}}}$ и $\lambda {{E}_{{in}}}$ были одного порядка, что позволяет не только уменьшить невязку по начальным условиям, но и в целом повысить точность. На рис. 5 показан результат совместной работы синусоидального отображения и подбора $\lambda $. Благодаря их использованию, PINN может воспроизводить области с высокочастотными колебаниями на границах окна, уменьшая на порядок средне-квадратичную ошибку решения.

3.6. Мета-PINN

Все вышеперечисленные PINN были обучены для одного сценария. Это означает, что при изменении любого параметра, например, энергии импульса, PINN необходимо переобучить. Этот процесс может занимать до двух часов на GPU для одного сценария. При такой скорости PINN не может превзойти SSFM, который занимает менее 1 ч на CPU (см. табл. 3). Однако PINN способен решать параметрические уравнения. Такой тип PINN часто называют мета-PINN [23]. В этом случае мета-PINN имеет дополнительные входы для параметров уравнения.

Рассмотрим задачу, в которой энергия входного импульса может изменяться в диапазоне от 1е-3 до 10е-3 нДж. Результаты для нескольких энергий приведены на рис. 6. С увеличением энергии импульса нелинейные эффекты становятся более выраженными, но это не мешает PINN предсказывать решение с приемлемой точностью. Отметим, что мета-PINN обладает хорошей обобщающей способностью. Хотя обучение проводилось только для энергий импульса k × 10^–3 нДж, где k – целое число, мета-PINN дает точные результаты для промежуточных значений энергии. Процедура обучения мета-PINN занимает примерно в 2–3 раза больше времени, чем PINN для одного сценария, в то время как предсказание одного сценария с помощью мета-PINN занимает лишь доли секунды. Это говорит о том, что мета-PINN может быть более эффективным, чем традиционные методы (в частности, метод разделения шагов), которые могут решать только один сценарий за раз.