CC BY
20
Оптимизационный метод расчета оптических элементов при протяженном источнике излучения
Е.В. Бызов1, С.В. Кравченко1, М.А. Моисеев1, Л.Л. Досколович 12 1 ИСОИ РАН - филиал ФНИЦ «Кристаллография и фотоника» РАН, 443001, Россия, г. Самара, ул. Молодогвардейская, д. 151, 2 Самарский национальный исследовательский университет имени академика С.П. Королёва, 443086, Россия, г. Самара, Московское шоссе, д. 34
Аннотация
Рассмотрен метод расчёта оптического элемента с двумя поверхностями свободной формы, формирующего заданное распределение освещенности для протяженного источника излучения. Метод основан на представлении поверхностей оптического элемента бикубическими сплайнами и последующей оптимизации их параметров с помощью квазиньютоновского метода, реализованного в пакете Matlab. Для расчета целевой функции предложен вариант метода трассировки лучей. С использованием предложенного метода рассчитан оптический элемент, имеющий рекордные характеристики: отношение высоты элемента к размеру источника излучения - 1,6; световая эффективность - 89,1 %; равномерность формируемого распределения (отношение минимальной и средней освещенности) в заданной квадратной области - 0,92.
Ключевые слова: поверхность свободной формы, трассировка лучей, оптимизация, оптика для светодиодов, распределение освещенности, неизображающая оптика, оптический дизайн.
Цитирование: Бызов, Е.В. Оптимизационный метод расчета оптических элементов при протяженном источнике излучения / Е.В. Бызов, С.В. Кравченко, М.А. Моисеев, Л.Л. Досколович // Компьютерная оптика. - 2020. - Т. 44, № 5. - С. 712-720. - DOI: 10.18287/2412-6179-CO-762.
Citation: Byzov EV, Kravchenko SV, Moiseev MA, Doskolovich LL. Optimization method for designing optical elements with an extended light source. Computer Optics 2020; 44(5): 712-720. DOI: 10.18287/2412-6179-CO-762.
Введение
В настоящее время в системах освещения широко используются светоизлучающие диоды (СИД). Для формирования заданных световых распределений в светодиодных системах освещения используется так называемая вторичная оптика (преломляющие оптические элементы с поверхностями «свободной формы», устанавливаемые непосредственно над излучающим элементом светодиода). В общем случае задача расчета преломляющего оптического элемента, формирующего заданное световое распределение, относится к классу обратных задач неизображающей оптики. Ряд эффективных методов для решения данной задачи был предложен только в приближении точечного источника излучения [1-10]. В настоящее время большой интерес представляет дизайн компактной вторичной оптики для мощных СИД, имеющих излучающий элемент относительно большого размера. В этом случае соотношение размеров оптического элемента к размеру светодиода может быть совсем небольшим (около двух или даже меньше). При этом решение, полученное в рамках приближения точечного источника, как правило, сильно отличается от оптимального решения.
Обратные задачи расчета оптических элементов при размерном источнике являются существенно более сложными. Для решения таких задач было пред-
ложено несколько подходов, включающих так называемые методы обратной связи (англ. feedback methods) [11, 12], методы «прямой оптимизации» [13 - 16] и методы «краевого луча» [17 - 23], основанные на использовании принципа краевого луча (англ. edge ray principle) для контроля формируемого светового распределения. Методы обратной связи используют в качестве начального приближения оптический элемент, рассчитанный для точечного источника излучения. Далее поверхность оптического элемента итерационно корректируется с учетом протяженного характера источника излучения [11, 12]. В методах прямой оптимизации поверхность оптического элемента задается в некотором параметрическом виде, например, в виде полинома [13] или сплайна [14, 16]. Параметры, задающие поверхность, рассматриваются как аргументы функции невязки, описывающей отличие заданного и формируемого распределений. Далее расчет оптической поверхности сводится к минимизации функции невязки с использованием какого-либо оптимизационного алгоритма. Отметим, что минимизация функции невязки требует многократного решения задачи расчета распределения освещенности, формируемого элементом при заданном наборе параметров поверхности. Данный расчет, как правило, осуществляется с использованием метода трассировки световых лучей.
Методы, использующие принцип краевого луча, «пришли» из методов расчета концентраторов солнечной энергии [17, 18]. В последние годы были предложены эффективные варианты методов краевого луча для расчета осесимметричных оптических элементов с двумя «рабочими» поверхностями, формирующих заданные радиально симметричные распределения освещенности [19-23]. В этих методах профили осесимметричного оптического элемента строятся последовательно (сегмент за сегментом) с использованием краевых лучей так, чтобы обеспечить некоторое требуемое соотношение между заданным световым потоком в текущем фрагменте формируемого светового распределения и светящейся областью поверхности оптического элемента, посылающей излучение в этот фрагмент. В недавней работе [24] в журнале "Optica" был рассмотрен метод «подстроенных» волновых фронтов (МПВФ, англ. wave-front tailoring method) для расчета оптических элементов, формирующих заданные световые распределения в дальней зоне, которые не обладают радиальной симметрией. МПВФ является развитием так называемого метода одновременного расчёта нескольких поверхностей (англ. simultaneous multiple surface design method) [1, 25, 26] и основан на преобразовании «входящих» волновых фронтов, испускаемых из краевых точек источника излучения, в «исходящие» волновые фронты, обеспечивающие формирование заданного светового распределения. Представленный в [24] пример расчета очень компактного оптического элемента c двумя поверхностями, формирующего почти идеальное световое распределение в квадратной области, демонстрирует высокие рабочие характеристики МПВФ. Поэтому работу [24] можно рассматривать как наиболее значимый результат в области расчета оптических элементов с размерным источником. В то же время, по мнению авторов настоящей статьи, МПВФ является сложным методом с точки зрения его понимания, программной реализации и применения. Кроме того, МПВФ имеет ряд существенных ограничений. Например, непонятно, как применить метод при неквадратной форме освещаемой области, в ближней зоне, при неламбертовском источнике излучения и т.п.
В настоящей работе мы предлагаем метод прямой оптимизации для решения обратных задач расчета преломляющих оптических элементов с двумя поверхностями. Методы прямой оптимизации обычно критикуются за высокую вычислительную сложность, связанную с необходимостью многократного решения прямой задачи (задачи расчёта, формируемого элементом распределения освещенности) [1, 19, 20]. Рассмотренный в работе метод трассировки световых лучей решает прямую задачу расчета освещенности в общем асимметричном случае за доли секунды на обычном персональном компьютере. Это, в совокупности c использованием квазиньютоновского
метода оптимизации, реализованного в пакете МайаЪ, позволило авторам рассчитать оптический элемент, превосходящий по своим характеристикам оптический элемент из работы [24]. Кроме того, предложенный метод, по сравнению с МПВФ, является более универсальным и не имеет ограничений, связанных с формой освещаемой области и выполнением условия дальней зоны.
1. Постановка задачи
Рассмотрим задачу расчёта оптического элемента. Пусть в области плоскости г = 0 находится протяженный источник, излучающий в полусферу г > 0. Центр источника (области 5) совпадает с началом координат. Излучение источника описывается функцией яркости Цф, у, х), хе5, где фе[0,2 л) и уе[0,л / 2) -полярный и азимутальный углы сферической системы координат, а х = (х, у) - декартовы координаты в плоскости источника. Отсчет угла у производится от оси 07. Излучение от источника падает на оптический элемент (показатель преломления материала п > 1) с двумя преломляющими поверхностями (рис. 1). Поверхности элемента задаются в сферических координатах функциями длины радиус-векторов Гт(ф, у), т = 1,2. Показатели преломления среды над и под элементом равны по = 1. Задача состоит в расчете поверхностей элемента гт(ф, у), т = 1,2 из условия формирования в плоскости г = / заданного распределения освещенности £/(и), иеД где и = (и,у) -декартовы координаты в данной плоскости.
z E(u,v) и=и(<р,у;хо)
2=/ f и
по
1 п Ыф.У)
1
source О Хо х
Рис. 1. Геометрия задачи
2. Решение прямой задачи
Как отмечено выше, оптимизационные методы решения обратной задачи расчета оптического элемента, формирующего заданное распределение освещенности, включают многократное решение прямой задачи, состоящей в расчете распределения освещенности -E(u), формируемого оптическим элементом при некоторых заданных поверхностях гт(ф, у), m = 1,2. Поэтому в данном пункте мы рассмотрим используемый метод решения прямой задачи. Для расчета £(u) рассмотрим сначала вычисление распреде-
ления освещенности Е(и; х0), создаваемого элементом ёБ протяженного источника с центром в точке х0е& Направления лучей, исходящих из точки х0, будем описывать сферическими координатами (ф, у). Будем считать, что исходящие лучи после преломления на поверхностях оптического элемента попадают в точки выходной плоскости г = /. В этом случае можно говорить, что оптический элемент задаёт некоторое «лучевое отображение» и = и(ф, у, х0), связывающее координаты (ф, у) исходящих лучей и координаты их прихода в выходную плоскость (рис. 1). Координата источника х0 в данном отображении рассматривается как параметр. Отображение и(ф, у, хо) определяет формируемое в плоскости г = / распределение освещенности Е(и; хо). Данное распределение получается из закона сохранения светового потока:
E (u )d 2u = T (ф, y; x0 )L (ф, y, x0 )d QdS±
(1)
где T (ф, y, x0) - коэффициент пропускания, описывающий френелевские потери при преломлении луча на поверхностях элемента, d Q = sin ydфdy - элемент телесного угла, dS± = cos ydS. Перепишем соотношение (1) в интегральном виде. Для этого воспользуемся фильтрующим свойством дельта-функции:
E (u; xo ) = flE(t;xo)s(t-u)d2t .
(2)
Выполним в последнем интеграле замену t = u^, y, x0, c) и, с учетом (1), получим
E (u; x0) = dSj T (ф, y; x0 )L (ф, y, x0)
®
x8(u - u (ф, y; x0 ))sin y cos yd фd y,
(3)
где © = {(ф, у)|фе[0,2я),уе[0,я /2)}. Для численных расчётов дельта-функцию в (3) нужно аппроксимировать функцией Гаусса [16]:
5П (u) = —exp (-u2/ ст2).
(4)
Такая аппроксимация позволяет получить усредненное распределение освещенности в выходной плоскости, причем радиус усреднения определяется параметром ст функции Гаусса.
Чтобы вычислить полное распределение освещенности, формируемое всем протяженным источником, необходимо проинтегрировать выражение (3) по области источника:
где
E (u ) = J E (u; x )d 2x =
S
= JJ f (x,ф,y)8- (u -u(ф,y;x))dфdydxdy,
S
f (x, ф, y) = T (ф, y; x )L (ф, y, x)sin y cos y .
(5)
(6)
В данной работе расчет распределения освещенности Е(и) основан на вычислении интеграла (5) методом Монте-Карло [27]. В этом случае генерируется N значений (реализаций) ^ = (хк, ук, фк, у к), к =1,..., N случайного вектора с независимыми компонентами, равномерно распределенного в области £*©. Каждую реализацию = (хк, ук, фк, ук) можно интерпретировать как луч, исходящий из случайной точки хк = (хк,ук,)еБ и имеющий случайное направление (фк, ук) е©. Далее, согласно методу Монте-Карло, за приближенное значение интеграла (5) принимают величину
En (u) =
N
= AV (u - u (, y k; x* ))f (x*, ф*, y k),
(7)
где ДУ=У/N, V = л2||£|| - «объём» области Бх©, ||Б|| -площадь области Б. Легко видеть, что (7) соответствует сумме элементарных световых потоков /(хк, фк, ук) ДУ, излученных случайно выбранными элементами источника в случайных направлениях. При этом данные световые потоки распределяются по окрестностям точек ик = и(фк, ук, хк) с гауссовским весом (4). Известно, что математическое ожидание случайной величины EN (и) совпадает с точным значением интеграла (5), а дисперсия обратно пропорциональна квадратному корню из числа лучей N (т.е. стремится к нулю с ростом К) [27].
Для расчета освещенности по формуле (7) необходимо выполнить трассировку лучей Эта трассировка требуется для вычисления координат ик = и(фк, ук, хк) точек пересечения излученных лучей с выходной плоскостью. Процедура трассировки лучей включает поиск точек пересечения падающих лучей с поверхностями оптического элемента, вычисление преломленных лучей и коэффициента пропускания (функция Т (ф, у, х) в (1), (6)). В данной работе для расчета точек пересечения лучей с поверхностями элемента выполнялась триангуляция оптических поверхностей гт (ф, у), т = 1,2. Поскольку поверхности Гт (ф, у) заданы в параметрическом виде, то триангуляция поверхности просто задается триангуляцией области определения ©. Число и размер треугольников в области © выбирались из условия представления функций Гт (ф, у), т = 1,2, с заданной точностью. Простейший алгоритм поиска пересечения луча с триангулированной поверхностью состоит в последовательной проверке условия пересечения луча с каждым треугольником (примитивом). Вычислительная сложность такого алгоритма составляет 0(МГ), где N - количество примитивов. Для уменьшения вычислительной сложности используется упорядочивание примитивов в пространстве. Такой подход позволяет эффективно исключать из поиска примитивы, через которые падающий луч не проходит. В данной работе использовалась структура бинарного
k =1
пространственного разбиения, называемая к^дерево [28], которая позволяет использовать алгоритм бинарного поиска [29, 30] для нахождения пересекаемого лучом примитива со средней вычислительной сложностью МГ). Направление преломленного луча в точке пересечения луча с поверхностью определяется законом Снеллиуса. Для расчета коэффициента пропускания Т (ф, у, х), который равен произведению коэффициентов пропускания для двух поверхностей элемента, использовались формулы Френеля для неполяризованного света.
Таким образом, метод решения прямой задачи включает трассировку N случайных лучей ^ = (хк, Ук, фк, у к) для вычисления координат щ = и(фк, ук, хк) и последующий расчет формируемого распределения освещенности по формуле (7).
3. Решение обратной задачи
Рассмотрим решение обратной задачи расчета поверхностей гт (ф, у), т = 1,2, оптического элемента, обеспечивающих формирование заданного распределения освещенности £Хи), ие Б. В настоящей работе предлагается определить функции Гт (ф, у), т = 1,2, в виде бикубических сплайнов [31], заданных на сетках Лт = (фт,,-,ут/), I =1,...,Мт, / = 1,...,Мт. Бикубический сплайн в каждой ячейке
Ст,/ = {(ф, у)|фт,/ < ф < фт.,+1, у»у < у < у/ } представляет собой многочлен вида:
.(ф, уЬЕЕ'
,(>■, j)
ф-ф»
p= 0 q=0 \ Р
фш,1+1 -ф»
у-у»
(8)
уm,j+1 уm,j у
Коэффициенты сплайнов а^^д могут быть выражены через значения длины радиус-векторов Гт,, = Гт(ф/, у/), значения первых производных
» -
дг
= (ф,, у j), j (ф,, у j)
дг
дф
ду
и значения смешанных производных
г (фу) =
д2 г w ' m
дфду
(ф/, у j)
в узлах сеток Лт, т = 1,2 [32]. Поэтому далее наборы значений гт1,/,Г,^-, г^/,Г^у] будем называть параметрами сплайнов. Стоит отметить, что количество точек, используемых при триангуляции поверхности, обычно на порядки превышает количество узлов Лт, т = 1,2.
(ф) Г(у) Г(фуЛ т = 12 -
Обозначим c =( г г»4/ ■ г■ г»'1'1''))
Ajruvx v.» \'m,,, j?'m,,, p ' m,,, j?'m, Ь
вектора из параметров сплайнов для первой и второй поверхности. Формируемое оптическим элементом распределение освещенности в выходной плоскости
г = / будет определяться параметрами сплайнов с = (С1, С2). Поэтому далее это распределение будем обозначать Е(и; с).
При определении поверхностей элемента в виде бикубических сплайнов задачу расчета оптического элемента можно сформулировать как задачу минимизации некоторой функции е(с), зависящей от параметров сплайнов с = (С1, С2) и представляющей отличие формируемого элементом распределения освещенности Е(и; с) от требуемого распределения Е^ (и), иеБ. В настоящей работе в качестве такой функции будем использовать среднеквадратичное отклонение (СКО) формируемого распределения освещённости от требуемого. В этом случае расчет оптического элемента сводится к следующей задаче оптимизации:
*) = JDJJ(£ (u)-E(u;c))2 d2u
^ min,
(9)
где ||D|| - площадь области D. Для решения задачи (9) был использован квазиньютоновский метод Бройдена-Флетчера-Гольдфарба-Шанно, реализованный в пакете Matlab в функции fminunc. Для использования указанной функции в простейшем случае требуется только написать подпрограмму, вычисляющую функцию (9) (при этом оценки градиента и гессиана функции e(c) вычисляются по конечно-разностным формулам). Вычисление e(c) сводится к расчёту распределения освещённости E(u; c), формируемого оптическим элементом при заданном векторе параметров с. Метод расчета освещённости E(u; c) подробно рассмотрен в предыдущем параграфе. Данный метод был также реализован в программной среде Matlab.
Отметим, что для представления поверхностей оптического элемента можно было выбрать и другие функции, например, полиномы [13] или иррациональные B-сплайны [14]. Авторы данной статьи решили использовать бикубические сплайны (8) в сферической системе координат, поскольку ранее такое представление позволило им эффективно решать различные задачи расчета оптических элементов при точечном источнике излучения [16, 33, 34].
4. Пример расчета
В качестве примера был рассчитан оптический элемент (показатель преломления материала n = 1,491), формирующий в плоскости z =/= 1000 мм «сглаженное равномерное» распределение освещенности в квадратной области при протяженном лам-бертовском источнике излучения 1*1 мм. Данное распределение является постоянным в квадратной области с размером стороны d =2/-tg(30°) «1155 мм и далее плавно спадает до нуля (рис. 2в). Сечение формируемого распределения показано пунктиром на рис. 3г. Такие вид и размеры распределения были выбраны для сравнения с результатами недавней статьи в журнале "Optica" [24]. В этой статье с использова-
X
нием МПВФ рассмотрен расчет оптического элемента, формирующего аналогичное распределение освещенности при таких же размерах источника излучения. Рассчитанный элемент является очень компактным, отношение высоты элемента к размеру диагонали источника равно 1,67. Это хороший показатель с учетом высокого качества, сформированного элементом распределения [24].
Рассмотрим применение предложенного в настоящей работе оптимизационного метода для решения рассмотренной задачи. В силу симметрии формируемого распределения освещенности, рассчитываемый оптический элемент имеет следующие плоскости симметрии: у = 0 (ф = 0), х = 0 (ф = ж / 2), у = ± х (ф = ж / 4 и ф = 3ж / 4). Поэтому поверхности оптического элемента достаточно задать только при фе[0, ж /4]. При этом «полные» поверхности оптического элемента (при фе[0,2 ж) получаются из поверхностей при фе[0, ж /4] последовательным отражением относительно плоскостей симметрии ф = ж / 4, ф = ж / 2 и ф = 0. Соответственно в качестве параметров оптимизации (параметров бикубических сплайнов, задающих радиус-векторы поверхностей элемента) достаточно рассматривать только параметры сплайнов при фе[0, ж /4]. Согласно пункту 4, параметры сплайнов определяются из решения задачи оптимизации (9). В качестве начального приближения для процесса оп-
тимизации мы использовали радиально-симметрич-ный элемент, рассчитанный для точечного источника из условия формирования в плоскости г = 1000 мм постоянного распределения освещенности в круге с радиусом Я = 500 мм. Расчет такого элемента подробно описан в [34] и сводится к интегрированию системы из трех обыкновенных дифференциальных уравнений. Поверхности рассчитанного элемента были аппроксимированы бикубическими сплайнами. Полученный в результате оптический элемент показан на рис. 2а. Оптимизируемые узлы сплайнов, представляющих внутреннюю и внешнюю поверхности элемента при фе[0, ж /4], отмечены на рис. 2а синими и красными точками. Заметим, что для того, чтобы полные поверхности оптического элемента не имели изломов, в узлах сплайновых сеток при ф = 0 и ф = ж /4 производные сплайнов т^ и были положены
равными нулю и исключены из процесса оптимизации. В узлах, расположенных в вершинах поверхностей (при у = 0), в качестве оптимизируемых параметров рассматривались только значения радиус-векторов гт (значения первых и вторых смешанных производных были положены равными нулю). При указанных условиях и выбранных сплайновых сетках число оптимизируемых параметров сплайнов равно 13 (внутренняя поверхность) и 85 (внешняя поверхность).
1200
1200
600
-600
-1200
б)
-1200
600 1200
-1200 -600
Рис. 2. Начальное приближение для оптического элемента (чёрными и серыми точками отмечены оптимизируемые параметры бикубических сплайнов) (а); нормированное распределение освещённости, формируемое оптическим элементом в плоскости г = 1000 тт при ламбертовском источнике излучения с размером 1*1 мм2 (б); требуемое сглаженное равномерное распределение освещённости в квадратной области (в)
На рис. 2б показано распределение освещенности, формируемое оптическим элементом на рис. 2а при протяженном ламбертовском источнике с размером 1x1 мм2. Данное распределение сильно отличается от требуемого распределения на рис. 2в, однако, как будет показано далее, оптимизация выбранных начальных поверхностей позволит получить хорошее решение задачи.
Расчет оптического элемента (минимизация функции (9)) осуществлялся с использованием функции йштипс из стандартной библиотеки МаЙаЬ. Для вычисления функции Е(и; с) в процессе оптимизации использовался метод трассировки лучей, рассмотренный в параграфе 3. При трассировке использовалось N = 500000 лучей. Такое число лучей обеспечивает приемлемый компромисс между временем и точностью вычислений (при увеличении числа лучей в 2 раза расчетное распределение освещенности Е(и; с) меняется не более чем 2-3 %). При числе лучей
N = 500000 среднее время одного вычисления функции e(c) на персональном компьютере Intel Core i7-3770 составляет 0,67 секунды. В процессе минимизации функция e(c) была вычислена 19296 раз, при этом полное время расчета элемента составило около 3,5 часов. Рассчитанный оптический элемент показан на рис. 3. Отношение высоты рассчитанного оптического элемента к размеру диагонали источника составило 1,6 (рис. 3б), что чуть лучше, чем для элемента из работы [24], формирующего аналогичное распределение. Важно отметить, что оптический элемент из работы [24], в отличие от рассчитанного элемента (рис. 3а, б), менее приспособлен для производства и практического применения. Действительно, оптический элемент из работы [24] не является «замкнутым» и состоит из двух отдельных фрагментов поверхностей, которые захватывают лишь часть излученного источником светового потока (см. рис. 5 в [24]).
4.17 мм
2,26 мм
1200
600
-600
-1200
б)
1,0
0,5
\0
z /
в) -1200 -600 0 600 1200 г) -1200 -600 0 600 и,у,мм
Рис. 3. Рассчитанный оптический элемент для формирования сглаженного равномерного распределения освещенности в квадратной области в плоскости г = 1000 тт при ламбертовском источнике с размером 1*1 мм2 (а); сечения оптического элемента и источника излучения в плоскостях ср= 0 и ср= ж/4 (б); нормированное распределение освещенности, формируемое оптическим элементом (распределение рассчитано в программе TracePro при числе лучей N = 1000000) (в); сечение формируемого распределения (г). Пунктиром показано требуемое распределение
Распределение освещенности, формируемое оптическим элементом, приведено на рис. 3в. Отметим, что для подтверждения правильности выполненных расчетов данное распределение было рассчитано в коммерческом программном обеспечении ТгасеРго [35]. Полученное распределение хорошо совпадает с требуемым распределением освещенности, сечение которого показано пунктирной линией на рис. 3г. Доля излученной источником энергии, попавшей в требуемую область, составила около 86 %. Следуя работе
[24], для характеристики равномерности формируемого светового распределения будем использовать отношение Ет„ / Еауегаёе минимальной и средней освещенности в заданной квадратной области с размером стороны d = 1155 мм. Для распределения освещенности на рис. 3в Ет„ / Еaverage = 0,92. Для оптического элемента из работы [24], формирующего аналогичное распределение освещенности, значение Ет„ / Еaverage немного хуже и составляет 0,9. Рассчитанный оптический элемент, по сравнению с оптиче-
ским элементом из работы [24], имеет большую световую эффективность. Под эффективностью понимается доля излученного источником светового потока, которая попадает в выходную плоскость. Для оптического элемента из работы [24] световая эффективность, которая была рассчитана без учета френелев-ских потерь, составила 80 %. Для оптического элемента на рис. 3 световая эффективность, рассчитанная в программе ТгасеРго с учетом френелевских потерь, составила 89,1 %. Это значение эффективности близко к теоретически достижимому пределу. Действительно, для плоской границы раздела между воздухом и материалом с показателем преломления п коэффициент пропускания Френеля при нормальном падении света равен Т\= 4п / (п +1)2. Поскольку оптический элемент имеет две преломляющие поверхности, то результирующий коэффициент пропускания может быть оценен как Т< Т12= [4п/(п +1)2]2. В этом случае теоретически достижимое значение световой эффективности равно ^ =1-Т2. В рассматриваемом примере п = 1,491 и значение теоретически достижимой эффективности составляет 92,4 %.
Представленный пример расчета оптического элемента подтверждает высокие рабочие характеристики предложенного оптимизационного метода. Важно отметить, что предложенный метод, по сравнению с методом работы [24], является более универсальным и не имеет ограничений, связанных с формой освещаемой области и выполнением условия дальней зоны. Применение метода для расчета оптических элементов, формирующих более сложные распределения, будет предметом дальнейших исследований.
Заключение
Рассмотрен оптимизационный метод расчёта оптического элемента с двумя поверхностями свободной формы, формирующего заданное распределение освещенности для протяженного источника излучения. В рамках метода поверхности оптического элемента задаются в виде бикубических сплайнов. Параметры сплайнов рассматриваются как аргументы целевой функции, описывающей среднеквадратиче-ское отклонение между формируемым и заданным распределениями. Для расчета целевой функции предложен вариант метода трассировки лучей. Минимизация целевой функции осуществляется с помощью квазиньютоновского метода оптимизации, реализованного в пакете МайаЪ.
С помощью представленного метода был рассчитан компактный оптический элемент, формирующий сглаженное равномерное распределение освещенности в квадратной области с угловым размером 60° при протяженном ламбертовском источнике излучения. Отношение высоты рассчитанного элемента к размеру источника составляет всего 1,6. Световая эффективность оптического элемента равна 89,1 %, а
равномерность формируемого распределения, определяемая как отношение средней и минимальной освещенности, составляет 0,92. По нашему мнению, характеристики рассчитанного элемента в настоящий момент являются рекордными.
Благодарности
Работа выполнена при поддержке Министерства науки и высшего образования в рамках выполнения работ по Государственному заданию ФНИЦ «Кристаллография и фотоника» РАН в части численного моделирования работы оптических элементов и Российского научного фонда (проект № 18-19-00326) в части разработки метода и расчёта поверхностей оптических элементов.
Литература
1. Wu, R. Design of freeform illumination optics / R. Wu, Z. Feng, Z. Zheng, R. Liang, P. Benitez, J.C. Minano, F. Duerr // Laser & Photonics Reviews. - 2018. - Vol. 12, Issue 7. - 1700310.
2. Wu, R. Freeform illumination design: a nonlinear boundary problem for the elliptic Monge-Ampere equation / R. Wu, L. Xu, P. Liu, Y. Zhang, Z. Zheng, H. Li, X. Liu // Optics Letters. - 2013. - Vol. 38, Issue 2. - P. 229-231.
3. Mao, X. Design of a smooth freeform illumination system for a point light source based on polar-type optimal transport mapping / X. Mao, S. Xu, X. Hu, Y. Xie // Applied Optics. - 2017. - Vol. 56, Issue 22. - P. 6324-6331.
4. Wu, R. Formulating the design of two freeform lens surfaces for point-like light sources / R. Wu, S. Chang, Z. Zheng, L. Zhao, X. Liu // Optics Letters. - 2018. - Vol. 43, Issue 7.
- P. 1619-1622.
5. Brix, K. Designing illumination lenses and mirrors by the numerical solution of Monge-Ampere equations / K. Brix, Y. Hafizogullari, A. Platen // Journal of the Optical Society of America A. - 2015. - Vol. 32, Issue 11. - P. 2227-2236.
6. Yadav, N.K. Computation of double freeform optical surfaces using a Monge-Ampere solver: Application to beam shaping / N.K. Yadav, J.H.M. ten Thije Boonkkamp, W.L. Ijzerman // Optics Communication. - 2019. - Vol. 439. - P. 251-259.
7. Schwartzburg, Y. High-contrast computational caustic design / Y. Schwartzburg, R. Testuz, A. Tagliasacchi, M. Pauly // ACM Transactions on Graphics. - 2014. -Vol. 33, Issue 4. - 74.
8. Oliker, V. Controlling light with freeform multifocal lens designed with supporting quadric method(SQM) / V. Oliker // Optics Express. - 2017. - Vol. 25, Issue 4. - P. A58-A72.
9. Doskolovich, L.L. Designing double freeform surfaces for collimated beam shaping with optimal mass transportation and linear assignment problems / L.L. Doskolovich,
D.A. Bykov, E.S. Andreev, E.A. Bezus, V. Oliker // Optics Express. - 2018. - Vol. 26, Issue 19. - P. 24602-24613.
10. Doskolovich, L.L. Optimal mass transportation and linear assignment problems in the design of freeform refractive optical elements generating far-field irradiance distributions / L.L. Doskolovich, D.A. Bykov, A.A. Mingazov,
E.A. Bezus // Optics Express. - 2019. - Vol. 27, Issue 9. -P. 13083-13097.
11. Luo, Y. Design of compact and smooth free-form optical system with uniform illuminance for LED source / Y. Luo, Z. Feng, Y. Han, H. Li // Optics Express. - 2010. - Vol. 18.
- P. 9055-9063.
12. Li, Z. Energy feedback freeform lenses for uniform illumination of extended light source LEDs / Z. Li, S. Yu, L. Lin, Y. Tang, X. Ding, W. Yuan, B. Yu // Applied Optics. -2016. - Vol. 55. - P. 10375-10381.
13. Liu, Z. Parametric optimization method for the design of high-efficiency free-form illumination system with a LED source / Z. Liu, P. Liu, and F. Yu, // Chinese Optics Letters. - 2012. - Vol. 10. - P. 112201-112201.
14. Fournier, F. Optimization of freeform lightpipes for light-emitting-diode projectors / F. Fournier, J. Rolland // Applied Optics. - 2008. - Vol. 47. - P. 957-966.
15. Zhao, S. Integral freeform illumination lens design of LED based pico-projector / S. Zhao, K. Wang, F. Chen, Z. Qin, S. Liu // Applied Optics. - 2013. - Vol. 52. - P. 2985-2993.
16. Moiseev, M.A. Design of refractive spline surface for generating required irradiance distribution with large angular dimension / M.A. Moiseev, L.L. Doskolovich // Journal of Modern Optics. - 2010. - Vol. 57, Issue 7. - P. 536-544. -DOI: 10.1080/09500341003764069.
17. Ries, H. Tailored edge-ray reflectors for illumination / H. Ries, R. Winston // Journal of the Optical Society of America A. - 1994. - Vol. 11. - P. 1260-1264.
18. Chaves, J. Introduction to nonimaging optics / J. Chaves. -2nd ed. - CRC Press, 2016.
19. Wu, R. Direct three-dimensional design of compact and ultra-efficient freeform lenses for extended light sources / R. Wu, C.Y. Huang, X. Zhu, H.-N. Cheng, R. Liang // Optica. - 2016. - Vol. 3. - P. 840-843.
20. Wu, R. Direct design of aspherical lenses for extended non-Lambertian sources in two-dimensional geometry / R. Wu, H. Hua, P. Benitez, J.C. Minano // Optics Letters. - 2015. -Vol. 40. - P. 3037-3040.
21. Hu, S. Ultra-compact LED lens with double freeform surfaces for uniform illumination / S. Hu, K. Du, T. Mei, L. Wan, N. Zhu // Optics Express. - 2015. - Vol. 23. -P. 20350-20355.
22. Li, X. Prescribed intensity in 3D rotational geometry for extended sources by using a conversion function in 2D design / X. Li, P. Ge, H. Wang // Applied Optics. - 2017. -Vol. 56. - P. 1795-1798.
23. Li, X. An efficient design method for LED surface sources in three-dimensional rotational geometry using projected angle difference / X. Li, P. Ge, H. Wang // Lighting Research & Technology. - 2019. - Vol. 51, Issue 3. - P. 457-464.
24. Sorgato, S. Design of illumination optics with extended sources based on wavefront tailoring / S. Sorgato, J. Chaves, H. Thienpont, F. Duerr // Optica. - 2019. - Vol. 6. - P. 966971.
25. Sorgato, S. Compact illumination optic with three freeform surfaces for improved beam control / S. Sorgato, R. Mohedano, J. Chaves, M. Hernández, J. Blen,
D. Grabovickic, P. Benítez, J. Miñano, H. Thienpont,
F. Duerr // Optics Express. - 2017. - Vol. 25. - P. 2962729641.
26. Dross, O. Review of SMS design methods and real world applications / O. Dross, R. Mohedano, P. Benítez, J.C. Miñano, J. Chaves, J. Blen, M. Hernández, F. Muñoz // Proceedings of SPIE. - 2004. - Vol. 5529. - P. 35-47.
27. Robert, C.P. Monte Carlo statistical methods / C.P. Robert,
G. Casella. - 2nd ed. - Springer, 2004. - ISBN: 978-1-44191939-7.
28. Bentley, J.L. Multidimensional binary search trees used for associative searching / J.L. Bentley // Communications of the ACM. - 1975. - Vol. 18, Issue 9. - P. 509-517.
29. Cormen, T.H. Introduction to algorithms / T.H. Cormen, C.E. Leiserson, R.L. Rivest. - MIT Press, 1989.
30. Sedgewick, R. Algorithms in C++. / R. Sedgewick. - Addi-son-Wesley, 1992.
31. Шикин, Е.В. Кривые и поверхности на экране компьютера. Руководство по сплайнам для пользователей /
E.В. Шикин, Л.И. Плис. - М.: ДИАЛОГ-МИФИ, 1996.
32. Бикубическая интерполяция [Электронный ресурс]. -URL: https://en.wikipedia.org/wiki/Bicubic_interpolation (дата обращения 23.07.2020 г.)
33. Moiseev, M.A. Design of efficient LED optics with two free-form surfaces / M.A. Moiseev, S.V. Kravchenko, L.L. Doskolovich // Optics Express. - 2014. - Vol. 22. -P. A1926-A1935. - DOI: 10.1364/0E.22.0A1926.
34. Kravchenko, S.V. Design of axisymmetric double-surface refractive optical elements generating required illuminance distributions / S.V. Kravchenko, M.A. Moiseev, E.V. Byzov, L.L. Doskolovich // Optics Communications. - 2020. -Vol. 459. - 124976. - DOI: 10.1016/j.optcom.2019.124976.
35. TracePro - программное обеспечение для дизайна и моделирования оптических систем. - URL: https://www.lambdares.com/tracepro/ (дата обращения 23.07.2020 г.).
Сведения об авторах
Бызов Егор Владимирович, 1988 года рождения. В 2014 году с отличием окончил магистратуру Самарского государственного аэрокосмического университета имени академика С.П. Королёва (СГАУ, ныне Самарский национальный исследовательский университет имени академика С.П. Королева) по направлению «Прикладные математика и физика». Область научных интересов: методы расчетов формирующей неизображающей оптики для светодиодов. E-mail: egor. byzov@gmail. com.
Кравченко Сергей Васильевич, 1991 года рождения. В 2014 году с отличием окончил магистратуру Самарского государственного аэрокосмического университета имени академика С.П. Королёва (ныне Самарский национальный исследовательский университет имени академика С.П. Королева) по направлению «Прикладные математика и физика». Кандидат физико-математических наук (2018 год). Является членом международного научного общества SPIE. В списке научных трудов С.В. Кравченко 18 публикаций. E-mail: sergey.kravchenko91@gmail.com .
Моисеев Михаил Александрович, 1986 года рождения, в 2008 году с отличием окончил Самарский государственный аэрокосмический университет имени академика С.П. Королёва (ныне Самарский национальный исследовательский университет имени академика С. П. Королева) по специальности 010600 «Прикладные математика и физика». Кандидат физико-математических наук (2011 год), работает научным сотрудником в Лабо-
ратории дифракционной оптики ИСОИ РАН - филиале ФНИЦ «Кристаллография и фотоника» РАН. Является членом международного научного общества SPIE. М.А. Моисеев - специалист в области решения обратных задач геометрической оптики. В списке научных работ М.А. Моисеева 35 статей и 7 авторских свидетельств. E-mail: mikhail@smr.ru .
Досколович Леонид Леонидович, в 1989 году с отличием окончил Куйбышевский авиационный институт (КуАИ, ныне - Самарский национальный исследовательский университет имени академика С.П. Королева) по специальности «Прикладная математика». Доктор физико-математических наук (2001 год), профессор, работает заведующим лабораторией дифракционной оптики ИСОИ РАН - филиала ФНИЦ «Кристаллография и фотоника» РАН, профессором кафедры технической кибернетики Самарского университета. Специалист в области дифракционной оптики, лазерных информационных технологий, нанофотоники. E-mail: leonid@ipsiras.ru .
ГРНТИ: 29.31.29
Поступила в редакцию 1 июня 2020 г. Окончательный вариант - 27 июля 2020 г.
Optimization method for designing optical elements with an extended light source
E. V. Byzov ', S. V. Kravchenko ', M.A. Moiseev', L.L. Doskolovich 1,2 'IPSIRAS - Branch of the FSRC "Crystallography and Photonics" RAS, 443001, Samara, Russia, Molodogvardeyskaya 151, 2Samara National Research University, 443086, Samara, Russia, Moskovskoye Shosse 34
Abstract
A method for designing an optical element with two free-form surfaces generating a prescribed illuminance distribution in the case of an extended light source is considered. The method is based on the representation of the optical element surfaces by bicubic splines and on the subsequent optimization of their parameters using a quasi-Newton method implemented in the Matlab software. To calculate the merit function, a version of the ray tracing method is proposed. Using the proposed method, an optical element with record characteristics was designed: the ratio of the element height to the source size is 1.6; luminous efficiency is 89.1 %; uniformity of the generated distribution (the ratio of the minimum and average illuminance) in a given square region is 0.92.
Keywords: freeform surface, ray tracing, optimization, LED optics, illuminance distribution, nonimaging optics, optical design.
Acknowledgements: This work was financially supported by the RF Ministry of Science and Higher Education within the State assignment to FSRC «Crystallography and Photonics» RAS under agreement 007-0/^3363/26 (Numerical simulation of the optical elements), Russian Science Foundation under project No. 18-19-00326 (Development of the method and design of the optical elements).
Citation: Byzov EV, Kravchenko SV, Moiseev MA, Doskolovich LL. Optimization method for designing optical elements with an extended light source. Computer Optics 2020; 44(5): 712-720. DOI: 10.18287/2412-6179-C0-762.
References
[1] Wu R, Feng Z, Zheng Z, Liang R, Benitez P, Minano JC, Duerr F. Design of freeform illumination optics. Laser Photon Rev 2018; 12(7): 1700310.
[2] R. Wu, L. Xu, P. Liu, Y. Zhang, Z. Zheng, H. Li, and X. Liu, "Freeform illumination design: a nonlinear boundary problem for the elliptic Monge-Ampere equation," Opt. Lett. 38(2), 229-231 (2013).
[3] Mao X, Xu S, Hu X, Xie Y. Design of a smooth freeform illumination system for a point light source based on polartype optimal transport mapping. Appl Opt 2017; 56(22): 6324-6331.
[4] Wu R, Chang S, Zheng Z, Zhao L, Liu X. Formulating the design of two freeform lens surfaces for point-like light sources. Opt Lett 2018; 43(7): 1619-1622.
[5] Brix K, Hafizogullari Y, Platen A. Designing illumination lenses and mirrors by the numerical solution of Monge-Ampere equations. J Opt Soc Am A 2015; 32(11): 2227-2236.
[6] Yadav NK, ten Thije Boonkkamp JHM, Ijzerman WL. Computation of double freeform optical surfaces using a Monge-Ampere solver: Application to beam shaping. Opt Commun 2019; 439: 251-259.
[7] Schwartzburg Y, Testuz R, Tagliasacchi A, Pauly M. High-contrast computational caustic design. ACM Trans Graph 2014; 33(4): 74.
[8] Oliker V. Controlling light with freeform multifocal lens designed with supporting quadric method(SQM). Opt Express 2017; 25(4): A58-A72.
[9] Doskolovich LL, Bykov DA, Andreev ES, Bezus EA, Oliker V. Designing double freeform surfaces for collimat-ed beam shaping with optimal mass transportation and linear assignment problems. Opt Express 2018; 26(19): 24602-24613.
[10] Doskolovich LL, Bykov DA, Mingazov AA, Bezus EA. Optimal mass transportation and linear assignment problems in the design of freeform refractive optical elements generating far-field irradiance distributions / Opt Express 2019; 27(9): 13083-13097.
[11] Luo Y, Feng Z, Han Y, Li H. Design of compact and smooth free-form optical system with uniform illuminance for LED source. Opt Express 2010; 18: 9055-9063.
[12] Li Z, Yu S, Lin L, Tang Y, Ding X, Yuan W, Yu B. Energy feedback freeform lenses for uniform illumination of extended light source LEDs. Appl Opt 2016; 55: 1037510381.
[13] Liu Z, Liu P, Yu F. Parametric optimization method for the design of high-efficiency free-form illumination system with a LED source. Chin Opt Lett 2012; 10: 112201112201.
[14] Fournier F, Rolland J. Optimization of freeform lightpipes for light-emitting-diode projectors. Appl Opt 2008; 47: 957-966.
[15] Zhao S, Wang K, Chen F, Qin Z, Liu S. Integral freeform illumination lens design of LED based pico-projector. Appl Opt 2013; 52: 2985-2993.
[16] Moiseev MA, Doskolovich LL. Design of refractive spline surface for generating required irradiance distribution with large angular dimension. J Mod Opt 2010; 57(7): 536-544. DOI: 10.1080/09500341003764069.
[17] Ries H, Winston R. Tailored edge-ray reflectors for illumination. J Opt Soc Am A 1994; 11: 1260-1264.
[18] Chaves J. Introduction to nonimaging optics. 2nd ed. CRC Press; 2016.
[19] Wu R, Huang CY, Zhu X, Cheng H-N, Liang R. Direct three-dimensional design of compact and ultra-efficient freeform lenses for extended light sources. Optica 2016; 3: 840-843.
[20] Wu R, Hua H, Benítez P, Miñano JC. Direct design of as-pherical lenses for extended non-Lambertian sources in two-dimensional geometry. Opt Lett 2015; 40: 3037-3040.
[21] Hu S, Du K, Mei T, Wan L, Zhu N. Ultra-compact LED lens with double freeform surfaces for uniform illumination. Opt Express 2015; 23: 20350-20355.
[22] Li X, Ge P, Wang H. Prescribed intensity in 3D rotational geometry for extended sources by using a conversion function in 2D design, Appl Opt 2017; 56: 1795-1798.
[23] Li X, Ge P, Wang H. An efficient design method for LED surface sources in three-dimensional rotational geometry using projected angle difference Light Res Technol 2019; 51(3): 457-464.
[24] Sorgato S, Chaves J, Thienpont H, Duerr F. Design of illumination optics with extended sources based on wave-front tailoring. Optica 2019; 6: 966-971.
[25] Sorgato S, Mohedano R, Chaves J, Hernández M, Blen J, Grabovickic D, Benítez P, Miñano J, Thienpont H, Duerr F. Compact illumination optic with three freeform surfaces for improved beam control. Opt Express 2017; 25: 2962729641.
[26] Dross O, Mohedano R, Benítez P, Miñano JC, Chaves J, Blen J, Hernández M, Muñoz F. Review of SMS design methods and real world applications. Proc SPIE 2004; 5529: 35-47.
[27] Robert CP, Casella G. Monte Carlo statistical methods. 2nd ed. Springer; 2004. ISBN: 978-1-4419-1939-7.
[28] Bentley JL. Multidimensional binary search trees used for associative searching. Commun ACM 1975; 18(9): 509517.
[29] Cormen TH, Leiserson CE, Rivest RL. Introduction to algorithms. MIT Press; 1989.
[30] Sedgewick R. Algorithms in C++. Addison-Wesley; 1992.
[31] Shikin EV, Plis LI. Curves and surfaces on a computer screen. Spline guide for users [In Russian]. Moscow: "DI-ALOG-MIFI" Publisher; 1996.
[32] Bicubic interpolation. Source: <https://en.wikipedia.org/wiki/Bicubic_interpolation>.
[33] Moiseev MA, Kravchenko SV, Doskolovich LL. Design of efficient LED optics with two free-form surfaces. Opt Express 2014; 22: A1926-A1935. DOI: 10.1364/0E.22.0A1926.
[34] Kravchenko SV, Moiseev MA, Byzov EV, Doskolovich LL. Design of axisymmetric double-surface refractive optical elements generating required illuminance distributions. Opt Commun 2020; 459: 124976.
[35] TracePro -software for design and analysis of illumination and optical systems. Source: <https://www.lambdares. com/tracepro/>.
Authors' information
Egor Vladimirovich Byzov (b. 1988) graduated with honors (2014) from Samara State Aerospace University named after S.P. Korolyov (now - Samara National Research University named after academician S.P. Korolyov), majoring in Applied Mathematics and Physics. Research interests: design methods of nonimaging optics for LEDs. E-mail: egor.byzov@gmail.com .
Sergey Vasilyevich Kravchenko (b. 1991) graduated with honors (2014) from Samara State Aerospace University named after S.P. Korolyov (now - Samara National Research University named after academician S.P. Korolyov), majoring in Applied Mathematics and Physics. He is a co-author of 18 scientific papers. He received his Candidate's Degree in Physics and Mathematics (2018) from Samara State Aerospace University. E-mail: sergey.kravchenko91@gmail.com .
Mikhail Alexandrovich Moiseev (b. 1986) graduated with honors (2008) from S.P. Korolyov Samara State Aerospace University (now - Samara National Research University named after academician S.P. Korolyov), majoring in Applied Mathematics and Physics. He received his Candidate's Degree in Physics and Mathematics (2011) from Samara State Aerospace University. He is a researcher at the Diffraction Optics laboratory of the IPSI RAS - Branch of the FSRC "Crystallography and Photonics" RAS. He is a SPIE-member. He is coauthor of 23 scientific papers and seven inventions and patents. His current research interests include nonimaging optics. http://www. ipsi. smr. ru/staff/MoiseevM. htm .
Leonid Leonidovich Doskolovich graduated with honors (1989) from S.P. Korolyov Kuibyshev Aviation Institute (presently, Samara National Research University), majoring in Applied Mathematics. He received his Doctor in Physics & Maths (2001) degree from Samara National Research University. Head of laboratory of the IPSI RAS - Branch of the FSRC "Crystallography and Photonics" RAS, professor at Samara University's Technical Cybernetics sub-department. Current research interests include diffractive optics, laser information technologies, nanophotonics. E-mail: leonid@ipsiras.ru .
Received June 1, 2020. The final version - Jule 27, 2020.
