diff --git a/mbtrack2/tracking/__init__.py b/mbtrack2/tracking/__init__.py index 2e53921b072f942afcde03e87a5e83a6eb086b41..6f72151ebe7cf8ce7a8e72a12395385fc6d08fc1 100644 --- a/mbtrack2/tracking/__init__.py +++ b/mbtrack2/tracking/__init__.py @@ -8,7 +8,9 @@ from mbtrack2.tracking.synchrotron import Synchrotron from mbtrack2.tracking.rf import (RFCavity, CavityResonator, ProportionalLoop, - TunerLoop) + TunerLoop, + ProportionalIntegralLoop, + DirectFeedback) from mbtrack2.tracking.parallel import Mpi from mbtrack2.tracking.element import (Element, LongitudinalMap, diff --git a/mbtrack2/tracking/rf.py b/mbtrack2/tracking/rf.py index ad397fff0b3697bda90e492f289a13994448b9b9..62a8b36d5311b99d743d52c6fb1c37d180c7a219 100644 --- a/mbtrack2/tracking/rf.py +++ b/mbtrack2/tracking/rf.py @@ -325,7 +325,7 @@ class CavityResonator(): # apply different kind of RF feedback for fb in self.feedback: - fb.track(self) + fb.track() self.nturn += 1 @@ -907,15 +907,15 @@ class ProportionalLoop(): """ Proportional feedback loop to control a CavityResonator amplitude and phase. - Feedback setpoints are cavity_resonator.Vc and cavity_resonator.theta. + Feedback setpoints are cav_res.Vc and cav_res.theta. The loop must be added to the CavityResonator object using: - cavity_resonator.feedback.append(loop) + cav_res.feedback.append(loop) Parameters ---------- ring : Synchrotron object - cavity_resonator : CavityResonator + cav_res : CavityResonator The CavityResonator which is to be controlled. gain_A : float Amplitude (voltage) gain of the feedback. @@ -925,15 +925,16 @@ class ProportionalLoop(): Feedback delay in unit of turns. """ - def __init__(self, ring, cavity_resonator, gain_A, gain_P, delay): + def __init__(self, ring, cav_res, gain_A, gain_P, delay): self.ring = ring + self.cav_res = cav_res self.gain_A = gain_A self.gain_P = gain_P self.delay = int(delay) - self.volt_delay = np.ones(self.delay)*cavity_resonator.Vc - self.phase_delay = np.ones(self.delay)*cavity_resonator.theta + self.volt_delay = np.ones(self.delay)*self.cav_res.Vc + self.phase_delay = np.ones(self.delay)*self.cav_res.theta - def track(self, cavity_resonator): + def track(self): """ Tracking method for the amplitude and phase loop. @@ -942,15 +943,15 @@ class ProportionalLoop(): None. """ - diff_A = self.volt_delay[-1] - cavity_resonator.Vc - diff_P = self.phase_delay[-1] - cavity_resonator.theta - cavity_resonator.Vg -= self.gain_A*diff_A - cavity_resonator.theta_g -= self.gain_P*diff_P - cavity_resonator.generator_phasor_record = np.ones(self.ring.h)*cavity_resonator.generator_phasor + diff_A = self.volt_delay[-1] - self.cav_res.Vc + diff_P = self.phase_delay[-1] - self.cav_res.theta + self.cav_res.Vg -= self.gain_A*diff_A + self.cav_res.theta_g -= self.gain_P*diff_P + self.cav_res.generator_phasor_record = np.ones(self.ring.h)*self.cav_res.generator_phasor self.volt_delay = np.roll(self.volt_delay, 1) self.phase_delay = np.roll(self.phase_delay, 1) - self.volt_delay[0] = cavity_resonator.cavity_voltage - self.phase_delay[0] = cavity_resonator.cavity_phase + self.volt_delay[0] = self.cav_res.cavity_voltage + self.phase_delay[0] = self.cav_res.cavity_phase class TunerLoop(): """ @@ -960,12 +961,12 @@ class TunerLoop(): Only a proportional controller is assumed. The loop must be added to the CavityResonator object using: - cavity_resonator.feedback.append(loop) + cav_res.feedback.append(loop) Parameters ---------- ring : Synchrotron object - cavity_resonator : CavityResonator + cav_res : CavityResonator The CavityResonator which is to be controlled. gain : float Proportional gain of the tuner loop. @@ -981,11 +982,12 @@ class TunerLoop(): Tuning offset in [rad]. """ - def __init__(self, ring, cavity_resonator, gain=0.01, avering_period=0, + def __init__(self, ring, cav_res, gain=0.01, avering_period=0, offset=0): self.ring = ring + self.cav_res = cav_res if avering_period == 0: - fs = self.ring.synchrotron_tune(cavity_resonator.Vc)*self.ring.f1/self.ring.h + fs = self.ring.synchrotron_tune(self.cav_res.Vc)*self.ring.f1/self.ring.h avering_period = 2/fs/self.ring.T0 self.Pgain = gain @@ -994,7 +996,7 @@ class TunerLoop(): self.diff = 0 self.count = 0 - def track(self, cavity_resonator): + def track(self): """ Tracking method for the tuner loop. @@ -1005,9 +1007,501 @@ class TunerLoop(): """ if self.count == self.avering_period: diff = self.diff/self.avering_period - self.offset - cavity_resonator.psi -= diff * self.Pgain + self.cav_res.psi -= diff * self.Pgain self.count = 0 self.diff = 0 else: - self.diff += cavity_resonator.cavity_phase - cavity_resonator.theta_g + cavity_resonator.psi + self.diff += self.cav_res.cavity_phase - self.cav_res.theta_g + self.cav_res.psi self.count += 1 + +class ProportionalIntegralLoop(): + """ + Proportional Integral (PI) loop to control a CavityResonator amplitude and + phase via generator current (ig) [1,2]. + + The ProportionalIntegralLoop should be initialized only after generator + parameters are set. + + The loop must be added to the CavityResonator object using: + cav_res.feedback.append(loop) + + Parameters + ---------- + ring : Synchrotron object + cav_res : CavityResonator object + CavityResonator in which the loop will be added. + gain : float or float list + Pgain and Igain of the feedback system + For IQ feedback, same gain set is used for I and Q. + In case of normal conducting cavity (QL~1e4), + the Pgain of ~1.0 and Igain of ~1e4(5) are usually used. + In case of super conducting cavity (QL > 1e6), + the Pgain of ~100 can be used. + In a "bad" parameter set, unstable oscillation of Vc is caused. + So, parameter scan of Vcloop should be made. + sample_num : int + Sample number to monitor Vc + The averaged Vc value in sample_num is monitored. + Units are in bucket numbers. + every : int + Interval to monitor and change Vc + Units are in bucket numbers. + delay : int + Loop delay of the PI feedback system. + Units are in bucket numbers. + Vref : array or complex, optional + Target (set point) value of the feedback. + If None, cav_res.Vc and cav_res.theta are set as the reference. + Default is None. + FF : bool, optional + True is recommended to prevent a Vc drop in the beginning of the tracking. + In case of small Pgain (QL ~ 1e4), Vc drop may cause baem loss due to static Robinson. + Default is True. + matrix : bool, optional + If True, Ig2Vg_matrix is used in Ig2Vg. + Ig2Vg_matrix is faster but init_Ig2Vg_matrix must be called after each + parameter change of cav_res. + If False, Ig2Vg_formula is used in Ig2Vg. + Default is False. + + Methods + ------- + track() + Tracking method for the Cavity PI control feedback. + Ig2Vg_formula() + Return Vg from Ig. + init_Ig2Vg_matrix() + Initialize matrix for Ig2Vg_matrix. + Ig2Vg_matrix() + Return Vg from Ig using matrix formalism. + Ig2Vg() + Go from Ig to Vg and apply values. + Vg2Ig(Vg) + Return Ig from Vg (assuming constant Vg). + + Notes + ----- + Assumption : delay >> every ~ sample_num + + Adjusting ig(Vg) parameter to keep the Vc constant (to target(Vref) values). + The "track" method is + 1) Monitoring the Vc phasor + Mean Vc value between specified bunch number (sample_num) is monitored + with specified interval period (every). + 2) Changing the ig phasor + The ig phasor is changed according the difference of the specified + reference values (Vref) with specified gain (gain). + By using ig instead of Vg, the cavity response can be taken account. + 3) ig changes are reflected to Vg after the specifed delay (delay) of the system + + Vc,IQ-->(rot)-->(-)-->(V->I,fac)-->PI-->Ig,IQ {--> Vg,IQ} + | + Ref + + Examples + -------- + PF; E0=2.5GeV, C=187m, frf=500MHz + QL=11800, fs0=23kHz + ==> gain=[0.5,1e4], sample_num=8, every=7(13ns), delay=500(1us) + is one reasonable parameter set. + + References + ---------- + [1] : https://en.wikipedia.org/wiki/PID_controller + [2] : Yamamoto, N., Takahashi, T., & Sakanaka, S. (2018). Reduction and + compensation of the transient beam loading effect in a double rf system + of synchrotron light sources. PRAB, 21(1), 012001. + + """ + def __init__(self, ring, cav_res, gain, sample_num, every, delay, + Vref=None, FF=True, matrix=False): + self.ring = ring + self.cav_res = cav_res + self.Ig2Vg_mat = np.zeros((self.ring.h, self.ring.h), dtype=complex) + self.gain = gain + self.FF = FF + if not self.FF: + self.FFconst = 0 + self.matrix = matrix + + if delay > 0: + self.delay = int(delay) + else: + self.delay = 1 + if every > 0: + self.every = int(every) + else: + self.every = 1 + record_size = int(np.ceil(self.delay/self.every)) + if record_size < 1: + raise ValueError("Bad parameter set : delay or sample_every") + self.sample_num = int(sample_num) + + # init lists for FB process + self.ig_phasor = np.ones(self.ring.h, dtype=complex)*self.Vg2Ig(self.cav_res.generator_phasor) + self.ig_phasor_record = self.ig_phasor + self.vc_previous = np.ones(self.sample_num)*self.cav_res.cavity_phasor + self.diff_record = np.zeros(record_size, dtype=complex) + self.I_record = 0+0j + + if Vref is None: + self.Vref = self.cav_res.Vc*np.exp(1j*self.cav_res.theta) + else: + self.Vref = Vref + + self.sample_list = range(0, self.ring.h, self.every) + + # Pre caclulation for Ig2Vg + if self.matrix: + self.init_Ig2Vg_matrix() + + def track(self, apply_changes=True): + """ + Tracking method for the Cavity PI control feedback. + + Returns + ------- + None. + + """ + vc_list = np.concatenate([self.vc_previous, self.cav_res.cavity_phasor_record]) #This line is slowing down the process. + self.ig_phasor.fill(self.ig_phasor[-1]) + + for index in self.sample_list: + # 2) updating Ig using last item of the list + diff = self.diff_record[-1] - self.FFconst + self.I_record += diff/self.ring.f1 + fb_value = self.gain[0]*diff + self.gain[1]*self.I_record + self.ig_phasor[index:] = self.Vg2Ig(fb_value) + self.FFconst + # Shift the record + self.diff_record = np.roll(self.diff_record, 1) + # 1) recording diff as a first item of the list + mean_vc = np.mean(vc_list[index:self.sample_num+index])*self.Rot + self.diff_record[0] = self.VrefRot - mean_vc + # update sample_list for next turn + self.sample_list = range(index+self.every-self.ring.h,self.ring.h,self.every) + # update vc_previous for next turn + self.vc_previous = self.cav_res.cavity_phasor_record[- self.sample_num:] + + self.ig_phasor_record = self.ig_phasor + + if apply_changes: + self.Ig2Vg() + + @property + def Vref(self): + """Return the target (set point) value of the feedback as a complex.""" + return self._Vref + + @Vref.setter + def Vref(self, value): + """ + Set the target (set point) value of the feedback. + + Vref can be specified as an IQ-complex style or array-like [Amp, Phase]. + If not specified, cav_res.Vc and cav_res.theta are used. + + """ + if isinstance(value,complex): + self._Vref = value + else: + self._Vref = value[0]*np.exp(1j*value[1]) + self.Rot = np.exp(-1j*np.angle(self._Vref)) + self.VrefRot = self._Vref*self.Rot + if self.FF: + self.FFconst = np.mean(self.ig_phasor) + + def Ig2Vg_formula(self): + """ + Return Vg from Ig. + + Eq.26 of ref [2]. + """ + generator_phasor_record = np.zeros_like(self.cav_res.generator_phasor_record, dtype=complex) + for index in range(self.ring.h): + generator_phasor_record[index] = (self.cav_res.generator_phasor_record[index-1] + *np.exp(-1/self.cav_res.filling_time*(1 - 1j*np.tan(self.cav_res.psi))*self.ring.T1) + + self.ig_phasor_record[index]*self.cav_res.loss_factor*self.ring.T1) + return generator_phasor_record + + def init_Ig2Vg_matrix(self): + """ + Initialize matrix for Ig2Vg_matrix. + + Shoud be called before first use of Ig2Vg_matrix and after each cavity + parameter change. + """ + k = np.arange(0, self.ring.h) + self.Ig2Vg_vec = np.exp(-1/self.cav_res.filling_time*(1 - 1j*np.tan(self.cav_res.psi))*self.ring.T1*(k+1)) + tempV = np.exp(-1/self.cav_res.filling_time*self.ring.T1*k*(1 - 1j*np.tan(self.cav_res.psi))) + for idx in np.arange(self.ring.h): + self.Ig2Vg_mat[idx:,idx] = tempV[:self.ring.h-idx] + + def Ig2Vg_matrix(self): + """ + Return Vg from Ig using matrix formalism. + Faster but self.init_Ig2Vg should be called after each CavityResonator + parameter change. + """ + generator_phasor_record = (self.Ig2Vg_vec*self.cav_res.generator_phasor_record[-1] + + self.Ig2Vg_mat.dot(self.ig_phasor_record)*self.cav_res.loss_factor*self.ring.T1) + return generator_phasor_record + + def Ig2Vg(self): + """ + Go from Ig to Vg. + + Apply new values to cav_res.generator_phasor_record, cav_res.Vg and + cav_res.theta_g from ig_phasor_record. + """ + if self.matrix: + self.cav_res.generator_phasor_record = self.Ig2Vg_matrix() + else: + self.cav_res.generator_phasor_record = self.Ig2Vg_formula() + self.cav_res.Vg = np.mean(np.abs(self.cav_res.generator_phasor_record)) + self.cav_res.theta_g = np.mean(np.angle(self.cav_res.generator_phasor_record)) + + def Vg2Ig(self, Vg): + """ + Return Ig from Vg (assuming constant Vg). + + Eq.25 of ref [2] assuming the dVg/dt = 0. + """ + fac = 1/self.cav_res.filling_time/self.cav_res.loss_factor # 1/R_L + ig = fac*Vg*( 1 - 1j*np.tan(self.cav_res.psi) ) + return ig + +class DirectFeedback(ProportionalIntegralLoop): + """ + Direct RF FB (DFB) via generator current (ig) using PI controller [1]. + + The DirectFeedback inherits from ProportionalIntegralLoop so all + mandatory parameters from ProportionalIntegralLoop should be passed as + kwargs when initializing a DirectFeedback object. + + To avoid cavity-beam unmatching (large synchrotron oscilation of beam), + the CavityResonator generator parameters should be set before + initialization. + + Parameters + ---------- + DFB_gain : float + Gain of the DFB. + DFB_phase_shift : float + Phase shift of the DFB. + DFB_sample_num : int, optional + Sample number to monitor Vc for the DFB. + The averaged Vc value in DFB_sample_num is monitored. + Units are in bucket numbers. + If None, value from the ProportionalIntegralLoop is used. + Default is None. + DFB_every : int, optional + Interval to monitor and change Vc for the DFB. + Units are in bucket numbers. + If None, value from the ProportionalIntegralLoop is used. + Default is None. + DFB_delay : int, optional + Loop delay of the DFB. + Units are in bucket numbers. + If None, value from the ProportionalIntegralLoop is used. + Default is None. + + Attributes + ---------- + phase_shift : float + Phase shift applied to DFB, defined as psi - DFB_phase_shift. + DFB_psi: float + Return detuning angle value with Direct RF feedback in [rad]. + DFB_alpha : float + Return the amplitude ratio alpha of the DFB. + DFB_gamma : float + Return the gain gamma of the DFB. + DFB_Rs : float + Return the shunt impedance of the DFB in [ohm]. + + Methods + ------- + DFB_parameter_set(DFB_gain, DFB_phase_shift) + Set DFB gain and phase shift parameters. + track() + Tracking method for the Direct RF feedback. + DFB_Vg() + Return the generator voltage main and DFB components in [V]. + DFB_fs() + Return the modified synchrotron frequency in [Hz]. + + References + ---------- + [1] : Akai, K. (2022). Stability analysis of rf accelerating mode with + feedback loops under heavy beam loading in SuperKEKB. PRAB, 25(10), + 102002. + + """ + def __init__(self, DFB_gain, DFB_phase_shift, DFB_sample_num=None, + DFB_every=None, DFB_delay=None, **kwargs): + super(DirectFeedback, self).__init__(**kwargs) + + if DFB_delay is not None: + self.DFB_delay = int(DFB_delay) + else: + self.DFB_delay = self.delay + + if DFB_sample_num is not None: + self.DFB_sample_num = int(DFB_sample_num) + else: + self.DFB_sample_num = self.sample_num + + if DFB_every is not None: + self.DFB_every = int(DFB_every) + else: + self.DFB_every = self.every + + record_size = int(np.ceil(self.DFB_delay/self.DFB_every)) + if record_size < 1: + raise ValueError("Bad parameter set : delay or sample_every") + + self.DFB_parameter_set(DFB_gain, DFB_phase_shift) + if np.sum(np.abs(self.cav_res.beam_phasor)) == 0: + cavity_phasor = self.cav_res.Vc*np.exp(1j*self.cav_res.theta) + else: + cavity_phasor = np.mean(self.cav_res.cavity_phasor_record) + self.DFB_VcRecord = np.ones(record_size, dtype=complex)*cavity_phasor + self.DFB_vc_previous = np.ones(self.DFB_sample_num, dtype=complex)*cavity_phasor + + self.DFB_sample_list = range(0, self.ring.h, self.DFB_every) + + @property + def DFB_phase_shift(self): + """Return DFB phase shift.""" + return self._DFB_phase_shift + + @DFB_phase_shift.setter + def DFB_phase_shift(self, value): + """Set DFB_phase_shift and phase_shift""" + self._DFB_phase_shift = value + self._phase_shift = self.cav_res.psi - value + + @property + def phase_shift(self): + """ + Return phase shift applied to DFB. + Defined as self.cav_res.psi - self.DFB_phase_shift. + """ + return self._phase_shift + + @property + def DFB_psi(self): + """ + Return detuning angle value with Direct RF feedback in [rad]. + + Fig.4 of ref [1]. + """ + return (np.angle(np.mean(self.cav_res.cavity_phasor_record)) - + np.angle(np.mean(self.ig_phasor_record))) + + @property + def DFB_alpha(self): + """ + Return the amplitude ratio alpha of the DFB. + + Near Eq. (13) of [1]. + """ + fac = np.abs(np.mean(self.DFB_ig_phasor)/np.mean(self.ig_phasor_record)) + return 20*np.log10(fac) + + @property + def DFB_gamma(self): + """ + Return the gain gamma of the DFB. + + Near Eq. (13) of [1]. + """ + fac = np.abs(np.mean(self.DFB_ig_phasor)/np.mean(self.ig_phasor_record)) + return fac/(1-fac) + + @property + def DFB_Rs(self): + """ + Return the shunt impedance of the DFB in [ohm]. + + Eq. (15) of [1]. + """ + return self.cav_res.Rs/(1+self.DFB_gamma*np.cos(self.DFB_psi)) + + def DFB_parameter_set(self, DFB_gain, DFB_phase_shift): + """ + Set DFB gain and phase shift parameters. + + Parameters + ---------- + DFB_gain : float + Gain of the DFB. + DFB_phase_shift : float + Phase shift of the DFB. + + """ + self.DFB_gain = DFB_gain + self.DFB_phase_shift = DFB_phase_shift + + if np.sum(np.abs(self.cav_res.beam_phasor)) == 0: + vc = np.ones(self.ring.h)*self.cav_res.Vc*np.exp(1j*self.cav_res.theta) + else: + vc = self.cav_res.cavity_phasor_record + vg_drf = self.DFB_gain*vc*np.exp(1j*self.phase_shift) + self.DFB_ig_phasor = self.Vg2Ig(vg_drf) + + self.ig_phasor = self.ig_phasor_record - self.DFB_ig_phasor + if self.FF: + self.FFconst = np.mean(self.ig_phasor) + + def track(self): + """ + Tracking method for the Direct RF feedback. + + Returns + ------- + None. + + """ + super(DirectFeedback, self).track(False) + + vc_list = np.concatenate([self.DFB_vc_previous, self.cav_res.cavity_phasor_record]) + self.DFB_ig_phasor = np.roll(self.DFB_ig_phasor, 1) + for index in self.DFB_sample_list: + # 2) updating Ig using last item of the list + vg_drf = self.DFB_gain*self.DFB_VcRecord[-1]*np.exp(1j*self.phase_shift) + self.DFB_ig_phasor[index:] = self.Vg2Ig(vg_drf) + # Shift the record + self.DFB_VcRecord = np.roll(self.DFB_VcRecord, 1) + # 1) recording Vc + mean_vc = np.mean(vc_list[index:self.DFB_sample_num+index]) + self.DFB_VcRecord[0] = mean_vc + # update sample_list for next turn + self.DFB_sample_list = range(index+self.DFB_every-self.ring.h, self.ring.h, self.DFB_every) + # update vc_previous for next turn + self.DFB_vc_previous = self.cav_res.cavity_phasor_record[- self.DFB_sample_num:] + + self.ig_phasor_record = self.ig_phasor + self.DFB_ig_phasor + + self.Ig2Vg() + + def DFB_Vg(self, vc=-1): + """Return the generator voltage main and DFB components in [V].""" + if vc == -1: + vc = np.mean(self.cav_res.cavity_phasor_record) + vg_drf=self.DFB_gain*vc*np.exp(1j*self.phase_shift) + vg_main=np.mean(self.cav_res.generator_phasor_record)-vg_drf + return vg_main, vg_drf + + def DFB_fs(self, vg_main=-1, vg_drf=-1): + """Return the modified synchrotron frequency in [Hz].""" + vc = np.mean(self.cav_res.cavity_phasor_record) + if vg_drf ==-1: + vg_drf=self.DFB_gain*vc*np.exp(1j*self.phase_shift) + if vg_main == -1: + vg_main=np.mean(self.cav_res.generator_phasor_record)-vg_drf + vg_sum = np.abs(vg_main)*np.sin(np.angle(vg_main))+np.abs(vg_drf)*np.sin(np.angle(vg_drf)) + omega_s = 0 + if (vg_sum) > 0.0: + omega_s=np.sqrt(self.ring.ac*self.ring.omega1*(vg_sum)/self.ring.E0/self.ring.T0) + return omega_s/2/np.pi