# Earthquake Control: An Emerging Application for Robust Control. Theory and Experimental Tests

Diego Gutiérrez-Oribio, Georgios Tzortzopoulos, Ioannis Stefanou, Franck Plestan

<figure><img src="image_1.png" /></figure>

## To cite this version:

Diego Gutiérrez-Oribio, Georgios Tzortzopoulos, Ioannis Stefanou, Franck Plestan. Earthquake Control: An Emerging Application for Robust Control. Theory and Experimental Tests. IEEE Transactions on Control Systems Technology, 2023, 31 (4), pp.1747-1761. 10.1109/tcst.2023.3242431. hal-05232690

HAL Id: hal-05232690

https://hal.science/hal-05232690v1

Submitted on 1 Sep 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# Earthquake Control: An Emerging Application for Robust Control. Theory and Experimental Tests

Diego Gutiérrez-Oribio, Georgios Tzortzopoulos, Ioannis Stefanou and Franck Plestan

**Abstract**— This paper addresses the possibility of using robust control theory for preventing earthquakes through fluid injections in the earth's crust. The designed robust controllers drive aseismically a fault system to a new equilibrium point of lower energy by tracking a slow reference signal. The control design is based on a reduced-order nonlinear model able to reproduce earthquake-like instabilities. Uncertainties related to the frictional and mechanical properties of the underlying physical process and external perturbations are considered. Two types of controllers are derived. The first one is based on sliding-mode theory and leads to local finite-time convergence of the tracking error and rejection of Lipschitz w.r.t. time perturbations. The second controller is based on LQR control and presents global exponential stability of the tracking error and rejection of Lipschitz w.r.t. states perturbations. Both controllers generate a continuous control signal, attenuating the chattering effect in the case of the sliding-mode algorithms. The developed controllers are tested extensively and compared on the basis of numerical simulations and experiments in the laboratory. The present work opens new perspectives for the application of robust nonlinear control theory to complex geosystems, earthquakes and sustainable energy production.

**Index Terms**— Controlling earthquakes, Stability of nonlinear systems, Robust control, Sliding-Mode Control.

## I. INTRODUCTION

Earthquakes are dynamic instabilities that occur in the earth's crust. 65% of the most catastrophic earthquakes happen at depth up to 12 km and are of natural causes [1]. However, earthquakes also occur due to anthropogenic causes. It is nowadays established that injecting fluids in the earth's crust can reactivate seismic faults, leading to important earthquake events (see [2], [3] and [4], to name a few examples).

In this paper, fluid injections are seen from another perspective. Instead of considering them as an earthquake triggering mechanism, they are seen as an input to a dynamical system, which can stabilize it and achieve tracking over a reference signal, if it is adequately designed. This dynamical system is the physical process leading to earthquake instabilities and it is characterized by important nonlinearities due to friction. Moreover, it can present many uncertainties and unmodelled

Diego Gutiérrez-Oribio, Georgios Tzortzopoulos and Ioannis Stefanou are in the École Centrale de Nantes - GeM, 1 Rue de la Noë, 44321 Nantes, France (e-mail: diego.gutiérrez-Oribio@ec-nantes.fr, georgios.tzortzopoulos@ec-nantes.fr, ioannis.stefanou@ec-nantes.fr).

Franck Plestan is in the École Centrale de Nantes - LS2N, UMR CNRS 6004, 1 Rue de la Noë, 44321 Nantes, France (e-mail: franck.plestan@ec-nantes.fr).

dynamics, that a successful control scheme has to compensate. Finally, unlike many existing applications of control theory that target decreasing the response time of the system, here the aim is the opposite, i.e., to slow down the system dynamics. The above mentioned characteristics of this system result in a challenging problem for control theory, with important applications in energy production (e.g., oil, gas, deep geothermal energy, and CO₂ sequestration) and earthquake prevention.

Existing strategies for earthquake control are very limited and ad-hoc. One can refer, for instance, to the field experiments in Rangely, Colorado, US, [5], where seismicity was reduced by turning off the pore pressure. In Dale, New York, US, where earthquakes of magnitude 1 to 1.4 were arrested by dropping the top hole pressure below 5 MPa [6]. More recent field experiments involve the well monitored tests by [7], [8]. However, as mentioned above, all these experiments were based on trial and error and they were not based on control theory.

Recently, an LQR control was designed to stabilize and perform tracking of an earthquake modelled by a MIMO system [9], whereas a double-scale asymptotic approach was employed to design a transfer function-based control in [10]. These first applications of control theory to this problem have shown that earthquakes could be controlled, at least from a mathematical point of view.

The objective of this paper is twofold. First, to evaluate the performance of different controllers and, second, to test them in the laboratory with a specially designed apparatus [11]. The design of the controllers is based on a reduced model for earthquakes. This reduced-order model establishes an average behaviour of a single earthquake fault (see [1], among others). It consists of a single mass that can slide on a rough surface under friction. The frictional interface is usually a complex structure (see [12], [13]), where various physico-chemical phenomena take place during seismic slip (see [1], [14]–[16]). As a result, the friction coefficient and its weakening, not only depend on the slip and the slip-rate, but also on the evolution of the microstructural network, the grain size, the presence and pressure of interstitial fluids, the temperature, time, the reactivation of chemical reactions and other multiphysics couplings (see [17]–[20]).

All these complex phenomena induce the presence of uncertainties and/or perturbations to the plant. As a consequence, they need to be compensated by a robust controller able to obtain a slow-aseismic response. A classic robust approach versus constant perturbations, is the integral action (see [21,
Chapter 12]). Among many robust control approaches, one can also cite the sliding-mode theory [22], [23]. This type of control is known for being insensitive to bounded and matched perturbations, leading to finite-time convergence. Sliding-mode controllers are also known for the simplicity of their gains tuning. The problem with these controllers though is the use of a discontinuous function, sign, that may lead to the so-called *chattering* effect, possibly damaging the actuators. Recently, in order to address the above mentioned drawback, the Continuous Higher-Order Sliding-Modes Algorithms (CHOSMA, see [24]–[26]) have been developed to keep the interesting features of the sliding-mode algorithms (compensation of Lipschitz w.r.t. time perturbations in finite-time), while using a continuous control signal.

In this paper, two control strategies are proposed to achieve a slow-seismic response: the first one is based on CHOSMA and the other based on a Linear Quadratic Regulator extended with integral action (e-LQR based on the original LQR control in [27]). The former is robust against Lipschitz w.r.t. time perturbations, while the latter is robust against Lipschitz w.r.t. state perturbations, presenting both continuous control signal. The designed controllers stabilize and reduce the natural response time of the system, making the energy dissipation to be many orders of magnitude slower compared to a real earthquake event. This is an uncommon paradigm in control theory, where usually the objective is to drive the states of the system to the origin as fast as possible. Finally, to test the feasibility of the presented algorithms, simulations and experimental confirmation using a real laboratory test benchmark, able to reproduce earthquake-like instabilities, are presented. This allows testing the controllers and comparing their performance.

The outline of this work is as follows. The description of the reduced-order model for reproducing earthquake-like instabilities, its instability condition, and the control objectives are given in Section II. The experimental setup is shown in the same section, while the design of the robust control strategies is detailed in Section III. Simulations and experimental results are shown in Section IV, and concluding remarks are made in Section V.

## A. Preliminaries

Throughout the text, the term "Lipschitz w.r.t. the time/state function" is used to call a function that fulfils a Lipschitz condition with respect to the time/state. Let us consider the time-varying differential equation

$$ \dot{x} = f(t, x(t)), \quad t \ge t_0, $$

where $x(t) \in \mathbb{R}^n$ is the state vector; $f: \mathbb{R}_{\ge 0} \times \mathbb{R}^n \to \mathbb{R}^n$ is a function that can be discontinuous, measurable with respect to $t$, and $f(t, 0) = 0$. The initial condition $x(t_0) \in \mathbb{R}^n$ at time instant $t_0 \in \mathbb{R}$ is denoted as $x_0$. The solution of the system is understood in the Filippov's sense (see [28]). Let $\Omega$ be an open subset of $\mathbb{R}^n$, such that $0 \in \Omega$.

**Definition 1:** [21], [29]. The origin, $x = 0$, of the latter system is said to be:

* Locally Stable (LS) if for any $\epsilon > 0$ there is $\delta = \delta(\epsilon, t_0) > 0$

such that if $||x_0|| \le \delta$ then $||x(t)|| \le \epsilon$ for any $x_0 \in \Omega$ and for all $t \ge t_0$.

* Locally Asymptotically Stable (LAS) if it is LS and there is $c = c(t_0) > 0$ such that if $||x_0|| \le c$ then $x(t) \to 0$ for any $x_0 \in \Omega$ and for all $t \ge t_0$.

* Locally Exponentially Stable (LES) if it is LAS and there are $k > 0$, $\lambda > 0$ such that $||x(t)|| \le k ||x_0|| e^{-\lambda(t-t_0)}$ for any $x_0 \in \Omega$ and for all $t \ge t_0$.

* Locally finite-time Stable (LFTS) if it is LAS and $x(t) = 0$ for all $t \ge T(t_0, x_0)$, where $T: \mathbb{R}_{\ge 0} \times \mathbb{R}^n \to \mathbb{R}_{\ge 0}$ is called the settling-time function.

If $\Omega = \mathbb{R}^n$, then $x = 0$ is said to be Globally Stable (GS), Globally Asymptotically Stable (GAS), Globally Exponentially Stable (GES), Globally finite-time Stable (GFTS), respectively.

The definition of weighted homogeneity is introduced to be used in the sequel.

**Definition 2:** [30], [31]. Consider the vector $x \in \mathbb{R}^n$. Its dilation operator is defined as $\Delta_{\epsilon}^r x := (\epsilon^{r_1}x_1, ..., \epsilon^{r_n}x_n)$, $\forall \epsilon > 0$, where $r_i > 0$ are the weights of the coordinates and $r = (r_1, ..., r_n)$ is the vector of weights. A function $V: \mathbb{R}^n \to \mathbb{R}$ (or a vector field $f: \mathbb{R}^n \to \mathbb{R}^n$, or vector-set $F(x) \subset \mathbb{R}^n$) is called $r$-homogeneous of degree $m \in \mathbb{R}$ if the identity $V(\Delta_{\epsilon}^r) = \epsilon^m V(x)$ holds (or $f(\Delta_{\epsilon}^r x) = \epsilon^m \Delta_{\epsilon}^r f(x)$, or $F(\Delta_{\epsilon}^r x) = \epsilon^m \Delta_{\epsilon}^r F(x)$).

The following result is well-known for continuous homogeneous functions (see [26], [32]), and can be extended to semi-continuous functions [33]:

**Lemma 1:** Let $\eta: \mathbb{R}^n \to \mathbb{R}$ and $\gamma: \mathbb{R}^n \to \mathbb{R}$ be two $r$-homogeneous and upper semi-continuous single-valued functions, with the same weights $r = (r_1, ..., r_n)$ and homogeneity degree $m > 0$. Suppose that $\gamma(x) \le 0$ in $\mathbb{R}^n$. If

$$ \{x \in \mathbb{R}^n \setminus \{0\} : \gamma(x) = 0\} \subseteq \{x \in \mathbb{R}^n \setminus \{0\} : \eta(x) < 0\}, $$

then there exists a real number $\lambda^*$ and a constant $c > 0$ so that, for all $\lambda \ge \lambda^*$ and for all $x \in \mathbb{R}^n \setminus \{0\}$ the following inequality is satisfied $\eta(x) + \lambda\gamma(x) \le -c||x|_{r,p}^m$.

Define the function $[\cdot]^{\gamma} := |\cdot|^{\gamma}\text{sign}(\cdot)$, for any $\gamma \in \mathbb{R}_{\ge 0}$ with $\text{sign}(x) = \begin{cases} 1 & x > 0 \\ [-1, 1] & x = 0 \\ -1 & x < 0 \end{cases}$.

The following Lemma is simple (just monotonicity) but useful

**Lemma 2:** [34] Consider the real variables $x, y$, it is always true that $\text{sign}([x+y]^{\beta} - [y]^{\beta}) = \text{sign}(x)$, $\beta > 0$.

## II. PROBLEM STATEMENT

Two scenarios for earthquake modelling used in this paper are presented in this Section. The first one is a reduced model of a real earthquake on which the control designs are performed and numerical simulations are made, whereas the second one, is a novel experimental setup designed to reproduce and then control earthquake-like instabilities in the laboratory.
## A. Reduced Model for Earthquakes

The dynamics of earthquakes can be represented, in average/energetical sense, with the spring-slider analogue system (see [1], [9], [10], [35], [36]) depicted in Fig. 1.

<figure><img src="image_3.png" /><figcaption>Fig. 1. Reduced mechanical model for reproducing earthquake-like instabilities.</figcaption></figure>

This mechanical system consists of a mass, m, which slides on a frictional interface (equivalent to a seismic fault). The mass is connected to a Kelvin-Voigt configuration composed by a spring with stiffness k (equivalent to the apparent elasticity of the host rock) and a dashpot with damping coefficient η (equivalent to the apparent viscosity of the host rock). At the other extremity of the Kelvin-Voigt configuration a constant velocity, v∞, is applied (equivalent to the far-field motion of the tectonic plates). It is assumed Coulomb friction with a friction coefficient μ(δ, δ̇) that depends on the slip δ (block's displacement) and the slip-rate δ̇ (block's velocity). According to Terzaghi's principle, the frictional stress τ takes the following form: τ = μ(δ, δ̇) (σ'_{n} - p), where σ'_{n} is the constant/reference average effective stress (e.g., the overburden due to the weight of the rocks and the interstitial fluid pressure) and p the fluid pressure developed due to injecting fluid. p is the input to the system for which the controllers will be designed and tested.

According to [1] and [35], approximately a rock mass of volume L<sub>ac</sub><sup>3</sup> is mobilized during an earthquake event, where L<sub>ac</sub> is equal to the length of the seismic fault. Therefore, the mobilized mass during an earthquake event is m ≈ ρL<sub>ac</sub><sup>3</sup>, where ρ is the density of the surrounding to a seismic fault rocks. The fault length can be calculated as L<sub>ac</sub> = G/¯k, where G is the shear-modulus of the host rock and ¯k = k/L<sub>ac</sub><sup>2</sup>, its apparent normalized elastic stiffness. The damping coefficient η is given by η = 2ζmω<sub>n</sub>, where ζ is the damping ratio and ω<sub>n</sub> = √k/m, the natural frequency of the reduced system.

Applying the force balance equation, system in Fig. 1 can be represented by the following mathematical model

$$m\ddot{\delta} = -\mu(\delta, \dot{\delta})A(\sigma'_{n} - p) + k(\delta_{\infty} - \delta) + \eta(\nu_{\infty} - \dot{\delta}) + \varphi_{e}(\delta, \dot{\delta}, t), \quad (1)$$

where A ≈ L<sub>ac</sub><sup>2</sup> is the effective contact area (fault rupture area), δ<sub>∞</sub> = v<sub>∞</sub>t the displacement at the extremity of the Kelvin-Voigt configuration, δ̇ the acceleration of the mobilized block, and ϕ<sub>e</sub>(δ, δ̇, t) is a perturbation affecting the system, e.g., an external perturbation or unmodelled dynamics due to the complex frictional phenomena.

In this paper, the friction coefficient μ(δ, δ̇) is assumed to fulfil

$$0 < \mu_{\text{res}} \le \mu(\delta, \dot{\delta}) \le \mu_{\text{max}}, \quad (2)$$

where the constants μ<sub>res</sub> and μ<sub>max</sub> are given. Such assumption is fulfilled by friction laws used in fault mechanics, like the slip-weakening friction law [35], the slip-rate weakening law [37], and the rate-and-state friction law [38], [39]. Notice that the exact frictional rheology is not known in reality, which needs the design of robust controllers, as will be done in Section III.

Additionally, according to [35], the seismic magnitude M<sub>w</sub> is defined as

$$M_w = \frac{2}{3} \log_{10} M_0 - 6.07, \quad M_0 = L_{ac}^3 \Delta\tau, \quad (3)$$

where M<sub>0</sub> is the seismic moment measured in [Nm] and Δτ = (μ<sub>max</sub> - μ<sub>res</sub>)σ'<sub>n</sub>.

## B. Shifted System and Instability Condition

Defining the state variables z<sub>1</sub> = δ and z<sub>2</sub> = δ̇, the state representation of system (1) is

$$\begin{aligned} \dot{z}_1 &= z_2, \\ \dot{z}_2 &= -\mu(z_1, z_2)\hat{N}(\sigma'_n - p) + \hat{k}(\delta_\infty - z_1) + \hat{\eta}(v_\infty - z_2) \\ &\quad + \hat{\varphi}_e(z_1, z_2, t), \end{aligned}$$

where N̂ = A/m, k̂ = k/m, η̂ = η/m and ϕ̂<sub>e</sub>(z<sub>1</sub>, z<sub>2</sub>, t) = ϕ<sub>e</sub>(z<sub>1</sub>, z<sub>2</sub>, t)/m.

The set of equilibrium points (z<sub>1</sub><sup>*</sup>, z<sub>2</sub><sup>*</sup>) of the above system in open loop and without perturbation, i.e. ϕ̂<sub>e</sub>(z<sub>1</sub>, z<sub>2</sub>, t) = 0, is described by

$$z_1^* = -\mu(z_1^*, z_2^*) \frac{\hat{N}}{\hat{k}} \sigma_n' + \delta_\infty + \frac{\hat{\eta}}{\hat{k}} v_\infty, \quad z_2^* = 0.$$

Note that the equilibrium (z<sub>1</sub><sup>*</sup>, z<sub>2</sub><sup>*</sup>) depends on the friction coefficient μ(z<sub>1</sub><sup>*</sup>, z<sub>2</sub><sup>*</sup>). In this paper, the controller design will be made when the system reaches the above equilibria, which can be unstable. The maximum value of the friction coefficient, μ(0, 0) = μ<sub>max</sub>, is considered for being on the verge of slip and the system is shifted as follows. Setting z<sub>1</sub><sup>*</sup> = 0 and the new state variables x<sub>1</sub> = z<sub>1</sub> - z<sub>1</sub><sup>*</sup> and x<sub>2</sub> = z<sub>2</sub> - z<sub>2</sub><sup>*</sup>, the shifted system reads as

$$\begin{aligned} \dot{x}_1 &= x_2, \\ \dot{x}_2 &= -[\mu(x_1, x_2) - \mu^*] \hat{N} \sigma'_n + \mu(x_1, x_2) \hat{N} p - \hat{k} x_1 - \hat{\eta} x_2 \\ &\quad + \hat{\varphi}_e(x_1, x_2, t), \end{aligned} \quad (4)$$

where μ<sup>*</sup> = μ(0, 0) = μ<sub>max</sub>. Note that if ϕ̂<sub>e</sub>(x<sub>1</sub>, x<sub>2</sub>, t) = 0, system (4) has an equilibrium point located at the origin x<sub>1</sub><sup>*</sup> = x<sub>2</sub><sup>*</sup> = 0 in open loop.

In order to analyse the stability of the origin of system (4) without the perturbation term ϕ̂<sub>e</sub>(x<sub>1</sub>, x<sub>2</sub>, t), consider its Jacobian matrix J(x<sub>1</sub>, x<sub>2</sub>) evaluated at the origin as

$$J(0, 0) = \begin{bmatrix} 0 & 1 \\ -\hat{k} - \hat{N} \sigma'_{n} \frac{\partial \mu}{\partial x_1} |_{(x_1, x_2)=(0,0)} & -\hat{\eta} - \frac{\partial \mu}{\partial x_2} |_{(x_1, x_2)=(0,0)} \\ -\hat{k} - \hat{N} \sigma'_{n} \frac{\partial \mu}{\partial x_1} |_{(x_1, x_2)=(0,0)} & 0 \end{bmatrix},$$

where the conditions to have an unstable origin are

$$k < -A\sigma'_{n} \frac{\partial \mu}{\partial x_1} |_{(x_1, x_2)=(0,0)}, \quad \eta < -m \frac{\partial \mu}{\partial x_2} |_{(x_1, x_2)=(0,0)}. \quad (5)$$
The above inequalities are in accordance with the nominal studies of [1], [9], [40] and they show that dynamic instability will take place when the elastic unloading of the springs or the apparent viscosity of the host rock cannot be counterbalanced by friction.

### C. Experimental setup

A novel experimental apparatus for reproducing and controlling earthquake-like instabilities in the laboratory, has been designed (see [11], [41], [42]). This experimental setup is depicted in Fig. 2.

Two loading systems are used in order to apply shear (displacement controlled by the vertical loading system) and normal (pressure controlled by the horizontal loading system) stresses to the sheared interfaces. The horizontal loading system consists of an inflatable rubber bag whose pressure can be controlled, in real-time, through a fast response electro-pneumatic pressure regulator. This system can simulate fluid injection/extraction into/from the fault interface by properly adjusting the (air) pressure in the rubber bag to desired levels resulting in variations of the effective normal stress in the sheared interfaces. The vertical loading system consists of a press compressing slowly a linear-elastic spring, which simulates the stored energy of the earthquake.

The specimen consists of three 3D printed samples of sand particles (see [43] for more details about this surrogate material and its characterization) and it is located below the spring. The middle block of the specimen simulates the mobilized mass of the rocks and its frictional interfaces with the adjacent blocks the seismic fault. No specific elements are used for damping, which is provided naturally by the various components of the experimental setup and the specimen itself. The friction coefficient between the blocks can be measured in the experimental setup (see [41], [43]). However, these tests are not presented here because the objective is to design robust controllers that are agnostic to the frictional rheology. The frictional rheology in a real fault might depend on the slip, slip-rate and time, but also in thermal, chemical and other processes. Moreover, due to heterogeneity, the frictional properties of faults involve uncertainties and contrary to a laboratory experiment, they cannot be inferred with high accuracy. Therefore, successful controllers must be able to compensate all these sources of uncertainties and to be able to obtain slow aseismic response with minimum information about the system.

For measuring the applied vertical load and the horizontal forces, load cells are used. For measuring the average slip of the middle block, two vertically placed Linear Variable Differential Transformers (LVDTs) are used. Their readings are averaged to eliminate parasitic measurements related to possible rotation of the middle block. Finally, the pump for controlling the pressure in the rubber bag simulates the fluid injection to the fault and can supply pressures up to 1 [MPa]. Note that in this setup, increasing the pressure in the rubber bag corresponds to decreasing the pressure, *p*, in the fault (see also eq. (1)). However in the following, the distinction between pressure in the bag and pressure in the fault is not

made and the results are presented in function of the pressure at the fault (i.e., $p = -p_{bag}$). All these sensors are connected to a data acquisition device processed using the LabVIEW software, which processes the data with a sampling rate of 1 [ms]. LabVIEW also allows the implementation of our robust controllers for controlling earthquake-like instabilities. See [41, Chapter 5] for more details about the experimental setup.

**Remark 1:** The electro-pneumatic pressure regulator has a response rate of 1 [ms]. It is controlled by a PID algorithm, which is faster than the response rate of the pump. The operating system of the computer in which the controllers were implemented is of 1 to 2 [ms]. The above characteristics times are much lower than the characteristic time of the system instability, which is 50 [ms]. Therefore, the actuator dynamics (the pump) was not taken into account in the system (1).

### D. Control Objective

As shown in system (4) and Fig. 1, the fluid pressure *p* is the only input acting on the dynamics of the mechanical system. In a real-scale scenario, fluid injections in the earth's crust change the fluid pore pressure over seismic faults [8]. As shown in [36] among others, this can destabilize the fault system and induce/trigger larger earthquakes. In order to illustrate this phenomenon, two numerical simulations of the dynamical system (4) are shown in Fig. 3. The first one has no input pressure (Natural earthquake using $p = 0$ [MPa]) and the second one has a constant pressure (Induced earthquake using $p = 5$ [MPa]). Both simulations were made starting from the origin and with an external perturbation $\hat{\phi}_e = 3.2 \times 10^{-15}$ [m/s²] applied at $t = 0$ [s] just to move the states out of the equilibrium point. Note how an earthquake event ($\delta$ evolving with high rate) has been triggered in both cases, but a larger and faster one was obtained on the second case where $p > 0$. The system parameters used for these simulations are shown in Table I (Real-fault column).

In the previous paragraph, the fact that open-loop fluid injections in the earth's crust might risk to stimulate large seismic events is highlighted. To prevent this, one could adequately adjust the fluid pressure (input *p*) by employing control techniques in order to stabilize the system (4) origin and/or track a reference input signal, releasing the stored elastic energy smoothly ($\delta$ evolving slowly) and not abruptly as shown above.

From a mathematical point of view, a reduced-order model for earthquakes like the system (1) captures well the low-frequency instabilities, which are of most importance, because they carry most of the seismic energy. For the experimental confirmation, the designed apparatus in Fig. 2 reproduces these low-frequency instabilities and successfully designed controllers must be able to dissipate them slowly. See [42], [44], for theoretical works for the control of more detailed models of the physical process of rupture of a single mature seismic fault.

Therefore, the objective in the Section III is to design a control law *p* driving $\delta$ and $\dot{\delta}$ in system (4) to follow some desired predefined references of slow slip rate, $r(t)$, $\dot{r}(t)$,
<figure><img src="image_4.png" /><figcaption>Fig. 2. Experimental apparatus for reproducing and controlling earthquake-like instabilities: (a)Schematic figure, (b)Real configuration</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Slip and slip-rate in terms of time for two different scenarios. Blue curve: Natural earthquake. Red Dashed Curve: Induced earthquake.</figcaption></figure>

resulting in a slow-aseismic response, despite the presence of $\hat{\phi}_e$ and parametric uncertainties. Furthermore, this control design will be tested in Section IV on simulations and on the experimental setup depicted in Fig. 2.

The designed controllers will reduce the natural response time of the system slowing its energy dissipation and eliminating bursts of kinetic energy (earthquake phenomenon). This is an uncommon paradigm in control theory where usually the objective is to drive the states of the system to the origin as fast as possible

### III. CONTROL DESIGN

Two control approaches are presented in this section: the first one is based on continuous sliding-mode control theory [25], [26] whereas the second one is based on LQR control

[27]. Both these controllers must be able to force the states of system (4) to track a predefined reference, using a continuous control signal, in spite of uncertainties/perturbations of the system (4).

It is worth emphasizing that due to the nature of the physical process of the earthquake phenomenon, the controllers have to be robust to uncertainties regarding the material properties and plant dynamics, as well as to disturbances of the input signals and measurements. This motivates the choice of sliding-modes control, which will be compared with a robust LQR control in order to slow down the system and make it follow a desired reference trajectory.

The desired reference for the output $y = x_1$ is a smooth function reading as

$$r(t) = d_{max}s^3(10 - 15s + 6s^2), \quad (6)$$

where $s = t/t_{op}$, $d_{max}$ the target displacement and $t_{op}$ the operational time of the tracking strategy. The constant $d_{max}$ is the distance the fault slides dynamically in order to reach its sequent stable equilibrium point. Notice that the parameter $t_{op}$ is free to be decided depending on the earthquake control scenario that one wants to apply. Nevertheless, $t_{op}$ has to be sufficiently high with respect to the characteristic time of the earthquake phenomenon, but low enough to achieve aseismic slip with higher velocity than the far-field velocity ($v_\infty$ in (1)), for the control scenario to make sense.

The choice of the reference output $y = x_1$ is motivated by the need to control the average slip over the fault. This average slip is directly connected with the magnitude of an earthquake through the seismic moment [35]. Therefore, by controlling the rate of the average slip, the system is forced to release its energy in a quasi-static way, i.e., aseismically. See [9], [42] for more details.

#### A. Sliding Mode based Control

To perform the tracking of the desired references $r(t)$, $\dot{r}(t)$, a sliding-mode-based control is designed. Defining the track-
ing error variables,

$$e_1 = x_1 - r, \quad e_2 = x_2 - \dot{r}, \qquad (7)$$

the error dynamics reads as

$$\begin{aligned} \dot{e}_1 &= e_2, \\ \dot{e}_2 &= -[\mu(e_1 + r, e_2 + \dot{r}) - \mu^*] \hat{N} \sigma'_n + \mu(e_1 + r) \hat{N} p \\ &\quad - \hat{k}(e_1 + r) - \hat{\eta}(e_2 + \dot{r}) + \hat{\varphi}_e(e_1 + r, e_2 + \dot{r}, t) - \ddot{r}. \end{aligned} \qquad (8)$$

If the exact knowledge of the system parameters and the system dynamics would be available, all the known dynamics in $\dot{e}_2$ can be compensated, in order to get the nominal error system,

$$\dot{e}_1 = e_2, \quad \dot{e}_2 = \nu, \qquad (9)$$

by designing the control p as

$$p = \frac{1}{\mu(e_1 + r, e_2 + \dot{r}) \hat{N}} \left\{ \nu + [\mu(e_1 + r) - \mu^*] \hat{N} \sigma'_n + \hat{k}(e_1 + r) + \hat{\eta}(e_2 + \dot{r}) - \hat{\varphi}_e(e_1 + r, e_2 + \dot{r}, t) + \ddot{r} \right\}, \qquad (10)$$

and with the new control input, $\nu$, designed to force $e_1, e_2$ towards zero. A solution could be the linear feedback control $\nu = -k_1 e_1 - k_2 e_2$ with any $k_1, k_2 > 0$. However, system (9) is valid only in the nominal case. If this is not the scenario, the application of the state-feedback (10) with uncertain parameters will not lead to (9). In this case, consider the feedback control

$$p = \frac{1}{\mu_0 \hat{N}_0} \nu, \qquad (11)$$

where the sub index ‘0’ indicates the nominal value of the real considered parameter. Notice that some additional nominal parameters could be used in (11), as $\hat{k}_0, \hat{\eta}_0$ or even the known term $\ddot{r}(t)$. However, the objective here is to design a controller requiring a limited amount of information.

Therefore, the closed-loop system obtained from (8) and (11) reads as

$$\dot{e}_1 = e_2, \quad \dot{e}_2 = \beta(t, e) [\nu + h(t, e)], \qquad (12)$$

where $e = [e_1, e_2]^T$, $\beta(t, e)$, is the uncertain control coefficient, and $h(t, e)$ is a matched perturbation affecting the system. These terms read as

$$\begin{aligned} \beta(t, e) &= \frac{\mu(e_1 + r, e_2 + \dot{r}) \hat{N}}{\mu_0 \hat{N}_0}, \\ h(t, e) &= \frac{1}{\beta(t, e)} \left\{ -[\mu(e_1 + r, e_2 + \dot{r}) - \mu^*] \hat{N} \sigma'_n - \hat{k}(e_1 + r) - \hat{\eta}(e_2 + \dot{r}) + \hat{\varphi}_e(e_1 + r, e_2 + \dot{r}, t) - \ddot{r} \right\}. \end{aligned} \qquad (13)$$

Both these terms are assumed to fulfil in the operating domain

$$0 < b_m \le \beta(t, e) \le b_M, \quad \left| \frac{dh(t, e)}{dt} \right| \le \bar{L}, \qquad (14)$$

with known constants $b_m, b_M, \bar{L}$.

**Remark 2:** The condition for $\beta(t, e)$ in (14) is satisfied because of the assumption of $\mu(e_1 + r, e_2 + \dot{r}) = \mu(x_1, x_2)$ in (2). The condition for $h(t, e)$ in (14) is satisfied (locally inside of a domain) because of the definition of $r(t)$ in (6) and if the external perturbation term $\hat{\varphi}_e(e_1 + r, e_2 + \dot{r}, t)$ is Lipschitz w.r.t. time. As a result, the tracking result obtained in the sequel is valid locally.

The design of the control input $\nu$ able to stabilize (12) at $e_1 = e_2 = 0$, despite the presence of $\beta(t, e), h(t, e)$, results in an aseismic motion of system (4). For this purpose, consider the Second-Order Continuous Twisting Algorithm (2-CTA) introduced in [25]

$$\begin{aligned} \nu &= -\lambda^{\frac{2}{3}} k_1 [e_1]^{\frac{1}{3}} - \lambda^{\frac{1}{2}} k_2 [e_2]^{\frac{1}{2}} + \xi_1, \\ \dot{\xi}_1 &= -\lambda k_3 [e_1]^0 - \lambda k_4 [e_2]^0, \end{aligned} \qquad (15)$$

and the Second-Order Discontinuous Integral Algorithm (2-DIA) introduced in [26]

$$\begin{aligned} \nu &= -\lambda^{\frac{1}{2}} k_{I2} \left[ [e_2]^{\frac{3}{2}} + \lambda^{\frac{1}{2}} k_{I1}^{\frac{3}{2}} e_1 \right]^{\frac{1}{3}} + \xi_1, \\ \dot{\xi}_1 &= -\lambda k_{I3} \left[ e_1 + \lambda^{-\frac{1}{2}} k_{I4} [e_2]^{\frac{3}{2}} \right]^{0}. \end{aligned} \qquad (16)$$

Both of these algorithms consist of a static homogeneous finite-time controller and a discontinuous integral action, aimed at estimating and compensating the uncertainties and perturbations. Notice that the presence of the discontinuous function, $[\cdot]^0$, in the integral action finally results in a continuous control signal.

**Theorem 1:** The origin of the system (12) is LFTS, despite the presence of the Lipschitz w.r.t. the time uncertainties/perturbations $h(t, e)$ and bounded uncertain coefficient $\beta(t, e)$ satisfying (14), if the control $\nu$ takes the form of (15) or (16), with gains appropriately chosen.

As a consequence of Theorem 1, the state variables $\delta, \dot{\delta}$ of system (1), are locally driven in finite-time to the desired references $r(t), \dot{r}(t)$ defined in (6).

**Proof:** The closed-loop system (12), with controller (15) reads as

$$\begin{aligned} \dot{e}_1 &= e_2, \\ \dot{e}_2 &= \beta(t, e) \left( -\lambda^{\frac{2}{3}} k_1 [e_1]^{\frac{1}{3}} - \lambda^{\frac{1}{2}} k_2 [e_2]^{\frac{1}{2}} + e_3 \right), \quad (17) \\ \dot{e}_3 &= -\lambda k_3 [e_1]^0 - \lambda k_4 [e_2]^0 + \dot{h}(t, e), \end{aligned}$$
with $e_3 = \xi_1 + h(t, e)$. If the controller (16) is used, one gets

$$\begin{aligned} \dot{e}_1 &= e_2, \\ \dot{e}_2 &= \beta(t, e) \left( -\lambda^{\frac{1}{2}} k_{I2} \left[ [e_2]^{\frac{3}{2}} + \lambda^{\frac{1}{3}} k_{I1}^{\frac{3}{2}} e_1 \right]^{\frac{1}{3}} + e_3 \right), \quad (18) \\ \dot{e}_3 &= -\lambda k_{I3} \left[ e_1 + \lambda^{-\frac{1}{2}} k_{I4} [e_2]^{\frac{3}{2}} \right]^{0} + \dot{h}(t, e). \end{aligned}$$

The solutions of both systems are understood in the Filippov's sense (see [28]). It will be shown that given bound parameters $b_m, b_M$ and a fixed $\bar{L}$, it is possible to find the values of both control gains such that $e = 0$ is LFTS, but for a not assignable value $L^*$ of the Lipschitz bound $\bar{L}$. To meet this value, one just has to scale the previously obtained gains using $\lambda$ such that $\lambda L^* \ge \bar{L}$.
The scaling with $\lambda > 0$ does not alter the stability of the system. This can be shown by first performing the linear change of variables $z = \lambda e$ in systems (17) and (18). The system in the variables $z$ has exactly the same form as (17) and (18), respectively, but with the gains scaled as in (15), i.e., $(k_1, k_2, k_3, k_4) \rightarrow (k_1\lambda^{\frac{2}{3}}, k_2\lambda^{\frac{1}{2}}, k_3\lambda, k_4\lambda)$ and as in (16), i.e., $(k_{I1}, k_{I2}, k_{I3}, k_{I4}) \rightarrow (k_{I1}\lambda^{\frac{1}{3}}, k_{I2}\lambda^{\frac{1}{2}}, k_{I3}\lambda, k_{I4}\lambda^{-\frac{1}{2}})$ (see [25], [45] for more details on each scaling). Therefore, both systems are equivalent in terms of stability.

The details of the proof for the closed-loop dynamics of system (18) with the controller (16), which follows closely the idea in [45], [46], is presented. The proof for the controller (15) is similar, and a brief comment below is provided. First, the gain $\lambda = 1$ is fixed. Then, a homogeneous and continuously differentiable Lyapunov function candidate is considered as

$$V(\phi) = \frac{3}{5}\gamma_1 |\phi_1|^{\frac{5}{3}} + \phi_1\phi_2 + \frac{2}{5}k_{I1}^{-\frac{3}{2}}|\phi_2|^{\frac{5}{2}} + \frac{1}{5}|\phi_3|^5,$$

where

$$\phi_1 = e_1 - [\phi_3]^3, \quad \phi_2 = e_2, \quad \phi_3 = k_{I1}^{-\frac{1}{2}}k_{I2}^{-1}e_3.$$

$V(\phi)$ is positive definite if $k_{I1} > 0$ and $\gamma_1 > 0$ is selected sufficiently large. Its derivative along the trajectories of the system (18) is

$$\dot{V} = -k_{I1}^{\frac{1}{2}}k_{I2}\beta(t, e) \left( \phi_1 + k_{I1}^{-\frac{3}{2}}[\phi_2]^{\frac{3}{2}} \right) \times \left[ \left[ \left( \phi_1 + k_{I1}^{-\frac{3}{2}}[\phi_2]^{\frac{3}{2}} \right) + [\phi_3]^3 \right]^{\frac{1}{3}} - \phi_3 \right] + \left( \gamma_1[\phi_1]^{\frac{2}{3}} + \phi_2 \right) \phi_2 - k_{I3}k_{I1}^{-\frac{1}{2}}k_{I2}^{-1}[\phi_3]^4 - 3|\phi_3|^2 \left( \gamma_1[\phi_1]^{\frac{2}{3}} + \phi_2 \right) \right] \left\{ \left[ \phi_1 + [\phi_3]^3 + k_{I4}[\phi_2]^{\frac{3}{2}} \right]^{0} - \frac{1}{k_{I3}}\dot{h}(t, e) \right\}.$$

Using the bounds for $\beta(t, e)$ and $h(t, e)$ in (14), the derivative reads as

$$\dot{V} \leq -\tilde{k}_{I2}b_m \left( \phi_1 + k_{I1}^{-\frac{3}{2}}[\phi_2]^{\frac{3}{2}} \right) \times \left[ \left[ \left( \phi_1 + k_{I1}^{-\frac{3}{2}}[\phi_2]^{\frac{3}{2}} \right) + [\phi_3]^3 \right]^{\frac{1}{3}} - \phi_3 \right] + \left( \gamma_1[\phi_1]^{\frac{2}{3}} + \phi_2 \right) \phi_2 - \tilde{k}_{I3}[\phi_3]^4 - 3|\phi_3|^2 \left( \gamma_1[\phi_1]^{\frac{2}{3}} + \phi_2 \right) \right] \left\{ \left[ \phi_1 + [\phi_3]^3 + k_{I4}[\phi_2]^{\frac{3}{2}} \right]^{0} - [-\tilde{L}^*, \tilde{L}^*] \right\},$$

where $\tilde{L}^* = \frac{L^*}{k_{I3}}$, $\tilde{k}_{I2} = k_{I1}^{\frac{1}{2}}k_{I2}$ and $\tilde{k}_{I3} = k_{I3}k_{I1}^{-\frac{1}{2}}k_{I2}^{-1}$. Note that the value of $L^*$ is not given, but it has to be found. Using Lemma 2, the first term can be proven to be negative semi-definite and it vanishes only on the set $S_1 = \{k_{I1}^{\frac{3}{2}}\phi_1 + [\phi_2]^{\frac{3}{2}} = 0\}$.

Evaluating $\dot{V}$ on this set results in

$$\dot{V}|_{S_1} \leq - (\gamma_1 k_{I1}^{-1} - 1) |\phi_2|^2 - \tilde{k}_{I3} |\phi_3|^2 \left[ 3 (\gamma_1 k_{I1}^{-1} - 1) \phi_2 + [\phi_3]^2 \right] \times \left\{ \left[ -k_{I1}^{-\frac{3}{2}} + k_{I4} \right] [\phi_2]^{\frac{3}{2}} + [\phi_3]^3 \right\}^0 - [-\tilde{L}^*, \tilde{L}^*]$$

If $\gamma_1 > k_{I1}$, the first term is negative semidefinite, and it is zero only on the set $S_2 = \{\phi_2 = 0\}$. Evaluating $\dot{V}|_{S_1}$ on $S_2$ reads as

$$\dot{V}|_{S_1 \cap S_2} \leq -\tilde{k}_{I3} [\phi_3]^4 \{[\phi_3]^0 - [-\tilde{L}^*, \tilde{L}^*]\}.$$

This is negative if $\tilde{L}^* = \frac{L^*}{k_{I3}} < 1$. By Lemma 1, $\dot{V}|_{S_1} < 0$ is proven by selecting $\tilde{k}_{I3} > 0$ small. Using Lemma 1, again, it is possible to make $\dot{V} < 0$ selecting $\tilde{k}_{I2} > 0$ sufficiently large.

For the closed-loop system (17), the smooth and homogeneous Lyapunov function

$$V(e) = \alpha_1 |e_1|^{\frac{5}{3}} + \alpha_2 e_1 e_2 + \alpha_3 |e_2|^{\frac{5}{2}} + \alpha_4 e_1 [e_3]^2 - \alpha_5 e_2 e_3^3 + \alpha_6 |e_3|^5,$$

can be selected. It has been shown in [25] that the origin of system (17) is GFTS for appropriate selected gains $k_1 > 0$, $k_2 > 0$, $k_3 > 0$, $k_4$. Although this result has been obtained in [25] using an SOS algorithm and assuming that the control coefficient is known and constant, the proof can be also extended to the actual case using similar arguments as those used above for the controller (16).

Therefore, the origin of systems (17) and (18) is asymptotically stable, but locally, as a consequence of Remark 2. Moreover, since the systems are homogeneous of negative degree, such origins are LFTS (see [47]).

## B. LQR-based Control

Following the sliding-mode based control design, this subsection presents an extended Linear Quadratic Regulator (e-LQR). The term extended is due to the integral action added to a standard LQR algorithm, in such a way that Lipschitz w.r.t. the states uncertainties are compensated with the resultant control.

Starting from (4), the plant is extended with a double integrator to improve the tracking of the reference trajectory. The stability of this augmented plant controlled by a LQR controller is being proved by a Lyapunov's approach. For this purpose, define

$$\dot{\xi}_1 = x_1 - r, \quad \dot{\xi}_2 = \xi_1. \qquad (19)$$
The matrix form of system (4) is

$$
\begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \end{bmatrix}_{\dot{x}} = \underbrace{\begin{bmatrix} 0 & 1 \\ -\hat{k} & -\hat{\eta} \end{bmatrix}}_{A(t)} \underbrace{\begin{bmatrix} x_1 \\ x_2 \end{bmatrix}}_{x} + \underbrace{\begin{bmatrix} 0 \\ \mu(x_1, x_2) \hat{N} \end{bmatrix}}_{B(t,x)} p + \underbrace{\begin{bmatrix} 0 \\ -[\mu(x_1, x_2) - \mu^*] \hat{N} \sigma'_n + \hat{\phi}_e(x_1, x_2, t) \end{bmatrix}}_{g(t,x)}, \quad (20)
$$

whereas the matrix form of $(\xi_1, \xi_2)$-system is

$$
\underbrace{\begin{bmatrix} \dot{\xi}_1 \\ \dot{\xi}_2 \end{bmatrix}}_{\dot{\xi}} = \underbrace{\begin{bmatrix} 0 & 0 \\ 1 & 0 \end{bmatrix}}_{C_{\xi}} \underbrace{\begin{bmatrix} \xi_1 \\ \xi_2 \end{bmatrix}}_{\xi} + \underbrace{\begin{bmatrix} 1 & 0 \\ 0 & 0 \end{bmatrix}}_{C_x} \underbrace{\begin{bmatrix} x_1 \\ x_2 \end{bmatrix}}_{x} + \underbrace{\begin{bmatrix} -r(t) \\ 0 \end{bmatrix}}_{r_{\xi}(t)}. \quad (21)
$$

Following a conventional integral control design (see [21, Chapter 12] for example) and supposing a constant reference $r(t) = r_0$, the augmented system composed by (20)-(21) reads as

$$
\underbrace{\begin{bmatrix} \dot{x} \\ \dot{\xi} \end{bmatrix}}_{\dot{x}_a} = \underbrace{\begin{bmatrix} A(t) & 0_{2 \times 2} \\ C_x & C_\xi \end{bmatrix}}_{A_a(t)} \underbrace{\begin{bmatrix} x \\ \xi \end{bmatrix}}_{x_a} + \underbrace{\begin{bmatrix} B(t, x) \\ 0_{2 \times 1} \end{bmatrix}}_{B_a(t, x_a)} p + \underbrace{\begin{bmatrix} g(t, x) \\ 0_{2 \times 1} \end{bmatrix}}_{g_a(t, x_a)}. \quad (22)
$$

Consider additive (matched) uncertainties to $A_a(t)$ and $B_a(t, x_a)$ such that,

$$A_a(t) = A_0 + \Delta A(t), \quad (23)$$

$$B_a(t, x_a) = B_0 + \Delta B(t, x_a), \quad (24)$$

where

$$A_0 = \begin{bmatrix} 0 & 1 & 0 & 0 \\ -\hat{k}_0 & -\hat{\eta}_0 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix}, \Delta B(t, x_a) = \begin{bmatrix} 0 \\ \Delta[\mu(x_1, x_2)\hat{N}] \\ 0 \\ 0 \end{bmatrix}$$

$$\Delta A(t) = \begin{bmatrix} 0 & 0 & 0 & 0 \\ -\Delta \hat{k} & -\Delta \hat{\eta} & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix}, B_0 = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}.$$

The sub-index '0' represents the nominal value, whereas the quantity with the prefix 'Δ' corresponds to the uncertainties of the respective variable. In addition, the assumption (2) was used for the friction coefficient $\mu(x_1, x_2)$, so the variation $\Delta B(t, x_a)$ is always positive semi-definite.

From (22), (23), and (24), one gets

$$\dot{x}_a = A_0 x_a + \Delta B(t, x_a) p + B_0 [p + h(t, x_a)], \quad (25)$$

where $h(t, x_a) = B_0^+ \Delta A(t) x_a + B_0^+ g(t, x_a)$, with $B_0^+$ is the Moore-Penrose inverse matrix of $B_0$, i.e.,

$$B_0^+ = \begin{bmatrix} 0 & \frac{1}{\mu_{res}\hat{N}} & 0 & 0 \end{bmatrix}. \quad (26)$$

Then, the nonlinear vector $h(t, x_a)$ can be written as

$$
\begin{align}
h(t, x_a) &= B_0^+ [\Delta A(t) x_a + g(t, x_a)] \nonumber \\
&= -\frac{1}{\mu_{res} \hat{N}} \left\{ \Delta \hat{\eta} x_2 + \Delta \hat{k} x_1 \right. \tag{27} \\
&\qquad \left. + [\mu(x_1, x_2) - \mu^*] \sigma'_n \hat{N} - \hat{\phi}_e(x_1, x_2, t) \right\}. \nonumber
\end{align}
$$

Assuming the external perturbation $\hat{\phi}_e(x_1, x_2, t)$ to be Lipschitz w.r.t. the states, i.e., $|\hat{\phi}_e(x_1, x_2, t)| \le |\hat{\phi}_{1e}x_1 + \hat{\phi}_{2e}x_2|$ for some known positive constants $\hat{\phi}_{1e}, \hat{\phi}_{2e}$, the norm of $h(t, x_a)$ reads as

$$
\begin{align*}
\|h(t, x_a)\| &= \left\| \frac{1}{\mu_{res} \hat{N}} \right. \\
&\qquad \left. \times \left[ \Delta \hat{k} + \frac{\mu(x_1, x_2) - \mu^*}{x_1} \sigma'_n \hat{N} + \hat{\phi}_{1e} \quad \Delta \hat{\eta} + \hat{\phi}_{2e} \quad 0 \quad 0 \right] x_a \right\| \\
&\le \left\| \frac{1}{\mu_{res} \hat{N}} \right. \\
&\qquad \left. \times \left[ \Delta \hat{k}_{\max} + \mu_{x_1}^{\max} \sigma'_n \Delta \hat{N}_{\max} + \hat{\phi}_{1e} \quad \Delta \hat{\eta}_{\max} + \hat{\phi}_{2e} \quad 0 \quad 0 \right] x_a \right\|,
\end{align*}
$$

where the bound $\mu_{x_1}^{\max} = \left| \frac{\partial \mu(x_1, x_2)}{\partial x_1} \right|_{\max}$ corresponds to the maximum absolute softening slope of the friction. The subscript 'max' denotes the maximum variation in absolute term from the respective nominal values.

Finally, one gets the bounds of the variation coefficient $\Delta B(t, x_a)$ and perturbation $h(t, x_a)$ in system (25) as,

$$0 \le \Delta B(t, x_a), \quad \|h(t, x_a)\| \le \|Gx_a\|, \quad (28)$$

with $G$ defined as,

$$G = \frac{1}{[\mu(x_1)\hat{N}]_{\min}} \times \left[ \Delta \hat{k}_{\max} + \mu_{\hat{x}_1}^{\max} \sigma'_n \Delta \hat{N}_{\max} + \hat{\phi}_{1e} \quad \Delta \hat{\eta}_{\max} + \hat{\phi}_{2e} \quad 0 \quad 0 \right]. \quad (29)$$

Inspired from the original LQR control proposed in [27], define the e-LQR control input $p$ designed for the augmented system (25) as

$$p = -R^{-1} B_0^T \Theta x_a = -[k_1 \quad k_2 \quad k_3 \quad k_4] [x_1 \quad x_2 \quad \xi_1 \quad \xi_2]^T, \quad (30)$$

where $R$ is a positive definite matrix to be chosen and $\Theta$ the positive-definite solution of the following Continuous Algebraic Riccati Equation (CARE),

$$A_0^T \Theta + \Theta A_0 - \Theta B_0 R^{-1} B_0^T \Theta = -Q, \quad (31)$$

with

$$Q = Q_0 + \|R^{1/2}\|^2 G^T G, \quad (32)$$

and $Q_0$ a positive definite matrix to be chosen.

The next Theorem states the main result of the e-LQR control.

**Theorem 2:** The origin of the augmented closed-loop system (25) is GES, with the presence of the Lipschitz w.r.t. the states perturbation $h(t, x_a)$ and positive semi-definite term $\Delta B(t, x_a)$ satisfying (28) and (29), when the control input $p$ takes the form of (30)-(32).

As a consequence of Theorem 2, the slip $\delta$ and slip-rate fault $\dot{\delta}$ of the original system (1) are driven globally and exponentially to a desired constant reference $r(t) = r_0$.

*Proof:* Select as Lyapunov candidate the positive definite and radially unbounded function $V(x_a) = x_a^T \Theta x_a$, where
$\Theta = \Theta^T > 0_{4 \times 4}$ as a positive definite matrix. The following inequalities are true for the Lyapunov candidate

$$ \lambda_{\min}(\Theta) \|x_a\|^2 \leq V(x_a) \leq \lambda_{\max}(\Theta) \|x_a\|^2, \quad (33) $$

$$ \left\| \frac{\partial V(x_a)}{\partial x_a} \right\| \leq 2\lambda_{\max}(\Theta) \|x_a\|, \quad (34) $$

where $\lambda_{\min}(\Theta)$ and $\lambda_{\max}(\Theta)$ are the minimum and maximum eigenvalues of the matrix $\Theta$, respectively.

The time derivative of the Lyapunov candidate $V(x_a)$ along the trajectories of closed loop system (25) and (30) reads as,

$$ \dot{V} = \dot{x}_a^T \Theta x_a + x_a^T \Theta \dot{x}_a = x_a^T (A_0^T \Theta + \Theta A_0 - \Theta B_0 R^{-1} B_0^T \Theta) x_a \quad (35a) $$

$$ -2x_a^T \Theta \frac{B_0 R^{-1} \Delta B^T + (B_0 R^{-1} \Delta B^T)^T}{2} \Theta x_a \quad (35b) $$

$$ -x_a^T \Theta B_0 R^{-1} B_0^T \Theta x_a + 2x_a^T \Theta B_0 h, \quad (35c) $$

where the dependencies of $h(t, x_a)$, $\Delta B(t, x_a)$ have been omitted for simpler notation.

The term (35a) is the CARE as defined in (31). Moreover, due to the fact that the nominal matrix $B_0$ has been selected to have the term $\Delta B(t, x_a)$ non-negative, the non-Hermitian matrix $B_0 R^{-1} \Delta B^T$ is always positive semi-definite and the term on (35b) is a non-positive scalar.

Therefore, letting $z = R^{-1/2} B_0^T \Theta x_a$, one gets

$$ \dot{V} \leq -x_a^T Q x_a - z^T z + 2z^T R^{1/2} h, \quad (36) $$

$$ \text{with } R^{-1} = R^{-1/2} R^{-1/2}. $$

Using the definition of $Q$ in (32) and the Lipschitz condition for $h$ in (28), it comes

$$ \begin{aligned} \dot{V} &\leq -x_a^T Q_0 x_a - \|R^{1/2}\|^2 x_a^T G^T G x_a - z^T z + 2z^T R^{1/2} h, \\ &\leq -x_a^T Q_0 x_a - \|R^{1/2}\|^2 \|G x_a\|^2 - \|z\|^2 \\ &\quad + 2 \|R^{1/2}\| \|z\| \|G x_a\|, \\ &\leq -x_a^T Q_0 x_a - (\|R^{1/2}\| \|G x_a\| - \|z\|)^2, \\ &\leq -x_a^T Q_0 x_a \leq -\lambda_{\min}(Q_0) \|x_a\|^2 < 0. \end{aligned} \quad (37) $$

As a conclusion, $V(x_a)$ is a Lyapunov function for system (25). Therefore, its origin is GES.

### C. Control Strategies Comparison

The designed control strategies, $p$, shown in (11), (15), (16) and (30)-(32) represent the fluid injected to the fault able to achieve aseismic response in the model (1) by tracking a slow reference, robustly and by using a continuous control signal. This is performed despite the presence of uncertainties and/or disturbances assumed to be as (12),(14) for the sliding mode algorithms, and to be as (25),(28),(29) for the e-LQR. For the earthquake application, this basically means that the friction coefficient must be Lipschitz w.r.t. the time and the states of the dynamical system. Such assumption is fulfilled by the most common friction laws used in fault mechanics (see [42] for a mathematical proof).

Now, some properties of each control strategy are discussed in the following, highlighting the differences between them. Sliding mode-based control:

* The origin of system (12), with bounded control coefficient $\beta(t, e)$ and Lipschitz w.r.t. the time perturbation $h(t, e)$ assumed as (14), is LFTS.

* Calculation of the gains for the 2-CTA (15) and 2-DIA (16) controllers, can be obtained using a Sum Of Square algorithm for the 2-CTA (see [25]) and by performing a maximization of homogeneous functions for the 2-DIA (see [45], [46]).

* The reference signal $r(t)$ that can be tracked, has to be chosen as $|r^{(3)}(t)| \leq \gamma$ with a positive constant $\gamma$, in order to fulfil the assumptions in (14). Note that the selected reference (6) fulfils this condition for all $t \in [0, t_{op}]$.

* Systems (17) and (18) are homogeneous vector-set of degree $d = -1$ and weights $(r_1, r_2, r_3) = (3, 2, 1)$. Due to homogeneity properties [47], the theoretical precision of the states after the transient are $|e_1| < \Delta_1 T_s^3$, $|e_2| < \Delta_2 T_s^2$ and $|e_3| < \Delta_3 T_s$, where $\Delta_i > 0$ with $i = \{1, ..., 3\}$ and $T_s$ the sampling time.

LQR-based control:

* The origin of system (25), with positive semi-definite variation coefficient $\Delta B(t, x_a)$ and Lipschitz w.r.t. the states $h(t, x_a)$ assumed as (28), is GES.

* Calculation of the gains for the e-LQR control (30) are obtained by solving the CARE (31)-(32).

* The classical version of an integral control will track a constant reference, i.e., $r(t) = r_0$ (see [21, Chapter 12], for example). According to the internal model principle (see [48]), the use of a double integrator (19) will be able to follow linear time references, i.e., $r(t) = \alpha_1 t + \alpha_2$, with some $\alpha_1, \alpha_2 \in \mathbb{R}$. The steady-state error tracking of the target reference (6) using the presented control will not be zero, but it will be improved by using this double integrator scheme. In addition, this error can become smaller increasing the e-LQR integral gains.

## IV. NUMERICAL SIMULATIONS AND EXPERIMENTAL CONFIRMATION

In order to illustrate the performance of both previous control algorithms, simulations and laboratory experiments have been performed based on the shifted system described by (4). The definition of the next parameters to be used are recalled

$$ \hat{N} = \frac{A}{m}, \quad \hat{k} = \frac{k}{m}, \quad \hat{\eta} = \frac{\eta}{m}, \quad A = L_{ac}^2, $$

$$ k = \frac{G}{L_{ac}}, \quad m = \rho L_{ac}^3, \quad \hat{\varphi}_e(x_1, x_2, t) = \frac{\varphi_e(x_1, x_2, t)}{m}, $$

where $\varphi_e(x_1, x_2, t)$ is an external perturbation and $\mu(x_1, x_2)$ the friction coefficient. In the following numerical examples, a friction coefficient of the form

$$ \mu(x_1) = \mu_{res} - \Delta\mu \cdot e^{-x_1/d_c}, \quad (38) $$

is considered with $\Delta\mu < 0$. Such function is defined as a slip-weakening friction law [35] and it evolves from an initial value
$\mu_{max} = \mu_{res} - \Delta\mu$ (static friction coefficient), to a residual one $\mu_{res}$ (kinetic friction coefficient) in a characteristic slip $d_c$, as shown by Fig. 4.

<figure><img src="image_6.png" /><figcaption>Fig. 4. Slip-weakening friction law used in the simulations.</figcaption></figure>

The choice of a friction law that depends only on the slip is argued by the difficulty of modelling the complex structure region of frictional interfaces. Nevertheless, the friction coefficient (38) fulfils the assumption imposed in (2) and it was used for calculating the controller gains.

Two cases have been studied in this paper. For the numerical simulations, the parameters of the reduced model for earthquakes (see eq. (4)) are chosen to correspond to an earthquake of magnitude $M_w = 5.8$, where simulations are being made using MATLAB Simulink. For the experiments, a different set of parameters was chosen corresponding to the laboratory-fault setup (Fig. 2), corresponding to an earthquake of magnitude $M_w = -5.2$. The mechanical and frictional parameters for these two cases are given in Table I.

**Remark 3:** Based on the Buckingham $\pi$ theorem [49], the lab-fault dynamics can be upscaled to obtain the real-fault dynamics under appropriate scaling laws (see [41, Appendix E] for further details).

Note how the seismic magnitude of the lab-fault allows to reproduce safely in the laboratory earthquake-like instabilities, but such experimental results can be upscaled to real earthquake events.

**Remark 4:** The controller gains used for the numerical simulations were obtained by upscaling the gains chosen for the experimental tests. This was made in order to be consistent with the upscale process commented above. See Appendix I for more details.

## A. Numerical Simulations

The numerical simulations on the shifted system (4) using the real-fault parameters from Table I, have been performed in MATLAB Simulink with Dormand-Prince's integration method of a fixed time step equal to $T_s^R = 184$ [ms].

The controller (11) with $\mu_0 = \mu_{res}^L, \hat{N}_0 = A^L/m^L, \nu$ chosen as (15) for the 2-CTA and as (16) for the 2-DIA, and the e-LQR control (30) have been implemented.<fnref n="1" /> The controllers gains have been selected as

* 2-CTA: $k_1 = 781.37, k_2 = 3.22 \times 10^3, k_3 = 4.51 \times 10^{-4}, k_4 = 2.15 \times 10^{-4}$ and $\lambda = 500$.

* 2-DIA: $k_{I1} = 5.21 \times 10^{-2}, k_{I2} = 3.23 \times 10^3, k_{I3} = 3.91 \times 10^{-4}, k_{I4} = 0$ and $\lambda = 500$.

* e-LQR: $k_1 = 1.88 \times 10^9, k_2 = 5.79 \times 10^8, k_3 = 1.02 \times 10^6$ and $k_4 = 18.52$.

The reference signal (6) is used with $d_{max} = 785$ [mm] and $t_{op} = 184.17$ [h]. Furthermore, an external perturbation $\hat{\phi}_e(x_1, x_2, t) = 3.2 \times 10^{-4} \sin(0.69t) + 3.2 \times 10^{-6} x_1$ [m/s²] has been added to the system, and the initial condition $x_1(0) = x_2(0) = 0$ was chosen for the three algorithms. It is important to notice that the external perturbation has a Lipschitz w.r.t. the time term and a Lipschitz w.r.t. the states term. Furthermore, the chosen frequency in the sinusoidal perturbation is equal to the natural frequency of the system, i.e., $\omega_n = \sqrt{k/m} \approx 0.69$ [rad/s]. The magnitude of both terms in such perturbation has been selected to appreciate a change in the system response and to check the robustness of the controllers versus different kinds of functions and resonance behaviour.

The numerical results are presented in Figs. 5-7. In the three simulations, the slip $x_1$ is followed the desired reference in order to dissipate slowly all the stored energy avoiding earthquake-like events. This can be seen by comparing Figs. 5 and 6, where the lack of a control input in the earthquake phenomenon results in the state $x_1$ evolving much faster than the controlled scenarios. Sliding-mode controllers present better results in the displacement tracking, but in terms of velocity, the e-LQR controller shows a more stable response without oscillations. Nevertheless, the three algorithms fulfil the tracking task despite the presence of the external perturbation and the use of nominal system parameters for the control design (robustness). In particular, the control signal generated by the e-LQR controller shows a smooth behaviour without demanding excess actuator response. On the other hand, the sliding-mode controllers (2-CTA and 2-DIA) show an oscillatory behaviour but not reaching high oscillations as the chattering effect (see zoom made in Fig. 7). This is because such sliding-mode algorithms make use of the discontinuous sign function, which provokes oscillations of high (in theory infinite) frequency. Nevertheless, the final control signal that the plant experience is continuous, thanks to the type of controllers employed. The above comment could explain the presence of more oscillations compared to the linear control, mainly in the velocity and in the control input.

<figure><img src="image_7.png" /><figcaption>Fig. 5. Earthquake phenomenon in the real-fault simulation.</figcaption></figure>

To make a further comparison between the three controllers, the mean integrated error (MISE), the average power (RMS), and the maximum errors have been calculated on the steady

<footnote n="1">¹Superscript *L* means lab-fault parameters from Table I.</footnote>
<table><caption>TABLE I
MECHANICAL AND FRICTIONAL PROPERTIES ADOPTED FOR THE REAL-SCALE SIMULATIONS (R) AND THE LABORATORY EXPERIMENTS (L)</caption><thead><tr><th>Parameter</th><th>Description</th><th>Simulations<br>Real-fault (R)</th><th>Experiments<sup>1</sup><br>Lab-fault (L)</th><th>Scaling factor<sup>2</sup><br>(R/L)</th></tr></thead><tbody><tr><td>ρ</td><td>density</td><td>2500 [kg/m<sup>3</sup>]</td><td>1385 [kg/m<sup>3</sup>]</td><td>λ<sub>ρ</sub> = 1.81</td></tr><tr><td>G</td><td>shear modulus</td><td>30 [GPa]</td><td>225.5 [kPa]</td><td>λ<sub>G</sub> = 1.33 × 10<sup>5</sup></td></tr><tr><td>η</td><td>damping coefficient</td><td>5 × 10<sup>14</sup> [kg/s]</td><td>408 [kg/s]</td><td>λ<sub>η</sub> = 1.23 × 10<sup>12</sup></td></tr><tr><td>L<sub>ac</sub></td><td>activated fault length</td><td>5 [km]</td><td>0.1 [m]</td><td>λ<sub>L<sub>ac</sub></sub> = 5 × 10<sup>4</sup></td></tr><tr><td>σ'<sub>n</sub></td><td>effective normal stress</td><td>50 [MPa]</td><td>0.1 [MPa]</td><td>λ<sub>p</sub> = 500</td></tr><tr><td>μ<sub>res</sub></td><td>residual friction</td><td>0.2353</td><td>0.4</td><td>λ<sub>μ</sub> = 0.5882</td></tr><tr><td>Δμ</td><td>friction drop</td><td>-0.1</td><td>-0.17</td><td>λ<sub>μ</sub></td></tr><tr><td>d<sub>c</sub></td><td>characteristic slip distance</td><td>276.35 [mm]</td><td>2.5 [mm]</td><td>λ<sub>δ</sub> = 110.54</td></tr><tr><td>d<sub>max</sub></td><td>maximum displacement</td><td>785 [mm]</td><td>7.1 [mm]</td><td>λ<sub>δ</sub></td></tr><tr><td>t<sub>op</sub></td><td>operation time</td><td>184.17 [h]</td><td>1 [h]</td><td>λ<sub>t</sub> = 184.17</td></tr><tr><td>T<sub>s</sub></td><td>sampling time</td><td>184 [ms]</td><td>1 [ms]</td><td>λ<sub>t</sub></td></tr><tr><td>M<sub>0</sub></td><td>seismic moment</td><td>6.25 × 10<sup>17</sup> [Nm]</td><td>17 [Nm]</td><td>λ<sub>M<sub>0</sub></sub> = 3.68 × 10<sup>16</sup></td></tr><tr><td>M<sub>w</sub></td><td>seismic magnitude</td><td>5.8</td><td>-5.2</td><td>λ<sub>M<sub>w</sub></sub> = 4.97</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Fig. 6. States in the real-fault simulations.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7. Control signal in the real-fault simulations.</figcaption></figure>

state (a time $t = t_{ss}$). All of these were normalized w.r.t. the e-LQR results. This can be seen in Fig. 8. The position results are smaller for the 2-DIA and the 2-CTA, but the e-LQR

velocity errors are smaller. This is consistent to the results shown in latter Figs. In addition, one can observe that the maximum position error for the e-LQR control is of the order of ~ 1 [mm], which can be considered negligible in a total slip of ~ 800 [mm]. On the other hand, the maximum velocity error obtained with the 2-CTA controller was of $3.5 \times 10^{-3}$ [m/s], which twenty times smaller than the maximum slip-rate of the original earthquake event (see Fig. 5). Furthermore, the three algorithms present the same RMS value, so they spent the same amount of energy (in average) to solve the tracking problem.

## B. Experimental Confirmation

In real applications, measurements of slip and slip-rate in system (1) can be obtained through microseismicity measurements and other geophysical methods (e.g., GPS, LIDAR and topographical measurements). For our case, the lab-fault system depicted in Fig. 2 measures the slip, $x_1$, using two LVDTs. In order to obtain the slip rate velocity $x_2$ required by the designed controllers, a robust exact filtering differentiator [50] has been implemented and reads as

$$ 
\begin{aligned}
\dot{w}_1 &= -5\lambda_d^{\frac{1}{5}}[w_1]^{\frac{4}{5}} + w_2, \\
\dot{w}_2 &= -10.03\lambda_d^{\frac{2}{5}}[w_1]^{\frac{3}{5}} + \hat{x}_1 - x_1, \\
\dot{\hat{x}}_1 &= -9.3\lambda_d^{\frac{3}{5}}[w_1]^{\frac{2}{5}} + \hat{x}_2, \\
\dot{\hat{x}}_2 &= -4.57\lambda_d^{\frac{4}{5}}[w_1]^{\frac{1}{5}} + \hat{x}_3, \\
\dot{\hat{x}}_3 &= -1.1\lambda_d^{0},
\end{aligned}
\quad (39) $$

with $λ_d = 1 \times 10^{-5}$. Therefore, all the designed control algorithms, $p$, have been implemented as functions of the estimated slip and slip rate, i.e., $p(\hat{x}_1, \hat{x}_2)$ instead of $p(x_1, x_2)$. This would reduce the possible noise from the estimations due to the LVDTs measurements.

**Remark 5:** Differentiator (39) provides the second derivative of the input $x_1(t)$ while filtering the signal with a second order filter, if $|x_1^{(3)}(t)| \le L_d$ and $λ_d > L_d$ (see proof

<footnote n="1">¹These values correspond to a double interface of an area of A = 100 [cm²] and a spring coefficient of k = 45.1 [N/mm] used for the experimental apparatus depicted in Fig. 2.</footnote>

<footnote n="2">²Due to the fact that seismic magnitude is based on a logarithm scale (see eq. (3)), the scaling factor between both faults is equal to $λ_{M_w} = \frac{2}{3} \log_{10} λ_{M_0} - 6.07$.</footnote>
<figure><img src="image_10.png" /><figcaption>Fig. 8. Errors and average power comparison in the simulations.</figcaption></figure>

in [50]). Furthermore, the authors of the present paper are fully aware of the separation principle problem with nonlinear systems. However, the purpose of the following results consists of experimentally showing the convergence of the complete closed loop-system (Plant, control and differentiator).

The designed controllers have been tested in the laboratory apparatus shown in Fig. 2. The experimental process is as follows. First, the spring is compressed by a press, taking the system close to its unstable equilibrium point. Without control, the stored elastic energy of the spring is released abruptly, causing an earthquake-like response (fast seismic slip), as shown in Fig. 9.

<figure><img src="image_11.png" /><figcaption>Fig. 9. Earthquake phenomenon in the lab-fault experiment.</figcaption></figure>

Then, the controller (11) with $\mu_0 = \mu_{res}^L$, $\hat{N}_0 = A^L/m^L$, $\nu$ chosen as (15) for the 2-CTA and (16) for the 2-DIA, and the e-LQR control defined by (30) have been implemented<fnref n="2" />. The lab-fault system parameters shown in Table I are considered as nominal values of the plant for calculating the controller gains as

* 2-CTA: $k_1 = 7.5$, $k_2 = 5$, $k_3 = 1.66 \times 10^{-4}$, $k_4 = 7.92 \times 10^{-5}$ and $\lambda = 500$.

* 2-DIA: $k_{I1} = 2$, $k_{I2} = 5$, $k_{I3} = 1.44 \times 10^{-4}$, $k_{I4} = 0$ and $\lambda = 500$.

* e-LQR: $k_1 = 4.17 \times 10^8$, $k_2 = 6.94 \times 10^5$, $k_3 = 4.17 \times 10^7$ and $k_4 = 1.39 \times 10^5$.

The reference signal to be tracked is chosen as (6), with $d_{max} = 7.5$ [mm] and $t_{op} = 1$ [h]. In other words, the controllers not only must avoid the unstable behaviour but to follow an aseismic slip equal to 7.5 [mm] in one hour, evolving according to the sigmoid function (6). The results are displayed in Figs. 10-11. The three controllers are able to fulfil the task of controlling an earthquake-like instability by increasing the response time of the system (from around 80 [ms] to 1 [h]). The three controllers achieve the steady state around 9 [min] and they keep the states around the reference signal. The e-LQR presents a large overshoot in the control signal at the beginning of the experiment (see Fig. 11) and the 2-DIA presents higher peaks on the velocity. Fig. 11 shows the three continuous control signals used for the tracking and how the measured friction coefficient evolves in time. Note how this latter is always bounded and always higher than the minimum value used for designing the controller ($\mu_{res} = 0.4$). It also shows high frequency oscillations compared with the nominal friction of Fig. (4), used in the numerical simulations. This behaviour is due to the grain size of the interfaces and other unmodelled dynamics, that the controllers successfully compensate.

Assuring a fair comparison between the three presented controllers is not straightforward, due to the presence of different uncertainties and disturbances in each experiment. In particular, the tested samples are not the same between different tests, the experiments do not initiate exactly in the same initial point, and better gains may exist to optimally tune the three different controllers. Furthermore, as discussed before, both control strategies present different theoretical properties. Nevertheless, the error comparison used in the

<footnote n="2">²Superscript *L* means lab-fault parameters from Table I.</footnote>
<figure><img src="image_12.png" /><figcaption>Fig. 10. States in the experimental confirmation.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11. Control signal and friction coefficient in the experimental confirmation.</figcaption></figure>

simulations has been made again for the experimental results. This is presented in Fig. 12. The lowest errors were obtained with the 2-CTA using the smallest average power, but the three algorithms present negligible absolute errors. A video of the experiment is available in [51].

V. CONCLUSIONS

In this paper, a challenging, emerging application of robust nonlinear control theory is presented for preventing natural and anthropogenic seismicity. The control is designed to perform tracking of a slow reference based on a reduced-order model for earthquakes. Two types of controllers are presented: the first one is based on sliding-mode theory and the

other on LQR control. The first one results in local finite-time convergence of the tracking error, while the second presents global exponential stability. Both controllers are designed to generate a continuous control signal and use integral action to compensate different kinds of perturbations. The algorithms have been tested and compared in numerical simulations over a real-fault system and in a specially designed experimental apparatus, showing that both types of controllers succeed in achieving tracking to a new stable equilibrium of lower energy. In the simulations, the best precision for the slip error was obtained with the sliding-mode algorithms, but the LQR control was better achieving a smaller error in the slip-rate error. With respect to the experimental tests, both controllers were successful in compensating unmodelled dynamics and parameter uncertainties present in the real systems, but the 2-CTA sliding-mode algorithm achieved the best results. As far as it concerns, the performance of the controllers could be upscaled to a real earthquake due to the existence of scaling laws between both faults. The design of controllers based on more detailed and complete models representing the earthquake phenomenon (e.g., coupled elastic and diffusion partial differential equations) remains as future work.

ACKNOWLEDGMENT

The authors would like to acknowledge the support of the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant agreement no. 757848 CoQuake).

APPENDIX I
GAIN SCALING

To show how the experimental results can be upscaled to a more realistic earthquake event, the lab-fault system parameters were scaled to obtain the real-fault parameters, according to Table I. Furthermore, the controller gains used in the experiments (IV-B) were manipulated in the same sense to obtain the set of gains (IV-A) for the numerical simulations.

Recalling (11), (15), (16) and (30), the controllers to be used in the lab-fault experiments are defined as

$$2\text{-CTA: } p^L = \left[-\lambda^{\frac{3}{2}} k_1 \left[e_1^L\right]^{\frac{1}{3}} - \lambda^{\frac{1}{2}} k_2 \left[e_2^L\right]^{\frac{1}{2}} + \xi_1^L\right] \frac{1}{\mu_0^L \hat{N}_0^L}, \\ \dot{\xi}_1^L = -\lambda k_3 \left[e_1^L\right]^0 - \lambda k_4 \left[e_2^L\right]^0,$$

$$2\text{-DIA: } p^L = - \left[ \lambda^{\frac{1}{2}} k_{I2} \left[ e_2^L \right]^{\frac{3}{2}} + \lambda^{\frac{1}{2}} k_{I1}^{\frac{3}{2}} e_1^L \right]^{\frac{1}{3}} + \xi_1^L \right] \frac{1}{\mu_0^L \hat{N}_0^L}, \\ \dot{\xi}_1^L = - \lambda k_{I3} \left[ e_1^L + \lambda^{-\frac{1}{2}} k_{I4} \left[ e_2^L \right]^{\frac{3}{2}} \right]^0,$$

$$\text{e-LQR: } p^L = -k_1 x_1^L - k_2 x_2^L - k_3 \xi_1^L - k_4 \xi_2^L, \\ \text{using the superscript } L \text{ as notation.}$$

According to Table I, the upscaled pressure, slip and slip-rate, error variables and integral terms can be obtained as $p^R = p^L \lambda_p$, $e_1^R = e_1^L \lambda_\delta$, $x_1^R = x_1^L \lambda_\delta$, $e_2^R = e_2^L \lambda_v$, $x_2^R = x_2^L \lambda_v$,
<figure><img src="image_14.png" /><figcaption>Fig. 12. Errors and average power comparison in the experiments.</figcaption></figure>

$\xi_1^R = \xi_1^L \lambda_\delta \lambda_t$, $\xi_2^R = \xi_2^L \lambda_\delta \lambda_t^2$, $\lambda_v = \lambda_\delta / \lambda_t$. Therefore, the latter controllers for the real-fault simulations read as

$$2\text{-CTA: } p^R = \left[-\lambda^{\frac{2}{3}} \bar{k}_1 \left[e_1^R\right]^{\frac{1}{3}} - \lambda^{\frac{1}{2}} \bar{k}_2 \left[e_2^R\right]^{\frac{1}{2}} + \xi_1^R\right] \frac{1}{\mu_0^L \bar{N}_0^L},$$
$$\dot{\xi}_1^R = -\lambda \bar{k}_3 \left[e_1^R\right]^0 - \lambda \bar{k}_4 \left[e_2^R\right]^0,$$

$$\text{where } (\bar{k}_1, \bar{k}_2, \bar{k}_3, \bar{k}_4) = \left( k_1 \frac{\lambda_p}{\lambda_\delta^{1/3}}, k_2 \frac{\lambda_p}{\lambda_v^{1/2}}, k_3 \frac{\lambda_p}{\lambda_t}, k_4 \frac{\lambda_p}{\lambda_v} \right),$$

$$2\text{-DIA: } p^R = - \left[ \lambda^{\frac{1}{2}} \bar{k}_{I2} \left[ \left[e_2^R\right]^{\frac{3}{2}} + \lambda^{\frac{1}{2}} \bar{k}_{I1}^{\frac{3}{2}} e_1^R \right]^{\frac{1}{3}} + \xi_1^R \right] \frac{1}{\mu_0^L \bar{N}_0^L},$$
$$\dot{\xi}_1^R = - \lambda \bar{k}_{I3} \left[ e_1^R + \lambda^{-\frac{1}{2}} \bar{k}_{I4} \left[e_2^R\right]^{\frac{3}{2}} \right]^0,$$

where

$$(\bar{k}_{I1}, \bar{k}_{I2}, \bar{k}_{I3}, \bar{k}_{I4}) = \left( k_{I1} \frac{\lambda_v}{\lambda_\delta^{2/3}}, k_{I2} \frac{\lambda_p}{\lambda_v^{1/2}}, k_{I3} \frac{\lambda_p}{\lambda_t}, k_{I4} \frac{\lambda_\delta}{\lambda_v^{3/2}} \right),$$

$$\text{e-LQR: } p^R = -\bar{k}_1 x_1^R - \bar{k}_2 x_2^R - \bar{k}_3 \xi_1^R - \bar{k}_4 \xi_2^R,$$

$$\text{where } (\bar{k}_1, \bar{k}_2, \bar{k}_3, \bar{k}_4) = \left( k_1 \frac{\lambda_p}{\lambda_\delta}, k_2 \frac{\lambda_p}{\lambda_v}, k_3 \frac{\lambda_p}{\lambda_\delta \lambda_t}, k_4 \frac{\lambda_p}{\lambda_\delta \lambda_v^2} \right),$$
and using the superscript *R* as notation.

## REFERENCES

[1] C. H. Scholz, *The Mechanics of Earthquakes and Faulting*. USA: Cambridge University Press, 2002.

[2] J. L. Rubinstein and A. B. Mahani, “Myths and facts on wastewater injection, hydraulic fracturing, enhanced oil recovery, and induced seismicity,” *Seismological Research Letters*, vol. 86, no. 4, pp. 1060–1067, 2015.

[3] K. M. Keranen, H. M. Savage, G. A. Abers, and E. S. Cochran, “Potentially induced earthquakes in Oklahoma, USA: Links between wastewater injection7 and the 2011 Mw 5.7 earthquake sequence,” *Geology*, vol. 41, no. 6, pp. 1060–1067, 2013.

[4] M. Zastrow, “South korea accepts geothermal plant probably caused destructive quake,” *Nature*, 2019.

[5] C. B. Raleigh, J. H. Healy, and J. D. Bredehoeft, “An Experiment in Earthquake Control at Rangely, Colorado,” *Science*, vol. 191, no. 4233, pp. 1230–1237, 1976.

[6] J. Fletcher and L. Sykes, “Earthquakes related to hydraulic mining and natural seismic activity in western New York State,” *Journal of Geophysical Research*, vol. 82, no. 26, pp. 3767–3780, 1977.

[7] Y. Guglielmi, F. Cappa, J.-P. Avouac, P. Henry, and D. Elsworth, “Seismicity triggered by fluid injection-induced aseismic slip,” *Science*, vol. 348, no. 6240, pp. 1224–1226, 2015.

[8] F. Cappa, M. M. Scuderi, C. Collettini, Y. Guglielmi, and J. P. Avouac, “Stabilization of fault slip by fluid injection in the laboratory and in situ,” *Science Advances*, vol. 5, no. 3, p. eaau4065, 2019.

[9] I. Stefanou, “Controlling anthropogenic and natural seismicity: Insights from active stabilization of the spring-slider model,” *Journal of Geophysical Research: Solid Earth*, vol. 124, no. 8, pp. 8786–8802, 2019.

[10] I. Stefanou, “Control instabilities and incite slow-slip in generalized burridge-knopoff models,” *arXiv:2008.03755*, 2020.

[11] I. Stefanou, G. Tzortzopoulos, and P. Braun, “Banc d'essai en cisaillement,” FR2109102, Patent Pending.

[12] Y. Ben-Zion and C. G. Sammis, “Characterization of fault zones,” *Pure and applied geophysics*, vol. 160, pp. 677–715, 2003.

[13] K. Brodie, D. Fettes, B. Harte, and R. Schmid, “Structural terms including fault rock terms,” tech. rep., Department of Earth Sciences, University of Manchester, UK, 2007.

[14] H. Rattez, I. Stefanou, and J. Sulem, “The importance of Thermo-Hydro-Mechanical couplings and microstructure to strain localization in 3D continua with application to seismic faults. Part I: Theory and linear stability analysis,” *Journal of the Mechanics and Physics of Solids*, vol. 115, pp. 54–76, 2018.

[15] H. Rattez, I. Stefanou, J. Sulem, E. Veveakis, and T. Poulet, “The importance of Thermo-Hydro-Mechanical couplings and microstructure to strain localization in 3D continua with application to seismic faults. Part II: Numerical implementation and post-bifurcation analysis,” *Journal of the Mechanics and Physics of Solids*, vol. 115, pp. 1–29, 2018.

[16] Z. Reches and D. A. Lockner, “Fault weakening and earthquake instability by powder lubrication,” *Nature*, vol. 467, no. 7314, pp. 452–455, 2010.

[17] N. Brantut and J. Sulem, “Strain Localization and Slip Instability in a Strain-Rate Hardening, Chemically Weakening Material,” *Journal of Applied Mechanics*, vol. 79, no. 3, p. 031004, 2012.

[18] E. Veveakis, S. Alevizos, and I. Vardoulakis, “Chemical reaction capping of thermal instabilities during shear of frictional faults,” *Journal of the Mechanics and Physics of Solids*, vol. 58, no. 9, pp. 1175–1194, 2010.

[19] E. Veveakis, I. Stefanou, and J. Sulem, “Failure in shear bands for granular materials: thermo-hydro-chemo-mechanical effects,” *Geotechnique Letters*, vol. 3, no. 2, pp. 31–36, 2013.

[20] N. A. Collins-Craft, I. Stefanou, J. Sulem, and I. Einav, “A Cosserat
Breakage Mechanics model for brittle granular media,” *Journal of the Mechanics and Physics of Solids*, vol. 141, p. 103975, 2020.

[21] H. Khalil, *Nonlinear Systems*. New Jersey, U.S.A.: Prentice Hall, 2002.

[22] V. Utkin, *Sliding modes in control and optimization*. Berlin, Germany: Springer Verlag, 1992.

[23] Y. Shtessel, C. Edwards, L. Fridman, and A. Levant, *Sliding Mode Control and Observation*. Intuitive theory of sliding mode control, New York, USA: Birkhauser, 2014.

[24] A. Chalanga, S. Kamal, and B. Bandyopadhyay, "Continuous integral sliding mode control: A chattering free approach," in *IEEE International Symposium on Industrial Electronics*, (Taipei, Taiwan), 2013.

[25] V. Torres-González, T. Sánchez, L. Fridman, and J. Moreno, "Design of continuous twisting algorithm," *Automatica*, vol. 80, pp. 119–126, 2017.

[26] J. Moreno, "Discontinuous integral control for mechanical systems," in *International Workshop on Variable Structure Systems*, (Nanjing, China), 2016.

[27] R. E. Kalman, "Contributions to the Theory of Optimal Control," *Boletín de la Sociedad Matemática Mexicana*, vol. 5, pp. 102–119, 1960.

[28] A. Filippov, *Differential Equations with Discontinuous Right-hand Sides*. Dordrecht, The Netherlands: Kluwer Academic Publishers, 1988.

[29] S. Bhat and D. Bernstein, "Finite-time stability of homogeneous systems," in *American Control Conference*, (New Mexico, USA), pp. 2513–2514, 1997.

[30] A. Baccioti and L. Rosier, *Lyapunov functions and stability in control theory*. New York: Springer-Verlag, 2nd ed., 2005.

[31] E. Bernuau, D. Efimov, W. Perruquetti, and A. Polyakov, "On homogeneity and its application in sliding mode control," *Journal of the Franklin Institute*, vol. 351, no. 4, pp. 1866–1901, 2014.

[32] V. Andrieu, L. Praly, and A. Astolfi, "Homogeneous approximation, recursive observer design and output feedback," *SIAM J. Control Optim.*, vol. 47, no. 4, pp. 1814–1850, 2008.

[33] E. Cruz-Zavala and J. Moreno, "Homogeneous high order sliding mode design: a Lyapunov approach," *Automatica*, vol. 80, pp. 232–238, 2017.

[34] A. Mercado-Uribe. and J. A. Moreno, "Output feedback discontinuous integral controller for siso nonlinear systems," in *15th International Workshop on Variable Structure Systems*, (Graz, Austria), 2018.

[35] H. Kanamori and E. E. Brodsky, "The physics of earthquakes," *Reports on Progress in Physics*, vol. 67, no. 8, pp. 1429–1496, 2004.

[36] G. Tzortzopoulos, P. Braun, and I. Stefanou, "Absorbent porous paper reveals how earthquakes could be mitigated," *Geophysical Research Letters*, vol. 48, no. 3, p. e2020GL090792, 2021.

[37] J. Huang and D. L. Turcotte, "Chaotic seismic faulting with a mass-spring model and velocity-weakening friction," *Pure and Applied Geophysics* PAGEOPH, vol. 138, no. 4, pp. 569–589, 1992.

[38] J. H. Dieterich, "Constitutive Properties of Faults With Simulated Gouge," in *Mechanical Behavior of Crustal Rocks*, vol. 24, pp. 103–120, 1981.

[39] A. Ruina, "Slip instability and state variable friction laws," *Journal of Geophysical Research: Solid Earth*, vol. 88, no. B12, pp. 10359–10370, 1983.

[40] J. H. Dieterich, "Modeling of rock friction: 1. experimental results and constitutive equations," *Journal of Geophysical Research: Solid Earth*, vol. 84, no. B5, pp. 2161–2168, 1979.

[41] G. Tzortzopoulos, *Controlling earthQuakes (CoQuake) in the laboratory using pertinent fault stimulating techniques*. PhD thesis, Ecole Centrale de Nantes, France, 2021.

[42] I. Stefanou and G. Tzortzopoulos, "Preventing instabilities and inducing controlled, slow-slip in frictionally unstable systems," *Journal of Geophysical Research: Solid Earth*, vol. 127, no. 7, p. e2021JB023410, 2022.

[43] P. Braun, G. Tzortzopoulos, and I. Stefanou, "Design of Sand-Based, 3-D-Printed Analog Faults With Controlled Frictional Properties," *Journal of Geophysical Research: Solid Earth*, vol. 126, p. e2020JB020520, 2021.

[44] D. Gutiérrez-Oribio, I. Stefanou, and F. Plestan, "Passivity-based control of underactuated mechanical systems with Coulomb friction: Application to earthquake prevention," *arXiv:2207.07181*, 2022.

[45] D. Gutiérrez-Oribio, J. A. Mercado-Uribe, J. A. Moreno, and L. Fridman, "Robust global stabilization of a class of underactuated mechanical systems of two degrees of freedom," *International Journal of Robust and Nonlinear Control*, vol. 31, no. 9, pp. 3439–3454, 2021.

[46] J. A. Mercado-Uribe and J. A. Moreno, "Discontinuous integral action for arbitrary relative degree in sliding-mode control," *Automatica*, vol. 118, p. 109018, 2020.

[47] A. Levant, "Homogeneity approach to high-order sliding mode design," *Automatica*, vol. 41, no. 5, pp. 823–830, 2005.

[48] B. A. Francis and W. M. Wonham, "The internal model principle of control theory," *Automatica*, vol. 12, no. 5, pp. 457–465, 1976.

[49] J. D. Logan, *Applied Mathematics*. USA: John Wiley & Sons, 2013.

[50] A. Levant and M. Livne, "Robust exact filtering differentiators," *European Journal of Control*, vol. 55, pp. 33–44, 2020.

[51] "Controlling earthquake-like instabilities in the laboratory: LIVE experiment presented during the PhD defense of Dr. George Tzortzopoulos," 2021. CoQuake website.

<figure><img src="image_15.png" /></figure>

Diego Gutiérrez-Oribio was born in México and received his M.S. degree and Bachelor degree in Mechatronic Engineering from the National Autonomous University of Mexico (UNAM), Mexico City, Mexico in 2016 and 2013, respectively. In 2021, he received his Ph.D. in Electrical Engineering - Automatic Control (with honours) at UNAM. He is currently a post-doctoral member in the GeM Laboratory in the Ecole Centrale, Nantes, France. His current research of interests include the control of mechanical systems, robust and nonlinear control, continuous higher-order sliding-mode control and control of earthquake phenomena.

<figure><img src="image_16.png" /></figure>

Georgios Tzortzopoulos does his PhD thesis on earthquake control using pertinent fault stimulating experimental techniques in the laboratory. He graduated in 2016 from the National Technical University of Athens (NTUA) holding a Diploma in Civil Engineering and a Master's Degree in Structural Engineering. He did his Master Thesis at the Laboratory of Structural Analysis and Antiseismic Research of NTUA. In the meantime, he fulfilled his military obligations in Greece. His current research interests are oriented towards fault mechanics and earthquake control using advanced experimental techniques.

<figure><img src="image_17.png" /></figure>

Ioannis Stefanou is Professor at the Ecole Centrale de Nantes (ECN), France. He studied civil engineering, mechanics and applied mathematics at the National Technical University of Athens. He then did his PhD thesis at the Laboratory of Geomaterials of the same institution. His main research topics are the mechanics of (geo-)materials, structural dynamics, geomechanics, fault reactivation and earthquake control. He is the PI of the ERC-StG project "Controlling earthQuakes - CoQuake" (www.coquake.eu) and of the Connect Talent project "Controlling Extreme EVents - CEEV" (www.blastructures.eu) awarded by the Pays de la Loire.

<figure><img src="image_18.png" /></figure>

Franck Plestan received the Ph.D. in Automatic Control from the Ecole Centrale de Nantes, France, in 1995. From September 1996 to August 2000, he was with Louis Pasteur University, Strasbourg, France. In September 2000, he joined the Ecole Centrale de Nantes, Nantes, France where he is currently Professor. His research interests include robust control (adaptive/higher order sliding mode, time-delay systems) and nonlinear observer design. He is also working in several application domains as pneumatic actuators, automotive, flying systems, renewable energy systems.