Цифровой двойник ускорителя

Задача

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

Цифровой двойник представляет собой программу, неотличимую от настоящей машины для всех остальных программ. Для проекта Accumulator требуется ферма двойников, то есть десятки виртуальных ускорителей, запущенных одновременно и отвечающих на те же команды, что и оборудование. На такой ферме отрабатываются алгоритмы AI/ML и параллельная оптимизация без какого-либо риска и в обстановке, приближенной к реальной смене.

Ферма собрана из четырёх частей. Физическим ядром служит Elegant, расчётный код динамики пучка. Оболочкой является SCAUT, библиотека, оркеструющая эксперименты на линаке инжектора СКИФ; её слои разобраны в главе «От скрипта к приложению». Внешним интерфейсом выступает EPICS IOC, а тиражирование собранной системы осуществляется контейнерами Docker.

Слайды к главе. Устройство двойника ACCUMULATOR, граница перевода единиц, сканирование с библиотекой SCAUT и ферма двойников изложены также в четвёртой части лекции «Разбор реальных кодов» с демонстрациями в терминале; слайды лекции доступны на сайте книги и в PDF.

Устройство одного двойника

Двойник реализован не как скрипт, запускающий Elegant, а как сервер, работающий по протоколу EPICS Channel Access с теми же именами каналов и в тех же единицах, что и реальная машина.

        ┌──────────────────────────────┐
        │  Клиент (Python / CS-Studio) │◄────────────┐
        └───────────────┬──────────────┘             │
                        │ caput / caget              │ readback PV
                        ▼                            │
        ┌──────────────────────────────┐             │
        │     EPICS Channel Access     │             │
        └───────────────┬──────────────┘             │
                        │ PV: ACC1:MG-LA1:QLF1:K1    │
                        ▼                            │
        ┌──────────────────────────────┐             │
        │      ElegantIOC (pcaspy)     │─────────────┘
        └───────────────┬──────────────┘◄────────────┐
                        │ eleput / eleget            │ результат
                        ▼                            │
        ┌──────────────────────────────┐             │
        │   Elegant (физическое ядро)  │─────────────┘
        └──────────────────────────────┘

Любой клиент, написанный для настоящей машины, работает с двойником без единой правки. Скрипт коррекции, операторский экран и оптимизатор вызывают caget и caput и не различают, с чем взаимодействуют: с оборудованием или с контейнером, запущенным на ноутбуке.

Сервер построен на библиотеке pcaspy. Точка входа accumulator/main.py разбирает файл структуры и передаёт полученный словарь каналов объекту SimpleServer вместе с префиксом экземпляра, взятым из переменной окружения. Затем запускается цикл, в котором запросы, принимаемые по сети, обрабатываются с частотой, заданной в IOC_REFRESH_RATE. Вся физика вынесена в драйвер ElegantIOC, унаследованный от pcaspy.Driver: он содержит два метода, write и read, вызываемые сервером при каждом обращении к каналу.

Структура вместо списка каналов

Список каналов нигде не задаётся вручную, а выводится из описания магнитной структуры. Файл config.lte описывает линак пятью секциями, первая из которых приведена ниже.

MG-LA1.QLF1: QUAD, L=0.115, K1=-10.115
MG-LA1.D250: DRIFT, L=0.250
MG-LA1.CL2:  KICKER
MG-LA1.D150: DRIFT, L=0.100
MG-LA1.QLD1: QUAD, L=0.115, K1=8.67
MG-LA1.D400: DRIFT, L=0.400
BI-LA1.PK3:  MONI
MG-LA1.D100: DRIFT, L=0.100
AS-LA1: LINE=(MG-LA1.QLF1, MG-LA1.D250, MG-LA1.CL2, MG-LA1.D150, MG-LA1.QLD1, MG-LA1.D400, BI-LA1.PK3, MG-LA1.D100)

Секции состоят из квадруполей, корректоров, пикапов и дрейфов; длину дрейфа задаёт параметр L, а не имя, и у MG-LA1.D150 в первой секции она равна 0,100 м. Импульс, набираемый пучком, растёт от 38 МэВ/c на первой секции до 182 МэВ/c на пятой, а ускорение обеспечивают четыре резонатора с напряжением 36 МВ и частотой 2856 МГц, записанные тем же синтаксисом. Дрейфы и линии, перечисленные в файле, каналов не порождают: канал создаётся только для элемента, тип которого найден в таблице обработчиков.

Рядом расположена таблица: для каждого типа элемента перечислены параметры, выставляемые наружу, и отмечено, какие из них доступны для записи, а какие только для чтения.

ELEMENT_HANDLERS = {
    'QUAD': {
        'K1':    {'pv_suffix': 'K1',    'is_read_only': False, 'precision': 6, 'scale': 1.0},
        'DX':    {'pv_suffix': 'DX',    'is_read_only': False, 'precision': 6, 'scale': PK_MONITOR_SCALE_FACTOR},
        'betax': {'pv_suffix': 'betax', 'is_read_only': True,  'precision': 6, 'scale': 1.0},
        'betay': {'pv_suffix': 'betay', 'is_read_only': True,  'precision': 6, 'scale': 1.0},
        ...
    },
    'MONI': {
        'Cx': {'pv_suffix': 'Cx', 'is_read_only': True, 'precision': 6, 'scale': PK_MONITOR_SCALE_FACTOR},
        'Cy': {'pv_suffix': 'Cy', 'is_read_only': True, 'precision': 6, 'scale': PK_MONITOR_SCALE_FACTOR},
        ...
    },
}

Построчный разбор для каждого распознанного элемента порождает каналы в соответствии с этой таблицей.

etype = tokens[0].upper()
if etype in ELEMENT_HANDLERS:
    for param_name, field_def in ELEMENT_HANDLERS[etype].items():
        pv_name = get_pv_name(name, field_def['pv_suffix'])
        initial_val_ele = parse_param_value(rest, param_name)
        ...
        pvdb[pv_name] = {'type': 'float', 'prec': field_def['precision'],
                         'value': initial_val_ele * scale}
        element_map[pv_name] = {'name': name, 'type': etype,
                                'param': param_name,
                                'read_only': field_def['is_read_only'],
                                'scale': scale}

Для линака из пяти секций получается 116 каналов: семь квадруполей по восемь каналов, пять корректоров по четыре, пять пикапов по четыре, четыре резонатора по два и двенадцать каналов на начальные условия пучка. При добавлении в структуру ещё одного квадруполя восемь новых каналов появляются автоматически, без правки кода. Таким образом, таблица описывает предметную область, а программа её интерпретирует.

Граница, на которой переводятся единицы

Физика и система управления оперируют разными величинами. Внутри расчёта квадруполь описан так, как принято в физике пучков: K1 означает градиент, нормированный на жёсткость. Величиной же, задаваемой с пульта, является ток источника питания в амперах. Перевод между ними осуществляется по формуле

$$ k_1 = 0{,}2998,\frac{G\ [\text{Тл/м}]}{p\ [\text{ГэВ}/c]}, $$

а коэффициент записан в конфигурации вместе с его выводом.

ENERGY_LA1 = 38e-3  # [GeV/c]
# 6 A ~ 6.3 [T/m] & k1 = 0.2998 * G [T/m] / P [GeV/c]
QL_QUAD_SCALE_FACTOR = 1 / 6 * 6.3 * 0.2998
PK_MONITOR_SCALE_FACTOR = 1e3  # m to mm

PV_MAPPING = {
    "MG-LA1.QLF1.K1": "MG-LA1:QLF1-Iadd:Set",
    "BI-LA1.PK3.Cx":  "BI-LA1:PK3-fastAvX:Mea",
}
SCALE_MAPPING = {
    "MG-LA1.QLF1.K1": ENERGY_LA1 / QL_QUAD_SCALE_FACTOR,
}

При подстановке чисел коэффициент оказывается равным 0,1207, и значение \( K_1 = -10{,}115 \) соответствует −1,22 А. В скане, записанном на работающей машине, тот же канал держит −1,3 А, то есть двойник и оборудование расходятся на 6 %. Для пересчёта «ток — градиент», построенного по одной опорной точке 6 А ~ 6,3 Тл/м, такая точность является ожидаемой, и на этом уровне сценарии смены уже можно проверять. Для остальных секций в той же конфигурации записаны импульсы 76, 114, 152 и 190 МэВ/c, то есть с шагом 38 МэВ/c, тогда как четыре резонатора по 36 МВ дают 74, 110, 146 и 182 МэВ/c; масштабы секций со второй по пятую расходятся с импульсом структуры на 2,7–4,4 %, и одно из двух значений необходимо привести к другому.

Здесь же решается вторая задача — задача имён. По умолчанию канал называется MG-LA1:QLF1:K1, однако таблица PV_MAPPING переименовывает его в MG-LA1:QLF1-Iadd:Set, то есть так, как этот источник питания назван на пульте. Пикап, измеряющий горизонтальное положение, отвечает на имя BI-LA1:PK3-fastAvX:Mea. Именно эти имена ожидает помощник оператора из следующей главы.

Обработка одного caput

Запись в канал не сохраняет переданное значение, а пересчитывает физику, заложенную в структуру. Метод write приведён ниже в сокращении: проверка того, что канал есть в element_map, запись в журнал перед eleput и обработка исключения опущены.

def write(self, reason, value):
    element_info = self.element_map[reason]
    if element_info['read_only']:
        logger.warning(f"Attempt to write to read-only monitor {reason}")
        return False
    sim_name = get_sim_name(element_info['name'], element_info['param'])
    scale = element_info.get('scale', 1.0)
    val_to_write = value / scale
    eleput(sim_name, val_to_write)
    self.setParam(reason, value)
    return True

Сервер находит элемент, указанный в имени канала, проверяет право на запись и делит значение на масштаб, переводя амперы обратно в \( K_1 \). Далее функция eleput, взятая из SCAUT, записывает параметр и заново запускает Elegant, трассирующий пучок через всю структуру. Один вызов caput влечёт полный пересчёт: за ответом стоит решение уравнений движения, а не заранее заготовленная таблица. Обратная функция eleget извлекает рассчитанные величины из файлов результатов, когда клиент читает канал-монитор.

Опущенная обработка исключения существенна. Если eleput завершается ошибкой, метод write записывает её в журнал и возвращает False, то есть запись отклоняется. Чтение устроено иначе: при исключении в eleget метод read записывает ошибку в журнал и возвращает прежнее значение канала, не выставляя ни статуса, ни тревоги, и клиент получает вместо отказа устаревшее число. В EPICS для такого случая предусмотрены статус и тревога канала, которые позволяют клиенту отличить устаревшее значение от рассчитанного, и двойник, неотличимый от машины, должен их выставлять.

Проверка описанного механизма выполняется двумя командами.

# положение пучка на третьем пикапе, миллиметры
caget ACC1:BI-LA1:PK3-fastAvX:Mea

# добавка к току первого квадруполя, амперы
caput ACC1:MG-LA1:QLF1-Iadd:Set -1.3

Ферма

Подгонка двойника представляет собой задачу оптимизации с дорого вычисляемой целевой функцией. Один шаг требует десятков запусков Elegant, и тысяча испытаний, выполненных последовательно, заняла бы сутки. Распараллелить сам Elegant внутри процесса невозможно: он однопоточный и обменивается данными через файлы, размещённые в рабочем каталоге. Однако можно запустить множество его копий, по контейнеру на экземпляр, каждый со своей файловой системой и своим префиксом каналов.

                     Локальная сеть
   ┌────────────────────────────────────────────┐
   │  Компьютер 1                               │
   │     ACC1   ACC2   ...   ACC30              │◄──────┐
   │                                            │       │
   │  Компьютер 2                               │       │   ┌───────────────┐
   │     ACC31  ACC32  ...   ACC60              │◄──────┼───│  Оркестратор  │
   │                                            │       │   │  Optuna / ГА  │
   │  Компьютер 3                               │       │   └───────────────┘
   │     ACC61  ACC62  ...   ACC90              │◄──────┘
   └────────────────────────────────────────────┘

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

accelerator-1:
  image: accumulator:latest
  container_name: acc-emulator-1
  env_file: accumulator.env
  environment:
    - 'ACCUMULATOR_PREFIX=ACC1:'
  ports: ['8501:8501', '5064:5064', '5065:5065/udp']

Описание из девяноста одинаковых блоков вручную не составляется. Скрипт generate_compose.py, вызванный с числом экземпляров и номером первого из них, распределяет порты с заданным шагом и записывает готовый файл для одного компьютера: python scripts/generate_compose.py --count 30 --start 31 для второго, и так далее. Веб-панель занимает порты 8501, 8502 и далее подряд, а порты Channel Access разнесены с шагом в тысячу, и отсчёт портов в каждом файле начинается заново. С --count 90 в одном файле порты Channel Access экземпляров с 62-го по 90-й вышли бы за 65535 (у девяностого 94064 и 94065), а проверки диапазона генератор не делает. Клиенту, которому необходимы все девяносто экземпляров одновременно, передаётся список адресов в переменной EPICS_CA_NAME_SERVERS: поиск канала широковещательной рассылкой заменяется прямым опросом перечисленных серверов.

Внутри каждого контейнера supervisor поддерживает два процесса, IOC и веб-панель, а образ содержит Elegant 2025.2.0 и набор утилит SDDS, установленных поверх Ubuntu 22.04. Панель, реализованная на Streamlit, читает тот же файл структуры, отображает каналы, сгруппированные по элементам, и позволяет править уставки из графического интерфейса без кода на стороне клиента.

Оркестратор распределяет испытания по свободным экземплярам через очередь.

INSTANCE_POOL = queue.Queue()
for _p in PREFIXES:            # ACC1:, ACC2:, ... ACC89:
    INSTANCE_POOL.put(_p)

def objective(trial: optuna.Trial) -> float:
    prefix = INSTANCE_POOL.get()      # ждём свободный экземпляр
    try:
        values = [trial.suggest_float(p, lo, hi)
                  for p, (lo, hi) in zip(actuator_pvs, bounds)]
        apply_actuators(prefix, actuator_pvs, values, strip=REAL_PREFIX)
        ...
        return float(np.mean(step_losses))
    finally:
        INSTANCE_POOL.put(prefix)     # возвращаем при любом исходе

Очередь здесь передаёт не данные, а право пользоваться ресурсом. Восемьдесят девять префиксов помещаются в пул, а девяностый экземпляр остаётся нетронутым: по нему впоследствии сверяется состояние двойника до подгонки. Число параллельных заданий равно размеру пула, а перебор ограничен тысячей испытаний и тысячей секунд.

Тысяча вариантов коррекции оптики: по осям — добавки к токам квадруполей

Цикл, ради которого всё затевалось

Двойник является звеном цикла, превращающего снятые измерения в уставки. Тетрадь optuna_tuner.ipynb проходит этот цикл целиком.

  1. Скан на реальной машине. Десять корректоров проходят заданную сетку, семнадцать каналов записываются как отклик, и результат сохраняется в JSON, снятый за 28 шагов.
  2. Подгонка двойника. Optuna перебирает 26 параметров, девятнадцать сдвигов и углов рассогласования и семь токов квадруполей, добиваясь совпадения отклика двойника с измеренным. Невязка вычисляется как средняя абсолютная разность, нормированная на шум пикапов, а испытание, прерванное ошибкой связи, получает штраф вместо остановки всего перебора; ошибку eleget на стороне сервера, как сказано выше, этот штраф не обнаруживает. Заведомо неудачные варианты отсекает MedianPruner, и на каждое испытание берётся десять шагов скана из двадцати восьми.
  3. Коррекция на двойнике. Подогнанной модели в качестве цели задаются бета-функции номинального расчёта, четырнадцать чисел на семи квадруполях, и подбираются токи в пределах ±0,2 А от текущих.
  4. Перенос на оборудование. Функция apply_to_accelerator сначала читает текущие значения, сохраняет прочитанное в файл снимка, печатает таблицу «было, станет, разница» и только после этого, и только при явно переданном dry_run=False, записывает значения в машину.
  5. Проверка. Новый скан подтверждает улучшение.

После второго шага двойник моделирует не «линак вообще», а конкретную машину в конкретный день, со всеми рассогласованиями, зафиксированными в скане. В этом заключается разница между моделью и цифровым двойником.

Настроенный экземпляр получает отдельный префикс ACCT:, номинальный расчёт размещён под префиксом ACCM:, а за ACC: находится настоящая машина. Бета-функции, снятые со всех трёх, сводятся в одну таблицу, в которой видна невязка как до подгонки, так и после неё.

Назначение фермы

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

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

Оптимизация бета-функции на виртуальном ускорителе

Обучение ML/AI-моделей. Девяносто экземпляров порождают размеченные данные в объёмах, недостижимых на живой установке: каждая пара «уставки — отклик» обходится в один прогон Elegant, а не в часть смены. Выборка, набранная на ферме, покрывает и те режимы, до которых на реальной машине очередь не доходит. На этих данных обучаются модели коррекции орбиты и помощник оператора Сава, которому посвящена следующая глава.

Коррекция орбиты моделью машинного обучения

Полезные ссылки

  • Accumulator, исходники двойника, разобранного в этой главе.
  • SCAUT, библиотека, оркеструющая эксперименты на ускорителе.
  • Elegant, код расчёта динамики пучка, положенный в основу двойника.
  • EPICS, система управления, протокол которой воспроизводит двойник.
  • pcaspy, библиотека, позволяющая написать сервер EPICS на Python.
  • Optuna, оптимизатор, с помощью которого двойник подгоняется под машину.