Location decision of low-altitude service station for transfer flight based on modified immune algorithm

Problem description

According to Guiding Opinions on Construction and Management of the General Aviation Flight Service Station System (CAAC, 2012), the LAFSS is divided into four classes which are nation, region and Class A and B, namely L = {1,2,3,4}. Since the national and regional administrative establishment has been completed, the grade L of LAFSS selected in this article is three or four. Each provincial administrative region meets the requirements of setting up one to three flight service stations of Class A in principle and several Class B as needed. Referring low altitude collaborative management, the categories of service demand in low-altitude flight is divided into two types, namely D = {1,2}, which represents demand in field area or transition corridor respectively. Each demand object has three service coverage modes, namely C = {0,1,2}, which represents independence, joint and repeat coverage respectively. The required service overlay is shown in Fig. 1.

Demand I is a collection of demands from U₁ city and H₂ county, but taken into account the provincial level of P₁, which adopts the regional 1^st level and separate coverage mode. Demand II relies on urban C₂, which adopts the Class A mode and separate coverage. Demand III from county H₃ and the navigable town t₃ belongs to the transfer corridor, which adopts the Class A and joint service coverage mode. Demand IV is part of local navigable town t₂, which is divided into Class A and B and joint coverage mode.

If the service distance or response time from any flight demand area to the candidate site is less than the safety requirements Ø, the demand area is covered by the candidate site. Since the flight service mode may be independence, joint or repeat coverage, the service coverage level function of low-altitude flight demand based on the coverage theory (Yao, 2021) is as follows in Eq. (1).

(1) $$C_L=\sum _{l\in L}{\rho }_l\sum _{c\in C}{\rho }_c+\mathrm{\partial }$$where ${\rho }_l$ is the weight of the adopted 1^st level flight service station. ${\rho }_c$ is the weight of the service coverage method. ∂ is the service random fluctuation coefficient.

Model assumptions

Physical space requirements of low altitude flight include in-field flight and transfer corridor flight. The former is a discrete combination optimization problem, while the latter needs to carry out task relay at each service station within the range of the transfer to achieve service demand coverage within the range of flight path, which belongs to continuous decision planning. Therefore, setting coverage and median quantification (Sharifi et al., 2023), this article establishes site selection mathematical model of multi-objective equilibrium optimization with minimum economic cost, fastest service response and minimum effective service demand distribution distance.

In order to facilitate the construction of mathematical model, the following assumptions are made.

• If service radius is reachable, a flight service station can cover multiple flight demand areas. And

• A flight demand area is served by only one station.

• Different levels of flight service stations correspond to different service radius, flight heights, coverage grade and decision weights.

• The same area of flight demand corresponds to different operational types, but the largest range is chosen.

• The candidate stations for low-altitude flight service are generated at the demand gathering area or nearby towns, without considering random installation.

Parameter setting

Based on the flight course in low-altitude airspace, the relative parameters include in low-altitude transit demand, candidate flight service sites, siting decisions, service response and signal coverage.

There are K transit corridor flights. All the main parameters of every transit corridor are described in set.

• ${\mathrm{\Omega }}_d=[{\omega }_1^d,{\omega }_2^d...{\omega }_k^d...{\omega }_K^d]$ is a set of transfer requirements, k ∈ K.

• ${\omega }_k^d=[P_s^k(x,y,z)...P_t^k(x,y,z)...P_e^k(x,y,z),h_t^k]$ is the position coordinates of any service requirement object, whose position contains start, any moment t and the end.

• S = {S₁, S₂…S_j…, S_N} is a set of candidate sites, j ∈ N. Because of low-altitude flight activity is a socially derived demand, the site selection usually relies on a navigable airport, takeoff-landing point, or a nearby town. Allowing layout located will only be at a given finite number of location points.

• Sj = {o_j, r_j, q_j, h_j, C_j} is a attributes set of any each candidate, which contains center coordinates o_j. service radius r_j, service supply q_j, service coverage height h_j, and service station cost C_j. Service station cost C_j consists of land cost C_j¹, facility cost C_j², operation and maintenance cost C_j³, and scalable set-aside cost C_j⁴.

• Three decision variables are x_j, $y_j^k$ and $g_j^{j+1}$. x_j is assigned to be value 1 if service candidate station selected or 0 otherwise. When the $y_j^k$ is assigned to be value 1 if demand point is covered individually, or 0 otherwise. $g_j^{j+1}$ is assigned to be value 1 if service overlap intersection, or 0 otherwise.

• Service response. $t_j^k$ is the time which low service candidate station S_j responds to demand point ${\omega }_k^d$. $t_j^k$depends on the service distance $l_j^k$ between the geometric center of supply and demand target $o_k^d$ and O_j, the signal propagation speed of the communication navigation facilities adopted by each candidate station $V_c$, the human-machine interaction response time $T_a$ and the signal coverage variance error $T_e$.

Considering the influence of terrain and geomorphic changes and blocking of communication and navigation signals by ground objects, the value of $t_j^i$ follows σ Poisson distribution. $t_j^k$ is taken according to Eq. (2).

(2) $$t_j^k=\sqrt{{[2(l_j^k÷V_c)]}^2+T_a^2+\sigma T_e^2}$$

The radio transmission speed is 3 × 10⁸ m/s (Sakai et al., 2021), while service distance is usually between 200 and 300 km, and the signal transmission V_C does not change much. However, due to the obstruction of the terrain, T_a and T_e are quite different. In this article, the signal coverage errors of LAFSS in plains, hills and plateaus are respectively 1 to 3 s, i.e., $T_e\in [1,3]$, and human-machine interaction response time $T_a\in [1,6]$.

Service coverage analysis

T_k is the total flight time of any transfer flight path ${\omega }_k^d$. The candidate service stations at the start are determined. With the transition flight, there may be multiple candidate service stations o_j and o_j′ and other service radii covering the track range. The service relay process is shown in Fig. 2.

From the start position $P_{\mathrm{s}}^k(x,y,z)$ to the final destination $P_e^k(x,y,z)$, each service station through the whole process can be likened to be completion of the task. Therefore, the location of LAFSS under such flight requirements is optimized. In this article, the flight time is taken as the continuous interval and the service range of the discrete candidate stations around the track is taken as the coverage space to determine the best location scheme and achieve continuous decision planning.

Coordinates conversion

Cartesian coordinates are used to describe the movement process of transition flight more visually and graphically, which are easy to realize the combination of number and shape. The longitude and latitude coordinates obtained by the GIS are set as (B,L,H) (Zhou et al., 2023), and the geographic coordinates are projected to the plane coordinate system. The coordinate conversion mode is shown in Fig. 3.

The transformation equation for the Gauss–Kruger projection is shown in Eq. (3).

(3) $$\begin{array}{c}X=\left(N+H\right)\ast \cos L\\ Y=\left(N+H\right)\ast \cos B\ast \sin L\\ Z=\left[N\ast \left(1-e^2\right)+H\right]\ast \sin B\end{array}\}$$where r_max and r_min is the longest and shortest radius. f is the ellipsoidal flatness. e is the first eccentricity which value is taken by $\sqrt{{\mathrm{r}}_{\mathrm{m}\mathrm{a}\mathrm{x}}^2-r_{min}^2}\mathrm{\setminus }{r}_{max}$. W is the first auxiliary coefficient i.e., calculated by $\sqrt{\left(1-e^2\right)B}$. N is the radius of curvature of the ellipsoidal surface given by $r_{max}\mathrm{\setminus }W$.

Determination of transit trajectory

If calculate roughly, simple aerial linear trajectory is given out. o-xyz coordinate system should be built as geodetic plane. The plane coordinates of the origin place is (0,0,0), V_t is the velocity at any time t_k in any flight path of the aircraft, and α is the angle between direction of travel and ground plane (Højlund et al., 2022). The coordinate system was created as shown in Fig. 4.

Based on the equation of motion with direct track, the coordinates of any position P^k_i accelerated at any time t_k are expressed as in Eq. (4).

(4) $$\begin{array}{rl}& x_t=(v_0{t}_k+{\displaystyle \frac{1}{2}a{t}_k^2+{\displaystyle \frac{1}{3}b{t}_k^3}})\cos \alpha \\ & y_t=\pm \mathrm{\Delta }{y}_{t_k}\\ & z_t=(v_0{t}_k+{\displaystyle \frac{1}{2}a{t}_k^2+{\displaystyle \frac{1}{3}b{t}_k^3}})\sin \alpha \end{array}$$

However, actual low-altitude flight path consists of a combination of straight lines and turning curves. Flight speed is a nonlinear continuous change process from small to large and then decreasing to 0 (Li & Peng, 2007). Therefore, in order to improve coverage accuracy of the service requirements, three-dimensional coordinates of improving hybrid motion trajectory are established in this article, as shown in Fig. 5.

β is the angle between the tangential velocity and the ground plane, V_tk is the instantaneous velocity of the aircraft at time t_k in the flight trajectory, V_t′ is the projection of the instantaneous velocity in the ground plane. γ is the angle between V_t′ and the positive direction of the y-axis. In the turning section, Ψ is defined to describe the angle between the plane and the horizontal plane where its circular route is located. Φ is the angle between the intersection line of the two planes and the y-axis.

Let the coordinates of initial position be (x₀,y₀,z₀). When projected trajectory on the o-x-y plane is a straight line, coordinates of any position P^k_i (x_t,y_t,z_t) at moment t_k is described in Eq. (5).

(5) $$\{\begin{array}{l}{x}_t=x_0+\left(v_{0}tk+\frac{1}{2}a{t}_k^2+\frac{1}{3}b{t}_k^3\right)\cos \beta \sin \gamma \\ y_t=y_0+\left(v_{0}tk+\frac{1}{2}a{t}_k^2+\frac{1}{3}b{t}_k^3\right)\cos \beta \cos \gamma \\ z_t=z_0+\left(v_0{t}_k+\frac{1}{2}a{t}_k^2+\frac{1}{3}b{t}_k^3\right)\sin \beta \end{array}$$

When the projection of trajectory on the o-x-y plane is a curve, the radius of the turn is r, the center of the circle O₁ coordinates (x₀₁,y_01,z₀₁), O₁′ is the O₁ projection, the radius O1C is parallel to AB, C′ is the projection of C.

Thus, the coordinates in the coordinate system at the moment t_k are described in Eq. (6).

(6) $$\{\begin{array}{l}{x}_{t_k}=x_{01}+X\cos \Phi =x_{01}+r\sin \left(\theta 1+\frac{V{t}_k}{r}\right)\cos \Psi \cos \Phi \\ y_{t_k}=y_{01}+Y\sin \Phi =y01+r\cos \left(\theta 1+\frac{V{t}_k}{r}\right)\sin \Phi \\ z_{t_k}=z_{01}+r\sin \theta 2\sin \Psi =z01+r\sin \left(\theta 1+\frac{V{t}_k}{r}\right)\sin \Psi \end{array}$$

The precise coverage range of service requirement for k^th low-altitude transition corridor in target object scope can be determined at any time t_k.

Objective function

We assume that each candidate service site is independent, each service station can provide serve for multiple demand areas, but one demand area is only served by one candidate site. It belongs to the hierarchical setting location problem of different service levels.

According to transfer flight requirements ${\mathrm{\Omega }}_d$, the most appropriate location is select to set up station S_j. At last the optimal site plan S* will be obtained which should ensure flight safety, pay out the minimum cost, obtain the maximum effective coverage area, decrease the minimum service response time and repeating coverage area. In a word, we achieve the maximum economic benefit. However, the optimization object S_j ∈ S is a discrete state in the solution space and belongs to discrete multi-objective combinatorial optimization. Main optimization objectives are shown in Eqs. (7)–(10).

If the candidate site S_j is selected, the sum of construction, operation, maintenance and subsequent scalable uncertain costs will be calculated in Eq. (7).

(7) $$f_1(x)=\sum _{j\in N}{x}_j(C_j^1+C_j^2+C_j^3+C_j^4)$$

Based on the response time priority principle, the distance between candidate sites and flight demand is determined. And then, taking the distance over speed of the propagation of the navigation facilities, the service response time is obtained. The objective function of response time is expressed in Eq. (8).

(8) $$f_2(x)=\sum _{j\in N}{x}_j\sum _{k\in K}{y}_j^k\sqrt{{[2(l_j^k÷V_c)]}^2+T_a^2+\sigma T_e^2}$$where $y_j^k$ means that there are two or more joint coverage, then we will choose the minimum distance by comparing function in Eq. (23).

Compared the distance between the geometric center of two adjacent low altitude flight service candidate sites with service radius, the overlap area of service coverage is judged. The repeating coverage area is calculated by using transformed right-angle plane coordinates. Taking account of all the distances between adjacent candidate sites within study area, the repetitive coverage objective function is expressed in Eq. (9).

(9) $$f_3(x)=\sum _{j\in N}{x}_j{g}_j^{j+1}\left({\displaystyle \frac{{\theta }_j\pi {(r_j^l)}^2+{\theta }_{j+1}\pi {(r_{j+1}^l)}^2}{360}-d_{{}_{j,j+1}}\mathrm{\Delta }{h}_{j,j+1}}\right)$$where d_j,_j+1 is the distance between the geometric centers of adjacent sites, the half-height of overlapping area is Δh_j,j+1, and the corresponding circle center angles of the intersecting areas are θ_j and θ_j+1.

The effective service coverage area not only depends on the number and rank of candidate sites, but also is influenced by the repeating coverage area. The corresponding objective function is expressed in Eq. (10).

(10) $$f_4(x)=\sum _{j\in N}{x}_j\sum _{k\in K}\sum _{j\in N}\sum _{l\in L}{y}_j^k\pi {\left(r_j^l\right)}^2-f_3(x)$$

The service efficiency is not only dependent on the quality and the location of selected service sites, but also is related to the service distribution scheme between sites and demand points closely.

Multi-objective model

The mathematical model of LAFSS location with the shortest total service response time, the lowest cost, the largest service coverage area and the lowest repeated coverage is established based on the multi-objective programming theory as follows in Eqs. (11)–(14).

(11) $$min{Z}_1=\sum _{l\in L}^L{\rho }_l\sum _{j=1}^N(C_j^1+C_j^2+C_j^3+C_j^4)x_j$$

(12) $$min{Z}_2=\sum C_L\sum _{j=1}^N\sum _{k\in K}{x}_i\sum _{j\in K}{y}_j^k\sqrt{{[2(l_j^k÷V_c)]}^2+T_a^2+\sigma T_e^2}$$

(13) $$min{Z}_3=\sum _{j=1}^N\sum _{l\in L}{x}_j{g}_j^{j+1}({\displaystyle \frac{{\theta }_j\pi {\left(r_j^l\right)}^2+{\theta }_{j+1}\pi {\left(r_{j+1}^l\right)}^2}{360}-d_{j,j+1}\mathrm{\Delta }{h}_{j,j+1})}$$

(14) $$max{Z}_4=\sum _{c\in C}{\rho }_c\sum _{j=1}^N{x}_j\sum _{k\in K}{\sum _{j\in N}{y}_j^k\pi \left(r_j^l\right)}^2-Z_3$$

Constraint s.t:

• Service response time is determined by the service coverage altitude of the flight service station and the demand altitude in Eq. (15).

  (15) $$h_k^d/V_c\le t_j^k\le [h_j/(V_c+T_a+T_e)]$$

• The service mode of the demand area is limited by the completed LAFSS in Eq. (16).

  (16) $$\mathrm{\forall }l\in L,c\in C,d\in D,k\in K,j\in N$$

• The number assigned service tasks does not exceed the number of selected service stations in Eq. (17).

  (17) $$Y_j^k\le X_j$$

• Each demand point is served by one and only one service station in Eq. (18).

  (18) $$\sum _{j\in N}\sum _{k\in K}{y}_j^k=1$$

• Decision variables for site selection, service and repeat coverage are list separately from Eqs. (19)–(21).

(19) $$x_j=\{0,1\}$$

(20) $$y_j^k=\{0,1\}$$

(21) $$g_j^{j+1}=\{0,1\}$$

Goal conflict analysis

Mathematical model for site selection of LAFSS is a MOOP. The level, number and location of stations determine the four goals. There are some conflicts between objectives, such as the response time f₂(x) is inversely proportional to the station cost f₁(x), the level and number of stations affecting the station cost f₁(x) are proportional to the service coverage area f₃(x) and the repeating coverage area f₄(x). However, the optimization objectives except f₄(x) are as small as possible. The conflict analysis among the four objective functions is shown in Fig. 6.

All objectives are interrelated and in conflict with each other, it is impossible to obtain a global optimal solution that makes multiple objectives reach the optimal. Usually, we can only make balance between multiple targets in a given area.

Algorithm design

Combined four objectives for service safety with efficiency constraints, this article proposes an improved two-layer response of immune optimization algorithm based on the operating mechanism of biological immune system to solve the MOOP of location decision of LAFSS. The IA has strong adaptability and robustness in solving MOOP (Jerne, 1974) and is a new intelligent optimization algorithm constructed artificially (Hajela & Lee, 1996). Compared with genetic algorithm, IA can produce a variety of specific antibodies, which has been proved to be an effective way to generate Pareto-Edge worth end faces (de Castro & Timmis, 2002).

The decision process is a kind of multi-process parallel optimization, which can obtain multiple sub-optimal solutions while seeking the optimal solution. In order to avoid the loss of the optimal equilibrium solution due to premature maturity, local adjustment of station rank and clonal reproduction of adjacent stations are used.

• Antibody encoding

  Besides the candidate station sites are selected or not in the process of LAFSS location planning, there are also corresponding service station grades. Otherwise, this article chooses natural number coding. The IA treats the problem to be solved as the antigen which is the demand area of low-altitude flight covered by the service. The solution of the problem is the antibody, which is the best location. Each antibody corresponds to a location scheme. The value of the loci in the antibody corresponds the candidate site selected or not and the grade.

Natural number coding is applied, and the encoding length of the solution is N i.e., the number of candidate sites. The antibody code is denoted by $A_m=(S_1^l,S_2^l...S_j^l...S_N^l)$, where $j\in N,l\in L,m\in M$. A_m represents the situation where the location scheme corresponding to the m^th antibody is selected for the j^th candidate site and the service station level l corresponding to the selected site. All the situations selected are as following in Eq. (22).

(22) $$S_j^l=\{\begin{array}{c}4,\mathrm{S}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{e}\mathrm{d}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\mathrm{e}\mathrm{d}\mathrm{a}\mathrm{s}4\\ 3,\mathrm{S}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{e}\mathrm{d}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\mathrm{e}\mathrm{d}\mathrm{a}\mathrm{s}3\\ 2,\mathrm{S}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{e}\mathrm{d}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\mathrm{e}\mathrm{d}\mathrm{a}\mathrm{s}2\\ 1,\mathrm{S}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{e}\mathrm{d}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\mathrm{e}\mathrm{d}\mathrm{a}\mathrm{s}4\\ 0,\mathrm{N}\mathrm{o}\mathrm{t}\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{e}\mathrm{d}\end{array}$$

• Generation of initial population

The generation of initial population plays an important role in optimization iteration of IA. In order to construct a feasible initial solution quickly, this article introduces the concept of clustering effect in economics, i.e., the centripetal force generated by the spatial concentration of various industries to attract economic activities to a certain region. According to the requirements of existing stations, initial candidate points of LAFSS are locked in the existing navigable airports, take-off and landing points or cities. The range of candidate points is further narrowed, the number of iterations is reduced, and the solution accuracy is increased.

Calculated track position $P_t^k(x,y,z)$ at any time t of the transition flight corridor k, the candidate station is classified as one of the initial population members when the distance $d_{t_j}^k$ is less than service radius r_j. The judgment function is shown in Eq. (23).

(23) $$d_{tj}^k={|\mathrm{O}}_j^{}(x,y,z)P_t^k(x,y,z)|\le r_j$$

• Colony diversity

The weight method is commonly used to transform the MOOP into a single objective optimization problem (Wang, 2017), affinity evaluation function between antibodies and antigens is shown in Eq. (24).

(24) $$f(S_j^l)={\tau }_{{}_1}{\displaystyle \frac{1}{f_1(x)}+{\tau }_{{}_2}{\displaystyle \frac{1}{f_2(x)}+{\tau }_{{}_3}{\displaystyle \frac{1}{f_3(x)}+{\tau }_{{}_4}{f}_4(x)}}}$$

The value of ${\tau }_{{}_1},{\tau }_{{}_2},{\tau }_{{}_3}\mathrm{a}\mathrm{n}\mathrm{d}{\tau }_{{}_4}$ is the weight coefficient of each objective function, and the sum is 1. This function can evaluate the merits of the candidate site solution. The value of the evaluation function $f(S_j^l)\in (0,f_4(x)]$, the larger the value, the better the candidate site solution.

A deformable R−bit continuous method is used to calculate affinity between antibodies in Eq. (25).

(25) $${\mathrm{S}}_{v,s}(i)=k_{{\mathrm{A}}_i{\mathrm{A}}_m}\mathrm{\setminus }{L}_{A_m}$$where L_Am is the total length of the antibody code. $k_{{\mathrm{A}}_i{\mathrm{A}}_m}$ is the same number of coding bits for both antibodies.

Antibody concentration represents the ratio of similar antibodies among antibody groups that also means the proportion of similar solutions in the site selection of LAFSS. It can characterize the diversity of antibody population (Cheng & Cui, 2020).

(26) $${\mathrm{C}}_{{\mathrm{A}}_{\mathrm{m}}}=\sum _{i=1}^M{B}_{AS}\mathrm{\setminus }M$$where M is the size of the population, and $B_{\mathrm{A}\mathrm{S}}$ indicates the characteristic value of the inter-antibody affinity.

When the inter-antibody affinity $S_{v,s}>\epsilon$ is present, the $B_{l,v}$ = 1. Otherwise $B_{l,v}=0$. where $\epsilon$ is a predefined threshold value. Excitation degree is the final evaluation result of antibody quality, which is usually obtained by mathematical operation of antibody affinity ${\mathrm{S}}_{v,s}(i)$ and antibody concentration C(A_m). The function of excitation degree is as follows in Eq. (27).

(27) $$Si{m}_{\mathrm{A}\mathrm{S}}(\mathrm{A}\mathrm{m})=\tau \times {\mathrm{S}}_{\mathrm{v},\mathrm{s}}(i)-\ell \times \mathrm{C}(\mathrm{A}\mathrm{m})$$where $\tau ,\ell$ indicates the affinity and concentration adjustment factors.

The matching between location scheme and optimization problem can stimulate the generation of new location schemes (Bin, Zhao & Yang, 2018). The probability that an individual is selected for a mutation operation is the expected reproduction probability in Eq. (28).

(28) $${\text{P}}_{\text{S}}({\text{A}}_m)=\Gamma \frac{f\left(S_j^l\right)}{{\displaystyle \sum f\left(S_j^l\right)}}+\left(1-\Gamma \right)\frac{C_{A_M}}{{\displaystyle \sum C_{A_M}}}$$where Г is the constant range from 0 to 1. The greater the affinity between antibody and antigen $f(S_j^l)$, the greater the expected reproduction probability P_S and the greater the possibility of being selected as a mutant. The higher the antibody concentration C(A_m), the lower the expected reproduction probability P_S and the lower the possibility of being selected as a mutant individual.

Improving operator

IA is usually insufficient ability to maintain population diversity and easy to fall into local optimization, which leads to reduce convergence speed and optimization performance greatly accompanied by site selection scale. In this article, clone variation, variable threshold selection operator and population update are adopted to preserve the self-regulation mechanism of the original immune optimization algorithm, which is helpful to increase the diversity of location schemes, improve the global optimization ability of the algorithm, and enhance the accuracy of the late solution of the algorithm to achieve rapid convergence.

Each antibody in the memory bank is cloned, the cloning scale is B_c, and each clone copy is mutated according to probability p_m. The clone reproduction function B(A_m) is calculated as follows in Eq. (29).

(29) $$\text{B}({\text{A}}_{\text{m}})=round\left[\zeta A_m\backslash \left(C_{A_M}{\displaystyle \sum _j^{U}f(S_j^l)}\right)\right]$$where $\zeta$ is the cloning factor, $f(S_j^l)$ is the value of the affinity function of the antibody, $C_{A_m}$ is the value of the antibody A_m concentration value, $\sum _j^{U}f(S_j^l)$ is the sum of antibody affinity values in the clone parent population, and U is the clone parent population size, round(·) is the integer function.

The immune selection operator can determine which antibodies can enter the clonal selection operation based on the expected selection probability (Xiong & Gong, 2019). Setting different thresholds is conducive to improving the solving efficiency according to different search periods (Yoo & Hajela, 1999). For this purpose, antibody concentration is changed in stages by a piecewise function, the equation is as follows in Eq. (30).

(30) $$T=\{\begin{array}{c}a,ift\le t_{ave}\\ a\sin \left({\displaystyle \frac{\pi }{2}\times {\displaystyle \frac{t_{max}-t}{t_{max}-t_{ave}}}}\right),ift>t_{ave}\end{array}$$where $\vartheta$ is the threshold selection factor. t indicates the current iteration number. t_ave is the average of evolutionary generations. t_max indicates the maximum iteration number.

The first antibodies N_p with better affinity and the re-initialized $Popsize-Np$ antibodies were selected to form a new generation population. According to the three modes of service coverage proposed in Section 1, some location points of the diversion track may be covered by two or more flight service stations at the same time, as shown in Fig. 7.

$P_t^k(x,y,z)$ indicates that the location belongs to joint coverage service. $d_{ti}^k$ and $d_{tj}^k$ indicate the distance between the location and two adjacent service sites, respectively. In this article, the signal space masking is not considered for the time being. Adjusting the service level of candidate sites, changing their service radius, and re-calculating their respective service distances, the affinity value of the IA is determined in Eq. (31).

(31) $${\mathrm{S}}_{v,s}^l(i)={\displaystyle \frac{r_i^l}{r_i^l+d_{ti}^k}}$$

Implementation flow and algorithm complexity

MIA in this article adopts the group search strategy and the optimal solution with high probability is got through iterative calculation (Wang, 2017). The implementation process of the MIA is shown in Fig. 8.

The MIA in this section consists of three main modules (Yang et al., 2019), which are antigen recognition and initial antibody generation, intrinsic immune evaluation and adaptive immune manipulation.

The solution steps are description from 1^st to 8^th step.

1^st Step: Antigen identification. According to the set of low-altitude transit flight regions ${\mathrm{\Omega }}_d=[{\omega }_1^d,{\omega }_2^d...{\omega }_k^d...{\omega }_K^d]$, establish the time and position of track coordinates for any transit corridor k respectively, convert the transit corridor flight into coordinate relay, construct the affinity function $d_t^k(ij)$ and various constraints.

2^nd Step: Generate the initial population. Based on the actual flight trajectory coordinates, determine all the candidate station solutions that satisfied $d_{ti}^k\le r_i^l$ around the t-trajectory at any moment $A_m=(S_1^l,S_2^l...S_j^l...S_N^l)$.

3^rd Step: Affinity evaluation. The affinity functions $f(S_j^l)$ and ${\mathrm{S}}_{v,s}(i)$ are called to evaluate the antibody-antibody affinity eigenvalues $B_{\mathrm{A}\mathrm{S}}$ in the population $S_{v,s}$.

4^th Step: Inherent immunity termination judgment. Judge whether the weight coefficients of each objective function ${\tau }_{{}_1},{\tau }_{{}_2},{\tau }_{{}_3}\mathrm{a}\mathrm{n}\mathrm{d}{\tau }_{{}_4}$ satisfy the termination condition, if so, the algorithm terminates, turn to 8^th step. otherwise, turn to 5^th step to continue the second layer of immune search.

5^th Step: Calculation of antibody concentration and excitation degree. According to the principle of lowest rank l configuration of candidate sites and minimum service distance priority $\mathrm{m}\mathrm{i}\mathrm{n}{d}_t^k(ij)$, based on the concentration of antibody C(A_m) and the degree of excitation $Si{m}_{AS}(A_m)$.

6^th Step: Adaptive immune processing. Including the immune selection threshold $T\left(m\right)$, clonal propagation B(A_m), variation probability P_s(A_m), to determine the antibody with high excitation is selected to enter the clonal selection operation (Chen & Zhang, 2023).

7^th Step: Population refresh. Using a change in site level, the population update function ${\mathrm{S}}_{v,s}^l(i)$ is called to replace the less motivated antibodies in the population with the generated new antibodies, and determine whether a new generation of antibodies is formed, and if so, go to 6^th step to continue the search for superiority. otherwise, go to next step.

8^th Step: Site selection decision scheme output. Output the location, number and grade of the final site, and call the objective calculation functions f₁(x), f₂(x), f₃(x) and f₄(x) respectively to get the indicators of the optimal solution.

Considering the population size assumed in Problem description, there are k transition flight corridors. Each flight corridor has P_t^k location points, and there are N candidate service sites. However, according to the service coverage requirements of each transition flight corridor k, only site N^k that $d_{ti}^k$ is less than or equal $r_{ti}^l$ can be included in the search scope.

According to the rules for calculating algorithm complexity (Qian, 2017), the total time complexity of executing one generation of search for the immune search algorithm for this site selection is shown in Eq. (32).

(32) $$T(n)=2O(N^k/2)+3O(N^k)+O(3N^k/2)+O((2N_{{}^c}^k){)}^2+O(2N_{{}^c}^k\log (2N_{{}^c}^k))$$where N^k denotes the number of candidate sites around the k^th transit flight path. $N^k{}_c\left(N^k{}_c\le N^k\right)$ means the number of candidate sites after immune siting and ranking.

Demand characteristics

This article take the service coverage of low-altitude transition corridor flight in Sichuan as the research object to solve and simulate the multi-objective double-layer immune location decision. The results are compared with the traditional dynamic programming and other intelligent search algorithms to verify the feasibility and optimization effect of the model and algorithm. The MIA for LAFSS is simulated by the MATLAB platform. The final visualization of location decision is carried out based on the Spring open source application framework of Java platform. Easyui is used to obtain Cesiumjs map data interface through the front-end development, and the site selection scheme is presented with the terrain data of the map of the world as the background.

Sichuan Province is one of autonomous flight areas to reform the coordinated management of low altitude airspace in China, and has explored and practiced the integrated operation mode of low altitude airspace based on visual autonomous flight. By the end of 2022, 385,378 general aviation flights and 114,606 h have been achieved (Hu, 2021). The Low Altitude Airspace Cooperative Operation Center of Sichuan Province which was completed in 2018 has formulated regulations on the coordinated management and use of low altitude airspace.

Reference to “Sichuan General Aviation Industry Development Plan (2019–2025)”, the implementation strategy of general aviation industry will create a core-two wings, multi-point N network (Sichuan Provincial Development and Reform Commission, 2017). According to the policy planning and development strategy, this article selects the east-west and north-south transfer flight corridors respectively. The schematic diagram of track layout is shown in Fig. 9.

Low altitude in Sichuan has set up five cooperative management airspace and five low altitude visual flight corridors, a total of more than 6,600 km² of pilot airspace (Sichuan Province National Defense Science Technology and Industry Office, 2018). The two transfer flights selected in this article span four regions which are with large differences in terrain change and economic development. The north-south transfer is mainly used for navigation training, industrial and agricultural operations and sightseeing business. The east-west transfer is a high-altitude flight, which undertakes the task of emergency rescue. The codes of main requirement identification points on the two transfer paths are shown in Table 1.

The coordinate conversion method described in Coordinates conversion section is used to convert geography to Cartesian. Based on the location $P_t^k(x,y,z)$ of every track at any time t, all candidate stations that meet the service distance corresponding are searched nearby the current station level.

Coding

According to the data preparation requirements of IA, all the candidate sites with class A and Class B in the vicinity of two diversion flight paths is coding. The encoding results are shown in Table 2. There are 84 candidate sites, which includes 13 existing airports and 71 planned airports. The candidate sites are classified five parts according to topographic changes. The Class B flight service stations that have been built are Zigong and Guangyuan respectively, and the initial allocation in coding is the same as that of other candidate stations, and the grade is assigned to be four.

By encoding, a scheme A_m is formed, i.e., a solution. The total encoding length N of the solution is 84, i.e., the number of candidate sites Fig. 10. The part of the solution ${\omega }_k^d$ represents the number of the requirement point that the solution satisfies. The m part of the solution represents the candidate site label j. The s part of the solution indicates whether or not the corresponding candidate site is selected and the level of selection.

${\omega }_k^d$ indicates that the solution meets the requirement. The candidate site NO.1 has a service station level of 4. Candidate site NO.2 was not selected. The candidate site NO.4 has a service station level of 3. The candidate site NO.73 has a service station level of 1. The candidate site NO.74 has a service station level of 2.

Parameter taking

The sites with level 3 and level 4 layout are mainly selected in the scope of the example in this article. The required infrastructure construction, production and operation, post-maintenance, and reservation costs are random numbers between 5 to 3, 3 to 2, 2 to 1, and 1 to 0.5 million, respectively. It is also assumed that costs other than infrastructure are proportional to operating time. The costs for each candidate site are shown in Table 3.

The four key parameters in the algorithm are memory capacity overbest, expected reproduction probability parameter ps, mutation probability pm, and concentration threshold. Suppose that the population size of the algorithm is sizepop = 50, the number of iterations is 600, and other parameters are shown in Table 4.

Table 4 shows four reasonable level values for the four key parameters of the experiment, with each pair of parameters being run 20 times.

Dynamic programming site selection based on service distance

The dynamic programming (DP) algorithm (Annear, Akhavan-Tabatabaei & Schmid, 2023) sets the general aviation airports that the transition flight track passes through as node $P_t^k(x,y,z)$, determines the distance between each node and the nearby candidate stations according to the screening condition $d_{ti}^k\le r_{si}^l$ of effective service distance, and obtains the dynamic programming path with service radius as the cost. Taking the north-south transition low-altitude general navigation flight as an example, the minimum cost dynamic programming graph as shown in Fig. 11 is obtained according to the distance of the candidate stations around the track.

The path from node P¹_i to P¹_j is unidirectional and bidirectional, and the path without arrow is bidirectional, which provides a variety of options. Starting from P¹_S, according to the 10 important nodes experienced by the flight trajectory, the service distance was calculated in stages, and the optimal location and service allocation scheme under this method was finally selected. The results of the dynamic site selection are shown in Fig. 12.

The number of sites selected by the dynamic programming scheme is eight, including one in the region and seven in class B. According to the principle of minimum cost and maximum flow, a regional LAFSS is used to cover Chengdu Plain, north Sichuan and east Sichuan.

However, its utilization seems to be extremely high, but the quality of service is affected in many areas due to signal blocking. Therefore, in this article, the inherent and adaptive characteristics of double-layer immunity, dynamic selection and clonal mutation, flexible adjustment of the level, find the optimal scheme.

MIA site selection optimization with two-layer response

Two demand corridors are set for transition flight, and each corridor involves 12 key turning points. The number of candidate sites is 84, which are two Class B flight service stations, 16 general aviation airports, and 68 cities and towns with planned general aviation airports. The level initialization of all candidate stations is set to four, and the service radius is set to 100 km. According to the safe flight criterion that the flight demand is covered by the service, the initial population is obtained as shown in Fig. 13.

The initial population consisted of 14 flight service stations, including three in Chengdu Plain, two in East Sichuan Hills, two in north Sichuan, three in south Sichuan and four in West Sichuan Plateau. Xindu in Chengdu Plain, weiyuan in southeast Sichuan and yingshan in northeast Sichuan have higher service efficiency. Except for qingchuan in the northwest of Sichuan, there are two service stations, and the remaining 10 stations all cover a navigable airport. Therefore, the repeated coverage areas are more in north Sichuan, Chengdu Plain and east Sichuan hills.

Innate immunity is a natural immune defense function formed during the development of the body. Zigong Fengming Airport S33, as an existing Class B flight service station, is compared to the non-specific defense function that has been available since birth. It is used as an immune body to search for germline development and evolution. In this article, according to the cost fluctuation of site selection, the auto-immune search was carried out according to the service radius of 120 and 150 km respectively. Finally, the innate immunity site selection schemes P¹_IA and P²_IA were obtained respectively, as shown in Table 5.

The comparison of the two site selection schemes is known, the quantity of service station is different by two Class B, the cost is increased by 22.22%, the maximum service demand ratio P²_IA is increased by 11.11%, and the equilibrium P¹_IA is increased by 11.43%. Based on Minitab Statistical Software platform, this article used the Taguchi experimental design method to discuss parameter setting and sensitivity analysis. The horizontal trend of the influence factors shown in Fig. 14 was obtained by using multivariate analysis of variance.

Therefore, other relevant parameters of the algorithm were set according to the horizontal trend in Fig. 15. Population size sizepop = 50, memory bank capacity overbest = 8, expected reproduction probability parameter ps = 0.9, mutation probability pm = 0.7, concentration threshold ε = 0.7.

Due to the initial effect of innate immunity with invasive antigen substances, the demand for flights in the western Sichuan region is low, so that the number of stations in this region for both schemes is only one. In order to avoid the disadvantages of poor compressive strength and no specific selectivity, adaptive immunity needs to be introduced. Thus, the advantages of randomness, parallelism, global convergence and population diversity are utilized to obtain an optimized location decision scheme.

The location of LAFSS in Sichuan Province can be considered as containing five main antigens from east, west, north, south and Chengdu Plain respectively, and multiple sub-antigens and antibody groups can be decomposed from each region through adaptive adjustment. In this way, more excellent antibodies can be provided for future high-level immunity, and multiple extremes or global equilibrium optimal search characteristics can be played. The change characteristics of Sv, s(i) and f(S^l_j) are obtained in Fig. 15.

In order to show the comparison of each index function more intuitively, this article adopts normalization processing. The term of the maximum value of the current objective function is set as one, and the rest is taken in this proportion.

The statistical characteristics of economic cost ratio are 1.49 in Chengdu Plain, 1.59 in north Sichuan, 1.3 in east Sichuan, 1.83 in south Sichuan and 2.29 in west Sichuan. The highest service response ratio is 2.83 in western Sichuan, which is 3.5 times higher than the lowest of 0.89 in Chengdu Plain. The highest proportion of repeated coverage is 2.6 in western Sichuan, which is 25 times higher than the lowest in Chengdu Plain. The highest effective coverage ratio is 1.95 in Chengdu Plain, which is 2.5 times higher than 0.76 in western Sichuan Plateau. However, the highest Sv, s(i) and f(S^l_j) belong to the south and east regions of Sichuan where all costs are more balanced, respectively.

According to the variation rules of Sv,s(i) and f(S^l_j), cloning will continue to be used to reproduce the above two regional good genes. After repeated coverage of mutations and selection, Chengdu Plain, as the center of southwest China, has been established as an excellent gene, four site selection indicators are well balanced, and total affinity is the best. The other two candidate sites for expansion and upgrading are “yuexi” in southern Sichuan and “aba” in western Sichuan. After 600 iterations, the mating probability is 0.6 and the mutation probability is 0.1. The decision scheme and service allocation as shown in Fig. 16 are obtained through the simulation platform with MIA.

According to the dynamic programming decision in Fig. 12 and the site selection results completed in this section, the index comparison of the three algorithms is obtained. The comparison results are shown in Table 6.

The optimization effect of the MIA on the location decision of LASS in Sichuan Province is as follows: the cost index Z1 is reduced by 26.46%, the service response is increased by 25%, the repeated coverage is reduced by 29.5%, and the effective service coverage area is increased by 17%. In terms of the time complexity of solving the problem, the iteration efficiency was increased by 36%, and the total immune affinity was increased by 9.5%.