ビームサイズを調整する例

機械学習環境を個人用PCにセットアップする に続いて、 より加速器に近い例を示します。評価関数をどのようにすべきか、など、色々な テストにも使えますので活用してください。なにかの参考になれば幸いです。

GPyTorch + Ocelot でビームサイズ最小化

下図に示すようにQM2台でビームサイズを調整する例を考えます。 初期条件（Twiss parameter; α,　β）は適当です。簡単のためディスパージョンηもゼロでスタートです。

実際の加速器ではQMの極性（QF/QD）が決まっていることが多いとは思いますが、 例としてバイポーラ電源がつながっていると想定します。

軌道計算

これくらいの例ならば自分で転送行列を書けば良いです。 ただし、今後の発展性を考えて、ここでは Ocelot を使います。 インストール方法はOcelotインストール方法を参照してください。

GPyOptによるビームサイズ最小化

コードをそのまま書きます。本当は適切にクラス化するべきですが、 この長さならベタ書きで分かると思うので global 変数を使ってQMのk値を関数側で変えてしまいます。

また、ビームサイズゼロは不可能なので、実際には ターゲットとなるビームサイズを定義する方が良いでしょう。

import GPy
import GPyOpt
import matplotlib.pyplot as plt
import numpy as np
import time

from ocelot import *
from ocelot.gui.accelerator import *

def func(x):
    global cnt, QM1, QM2, tws, lat, tws0
    x0=x[:,0][0]
    x1=x[:,1][0]

    # update optics
    QM1.k1 = x0
    QM2.k1 = x1
    lat.update_transfer_maps()
    tws = twiss(lat, tws0, nPoints=1000)

    idx = -1  # last point of the lattice
    sx = np.sqrt(tws[idx].beta_x)
    sy = np.sqrt(tws[idx].beta_y)

    # Evaluation Function Example
    val = np.abs(sx) + np.abs(sy)
    # val = np.abs(sx*sy)
    # val = np.log(np.abs((sx+10)*(sy+10)))
    # val = -1.0/np.abs((sx+10)*(sy+10))

    print("%d: %+8.3f,%+8.3f,%+8.3f,%+8.3f,%+8.3f"%(cnt, x0, x1, sx, sy, val))
    cnt = cnt + 1
    return val


# ==== Main =====

# Define Optics
## Drift
L20  = Drift(l=0.2, eid='L20')
L50  = Drift(l=0.5, eid='L50')
## Quadrupoles
QM1 = Quadrupole(l=0.2, k1=3.9, eid='QM1')
QM2 = Quadrupole(l=0.2, k1=-3.25, eid='QM2')
## Lattice
cell = (L50, QM1, L20, QM2, L50, L50)
lat = MagneticLattice(cell, stop=None)
print("length of the cell: ", lat.totalLen, " m")
tws0 = Twiss()
tws0.beta_x = 20.0
tws0.beta_y = 30.0
tws0.alpha_x = -10.0
tws0.alpha_y = -5.0
tws0.emit_x = 1.0
tws0.emit_y = 1.0
tws0.E = 1.0

tws = twiss(lat, tws0, nPoints=1000)

# GP Optimization
print("==== start ====")
print("# cnt, x0, x1, sx, sy, eval_val")

cnt = 0
bounds = [{'name':'x0', 'type':'continuous', 'domain':(-20, 20)},
          {'name':'x1', 'type':'continuous', 'domain':(-20, 20)}]

myBopt = GPyOpt.methods.BayesianOptimization(f=func,domain=bounds,initial_design_numdata=10,acquisition_type='LCB')

myBopt.run_optimization(max_iter=50)
myBopt.plot_acquisition()
myBopt.plot_convergence()
print("Best     = ", myBopt.x_opt)
print("eval_val = ", func(np.array([myBopt.x_opt])))

# Plot results
QM1.k1 = myBopt.x_opt[0]
QM2.k1 = myBopt.x_opt[1]
lat.update_transfer_maps()
tws = twiss(lat, tws0, nPoints=1000)

s = [p.s for p in tws]
beta_x = [p.beta_x for p in tws]
beta_y = [p.beta_y for p in tws]

plt.plot(s, beta_x, 'b-');
plt.plot(s, beta_y, 'r-');

plt.legend(['Horiz', 'Vert'])
plt.xlabel('s [m]')
plt.ylabel('beta [m]')
plt.grid(True)
plt.show()

計算結果例1 : QF+QD に落ち着いた場合

計算結果例2 : スキャンする範囲を QD+QF に限定した場合

次のステップ

適当に色々と変更して試すと良い

• ビームの初期パラメータを変える：簡単のためには円形ビーム（βx = βy ）にして、αx=αy=0 にすると、もっと単純
• 解析解との比較
• 評価関数を色々変える：sx+sy, sx*sy, log(sx+sy+), -1/(sx+sy), ....ゼロ割に気をつけて適当なオフセットを入れると-1/((sx+10)*(sy+10)) とか....
  • それぞれで最小値を見つける速度が変わる
  • 単なるビームサイズだけではなく、他の条件（個々のQM k値を出来るだけ小さくするとか、途中のビームサイズを小さくするとか(betatron関数を小さくする））などを入れてみる
• Optics的には Bendを入れてdispersionが入った状態にするとか、空間電荷効果を入れるとか、より実践的な形に
• 測定したビームサイズに意図的に「ノイズ」を入れて、GPの応答がどうなるかを確認する（重要！）
• GPのカーネルを変える
• 他のGPツール との比較
• 他の最適化手法（Nelder–Mead(downhill simplex)とか、最急降下法とか、遺伝的アルゴリズムとか）との比較、あるいは組み合わせを考える