Application of multi-objective optimization to match turn pattern measurements for cyclotrons

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 optimization 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.5 mm over all 182 turns.


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 mA, 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 [1]. A first step towards a realistic numerical model of the PSI Ring cyclotron using OPAL was taken in [2], but only the last few turns before extraction and one of the 18 trim coils were included. Another preceding study [3] 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 usage of multi-objective optimisation or a scan of the full parameter space. Several papers such as [4][5][6][7][8][9][10] already showed the successful application of evolutionary algorithms like particle swarm, differential evolution or NSGA-II [11] 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 [12]. 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.5 mm. With the exception of [3,13] the authors aren't aware of any paper that tries to match the measured turn pattern with simulation in context of cyclotrons. Especially the usage of MOGAs to achieve this objective is unique.
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. 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 mA (max. 2.4 mA) proton beams with the four main cavities (Cav. 1 -4) and one flat top cavity (Cav. 5) at 50.65 MHz from 72 MeV to 590 MeV. A beam is injected at an azimuth of 110 • and a radius around 2 m. After typically 182 turns it is extracted and guided to several targets to produce either muons (through pion decay) or neutrons.

A. Radial probes
The cyclotron is equipped with a total of 5 wire probes that enable beam profiles [14]. However, here we use only data of the probes RRI2 (turns 1 to 16) and RRL (turns 9 to 182).
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.1 mm and 0.5 mm 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 (normalised) 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) with azimuth ϕ and s ∈ [s 1 , s 2 ] and offset to the origin where a ∈ R.  [14]. Injection elements are marked in green. Figure 2. Mathematical description of probe positioning. The probe orientation is given by s with begin s1 and end s2. The offset of the probe from the machine center is indicated by the dashed line with symbol a. [14] B. Peak detection of probe measurement To determine the radial beam position at each turn the radial profile from the wire probe needs to be analysed. 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   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.

C. Measurement of centred 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-centred 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 centred accurately enough that individual turns are well separated. Only then, it is possible to accurately count the number of turns [15,16].
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 where ν r is the radial tune, A is the betatron amplitude and φ the betatron phase. The radius gain per turn dr dn 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 centred beam the currents of AND1 and AND2 have to be chosen such that A ≈ 0. A straightforward method, used also at PSI, is to measure the two-dimensional maps A(I 1 , I 2 ) and φ(I 1 , I 2 ), where I 1 is the current in AND1 and I 2 the current in AND2 respectively. Then A and φ can be interpolated.
The probe measurements used in this paper are performed with a beam intensity of 88 µA. The beam profiles are given in Fig. 3 and Fig. 4. The corresponding turn separation, i.e. the distance between neighbouring turns, at the probes is shown in Fig. 6 and Fig. 7. At injection the turn separation is at maximum 27 mm which shrinks to 6.1 mm for the last 20 turns.

D. Trim coils
The PSI Ring cyclotron is equipped with 18 trim coils which allow compensation of field errors and manipulation and optimisation 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.

Measurement fitting
The presented trim coil model is based on measurements of ∆ sin(ϕ) [17] in Eq. (4) as depicted in Fig. 8 with beam phase ϕ. As stated in [18], the beam phase relates to ∆B k , the magnetic field change due to trim coil k, by with radius r, magnetic field B(r), energy gain V (r), charge q, kinetic energy E(r) and relativistic factor γ(r).
In the development of the trim coil model the simplified relation 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 normalised fields as discussed later. Each trim coil phase data was approximated by a rational function, i.e.
with m, n ∈ N 0 and m > n. The coefficients were computed by a Python script using the non-linear leastsquares method of SciPy [19] 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 with f (r) ≡ df (r)/dr. The normalised magnetic field and its derivative of each trim coil are depicted in Fig. 10.  (ϕ)  TC1   TC2   TC3   TC4   TC5   TC6 TC7 TC8  TC9 TC10   TC11   TC12   TC13   TC14  TC15 TC16 TC17 TC18 Figure 8. Measurement of the beam phase ϕ shift due to trim coils [17].

III. OPAL
The open source library OPAL [12] is a parallel electrostatic Particle-In-Cell (ES-PIC) framework for largescale particle accelerator simulations. In the sequel the OPAL-cycl flavour 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  with charge q, mass m 0 and the six-dimensional parti- ) and B ≡ B(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) [20] 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.

A. Probe element and peak detection
In the cyclotron flavour 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 localisation 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.

B. Multi-objective optimisation
Since release version 2.0.0, OPAL is equipped with a multi-objective genetic algorithm NSGA-II (non-  dominated sorting genetic algorithm) implementation [21]. The new OPAL feature was already applied in [22]. 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 userdefined 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 ∞ -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.

C. New trim coil model
The existing trim coil model in OPAL [2] 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 • . The new model uses a more general description by rational functions as described in Sec. II D 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 [2].
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 userdefined threshold that is the lower limit to apply the additional fields. The implementation of the trim coils assumes normalised 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. TC(r) = B max n i=0 a i r i m j=0 b j r j with n, m ∈ N 0 ∧ r ∈ [r min , r max ]. 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.
where N is the number of turns and r m i and r s i are the i-th turn radii.
An iterative process to get a model that is in good agreement with measurements applied multi-objective optimisation and local search. Furthermore, the input parameter space between different optimisations 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 C 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 discretisation. 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.5 mm, mean absolute error (MAE) 0.6 mm and mean squared error (MSE) 0.4 mm 2 . The reference simulation is selected based on the observation that the turn radii differ only in the order of O(0.2) mm at the probes compared to 11 520 steps per turn. A significant improvement is only achieved with an integrator of higher order.

A. 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 optimisation 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. V 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, ϕ of Eq. (1) while keeping the length of each probe s ∈ [s 0 + t, s 1 + t] with t ∈ R fixed.

B. 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 optimisation approach as explained in the following.

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 σ [l,u]  with N = u − l + 1 the number of aggregated turns or the l ∞ -error norm where r m i and r s i are the i-th turn radii of the measurement and simulation respectively. The l ∞ -error norm turned out to suit the optimisation better.

Reduction of trim coil support
Due to the field overlap of neighbouring 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.

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 coild 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 D).

Single particle tracking
The radial profiles are measured using a low intensity beam, i.e. 88 µ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 s on a single core.

C. Multi-objective optimisation
The dimension of the design variable space required a rather large number of individuals per generation in order to sample the space sufficiently. Due to the hardware architecture of Piz Daint where a node is equipped with 2 Intel Xeon E5-2695 v4 @ 2.1 GHz (2 × 18 cores, 64/128 GB RAM) processors, a total number of 8062 individuals was selected in which two cores of the 224 nodes (36 · 224 − 2 = 8062) were reserved for individual post-processing and bookkeeping. Since there are several probes and, thus, several files, the objectives are split into RRI2 and RRL objectives. Furthermore, for the optimisation it was advantageous to split RRL into more than one objective. The measurement was therefore divided into five objectives such that all have approximately the same amount of turns and that each trim coil should influence at most one objective at the probes.
The optimisation 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 σ j , i.e. We only show the evolution of the last optimisation 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 summarised in Tab. I. According to Fig. 6 the smallest turn separation at RRI2 is 18 mm σ [1,16] = 6.4 mm, where the symbol σ [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). σ [106,148] σ [149,182] σ [1,16] σ [32,61] σ [62,105] σ [9,31] Table I. Result of best individual obtained by optimisation using the l∞-error norm for each objective. The label σ [l,u] indicates an objective for the turns in the range [l, u].

D. Local search
Once the best individual from the genetic algorithm was selected, a local search around this individual was done. The chosen local search involved changing a single parameter value iteratively. While this approach is somewhat simplistic the turn pattern matching improved 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 ∞ -error), the second largest error, the third largest error, and a weighted l 2error 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 mm for the l ∞ -error). In Fig. 18 the l ∞ -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 mm to less than 5 mm. It can also be seen that the error occasionally increases which avoids the local optima. Once no improvement in the l ∞ -error over 3000 iterations was observed the local search was stopped. In Fig. 19 and Fig. 20 the effect on the l ∞ -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 ∞ -error. The explanation for the latter is that the largest mismatch during the scan was often in one of the first turns .   tc03mb  tc02mb  tc01mb  phimain2  vmaincav3  phimain4  rftshift  phimain1  phimain3  phirfft  tc04mb  pdisft  rinit  vftcav  rri2a  phift  tc05mb  tc06mb  rrlphi  rrla  tc07mb  tc11mb  tc13mb  tc08mb  tc12mb  tc14mb  tc15mb  tc09mb  tc16mb  tc10mb  vchange  vmaincav4  vmaincav2  phiinit  vmaincav1  prinit  rri2phi

V. DISCUSSION
The starting point of the local search was the best individual of the MOGA as explained in Sec. IV D. This additional step could reduce the error spread and maximum absolute error (at turn 2) compared to the mea- tc06mb  tc08mb  tc09mb  tc07mb  tc10mb  tc12mb  tc13mb  tc11mb  tc14mb  tc15mb  tc16mb  vmaincav1  phiinit  vmaincav4  vchange  vmaincav2  rinit  prinit  tc05mb  benergy  rrla  rri2a  rri2phi  tc01mb  phift  tc02mb  phirfft  rrlphi  phimain1  vftcav  phimain2  phimain3  phimain4  rftshift  pdisft  tc03mb  tc04mb   surement (cf. Tab. II). 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 µ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 mm in comparison to the single particle simulation. The mean absolute error (MAE) and mean squared error rise only by +0.1 mm and +0.2 mm 2 respectively.
A multi particle tracking simulation with space charge doesn't change the turn pattern perceptibly (cf. Tab. III). The l ∞ -error between both multi particle simulations differs by 0.05 mm and MAE 0.00 mm. 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. IV. The systematic error is 3.9 mm which is comparable to 4.5 mm of the local search (cf. Tab. II). The MAE also differs only by 0.3 mm. The difference of the MSE is, however, 2.5 mm.

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, Table III. Maximum absolute error (l∞-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.  Figure 21. Error of the turn radius at RRI2 between measurement and simulation of the best individual obtained by multi-objective optimisation and local search.
radial probes as well as the time discretisation 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 [2] 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 centred beam in the PSI Ring cyclotron. Due to 16 trim coils and 32 other design variables, the usage of a multi-objective optimisation genetic algorithm (MOGA) was indispensable to match all 182 turns (i.e. objectives) of the machine. 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.5 mm. Despite several simplifications on the optimisation 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.
The proposed approach combining MOGA with a local search is unique to this kind of problem and might be used as a guideline for future projects. One of them is the DAEδALUS project [23,24], a proposed search for CP-  violation in the neutrino sector. The DAEδ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.

ACKNOWLEDGMENTS
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δALUS is gratefully acknowledged as well. Furthermore, we appreciate the expertise of M. Kranjčevic regarding the multi-objective optimisation using evolutionary algorithms. All optimisations 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.