CC BY
0
УДК 519.853.62+519.632.4+535.421 БСХ: 10.14529/mmp240205
0 ВЫЧИСЛЕНИИ ФУНКЦИОНАЛЬНОЙ ПРОИЗВОДНОЙ
ДЛЯ ОДНОЙ ЗАДАЧИ ОПТИМАЛЬНОГО УПРАВЛЕНИЯ
М.О. Корпусов1'2, М.В. Артемьева2
1 Российский университет дружбы народов, г. Москва, Российская Федерация
2Московский государственный университет, г. Москва, Российская Федерация
Задача синтеза многослойной дифракционной решетки формулируется как задача оптимального управления и заключается в минимизации целевого функционала, зависящего от геометрических параметров профиля решетки. Градиентный метод является наиболее надежным и стабильным методом решения этой задачи. В статье представлен метод вычисления функциональной производной (градиента) целевого функционала, который выполняется путем решения сопряженной задачи со специальными граничными условиями. Кроме того, в статье обсуждается численная реализация этого решения и расчет градиента. Также представлены результаты вычислительного эксперимента.
Ключевые слова: производная функционала; градиент; сопряженная задача; задача оптимального управления; задача синтеза; дифракционные решетки.
Введение
В современной лазерной технике, системах коммуникации, космических исследованиях и многих других областях науки и техники для управления электромагнитным излучением широко применяются многослойные дифракционные решетки [1, 2]. При этом в различных оптических системах к дифракционным решеткам предъявляются совершенно разные требования. Однако в большинстве случаев решетка должна обладать максимально возможной дифракционной эффективностью в определенных порядках дифракции. Получить требуемую характеристику можно за счет подбора формы профиля решетки, поскольку именно геометрия профиля в основном определяет эффективность, с которой свет дифрагирует в каждый из порядков [2]. Таким образом, чтобы найти оптимальные параметры дифракционной решетки перед ее непосредственным изготовлением, необходимо решить задачу синтеза дифракционной решетки.
С математической точки зрения задачи синтеза дифракционных решеток являются задачами оптимального управления и формулируются как задачи минимизации целевого функционала. Целевой функционал строится таким образом, чтобы его минимум соответствовал максимальному значению дифракционной эффективности при требуемом порядке дифракции на отражение или пропускание и при этом сам функционал зависел от геометрических параметров решетки - параметров управления. Для минимизации целевого функционала применяется градиентный метод, при котором градиент вычисляется с помощью решения сопряженной задачи [3]. Это наиболее устойчивый метод с точки зрения увеличения числа управляющих параметров, и его сходимость к оптимальному решению математически обоснована [4].
В данной работе обсуждается метод вычисления градиента целевого функционала при помощи решения сопряженной задачи с выбранными специальным образом граничными условиями.
1. Постановка задачи синтеза дифракционной решетки
Рассмотрим одномерную многослойную дифракционную решетку, состоящую непосредственно из решетки с прямоугольным профилем штриха, подложки и системы однородных слоев, расположенной между ними (последнее, вообще говоря, не является обязательным). Решетка считается бесконечной и периодической вдоль оси х с периодом й, причем ось г перпендикулярна границам слоев структуры. Период профиля разбивается на К равных между собой подпериодов = й/К, как показано на рис. 1. Каждому к-му подпериоду соответствует свой фактор заполнения /.
fidi
fkdk
x
z
d
d
k
h
n
"III
Рис. 1. Многослойная бинарная дифракционная решетка
Математическая постановка задачи синтеза формулируется как задача оптимального управления и заключается в минимизации целевого функционала
шш ^ (/), где / = (/ / ... Д) (1)
- набор управляющих параметров. Для рассматриваемой в работе многослойной дифракционной решетки набор управляющих параметров составляют факторы заполнения {/}, которые оптимизируются таким образом, чтобы решетка обладала максимально близкой к 100% дифракционной эффективностью в п-ом дифракционном порядке (в спектре отражения для отражающей решетки или пропускания для пропускающей) в некотором заданном диапазоне длин волн.
В общем виде целевой функционал представляет собой интегральную характеристику дифракционной решетки в отношении квадратичного отклонения дифракционной эффективности в заданном диапазоне длин волн от желаемого значения (100%) и зависит от управляющих параметров решетки. В численных расчетах интеграл заменяется суммой по некоторому набору дискретных длин волн Аг, при этом функционал имеет следующий вид
ь
Р(/) = Е(1 - °Еп(/,А))2, (2)
1=1
где ПЕп(/, А¡) - дифракционная эффективность многослойной решетки в п-м порядке в спектре отражения или пропускания на некоторой длине волны А1.
Минимизация функционала осуществляется градиентными методами, и градиент вычисляется при помощи сопряженной задачи с поставленными специальным образом граничными условиями.
2. Вычисление функциональной производной
Перейдем к непосредственно обсуждению метода вычисления градиента и постановки прямой и сопряженной задач. В общем случае идея заключается в следующем
[3]. Функция и и параметр f связаны операторным уравнением
Д(и,/ ) = 0, (3)
которое предполагается разрешимым относительно И при заданном f и является обыкновенным дифференциальным уравнением с граничными условиями.
Формулу для вариации 8^ функционала Е(f) = Ф(И,f) можно представить через вариации 8f и 8И
^ = = ФиШ + Ф^/. (4)
Затем слагаемое Фи8И заменяется на функционал от 8f при помощи уравнения в вариациях
Яи 8И + Яг 8f = 0 (5)
и тождества Лагранжа
(Ф,Яи8И) = (ЯИ Ф,8И). (6)
Далее, функция Ф подбирается как решение уравнения ЯИФ = -Фи, а из уравнения в вариациях (5) выражается Яи8И = —Яf8f, что приводит к следующей цепочке равенств
Фи8И = -(ЯИФ, 8И) = -(Ф, Яи8И) = (Ф, Яf8f) = (Я}Ф, 8f). (7)
Таким образом, из выражения для вариации функционала (4) получаем общую формулу для функциональной производной (или градиента)
дЕ
^ = Ф/ + Д}Ф. (8)
Рассмотрим применение общей схемы непосредственно к нашей задаче и начнем с получения конкретного вида операторного уравнения.
2.1. Постановка прямой задачи дифракции плоской волны на дифракционной решетке
Для определенности рассмотрим волны с ТЕ-поляризацией. В этом случае уравнение Гельмгольца для у-компоненты электрического поля
§ + § + = 0 (9)
с ограничивающими расчетную область по оси г парциальными условиями излучения в интегральной форме
а
+ Уп{х)йх = 2гк% 5п,0, (10)
ъх,п У
\ иг
о
а
^-гк^Е
дг ЕУ
о
п
г=г„
Уп(х)^х = 0, (11)
и условиями непрерывности тангенциальных компонент поля Еу (условия сопряжения) на границах 7т между слоями с номерами т и т +1 решетки
[ЕУ ]l7m =0,
дЕу
dz
0, m = 0,1, 2,... ,M, (12)
Ym
при помощи применения неполного метода Галеркина [5] сводится к матрично-векторному обыкновенному дифференциальному уравнению второго порядка с граничными условиями третьего рода следующим образом. Приближенное решение EN (x, z) уравнения Гельмгольца ищется в виде конечного разложения по первым
N функциям базиса Yn(x) = ^J~^elkx'nX, соответствующим дифракционному порядку с номером n = 0, ±1, ±2,...
EN (x,z) = Y(x)U(z), (13)
где Y(x) = ÍY\(x), Y2(x),... , YN(x)) - вектор-строка первых N функций базиса,
U(z) = (Ui(z), U2(z),..., UN(z))T - вектор-столбец коэффициентов разложения, зависящих от переменной z.
Таким образом прямая задача для вектора коэффициентов U(z) имеет вид
U''(z)+ £ (z,f )U(z) = 0,
U'(0) + гГ(0) U(0) = 2гГ(0) Ao, (14)
U'(zs) - ir(M+1)U(zs) = 0,
где матрица
Г d
£(z,f ) = k2 / Y*(x)e(x, z)Y(x)dx - M2. (15)
0
Здесь M - диагональная матрица, на главной диагонали которой расположены постоянные распространения kx,n. Матрицы Г(0'м+1) представляют собой диагональные матрицы, на главных диагоналях которых расположены постоянные распространения ^ПП) дифракционных порядков во внешней среде и в подложке соответственно и A0 = (1, 0, 0,..., 0) - амплитудный вектор падающий на решетку волны со стороны внешней среды.
Решение полученной задачи (14) в m-м слое формально можно выписать аналитически
Um(z) Qme ( ) Am(zm-1) + Qme ( )Bm(zm-(16)
где Qm - матрица собственных векторов, а rm - диагональная матрица собственных значений матрицы £(z, f) в m-м слое, Am - амплитудные векторы облучающих структуру волн, Bm - амплитудные векторы волн, рассеянных структурой.
Численно задача (14) решается методом матриц переноса, применение которого подробно описано в работе [5]. В силу постоянства коэффициентов уравнения внутри слоя толщины hm для всей многослойной структуры, состоящей из M слоев с толщинами h1 ,h2,..., hM, матрица переноса P получается перемножением матриц переноса для каждого слоя
м Г /0 I\ 1
P = PhM PhM-1 ... Phl = nPhm , где Phm =eXp| i 0J , (17)
I - единичная матрица размерности N х N, и связывает векторы И(г) и V(z) -¿и'(г) на границах структуры г = 0 и г =
и(г8)
V(zs)
р Ги(0)\
рЫ0))
(18)
При этом граничные условия имеют следующий вид:
V(0) + Г(0) и(0) = 2Г(0) Ао, V(zs) - Г(м+1)И(г3) = 0.
(19)
Таким образом, уравнение (18) с граничными условиями (19) составляют систему линейных алгебраических уравнений для векторов И (г) и V(z) на границах многослойной структуры. Отметим, что эта система справедлива как для спектра отражения, так и для спектра пропускания, и ее решение может быть получено аналитически.
Вычислим аналитически элементы матрицы £ (г, f). Поскольку в однородных слоях величина диэлектрической проницаемости е(х, г) постоянна, матрица £ (г, f) имеет диагональный вид
£(z,f ) = коМх,г)1 - М2, (20)
где I - единичная матрица. В слое с решеткой распределение диэлектрической проницаемости е(х, г) является кусочно-постоянной функцией, которая зависит от фактора заполнения f (который является параметром управления) и имеет следующий вид:
е(х, г) =
е1, х €
х €
к-1 к-1
^2 ¿р^ ^2 ¿р + ¿ки р=0 р=0 к-1 к
¿р + ¿к fk, Е ¿р .р=0 р=0
(21)
причем к = 1, К, р = 0,к л «нулевой> подпериод равен (¿о = 0. Подставляя е(х, г) в
интеграл и проводя интегрирование по периоду решетки, равному d
(£ (z,f)),
к
2пг
к-1
к Со £ Лр
Ее
А2 (т -П2) к=1
к
р=0
[(ех -1)
' П1П2
Со/к Лк
- е1+е
Со Л к
-М
П1П2 '
Е4[(£1-1)/, + 1]-(М2)
. 2п
П1П2 '
¿к, получим
к=1
П1 = п2, (22) П1 = П2,
где обозначено Со = г—(п\ — пг).
d
Далее задачу (14) будем рассматривать в качестве операторного уравнения.
2.2. Вариация функционала
Перейдем к рассмотрению введенного ранее целевого функционала (2)
^(f ) = Ф(И^) = £ (1 - DEn(f, Аг))
(23)
1=1
1
1
где ДЕп(/, Л г) - дифракционная эффективность многослойной решетки в п-ом порядке в спектре отражения или пропускания на некоторой длине волны Лг. Вводя векторы Во - амплитудный вектор отраженной решеткой волны и Ам+1 - амплитудный вектор волны, прошедшей в подложку, запишем явное выражение для дифракционной эффективности в спектре отражения
| B0,n '2
bz,0
и в спектре пропускания
DEr,n(f,\l) = ^\B0,n\2, (24)
к(ш)
£>Я4>п(/, Аг) = \Ам+1,п\ ■ (25)
Для вычисления вариации функционала амплитудный коэффициент отраженной волны Во п в п-ом порядке удобно записать как
Во,п = аВо, (26)
где а = (0 • • • 1 • • • 0) - вектор-строка, который состоит из нулевых элементов и единицы, находящейся на п-ой позиции. Тогда квадрат модуля равен
|Во,п|2 = В*>пВо,п = Вд Нп Во, (27)
где
/° ••• 0\
Нп = ат а = I . 1 . I (28)
••• 0)
- матрица, у которой единица стоит на пересечении п-ой строки и п-го столбца. Аналогично для амплитудного коэффициента
прошедшей волны Ам+1,п в п-ом порядке
|АМ +1,п1 = АМ+1 Нп Ам+1. (29)
Получим вариацию функционала в виде (4)
ь
= -2 ^ [1 - ОЕп(/, Лг)] [6(ДЕп)и + 6(ЯЕп)^ . (30)
1=1
Поскольку ДЕп(/, Лг) явно от / не зависит, то 6(ДЕп)f = 0. Учитывая условия связи амплитудных векторов падающей, прошедшей и отраженной волн
Ао + Во = и(0), Ам+1 = ), (31)
вычислим вариацию 6(ДЕп)и. Рассмотрим сначала случай спектра отражения
8(ОЕг>п)и = 8 (4уВо^В0
V kz,0 / и
= (и(°) - А°) + (и(°) - А°Г ^¿и(0)). (32)
, о
к№
Для удобства дальнейших вычислений обозначим Х(0) = —щ- Нга (и(0) — Ао). Тогда
к(,0
выражение для Й(БЕп)и примет вид
6(ВЕг,п)и = [Х*(0) ЙИ(0)] * + Х*(0) ЙИ(0). (33)
Аналогично для случая спектра пропускания
( к{Ш)
,га)и — й ( —Щ- Ам+1 ^га А-М+1
V к(, 0 / и
к(Ш)
+ и*(г,) Кгайи(^))
к( ) к(,0
[X*(г3) 6И(г3)] * + Х*(га) ЙИ(га), (34)
к(™)
где обозначено Z(zs) = —
к( ) к(,0
И, таким образом, окончательно получаем формулы для вариация функционала в случае спектра отражения
ь
ЙЕ = -2 ^ [1 - БЕ^, Аг)] [(Х*(0) ЙИ(0))* + Х*(0) ЙИ(0)] (35)
1=1
и в случае спектра пропускания
ь
ЙЕ = -2 [1 - БЕ^, А*)] [(Х*(г3) ЙИ(га)) * + X* {га) 6И(г3)]. (36)
1=1
2.3. Задача в вариациях
Операторное уравнение для (14) в виде (3)
Я(И^): И''(г)+ £ (z,f )И(г) = 0,
И'(0) + ¿Г(0) И(0) = 2гГ(0) А0, (37)
И'(га) - ¿Г(м+1)И(га) = 0.
Для того, чтобы получить уравнение в вариациях в форме (5), проварьируем Я(И, f)
ЙЯ(И, f) = ЯиЙИ + Я/6f =
= ЙИ''(г) + £ (г, f )ЙИ(г) + £/ (г, f )И(г^ = 0, (38)
где Яи ЙИ = ЙИ'' (г) + £ (г, f )ЙИ(г) и Я/ 8f = £/ (г, f )И(z)Йf. Выполним варьирование граничных условий (14). На левой границе
(И + ЙИ)'(0) + ¿Г(0) (И + ЙИ)(0) = 2гГ(0)А0, (39)
И'(0) + ЙИ'(0) + ¿Г(0) И(0) + ¿Г(0) ЙИ(0) = 2гГ(0)А0. ( )
С учетом (14) приходим к следующему граничному условию в нуле для 6И(0)
6И'(0) + гГ(о) 6И(0) = 0. (40)
Аналогично на правой границе
6И (г,) — гГ(М+1) 6И( г,) = 0. (41)
Итак, в ходе проделанных выше преобразований получена задача в вариациях
6И"(г) + £ (г, / )6И(г) + £ (г, / )И(г)6/ = 0, 6И'(0) + гГ(о) 6И(0) = 0, (42)
гГ(М+1) 6И( г,) = 0.
2.4. Постановка сопряженной задачи
Чтобы получить сопряженное уравнение для функции Ф(г), воспользуемся тождеством Лагранжа (6) и выполним интегрирование по частям с учетом граничных условий (42) для 6И(г)
¿а
(Ф,Ди 6И) = У Ф*(г) (6И''(г)+ £ (г,/ )6И(г)) ¿г =
о
¿а
= У (Ф*''(г) + Ф*(г)£(г, /))6И(г)^ — [Ф*'(г,) — Ф*(г,)гГ(М+1)]6И(г,) + о
+ [Ф*'(0) + Ф*(0)гГ(о)]6И(0) = (ДИФ, 6И), (43)
*
где под операцией понимается комплексное сопряжение и транспонирование.
Здесь и далее функция Ф*(г) представляет собой вектор-строку комплексно сопряженных к функции U(z) коэффициентов разложения. Будем выбирать ее так, чтобы
Ф*"(г) + Ф* (z)£ (z,f) = 0. (44)
Это уравнение является общим для двух спектров. Однако граничные условия для спектра отражения будут несколько отличаться от граничных условий для спектра пропускания. Рассмотрим каждый спектр отдельно.
2.4.1. Спектр отражения
Используя явное выражение для вариации функционала в случае спектра отражения (35), на левой границе при z = 0 положим
[Ф*'(0) + Ф* (0)гГ(0)]Ш(0) = -Х*(0)Ш(0), (45)
а на правой при z = zs
-[Ф*'Ы - Ф*^)гГ(м+1)]¿U(zs) = 0. (46)
Тогда для функции Ф*^) имеем следующую постановку сопряженной задачи для случая спектра отражения
Ф*''ф + Ф* (z)£ (z,f) = 0,
Ф*'(0) + Ф* (0)ir(0) = -X* (0), (47)
Ф*'^) - Ф*(zs)гГ(м+1) = 0.
Однако полученную сопряженную задачу удобнее переписать для случая, когда функция Ф(г) - вектор-столбец и (Г)* - комплексно сопряженная матрица. Для этого применим к операцию комплексного сопряжения и транспонирования
Ф''(г)+ £*(г, /)Ф(г) = 0,
Ф'(0) - г (Г(0))* Ф(0) = -Х(0), (48)
Ф'(га) + г (Г(м+1))* Ф(га) = 0.
Таким образом, задача (48) является сопряженной к задаче (14) в случае спектра отражения.
2.4.2. Спектр пропускания
В случае спектра пропускания на левой границе при г = 0 положим
[Ф*' (0) + Ф*(0)гГ(0) ]ЙИ(0) = 0, (49)
а на правой границе, при г = га, с учетом формулы для вариации функционала в случае спектра пропускания (36)
-[Ф*'(га) - Ф*(га)гГ(м+1)]ЙИ(га) = -Х{г8). (50)
Тогда для вектора-строки Ф*(г) имеем
Ф*''(г) + Ф*(г)£ (г,/) = 0,
Ф*' (0) + Ф*(0)гГ(0) = 0, (51)
-Ф*' (га) + Ф*(га )гГ(м+1) = -X* {г3).
и, также применяя операцию комплексного сопряжения и транспонирования, перейдем к сопряженной задачи для вектора-столбца Ф(г)
Ф''(г)+ £*(г, /)Ф(г) = 0,
Ф'(0) - г (Г(0))* Ф(0) = 0, (52)
-Ф'(га) - г (Г(м+1))* Ф(га) = ^).
Задача (48) представляет собой сопряженную к (14) задачу в случае спектра отражения.
Решение сопряженных задач (48) и (52) в т-м слое будет аналогично решению (16) прямой задачи
Фт(г) = дтегГт((-(т-1)Ст(гт-1) + дтв-гГт((-(т-1)Бт(гт-1). (53)
Здесь Qm
и (Гт)2 матрицы собственных векторов и собственных значений матрицы £*(г, f) в т-м слое соответственно, причем матрица (Г^)2 - диагональная, а векторы С0 и Б0 соответствуют векторам А0 и В0.
2.5. Формула для градиента
С учетом предположений, сделанных в предыдущем разделе, интеграл (43) в случае спектра отражения сводится к следующему виду
(Ф, ЯийИ) = -Х*(0)ЙИ(0) = (ЯИФ, ЙИ). (54)
Выражая из (38) Д/¿/ и применяя тождество Лагранжа, преобразуем Х*(0)$И(0)
—X*(0)£И(0) = (ДИФ, ¿И) = (Ф, Ди^И) = —(Ф, Д/¿/) = — (Д*Ф, ¿/). (55) То же проделаем и для [Х*(0)£И(0)]*
— Х*(0)5И(0)]* = (ДиФ,^И)* = (Ф, Ди^И)* = —(Ф, Д/¿/)* = —(Д*Ф,/)*. (56)
Подставляя эти преобразованные величины в (35), находим выражение для вариации функционала с случае спектра отражения в операторном виде
ь
¿Р = — 2 £ [1 — Я£р,„(/, А,)] [(Д*Ф, ¿/)* + (Д*Ф, ¿/Я , (57)
1=1
из которого легко получается формула для градиента в случае спектра отражения
ЯР Ь
_ = -2 £ [1 - ЯДР>П(/, А,)] [(Д}Ф)* + Д}ф] = (58)
= —4 ЕЬ=1 [1 — Я£р>га(/, А,)] Яе (Д*Ф).
С учетом преобразований
) = Д* Ф, 2* (г,) = (Д* Ф)* (59)
в случае спектра пропускания получаются аналогичные выражения для вариации функционала
dF . „ „
¿Р = —^ [1 — ¿>Е,п(/, А,)] [(Д*Ф, ¿/)* + (Д*Ф, ¿/)] (60)
1=1
и для градиента
= -4 ^ [1 - Лг)] 11е(Д}Ф). (61)
^ 1=1
Перейдем от операторного представления к интегральному. Воспользовавшись тождеством Лагранжа (6) и выражением (38), определим интеграл Д* Ф
¿а
Д* Ф = у Ф*(г)£/ (г,/ )И(г)^, (62)
о
где элементы матрицы £/(г,/) вычисляются аналитически (см. формулу (73)), и получим выражение для функциональной производной в интегральном виде
= -ЯЗДЛ)] Ке . (63)
Формула (63) справедлива как для спектра отражения, так и для спектра пропускания, где ДЕга(/, А,) - дифракционная эффективность в п-м порядке спектра отражения или пропускания. Здесь И(г) и Ф*(г) - решения прямой и сопряженной задач, получаемые при помощи метода матриц переноса.
ПР ^, Р = ехР<{ ¿(^ А кт У (65)
3. Численная реализация решения сопряженной задачи и вычисления градиента
Алгоритм решения прямой задачи дифракции плоской волны на дифракционной решетки методом матриц переноса подробно изложен в работе [5]. В данном разделе кратко рассмотрим применение метода матриц переноса для решения сопряженной задачи.
Интеграл в (63) вычисляется по формуле численного интегрирования прямоугольниками, для чего слой с решеткой разбивается на М^ подслоев, причем матрица £(г, /) одинаковая для каждого подслоя.
3.1. Метод матриц переноса для решения сопряженной задачи
Матрицы переноса для каждого слоя и для всей структуры в случае сопряженной задачи имеет вид, аналогичный соответствующим матрицам в случае прямой задачи (17). То есть,
ФЫ\ Р( Ф(°Л (64)
где матрица Р, связывающая вектора Ф(г) и W(z) = —гФ'(г) на границах многослойной структуры, представляется через матрицы переноса слоя
М , ,
О I
Р^т, Р^т =ехР<; И
т=1 ^ ^
Для случая спектра отражения уравнение (64) с граничными условиями
W(0) — (Г(0))* Ф(О) = ¿Х(0),
W(zs) + (Г(М+1)Ф(*в) = 0, 1 ;
представляют собой систему сопряженных к системе (18)-(19) линейных алгебраических уравнений для векторов Ф(г) и W(z), решение которой так же можно получить аналитически.
Аналогично в случае спектра пропускания уравнение (64) и граничные условия
W(0) — (Г(0))* Ф(0) = О,
W(zs) + (Г(М+1))* Ф(*в) = — ), 1 ;
являются системой сопряженных к системе (18)—(19) линейных алгебраических уравнений для векторов Ф(г) и W(z), решение которой можно получить аналитически.
3.2. Вычисление градиента
Будем вычислять интеграл в (63) по численной формуле прямоугольников. Поскольку в однородных слоях матрица (г, /) = 0, будем вести расчеты на отрезке [О, где Л — толщина 1-го слоя с решеткой. Разобьем отрезок [0, на М^ равных Н\
частей АЛ,1 = ——. Решение прямой и сопряженной задач в //-й части разбиения решетки можно представить через матрицы РМ(АЛ1) и РМ(АЛ1) соответственно
'иДМ^ = ) /У(°)\ ( ф,(АЛ1) \ = Р,(АЛ1) ^(°)\ (68) У^АЛ^ = Р (АЛ1 ^и(0^ , ^(АЛ)) = Р (АЛ1^ Ф(°^ , (68)
где ¡л = 1, МН1. Поскольку значения функций при г = 0 известны из граничных условий, то при г = ДЛ^ для вектора и^ДЛ^) получим следующее выражение
) = Р?1(ДМи(0) + Р?2(Д^1 )У(0). (69)
Аналогично для вектора Ф* (ДЛг)
Ф*(ДМ = Р?1 (Д^1)Ф(0) - Р^2(Д^1 )^Х0). (70)
Здесь матрицы РМ(Д^1) и РМ(Д^1) в слоях решетки вычисляются аналогично Р-матрице всей структуры (17).
Тогда, учитывая полученные выражения, интеграл в (63) можно вычислить по формуле численного интегрирования прямоугольниками
[ Ф*(г)£(г, /)и(г)^ = ДЛ^ ФДД^ )£(Ль /^(Д^). (71)
Таким образом, численная формула для градиента имеет вид
Mhl
яр L
= _4^[l -DEn(fAi)]^e f 1=1
M=1
(72)
Матрица ^(Л1, /) для к-го подпериода решетки является трехмерной, ее элементы одинаковы для каждого слоя ДЛ1 разбиения решетки и вычисляется аналитически один раз сразу для всей решетки по следующей формуле:
© (z,f))
k
UIU2
(k-l
Л
2п\ 2 4, ч Co £dP+dk fk
t) J(£l~ )e ' Щ^П2>
2n\2 dk t .
T -J(£l~lh П1 = П2>
(73)
,2n
где при fc = 0 «нулевой> подпериод равен do = 0, и обозначено Со = г—(ni — П2).
d
4. Результаты вычислительного эксперимента
Предложенный в данной работе градиентный метод решения задач оптимального управления был реализован в виде программного кода на базе пакета программ Матьав. С помощью реализованного алгоритма решены задачи синтеза отражающей и пропускающей дифракционных решеток, параметры которых были взяты из работ [6] и [7] соответственно. Рассматривались задачи синтеза решеток с прямоугольным профилем, определяемым одним и десятью параметрами управления. Для проверки и отладки работы алгоритма градиент, вычисленный по формуле (72), сравнивался с градиентом, вычисленным конечно-разностным методом. Во всех расчетах удерживалось N = 9 дифракционных порядков, при интегрировании слой с решеткой разбивался на М^ = 15 равных частей.
4.1. Отражающая решетка
Синтезируется многослойная дифракционная решетка для использования в качестве зеркала внешнего резонатора полупроводникового лазера, обеспечивающая стабилизацию длины волны излучения и подстройку лазера в диапазоне длин волн от 485 до 515 нм (зелено-голубая область видимого спектра) [6]. Такая решетка должна обладать максимальной дифракционной эффективностью в конфигурации Литтроу в минус первом порядке в заданном диапазоне длин волн в спектре отражения.
Структура синтезируемой решетки представлена на рис. 2. В качестве однородной многослойной структуры выбрано 20-слойное четвертьволновое зеркало, с показателями преломления материалов слоев пя = 2,375 - оксид ниобия(У) ^Ь2О5) и п^ = 1,46 - кварц (БЮ2), настроенное на длину волны А = 500 нм, то есть на центральною длину волны заданного диапазона. Решетка выполняется из кварца, а подложка - из оксида ниобия с показателями преломления пёТ = 1,46 и пц\ = 2,375 соответственно.
->■ X
/а
Щ
X
а
н
п
ёг
пн
пь
пш
Рис. 2. Синтезируемая многослойная отражающая дифракционная решетка, определяемая одним параметром управления f. Параметры профиля решетки: d = 384,8 нм - период решетки, к = 438,6 нм - глубина профиля, к0 = 21,1 нм -толщина дополнительного однородного слоя из материала решетки
На первом этапе для того, чтобы убедиться в правильности полученных выражений, вычисленный при помощи сопряженной задачи градиент (72) сравнивается с градиентом, полученным конечно-разностным методом. Как видно на графике зависимости градиента от управляющего параметра, показном на рис. 3б, вычисленные разными методами значения градиента совпадают в пределах допустимой погрешности.
В ходе вычислений были получены два оптимальных значения фактора заполнения: при начальном приближении Д = 0,3 оптимальный фактор заполнения равен f = 0,1974, а при Д = 0,6 - f = 0,5074. Объяснение этому можно найти, обратившись к графику зависимости функционала от фактора заполнения (см. рис. 3а), на котором видны две точки минимума. В зависимости от выбора начального приближения по направлению антиградиента можно <прийти> в тот или иной минимум. И, как показано на рис. 4, эти две точки действительно являются оптимальным решением, поскольку значения дифракционной эффективности для них близки к единице, что и требовалось в постановке задачи.
Функционал
30 20 10 0
и. -10
тз
™ -20 -30 -40 -50
Градиент функционала
- конечно-разностный
- сопряженная задача
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 f
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 f
а)
б)
Рис. 3. а) Зависимость функционала от фактора заполнения; б) Зависимость градиента, вычисленного по формуле (72) и конечно-разностным методом, от фактора заполнения
Дифракционная эффективность
V 0.94
450 460 470 480 490 500 510 520 530 540 550 Л , нм
Рис. 4. Зависимость дифракционной эффективности в -1-м порядке от длины волны
4.2. Пропускающая решетка
Перейдем к рассмотрению задачи оптимизации бинарной пропускающей решетки с большим числом параметров, представленной в работе [7]. Данная решетка проектируется из фторида бария (БаЕ2) с показателем преломления пёТ = 1,396, чтобы заменить треугольную решетку, обеспечивающую пропускание энергии нормально падающего излучения с длиной волны А = 10 600 нм в первый дифракционный порядок с максимальной эффективностью. Для этого каждый период прямоугольной решетки, с d = 40 955,3 нм и к = 24 646,9 нм, разбивается на 10 равных частей, каждой из которых соответствует свой фактор заполнения (см. рис. 5). Такие решетки могут быть использованы в качестве оптического ограничителя для предотвращения повреждения работающих в ИК-диапазоне элементов наведения управляемых ракет ослепляющим лазерным оружием.
В ходе решения задачи оптимизации методом были найдены оптимальные значения факторов заполнения, представленные в таблице, при которых дифракционная эффективность в первом порядке составляет 90,34%.
6
5
и- 3
0
f1 f2 ' ' ' f10
Рис. 5. Схематическое изображение замены треугольной пропускающей решетки на бинарную пропускающую решетку
Таблица
Оптимальные значения факторов заполнения
h /2 /3 /4 /5 /б /7 /8 /9 /ю
0 0 0,0136 0,1537 0,2470 0,3186 0,3988 0,4492 0,5457 1
Кроме того, время работы минимизирующего алгоритма с градиентом, вычисляемым при помощи решения сопряженной задачи, оказывается в 5 меньше, чем время работы аналогичного алгоритма с вычислением градиента конечно-разностным методом. Это связано с тем, что на каждой итерации в первом случае целевой функционал вычисляется 2 раза (решение прямой и сопряженной задач), в то время как во втором случае - 11 раз (решение прямой задачи 1 раз для вычисления функционала и 10 раз для вычисления градиента).
Заключение
В работе представлена методика вычисления градиента целевого функционала для применения к решению задачи синтеза одномерных многослойных отражающих и пропускающих дифракционных решеток. Кроме того, приведена строгая математическая постановка задачи оптимизации и полная формулировка численного алгоритма ее решения. Результаты численных экспериментов показывают надежность градиентного метода и его устойчивость относительно увеличения количества параметров управления. Это позволяет в дальнейшем перейти к рассмотрению общей постановки задачи синтеза, при которой профиль решетки определяется произвольной функцией без каких-либо специальных априорных предположений (кроме практической реализуемости), являющейся результатом непосредственного решения задачи.
Работа выполнена при финансовой поддержке Российского научного фонда (проект № 23-11-00056).
Литература
1. Palmer, C. Diffraction Grating Handbook j C. Palmer, E. Loewen. - New York: Newport Corporation, 2005.
2. Loewen, E.G. Diffraction Gratings and Application j E.G. Loewen, E. Popov. - New York: Marcel Dekker, 1997.
3. Федоренко, Р.П. Приближенное решение задач оптимального управления j Р.П. Федо-ренко. - М.: Наука, 1978.
4. Артемьева, М.В. Решение задач синтеза дифракционных решеток для практических приложений j М.В. Артемьева, А.Н. Боголюбов, А.А. Петухов j j Физические основы приборостроения. - 2020. - Т. 9, № 3. - С. 4-13.
5. Артемьева, М.В. Метод матриц переноса для решения задачи дифракции плоской волны с ТЕ-поляризацией на одномерной бинарной дифракционной решетке /М.В. Артемьева, А.Н. Боголюбов, А.А. Петухов // Физические основы приборостроения. - 2022. - Т. 11, № 2. - С. 40-48.
6. Петухов, А.А. Математическое моделирование многослойных дифракционных решеток / А.А. Петухов, А.Н. Боголюбов, М.К. Трубецков // Физические основы приборостроения. - 2014. - Т. 3, № 4. - С. 20-27.
7. Nan, Lu. Design of Transmission Blazed Binary Gratings for Optical Limiting with the Form-Birefringence Theory / Lu Nan, Dengfeng Kuang, Guoguang Mu // Applied Optics. - 2008. -V. 47, № 21. - P. 3743-3750.
Максим Олегович Корпусов, доктор физико-математических наук, профессор, кафедра математики, Московский государственный университет; математический институт, Российский университет дружбы народов (г. Москва, Российская Федерация), korpusov@gmail.com.
Маргарита Витальевна Артемьева, аспирант, кафедра математики, Московский государственный университет (г. Москва, Российская Федерация), artemeva.mv14@physics.msu.ru.
Поступила в редакцию 17 апреля 2024 г-
MSC 35Q60, 65L12, 74H35 DOI: 10.14529/mmp240205
ON CALCULATING OF FUNCTIONAL DERIVATIVE FOR AN OPTIMAL CONTROL PROBLEM
M.O. Korpusov1'2, M.V. Artemeva2
1People's Friendship University of Russia, Moscow, Russian Federation 2Moscow State University, Moscow, Russian Federation E-mail: korpusov@gmail.com, artemeva.mv14@physics.msu.ru
The synthesis problem of a multilayer diffraction grating is formulated as an optimal control problem and consists in minimizing the cost functional depending on the geometric parameters of the grating profile. The gradient method is the most reliable and stable method for solving this problem. The paper deals with a method for calculating the gradient of the cost functional, which is done by solving a cojugate problem with special boundary conditions. Additionally, the paper discusses the numerical implementation of this solution and the calculation of the gradient. The results from a computational experiment are also presented.
Keywords: functional derivative; gradient; conjugate problem; optimal control problem; synthesis problem; diffraction gratings.
References
1. Palmer C., Loewen E. Diffraction Grating Handbook. New York, Newport Corporation, 2005.
2. Loewen E.G., Popov E. Diffraction Gratings and Application. New York, Marcel Dekker, 1997.
3. Fedorenko R.P. Approximate Solution of Optimal Control Problems. Moscow, Nauka, 1978.
4. Artemeva M.V., Bogolyubov A.N., Petukhov A.A. Solving Diffraction Grating Synthesis Problems for Practical Applications. Physical Bases of Instrumentation, 2020, vol. 9, no. 3, pp. 4-13. DOI: 10.25210/jfop-2003
5. Artemeva M.V., Bogolyubov A.N., Petukhov A.A. Transfer Matrix Method for Solving the Problem of Diffraction of a Plane Wave with Te-Polarization on a One-Dimensional Binary Diffraction Grating. Physical Bases of Instrumentation, 2022, vol. 11, no. 2, pp. 40-48. DOI: 10.25210/jfop-2202-040048
6. Petukhov A.A., Bogolyubov A.N., Trubetskov M.K. Mathematical Modeling of Multilayer Diffraction Gratings. Physical Bases of Instrumentation, 2014, vol. 3, no. 4, pp. 20-27. DOI: 10.25210/jfop-1404
7. Lu Nan, Kuang Dengfeng, Mu Guoguang. Design of Transmission Blazed Binary Gratings for Optical Limiting with the Form-Birefringence Theory. Applied Optics, 2008, vol. 47, no. 21, pp. 3743-3750. DOI: 10.1364/ao.47.003743
Received April 17, 2024
