Application of particle swarm optimization in optimal placement of tsunami sensors

The tsunami sensor optimal placement problem

Consider the domain with water areas Ω, and the subduction zone P. Let D be the part of the domain Ω where the sensors can be stationed. A sample illustration is shown in Fig. 1. Our goal is to determine the optimal placement of L sensors for the earliest detection of tsunami waves, arising from any source point in P. We denote {p_j}^P_{j = 1} as the points located in the subduction zone P, where p_j = (x_j,y_j). Moreover, let {q_i}^L_{i = 1} denote the sensors to be placed in D, where q_i = (x_i,y_i). We also let the configuration Q = {q₁,…,q_L} of L sensors represent a potential solution to the problem of interest.

Assume a source point p_j ∈ P generates some perturbation. The wave caused by such event will move away from the source, arriving at some water point x ∈ D. We denote (1) $$\tau (p_j,x):=\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{v}\mathrm{e}\mathrm{l}\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}\mathrm{f}\mathrm{r}\mathrm{o}\mathrm{m}{p}_j\mathrm{t}\mathrm{o}x$$

To calculate the travel time τ in (1) from p_j to an arbitrary point x ∈ D, we solve the 2D nonlinear SWE by initiating the wave from p_j. The time it takes for the wave to reach x will be τ(p_j, x). The numerical solution of the 2D nonlinear SWE is discussed in the next section. The time it takes for Q to detect a tsunami from p_j is computed as $$t(p_j,\mathbf{Q})=\underset{1\le i\le L}{min}\tau (p_j,q_i)$$Here, we take the minimum travel time over all L sensors since we want detection by one sensor to be as early as possible. Next, we need to take into account all possible source points p_j to obtain the tsunami registration time T(Q) by configuration Q, that is, we set $$T(\mathbf{Q})=\underset{1\le j\le P}{max}t(p_j,\mathbf{Q})$$where the maximum of the earlier expression is taken over all source points p_j. Thus, we formulate the optimization problem as follows:

Given L sensors, find Q* = {q₁,…,q_L} that minimizes T(Q) over D, or equivalently, (2) $${\mathbf{Q}}^{\star }=\underset{\mathbf{Q}\in D}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}T(\mathbf{Q}).$$The minimization problem above is based on the formulation in Astrakova et al. (2009).

The shallow water equations

Shallow water equations are a set of hyperbolic partial differential equations that describe fluid flow (Sadaka, 2012). They are often used to describe geophysical flows (such as rivers, lakes, coastal areas, etc.), rainfall runoff, tsunami propagation, or even atmospheric flows, given the appropriate initial conditions and source terms (Backhaus, 1983; Cea & Bladé, 2015; Rahman & Mandokhail, 2018; Galewsky, Scott & Polvani, 2004; Fu et al., 2017; Liu et al., 2008, 1995). SWE are derived from the Navier–Stokes equations, with the main assumption that the horizontal length scale is much greater than the vertical length scale. The 2D nonlinear shallow water equations are given by $$h_t+(hu{)}_x+(hv{)}_y=0$$ $$(hu{)}_t+{\left(h{u}^2+{\displaystyle \frac{1}{2}g{h}^2}\right)}_x+(huv{)}_y=gh{z}_x$$ $$(hv{)}_t+(huv{)}_x+{\left(h{v}^2+{\displaystyle \frac{1}{2}g{h}^2}\right)}_y=gh{z}_y$$Here, h is the local water depth, u,v are the velocities in the x and y directions, respectively, g = 9.81 m/s² is the gravitational acceleration, and ∇z is the bed profile. An equivalent form of this system is given by (3) $$h_t+(hu{)}_x+(hv{)}_y=0$$ (4) $$u_t+u{u}_x+v{u}_y+g{\eta }_x=0$$ (5) $$v_t+u{v}_x+v{v}_y+g{\eta }_y=0$$derived by setting h = η + z where η is the surface elevation. Note that we can rewrite uu_x, vu_y, uv_x, and vv_y in terms of _uq = hu, _vq = hv, u, v, and h, that is (6) $$u{u}_x+v{u}_y={\displaystyle \frac{1}{h}((uqu{)}_x-u{q}_{x}u)+{\displaystyle \frac{1}{h}((vqu{)}_y-v{q}_{y}u)}}$$

(7) $$u{v}_x+v{v}_y={\displaystyle \frac{1}{h}((uqv{)}_x-u{q}_{x}v)+{\displaystyle \frac{1}{h}((vqv{)}_y-v{q}_{y}v).}}$$

A compilation of the analytical solutions of the shallow water equations (both 1D and 2D) for different cases are presented in Delestre et al. (2013). In this article, we solve the 2D nonlinear SWE numerically by implementing the finite volume method in a staggered grid with a momentum conservative scheme (Pudjaprasetya & Magdalena, 2014; Stelling & Duinmeijer, 2003), which will be discussed below. The solution of the SWE provides us the value of h at any mesh point in the domain in a given time. Hence, we can determine the time it takes for a disturbance to reach any mesh point in the domain for the first time (i.e., minimal time), which is the unknown in (1). We will use linear interpolation to determine wave travel time at any given point, using the three closest mesh points to this point.

The Finite Volume Method (FVM) is a method in solving partial differential equations (PDEs) by evaluating the conservative variables across the volume (LeVeque, 2002). The idea is to divide the domain into parts (see Fig. 2), which we will refer to as cells/control volumes, then integrate the equation over that volume. We then use Gauss’ Theorem to transform the volume integral into a surface integral. Evaluating this integral will give us the FVM discretization of a PDE.

Finite Volume Method has been extensively used in various fields such as heat and mass transfer, gas dynamics and fluid flow problems (Guedri et al., 2009, Lukáčová-Medvid’ová, Morton & Warnecke, 2002; Demirdžić & Perić, 1990). Its popularity as a discretization method stems from its high flexibility. The discretization is carried out directly in the physical domain without the need of transformation between the physical and computational system. Furthermore, it is easier to implement in unstructured meshes in comparison to other methods. Another important feature of FVM is that it preserves conservation, making it quite attractive when modeling problems where flux is of importance (Moukalled, Mangani & Darwish, 2016).

We will use a staggered grid for the finite volume method. In a staggered grid, the scalar variables (mass, pressure, density, etc.) are defined at the cell centers, while the velocity or momentum variables are defined at the center of the cell faces. This is different from a collocated grid where all variables are stored in the same position (Harlow & Welch, 1965). Figure 3 shows the 2D staggered grid.

The implementation of FVM on a staggered grid on SWE is described as follows. Consider a rectangular spatial domain (0, L_x) × (0, L_y) with hard wall boundary conditions, that is, u(0, y, t) = u(L_x, y, t) = 0 and v(x, 0, t) = v(x, L_y, t) = 0. The domain is meshed with a grid of M × N cells. The center of each cell is denoted by x_i,j, while the center of the bottom edge, top edge, left edge, and right edge of each cell are denoted by $x_{i,j-{\displaystyle \frac{1}{2}}}$, $x_{i,j+{\displaystyle \frac{1}{2}}}$, $x_{i-{\displaystyle \frac{1}{2},j}}$, and $x_{i+{\displaystyle \frac{1}{2},j}}$, respectively. Note that the sizes of h, u, and v are M × N, M × (N + 1), and (M + 1) × N, respectively.

The approximation of the mass Eq. (3) at the cell centered at x_i,j is $${\displaystyle \frac{d{h}_{i,j}}{dt}+{\displaystyle \frac{\ast h_{i+{\displaystyle \frac{1}{2},j}}{u}_{i+{\displaystyle \frac{1}{2},j}}-\ast h_{i-{\displaystyle \frac{1}{2},j}}{u}_{i-{\displaystyle \frac{1}{2},j}}}{\mathrm{\Delta }x}+{\displaystyle \frac{\ast h_{i,j+{\displaystyle \frac{1}{2}}}{v}_{i,j+{\displaystyle \frac{1}{2}}}-\ast h_{i,j-{\displaystyle \frac{1}{2}}}{v}_{i,j-{\displaystyle \frac{1}{2}}}}{\mathrm{\Delta }y}=0}}}$$

with upwind approximations $$\ast h_{i+{\displaystyle \frac{1}{2},j}}=\{\begin{array}{cc}{h}_{i,j},\hfill & if{u}_{i+{\displaystyle \frac{1}{2},j}}\ge 0\hfill \\ h_{i+1,j},\hfill & if{u}_{i+{\displaystyle \frac{1}{2},j}}<0,\hfill \\ \hfill & \hfill \end{array}\hspace{1em}\ast h_{i,j+{\displaystyle \frac{1}{2}}}=\{\begin{array}{cc}{h}_{i,j},\hfill & if{v}_{i,j+{\displaystyle \frac{1}{2}}}\ge 0\hfill \\ h_{i,j+1},\hfill & if{v}_{i,j+{\displaystyle \frac{1}{2}}}<0.\hfill \end{array}$$The approximation of the momentum Eq. (4) is implemented at the cell centered at $x_{i+{\displaystyle \frac{1}{2},j}}$, and (5) at the cell centered at $x_{i,j+{\displaystyle \frac{1}{2}}}$. We use the relations (6) and (7) for the momentum conservative approximation of the advection terms. For positive flow directions u > 0, v > 0, we have $${\displaystyle \frac{d{u}_{i+{\displaystyle \frac{1}{2},j}}}{dt}+{\displaystyle \frac{u{\overline{q}}_{i,j}}{{\overline{h}}_{i+{\displaystyle \frac{1}{2},j}}}\left({\displaystyle \frac{u_{i+{\displaystyle \frac{1}{2},j}}-u_{i-{\displaystyle \frac{1}{2},j}}}{\mathrm{\Delta }x}}\right)+{\displaystyle \frac{v{\overline{q}}_{i,j-{\displaystyle \frac{1}{2}}}}{{\overline{h}}_{i+{\displaystyle \frac{1}{2},j}}}\left({\displaystyle \frac{u_{i+{\displaystyle \frac{1}{2},j}}-u_{i+{\displaystyle \frac{1}{2},j-1}}}{\mathrm{\Delta }y}}\right)+g{\displaystyle \frac{{\eta }_{i+1,j}-{\eta }_{i,j}}{\mathrm{\Delta }x}=0}}}}$$ $${\displaystyle \frac{d{v}_{i,j+{\displaystyle \frac{1}{2}}}}{dt}+{\displaystyle \frac{u{\overline{q}}_{i-{\displaystyle \frac{1}{2},j}}}{{\overline{h}}_{i,j+{\displaystyle \frac{1}{2}}}}\left({\displaystyle \frac{v_{i,j+{\displaystyle \frac{1}{2}}}-v_{i-1,j+{\displaystyle \frac{1}{2}}}}{\mathrm{\Delta }x}}\right)+{\displaystyle \frac{v{\overline{q}}_{i,j}}{{\overline{h}}_{i,j+{\displaystyle \frac{1}{2}}}}\left({\displaystyle \frac{v_{i,j+{\displaystyle \frac{1}{2}}}-v_{i,j-{\displaystyle \frac{1}{2}}}}{\mathrm{\Delta }y}}\right)+g{\displaystyle \frac{{\eta }_{i,j+1}-{\eta }_{i,j}}{\mathrm{\Delta }y}=0,}}}}$$where $${\overline{h}}_{i+{\displaystyle \frac{1}{2},j}}={\displaystyle \frac{1}{2}(h_{i+1,j}+h_{i,j}),\hspace{1em}\hspace{1em}\hspace{1em}{\overline{h}}_{i,j+{\displaystyle \frac{1}{2}}}={\displaystyle \frac{1}{2}(h_{i,j+1}+h_{i,j}),}}$$ $$u{\overline{q}}_{i,j}={\displaystyle \frac{1}{2}\left(u{\overline{q}}_{i+{\displaystyle \frac{1}{2},j}}+u{\overline{q}}_{i-{\displaystyle \frac{1}{2},j}}\right),\hspace{1em}u{\overline{q}}_{i+{\displaystyle \frac{1}{2},j}}=\ast h_{i+{\displaystyle \frac{1}{2},j}}{u}_{i+{\displaystyle \frac{1}{2},j}},}$$

$$v{\overline{q}}_{i,j}={\displaystyle \frac{1}{2}\left(v{\overline{q}}_{i,j+{\displaystyle \frac{1}{2}}}+v{\overline{q}}_{i,j-{\displaystyle \frac{1}{2}}}\right),\hspace{1em}v{\overline{q}}_{i,j+{\displaystyle \frac{1}{2}}}=\ast h_{i,j+{\displaystyle \frac{1}{2}}}{v}_{i,j+{\displaystyle \frac{1}{2}}},}$$ $$v{\overline{q}}_{i,j-{\displaystyle \frac{1}{2}}}={\overline{h}}_{i,j-{\displaystyle \frac{1}{2}}}{\overline{v}}_{i,j-{\displaystyle \frac{1}{2}}},\hspace{1em}{\overline{v}}_{i,j-{\displaystyle \frac{1}{2}}}={\displaystyle \frac{1}{2}\left(v_{i+1,j-{\displaystyle \frac{1}{2}}}+v_{i-1,j-{\displaystyle \frac{1}{2}}}\right)}$$ $$u{\overline{q}}_{i-{\displaystyle \frac{1}{2},j}}={\overline{h}}_{i-{\displaystyle \frac{1}{2},j}}{\overline{u}}_{i-{\displaystyle \frac{1}{2},j}},\hspace{1em}{\overline{u}}_{i-{\displaystyle \frac{1}{2},j}}={\displaystyle \frac{1}{2}\left(u_{i-{\displaystyle \frac{1}{2},j+1}}+u_{i-{\displaystyle \frac{1}{2},j-1}}\right)}$$

Explicit schemes require a careful selection of the time step to fulfill the stability of the equation. One of the observations of Courant, Friedrichs & Lewy (1967) is that in order for solutions of a difference equation to converge to the solution of a partial differential equation, the difference scheme should use all the information contained in the initial data that influence the solution. This condition is called the CFL condition. The CFL number (also known as the Courant number) is given by $$c={\displaystyle \frac{\mathrm{\Delta }t\mathrm{m}\mathrm{a}\mathrm{x}({\lambda }_{\mathrm{x}},{\lambda }_{\mathrm{y}})}{min(\mathrm{\Delta }x,\mathrm{\Delta }y)}}$$To satisfy the CFL condition, the grid speed, given by ${\displaystyle \frac{\mathrm{\Delta }x}{\mathrm{\Delta }t}}$ and ${\displaystyle \frac{\mathrm{\Delta }y}{\mathrm{\Delta }t}}$, must be at least as large as the propagation speed in x and y, given by λ_x and λ_y, respectively (Toro, 1999; Lax, 2013), that is, 0 ≤ c ≤ 1 holds. Thus, for the stability of SWE using FVM scheme on a staggered grid (Gunawan, 2015), we define the time step as follows: $$\mathrm{\Delta }t=c{\displaystyle \frac{min(\mathrm{\Delta }x,\mathrm{\Delta }y)}{\underset{i,j}{max}\left({\displaystyle \frac{1}{2h_{i,j}}\left|u{q}_{i+{\displaystyle \frac{1}{2},j}}+u{q}_{i-{\displaystyle \frac{1}{2},j}}\right|+\sqrt{g{h}_{i,j}},{\displaystyle \frac{1}{2h_{i,j}}\left|v{q}_{i,j+{\displaystyle \frac{1}{2}}}+v{q}_{i,j-{\displaystyle \frac{1}{2}}}\right|+\sqrt{g{h}_{i,j}}}}\right)}.}$$The above numerical scheme has been validated through several benchmark tests in 1D and 2D (e.g., dam breaks, transcritical flows, wave shoaling), where the obtained numerical results are in very good agreement with the analytical solutions (Pudjaprasetya & Magdalena, 2014; Magdalena, Iryanto & Reeve, 2018; Magdalena, Erwina & Pudjaprasetya, 2015; Magdalena & Pudjaprasetya, 2014; Magdalena, Pudjaprasetya & Wiryanto, 2014).

The particle swarm optimization algorithm

It was illustrated in Ferrolino, Lope & Mendoza (2020) that (2) is neither convex nor differentiable. Thus, gradient-based and local search algorithms may not be suitable to solve the problem. This justifies the use of population-based algorithms in determining the optimal locations.

For this work, we explore the use of the PSO algorithm (Kennedy & Eberhart, 1995; Shi & Eberhart, 1998) to solve the optimization problem presented in (2). This algorithm has been used to solve several optimization problems because of its speed and efficacy. It has been successfully used to solve unconstrained and constrained problems in a variety of fields such as engineering, communications theory, and operations research (Zhang, Wang & Ji, 2015). Other applications are in solving problems in clustering analysis (Chen & Ye, 2004), scheduling (Yu, Yuan & Wang, 2007; Zhang et al., 2014), facility location (Hajipour & Pasandideh, 2012) and vehicle routing (Belmecheri et al., 2013).

Particle Swarm Optimization is a population-based algorithm having a group of individuals (known as particles) moving in a search space to find the best solutions. The new positions of the particles are determined by their velocities, which are updated based on their own experience and the experience of neighboring particles. Through this, PSO can balance exploitation and exploration (He & Wang, 2007).

Algorithm 1 summarizes the steps of PSO, with some alterations proposed in Mezura-Montes & Coello (2011); Pedersen (2010). We use the MATLAB built-in command particleswarm in our implementations.

Semicircle subduction zone with flat bottom

Consider the same domain shown in Fig. 1. Here, D is the water surface with constant depth and the rest is dry land. The optimal location of L = 1, …, 4 sensors are shown in Fig. 4. Note that for the case when L = 1, the optimal location of the sensor converged to the center of the circle. Moreover, as we employ more sensors, they tend to move closer to the subduction zone, and mimic its shape.

Now, we study how tsunami detection time will be influenced by changes in the number of sensors. Figure 5 plots the relationship between these two variables. We can see that detection time decreases as the number of sensors increases. Moreover, we can see from Table 1 that when we allocate additional sensors from 4, there is no significant difference compared to other cases. Hence, L = 4 sensors may be adequate in giving us fast time registration, without introducing too much costs.

Semicircle subduction zone with inverse tangent bottom

In this subsection, we examine how a different bathymetry will affect the optimal location of tsunami sensors. Specifically, we wish to investigate how the locations will change if there will be variations in the water depth. For this purpose, we consider an inverse tangent bottom topography as shown in Fig. 6. The optimal location of L = 1 and 2 sensors acquired using the inverse tangent bottom in comparison to the optimal locations attained using a flat bottom are presented in Fig. 7.

Observe in Fig. 7A that the sensor moves towards the shallower area. This phenomenon happens since waves travel slower in shallower areas, caused by the force exerted on them by the seabed. A similar trend can be observed for the two sensors as shown in Fig. 7B. Observe that the sensor in the right moves upwards and towards the shallower area. However, this will make tsunami detection time from the source points near the center longer. Thus, to accommodate these source points, the sensor on the left move rightwards and downwards. Note that this will not greatly affect tsunami detection time from the source points in deeper areas since waves travel much faster there.

Optimal placement of sensors in Cotabato trench

The Cotabato trench is considered as one of the most dangerous major fault zones in the Philippines since it has high tsunamigenic potential. It is located near a southwestern coast of Mindanao, Philippines. Some examples of tsunamigenic earthquakes that occured in this trench are the 1918 Celebes Sea earthquake (M 8.0), the 1976 Moro Gulf earthquake (M 8.1), and the 2002 Mindanao earthquake (M 7.5) (Stewart & Cohn, 1979; Bautista et al., 2012).

Consider a portion of the Cotabato Trench as shown in Fig. 8. The top left corner has latitude 6.3142° and longitude 124.1083°, while the bottom right corner has latitude 5.6158° and longitude 124.8067°. The subduction zone P is illustrated by the red line, and D is the water surface above the subduction zone. The corresponding bathymetric profile of this trench is shown in Fig. 9. We obtained the bathymetric profile of Cotabato trench from the General Bathymetic Charts of the Oceans (GEBCO), https://www.gebco.net/.

The placements of L = 1,…,6 sensor/s with minimum guaranteed time are illustrated in Fig. 10. Note here that as we employ more sensors, they become distributed along the subduction zone. The plot describing the relationship between detection time and the sensor’s number is shown in Fig. 11, while their corresponding values are presented in Table 2. Observe that increasing the number of sensors from 6 only shows little improvement in detection time. So, L = 6 sensors may be adequate in giving us fast time registration, without introducing additional costs.