# Advances in sliding mode control of earthquakes via boundary tracking of wave and diffusion PDEs

Diego Gutierrez-Oribio, Yury Orlov, Ioannis Stefanou, Franck Plestan

► To cite this version:

Diego Gutierrez-Oribio, Yury Orlov, Ioannis Stefanou, Franck Plestan. Advances in sliding mode control of earthquakes via boundary tracking of wave and diffusion PDEs. 16th International Workshop on Variable Structure Systems and Sliding Mode Control, Sep 2022, Rio de Janeiro, Brazil. 10.1109/vss57184.2022.9902111. hal-03774399

HAL Id: hal-03774399

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

Submitted on 10 Sep 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Advances in Sliding Mode Control of Earthquakes via Boundary Tracking of Wave and Diffusion PDEs

Diego Gutiérrez-Oribio<fnref n="1" />, Yury Orlov<fnref n="2" />, Ioannis Stefanou<fnref n="1" /> and Franck Plestan<fnref n="3" />

**Abstract**— Two recent results on earthquake control are summarized. A simplified model of an earthquake phenomenon is addressed by means of a cascade system of a 1D wave equation, representing the fault slip and wave propagation, and a 1D diffusion equation, representing the actuator dynamics as a diffusion process. In order to avoid a fast slip (earthquake-like behaviour), the control is designed to follow a slow reference. The control strategies are presented separately for both Partial Differential Equations (PDEs). For the wave PDE, a homogeneous boundary tracking control is developed to achieve exponential Input to State Stability (eISS) of the error closed-loop dynamics. For the diffusion PDE, Proportional Integral (PI) and a discontinuous integral term are coupled to exponentially stabilize the error origin despite model uncertainties and perturbations. Simulations are additionally conducted to support the robustness and stability properties of the proposed control algorithms, by separately, obtaining critical remarks that will lead to the design of the single control for the cascade system in a future stage.

## I. INTRODUCTION

Earthquakes are dynamic instabilities that occur in the earth's crust [1]. Most of earthquakes are of natural causes. 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). Such earthquake instabilities are characterized by important nonlinearities due to friction. Such frictional term can present many uncertainties and unmodelled dynamics to be compensated by an appropriate control action. There are few works dealing with the control of earthquake phenomena (see [5], [6]) using a mathematical model of reduced order that can establish the average behaviour of a single earthquake fault. Nevertheless, the wave PDE model describes more precisely a fault slip in one dimension due to the inclusion of the wave propagation phenomenon, whereas the diffusion PDE model describes the evolution of the input pressure applied at the fault through diffusion process (actuator dynamics).

The wave equation is a second-order hyperbolic PDE used to model the propagation of waves in elastic medium. This may concern, for instance, acoustic, electromagnetic and mechanical waves. The boundary control design of such

systems has been addressed with the backstepping approach [7], [8], [9] and $H_\infty$ control [10], for example. Reference tracking has been not often addressed. It has been solved in [11], [12] only, using flatness-based control and adaptive control, respectively, and that calls for further investigation.

The diffusion equation is a first-order parabolic PDE used to describe heat and other diffusion processes such as chemical reactions, population growth, market price fluctuations, and fluids flow to name a few. The motion planning has been addressed for the diffusion equation using a flat output capable of parameterising the system state evolution [8], [13], [14], [15]. Such approach provides an open-loop control able to track the desired reference but only for the nominal system with no uncertainties and disturbances and starting at some specific initial condition. In order to robustify such a tracking control, one may invoke a feedback control. Boundary tracking control of the diffusion equation has been based on output feedback [16], [17] and backstepping [8], [15], to name a few.

It is well-known that real systems present uncertainties/disturbances in their dynamics, resulting in the necessity of robust control which is able of attenuating such phenomena. Sliding-mode control has been recognized since a long time as a powerful control method to counteract non-vanishing external disturbances and unmodelled dynamics, even for infinite dimension systems (see, e.g., [18] and a very recent monograph [19]). In [20] a dynamic sliding-mode control was designed to achieve the uncertain wave equation stabilization; in [21] sliding-mode techniques were used to perform the tracking of the wave equation, but both with distributed control input. With respect to the sliding-modes control in the diffusion equation, the design of robust controllers was made in [22] and the tracking task was solved in [21] compensating Lipschitz-in-time disturbances by using distributed control.

This paper reviews recent advances in the earthquake control task based on two recent results [23] and [24]. A simplified earthquake model is considered as a coupling between a 1D wave equation and a 1D diffusion equation. The control design is addressed separately, designing a tracking boundary control for each PDE. For the wave equation case, a homogeneous algorithm achieves exponential eISS of the error closed-loop system. A PI feedback plus a discontinuous component for the diffusion equation obtains global exponential stability of the error closed-loop system despite the presence of uncertainties and/or disturbances. Both tracking control strategies are tested in simulations, supporting the reliability and robustness of the controllers.

<footnote n="1">¹ Nantes Université, École Centrale Nantes, France
CNRS, UMR 6183, F-44000 Nantes, (diego.gutierrez-oribio@ec-nantes.fr and ioannis.stefanou@ec-nantes.fr)</footnote>

<footnote n="2">² CICESE Research Center, 22860 Ensenada, México (yorlov@cicese.mx)</footnote>

<footnote n="3">³ Nantes Université, École Centrale Nantes, CNRS, LS2N, UMR 6004, F-44000 Nantes, France (franck.plestan@ec-nantes.fr)</footnote>
This will lead to the design of the single control for the cascade system in future works.

The outline of this work is as follows. The simplified earthquake model and the control objective are introduced in Section II. The design of the boundary tracking control for the wave equation and its simulations are given in Section III, whereas the tracking boundary control for the diffusion equation is accompanied with simulations in Section IV. The proposed tracking strategies are discussed in Section V, and concluding remarks are collected in Section VI.

**Notation:** Define the function $[\cdot]^{\gamma} := |\cdot|^{\gamma}\text{sign}(\cdot)$, for any $\gamma \in \mathbb{R}_{\ge 0}$. $H^l(a, b)$ with $a \le b$ and $l = \{0, 1, 2, ...\}$ denotes the Sobolev space of absolutely continuous scalar functions $u(x)$ on $(a, b)$ with square integrable derivatives $u^{(i)}(x)$ up to the order $l$. $H^l$-norm reads as $||u(\cdot)||_{H^l(a,b)} = \sqrt{\int_a^b \sum_{i=0}^l [u^{(i)}(x)]^2 dx}$. The symbols $||u(\cdot)||_{H^0(a,b)} = ||u(\cdot)||_{L^2(a,b)} = ||u(\cdot)||$ are also used as throughout the paper. Finally, $u_x = \partial u/\partial x$ and $u_{xx} = \partial^2 u/\partial x^2$ are in play as well.

## II. PROBLEM STATEMENT

The physical process of an earthquake is modelled by a three dimensional coupled PDE (see Biot theory, [25]). Here we simplify the problem to 1D by assuming invariance in the spatial variables belonging to the fault plane. Moreover, we consider that fluid injection rates are low. As a result, the physical process can be described by a 1D wave equation and a 1D diffusion equation, which are coupled through their BC on the seismic fault. This coupling is owed to the dependency of the fault friction (1D wave equation Neumann BC) to the fluid pressure (1D diffusion equation Dirichlet BC) at $x = 0$, as seen in Fig. 1 and system (1).

$$ 
\begin{aligned}
u_{tt}(x,t) &= c^2 u_{xx}(x,t), \\
u_x(0,t) &= \frac{\mu(u(0,t), u_t(0,t), t) - \mu_0}{G} \sigma' \\
&\quad - \frac{\mu(u(0,t), u_t(0,t), t)}{G} p(0,t), \quad u(D,t) = 0, \\
p_t(x,t) &= dp_{xx}(x,t), \\
p_x(0,t) &= 0, \quad p_x(D_2,t) = q + \varphi(t).
\end{aligned} 
\quad (1) $$

The first three equations of (1) (1D wave equation) describe fault slip on the boundary $x = 0$ and the wave propagation due to deformation propagation along $x$ until the attached point $x = D > 0$. $[u(x,t), u_t(x,t)]$ is its state vector evolving in the space $H^1(0, D) \times H^0(0, D)$ (displacement and velocity, respectively). $x \in [0, D]$ is the space variable, $t \ge 0$ is the time variable, $c = \sqrt{G/\rho}$ is the shear velocity, where $G$ is the shear modulus and $\rho$ the density of the surrounding to the fault rocks. $\mu(u(0,t), u_t(0,t), t)$ is the friction coefficient and $\mu_0$ is its initial constant value<fnref n="1" />. On the frictional interface (fault), the difference between the constant effective “overburden” normal stress, $\sigma'$, and the pressure of the injected or extracted fluid in/from the

fault-like interface, $p(0, t)$, is applied (Terzaghi's principle of effective stress [26]), which justifies the boundary condition at $x = 0$. More details of this system are given in [5], [27], [28].

The left-hand boundary represents the coupling with the next three equations (1D diffusion equation) due to the presence of an actuator dynamics of the input pressure. $p(x, t)$ is its state vector (pressure) evolving in the space $H^1(0, D_2)$, $x \in [0, D_2]$ is the space variable, $t \ge 0$ is the time variable and $d$ is the thermal diffusivity. A well is located at the boundary $x = D_2 > 0$ where the flux $q$ (control input) is applied. The term $\varphi$ represents external disturbances acting on the system. $\varphi$ is supposed to be sufficiently smooth and to fulfil $|\dot{\varphi}| \le L$, with some a priori known constant $L$.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Simplified earthquake model, a cascade system of a 1D wave equation (blue) and a 1D diffusion equation (red). $u(x, t)$ denotes displacement (wave) and $p(x, t)$ denotes fluid pressure.</figcaption></figure>

It is well-known (see, e.g., [19, Chapter 3]) that for arbitrary initial conditions $u(x, 0)$, $u_t(x, 0)$ of class $H^1(0, D) \times H^0(0, D)$ and $p(x, 0)$ of class $H^1(0, D_2)$, there exists a mild solution of the open-loop boundary value problem (BVP) (1). Throughout the paper, only mild solutions of the corresponding BVP are in play.

Without control, system (1) presents a fast slip behaviour (earthquake). To prevent this, one should adequately adjust the well flux, $q$, to obtain a desired fluid pressure at the fault, $p(0, t)$, by employing control techniques in order to stabilize the system (1) origin and/or track a reference input signal, releasing the stored elastic energy smoothly and not abruptly. Therefore, the objective of this work is to design a control input $q$ of the diffusion equation, which drives the boundary states $u(0, t)$ and $u_t(0, t)$ of the wave equation, to a vicinity of desired references $r(t)$, $\dot{r}(t)$, respectively, despite the presence of uncertainties in the system parameters and/or disturbances $\varphi$. For this purpose, the tracking control design is performed separately for the wave and diffusion equations and supported with simulation evidences for each case.

## III. TRACKING FOR THE WAVE EQUATION

In order to drive the state $u(0, t)$ to a desired reference signal $r(t)$ on the wave equation

$$ 
\begin{aligned}
u_{tt}(x,t) &= c^2 u_{xx}(x,t), \\
u_x(0,t) &= \frac{\mu - \mu_0}{G} \sigma' - \frac{\mu}{G} p, \quad u(D,t) = 0,
\end{aligned} 
\quad (2) $$

<footnote n="1">¹From now on, the dependency on $(u(0, t), u_t(0, t), t)$ will be neglected for ease of reference.</footnote>
the error $w(x, t)$ is defined as

$$w(x,t) = u(x,t) - \frac{D-x}{D}r(t). \quad (3)$$

Note that if $w(x,t) = 0$, then the main tracking goal, $u(0,t) = r(t)$, will be solved. Moreover, this error is also consistent with the boundary condition $u(D,t) = 0$ of the system (2).

Then, the error dynamics of the wave equation reads as

$$\begin{aligned} w_{tt}(x,t) &= c^2 w_{xx}(x,t) - \frac{D-x}{D} \ddot{r}(t), \\ w_x(0,t) &= \frac{\mu - \mu_0}{G} \sigma' - \frac{\mu}{G} p + \frac{1}{D} r(t), \quad w(D,t) = 0. \end{aligned} \quad (4)$$

### A. Control Design

Define the control $p$ as

$$\begin{aligned} p &= -\frac{G}{\mu} \left( \nu - \frac{1}{D} r(t) - \frac{\mu - \mu_0}{G} \sigma' \right), \\ \nu &= \lambda_1 [w(0,t)]^a + \lambda_2 [w_t(0,t)]^a. \end{aligned} \quad (5)$$

with $a \in [0, 1]$ and $\lambda_1, \lambda_2$ positive gains to be chosen.

The proposed control (5) is composed by feedback linearisation terms and homogeneous algorithm $\nu$ inspired by [29]. When $a = 0$ this algorithm is nothing else than the well-known Twisting Controller [30] and is governed by discontinuous (multi-valued) right-hand side in the two manifolds $w(0,t) = 0$ and $w_t(0,t) = 0$. Solutions of the distributed parameter system (4) driven by this discontinuous controller are understood in the Filippov's sense [31]. Extension of the Filippov concept towards the infinite-dimensional setting may be found in [19]. The present paper focuses on the tracking synthesis whereas the well-posedness analysis of the closed-loop system (4), (5) is similar to that of [21] and it remains beyond the scope of the paper. Thus, for the closed-loop system in question, it is simply assumed that it possesses a unique Filippov's solution.

The stability of the error dynamics (4) is then in force.

**Theorem 1:** System (4) fulfils an eISS condition<fnref n="2" /> w.r.t. the $H^1(0, D) \times H^0(0, D)$ norm

$$\Gamma_w(t) = \left( |w(0,t)|^{1+a} + ||w_x(x,t)||^2 + \frac{1}{c^2} ||w_t(x,t)||^2 \right)^{1/2}$$

if the controller gains on (5) are tuned as $\lambda_1 > \frac{\lambda_2}{1+a}$, $\lambda_2 > 0$.

**Remark 1:** The proof of Theorem 1 and more details about it, can be found in [23]. Furthermore, the uncertain system control design, i.e. the system (2) with uncertainties in the parameters and perturbations matched to the control, is under preparation.

### B. Numerical Simulations

The friction coefficient chosen to perform numerical simulations is a slip-weakening friction law (see [33]) $\mu = \mu_{res} - \Delta\mu \cdot e^{-u(0,t)/d_c}$, with $\Delta\mu < 0$. A fault is reactivated when it reaches its static coefficient of friction, $\mu_{max} = \mu_{res} - \Delta\mu$,

following by a stress drop $\Delta\mu$ in a characteristic distance $d_c$, leaving the residual value $\mu_{res}$ as friction coefficient.

The desired reference is a smooth function chosen as $r(t) = d_{max}s^3(10 - 15s + 6s^2)$, where $s = t/t_{op}$, $d_{max}$ is the target displacement and $t_{op}$ is the operational time of the controller. The later is chosen in order to modify the speed and acceleration of the reference. The reference parameters were chosen as $d_{max} = 0.823$ [m] and $t_{op} = 600$ [s].

The system parameters are $\rho = 2500$ [kg/m³] (surrounding rock density), $G = 30$ [GPa] (shear modulus), $\sigma' = 50$ [MPa] (effective normal stress), $\mu_{res} = 0.5$ (residual friction), $\Delta\mu = -0.1$ (friction drop), $d_c = 10$ [mm] (characteristic slip distance) and $D = 5$ [km] (fault length).

The control (3), (5) and the system (2) have been implemented in Matlab Simulink with the Runge-Kutta's integration method of fixed step and a sampling time equal to 1 [ms]. The wave equation is implemented using the finite-differences approximation technique, discretizing the spatial domain $x \in [0, D]$ into 51 ordinary differential equations (ODE).

Three different values of $a$ in (5) were chosen to compare the results. The respective gains $\lambda_1, \lambda_2$ were selected according to condition in Theorem 1:

* $a = 0$ with $\lambda_1 = 1.1 \times 10^{-8}$ and $\lambda_2 = 1 \times 10^{-8}$.

* $a = 0.5$ with $\lambda_1 = 7 \times 10^{-9}$ and $\lambda_2 = 1 \times 10^{-8}$.

* $a = 1$ with $\lambda_1 = 5.1 \times 10^{-9}$ and $\lambda_2 = 1 \times 10^{-8}$.

In order to evaluate the importance of the control implementation, a scenario without control (like a "standard" earthquake) is also presented. The results can be seen in Figures 2-4. The three cases drive the states $u(0,t)$, $u_t(0,t)$ close to the desired reference, dissipating the earthquake energy slowly in comparison with the control-free scenario. The case $a = 0$ achieves a much better precision for $w(0,t)$ but not for $w_t(0,t)$. This is due to the presence of high-frequency oscillations in $w_t(0,t)$ due to the use of discontinuous control when $a = 0$. Nevertheless, the control signal $p$ (Fig. 4) is not affected by a chattering effect due to its small gains.

## IV. TRACKING FOR THE DIFFUSION EQUATION

Consider the diffusion equation

$$\begin{aligned} p_t(x,t) &= dp_{xx}(x,t), \\ p_x(0,t) &= 0, \quad p_x(D_2,t) = q + \varphi(t), \end{aligned} \quad (6)$$

where the parameters $d, D_2$ are assumed to be uncertain, but bounded as $0 < d_m \le d \le d_M$, $0 < D_m \le D_2 \le D_M$, by some known constant $d_m, d_M, D_m, D_M$.

The objective here is to design a control input $q$ able to drive the output $y(t) = p(0,t)$ of the underlying BVP (6) to a desired reference $r_p(t)$ despite the presence of uncertainties and/or disturbances. Furthermore, the measurement of the state at $x = D_2$ is assumed to be available for the control design, i.e., $p(D_2,t)$. Yet, the exact value of the parameter $D_2$ is uncertain and only a nominal value, $D_n$, is on the use. This is because the exact fault location in the cascade system is unknown (the point $x = 0$ in Fig. 1).

<footnote n="2">²See [32] for more details on ISS in infinite-dimensional setting.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 2. Slip and slip-rate at the fault ($x = 0$) versus time (sec).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Error variables at the fault ($x = 0$) versus time (sec).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Control signal $p$ versus time (sec).</figcaption></figure>

### A. Motion Planning

Following [13],[8, Chapter 12], the reference trajectory generation for the diffusion equation (6) can be obtained through its flat output $y(t) = p(0, t)$. The trajectory generation is further performed for the unperturbed system with

$$ \varphi(t) \equiv 0. $$

Thus, the reference state trajectory is specified to

$$ \bar{p}(x, t) = \sum_{i=0}^{\infty} r_{p}^{(i)}(t) \frac{x^{2i}}{d^{i}(2i)!}, \quad (7) $$

and the nominal input signal from the BC at $x = D$ is defined as

$$ \bar{q} = \sum_{i=1}^{\infty} r_{p}^{(i)}(t) \frac{D_{2}^{2i-1}}{d^{i}(2i-1)!}. \quad (8) $$

In order to guarantee that $p(x, t) \to \bar{p}(x, t)$ and $y(t) \to r_p(t)$, the convergence of (7) has to be guaranteed. The next theorem states the conditions to be imposed on the to-be-tracked reference in order to achieve convergence.

**Theorem 2:** The series (7) is absolutely convergent provided that the output reference $r_p(t)$ fulfils

$$ \sup_{t \ge 0} |r_p^{(i+1)}(t)| \le \frac{2d}{D_2^2} (i+1) \sup_{t \ge 0} |r_p^{(i)}(t)| \quad \forall \ i \in \mathbb{Z}_{\ge 0}. \quad (9) $$

### B. Control Design

The output $y(t) = p(0, t)$ to-be-tracked matches the reference trajectory $r_p(t)$ only if the state trajectory satisfies the same initial condition as that of the plant, i.e., if $p(x, 0) = \bar{p}(x, 0)$. Furthermore, the nominal control (8) is not able to compensate any non-trivial disturbance $\varphi(t) \ne 0$. Next, the tracking problem of interest is addressed for the arbitrarily initialized state to-be-tracked trajectory in the presence of matched disturbances. Setting the state deviation

$$ \tilde{p}(x, t) = p(x, t) - \bar{p}(x, t), \quad (10) $$

from the reference trajectory $\bar{p}(x, t)$ (7), the error dynamics (10) is governed by

$$ \begin{aligned} \tilde{p}_t(x, t) &= d\tilde{p}_{xx}(x, t), \\ \tilde{p}_x(0, t) &= 0, \quad \tilde{p}_x(D_2, t) = \tilde{q} + \varphi(t), \end{aligned} \quad (11) $$

where the control input $q = \tilde{q} + \bar{q}$ is precomposed in terms of the reference input $\bar{q}$, determined by (8), and a new virtual component $\tilde{q}$. Now the tracking objective for the nominal reference (7) is reduced to the design of the virtual input $\tilde{q}$ exponentially stabilizing the error (10) at the origin.

Consider the following control law $\tilde{q}$ defined as

$$ \begin{aligned} \tilde{q} &= -\frac{\lambda_1}{D_n} \tilde{p}(D_2, t) + \nu, \\ \dot{\nu} &= -\lambda_2 \tilde{p}(D_2, t) - \lambda_3 [\tilde{p}(D_2, t)]^0, \end{aligned} \quad (12) $$

where $\lambda_1, \lambda_2, \lambda_3$ are the gains that must be tuned and $D_n$ is the nominal value of $D_2$.

The proposed feedback law (12) is composed by a PI control and a discontinuous term passing through an integrator. It generates a continuous control signal despite a discontinuous (multi-valued) right-hand side in the manifold $\tilde{p}(D_2, t) = 0$. The closed-loop system in question is assumed to admit a Filippov solution as for the wave equation case.
The closed-loop system (11), driven by (12), reads as

$$
\begin{align}
\tilde{p}_t(x,t) &= d\tilde{p}_{xx}(x,t), \nonumber \\
\tilde{p}_x(0,t) &= 0, \quad \tilde{p}_x(D_2,t) = -\frac{\lambda_1}{D_n}\tilde{p}(D_2,t) + \delta, \tag{13} \\
\dot{\delta} &= -\lambda_2\tilde{p}(D_2,t) - \lambda_3[\tilde{p}(D_2,t)]^0 + \dot{\varphi}(t), \nonumber
\end{align}
$$

with $\delta = \nu + \varphi$, being substituted into (12) for $\nu = \delta - \varphi$
for deriving $\delta$-dynamics (13). The following result is then in
force.

**Theorem 3:** Let the error dynamics (11) be driven by the feedback $\tilde{q}$, governed by (12) and tuned in accordance with $\lambda_1 > \frac{D_n}{D_m}$, $\lambda_2 > d_M$, $\lambda_3 > L$. Then, the closed-loop error dynamics (13) are globally exponentially stable despite the presence of any boundary disturbance $\varphi$ satisfying $|\dot{\varphi}| \le L$.

**Remark 2:** The proof of Theorems 2, 3 and more details about them, can be found in [24].

C. Numerical Simulations

The proposed control strategy $q = \tilde{q} + \bar{q}$ with (8),(12) has been implemented in the system (6) using Matlab Simulink with Euler's integration method of fixed step and a sampling time equal to 50 [ms]. The diffusion equation was implemented using $d = 0.05$, $D_2 = 10$ and $\varphi(t) = 0.01t + 2\sin(t)$ and the finite-differences approximation technique, discretizing the spatial domain $x \in [0, D_2]$ into 51 ordinary differential equations (ODE). The initial condition was set in $p(x, 0) = 0$.

The tracking reference $r_p(t)$ and the corresponding nom-
inal reference $\tilde{p}(x, t)$ and control $\tilde{q}$, are defined depending
on the simulation time and expressions (7)-(8):

1) $t \in [0, 1 \times 10^4)$ [s]: $r_p(t) = At$, $A = 1 \times 10^{-3}$.

2) $t \in [1 \times 10^4, 2 \times 10^4)$ [s]: $r_p(t) = 20$.

3) $t \in [2 \times 10^4, 3 \times 10^4)$ [s]:
$$r_p(t) = Ae^{[-\beta(t-2 \times 10^4)]}, \quad A = 20, \beta = 1 \times 10^{-3}.$$

4) $t \in [3 \times 10^4, 4 \times 10^4]$ [s]: $r_p(t) = A \sin(\omega t)$, $A = 10$, $\omega = 1 \times 10^{-3}$.

All of the proposed output references fulfil the condition (9). The respective gains of the error control $\tilde{q}$ have been selected according to Theorem 3 as $\lambda_1 = 1$, $\lambda_2 = 10$, $\lambda_3 = 2.5$ and $D_n = 9$, $D_m = 9$, $d_M = 0.06$. The results are displayed in Figs. 5-7. The main objective of stabilizing the output $y = p(0, t)$ over the four successive references $r_p(t)$ is achieved despite the presence of the unbounded but Lipschitz perturbation $\varphi(t)$ (see Fig. 5). Furthermore, the designed control strategy is able to stabilize the norm of the error despite the abrupt change between references, as seen in Fig. 6. The boundary control $q$ shown in Fig. 7 clarifies how such disturbance compensation is performed. This is due to the presence of the discontinuous term on the control design. Nevertheless, the control signal generated is continuous throughout the tracking task.

V. DESIGN DISCUSSION

Attractive design features of the proposed controllers
and their properties, revealed in simulations, are further

<figure><img src="image_6.png" /><figcaption>Fig. 5. Tracking of the output $y = p(0, t)$ over the successive references $r_p(t)$ versus time (sec).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Norm of the error versus time (sec).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Control signal $q$ versus time (sec).</figcaption></figure>

discussed.

• An eISS result is obtained instead of asymptotic stability due to the presence of the unmatched disturbance $\frac{D-x}{D}\ddot{r}(t)$ in (4). Nevertheless, the ultimate bound can be reduced either incrementing the control gains $\lambda_1, \lambda_2$ or choosing a slow tracking reference $r(t)$ making $L_{\dot{r}}$ smaller.

* Both proposed boundary controllers (5) and (12) require the only information of the wave equation at $x = 0$, i.e., $[u(0, t), u_t(0, t)]$, or the boundary diffusion equation information $p(D_2, t)$. This greatly minimizes the amount of information needed to perform the tracking task in both PDE's.

* Due to the exponential decay of the error dynamics (13) and by virtue of $\delta = \nu + \varphi$, the virtual input $\nu$ of the control (12) approaches the negative disturbance value, i.e., $\nu(t) \to -\varphi(t)$ as $t \to \infty$.

* Given the overall cascade system (1), the tracking for the diffusion equation calls for future investigation with the motion planing $r_p(t) = p(0, t)$, where $p(0, t)$ takes the place of the designed tracking control for the wave equation (5). Remarkably, whenever condition (9) is to be fulfilled, the
only parameter choice in the proposed control (5) is $a = 1$ (the linear case).

## VI. CONCLUSIONS

The present work summarizes recent advances in earthquake control. A cascade system, composed of a 1D wave equation coupled to a 1D diffusion equation, is presented as a simplified earthquake model. In order to ensure slow aseismic response of the system, a tracking control is designed in two steps. It is shown that a homogeneous boundary control is able to drive the wave equation states close to a bounded reference signal, resulting in eISS of the closed-loop system, whereas a PI plus discontinuous integral term control presented global exponential stability of the error system, despite the presence of uncertainties and/or disturbances in the model. Simulations show that the discontinuous version of the homogeneous algorithm achieves the best precision and robustness for the price of using high-frequency control signals. In turn, for the diffusion equation, simulations show robust tracking for different kinds of references, using continuous control signals and therefore diminishing the chattering effect in the plant dynamics. The presented control strategies give a deeper insight on the control of the coupled system to be designed, resulting in performing the motion planing of the diffusion equation choosing the linear case of the homogeneous control of the wave equation. The further development of such a cascade strategy calls for future work, as well as the 3D case of the earthquake phenomenon does.

## ACKNOWLEDGEMENTS

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). Prof. Yury Orlov work has been partially supported by Atlanstic2020, a research program of Région Pays de la Loire in France, and by CONACYT grant A1-S-9270.

## 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] 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.

[6] ——, “Control instabilities and incite slow-slip in generalized Burridge-Knopoff models,” *arXiv:2008.03755*, 2020.

[7] B. d'Andréa Novel and J. M. Coron, “Exponential stabilization of an overhead crane with flexible cable via a backstepping approach,” *Automatica*, vol. 36, no. 4, p. 587–593, 2000.

[8] M. Krstic and A. Smyshlyaev, *Boundary Control of PDEs*. SIAM Advances in Design and Control, 2008.

[9] I. Karafyllis and M. Krstic, *Input-to-state stability for PDEs*. Cham, Switzerland: Springer International Publishing, 2019.

[10] E. Fridman and Y. Orlov, “An LMI approach to $H_\infty$ boundary control of semilinear parabolic and hyperbolic systems,” *Automatica*, vol. 45, pp. 2060–2066, 2009.

[11] T. Meurer and A. Kugi, “Tracking control design for a wave equation with dynamic boundary conditions modeling a piezoelectric stack actuator,” *International Journal of Robust and Nonlinear Control*, vol. 21, pp. 542–562, 2011.

[12] W. Guo and B. Z. Guo, “Performance output tracking for a wave equation subject to unmatched general boundary harmonic disturbance,” *Automatica*, vol. 68, pp. 194–202, 2016.

[13] B. Laroche, P. Martin, and P. Rouchon, “Motion planning for the heat equation,” *International Journal of Robust and Nonlinear Control*, vol. 10, pp. 629–643, 2000.

[14] A. F. Lynch and J. Rudolph, “Flatness-based boundary control of a class of quasilinear parabolic distributed parameter systems,” *International Journal of Control*, vol. 75, no. 15, pp. 1219–1230, 2002.

[15] T. Meurer and A. Kugi, “Tracking control for boundary controlled parabolic PDEs with varying parameters: Combining backstepping and differential flatness,” *Automatica*, vol. 45, p. 1182–1194, 2009.

[16] F.-F. Jin and B.-Z. Guo, “Performance boundary output tracking for one-dimensional heat equation with boundary unmatched disturbance,” *Automatica*, vol. 96, pp. 1–10, 2018.

[17] X.-H. Wu and H. Feng, “Output tracking for a 1-d heat equation with non-collocated configurations,” *Journal of the Franklin Institute*, vol. 357, p. 3299–3315, 2020.

[18] Y. Orlov and V. I. Utkin, “Sliding mode control in indefinite-dimensional systems,” *Automatica*, vol. 23, no. 6, p. 753–757, 1987.

[19] Y. Orlov, *Nonsmooth Lyapunov Analysis in Finite and Infinite Dimensions*. Cham, Switzerland: Springer International Publishing, 2020.

[20] Y. Orlov, A. Pisano, and E. Usai, “Exponential stabilization of the uncertain wave equation via distributed dynamic input extension,” *IEEE Transactions on Automatic Control*, vol. 56, no. 1, pp. 212–217, 2011.

[21] A. Pisano, Y. Orlov, and E. Usai, “Tracking control of the uncertain heat and wave equation via power-fractional and sliding-mode techniques,” *SIAM J. Control Optim.*, vol. 49, no. 2, pp. 363–382, 2011.

[22] A. Pisano and Y. Orlov, “Boundary second-order sliding-mode control of an uncertain heat process with unbounded matched perturbation,” *Automatica*, vol. 48, pp. 1768–1775, 2012.

[23] D. Gutiérrez-Oribio, Y. Orlov, I. Stefanou, and F. Plestan, “Tracking for a wave equation using homogeneous boundary control,” in *20th European Control Conference*, London, UK, 2022.

[24] ——, “Robust tracking for the diffusion equation using boundary control,” in *IEEE Conference on Decision and Control (Accepted)*, Cancún, México, 2022.

[25] M. A. Biot, “General theory of three-dimensional consolidation,” *Journal of Applied Physics*, vol. 12, no. 155, pp. 155–164, 1941.

[26] K. Terzaghi, *Theoretical Soil Mechanics*. John Wiley & Sons, Inc., 1943.

[27] 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.

[28] 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.

[29] Y. Orlov, Y. Aoustin, and C. Chevallereau, “Finite time stabilization of a perturbed double integrator—part I: Continuous sliding mode-based output feedback synthesis,” *IEEE Transactions on Automatic Control*, vol. 56, no. 3, pp. 614–618, 2011.

[30] S. Emelyanov, S. Korovin, and L. Levantovsky, “Higher-order sliding regimes in the binary control systems,” *Soviet Physics Doklady*, vol. 31, no. 4, pp. 291–293, 1986.

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

[32] S. Dashkovskiy and A. Mironchenko, “Input-to-state stability of infinite-dimensional control systems,” *Math. Control Signals Syst.*, vol. 25, pp. 1–35, 2013.

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