Example 002: Fitting SPS-BWS simulation data

import numpy as np
import matplotlib.pyplot as plt
import iddefix
from scipy.constants import c

Simulation data example

This example features simulated wakefield data of the “Beam Wire Scanner” device.

Fitting on fully decayed wakefield

Import data

# Importing impedance data
data = np.loadtxt('examples/data/003_beam_wire_scanner.txt', comments='#', delimiter='\t')

# Extracting frequency and impedance

frequency = data[:,0]*1e9 # Convert to GHz
real_impedance = data[:,1]
imag_impedance = data[:,2]

impedance = real_impedance + 1j*imag_impedance

plt.plot(frequency, impedance.real, label='Real')
plt.plot(frequency, impedance.imag, label='Imag')
plt.plot(frequency, np.abs(impedance), label='Abs')

plt.xlabel('Frequency [Hz]')
plt.ylabel('Impedance [Ohm]')
plt.legend()

<matplotlib.legend.Legend at 0x7f8631cacad0>

For the beam wire scanner, the power loss due to the real part of the impedance is crucial. Thus this example will demonstrate fitting upon the real part, as opposed to the prior examples which fitted upon the absolute impedance magnitude

# Bounds on resonators parameters
""" In this example, we will use the SmartBoundDetermination class to determine the bounds on the resonators parameters. 
By using SmartBoundDetermination ther is no need to specify the bounds manually."""

SBD_parameterBound = iddefix.SmartBoundDetermination(frequency, impedance.real, minimum_peak_height=2)

parameterBounds = SBD_parameterBound.find()

N_resonators = SBD_parameterBound.N_resonators

SBD_parameterBound.inspect()

Running IDDEFIX DE

Running the DE algorithm with IDDEFIX and chosen parameters.

%%time
DE_model = iddefix.EvolutionaryAlgorithm(frequency, 
                                         impedance.real,
                                         N_resonators=N_resonators, 
                                         parameterBounds=SBD_parameterBound.find(),
                                         plane="longitudinal", 
                                         objectiveFunction=iddefix.ObjectiveFunctions.sumOfSquaredErrorReal)

DE_model.run_differential_evolution(maxiter=30000,
                                    popsize=90,
                                    tol=0.0001,
                                    mutation=(0.3, 0.8),
                                    crossover_rate=0.5)
print(DE_model.warning)

Optimization Progress %: 104.20929085863378it [02:17,  1.32s/it]                                

----------------------------------------------------------------------
Resonator |   Rs [Ohm/m or Ohm]    |        Q         |    fres [Hz]     
----------------------------------------------------------------------
    1     |        8.83e+00        |      36.41       |    1.483e+08     
    2     |        2.78e+00        |      18.71       |    3.618e+08     
    3     |        1.52e+01        |      10.77       |    6.504e+08     
    4     |        1.28e+02        |      10.97       |    7.661e+08     
    5     |        1.13e+02        |      27.00       |    8.155e+08     
    6     |        4.30e+01        |      35.81       |    1.001e+09     
----------------------------------------------------------------------
callback function requested stop early
CPU times: user 1min 29s, sys: 6.6 s, total: 1min 35s
Wall time: 2min 17s

Minimization step

To further refine the solution obtained by the DE algorithm, a second optimization step is applied using the Nelder-Mead minimization algorithm.

DE_model.run_minimization_algorithm(margin=0.5)

Method for minimization : Nelder-Mead


----------------------------------------------------------------------
Resonator |   Rs [Ohm/m or Ohm]    |        Q         |    fres [Hz]     
----------------------------------------------------------------------
    1     |        7.09e+00        |      22.04       |    1.483e+08     
    2     |        2.90e+00        |      25.50       |    3.618e+08     
    3     |        7.62e+00        |       6.66       |    6.292e+08     
    4     |        1.30e+02        |       9.85       |    7.687e+08     
    5     |        9.22e+01        |      28.52       |    8.164e+08     
    6     |        4.30e+01        |      38.29       |    1.001e+09     
----------------------------------------------------------------------

plt.figure(figsize=(10, 6))
result_DE = DE_model.get_impedance(use_minimization=False).real
result_DE_MIN = DE_model.get_impedance().real

plt.plot(frequency, impedance.real, lw=5, label='Real', color='black')
plt.plot(frequency, result_DE,  lw=2, label='DE - Res. fitting')
plt.plot(frequency, result_DE_MIN,  lw=2, label='DE + Min. - Res. fitting')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Impedance [Ohm]')
plt.legend()

<matplotlib.legend.Legend at 0x7f8631adf690>

And the wake too:

plt.figure(figsize=(10, 6))

timeframe = np.linspace(-1, 20, 1000)/c

result_DE = DE_model.get_wake_potential(timeframe, 1e-10, use_minimization=False)/1e12
result_DE_MIN = DE_model.get_wake_potential(timeframe, sigma=1e-10)/1e12

plt.plot(timeframe, result_DE,  lw=2, label='DE - Res. fitting')
plt.plot(timeframe, result_DE_MIN,  lw=2, label='DE + Min. - Res. fitting')
plt.xlabel('Wake length [m]')
plt.ylabel('Wake potential [V/pC]')
plt.legend()

<matplotlib.legend.Legend at 0x7f8631b67750>