Example 002: Extrapolation of simulation data

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

Simulation data example

This example features CST simulated wakefield data of the “Accelerator Cavity” introductory example.

Fitting on fully decayed wakefield

Import data

# Importing impedance data
data_fully_decayed = np.loadtxt('examples/data/002_impedance_acceleratorCavity_fully_decayed.txt', comments='#', delimiter='\t')

# Extracting frequency and impedance

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

impedance_fd = real_impedance + 1j*imag_impedance

Visually inspecting the impedance

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

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

<matplotlib.legend.Legend at 0x7f797cbb0f50>

../_images/9fcc983585b0d946f564a2a8ffb37e4109d7b0c0a986cfaeb53cf49b717164a4.png

Fitting resonators with IDDEFIX on the absolute magnitude of the 3 resonator impdance spectrum

# Setting amount of resonators to fit

N_resonators = 2

# 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, np.abs(impedance_fd), minimum_peak_height=1000)

SBD_parameterBound.inspect()

../_images/ac2ade0a57d8ff318d5f7589a068b53bc974a1e81fef51790b7cef896b2567f8.png

Running IDDEFIX DE

Running the DE algorithm with IDDEFIX and chosen parameters.

%%time
DE_model_fd = iddefix.EvolutionaryAlgorithm(frequency, 
                                         impedance_fd,
                                         N_resonators=N_resonators, 
                                         parameterBounds=SBD_parameterBound.find(),
                                         plane="longitudinal", 
                                         objectiveFunction=iddefix.ObjectiveFunctions.sumOfSquaredError)

DE_model_fd.run_differential_evolution(maxiter=30000,
                                    popsize=30,
                                    tol=0.01,
                                    mutation=(0.4, 1.0),
                                    crossover_rate=0.7)
print(DE_model_fd.warning)

Optimization Progress %:   0%|          | 0/100 [00:00<?, ?it/s]

Optimization Progress %: 112.70799498485952it [00:07, 14.96it/s]

----------------------------------------------------------------------
Resonator |   Rs [Ohm/m or Ohm]    |        Q         |    fres [Hz]     
----------------------------------------------------------------------
    1     |        3.04e+03        |      69.97       |    5.453e+08     
    2     |        2.37e+03        |      64.96       |    8.123e+08     
----------------------------------------------------------------------
callback function requested stop early
CPU times: user 3.44 s, sys: 2.13 s, total: 5.57 s
Wall time: 7.58 s

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_fd.run_minimization_algorithm()

Method for minimization : Nelder-Mead


----------------------------------------------------------------------
Resonator |   Rs [Ohm/m or Ohm]    |        Q         |    fres [Hz]     
----------------------------------------------------------------------
    1     |        3.03e+03        |      69.88       |    5.453e+08     
    2     |        2.36e+03        |      64.55       |    8.123e+08     
----------------------------------------------------------------------

Assesing the fitting visually

plt.figure(figsize=(10, 6))
result_DE = np.abs(DE_model_fd.get_impedance(use_minimization=False))
result_DE_MIN = np.abs(DE_model_fd.get_impedance())

plt.plot(frequency, np.abs(impedance_fd), 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 0x7f79a89217d0>

../_images/581b22527ac6bdbc776fcb84541596fd30f334efadb326e53b0190f848fddff9.png

Fitting on partially decayed wakefield