Коррекция равновесной орбиты с применением матрицы отклика и нейронных сетей

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

Постановка задачи

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

Для устранения этих искажений в кольце установлены два набора устройств.

  • Мониторы положения пучка (BPM, пикапы) — \(N_{bpm}\) штук, измеряющие отклонение орбиты \(x_i\) в своих точках;
  • Дипольные корректоры — \(N_{cor}\) штук, слабые магниты, каждый из которых поворачивает траекторию на угол \(\theta_j\).

Задача коррекции, стоящая перед оператором каждую смену, состоит в следующем: по измеренному вектору отклонений \(\vec{x}\) найти такие углы корректоров \(\vec{\theta}\), при которых орбита окажется как можно ближе к расчётной.

Матрица отклика

Пока отклонения малы, ускоритель остаётся линейной системой. Реакция орбиты на \(j\)-й корректор, включённый в одиночку, описывается матрицей отклика (response matrix).

$$ x_i = \sum_{j} R_{ij},\theta_j, \qquad R_{ij} = \frac{\partial x_i}{\partial \theta_j}. $$

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

$$ R_{ij} = \frac{\sqrt{\beta_i \beta_j}}{2\sin(\pi\nu)}\cos\left(|\varphi_i - \varphi_j| - \pi\nu\right), $$

где \(\beta_i\), \(\beta_j\) — значения бета-функции в точках монитора и корректора, \(\varphi_i\), \(\varphi_j\) — набеги бетатронной фазы, \(\nu\) — бетатронная частота (число колебаний на оборот). Матрицу можно рассчитать по модели оптики, а можно измерить, отклоняя каждый корректор по очереди и записывая отклик всех мониторов. Обычно поступают именно так, одновременно проверяя модель.

Построим модельное кольцо на Python, устроенное по этой же формуле.

import numpy as np

rng = np.random.default_rng(42)

n_bpm, n_cor = 64, 48       # мониторы и корректоры
nu = 8.42                   # бетатронная частота кольца

# расставим элементы по кольцу и разыграем бета-функции
phi_bpm = np.sort(rng.uniform(0, 2 * np.pi * nu, n_bpm))   # фазы мониторов
phi_cor = np.sort(rng.uniform(0, 2 * np.pi * nu, n_cor))   # фазы корректоров
beta_bpm = rng.uniform(4.0, 24.0, n_bpm)                   # бета в мониторах, м
beta_cor = rng.uniform(4.0, 24.0, n_cor)                   # бета в корректорах, м

def response_matrix(beta_b, phi_b, beta_c, phi_c, nu):
    """Матрица отклика замкнутой орбиты кольца, мм/мрад."""
    dphi = np.abs(phi_b[:, None] - phi_c[None, :])
    return (np.sqrt(beta_b[:, None] * beta_c[None, :])
            / (2 * np.sin(np.pi * nu))
            * np.cos(dphi - np.pi * nu))

R = response_matrix(beta_bpm, phi_bpm, beta_cor, phi_cor, nu)
print(R.shape)   # (64, 48)

Внесём в машину нормально распределённые случайные ошибки и рассмотрим полученную орбиту.

# истинные ошибки — как будто от смещений квадруполей
theta_err = rng.normal(0, 0.05, n_cor)          # мрад
noise = rng.normal(0, 0.02, n_bpm)              # шум мониторов, мм
x_measured = R @ theta_err + noise              # измеренная орбита

print(f"СКО орбиты до коррекции: {x_measured.std():.3f} мм")

СКО орбиты до коррекции: 1.539 мм

Коррекция через SVD

Нахождение токов корректоров сводится к решению обратной задачи \(R\vec{\theta} \approx -\vec{x}\). Матрица прямоугольная (\(64 \times 48\)), система переопределена, поэтому решение ищется в смысле наименьших квадратов. Стандартным инструментом в физике ускорителей является сингулярное разложение (SVD).

$$ R = U,\Sigma,V^T, \qquad \vec{\theta} = -V,\Sigma^{-1}U^T \vec{x}. $$

Часть сингулярных чисел \(\sigma_k\) близка к нулю. Соответствующие комбинации корректоров почти не влияют на орбиту в мониторах, и при обращении \(1/\sigma_k\) для таких мод неограниченно растёт, вследствие чего шум мониторов превращается в огромные токи корректоров. Это классическая некорректная обратная задача, и решается она отсечением малых сингулярных чисел.

def svd_correction(R, x, n_modes):
    """Коррекция орбиты с обрезанием SVD-спектра до n_modes мод."""
    U, s, Vt = np.linalg.svd(R, full_matrices=False)
    s_inv = np.zeros_like(s)
    s_inv[:n_modes] = 1.0 / s[:n_modes]     # оставляем только сильные моды
    return -(Vt.T * s_inv) @ (U.T @ x)

for n_modes in [10, 30, 48]:
    theta_corr = svd_correction(R, x_measured, n_modes)
    x_after = R @ (theta_err + theta_corr) + rng.normal(0, 0.02, n_bpm)
    print(f"мод: {n_modes:2d} | СКО орбиты: {x_after.std():6.3f} мм"
          f" | максимальный угол: {np.abs(theta_corr).max():6.3f} мрад")

мод: 10 | СКО орбиты:  0.297 мм | максимальный угол:  0.061 мрад
мод: 30 | СКО орбиты:  0.028 мм | максимальный угол:  0.126 мрад
мод: 48 | СКО орбиты:  0.030 мм | максимальный угол: 1343310628796.820 мрад

Десять мод уменьшают орбиту в пять раз, тридцать — в пятьдесят пять, причём токи корректоров остаются приемлемыми. Все сорок восемь мод орбиту не улучшают (0.030 против 0.028 мм), однако требуемые углы возрастают до \(10^{12}\) мрад.

Происхождение величины \(10^{12}\) видно из спектра. Наибольшее сингулярное число здесь равно 183.6, сорок первое — 0.057, а последние семь лежат в районе \(10^{-14}\), то есть на уровне ошибки округления. Матрица модели вырождена: её ранг равен 41, а не 48. Семь комбинаций корректоров не влияют на орбиту, и деление на малое число оказывается делением на ноль, размытый арифметикой с плавающей точкой.

На реальном кольце измеренная матрица отклика обусловлена плохо, но не вырождена: отношение наибольшего сингулярного числа к наименьшему там составляет \(10^2 - 10^3\), а не \(10^{16}\), как здесь. Усекать спектр всё равно приходится, однако по другой причине.

Число удерживаемых мод является компромиссом между качеством коррекции и величиной токов. Родственный подход даёт регуляризация Тихонова, где вместо жёсткого усечения минимизируется \(|R\theta + x|^2 + \lambda|\theta|^2\).

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

Причём здесь нейронные сети

Линейный подход имеет границы применимости.

  • при больших искажениях и в сильно нелинейных элементах (секступолях) отклик, предсказанный линейной теорией, расходится с действительным;
  • реальная оптика дрейфует, и матрица, измеренная утром, к вечеру уже несколько отличается;
  • измерение полной матрицы отклика является длительной процедурой, отнимающей пучковое время.

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

Построим такой корректор на многослойном перцептроне.

from sklearn.neural_network import MLPRegressor
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline

# обучающая выборка: случайные ошибки -> орбиты (с шумом мониторов)
n_samples = 20_000
thetas = rng.normal(0, 0.05, (n_samples, n_cor))
orbits = thetas @ R.T + rng.normal(0, 0.02, (n_samples, n_bpm))

model = make_pipeline(
    StandardScaler(),
    MLPRegressor(hidden_layer_sizes=(128, 128), activation="relu",
                 max_iter=200, early_stopping=True, random_state=0),
)
model.fit(orbits, thetas)                      # учимся восстанавливать ошибки

theta_pred = model.predict(x_measured[None, :])[0]
x_after_nn = R @ (theta_err - theta_pred) + rng.normal(0, 0.02, n_bpm)
print(f"СКО орбиты после нейросетевой коррекции: {x_after_nn.std():.3f} мм")

СКО орбиты после нейросетевой коррекции: 0.151 мм

Значение 0.151 мм в десять раз лучше, чем до коррекции, но в пять раз хуже, чем даёт SVD с тридцатью модами. Это ожидаемо. Модельная машина линейна по построению, а для линейной задачи SVD является точным решением в смысле наименьших квадратов; превзойти его невозможно, и сеть, обученная на тех же данных, может лишь приблизиться к нему, затратив двадцать тысяч примеров и минуту обучения.

Такой тест обязателен. На линейной задаче нейросеть не нужна. Если бы сеть превзошла SVD, это означало бы ошибку в постановке эксперимента. Машинное обучение необходимо там, где линейная модель перестаёт быть точной.

Гибридный корректор. Сеть учится на остатках

Реальная машина линейна лишь приближённо. Корректоры насыщаются, секступоли вносят нелинейный отклик, железо магнитов имеет гистерезис, и после SVD-коррекции в орбите остаётся систематический остаток, не описываемый линейной моделью. Именно он передаётся сети.

Такая конструкция используется на действующих установках; проверим её на модели.

$$ \Delta\theta = \underbrace{\Delta\theta_{\text{SVD}}}{\text{линейная часть}} + \underbrace{\Delta\theta{\text{НС}}}_{\text{поправка на нелинейность}}. $$

Сеть обучается не отображению «орбита → токи» целиком, а только разности между истинной ошибкой и ответом, предложенным SVD. Задача на порядок проще: вместо всей линейной алгебры кольца сеть уточняет почти правильный ответ линейной части. Сети достаточно двух скрытых слоёв 64–32.

Усилим нелинейность: примем предел корректора в несколько раз меньше типичного требуемого угла (ошибки разыгрываются с \(\sigma = 0.05\) мрад), чтобы линейная модель заведомо перестала быть точной. Кольцо, матрица отклика R и число удерживаемых мод остаются прежними, заданными в начале главы.

THETA_MAX = 0.01          # предел корректора, мрад
sigma_bpm = 0.02          # шум мониторов, мм
n_modes = 30              # удерживаемых SVD-мод, как и раньше

def machine(theta):
    # кольцо с насыщающимися корректорами: отклик уже не линеен
    return R @ (THETA_MAX * np.tanh(theta / THETA_MAX))

Псевдообратную матрицу из усечённого спектра рассчитаем один раз заранее: коррекция сводится к умножению на неё, и разложение \(R\) на каждом шаге не требуется.

U, s, Vt = np.linalg.svd(R, full_matrices=False)
s_inv = np.zeros_like(s)
s_inv[:n_modes] = 1.0 / s[:n_modes]
svd_matrix = -(Vt.T * s_inv) @ U.T        # то же, что svd_correction, но матрицей

Обучающая выборка строится так же, как раньше, однако целью служит остаток, не устранённый SVD.

n_samples = 20_000
thetas = rng.normal(0, 0.05, (n_samples, n_cor))
orbits = (THETA_MAX * np.tanh(thetas / THETA_MAX)) @ R.T \
         + rng.normal(0, sigma_bpm, (n_samples, n_bpm))

svd_answer = orbits @ svd_matrix.T        # что предложил бы SVD
residuals = -thetas - svd_answer          # чего ему не хватило

Сеть обучается предсказывать этот остаток по орбите, а шаг коррекции суммирует оба вклада.

net = make_pipeline(
    StandardScaler(),
    MLPRegressor(hidden_layer_sizes=(64, 32), max_iter=300,
                 early_stopping=True, random_state=0),
)
net.fit(orbits, residuals)

def hybrid_step(x):
    # один шаг коррекции: линейное решение плюс поправка сети
    return svd_matrix @ x + net.predict(x[None, :])[0]

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

Проведём 100 испытаний со случайными ошибками, генерируемыми заново для каждого, и сравним чистый итеративный SVD с гибридом. В таблице приведены среднее по испытаниям и среднеквадратичный разброс.

def run(step, n_iter=6, n_trials=100):
    out = np.zeros((n_trials, n_iter + 1))
    for t in range(n_trials):
        theta = rng.normal(0, 0.05, n_cor)
        for k in range(n_iter + 1):
            x = machine(theta) + rng.normal(0, sigma_bpm, n_bpm)
            out[t, k] = x.std()
            if k < n_iter:
                theta = theta + step(x)
    return out.mean(0), out.std(0)

svd_mean, svd_std = run(lambda x: svd_matrix @ x)   # чистый итеративный SVD
nn_mean, nn_std = run(hybrid_step)                  # гибрид SVD плюс сеть

print(f"до коррекции {svd_mean[0]:.3f} ± {svd_std[0]:.3f} мм")
print()
print("итер |             SVD |   гибрид SVD+НС | выигрыш")
for k in range(1, len(svd_mean)):
    gain = 100 * (svd_mean[k] - nn_mean[k]) / svd_mean[k]
    print(f"{k:4d} | {svd_mean[k]:7.3f} ± {svd_std[k]:5.3f} |"
          f" {nn_mean[k]:7.3f} ± {nn_std[k]:5.3f} | {gain:+5.1f}%")

до коррекции 0.306 ± 0.099 мм

итер |             SVD |   гибрид SVD+НС | выигрыш
   1 |   0.263 ± 0.093 |   0.191 ± 0.063 | +27.1%
   2 |   0.226 ± 0.079 |   0.198 ± 0.068 | +12.2%
   3 |   0.196 ± 0.063 |   0.214 ± 0.078 |  -9.1%
   4 |   0.160 ± 0.047 |   0.223 ± 0.081 | -38.9%
   5 |   0.138 ± 0.045 |   0.232 ± 0.084 | -68.0%
   6 |   0.117 ± 0.039 |   0.242 ± 0.077 | -107.1%

Из таблицы следуют три вывода.

Сеть выигрывает на первых итерациях. На первом шаге гибрид даёт орбиту более чем на четверть лучше, 0.191 против 0.263 мм. Это вклад нелинейности, которую SVD не учитывает, а сеть выучила. На реальной машине выигрыш измеряется не в миллиметрах, а в пучковом времени: каждая итерация коррекции представляет собой отдельный цикл «выставить токи, дождаться стабилизации, измерить орбиту», и сокращение их числа ценнее, чем выигрыш последних микрон.

Начиная с третьей итерации гибрид проигрывает, и разрыв растёт. Чистый SVD к шестой итерации доводит орбиту до 0.117 мм, а гибрид выходит на плато и медленно растёт к 0.242 мм. Плато возникает из-за систематической ошибки сети. Обученная на конечной выборке, она предсказывает остаток с погрешностью, и на каждой итерации эта погрешность вносится в орбиту. Сеть сама становится источником ошибки.

Итеративность решает ту же задачу, что и сеть, и решает её лучше. Два улучшения направлены на одну проблему, поэтому суммировать их бессмысленно. Если сравнить гибрид после трёх итераций (0.214 мм) с однократным SVD (0.263 мм), получилось бы улучшение на 19 %, которое было бы приписано сети, тогда как принадлежит оно итерациям. Вклад каждого улучшения необходимо разделять. В противном случае легко опубликовать результат, который не выдержит попытки воспроизведения.

У гибридного корректора существует оптимальное число итераций, и его необходимо измерить. Здесь оно равно единице: уже на второй итерации гибрид даёт 0.198 мм против 0.191 мм на первой, а его преимущество над чистым SVD сокращается с +27.1 % до +12.2 %, чтобы на третьей смениться проигрышем. Останавливаться следует там, где поправка, предсказанная сетью, перестаёт превышать её собственную погрешность.

Величина выигрыша зависит от кольца: на семи разыгранных наборах бета-функций и фаз первая итерация давала от 2 % до 27 % улучшения. Форма кривой — выигрыш на первых шагах, плато, проигрыш далее — повторяется всегда. Поэтому в таблице приведён разброс, а не одно среднее; без него из таких данных можно «доказать» что угодно.

Стоимость одного шага также имеет значение. Один шаг SVD сводится к умножению матрицы \(48\times64\) на вектор и занимает 1.2 мкс на Apple M4; один вызов сети, обёрнутой в scikit-learn, занимает 110 мкс, почти в сто раз дольше. Для медленной обратной связи с частотой в единицы герц это несущественно, а для быстрой (килогерцы) уже существенно, и сеть потребуется вынести из Python в скомпилированный код.

Применимость на реальной установке

На линейном ускорителе инжектора ЦКП «СКИФ» коррекция орбиты выполняется по измеренной матрице отклика, для чего служат 7 датчиков положения пучка, 6 больших и 8 малых рамочных двухкоординатных корректоров. Матрица там измеряется, а не берётся из модели. Каждый корректор по очереди отклоняется на известную величину, и записывается отклик всех пикапов. Процедура длительная, однако одновременно она проверяет модель оптики. Коррекция выполняется итеративно, с повторным измерением матрицы отклика между итерациями. Без неё транспортировка пучка по ускорителю и каналу сопровождается значительными потерями, что обнаруживается при запуске каждой новой машины.

Ниже приведены три обстоятельства, не заметные на модели, но определяющие успех на пульте.

  • Выбор числа сингулярных чисел остаётся главным настроечным параметром. При слишком большом числе в решение попадают шумовые компоненты, усиленные делением на малые сингулярные числа, и коррекция ухудшается; при слишком малом коррекция оказывается недостаточной. Число выбирается по минимуму невязки, перебором, как показано выше.
  • Ограничения на токи корректоров входят в алгоритм, а не в проверку после него. Источник питания корректора имеет предел в несколько ампер; если ограничение не заложено в постановку задачи, оптимизатор выдаёт неисполнимый ответ, подобный варианту со всеми сорока восемью модами выше.
  • Сеть, если она добавляется, инициализируется SVD-решением. Это та же работа по остаткам, что и выше на модели: линейная часть даёт разумное приближение, сеть отвечает за поправку. Декомпозиция обеспечивает и предсказуемую деградацию: если сеть выдаст некорректный результат, останется рабочее SVD-решение. Гибридная схема с ограничениями на токи и SVD-инициализацией проверялась, в частности, на модели накопителя ВЭПП-5, где стоят 16 датчиков положения и 16 корректоров с пределом \(\pm 7\) А.

Нейросетевой корректор всегда работает рядом с классическим, а не вместо него. SVD-коррекция с усечением мод является прозрачным и предсказуемым алгоритмом, сеть — непрозрачным дополнением для режимов, где линейная модель перестаёт работать. Предложенные сетью токи проходят те же ограничения и проверки, что и любые другие.

Задания для самостоятельной работы

  1. На основе приведённого выше кода построить зависимость СКО орбиты и максимального тока корректора от числа удерживаемых SVD-мод, получив классическую L-кривую регуляризации.
  2. Усилить нелинейность (снизить \(\theta_{max}\) вдвое) и повторить таблицу из 100 испытаний, определив величину насыщения, при которой преимущество гибрида перестаёт исчезать к третьей итерации.
  3. Испортить один монитор (умножить его показания на 3) и определить, какой из методов устойчивее к неисправному датчику; обучить сеть на данных, испорченных такими же неисправностями.
  4. Заложить в коррекцию предел корректора, ограничив предлагаемые углы по \(\pm\theta_{max}\) до подачи на машину, и оценить, насколько это меняет число итераций до выхода на плато.

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