Оптимизация огибающей пучка заряженных частиц с помощью генетического алгоритма
Sometimes it’s useful to use a computer to adjust the envelope.
Genetic algorithms work similarly to natural selection in nature. The basis of the genetic algorithm is the individual, chromosomes (genes), population (set of individuals), the functions of adaptability, selection, crossing and mutation. In practice, everything happens like this: the creation of the first generation, assessment, selection, crossing, mutation, new generation, evaluation, etc. until we get the desired result.
DEAP is a framework for working with genetic algorithms, containing many ready-made tools, the main thing is to know how to use them.
Create a beam and accelerator.
import kenv as kv
beam = kv.Beam(energy=2,
current=2e3,
radius=50e-3,
rp=50e-3,
normalized_emittance=1000e-6)
accelerator = kv.Accelerator(0, 5, 0.01)
Solenoids = [
[ 0.5000, 0.02, 'Bz.dat', 'Sol. 1'],
[ 1.0000, 0.02, 'Bz.dat', 'Sol. 2'],
[ 1.5000, 0.02, 'Bz.dat', 'Sol. 3'],
[ 2.0000, 0.02, 'Bz.dat', 'Sol. 4'],
[ 2.5000, 0.02, 'Bz.dat', 'Sol. 5'],
[ 3.0000, 0.02, 'Bz.dat', 'Sol. 6'],
[ 3.5000, 0.02, 'Bz.dat', 'Sol. 7'],
[ 4.0000, 0.02, 'Bz.dat', 'Sol. 8'],
[ 4.5000, 0.02, 'Bz.dat', 'Sol. 9'],
]
for z0, B0, filename, name in Solenoids:
accelerator.Bz_beamline[name] = kv.Element(z0, B0, filename, name)
accelerator.compile()
simulation = kv.Simulation(beam, accelerator)
simulation.track()
Graphic
matplotlib
import holoviews as hv
hv.extension('matplotlib')
%opts Layout [tight=True]
%output size=150 backend='matplotlib' fig='svg'
%opts Area Curve [aspect=3 show_grid=True]
%opts Area (alpha=0.25)
%opts Curve (alpha=0.5)
%opts Area.Beam [aspect=3 show_grid=True] (color='red' alpha=0.3)
import warnings
warnings.filterwarnings('ignore')
dim_z = hv.Dimension('z', unit='m', range=(accelerator.start, accelerator.stop))
dim_Bz = hv.Dimension('Bz', unit='T', label='Bz', range=(0, 0.1))
dim_Ez = hv.Dimension('Ez', unit='MV/m', label='Ez')
dim_r = hv.Dimension('r', label="Beam r", unit='mm', range=(0, 150))
z_Bz= hv.Area((accelerator.parameter,accelerator.Bz(accelerator.parameter)), kdims=[dim_z], vdims=[dim_Bz])
z_r = hv.Area(((accelerator.parameter,simulation.envelope_x(accelerator.parameter)*1e3)), kdims=[dim_z], vdims=[dim_r], group='Beam')
(z_r+z_Bz).cols(1)
As you can see, the envelope is far from perfect.
Apply the genetic algorithm.
Genetic algorithm
First we construct a model envelope.
import numpy as np
coefficients = np.polyfit([accelerator.start, (accelerator.start+accelerator.stop)/2, accelerator.stop], [0.15, 0.05, 0.03], 2)
envelope_mod = coefficients[0]*accelerator.parameter**2 + coefficients[1]*accelerator.parameter+ coefficients[2]
z_env = hv.Curve((accelerator.parameter, envelope_mod*1e3), kdims=[dim_z], vdims=[dim_r], label='Envelope_mod', group='Beam')
z_env
Connect the necessary libraries.
import random
from deap import creator, base, tools, algorithms
The first step is to create the necessary types. Usually it is fitness and the individual.
creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
creator.create("Individual", list, fitness=creator.FitnessMin)
Next, you need to create a toolkit. To do this, you need to set the basic parameters for our algorithm.
Sol_B = 0.02 # [T] average Bz in the solenoids
Sol_B_max = 0.05 # [T] max Bz
Sol_B_min = 0.005 # [T] min Bz
Sol_Num = 9 # quantity
CXPB = 0.4 # cross chance
MUTPB = 0.6 # Mutation probability
NGEN = 100 # Number of generations
POP = 100 # Number of individuals
toolbox = base.Toolbox()
toolbox.register("attr_float", random.gauss, Sol_B, ((Sol_B_max - Sol_B_min)/2)**2)
toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attr_float, n=Sol_Num)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
Create a fitness function.
def evalution_envelope(individual):
for x in range(Sol_Num):
accelerator.Bz_beamline['Sol. %.d'%(x+1)].max_field = individual[x]
accelerator.compile()
simulation = kv.Simulation(beam, accelerator)
simulation.track()
tuple(envelope_mod)
tuple(simulation.envelope_x(accelerator.parameter))
sqerrors = np.sqrt((envelope_mod-simulation.envelope_x(accelerator.parameter))**2)
return sum(sqerrors)/len(accelerator.parameter),
toolbox.register("evaluate", evalution_envelope)
toolbox.register("mate", tools.cxUniform, indpb=MUTPB )
toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=((Sol_B_max - Sol_B_min)/2)**2, indpb=MUTPB)
toolbox.register("select", tools.selTournament, tournsize=NGEN//3)
Logbook.
stats_fit = tools.Statistics(lambda ind: ind.fitness.values)
mstats = tools.MultiStatistics(fitness=stats_fit)
mstats.register("avg", np.mean)
mstats.register("std", np.std)
mstats.register("min", np.min)
mstats.register("max", np.max)
logbook = tools.Logbook()
population = toolbox.population(n=POP)
pop, logbook = algorithms.eaSimple(population, toolbox, cxpb=CXPB, mutpb=MUTPB, ngen=NGEN, stats=mstats, verbose=True)
fitness
--------------------------------------------------------------------------
gen nevals avg gen max min nevals std
0 100 0.0373577 0 0.0408331 0.0337452 100 0.00113532
1 82 0.0348016 1 0.0379858 0.0299872 82 0.00134257
2 69 0.0317311 2 0.0355503 0.0268419 69 0.00158446
3 70 0.0284586 3 0.0310798 0.0258454 70 0.00116153
4 77 0.0262835 4 0.028602 0.0245999 77 0.00073413
5 72 0.0248689 5 0.0268379 0.0234427 72 0.000541174
6 84 0.0237668 6 0.0256165 0.0221611 84 0.00067109
7 73 0.0225147 7 0.0241554 0.0203067 73 0.000610466
8 88 0.0210856 8 0.0235622 0.0194378 88 0.000755119
9 80 0.0197334 9 0.0209437 0.0190342 80 0.000337098
10 81 0.0191692 10 0.0204253 0.0183268 81 0.00035722
11 75 0.0185557 11 0.019661 0.0178986 75 0.000283725
12 78 0.0181075 12 0.0190688 0.0175483 78 0.00029574
13 72 0.0177022 13 0.018703 0.0172446 72 0.000248748
14 75 0.0174452 14 0.0183427 0.0169873 75 0.000190804
15 70 0.017232 15 0.0178082 0.0169012 70 0.000181106
16 79 0.0170371 16 0.0180127 0.0167734 79 0.000208848
17 80 0.0168772 17 0.0176019 0.0165999 80 0.000167511
18 80 0.0167085 18 0.0175836 0.0164198 80 0.000174468
19 77 0.0165316 19 0.0171106 0.0162831 77 0.000152431
20 76 0.0163782 20 0.0168837 0.0160568 76 0.000171506
21 69 0.0161959 21 0.0170411 0.015828 69 0.000194345
22 77 0.0159477 22 0.0166562 0.0156399 77 0.000157436
23 79 0.0158229 23 0.0184114 0.0155106 79 0.000314105
24 80 0.0156802 24 0.0167603 0.0153466 80 0.000234712
25 72 0.0154841 25 0.0165787 0.0152393 72 0.000184083
26 74 0.0153536 26 0.0160787 0.0150713 74 0.000201229
27 83 0.0152175 27 0.0160113 0.0147709 83 0.000237473
28 79 0.0149943 28 0.0159486 0.014535 79 0.00027873
29 80 0.0147364 29 0.0155165 0.0143068 80 0.000219007
30 78 0.0145096 30 0.0155232 0.01421 78 0.000219074
31 76 0.0143674 31 0.0162176 0.0141153 76 0.000292944
32 77 0.0142465 32 0.0154183 0.0139755 77 0.000228094
33 74 0.0141755 33 0.016872 0.0137429 74 0.000370584
34 65 0.0140263 34 0.0165405 0.0135535 65 0.000369079
35 79 0.0138775 35 0.0165023 0.0134757 79 0.000469816
36 78 0.0136499 36 0.0165048 0.0132913 78 0.00036574
37 77 0.0135281 37 0.0162187 0.0131399 77 0.000370037
38 72 0.0133894 38 0.0161935 0.013018 72 0.000407104
39 74 0.0132855 39 0.0151958 0.0129085 74 0.000356725
40 73 0.0132202 40 0.0158132 0.0127831 73 0.000502361
41 78 0.0131469 41 0.0154469 0.0126621 78 0.000549534
42 72 0.0128525 42 0.0153473 0.0125449 72 0.000371909
43 81 0.0128408 43 0.0143709 0.0124426 81 0.00037498
44 80 0.0128018 44 0.0143546 0.0122802 80 0.00042712
45 76 0.0126884 45 0.0149749 0.0122105 76 0.000484442
46 74 0.0124913 46 0.0141759 0.0121807 74 0.000374844
47 84 0.0124814 47 0.0150409 0.012149 84 0.000478495
48 65 0.012378 48 0.0143345 0.0120466 65 0.000396209
49 77 0.0123797 49 0.0141054 0.0119723 77 0.000465145
50 66 0.0122275 50 0.0136186 0.0119299 66 0.000324548
51 77 0.0122757 51 0.0142704 0.0118622 77 0.000472376
52 74 0.0121928 52 0.0138265 0.0118082 74 0.000460488
53 73 0.0121423 53 0.0145869 0.0117309 73 0.000487374
54 70 0.0120654 54 0.0132641 0.0116092 70 0.000405956
55 77 0.012012 55 0.0137579 0.0115895 77 0.00050723
56 73 0.0118692 56 0.0141939 0.0115563 73 0.000437811
57 83 0.0118798 57 0.0139667 0.0114851 83 0.000426248
58 71 0.0117967 58 0.01344 0.0114796 71 0.000419319
59 72 0.0117301 59 0.0145933 0.0114334 72 0.000445674
60 78 0.0116524 60 0.0135045 0.0113808 78 0.000333071
61 78 0.0116611 61 0.0132638 0.0113299 78 0.000437407
62 78 0.0116445 62 0.0140102 0.0113001 78 0.000467812
63 80 0.0116662 63 0.0138093 0.0112541 80 0.000496356
64 78 0.0115682 64 0.0126411 0.0112403 78 0.000360796
65 73 0.0115044 65 0.0133941 0.0112293 73 0.000428872
66 81 0.0114277 66 0.0134346 0.0112152 81 0.000315043
67 72 0.0114574 67 0.0132544 0.0111495 72 0.000390321
68 81 0.0114764 68 0.0135241 0.0110932 81 0.000486222
69 75 0.011434 69 0.0132443 0.0110874 75 0.000469529
70 74 0.0113903 70 0.0129446 0.0109944 74 0.000420948
71 80 0.0112841 71 0.0132937 0.0109944 80 0.000408004
72 75 0.0112891 72 0.0133998 0.010978 75 0.000465794
73 79 0.0112265 73 0.0122116 0.0109464 79 0.00029769
74 80 0.0112526 74 0.0129213 0.0108979 80 0.000417258
75 71 0.0112418 75 0.0128844 0.0108979 71 0.000447559
76 73 0.0111925 76 0.0138475 0.0108402 73 0.000519714
77 78 0.0111465 77 0.0132176 0.010835 78 0.000422602
78 79 0.0112091 78 0.0133698 0.0107964 79 0.000558782
79 73 0.0110862 79 0.0131931 0.0107964 73 0.000465272
80 67 0.011036 80 0.0127483 0.0107964 67 0.000386788
top = tools.selBest(pop, k=10)
gen = np.array(logbook.select("gen"))
fit_min = np.array(logbook.chapters["fitness"].select("min"))
fit_max = np.array(logbook.chapters["fitness"].select("max"))
fit_avg = np.array(logbook.chapters["fitness"].select("avg"))
for x in range(Sol_Num):
accelerator.Bz_beamline['Sol. %.d'%(x+1)].max_field = top[0][x]
accelerator.compile()
simulation = kv.Simulation(beam, accelerator)
simulation.track()
z_Bz_gen= hv.Area((accelerator.parameter,accelerator.Bz(accelerator.parameter)), kdims=[dim_z], vdims=[dim_Bz])
z_r_gen = hv.Area(((accelerator.parameter,simulation.envelope_x(accelerator.parameter)*1e3)), kdims=[dim_z], vdims=[dim_r], group='Beam')
(z_r_gen*z_env+z_Bz_gen).cols(1)
