Trajectory planning for uncrewed combat aircraft penetration with dynamic radar cross section constraint

Kinematic model of UCAV

The kinematic model of the mass based on the inertial coordinate system can be described as:

(1) $$\{\begin{array}{rl}\dot{x}(t)& =v(t)cos\theta (t)cos\psi (t)\\ \dot{y}(t)& =v(t)cos\theta (t)sin\psi (t)\\ \dot{h}(t)& =v(t)sin\theta (t)\end{array},$$where $(x(t),y(t),h(t))$ are the coordinates of the position at time $t$, $v(t)$ is the speed of the UCAV. $\theta (t)$ is the track inclination, which is the angle between the velocity and the horizontal plane; $\psi (t)$ is the yaw angle, which is the angle between the projection of the velocity in the horizontal plane and the horizontal axis.

The simplified three-degree-of-freedom model of the UCAV can be expressed as:

(2) $$\{\begin{array}{rl}\dot{v}(t)& =\frac{Fcos\alpha (t)-D}{m}-gsin\theta (t)\\ \dot{\psi }(t)& =\frac{L+Fsin\alpha (t)sin\gamma (t)}{mv(t)cos\theta (t)}\\ \dot{\theta }(t)& =\frac{(L+Fsin\alpha (t))cos\gamma (t)}{mv(t)}-\frac{gcos\theta (t)}{v(t)}\end{array},$$where $g$ is the acceleration of gravity, $\gamma$ is the roll angle, $\alpha$ is the angle of attack, $m$ is the weight of the UCAV, F is the engine thrust, D is the aerodynamic drag, and L is the aerodynamic lift.

Describing the dynamic behavior of UCAVs during real-time trajectory planning based on the three-degree-of-freedom (3DOF) prime assumption:

(3) $$\sum =[x(t),y(t),h(t),v(t),\theta (t),\psi (t),\gamma (t),\alpha (t)].$$

The 3DOF model and the angle relationships are shown in Fig. 1. The drag and lift of the UCAV during flight is modeled as follows:

(4) $$\{\begin{array}{rl}D(t)& =0.5\rho v(t{)}^{2}S{C}_D(t)\\ L(t)& =0.5\rho v(t{)}^{2}S{C}_L(t)\end{array},$$where $\rho$ is the air density, S is the equivalent cross-sectional area, $C_L$ is the aerodynamic lift coefficient, and $C_D$ is the aerodynamic drag coefficient.

Remark 1: The following assumptions are made when establishing 3DOF for the UCAV: The UCAV is a rigid body, with no elastic deformation in the process of movement; the additional effects of the Earth’s rotation, revolution and curvature of the Earth are not taken into account; the tactics used are low-altitude surprise defense, with a relatively small range of change in flight altitude, and the acceleration of gravity and the density of the air are constant within the UAV’s flight altitude; that is to say, the mass of the UCAV is constant in the process of flight.

Dynamic RCS model

The RCS exhibits a direct positive relationship with its detectability by radar: a larger RCS corresponds to a stronger reflected signal captured by the radar, thereby enhancing the likelihood of target detection and recognition. When analyzing radar threats, the relative attitude and position of the UCAV and the radar need to be modeled to evaluate the RCS value (See Fig. 2; Hu et al., 2022).

(5) $$\{\begin{array}{rl}{x}^{\mathrm{\prime }}& =x-x_T\\ y^{\mathrm{\prime }}& =y-y_T\\ h^{\mathrm{\prime }}& =h-h_T\\ \varphi & =arctan(y^{\mathrm{\prime }}/x^{\mathrm{\prime }})\\ \lambda & =\varphi -\psi +\pi \\ \vartheta & =arctan(h^{\mathrm{\prime }}/\sqrt{{(x^{\mathrm{\prime }})}^2+{(y^{\mathrm{\prime }})}^2})\\ r& =\sqrt{{(x^{\mathrm{\prime }})}^2+{(y^{\mathrm{\prime }})}^2+{(h^{\mathrm{\prime }})}^2}\end{array},$$where $(x_T,y_T,h_T)$ is the position of the radar, $\vartheta$ is the height angle between the enemy radar and the UCAV, and $\varphi$ is the azimuth angle between them.

The RCS of an UCAV exhibits dynamic variations throughout its flight, as the real-time angle between the fuselage and the incident radar wave continuously shifts. Its directional scattering behaviors can be categorized into forward, tail, and vertical scattering patterns: Forward Scattering: When radar waves approach directly from the front of the fuselage, the UCAV demonstrates relatively low RCS values due to the alignment that minimizes wave reflection. Tail Scattering: A significant increase in RCS is observed when radar waves impinge from the rear direction, as the structural geometry at the tail section tends to enhance electromagnetic wave reflection. Vertical Scattering: Incidence of radar waves from the top or bottom of the fuselage results in comparatively higher RCS magnitudes, primarily influenced by the vertical surface configurations that offer larger effective reflection areas. The relationship between dynamic RCS $\varsigma$ and UCAV attitude $(\lambda ,\vartheta ,\gamma )$ is designed accordingly:

(6) $$\varsigma =\frac{\pi k_1^2{k}_2^2{k}_3^2}{{(k_1^{2}si{n}^2\delta co{s}^2\eta +k_2^{2}si{n}^2\delta si{n}^2\eta +k_3^{2}co{s}^2\delta )}^2},$$where $k_1=0.3172,k_2=0.1784,k_3=1.003$, and $\delta$ and $\eta$ are process variables.

(7) $$\{\begin{array}{c}\delta =arccos(cos\vartheta cos\lambda )\\ \eta =\gamma -arctan\left(\frac{tan\vartheta }{tan\lambda }\right)\end{array}.$$

The model will be used as one of the constraints.

Radar threat assessment

After the UCAV position and attitude are determined, the dynamic RCS model is used to calculate a. The probability of detection of the UCAV is calculated based on the characteristics of the radar and the distance. Solve for the signal power S according to the following equation.

(8) $$S=\frac{P_T{G}_T{G}_R{\lambda }^{\mathrm{\prime }2}\varsigma }{{(4\pi )}^3{r}_T^2{r}_R^2},$$where $P_T$ is the transmitting power of the radar, $G_T$ is the transmitting antenna gain, $G_R$ is the receiving antenna gain, ${\lambda }^{\mathrm{\prime }}$ is the wavelength, $r_T$ is the distance between the transmitting radar and the target, and $r_R$ is the distance between the receiving radar and the target. When the radar is a single antenna radar, $G_T=G_R$, $r_T=r_R$, then $S=\left(P_T{G}_T^2{\lambda }^{\mathrm{\prime }2}\varsigma \right)/(4\pi {)}^3{r}_T^4$.

The number of false detections per unit time produced by the radar system in a noisy background is Marcum.

(9) $$nfa=\frac{n_{p}lg2}{P_{fa}}.$$

In this expression, the variable $n_p$ represents the number of accumulated pulses, which has been set to a value of 10. Typically, 10 pulses accumulated by a detection radar are sufficient to achieve accurate detection of a target. The variable $P_{fa}$ represents the false alarm probability, which has been set to ${10}^{-11}$.

A threshold is preset as a initialized detection threshold $V_{{\mathrm{\Gamma }}_0}$ to distinguish between the target signal and noise.

(10) $$V_{{\tau }_0}=n_p-\sqrt{n_p}+2.3\sqrt{-lg{P}_{fa}}(\sqrt{-lg{P}_{fa}}+\sqrt{n_p}-1).$$

Update $V_{\tau }$ by iterating until convergence:

(11) $$V_{{\tau }_n}=V_{{\tau }_{n-1}}+\frac{0.5\frac{n_p}{nfa}-{\mathrm{\Gamma }}_I(n_p,V_{{\tau }_{n-1}})}{e^{(n_p-1)ln{V}_{{\tau }_{n-1}}-V_{{\tau }_{n-1}}-ln(n_p-1)!}},$$where ${\mathrm{\Gamma }}_I(a,b)={\int }_0^a\frac{e^{-x}{x}^{b-1}}{(b-1)!}dx$ is the incomplete gamma function. And check the convergence of $V_{\tau }$.

Then the probability that the UCAV is detected by the radar is solved as follows:

(12) $$P=\{\begin{array}{l}{e}^{-\frac{V_{\tau }}{1+n_{p}SNR}},n_p=1\\ 1-{\mathrm{\Gamma }}_I(V_{\tau },n_p-1)+{\left(1+\frac{n_p}{SNR}\right)}^{n_p-1}{\mathrm{\Gamma }}_I\left(\frac{n_{p}SNR{V}_{\tau }}{1+n_{p}SNR},n_p-1\right)e^{-\frac{V_{\tau }}{1+n_{p}SNR}},n_p>1\end{array},$$where SNR represents the signal-to-noise ratio and is related to S.

Objective function establishment

Trajectory planning is essentially a class of constrained optimization problems. During the penetration mission of an UCAV, trajectory planning and radar threat assessment based on dynamic RCS characteristics are paramount considerations for ensuring survivability and mission success. This dual focus on geometric path constraints and RCS-dependent observability metrics establishes a critical framework for minimizing exposure to hostile radar systems.

The Euclidean distance between the candidate track point and the target point is used as a cost indicator $C_d$ to measure the track length to guide the UCAV to fly in the target direction.

(13) $$C_d=\sum _{i=1}^{n_U}\Vert X_{i,t+\mathrm{\Delta }t}-T_i\Vert ,$$where, $X_{i,t+\mathrm{\Delta }t}$ is the candidate position of the i-th UCAV at the time $t$, and $T_i$ is the task target position of the i-th UCAV.

The UCAV based on dynamic RCS is subject to the radar threats of each radar detection at $C_r$.

(14) $$C_r=\sum _{i=1}^{n_U}\sum _{j=1}^{n_R}{P}_{ij}.$$

At $t+\mathrm{\Delta }t$ time, the detection threat coefficient of the i-th UCAV is $C_r$.

The overall objective function of the design is as follows:

(15) $$\begin{array}{c}minf={\xi }_1+{\xi }_2{C}_d+{\xi }_3{C}_r\\ s.t.lb\le [{\alpha }_i(t),{\gamma }_i(t),F_i(t)]\le ub\end{array},$$where, ${\xi }_1,{\xi }_2,{\xi }_3$ are the weight factors. They are used to balance the costs in the objective function $f$. $[{\alpha }_i(t),{\gamma }_i(t),F_i(t)]$ is a set of attack angle, rolling angle and thrust variable. $ub$ is the upper limit and $lb$ is the lower limit.

Constraints design

In addition to avoiding constraints such as enemy radar detection areas during penetration, UCAV also has many constraints. Affected by limited fuel and limited range, the range constraints of reasonable planning of the flight track and optimized flight speed and route length; flight performance constraints subject to the maximum flight speed and other flight performance limitations; mission time constraints that require the completion of penetration and corresponding tasks within the specified time; and altitude cruise conditions.

Consider the performance parameters of UCAV and set flight speed constraint:

(16) $$g_1=max\{v_{min}-v_t,v_t-v_{max}\}\le 0,$$where $v_t$ is the speed of the current moment, $v_{min}$ is the minimum flight speed allowed, and $v_{max}$ is the maximum flight speed.

Constrained by UCAV performance conditions, there is an upper limit in flight, which represents the highest flight altitude that UCAV can achieve, and is usually defined by the vertical distance relative to sea level. At the same time, from a safety perspective, in order to effectively reduce the risk of UCAV colliding with ground obstacles, its minimum flight altitude must be limited. Designing the flight altitude constraints of UCAV:

(17) $$g_2=max\{h_{min}-h_t,h_t-h_{max}\}\le 0,$$where $h_t$ is the height of the current moment, $h_{min}$ is the minimum height, and $h_{max}$ is the maximum height. Normally, $h_{min}=H+H_{safe}$. H is the terrain height corresponding to the track point, and $H_{safe}$ is the safe height of UCAV.

In addition, we do not use flight time as a constraint, but as an evaluation criterion for subsequent algorithms.

Dynamic selection scoring mechanism

In traditional BKA, the sequential execution of attacks and migration behaviors simplifies the algorithm design but sacrifices flexibility, speed, and resource efficiency. In practical applications, a dynamic adaptive strategy needs to replace the fixed stage division in order to better balance the global exploration and local exploitation and improve the algorithm performance.

Fusion of dynamic selection scoring mechanism before migration and attack behavior triggers. Dynamically adjusting attack scores and migration scores with the optimization process. The score for setting up an attack behavior is calculated as follows:

(18) $$Scor{e}_A={\overline{f}}_t-min(f_t),$$where $f_t$ is the overall objective function at time $t$, and ${\overline{f}}_t$ means the mean of $f_t$.

The attack behavior score is calculated as follows:

(19) $$Scor{e}_M=\sigma (f_t),$$where $\sigma$ means the standard deviation.

When the attacking behavior score is greater than the migratory behavior score, the attacking behavior is selected. When the attacking behavior score is less than the migrating behavior score, select the migrating behavior.

Random cross-up update strategy

In traditional BKA, if the fitness value of the current population is less than the fitness value of the random population, then the leader will give up leadership and join the migratory population. However, the location of the migrating population is not optimal, which can easily lead to the population falling into the local optimization trap.

Therefore, when $f_i$ is smaller than $f_{irand}$, the migratory population-based location update strategy is improved to a cross-exploration update strategy. This enhancement scheme is based on the crossover principle of genetic algorithms, promoting population diversity by recombining information from multiple parents (Goldberg, 1989). By employing random perturbation factors and six reference solutions to construct multiple differential vectors, it simulates the multi-parent crossover concept in genetic algorithms (Eiben & Smith, 2015). These vectors facilitate structured yet randomized exploration of the search space, significantly enhancing global search capabilities (Storn & Price, 1997). This improved mechanism contains a randomized perturbation factor G and 6 reference solutions for enhanced exploration capabilities. The cross-exploration strategy is as follows:

(20) $$X_{i,new}=X_{i,a}+G{\mathrm{\Delta }}_1+G({\mathrm{\Delta }}_1-{\mathrm{\Delta }}_2)+G({\mathrm{\Delta }}_2-{\mathrm{\Delta }}_3),$$where $X_{i,new}$ is the new position of the i-th individual, $G\sim U(-1,1)$ is the random perturbation factor, ${\mathrm{\Delta }}_1,{\mathrm{\Delta }}_2,{\mathrm{\Delta }}_3$ are difference vectors generated by randomizing individual positions for enhanced exploration, $X_{i,a}$ is a randomized individual location.

${\mathrm{\Delta }}_1,{\mathrm{\Delta }}_2,{\mathrm{\Delta }}_3$ introduce difference terms to the update strategy.

(21) $$\begin{array}{rl}{\mathrm{\Delta }}_1& =X_{i,a}-X_{i,b}\\ {\mathrm{\Delta }}_2& =X_{i,c}-X_{i,d}\\ {\mathrm{\Delta }}_3& =X_{i,e}-X_{i,f}\end{array},$$where $a,b,c,d,e$ and $f$ are indexes of six different individuals randomly selected from the population.

In this design, the complexity and diversity of the search paths are enhanced by generating nonlinear perturbation terms through the introduction of multi-order differences ${\mathrm{\Delta }}_1,{\mathrm{\Delta }}_2$ and ${\mathrm{\Delta }}_3$. In addition, the stochastic factor G is defined as a dynamic parameter in the interval $(-1,1)$, whose positive and negative values allow the solution vector to perform forward updating or reverse exploration in iterations, thus breaking the singularity of the search direction. At the same time, the algorithm constructs a differential perturbation pattern by introducing six randomly selected reference solutions, which form a multidirectional guiding effect in the solution space, avoiding the algorithm from falling into local stagnation due to over-reliance on a single optimal solution. This strategy, which incorporates multi-order differential perturbation and bi-directional stochastic factors, significantly improves the global search capability and local escape efficiency.

Multi-structure fusion

PO-MIBKA dynamically balances global exploration and local exploitation through a multi-strategy fusion mechanism to significantly improve the solution efficiency of complex optimization problems. In the global exploration stage, the Lévy flight mechanism is introduced, which utilizes its property of combining long step jumps with short step fine search to enable the algorithm to perform non-uniform and efficient exploration in the solution space. Meanwhile, a population diversity maintenance strategy based on individual similarity is designed to dynamically identify the degree of population aggregation and avoid premature convergence by calculating the similarity between individuals. In the local development stage, an adaptive feedback adjustment strategy is used to dynamically adjust the search intensity according to the convergence state of the current iteration.

Lévy flight mechanism: When the algorithm continues to iterate during the optimization process, if the average value of the historical optimal fitness is monitored to remain unchanged for five consecutive generations, the algorithm is determined to enter a stagnant state. To break this local convergence dilemma, the system will automatically activate a random wandering mechanism based on Lévy flights to redistribute the positional coordinates of individuals in the population (Zhang et al., 2024a). This strategy, through the introduced stochastic step size, can realize long-distance exploratory jumps as well as maintain fine search in local areas, thus effectively enhancing the diversity level of the population. The strategy model is designed as follows:

(22) $$X_{i,t+1}=X_{best,t}+R\ast Levy(X_{i,t})+R\ast |X_{best,t}-X_{i,t}|$$

(23) $$Levy(X_{i,t})=\omega s(X_{best,t}-X_{i,t}),$$where R is a random quantity obeying a normal distribution and $\omega \in [-1,1]$ is a scale factor, $s$ is the random wandering step size, which is given by:

(24) $$s=\frac{p}{|q{|}^{\frac{1}{\zeta }}},$$where, both $p\sim N(0,{\sigma }_p^2)$ and $q\sim N(0,{\sigma }_q^2)$ are parameters that follow a normal distribution, and $\zeta \in (0,2)$ is the exponential parameter. ${\sigma }_p$ and ${\sigma }_q$ are as follows:

(25) $${\sigma }_p={\left(\frac{\mathrm{\Gamma }(1+\zeta )sin(\pi \zeta /2)}{\mathrm{\Gamma }((1+\zeta )/2)\zeta 2^{(\zeta -1)/2}}\right)}^{\frac{1}{\zeta }}$$

(26) $${\sigma }_q=1.$$

Individual similarity strategy: To avoid large numbers of black-winged kites choosing similar migration directions during the migration behavior phase, cosine similarity was used to calculate the degree of crowding among individuals.

Construct vectors P and Q:

(27) $$\begin{array}{rl}P& =X_{i,t}-X_{i,best}\\ Q& =X_{j,t}-X_{j,best},(i\ne j)\end{array}.$$

The cosine similarity is calculated as follows:

(28) $$C=\frac{PQ}{|P||Q|}.$$

Other individuals are compared and the one with the smallest cosine similarity is selected as the update direction $X_{new,t+1}$. Calculate the new candidate position.

(29) $$X_{i,t+1}=X_{best,t+1}+r\times (X_{new,t+1}-K{X}_{i,t}),$$where $r\in [0,1]$ is a vector of random elements and K is a random number of 1 or 2.

Feedback mechanism: Feedback regulation is a method of control in which the output of a system is fed into the system as information. In the context of engineering control, Proportional-Integral-Derivative (PID) controller is one of the most frequently used types of automatic controllers, which is designed based on the concept of feedback regulation. Among the many PID control strategies, incremental PID control is a representative recursive algorithm (Li, Ang & Chong, 2006). It is based on the output value of the controlled object to implement the adjustment operation, in order to ensure the stable operation of the system. This control method has a wide range of potential applications and can be introduced into the optimization algorithm (Chen et al., 2024). Specifically, in the optimization algorithm, the search process can be regarded as a system, and the best fitness value obtained at each iteration can be taken as the output value, and incremental PID feedback adjustment can be applied.

Incremental PID output results in $\mathrm{\Delta }u$.

(30) $$\mathrm{\Delta }u=k_p{r}_1(e_k-e_{k-1})+k_i{r}_2{e}_k+k_d{r}_3(e_k-2e_{k-1}+e_{k-2}),$$where $k_p$ is the constant of the proportional term, $k_i$ is the constant of the integral term, $k_d$ is the constant of the differential term, $r_1,r_2,r_3\in [0,1]$ are the random numbers, and $e_k,e_{k-1},e_{k-2}$ are iterative biases.

When $t>1$, the current iteration bias $e_k$ is calculated as:

(31) $$e_k=X_{best,t}-X_{i,t}.$$

The deviation $e_{k-1}$ of the previous iteration is calculated as:

(32) $$e_{k-1}=e_k+\mathrm{\Delta }{X}_{best},$$where $\mathrm{\Delta }{X}_{best}=X_{best,t}-X_{best,t-1}$ is the amount of change in the target position. For ease of calculation, $e_{k-2}$ is calculated as:

(33) $$e_{k-2}=e_{k-1}.$$

When $t=1$, we assume $e_{k-2}=e_{k-1}=e_k$ that for ease of calculation.

Update the position after the feedback mechanism $X_{i,t+1}$:

(34) $$\begin{array}{rl}{X}_{i,t+1}& =X_{i,t}+k_f\mathrm{\Delta }{u}_t+(1-k_f)H_t\\ k_f& =r_{4}cos(t/T)\\ H_t& =(cos(1\text{-}\mathit{t}\mathit{/}\mathit{T})+\rho {\mathit{r}}_{\mathit{5}}){\mathit{e}}_{\mathit{k}}\end{array},$$where $r_4$ and $r_5$ are random elements, $\rho =(ln(T-t+2)/lnT{)}^2$ is regulatory factor. T is the total number of iterations.

Design of PO-MIBKA

Using the proposed improved PO-MIBKA, suitable flight path is found in the framework of the designed constrained optimization problem. In PO-MIBKA, there exist three key and major steps in the algorithm operation, which are the initialization of the population location setting, the execution of the population migration behavior, and the simulation of the population attack behavior. The full process is schematized in Fig. 3.

Initialization of PO-MIBKA: The first step in initializing the population is to create a set of random solutions, where the position of each individual black-winged kite is represented by a matrix. Assuming a population size of $pop$ and a problem dimension size of $dim$, the $j$-th dimensional position of the $i$-th black-winged kite can be represented as ${\mathrm{\Omega }}_{ij}$.

(35) $$\mathrm{\Omega }=\left[\begin{array}{c}{X}_1\\ X_2\\ ...\\ X_{pop}\end{array}\right]=\left[\begin{array}{cccc}{x}_{1,1}& x_{1,2}& ...& x_{1,dim}\\ x_{2,1}& x_{2,2}& ...& x_{2,dim}\\ & & ...& \\ x_{pop,1}& x_{pop,2}& ...& x_{pop,dim}\end{array}\right].$$

The location of each black-winged kite is usually distributed evenly over a range.

(36) $$x_i=x_i^{lb}+rand(x_i^{ub}-x_i^{lb}),$$where $i\in \{1,2,...,pop\}$, $x_i^{lb}$ and $x_i^{ub}$ are the lower and upper bounds of the $j$-th dimension, respectively, and rand is a randomly chosen value between 0 and 1.

After initializing the position, the fitness value of each individual is calculated and the individual with the smallest value is taken as the leader.

(37) $$X_{best,t}=X_{f_{min},t},$$

$X_{best,t}$ is the optimal location for this population at the time $t$.

Attack behavior: The capture of prey by black-winged kites is usually divided into two phases: circling flight in search of prey and swooping flight to feed on prey.

The formula for calculating black-winged kites circling for prey in the traditional BKA causes the update position to appear only in the upper right of the most advantageous point, and does not achieve a circling effect. To solve the above problem, we introduce a spiral factor to update the position. Assuming an iteration period of $T=100$, the possible update positions in the hovering behavior are shown in Fig 4.

Assume $\kappa \in (0,1)$ is a random number. If $\kappa >0.9$, the position is updated when hovering as follows:

(38) $$X_{i,t+1}=X_{i,t}+radiu{s}_{i,t}{v}_{i,t}$$

(39) $$radiu{s}_{i,t}=0.05e^{-2{\left(\frac{t}{T}\right)}^2}(1+sin(\kappa )),$$where $radiu{s}_{i,t}$ is the dynamic radius and $v_{i,t}$ is the spiral factor.

(40) $$v_{i,t}=\frac{1}{\sqrt{n}}\left[\begin{array}{c}cos(k_{1}t)\\ sin(k_{2}t)\\ cos(k_{3}t)\\ ...\end{array}\right],$$where $k_i$ are the frequency parameters of each dimension that control the spiral shape.

If $\kappa \le 0.9$, the black-winged kite makes a swooping attack. This action causes the status to be updated to:

(41) $$X_{i,t+1}=X_{i,t}+0.05e^{-2{\left(\frac{t}{T}\right)}^2}(2\kappa -1)X_{i,t}.$$

Migration behavior: In order to better adapt to change seasons, many birds migrate from north to south during winter to obtain better living conditions and resources. During the migration process, a leader will usually lead the group. The navigational ability of the leader is crucial to the success of the entire group in completing the migration. When the fitness value of the current population exceeds that of the randomly generated population, the leading individual will guide the population to move forward until reaching the target destination. Conversely, when the fitness value of the current population is lower than that of the random population, the leader abdicates its guiding role, and the population migration process is regulated by the proposed stochastic cross-updating mechanism. The model for population migration is described as follows:

(42) $$\begin{array}{rl}{X}_{i,t+1}& =X_{i,a}+G{\mathrm{\Delta }}_1+G({\mathrm{\Delta }}_1-{\mathrm{\Delta }}_2)+G({\mathrm{\Delta }}_2-{\mathrm{\Delta }}_3),f_i<f_{rand}\\ X_{i,t+1}& =X_{i,t}+C\times (X_{best,t}-m{X}_{i,t}),else\\ m& =2sin(\kappa +\frac{\pi }{2})\end{array}.$$

When $f_i<f_{rand}$, the model is updated in a manner consistent with Eq. (20). C is the Corsi mutation.

(43) $$C=\frac{\delta }{\pi ({\delta }^2+{(x\text{-}\mu )}^2)},$$where $\delta =1$ and $\mu =0$.

The above are the three main steps of PO-MIBKA, the pseudo code is shown in Algorithm 1.

Function test on PO-MIBKA

A comparison is made between PO-MIBKA and traditional BKA on classical test functions. The number of iterations for optimization of each function is 1,000 and optimal fitness curves are shown in Fig. 5

The improved PO-MIBKA achieves a faster optimization search speed than the traditional BKA algorithm and significantly improves its ability to find the global optimal solution or a better solution. The comparison of the optimal fitness values of the test functions at 1,000 iterations is shown in Table 1.

Complexity analysis

Based on the position update and fitness calculation mechanism of the PO-MIBKA algorithm, assume the complexity of the fitness function Fit( $\cdot$) is $O(b)$, the search space dimension is $d$, the population size is $n$, and the maximum number of iterations is T. In the initialization phase, the algorithm generates initial solutions and calculates their fitness, with a complexity of $O(nd+nb)$. In the main loop, each iteration includes the Lévy flight mechanism, score calculations, position updates under conditional judgment, and the best solution update. Among these, the Lévy flight and position update operations involve vector computations with a complexity of $O(nd)$; fitness calculation and score calculations require traversing the population, resulting in a complexity of $O(nb+n)$. Since the operations in the conditional branches are all related to the population size and dimension, the worst-case complexity per iteration is $O(nd+nb)$. Therefore, the overall time complexity is $O(Tn(d+b))$.