Particle swarming of sensor correction filters

Figure 1 shows a functional schematic of the SC signal path. Each block represents a type of sensor installed in the LIGO observatories and shows the contributions to their output signals in terms of inertial translation, tilt (or rotation), and self-noise. All variables carry a tilde to indicate the Fourier transform. From this block diagram, the sensor-corrected capacitive position sensor (CPS) signal, ã_sc, can be calculated based on the signals from the seismometer and rotation sensor on the ground

$$\begin{equation}{\tilde {\theta }}_{\mathrm{r}}={\tilde {\theta }}_{\mathrm{g}}-\mathrm{T}\mathrm{C}\left({\tilde {\theta }}_{\mathrm{g}}+{\tilde{n}}_{\text{BRS}}\right)\end{equation} \tag{ 1 }$$

$$\begin{equation}{\tilde{a}}_{\text{sc}}={\tilde {\dot {x}}}_{\mathrm{p}}-{\tilde {\dot {x}}}_{\mathrm{g}}+{\tilde{n}}_{\text{CPS}}+\mathrm{S}\mathrm{C}\left({\tilde {\dot {x}}}_{\mathrm{g}}-{\tilde {\theta }}_{\mathrm{r}}\frac{g}{\omega }\right),\end{equation} \tag{ 2 }$$

$$\begin{equation}{\tilde{a}}_{\text{sc}}=\underset{\text{platform}\;\text{motion}}{\underbrace{{\tilde {\dot {x}}}_{\mathrm{p}}}}+{\tilde{n}}_{\text{CPS}}-\underset{\text{ground}\;\text{injection}}{\underbrace{\left(1-\mathrm{S}\mathrm{C}\right){\tilde {\dot {x}}}_{\mathrm{g}}}}-\underset{\text{tilt}\;\text{injection}}{\underbrace{\mathrm{S}\mathrm{C}\left({\tilde {\theta }}_{\mathrm{r}}\frac{g}{\omega }+{\tilde{n}}_{\text{STS}2}\right)}}.\end{equation} \tag{ 3 }$$

Here $\tilde {\dot {x}}$ denotes translational motion, $\tilde {\theta }$ tilt motion, while subscripts g and p show whether this is associated with ground or platform motion respectively. The complex frequency responses of the sensor-correction and tilt-correction filters are denoted by SC and TC respectively, ω is the angular frequency, and g is acceleration due to gravity. The 'residual tilt', ${\tilde {\theta }}_{\mathrm{r}}$, is equal to the ground tilt in degrees of freedom where there is no BRS, as described below. The ñ terms denote the self-noise for the instrument indicated in the subscript. Table 1 describes the role of each instrument and shows the contributions that make up the output signals.

Table 1. The instruments used in the construction of the SC cost. The names of the sensors are: BRS—beam rotation sensors, STS2—Streckeisen STS-2 force-feedback seismometer, CPS—capacitive position sensor, T240—trillium T240 force-feedback seismometer, and GS13—Geotech Instruments GS-13 short-period seismometer.

Sensor Name	Measures	Contributions
BRS	Ground inertial rotation	${\tilde {\theta }}_{\mathrm{g}}+{\tilde{n}}_{\text{BRS}}$
STS2	Ground inertial translation	${\tilde {\dot {x}}}_{\mathrm{g}}+{\tilde{n}}_{\text{STS}2}+{\tilde {\theta }}_{\mathrm{g}}\frac{g}{\omega }$
CPS	Relative platform velocity	${\tilde {\dot {x}}}_{\mathrm{p}}-{\tilde {\dot {x}}}_{\mathrm{g}}+{\tilde{n}}_{\text{CPS}}$
T240	BSC-platform inertial translation	${\tilde {\dot {x}}}_{\mathrm{p}}+{\tilde{n}}_{\text{T}240}+{\tilde {\theta }}_{\mathrm{p}}\frac{g}{\omega }$
GS13	HAM-platform inertial translation	${\tilde {\dot {x}}}_{\mathrm{p}}+{\tilde{n}}_{\text{GS}13}+{\tilde {\theta }}_{\mathrm{p}}\frac{g}{\omega }$

The ground and platform inertial motions are measured by seismometers. The ground rotation can only be measured where there is a beam rotation sensor (BRS) [29]. At the time of writing, there were four BRS instruments at LIGO Livingston Observatory (LLO), measuring the rotation that affects degrees of freedom parallel to the beam-line direction at the corner-station and at both end stations. At LIGO Hanford Observatory (LHO) there were two BRS instruments, one in the beam-direction at each end station.

In the ideal case, the sensor-corrected CPS signal only includes the platform motion, which can subsequently be suppressed by feedback control. By constructing the sensor-corrected CPS signal offline from 2 h stretches of data, many different sensor-correction filters can be implemented and tested without requiring valuable operating time of the LIGO observatories.

The optimisation problem of the sensor-correction filter is now clear: we must add as much of the inertial ground-motion term as possible to the CPS signal, while minimising the injection of tilt and inertial-sensor self-noise. The compromise between the two terms in equation (3) is complicated further by the small frequency separation of important dynamics in their spectra. The dominant contribution to the RMS ground translation occurs between 0.1 and 0.3 Hz (the secondary micro-seismic peak), and the tilt-injection begins to dominate the seismometer signals below approximately 0.05 Hz [30–32].

The SC filter can be thought of as one part of a complementary pair of 'blending' filters [33], an example of SC filters can be found in figure 6. The sensor-correction filter attenuates tilt and IS noise, while the ground injection is reduced by a factor of one minus the sensor-correction filter as can be seen from equation (3). When installed into the detector only the SC filter is installed, the filter complement is used to show the ground injection suppression.

The optimisation is restricted to a band between 4 mHz and 4 Hz. Below 4 mHz, the effect of all ISs is negligible on the platform RMS velocity [32]. Above 4 Hz, the CPS contribution to the platform feedback signal becomes negligible [15, 16], due to the suppression provided by the high and low pass blending filters respectively.

The design of a cost function often encodes much of the complexity of an optimisation problem. To aid this process, we created a cost that was the integral over frequency of the tilt-injection and ground-injection terms identified in equation (3). The ability to directly compare these spectral integrands with spectra of the input signals was of great practical benefit during debugging, and helped to shape both the cost function and performance penalties described below. The final cost is also approximately equal to the residual RMS velocity of the platform.

The RMS velocity of the platform is a useful figure of merit (FOM) because it is correlated with the scattered light performance of the interferometer [34]. It also balances the need to control drift at low frequencies (better evaluated by the RMS position) and limiting control forces at high-frequencies (better evaluated by RMS acceleration).

At frequencies below approximately 50 mHz, the ground seismometer signal becomes completely dominated by tilt caused by wind pressure acting on the buildings [32]. Furthermore, the STS2 is AC-coupled at 8 mHz, and its output is no longer an inertial measure of the ground velocity. Therefore the STS2 signal below 50 mHz is considered to be part of the tilt-injection spectrum. Above this frequency, a simple power-law extrapolation of the tilt spectrum is made, as shown in figure 2. The small region from 50–100 mHz, which contains the primary micro-seismic peak, contains little spectral power and is de-weighted by the ground weighting function shown in figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The contribution to the total cost function is therefore the integral of the tilt-injection spectrum multiplied by the frequency response of the sensor-correction filter.

With no sensor-correction filter applied, the inertial ground motion can be completely transferred to the platform through the CPS sensor. It is therefore important to reduce this term as much as possible, especially at the secondary micro-seismic peak, seen between 0.1 and 0.3 Hz in the plots below.

Each isolation platform has high Q, multi-stage suspensions that further reduce vibration above approximately 1 Hz. The lowest frequency pendulum-modes of the 'quadruple' suspensions, at 0.45 Hz and 1 Hz, are responsible for a significant fraction of the residual motion in LIGO, and it was important to include these dynamics in our cost. The weighting shown in figure 3 includes a simplified and broadened fit of the suspension response [35].

At frequencies below 0.1 Hz the ground injection cost is linearly de-weighted to zero. There are two reasons for this: first, the spectrum of the ground seismometer is increasingly dominated by tilt, and second, actual inertial translation is almost completely common-mode, even down the 4 km long arms of LIGO. A vanishingly small actuation force is needed to combat the differential component. At frequencies above 1 Hz, the CPS sensor plays very little role in the motion of the platform, and the cost is de-weighted to improve the numerical precision of the optimisation in the critical 0.1–1 Hz region.

There are therefore two components to the final cost of the swarming process, the ground and tilt injection. The ground injection signal is the previously calculated sensor corrected CPS signal minus the inertial platform motion, ${\tilde {\dot {x}}}_{\text{inj}}$. This term is calculated in the time domain to account for the finite coherences between the ground seismometer, CPS, and IS on the platform. The tilt injection term is calculated in the frequency domain due to the estimation processed used to calculate the tilt integrand, $\tilde {b}$ above 50 mHz. This is multiplied by the frequency response of the SC filter producing the tilt injection term.

The ground injection term, ${\tilde {\dot {x}}}_{\text{inj}}$, is given by equation (4) and is defined here in the frequency domain. ã_sc is the previously calculated sensor corrected CPS, ${\tilde {\dot {x}}}_{\mathrm{p}}$ is the platform motion, ñ_IS is the noise of the IS, this is either a GS13 or a T240 and ${\tilde {\theta }}_{\mathrm{p}}$ is the platform tilt. In reality, the swarming process calculates this in the time domain to ensure the swarmed filter is stable. The tilt injection term, ${\tilde {\theta }}_{\text{inj}}$, is given by,

$$\begin{equation}{\tilde {\dot {x}}}_{\text{inj}}={\tilde {a}}_{\text{sc}}-{\dot {\tilde {x}}}_{\mathrm{p}}+{\tilde{n}}_{\text{IS}}+{\tilde {\theta }}_{\mathrm{p}}\frac{g}{\omega },\end{equation} \tag{ 4 }$$

$$\begin{equation}{\tilde {\theta }}_{\text{inj}}=\mathrm{S}\mathrm{C}\enspace \tilde {b},\end{equation} \tag{ 5 }$$

and is calculated in the frequency domain. Here the tilt fit curve, $\tilde {b}$, shown in figure 2, is multiplied by the SC filter and the ground injection term is multiplied by the ground weighting term, W_g, shown in figure 3. The tilt and ground injection terms are both integrated over the entire frequency band and summed in quadrature producing the final cost, given by

$$\begin{equation}{\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{t}}^{2}={\int }_{0.004\enspace \mathrm{H}\mathrm{z}}^{4\enspace \mathrm{H}\mathrm{z}}\vert {W}_{\mathrm{g}}{\tilde {x}}_{\text{inj}}{\vert }^{2}\enspace \mathrm{d}f+{\int }_{0.004\enspace \mathrm{H}\mathrm{z}}^{4\enspace \mathrm{H}\mathrm{z}}\vert {\tilde {\theta }}_{\text{inj}}{\vert }^{2}\enspace \mathrm{d}f.\end{equation} \tag{ 6 }$$

Using on-site data provides many advantages towards filter design. With the SC filter, the output of the platform seismometer acts as an in situ test of the quality of the filter. The goal of SC is to remove the ground motion in the CPS, such that the CPS follows the platform motion creating a quasi-IS. The T240 measures this and can therefore be used as a measure of quality, when factoring in the noise and tilt as measured by this sensor.

The designed filter can then be evaluated using other sets of input motion to verify their performance in a range of environmental conditions. Ensuring the cost function outlined above is physically motivated allows for quick comparisons between the performance of multiple filters and filters across different sets of environmental conditions.

Particle swarm routines aim to simulate social behaviour to explore a parameter space. A number of 'particles' are generated with randomised parameters as part of the initial population. In the context of this work, a particle contain the information necessary to build a zero pole gain (ZPK) system that defines the sensor-correction filter. A ZPK function can be defined simply as a mathematical function of the following form

$$\begin{equation}\frac{{X}_{\text{out}}}{{X}_{\text{in}}}=k\frac{\prod _{i=1}^{I}\left(s-{z}_{\mathrm{i}}\right)}{\prod _{j=1}^{J}\left(s-{p}_{\mathrm{j}}\right)},\end{equation} \tag{ 7 }$$

where X_out is the output, X_in is the input, z_i are the zeros, p_j are the poles and k is the system gain. The mapping from the degrees of freedom of the particles and the ZPK system is described in table 2. The particles are all then put through the cost function by the particle swarm routine. The cost function builds the ZPK of the filter from the particle's parameters and then a cost for this particle is calculated using the relevant equations discussed in section 3.4. Particles then determine how to change their parameters for the next generation based on the overall global minimum cost, the local minimum, and a random 'kick'. The process is then repeated until an exit condition of the swarm is met. During initial testing a swarm size of 500–1500 particles was able to sufficiently explore the parameter space whilst still enabling many iterations to be performed in a runtime less than 12 h on a typical laptop. Run time is linearly related to the number of particles as the run time is dominated by the number of amplitude spectral densities it must calculate for each particle.

To prevent overfitting to the training data the swarmed SC filter can only contain four complex poles and three complex zeros. The filters are seeded with a single real pole and three real zeros at 0 Hz to shape the filter and guarantee the required roll off at low frequency.

The use of a particle swarm allows for an unguided search of the created parameter space. It is forced to construct the initial generation of filters with random parameters for the frequency and Q (for complex poles and zeros) so that the parameter space for any potential filters is fully explored.

How the parameter space is explored by the optimisation routine was the first area considered when refining the filter generation tool. By choosing the space in which the roots are generated, the weighting of filter construction in early generations can be skewed to allow for faster and more reliable convergence. The positions of roots in an 'ideal' filter will typically not be distributed uniformly. If the filters were generated in a linear frequency space, where the value generated is the position of the root, then shaping the low frequency response in a wide data range would take a much greater optimisation time. Instead, to account for various physical factors in certain frequency bands, roots will typically be spread over a logarithmic scale with a high density at crossover points to allow for better shaping.

The degeneracy of the space is the largest problem with a natural frequency creation, (where the value of the particle is the frequency of the root). Any two roots of a polynomial can be interchanged and the resulting transfer function will remain the same, displayed by figure 4 (top). This artificially inflates the parameter space by a factor proportional to the number of roots factorial. Larger parameter spaces will take longer for the optimisation routine to fully explore. To solve this problem, roots were ordered as shown in figure 4 (bottom). Here the lowest frequency root is defined in natural frequency space and the subsequent roots defined in terms of the difference to the last one generated. This orders the roots, solving the degeneracy problem. Ordering of the roots presents different problems to the filter construction in the form of imposing boundaries on root placement. To maintain the appropriate slopes outside of the studied region, no poles or zeros can be placed there. When building a filter, the range of acceptable poles and zeros should therefore be bound to the range of frequencies that the data covers. However, by defining the relations between roots, no boundary conditions can be imposed as shown by figure 4, where each jump is less than the range of frequencies allowed, but still the final root lies outside the acceptable range. This wastes a large amount of the available computing power on filters which must be discarded.

Zoom In Zoom Out Reset image size

With effort on two fronts this problem was overcome. The first thing was to add a large scaling cost when an invalid filter was produced. This extra cost was based on a multiplicative term of how many roots and how far each root was outside the acceptable range. This created a 'bucket' in the parameter space which allowed particles to drift back into the acceptable range by imposing a gradient towards acceptability. The second tactic was a careful initialisation of the swarm matrix. Before the swarm was run, a first generation of filters was created in natural frequency space and then ordered from highest to lowest. The difference between these was then calculated and used as the parameters for generation in the swarm. This meant initial generation of particles were all viable filters and any particle which drifted out of range was pushed back by the aforementioned bucket. A graph showing the convergence of the different methods discussed is shown in figure 5, showing that after 120 generations with 500 particles the ordered generation initialised jumps achieved a significant improvement in best cost. It showed that both steps were necessary to see an improvement.

Zoom In Zoom Out Reset image size

The parameters used by the swarm to make a control filter are shown in table 2.

The data used for filter generation was collected during the commissioning break between the second and third observing runs, when the detector was in a damped only state (the isolation forward control loops were not in in operation). This simplifies the sample data as the many control techniques did not have to be considered in the processing of the input data, nor was it necessary to account for closed loop stability requirements.

The swarmed filters were created using data collected during times of above average environmental disturbance (500 nm s⁻¹ bandwidth-limited RMS-velocity between 0.1 Hz and 0.3 Hz motion, and 8–9 m s⁻¹ windspeed). Filters trained and built during different environmental conditions are discussed in chapter 6 of Cooper's thesis [28]. Since the detectors have the ability to switch between different filters, it would be possible to generate filters to account for seasonal variations of ground motion and wind speed, though this is outside the scope of the paper.

In the corner stations of both LIGO observatories one sensor-correction filter is applied to the STS-2 signal in each degree of freedom, and distributed to all chambers. This ensures that the residual tilt injection that couples through the SC path is the same for each of the isolation platforms, reducing differential platform translational motion [28] (chapter 6). Therefore when optimising the sensor-correction filters for the corner station, a filter must be produced that is effective in all chambers. To produce the required input data for the swarming process, signals from each chamber on the 'beam axis' were averaged together. For each of the horizontal translational degrees of freedom (X and Y) we combine the signals from all the relevant chambers in the main building.

The filters produced by the swarm were tested at the Hanford site between the 11th of October 2019 and 29th October 2019. The observatory had several other events occurring that resulted in the Internal Seismic Isolation (ISI) platforms operating at atmospheric pressure leading to overall higher noise. On site operators noted that the filter injected more ground tilt into the ISI, this is due to the filter having a factor 6 higher gain at low frequency compared to the previous SC filter.

The direct comparison of two SC filters that were active during two different stretches of time is not trivial due to the differences in the input ground and wind conditions. As a result, the cost function that was used to design the filters would not be appropriate. Moreover, the GS-13, used on the Horizontal Access Module (HAM) ISIs as a witness sensor is limited by its own self noise below approximately 0.05 Hz, is susceptible to platform tilt, and thus would not be a good witness sensor below this frequencies. To make a fair comparison, filter performance was evaluated during periods with similar environmental conditions. Both times had an average level of microseismic motion and low wind speed. To compensate for the remaining differences in input ground motion, the FOM term is weighted by the ratio of the amplitude spectral densities of the ground seismometer signals from the two testing periods.

As such, below approximately 0.1 Hz the performance of the filters was determined by the sensor corrected CPS signal and by the GS-13 or T-240 above this frequency. The exact form of the FOM is given by equation (9). The sensor corrected CPS provides information on the relative tilt injection and partially the suppression of the microseismic peak while the on-platform IS measures the suppression of the microseismic peak and the isolation performance above 0.1 Hz. It is important to quantify the FOM in terms of RMS velocity, as this is strongly correlated with scattered light in the interferometer moreover the velocity is a compromise between the need to control drift at low frequency (best evaluated by the RMS position) and limiting control forces at high-frequencies (best evaluated by the RMS acceleration).

To quantify the performance of a SC filter, the amplitude spectral density of the signals from each sensor described in table 1 was recorded. The spectra were standardised by converting them into velocity and normalising them to account for the frequency-dependence of their sensitivity, as appropriate. The spectra are then combined such that the CPS is used between 1 mHz and 0.1 Hz and the on-platform IS, is used between 0.1 Hz and 1.5 Hz. The IS is either a GS13 for the HAM chambers or a T240 for the BSC chambers. The resulting spectrum is defined as

$$\begin{equation}{\tilde {\dot {x}}}_{m}=\begin{cases}\mathrm{C}\mathrm{P}\mathrm{S},\quad \hfill & 0.001{\leqslant}f{\leqslant}0.1\hfill \\ \mathrm{I}\mathrm{S},\quad \hfill & 0.1{< }f{\leqslant}1.5.\hfill \end{cases}\end{equation} \tag{ 8 }$$

This 'combined' signal is then integrated over the band as shown producing the FOM,

$$\begin{equation}\mathrm{F}\mathrm{O}\mathrm{M}={\int }_{0.001}^{1.5}\vert {\tilde {\dot {x}}}_{m}{\vert }^{2}\enspace \mathrm{d}f.\end{equation} \tag{ 9 }$$

A comparison of the swarmed SC filter and the current filter is shown in figure 6. The swarmed filter has significantly less gain peaking around 50 mHz resulting in more tilt injection below 10 mHz. Figure 7 shows a comparison between the two SC filters when evaluated on HAM2 in the X degree of freedom. The RMS of each of the traces are calculated from 1.5 Hz to prevent the large spikes at 1.7 and 2.2 Hz due to the suspension resonances from saturating the RMS. Overall, the swarmed SC filter results in 27% lower RMS velocity of the ISI compared to the current filter when the input ground motion has been normalised. While the swarmed filter causes a factor 10 more tilt injection at 1 mHz this does little to affect the overall RMS of the platform. The swarmed filter makes most of its improvement between 0.4 and 0.9 Hz, of particular interest is the performance difference at 0.45 Hz, one of the suspension resonances, where the swarmed filter results in 33% less motion than the current filter.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 8 shows a comparison between the swarmed and current SC filter on ITMX in the X degree of freedom. The swarmed filter results in a more modest 18% reduction in RMS velocity of the ISI. Again, the swarmed filter causes more tilt injection between 30 and 200 mHz, though this has little effect on the overall velocity RMS. The biggest difference between the two filters comes between 0.4 and 0.9 Hz. At 0.45 Hz the swarmed filter results in 68% less motion than the current SC filter.

Zoom In Zoom Out Reset image size