=Paper=
{{Paper
|id=Vol-2022/paper18
|storemode=property
|title=
Обратные задачи моделирования на основе регуляризации и распределенных вычислений в среде Everest
(Inverse Problem in the Modeling on the Basis of Regularization and Distributed Computing in the Everest Environment)
|pdfUrl=https://ceur-ws.org/Vol-2022/paper18.pdf
|volume=Vol-2022
|authors=Alexander Afanasiev,Vladimir Voloshinov,Alexander Sokolov
|dblpUrl=https://dblp.org/rec/conf/rcdl/AfanasievVS17
}}
==
Обратные задачи моделирования на основе регуляризации и распределенных вычислений в среде Everest
(Inverse Problem in the Modeling on the Basis of Regularization and Distributed Computing in the Everest Environment)
==
https://ceur-ws.org/Vol-2022/paper18.pdf
Обратные задачи моделирования на основе
регуляризации и распределенных вычислений
в среде Everest
© А.П. Афанасьев © В.В. Волошинов © А.В. Соколов
Институт проблем передачи информации им. А.А. Харкевича РАН,
Москва, Россия
alexander.afanasyev@gmail.com vladimir.voloshinov@gmail.com
alexander.v.sokolov@gmail.com
Аннотация. Изложена методика оценки математических моделей физических явлений, проис-
ходящих в некоторой пространственной среде, на основе рядов экспериментальных данных. Целевая
функция в обратной оптимизационной задаче идентификации параметров модели включает регуля-
ризирующее слагаемое с неизвестными весовыми коэффициентами при вторых производных функ-
ций, описывающими исследуемое явление. Для выбора этих весовых коэффициентов применена
процедура перекрестной (взаимной) верификации, когда часть исходных экспериментальных данных
используется для «восстановления» остальных. Чем точнее результаты «взаимопроверки» для доста-
точно широкого набора проверочных тестов, тем «лучше» набор весовых коэффициентов. При
вы-боре направления их улучшения требуется решить большой набор вспомогательных подзадач
мате-матического программирования, для чего предложено использовать распределенную систему
серви-сов оптимизации на платформе Everest.
Ключевые слова: обратные задачи, регуляризация, распределенные вычисления, REST-
сервисы, платформа Everest.
Inverse Problem in the Modeling on the Basis of Regulariza-
tion and Distributed Computing in the Everest Environment
© A.P. Afanasiev © V.V. Voloshinov © A.V. Sokolov
Institute for Information Transmission Problems RAS (Kharkevich institute),
Moscow, Russia
alexander.afanasyev@gmail.com vladimir.voloshinov@gmail.com
alexander.v.sokolov@gmail.com
Abstract. A method for estimating mathematical models of physical spatial phenomena is presented.
Estimating is based on the series of experimental data. The objective function in the inverse optimization
problem of identification of model parameters includes a regularizing term with unknown weight coeffi-
cients for the 2nd derivatives of the spatial function describing the phenomenon. Successive cross-validation
procedure is used to choose values of weight coefficients. This cross-validation consists in approximation of
one subset of experimental data by processing of a complementary subset. The better accuracy of the “cross-
approximation”, the better set of weight coefficients. Choosing direction of the possible improvement re-
quires solving a number of subsidiary optimization problems. For that it is proposed to use distributed com-
puting environment of optimization services deployed via Everest toolkit.
Keywords: inverse problems, regularized (ridge) approximation, distributed computing, Everest plat-
form.
Введение пространственной среде, где для искомых зависимо-
стей имеет смысл понятие «гладкости» по перемен-
В работе предложен нестандартный подход к по- ным, описывающим состояние модели. Речь может
следовательному уточнению математических моде- идти о физических процессах, которые моделируют-
лей, описывающих физические явления в некоторой ся дважды дифференцируемыми функциями на пря-
мой, плоскости или в трехмерном пространстве.
Труды XIX Международной конференции «Ана- Данный подход является развитием методов оп-
литика и управление данными в областях с ин- тимизационной параметрической идентификации на
тенсивным использованием данных» (DAMDID/ основе экспериментальных данных с учетом воз-
RCDL’2017), Москва, Россия, 10–13 октября 2017 можных (и неизвестных) неточностей в этих дан-
года ных. Чем шире доступный набор измерений (экспе-
100
риментальных данных) и выше их точность, тем
более точную и подробную модель можно постро- (
∆ F f (⋅ xy ), K = ) K1 ∑ (zk − f ( xk , yk ))2 ; (3)
k∈K
ить. Иногда недостаток количественной информа-
ции можно компенсировать дополнительными зна- - характеристику негладкости
ниями о закономерностях исследуемого процесса,
уравнений его описывающих и т. д. (
∆ S f (⋅ xy ), β x , β y = )
Аналогом является обработка статистических ∂2 f
2
∂2 f
2
данных на основе сплайн-аппроксимации («сгла-
женных» кубических сплайнов), где коэффициент
β x2 ∫∫
XY
∂x
2
dx dy + 2 β x β y
XY
∫∫
∂x∂y dxdy + (4)
штрафа за интеграл квадрата 2-й производной ап- 2
проксимирующей зависимости играет роль пара- ∂2 f
метра регуляризации [3, 4, 9, 12, 13]. Особенность β y2 ∫∫
∂y 2 dxdy
предлагаемого подхода – минимизация суммы от- XY
клонения от экспериментальных данных и штрафа
(для сеточной функции вторые производные заме-
за «негладкость» при дополнительных ограничениях
няются разностными выражениями);
(соотношениях исследуемой модели).
- компромиссный критерий (функционал Тихонова
По сути предлагается схема постепенного уточ-
нения математической модели наблюдаемого физи- [10]), зависящий от параметров βx, βy и множества
ческого явления с количественной оценкой качества измерений, характеризуемого набором индексов K:
такого уточнения. Если предсказательная точность
∆ ( f (⋅xy ), β x , β y , K ) = ∆ F ( f (⋅xy ), K ) + ∆ S ( f (⋅xy ), β x , β y ). (5)
новой модели повысилась, мы на верном пути.
Будем искать функцию f(x,y) в классе дважды
1 Описание метода
непрерывно-дифференцируемых функций, решая
Опишем общую схему метода на примере по- следующую задачу минимизации, где переменными
строения модели, описывающей исследуемое явле- являются либо значения функции «на сетке» по пе-
ние с помощью переменных x, y, z. Здесь z является ременным x,y, либо параметры f(x,y):
измеряемой характеристикой, зависящей от x и y.
Требуется построить модель в форме явной и неяв- ∆ ( f (⋅ xy ), β x , β y , K ) :
ной зависимостей между указанными переменными:
f (⋅ xy ) = Arg min M (x, y, f ( x, y ) ) = 0,
* (6)
z = f ( x, y ), x ∈ X , y ∈ Y ; M ( x, y, z ) = 0. (1)
f ( ⋅ xy )
x ∈ X , y ∈ Y .
Здесь f(x,y) – неизвестная функция, которую и тре-
Решение задачи (6) для различных значений βx,
буется «восстановить» при дополнительных пред-
положениях о «физике» наблюдаемого явления, X, Y
βy дает функции f(x,y), соответствующие различным
соотношениям между «гладкостью» и «точностью»
– множества (интервалы) допустимых значений со-
(восстановления экспериментальных данных). Для
ответствующих переменных. Предполагаемая и
поиска «разумного» баланса, выраженного в значе-
подлежащая верификации физическая модель пред-
ставлена вторым уравнением (1). Неявная зависи- ниях βx, βy, воспользуемся процедурой перекрестной
мость M определяет дополнительные связи между верификации [3, 4]. Для этого разобьем множество
переменными. Их может быть несколько, т. е. M – индексов экспериментальных данных на набор не-
вектор-функция. Будем искать функцию f(x,y) либо в пересекающихся подмножеств
виде значений на достаточно «мелкой» сетке по пе-
ременным x,y, либо в классе функций, зависящих от
K= K , K K = ∅, i ≠ j .
i∈I i i j (7)
параметров, значения которых и нужно определить, Уберем из множества K одно из подмножеств Ki.
решив обратную задачу. Решим задачу минимизации (6) на оставшемся
Пусть у исследователя имеется набор экспери- наборе данных K∖Ki. Пусть ее решение f K*i ( x, y ) :
ментальных данных (измерений) следующего вида:
∆ ( f (⋅xy ), β x , β y , K \ K i ) :
{zk , xk , yk }, k ∈ K, K = 1, k max , (2) *
Ki
f (⋅xy ) = Arg min M ( x, y, f ( x, y ) ) = 0,
(8)
где значения zk измерены с некоторыми неизвест-
f ( ⋅ xy )
x ∈ X , y ∈Y .
ными погрешностями. Задача восстановления функ-
ции f часто является некорректной и требует регуля- Определим отклонение f K*i ( x, y ) от измерений
ризации. Будем искать зависимость f(x,y) (в пара-
метризованном или «сеточном» виде) методом регу-
∑ (z − f ( x , y )) . По-
2
для k ∈ K i по формуле k
*
Ki k k
ляризованной идентификации SvF, Smoothness-vs- k∈K i
Fitting, состоящем в поиске компромисса между
точностью совпадения с экспериментальными дан- вторив эту процедуру для всех подмножеств Ki, i∈I,
ными и гладкостью искомой функции f(x,y). получим «перекрестную» оценку точности модели
для заданных βx, βy:
Для формализации такого подхода введем:
- характеристику совпадения с измерениями
101
жить M(x,y,z) тождественно равной нулю). Для
F (β x , β y ) = (σ SvF (β x , β y )) =
1
∑ ∑
K i∈I k∈K i
2 2
(
zk − f K*i ( xk , yk ) . (9) ) функции одной переменной x требуется лишь один
Для получения минимальной погрешности моде- скалярный параметр β. При β→0 задача (6) стано-
лировании и выбора «оптимального» соотношения вится задачей сплайн-интерполяции, решением ко-
торой является т. н. кубический сплайн, т. е. функ-
близость–сложность будем искать βx, βy, которые
ция, имеющая минимальный (см. (4)) интеграл
(
минимизируют величину σ SvF β x , β y . Для найден- ) квадрата 2-ой производной и проходящая через все
ных в результате минимизации значений βx, βy ис- «экспериментальные» точки (2), Рис. 1b).
комая функция f*(x,y) определяется в результате ре-
шения основной задачи (6). Итоговая погрешность
аппроксимации (невязка) определяется по формуле
(σ ) = K1 ∑ (z − f ( x , y )) .
* 2
k
*
k k
2
(10)
k ∈K
Рисунок 1b Кубический сплайн, β→0, и значение
β*, полученное методом SvF
При β→∞ мы получим линейную функцию, ми-
Рисунок 1a Исходные данные и результат метода нимизирующую сумму квадратов отклонения от
наименьших квадратов β→∞ исходных данных, Рис. 1a. «Компромиссный» ре-
Таким образом, процедура расчетов соответству- зультат сплайн-аппроксимации со штрафом (за не-
ет двухуровневой задаче оптимизации: гладкость), определенным методом SvF, представ-
лен на Рис. 1b снизу.
( ) ∑ ∑( )
1 2
Φ βx,βy = z k − f K* i ( x k , y k ) → min , (11) Опишем процедуру решения задачи верхнего
K i∈I k∈K β x ,β y ≥0
i уровня (11). Введем обозначение P(a,β ) для произ-
вольного многочлена 2-го порядка вектора перемен-
где функции f K*i ( x, y ) определяются как решения ных β , зависящего от вектора коэффициентов a
независимых задач оптимизации (8).
P(a, β) = axx β x2 + axy β x β y + a yy β y2 + ax β x + a y β y + a0 , (12)
Результат метода SvF находится, в некотором
смысле, «между» результатами применения хорошо ( ) ( )
где a = a xx , a xy , a yy , a x , a0 , β = β x , β y .
известных методов наименьших квадратов и куби-
ческой сплайн-интерполяции. На Рис. 1 метод SvF Далее алгоритм строит последовательность зна-
продемонстрирован на примере восстановления чений βν, ν=1,2,… . Пусть, на N-ом шаге построены
функции одного переменного по набору исходных значения βν, ν=1:N, для которых вычислены значе-
данных (Рис. 1a) при отсутствии дополнительных ( )
ния Φν = Φ βν . Без ограничения общности (воз-
модельных соотношений (формально можно поло- можно, после перенумерации) можно считать, что
102
Φ N – минимальное из полученных значений. Будем мер, в экспериментах наблюдается «быстрый нело-
трактовать {
Φν , βν ν =1
N
}
как набор точек в R3. Рас-
кальный перенос тепла» – практически мгновенное
(по сравнению с «классической» тепловой диффузи-
смотрим задачу аппроксимации этих точек графи- ей) повышение температуры в центре плазменного
ком многочлена 2-го порядка (12), причем чем век- шнура после охлаждения его периферии или обрат-
тор βν ближе к «наилучшему» β N, тем с большим ный процесс – мгновенное понижение температуры
весом он будет учитываться. Кроме того, будем в центре при быстром нагреве периферии плазмы.
штрафовать за кривизну построенной функции с Здесь мы не будем обсуждать неясную физику
коэффициентом m (подлежащим выбору): этого явления. Ограничимся демонстрацией приме-
нения метода к предварительной обработке экспе-
∑e
− βn −β N
(Φ − P(a, β )) + m (a + a + a ) → min . (13)
n n 2 2
xx
2
xy
2
yy
a
риментальных данных из статьи [8], где обсуждает-
n=1:N
ся явление «быстрого» охлаждения плазмы в стел-
Пусть a*(m) оптимальное значение вектора перемен- лараторе LHD, http://www.lhd.nifs.ac.jp/en
ных в этой задаче выпуклого программирования. Исходными данными являются графики зависи-
Выберем значение m, минимизируя отклонение мостей температуры плазмы от расстояния (ρ) до
значения аппроксимирующего полинома от центра тороидальной камеры в различные моменты
времени (t), Рис. 2.
наилучшего значения Φ N . Тем самым для аппрок-
симации получим вновь двухуровневую задачу:
( )
P a* ( m ), β N − Φ N → min ,
m ≥0
(14)
где a*(m) – решение задачи «нижнего уровня» (13).
Здесь в задаче верхнего уровня нужно выбрать
единственную переменную m, а задача нижнего
уровня эффективно разрешима благодаря ее выпук-
лости. Поэтому для поиска оптимального штрафно-
го коэффициента m* применим тот или иной алго-
ритм минимизации функции одного переменного.
После аппроксимации зависимости Φ (β ) квад-
ратичной функцией P a* ( µ * ), β новое значение ( )
вектора β находится в результате решения следую-
Рисунок 2 Температура на разных расстояниях от
щей вспомогательной задачи, которая напоминает
центра тороидальной камеры T(ρ,t), [8]
метод линеаризации Пшеничного–Данилина [11],
применяемый к минимизации функции P a* ( µ * ), β ( ) Множеством K экспериментальных данных о
по β: температуре, см. (2), является набор пар индексов iρ,
it, значений расстояния и моментов времени. Разде-
(∇ P(a ( m ), β )) (β − β )+ 12 β − β
β
* * N T N N
2
→ min,
β ≥0
лим множество измерений температуры на 8 частей,
определяемых значениями ρ (горизонтальные груп-
где ∇β P () обозначает градиент многочлена (12) пы точек на Рис. 2). Для перекрестной верификации
по переменным βx, βy. Заметим, что решение этой будем использовать 6 множеств (см. (7)):
выпуклой задачи квадратичного программирования {( ) }
K i = iρ , it : it ∈ I t , iρ = 2 : 7. (15)
существует и единственно. Вычислим значение
Крайние значения (ρ1=0.02 и ρ9=1) в перекрестном
( )
Φ N +1 = Φ β N +1 , решив набор задач нижнего уровня оценивании не учитывались, т. к. в этой задаче важ-
(8). Если величина Φ β N +1 − Φ β N ( ) ( ) меньше неко- на интерполяция значений температуры плазмы.
Ниже приведен ряд моделей наблюдаемого явле-
торого порогового значения, работа алгоритма пре-
ния вида (1). Они отличаются друг от друга предпо-
кращается. Если это не так, то схема расчетов по-
ложениями о процессе распространения тепла в
вторяется для расширенного набора Φν , βν ν =1 . { } N +1
плазме в форме второй группы соотношений (1).
Цель – выбрать модель, которая бы максимально
3 Демонстрация метода при моделирова- точно описывала функцию T(ρ,t) – зависимость тем-
нии распространения тепла в высоко- пературы плазменного пучка от ρ и t. Можно ожи-
температурной плазме дать, что при «правильных» уточнениях предсказа-
тельная точность моделирования (по отклонению
Физические явления в горячей плазме, удержи-
функции T(ρ,t) от значений на Рис. 2) должна улуч-
ваемой сильным магнитным полем в тороидальных
шиться. Результаты расчетов по всем моделям при-
вакуумных камерах установок термоядерного синте-
ведены ниже в Таблице 1. Графики, иллюстрирую-
за, таких, как ТОКАМАКи и стеллараторы, важны
щие расчеты по моделям, доступны на странице
для перспектив термоядерной энергетики. Напри-
http://distcomp.ru/~vladimirv/damdid2017.
103
3.1 Модель 1 (сплайн-интерполяция без SvF) где T0(ρ) – функция температуры пучка в начальный
момент времени предполагается известной. В такой
Пусть у нас нет никаких предположений о физи-
постановке критерий (5) имеет вид (17), но, кроме
ческих причинах наблюдаемого явления. Будем ис-
коэффициентов βTρ , βTt, βπρ и βπt, минимизация про-
кать функцию T(ρ,t), проходящую через все точки
используемого набора измерений и имеющую ми- водится еще по переменной χ. Результаты модели-
нимальную кривизну в смысле формулы (4), а точ- рования показывают, что учет только диффузионно-
ность определять процедурой перекрестной верифи- го члена (без дополнительного «источника») не яв-
кации на подмножествах (15). Формально этот при- ляется удачным: точность моделирования падает.
ем соответствует случаю, β→0, Рис. 1b. 3.5 Модель 5 (Метод SvF, тепловая диффузия и
Получили задачу сплайн-интерполяции. При до- «неизвестный источник»)
статочно «необременительных» предположениях о Предположим, что наряду с «медленной» тепло-
расположении узлов интерполяции её решение су- вой диффузией в плазме действует некоторый до-
ществует и единственно [13]. полнительный механизм переноса тепла, который в
3.2 Модель 2 (Сплайн-аппроксимация) следующей формуле обозначен S(ρ,t):
Не имея, как и выше, никаких предположений о
«физике» явления, применим метод SvF, когда
∂T
(ρ , t ) = χ ∂ ρ ∂ (T (ρ , t ) − T0 (ρ )) + S ( ρ , t ). (19)
∂t ρ∂ρ ∂ρ
функция M(x,y,z) тождественно равна нулю (см. (1)
и Рис. 1b, снизу). Эта задача является задачей Здесь функция критерия (5) имеет вид, аналогичный
сплайн-аппроксимации с выбором штрафа за не- (17), но зависит от четырех β-коэффициентов для
гладкость на основе перекрестного оценивания. Как функций T(ρ,t) и S(ρ,t):
и для Модели 1, её решение существует и един-
ственно [13]. Несмотря на то, что построенная ∆ ({T (⋅ρt ), S (⋅ρt )}, βTρ , βTt , β Sρ , β St , K ) =
функция температуры уже не проходит через узлы, = ∆ F (T (⋅ρt ), K ) + ∆ S (T (⋅ρt ), βTρ , βTt ) + (20)
оценка точности заметно улучшилась (см. Табл. 1).
+ ∆ S (S (⋅ρt ), β Sρ , β St ).
3.3 Модель 3 (Метод SvF и простое дифференци-
альное уравнение) В итоге определяются коэффициенты βT,ρ , βT,t, βS,ρ и
βS,t, и коэффициент «теплопроводности» χ. Точность
Предположим, что процесс описывается про- моделирования заметно возросла (Табл. 1).
стейшим дифференциальным уравнением:
3.6 Сравнение результатов расчетов
∂T
(ρ , t ) = π (ρ , t ) (16) Как уже было объявлено, результаты расчетов
∂t
сведены в Таблицу 1. Графические изображения
Здесь неизвестными являются две функции T(ρ,t) и построенных зависимостей от ρ и t приведены на
π(ρ,t) от двух переменных. Поэтому в схеме SvF рисунках в http://distcomp.ru/~vladimirv/damdid2017.
проводится оптимизация по четырем коэффициен- Таблица 1 Результаты моделирования
там штрафа «за негладкость» обеих функций: βTρ, βTt Погрешность
и βπρ, βπ,t. Поскольку для неизвестной функции Модель Кросс- Аппроксимации,
π(ρ,t), правой части дифференциального уравнения валидации, % СКО 1, %
2
(16), нет экспериментальных данных, то основной 1 11.94 0
компромиссный «критерий» (5) имеет вид 2 9.25 3.55
({ } )
∆ T (⋅ ρt ), π (⋅ ρt ) , β Tρ , β Tt , β πρ , β πt , K =
3
4
9.19
10.00
3.55
8.53
( ) (
∆ F T (⋅ ρt ), K + ∆ S T (⋅ ρt ), β Tρ , β Tt + ) (17) 5 (χ=0.21) 1.85 0.59
∆ S (π (⋅ ρt ), β πρ , β πt )
Заметим, что переход от «простых» сплайнов к
Сравнение с предыдущей моделью (см. Табл. 1) более содержательным моделям 2–4 дал незначи-
не выявило принципиальных изменений. тельное улучшение точности перекрестной провер-
ки. Но расчет по Модели 5 дал многократное улуч-
3.4 Модель 4 (Метод SvF и тепловая диффузия) шение показателей точности, т. е. выделение неиз-
Здесь предполагается, что распространение тепла вестного «источника» (переносчика тепловой энер-
описывается классическим дифференциальным гии) является, по-видимому, верным уточнением
уравнением тепловой диффузии в цилиндрически модели. Приведенный пример демонстрирует при-
симметричной среде (как приближении тороидаль- менение метода для количественной проверки «ка-
ной камеры) с заранее неизвестным коэффициентом чества» различных математических моделей на
«теплопроводности» χ. имеющемся наборе экспериментальных данных.
∂T
(ρ , t ) = χ ∂ ρ ∂ (T (ρ , t ) − T0 (ρ )), (18)
∂t ρ∂ρ ∂ρ
1
Среднеквадратичное отклонение
2
Сплайн-интерполяция
104
гранжа при ограничениях). Формат этого файла со-
4 Возможности программной реализации ответствуют стандарту применяемого AML, и AML-
в среде Everest транслятор может его импортировать.
Предлагаемая методика основана на решении за- Также эти языки позволяют описывать сложно-
дач математического программирования. Для ее составные сценарии расчетов по оптимизационным
практического применения нужен набор программ- моделям: содержащие условные переходы, циклы,
ных инструментов для: 1) описания указанных за- динамическое формирование новых задач на основе
дач; 2) формирования структур данных, соответ- результатов решения предыдущих, создание набо-
ствующих отдельным экземплярам таких задач; ров задач для различных наборов значений парамет-
3) отправки этих данных пакетам численных мето- ров и т. п. Являясь по назначению языками про-
дов (решателям), для поиска решения; 4) обработки граммирования высокого уровня, AMPL и GAMS не
результатов работы решателей, например, для изме- отвечают требованиям, предъявляемым даже к про-
нения метода расчетов. Сложившаяся практика цедурным языкам (надо иметь ввиду «почтенный
применения оптимизационных моделей предусмат- возраст» языков, AMPL и GAMS появились в конце
ривает два способа организации расчетов. 1970-х годов). Например, в них нет понятия проце-
Первый, «низкоуровневый», на основе открытого дуры-функции, все переменные (кроме внутренних
программного интерфейса (API) решателя для неко- индексов циклов или операторов «итерирования»)
торого языка программирования (C/C++, C#, Java, являются глобальными и т. п.
Python и т. п.). Подготовка данных для отправки В связи с этим большой интерес вызывает систе-
решателю и обработка результатов производятся ма оптимизационного моделирования Pyomo (PY-
обычно на языке API решателя. Такой подход, хотя thon Optimization Modeling Objects) [2]
и может привести к созданию высокопроизводи- http://pyomo.org, основанная на популярном объект-
тельной системы расчетов, является достаточно но-ориентированном языке программирования Py-
трудоемким и требует привлечения высококвалифи- thon (Pyomo представляет собой специализирован-
цированных программистов. Кроме того, изменения ный Python-пакет). Четыре года лет назад Pyomo
в схеме расчетов или структуре применяемой опти- стал совместимым со стандартом AMPL. Это про-
мизационной модели требуют переписывания зна- изошло благодаря тому, что авторы AMPL 12 лет
чительных фрагментов программного кода. назад «раскрыли» внутренний формат AMPL-стаба.
Второй, «высокоуровневый», подход использует Принцип расчетов в системе Pyomo повторяет
алгебраические (декларативные) языки оптимизаци- схему применения языка AMPL: модель (в форме
онного моделирования, что гораздо удобнее, осо- набора Python-объектов) вместе с исходными дан-
бенно для поисковых исследований. Развитие таких ными (либо в виде Python-объектов, либо в формате
языков (AML, Algebraic Modelling Language – в ан- файлов с данными AMPL-формате) преобразуются в
глоязычной литературе) ведется уже более 30 лет. AMPL-стаб, передаваемый AMPL-решателю. Файл с
До настоящего времени наиболее популярными яв- решением можно считать специальной процедурой
ляются AMPL, GAMS. Основными составляющими пакета Pyomo, для оформления результата и/или
AML-систем являются: собственно язык для описа- подготовки исходных данных новых задач матема-
ния оптимизационных моделей, средства автомати- тического программирования.
ческого дифференцирования и унифицированный
4.1 Сведения о платформе Everest
интерфейс взаимодействия с пакетами.
AML-языки позволяют записать соотношения Разработка программного обеспечения для со-
оптимизационных задач (параметры, переменные, здания систем на основе REST-сервисов ведутся в
целевую функцию, ограничения, индексы парамет- нашем коллективе около шести лет. Первоначальная
ров, переменных и ограничений и т. п.), разделив и последующая стабильная версии ПО имели назва-
«символьное описание» задачи (т. н. модельное ние MathCloud, mathcloud.org. Последние три года
представление) и «конкретные данные» (наборы велась активная работа по переходу на новую вер-
индексов, значения числовых параметров и т. п.). сию программного инструментария, т. н. Everest [6,
Символьная модель и конкретные данные, пред- 7], http://everest.distcomp.org.
ставленные обычно текстовыми файлами, обрабаты- Программный инструментарий Everest является
ваются специальным транслятором. На выходе – системой с открытым кодом, свободно доступным
специальная структура данных в виде т. н. стаб- на популярном портале gitlab.com. Семантика Ever-
файла (stub, в терминологии AMPL), готового для est базируется на следующей иерархии понятий:
передачи решателям, совместимым с языком моде- • Приложение Everest – REST-сервис с REST-
лирования. Для нелинейных задач стаб также со- интерфейсом в формате JSON; приложение Ever-
держит правила вычисления первых и вторых про- est, вообще говоря, является абстракцией, для
изводных всех функций задачи математического которой не выделено никакого реального вычис-
программирования. Если численный метод находит лительного ресурса (его выбор и подключение
решение, то AML-совместимый решатель возвраща- происходят непосредственно перед вызовом);
ет файл, содержащий значения всех «прямых» и
двойственных переменных задачи (множителей Ла- • Вычислительный ресурс Everest – реальное вы-
числительное устройство или инфраструктура
105
(сервер, кластер, грид, облако), где производится 3. Подсистема балансировки вычислительной
обработка данных Everest-приложениями; нагрузки управляет выполнением заданий, распре-
• Агент доступа к вычислительным ресурсам Ev- деляя их между вычислительными ресурсами, под-
erest – программный модуль (на Python), обеспе- ключенными к одному приложению. Пользователь
чивающий подключение ресурса к системе Ever- может не знать, какие именно ресурсы обрабатыва-
est (некоторые основные типы ресурсов пред- ют его данные. Эта подсистема постоянно совер-
ставлены на Рис. 3); шенствуется разработчиками Everest, в частности,
ожидается возможность выбора различных политик
• сервер Everest (контейнер приложений) – цен- распределения заданий между ресурсами.
тральный сервер для: сохранения дескрипторов
всех приложений; регистрации пользователей; 4. Система контроля доступа к приложениям и
управления правами доступа к созданным при- защиты данных в Everest использует две техноло-
ложениям; управление очередями заданий; гии: защищенный обмен данными между Everest-
сервером и агентами по протоколу HTTPS/SSH;
• веб-интерфейс работы с сервером Everest, вклю- специальные «временные» ключи, т. н. токены, вы-
чающий средства создания приложений, запуска даваемые зарегистрированным пользователям Ever-
и контроля за ходом выполнения заданий. est с ограниченным сроком действия (7 дней).
Предъявление токенов обязательно и для работы с
веб-интерфейсом сервера, и при вызове приложений
через Everest API.
4.2 Сервис оптимизации
Базовым сервисом решения задач оптимизации в
Everest является сервис solve-ampl-stub решения
задач математического программирования, пред-
ставленных в виде AMPL-стаб-файла [1], пакетом
численных методов (решателем), указанным при
обращении к сервису. Сейчас сервис обеспечивает
унифицированный доступ к следующим пакетам,
позволяющим решать основные типы задач матема-
тического программирования
Рисунок 3 Архитектура программного инструмен- (LP/MILP/NLP/MINLP):
тария Everest • Ipopt (Coin-OR Interior Point Optimizer, NLP),
https://projects.coin-or.org/Ipopt;
Перечислим ряд особенностей Everest, важных с
точки зрения практического развертывания и при- • CBC (Coin-OR Branch-and-Cut, LP, MILP),
менения систем оптимизационного моделирования в https://projects.coin-or.org/Cbc;
распределенной вычислительной инфраструктуре. • SCIP (Solving Constraint Integer Programs, LP,
1. Автор приложения может «незаметно» для MILP, MINLP (билинейные невыпуклые)),
пользователей повышать/понижать вычислительную http://scip.zib.de
«мощность» приложений (сервисов), изменив спи- • Bonmin (COIN-OR Basic Open-source Nonlinear
сок ресурсов (фактически агентов доступа к ресур- (convex) Mixed Integer programming, MINLP),
сам), ассоциированных с данным приложением. https://projects.coin-or.org/Bonmin
Например, если комплект решателей будет установ- Данным приложением можно пользоваться как
лен на новом вычислительном сервере вместе с через его веб-интерфейс, так и через программный
агентом Everest, то производительность Everest- интерфейс Everest Python API.
приложения для решения задач оптимизации повы-
сится. 4.3 Сведения о программной архитектуре Pyomo-
2. Платформа Everest предлагает унифицирован- Everest
ный программный интерфейс (Everest Python API), Основным требованием к системе было: обеспе-
gitlab.com/everest/python-api. Он позволяет клиент- чить возможность выполнения любых программ
ским модулям на Python взаимодействовать с серви- (сценариев расчетов) на языке Python с использова-
сами Everest по модели асинхронных вызовов уда- нием пакета Pyomo и AMPL-совместимых решате-
ленных объектов и программировать вычислитель- лей в среде сервисов оптимизации Everest, возмож-
ные сценарии координированной обработки данных но, после некоторой модификации самой программы
несколькими приложениями Everest. При этом неза- согласно определенному набору правил. Общий
висимые задания будут выполняться одновременно принцип работы PyomoEverest (https://github.com/
несколькими приложениями (или одним приложе- distcomp/pyomo-everest) аналогичен разработанной
нием, но на разных вычислительных ресурсах, под- нами ранее системе AMPLx [5].
ключенных к этому приложению). Задача баланси- PyomoEverest состоит из двух элементов: 1) моду-
ровки вычислительной нагрузки между ресурсами ля на Python, который посредством Everest Python
возложена на сервер Everest. API обеспечивает взаимодействие с сервисом solve-
106
ampl-stub (см. выше); 2) пула вычислительных ре- ускорена за счет одновременного решения указан-
сурсов, подключенных к сервису solve-ampl-stub ных задач пулом решателей, установленных в рас-
посредством агентов доступа Everest. пределенной вычислительной среде.
Эта вычислительная схема характеризуется пе-
риодическим решением относительно небольшого
набора (десятки) независимых и относительно
сложных задач математического программирования
(время решения каждой современным решателем на
современном сервере – не менее пары минут). Для
ее реализации в режиме распределенных вычисле-
ний предлагается использовать систему Pyomo-
Everest. Решение всех задач выполняется сервисом
оптимизации Everest, причем независимые подзада-
чи могут решаться параллельно под управлением
службы балансировки вычислительной нагрузки
Everest на ресурсы, подключенных к сервису.
Предложенный подход указывает перспективное
Рисунок 4 Модификация Python/Pyomo кода по направление применения распределенных систем на
«шаблону» PyomoEverest основе высокоуровневых средств оптимизационного
моделирования.
Типовой прием «распараллеливания» фрагмента
Благодарности
алгоритма расчетов, записанного на Pyomo, пред-
ставлен на Рис. 4. Вверху, в рамке, приведены фраг- Работа поддержана Российским научным фондом
менты кода для выбора решателя и цикла for, где (грант № 16-11-10352).
решается набор независимых подзадач. Внизу нахо-
дится фрагмент модифицированного кода. Правило Литература
модификации в том, чтобы каждый цикл for или [1] Fourer, R., Gay, D. M., Kernighan, B. W.: AMPL:
while нужно заменить тремя группами операторов: A Modeling Language for Mathematical Pro-
1. цикл формирования набора подзадач в форме gramming, 2nd edition.: Duxbury Press (2002)
AMPL-стабов; [2] Hart, W. E., Laird, C., Watson, J.-P., Wood-
2. параллельное решение задач, представленных ruff, D. L.: Pyomo-optimization modeling in py-
своими стаб-файлами, приложением Everest thon, 67, 238 p. Springer (2012)
(сейчас – solve-ampl-stub) с подключенным к [3] Hastie, T. J., Tibshirani, R. J.: Generalized additive
нему пулом AMPL-совместимых решателей; models, 43, CRC Press (1990)
3. цикл обработки результатов, доставленных в [4] Hastie, T. J., Tibshirani, R. J., Friedman, J.: Unsu-
виде набора файлов *.sol с решениями подзадач. pervised Learning. The Elements of Statistical
Модифицированный фрагмент приведен в нижней Learning: Springer, pp. 485-585 (2009)
рамке. Здесь важно использование Python-класса [5] Smirnov, S., Voloshinov, V., Sukhosroslov, O.:
pyomo.core.base.SymbolMap. Экземпляр этого класса Distributed Optimization on the Base of AMPL
(структура symbol_map) создается при записи стаб- Modeling Language and Everest Platform. Proce-
файла подзадачи в первом цикле, сохраняется (здесь dia Computer Science, 101, pp. 313-322 (2016)
– в массиве _symbMaps) и применяется во втором [6] Sukhoroslov, O., Rubtsov, A., Volkov, S.: Devel-
цикле при чтении решений из *.sol-файлов для кор- opment of Distributed Computing Applications
ректного сопоставления их содержимого структуре and Services with Everest Cloud Platform. Com-
соответствующей оптимизационной подзадачи. puter Research and Modeling, 7 (3), pp. 593-599
Заключение (2015)
[7] Sukhoroslov, O., Volkov, S., Afanasiev, A. A.:
Приведенные результаты моделирования дина- Web-Based Platform for Publication and Distrib-
мики температуры плазмы подтверждают эффек- uted Execution of Computing Applications. Paral-
тивность предложенного метода «гладкой» регуля- lel and Distributed Computing, 14th Int. Symposi-
ризации для обработки экспериментальных данных. um on IEEE, pp. 175-184 (2015)
Практическое применение метода основано на [8] Tamura, N., et al.: Impact of Nonlocal Electron
«перекрестной» (взаимной) верификации – проце- Heat Transport on the High Temperature Plasmas
дуре определения «пропущенной» части экспери- of LHD. Nuclear Fusion, 47 (5), pp. 449 (2007)
ментальных данных по остальным измерениям. Это [9] Линник, В. Г., Соколов, А. В., Миронен-
известный байесовский подход к «настройке» пара- ко, И. В.: Паттерны 137CS и их трансформация
метров некоторой регрессии по статистическим в ландшафтах ополья Брянской области. Со-
данным [3, 4]. Поскольку здесь требуется решать временные тенденции развития биогеохимии.
наборы независимых задач математического про- М.: ГЕОХИ РАН, c. 423-434 (2016)
граммирования, то работа алгоритма может быть
107
[10] Морозов, В. А.: Регулярные методы решения онной сплайн-аппроксимации. Новосибирск:
некорректно поставленных задач. М.: Наука Изд-во ИВМиМГ СО РАН (2005)
(1987) [13] Тихонов, А. И.: О математических методах ав-
[11] Пшеничный, Б. Н.: Метод линеаризации. М.: томатизации обработки наблюдений. Пробле-
Наука (1983) мы вычислительной математики. М.: МГУ, c. 3-
[12] Роженко, А. И.: Теория и алгоритмы вариаци- 17 (1980)
108