Example 001: Analytical Resonator fitting
import numpy as np
import matplotlib.pyplot as plt
import iddefix
Modelling a three resonator impedance spectrum by three arbitrary resonators
Any arbitrary number of resonators can be initialized. This is example will demonstrate the creation of three resonators with the resonator Formulas functions inside iddeficx and later fitting.
Three random resonators (\(R_s\), Q, \(f_r\)):
400 \(\Omega\), 30, 0.2 GHz
1000 \(\Omega\), 10, 1 GHz
500 \(\Omega\), 20, 1.75 GHz
# Assigning the resonator parameters
parameters ={ '1': [400, 30, 0.2e9],
'2': [1000, 10, 1e9],
'3': [500, 20, 1.75e9],}
# Computing the impedance spectrum for the resonators
frequency = np.linspace(0, 2e9, 1000)
impedance = iddefix.Impedances.n_Resonator_longitudinal_imp(frequency, parameters)
Plotting the impedance spectrum
Plots of both the real- and imaginary part of the impedance, as well as the absolute of the complex 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 0x7f43de3cd590>
Fitting resonators with IDDEFIX on the absolute magnitude of the 3 resonator impdance spectrum
# Setting amount of resonators to fit
N_resonators = 3
# Bounds on resonators parameters
""" Bounds have this format [(Rs_min, Rs_max), (Q_min, Q_max), (fres_min, fres_max)].
ParameterBounds allows us to manually add a resonator with desired parameters """
parameterBounds = [(0, 2000), (1, 1e3), (0.1e9, 2e9),
(0, 2000), (1, 1e3), (0.1e9, 2e9),
(0, 2000), (1, 1e3), (0.1e9, 2e9)]
Running IDDEFIX DE
Running the DE algorithm with IDDEFIX and chosen parameters. A good rule-of-thumb is to have the population set at N = 5 * Nres. Mutation parameters and crossover_rate can be changed to adjust level of exploration/exploitation.
%%time
DE_model = iddefix.EvolutionaryAlgorithm(frequency,
impedance,
N_resonators=N_resonators,
parameterBounds=parameterBounds,
plane="longitudinal",
objectiveFunction=iddefix.ObjectiveFunctions.sumOfSquaredError)
DE_model.run_differential_evolution(maxiter=400,
popsize=45,
tol=0.01,
mutation=(0.4, 1.0),
crossover_rate=0.7)
print(DE_model.warning)
Optimization Progress %: 3%|▎ | 2.6286408213164307/100 [00:23<14:41, 9.05s/it]
----------------------------------------------------------------------
Resonator | Rs [Ohm/m or Ohm] | Q | fres [Hz]
----------------------------------------------------------------------
1 | 1.00e+03 | 9.84 | 1.000e+09
2 | 5.03e+02 | 19.97 | 1.750e+09
3 | 4.01e+02 | 30.64 | 1.998e+08
----------------------------------------------------------------------
Maximum number of iterations has been exceeded.
CPU times: user 15.7 s, sys: 7.29 s, total: 23 s
Wall time: 23.8 s
The found resonator parameters are close to being exactly correct. We can fit better by running for more generations or by doing the mimimization step:
Minimization step
To further refine the solution obtained by the DE algorithm, a second optimization step is applied using the Nelder-Mead minimization algorithm. This additional step starts with the results from the DE algorithm as the initial guess and iteratively adjusts the parameters within a predefined range of 10% above or below their original values. By doing so, the Nelder-Mead algorithm fine-tunes the solution to reduce the error further, leveraging its capability to explore the local parameter space efficiently. This two-step optimization ap- proach ensures a more precise fit by combining the global search power of DE with the local refinement capabilities of Nelder-Mead.
DE_model.run_minimization_algorithm()
Method for minimization : Nelder-Mead
----------------------------------------------------------------------
Resonator | Rs [Ohm/m or Ohm] | Q | fres [Hz]
----------------------------------------------------------------------
1 | 1.00e+03 | 10.00 | 1.000e+09
2 | 5.00e+02 | 20.00 | 1.750e+09
3 | 4.00e+02 | 30.00 | 2.000e+08
----------------------------------------------------------------------
Both the minimization and the differential evolution routines return the parameters in a list size [3*Nres]
print('Minimization parameters:')
print(DE_model.minimizationParameters)
print('Diferential evolution parameters:')
print(DE_model.evolutionParameters)
Minimization parameters:
[9.99999993e+02 9.99999983e+00 1.00000000e+09 4.99999992e+02
1.99999994e+01 1.75000000e+09 4.00000007e+02 3.00000019e+01
2.00000000e+08]
Diferential evolution parameters:
[9.99551882e+02 9.84124038e+00 9.99983157e+08 5.03335415e+02
1.99657188e+01 1.75026584e+09 4.01449425e+02 3.06411196e+01
1.99815069e+08]
There’s also a method to print them in Markdown syntax
print('For terminal:')
DE_model.display_resonator_parameters(DE_model.minimizationParameters)
print('\nFor Markdown:')
DE_model.display_resonator_parameters(DE_model.minimizationParameters, to_markdown=True)
For terminal:
----------------------------------------------------------------------
Resonator | Rs [Ohm/m or Ohm] | Q | fres [Hz]
----------------------------------------------------------------------
1 | 1.00e+03 | 10.00 | 1.000e+09
2 | 5.00e+02 | 20.00 | 1.750e+09
3 | 4.00e+02 | 30.00 | 2.000e+08
----------------------------------------------------------------------
For Markdown:
| Resonator | Rs [Ohm/m or Ohm] | Q | fres [Hz] |
|-----------|------------------|---|-----------|
| 1 | 1000 | 10 | 1e+09 |
| 2 | 500 | 20 | 1.75e+09 |
| 3 | 400 | 30 | 2e+08 |
Assesing the fitting visually
plt.figure(figsize=(10, 6))
result_DE = np.abs(DE_model.get_impedance(use_minimization=False))
result_DE_MIN = np.abs(DE_model.get_impedance())
plt.plot(frequency, np.abs(impedance), 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 0x7f43df0ef4d0>
Retrieve the wake function and wake potential
Using iddefix Resonator formalism formulas, derived by S. Joly, one can reconstruct the wake function and the wake potential using the analytical expressions available inside resonatorFormalism.py and the parameters
plt.figure(figsize=(10, 6))
time = np.linspace(-1e-9, 20e-9, 100000)
wake = DE_model.get_wake(time_data=time)
wake_potential = DE_model.get_wake_potential(time_data=time, sigma=1e-10)
plt.plot(time*1e9, wake, lw=2, label='Wake function', color='blue')
plt.plot(time*1e9, wake_potential, lw=2, label='Wake potential', color='red')
plt.xlabel('Time [ns]')
plt.ylabel('Wake potential & wake function')
plt.legend()
<matplotlib.legend.Legend at 0x7f441fe79390>
