Simulations of real-time system identification for superconducting cavities with a recursive least-squares algorithm

We explore the performance of a recursive least-squares algorithm to determine the bandwidth $\omega_{12}$ and the detuning $\Delta\omega$ of a superconducting cavity. We base the simulations on parameters of the ESS double-spoke cavities. Expressions for the signal-to-noise ratio of derived parameters are given to explore the applicability of the algorithm to other configurations.


I. INTRODUCTION
Superconducting accelerating cavities are used to accelerate protons [1,2], electrons [3][4][5], and heavy ions [6][7][8], both with pulsed [1,4] and with continuous beams [3,9].Owing to the low losses, the cavities have a very narrow bandwidth on the order of Hz for bare cavities and a few 100 Hz for cavities equipped with high-power couplers.In order to efficiently cool these cavities with liquid helium they are made of rather thin material, which makes them easily deformable and this changes their resonance frequency, often by an amount comparable to their bandwidth.In pulsed operation, the dominant deformation comes from the electro-magnetic pressure of the field inside the cavity, the Lorentz-force detuning [10,11], while cavities operating continuously are perturbed by so-called microphonics [12,13], caused by pressure variations of the liquid helium bath or mechanical perturbations, for example, by reciprocating pumps or by malfunctioning equipment.As a consequence of these perturbations, the cavities are detuned and force the power generators to increase their output to maintain fields necessary for stable operation of the beams.This reduces the efficiency of the system and requires an, often substantial, overhead of the power generation, forcing it to operate at a less than optimal working point.To avoid this sub-optimal mode of operation and to compensate the detuning, many accelerators employ active tuning systems that use stepper motors and piezo-actuators [14] to squeeze the cavities back in tune, which requires diagnostic systems to measure the detuning.These measurements are usually based arXiv:2309.01653v2[physics.acc-ph]20 Nov 2023 on comparing the phase of the signal that excites the cavity, measured with a directional coupler just upstream of the input coupler, to the phase of the field inside the cavity, measured by a field probe or antenna inside the cavity.Both analog [12,15] and digital [16,17] signal processing systems are used; often as part of the low-level radio-frequency (LLRF) feedback system that stabilizes the fields in the cavity.Even more elaborate systems, based on various system identification algorithms, are used or planned [18][19][20][21].All these algorithms normally rely on low-pass filtering the often noisy signals from the directional couplers and antennas in order to provide a reliable estimate of the cavity detuning and the bandwidth.
In this report, we focus on a complementary algorithm that continuously improves the estimated fit parameters by increasing the size of a system of equations.Instead of solving this rapidly growing system directly, we employ a recursive least-squares (RLS) algorithm [22,23], which only requires moderate numerical expenditure in each time step.Remarkably, asymptotically the difference between the continuously improving estimates of the fit parameters and the "true" values-the so-called estimation error-approaches zero [24] albeit at the expense of a limited ability to resolve changing parameters.We therefore introduce a finite memory when solving the system, which downgrades old measurements in favor of new ones.This allows us to handle even changing parameters at the expense of an increased noise level of the fit parameters.
In the following sections, we first introduce the model of the cavity and transform the continuous-time model to discrete time.In Section III we develop the RLS algorithm to identify the cavity parameters.In Section IV we explore the capabilities of the algorithm in simulations before calculating the signal-to-noise ratio in Section V and the conclusions.

II. MODEL
Accelerating cavities can be described by an equivalent circuit composed of a resistor R, an inductance L, and a capacitor C, all connected in parallel.This circuit is then excited by a current I and responds by a building up a voltage V = V r +iV i across the components.This voltage is decomposed into real (in-phase, I) and imaginary (out-of-phase, Q) components.
After averaging over the fast oscillations the evolution of the real and imaginary parts of the voltage envelope is given by the following state-space representation [16] of the system that describes the dynamics of the cavity voltage powered by a generator that provides the currents.The directional couplers used to measure the input signal, however, measure the forward component of the current ⃗ I + rather than the total current ⃗ I = ⃗ I + + ⃗ I − .
Close to resonance, it is straightforward to show that the measured forward current ⃗ I + , which proportional to the signal from the directional coupler, is related to the total current ⃗ I by with the coupling factor β = Q 0 /Q E given by the ratio of the intrinsic quality factor of the cavity Q 0 and the external quality factor defines the loaded quality factor Q L .Replacing the currents on the right-hand side of Equation 1 with the help of Equation 2 then leads to with ω E = ω/Q E .We also introduce ω 12 = ω/2Q L and the cavity resonance frequency ω.
Furthermore, we assume that magnitude and phase of all currents and voltages can be reliably measured after the hardware (antennas, cables, and amplifiers) is properly calibrated.
Equation 3 is in the standard form of a linear dynamical system ⃗ V = Ā⃗ V + B⃗ I + where ⃗ V is the column vector with real and imaginary part of the voltages and ⃗ I + that of the forward currents.The matrices Ā and B correspond to those in Equation 3 and are given by For the simulations we will convert the continuous-time system from Equation 3 to discrete time with time step ∆t, which corresponds to the sampling time if the system is implemented digitally.By replacing the derivatives of the voltages by finite differences where we label the time steps by t, Equation 3 becomes B = ω E ∆tR1, and the process noise ⃗ w t .We assume that the noise is uncorrelated and has magnitude σ p .It is thus characterized by its expectation value in the system identification process.We assume it is uncorrelated, has magnitude σ m , and

III. SYSTEM IDENTIFICATION
Now we turn to the task of extracting ω 12 ∆t and ∆ω∆t from continuously measured voltages ⃗ V ′ t and currents ⃗ I + t .In order to isolate the sought parameters, we rewrite Equation 6 in the form and B = ω E ∆tR1.After reorganizing this equation to we rewrite F ⃗ V ′ t on the right-hand side as We now introduce the abbreviations and stack Equation 9for consecutive times on top of each other.In this way, we obtain a growing system of equations to determine ω 12 ∆t and ∆ω∆t that we solve in the least-squares sense with the Moore-Penrose pseudo-inverse [26] ⃗ Here we introduce the abbreviation ⃗ q T to denote the estimated parameters at time step T .
We can avoid lengthy evaluations by calculating Equation 13recursively.With the definition P −1 T = U ⊤ T U T , its initial value P 0 = p 0 1, and the definition of U T from Equation 12we express P T +1 through P T in the following way We note that for all time steps t is proportional to the unit matrix 1.This renders the fit into two orthogonal and independent parts; one for each of the fit parameters.To proceed, we introduce the scalar quantity p T with P T = p T 1 and find that it obeys Taking the reciprocal leads to Note that we need to initialize this recursion with a non-zero value and set p 0 = 1 in the simulations.Despite being numerically unity, we carry p 0 through all equations, because it carries the inverse units of ⃗ V 2 T .We now turn to finding ⃗ q T +1 by writing Equation 13 for T + 1 = 1 Equations 17 and 18 constitute the algorithm to continuously update estimates for the two components of ⃗ q, the bandwidth q(1) = ω 12 ∆t and the detuning q(2) = ∆ω∆t, as new voltage and current measurements-both enter in G T +1 and ⃗ y T +2 -become available.We refer to the MATLAB [27] code on github [28] for the details of the implementation.
In Equations 17 and 18 new information from measurements are used to continuously improve the estimate of the fit parameters, but in situations where they change, we have to introduce a way to forget old information.Therefore, in order to emphasize newly added information we follow [22,29] and introduce a "forgetting factor" α = 1 − 1/N f where N f is the time horizon over which old information is downgraded in the last equality of Equation 14, which now reads We see that we only have to replace P T by P T /α, or equivalently p T by p T /α, in the derivation of Equations 17 and 18 and find for the update of p T and for the update of the estimated parameters ⃗ q T ⃗ q T +1 = 1 that are capable of following time-dependent system parameters.These expressions can be evaluated very efficiently.We find that the calculations in Equation 20 involve four multiplications and one inverse whereas the calculations in Equation 21involve ten multiplications if we reuse the expression in the square bracket.Thus, in total fourteen multiplication and one, computationally more expensive, inverse are required.This is about ten times the computational effort needed for a PI controller that typically requires three multiplications.The processing delay of the system identification algorithm should therefore be correspondingly longer.The details of the timing depend of course on the hardware used to implement these algorithms.In particular, on a field-programmable gate array, many operations can be done in parallel.

IV. SIMULATIONS
We base our simulations on parameters for the prototype spoke-cavity module [30] for the ESS [2], which operates at 352 MHz, has an external Q E [31] in the range 1.75 × 10 5 to 2.85 × 10 5 .One of the measured cavities exhibited a loaded-Q of Q L = 1.8 × 10 5 [32] while it was operating at a high gradient.The resulting bandwidth is f 12 = ω 12 /2π ≈ 1000 Hz.
The cavity showed Lorentz-force detuning on the order of a few hundred Hz [32][33][34]; for our simulations we typically use ∆f = ∆ω/2π = 500 Hz.Moreover, we use a process noise level of σ p = 10 −4 × V max and a measurement noise level of σ m = 10 −3 × V max , where V max is the peak voltage inside the cavity.We report the voltages and currents normalized to the values without detuning and denote them by v r , v i and i r , i i respectively.The peak voltage and current in those conditions then becomes unity.Furthermore, we assume that the data-acquisition system operates at a rate of 10 Msamples/s, resulting in ∆t = 100 ns.
We found that the forgetting horizon N f scales with the relative noise levels σ p and σ m .We use N f = 100, unless explicitely specified, because it gave good results.
The left-hand side in Figure 1 shows the normalized currents and voltages over the first 1000 iterations (100 µs), where the currents are turned on after 100 iterations.We observe that the real part of the current (black line) assumes its new value at that point whereas We can understand this behavior by noting that p T is proportional to the diagonal element of the empirical covariance matrix of the least-squares fit in Equation 13.
Therefore the square root of p T is proportional to the error bars of the fit parameter.Figure 2 shows p T for a simulation with N f = 100 (black solid) and N f = 10 (red dashes) for 5000 iterations.We observe that both curves initially increase during the period that the fit is noisy but then approach a constant value that determines the achievable error bars of the fit parameters.This value can be derived from Equation 20 by setting p T +1 = p T = p ∞ and solving for p Here ⃗ V ′ ∞ is the voltage inside the cavity.For the error bars of both components of ⃗ q we thus find σ m / N f ⃗ V ′2 ∞ ), a value that corresponds to the rms deviations of the fit parameters, shown, for example on the second half in Figure 1.Furthermore by construction, the off-diagonal elements of the matrix P T = p T 1 are zero, which indicates that the fit of the bandwidth and the detuning are orthogonal and that makes the algorithm very robust.Moreover, we found that instead of operating open-loop, using a PI-controller to control the cavity voltage does not significantly alter the performance of the system identification process.
We now explore the algorithm's ability to identify parameter changes during steady state operation.The left-hand side in Figure 3 illustrates the effect of microphonics on the the currents and voltages.We simulate this by an oscillation of ∆f with amplitude of f 12 /2 and frequency 1 kHz.Especially v i reveals this oscillation, though also v r oscillates.The right-hand side of Figure 3 shows how the algorithm correctly identifies f 12 and both the amplitude and oscillation frequency of ∆f .
Increasing the oscillation frequency to 20 kHz results in Figure 4 where we have reduced the duration of the simulation to 10 4 iterations in order to improve the visibility of oscillations on the plot.We see that the oscillations are still resolved, albeit at a lower amplitude, which is a consequence of the forgetting horizon N f = 100.It implicitly introduces averaging over N f iterations and thus behaves like a low-pass filter with a time constant of N f ∆t = 10 µs or a cutoff frequency on the order of 100 kHz that already causes some attenuation of the The oscillations are still seen, but the amplitude is significantly reduced.This can be partially alleviated by decreasing N f , albeit at the expense of an increased noise level.

kHz oscillation.
In Figure 5 we explore a rapid increase of the bandwidth, for example, due to a quench.
In the simulation, we simply double the value of ω 12 after 5000 iterations.The plots in the top-left of Figure 5 show the currents and voltages and on the top-right the fit parameters.
We find that the fitted bandwidth (black) is indeed doubled and that the reconstruction of the detuning is unaffected.The plot on the bottom left shows an enlarged view of the fit parameters around the time of the step.It shows that the doubled value is approached within about 2 × N f = 200 iterations.If we run the same simulation with a ten times reduced value of N f = 10, we obtain the plot on the bottom right.We find that the changed value is approached within a few tens of iterations, albeit at the expense of an increased noise level, which is consistent with the discussion regarding Figure 2. Balancing the noise level and the response is just a matter of adjusting the value of N f , the topic of the following section.

V. SIGNAL TO NOISE
In Section IV we already found that the asymptotic noise level N for constant parameters is given by where we denote the magnitude of ⃗ V ′ ∞ by V ′ ∞ .We now consider a situation where the system has reached a quasi-stationary state and that perturbations of the ω 12 and ∆ω are so small that they affect V ′ ∞ very little.We can therefore also use it to write p ∞ = 1/N f V ′2 ∞ despite temporally varying ω 12 and ∆ω.Replacing p T by p ∞ in Equation 21 then leads to Using Equation 9and 10 we rewrite ⃗ y T +2 as where the vector on the right-hand side with the subscript hw are the "true" values of the hardware.Combining these equations, utilizing Equation 15, and replacing V ′ T by V ′ ∞ we arrive at In the next step we use α = 1 − 1/N f and reshuffle terms to obtain Introducing τ f = N f ∆t, replacing the finite difference by a differential, and Laplacetransforming the resulting equation we find where s is the Laplace variable and we denote the Laplace transform of a variable by a tilde.
We obtain the the time dependence by replacing s = iω = 2πif (28) and find that the reconstructed system parameters ⃗ q are given by the hardware parameters passed through a low-pass filter with time constant τ f .Of particular interest is the absolute value of the amplitude of the detuning ∆ω at frequency ω, which is given by This constitutes the signal we strive to measure.For the signal-to-noise ratio S/N we then find where all parameters are explicitely written out in order to explore the trade-off among them.Apparently it depends on the magnitude (amplitude) of the detuning ∆ω hw and the relative accuracy of the voltage measurement σ m /V ′ ∞ , but also on the attenuation of an oscillation due to the forgetting time horizon N f .As long as S/N is sufficiently large, say 5 or so, the oscillation is discernible.

VI. CONCLUSIONS
We worked out an algorithm to determine the cavity bandwidth f 12 and the detuning ∆f by correlating signal from a directional coupler before the cavity and the voltages inside the cavity.The calculations are very efficient and given by Equations 17 and 18 for static parameters and by Equations 20 and 21 for time-varying parameters.These recursion equations are very compact and require only moderate resources, for example, on a field-programmable gate array.
Despite the absence of low-pass filtering, the RLS algorithm is resilient to noise of the measured voltages, because the forgetting horizon implicitly introduces a low-pass filter whose time constant is τ f = N f ∆t.We can taylor the performance by selecting a large value of N f , which reduces the noise of the reconstructed parameters, whereas smaller values of N f make the algorithm more responsive to parameter changes on faster times scales.The trade-off between achievable frequency resolution, N f , and measurement noise σ m can be explored with the help of Equation 30.

FIG. 1 :
FIG. 1: Left: the normalized currents (bottom) and the voltages (top) after starting to fill the cavity for 1000 iterations (100 µs).Right: the reconstructed fit parameters, the bandwidth f 12 (black) and the detuning ∆f (red).Note that the parameters are found despite the noise level (σ p = 10 −4 and σ m = 10 −3 of peak voltage) used in the simulation.

FIG. 2 :
FIG.2:The variable p T as a function of the iterations for N f = 100 (solid) and N f = 10 (dashes).

FIG. 3 :
FIG. 3: The currents and voltages (left) and the fit parameters (right) for 10 5 iterations (10 ms) while the detuning ∆f oscillates with an amplitude of 500 Hz and with a mechanical-mode frequency of 1 kHz.The oscillations are clearly visible on both phases of the voltage and the correctly reconstructed fit parameters.

FIG. 4 :
FIG. 4: The reconstructed fit parameters for a 20 kHz mechanical oscillation of the detuning ∆f .

FIG. 5 :
FIG. 5: The currents and voltages (top left) and fit parameters (top right) for 10 4 iterations (1 ms) as the bandwidth f 12 is doubled at iteration 5000.The bottom row shows an enlarged view of fit parameters around the time of the change.On the left, we use N f = 100 and on the right we use N f = 10.