Applications of mathematics to maritime search

The issue of searching missing aircraft is valuable after the event of MH370. This paper provides a global optimal model to foster the efficiency of maritime search. Firstly, the limited scope, a circle whose center is the last known position of the aircraft, should be estimated based on the historical data recorded before the disappearance of the aircraft. And Bayes' theorem is applied to calculate the probability that the plane falling in the region can be found. Secondly, the drift of aircraft debris under the influence of wind and current is considered via Finite Volume Community Ocean Model(FVCOM) and Monte Carlo Method(MC), which make the theory more reasonable. Finally, a global optimal model about vessel and aircraft quantitative constraints is established, which fully considers factors including the area of sea region to be searched, the maximum speed, search capabilities, initial distance of the vessels by introducing 0-1 decision variables.

Abstract. The issue of searching missing aircraft is valuable after the event of MH370. This paper provides a global optimal model to foster the efficiency of maritime search. Firstly, the limited scope, a circle whose center is the last known position of the aircraft, should be estimated based on the historical data recorded before the disappearance of the aircraft. And Bayes' theorem is applied to calculate the probability that the plane falling in the region can be found. Secondly, the drift of aircraft debris under the influence of wind and current is considered via Finite Volume Community Ocean Model(FVCOM) and Monte Carlo Method(MC), which make the theory more reasonable. Finally, a global optimal model about vessel and aircraft quantitative constraints is established, which fully considers factors including the area of sea region to be searched, the maximum speed, search capabilities, initial distance of the vessels by introducing 0-1 decision variables.
1. According to the historical data before the disappearance of the aircraft, how to estimate the scope used to search the plane?
The first step is determine whether the search scope is on the land or in the sea. We take the last known position (LKP) of distress target [13] as the center of the circle, the most distance between the last known time and known or estimated distress time as the radius of the circle, which contains the location of all possible distress target. Then we conduct the distress inference that based on the assumption of known facts and over in one's mind, establishing searching reference (Datum). Finally, we calculate the position of total probable error to obtain the reasonable scope to search, determining the search area of the process. 2. In the search range, whether the probability of finding the missing aircraft in each region is the same? If no, how to determine the probability?
We should analyze the influence that wind and current have on the probability of finding the airplane on the sea range to determine the probability distribution based on the Bayesian theory combined with FVCOM [4] model (current model).Then we calculate the specific influence results by Monte Carlo method [15], which will be used to give a feedback on the initial probability distribution and update the real-time probability distribution. 3. How to search for a lost plane?
We abstract problem as a global optimization problem with the number constraints of searching ships and aircrafts to establish a three-dimensional maritime searching and global optimization model [19]. The model takes into account the searching area, the maximum speed, the ability of ships and aircrafts to search, the initial distance, the maximum endurance time and other factors of searching ships and aircraft, and comprehensively analyze the relation between the search coverage time and search power quantity (cost), then we obtain the optimization (economic and feasible). The probability distribution in question 2 will be updated immediately.
This paper contains six parts-Introduction, Determine the scope of the search area limits initially, Based on Bayesian probability distribution model, Drift model, Research of global optimization model and Conclusions, solving problems about maritime search. The first part introduces how to analyze the problem about air crash and the outline about this paper, besides the background information. Then the second is written about how to make use of the historical data of flight and the last known position to estimate the possible range of wreckage. The third is using Bayesian probability theory and distribution model to narrow the possible range. The ocean current is considered in the fourth part to compute the location of aircraft accurately, after which global optimization model is established to arrange rescue flight 2. Determine the scope of the search area limits initially. The first step of the plan is to determine the scope of the search, which is used to guide the search effectively in the actual search operation. Usually we take the last known position(LKP) of the distress target as the center of the circle, and the longest distance calculated from the distance data of the last known time and the data of known or estimated distress time as the radius to get the possible position limit range, shown as follow Figure 1.
Searching in such a large range limit in many cases is impractical. At possible area, the probability of target location probability distribution is the key factor in the search plan, because it affects how to use the available search facilities.
So we use Bayes' theorem [12] to estimate the range of search probability distribution.
3. Based on Bayesian probability distribution model.

3.1.
Bayes' theorem. Typically, there are differences between the probability of an event A under the conditions of the event B and the probability of the event B under the conditions of the event A, however, both of which have determined relationship, Bayes' theorem states this relations. Its contents are: B appear under the premise, the probability of A appear equal probability of appearing under the premise of A multiplied by the probability of A B appears divided by the probability that B occurs again arise. Through contact A and B, it is calculated from the probability of occurrence of an event of another event, which dates back from the original results.
Bayes' theorem is stated mathematically as the following equation: where A and B are events. P (A) and P (B) are the probabilities of A and B, P (A|B), the conditional probability, is the probability of A given that B is true.
Bayes' technical principles divided into prior distribution and posterior distribution.
Prior distribution is a probability distribution of the overall distribution parameter θ. Bayesian believed that in any statistical inference questions about θ, in addition to use the information provided in the sample X, you also need to set a prior distribution of θ, which is an indispensable factor during the inference.
Bayesian regarded the prior distribution as the probability of a priori information about θ before sampling. Prior distributions do not have an objective, and it may be partially or entirely based on the belief in charge.
Posterior distribution is the calculation of the conditional distribution (θ|x) of θ under the condition of known X = x on the premise of the distribution P (θ) of sample X and the prior distribution π(θ) of the θ, based on the theory of probability and conditional probability distribution. Because of this distribution is got after sampling, it is called the posterior distribution. Bayesian believed that: this distribution combined information that provided by sample X and prior distribution π(θ).

3.2.
Probability distribution model. Based on the Bayes' theorem, the whole area is divided into a number of small grid, each of which has two probability value p and q. p is the probability that the missing aircraft is in the grid, and q is the probability to find the aircraft assumed that it is in the grid, if it turns out to be unsuccessful in finding the missing airplane after searching, the probability should be updated to Besides, if the probability of its original missing airplane in its place is r, then this will be updated as a probability Besides, with more new evidence found, we could constantly adjust it may be the probability of location. 4. Drift model. After determining the site where flight crashed on the ocean, we need to consider how the ocean impacts the bodies of the flight. As the time passed by, the bodies of the flight will be drifted by wind, waves and current. By quantifying these uncertainties, we can determine the optimal search area. But there are too many uncertain factors in ocean, and in order to establish a simplified model, we can define some limited conditions: 1. In this model, we ignore the effect of waves. The effect of waves to the bodies of the flight is very complex and hard to evaluate due to the difficulties of data forecast and collection, and it is also negligible comparing with the effect of wind or current. 2. Regardless of the wind, we assume that the velocity of drift objects is equal to surface water flow movement. When the water flow velocity changes, the movement of objects also changes. 3. Consider of only the wind drift error, wind induced drift error and initial position error. Now we consider wind and current respectively on the influence of drift.

External forces analysis.
4.1.1. Wind-induced drift. The wind will impact the bodies of the flight that on the surface, making these parts move in the direction of the wind, and the ocean will drag the parts of the bodies that under water. These two forces make the bodies of the flight drift. Because of the irregularity of the bodies, there is an angle will be exist between the direction of the drift trajectory and the direction of the wind, According to a prior, there is a linear relationship between the drift velocity of the floats and the wind velocity that ten miles above the ocean. So we define the wind impact as the average of ten minutes wind speed that ten miles above the surface W 10 . The wind direction vector can be divided into two parts: downwind leeway component(DWL) and crosswind leeway component(CWL). And it is shown as above Figure 2: According to the different species floater of the experimental research in the wind-induced drift conducted by Allen et al., DWL and CWL are supposed to be considered in a linear relationship of wind speed. The relation expressions of downwind leeway component L d and crosswind leeway component L c are We use the code of the FVCOM from Internet. For instance, the ocean surface current data of the day that FA447 crashed is shown as Graph 3 as above.

4.2.
Drift model under Monte Carlo method. According to the conditions we give above and the impact of wind and current we have quantified [3], the drift model can be simply regarded as the sum of flow induced drift and wind induced drift:  We then use second order Runge-Kutta to calculate, where the flow field data is the simulated result from the FVCOM model.
We assume that the positions of the bodies are expressed as a Markov process.
That is to say, the probability density function of the future state is only depends on the current status, and has no relationship with special process of the model system that achieve the status. On account of this assumption, the drift trajectory is satisfied with dx = V (x, t)dt + dε (7) where x is the position of the bodies ; V is the velocity function that the bodies under the action of external forces;dε is a random perturbation terms that consists of zero mean and known variance, which represents the uncertainties of external forces. Considering the impact of the flow on the determination of the search area is negligible, we just take the wind induce drift as dε into account.

4.3.
The error analysis of drift model. For the drift direction of the bodies is average distribution, we assume that during the drift, whether the initial direction is left or right, the drift direction of the bodies is fixed. As mentioned above, there has a linear relationship between wind-induced drift and wind velocity, and because of the existence of accidental factors [7], the slope and the intercept of the linear relationship should be added a random perturbation term respectively that is consistent with normal distribution N (0, σ), a n = a + ε n /20 b n = b + ε n /2 When W 10 = 10 m/s, make the slope and intercept in the same disturbance, then L n = (a + ε n /20) W 10 + (b + ε n /2) When the wind velocity is low, the Disturbance is mainly induced by the intercept (b + ε/2), and when the wind velocity is high, it is induced by slope (a + ε n /20).
Assuming that the disturbance wind field satisfy with normal distribution The disturbance wind velocity can be expressed as Then the entity of the wind induced drift function is

5.
Research of global optimization model.

5.1.
Model of the proposed. The joint aeronautical and maritime search is the searching activity coordinated by aircraft and ships, also known as the most efficient means of finding distress targets in the sea. The problem focus on determining the number of searching planes and searching ships so they can work cooperatively to implement prompt and effective search in the delimited area within the shortest time. We establish the global optimization mathematical models that can help to make an alternative scheme to realize the optimal usage of maritime three-dimensional search power by solving these mathematical models.
We assume the circumstances that several available rescue ships and rescue airplanes are around search sea area. Besides, it's possible that there are several passing ships inside and outside the search sea area. The power of the joint aeronautical and maritime search comes from all the available ships and planes stated above. Among them, the primary distance, optimal navigational speed and search ability can be differentiated between specialized rescue ships and passing ships. These indexes can also be differentiated among specialized rescue ships. In practice, the difficulty is that the feasible and effective search scheme is supposed to be achieved, which is exactly the question this paper tries to solve by establishing following mathematical models. Assume , quantity ceiling of search airplanes is Q a (1 ≤ Q a ≤ N ); 9. The whole search activity takes T (h); Then: to reach search sea area at full speed; 2. i ship takesT v i = TT v i to search the whole search sea area; 3. j airplane takesT a j = 2D a j /V a j to go back and forth to the search sea area at full speed; 4. j airplane takes −T a j = T − T 1 j to search the whole search sea area; 5. Sorties of j airplane is L j = T /T 1 j . These mathematical models aim to determine the optimal quantity of search ships and search airplanes to work coordinately, minimizing the time consumed for searching the sea area on the premises of the following decision variables introduced [6]: Assuming t s as the beginning time and t e as the ending time, the total search time is T = t e − t s .
In every search, the time of the rescue aircraft is equal to the maximum time of the best use of fuel. The time T L j consists of three parts. 1. The time of the plane a j to the area to search is T a j1 ; 2. The time of the plane a j to search on the area of searching is T a j2 ; 3. The time of the plane a j to turn back T a j3 . This paper assumes the time those flights go forth to and back from the search area is equal regardless of the time of plane supplies. So Based on the above analysis, meet full coverage of the need to meet the follow formula: It also means: Based on above formula, we can get the follow model (P 1 ): The model is to optimize the time T for search with two conditions. The first is the number of the ship, the second is the number of the plane for searching.

5.2.
Analysis model for solving complex. After establishing the model, the most urgent problem is how to settle it out. The model (P 1 ) established in the paper only has limited decision-making scheme (feasible solution), in which its decisionmaking variables x i and y j can only be integer 0 or 1 to ensure the existence of optimal solution, and exhaustive method can be taken in the theory (Enumerate all the feasible solution, and then compare the objective function value of each feasible solution) to find the most optimal. However, the optimal solution of model (P 1 ) by obtained from exhaustive method is not feasible due to the following reasons.
Provided the total number of available rescuing power [14] is n (n is the size of the problem referred to S = 2 n ), in the model (P 1 ), n = M + N , and the total number of programs is S = 2 n without regard to the total amount of all constraints. If computer completes the calculation of a program at the time of 1ns (1ns= 10 −9 s), then exhaustive method should be used to calculate the required time in the scale of the problem in 30, 50 and 100 cases, as shown in Table 1.
If we take all the constrains into consideration, the result of the selection scheme for the sea power is the combinatorial number N Q v (M is number of ships, Q v is total of all constrains), and for the air power is N Q a (N is number of planes, Q a is total of all constrains), so the total number of three-dimensional search scheme is   Table 2.
It can be seen from data in Table 1 and 2 that exhaustive method in larger scale cannot get the globally optimal solution [10] in a reasonable time. In maritime searching practice, searching personnel can always process a great deal of searching power information through various technical means, at the same time, the applying exhaustive method to obtain the global optimal solution is not feasible. So, the effective algorithm of the model (P 1 ) meeting the need of maritime search and practice is a valuable question.

Model algorithm.
The global optimized maritime three-dimensional search models (P 1 ) established in this paper have special objective function (a linear objective function with decision variable of numerator and denominator only taking 0 or 1) and constraint conditions (only two linear constraint condition). Like knapsack problem in inter programing, these problems belonging to a special fractional programming problem, called as 'fractional knapsack problem' in document, can be solved through Dinkelbach algorithm [9] in polynomial time to meet real needs.
In order to illustrate model (P 1 ) characteristics and solutions more clearly, we can implement equivalence transformation to objective function as follows: Assume: . . , M . According to the known conditions: So the model(P 1 ) can be equivalently expressed as the following model (P 2 ): Through previous analysis about the categorization of optimizing model [11], we can deduce that model (P 2 ) belongs to 0-1 programming model in integer programming as well as nonlinear programming model because its objective function is nonlinear function. The special point of the model is that the coefficient of decision variable x i (i = 1, 2, . . . , M ) is always 1, called as 'restraint of aggregate' in this paper.
Model (P 2 )'s objective function is a fractional expression. The coefficient p i of decision variable x i in numerator is positive, coefficient q i of decision variable x i in denominator is non-negative number, and constant term q 0 is positive. Now we study the value range of p 0 and approach of determining decision variable y j : According to expression 1. The ability of the airplane a j search area A a j , there is A a j > 0; 2. The time for plane to go back and forth to the search sea area and flying base for airplane a j isT a j , there isT L j = 2D a j /V a j > 0; 3. Maximum duration of flight of airplane a j is T L j , there is T L j > 0. In order to maximize the objective function T of model (P 2 ), p 0 = N j=1ā j y j should also be maximized, which is have all y j = 1, j = 1, 2, . . . ,N . Considering the restrict of aggregate of airplane, value range of decision variable y j of model (P 1 ) is as follows: 1. For T L j ≤T a j , j ∈ {1, 2, . . . , N }, y j ≡ 0; 2. For T L j >T a j , j = 1, 2, . . . ,N ,, if j ≤ Q a j , y i = 1, if j > Q a j , y j = 0. After determining the value of all decision variable y j , p 0 is a constant term and p 0 > 0, and the model of (P 2 ) can be solved based on Dinkelbach algorithm.

6.
Conclusions. This global optimal model is economic and feasible on the basis of analysis results under different vessel and aircraft quantity considering the relationship between searching coverage time and number of search cost. The probability of the existence of accidental factors, an error term that in the slope and the intercept are all added in the models to guarantee that optimization of search effort selection can be judged directly and effectively.
In this paper, the wind-induced drift is divided into downwind leeway component and crosswind leeway component to qualify the impact of the wind to the bodies of airplane. As for the impact of the current, Finite Volume Community Ocean Model(FVCOM) is used to stimulate the filed. To develop the drift model, Monte Carlo Method is introduced to calculate the movement of the aircraft debris, figuring out the ultimate position of the aircraft debris. These ways completely boost the accuracy and rationality of results.
Although this paper provides the search for missing airplanes with the optimization algorithm, there is no way to test the model due to lacking of actual data. In the future, the searcher would find the plane through this method rapidly to raise the possibility of passengers survived, if there is a similar event happen and having tremendous data.