CC BY
63
УДК 519.245
Вестник СПбГУ. Сер. 10. 2014. Вып. 1
Л. В. Владимирова
МЕТОДЫ МОНТЕ-КАРЛО В ЗАДАЧЕ ОПТИМИЗАЦИИ ДИНАМИКИ ПУЧКОВ
Санкт-Петербургский государственный университет, 199034, Санкт-Петербург, Российская Федерация
В работе рассматриваются задачи оптимизации динамики заряженных пучков, которые формулируются как задачи программного управления ансамблем траекторий соответствующей динамической системы; при этом критерии качества представляют собой функционалы, заданные на траекториях пучка. Исследуется вопрос численного расчета интегральных критериев специального вида по методу Монте-Карло. Разработан простой вычислительный алгоритм, получена формула для приближенного вычисления интеграла. В качестве примера рассмотрена задача оптимизации продольной динамики пучка в линейном волно-водном ускорителе. Предложены интегральные критерии качества и проведена многокритериальная оптимизация, в ходе которой расчет значений критериев осуществлялся по методу Монте-Карло. Представленные результаты свидетельствуют о повышении качества пучка. Библиогр. 13 назв. Ил. 1.
Ключевые слова: пучок заряженных частиц, оптимизация динамики пучка, интегральные функционалы, метод Монте-Карло, линейный ускоритель, многокритериальная оптимизация.
Vladimirova L. V. Monte Carlo methods in beam dynamics optimization problem // Vestnik of St. Petersburg University. Ser. 10. Applied mathematics, computer science, control processes. 2014. Issue 1. P. 31—39.
In this paper, the problems of charged beams dynamics optimization are considered, when formulated to be trajectory ensemble control problems for specific dynamical systems; quality criteria are presented to be the functionals defined on beam trajectories. The problem under investigation is numerical calculation of special formed integral criteria values by Monte Carlo method. Simple computational algorithm is developed and the formula for integral criterion approximation is obtained. The problem of optimization of longitudinal beam dynamics in waveguide linear accelerator is discussed as an example. Integral quality criteria are suggested, multicriteria optimization is carried out. In the course of optimization criteria values calculation is carried out on the basis of Monte Carlo method. The results presented demonstrate beam quality increment. Bibliogr. 13. Il. 1.
Keywords: charged particle beam, beam dynamics optimization, integral functionals, Monte Carlo method, linear accelerator, multicriteria optimization.
Введение. Задачи оптимизации параметров систем формирования и ускорения заряженных пучков представляют интерес как для создания и эксплуатации соответствующих приборов, так и с математической точки зрения. Д. А. Овсянниковым предложен подход, в соответствии с которым указанные задачи формулируются как задачи управления ансамблями траекторий соответствующих динамических систем. При этом критериями качества управляемого процесса (динамики частиц) являются заданные на траекториях пучка функционалы. Следует отметить, что удачный выбор функционала качества - необходимое условие успешной оптимизации. В монографиях [1, 2] рассмотрен ряд математических моделей управления пучками заряженных частиц, даны постановки задач оптимального управления, развиты математические методы управления пучками и представлены результаты численного моделирования и оптимизации динамики частиц в различных структурах.
© Л. В. Владимирова, 2014
В настоящей работе рассматриваются задачи управления пучками в системах, описываемых дифференциальными и интегродифференциальными уравнениями. В последнем случае динамика частиц описывается с учетом их взаимодействия. Используются интегральные критерии качества управляемого процесса, заданные на траекториях пучка с учетом плотности их распределения. Указанные критерии оцениваются на основе метода Монте-Карло.
Задача оптимального управления ансамблем траекторий. Рассмотрим задачу управления пучком заряженных частиц в линейном волноводном ускорителе. В рамках подхода, предложенного Д. А. Овсянниковым [1,2], введем следующую математическую модель управления. Пусть эволюция пучка частиц описывается диффе-
ренциальными уравнениями
§=/(*,*,«), (1)
^ + х, и) + p(t, x)divxf(t, X, и) = 0 (2)
с начальными условиями
x(0) = xo е Mo, (3)
p(0,x) = po(x). (4)
Здесь t е [0,T] - независимая переменная (число T фиксировано); x е Q С En -n-мерный вектор фазовых координат частицы, Q - открытое множество; u = u(t) - r-мерная вектор-функция управления; f (t, x, u) - n-мерная вектор-функция; Mo - откры-
тое ограниченное множество ненулевой меры; ро(х) - непрерывно дифференцируемая неотрицательная функция, / ро(хо) ¿хо = 1. Полагаем, что начальное положение ча-
Мо
стицы хо есть значение п-мерной случайной величины Хо, распределенной в Мо с плотностью распределения вероятностей ро(хо). Фазовое состояние частицы х(Ь, хо, и) будем рассматривать как значение п-мерной случайной величины Хь с плотностью распределения вероятностей р(£, х). Взаимосвязь плотностей распределения вероятностей ро(хо) и р(Ь,х) определяется уравнением (2) с начальным условием (4).
Считаем, что /, &ух/ непрерывны по совокупности своих аргументов и имеют непрерывные частные производные по х. Класс допустимых управлений Б состоит из кусочно-непрерывных вектор-функций и(Ь), принимающих значения в компакте и С Ег.
При сделанных предположениях для любой точки (Ь,х) € [0,Т] х О существует единственное решение задачи Коши для системы (1) с начальным условием х(£) = х при произвольном допустимом управлении. Будем предполагать, что при любом управлении и € Б решение х = х(Ь,хо,и), соответствующее начальному условию х(0) = хо, определено на всем интервале [О, Т] при любом хо € Мо (здесь Мо - замыкание множества Мо).
Пучок траекторий. В силу сделанных предположений каждому управлению и = и(Ь) можно сопоставить согласно системе (1), (2) семейство траекторий
х = х(Ь, хо, и), (5)
отвечающих всевозможным значениям случайной величины Хо. Семейство траекторий (5) называется пучком (или ансамблем) траекторий, исходящих из множества Мо. Образ множества Мо в силу системы (1), (2) при управлении и = и(Ь) в момент £
Мг,п = {хь = х(г, хо, и) : хо € Мо} , (6)
т. е. множество возможных значений случайной величины Х4, будем называть сечением пучка траекторий, исходящего из М0.
Постановка задачи. Введем функционал
I(и) = J ! Ф(1,Х1)р(1,Х1) ¿XI ¿Ь + J д(хт)р(Т,хт) dxт
(7)
о мм
Мт,,
где х), д(х) - непрерывно дифференцируемые неотрицательные функции. Функционал (7) в зависимости от вида функций Ф(£, х) и д(х) может иметь разный смысл.
Рассмотрим задачу минимизации функционала (7) по управлениям и € О. Это частный случай задачи программного управления ансамблем траекторий системы (1) с учетом плотности их распределения [1,2]. Мы интерпретируем функцию р^,х) как плотность распределения вероятностей случайной величины Х4, характеризующей фазовое состояние частицы при значении £ независимой переменной. Тогда интегралы § Ф(Ь,х1)р(Ь,х1) ¿х^, / д(хт)р(Т,хт) ¿хт могут интерпретироваться как математические ожидания случайных величин Ф(1,Х^ и д(Хт) соответственно. Физический же смысл указанных величин может быть различным. Управление и(0), доставляющее минимум функционалу (7), будем называть оптимальным управлением по отношению к функционалу (7).
Метод Монте-Карло для оценки интегралов. Метод Монте-Карло является по существу методом приближенного интегрирования [3]. Пусть случайная величина X, принимающая значения в некоторой заданной области О С Еп, имеет плотность распределения вероятностей р(х). Введем некоторую случайную величину п = Ф(Х). Ее математическое ожидание равно
МП = 1Ф(х)р(х)
По закону больших чисел [3, 4] с вероятностью 1 имеем
1
N
Ит — > 'ФГя^
NN '
к=1
Ф(х)р(х) ¿х,
(8)
где х^к\ к = 1, Ж, являются независимыми реализациями случайной величины X, распределенной в области О с плотностью вероятностей р(х); N - число реализаций х(к).
Пусть существует конечный второй момент случайной величины п = Ф(Х). Тогда в качестве частного случая центральной предельной теоремы [4] при тех же предположениях относительно х(к), что и выше, находим
Ит Р
N^00
л/К
1 N Г
— / Ф(х)р(х)(1х
к=1
< г
¿и,
и
е
а
здесь а2 = / Ф2(х)р(х) ¿х — / Ф(х)р(х) ¿х I - дисперсия случайной величины п
о \о )
Ф(Х) или с заданной вероятностью 7 неравенство
1 "
^Ф(х(к)) — Ф(х)р(х) ¿х
и_1 ^
к=1 о
га , .
< -!=■ 9
л/П
N
Значение г определяется по заданной доверительной вероятности 7
/ е~"2 дм = \ —7.
Следовательно, если интеграл f Ф(х)р(х) ¿х оценивается с помощью среднего
о
1 N
арифметического — Ф(ж^), то погрешность при заданной вероятности 7 убывает
N к=1
как ° (тм
Формулы (8), (9) служат обоснованием использования среднего арифметического 1 N
— Ф{хЮ) в качестве приближенного значения интеграла / Ф(х)р(х) <1х. N к=1 О
Таким образом, получаем формулу для приближенного вычисления интеграла
1 IV
/ Ф(х)р(х)сЬ « — (10)
О к=1
Подобная вычислительная схема очень часто используется в приложениях. Величина дисперсии, через которую определяется погрешность, может быть тоже оценена. Таким образом, имеется возможность определить число N, гарантирующее необходимую точность с заданной вероятностью 7 (надежностью), и оценить погрешность. Это важное свойство, так как для детерминированных формул приближенного вычисления интеграла не всегда удается даже при небольшой кратности интеграла найти погрешность приближенного вычисления.
Метод Монте-Карло выделяется своей простотой и общностью. Сходимость метода не зависит от кратности интеграла, но медленная сходимость является значительным недостатком метода. Существуют модификации метода Монте-Карло, которые могут обеспечить повышение порядка сходимости [3]. Правда, вычислительная процедура при этом усложняется.
Сходимость метода Монте-Карло есть сходимость по вероятности, однако вероятностные методы в достаточной мере оправдывают себя в практических приложениях.
Применим формулу (10) к оценке критерия качества (7). Оценим сначала второе слагаемое в (7)
Г 1 N
/ д{хт)р{Т,хт)<1хт «.Тт^Э^Р), (11)
^ 1,_1
Ыт. к=1
где к = 1, N, - значения случайной величины Хт, распределенной с плотностью р(Т,хт) на множестве Мти (6).
Заметим, что в задачах оптимизации динамики заряженных пучков величины хТ^, к = 1, -/V, могут быть получены как фазовые состояния N модельных частиц при значении Т независимой переменной. При этом предполагается, что начальные состояния частиц к = выбраны как значения п-мерной случайной величины Хо, рас-
пределенной в области Мо с плотностью ро(х).
Для оценки первого слагаемого в критерии (7) применим к интегралу по Ь формулу Симпсона. Имеем
т
/ / Х()р(Ь, Х() йхг СЬ «
0 Мь,и
~ Т7Т7 ( / Ф(0,х0)р(хо)(1хо+ / Ф{Т,хт)р{Т,хт)<1хт + (12)
2М \Мо,и Мт,и
М-1 N
+ 2Е I ф(ьт ,хт )р(гт,хгт) ¿хт) .
т=1 МЬт,и )
Здесь М(т1и, то = О, М, - сечения пучка траекторий (6), соответствующего управлению и, при значениях Ьт = тТ/М независимой переменной; для каждого т вектор х1т есть значение случайной величины Х1т. Отметим, что Ьо = 0, Ьм = Т.
Будем использовать именно формулу Симпсона, так как в этом случае сумма в правой части (12) включает слагаемые по начальному и конечному сечениям пучка траекторий. Интегралы в (12) по сечениям пучка траекторий М4гг11„, то = 0, М, оценим по формуле (10). В результате получим следующую оценку первого интеграла в (7):
т
/ / Ф(Ь, х^)р(Ь, х^) сМ «
0 Мь,и
(13)
Т / N п\ N п\ М-1 N ,,, \
ЕФМч)+Е^,4ч) + 2 Е ЕФ^хЙ') ■
¿М1\ \к=1 к=1 т=1 к=1 /
Суммируя правые части (11) и (13), находим оценку функционала качества (7).
Замечание. Система (1) с начальным условием (3) обычно интегрируется численно методом Рунге-Кутты, в результате всегда можем получить значения фазовых векторов. Соответствующие случайные величины Х1т распределены с плотностью Р^т, хгт), т = 0, М. Таким образом, задача оценки интегралов в функционале (7) облегчается тем, что не нужно специально моделировать значения случайных величин Х4т, т = (Щ.
Численные результаты. Рассмотрим для примера ускоритель на бегущей волне, а именно ЛУЭ-15-М. Основные характеристики прибора: начальная энергия ^о = 40 кэВ, длина Ь = 0.78 м, А = 0.10 м - длина ускоряющей волны [5-8].
Пусть 2 - продольная координата частицы, £ = — - приведенная продольная коор-
А
дината, 7 - приведенная энергия, ¡3 - приведенная скорость (/З7 = — 1), - фаза частицы. Продольное движение заряженной частицы описывается обыкновенными дифференциальными уравнениями [1, 2]
др
= 2п
В данном случае независимой переменной является £ и фазовый вектор х =
В качестве управляющих функций (компонент двумерного вектора и(£)) используются еЕо(ХрХ
и 1(£) = -т.—, где Ьо{л^) - амплитуда ускоряющей волны, то, е - соответственно
шов2
/е\ л, масса покоя и заряд электрона, гЫС) =--приведенная фазовая скорость уско-
в
ряющей волны.
Множество начальных данных для системы (14)
Мо= (1 + - Д7, 1 + + Д7) X (-7Г,7Г).
у шов2 шов2 у
Здесь 2А7 - начальный разброс энергий электронов на входе в ускоритель.
Осуществляем численное моделирование продольной динамики N модельных частиц. При этом начальные фазовые состояния частиц х^ = (70^, ^о^) ^ Мо, к = 1, М, представляют собой значения случайной величины Хо, распределенной с плотностью
вероятностей ро(хо) = —т— (т. е. равномерно) на множестве Мо-АпА'у
Задача оптимизации продольного движения частиц. Требуется подобрать компоненты управляющего вектора и(£) таким образом, чтобы достичь на выходе ускорителя при £ = Ь (Ь = Ь/А) заданной энергии 7 при максимальном захвате частиц в ускорение.
Задача сводится к минимизации следующих функционалов [5-8]: К\(и) = / Ф^-уь)р(Ь,чь,рь) д^ь дрь,
Мь,и
ь (15)
К2(и) = I I ф2(р$,рц) ди
0 ые,и
где
ФЛ!) = (1 - 1?,
{(р + п)2, р < -п, 0, -п < р <п, (Р - п)2, Р > п-
Функция $1(7) в первом критерии К\(и) имеет вид штрафной функции. Чем меньше разница между энергией на выходе ускорителя 7ь и требуемой энергии 7, тем меньше значение К1(и).
Второй критерий К2 (и) характеризует захват частиц в режим ускорения. Действительно, равенство интеграла нулю означает, что фазы частиц вдоль ускорителя остаются в интервале длиной 2п. Это означает, что частицы не покидают «своего» сгустка.
При решении этой задачи компоненты вектора управления были взяты в виде параметризованных функций
а(0 = и1(£, а) = а1 + а3) +а^, @ф(0 = и^, а) = а5 + а7) +а8 (16)
с вектором параметров а = (а1, а2,..., а%).
Для нахождения оптимизированных управлений проводилась многокритериальная оптимизация [9-13].
Алгоритм оптимизации.
1). Вектор параметров выбирается как реализация равномерно распределенного случайного вектора в параллелепипеде
П = {а : в* < <ц < в**, г = 178}. (17)
Здесь заданные величины в*, в**, г = 1,8, определяют разброс координат вектора параметров а = (а1 ,а2,...,а%). Проводятся . испытаний и находятся . пробных точек
3 = 1, ] 1 в пространстве параметров.
2). Для каждого вектора параметров определяются управления в соответствии с формулами (16).
3). Для каждого управления рассчитываются фазовые состояния N модельных частиц в точках = тЬ/М, то = О, М. Для этого в области Мо выбираются начальные
положения частиц к = численно решаются N задач Коши для системы (14)
(к) - -
и вычисляются фазовые векторы х^ , к = 1^, т = 0,М. При использовании полученных результатов осуществляется оценка критериев качества (15) по методу Монте-Карло по формулам (11), (13). Получаем . точек в пространстве критериев.
4). В пространстве критериев (К.,К2) находим приближенную компромиссную кривую. Соответствующие точки в пространстве параметров определят эффективные управления.
Данные для численного эксперимента.
Число модельных параметров . = 100;
требуемая величина выходной энергии 7 = 10.5;
значения векторов, определяющие параллелепипед (17):
в* = (0.700, 2.935, -3.546,1.294, 0.217,1.993, -2.792,0.657), в** = (0.718, 2.953, -3.528,1.312, 0.236, 2.011, -2.774, 0.675).
Результаты численного эсперимента. Получены параметры, определяющие эффективные управления (16), - это компоненты вектора
аэфф = (0.701, 2.94, -3.55,1.30, 0.218,1.99, -2.78,0.675).
На рисунке представлены фазоэнергетические распределения частиц для некоторого исходного управления, выбранного из физических соображений, и эффективного управления. На нем видно, что при оптимизации уменьшился разброс по фазам и средняя приведенная энергия частиц приблизилась к требуемой. Захват частиц в ускорение увеличился с 0.87 до 0.94.
Заключение. При исследовании многих задач оптимизации динамики заряженных пучков удобно использовать для оценки качества управляемого процесса интегральные функционалы специального вида. А именно, это интегралы по пучку траекторий или его сечению, причем в подынтегральную функцию входит как множитель плотность распределения частиц в фазовом пространстве.
В настоящей работе обоснована возможность и показана эффективность расчета таких интегралов по методу Монте-Карло. При этом указанная плотность трактуется как плотность распределения вероятностей случайной величины, описывающей фазовое состояние частицы.
Потребность в простом и эффективном методе оценки критериев качества вызвана необходимостью многократного их расчета в процессе оптимизации. Благодаря простоте алгоритма метода Монте-Карло вычисление значений критериев не требует значительных затрат времени.
Метод Монте-Карло особенно эффективен при большой кратности интегралов. Достоинством метода является возможность оценки погрешности на основе закона больших чисел. Для детерминированных методов такая оценка затруднительна даже для интегралов малой кратности.
В работе рассмотрена задача оптимизации продольной динамики пучка в линейном волноводном ускорителе. Предложены интегральные критерии качества и проведена многокритериальная оптимизация, в ходе которой расчет значений критериев осуществлялся по методу Монте-Карло. Представленные результаты свидетельствуют о повышении качества пучка.
Литература
1. Овсянников Д. А. Математические методы управления пучками. Л.: Изд-во Ленингр. ун-та, 1980. 276 с.
2. Овсянников Д. А. Моделирование и оптимизация динамики пучков заряженных частиц. Л.: Изд-во Ленингр. ун-та, 1990. 311 с.
3. Ермаков С. М. Метод Монте-Карло и смежные вопросы. М.: Наука, 1975. 472 с.
4. Феллер В. Введение в теорию вероятностей и ее приложения: в 2 т. /пер. с англ. Р. Л. Доб-рушина; под ред. Б. Б. Дынкина. М.: Мир, 1984. Т. 1. 528 с. (Feller W. An introduction to probability theory and its applications.)
5. Владимирова Л. В., Овсянников Д. А., Рубцова И. Д. и др. Оптимизация динамики частиц в ускорительной установке с учетом возбуждения резонатора пучком // Вопросы атомной науки и техники. Сер. Техника физич. эксперимента. 1987. Вып. 4(36). C. 64—67.
6. Владимирова Л. В. Оптимизация динамики пучков взаимодействующих частиц в линейном ускорителе // Изв. РАН. Теория и системы управления. 1995. № 6. С. 178—183.
7. Olemskoy I. V., Rubtsova I. D. Modeling and Optimization of Beam Dynamics in Resonance Bunching and Decelerating Systems // Proc. of the First Intern. Workshop: Beam Dynamics & Optimization. St. Peterburg, 1995. P. 143-153.
8. Ovsyannikov D. A., Ovsyannikov A. D., Vorogushin M. F., Svistunov Yu. A., Durkin A. P. Beam Dynamics Optimization: Models, Methods and Applications // Nuclear Instruments and Methods in Physics Research. Section A: Accelerators, Spectrometers, Detectors and Associated Equipment. 2006. Vol. 558, N 1. P. 11-19.
9. Соболь И. М., Статников Р. Б. Выбор оптимальных параметров в задачах со многими критериями. М.: Наука, 1981. 112 с.
10. Подиновский В. В., Ногин В. Д. Парето-оптимальные решения многокритериальных задач. М.: Наука, 1982. 256 c.
11. Губанов В. А., Захаров В. В., Коваленко А. Н. Введение в системный анализ. Л.: Изд-во Ленингр. ун-та, 1988. 227 с.
12. Vladimirova L. V. On Application of Multicriterial Optimization in Partical Beam Control Problem // Proc. of the 11th Intern. IFAC Workshop: Control Applications of Optimization. (July 3-6, 2000, St. Petersburg, Russia) / ed. V. Zakharov. Oxford, UK: Pergamon Press, 2001. Vol. 1. P. 363-367.
13. Владимирова Л. В. Многокритериальная оптимизация в прикладных задачах: учеб. пособие (для студентов и аспирантов физ.-мат. специальностей). СПб.: Изд-во С.-Петерб. ун-та, 2002. 24 c.
Статья рекомендована к печати проф. Д. А. Овсянниковым. Статья поступила в редакцию 31 октября 2013 г.
Контактная информация
Владимирова Людмила Васильевна - кандидат физико-математических наук, доцент; e-mail: [email protected]
Vladimirova Ludmila Vasilievna - candidate of physical and mathematical sciences, reader, St. Petersburg State University, 199034, St. Petersburg, Russian Federation; e-mail: [email protected]
