Application of particle swarm optimization in optimal placement of tsunami sensors
 Published
 Accepted
 Received
 Academic Editor
 Marcin Woźniak
 Subject Areas
 Optimization Theory and Computation, Scientific Computing and Simulation
 Keywords
 Particle swarm optimization, Nonlinear shallow water equations, Tsunami sensors, Tsunami early warning system, Heuristic algorithm, Finite volume method
 Copyright
 © 2020 Ferrolino et al.
 Licence
 This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, reproduction and adaptation in any medium and for any purpose provided that it is properly attributed. For attribution, the original author(s), title, publication source (PeerJ Computer Science) and either DOI or URL of the article must be cited.
 Cite this article
 2020. Application of particle swarm optimization in optimal placement of tsunami sensors. PeerJ Computer Science 6:e333 https://doi.org/10.7717/peerjcs.333
Abstract
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 methods focused on these in literature (Goda, Mori & Yasuda, 2019; Park et al., 2019; Goda & Risi, 2017). In this research, our goal is to determine the optimal locations of sensors for the earliest detection of tsunami waves. Through this, we can enable early and effective response to reduce casualties.
Traditionally, tide gauges, typically placed in a pier, are used for tsunami detection. They are employed to determine trends in the mean sea level, tidal computation, and harbor operations and navigation. However, since they are located in areas where tsunamis still have very high energy, they are frequently destroyed and fail to report the incoming tsunami and its characteristics. Furthermore, as they are near the coast, tsunami confirmation may arrive too late for timely evacuation measures which may lead to an ineffective local response (Barrick, 1979). Hence, we instead consider seafloormounted sensors, particularly the bottom pressure recorders (BPRs). These sensors are capable of detecting tsunamis, even in open oceans, by measuring variations in hydrostatic bottom pressure. They use an acoustic coupling to transmit data from the seafloor to the surface buoys, which then relay via satellite the information to a landbased station (Eble & Gonzalez, 1990). Further details on these tsunami sensor networks can be found in Rabinovich & Eblé (2015) and Nagai et al. (2007).
Studies on where to optimally place sensors of tsunami observation networks are limited. Some considered various technical aspects supported by expert judgments such as financial limitations (Araki et al., 2008) and legal aspects on geographical boundaries (Abe & Iwamura, 2013). Attempts to maximize the accuracy of prediction of some tsunami parameters (e.g., source, amplitude, etc.) were proposed in Mulia, Gusman & Satake (2017), Meza, Catalán & Tsushima (2020), Mulia et al. (2017), Saunders & Haase (2018). There are also researches that incorporated the effectiveness of tsunami warnings. In Omira et al. (2009), Schindelé, Loevenbruck & Hébert (2008), possible locations of sensors while considering installation constraints were investigated. However, the researchers did not apply optimization algorithms to solve the problem. In Braddock & Carmody (2001) and Groen, Botten & Blazek (2010), the optimal location of sensors among the candidate detection sites were proposed. In their investigations, they did not take into account bathymetric data since a fixed wave travel speed was used.
Astrakova et al. (2009) developed a problem of placing sensors optimally to detect tsunami waves as early as possible. However, the algorithm of computation for the wave travel time (from a source point to a point of interest) is based on the Huygen’s principle, which is computationally expensive. In Ferrolino, Lope & Mendoza (2020), wave velocity approximation was based on wave front kinematics. In this article, we compute for the wave travel time using the 2D nonlinear shallow water equations (SWE). This approach allows us to more accurately compute the tsunami detection time.
The 2D nonlinear SWE is generally used in the numerical simulation of tsunami propagation from the open ocean to the coast (Ulutas, 2012). They are accurate for solving long wave propagation (i.e., waves whose vertical length scale is much greater than their horizontal length scale), and runup and inundation problems due to the introduction of nonlinear terms, which are essential in tsunami modeling. Moreover, 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 Pudjaprasetya & Magdalena (2014), Stelling & Duinmeijer (2003), Magdalena, Rif’atin & Reeve (2020), Magdalena, Iryanto & Reeve (2018), Andadari & Magdalena (2019), Magdalena, Pudjaprasetya & Wiryanto (2014), and Magdalena & Rif’atin (2020). 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 populationbased 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 (1) $$\tau ({p}_{j},x):=\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{v}\mathrm{e}\mathrm{l}\phantom{\rule{thickmathspace}{0ex}}\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}\phantom{\rule{thickmathspace}{0ex}}\mathrm{f}\mathrm{r}\mathrm{o}\mathrm{m}\phantom{\rule{thickmathspace}{0ex}}{p}_{j\phantom{\rule{thickmathspace}{0ex}}}\mathrm{t}\mathrm{o}\phantom{\rule{thickmathspace}{0ex}}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}\phantom{\rule{thickmathspace}{0ex}}\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}\phantom{\rule{thickmathspace}{0ex}}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_{1},…,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}\phantom{\rule{1pt}{0ex}}\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^{2} 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 _{u}q = hu, _{v}q = 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 $$\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 & \phantom{\rule{1em}{0ex}}if{u}_{i+{\displaystyle \frac{1}{2},j}}\ge 0\hfill \\ {h}_{i+1,j},\hfill & \phantom{\rule{1em}{0ex}}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 & \phantom{\rule{1em}{0ex}}if{v}_{i,j+{\displaystyle \frac{1}{2}}}\ge 0\hfill \\ {h}_{i,j+1},\hfill & \phantom{\rule{1em}{0ex}}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 $$\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},j1}}}{\mathrm{\Delta}y}}\right)+g{\displaystyle \frac{{\eta}_{i+1,j}{\eta}_{i,j}}{\mathrm{\Delta}x}=0}}}$$ $$\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}_{i1,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}\phantom{\rule{thickmathspace}{0ex}}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}\phantom{\rule{thickmathspace}{0ex}}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}_{i1,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},j1}}\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\phantom{\rule{thickmathspace}{0ex}}\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 $\frac{\mathrm{\Delta}x}{\mathrm{\Delta}t}$ and $\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\phantom{\rule{1pt}{0ex}}{\displaystyle \frac{min(\mathrm{\Delta}x,\mathrm{\Delta}y)}{\underset{i,j}{max}\left({\displaystyle \frac{1}{2{h}_{i,j}}\leftu{q}_{i+{\displaystyle \frac{1}{2},j}}+u{q}_{i{\displaystyle \frac{1}{2},j}}\right+\sqrt{g{h}_{i,j}},\phantom{\rule{thickmathspace}{0ex}}{\displaystyle \frac{1}{2{h}_{i,j}}\leftv{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, gradientbased and local search algorithms may not be suitable to solve the problem. This justifies the use of populationbased 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 populationbased 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 MezuraMontes & Coello (2011); Pedersen (2010). We use the MATLAB builtin command particleswarm in our implementations.
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: 




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. 
6: Update N. If val = 0: 



If val = 1: 


7: Terminate if the stopping criterion has been satisfied. Else, go back to 3. 
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 realworld problem. Because PSO is heuristic, ten independent runs were 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.
L  Detection time (s)  Improvement in time due to an additional sensor 

1  3.84183  – 
2  2.51604  1.32579 
3  1.46286  1.05318 
4  0.88702  0.57584 
5  0.87654  0.01048 
6  0.87565  0.00089 
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.
L  Detection time (s)  Improvement in time due to an additional sensor 

1  224.9557  – 
2  90.0057  134.9500 
3  69.2824  20.7233 
4  57.2515  12.0309 
5  36.5071  20.7444 
6  24.1509  12.3562 
7  20.3382  3.8127 
8  18.0073  2.331 
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.