Application of particle swarm optimization in optimal placement of tsunami sensors

Rapid detection and early warning systems demonstrate crucial significance in tsunami risk reduction measures. So far, several tsunami observation networks have been deployed in tsunamigenic regions to issue effective local response. However, guidance on where to station these sensors are limited. In this article, we address the problem of determining the placement of tsunami sensors with the least possible tsunami detection time. We use the solutions of the 2D nonlinear shallow water equations to compute the wave travel time. The optimization problem is solved by implementing the particle swarm optimization algorithm. We apply our model to a simple test problem with varying depths. We also use our proposed method to determine the placement of sensors for early tsunami detection in Cotabato Trench, Philippines.


INTRODUCTION
While not the most prevalent among all natural disasters, tsunamis rank higher in scale compared to any others because of its destructive potential. Tsunamis are a series of ocean waves prompted by the displacement of a large volume of water. They can be generated by earthquakes, landslides, volcanic eruptions and even meteor impacts, although they mostly take place in subduction zones caused by underwater earthquakes. The sudden motion created by the said disturbances gives an enormous shove to the overlying water, which then travels away from the source.
Tsunamis usually have small initial amplitudes, so they can go unnoticed by sailors. However, they have far longer wavelengths unlike normal sea waves. Since the rate at which a wave loses its energy is inversely proportional to its wavelength, tsunami lose little energy as they propagate. Therefore, out in the depths of the ocean, tsunamis can travel at high speeds and cross great distances with only little energy loss (Liu, 2009). As they approach shallower waters, they slow down, and begin to grow in height. Once they crash ashore, they can cause widespread property damage, destruction of natural resources, and human injury or death.
Given the devastating potential of tsunamis, it is of great importance to develop techniques regarding risk assessment and damage mitigation. There are a number of they can take into account many types of tsunami generation (Zergani, Aziz & Viswanathan, 2015;LeVeque, George & Berger, 2011). In solving the 2D nonlinear SWE, we use the finite volume method in a staggered grid, with a momentum conservative scheme, as introduced in , Stelling & Duinmeijer (2003), Magdalena, Rif'atin & Reeve (2020), Magdalena, Iryanto & Reeve (2018), Andadari & Magdalena (2019), Magdalena, Pudjaprasetya & Wiryanto (2014), and . This method is very effective and robust because it is wellbalanced, conservative, and capable of handling complex bathymetry and topography.
We apply the Particle Swarm Optimization (PSO) in solving the optimal sensor location problem. This is a population-based algorithm inspired from the intelligent collective behavior of animal groups such as flocks of birds. Studies show that these animals are capable of sharing information among their group, and such ability grants them a survival advantage (Kennedy & Eberhart, 1995;Shi & Eberhart, 1998). PSO has gained much attention in solving optimization problems because it is easy to implement, computationally inexpensive, robust, and quick in convergence.
We first test our proposed approach on benchmark problems before using it to determine the optimal placement of sensors in the Cotabato trench, located in Mindanao, Philippines. The Philippines is vulnerable to tsunamis due to the presence of offshore faults and trenches (Valenzuela et al., 2020). To date, a total of 41 tidal waves classified as tsunamis have struck the country since 1589 (Bautista et al., 2012). Moreover, compared to neighboring countries, it has received less attention in tsunami research despite having regions with high seismic activity (Løvholt et al., 2012).
The remainder of this article is structured as follows: the methods implemented to solve the problem of optimal sensor placement are discussed in the "Methodology". "Numerical Results" obtained on different domains are then presented. Finally, we conclude our work and provide suggestions for future research.

METHODOLOGY
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 1 ,…,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 τðp j ; xÞ :¼ travel time from p j to x: (1) 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 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 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 1 ,…,q L } that minimizes T(Q) over D, or equivalently, Tð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., 2008Liu et al., , 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 Here, h is the local water depth, u,v are the velocities in the x and y directions, respectively, g = 9.81 m/s 2 is the gravitational acceleration, and ∇z is the bed profile. An equivalent form of this system is given by 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 u q = hu, v q = hv, u, v, and h, that is ðð v qvÞ 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 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, 2002Demirdž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À 1 2 , x i;jþ 1 2 , x iÀ 1 2 ; j , and x iþ 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 The approximation of the momentum Eq. (4) is implemented at the cell centered at x iþ 1 2 ; j , and (5) at the cell centered at x i;jþ 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 ; 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 To satisfy the CFL condition, the grid speed, given by Dx Dt and Dy Dt , 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:

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.

NUMERICAL RESULTS
We now present the results of our numerical simulations. The proposed method was first tested on a rectangular domain with a semicircular subduction zone, and then later applied to a real-world problem. Because PSO is heuristic, ten independent runs were Algorithm 1 Particle swarm optimization algorithm. Input: npop: swarm size, v: initial particle velocities within the range [−r,r], nvars: number of variables Output: best particle location 1: Initialization. Create particles at random of size npop within bounds. Also, initialize inertia ω and N to nbor = max(2,npop×frac), where nbor is the minimum neighborhood size and frac is the minimum neighbors fraction. Set the stall counter sc = 0.
2: Evaluate the cost function f for all particles. Set best = min (f(p i )) where p i is the position of particle i. In the latter iterations, p i will be the position of the best objective function that particle i has found. Let loc be the location such that best = f(loc).
3: For each particle i in the swarm at position x i , do the following: Select a subset S of N particles at random. Determine the best cost function value cfbest(S), and the position xbest(S) of cfbest(S).

Update the particle's velocity to
where v is the previous velocity, u 1 ,u 2 are uniformly (0,1) distributed random vectors of size nvars, and w 1 ,w 2 are the weights for exploitation and exploration, respectively.
Update the particle's position to x = x + v (implement thee bounds if necessary).

Calculate the cost function
4: Let F = min F(k) out of all k particles in the swarm. If F < best, set best = F and loc = x.

5:
If the best cost function value is lowered, set val = 0. Otherwise, val = 1.
Set N to nbor.
7: Terminate if the stopping criterion has been satisfied. Else, go back to 3.
implemented for each scenario that was considered, and we present here the average values of the results of these runs. Moreover, we set the population to 50 * L, where L is the number of sensors, and the number of function evaluations to 10,000.

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  (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.

CONCLUSION AND RECOMMENDATIONS
The problem of placing sensors optimally is considered to provide the earliest detection of tsunami waves, which leads to effective warning and response by local communities. A more accurate calculation of wave travel time is proposed, which is done by solving the 2D nonlinear shallow water equations. Moreover, the PSO algorithm was implemented to solve the minimization problem. The proposed model was first tested to simple domains with varying bathymetries, comprising of a semicircle subduction zone. The calculated sensor locations for the benchmark problems are geometrically sensible. Afterwards, the model is applied to calculate the optimal placement of tsunami sensors in Cotabato Trench. The actual bathymetry profile of this trench is used to acquire a more  realistic estimation of the wave travel time. We note here that our model can be used to any domain with any form of subduction zone or bathymetry. Moreover, it can consider other types of tsunami generation by changing the conditions in the 2D nonlinear shallow water equations.
Numerical results have shown an inverse relationship between the number of sensors and tsunami detection time. However, employing additional sensors implies more cost. Imposing constraints in cost and installation (e.g., depth, location) may be done in future works. Another recommendation is to consider maximizing tsunami warning efficacy and accuracy of estimation of tsunami parameters as additional objective functions.

ADDITIONAL INFORMATION AND DECLARATIONS Funding
This work was funded by the UP System Enhanced Creative Work and Research Grant (ECWRG-2019-2-11-R) and the Research Grant from Institut Teknologi Bandung. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.