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

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

Д.В. Сошников12, А.А. Мингазов1, Л.Л. Досколович1,2 1 Институт систем обработки изображений, НИЦ «Курчатовский институт», 443001, Россия, г. Самара, ул. Молодогвардейская, д. 151; 2 Самарский национальный исследовательский университет имени академика С.П. Королёва, 443086, Россия, г. Самара, Московское шоссе, д. 34

Аннотация

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

Ключевые слова: оптический дизайн, геометрическая оптика, задача Монжа -Канторовича, линейная вариационная задача, линейное программирование.

Цитирование: Сошников, Д.В. Многомасштабный подход в задаче линейного программирования для расчета преломляющего оптического элемента, формирующего заданное распределение освещенности в дальней зоне / Д.В. Сошников, А.А. Мингазов, Л.Л. Досколович // Компьютерная оптика. - 2024. - Т. 48, № 3. - С. 356-362. - DOI: 10.18287/2412-6179-CO-1379.

Citation: Soshnikov DV, Mingazov AA, Doskolovich LL. A multiscale approach to a linear programming problem for calculating a refractive optical element forming a given far-field illumination distribution. Computer Optics 2024; 48(3): 356-362. DOI: 10.18287/2412-6179-CO-1379.

Введение

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

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

Наиболее эффективные методы расчета [1 - 11] так или иначе основаны на данной взаимосвязи. Далеко не все задачи неизображающей оптики могут быть сформулированы в виде задачи перемещения масс. В частности, такой переформулировки нет для задачи формирования заданного распределения освещенности в ближней зоне. Тем не менее некоторые идеи, связанные с задачей перемещения масс, работают и в этих случаях [12, 13]. В [9, 11] описаны известные типы оптических элементов, для которых задача расчета сводится к задаче перемещения масс с некоторой функцией стоимости.

Мы рассматриваем задачу расчета преломляющего оптического элемента, формирующего заданное распределение освещенности в дальней зоне при параллельном освещающем пучке лучей. Эта задача может быть сформулирована как задача перемещения масс Монжа-Канторовича. Но в большинстве публикаций, например, в [9, 11], рассматривается только формулировка задачи перемещения масс для задачи формирования заданного распределения интенсивности, которая эквивалентна данной. В данной статье мы явно приводим задачу перемещения масс для распределения освещенности в дальней зоне. Также мы описываем связь функционала Лагранжа и вариаци-

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

Ранее делались попытки [14, 15] использования линейной вариационной задачи для расчета отражающего элемента, формирующего заданное распределение интенсивности. Их нельзя назвать удачными, поскольку использовалось решение задачи линейного программирования, являющейся прямой аппроксимацией линейной вариационной задачи. Соответственно, из-за высокой вычислительной сложности, которую мы отмечали ранее, авторам [14, 15] удалось получить решения только для изображений, аппроксимируемых сеткой из 64 точек.

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

Постановка задачи расчета оптического элемента

Здесь мы подробно опишем математическую формулировку задачи расчета преломляющего оптического элемента. По сравнению с реальной задачей вводятся два существенных упрощения. Во-первых, на данном этапе мы игнорируем факт наличия френе-левских потерь, считая, что весь световой поток достигает экрана. При этом в дальнейшем френелевские потери можно учесть, выбрав решение, минимизирующее отклонение лучей (минимум задачи перемещения масс) и, соответственно, близкое к оптимуму по френелевским потерям. Во-вторых, рассматривается приближение «дальней зоны», то есть мы пренебрегаем размерами оптического элемента. В этом случае точка экрана, в которую попадает луч, определяется только его направлением. Это достаточно распространенное приближение в светотехнических расчетах. Например, эквивалентная задача формирования заданного распределения интенсивности рассматривается в публикациях [1, 2, 9, 11, 17]. В отличие от этих статей мы рассмотрим задачу перемещения масс непосредственно для формирования распределения освещенности в удаленной плоскости, в результате функция стоимости будет иметь несколько иной вид, чем в указанных статьях.

Пусть Е3 - трехмерное пространство с координатами (иь и2, ¿). Рассмотрим параллельный световой пучок, распространяющийся в направлении е = (0,0,1) и имеющий в плоскости г = 0 распределение освещённости /(и), иеО. На удаленной плоскости П, задаваемой уравнением г =/, расположен экран.

Предположим также, что задана преломляющая поверхность у=Я(и), при этом показатель преломления при ъ< Я(х) равен щ, а при ъ< Я(и) равен н2. Она задает лучевое отображение у: О ^ П следующим образом. Луч, выходящий из точки плоскости г = 0 с координатами (и1, и2) и направлением е, преломляясь на поверхности г=Я(и), приобретает направление х(и) = (х;(и), т2(и), т3(и)). Поскольку расстояние / до экрана очень велико по сравнению с размером оптического элемента, можем пренебречь размерами оптического элемента и считать, что точка пересечения преломленного луча с плоскостью г =/ определяется только направлением х(и). Соответственно, отображение у в координатах задается формулой

у(и) =

/ ^(и) / Т1(и)

Тз(и) ' Тз(и)

В результате на плоскости г =/ формируется распределение освещённости Ь(х) в некоторой области Б, определяемое условием сохранения светового потока

|Ь(x)dx = | /{и)йи,

(1)

у-1(ю)

для любого борелевского ю с Б.

В случае дифференцируемой поверхности рефрактора соотношение (1) можно записать в виде уравнения Монжа - Ампера:

/ (и) = Ь(у(и)) • 3., (и),

(2)

где 3у(и) - якобиан отображения у, вычисленный в точке и.

Теперь, если поверхность Я неизвестна, но заданы функции /(и) и Ь(х) с носителями О и Б соответственно, возникает задача расчета преломляющей поверхности из условия формирования заданного распределение освещенности Ь на удаленной плоскости г =/ То есть требуется рассчитать поверхность Я, которая индуцирует биективное отображение у (описанным выше способом), для которого выполнено условие (2). Итоговая геометрия рассматриваемой нами задачи приведена на рис. 1.

Функция стоимости и задача перемещения масс

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

Утверждение 1. Плоскость, фокусирующая весь световой поток, распространяющийся в направлении е = (0,0,1), в точку (х1, х2,/) плоскости г =/ задается уравнением вида

г = К(и, х) - £, (3)

где £ - константа, а

со

K(u, x) =

n2 xu + n2 x2u2 nls] Xi2 + x| + f2 - n2 f

(4)

Рис. 1. Геометрия задачи

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

( , А

Я =

f

VxfTXFT

f2 ' Vx2 + x22 +

f2 ' Vxi2 + x2 + f2

Направление входящего луча е = (0,0,1), а единичную нормаль обозначим N = (Л^;,^2,^3). Поскольку векторы N, е, Я лежат в одной плоскости, можем записать

N = а • е + Ь • Я.

Векторно умножив это соотношение на N, получим

0 = а • е х N + Ь • Я х N.

Поскольку векторы единичные, то е х N = 8Шф,пс, Я х N = 5Шфге/ . По закону Снелла

п1 • е х N = п2 • Я х N.

Значит,

а=_ я х N

1

e х N

П2

То есть N || n1 • e -n2 • R, где

(

n1 • e - n2 • R =

n2 x.

n2 x2

tJx2 + x2 + f2 '

niyjx2 + x2 + f2 -n2f ^

yjx? + x2 + f2' Jxf+xfTJ2

Обозначим через Nz = (Nz1,Nz2, Nz3) вектор, параллельный n1 • e - n2 • R, для которого Nz3 =1. То есть

Nz =

n2 x,

n2 x2

n1y] x,2 + x| + f2 - n2 f n1 ■Jxxf+xf+f2 - n2 f

Уравнение плоскости с нормалью N2 имеет вид (3). □ Поверхность 2 = Я(и) является огибающей семейства поверхностей вида (3), если для некоторой дифференцируемой функции £(х) выполнено

R(u) = K(u, x) - S (x),

д

dxt

\K(u, x) - S(x)] = 0, i = 1,2.

(5)

Как и в статьях [1 - 5, 7, 8], мы будем рассматривать частный случай огибающих, когда во втором уравнении (5) имеет место не произвольная критическая точка, а точка минимума. В этом случае условие (5) можно записать

R(u) = min(^(u, x) - S (x)).

xeD

(6)

Лучевое соответствие при этом может быть задано формулой

у(и) = arg min(X(w, x) - S (x)).

xeD

Определение 2. Отображение T: G ^ D называется интегрируемым, если для некоторой функции Se C(D) имеет место представление

T (u) = argmin(^(u, x) - S (x)).

xeD

Соотношение (7) полностью определяет отображение T при известной функции S. Будем обозначать такое отображение 7®.

Определение 3. Будем называть отображение T: G ^ D удовлетворяющим условию сохранения светового потока, если для любого борелевского подмножества ю с D выполнено

J I(u)du = jz(x)dx.

Т-!(т) ш

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

Лемма 4. Условие сохранения потока для обратимого отображения Т: О ^ Б равносильно любому из условий:

1) для любого борелевского ш с Б выполнено:

J I(u)du = jz(x)dx;

T _1(m) rn

2) для любой h e C(D) выполнено: Jh(x) L( x)dx = Jh(T (u ))I (u )du;

D G

3) I(u) = L(T(u))-Jt (u), где Jt (u) - якобиан отображения T, вычисленный в точке u. Доказательство. Доказано в [4, lemma 4.11]. □ Задачу расчета оптического элемента можно

сформулировать следующим образом. Для заданных распределений I(u) и L(x) требуется найти лу-

чевое соответствие Т: О ^ Б, которое удовлетворяет условию сохранения светового потока (8) и условию интегрируемости (7), а также восстановить прямой и двойственный параметры оптического элемента.

Данная задача может быть сформулирована как задача перемещения масс Монжа-Канторовича.

Теорема 5. Рассмотрим функционал

Т(Т) = ]Ди, Т(и))/(и^и :

(9)

заданный на пространстве отображений Т: О ^ Б, удовлетворяющих условию сохранения светового потока (8). Минимум функционала Т является интегрируемым отображением.

Доказательство. Скажем несколько слов о доказательстве, хотя аналогичные рассуждения приводились в [6, 8, 11]. Рассмотрим функционал Лагранжа для рассматриваемой вариационной задачи, взяв в качестве ограничений эквивалентное условие 3 из леммы 4. Возьмем произвольную цеС(О), и пусть

Н(Т, ц) = \[К(и, Т (и)) / (и) + Х(Т (и)) >

О

х( Ь(Т (и)) • 3Т (и) - / (и ))^и.

(10)

Преобразуем, пользуясь эквивалентным условием сохранения светового потока 2 из леммы 4

Н(Т,ц) = |[Ди,Т(и)) -ЦТ(и))]/(и^и -

О

+|Х( х) Ь( x)dx.

(11)

Условный минимум функционала (9) при ограничениях (8) соответствует критической точке функционала (10). Вариация функционала Н легко вычисляется. Условие 5лН = 0, как нетрудно понять из формулы (10), равносильно условию сохранения светового потока. Условие равенства нулю вариации 5тН = 0 равносильно условию интегрируемости (5).

Функционал Лагранжа (11) можно модифицировать таким образом, что он будет зависеть только от функции X. Заметим, что в критической точке функционала Н отображение Т и X связаны соотношением (7), то есть Т = Т х. Соответственно, при рассмотрении функционала мы можем заранее ограничиться парами (Т х, X). Соответственно, мы можем рассматривать функционал

Н(Х) = |[Ди, Тх (и)) - Х(Тх (и))]/(и^и -

О

+|Х( х) Ь( х^х,

(12)

или можно записать иначе

Н(Х) = |гпт[Ди, х) — Х( х)]/(и^и +|Х( х) Ь( х)ах. (13)

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

Функционал Н связан с линейным функционалом Я, который мы будем рассматривать далее. Определим Я на пространстве пар функций С(О)*С(Б) следующей формулой

Я(ц, X) = |ц(и)/(и)du + |Х(х)Ь(x)dx.

О Б

Заметим, что Н(Х) = Я (цх, X) для функции

^ (и) = шт[Ди, х) - X(х)].

(14)

(15)

Утверждение 6. Рассмотрим пару непрерывных функций (ц, X)eC(О)xC(Б). Тогда

Я (ц, X) < Я (^, X).

Доказательство. Аналогичные рассуждение приведены в доказательстве теоремы 1 в [18]. □

Кроме того, заметим, что для функции X и функции ц, определенной соотношением (15), выполнено неравенство

X(x) + ц(и) = X(х) + тш[^(и, х) - X(х)] <

хеБ

< X(х) + К(и, х) - X(х) = К(и, х).

(16)

Теперь мы можем рассмотреть линейный функционал О(ц, X) на пространстве пар функций с ограничениями (16):

П = {(ц, X) е С (О) х С (Б) | ц(и) + X( х) <

< К(и,х) Уи е О,Ух е Б}.

Утверждение 7. Максимумы функционалов Я и Н совпадают. Если максимум Н достигается в X0, то максимум Я достигается в точке (^ , X0).

Доказательство. Из утверждения 6 следует: для каждого значения Н(к) существует равное значение О(ц\ X). Наоборот, для пары (ц, X)еQ невозможно сопоставить значение О(ц, X). Тем не менее, эта пара не является экстремальной (если не имеет вид (ц\ X)). □

Многомасштабный подход

Мы будем рассматривать численное решение максимизации функционала

Я (ц, X) = |ц(и) / (и )du + ]Х( х) Ь( x)dx

^ тах

на пространстве

П = {(ц, X) е С (О) х С (Б) | X( х) + ц(и) < < К(и, х) У и е О, Ух е Б}.

(17)

(18)

Разобьем области О и Б на непересекающиеся области О г и Б, соответственно,

N М

О = ЦО,, Б = ЦБ,.

И пусть в каждой области О, выбраны точки и,, а в каждой Б, выбрана точка х,. Функция ц(и) заменяется набором значений ц, в точках и,, а X(x) заменяется набором значений X, в точках х,. Аналогично распределения освещенности описываются наборами чисел Ь и Ь,. Обозначим также К, = К(и,, х,). В результате мы получим задачу линейного программирования

Я(M, X) = XM1 ' + XXjLj ^ max,

i=i j=i

Mi +X^ <K,i = 17N, j = W.

(19)

(20)

Дискретная задача определяется набором данных {(G,u e Gj, {(Dj,xj e Dj)}=. Такая дискретная схема (19), (20) используется в статьях [14, 15]. В силу огромного количество ограничений (20), даже в случае изображений небольшого размера, использовать такую дискретную схему удается лишь в очень ограниченных случаях. Скажем, в упомянутых статьях рассматривалась задача формирования изображений, имеющих размеры 5 на 5 точек. Для преодоления этой трудности можно использовать многомасштабный подход, различные варианты которого возникают и в других вариантах задачи перемещения масс, например, [16, 17]. Многомасштабный подход предполагает последовательное решение задачи на сгущающихся сетках, причем решение на более грубой сетке выступает либо начальным приближением для решения на грубой сетке, как [16, 17], либо упрощает задачу, как в рассматриваемом случае. В нашей ситуации мы используем решение на грубой сетке, чтобы избавиться от большого числа несущественных ограничений.

Предположим, что нам известны наборы чисел (Mi), (Xj), являющиеся решением задачи (21), (22) для набора данных дискретизации {(G°ld,u°ld)}-=°и , {(Dfd, xold)jf .

Заметим, что условие того, что лучевое отображение у переводит точку u в точку x, формулируется в терминах функций м, X как

M(u) + X( x) = K(u, x).

Для нас будут особенно важны пары точек, которые близки к лучевой паре (u, y(u)) для данного приближения. Более точно, пусть выбрано некоторое число е > 0. Будем называть пару точек (ui, xj) квазилучевой, если выполнено

| Mi +Xj - Kj |< е.

Допустим, у нас есть более тонкая сетка new, определяемая набором данных {(G"ew, u"ew) ,

{(D']ew, xjw )}М="™ , которая является подразбиением сетки old. Последнее означает, что для каждого i имеется набор индексов а1, ..., anj, для которых выполнено

ni

Goid = ,

1=1

а также для каждого j имеется набор индексов Рь ..., Р". , для которых выполнено

Df =

new

Р1 •

Соответственно, для точки ui грубой сетки old точки u^,, ..., uan. тонкой сетки new будем называть соседями и обозначать их множество Neigh(ui), аналогично для xj грубой сетки old соседями будем называть точки xp1, ..., xpm, тонкой сетки new и обозначать их множество Neigh(xj).

Если пара (ui, xj) сетки old является квазилучевой, пары точек тонкой сетки Neigh(u1)^-Neigh(xj) будем называть потенциально лучевыми. Теперь мы можем записать упрощенную задачу на тонкой сетке.

Я (М, X) = Y^Milr + XX jLj

^ max,

j=1

Mi +Xj < Kj

new j Kij

(21) (22)

для всех потенциально лучевых пар (и, х,).

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

Теперь, имея последовательность измельчающихся сеток, представленных наборами данных {(О,да,е О,®)}^ , {(Б,к),хк е Б™, , цепочка решений зависит только от набора е(к). Выбор последовательности е(к) имеет большое значение.

Расчетные примеры

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

В обозначениях, введенных ранее, функция освещенности источника на плоскости г = 0 имеет вид

I (u ) =

1,| u1 |< 1,| u21< 1; 0, иначе.

m

Расстояние до плоскости f= 1000, показатели преломления щ = 1,49 (полиметилметакрилат), n2 = 1. Мы начинаем с аппроксимации требуемого распределения освещенности сеткой из 16 на 16 точек. На следующей итерации мы измельчаем обе сетки в 4 раза (вдвое равномерно уплотняем сетку по горизонтали и вертикали). Нам потребовалось три итерации, то есть в итоге изображение аппроксимируется дискретной сеткой в 64 на 64 точки. Использовались следующие значения параметров для итераций: е(1) = 10 -5, е(2) = 10 -7. Время расчета составило несколько часов, притом что без применения многомасштабного подхода не удавалось решить задачу даже на сетке 32 на 32 точки. Все получаемые в процессе задачи линейного программирования решаются симплекс-методом. Дискретные значения, получаемые в результате решения задачи линейного программирования, аппроксимировались сплайном в программе Rhinoceros 3D. Полученный таким образом оптический элемент и результаты моделирования в программе TracePro приведены на рис. 2.

О

О 100 во

х, мм

Рис. 2. Преломляющая поверхность, формирующая заданное распределение освещенности в виде кольца (а);

результат моделирования работы преломляющей поверхности в программе для светотехнических расчетов ТгасеРго (б)

Нормированное (на !2-норму требуемого распределения) среднеквадратичное отклонение полученного распределения от постоянного распределения на кольце - 8,4 %, световая эффективность оптического элемента - 97,2 %.

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

I (u) =

|1,

10, иначе.

Дискретные сетки, число итераций и пороги е те же, что и в предыдущем примере. Рассчитанный оптический элемент и результаты моделирования в программе ТгасеРго приведены на рис. 3.

Нормированное среднеквадратичное отклонение полученного распределения от постоянного распределения на кольце - 9,3 %, световая эффективность оптического элемента - 98,3 %.

+

а) (6) „мм ±ии ""

Рис. 3. Преломляющая поверхность, формирующая заданное распределение освещенности в виде креста (а);

результат моделирования работы преломляющей поверхности в программе для светотехнических расчетов TracePro (б)

Заключение

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

Благодарности

Работа выполнена за счет средств государственного задания в сфере научной деятельности (проект FSSS-2024-0014) в части разработки и реализации многомасштабного метода и государственного задания НИЦ «Курчатовский институт» в части установления взаимосвязи линейного функционала с функционалом стоимости.

References

[1] Glimm T, Oliker VI. Optical design of single reflector systems and the Monge-Kantorovich mass transfer problem. J Math Sci 2003; 117: 4096-4108. DOI: 10.1023/A:1024856201493.

[2] Wang XJ. On design of a reflector antenna II. Calc Var Partial Differ Equ 2004; 20: 329-341. DOI: 10.1007/s00526-003-0239-4.

[3] Glimm T, Oliker VI. Optical design of two-reflector systems, the Monge-Kantorovich mass transfer problem and Fermat's principle. Indiana Univ Math J 2004; 53(5): 1255-1277. DOI: 10.1512/iumj.2004.53.2455.

[4] Rubinstein J, Wolansky G. Intensity control with a freeform lens. J Opt Soc Am A 2007; 24(2): 463-469. DOI: 10.13 64/JOSAA.24.000463.

[5] Doskolovich LL, Mingazov AA, Bykov DA, Andreev ES, Bezus EA. Variational approach to calculation of light field eikonal function for illuminating a prescribed region. Opt Express 2017; 25(22): 26378-26392. DOI: 10.1364ЮЕ.25.026378.

[6] Doskolovich LL, Mingazov AA, Bykov DA, Andreev ES. Variational approach to eikonal function computation. Computer Optics 2018; 42(4): 557-567. DOI: 10.18287/2412-6179-2018-42-4-557-567.

[7] Doskolovich LL, Bykov DA, Andreev ES, Bezus EA, Oliker V. Designing double freeform surfaces for collimat-ed beam shaping with optimal mass transportation and lin-

ear assignment problems. Opt Express 2018; 26(19): 24602-24613. DOI: 10.1364/OE.26.024602.

[8] Bykov DA, Doskolovich LL, Mingazov AA, Andreev ES, Kazanskiy NL. Linear assignment problem in the design of freeform refractive optical elements generating prescribed irradiance distributions. Opt Express 2018; 26(21): 2781227825. DOI: 10.1364/OE.26.027812.

[9] Yadav NK. Monge-Ampere problems with non-quadratic cost function: application to freeform optics. Eindhoven: Technische Universiteit Eindhoven; 2018.

[10] Gutierrez CE, Huang Q. The near field refractor. Annales de l'Institut Henri Poincaré C, Analyse non linéaire 2014; 31(4): 655-684. DOI: 10.1016/j.anihpc.2013.07.001.

[11] Mingazov AA, Doskolovich LL, Bykov DA, Byzov EV. Support quadric method in non-imaging optics problems that can be reformulated as a mass transfer problem. Computer Optics 2022; 46(3): 353-365. DOI: 10.18287/2412-6179-CO-1055.

[12] Graf T, Oliker V. An optimal transport approach to near-field reflector problem in optical design. Inv Problems 2012; 28(2): 025001. DOI: 10.1088/02665611/28/2/025001.

[13] Schwartzburg Y, Testuz R, Tagliasacchi A, Pauly M. High-contrast computational caustic design. ACM Trans on Graph 2014; 33(4). DOI: 10.1145/2601097.2601200.

[14] Canavesi ^ Cassarly WJ, Rolland JP. Direct calculation algorithm for two-dimensional reflector design. Opt Lett 2012; 37(18): 3852-3854. DOI: 10.1364/0L.37.003852.

[15] Canavesi Q Cassarly WJ, Rolland JP. Observations on the linear programming formulation of the single reflector design problem. Opt Express 2012; 20(4): 4050-4055. DOI: 10.1364/OE.20.004050.

[16] Merigot Q. A multiscale approach to optimal transport. Comp Graph Forum 2011; 30(5): 1583-1592. DOI: 10.1111/j.1467-8659.2011.02032.x.

[17] Bykov DA, Doskolovich LL, Bezus EA. Multiscale approach and linear assignment problem in designing mirrors generating far-field irradiance distributions. Opt Lett 2020; 45(13): 3549-3552. DOI: 10.1364/OL.393895.

[18] Doskolovich LL, Mingazov AA, Bykov DA, Andreev ES. Variational approach to eikonal function computation. Computer Optics 2018; 42(4): 557-567. DOI: 10.18287/2412-6179-2018-42-4-557-567.

Сведения об авторах

Сошников Даниил Владимирович, в 2022 году с отличием окончил бакалавриат Самарского национального исследовательского университета имени академика С.П. Королёва по специальности «Прикладная математика и информатика». E-mail: soshnikov.d.v@mail.ru

Мингазов Альберт Айдарович, в 2010 году окончил Самарский государственный университет (СамГУ, ныне - Самарский национальный исследовательский университет имени академика С.П. Королёва) по специальности «Математика». В 2012 году окончил магистратуру НИУ ВШЭ по направлению «Математика», в 2015 окончил аспирантуру Санкт-Петербургского отделения математического института имени В. А. Стеклова. Кандидат физико-математических наук (2015), научный сотрудник лаборатории дифракционной оптики Института систем обработки изображений, НИЦ «Курчатовский институт». E-mail: mingazov88@gmail. com

Досколович Леонид Леонидович, в 1989 году с отличием окончил Куйбышевский авиационный институт (КуАИ, ныне - Самарский национальный исследовательский университет имени академика С.П. Королёва) по специальности «Прикладная математика». Доктор физико-математических наук (2001 год), профессор, главный научный сотрудник лаборатории дифракционной оптики Института систем обработки изображений, НИЦ «Курчатовский институт», профессор кафедры технической кибернетики Самарского университета и ведущий научный сотрудник научно-исследовательской лаборатории прорывных технологий дистанционного зондирования Земли Самарского университета. Специалист в области дифракционной оптики, лазерных информационных технологий, нанофотоники. E-mail: leonid@ipsiras.ru

ГРНТИ: 29.31.29

Поступила в редакцию 26 июня 2023 г. Окончательный вариант - 16 сентября 2023 г.

A multiscale approach to a linear programming problem for calculating a refractive optical element forming a given far-field illumination distribution

D.V. SoshnikovA.A. Mingazov1, L.L. Doskolovich1,2 1 Image Processing Systems Institute, NRC "Kurchatov Institute";

443001, Samara, Russia, Molodogvardeyskaya 151;

2 Samara National Research University, 443086, Samara, Russia, Moskovskoye Shosse 34

Abstract

We consider an approach to the calculation of a refractive optical element generating a prescribed irradiance distribution in the far field for a plane incident beam based on a certain variational problem. We consider an explicit formulation of this problem in the form of the Monge-Kantorovich mass transfer problem. We also demonstrate a connection between the mass transfer problem and the dual linear variational problem. A numerical solution of the linear variational problem is also considered. The direct solution of this type of problem presents a huge computational complexity. To overcome this difficulty we use the so-called multiscale approach based on constructing a chain of approximations that are solutions on refining grids.

Keywords: optical design, geometric optics, Monge-Kantorovich problem, linear variational problem, linear programming.

Citation: Soshnikov DV, Mingazov AA, Doskolovich LL. A multiscale approach to a linear programming problem for calculating a refractive optical element forming a given far-field illumination distribution. Computer Optics 2024; 48(3): 356-362. DOI: 10.18287/2412-6179-CO-1379.

Acknowledgements: This work was performed within the State assignment for scientific research to Samara University (project FSSS-2024-0014) in part of the development and implementation of the multiscale approach and the State assignment of the NRC "Kurchatov Institute" in part of establishing the relationship between the linear functional and the cost functional.

Authors' information

Daniil Vadimovich Soshnikov, graduated with honors (2022) from Samara National Research University with a major in Applied Mathematics and Computer Science. E-mail: soshnikov. d.v@mail.ru

Albert Aidarovich Mingazov graduated from Samara State University (presently, Samara National Reseach University) in 2010, majoring in Mathematics. In 2012 he got the master degree in Mathematics at National Research University Higher School of Economics, in 2015 he graduated from the post-graduate course of the St. Peterburg department of Steklov Mathematical Institute of Russian Academy of Sciences. Candidate of Phisics and Mathematics (2015). Currently he is a researcher at Diffractive Optics laboratory of the Image Processing Systems Institute, NRC "Kurchatov Institute". E-mail: mingazov88@gmail.com

Leonid Leonidovich Doskolovich graduated with honours (1989) from S.P. Korolyov Kuibyshev Aviation Institute (presently, Samara National Reseach University), majoring in Applied Mathematics. He received his Doctor in Physics & Maths (2001) degree from Samara State Aerospace University. Main reseacher of the Image Processing Systems Institute, NRC "Kurchatov Institute", professor at Technical Cybernetics department of National Research University, senior researcher at the Breakthrough Technologies for Earth's Remote Sensing laboratory at SSAU. His leading research interests include diffractive optics, laser information technologies, nanophotonics. E-mail: leonid@ipsiras.ru

Received June 26, 2023. The final version - September 16, 2023.