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

Задача

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

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

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

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

Двойник реализован не как скрипт, запускающий 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.150
MG-LA1.QLD1: QUAD, L=0.115, K1=8.67
MG-LA1.D400: DRIFT, L=0.400
BI-LA1.PK3:  MONI

Секции состоят из квадруполей, корректоров, пикапов и дрейфов. Импульс, набираемый пучком, растёт от 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 Тл/м, такая точность является ожидаемой, и на этом уровне сценарии смены уже можно проверять.

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

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

Запись в канал не сохраняет переданное значение, а пересчитывает физику, заложенную в структуру.

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 извлекает рассчитанные величины из файлов результатов, когда клиент читает канал-монитор.

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

# положение пучка на третьем пикапе, миллиметры
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 90. Веб-панель занимает порты 8501, 8502 и далее подряд, а порты Channel Access разнесены с шагом в тысячу. Клиенту, которому необходимы все девяносто экземпляров одновременно, передаётся список адресов в переменной 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 параметров, девятнадцать сдвигов и углов рассогласования и семь токов квадруполей, добиваясь совпадения отклика двойника с измеренным. Невязка вычисляется как средняя абсолютная разность, нормированная на шум пикапов, а испытание, прерванное ошибкой связи, получает штраф вместо остановки всего перебора. Заведомо неудачные варианты отсекает 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, оптимизатор, с помощью которого двойник подгоняется под машину.