Genetic algorithm¶
Due to the speed of KENV, we can apply a genetic algorithm to reconstruct the envelope from the experimental data.
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.
Genetic algorithm for building an envelope¶
As an example, let’s consider a simple problem. Sometimes it’s useful to use a computer to adjust the envelope.
Create a beam and accelerator.¶
[1]:
import kenv as kv
[2]:
beam = kv.Beam(energy=2,
current=2e3,
radius=50e-3,
rp=50e-3,
normalized_emittance=1000e-6)
[3]:
accelerator = kv.Accelerator(0, 5, 0.01)
[25]:
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'],
]
[26]:
for z0, B0, filename, name in Solenoids:
accelerator.Bz_beamline[name] = kv.Element(z0, B0, filename, name)
[27]:
accelerator.compile()
[28]:
simulation = kv.Simulation(beam, accelerator)
[29]:
simulation.track()
Graphics¶
[30]:
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')
[31]:
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))
[32]:
z_Bz= hv.Area((accelerator.z,accelerator.Bz(accelerator.z)), kdims=[dim_z], vdims=[dim_Bz])
z_r = hv.Area(((accelerator.z,simulation.envelope_x(accelerator.z)*1e3)), kdims=[dim_z],
vdims=[dim_r], group='Beam')
(z_r+z_Bz).cols(1)
[32]:
As you can see, the envelope is far from perfect.
Apply the genetic algorithm.
Deap¶
First we construct a model envelope.
[12]:
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')
[13]:
z_env
[13]:
Connect the necessary libraries.
[14]:
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.
[15]:
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.
[16]:
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 = 50 # Number of generations
POP = 50 # Number of individuals
[17]:
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.
[18]:
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),
[19]:
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.
[20]:
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()
[21]:
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 50 0.0373014 0 0.0411582 0.0344921 50 0.00155064
1 38 0.0348281 1 0.0375636 0.0327327 38 0.000898671
2 43 0.0334219 2 0.0354142 0.0296962 43 0.000987574
3 40 0.0317281 3 0.0346163 0.0293286 40 0.00120335
4 42 0.0295585 4 0.0319368 0.0273093 42 0.000888509
5 35 0.0279077 5 0.0309627 0.0251781 35 0.000986704
6 35 0.0262448 6 0.0288516 0.0246691 35 0.000886848
7 38 0.0248704 7 0.0267406 0.0239342 38 0.000504496
8 37 0.0241506 8 0.0255704 0.023129 37 0.000552539
9 38 0.0233159 9 0.0250388 0.0222832 38 0.000536844
10 39 0.0225669 10 0.023799 0.0214839 39 0.000552796
11 40 0.0215898 11 0.0230274 0.020334 40 0.000471506
12 35 0.0207404 12 0.0225252 0.0198564 35 0.000507565
13 38 0.0202193 13 0.0214259 0.0194126 38 0.000458111
14 41 0.0196056 14 0.0206291 0.0185125 41 0.000421382
15 38 0.0189491 15 0.0205539 0.0171812 38 0.000492162
16 35 0.0180715 16 0.0194843 0.0170724 35 0.000519289
17 39 0.0171876 17 0.0179657 0.0168041 39 0.000176087
18 42 0.0169952 18 0.0175188 0.0166024 42 0.000203452
19 35 0.0167067 19 0.0169669 0.0164402 35 0.000117499
20 37 0.0165706 20 0.0171469 0.0162325 37 0.000148736
21 37 0.0164005 21 0.0172883 0.0160788 37 0.000218665
22 39 0.0162667 22 0.0172115 0.0159846 39 0.00020546
23 38 0.0161243 23 0.0165658 0.0158296 38 0.000175362
24 34 0.01599 24 0.016868 0.0154732 34 0.00021438
25 35 0.0157555 25 0.0162781 0.0152689 35 0.000235593
26 35 0.0155353 26 0.0168982 0.0152536 35 0.000299995
27 38 0.0153809 27 0.0159431 0.0149969 38 0.000197816
28 39 0.01522 28 0.01601 0.0148777 39 0.000250282
29 40 0.0150565 29 0.0158106 0.014793 40 0.000202265
30 37 0.0149094 30 0.0154494 0.0146449 37 0.000170456
31 40 0.0148635 31 0.0159437 0.0146449 40 0.000212915
32 38 0.0146968 32 0.0151042 0.0143436 38 0.000157724
33 34 0.0146054 33 0.0153208 0.0142456 34 0.000255638
34 42 0.0143601 34 0.0151883 0.0140128 42 0.000229098
35 46 0.0142238 35 0.0148054 0.0139793 46 0.000192102
36 38 0.0141098 36 0.0146636 0.0138423 38 0.000158518
37 35 0.013988 37 0.0163101 0.013766 35 0.000366857
38 41 0.0139639 38 0.0149201 0.0135636 41 0.000268292
39 38 0.0138088 39 0.0147568 0.0135602 38 0.000252177
40 34 0.0137046 40 0.0154014 0.0134745 34 0.000319706
41 41 0.0136914 41 0.0154436 0.0134605 41 0.000330557
42 36 0.0136462 42 0.0151335 0.0132983 36 0.000335113
43 35 0.0135327 43 0.0150908 0.013235 35 0.000372443
44 41 0.0133991 44 0.0153258 0.0130282 41 0.000328673
45 42 0.0133896 45 0.0154806 0.0129309 42 0.00045044
46 38 0.0131798 46 0.014134 0.0128621 38 0.000253342
47 38 0.0131567 47 0.0144235 0.0127197 38 0.000373819
48 38 0.0130646 48 0.0146853 0.0126411 38 0.000410193
49 36 0.0129632 49 0.0155121 0.0125591 36 0.000601679
50 36 0.0128911 50 0.0145057 0.0125462 36 0.000479437
[22]:
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"))
[23]:
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)
[23]:
[ ]: