Matching of turn pattern measurements for cyclotrons using multi-objective optimization

TL;DR

This paper introduces a multi-objective optimization approach using evolutionary algorithms to accurately match turn patterns in cyclotrons, specifically addressing uncertainties in simulation inputs and incorporating a realistic trim coil model within the OPAL framework.

Contribution

A novel multi-objective optimization method with a realistic trim coil model for better turn pattern matching in cyclotrons is presented.

Findings

Reduced maximum absolute error to 4.54 mm across 182 turns

Successfully integrated a realistic trim coil model into OPAL

Demonstrated effectiveness of evolutionary optimization in complex parameter spaces

Abstract

The usage of numerical models to study the evolution of particle beams is an essential step in the design process of particle accelerators However, uncertainties of input quantities such as beam energy and magnetic field lead to simulation results that do not fully agree with measurements, hence the final machine will behave slightly differently than the simulations In case of cyclotrons such discrepancies affect the overall turn pattern or may even alter the number of turns in the machine Inaccuracies at the PSI Ring cyclotron facility that may harm the isochronism are compensated by additional magnetic fields provided by 18 trim coils These are often absent from simulations or their implementation is very simplistic In this paper a newly developed realistic trim coil model within the particle accelerator framework OPAL is presented that was used to match the turn pattern of the PSI Ring cyclotron Due to the high-dimensional search space consisting of 48 design variables (simulation input parameters) and 182 objectives (i.e turns) simulation and measurement cannot be matched in a straightforward manner Instead, an evolutionary multi-objective optimisation with a population size of more than 8000 individuals per generation together with a local search approach were applied that reduced the maximum absolute error to 4.54 mm over all 182 turns.

Figures

Tables

Table 1. Table 1: Result of best individual obtained by optimization using the l ∞ subscript 𝑙 l_{\infty} -error norm for each objective. The label σ [ l , u ] subscript 𝜎 𝑙 𝑢 \sigma_{[l,u]} indicates an objective for the turns in the range [ l , u ] 𝑙 𝑢 [l,u] .

Objective	$l_{\infty }$-error	Probe
${\sigma }_{[l,u]}$	($\mathrm{mm}$)
${\sigma }_{[1,16]}$	$6.38$	RRI2
${\sigma }_{[9,31]}$	$3.76$	RRL
${\sigma }_{[32,61]}$	$6.34$	RRL
${\sigma }_{[62,105]}$	$4.39$	RRL
${\sigma }_{[106,148]}$	$2.91$	RRL
${\sigma }_{[149,182]}$	$3.27$	RRL

Table 2. Table 2: Maximum absolute error ( l ∞ subscript 𝑙 l_{\infty} -norm), mean absolute error (MAE) and the mean squared error (MSE) of the best individual of the optimizer and local search compared to the measurement. In both cases the maximum error is at turn 2 2 2 .

Method	$l_{\infty }$-norm	MAE	MSE
	($\mathrm{mm}$)	($\mathrm{mm}$)	($\mathrm{mm}{}^2$)
optimizer	$6.38$	$\mathrm{2.029\hspace{0.17em}052\hspace{0.17em}631\hspace{0.17em}58}$	$\mathrm{6.301\hspace{0.17em}941\hspace{0.17em}052\hspace{0.17em}63}$
local search	$4.54$	$\mathrm{1.400\hspace{0.17em}578\hspace{0.17em}947\hspace{0.17em}37}$	$\mathrm{3.404\hspace{0.17em}495\hspace{0.17em}263\hspace{0.17em}16}$

Table 3. Table 3: Maximum absolute error ( l ∞ subscript 𝑙 l_{\infty} -norm), mean absolute error (MAE) and the mean squared error (MSE) of the measurement or multi particle tracking simulation including space charge to the multi particle tracking simulation neglecting space charge.

Comparison to	$l_{\infty }$-norm	MAE	MSE
	($\mathrm{mm}$)	($\mathrm{mm}$)	($\mathrm{mm}{}^2$)
measurement	$4.64$	$\mathrm{1.462\hspace{0.17em}368\hspace{0.17em}421\hspace{0.17em}05}$	$\mathrm{3.585\hspace{0.17em}795\hspace{0.17em}263\hspace{0.17em}16}$
space charge	$\mathrm{0.050\hspace{0.17em}000\hspace{0.17em}000\hspace{0.17em}000\hspace{0.17em}2}$	$\mathrm{0.004\hspace{0.17em}157\hspace{0.17em}894\hspace{0.17em}736\hspace{0.17em}85}$	$\mathrm{0.000\hspace{0.17em}056\hspace{0.17em}315\hspace{0.17em}789\hspace{0.17em}474\hspace{0.17em}1}$

Table 4. Table 4: Estimation of the lower bound of the error due to model simplifications and measurement inaccuracies.

Error source	$l_{\infty }$-norm	MAE	MSE
	($\mathrm{mm}$)	($\mathrm{mm}$)	($\mathrm{mm}{}^2$)
measurement (cf. Fig. 5)	$2.2$	$\mathrm{0.397\hspace{0.17em}988\hspace{0.17em}505\hspace{0.17em}747}$	$\mathrm{0.302\hspace{0.17em}040\hspace{0.17em}229\hspace{0.17em}885}$
multi particle (no space charge)	$0.1$	$0.1$	$0.2$
space charge effect	$\mathrm{0.050\hspace{0.17em}000\hspace{0.17em}000\hspace{0.17em}000\hspace{0.17em}2}$	$\mathrm{0.004\hspace{0.17em}157\hspace{0.17em}894\hspace{0.17em}736\hspace{0.17em}85}$	$\mathrm{0.000\hspace{0.17em}056\hspace{0.17em}315\hspace{0.17em}789\hspace{0.17em}474\hspace{0.17em}1}$
step size (cf. Sec. IV)	$1.48$	$\mathrm{0.560\hspace{0.17em}344\hspace{0.17em}827\hspace{0.17em}586}$	$\mathrm{0.449\hspace{0.17em}387\hspace{0.17em}356\hspace{0.17em}322}$
sum	$3.9$	$1.1$	$0.9$

Table 5. Table 5: Design variable abbreviations and their meaning.

Design variable	Unit	Meaning
benergy	$\mathrm{GeV}$	injection beam energy
pdisft	$\mathrm{mm}$	displacement of flat top’s axis from global center
pdismain1 - pdismain4	$\mathrm{mm}$	displacement of main cavity’s axis from global center
phift	$\mathrm{deg}$	flat top cavity angle w.r.t. global coordinate system
phiinit	$\mathrm{deg}$	injection angle of beam
phimain1 - phimain4	$\mathrm{deg}$	main cavity’s angle w.r.t. the center line of sector magnet 1
phirfft	$\mathrm{deg}$	phase of flat top
prinit	$\beta \gamma$	injection radial momentum
rftshift	$\mathrm{mm}$	flat top cavity displacement in radial direction
rinit	$\mathrm{mm}$	injection radius w.r.t. the global coordinate system
rmainshift1 - rmainshift4	$\mathrm{mm}$	main RF cavity displacement in radial direction
rri2a	$\mathrm{mm}$	$a$ of Eq. (1) for RRI2
rri2phi	$\mathrm{deg}$	$\phi$ of Eq. (1) for RRI2
rri2shift	$\mathrm{mm}$	start position of RRI2 in radial direction ($>0:$ outwards)
rrla	$\mathrm{mm}$	$a$ of Eq. (1) for RRL
rrlphi	$\mathrm{deg}$	$\phi$ of Eq. (1) for RRL
rrlshift	$\mathrm{mm}$	start position of RRL in radial direction ($>0:$ outwards)
tc01mb - tc16mb	$T$	trim coil maximum magnetic field
vchange	$\mathrm{MV}$	extra RF voltage change of main cavities (in total)
vftcav	$\mathrm{MV}$	RF voltage on flat top cavity
vmaincav1 - vmaincav4	$\mathrm{MV}$	RF voltage on main cavity 1 - 4

Equations

$$\left(\begin{array}{c}x\\ y\end{array}\right)=s\cdot \left(\begin{array}{c}[r]\cos (\phi )\\ \sin (\phi )\end{array}\right)+\left(\begin{array}{c}{x}_0\\ y_0\end{array}\right)$$

$$\left(\begin{array}{c}{x}_0\\ y_0\end{array}\right)=a\cdot \left(\begin{array}{c}[r]\sin (\phi )\\ -\cos (\phi )\end{array}\right)$$

$$r_n=r_0+⟨\frac{dr}{dn}⟩n+A\sin (2\pi {\nu }_{r}n+\varphi ),$$

$$\Delta B_k\sim -\frac{qB(r)V(r)r}{E(r)\gamma (r)(\gamma (r)+1)}\frac{d\sin (\phi )}{dr}$$

$$\Delta B_k\sim -\frac{d\sin (\phi )}{dr}$$

$$[\Delta \sin (\phi )](r)\approx \frac{f(r)}{g(r)}=\frac{{\sum }_{i=0}^n{a}_i{r}^i}{{\sum }_{j=0}^m{b}_j{r}^j}$$

$$B(r)\sim \frac{d}{dr}\frac{f(r)}{g(r)}=\frac{f(r)g^{\prime }(r)-f^{\prime }(r)g(r)}{g^2(r)},$$

$$\frac{df}{dt}=\frac{\partial f}{\partial t}+\mathbf{v}\cdot {\nabla }_{\mathbf{x}}f+\frac{q}{m_0}\left(\mathbf{E}+\mathbf{v}\times \mathbf{B}\right)\cdot {\nabla }_{\mathbf{v}}f=0,$$

$$\text{mbox}TC(r)=B_{\max }\frac{{\sum }_{i=0}^n{a}_i{r}^i}{{\sum }_{j=0}^m{b}_j{r}^j}$$

$$\min \underset{i=1\dots N}{\max }\left|r_i^m-r_i^s\right|,$$

$${\sigma }_{[l,u]}=\frac{1}{N}\sqrt{\sum _{i=l}^u(r_i^m-r_i^s{)}^2}$$

$${\sigma }_{[l,u]}=\underset{i=l\dots u}{\max }\left|r_i^m-r_i^s\right|$$

$$\underset{i=1,...,N}{\min }{\left(\sum _{j=1}^M{\sigma }_j\right)}_i.$$

Full text

I Introduction

The PSI Ring cyclotron was commissioned in the mid-1970s and has been in user operation since. It has a long history of upgrades and improvements that made it possible to operate the machine with currents of up to $2.4\text{\,}\mathrm{m}\mathrm{A}$, a figure exceeding the design specification by a factor of $24$. But operation and development of accelerators over several decades is challenging in many ways. Keeping documentation up-to-date has been proven very challenging and in some cases even impossible. Beamlines or insertion devices have been modified or replaced, sputtering processes change the form of collimators and apertures. However, a post commissioning survey of possibly activated accelerator components is difficult and risky as it requires a partial or total disassembly, hence components may not be accessible with reasonable effort. Thus, it is no surprise that the work to improve, optimise, replace or test models of accelerators continues even after decades of successful operation.

Here we report our efforts to model and fit the turn pattern of the PSI Ring cyclotron, for which accurate magnetic field data of the trim coil fields are not available. Instead the average field profiles of the trim coils are derived from measurements of beam phase shifts. If a local field change by some inner trim coil causes a local radial shift of some turn, this subsequently leads, in combination with a slight difference in betatron tune, -phase and -amplitude, to a significant difference in the overall turn pattern. This becomes more and more significant from turn to turn. Hence, the overall fit should be most sensitive to the innermost trim coils.

Besides the contributions of the trim coils to the total magnetic field, the accuracy of the voltage profiles of the RF resonators plays a crucial role for the exact form of the turn pattern too. The PSI Ring cyclotron is equipped with 5 RF cavities, i.e. 4 main cavities and a third harmonic flattop cavity (see Fig. 1). While all main cavities could, in principle, have the same field profiles, they are not always operated at the exact same voltage. Furthermore, due to the sheer size of the Ring cyclotron, the exact position of the cavities might slightly differ from one to another. Hence, a fit of the turn pattern in the Ring must very likely allow for (small) variations of the positions and voltages of the cavities.

Various computer codes are able to compute turn patterns of cyclotrons. A survey of the most common cyclotron codes is given in Smirnov (2017). A first step towards a realistic numerical model of the PSI Ring cyclotron using OPAL was taken in Bi et al. (2011), but only the last few turns before extraction and one of the 18 trim coils were included. Another preceding study Pogue et al. (2017) has shown, however, that it is a challenging task to match all $182$ turns and likely requires more free parameters with at least all trim coil amplitudes, but also cavity voltages and possibly cavity alignment errors.

Due to the large design and objective spaces a simple parameter tweaking by hand is infeasible. A remedy is the use of global optimization techniques. In this paper an available framework of an evolutionary algorithm is applied that is complemented with a local search. Several papers such as Huang and Safranek (2014); Pang and Rybarcyk (2014); Yang et al. (2009); Husain and Ghodke (2018); Bazarov and Sinclair (2005); Yang et al. (2011); Hofler et al. (2013) already showed the successful application of evolutionary algorithms like particle swarm, differential evolution or NSGA-II Deb et al. (2002) in connection with particle accelerator modelling.

In this paper a general trim coil model is presented that is integrated into the particle accelerator framework OPAL Adelmann et al. (2019). It allows a more realistic description of the magnetic fields based on measured data. Together with the built-in multi-objective genetic algorithm (MOGA) and a local search it was possible to match all turns of the PSI Ring cyclotron to a maximum absolute error of $4.54\text{\,}\mathrm{m}\mathrm{m}$. With the exception of Baumgarten (2017); Pogue et al. (2017) the authors aren’t aware of any paper that tries to match the measured turn pattern with simulation in context of cyclotrons.

The paper is structured as follows: In Sec. II the aforementioned cyclotron is described and the new trim coil model is explained before the next chapter discusses its modelling and implementation within OPAL. The results of both approaches, i.e. local search and MOGA, are shown in Sec. IV with a closer discussion in Sec. V. Final remarks and a conclusion are gathered in the last section.

II PSI Ring Cyclotron

Fig. 1 shows the eight sector (SM1 - SM8) Ring cyclotron at PSI that is the last accelerating stage of the HIPA (High Intensity Proton Accelerator) facility accelerating routinely $2.2\text{\,}\mathrm{m}\mathrm{A}$ (max. $2.4\text{\,}\mathrm{m}\mathrm{A}$) proton beams with the four main cavities (Cav. $1$ - $4$) and one flat top cavity (Cav. 5) at $50.65\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$ from $72\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}590\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$. A beam is injected at an azimuth of $110\text{\,}\mathrm{\SIUnitSymbolDegree}$ and a radius around $2\text{\,}\mathrm{m}$. After typically $182$ turns it is extracted and guided to several targets to produce either muons (through pion decay) or neutrons.

II.1 Radial probes

The cyclotron is equipped with a total of $5$ wire probes that enable beam profiles Dölling . However, here we use only data of the probes RRI2 (turns $116$) and RRL (turns $9182$).

The beam profiles are obtained by measuring the current of the wire while it moves radially through the median plane and crosses subsequently multiple turns. The step width of the RRI2 and RRL probes are $0.05\text{\,}\mathrm{m}\mathrm{m}$ and $0.5\text{\,}\mathrm{m}\mathrm{m}$ respectively. The RRL probe has a single vertical carbon wire, while the RRI2 probe has $3$ carbon wires, $2$ crossed and $1$ vertical, the combination of which enables to obtain information about shape and vertical position of the beam. Here we use exclusively information of the vertical wires. Examples of the (normalized) profile measurement for the RRI2 and RRL probes are given in Fig. 3 and Fig. 4 respectively. The wire position, relative to the machine center at $(0,0)$, is described by (cf. Fig. 2)

[TABLE]

with azimuth $\varphi$ and $s\in\left[s_{1},s_{2}\right]$ and offset to the origin

[TABLE]

where $a\in\mathbb{R}$.

II.2 Peak detection of probe measurement

To determine the radial beam position at each turn the radial profile from the wire probe needs to be analyzed. This is done with a robust and straightforward peak detection algorithm that searches for downward zero-crossings in the smoothed first derivative with thresholds on minimum peak value, area and slope. The identified peaks of the measurements are indicated in Fig. 3 and Fig. 4 with red dots.

In order to estimate the error of the measurements the reference measurement of Fig. 4 was compared to measurements with a lower and higher beam current with the same machine condition. The histogram of the changes in peak positions is shown in Fig. 5. As seen in the figure a change in the beam current does not influence the peak positions significantly.

II.3 Measurement of centered beam

The beam is extracted from the PSI Ring cyclotron using an electrostatic extractor, the septum that is located in the gap between the last two turns. The standard production setup of the PSI Ring cyclotron makes use of a non-centered beam such that the beam gap for the septum is enlarged by the beam precession. In order to obtain a proper scan of the turn pattern with a long radial probe, the beam has first to be centered accurately enough that individual turns are well separated. Only then, it is possible to accurately count the number of turns Joho (1979); Stammbach (1993).

The beam centering of the PSI Ring cyclotron is determined by beam energy, radius and angle. The former is fixed by the extracted beam energy from Injector 2 but the latter can be manipulated by the last two injection magnets AND1 and AND2 and the voltage of the electrostatic injection channel EIC (cf. Fig. 1). The centering of the beam is quantified by a numerical analysis of the data of a radial injection probe (RRI2). The radial positions $r_{n}$ of turn number $n$ can approximately be described by

[TABLE]

where $\nu_{r}$ is the radial tune, $A$ is the betatron amplitude and $\phi$ the betatron phase. The radius gain per turn $\left\langle{\frac{dr}{dn}}\right\rangle$ can be assumed to be approximately constant over a small range of turns where adjacent turns do not overlap if $2A$ is smaller than the radius gain. For a centered beam the currents of AND1 and AND2 have to be chosen such that $A\approx 0$. A straightforward method, used also at PSI, is to measure the two-dimensional maps $A(I_{1},I_{2})$ and $\phi(I_{1},I_{2})$, where $I_{1}$ is the current in AND1 and $I_{2}$ the current in AND2 respectively. Then $A$ and $\phi$ can be interpolated.

The probe measurements used in this paper are performed with a beam intensity of $88\text{\,}\mathrm{\SIUnitSymbolMicro A}$. The beam profiles are given in Fig. 3 and Fig. 4. The corresponding turn separation, i.e. the distance between neighboring turns, at the probes is shown in Fig. 6 and Fig. 7. At injection the turn separation is at maximum $27\text{\,}\mathrm{m}\mathrm{m}$ which shrinks to $6.094\,210\,526\,32\text{\,}\mathrm{m}\mathrm{m}$ for the last $20$ turns.

II.4 Trim coils

The PSI Ring cyclotron is equipped with 18 trim coils which allow compensation of field errors and manipulation and optimization of the beam phase and isochronism. The trim coils are referred to as TC1 to TC18, from small to large radius. The trim coils TC6 to TC14 are only positioned on top of the odd numbered sector magnets, i.e. SM1, SM3 etc. The other trim coils are on all magnets.

The trim coils allow to change the field and to shift the beam phase. This can also alter the energy gain per turn. Furthermore the trim coils enable, within certain limits, manipulation of the tune diagram of the ring cyclotron and thus to either avoid resonances or change their position and influence on the beam.

As mentioned in the introduction, the lack of magnetic field data of the trim coils makes it currently impossible to model the fields accurately. Therefore, the simulation model developed in this paper is based on the average field profiles obtained by measurements of the beam phase shifts as subsequently explained.

II.4.1 Measurement fitting

The presented trim coil model is based on measurements of $\Delta\sin(\varphi)$ Adam and Joho (1974) in Eq. (4) as depicted in Fig. 8 with beam phase $\varphi$. As stated in Parfenova et al. (2016), the beam phase relates to $\Delta B_{k}$, the magnetic field change due to trim coil $k$, by

[TABLE]

with radius $r$, magnetic field $B(r)$, energy gain $V(r)$, charge $q$, kinetic energy $E(r)$ and relativistic factor $\gamma(r)$. In the development of the trim coil model the simplified relation

[TABLE]

was used instead. Since the neglected factor of Eq. (2) varies little over the radial range of a single trim coil, the negligence in Eq. (3) doesn’t deteriorate the model. In order to obtain the magnetic field magnitude as an additional degree of freedom the numerical model relies on normalized fields as discussed later. Each trim coil phase data was approximated by a rational function, i.e.

[TABLE]

with $m,n\in\mathbb{N}_{0}$ and $m>n$. The coefficients were computed by a Python script using the non-linear least-squares method of SciPy Jones et al. (01) where $n=2$, $m=4$ for trim coils TC2 - TC15 and $n=4$, $m=5$ for TC1 and TC16 - TC18 respectively. The fits of the data are shown in Fig. 9. The selection of the parameters $n$ and $m$ was done empirically trying to keep the polynomial degree small. As a result of Eq. (4) the corresponding magnetic field is therefore simply given by

[TABLE]

with $f^{\prime}(r)\equiv df(r)/dr$. The normalized magnetic field and its derivative of each trim coil are depicted in Fig. 10.

III OPAL

The open source library OPAL Adelmann et al. (2019) is a parallel electrostatic Particle-In-Cell (ES-PIC) framework for large-scale particle accelerator simulations. In the sequel the OPAL-cycl flavor that is used for all simulations in this study is also referred as OPAL. The particles are evolved in time $t$ by either a fourth order Runge-Kutta or a second order Leapfrog according to the collisionless Boltzmann (or Vlasov-Poisson) equation

[TABLE]

with charge $q$, mass $m_{0}$ and the six-dimensional particle density function $f(\mathbf{x},\mathbf{v},t)$ with $(\mathbf{x},\mathbf{v})\in\mathbb{R}^{3\times 3}$. The electromagnetic fields $\mathbf{E}\equiv\mathbf{E}(\mathbf{x},t)$ and $\mathbf{B}\equiv\mathbf{B}(\mathbf{x},t)$ consist of a bunch internal and external contribution. The bunch self-field is obtained in the beam rest frame by either a FFT Poisson solver or a Smoothed Aggregation Algebraic Multigrid (SAAMG) Adelmann et al. (2010) solver that is able to handle arbitrary accelerator geometries.

The following subsections highlight three features of OPAL that were used as well as extended for the purpose of this study.

III.1 Probe element and peak detection

In the cyclotron flavor of OPAL the probe is a special element placed on the midplane in Cartesian coordinates to record particles in simulations. The origin is the machine center, therefore, the positioning according to Eq. (1) with Fig. 2 is directly applicable. The OPAL syntax is shown in Fig. 11. Since OPAL has fixed time steps, particle position recording at the probe is done by a linear extrapolation of the particle direction from the closest tracking point towards the probe axis. In order to compare measurement and simulation the original description in OPAL was extended to write a particle histogram and a file collecting the peak locations. In single particle tracking the peak and thus turn detection is trivial. The localization of a turn in a multi particle simulation is achieved by summing up all radii where particles hit the probe. The mean is then determined as the orbit radius.

III.2 Multi-objective optimization

Since release version 2.0.0, OPAL is equipped with a multi-objective genetic algorithm NSGA-II (non-dominated sorting genetic algorithm) implementation Ineichen (2013). The new OPAL feature was already applied in Neveu et al. (2018). Initially a population of $n$ randomly spawned individuals within a predefined hyperspace of design variables is created after which a new set of individuals for the next population within the same bounds is selected by mixing the current $k$-fittest individuals based on a crossover pattern as well as a gene mutation algorithm. The fitness of an individual is determined by the evaluation of user-defined functions on the objectives that in this paper are the peak differences between measurement and simulation. In this regard the framework was extended to be able to read peak files and to evaluate the $l_{\infty}$-error and $l_{2}$-error norms on the parsed peak data. Furthermore, an expression to check the number of turns in simulations was added that served as an individual constraint. One of the main benefits of MOGA is its ability to search in parallel as it can evaluate all individuals of a single generation at the same time.

III.3 New trim coil model

The existing trim coil model in OPAL Bi et al. (2011) was especially designed to fit the shape of TC15 of the PSI Ring cyclotron because its interest focused on the turns close to extraction. Furthermore, the field was contributed not only local to the sector magnets but continuous smeared out on $360\text{\,}\mathrm{\SIUnitSymbolDegree}$. The new model uses a more general description by rational functions as described in Sec. II.4.1. This representation of the field allows a simple analytical differentiation to obtain the necessary derivative for the magnetic field interpolation to the position of each particle. That way the model is not restricted to the specific shape of TC15 in Bi et al. (2011).

The new trim coil model does not support an azimuthally limited field definition. However, the trim coil fields are restricted to the sector magnets by a user-defined threshold that is the lower limit to apply the additional fields. The implementation of the trim coils assumes normalized polynomial coefficients such that the maximum value of the field is $1.0$, thus, the maximum field strength $B_{\max}$ is an auxiliary tuning parameter, i.e.

[TABLE]

with $n,m\in\mathbb{N}_{0}\wedge r\in\left[r_{\min},r_{\max}\right]$. Secondly, the trim coil field is restricted in the radial direction by two extra parameters $r_{min}$ and $r_{max}$ to allow more flexibility. Nevertheless, the bounds have to be selected carefully to avoid a discontinuity in the magnetic field. In azimuthal direction the implementation uses a linear decaying field to prevent the previously mentioned issue. Since the functions $f(r)$ and $g(r)$ in Eq. (4) are polynomials in radius $r$, the derivative Eq. (5) is a rational function again. The model can therefore accept either the phase or magnetic field as input.

A template of a trim coil definition in an OPAL input file is given in Fig. 12. The parameter TYPE specifies if the polynomial represents the phase PSI-PHASE or the magnetic field PSI-BFIELD. In order to be applied the trim coil elements have to be appended to a list in the cyclotron command as depicted in Fig. 13.

IV Turn Pattern Matching

A measure for the quality of the pattern matching is the maximal peak difference between measurement $m$ and simulation $s$, i.e.

[TABLE]

where $N$ is the number of turns and $r_{i}^{m}$ and $r_{i}^{s}$ are the $i$-th turn radii.

An iterative process to get a model that is in good agreement with measurements applied multi-objective optimization and local search. Furthermore, the input parameter space between different optimizations was flexible, i.e. design variables (DVARs) were added and removed. The selection of the design variables is described in detail in the subsequent section.

Due to the decrease of the turn separation as discussed in Sec. II.3 and the increase of the circumference of the machine the choice of the number of steps per turn in simulation needs to be chosen carefully to obtain reliable results. Furthermore, the extrapolation method that is used to get the point where the particle hits the probe depends also on this time discretization. After a comparison between different number of steps per turn and a reference simulation with $23\,040$ integration steps per turn (cf. Fig. 15 and Fig. 16), the optimal number of steps per turn for the fourth order Runge-Kutta integrator w.r.t. accuracy and runtime turned out to be $2880$. It compares to the reference simulation with a maximum absolute error of $1.48\text{\,}\mathrm{m}\mathrm{m}$, mean absolute error (MAE) $0.560\,344\,827\,586\text{\,}\mathrm{m}\mathrm{m}$ and mean squared error (MSE) $0.449\,387\,356\,322\text{\,}\mathrm{m}\mathrm{m}\mathrm{{}^{2}}$. The reference simulation is selected based on the observation that the turn radii differ only in the order of $\mathcal{O}($0.19$)\$\mathrm{m}\mathrm{m}$$ at the probes compared to $11\,520$ steps per turn. A significant improvement is only achieved with an integrator of higher order.

IV.1 Design variable selection

As previously mentioned the selection of the design variables wasn’t obvious at first. Initially, the beam injection parameters and the peak magnetic field of the first trim coils were considered in order to match the first turns at injection. However, the idea of matching basically turn by turn starting at injection failed soon since the turn pattern difference between measurement and simulation started to diverge at later turns due to the wrong energy gain per turn. As a consequence the RF cavity parameters together with a constraint on the number of completed turns were added to guide the optimization towards solutions with the right number of turns. For RF cavities their voltages and positions were varied where a position encompasses the angle, radial position and displacement from the global center. For the flat top cavity also the phase angle was added as a parameter. The angle between RRI2 and RRL and their radial position were varied in order to smooth the transition between the probes.

Since TC18 is turned off, it was not added as design variable. Also TC17 that influences the last few turns was discarded since the extraction channel is not simulated. However, these turns are still corrected in simulation by TC16.

The final list of $48$ DVARs is given in Tab. 5 of the appendix. As a further clarification they are also depicted in the drawing of the cyclotron in Fig. 14. The angle between the probes RRI2 and RRL as indicated by ⑦ in the plot is adjusted using the variables $a$, $\varphi$ of Eq. (1) while keeping the length of each probe $s\in\left[s_{0}+t,s_{1}+t\right]$ with $t\in\mathbb{R}$ fixed.

IV.2 Model simplifications

Beside numerical approximations on the design of the new trim coil model the large number of DVARs (i.e. $48$) and objectives (i.e. $182$ turns) required further simplifications on the optimization approach as explained in the following.

IV.2.1 Aggregation of turns

In case of the PSI Ring cyclotron the number of objectives (i.e. turns) is $182$. In order to reduce this space multiple turns were clustered to single objectives $\sigma_{[l,u]}$ with turns in the range $[l,u]\in[1,182]$ by either the $l_{2}$-error norm

[TABLE]

with $N=u-l+1$ the number of aggregated turns or the $l_{\infty}$-error norm

[TABLE]

where $r_{i}^{m}$ and $r_{i}^{s}$ are the $i$-th turn radii of the measurement and simulation respectively. The $l_{\infty}$-norm suits our definition of the measure for the pattern matching quality, i.e. Eq. (6).

IV.2.2 Reduction of trim coil support

Due to the field overlap of neighboring trim coils a valid assumption is the partial cancellation of the field tails. That’s why the model uses a reduced radial support. Only trim coil TC1 uses the full range on the lower half.

IV.2.3 Location of trim coils

In order to limit the trim coil field in the azimuthal direction the user provides a lower bound of the magnetic field by the attribute TRIMCOILTHRESHOLD (cf. Fig. 13) above which the trim coil field is applied. This is a limitation of the new model since the real machine provides the field of all trim coils only on specific sector magnets (cf. Sec. II.4).

IV.2.4 Single particle tracking

The radial profiles are measured using a low intensity beam, i.e. $88\text{\,}\mathrm{\SIUnitSymbolMicro A}$. The negligence of space charge in order to lower the time to solution of a single simulation is therefore a reasonable assumption. A further simplification to single particle tracking is motivated by the observation that peaks are detected at the centroid of the beam (cf. Fig. 3 and Fig. 4). This reduces the time to model the full machine to approximately $2\text{\,}\mathrm{s}$ on a single core.

IV.3 Multi-objective optimization

The dimension of the design variable space required a rather large number of individuals per generation in order to sample the space sufficiently. All optimizations were performed on Piz Daint CSCS (2019), a supercomputer of the Swiss National Supercomputing Centre (CSCS). Due to its hardware architecture where a node is equipped with $2$ Intel Xeon E5-2695 v4 $@$ $2.10\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ ($2\times 18$ cores, $64\text{/}128\text{\,}\mathrm{G}\mathrm{B}$ RAM) processors, a total number of $8062$ individuals was selected in which two cores of the $224$ nodes ($$36$\cdot$224$-2=8062$) were reserved for individual post-processing and bookkeeping. Since a single objective according to Sec. IV.2.1 didn’t perform well, the turns were grouped into a total of finally six objectives. A reason might be local optima from which MOGAs are more likely to escape as discussed in Ishibuchi et al. (2006). While RRI2 was kept as a single objective, RRL was split into five objectives where each had approximately the same amount of turns and was influenced by a single trim coil only.

The optimization consisted of several independent runs with initially large bounds for each design variable. These bounds were narrowed according to the best individual. A best individual per generation is defined as the smallest sum of all $M$ objectives $\sigma_{j}$, i.e.

[TABLE]

We only show the evolution of the last optimization in Fig. 17 that was stopped after $79$ generations since after $26$ generations no significant error reduction was observed.

The objective values of the best individual over all generations are summarized in Tab. 1. According to Fig. 6 the smallest turn separation at RRI2 is 18\text{\,}\mathrm{m}\mathrm{m}$\gg\sigma_{[1,16]}=$6.38\text{\,}\mathrm{m}\mathrm{m}, where the symbol $\sigma_{[l,u]}$ indicates a single objective for the turns $l$ to $u$. Therefore, the deviation to the measurement is less than half a turn for the maximum absolute error. In case of RRL the difference is also always below the turn separation (cf. Fig. 7).

IV.4 Local search

While the genetic algorithm can in principle search a large variable space effectively, it was not able to find a nearby better solution in a reasonable amount of time. We suspect this is due to high sensitivity of the design variables near the optimum, the heuristics of genetic algorithms and the large dimensionality of the design variable space. Therefore, once the best individual from the genetic algorithm was selected, a local search around this individual was done to find the optimum. The chosen local search involved changing a single parameter value iteratively. This approach reduced the turn pattern error significantly.

Defining a good metric for the search was crucial since the iterative search is likely to stop in a local optimum. To avoid local optima several norms were used simultaneously, namely the maximal error ($l_{\infty}$-error), the second largest error, the third largest error, and a weighted $l_{2}$-error where the RRI2 turns were weighted equally to the RRL turns. A parameter was allowed to change when there was an improvement in any of the norms while not worsening the other norms significantly ($0.01\text{\,}\mathrm{m}\mathrm{m}$ for the $l_{\infty}$-error). This is equivalent to a multi-objective optimization. In Fig. 18 the $l_{\infty}$-error is shown during the iterative search. It can be seen that there was a significant improvement in the beginning, reducing the error from more than $6\text{\,}\mathrm{m}\mathrm{m}$ to less than $5\text{\,}\mathrm{m}\mathrm{m}$. It can also be seen that the error occasionally increases which avoids the local optima. Once no improvement in the $l_{\infty}$-error over $30\,000$ iterations was observed the local search was stopped.

In Fig. 19 and Fig. 20 the effect on the $l_{\infty}$-error and $l_{2}$-error per design variable is shown. It can be seen that as expected the trim coils generally improve the $l_{2}$-error, while the first trim coils and RF improve the $l_{\infty}$-error. The explanation for the latter is that the largest mismatch during the scan was often in one of the first turns.

V Discussion

The starting point of the local search was the best individual of the MOGA as explained in Sec. IV.4. This additional step could reduce the error spread and maximum absolute error (at turn $2$) compared to the measurement (cf. Tab. 2). The error of the turn radius of both methods is shown in Fig. 21 and Fig. 22.

The result of the single particle local search is verified with two multi particle tracking simulations at $88\text{\,}\mathrm{\SIUnitSymbolMicro A}$ having $360\,000$ macro particles and either space charge (i.e. FFT Poisson solver) switched on or off. The multi particle tracking (no space charge) changes the error compared to the measurement only slightly. The maximum absolute error is increased by $0.1\text{\,}\mathrm{m}\mathrm{m}$ in comparison to the single particle simulation. The mean absolute error (MAE) and mean squared error rise only by $+$0.061,789,474\text{,}\mathrm{m}\mathrm{m} and $+$0.1813\text{\,}\mathrm{m}\mathrm{m}\mathrm{{}^{2}} respectively.

A multi particle tracking simulation with space charge doesn’t change the turn pattern perceptibly (cf. Tab. 3). The $l_{\infty}$-error between both multi particle simulations differs by $0.050\,000\,000\,000\,2\text{\,}\mathrm{m}\mathrm{m}$ and MAE $0.004\,157\,894\,736\,85\text{\,}\mathrm{m}\mathrm{m}$. These observations confirm the model assumptions to neglect space charge and to use a single particle only in order to match the turn pattern.

An estimation of the error due to the measurement and model simplifications is given in Tab. 4. The systematic error is $3.9\text{\,}\mathrm{m}\mathrm{m}$ which is comparable to $4.54\text{\,}\mathrm{m}\mathrm{m}$ of the local search (cf. Tab. 2). The MAE also differs only by $0.3\text{\,}\mathrm{m}\mathrm{m}$. The difference of the MSE is, however, $2.5\text{\,}\mathrm{m}\mathrm{m}$.

VI Conclusion

A realistic simulation of existing cyclotrons heavily depends on measured data and the accurate parameter specification of the machine. Furthermore, the precision of the numerical models of devices such as RF cavities, radial probes as well as the time discretization have an impact on the result. While the latter is improved by a higher resolution at the expense of longer time to solutions of the simulation or a more accurate time integrator, numerical models of devices are only enhanced by better methods.

A new, more realistic trim coil model in the beam dynamics code OPAL that supersedes the model developed in Bi et al. (2011) was presented. The model uses rational functions to describe the shape of the trim coil field in radial direction. Although the model was applied to the PSI Ring cyclotron, it is applicable to any circular type of machine.

Thanks to the flexibility of the model, it could be used to match the turn pattern of the simulation with measurements of a centered beam in the PSI Ring cyclotron. In order to match all 182 turns a multi-objective optimization was applied with a parameter space including 16 trim coils and 32 other design variables, such as beam injection energy, RF cavity voltages and various element positions. The full list of design variables is given in Tab. 5. This process was complemented with a local search starting from the best individual of the MOGA. That way the absolute error between simulation and measurement could be reduced to at most $4.54\text{\,}\mathrm{m}\mathrm{m}$. Despite several simplifications on the optimization procedure, the multi particle tracking without space charge verified the matching of the single particle tracking. Nevertheless, the numerical model can be improved further, especially the azimuthal location of the trim coil field could be enhanced and, if possible, a 3-dimensional representation is aimed. In addition future work will include the matching of a non-centered beam and the beam profile.

The proposed approach of multi-objective optimization is unique to this kind of problem and might be used as a guideline for future projects. One of them is the DAE$\delta$ALUS project Alonso and DAE$\delta$ALUS Collaboration (2013); Adelmann and DAE$\delta$ALUS Collaboration (2013), a proposed search for CP-violation in the neutrino sector. The DAE$\delta$ALUS Superconducting Ring Cyclotron (DSRC) shares many similarities with the PSI Ring cyclotron. Future design studies will benefit greatly from the newly developed methods that were presented here.

Acknowledgements.

The authors thank W. Joho, M. Humbel, R. Dölling and H. Zhang for their helpful discussions about trim coils and probe measurements. The contribution of D. Winklehner regarding the DAE$\delta$ALUS is gratefully acknowledged as well. Furthermore, we appreciate the expertise of M. Kranjčevic regarding the multi-objective optimization using evolutionary algorithms. All optimizations were performed using the HPC resources provided by the Swiss National Supercomputing Centre (CSCS). This project is funded by the Swiss National Science Foundation (SNSF) under contract number 200021_159936.

Appendix A Design Variables