A08026022 Liu Zongchen Li Bohao Yang Hongjun 1-19

Precise positioning of multiple rocket debris

Abstract

This paper starts with the establishment of an optimization model; proposes an algebraic model and an optimization model based on the sonic boom time received by monitoring equipment and the three-dimensional coordinate information of the equipment. It enables faster and more accurate positioning of the landing site of rocket debris. The model precisely locates the three-dimensional position and time of the sonic boom by minimizing the sum of absolute differences between predicted and observed sonic boom times. Additionally, the Particle Swarm Optimization (PSO) algorithm is introduced to solve this optimization problem, providing an effective method for quickly and accurately locating the landing site of rocket debris, and offering a reference for the safe recovery and reuse of rocket debris.

For Problem 1: First, latitude and longitude are converted into three-dimensional coordinates. Then, a system of nonlinear equations is established based on the Euclidean distance formula (distance formula between two points). The Particle Swarm Optimization (PSO) algorithm is used to solve the nonlinear equations. The global search capability and parallel computing features of the PSO algorithm help reduce the complexity of the algorithm, improve computational efficiency, and decrease the number of function iterations. Moreover, it enhances the accuracy of the results. The results show that at least five monitoring devices are required, the coordinates of the sonic boom of the debris are (36.0991, -6.0368, 2.000), and the initial time of the sonic boom is 1 s.

For Problem 2: Based on the model established in Problem 1, each pair of monitoring devices and the sonic boom location can form a triangle. According to the triangle inequality theorem, the sum of any two sides must be greater than the third side, and the difference must be less than the third side. By classifying and filtering out unsuitable data through this method iteratively, we can determine which sonic boom location the vibration waves received by the monitoring devices originated from. To determine the positions and times of the four debris when sonic booms occurred in the air, at least six monitoring devices are required.

For Problem 3: Based on the triangulation model established in Problem 2, input the given data into the data recognition model to obtain the classification results of the four sets of data corresponding to the first batch of six devices. Then, verify and identify the device data through testing, and finally use the model from Problem 1 to solve for the sonic boom locations and times.

For Problem 4: Error must be taken into account, which requires modifying the algorithm from Problem 3. The modified algorithm is used to identify data with errors, thereby solving for the corresponding positions and times. The obtained positions will differ from the actual positions, allowing the construction of a position correction model to determine the true positions and times.

Keywords: Euclidean distance, nonlinear equations, particle swarm optimization, triangle inequality theorem, data recognition model

1. Problem Restatement

1.1 Background Knowledge

Rocket debris precise positioning technology is a critical technique that plays a vital role in ensuring space launch safety, improving debris recovery efficiency, and protecting ground safety. In the modern aerospace field, rockets serve as key vehicles for space exploration. However, due to the high cost and complexity of rocket launches, coupled with the significant research value of rocket debris, the recovery and reuse of rocket debris have become crucial. Multi-stage rocket launches are one of the common methods, involving the separation and return of booster stages to Earth after mission completion. These fallen components may contain valuable materials and equipment. Accurately locating rocket debris aids in the rapid search and recovery of these valuable parts, thereby reducing the cost of space missions. Moreover, these fallen components may pose threats to ground facilities and safety, or even fall into natural environments such as oceans, land, or deserts, potentially impacting the ecological environment. Therefore, the precise positioning and swift recovery of such rocket debris are of utmost importance.

1.2 Problems to Be Solved

The vast majority of rockets are multi-stage rockets. After completing their designated tasks, the lower-stage rockets or boosters separate via interstage separation devices and fall. During their descent to the ground, the debris generates transonic sonic booms. To quickly recover rocket debris, multiple shockwave monitoring devices are deployed within the theoretical debris impact zone to capture the transonic sonic booms transmitted from different rocket debris in the air. Based on the arrival times of these sonic booms, the positions of the debris when the sonic booms occur are located, and ballistic extrapolation is then used to achieve rapid and precise positioning of the debris landing points.

Problem 1: By establishing a mathematical model, determine the position coordinates and start time of a single piece of debris in the air when a sonic boom occurs, as well as the number of monitoring devices required to detect the sonic booms.

Problem 2: By establishing a data identification model, analyze and determine which four rocket debris the four sets of shockwaves received by the monitoring devices originate from. Then, under the model from Problem 1, solve for the minimum number of monitoring devices required.

Problem 3: With an increase in monitoring points, based on Problem 2, use the enumeration method to provide classification results for all data. Then, run the model from Problem 2 separately to find applicable solutions, and finally determine the positions and times of the sonic booms generated by the four pieces of debris in the air.

Question 4: Considering measurement errors, reconstruct the model from Question 2, provide a corrected model example, establish an error interval for each monitoring point, analyze result deviations, and propose a solution for precise aerial localization of wreckage.

2. Problem Analysis

2.1 Analysis of Question 1

To determine the coordinates and time of a single wreckage's sonic boom, we first clarify the output type. The position can be treated as a 3D coordinate solution, while time should align with the observation system. By framing it as 3D coordinates, we can use the Euclidean distance formula (the standard coordinate distance between two points) to calculate the distance between the sonic boom location and monitoring equipment. With known speed of sound and unknown time variables, we can formulate equations. However, the degree of unknowns

> 1

leads to a system of nonlinear equations.

Subsequently, we need to solve this system of nonlinear equations to determine the location and time of the rocket debris' sonic boom. For solving nonlinear equations, what we seek are approximate solutions, and our goal is to make the results more precise, which requires more accurate initial values. To enhance precision, we decided to use the Particle Swarm Optimization (PSO) algorithm to solve the nonlinear equations. The quality of the solution is evaluated through fitness, making it easy to implement, highly accurate, and quick to converge. The optimal solution is found through iteration, where in each iteration, particles move based on two extremes: the first being the best solution found by the particle itself, and the second being the best solution found by the entire population. This process determines the sonic boom's location and time. Finally, to verify the accuracy of the solution, we use the remaining data for validation. By comparing the differences between the actual monitored sonic boom locations and those calculated by the model, we assess the model's accuracy and reliability.

2.2 Analysis of Problem Two

The objective of this problem is to classify the four sets of shockwaves and determine the locations and times of the four debris events. First, coordinate transformation is performed to ensure all device position data can be compared and calculated under the same coordinate system. Secondly, as stated in the problem, the data from each monitoring device can only be selected once. The sonic boom location of the debris in space and any two monitoring devices can form a triangle, which must adhere to the triangle rule. The enumeration method is used to identify all possible solutions, and as the number of monitoring points increases, the number of possible solutions also rises. The possible

Substitute the possible solutions into the inequality, i.e., the triangle inequality theorem (the sum of any two sides is greater than the third side, and the difference is less than the third side), to eliminate all non-conforming solutions, thereby reducing the range of possible solutions. Finally, we can use the nonlinear equation model established in Problem One to solve the remaining solutions, thereby obtaining the precise location of the monitoring points and the time of the sonic boom.

2.3 Analysis of Problem Three

Based on the data identification model established in Problem Two, the data from Problem Three is substituted. Using the model from Problem One, the enumeration method is applied to provide classification results for all data. As the number of monitoring points increases, the model from Problem Two is run separately to find applicable solutions. By continuously adjusting the number of monitoring devices, the accuracy of the results is improved. Through iterative solving, the classification results for each set of data are obtained. Subsequently, the nonlinear equation system established in Problem One is used to solve the classified data, determining the locations and times of the sonic booms corresponding to the rocket debris in each dataset.

2.4 Analysis of Problem 4

According to Problem 4, the device recording time has a random error of 0.5 s. Considering the impact of this error on the triangle method, we reconstruct the triangle method for constraint. Based on the monitoring device's accuracy and error range, an error interval is constructed for each monitoring point. These error intervals are then combined with the triangle theorem to build a system of inequalities. All inequalities must be satisfied; if any fails, the corresponding solution must be eliminated. After filtering out invalid solutions, the precise sonic boom location is determined using the data identification model from Problem 2.

3. Symbol Description

The main symbolic variables and their meanings used in this paper are shown in Table 3-1.

Symbol	Explanation says 勗
$χ_{B i}$	Represents the horizontal coordinate of monitoring point $B_{i}$
$y_{B i}$	Represents the vertical coordinate of monitoring point $B_{i}$ ;
$h_{i}$	Indicates the elevation of the $i$ th monitoring device
$n$	Indicates the number of monitoring devices
$B_{i} (x_{i}, y_{i}, h_{i})$	Indicates the position coordinates of the $^{a}$ th monitoring device
$V$	Indicates the propagation speed of the shock wave is $340 m / s$
$t_{i}$	Indicates the time of receiving signals from each station
$t$	Indicates the time elapsed since the rocket's sonic boom
$X = (x_{1}, x_{2}, \dots . x_{i}, \dots, x_{m})$	Represents a population composed of particles
$v_{i} = {(v_{i 1}, v_{i 2}, \dots . . v_{i n})}^{T}$	Represents the velocity of the particle population
$p_{i} = {(p_{i 1}, p_{i 2}, \dots . . p_{in})}^{T}$	Personal best value of particle $^{i}$
$B, C, D, E$	Decision coefficient matrix for 4 classes

Table 3-1

4. Model assumptions

(1) Ignore atmospheric turbulence and other complex atmospheric phenomena, and neglect friction between rocket debris and the atmosphere.
  (2) Assume the data error is negligible and approximates 0.
  (3) Assume all error values are approximately equal.
  (4) Assume the mass of rocket debris is equal or has no impact on experimental results.
(5) Assume the error in sensor data is known and can be processed using statistical methods.
  (6) Assume the error in collecting small signals is the same each time.

5. Model Establishment and Solution

5.1 Model Establishment and Solution for Problem 1

Let one of the detection points A be the coordinate origin, with the longitude direction as the x-axis and the latitude direction as the y-axis, thus constructing a plane rectangular coordinate system. If the longitude and latitude coordinates of the origin are

A (0, 0)

, and the longitude and latitude coordinates of other detection points are

B_{i} (x_{经度， B i}, y_{线度， B i})

, based on the approximate distance value of 111.263 km per degree of latitude and 97.304 km per degree of longitude, it can be concluded that

{\begin{cases} x_{B i} = 97.304 (x_{经度, B i} - x_{经度}) \\ y_{B i} = 111.263 (y_{经度, B i} - y_{线度}) \end{cases}

Where:

x_{B i}

represents the plane abscissa of detection point

B_{i}

;

y_{B i}

represents the plane longitudinal coordinate of detection point

B_{i}

;
There are

n

monitoring devices, and the position coordinate of the

i

th monitoring device is

B_{i} (x_{i}, y_{i}, h_{i})

. Here,

x_{i}

is the horizontal projection abscissa of the accuracy of the

i

th monitoring device,

y_{i}

is the horizontal projection ordinate of the

i

th monitoring device,

h_{i}

represents the elevation of the

i

th monitoring device. If

C

is the sonic boom position

(x, y, h)

, then the Euclidean distance can be obtained as follows:

d_{i} = \sqrt{{(x_{i} - x)}^{2} + {(y_{i} - y)}^{2} + {(h_{i} - h)}^{2}} (i = 1, 2, \dots, n)

While

d_{i} = V (t_{i} - t)

Where:

V

indicates the propagation speed of the shock wave is

340 m / s

t_{i}

represents the time of signal reception from each station.

t

represents the time elapsed since the rocket's sonic boom.
A system of nonlinear equations is established from the above equation.

{\begin{matrix} (\sqrt{{(x_{1} - x)}^{2} + {(y_{1} - y)}^{2} + {(h_{1} - h)}^{2}}) = V (t_{1} - t) \\ (\sqrt{{(x_{2} - x)}^{2} + {(y_{2} - y)}^{2} + {(h_{2} - h)}^{2}}) = V (t_{2} - t) \\ \dots \dots . \\ (\sqrt{{(x_{n} - x)}^{2} + {(y_{n} - y)}^{2} + {(h_{n} - h)}^{2}}) = V (t_{n} - t) \end{matrix}

Let

f_{i} (x, y, h, t) = | {({(x_{i} - x)}^{2} + {(y_{i} - y)}^{2} + {(h_{i} - h)}^{2})}^{\frac{1}{2}} - V (t_{i} - t) | (i = 1, 2 \dots . n)

There exists

f_{i} (x, y, h, t) = 0 (i = 1, 2 \dots n)

Min f = \sum_{i = 1}^{i = n} | {({(x_{i} - x)}^{2} + {(y_{i} - y)}^{2} + {(h_{i} - h)}^{2})}^{\frac{1}{2}} - V (t_{i} - t) |

By finding the extremum of this equation, the obtained result is an approximate solution to equation (5.1-4). This paper will employ the Particle Swarm Optimization algorithm to solve it.

Particle Swarm Optimization

In a

n

-dimensional search space, a population consisting of m particles is denoted as

X = (x_{1}, x_{2}, \dots . x_{i}, \dots, x_{m})

. Here, the position of the

i

-th particle is

x_{i} = {(x_{i 1}, x_{i 2}, \dots . x_{i}, \dots ., x_{i n})}^{T}

, and its velocity is

v_{i} = {(v_{i 1}, v_{i 2}, \dots . v_{i n})}^{T}

. The individual extremum of particle

i

p_{i} = {(p_{i 1}, p_{i 2}, \dots . p_{i n})}^{T}

, and the global extremum of the population is

p_{g} = {(p_{g_{1}}, p_{g_{2}}, \dots . . p_{g n})}^{T}

. During the search process, particles update their velocity and position by tracking two target values: one is the best solution found by the particle itself, known as the individual extremum; the other is the best solution found by the entire population, known as the global extremum. The iterative calculation formula is:

\begin{matrix} {\begin{matrix} v_{i d} = w v_{i d} + c_{1} ξ (p_{i d} - x_{i d}) + c_{2} η (p_{g d} - x_{i d}) \\ x_{i d} = x_{i d} + v_{i d} \end{matrix} \\ {\begin{matrix} v_{i d}^{k + 1} = w v_{i d}^{k} + c_{1} ξ (p_{i d}^{k} - x_{i d}^{k}) + c_{2} η (p_{g d}^{k} - x_{i d}^{k}) \\ x_{i d}^{k + 1} = x_{i d}^{k} + v_{i d}^{k + 1} \end{matrix} \end{matrix}

In the equation,

v_{i}^{k}

represents the velocity of particle

i

at the

k

-th iteration,

v_{i d}^{k}

is the velocity component of

v_{i}^{k}

in the

d

-th dimension;

x_{i}^{k}

is the position of particle

i

at the

k

-th iteration,

x_{i d}^{k}

is the position component of

x_{i}^{k}

in the

d

-th dimension;

p_{i}^{k}

is the individual extremum of particle

i

at the

k

-th iteration,

p_{i d}^{k}

is the individual extremum component of

p_{i}^{k}

in the

d

-th dimension;

p_{g}^{k}

is the global extremum of the particle swarm at the

k

-th iteration.

The global extremum,

p_{g d}^{k}

is the global extremum component of the

p_{g}^{k}

-th dimension in the

d

-th dimension;

w

is called the inertia weight,

c_{1}, c_{2}

is called the learning factor,

ξ, n

is a uniformly distributed random number within the interval

[0, 1]

, i.e.,

ξ, n \in U [0, 1]

Set the learning factor

c_{1}, c_{2}

, inertia weight

w

, and maximum evolution generation

K_{max}

. Initialize the current evolution generation to

k = 1

. Randomly generate the particle swarm

X (k)

within the defined domain, including the initial positions

x (k)

and velocities

v (k)

of the particles, and set their initial positions as the individual historical best positions, i.e., the individual extremum.

Construct the fitness function

f = \sum_{i = 1}^{i = n} | {({(x_{i} - x)}^{2} + {(y_{i} - y)}^{2} + {(h_{i} - h)}^{2})}^{\frac{1}{2}} - V (t_{i} - t) | (5.1 - 9)

.
Based on the current position of each particle, calculate the fitness value of the particles using the fitness function (5.1-9). Then, select the position of the particle with the highest fitness as the historical best position of the particle swarm, i.e., the global extremum.

Recalculate the fitness value of each particle's current position and compare it with its individual extremum. If the current position is better, update the particle's historical best position to its current position.

Compare each particle's personal best with the global best. If it is better, update the particle's personal best to the current global best, meaning the particle's current position becomes the historical best position of the swarm.

The following coordinates can then be substituted

Device	$x (km)$	$y (km)$	$h (km)$	Sonic boom arrival time (s)
A	0	0	0.824	100. 767
B	52. 446856	28.03752	0.727	112. 22
C	45.830184	64.64206	0.742	188.02
D	0.97304	69.09246	0.85	258.985
E	27.537032	45.95038	0.786	118.443
F	21.990704	79.77342	0.678	266.871

-18.876976

-9.23458

0.575

163.024

Table 1: Plane coordinates of each monitoring station
Use the particle swarm algorithm to solve the position and start time of the rocket sonic boom. The solution results are as follows.

Number of detection devices	$x (km)$	y (km)	$h (km)$	Start time of the sonic boom (s)	Target value
2	0.0014	0.0032	0.0010	1	56.464
3	0.0017	0.0016	0.0010	1	71.780
4	0.0037	0.0034	0.0010	1	90.049
5	21.9398	5.7285	2	1	78.906
6	33.5709	-4.7046	2	1	64.392
7	36.0991	-6.0368	2	1	31.436

Table 2: Solution results using different numbers of detection devices
At least 5 monitoring devices need to be deployed. The coordinates of the sonic boom's location are (36.0991, -6.0368, 2.000), and the starting time of the sonic boom is 1 s.

5.2 Model Establishment and Solution for Problem Two

(1) Determine the minimum number of monitoring devices required to locate the positions and times of the 4 wreckage occurrences
During the classification process, each data point can only be selected once. Let the time when n monitoring devices collect 4 sets of signals be A.

A = {[\begin{array}{cccc} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ ⋮ & ⋮ & ⋮ & ⋮ \\ a_{n 1} & a_{n 2} & a_{n 3} & a_{n 4} \end{array}]}_{n * 4}

On this basis, let the decision coefficient matrices for the 4 classes be B, C, D, E

B = {[b_{i j}]}_{4 * n}

(5.2-2), where

b_{i j} = {\begin{matrix} 0 \\ 1 \end{matrix}

b_{i j} = 1

indicates that the time set received by the device belongs to category B;

b_{i j} = 0

indicates that the time set received by the device does not belong to category B
Since there is only one set of data in each monitoring device that belongs to this category,

\sum_{i = 1}^{i = 4} b_{i j} = 1 (j = 1, 2, \dots, n)

It can be determined that the data in this category is

A_{1}

A_{1} = A \cdot B

Similarly,

C = {[c_{i j}]}_{4 \times n (5.2 - 6),} D = {[d_{i j}]}_{4 \times n (5.2 - 7)}, E = {[e_{i j}]}_{4 \times n (5.2 - 8)}

Among them,

^{C_{i j}} = {\begin{matrix} 0 \\ 1 & (5.2 - 9) \end{matrix} d_{i j} = {\begin{matrix} 0 \\ 1 & (5.2 - 10) \end{matrix}, e_{i j} = {\begin{matrix} 0 \\ 1 \end{matrix} (5.2 - 11)

\sum_{i = 1}^{4} c_{i j} = 1_{(5.2 - 12)}, \sum_{i = 1}^{4} d_{i j} = 1 (5.2 - 13), \sum_{i = 1}^{4} e_{i j} = 1

and all data can only be selected once, i.e., there exists

F

F = B + C + D + E (5.2 - 15)

where:

F = {[\begin{array}{lllll} 1 & 1 & 1 & \dots & 1 \\ 1 & 1 & 1 & \dots & 1 \\ 1 & 1 & 1 & \dots & 1 \\ 1 & 1 & 1 & \dots & 1 \end{array}]}_{4 * n} (5.2 - 16)

  Matrix multiplication forms a system of homogeneous linear equations
Theorem: When the number of equations is far less than the number of unknowns in the system, there will be infinitely many solutions
  So we need to impose constraints
  Using the triangle rule
Let the sonic boom position be

G

, the monitoring point position be

H, K

, forming

△ G H K

Point H received the sonic boom at time

t_{H}, K

, point received the sonic boom at time

t_{K}

, with speed v

| H G - H K | < K G < H K + H G (5.2 - 17)

Divide by v simultaneously, then the time for point

k

to experience the same sonic boom is

t_{k}

, yielding

| t_{H} - (\frac{H K}{V}) | < t_{K} < (\frac{H K}{V}) + t_{H}

(2) Classify the four sets of shock waves

First, generate all possible solutions using the enumeration method, set the initial value

n = 1

, and randomly classify the data of the nth monitoring device.

Then execute

n = n + 1

to calculate all permutations for the data of the

n

th monitoring device and the data of the

n

th monitoring device,

Obtain the solution set

M

.
By establishing inequality relationships between the first

n - 1

monitoring devices and the

n

-th monitoring device relative to the sonic boom location among

n

devices, eliminate non-compliant solutions from set

M

to derive solution set

M^{'}

Finally, as n increases, more constraints are introduced, leading to the optimal solution. Using the particle swarm algorithm from Problem 1, substitute

M^{'}

into (5.1–6) to solve for the coordinates of the sonic boom location and its occurrence time for each dataset.

Determine the minimum number of monitoring devices required.

The number of equations is

4 n + 4 n = 8 n (\sum_{i = 1}^{i = 4} b_{i j} = 1 (j = 1, 2, \dots, n)

\sum_{i = 1}^{i = 4} c_{i j} = 1 (j = 1, 2, \dots, n), \sum_{i = 1}^{i = 4} d_{i j} = 1 (j = 1, 2, \dots, n)

\sum_{i = 1}^{i = 4} e_{i j} = 1 (j = 1, 2, \dots, n)

F = {[\begin{array}{ccccc} 1 & 1 & 1 & \dots & 1 \\ 1 & 1 & 1 & \dots & 1 \\ 1 & 1 & 1 & \dots & 1 \\ 1 & 1 & 1 & \dots & 1 \end{array}]}_{4 n}

The inequality constraint is

4 C_{n}^{2}

units

\begin{aligned} 令^{x_{n}} & = 16 n - 8 n - 4 C_{n}^{2} \\ = 8 n - 4 C_{n}^{2} \end{aligned}

When

x_{n} > 0

, multiple solutions may exist
When

x_{n} < 0

, a unique solution may exist (Theorem: A unique solution is only possible when the number of unknowns is far fewer than the number of equations)
When

n \in N^{*}, n > 5

x_{n} < 0

, a unique solution may only exist if there are at least

n = 6

units

5.3 Model Establishment and Solution for Problem Three

First, based on the model from Problem Two, we can arbitrarily select data from 6 devices and sequentially input them into the recognition algorithm established in Problem 2. The classification results for the corresponding 4 sets of data from each device can then be obtained through solving.

Table 3: Classification Results of the 4 Sets of Data Corresponding to the First Batch of 6 Devices

Device	Group 1	Group 2	Group 3	Group 4
A	100.767	164.229	214.85	270.065

B	112.22	92.453	196.583	169.362
C	110.696	75.56	188.02	156.936
E	118.443	78.6	126.669	86.216
F	67.274	175.482	266.87	166.27
G	163.024	103.738	210.306	206.789

Subsequently, reselect 6 given data points and input them into the recognition algorithm established in Problem 2. Obtain the classification results of the corresponding 4 sets of data for each device through computation.

Table 4: Classification results of the 4 sets of data corresponding to the 6 devices in the second batch.

Device	Group 1	Group 2	Group 3	Group 4
A	100.767	164.229	214.85	270.065
B	112.22	92.453	196.583	169.362
C	110.696	75.56	188.02	156.936
D	94.653	141.409	196.517	258.985
F	67.274	175.482	266.87	166.27
G	163.024	103.738	210.306	206.789

The comparison shows that this data recognition algorithm can effectively identify device data with accurate results.

Finally, substitute the data from the above table into the model established in Problem 1, and use the particle swarm algorithm to solve for the location and time of the rocket's sonic boom.

Classification Group	$x (km)$	y (km)	$h (km)$	Start time of the sonic boom (s)
Group 1	18.216	69.957	4. 000	13. 323
Group 2	24.478	49.801	0.001	1
Group 3	12.637	48.772	0.001	66.661
Group 4	29.036	48.455	0.001	77.324

Table 5: Monitoring Equipment Data Solution Results

Figure 5.3 Particle Swarm Optimization (PSO) Algorithm Iteration Diagram

According to the iteration graph of the Particle Swarm Optimization (PSO) algorithm, the convergence speed is satisfactory throughout the iteration process. Within 500 iterations, the fitness functions corresponding to each group of data tend to stabilize. Within these 500 iterations, the PSO algorithm can perform multi-point searches in parallel, thereby expanding the search space and reducing the likelihood of local convergence. The parallelism of the PSO algorithm is manifested in two aspects: first, the inherent parallelism of the algorithm allows thousands of computers to independently perform population evolution calculations without affecting computational efficiency, enabling rapid global convergence; second, the embedded parallelism organizes searches through populations, allowing simultaneous exploration of multiple regions in the solution space, thus achieving higher efficiency with fewer computations.

5.4 Model Establishment and Solution for Problem Four

Due to various factors, there are errors between monitoring devices, and this error is resolved according to inequalities
The sonic boom arrival time at point

H

t_{H}, H

, and the actual time at point

t_{H}

is actual, then

t_{H 实际} - t_{H} \leq t_{H} \leq t_{H 实际} + Δ t_{H}

According to the triangle theorem:

| H G - H k | < k G < H k + H G

Let k receive the same sonic boom time as

^{t_{k} 实际}

| \frac{H k}{V} + t_{实际 H} | < t_{实际 k} < \frac{H k}{V} + t_{实际 H}

\begin{aligned} N_{1} = {t_{实际 k | \frac{H k}{v} - t_{H} | - max (Δ t) < t_{实际} < \frac{H k}{v} + t_{H} + max (Δ t)}}^{N_{2} = {t_{实际 k} | t_{k} - max (Δ t) \leq t_{实际 k} \leq t_{k} + max (Δ t) |}} (5.4-5) \end{aligned}

N_{1} \cap N_{2} \neq ϕ

, then possible solutions may exist
If

N_{1} \cap N_{2} = ϕ

, then no possible solutions exist
Due to a 0.5-second error, the model was used for identification based on the table provided in Question 3,
Table 6: Classification results of 4 data groups corresponding to 7 devices

Device	Group 1	Group 2	Group 3	Group 4
A	100.78	164.26	214.93	270.07
B	112.24	92.46	196.68	169.40
C	110.73	75.59	188.11	156.99
D	94.70	141.45	196.6	258.08
E	118.45	78.66	126.76	86.31
F	67.36	175.57	266.96	166.30
G	163.04	103.78	210.36	206.82

Select partial data from above to solve for the sonic boom location, using the model from Problem 1 to obtain the position and time of the sonic boom

Table 7: Location and time of sonic boom for each group of data

Classification group	Actual Position				Theoretical Position				Deviation (km)
Classification group	X (km)	Y (km)	H (km)	Sonic boom initiation time ( $s$ )	X (km)	Y (km)	H (km)	Sonic boom initiation time (s)	Deviation (km)
Group 1	18.216	69.957	4	13.323	18.192	69.915	0.001	13.780	4
Group 2	24.478	49.801	0.001	1	24.309	49.801	4	1	4
Group 3	12.637	48.772	0.001	66.661	18.614	38.203	0.001	91.929	12. 14
Group 4	29.036	48. 455	0. 001	77.324	27. 862	48.508	4	74. 200	4. 17

It is necessary to reduce the time error, otherwise the deviation value will be too large. Therefore, by adding the time error, a new Euclidean distance is obtained.
Formula

V (t_{i j 理论} - t_{i}) = \sqrt{{(x_{j} - x_{i})}^{2} + {(y_{j} - y_{i})}^{2} + {(h_{j} - h_{i})}^{2}}

(5.4–6)
Where:

(x_{j}, y_{j}, h_{j})

represents the coordinates of the

j

th monitoring device.

From this, the time error can be obtained:

Δ t_{i j} = t_{i j 理论} - t_{i j 实际} (i - 1, 2, 3 \dots n, j = 1, 2, 3 \dots m)

Based on this, the following algorithm can be proposed:
Given the initial value

A_{i j}^{n}, t_{i j 实际}^{n = 1)}

, the positions

(x_{j}, y_{j}, h_{j})

of each monitoring device, then use Model 1 to solve for the

i

th set of impact positions and time

(x_{i}^{n}, y_{i}^{n}, h_{i}^{n}), t_{i}^{n}

, determine the error. If

D = {d_{j} ∣ d_{j} > 1}

and

D \neq ϕ

, then the time when each monitoring point receives the information can be calculated.

\begin{aligned} t_{i j 理论}^{(n)} = t_{i}^{(n)} + \frac{1}{v} \sqrt{{(x_{j} - x_{i}^{(n)})}^{2} + {(y_{j} - y_{i}^{(n)})}^{2} + {(h_{j} - h_{i}^{(n)})}^{2}} \\ i = 1, 2, \dots ., N; j = 1, 2, \dots, M \end{aligned}

Calculate the average time error.

Δ t_{j}^{(n)} = \frac{1}{N} \sum_{j = 1}^{j = M} Δ t_{i j}^{(n)} (i = 1, 2, \dots, N)

Among them:

Δ t_{i j}^{(n)} = t_{i j}^{(n)} 理论 - t_{i j}^{(n)}

Update the monitoring point time for the

n + 1

th time

t_{i j 实际}^{(n + 1)} = t_{i j}^{(n)} + Δ t_{j}^{(n)} (i = 1, 2 \dots, N; j = 1, 2, \dots, M)

Update data

A_{i j}^{(n + 1)}

, then restart error judgment from the new data;

Table 8: Location and time of sonic booms for each group of data

Classification Group	Actual Position				Adjusted Theoretical Position				Deviation (km)
Classification Group	X (km)	Y (km)	H (km)	Sonic boom initiation time (s)	X (km)	Y (km)	H (km)	Sonic boom initiation time (s)	Deviation (km)
Group 1	18. 216	69.957	4	13. 323	18. 192	69.915	3.99	13.780	0.05
Group 2	24. 478	49. 801	0.001	1	24. 309	49. 801	0.9	1.000	0.91
Group 3	12.637	48.772	0.001	66. 661	12. 963	48.701	0.001	65. 900	0.33
Group 4	29.036	48.455	0.001	77.324	29.96	48.508	0.001	77. 200	0.93

6. Model Evaluation and Promotion

6.1 Advantages of the Model

(1) To reduce cumbersome calculations and iteration times, this model does not employ conventional methods typically used for nonlinear equation systems such as quasi-Newton or Newton methods, but instead utilizes the particle swarm optimization algorithm to solve the model.
(2) This model can efficiently process large amounts of data and quickly generate positioning results, reducing computational complexity.
(3) The constraints are easy to identify, as it only employs the triangle rule and eliminates unsuitable data through enumeration, without using overly complex methods.

6.2 Model Shortcomings and Improvement Directions

(1) We have ignored certain factors, which could actually have some impact on the results. If sensor data contains errors or is missing, the model's positioning results may be affected.
(2) The landing point of rocket debris is influenced by various factors such as wind speed, wind direction, terrain, etc. Therefore, these environmental factors must be fully considered during the modeling process, along with their potential dynamic changes over time. However, due to our team's limited capabilities, we are unable to achieve this.
(3) In complex environments, the accuracy of positioning may decrease. By integrating multiple data sources such as satellite remote sensing, ground radar, and drones, we can achieve the fusion and complementarity of multi-source information, thereby enhancing positioning accuracy.

6.3 Model Generalization

Not only can rocket debris be located, but all flying devices can also be positioned. It can also be used for GPS positioning.

7. References

[1] Jiang Qiyuan. Mathematical Model [M]. Beijing: Higher Education Press, 1987.
[2] Liang Jun, Cheng Can. Improved Particle Swarm Optimization Algorithm [J]. Computer Engineering and Design, 2008, 29(11): 2893-2896.
[3] Yang Zhaoxia, Fang Jianwen, Li Jiarong, et al. The Role of Particle Swarm Optimization Algorithm in Multi-parameter Fitting [J]. Journal of Zhejiang Normal University, 2008, 31(2): 173-177
[4] Xue Ting. Research and Improvement of Particle Swarm Optimization Algorithm [D]. Dalian: Dalian Maritime University, 2008.
[5] Du Yuping. Research on the Improvement of Particle Swarm Algorithm [D]. Xi'an: Northwest University, 2008.
[6] Chen Xiaoqian, Luo Shibin, Wang Zhenguo, et al. Research on Pre- and Post-processing in the Application of BP Neural Networks [J]. Systems Engineering - Theory & Practice, 2002, 22(1): 65-70.
[7] Jiang Baochuan, Hu Junming. Research on Improved Particle Swarm Algorithm for Solving Multimodal Functions [J]. Journal of Ningbo University, 2008, 21(2): 150-154.
[8] Zhang Xuanping, Du Yuping, Qin Guoqiang. An Adaptive Particle Swarm Optimization Algorithm with Dynamically Changing Inertia Weight [J]. Journal of Xi'an Jiaotong University, 2005, 39(10): 1039-1042.
[9] Feng Xiang, Chen Guolong, Guo Wenzhong. Setting and Experimental Analysis of Acceleration Factors in Particle Swarm Optimization Algorithm [J]. Journal of Jimei University, 2006, 11(2): 146-151.

Appendix

Problem 1 Program 1

%% 问题一数据转化
A= [
110.241 27.204 824 100.767
110.780 27.456 727 112.220
110.712 27.785 742 188.020
110.251 27.825 850 258.985
110.524 27.617 786 118.443
110.467 27.921 678 266.871
110.047 27.121 575 163.024];
A1=[110.241 27.204]; % 坐标原点所在的经纬度
A2=A(:,1:1:2) - [ones(size(A,1),1)*A1(1) ones(size(A,1),1)*A1(2)];
A3=[A2(:, 1)*97.304 A2 (:, 2)*111.260 A(:,3)/1000 A(:,4)]

Question 1 Program 2 Main Function

%% 基于粒子群算法求解
clc;clear all;close all;
A= [
110.241 27.204 824 100.767
110.780 27.456 727 112.220

110.712 27.785 742 188.020
110.251 27.825 850 258.985
110.524 27.617 786 118.443
110.467 27.921 678 266.871
110.047 27.121 575 163.024];
A1=[110.241 27.204]; % 坐标原点所在的经纬度
A2=A(:,1:1:2) - [ones(size(A,1),1)*A1(1) ones(size(A,1),1)*A1(2)];
A4=[A2(:,1)*97.304 A2(:,2)*111.260 A(:, 3)/1000 A(:,4)];
h= [ ]
W= [ ]
for k=2:1:7
% 构建数据集 A3
N1=A4(1:1:k,:);
N2=zeros(7-k,4);
A3 = [N1
    N2];
%粒子群参数设定
pop = 500;%种群数量
dim = 4;%变量维度 cl
ub = [500,500,2,100];%粒子上边界信息
lb = [-500 ,-500 ,0.001,1];%粒子下边界信息
vmax = [2,2,2,2];%粒子的速度上边界
vmin = [-2,-2,-2,-2]; %粒子的速度下边界
maxIter = 500;%最大迭代次数
fobj = @(x)fun(x,A3); %设置适应度函数为 fun(x);
%粒子群求解问题
[Best_Pos,Best_fitness,IterCurve] =
pso(pop,dim,ub,lb, fobj,vmax,vmin,maxIter);

%绘制迭代曲线
figure
plot(IterCurve,'r-','linewidth',1.5);
grid on;%网格开
title('粒子群迭代曲线')
xlabel('迭代次数')
ylabel('适应度值')
disp(['求解得到的 x,y,h,t 为:', num2str(Best_Pos(1)),'
',num2str(Best_Pos(2)),' ',num2str(Best_Pos(3)),'
',num2str(Best_Pos(4))]);
disp(['最优解对应的函数值为: ', num2str(Best_fitness)]);
h= [h
    Best_Pos];
W= [W
    Best_fitness];
end
W
h

Question 1 Particle Swarm Main Program

%%适应度函数
function fitness = fun(x,A3)
    x1 = x(1); %TS
    x2 = x(2); %Th
    x3 = x(3); %R
    x4 = x(4); %L
for i=1:1:7
F(i)=((x1-A3(i,1)).^2+(x2-A3(i,2)).^2+(x3-A3(i, 3)).^2).^0.5-((A3(i, 4)