Регуляризация обратной задачи определения параметров нелокального переноса энергии в термоядерной плазме и ее решение методами оптимизации А. В. Соколовa, В. В. Волошиновб Институт проблем передачи информации им. А.А. Харкевича Российской академии наук, 127051, г. Москва, Большой Каретный переулок, д.19 стр. 1 E-mail: a alexander.v.sokolov@gmail.com, б vladimir.voloshinov@gmail.com Рассматривается задача построения параметризованной модели нелокального переноса тепла в вы- сокотемпературной плазме в установках термоядерного синтеза (токамаки, стеллараторы). Эта обратная задача является некорректной и требует регуляризации. Предлагаемый метод ее решения основан на ре- гуляризованной аппроксимации исходных данных вместе с оптимизационной идентификацией парамет- ров модели нелокального теплообмена. Критерием оптимизации является точность выполнения соотно- шений модели плюс мера совпадения с экспериментальными данными и гладкость зависимостей темпе- ратур и определяемых коэффициентов поглощения, «взвешенные» с коэффициентами регуляризации (заранее неизвестными). Для наилучшего выбора указанных весовых коэффициентов применяется про- цедура т.н. перекрестной (взаимной) верификации, в ходе которой часть рядов исходных данных исполь- зуется для «восстановления» остальных. Чем точнее результаты «взаимопроверки» для достаточно ши- рокого набора проверочных тестов, тем «лучше» набор весовых коэффициентов. При переходе к сеточ- ным функциям, заданным в узлах пространственной сетки, каждый этап предлагаемого метода сводится к двухуровневой конечномерной нелинейной задаче математического программирования. Поскольку на «нижнем» уровне требуется решать набор независимых подзадач взаимной верификации, при парал- лельном решении этих подзадач в распределенной среде сервисов оптимизации эта фаза расчетов может быть значительно ускорена. Ключевые слова: обратные задачи, оптимизационная идентификация, регуляризованная аппрокси- мация, двухуровневая оптимизация, распределенные вычисления. Работа выполнена при финансовой поддержке Российского научного фонда (проект № 16-11-10352). © 2016 Александр Витальевич Соколов, Владимир Владимирович Волошинов 441 Введение В данной работе предлагается нестандартный подход к последовательному уточнению фи- зической модели описывающей явления в некоторой пространственной среде. Данный подход является развитием обычных методов оптимизационной идентификации параметров моделей на основе экспериментальных данных с точки зрения учета возможных (часто, заранее неизвестных) неточностей в этих данных. Чем шире доступный набор экспери- ментальных данных и чем выше их точность, тем более точную и подробную модель можно построить. Иногда недостаток количественной информации может компенсироваться наличием значительного объема качественных знаний о закономерностях исследуемого процесса, урав- нений его описывающих и т.д. В данной работе описывается вычислительная схема последовательного уточнения физи- ческой модели, методом оптимизационной идентификации ее параметров с одновременной ре- гуляризацией имеющихся наборов экспериментальных данных. В некотором смысле, предлага- ется схема постепенного усложнения математической модели наблюдаемого физического яв- ления с количественной оценкой качества сделанного усложнения. По принципу - если предсказательная точность новой модели по сравнению с предыдущей повысилась - мы на вер- ном пути. Основной идеей, лежащей в основе вычислительной схемы является «параметризо- ванный» компромисс между «простотой» (характеризуемой гладкостью определяемых функ- ций) и «точностью» (аппроксимации экспериментальных данных). Поэтому Соколов А.В. предложил называть его «метод SvF» (Simplicity vs Fitting). Постановка оптимизационной задачи. Метод SvF. Рассматривается модель нелокального переноса тепла в высокотемпературной плазме в тороидальных камерах установок термоядерного синтеза [Kukushkin, Sdvizhenskii, …, 2015; Kukushkin, Kulichenko, …, 2016]. В упрощенном виде (без учета излучения) модель имеет вид следующего интегро-дифференциального уравнения, описывающего динамику температуры плазмы как функцию расстояния от центра (плазменного пучка) ρ ϵ [0,0.4] (доля диаметра ка- меры) и времени наблюдений t ϵ [2.8, 2.815] (сек):  k  T  r,t  T  r ,t ne  r ,t rdr  k T0  r  T0  r ne  r rdr 1 1 0 k T   ,t   0 1  k T0 (  )  0 1 0 k  T  r,t  ne  r ,t rdr 0 k T0  r  ne  r rdr 0 (1)     T   0  rel   T   ,t T0       3   ,t  .     t Рассмотрим случай, когда электронные плотности ne и ne0, распределение температуры в начальный момент времени Т0 и константа χ0 считаются известными. Для поиска неизвестной функции поглощения k(T), «подкорректированной» функции температуры T(ρ,t) и константы χrel используются экспериментальные данные измерений температуры T(ρ,t) в заданных точках по радиусу и времени: Tij , ,t , i  I , j  J , I  1,N , J  1,N , i j I J (2) где NI = 8, ρi=0.02, 0.19, 0.28, 0.36, 0.43, 0.58, 0.64, 0.7 NJ = 16, tj=2.800+0.001*(j-1) Традиционный подход к решению задач идентификации состоит в выборе неизвестных параметров (функций T(ρ,t) и k(T), числа χrel), которые обеспечивают максимальную (в выбран- ной метрике) близость решения к экспериментальным значениям. Такой подход часто приводит к некорректным постановкам, так как рассматриваемые задачи по своей природе являются об- 442 ратными. Регуляризация их решения может быть проведена с использованием функционалов сложности [Тихонов, 1980] математических моделей, конкретный вид которых зависит как от специфики объекта исследований, так и от применяемого типа моделей, а также объема и каче- ства экспериментальных данных. Для регуляризации задачи поиска неизвестных функций модели (1) по данным (2) выберем в качестве меры близости температуры T(ρ,t) к экспериментальным данным из множества S=IJ следующую функцию   2 1 F (T (  , t ), S )   T T ( i ,t j ) , (3) |S | i , jS ij а в качестве меры сложности функционал, зависящий от вторых производных T(ρ,t) и k(T) 2 2  d 2k    2T  C  k (T ),T (  ,t ),  k ,   ,  t      2 2  dT   2     2  d  dt D Dt    k DT  dT  2 2 (4)   2T    2T   2  t       d  dt t D Dt  t 2  2  d dt D Dt  t  Тогда минимизация целевой функции    G k (T ),T (  ,t ), k ,  ,t ,S  F T (  ,t ),S   C k (T ),T (  ,t ),  k ,   , t  (5) при условии выполнения уравнения (1) для набора данных (2) определяет для различных βk, βρ, βt ≥0 семейство решений с различными соотношениями сложность-близость. Для выбора «разумного» решения из этого семейства (нахождения βk, βρ, βt) воспользу- емся процедурой перекрестного оценивания. Разобьём множество (2) измерений температуры S  I  J на NI (сейчас, 8) частей, соответствующих различным значениями ρi S  (i, j ): jJ  для i  I . (6) i Уберем из множества S одно из подмножеств Si. Решим задачу минимизации (5) на со- кращенном наборе данных S \ Si и проверим, что получилось (верификация) - подсчитаем от- клонение от измерений для подмножества Si. Повторим процедуру для всех |I| подмножеств (6). В итоге получим оценку точности моделирования для заданных βk, βρ, βt. Для получения окончательной оценки погрешности моделировании (и выбора «опти- мального» соотношения сложность-близость) найдем βk, βρ, βt, которые минимизируют оценку погрешности моделирования:  SvF   |S1|  ,min    Tij Ti ( i ,t j )  , 2 2 *  , iI jS k  t i (7) где T * (  , t ) , аппроксимация для массива измерений с удаленным Si, определяются из i min Ti (  ,t ), ki (T )  G ki (T ),Ti (  ,t ),k ,  ,t ,S \ Si  (8) Итак, для заданного разбиения множества измерений (2) на подмножества (6) находим βk, βρ, βt из (7), где функционалы определяются в (3)-(5), а Ti(ρ,t) находятся из (8) с учетом урав- нения модели (1). 443 Особенности программной реализации После замены «непрерывной» функции температуры T(ρ,t) на сеточную и представления искомой функции поглощения излучения K(T) на полином с неизвестными коэффициентами процедура идентификации параметров на базе метода SvF представляет собой задачу двух- уровневой оптимизации – на нижнем уровне для заданного набора коэффициентов β решается NI (8) задач поиска функций Ki (T) и Ti(ρ,t), а на верхнем уровне осуществляется поиск β, мини- мизирующих оценку погрешности моделирования (7). В настоящее время вычислительная схема SvF реализована на языке Python в системе на базе специализированного пакета поддержки расчетов по оптимизационным моделям Pyomo (PYthon Optimization Modeling Objects) [Hart, Laird, …, 2012]. Три года назад, благодаря тому, что авторы популярного алгебраического языка оптимизационного моделирования AMPL в 2005 году «раскрыли» внутренний формат AMPL-стаба, [Gay, 2005], Pyomo стал совместимым со стандартом AMPL. Фактически, принцип расчетов в системе Pyomo полностью повторяет схему применения языка AMPL: модель (в форме набора Python-объектов) вместе с исходными данными (либо в виде Python-объектов, либо в формате файлов с данными AMPL-формате) преобразуются в AMPL-стаб, который обрабатывается AMPL-совместимым решателем. Файл с решением считывается специальной процедурой пакета Pyomo, чтобы обработать результат и/или подготовить данные для новых задач математического программирования. Для решения нелинейных задач математического программирования большой размерно- сти (5) (при условии выполнения «модельного» уравнения (1)) применяется пакет Ipopt, дос- тупный в исходных кодах https://projects.coin-or.org/Ipopt. Подготовка данных для пакета про- изводится пакетом Pyomo. Для решения задачи верхнего уровня (7) относительно трех пере- менных (βk, βρ, βt), пока применяется некоторый эвристический алгоритм. Он состоит в последовательном улучшении коэффициентов вдоль некоторого направления их изменения. Выбор этого направления проводится в результате построения сплайн-интерполяции значений SvF (см. (7)), полученных для предыдущих наборов значений βk, βρ, βt . Поскольку на «нижнем уровне» (взаимной верификации) метода SvF, требуется решить набор независимых задач, соответствующих критерию (8), то на этом этапе вычисления могут быть распараллелены. К настоящему времени создан прототип «адаптера», позволяющего ре- шать наборы независимых задач математического программирования в распределенной среде Everest [Sukhoroslov, Volkov, …, 2015], где развернута система сервисов на основе пакетов чис- ленных методов оптимизации, включая Ipopt. Рис. 1. Функция T(ρ,t) Рис. 2. Функция K(T) 444 Некоторые результаты Решение численного аналога задачи (1)-(8) проводилось для шага по радиусу hρ = 0.02 и шага по времени ht = 0.02. На Рис. 1, 2 приведены графики найденных функций T(ρ,t) и K(T). При этом наблюдался эффект существенного повышения точности модели (1) при очень не- большой (на несколько процентов) коррекции экспериментальных данных (в данном случае – значений температуры Tij, см. (2)). Выводы Проведенные расчеты показали эффективность предложенного метода обработки наблю- дений физического эксперимента динамики температура плазмы и реализации соответствую- щих алгоритмов в среде приложений Everest. Использование метода SvF, основанного на фор- мализации функционалов сложности и применении процедуры перекрестной проверки, позво- лило не только идентифицировать неизвестные (сеточные и полиномиальные) функции, но и оценить достигнутую точность моделирования, с учетом возможной погрешности в предостав- ленных экспериментальных данных. Список литературы Kukushkin A.B., Sdvizhenskii P.A., Voloshinov V.V., Kulichenko A.A. A Model of Recovering the Fast Nonlocal Transport Parameters in Magnetic Fusion Plasmas // EPS 2015: 42nd European Physi- cal Society Conference on Plasma Physics. Mulhouse: European Physical Society. — 2015. — Vol. 39E. — P. 5.182. ISBN 2-914771-98-3. URL: http://ocs.ciemat.es/EPS2015PAP/pdf/P5.182.pdf. Kukushkin A.B., Kulichenko A.A., Sdvizhenskii P.A., Sokolov A.V., Voloshinov V.V. Inverse Problem for Fast Nonlocal Heat Transport Events in Magnetic Fusion Plasmas // EPS 2016: 43rd Europe- an Physical Society Conference on Plasma Physics. Mulhouse: European Physical Society. — 2016. — Vol. 40A. — P. 2.028. ISBN 2-914771-99-1. URL: http://ocs.ciemat.es/EPS2016PAP/pdf/P2.028.pdf. Тихонов А.Н. О математических методах автоматизации обработки наблюдений // В сб. Про- блемы вычислительной математики.— 1980. — С. 3–17. Tikhonov A.N. O matematicheskih metodah avtomatizacii obrabotki nabljudenij [About mathematical methods of au- tomating observation processing] // Sbornik Problemy vychislitel'noj matematiki [Collection Problems of computa- tional mathematics]. — 1980. — P. 3–17. Hart W.E., Laird C., Watson J.-P., Woodruff D.L. Pyomo–optimization modeling in python. Springer. — 2012. — Vol. 67. –– 238 p. Gay D. Writing. nl files. Sandia National Laboratories Technical Report No. 2005-7907P. Technical Report No. 2005-7907P. — 2005. — 14 p. Sukhoroslov O., Volkov S., Afanasiev A. A Web-Based Platform for Publication and Distributed Exe- cution of Computing Applications // Parallel and Distributed Computing (ISPDC), 2015 14th In- ternational Symposium on IEEE. — 2015. — P. 175–184. 445 Optimization and regularization methods in the inverse problem of parametric modeling of a non-local heat transfer in the fusion plasma A. V. Sokolova, V. V. Voloshinovb Institute for Information Transmission Problems of the Russian Academy of Sciences, 19, build.1, Bolshoy Karetny per., Moscow, 127051, Russia E-mail: a alexander.v.sokolov@gmail.com, b vladimir.voloshinov@gmail.com We consider parametric model of the phenomenon of so called non-local heat transport in high temperature magnetic fusion plasmas inside thermonuclear helical devices (tokamaks, stellarators). Underlying mathematical problem, actually - inverse problems, is ill-posed and some regularization is required. We propose a method based on regularized approximation of the given experimental data in conjunction with identification of un- known model parameters via optimization problem. Optimization criterion is to minimize model relation error plus experimental data approximation error and smoothness of spatial temperature functions and energy absorp- tion function which are weighted with some coefficients (unknown in advance). A procedure of, so called, cross- validation is used to find the best values of these coefficients. Actually, this cross-validation consists in restoring of one subset of experimental data via a complementary subset. The better accuracy of cross-validation, the bet- ter choice of weight coefficients. Upon transition to grid-functions, defined in the nodes of spatial mesh, every stage of the method becomes two-level, nonlinear mathematical programming problem (finite dimensional). At the ‘‘low’’ level we need to solve a number of independent optimization problems (cross-validation), so parallel solving may dramatically speeds-up this phase of the algorithm. Keywords: inverse problems, optimization identification, regularized (ridge) approximation, bilevel opti- mization, distributed computing. The work was supported by Russian Science Foundation (project #16-11-10352). © 2016 Alexander V. Sokolov, Vladimir V. Voloshinov 446