# H-infinity Optimization for Active Seismic Isolation Systems in
Gravitational-Wave Detectors

Terrence Tsang<EMAIL_ADDRESS>Cardiff University, Cardiff CF24 3AA,
United Kindom  Fabi$\acute{\mathrm{a}}$n Erasmo Pe$\tilde{\mathrm{n}}$a
Arellano  Universidad de Guadalajara, 44100 Guadalajara, Jalisco, Mexico
Takafumi Ushiba  Institute for Cosmic Ray Research, KAGRA Observatory, The
University of Tokyo, 238 Higashi-Mozumi, Kamioka-cho, Hida City, Gifu
506-1205, Japan  Ryutaro Takahashi  Gravitational Wave Science Project,
National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka City, Tokyo
181-8588, Japan  Yoichi Aso  Gravitational Wave Science Project, National
Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka City, Tokyo 181-8588,
Japan  Katherine Dooley  Cardiff University, Cardiff CF24 3AA, United Kindom

###### Abstract

Low-frequency seismic noise poses a significant challenge in gravitational-
wave detectors, prompting the use of active seismic isolation systems. These
systems employ advanced control configurations, including feedback control,
sensor correction, and sensor fusion, to mitigate seismic disturbances.
Optimal performance of these systems can be achieved through the application
of the H-infinity method, which optimizes control filters. As a demonstration,
the H-infinity method was employed to optimize a sensor correction filter and
complementary filters for sensor fusion for a KAGRA Type-B suspension,
resulting in the development of an instrumentation-optimal “super sensor”. The
super sensor exhibits overall displacement noise RMS performance comparable to
a relative sensor, yet offers significantly better band-limited RMS
performance around the secondary microseism frequency, making it highly robust
to seismic noise. This highlights the efficacy of the H-infinity method in
generating seismic isolation filters tailored to instrument limitations,
thereby transforming the system into a fully functional active seismic
isolation system.

††preprint: APS/123-QED

##  I Introduction

###  I.1 Ground-Based gravitational-wave detector

Gravitational-wave observatories, such as Advanced LIGO [18], Advanced Virgo
[11], and KAGRA [14], are expansive laser interferometers spanning kilometers
in length. Their primary purpose is to detect gravitational radiation
emanating from astronomical events involving massive objects undergoing
significant accelerations.

To date, the gravitational waves captured by these observatories have
originated from compact binary inspirals. These include binary neutron stars,
binary black holes, and neutron star-black hole binaries [5, 7, 9, 8].
Additionally, ongoing searches target different types of gravitational-wave
signals. These include continuous gravitational waves emitted by asymmetric
rotating neutron stars, burst gravitational waves stemming from core collapse
supernovae, and a stochastic gravitational-wave background [4, 22].

###  I.2 Noise sources

Gravitational waves are minuscule distortions in the fabric of spacetime,
necessitating the use of incredibly sensitive laser interferometers for their
detection. Modern gravitational-wave detectors achieve remarkable sensitivity,
capable of detecting strain at levels lower than $10^{-21}\
\mathrm{m}/\mathrm{m}/\sqrt{\mathrm{Hz}}$ at 100 Hz [3]. However, achieving
such high sensitivity requires effective suppression of various sources of
noise.

Three primary sources of noise in gravitational-wave detectors are quantum
noise, thermal noise, scattered light noise, control noise, and seismic noise
[27]. Seismic noise in gravitational-wave detectors arises from the random
motion of the ground, leading to decoherent motion between the mirrors and
causing random fluctuations in arm lengths. Seismic noise impacts
gravitational-wave detectors in two main ways. In the detection band,
typically around 10 to few kHz, seismic noise induces strain noise orders of
magnitude higher than gravitational-wave strain itself. At lower frequencies,
seismic noise can disrupt the control systems that sustain the optical
cavities via saturation of digital control signals or misalignment of optics.
To counteract these effects, the main optics of the interferometer are
suspended using multiple pendulum suspensions, while active seismic isolation
systems help mitigate the disruptive effects of seismic noise [23, 12, 15].

###  I.3 Seismic isolation

Figure 1: Amplitude spectral densities of the example seismic and sensor
noises taken from the KAGRA detector.

As mentioned previously, the ground motion, also known as seismic noise,
impacts gravitational-wave detectors in two main ways:

  * •

Small differential motions between the test mass mirrors above 10 Hz lead to
arm-length fluctuations, which are quantified as one of the strain noises of
the gravitational-wave detector.

  * •

Larger motions at lower frequencies can cause the optical cavities formed by
the mirrors to deviate from their operating point, potentially disabling the
detector’s observation mode.

Addressing the former issue ensures that the gravitational-wave detector
maintains proper sensitivity, while addressing the latter improves the
detector’s duty cycle. The level of seismic noise varies depending on the
location, and the seismic noise data presented in Fig. 1 serves as an
illustrative example.

To detect gravitational waves corresponding to strains on the order of
$10^{-21}\ \mathrm{m}/\mathrm{m}/\sqrt{\mathrm{Hz}}$, the differential motion
between the test mass mirrors must be lower than $10^{-18}\
\mathrm{m}/\sqrt{\mathrm{Hz}}$, assuming an interferometer arm length on the
order of $10^{3}$ meters. The seismic noise at 10 Hz is approximately
$10^{-11}\ \mathrm{m}/\sqrt{\mathrm{Hz}}$, as depicted in Fig. 1.
Consequently, a suppression of more than $10^{7}$ times is necessary. This
suppression can be achieved by suspending the test masses with multiple
pendulums, which act as stacked spring-mass systems under the assumption of
small angles. The suppression ratio of a multiple pendulum increases
exponentially with the number of pendulums above the resonance frequencies and
is proportional to $f^{2n}$, where $f$ is frequency and $n$ is the number of
pendulum stages, For instance, a quadruple pendulum resonating at 1 Hz can
achieve a seismic noise suppression of $10^{8}$ at 10 Hz. Hence, the primary
optics in gravitational-wave detectors are typically suspended by multiple
pendulums.

The motion of the ground has most energy concentrated around 0.1 and 0.5 Hz,
as shown as a peak in Fig. 1. This is known as the secondary microseism, which
is a consistent hum caused by the interaction between ocean waves traveling in
opposite directions. In the example shown, the seismic noise peaks at around
$10^{-6}\ \mathrm{m}/\sqrt{\mathrm{Hz}}$ and has an RMS of around $2\times
10^{-7}\ \mathrm{m}$.

Gravitational-wave detectors utilize high-finesse kilometer-scale optical
cavities to improve gravitational wave detection. Maintaining the stability of
such optical cavities requires that the separation fluctuation of the test
masses remains within a fraction of the wavelength of the resonating laser
field, around $10^{-14}\ \mathrm{m}$ [6]. However, secondary microseisms can
often induce RMS motion levels orders of magnitude higher than this threshold.
For instance, at LIGO Livingston, secondary microseisms contribute to an RMS
displacement of approximately $10^{-6}$ meters, which is comparable to the
wavelength of the laser (1064 nm). Consequently, attenuation of seismic noise
at lower frequencies is essential for optimal detector performance.

Within gravitational-wave detectors, numerous layers of control systems
collaborate to minimize the residual motion of the optics to a level conducive
to sustaining optical cavities. Seismic isolation serves as the foundational
layer that mitigates external disturbances, thereby enabling higher-level
control systems, such as interferometric control systems, to operate
effectively. Its primary responsibility is to sufficiently suppress the motion
of the all optics, allowing for the engagement of higher-level control
systems.

To achieve seismic isolation at lower frequencies, optics in gravitational-
wave detectors are mounted on isolation platforms, typically together with the
entire multiple-pendulum suspension. The inertial motion of an isolation
platform is regulated via a feedback mechanism, actively suppressing the
seismic noise-induced platform motion.

The feedback control performance is limited by the ability to accurately
measure the inertial motion of the isolation platform. To ensure accurate
measurement, multiple sensors, namely a relative sensor, a feedback inertial
sensor (referred to simply as an inertial sensor hereafter), and a ground
seismometer (referred to simply as a seismometer hereafter), are commonly used
in conjunction with two control configurations: sensor correction [23] and
sensor fusion [30].

The sensor correction method aims to combine the readouts of a relative sensor
and a seismometer, producing a corrected relative sensor that reads inertial
motion by minimizing the ground coupling in the relative sensor. The sensor
fusion (or sensor blending) method combines two different sensors, in this
case, the sensor-corrected relative sensor and the inertial sensor, both
measuring the inertial motion of the isolation platform. This results in a
“super sensor” that leverages the superior aspects of both sensors, achieving
superior noise performance.

The operation of the feedback control, sensor correction, and sensor fusion
systems requires proper designs of infinite impulse response (IIR) filters,
namely the feedback controller, the sensor correction filter, and the
complementary filters, respectively. Conventionally, these filters are
designed by manually shaping the frequency responses to satisfy gain and noise
suppression requirements. However, designing such a filter manually, by
placing zeros and poles and fine-tuning its shape, can be exceedingly
challenging, labor-intensive, and inconsistent. Therefore, an optimization
approach is sought instead, aiming to robustly optimize the IIR filters.

Quantifying the required amount of seismic noise suppression poses a challenge
for active isolation systems due to the complexity of the couplings between
seismic disturbances and optical cavities. One such coupling is length-to-
angular coupling, which tilts the beam axis of the optical cavities as a
result of the mechanical nature of multiple pendulum suspensions.

Rather than setting arbitrary performance goals for active seismic isolation,
this paper concentrates on optimizing seismic control systems to achieve
instrumentational optimality for the particular control problem. The objective
of the optimization problem is to devise a methodology for optimizing control
filters that yield control performance closest to the instrumentation limit.
The instrumentation limit of a control problem is dictated by the quality of
the instrumentation, including sensors, whose noise performances can be
quantified as spectral density in the frequency domain, as illustrated by the
sensor noises depicted in Fig. 1.

The H-infinity ($\mathcal{H}_{\infty}$) control method has been proven to be
well-suited for optimizing control filters in seismic isolation systems [29,
30, 31]. In a spectral density plot, the ideal control performance should
closely approach the instrumentation limit. Consequently, the deviation of the
actual control performance from this limit serves as a valuable measure of
merit, or a cost function, which can be minimized using an optimization
algorithm. However, this deviation is a frequency-dependent quantity that
requires minimization across the entire frequency domain and cannot be easily
condensed into a single value. The H-infinity norm serves as a metric for
measuring the maximum deviation, and minimizing it is equivalent to specifying
a performance goal equivalent to the instrumentation limit. With the
capability to specify frequency-dependent control objectives, the H-infinity
method emerges as an ideal tool for optimizing seismic control systems.

This paper is organized as follows: Section II introduces active seismic
isolation systems and common seismic control techniques employed in
gravitational-wave detectors. This section delves into the three design
freedoms within active isolation systems, namely the feedback controller,
sensor correction filters, and complementary filters in a sensor fusion
configuration. Section III presents a solution to the three control problems.
It addresses the complementary filter problem using the H-infinity method and
generalizes both the sensor correction problem and the feedback control
problem into a complementary filter problem, thereby making them solvable
using a unified solution. Additionally, this section discusses practical
issues encountered during the implementation of the H-infinity method, such as
noise modeling issues, integration drifts, and numerical challenges. In
Section IV, the paper exemplifies the application of H-infinity sensor
correction and sensor fusion methods, reporting results utilizing sensors from
KAGRA suspensions. Section V explores potential future enhancements for the
H-infinity method and for active seismic isolation systems in general.
Finally, Section VI provides concluding remarks.

##  II Active seismic isolation

Figure 2: A conceptual illustration of an active isolation system suspending a
mirror in gravitational-wave detector.

Seismic isolation systems in gravitational-wave detectors rely on self-
correcting (feedback) and active cancellation (feedforward) control systems,
integrating sensors and actuators to achieve active isolation. Feedback
control hinges on accurately measuring the inertial motion of the isolation
platform and shaping the feedback controller in the frequency domain to
minimize disturbance and sensing noise simultaneously. To enhance inertial
sensing without upgrading sensors, gravitational-wave detectors commonly use
sensor correction and sensor fusion techniques [23, 17, 21, 31]. Sensor
correction removes seismic noise from relative displacement sensors using a
ground witness, while sensor fusion combines multiple sensors to achieve
superior noise performance. Both techniques require frequency-dependent filter
shaping for optimal sensing. Depicted in Fig. 2 is a conceptual representation
of a generic active isolation system, serving as a model for various systems,
such as the HAM-ISI in LIGO, the Superattenuator in Virgo, and the Type-A
suspension in KAGRA. The system features an isolation table coupled to the
ground and acts as the suspension point for primary optics. Two sensors, a
relative sensor and an inertial sensor, measure the motion of the isolation
table. The seismic noise-coupled relative sensor is corrected using sensor
correction, and the sensor-corrected relative sensor and the inertial sensor
are fused to improve noise performance. The resulting super sensor readout is
fed into a feedback loop to actively adjust the position of the isolation
table, achieving active isolation.

To facilitate further discussions on these control schemes, we outline several
assumptions. In reality, the isolation table is suspended in more than one
degree of freedom. However, for simplicity, we assume that control matrices
such as sensing and actuation matrices are properly implemented to achieve
coordinate transformations between the sensing (and actuation) and control
bases. Additionally, we assume that mechanical degrees of freedom are
adequately decoupled into orthogonal bases, enabling the treatment of full
multiple-degrees-of-freedom active isolation systems as multiple single-
degree-of-freedom systems. Consequently, the ensuing discussion focuses on a
single-degree-of-freedom active isolation, akin to the one depicted in Fig. 2.
Finally, we presume that the sensing elements, including the relative sensor,
inertial sensor, and seismometer, are calibrated and inter-calibrated into
displacement units. In practice, this is addressed through the implementation
of calibration filters that invert the frequency responses of the sensors
[31].

###  II.1 Feedback control systems

Figure 3: Feedback control system with reference tracking, disturbance and
sensing noise. $R(s)$ is the reference setpoint, $K(s)$ is the feedback
controller, $P(s)$ is the physical plant, $D(s)$ is the disturbance, $N(s)$ is
the sensing noise, and $X(s)$ is the displacement output.

Seismic isolation systems within gravitational-wave detectors heavily rely on
feedback control mechanisms to achieve isolation, especially targeting
frequencies around the secondary microseism, typically between 0.1 and 0.5 Hz.
Represented by a block diagram shown in Fig. 3, the feedback control system
involves inertial displacement $X(s)$, disturbances $D(s)$, sensing noise
$N(s)$, reference signal $R(s)$, physical system $P(s)$, and controller
$K(s)$. Here, $s$ represents the complex variable in the Laplace domain. At
non-zero frequencies, the displacement reads

$X(s)=\frac{1}{1+L(s)}D(s)-\frac{L(s)}{1+L(s)}N(s)\,,$ (1)

where $L(s)=K(s)P(s)$ is the open-loop transfer function. Note that the
reference signal $R(s)$ is typically a constant. The primary objective is to
minimize or satisfy certain frequency-dependent goals for displacement $X(s)$
through the design of controller $K(s)$.

###  II.2 Sensor correction

The mechanism enabling the removal of seismic noise from relative sensors and
effectively turn them into inertial sensors is referred to as sensor
correction. To implement sensor correction, a seismometer is placed near the
isolation platform. This seismometer, measuring ground motion $X_{g}(s)$,
serves as a witness sensor in a feed-forward noise cancellation configuration
to eliminate seismic noise coupling from the relative sensor readout. In an
ideal scenario, summing the relative readout with the inertial sensor readout
would yield a corrected readout devoid of seismic noise coupling. However,
inertial sensors typically exhibit a frequency-dependent response resembling a
high-order high-pass filter, necessitating a high-order integrator to invert
their frequency response. Consequently, readout noises, including quantization
noise, are amplified by the inverse filter, rendering inertial sensors
inherently poor at low frequency compared to relative sensors. Therefore,
direct summation of the relative and inertial readouts yields a seismic noise-
free readout containing low-frequency noise unsuitable for the precision
required for coarse interferometer alignment. Instead, the inertial signal
must be filtered before summation with the relative readout. The block diagram
of this sensor correction configuration is depicted in Fig. 4.

Figure 4: Sensor correction for relative sensors. $X(s)$ is the inertial
displacement of the isolation platform, $X_{g}(s)$ is the displacement of the
ground, $N_{\mathrm{rel}}$ is the intrinsic noise of the relative sensor,
$N_{\mathrm{seis}}(s)$ is the seismometer noise, $H_{\mathrm{sc}}$ is the
sensor correction filter, and $Y_{\mathrm{rel,sc}}(s)$ is the readout of the
sensor-corrected relative sensor.

In Fig. 4, two forward paths are depicted: the top path illustrates the
relative sensor readout, while the bottom path corresponds to the a
seismometer signal being filtered by a sensor correction filter
$H_{\mathrm{sc}}$. The two signals are then summed to form a sensor-corrected
relative readout $Y_{\mathrm{rel,sc}}(s)$, expressed as

$Y_{\mathrm{rel,sc}}(s)=X(s)+N_{\mathrm{rel,sc}}(s)\,,$ (2)

where the noise of the sensor-corrected relative sensor,
$N_{\mathrm{rel,sc}}(s)$, is given by

$N_{\mathrm{rel,sc}}(s)=N_{\mathrm{rel}}(s)+H_{\mathrm{sc}}(s)N_{\mathrm{seis}}(s)\\\
-\left[1-H_{\mathrm{sc}}(s)\right]X_{g}(s)\,,$ (3)

and $N_{\mathrm{seis}}(s)$ represents the seismometer noise. With the
corrected relative sensor serving as the feedback sensor, the noise
$N_{\mathrm{rel,sc}}(s)$ becomes another constraint alongside the disturbance
$D(s)$. However, the sensing noise in this configuration depends on the filter
$H_{\mathrm{sc}}(s)$ and can be manipulated to achieve a desirable shape in
the frequency domain.

Similar to the controller $K(s)$ in a feedback system, the sensor correction
filter $H_{\mathrm{sc}}(s)$ can be tailored to accomplish specific control
objectives, such as minimizing the sensing noise $N_{\mathrm{rel,sc}}(s)$.
Eqn. (3) delineates the corrected relative sensing noise
$N_{\mathrm{rel,sc}}(s)$ as a sum of three components. The first term
represents the intrinsic relative sensor noise $N_{\mathrm{rel}}(s)$, which
can be perceived as ambient noise unaffected by the sensor correction filter
$H_{sc}(s)$. The latter two components correspond to the filtered seismometer
noise $H_{\mathrm{sc}}(s)N_{\mathrm{seis}}(s)$ and the subtracted ground
motion $[1-H_{\mathrm{sc}}(s)]X_{g}(s)$. Both components are contingent on the
sensor correction filter $H_{sc}(s)$. However, akin to disturbance and sensing
noise, they cannot be simultaneously minimized since $H_{\mathrm{sc}}(s)$ and
its complement, $1-H_{\mathrm{sc}}(s)$, cannot be small simultaneously. This
predicament resembles the challenge encountered in the feedback control system
described in Sec. II.1, wherein designing a control filter to suppress two
noises with distinct frequency contents poses a similar difficulty.
Nonetheless, as we will elaborate in subsequent sections, the sensor
correction problem differs due to the presence of ambient noise
$N_{\mathrm{rel}}(s)$. Consequently, noises attenuated below the ambient noise
are over-suppressed, given that the sensing noise will be dominated by the
ambient.

###  II.3 Sensor fusion

Inertial sensors prove valuable in measuring ground motion for sensor
correction and can also gauge the motion of the isolation platform for
feedback control. Generally, they exhibit superior noise performance at high
frequencies compared to relative sensors. Consequently, when utilized in the
feedback configuration, they can offer better seismic isolation and introduce
less noise at high frequencies. However, as previously noted, inertial sensors
suffer from inferior low-frequency performance, which could significantly
constrain the control bandwidth. To address this limitation, a solution
involves merging the relative displacement sensor and the inertial sensor into
a so-called “super sensor”, combining the low-frequency reliability of a
(sensor-corrected) relative sensor with the high-frequency noise performance
of an inertial sensor. This integration can complement sensor correction, as
described in Sec. II.2, to achieve optimal sensing performance and
consequently enhance seismic isolation performance.

Figure 5: A sensor fusion configuration using complementary filters. $X(s)$ is
the displacement being measured, $N_{1}(s)$ and $N_{2}(s)$ are the sensing
noises of the sensors being fused, $H_{1}(s)$ and $H_{2}(s)$ are the
complementary filters, and $Y_{\mathrm{super}}(s)$ is the readout of the super
sensor.

In gravitational-wave detectors, blending two sensors is commonly achieved
through sensor fusion using complementary filters. The block diagram
illustrating a sensor fusion system with two sensors is depicted in Fig. 5.
Here, subscripts 1 and 2 denote two different sensors without loss of
generality. Both sensors measure a common displacement signal $X(s)$, each
with its respective sensor noise $N_{1}(s)$ and $N_{2}(s)$. Note that the
noises $N_{1}(s)$ are $N_{2}(s)$ refer to the deviation from the the actual
signal $X(s)$ and can be composed of a superposition of difference noise
sources. For example, a typical choice for $N_{1}(s)$ is the noise of the
sensor-corrected relative sensor, which is a combination of the relative
sensor noise, the subtracted ground motion, and the filtered seismometer
noise, as depicted in Eqn. (3). The measured signals undergo filtering by
complementary filters $H_{1}(s)$ and $H_{2}(s)$, respectively. Subsequently,
they are combined to form a super sensor reading $Y_{\mathrm{super}}(s)$:

$Y_{\mathrm{super}}(s)=\left[H_{1}(s)+H_{2}(s)\right]X(s)+N_{\mathrm{super}}(s)\,,$
(4)

where the sensing noise of the super sensor $N_{\mathrm{super}}(s)$ is given
by

$N_{\mathrm{super}}(s)=H_{1}(s)N_{1}(s)+H_{2}(s)N_{2}(s)\,.$ (5)

For the super sensor to accurately represent an undistorted displacement
signal X(s), the complementary filters must satisfy

$H_{1}(s)+H_{2}(s)=1\,,$ (6)

hence termed complementary filters. This constraint implies that one filter is
fully determined by the other, e.g., $H_{2}(s)=1-H_{1}(s)$. Consequently, only
one filter needs to be designed, and the other is automatically constrained.

The expression for the super sensor noise $N_{\mathrm{super}}(s)$ is
determined by Eqn. (5). Minimizing $N_{\mathrm{super}}(s)$ is crucial as it
dominates the output of the feedback control system at high loop gains.
Substituting Eqn. (6) into Eqn. (5) yields
$N_{\mathrm{super}}(s)=H_{1}(s)N_{1}(s)+[1-H_{1}(s)]N_{2}(s)$, indicating that
the super sensor noise depends solely on filter $H_{1}(s)$. Achieving minimum
$N_{\mathrm{super}}(s)$ entails ensuring that the spectrum $\lvert
N_{\mathrm{super}}(j\omega)\rvert$ converges to the lower of $\lvert
N_{1}(j\omega)\rvert$ or $|N_{2}(j\omega)\rvert$, given the wide dynamic range
between $N_{1}(s)$ and $N_{2}(s)$. In Sec. III.1, we will delve into
approaches for designing the ideal filter shape to meet this objective.

##  III Solving the seismic isolation control problems

###  III.1 Solving the complementary filter problem

The complementary filter problem outlined in Sec. II.3 addresses the
frequency-dependent specification of the sensing noise of the super sensor,
denoted as $N_{\mathrm{super}}(s)$. The objective is to shape a complementary
filter, either $H_{1}(s)$ or $H_{2}(s)$, in order to closely match the
frequency spectrum of $N_{\mathrm{super}}(s)$ to the minimum of the two sensor
noises, $\min(\lvert N_{1}(j\omega)\rvert,\,\lvert N_{2}(j\omega)\rvert)$.

The H-infinity method is commonly employed in a mixed-sensitivity setting to
optimize feedback controllers, ensuring that the closed-loop sensitivity
functions satisfy specific frequency-dependent weightings [28]. This method
proves to be ideal for optimizing complementary filters for active seismic
isolation systems in gravitational-wave detectors [30].

This section sets the stage for the complementary filter problem within an
H-infinity optimization framework and demonstrates that the resulting
H-infinity optimal complementary filters effectively address the filter
shaping challenge. In conventional filter shaping scenarios, control filters
are typically designed based on specified requirements such as target
amplitude or attenuation. Previous research has shown that the H-infinity
method can shape complementary filters to meet certain attenuation
requirements [29]. However, a limitation of this approach is that the
attenuation specifications, being user-defined preferences, do not inherently
optimize noise performance. To overcome this limitation, the method was
further refined in subsequent research to incorporate noise spectrum
specifications into the cost function, facilitating noise minimization [30].
In this section, we revisit the H-infinity optimization method outlined in
previous studies, and in Sec. III.2, we extend this approach to address sensor
correction and feedback control problems.

Figure 6: Generalized plant representation of the complementary filter
problem. $\mathbf{P}(f)$ is the generalized plant, $\mathbf{K}$ is the
feedback controller of the generalized plant, $w_{1}$ and $w_{2}$ are input
signals, $W_{1}(s)$ and $W_{2}(s)$ are weighting functions, $\hat{N}_{1}(s)$
and $\hat{N}_{2}(s)$ are noise models, $H_{1}(s)$ is a complementary filter,
and $z$ is the output error signal.

The H-infinity method serves as an optimal control technique for minimizing a
cost function, specifically the H-infinity norm, through H-infinity synthesis.
To utilize the H-infinity method effectively, the cost function must be
formulated in the form of a generalized plant representation. Fig. 6
illustrates the generalized plant corresponding to the complementary filter
problem, delineating the interconnections of the generalized plant
$\mathbf{P}(s)$ and its feedback controller $\mathbf{K}(s)$. The generalized
plant $\mathbf{P}(s)$ represents a system with two inputs, $w_{1}$ and
$w_{2}$, and one output, $z$, while the controller $\mathbf{K}(s)$ is a system
with one input and one output, represented as $\mathbf{K}(s)=[H_{1}(s)]$,
where $H_{1}(s)$ denotes the complementary filter subject to optimization.
It’s important to note that the generalized plant $\mathbf{P}(s)$ and the
controller $\mathbf{K}(s)$ in this context should not be confused with the
physical plant $P(s)$ and the feedback controller $K(s)$ in the feedback
control system described in Sec. II.1.

The generalized plant $\mathbf{P}(s)$ can be represented by a transfer matrix

$\mathbf{P}(s)=\begin{bmatrix}0&W_{2}(s)\hat{N}_{2}(s)&1\\\
W_{1}(s)\hat{N}_{1}(s)&-W_{2}(s)\hat{N}_{2}(s)&0\end{bmatrix}\,.$ (7)

Here, $W_{1}(s)$ and $W_{2}(s)$ represent weighting functions, to be
elaborated on later in this section, and $\hat{N}_{1}(s)$ and $\hat{N}_{2}(s)$
denote noise models. The magnitude responses of the noise models
$\hat{N}_{1}(s)$ and $\hat{N}_{2}(s)$ are tailored to reflect the amplitudes
of the sensor noises $N_{1}(s)$ and $N_{2}(s)$ in the frequency domain, such
as amplitude spectral densities. Various methods can achieve such shaping, one
of which is detailed in Sec.III.3.

By setting unitary weighting functions $W_{1}(s)=W_{2}(s)=1$ and unity inputs
$w_{1}=w_{2}=1$, the output becomes
$z=H_{1}(s)\hat{N}_{1}(s)+[1-H_{1}(s)]\hat{N}_{2}(s)$, resembling the
expression for the super sensor noise in Eqn. (5). This weighting
specification holds significance when minimizing the super sensor noise in the
root mean square (RMS) sense. Alternatively, the same generalized plant can be
specified as an $\mathcal{H}_{2}$ problem, which aims to minimize the RMS
noise value of the super sensor [28]. However, this approach may not be as
meaningful due to the large dynamic range of the noises, necessitating
suppression across a wide range of magnitudes. While weighting functions can
specify the relative importance of noise suppression at different frequencies,
such specification may render the cost function arbitrary and compromise
optimality. In contrast, in the H-infinity case discussed below, the weighting
functions hold more significance.

The H-infinity synthesis yields an H-infinity optimal controller that
minimizes the H-infinity norm of the closed-loop transfer matrix G(s), given
by:

$\mathbf{G}(s)=\\\
\begin{bmatrix}H_{1}(s)\hat{N}_{1}(s)W_{1}(s)&\left[1-H_{1}(s)\right]\hat{N}_{2}(s)W_{2}(s)\end{bmatrix}\,.$
(8)

The H-infinity norm is defined as:

$\begin{split}\left\lVert\mathbf{G}(s)\right\rVert_{\infty}&=\sup_{\omega}\bar{\sigma}(\mathbf{G}(j\omega))\\\
&=\sup_{\omega}\sqrt{\lambda(\mathbf{G}(j\omega)\mathbf{G}(j\omega)^{H})}\,,\end{split}$
(9)

where $\sup(\cdot)$ denotes the supremum, $\bar{\sigma}(\cdot)$ denotes the
maximum singular value, $G(j\omega)^{H}$ is the conjugate transpose of
$G(j\omega)$, and $\lambda(\cdot)$ denotes the eigenvalue. The eigenvalue
$\lambda(\mathbf{G}(j\omega)\mathbf{G}(j\omega)^{H})$ can be evaluated as a
sum of squares:

$\lambda(\mathbf{G}(j\omega)\mathbf{G}(j\omega)^{H})=\\\ \left\lvert
H_{1}(j\omega)\hat{N}_{1}(j\omega)W_{1}(j\omega)\right\rvert^{2}\\\
+\left\lvert\left[1-H_{1}(j\omega)\right]\hat{N}_{2}(j\omega)W_{2}(j\omega)\right\rvert^{2}\,.$
(10)

Substituting Eqn. (10) into the H-infinity norm (9) reveals that it is the
supremum of a quadrature sum, in the form of $\sqrt{a^{2}+b^{2}}$, where $a$
and $b$ are real values. When $a>b$, the quadrature sum can be approximated as
$\sqrt{a^{2}+b^{2}}=a\sqrt{1+b^{2}/a^{2}}\approx a$. Also, the amplitude
responses of the transfer functions typically differ and have a large dynamic
range. Therefore, when the first term in Eqn. (10) is far greater than the
second term, the H-infinity norm can be approximated as

$\left\lVert\mathbf{G}(s)\right\rVert_{\infty}\approx\sup_{\omega}\left\lvert
H_{1}(j\omega)\hat{N}_{1}(j\omega)W_{1}(j\omega)\right\rvert\,,$ (11)

with at most a multiplicative error of $\sqrt{2}$, when the two terms are
equal, which is negligible in logarithmic scale.

Upon H-infinity optimization, if the value of the H-infinity norm approaches a
value $\gamma$ such that $\lVert\mathbf{G}(s)\rVert_{\infty}=\gamma$, then
Eqn. (11) can be rewritten into an inequality:

$\left\lvert
H_{1}(j\omega)\hat{N}_{1}(j\omega)\right\rvert\leq\gamma\left\lvert
W_{1}(j\omega)\right\rvert^{-1}\,.$ (12)

A similar argument can be made at other frequencies, yielding another
inequality:

$\left\lvert(1-H_{1}(j\omega))\hat{N}_{2}(j\omega)\right\rvert\leq\gamma\left\lvert
W_{2}(j\omega)\right\rvert^{-1}\,.$ (13)

Both inequalities represent the upper bounds of the filtered sensor noises,
and since they hold true at all frequencies, the reciprocals of the weighting
functions $W_{1}(s)$ and $W_{2}(s)$ denote the frequency-dependent
specifications for the filtered sensor noises $H_{1}(s)N_{1}(s)$ and
$H_{2}(s)N_{2}(s)$, respectively.

The H-infinity norm, which is being minimized, represents the deviation from
the target specifications. In other words, H-infinity optimization ensures
that the noise of the super sensor has the least deviation from the target
noise spectrum.

This method allows for choosing the target noise spectrum of the super sensor.
While this can be useful when a specific spectrum can be specified, it may
still not be noise-optimal. Ideally, the super sensor should have sensing
noise as close to the minimum of $\min(\lvert N_{1}(j\omega)\rvert,\,\lvert
N_{2}(j\omega)\rvert)$ as possible. To achieve noise-optimal sensor fusion,
the weighting functions $W_{1}(s)$ and $W_{2}(s)$ should be shaped such that
$\lvert W_{1}(j\omega)\rvert^{-1}=\lvert W_{2}(j\omega)\rvert^{-1}=\min(\lvert
N_{1}(j\omega)\rvert,\,\lvert N_{2}(j\omega)\rvert)$.

However, there’s a more elegant way to achieve the equivalent weighting
without shaping additional transfer functions. When one sensor noise is
greater than the other, say,$\lvert N_{1}(j\omega)\rvert>\lvert
N_{2}(j\omega)\rvert$, the filtered sensor noise dominates the super sensor
noise, i.e., $\lvert N_{\mathrm{super}}(j\omega)\rvert\approx\lvert
H_{1}(j\omega)N_{1}(j\omega)\rvert$. The other sensor noise becomes the
minimum, i.e., $\min(\lvert N_{1}(j\omega)\rvert,\,\lvert
N_{2}(j\omega)\rvert)=\lvert N_{2}(j\omega)\rvert$. In this case, the super
sensor is noise-optimal when the filtered noise $\lvert
H_{1}(j\omega)N_{1}(j\omega)\rvert$ is required to meet the minimum, $\lvert
N_{2}(j\omega)\rvert$. Consequently, the weighting functions can be specified
as:

$W_{1}(s)=\frac{1}{\hat{N}_{2}(s)}\,$ (14)

and

$W_{2}(s)=\frac{1}{\hat{N}_{1}(s)}\,.$ (15)

With these weight specifications, the minimum value of the H-infinity norm is
$\gamma=1$, indicating that the H-infinity optimal complementary filters give
the ideal super sensor with minimum sensor noise. As later shown in Sec. IV.3,
the H-infinity norms obtained are $\gamma=1.033$ and $\gamma=3.5$ for the
sensor correction and complementary filter problems, respectively.

Once the weighting functions are specified, the generalized plant
$\mathbf{P}(s)$, in the form of the transfer matrix in Eqn. (7), is fully
defined. The H-infinity problem can then be solved using available control
system packages such as Python Control [1] and Matlab robust control toolbox.
The `control.hsynthesis()` function wraps around the `SB10AD` Fortran
subroutine, which is a implementation of the solutions in Refs. [19, 20, 26].

In this work, an open-source Python library called Kontrol [32] was written,
which wraps around Python Control. The Kontrol library provides a class called
`kontrol.ComplementaryFilter` with attributes `noise1`, `noise2`, `weight1`,
and `weight2`, which correspond to $\hat{N}_{1}(s)$, $\hat{N}_{2}(s)$,
$W_{1}(s)$, and $W_{2}(s)$, respectively. The method `hinfsynthesis()` in the
Kontrol library wraps around `control.hinfsyn()` from Python Control and can
be used to solve the H-infinity problem, outputting two H-infinity optimal
complementary filters $H_{1}(s)$ and $H_{2}(s)=1-H_{1}(s)$.

###  III.2 Solving the sensor correction and feedback control problems

The feedback control problem is analogous to the complementary filter problem.
Both are spectrum minimization problems: one minimizes the displacement of the
isolation platform, and the other minimizes the measurement noise of the super
sensor.

Comparing Eqn. (1) and Eqn. (6), we see that the two problems are equivalent
by making the substitutions $N_{1}(s)\to D(s)$ and $N_{2}(s)\to-N(s)$. The
weighting functions can be defined similarly as the inverses of the
disturbance and sensing noise transfer function models. Solving the H-infinity
problem with such substitutions yields the optimal suppressions $1/(1+L(s))$
and $L(s)/(1+L(s))$. Then, the H-infinity optimal feedback controller can be
derived using:

$K(s)=\frac{1}{P(s)}\frac{1-S(s)}{S(s)}\,.$ (16)

Here, $S(s)$ is the sensitivity function $1/(1+L(s))$ obtained from the
H-infinity synthesis (analogous to the optimal complementary filter $H_{1}(s)$
in Sec. III.1), and $P(s)$ is the transfer function from the actuation signal
to the displacement of the isolation platform.

Alternatively, the feedback control problem can be solved using a signal-based
mixed-sensitivity configuration, which outputs the optimal controller $K(s)$
directly without needing further conversion. This would require users to
define a slightly different generalized plant, with the same weighting
functions and noise models, while the resulting optimal controller is almost
identical to the aforementioned method. The alternative method is provided in
Ref. [31] and is not detailed here.

The sensor correction problem is slightly different from the feedback control
problem and the complementary filter problem due to the presence of an ambient
term $N_{\mathrm{rel}}(s)$ in the error function (3). This ambient term
defines the minimum noise of the sensor-corrected relative displacement sensor
differently, making the weighting functions in Eqn. (14) and Eqn. (15)
unsuitable for optimizing sensor correction filters. Using a similar
specification for sensor correction would result in suboptimal sensor
correction filters that oversuppress the seismometer noise and the seismic
noise, causing unnecessary noise amplification while leaving the ambient
relative sensor noise dominant.

Recall that the specification in Eqn. (14) and (15) is an alternative way to
specify a lower boundary for the super sensor noise, i.e., $\lvert
W_{1}(j\omega)\rvert^{-1}=\lvert W_{2}(j\omega)\rvert^{-1}=\min(\lvert
N_{1}(j\omega)\rvert,\lvert N_{2}(j\omega)\rvert)$. The sensor correction
problem can be solved by identifying the correct lower boundary. In this case,
it would be $\max[\lvert N_{\mathrm{rel}}(j\omega)\rvert,\min(\lvert
N_{\mathrm{seis}}(j\omega)\rvert,\lvert X_{g}(j\omega)\rvert)]$, which is
either the relative sensor noise or the minimum between the seismometer noise
and the seismic noise, whichever is higher. Once the correct minimum spectrum
is identified, the sensor correction problem can be solved by setting the
weightings $W_{1}(s)$ and $W_{2}(s)$ to the reciprocal of the model of that
minimum, and by substituting $\hat{N}_{1}(s)$ and $\hat{N}_{2}(s)$ with the
corresponding models of the seismometer noise and the seismic noise spectra.

###  III.3 Shaping noise models

To utilize the H-infinity method described in Sec. III.1, the frequency
spectrum of the sensor noises$N_{1}(s)$ and $N_{2}(s)$ must be modeled as the
magnitude responses of the transfer functions $\hat{N}_{1}(s)$ and
$\hat{N}_{2}(s)$, respectively. This can be achieved by fitting a transfer
function to the frequency domain data through a cost function minimization,
similar to any fitting problems. An example cost function could be:

$J(\bm{\theta})=\frac{1}{M}\sum_{\omega}\left[\log\lvert
N(j\omega)\rvert-\log\lvert\hat{N}(j\omega;\bm{\theta})\rvert\right]^{n}\,.$
(17)

Here, $\bm{\theta}$ is an array of parameters that defines the noise model
$\hat{N}(s,\bm{\theta})$, $N(j\omega)$ is the frequency spectrum to be fitted,
$M$ is the number of frequency bins, and $n$ is a positive integer.

To solve the minimization problem, the transfer function must be parameterized
so that an optimization algorithm can iteratively evaluate the cost function
with different parameters. There are two ways to parameterize a fixed-order
transfer function: zero-pole-gain (ZPK) model and ratio of polynomials. The
ZPK model is represented as:

$T_{\mathrm{zpk}}(s)=k\frac{\prod_{i}s-z_{i}}{\prod_{i}s-p_{i}}\,.$ (18)

Here, $k$ is a real constant, and $z_{i}$’s and $p_{i}$’s are the zeros and
poles of the transfer function, respectively. On the other hand, the ratio of
polynomials representation is:

$T_{\mathrm{poly}}(s)=\frac{\sum_{i}b_{i}s^{i}}{\sum_{i}a_{i}s^{i}}\,,$ (19)

and $b_{i}$’s and $a_{i}$’s are the coeffifients of the numerator and
denominator polynomials, respectlively.

Both models can represent all possible transfer functions, but the ZPK model
requires a predetermined structure, i.e., a fixed number of simple and
imaginary poles, since physical transfer functions require complex zeros and
poles to come in conjugate pairs. This puts a limit on the parameter space,
making the ZPK model suboptimal for fitting.

In contrast, the polynomial representation allows the optimization process to
scan the full parameter space. However, this can be numerically challenging as
the cost function may not be bounded, e.g., transfer functions having right-
half-plane poles. Therefore, polynomial models cannot be used to fit frequency
domain data directly. To overcome these problems, we propose a 2-step method
(with an optional step) to acquire a transfer function model.

It is important to note that the H-infinity method requires the plant to be
proper and have a bounded magnitude response. This means that the weighting
functions and the noise models cannot contain integrators, i.e., no poles at 0
Hz. If Eqn. (14) and Eqn. (15) are used to specify the weighting functions,
this further restricts the transfer function models to contain no zeros at 0
and have flat responses at both ends of the spectrum. In other words, the
transfer function models are required to have the same number of zeros and
poles while having no zeros and poles at 0 Hz. To enforce this, the frequency
domain data needs to be padded with edge values at both ends of the spectrum.

As an optional intermediate step, it can be useful to model the frequency
spectrum empirically prior to fitting the transfer function. There are several
advantages to using an empirical model. Firstly, the empirical model may have
fewer parameters compared to the transfer function models due to fractional-
order frequency dependency. For example, the amplitude spectral densities of
the relative sensor noise and inertial sensor noises can be well approximated
by the model:

$N_{\mathrm{model}}(f)=\left[\left(\frac{n_{a}}{f^{a}}\right)^{2}+\left(\frac{n_{b}}{f^{b}}\right)^{2}\right]^{\frac{1}{2}}\,.$
(20)

This model contains only 4 real-valued parameters ($n_{a}$, $n_{b}$, $a$, and
$b$), while the required transfer functions approximating these noises
typically contain 6 or more parameters. This makes fitting the frequency
spectrum measurement, which can contain a large amount of data points,
relatively easy.

Secondly, the empirical model can be treated as a denoised spectrum to remove
unwanted features in the data. Also, spectra are often measured with linear
frequency spacing. When the frequency axis is evenly spaced, the data clusters
at the high-frequency space in a logarithmic scale, reducing the relative
importance of low-frequency data points. With the empirical model, users have
the option to evaluate the model with a logarithmic-spaced frequency axis,
which normalizes the relative weighting between low and high-frequency data.
Moreover, the spectrum can be represented by a few number of data points,
which can largely reduce the cost function evaluation time. All of this helps
in acquiring a more accurate transfer function with a magnitude response
representing a frequency spectrum easily.

While fitting a transfer function in polynomial form is not numerically
stable, this is not the case with the zero-pole-gain (ZPK) model. The simple
zeros and simple poles in a ZPK model represent the corner frequency at which
the magnitude response changes slope. The minimization process can be made
numerically stable if the zeros and poles are bounded to the left-hand side of
the complex plane. And, since the spectrum is flat at both ends of the
frequency axis, the zeros and poles must locate within the frequency band of
interest, which defines the parameter boundary for the zeros and poles. The
gain in the model is also closely related to the DC value of the spectrum, and
that can be relatively easy to be bounded or obtained directly. This allows
the fitting process to be done by minimizing the cost function Eqn. (17) via a
global optimization algorithm, such as differential evolution. Upon
optimization, a transfer function with only real zeros and poles is obtained,
and its magnitude response should coarsely follow the spectrum to be modeled.

As previously mentioned, the ZPK model does not allow the optimization process
to search within the full parameter space, and the obtained ZPK model is not
necessarily the best fit. Instead of using the ZPK model as a noise model, it
can be expanded into a ratio of polynomials. These coefficients correspond to
a transfer function model that has a magnitude response already close to
optimal, which makes them perfect as the initial parameters for local
minimization. In this case, the local minimization should result in a transfer
function model that has a better fit than the initial ZPK model as the
parameter space gets expanded. Also, the local optimization serves the purpose
of fine-tuning the coefficients without largely altering them, reducing the
risk of reaching numerical instability. This way, it becomes possible to fit a
transfer function to frequency domain data using common minimization
algorithms without worrying about numerical stability. While a transfer
function obtained this way may contain right-half-plane zeros, they can be
eliminated by flipping them along the imaginary axis, which only changes the
phase response but not the magnitude response of the transfer function,
ensuring it remains the best fit to the data.

The procedures outlined above allow one to determine the order, i.e., the
number of zeros and poles, of the transfer function model to be fitted.
Initially, the selection of the system order may seem arbitrary, but it can be
iteratively determined through the fitting process. For signals with
fractional-order frequency dependency in the amplitude spectrum, it may
require infinitely many zeros and poles to fully represent the signal in
frequency space. However, adjusting the target spectrum with flat ends can
help practically represent the fractional-order signal with a relatively low
amount of zeros and poles.

To approach the optimal system order, one can repeat the transfer function
fitting process with increasing system order. The minimum order can be
determined from the maximum slope of the spectrum. For example, if the
spectrum has an $f^{-3.5}$ frequency dependency, then rounding up, it follows
that a minimum order of 4, corresponding to 4 zeros and 4 poles, is required
to achieve that slope. From here, one can repeat the fitting process with
increasing order. This iteration does not necessarily continue indefinitely
since redundant zero-pole pairs get placed out of the frequency band of
interest and potentially cancel each other out. Placing zeros and poles
outside the frequency band has minimal effect on the cost function and are
therefore unnecessary. The presence of out-of-bounds zeros and poles can be
used as a termination signal for the iteration, and this way, the optimal
system order can be determined.

There are a few techniques that can be employed to assist the fitting process
and make it more efficient. When fitting with the ZPK model, which relies on a
global optimization scheme, changing the random seed can be useful if the
minimization does not converge into a global minimum. This can be more
effective than directly increasing the system order, as it reduces the
complexity of the final H-infinity filter. Another useful technique is
parameter scaling. For instance, instead of directly optimizing a model
parameter $x$, another parameter $y=\log x$ can be optimized, and the model
parameter is obtained via $x=10^{y}$. This is particularly useful for
parameters that scale logarithmically or have a large dynamic range, such as
the coefficients of the polynomial model and the value of zeros and poles.
Parameter scaling helps evenly space the initial agents, especially in the
case of differential evolution, in the correctly scaled parameter space and
achieves faster convergence as fewer minimization iterations are needed.

It’s worth mentioning that the procedures outlined above are generally useful
for applications that require transfer functions to be optimized via a
mathematical optimization process and are not limited to spectrum fitting.
These procedures are available in the `Kontrol` package via the function
`kontrol.curvefit.spectrum_fit()`. The function takes four mandatory
arguments: the independent variable `xdata`, the dependent variable `ydata`,
and the numbers of zeros and poles `nzero` and `npole`. For the purpose of
H-infinity synthesis, the independent variable and the dependent variable are
simply the frequency axis and the amplitude spectrum. The numbers of zeros and
poles are required to be the same, which determines the system order.

###  III.4 Making practical H-infinity optimal filters

####  III.4.1 Integration drift

The H-infinity optimal filters can sometimes be impractical and unsuitable for
real-time control systems. As discussed in Sec. III.3, each weighting function
and noise model specified in the generalized plant has an equal number of
zeros and poles, resulting in magnitude responses with definite values at both
zero and infinite frequencies, making them flat at both ends. While H-infinity
optimal controllers match the number of states of the specified generalized
plant, they must also exhibit flat responses at both ends of the spectrum.
Consequently, H-infinity filters lack true roll-off at high and low
frequencies, which can be problematic, especially when filtering inertial
sensors prone to integration drift.

Integration drift in inertial sensors leads to progressively larger errors in
displacement due to small errors in measured velocity or acceleration being
integrated over time. With position control activated, the isolation platform
drifts over time, and the integrated signal in the digital system may
saturate, leading to numerical issues. Thus, inertial sensors necessitate true
low-frequency roll-off beyond the frequency band synthesized for by H-infinity
filters.

However, incorporating additional filters along the sensing path may introduce
complications for feedback control, particularly when used in conjunction with
sensor fusion. If the prefilter $H_{\mathrm{prefilt}}(s)$ is employed to
prefilter the inertial sensor in a sensor fusion configuration alongside
complementary filters $H_{1}(s)$ and $H_{2}(s)$, the super sensor reads
$(H_{1}(s)+H_{\mathrm{prefilt}}(s)H_{2}(s))X(s)$, where
$H_{1}(s)+H_{\mathrm{prefilt}}(s)H_{2}(s)\neq 1$ and $X(s)$ is the desired
displacement signal. Consequently, the super sensor exhibits a non-unity
frequency response, distorting the displacement signal. Failure to compensate
for this sensor response can result in an unexpected plant, potentially
leading to excessive gain peaking or instability when the super sensor is used
for feedback control.

Various methods exist to normalize the frequency response of the super sensor.
One approach is to invert the unwanted response from the super sensor by
applying a compensation filter
$H_{\mathrm{comp}}(s)=(H_{1}(s)+H_{\mathrm{prefilt}}(s)H_{2}(s))^{-1}$ to the
super sensor readout, thereby recovering the displacement signal $X(s)$.
Alternatively, the complementary filters can be modified to address the issue.
For instance, one can redefine the complementary filters as

$H_{1}^{\prime}(s)=\frac{H_{1}(s)}{H_{1}(s)+H_{\mathrm{prefilt}}(s)H_{2}(s)}\,,$
(21)

and

$H_{2}^{\prime}(s)=\frac{H_{\mathrm{prefilt}}(s)H_{2}(s)}{H_{1}(s)+H_{\mathrm{prefilt}}(s)H_{2}(s)}\,.$
(22)

as proposed in Ref. [34]. Alternatively, a simpler approach involves
redefining the other complementary filter set as

$H_{2}^{\prime\prime}(s)=H_{\mathrm{prefilt}}(s)H_{2}(s)\,,$ (23)

and

$H_{1}^{\prime\prime}(s)=1-H_{2}^{\prime\prime}(s)\,.$ (24)

Both sets of modified complementary filters can address the signal distortion
issue but may exacerbate noise performance due to the inevitable waterbed
effect. The set described in Eqn. (21) and (22) aims to counteract the
distortion arising from prefiltering inertial sensors. Typically, the
distortion caused by prefiltering manifests as a notch around the blending
frequency. Consequently, the compensation filter exhibits an anti-notch
behavior, resulting in a resonance-like magnitude response. When used as a
compensation filter for the complementary filters, this anti-notch filter
locally amplifies noise around the blending frequency, potentially worsening
the noise RMS performance.

The second approach, utilizing modified complementary filters as depicted in
Eqn. (23) and (24), could offer a more feasible solution. The new high-pass
complementary filter $H_{2}^{\prime\prime}(s)$ maintains similar magnitude
properties to the H-infinity optimal filter $H_{2}(s)$ within the frequency
band of interest, continuing the roll-off beyond frequencies of interest where
the H-infinity filter becomes flat. Consequently, $H_{2}^{\prime\prime}(s)$
remains H-infinity optimal within the frequency band of interest. The
magnitude response of the low-pass complementary filter
$H_{1}^{\prime\prime}(s)\equiv 1-H_{2}^{\prime\prime}(s)$ remains largely
unchanged below the blending frequency. However, the waterbed effect becomes
relevant at higher frequencies where the magnitude is already significantly
lower than unity, resulting in minimal deterioration in noise RMS performance.
Nonetheless, $H_{1}^{\prime\prime}(s)$ is not H-infinity optimal, and the
filtered noise $H_{1}^{\prime\prime}(s)N_{1}(s)$ may exhibit slightly higher
levels compared to the H-infinity optimal case in certain frequency regions.
Consequently, this approach is adopted in this study.

The extent of noise performance deterioration is contingent upon how the
frequency band of interest is defined, as it determines the cutoff frequency
of the high-pass prefilter. Generally, noise amplification is less pronounced
when the cutoff frequency of the prefilter is further away from the blending
frequency. However, extending the frequency band of interest to much lower
frequencies poses numerical challenges for noise modeling and H-infinity
synthesis. Therefore, noise performance degradation resulting from
prefiltering is unavoidable and necessitates experimentation to achieve
satisfactory results. In this paper, noise spectra are modeled between
$3\times 10^{-3}$ and $10^{2}$ Hz.

####  III.4.2 Numerical issues

The H-infinity problem can only be solved numerically, and certain numerical
issues may render H-infinity controllers impractical for implementation. The
Python function `control.hinfsyn()` might generate controllers with
coefficients that are insignificantly small, often several orders of magnitude
smaller than other coefficients. Consequently, such controllers may contain
zeros or poles positioned at extremely high frequencies, typically beyond the
Nyquist frequency of the digital control system. While this issue is
frequently overlooked during offline noise performance analysis in Python, it
becomes apparent when implemented in real-time digital systems with finite
sampling rates, leading to misinterpretation and rendering the filter unusable
[31]. One solution to this problem involves zeroing out the negligibly small
coefficients. By doing so, outlier zeros and poles are automatically removed
without altering the frequency response below the Nyquist frequency. This
cleanup process can be performed manually or by utilizing Kontrol’s
`ComplementaryFilter.hinfsynthesis()` method, which internally compares
coefficients and automatically removes outliers.

##  IV Results: KAGRA Type-B seismic attenuation system

Figure 7: KAGRA Type-B Seismic Attenuation System.

In this section, the H-infinity method is applied to optimize sensor
correction filters and complementary filters for an active seismic isolation
system in KAGRA. Fig. 7 illustrates a Type-B suspension utilized to suspend
the beamsplitter (BS) and the signal recycling mirrors (SR2, SR3, and SRM) of
the KAGRA interferometer [15]. This suspension comprises a quintuple pendulum
arrangement, with the first pendulum stage serving as the isolation table,
functioning as an inverted pendulum, and the suspended mirror acting as the
fifth pendulum stage. Equipped with 3 L-4C geophones from Sercel, they measure
the inertial motion of the isolation table, while 3 custom linear variable
differential transformers (LVDTs) measure the relative displacement between
the isolation table and the ground-attached base. Additionally, a Trillium
Compact seismometer placed nearby monitors ground motion for sensor correction
purposes.

To apply the H-infinity method, the seismic noise and sensor noises from the
inertial sensor, relative sensor, and seismometer must be modeled. These noise
models are derived using the procedures outlined in Sec. III.3, utilizing a
dataset obtained at 9:00 a.m. UTC on April 26th, 2023. The dataset comprises
free spectra of the isolation table and simultaneously measured ground
displacement spectra. Further details regarding the modeling procedures are
elaborated in subsequent sections.

Initially, the noise models are employed to optimize a sensor correction
filter, as described in Sec. IV.1, yielding a sensor-corrected relative
sensor. Subsequently, in Sec. IV.2, the sensing noise of the sensor-corrected
relative sensor is modeled and utilized to optimize complementary filters,
integrating the sensor-corrected relative sensor and the inertial sensor on
the isolation table. The combination of the optimized sensor correction and
sensor fusion configurations results in a super sensor exhibiting superior
noise properties.

Furthermore, the optimized sensor correction and complementary filters
underwent testing in a real-time system at 5:00 a.m. UTC on September 12th,
2023, under different seismic conditions, serving as validation. The results
of these validation tests are presented in Sec IV.3.

###  IV.1 Sensor correction of a linear variable differential transformer

####  IV.1.1 Seismic noise, seismometer noise, and relative sensor noise
models

In Fig. 8, the seismometer readout (blue solid line), the relative sensor
readout (orange solid line), and the amplitude responses of the noise models
for the relative sensor noise (green dashed line), seismometer noise (red
dash-dot line), and seismic noise (purple dotted line) are depicted. The
seismometer captured the secondary microseism, peaking at around 1 $\mu$m
between 0.2 and 0.3 Hz. This seismic noise is also evident in the relative
sensor readout. This is the primary noise component targeted for removal by
sensor correction. Additionally, a peak is observed around 0.08 Hz,
corresponding to the resonance frequency of the isolation table due to its
inverted pendulum configuration.

Figure 8: Amplitude spectral densities of the seismometer readout, the
relative sensor readout, and the noise models taken at KAGRA along the Y-arm
(horizontal) direction.

The seismometer noise dominates the readout at low frequencies. To model
seismometer noise, only data below 0.06 Hz is utilized. To obtain the seismic
noise model, the seismometer spectrum is extended below 0.08 Hz with the local
minima value to remove seismometer noise and is extended with above 0.2 to 0.3
Hz with local maxima value. The latter extension ensures seismic noise
reduction above the secondary microseism, effectively transforming the sensor
correction filter into a high-pass filter, consistent with its traditional
design [23]. Although suboptimal, this approach is necessary because seismic
noise is lower than the relative sensor noise at higher frequencies.
Amplifying seismic noise at higher frequencies could lead to degradation of
noise RMS performance, particularly in anthropogenic and acoustic bands where
seismic noise varies significantly throughout the day. However, since the
sensor correction filter is static, it should be designed in a way such that
environmental changes can be handled sufficiently. Therefore, padding the
spectrum at frequencies higher than the microseismic peak is a reasonable
augmentation to maintain high seismic attenuation. As for the relative sensor,
the noise is modeled below 0.03 Hz and above 0.8 Hz, neglecting resonances and
ground coupling.

Following the empirical modeling, ZPK modeling, and transfer function modeling
procedures outlined in Sec. III.3, transfer function models for relative
sensor noise, seismometer noise, and seismic noise are obtained. The amplitude
responses of these models are depicted as green dashed, red dash-dot, and
purple dotted lines in Fig. 8, respectively. The empirical models, which
remain unflattened outside the defined frequency band of interest, spanning
from 0.003 Hz to 100 Hz, are useful for estimating the performance of the
optimized sensor correction filter later on.

####  IV.1.2 H-infinity optimal sensor correction filter

Figure 9: Magnitude responses of the H-infinity sensor correction filters and
its complement, which denotes the seismic transmissivities.

Using the transfer function models depicted in Fig. 8, the sensor correction
filter can be optimized utilizing the H-infinity method outlined in Sections
III.1 and III.2. The seismometer noise model and the seismic noise model are
specified as the noise models $\hat{N}_{1}(s)$ and $\hat{N}_{2}(s)$,
respectively. Weighting functions are employed to denote the target noise
levels of the seismometer noise and the seismic noise. These weighting
functions are defined as the maximum between the relative sensor noise and the
minimum between the seismic and seismometer noises. Given that the minimum
noise level can be approximated as simply the relative sensor noise,
$\max[|N_{\mathrm{rel}}(j\omega)|,\min(|N_{\mathrm{seis}}(j\omega)|,|X_{g}(j\omega)|)]\approx|N_{\mathrm{rel}}(j\omega)|$,
as it is higher than the lower bound between the seismic and seismometer
noises, the weighting functions (14) and (15) are set to be the reciprocal of
the relative sensor noise model.

With the noise model and weighting function specifications in place, the
sensor correction filter is optimized using the `hinfsynthesis()` method in
Kontrol. Fig. 9 illustrates the magnitude responses of the optimized sensor
correction filter and its complement, i.e., the seismic transmissivity. In the
figure, the solid blue and orange lines represent the magnitude responses of
the H-infinity optimal sensor correction filter and the seismic
transmissivity, respectively. However, the sensor correction filter, optimized
with respect to the specified noise models and weights, is only well-defined
within the frequency range of 0.003 Hz to 100 Hz and does not roll off below
0.003 Hz, as indicated in the figure.

To address this, an additional 3rd-order high-pass filter with a cutoff
frequency at 0.003 Hz is introduced as a prefilter to achieve continuous roll-
off at low frequencies for the sensor correction filter. Consequently, the
sensor correction filter with prefiltering now exhibits roll-off at low
frequencies, thereby providing appropriate seismometer noise attenuation. As a
tradeoff, the seismic transmissivity increases by a negligible amount,
particularly across 0.2 and 1 Hz.

####  IV.1.3 Estimated sensor correction performance

Figure 10: Amplitude spectral densities of the seismometer readout,
uncorrected relative sensor noise, and estimated sensing noise of the sensor-
corrected relative sensor.

Fig. 10 illustrates the estimated noise level of the sensor-corrected relative
sensor, along with the uncorrected relative sensor noise and the seismometer
readout. Here, the uncorrected relative sensor measures the relative motion
between the isolation table and the ground. Therefore, the ground motion is
included as a sensing noise of the uncorrected relative sensor. Both the
uncorrected and corrected relative sensor noises utilize noise models of the
relative sensor and the seismic noise that reflects the actual frequency
content of the seismic noise.

The minimum H-infinity norm of this problem is $\gamma=1.033$, indicating that
the optimized sensor correction filter suppresses seismic and seismometer
noises to at most 1.033 times the amplitude of the relative sensor noise.
Consequently, the sensing noise of the sensor-corrected relative sensor should
be almost indistinguishable from the intrinsic noise of the relative sensor,
as demonstrated in Fig. 10.

To assess the sensor correction performance, the noise performances of the
uncorrected and corrected relative sensors are evaluated in band-limited RMS
values in displacement and velocity around 0.1 and 0.5 Hz. The band-limited
noise RMS values are 0.160 $\mu\mathrm{m}$ and 0.0367 $\mu\mathrm{m}$ in
displacement for the uncorrected and corrected relative sensors, respectively.
Similarly, the band-limited noise RMS values in velocity are 0.265
$\mu\mathrm{m}$ and 0.0547 $\mu\mathrm{m}$ for the uncorrected and corrected
relative sensors, respectively.

It is evident that the H-infinity optimal sensor correction filter provides a
4.36 times and 4.84 times reduction in RMS displacement and velocity,
respectively, around the secondary microseism peak. Additionally, there is
approximately a 7 times peak attenuation at 0.265 Hz, where the secondary
microseism peak resides.

###  IV.2 Sensor fusion of a sensor-corrected linear variable differential
transformer and a geophone

The sensor-corrected relative sensor introduced in Sec. IV A is already
suitable for active damping and seismic isolation. When paired with a high-
gain controller, it offers control performance up to the noise level of the
corrected relative sensor, which is approximately
$10^{-1}\,\mu\mathrm{m}/\sqrt{\mathrm{Hz}}$ at 0.2 Hz, as indicated in Fig.
10. In contrast, inertial sensors typically exhibit better noise performance
above 0.1 Hz.

Therefore, combining the inertial sensor with the corrected relative sensor
has the potential to enhance control performance above 0.1 Hz, leading to
improved active damping and seismic isolation. By leveraging the superior
noise performance of the inertial sensor in higher frequency ranges, the
overall control system can achieve enhanced effectiveness in mitigating
disturbances and maintaining system stability.

####  IV.2.1 Sensor-corrected relative sensor and inertial sensor noise
models.

Figure 11: Amplitude spectral densities of the corrected relative sensor
noise, inertial sensor noise, and their models.

In Fig. 11, the estimated noise of the corrected relative sensor (blue solid)
and the inertial sensor readout (orange solid) are shown, accompanied by the
magnitude responses of their corresponding transfer function noise models
(green dashed and red dash-dot lines). The estimated corrected relative sensor
noise is directly derived from Sec. IV.1.

The inertial sensor reveals features not observable at 0.08 Hz, 0.4 Hz, and
0.65 Hz, corresponding to pendulum modes of the multiple pendulum chain.
Motion induced by microseismic activity on the isolation platform is
observable around 0.2 Hz.

By leveraging the capability of the inertial sensor to detect these modes, it
becomes feasible to actively damp them, thereby reducing the residual motion
of the suspended optics. This is particularly beneficial at higher stages
where control noise is more effectively reduced via passive isolation. In
contrast, control noise due to damping control loops at the lower pendulum
stages are not as profoundly reduced by passive isolation and are subjected
tighter gain and noise requirements.

####  IV.2.2 H-infinity optimal complementary filters

Figure 12: Magnitude responses of the H-infinity complementary filters.

The two noise models in Fig. 11 are the noise models in the generalized plant.
In this case, the weighting functions are specified according to Eqn. (14) and
Eqn. (15), representing the inverse of the noise models. With these
specifications, H-infinity synthesis yields optimal complementary filters, as
depicted in Fig. 12.

In Fig. 12, the blue and orange solid lines represent the magnitude responses
of the optimized low-pass and high-pass filters, respectively. Additionally,
the prefilter-practicalized low-pass and high-pass filters, obtained from Eqn
(21) and Eqn. (22), are shown as green dashed and red dot-dash lines in the
figure. Here, the prefilter is a 4th-order high-pass filter with a cutoff
frequency at 0.003 Hz. Again, the prefilter is required here to convert the
H-infinity high-pass filter into a one that has the required roll-off at below
0.003 Hz.

It is evident from the figure that with the prefilter, the high-pass filter
continues to roll off at frequencies below the frequency band of interest.
However, the low-pass filter experiences noticeable amplification between 0.2
to 100 Hz, as a tradeoff. This amplification is more pronounced compared to
the sensor correction case, primarily due to the higher order and closer
cutoff frequency of the required prefilter relative to the blending frequency.

####  IV.2.3 Estimated sensor fusion performance

Figure 13: Amplitude spectral densities of the estimated sensor-corrected
relative sensor noise, the inertial sensor noise, and the estimated sensing
noise of the super sensor.

Using the estimated sensor-corrected relative sensor noise and the empirical
model of the inertial sensor noise, the predicted noise of the super sensor is
depicted in Fig. 13. Without the prefilter, the complementary filters were
shaped perfectly according to the given noise information, resulting in the
super sensor noise maintaining a similar but minimized deviation from the
lower noise limit at all frequencies. The H-infinity norm corresponding to
this problem is $\gamma=3.5$, indicating that the deviation is also around 3.5
times the lower limit, and the super sensor noise hovers around that amount
above the lower limit. A higher $\gamma$ value indicates stricter
specifications, which is applicable in this case where the relative sensor
noise and inertial sensor noise significantly differ in frequency dependency,
necessitating high-order roll-offs.

With the prefilter, the low-frequency inertial sensor noise gets further
suppressed, notably below 0.01 Hz, while the relative sensor noise gets
amplified around 0.2 to 100 Hz. The super sensor noise, with prefilter, still
generally follows the shape of the lower noise limit but has some parts lower
and others higher than the original super sensor noise.

The overall displacement noise RMS values of the super sensor with and without
prefilter are 0.199 $\mu\mathrm{m}$ and 0.275 $\mu\mathrm{m}$ , respectively.
Additionally, the overall velocity noise RMS values of the super sensor with
and without prefilter are 0.044 $\mu\mathrm{m/s}$ and 0.0429
$\mu\mathrm{m/s}$, respectively. The prefilter reduces the displacement noise
by about 28% with the cost of increasing velocity noise by about 2.5%. This
demonstrates that the prefilter offers a good tradeoff, even with the more
pronounced amplification seen in the low-pass filter around 0.2 to 100 Hz.

The overall displacement RMS value of the sensor-corrected relative sensor
noise is 0.120 $\mu\mathrm{m}$, whereas that for the super sensor noise is
0.199 $\mu\mathrm{m}$. This higher displacement RMS is attributed to the high
inertial sensing noise below 0.1 Hz. Moreover, the band-limited displacement
and velocity RMS values are 0.0188 $\mu\mathrm{m}$ and 0.0150
$\mu\mathrm{m/s}$ for the super sensor noise, offering a 1.95 times and 3.78
times reduction in band-limited displacement and velocity RMS, respectively,
compared to the sensor-corrected relative sensor, which is 0.0367
$\mu\mathrm{m}$ and 0.0547 $\mu\mathrm{m/s}$ as mentioned in Sec. IV.1.3.
Additionally, the overall velocity RMS of the super sensor noise is 0.0440
$\mu\mathrm{m/s}$, much lower than the 8.45 $\mu\mathrm{m/s}$ of the sensor-
corrected relative sensor noise. This indicates that the super sensor offers
superior noise performance, particularly around the secondary microseism, and
a better velocity noise performance, which is crucial for improving the
chances of lock-acquisition of the interferometer.

###  IV.3 Implementation with a KAGRA seismic isolation system

As a validation of the control filters, the H-infinity optimized sensor
correction filter and complementary filters in Sec. IV.1 and Sec. IV.2 were
tested with a KAGRA Type-B seismic attenuation system. With the isolation
platform uncontrolled, the readouts from the sensor-corrected relative sensor
and the super sensor are compared with the seismic noise and the original
relative sensor. These sensors measure the common motion of the table but have
different noise characteristics that can be compared. In this section, the
overall RMS values, the band-limited RMS values around the microseism, and the
amplitude spectral densities of the readouts at 10 Hz are compared.

Figure 14: Amplitude spectral densities of the measured horizontal isolation
platform displacement spectrums and the sensor noise empirical models.

Fig. 14 shows the open-loop displacement of the isolation platform measured by
the uncorrected relative sensor (blue solid), the sensor-corrected relative
sensor (orange solid), and the super sensor (green solid). The seismometer
readout, representing the ground displacement, is also shown in gray.
Additionally, the empirical models of the relative sensor (red dashed) and the
inertial sensor (purple dash-dot) noises are plotted for reference. The
secondary microseism peaks at around $0.2\,\mu\mathrm{m}/\sqrt{\mathrm{Hz}}$,
significantly lower than the one used to optimize the sensor correction and
complementary filters. Therefore, this can be used as a validation test for
the optimized sensor correction and complementary filters under a suboptimal
condition.

Qualitatively, the actual performance of the sensor correction and the sensor
fusion matched the expectation set in Sec. IV.1 and Sec. IV.2. The secondary
microseism around 0.2 Hz is well suppressed by the sensor correction method.
The secondary microseismic peak is not observable in the sensor-corrected
relative sensor readout, and the injected seismometer noise is minimal at
lower frequencies. The readout of the relative sensor is coupled to the motion
of the isolation platform and the seismic noise, with an RMS value of 0.0982
$\mu\mathrm{m}$. With the sensor correction engaged, the readout of the
corrected relative sensor has an RMS value of 0.0879 $\mu\mathrm{m}$. However,
the RMS value of the super sensor readout is 0.117 $\mu\mathrm{m}$, slightly
higher but comparable to the relative sensor readouts, corrected and
uncorrected. This is a result of complementary filters working in a suboptimal
seismic condition, where low-frequency inertial sensor noise was injected as a
tradeoff for excessive seismic attenuation.

To better assess the seismic isolation performance, the band-limited RMS
values of the readouts, around the secondary microseism, from 0.1 to 0.5 Hz
are compared. These values provide insight into how effectively each sensor
mitigates the effects of the microseism. The band-limited RMS values in
displacement units for the uncorrected relative sensor, sensor-corrected
relative sensor, and the super sensor are 0.0500 $\mu\mathrm{m}$, 0.0261
$\mu\mathrm{m}$, and 0.0150 $\mu\mathrm{m}$, respectively. Similarly, the
band-limited RMS values in velocity units show a similar trend, with values of
0.0773 $\mu\mathrm{m/s}$, 0.0403 $\mu\mathrm{m/s}$, and 0.0133
$\mu\mathrm{m/s}$ for the uncorrected relative sensor, sensor-corrected
relative sensor, and the super sensor readouts, respectively. Although these
RMS values do not provide absolute performance measures, they offer a relative
comparison of each sensor’s susceptibility to the secondary microseism. The
super sensor demonstrates the lowest band-limited RMS value, indicating its
superior robustness and suitability for long-term operation, especially in
environments where seismic conditions change gradually. This highlights the
utility of the H-infinity method in optimizing sensor correction and
complementary filters, enabling functional and practical configurations for
active isolation systems in gravitational-wave detectors. Nevertheless, the
successful implementation of the sensor correction and sensor fusion controls
in active isolation systems contributes to the improvement of duty cycle of
KAGRA from 53% during the O3GK observation run [10] to 80% during the O4a
observation run [33].

##  V Discussions and Future Work

###  V.1 Elevated noise RMS values of the super sensor

As discussed in sections IV.2 and IV.3, the super sensor exhibits poorer noise
performance at low frequencies due to high inertial sensing noise, a common
challenge in gravitational-wave detectors’ sensor fusion systems, unrelated to
the H-infinity method. However, despite this elevated low-frequency noise, the
super sensor offers superior attenuation in the secondary microseism frequency
band, where seismic noise is non-stationary. Essentially, the super sensor
demonstrates greater robustness against variations in secondary microseism,
with less coupling to this phenomenon compared to relative sensors.
Consequently, the advantages of employing a super sensor become evident when
secondary microseism levels are heightened.

###  V.2 Problems with cascaded systems

The sensor fusion system involving a sensor-corrected relative sensor operates
as a cascaded control system, where the performance of each control system
depends on the other. However, optimizing individual control systems using the
H-infinity method does not ensure overall instrumentation optimality
throughout the frequency space. This is evident in the sensor fusion setup
described in Sec. IV.2, where the super sensor noise deviates multiplicatively
from the relative sensor noise by the product of their respective H-infinity
norms, 1.033 and 3.5, at low frequencies but only deviates by 3.5 from the
inertial sensor noise at high frequencies. To address this, reformulating the
generalized plant to co-optimize the necessary control filters for the
cascaded system can ensure true H-infinity optimality, potentially mitigating
the elevated noise RMS of the super sensor.

###  V.3 Prefiltering problem

The necessity of prefiltering inertial sensors, combined with the
complementary condition described in Eqn. (6), presents a notable challenge
for the H-infinity method. As outlined in Sections IV.1 and IV.2, the process
of reconfiguring the complementary low-pass filter from a complementary high-
pass filter, with the inclusion of a prefilter, introduces amplification in
the low-pass filter, thereby compromising the optimality of the complementary
filters in the H-infinity framework. Consequently, the low-pass filters
discussed in the preceding sections demonstrate a slight degradation in noise
performance within the frequency range of 0.2 to 100 Hz.

This challenge stems from the inability to fully define the spectrum of
inertial sensor noise, particularly at lower frequencies. This issue mirrors a
well-known research problem in H-infinity methods, where integral actions are
desired in H-infinity controllers. Various solutions documented in
contemporary control theory literature address the optimization of controllers
with integral actions, such as the “Shift with approximate inverse method”
[24]. By adapting these solutions to address sensor fusion problems, the
prefiltering challenge can be mitigated effectively.

###  V.4 Adaptive seismic control

The H-infinity filters optimized for specific seismic conditions, as
demonstrated in Sec. IV.2, may not always remain optimal. Particularly, the
excessive suppression around the secondary microseismic frequency, unnecessary
during periods of low secondary microseism magnitude, led to noise injection
at lower frequencies, resulting in higher noise RMS. However, this is not a
drawback of the H-infinity method but rather a limitation of the current
digital control systems in gravitational-wave detectors, which rely on static
controllers. The assumption of stationary signals in the H-infinity method
conflicts with the non-stationary nature of seismic noise, causing the optimal
filters to appear suboptimal and provide worse performance during seismic
condition changes.

Seismic noise, predominantly RMS-dominated by the primary and secondary
microseismic peaks, exhibits significant variation. While the primary
microseism typically remains unobservable, the secondary microseism,
characterized by its background hum, fluctuates in magnitude and frequency due
to weather and seasonal changes. These fluctuations, coupled with the
transient seismic waves from events like distant earthquakes, challenge the
seismic isolation in gravitational-wave detectors, limiting their duty cycle.
To address this challenge, two solutions, short-term and long-term solutions,
are briefly discussed here.

The short-term solution involves switching control filters to achieve semi-
adaptive control without overhauling the digital control system. Conversely,
the long-term solution necessitates reprogramming the system to implement
time-domain adaptive control methods. The H-infinity method opens up the
opportunity for real-time control filter optimization, but automation
challenges persist, particularly in real-time modeling of the seismometer
noise for sensor correction. Alternatively, a semi-adaptive control approach
could involve regularly switching, similar to a gain scheduling approach,
seismic control filters optimized for different seismic environments,
selecting filters based on a custom criteria or cost function evaluated using
real-time seismic noise measurements.

A permanent adaptive isolation control solution requires reprogramming digital
control systems to implement advanced methods like model predictive control
(MPC) [25]. In MPC, seismic disturbances are forecasted in real-time and
control signals are optimized accordingly, offering adaptiveness within short
time horizons. MPC can complement the existing LTI control system, forming a
cascaded control system that maximizes robustness against seismic
disturbances. Ultimately, the MPC framework represents a long-term development
goal for active isolation systems, with the aim to eliminate observation
downtime in gravitational-wave detectors.

##  VI Conclusions

Gravitational-wave detectors, being large-scale interferometers with precisely
aligned optics, are vulnerable to environmental disturbances like earthquakes.
To counteract this, current detectors use active isolation systems employing
control techniques such as feedback control, sensor correction, and sensor
fusion. Optimizing these control configurations is crucial, and the H-infinity
method proves useful in creating optimal control filters. These filters
provide control performances close to the instrumentation limit.

The H-infinity method has shown effectiveness in optimizing sensor correction
and fusion configurations for suspension systems like the KAGRA Type-B,
resulting in a highly functional super sensor with superior isolation
performance, particularly around secondary microseism frequencies. However,
challenges remain, such as elevated noise RMS in the super sensor, cascaded
system issues, and prefiltering challenges. Improvements to the H-infinity
method are necessary to address these challenges comprehensively.

Seismic noise variation, particularly secondary microseisms, presents a hurdle
for static control systems, which may struggle to robustly handle such
changes. Nevertheless, the H-infinity method can generate sets of control
filters optimized for different seismic conditions. These sets can be employed
in a semi-adaptive control system, where filters are regularly adjusted, akin
to gain scheduling. While alternative approaches like Model Predictive Control
offer full adaptiveness to dynamic changes in seismic noise, their
implementation requires extensive development over the long term.

###### Acknowledgements.

This material is based upon work supported by NSF’s LIGO Laboratory which is a
major facility fully funded by the National Science Foundation. This work was
supported by Ministry of Education, Culture, Sports, Science and Technology
(MEXT), Japan Society for the Promotion of Science (JSPS) in Japan; National
Research Foundation (NRF) and Ministry of Science and ICT (MSIT) in Korea;
Academia Sinica (AS) and National Science and Technology Council (NSTC) in
Taiwan. This work was supported by the Science and Technology Facilities
Council grant number ST/V005618/1.

## References

  * [1] Control systems library for python.  http://github.com/python-control/python-control, 2017.
  * [2] J. Aasi et al.  Enhancing the sensitivity of the LIGO gravitational wave detector by using squeezed states of light.  Nature Photon., 7:613–619, 2013.
  * [3] B. P. Abbott et al.  Observation of gravitational waves from a binary black hole merger.  Phys. Rev. Lett., 116:061102, Feb 2016.
  * [4] B. P. Abbott et al.  Prospects for observing and localizing gravitational-wave transients with Advanced LIGO, Advanced Virgo and KAGRA.  Living Rev. Rel., 21(1):3, 2018.
  * [5] B. P. Abbott et al.  Gwtc-1: A gravitational-wave transient catalog of compact binary mergers observed by ligo and virgo during the first and second observing runs.  Phys. Rev. X, 9:031040, Sep 2019.
  * [6] Benjamin P. Abbott et al.  Sensitivity of the Advanced LIGO detectors at the beginning of gravitational wave astronomy.  Phys. Rev. D, 93(11):112004, 2016.  [Addendum: Phys.Rev.D 97, 059901 (2018)].
  * [7] R. Abbott et al.  Gwtc-2: Compact binary coalescences observed by ligo and virgo during the first half of the third observing run.  Phys. Rev. X, 11:021053, Jun 2021.
  * [8] R. Abbott et al.  Gwtc-3: Compact binary coalescences observed by ligo and virgo during the second part of the third observing run.  Phys. Rev. X, 13:041039, Dec 2023.
  * [9] R. Abbott et al.  Gwtc-2.1: Deep extended catalog of compact binary coalescences observed by ligo and virgo during the first half of the third observing run.  Phys. Rev. D, 109:022001, Jan 2024.
  * [10] H. Abe et al.  Performance of the KAGRA detector during the first joint observation with GEO 600 (O3GK).  PTEP, 2023(10):10A101, 2023.
  * [11] F Acernese et al.  Advanced virgo: a second-generation interferometric gravitational wave detector.  Classical and Quantum Gravity, 32(2):024001, dec 2014.
  * [12] Y. Akiyama et al.  Vibration isolation system with a compact damping system for power recycling mirrors of KAGRA.  Class. Quant. Grav., 36(9):095015, 2019.
  * [13] T. Akutsu et al.  First cryogenic test operation of underground km-scale gravitational-wave observatory KAGRA.  Class. Quant. Grav., 36(16):165008, 2019.
  * [14] T Akutsu et al.  Overview of KAGRA: Detector design and construction history.  Progress of Theoretical and Experimental Physics, 2021(5):05A101, 08 2020.
  * [15] T Akutsu et al.  Vibration isolation systems for the beam splitter and signal recycling mirrors of the kagra gravitational wave detector.  Classical and Quantum Gravity, 38(6):065011, mar 2021.
  * [16] G.J. Balas, J.C. Doyle, K. Glover, A. Packard, and R. Smith.  $\mu$-Analysis and Synthesis Toolbox.  The MathWorks Inc., Natick, MA, 1995.
  * [17] Jonathan J Carter, Samuel J Cooper, Edward Thrift, Joseph Briggs, Jim Warner, Michael P Ross, and Conor M Mow-Lowry.  Particle swarming of sensor correction filters.  Classical and Quantum Gravity, 37(20):205009, sep 2020.
  * [18] The LIGO Scientific Collaboration.  Advanced ligo.  Classical and Quantum Gravity, 32(7):074001, mar 2015.
  * [19] J.C. Doyle, K. Glover, P.P. Khargonekar, and B.A. Francis.  State-space solutions to standard $\mathcal{H}_{2}$ and $\mathcal{H}_{\infty}$ control problems.  IEEE Transactions on Automatic Control, 34(8):831–847, 1989.
  * [20] Keith Glover and John C. Doyle.  State-space formulae for all stabilizing controllers that satisfy and $H_{\infty}$ norm bound and relations to risk sensitivity.  Systems & Control Letters, 11:167–172, 1988.
  * [21] Wensheng Hua, Ranav Adhikari, Dan Debra, Joseph Giaime, Giles Hammond, C Hardham, Michaela Hennessy, Jonathan How, Brian Lantz, M Macinnis, R Mittleman, Sam Richman, Norna Robertson, Jimmy Rollins, David Shoemaker, Robin Stebbins, and Stanford Univ.  Low frequency active vibration isolation for advanced ligo.  Proceedings of SPIE - The International Society for Optical Engineering, 5500, 09 2004.
  * [22] Michele Maggiore, Chris Van Den Broeck, Nicola Bartolo, Enis Belgacem, Daniele Bertacca, Marie Anne Bizouard, Marica Branchesi, Sebastien Clesse, Stefano Foffa, Juan García-Bellido, Stefan Grimm, Jan Harms, Tanja Hinderer, Sabino Matarrese, Cristiano Palomba, Marco Peloso, Angelo Ricciardone, and Mairi Sakellariadou.  Science case for the einstein telescope.  Journal of Cosmology and Astroparticle Physics, 2020(03):050, mar 2020.
  * [23] F. Matichard et al.  Seismic isolation of Advanced LIGO: Review of strategy, instrumentation and performance.  Class. Quant. Grav., 32(18):185003, 2015.
  * [24] Marcio Miranda, Ricardo Takahashi, and Fábio Jota.  H-infinity controllers with integral action: an experimental evaluation.  Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering, 224:21–29, 01 2010.
  * [25] Michael Nikolaou.  Model predictive controllers: A critical synthesis of theory and industrial needs.  Advances in Chemical Engineering, 26:131–204, 2001.
  * [26] P.Hr. Petkov, D.W. Gu, and M.M. Konstantinow.  Fortran 77 routines for Hinf and H2 design of continuous-time linear control systems.  1998\.
  * [27] Peter R Saulson.  Fundamentals of Interferometric Gravitational Wave Detectors.  WORLD SCIENTIFIC, 2nd edition, 2017.
  * [28] Sigurd Skogestad and Ian Postlethwaite.  Multivariable Feedback Control: Analysis and Design.  John Wiley & Sons, Inc., Hoboken, NJ, USA, 2005.
  * [29] Dehaeze Thomas, Mohit Verma, and Collette Christophe.  Complementary filters shaping using $\mathcal{H}_{\infty}$ synthesis.  pages 459–464, 11 2019.
  * [30] T. T. L. Tsang, T. G. F. Li, T. Dehaeze, and C. Collette.  Optimal sensor fusion method for active vibration isolation systems in ground-based gravitational-wave detectors.  Class. Quant. Grav., 39(18):185007, 2022.
  * [31] Terrence Tak Lun Tsang.  Optimizing active vibration isolation systems in ground-based interferometric gravitational-wave detectors, 2022.
  * [32] T.T.L. Tsang.  Kontrol.  https://kontrol.readthedocs.io/, 2020.
  * [33] T Ushiba et al.  Status of the underground gravitational wave detector kagra, 2023.  https://indico.cern.ch/event/1199289/contributions/5449057/attachments/2703561/4692881/taup2023_20230830.pdf.
  * [34] J.V. van Heijningen.  Turn up the bass!: Low-frequency performance improvement of seismic attenuation systems and vibration sensors for next generation gravitational wave detectors.  PhD thesis, Vrije Universiteit Amsterdam, 2018.