Particle swarming of sensor correction filters

Reducing the impact of seismic activity on the motion of suspended optics is essential for the operation of ground-based gravitational wave detectors. During periods of increased seismic activity, low-frequency ground translation and tilt cause the Advanced LIGO observatories to lose `lock', reducing their duty cycles. This paper applies modern global-optimisation algorithms to aid in the design of the `sensor correction' filter, used in the control of the active platforms. It is shown that a particle swarm algorithm that minimises a cost-function approximating the differential RMS velocity between platforms can produce control filters that perform better across most frequencies in the control bandwidth than those currently installed. These tests were conducted using training data from the LIGO Hanford Observatory seismic instruments and simulations of the HAM-ISI (Horizontal Access Module Internal Seismic Isolation) platforms. These results show that new methods of producing control filters are ready for use at LIGO. The filters were implemented at LIGO's Hanford Observatory, and use the resulting data to refine the cost function.


I. INTRODUCTION
The first observation of gravitational waves from a binary black hole [1] was just a few years ago, but there has since been 11 confirmed astrophysical signals during the O1 and O2 observing periods [2], and many more event candidates during O3 [3]. For the Advanced LIGO detectors [4,5], the increased frequency of detections has placed a greater emphasis on the stability and consistency of the observatories.
The extraordinary sensitivity of the detectors places stringent requirements on the residual vibration in the measurement band [6]. Additionally, for the observatories to enter into an operational state, the mirrors must be 'locked' to an ideal operating point [7][8][9]. Low frequency ground motion can cause the mirrors to move out of this operating point, a process that is known as loosing lock, if this motion is not sufficiently suppressed. To reduce low-frequency vibrations, a complex multi-stage active isolation system is employed [10,11]. The active control employs an array of differential and inertial sensors on each stage. At frequencies between 0.3 and 15 Hz, feedback control provides most of the vibration suppression; near the microseismic peak (50-300 mHz) [12] it is done by feed-forward control.
'Sensor correction' (SC) is a control technique where signals from ground seismometers are filtered and summed with the output of differential displacement sensors that measure the position of LIGO's isolated platforms, as shown in Figure 1. This removes the ground contribution of the displacement sensor to yield the inertial platform motion. The aim is to use the excellent noise performance of vault-grade seismometers at low frequencies. However, there is a cut-off below 0.1 Hz where ground seismometers are dominated by their tilt susceptibility [13]. This noise is avoided by using the displacement sensors as the feedback signal. A good sensor correction filter aggressively tackles the translational ground motion without including frequencies where noise dominates, in particular, tilting contamination. This paper presents a novel approach for improving sensor-correction filter design [14,15] based on particle swarm optimisation. Improved filters are found with an unguided search using real data, for a variety of environmental conditions. By using real data-streams from all three types of installed sensors, and applying IIR filters exactly as LIGO's real-time control system does, the tool becomes immune to a slew of systematic issues including finite coherence, differing transfer-functions, and numerical implementation effects. A cost function was designed to mimic the physical quantity we believe is most important for LIGO's performance, the RMS inertial velocity. The velocity spectrum was shaped in an a-priori fashion to account for known physical effects of common-mode rejection and resonances of the suspended optics. Finally, the search variables were reparameterised to remove degeneracy in the search space.
In section 2, comparisons between current filter designs methods and this technique are discussed. Section 3 details the implementation of the particle swarm. Section 4 reviews the filters and their performance when installed at the Hanford site.

II. NUMERICALLY OPTIMISED FILTER DESIGN
Filter design is a challenge in many fields; each of which have their own relevant filter design methods. The SC filters used at LIGO are hand-tuned by the operators and commissioners.
While this can be an effective design method, it requires substantial experience and an understanding of the performance requirements. Until now there has not been a quantitative performance metric, and it has been difficult to evaluate the relative quality of different filters.
Furthermore, there exists the possibility that the designs are limited by human biases. Additionally, since this design and testing process is time consuming, the filters are only replaced if there are obvious performance issues.
Other fields have developed a wide array of tools for generating of control filters to suit their specific needs. An example often used for motion suppression in car suspensions is reinforced learning automata, such as described by Howell et al. [16] and other later groups [17]. Practical limitations prevent this method from direct application to gravitational wave detectors. It would require in operation testing so that the direct output of the filter can be evaluated. In operation measurements are not feasible for a gravitational-wave detector, due to the unacceptable loss of observation time of the testing process; this strongly motivates the development of rigorous off-line testing procedures.
Storn et al. [18] pioneered a method for when there is a priori knowledge of the desired filter response. They use a differential evolution optimiser to fit poles and zeros to a prescribed phase and magnitude response, whilst employing cost 'penalties' in order to meet stringent stability conditions. This method is not well suited for designing a sensor-correction filter as the desired response is unknown due to the complex noise dynamics and limited coherence of the sensors in the isolation system. Other groups have used particle swarms to design control filters based upon Storn's work [19,20]. They show that the particle swarm optimisation algorithm can effectively design filters due to its ability to sample large parameter spaces. A direct comparison between genetic algorithms and particle swarms is presented in [19]. The particle swarm reached a lower minimum more consistently, though there was little difference between the two methods; a similar conclusion was reached during the development of the optimisation method presented here. Modifications to the particle-swarm optimiser are described in [20] that made it similar to the "branch and prune" optimisation technique. The changes significantly improved the optimisation performance, and this illustrates the relative flexibility and adaptability of the particle-swarm technique.
A numerical optimiser enables a conceptually simple way to design SC filters. The only bound on the solver are the asymptotes outside the band of interest (4 mHz-4 Hz), which are entirely controlled by the number of poles and zeros in the filter. These asymptotes have strong physical motivation, however. Within the entire frequency band of interest a completely unguided and randomly seeded search is made. This allows for non-obvious, but beneficial features to be incorporated in the filter. Although conceptually simple, implementing such an optimiser proved to be a non trivial task.
Data collected from the wide array of on-site sensors was used for the calculation of filter performance during generation. Calculating the overall ground injection using this numerical, time domain approach allows for the inclusion of finite coherence between multiple sensors preventing the over-estimation of the sensor correction filter. Furthermore, many real world effects of implementing filters such as differing transfer functions between sensors are already incorporated in the resultant filters. Because the data processed is in the time domain, only real coherence can be subtracted. Finally, doing so with time domain data is the numerically closest operation to on site action. This means an unstable filter will be inherently penalised as it will cause a large increase motion in the system. While this work focuses on sensor correction filters, the tools modular design means that it can be adapted to design various control filters within the seismic isolation system. The problem for filter generation then becomes purely a matter of finding an appropriate cost function.

III. GUIDING THE SWARM
All numerical optimisation processes require a cost function. Using a particle swarm optimisation routine allows for a computationally intensive one. The cost function mimics the implementation of the sensor-correction filter performed by the LIGO's real-time control system. From there, the filter's performance was analysed to estimate the disturbances injected into LIGO's isolated platforms. Additional weighting was applied to account for key features such as suspension resonances.

Sensor Name
Measures Contributions BRS Ground inertial rotationθ g +ñ BRS

CPS
Relative platform velocityx p −x g +ñ CPS

T240
BSC-platform inertial translationx p +ñ T240 +θ p g ω GS13 HAM-platform inertial translationx p +ñ GS13 +θ p g ω Table I. The instruments used in the construction of the sensor correction 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. (1) Herex denotes translational motion,θ 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',θ 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 I  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 long sections of data, many different sensor-correction filters can be implemented and tested without requiring valuable operating time of the LIGO observatories.
Signals are processed to flatten the frequency response in band and are converted to velocity and radians for translation and rotational motion respectively.
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 [23][24][25].
The sensor correction filter can be thought of as one part of a complementary pair of 'blending' filters [26], an example of sensor correction filters can be found in figure 6. The sensorcorrection filter attenuates tilt and inertial sensor noise, while the ground injection is reduced by a factor of one minus the sensor-correction filter. When installed into the detector only the sensor correction 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 inertial sensors is negligible on the platform RMS velocity [25]. Above 4 Hz, the CPS contribution to the platform feedback signal becomes negligible [10,11], due to the suppression provided by the high and low pass blending filters respectively.

B. A physically motivated cost function
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 tiltinjection 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 because it is correlated with the scattered light performance of the interferometer. 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).

C. Tilt injection
At frequencies below approximately 50 mHz, the ground seismometer signal becomes completely dominated by tilt caused by wind pressure acting on the buildings [25]. 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 microseismic peak, contains little spectral power and is de-weighted by the ground weighting function shown in figure 3.
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. A linear fit is used in this region to estimate the tilt coupling for the rest of the frequency band.

D. Ground injection
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.
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,x inj . This term is calculated in the time domain to account for the finite coherences between the ground seismometer, CPS, and inertial sensor on the platform. The tilt injection term is calculated in the frequency domain due to the estimation processed used to calculate the tilt integrand,b above 50 mHz. This is multiplied by the frequency response of the sensor correction filter producing the tilt injection term.
The ground injection term,x inj , is given by Equation 4 and is defined here in the frequency domain.ã sc is the previously calculated sensor corrected CPS,x p is the platform motion,ñ IS is the noise of the inertial sensor, this is either a GS13 or a T240 andθ 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,θ inj , is given by, and is calculated in the frequency domain. Here the tilt fit curve,b, shown in Figure 2, is multiplied by the sensor correction filter and the ground injection term is multiplied by the ground weighting term, W g , shown in Using on-site data provides many advantages towards filter design. With the sensor correction filter, the output of the platform seismometer acts as an in-situ test of the quality of the filter.
The goal of sensor correction is to remove the ground motion in the CPS, such that the CPS follows the platform motion creating a quasi-inertial sensor. 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.

E. Particle Swarm Optimisation
Particle swarm routines aim to simulate social behavior to explore a parameter space. A number of 'particles' are generated with randomised parameters. In the context of this work, these contain the information necessary to build a control filter Zero Pole Gain (ZPK) function.
The routine then calculates a cost for each of the particles using the cost function. In this case it must also build the ZPK of the filter from the generated parameters and then calculate the cost using the relevant equations discussed in section III A. Particles then determine how to change their parameters based on the overall global minimum, 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 reasonable time frame on a typical laptop.
For use in LIGO, each designed sensor correction filter should use no more than 20 second order sections, this means the swarmed sensor correction filter can only contain 4 complex poles and 3 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.

F. Reparameterisation
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 a "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.
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 Two nodes are interchangeable  Table II. A list of the parameters sent to the swarm and how they are converted into a filter function.
Each frequency and Q is used to create one complementary pair of complex roots. The range of Q values allow the poles and zeroes to become slightly over-damped, making a pair of real poles.
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.
The parameters used by the swarm to make a control filter are shown in Table II. the effect of swarm initialisation. All data is measured relative to the best filter. This corresponds to a 13% improvement in frequency space. The tests were run several times and those with similar cost in initial generation displayed.

IV. IMPLEMENTATION AT LIGO HANFORD OBSERVATORY
A. Choice of training data 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, it is assumed that if the filter was designed during these times it would be more robust in a wider range of environmental conditions. 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.

B. Optimising Multiple Chambers
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 sensor correction path is the same for each of the isolation platforms, reducing differential platform translational motion [21](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.

C. Comparison with current filter
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 ISI (Internal Seismic Isolation) 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 sensor correction filter.
The direct comparison of two sensor correction filters that were active during two different stretches of time is complex 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 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.
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 figure of merit is given by Equation 8. The sensor corrected CPS provides information on the relative tilt injection and partially the suppression of the microseismic peak while the on-platform inertial sensor measures the suppression of the microseismic peak and the isolation performance above 0.1 Hz.
To quantify the performance of each sensor correction filter, the amplitude spectral density was taken of each sensors in table I, whereby the signals are converted into displacement and are plant inverted where appropriate. The signal from the CPS from 1 mHz is stitched with the on-platform inertial sensor, IS, either the GS13 or T240 as appropriate, at 0.1 Hz as shown by, . This 'stitched' signal is then integrated over the band as shown producing the figure of merit (FOM), .
A comparison of the swarmed sensor correction 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 sensor correction filters when evaluated on HAM2 in the X degree of freedom. To compensate for the differences in input ground motion when the filters were installed, the figure of merit term, is weighted by the ratio of the signal from the ground seismometers during the two times. 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 sensor correction 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.

V. CONCLUSIONS
Many strategies are being developed to decrease lock loss due to ground motions at the LIGO facilities. Improving the control filter design presents an opportunity to make significant gains with simple and affordable changes. A robust working tool to design sensor correction filters has been developed, and filters created by this method have been deployed on-site. This work focused on one specific filter although its design principles can be adapted for other filters. The blend filters present in all the isolation systems are excellent candidates for testing. Designing control filters in this manner presented some challenges that had to be overcome, but techniques were developed, trialled, and tested throughout have solved many of the issues and may prove relevant to other applications. With this tool, better seismic isolation can be achieved leading to more observing time, more source detections and, ultimately, a greater understanding of some of the most fascinating phenomena in the universe.