O

# Introduction

Large-scale multi- objective design optimization enables the automated identification of working points in high dimensional search and decision spaces. The design, commission, and operation of an accelerator facility is an excellent example of such a non-trivial problem.

Such a task involves significant computer modeling using accelerator simulation codes such as PSI’s OPAL (Object Oriented Parallel Accelerator Library) framework. Despite the fact that these codes are parallel, typical simulation parameters (e.g. a small number of macroparticles and mesh size) for individual runs, limit their scalability to several hundred or a few thousand processors. This represents a strong impediment in view of the petascale regime, therefore parallelisation on multiple levels, e.g. running multiple parallel simulations in parallel, becomes a necessity. On the other hand, introducing a low-dimensional scalable model enables multi-resolution simulation runs.

optPilot is a general-purpose framework for simulation-based multi-objective optimization methods that allows the automatic investigation of optimal sets of parameters.

The implementation is based on a Master/Slave paradigm, employing several masters and groups of workers to prevent communication hotspots at master processes. In addition, we exploit information about the underlying network topology when placing master processes and assigning roles. Solution states are exchanged between masters in a "rumor routing" fashion on a social network graph using one-sided communication. Using evolutionary algorithms and OPAL simulations as optimizer and forward solver in our framework, we demonstrate the feasibility and scalability of our approach on real applications in the domain of particle accelerators.

# Running

• set ENV variables (only required for OPAL):

• FIELDMAPS: folder containing fieldmaps
• TEMPLATES: directory containing the template input file
• SIMTMPDIR: directory for temporary creation of simulation files (must exist)
• in case of more than one optimiser

• --sol-synch: after how many generations solutions between the optimisers are exchanged
• set program arguments:

• --inputfile=fname: input file containing optimization problem
• --outfile=fname: name used in output file generation
• --outdir=dirname: name of directory used to store generation output files (generated if non-existing)
• --initialPopulation=num: size of the initial population
• --num-masters=num: number of master nodes
• --num-coworkers=num: number processors per worker
• --selector=path: path of the selector (PISA only)
• --dump-dat=freq: dump old generation data format with frequency (PISA only)
• --num-ind-gen=num: number of individuals in a generation (PISA only)
• and for convergence we support:

• --maxGenerations=num: number of generations to run
• --epslion=num: tolerance of hypervolume criteria
• --expected-hypervol=num: the reference hypervolume
• --conv-hvol-prog=num: converge if change in hypervolume is smaller
• run, stand-alone example:

  mpirun -np 4 opt-pilot.exe \
--inputfile=FiPha3Opt1.tmpl --outfile=results.dat \
--maxGenerations=5 --initialPopulation=10 --num-masters=1
• run, OPAL equivalent (note the additional duplicated inputfile):

  mpirun -np 4 opal \
--inputfile=FiPha3Opt1.tmpl --outfile=results.dat \
--maxGenerations=5 --initialPopulation=10 --num-masters=1 \
FiPha3Opt1.tmpl

# Basic Syntax

One needs to define the design variables, objectives and constraints one by one:

d1: DVAR, VARIABLE="x1", LOWERBOUND="-1.0", UPPERBOUND="1.0";
d2: DVAR, VARIABLE="x2", LOWERBOUND="-1.0", UPPERBOUND="1.0";
d3: DVAR, VARIABLE="x3", LOWERBOUND="-1.0", UPPERBOUND="1.0";

This defines three design variables named d1, d2 and d3. The variable names x1, x2 and x3 can be used in objectives and constraints. Bounds for design variables should always be given. Note the string notation.

obj1: OBJECTIVE, EXPR="-energy";
obj2: OBJECTIVE, EXPR="emit_x";

This defines two objectives named obj1 and obj2. The optimiser knows several mathematical functions and methods to access output files (see below). Without the use of these special methods design variables or the variables from the .stat output file at the end of the simulation can be used for optimisation. Note that objectives are always minimised, so in this case a solution where the final energy is maximal and the final horizontal emittance is minimal is looked for.

con1: CONSTRAINT, EXPR="x1+x2<1.0";
con2: CONSTRAINT, EXPR="numParticles>1000";

This defines two constraints. The syntax is similar to the OBJECTIVE syntax. The first constraint consists of only design variables and will be evaluated before the simulation. The second constraint will be evaluated after the simulation.

After defining the individual elements, one needs to define the sets of design variables, objectives and constraints:

dvars:   DVARS       = (d1, d2, d3);
objs:    OBJECTIVES  = (obj1, obj2);
constrs: CONSTRAINTS = (con1, con2);

Zero constraints can be defined as:

constrs: CONSTRAINTS = ();

Finally, the optimiser command should be added:

opt: OPTIMIZE, OBJECTIVES=objs, DVARS=dvars, CONSTRAINTS=constrs;

# Expression syntax

The EXPR The optimiser parser knows the following mathematical functions (which are mapped directly to the STL functions with the same name):

• sqrt(x) : square root of x
• pow(x,k) : x to the power k
• exp(x) : e to the power x
• log(x) : natural logarithm of x
• ceil(x) : round x upward to the smallest integral value that is not less than x
• floor(x) : round x downward to smallest integral value that is not greater than x
• fabs(x) : absolute value of x
• fmod(x,y): floating point remainder of x/y
• sin(x) : sine of angle x (in radians)
• asin(x) : the arcsin of x (return value in radians)
• cos(x) : cosine of angle x (in radians)
• acos(x) : the arc cosine of x (return value in radians)
• tan(x) : tangent of angle x (in radians)
• atan(x) : the arc tangent of x (return value in radians)

In addition the optimiser parser has one non-STL function:

• sq(x) : square of x

There are several methods to access output data:

• fromFile(filename) : reads vector data from file. If the file contains more than one value the sum is returned.
• sddsVariableAt(variable,s,filename) : interpolates a variable (first argument) from SDDS file (third argument) near a specific s position (second argument)
• sameSDDSVariableAt(variable,s) : same as sddsVariableAt for the main .stat output file.
• sumErrSq(measfilename,variable,filename) : computes the square root of the sum of all measurement errors (given as first and third argument) for a variable (second argument) according to: result = \frac{1}{n} * \sqrt{\sum_{i=0}^n (measurement_i - value_i)^2}.
• radialPeak(filename,n) : returns n-th peak of a radial probe output file.
• sumErrSqRadialPeak(measfilename,filename,start,end) : computes the square root of the sum of all peak errors (given as first and second argument) for a range of peaks (third argument and fourth argument) according to: result = \frac{1}{n} * \sqrt{\sum_{i=start}^{end} (measurement_i - value_i)^2}.
• probVariableWithID(variable,id,probelossfile) : returns the value of the variable (first argument) with a certain ID (second argument) from probe loss file (third argument).

# FODO Example

export PATH=~adelmann/amas/bin:\$PATH

The set of input files are located here: fodo.tgz

This is a very simple example of a fodo cell. We want to achieve in both the x and y plane a desired beam size.

OPTION, VERSION=10600;
REAL SOL  = 2.9979246E8;
REAL Pcen = 100.0E6;
REAL BRho = Pcen/SOL;
REAL QK1  = 6.2519;
REAL QSTR = QK1*BRho/2.0;
REAL qb=77.0e-12;
REAL bfreq=1300.0E6;
REAL bcurrent=qb*bfreq;

QDX1: QUADRUPOLE, ELEMEDGE=0.0,  L=0.10, K1=17.7 * BRho / 2;
QFX1: QUADRUPOLE, ELEMEDGE=0.91, L=0.20, K1=-17.5 * BRho / 2;
QDX2: QUADRUPOLE, ELEMEDGE=1.92, L=0.10, K1=17.5 * BRho / 2;

FODO_Full: LINE = (QDX1, QFX1, QDX2);

// SC calculations on:
Fs1:FIELDSOLVER, NONE = FFT, MX = 8, MY = 8, MT = 8,
PARFFTX = true, PARFFTY = true, PARFFTT = true,
BCFFTX = open, BCFFTY = open, BCFFTT = open,
BBOXINCR = 1, GREENSF = INTEGRATED;

Dist1:DISTRIBUTION, DISTRIBUTION=GAUSS,
sigmax=3E-3, sigmapx=1E-5,
sigmay=3E-3, sigmapy=1E-5,
sigmat=3E-3, sigmapt=1E-5, CORRX=0, CORRY=0, CORRT=0;

beam1: BEAM, PARTICLE = ELECTRON, pc = P0, NPART = 5000, BFREQ = bfreq,
BCURRENT = bcurrent, CHARGE = -1;

SELECT, LINE=FODO_Full;

TRACK, LINE=FODO_Full, BEAM=beam1, MAXSTEPS={6e5}, DT={1e-12}, ZSTOP={2.02};
RUN, METHOD = "PARALLEL-T", BEAM = beam1, FIELDSOLVER = Fs1, DISTRIBUTION = Dist1;
ENDTRACK;
QUIT;

With some random initial conditions we get the following solution (with the python script from pyOpalTools)

visualize_pf.py --objectives=drmsx,drmsy,QDX1_K1 --path=results-fodo-new --generation=500

// ---- OPTIMIZER SECTION -------
dv0: DVAR, VARIABLE="QDX1_K1", LOWERBOUND="0", UPPERBOUND="35";
dv1: DVAR, VARIABLE="QDX2_K1", LOWERBOUND="0", UPPERBOUND="35";
dv2: DVAR, VARIABLE="QFX1_K1", LOWERBOUND="-35", UPPERBOUND="0";
dvars: DVARS=(dv0,dv1,dv2);

drmsx:OBJECTIVE,EXPR="fabs(sameSDDSVariableAt("rms_x",2.00) - 0.0000977)";
drmsy:OBJECTIVE,EXPR="fabs(sameSDDSVariableAt("rms_y",2.00) - 0.0001833)";

objs: OBJECTIVES=(drmsx,drmsy);
constrs: CONSTRAINTS = ();
opt: OPTIMIZE, OBJECTIVES = objs, DVARS = dvars, CONSTRAINTS = constrs;

Now look at the initial pareto front:

And the final Pareto front after 1000 generations:

The solutions for each generation will be saved in either a plain ASCII or JSON format. The Pareto front of all individuals is in a JSON file, ParetoFront_.json. There are three log files opt.trace.0, opt-progress.0 and pilot.trace.0 that log the job management. E.g. to count the total number of dispatched simulations:
cat opt.trace.0 | grep dispatched | wc -l
cat opt.trace.0 | grep invalid | wc -l