A plague epidemic disease optimizing approach

Abstract

In order to solve some optimization problems with multi-local optimal solutions, a plague infectious disease optimization (PIDO) algorithm is proposed by the dynamic model of plague infectious disease with pulse vaccination and time delay. In this algorithm, it is assumed that there are several villagers living in a village, each villager is characterized by some characteristics. The plague virus is prevalent in the village, and the villagers contract the infectious disease through effective contact with sick rats. The plague virus attacks is the few characteristics of the human body. Under the action of the plague virus, the growth status of each villager will be randomly transformed among 4 states of susceptibility, exposure, morbidity and recovery, thus a random search is achieved for the global optimal solution. The physical strength degree of villagers is described by the human health index (HHI). The higher the villager’s HHI index, the stronger the physique and the higher the surviving likelihood. 9 operators (S_S, S_E, E_E, E_I, E_R, I_I, I_R, R_R, R_S) are designed in the PIDO algorithm, and each operator only deals with the 1/1000∼1/100 of the total number of variables each time. The case study results show that PIDO algorithm has the characteristics of fast search speed and global convergence, and it is suitable for solving global optimization problems with higher dimensions.

Keywords

Swarm intelligence optimization global optimization plague transmission dynamic model random search

1 Introduction

When solving global optimal solutions with a large number of local optimal solution optimization problems, the swarm intelligence optimization algorithm has unique advantages [1]. When an optimization problem is solved in this kind of algorithm, many trial solutions are started at the same time, and iterative calculations are implemented. Some very difficult multi-local optimization problems can be solved by this “group attack” method. Unfortunately, there is no swarm intelligence optimization algorithm that can solve all types of optimization problems [2]. People have proposed a lot of swarm intelligence optimization algorithms, such as genetic algorithm [3], ant colony algorithm [4, 5], particle swarm algorithm [6 –8], biogeographic algorithm [9], differential evolution algorithm [10, 11], bee colony algorithm [12].

There are a large number of local optimal solution optimization problems to solve global optimal solutions. Multi Access Interference (MAI) is the main source limiting the capacity and quality of the Multiple-Input Multiple-Output Ortho-gonal Frequency Division Multiplexing (MIMO-OFDM) system which fulfills the demand of high-speed transmission rate and high quality of service for future underwater acoustic (UWA) communication. Multi User Detection (MUD) is needed to overcome the per-formance degradation caused by MAI. Both local and global optimal solutions are obtained in Bionic Binary Spotted Hyena Optimizer (BBSHO) algorithm using the Position Coordinate Vectors (PCVs) of the social behavior of spotted hyenas to achieve MUD [13]. Further, Extremal Optimization (EO) is introduced in BBSHO algorithm to improve the local search ability within the search space. With the measurement of the Earth’s magnetic field, magnetic compass can provide high frequency heading information. However, it suffers from local magnetic interference. An intelligent ellipsoid calibration method based on the grey wolf is proposed to generate optimal parameters for magnetic compass to generate high performance heading information [14]. With the analysis of the projection relationship among the navigation coordinate frame, the body frame and the local horizontal frame, the heading ellipsoid equation is constructed. Time evolutionary optimization (TEO) is a novel population-based meta-heuristic optimization algorithm, inspired by natural selection and evolution of creatures over time. Time and the environment are two main factors of evolution at TEO. Enhanced time evolutionary optimization (ETEO) is presented [15]. ETEO is the new version of TEO which modifies time evolutionary factor and applied population clustering. Population clustering amplified environmental factor to increase the efficiency of ETEO. For this purpose, a memory is used to save some best designs and ETEO can escape from local optimal points. The algorithm was validated by solving several constraint benchmarks and engineering design problems. The material 20MnMoNi55 reactor pressure vessel (RPV) steel has safety-critical applications and may experience cyclic loading during atypical operating conditions. Simulation of the material behaviour subjected to cyclic loading requires suitable kinematic hardening rules and accurate material parameters. The determination of Chaboche kinematic hardening parameters from experimental data is not well defined and the same is true for the cyclic yield stress. These parameters are interdependent and need optimization, maintaining the same hierarchy of dependence. Evolutionary search algorithms are implemented in nested applications to find the material parameters phenomenologically [16]. The optimum material parameters offer good simulation of the actual material behaviour in cyclic plastic loading.

Many-objective optimization problems (MaOPs) contain four or more conflicting objectives to be optimized. A number of efficient decomposition-based evolutionary algorithms have been developed in recent years to solve them. However, computationally expensive MaOPs have been scarcely investigated. Typically, surrogate-assisted methods have been used in the literature to tackle computationally expensive problems, but such studies have largely focused on problems with 1–3 objectives. An approach called hybrid surrogate-assisted many-objective evolutionary algorithm is presented to solve computationally expensive MaOPs [17]. The key features of the approach include: (1) the use of multiple surrogates to effectively approximate a wide range of objective functions; (2) use of two sets of reference vectors for improved performance on irregular Pareto fronts (PFs); (3) effective use of archive solutions during offspring generation; and (4) a local improvement scheme for generating high-quality infill solutions.

In the swarm intelligence optimization algorithm, each trial solution is compared to an individual with biological characteristics. Therefore, some special biological activity scenarios are used to construct swarm intelligence optimization algorithms. However, the biological activity scenarios that have been proposed for swarm intelligence optimization algorithms are very simple and have weak mathematical foundations, such as genetic algorithm, ant colony algorithm, biogeographic algorithm, and bee colony algorithm. Some swarm intelligence optimization algorithms do not even have biological activity scenarios, such as particle swarm algorithm and differential evolution algorithm. If a swarm intelligence optimization algorithm has a clear and vivid biological activity scene, and the scene can be accurately described by mathematical theory, then this scene is extremely conducive to constructing a swarm intelligence optimization algorithm. Because the swarm intelligence optimization algorithm is based on good mathematics, its performance is easy to be analyzed. The plague infectious disease optimization (PIDO) in this paper is based on such a biological activity scenario.

Plague is a powerful infectious disease transmitted by Yersinia pestis through rodent fleas[18], it has a very high fatality rate. Plague has been three pandemics in human history. It is a natural foci disease that is widespread among wild rodents. In order to control the plague’s hazards, tremendous efforts has been spent in the world, and a large number of mathematical models are built to study the transmission mechanism of the plague. It is mainly manifested in the following four aspects:

Differential equations are used to establish the spreading dynamics model of the plague between the rats and human, the model is used to describe the epidemic characteristics of the plague and analyze the relationship between the various parameters in the model and the development of the epidemic [19 –22].

Statistical mathematical models are used to study the correlation between natural factors (such as climate, temperature, rainfall.) of the plague foci and rodent density, vector density, plague infection rate and other parameters, and their hysteresis characteristics [23, 24].

The biological characteristics of the plague rat is studied from the aspects of ecological adaptability, reproduction, growth and development, blood biochemistry, hibernation and its own pathogen spectrum, and a plague animal model is established, basis is provided to reveal the pathogenic mechanism, evaluate the vaccine effect, and treat the plague for program selection [25, 26].

Screen out the monitoring indicators that are most relevant to the epidemic of animal plague, construct a forecast and early warning model of plague epidemic in plague foci, these are used to control the large-scale outbreak of plague [27, 28].

When solving the optimization problem, a heuristic search strategy is used in the swarm intelligence optimization algorithm, and a large number of individuals are started to search at the same time. The behavior of individuals and the synergy between individuals have become very important. Individual behavior can be described by operators, and the synergy between individuals is the logical structure of the swarm intelligence optimization algorithm. The PIDO algorithm is based on a time-lag plague infectious disease model with impulse vaccination [18], where the biological meaning of an individual is a human individual in a certain ecosystem. The synergy between individuals is reflected by the transmission mechanism of the plague infectious disease, and the individual’s behavior is reflected by the individual’s random transition between the four states of susceptibility, exposure, disease, and cure.

The PIDO algorithm adopts a different design idea from the existing swarm intelligence algorithm, and the time-delay plague epidemic model of impulse vaccination is transformed into a method that can solve a large number of local optimal solution optimization problems, and the constructed operator can be sufficient reflecting the interaction relationship of the time-delayed plague epidemic model of impulse vaccination, thus the basic idea of the plague epidemic dynamics theory is reflected, the algorithm has global convergence.

2 Optimal algorithm design method for plague infectious disease

Suppose the optimization problem to be solved is Equation (1): ${\begin{matrix} min f (X) \\ X \in H \subset R^{n} \end{matrix}$ (1)

In the Equation (1), Rⁿ is the n-dimensional Euclidean space, n is the total number of real variables which are included in the optimization problem; X = (x₁,x₂,...,x_n) is an n-dimensional solution vector; f(X) is the objective function, The mathematical expression of f(X) has no restrictions; H is the search space.

2.1 PIDO algorithm scene design

Suppose there are N villagers living in a village Z, and each villager is represented by a number as 1,2,...,N; each villager is represented by n features, and one feature is equivalent to a certain organ of the villager, namely for the i-th villager, its characteristic is expressed as (x_i,1,x_i,2,...,x_i,n), i = 1∼N; the plague is prevalent in the village Z, and villagers make effective contact with rats carrying the plague virus, such as being bitten by a rat, eating rat meat, or accidentally eating food contaminated by rat feces, they will be infected with the plague, which will spread among the villagers; the plague virus attacks certain organs of the villagers. Villagers who are infected with the plague are called susceptible. After they are infected with the plague, they will not immediately become ill, and the virus in their bodies enters the incubation period; the villagers who have the plague virus in the incubation period are called as the exposed persons, and they will also pass the virus to other villagers effectively in contact; those exposed after the incubation period will become ill. These villagers are called the sick, and they will pass the virus in their bodies to other villagers who are in effective contact with them; the sick and exposed people can be cured through medical treatment. Those who are cured and exposed are called cured; in order to prevent the plague virus from harming the villagers, the villagers will be vaccinated every once in a while. After vaccination, the villagers will not be 100% successfully immunized; the villagers who are successfully immunized will not be infected with the plague virus for a period of time; the villagers who have not been successfully immunized are still susceptible; the immunity of the successfully immunized villagers will automatically disappear and lose immunity after a period of time; villagers who have not been immunized or have lost immunity will be infected with the plague virus again. Under the influence of the plague virus, the growth status of each villager in the village will switch between the four states of susceptibility, exposure, disease, and cure. This conversion is mapped to the H of the optimization problem Equation (1), which is equivalent to each trial solution jumping from one position to another in H, so the purpose of random search for H is achieved. The physical strength of each villager is determined by the characteristics of the villager. A villager with a strong physique can continue to grow, while a villager with a weak physique stops growing.

The above scene is projected into the search for the global optimal solution of the optimization problem Equation (1), the meaning is as follows:

The H of the optimization problem corresponds to the village Z. A villager in this village corresponds to a tentative solution of the optimization problem Equation (1). The set of tentative solutions corresponding to N villagers is X = {X₁,X₂,...,X_N}, X_i = (x_i,1, x_i,2 ... ,x_i,n), i = 1∼N. A feature j of villager i corresponds to a variable x_i,j of solution X_i. Therefore, the villager i is equivalent to the trial solution X_i. The physical strength of the villagers is represented by the HHI (human health index, HHI) index, and the HHI index of the villager i is calculated by the following formula (2): $HHI (X_{i}) = {\begin{matrix} \frac{1}{1 + f (X_{i})}, if f (X_{i}) > 0 \\ 1 + | f (X_{i}) |, if f (X_{i}) ⩽ 0 \end{matrix}$ (2)

HHI is an index for evaluating the quality of the solution.

2.2 Delayed Plague Infectious Disease Model with Pulse Vaccination

For the convenience of discussion, the villagers living in a village are divided into 4 categories:

Susceptible (category S): all villagers who are not infected with the plague. If the S category villagers effectively contact the villagers with the plague virus, they may be infected with the plague.

Exposed persons (Category E): namely, villagers who have been infected with the plague virus but have not yet developed the disease, and they are potential patients.

Those who are sick (Class I): Villagers who have been infected with the plague virus and have become ill. Class I villagers who have effective contact with Class S villagers may transmit the virus to Class S villagers.

Cure (category R): the healer of E or I villagers. This group of villagers will not get the disease temporarily, but after a period of time, if they effectively contact the E or I disease villagers, they will return to sick.

The delayed plague infectious disease model with impulse vaccination is shown in formula (3) [18].

${\begin{matrix} \frac{dS (t)}{dt} = μ - β I (t) (1 + vI (t) + μ) S (t) - q μ I (t), t \neq kT \\ \frac{dI (t)}{dt} = - β e^{- μ τ} I (t - τ) (1 + vI (t - τ) S (t - τ)) + {qe}^{- μ τ} I (t - τ) - (r - μ) I (t), t \neq kT \\ \frac{dR (t)}{dt} = r I (t) - μ R (t), t \neq kT \\ E (t) = 1 - S (t) - 1 - I (t) - R (t), t \neq kT \\ {S (t}^{+}) = S (t) - b μ (S (t) + E (t) + R (t)) - bp μ I (t), t = kT \\ {E (t}^{+}) = E (t), t = kT \\ {I (t}^{+}) = I (t), t = kT \\ R {(t}^{+}) = R (t) - b μ (S (t) + E (t) + R (t)) - bp μ I (t), t = kT \end{matrix}$ (3)

S(t) is the proportion of susceptible people at time t; E(t) is the proportion of exposed persons at time t; I(t) is the proportion of sick people at time t; R(t) is the proportion of cure people at time t. In formula (3), t is the period; R(t), I(t), E(t), and S(t) are the proportions of people belonging to R, I, E, and S respectively in period t, and R (t)≤0, I(t)≤0, E(t)≤0, S(t)≤0, R(t) + I(t) + E(t) + S(t) = 1. Others’ parameters are shown in Table 1 [18]. In the table, Rand(A,B) means to generate a uniformly distributed random number in the interval [A,B]; INT(w) means to round the real number w to integer.

Table 1

Meanings of parameters and method of taking value

Parameter name	Parameter meaning	Value restrictions	Value method
μ	Birth rate	μ> 0	μ^t = Rand(0.15, 0.35)
β	Linear infection rate	β> 0	β^t = Rand(0.15, 0.35)
ν	Nonlinear infection rate	ν≤0	ν^t = Rand(0.015, 0.085)
q	Vertical transmission rate	q > 0	q^t = Rand(0.15, 0.35)
p	Horizontal infection rate	p > 0	p^t = Rand(0.15, 0.35)
r	Cure rate	r > 0	r^t = Rand(0.15, 0.35)
b	Immunization success rate	b > 0	b^t = Rand(0.15, 0.35)
τ	Incubation period	τ> 0	τ_i = INT(Rand(1, 4)), i = 1∼N
V	Immunization period	V > 0	V_i = INT(Rand(1, 4))–1, i = 1∼N
T	Vaccination cycle	T is a positive

At period t, because any person can only be in one of R, I, E, and S, R(t), I(t), E(t), S(t) are actually the probability that a person belongs to R, I, E, or S; when a person belongs to R, I, E, or S, it means that a person is in R state, I state, E state or S state. The following Equations (4–16) are contributed in this paper.

Therefore, Equation (3) can be applied to any person in the crowd, namely:

${\begin{matrix} \frac{{dS}_{i} (t)}{dt} = μ - β I_{i} (t) (1 + {vI}_{i} (t) + μ) S_{i} (t) - q μ I_{i} (t), t \neq kT \\ \frac{{dI}_{i} (t)}{dt} = - β e^{- μ τ i} I_{i} (t - τ_{i}) (1 + {vI}_{i} (t - τ_{i}) S_{i} (t - τ_{i})) + {qe}^{- μ τ_{i}} I_{i} (t - τ_{i}) - (r - μ) I_{i} (t), t \neq kT \\ \frac{{dR}_{i} (t)}{dt} = r I_{i} (t) - μ R_{i} (t), t \neq kT \\ E_{i} (t) = 1 - S_{i} (t) - I_{i} (t) - R_{i} (t), t \neq kT \\ S_{i} (t^{+}) = S_{i} (t) - b μ (S_{i} (t) + E_{i} (t) + R_{i} (t)) - bp μ I_{i} (t), t = kT \\ E_{i} (t^{+}) = E_{i} (t), t = kT \\ I_{i} (t^{+}) = I_{i} (t), t = kT \\ R_{i} (t^{+}) = R_{i} (t) - b μ (S_{i} (t) + E_{i} (t) + R_{i} (t)) - bp μ I_{i} (t), t = kT \\ i = 1 \sim N \end{matrix}$ (4)

Equation (4) is used to calculate the probability that each person is in the R state, I state, E state and S state during the period t.

The values of the t parameters b, r, p, q, v, β, and μ are respectively b^t, r^t, p^t, q^t, v^t, β^t, μ ^t ; for the convenience of calculation, the formula (4) is changed to recursive discrete expression (5):

${\begin{matrix} S_{i} (t + 1) = S_{i} (t) + μ^{t} - β^{t} I_{i} (t) (1 + v^{t} I_{i} (t) + μ^{t}) S_{i} (t) - q^{t} μ^{t} I_{i} (t) \\ I_{i} (t + 1) = I_{i} (t) - β^{t} e^{- μ^{t} τ_{i}} I_{i} (t - τ_{i}) (1 + v^{t} I_{i} (t - τ_{i}) S_{i} (t - τ_{i})) + q^{t} e^{- μ^{t} τ_{i}} I_{i} (t - τ_{i}) - (r^{t} - μ^{t}) I_{i} (t) \\ R_{i} (t + 1) = R_{i} (t) + r^{t} I_{i} (t) - μ^{t} R_{i} (t) \\ E_{i} (t) = 1 - S_{i} (t) - I_{i} (t) - R_{i} (t) \\ i = 1 \sim N, t \neq kT \end{matrix}$ (5) ${\begin{matrix} S_{i} (t + 1) = S_{i} (t) - b^{t} μ^{t} (S_{i} (t) + E_{i} (t) + R_{i} (t)) - b^{t} p^{t} μ^{t} I_{i} (t) \\ E_{i} (t + 1) = E_{i} (t) \\ I_{i} (t + 1) = I_{i} (t) \\ R_{i} (t + 1) = R_{i} (t) + b^{t} μ^{t} (S_{i} (t) + E_{i} (t) + R_{i} (t)) + b^{t} p^{t} μ^{t} I_{i} (t) \\ i = 1 \sim N, t = kT \end{matrix}$ (6)

In Equations (5) and (6), R_i(t), I_i(t), E_i(t), and S_i(t) are the probability that individual i belongs to the R, I, E, and S groups, respectively, at period t. R_i(t)≤0, I_i(t)≤0, E_i(t)≤0, S_i(t)≤0, R_i(t) + I_i(t) + E_i(t) + S_i(t) = 1; The values of parameters b^t, r^t, p^t, q^t, v^t, β^t, μ^t, τ _i are shown in Table 1.

At period t, randomly generated villagers’ b^t, r^t, p^t, q^t, v^t, β^t, μ^t, τ_i, and formulas (5) and (6) are used to calculate villagers i’s R_i(t), I_i(t), E_i(t) and S_i(t). In the period t, the villager i is in R state, I state, E state and S state. Which of the four states is he in? This is determined by the probability of R_i(t), I_i(t), E_i(t) and S_i(t), that is, which of R_i(t), I_i(t), E_i(t) and S_i(t) is larger, the probability that the corresponding state is selected is also greater. Table 2 shows the spread of plague virus among villagers.

Table 2

Transformation of legal states of plague virus transmission

Period t-1	Period t	Legal state transition	Operator name
S	S	S⟶S	S_S
	E	S⟶E	S_E
E	E	E⟶E	E_E
	I	E⟶I	E_I
	R	E⟶R	E_R
I	I	I⟶I	I_I
	R	I⟶R	I_R
R	R	R⟶R	R_R
	S	R⟶S	R_S

The nine legal state transitions shown in Table 2 can be represented in Fig. 1, where each arc in the figure corresponds to an operator. The 9 legal state transitions shown in Table 2 can be described by 9 operators: S_S, S_E, E_E, E_I, E_R, I_I, I_R, R_R, R_S.

2.3 Feature set generation method

The generation method of the dominant villager set PS^s: L villagers are randomly selected from the villagers in state s. These villagers have a higher HHI index than the current villager i, forming the dominant villager set PS^s, s∈{R,I,E,S}; L is called the number of villagers exerting influence.

The generation method of category villager set CS^s: L villagers are randomly selected from the villagers in state s to form category villager set CS^s, s∈{R,I,E,S}.

Fig. 1

9 legal state transitions.

2.4 PIDO evolution operator design method

S_S operator. The S_S operator is described as the situation where the villagers in the S state at time t have not contracted infectious diseases at time t + 1. For villager i in the susceptible state S, there is formula (7): $v_{i, j} (t + 1) = {\begin{matrix} \sum_{s \in {PS}^{S}} α_{s} x_{s, j} (t), | {PS}^{S} | > 0 \\ x_{i, j} (t), | {PS}^{S} | = 0 \end{matrix}$ (7) In formula (7), v_i,j(t + 1), x_i,j(t) are the state values of the characteristic j of villager _i in period t + 1 and period t respectively; α _s is a constant, α _s = Rand(0.3,0.7).

S_E operator. The S_E operator is described as the situation where the villagers who were in the S state at the time t + 1 were infected with the plague virus through effective contact with the villagers who were already in the E or I state at the time t + 1. Let L villagers in E state or I state randomly select the state value of feature j, after processing, pass to the corresponding feature j of villager i in S state, so that they will be infected with the plague virus and become exposed. That is, for villager i in S state, there is formula (8). $v_{i, j} (t + 1) = {\begin{matrix} \sum_{s \in {CS}^{E} \cup {CS}^{I}} α_{s} x_{s, j} (t), | {CS}^{E} \cup {CS}^{I} | > 0 \\ x_{i, j} (t), | {CS}^{E} \cup {CS}^{I} | = 0 \end{matrix}$ (8)

E_E operator. The E_E operator is described as a villager who has been in E state at time t is still in the incubation period at time t + 1 because the incubation period has not expired. Let the state value of the characteristic j of the villager in the E state but whose HHI index is higher than that of the current villager i be processed and passed to the corresponding characteristic j of the villager i in the E state, the HHI index is increased. That is, for villager i in E state, there is formula (9). $v_{i, j} (t + 1) = {\begin{matrix} \sum_{s \in {PS}^{E}} α_{s} x_{s, j} (t), | {PS}^{E} | > 0 \\ x_{i, j} (t), | {PS}^{E} | = 0 \end{matrix}$ (9)

E_I operator. The E_I operator is described as a villager who has been in E state at time t, at time t + 1, the incubation period has expired and the plague virus begins to attack. After the state value of the feature j of the L villagers in the I state is processed, it is transmitted to the corresponding feature j of the villager i in the E state to make them sick. That is, for villager i in E state, there is formula (10). $v_{i, j} (t + 1) = {\begin{matrix} \sum_{s \in {CS}^{I}} α_{s} x_{s, j} (t), | {CS}^{I} | > 0 \\ x_{i, j} (t), | {CS}^{I} | = 0 \end{matrix}$ (10)

E_R operator. The E_R operator is described as a villager who has been in the E state at time t, and is cured by treatment at time t + 1. Let the state value of the characteristic j of the L villagers in the R state be processed, and pass to the corresponding characteristic j of the villager i in the E state to make them healed. That is, for villager i in E state, there is formula (11). $v_{i, j} (t + 1) = {\begin{matrix} \sum_{s \in {CS}^{R}} α_{s} x_{s, j} (t), | {CS}^{R} | > 0 \\ x_{i, j} (t), | {CS}^{R} | = 0 \end{matrix}$ (11)

I_I operator. The I_I operator is described as a villager who has been in the I state at time t and is still in the I state at time t + 1. Let the state value of the characteristic j of the villager in the I state but whose HHI index is higher than that of the current villager i be processed and transmitted to the corresponding characteristic j of the villager i in the I state to increase the HHI index. That is, for villager i in I state, there is formula (12). $v_{i, j} (t + 1) = {\begin{matrix} \sum_{s \in {PS}^{I}} α_{s} x_{s, j} (t), | {PS}^{I} | > 0 \\ x_{i, j} (t), | {PS}^{I} | = 0 \end{matrix}$ (12)

I_R operator. The I_R operator is described as a villager in the I state at time t, who is cured by treatment at time t + 1. Let the state value of the characteristic j of the L villagers in the R state be processed and transmitted to the corresponding characteristic j of the current villager i to make them healed. That is, for villager i in I state, there is formula (13). $v_{i, j} (t + 1) = {\begin{matrix} \sum_{s \in {CS}^{R}} α_{s} x_{s, j} (t), | {CS}^{R} | > 0 \\ x_{i, j} (t), | {CS}^{R} | = 0 \end{matrix}$ (13)

R_R operator. The R_R operator is described as a villager who has been in the R state at the time t, and is still in the R state at the time t + 1. Let the state value of the characteristic j of the villager in the R state but whose HHI index is higher than that of the current villager i be processed and transmitted to the corresponding characteristic j of the villager i in the R state to increase the HHI index. That is, for villager i in R state, there is formula (14). $v_{i, j} (t + 1) = {\begin{matrix} \sum_{s \in {PS}^{R}} α_{s} x_{s, j} (t), | {PS}^{R} | > 0 \\ x_{i, j} (t), | {PS}^{R} | = 0 \end{matrix}$ (14)

R_S operator. The R_S operator is described as a villager who was in the R state at the time t and turned to the S state due to the loss of immunity at the time t + 1. Let the state value of the feature j of the L villagers in the S state be processed and sent to the corresponding feature j of the current villager i to turn it into the S state. That is, for villager i in R state, there is formula (15). $v_{i, j} (t + 1) = {\begin{matrix} \sum_{s \in {CS}^{S}} α_{s} x_{s, j} (t), | {CS}^{S} | > 0 \\ x_{i, j} (t), | {CS}^{S} | = 0 \end{matrix}$ (15)Growth operator. This operator describes the growth of villagers, there is formula (16).

$X_{i} (t + 1) = {\begin{matrix} V_{i} (t + 1), ifHHI (V_{i} (t + 1)) > HHI (X_{i} (t)) \\ X_{i} (t), Others \\ i = 1 \sim N \end{matrix}$ (16)

In formula (16), X_i(t) = (x_i,1(t),x_i,2(t),...,x_i,n(t)), V_i(t + 1) = (v_i,1(t + 1), v_i,2(t + 1),...,v_i,n(t + 1)). HHI(X_i(t)) and HHI(V_i(t + 1)) are calculated according to formula (2).

Equations (7)–(15) are the four state transition metrics, that is, the process of each step of the algorithm. Equation (16) is to ensure that each step of the algorithm solution will not be worse than the previous step.

2.5 Construction of PIDO algorithm

Step 1 (initialization):

Let t = 0. The number of evolutionary periods G = 8000 ∼60000, the error requirement ɛ= 10^–5∼10^–10, N = 50∼500, the probability of successful vaccination of villagers Q₀ = 0.5∼0.9, the probability of villagers being affected E₀ = 1/1000∼1 /100, L = 3∼6, T = 3∼10.

Randomly select the trial solution [X₁(0), X₂(0),..., X_N(0)] corresponding to N villagers in H.

Let V(i) = 0, i = 1∼N. V(i) = 0 means that the vaccination of villager i is unsuccessful or not; V(i) > 0 means that the vaccination of villager i is successful.

Randomly select the global optimal solution Y* initial value.

Step 2: $a_{k}^{i}$ = Rand(0,1), k = 1∼4; $S_{0} = \sum_{k = 1}^{4} a_{k}^{i}$ , $S_{i} (0) = a_{1}^{i} / - S_{0}$ , $E_{i} (0) = a_{2}^{i} / - S_{0}$ , $I_{i} (0) = a_{3}^{i} / - S_{0}$ , R_i (0) = 1 - S_i (0) - E_i (0) - I_i (0), i = 1∼N.

Step 3: Calculate the RIES status of villager i, RIES_i(0) = RIES(S_i(0), E_i(0),I_i(0), R_i(0)), i = 1∼N. Among them, RIES_i(0) is the state of villager i in period 0; function RIES(S_i(0), E_i(0), I_i(0), R_i(0)) is used to determine which state of villager i will be.

For t = 0 to G

Calculate the parameters b^t, r^t, p^t, q^t,ν^t,β^t,μ ^t according to the method given in Table 1;

For all s∈{R, I, E, S}, generate feature sets PS^s and CS^s;

For i = 1 to N

τ _i = INT(Rand(1, 4)), V_i = INT(Rand(1, 4)) –1;

If T cannot divide t then

Calculate R_i(t + 1), I_i(t + 1), E_i(t + 1) and S_i(t + 1) according to formula (5);

Else

q₀ = Rand(0,1)

If q₀≤Q₀ then

Calculate R_i(t + 1), I_i(t + 1), E_i(t + 1) and S_i(t + 1) according to formula (6);

V(i) = t + 1;

Else

Calculate R_i(t + 1), I_i(t + 1), E_i(t + 1) and S_i(t + 1) according to formula (5);

End if

For j = 1 to n

p₀ = Rand(0,1)

If p₀≤E₀ then

If RIES_i(t) = S then

If RIES_i(t + 1) = S then

Execute the S_S operator according to formula (7) to obtain v_i,j(t + 1);

Else if RIES_i(t + 1) = E and V(i) = 0 then

LP(i) = t + 1;

Execute the S_E operator according to formula (8) to obtain v_i,j(t + 1);

Else

v_i,j(t + 1) = x_i,j(t), RIES_i(t + 1) = RIES_i(t);

End if

Else

if RIES_i(t) = E then

If RIES_i(t + 1) = E then

Execute the E_E operator according to formula (9) to obtain v_i,j(t + 1);

Else if RIES_i(t + 1) = I and (t + 1-LP(i)) > τ _i then

Execute the E_I operator according to formula (10) to obtain v_i,j(t + 1);

Else if RIES_i(t + 1) = R then

Execute the E_R operator according to formula (11) to obtain v_i,j(t + 1);

Else

v_i,j(t + 1) = x_i,j(t), RIES_i(t + 1) = RIES_i(t);

End if

Else if RIES_i(t) = I then

If RIES_i(t + 1) = I then

Execute the I_I operator according to formula (12) to obtain v_i,j(t + 1);

Else if RIES_i(t + 1) = R then

Execute the I_R operator according to formula (13) to obtain v_i,j(t + 1);

Else

v_i,j(t + 1) = x_i,j(t), RIES_i(t + 1) = RIES_i(t);

End if

Else if RIES_i(t) = R then

If RIES_i(t + 1) = R then

Execute R_R operator according to formula (14) to obtain v_i,j(t + 1);

Else if RIES_i(t + 1) = S and (t + 1-V(i)) > V_i then

V(i) = 0;

Execute R_S operator according to formula (15) to obtain v_i,j(t + 1);

Else

V_i,j(t + 1) = x_i,j(t), RIES_i(t + 1) = RIES_i(t);

End if

Else

V_i,j(t + 1) = x_i,j(t), RIES_i(t + 1) = RIES_i(t);

End if

End for

Execute the growth operator according to formula (16) to obtain X_i(t + 1);

End for

If |X^*t +1-Y^*|< ɛ then

Go to step 4;

End if

Y* = X*^t +1;

End for

Step 4 (End):

The function RIES(p_S, p_E, p_I, p_R) is defined as follows:

Function RIES(p_S, p_E, p_I, p_R) // p_S, p_E, p_I, p_R are the probability of occurrence of S, E, I, R state

respectively

w = Rand(0,1);

If w≤p_S then Return S;

Else if p_S < w≤p_S + p_E then Return E;

Else if p_S + p_E < w≤p_S + p_E + p_I then Return I;

Else if p_S + p_E + p_I < w≤p_S + p_E + p_I + p_R then Return R;

End if

End function

2.6 Features of PIDO algorithm

Time complexity. The time complexity calculation process of PIDO algorithm is shown in Table 3.

Markov characteristics. From Equations (7) to (15) defined by the operators of S_S, S_E, E_E, E_I, E_R, I_I, I_R, R_R, R_S, etc., any new trial solution X(t + 1) is only related to its state X(t) in period t, and it has nothing to do with its previous evolution to the current state, indicating that the evolution process of PIDO algorithm has Markov characteristics.

From the definition of growth operator (16), it is known that the HHI index of any villager in period t + 1 will never be lower than its HHI index in period t, it shows that the evolution process of PIDO algorithm has the characteristic of “not bad step by step".

Table 3
Time complexity table for PIDO algorithm

Operating Time complexity Maximum number of cycles

Initialization O(4n + 3nN + n(5 + 7 N) + 1

3 N(n² + 4n + 8))

Operators: S_S, S_E, E_E, E_I, (N + 8)(G + 6)

E_R, I_I, I_R, R_R, R_S O((N + 4L + 6)nE₀/12)

State preservation 4 + 3n (N + 8)(G + 6)

Objective function O(n²) (N + 8)(G + 6)

Growth operator O(4n) (N + 8)(G + 6)

Result output O(n) 1

Operating	Time complexity	Maximum number of cycles
Initialization	O(4n + 3nN + n(5 + 7 N) +	1
	3 N(n² + 4n + 8))
Operators: S_S, S_E, E_E, E_I,	(N + 8)(G + 6)
E_R, I_I, I_R, R_R, R_S	O((N + 4L + 6)nE₀/12)
State preservation	4 + 3n	(N + 8)(G + 6)
Objective function	O(n²)	(N + 8)(G + 6)
Growth operator	O(4n)	(N + 8)(G + 6)
Result output	O(n)	1

3 Proof of global convergence for PIDO algorithm

According to PIDO algorithm, the village Z is equivalent to the search space H, and the N villagers {X₁ (t) , X₂ (t) , . . . , X_N (t)} in the village are arranged into a sequence X^t₁, X^t₂,..., X^t_N. The number of villagers is N, a villager is equivalent to a trial solution, and its objective function value is f(X^t_i), and the set of the states of all villagers is F: $F = {f (X_{i}^{t}) | X_{i}^{t} \in H}$ . Further order: $F = {F_{1}, F_{2}, \dots, F_{N}}, F_{1} ⩽ F_{2} ⩽ \dots ⩽ F_{N}$ (17)

Without loss of generality, F₁ is also desired to make the global optimal solution. From the subscript of formula (17), the set U = {1,2,...,N} can be obtained:

∀i∈U, i is the state that villager i may be in when performing random search. Suppose that the best objective function value found by villager i in period t is F_i, and its corresponding state is i. Obviously, from Equation (17), when the villager i searches at time t + 1, if it transitions to a better state k, it should satisfy k < i; on the contrary, if it transitions to a worse state k, it should satisfy k > i, as shown in Fig. 2.

Fig. 2

State transition of PIDO algorithm in random search.

∀X∈H, F₁≤F(X)≤F_N, H is divided into subsets: $X_{S}^{i} = X | X \in H$ and F (X) = F_i, i = 1, 2, . . . , N. Obviously, State switching in $X_{S}^{i}$ cannot change the objective function value. In addition, there are obviously:

$\begin{matrix} \sum_{i = 1}^{n} | X_{S}^{i} | = N; \forall i \in {1, 2, \dots, N}, X_{S}^{i} \neq φ; \\ \forall i \neq j, X_{S}^{i} \cap X_{S}^{j} = φ, \cup_{i = 1}^{N} X {S^{i}}_{S}^{i} = H \end{matrix}$

Let X_i,j(i = 1,2,...,N,j = 1,2,...,|Xⁱ_S|) denote the j-th state in Xⁱ_S. During the evolution of the villagers, the transition from one state (i,j) to another state (k,l) can be expressed as X_i,j⟶X_k,l. Assumptions: the transition probability from X_i,j to X_k,l is p_ij,kl, the transition probability from X_i,j to any state in X^k_S is p_ij,k. From any state in Xⁱ_S to any one of X^k_S, the transition probability of the state is p_i,k, then the following formula: $\begin{matrix} p_{ij, k} = \sum_{l = 1}^{| X_{S}^{k} |} p_{ij, kl}, \sum_{k = 1}^{N} p_{ij, k} = 1, p_{i, k} ⩾ p_{ij, k} \\ p_{i, k} ⩾ p_{ij, k}, \Rightarrow \sum_{k = 1}^{N} p_{i, k} ⩾ \sum_{k = 1}^{N} p_{ij, k} = 1 and 0 ⩽ \sum_{k = 1}^{N} p_{i, k} ⩽ 1 \Rightarrow \\ \sum_{k = 1}^{N} p_{i, k} = 1 \end{matrix}$ (18)

Lemma 1: In the PIDO algorithm, ∀X_i,j∈Xⁱ_S,i = 1,2,..., N,j = 1,2,...,| Xⁱ_S |, satisfying: $\forall k > i, p_{i, k} = 0$ (19) $\exists k < {i, p}_{i, k} > 0$ (20)

Proof:

The proof of lemma (19): suppose state i is the state of villager i in period t, and its corresponding position in H is X^t. From Section 2.6, the random search process of PIDO algorithm has not bad characteristic of step by step. Therefore, in period t + 1, villager i will not transfer to any worse state, as shown in Fig. 2, with the following formula: $\begin{matrix} f (X^{t + 1}) ⩽ f (X^{t}) \Rightarrow \forall k > i, p_{ij, kl} = 0 \\ \Rightarrow \forall k > i, p_{ij, k} = \sum_{l = 1}^{| X_{S}^{k} |} p_{ij, kl} = 0 \Rightarrow \forall k > i, p_{i, k} = 0 \end{matrix}$

The proof of lemma (20): Suppose state i is the state of villager i at period t. In period t + 1, villager i randomly chooses S_S, S_E, E_E, E_I, E_R, I_I, I_R, R_R, R_S operator to evolve, it is moved to a better state k. At this time, there are two situations as follows:

If state i is the global optimal state, i.e. i = 1, the next transition will remain in the original state, i.e. k = i = 1. This is because from Section 2.6, villager i will not move to another state worse than the original state, so they must stay in the original state i with the probability p_1,1 = 1. Because p_1,1 = 1 > 0, the proposition is proved.

If state i is not the global optimal state, there must be at least one intermediate state k between state 1 and current state i, as shown in Fig. 2, making F₁≤F_k < F_i, that is, 1≤k < i, At that time, villager i can move from the current state i to a new and better state k, that is, p_i,k > 0. The proposition is proved.

In summary, available ∃k < i, p_i,k > 0.

Theorem 1: Let P ′ be a reducible random matrix of order n, that is, $P^{'} = [\begin{matrix} C & \dots & {0 \end{matrix} R \dots {T]$ can be obtained after the same row and column transformation, where C is a primitive random matrix of order m, and T ¬= 0, R ¬= 0, then: $\begin{matrix} P^{' \infty} = lim_{k \to \infty} P^{' k} = lim_{k \to \infty} \\ [\begin{matrix} C^{k} & \dots & 0 \\ \sum_{i = 1}^{k - 1} T^{i} {RC}^{k - i} & \dots & T^{k} \end{matrix}] = [\begin{matrix} C^{\infty} & \dots & 0 \\ R^{\infty} & \dots & 0 \end{matrix}] \end{matrix}$

The above matrix is a stable random matrix, P ^∞ = 1′P ^∞, P ^′∞ = P ^′0 P ^′∞ is uniquely determined and has nothing to do with the initial distribution, P^′∞ satisfies the condition: $P^{' \infty} = {[p_{ij}]}_{n \times n}, {\begin{matrix} p_{ij} > 0, 1 ⩽ i ⩽ n, 1 ⩽ j ⩽ m \\ p_{ij} = 0, 1 ⩽ i ⩽ n, m < j ⩽ n \end{matrix}$

The proof process of Theorem 1 can be found in literature [29].

Theorem 2: PIDO algorithm has global convergence.

Proof: From Section 2.6, the search process of PIDO algorithm has Markov characteristics. For each Xⁱ_S,i = 1,2,...,N can be regarded as a state on a finite Markov chain. According to the conclusion of formula (19) in the lemma, the state transition matrix of the Markov chain can be obtained as the following formula: $P^{'} = [\begin{matrix} p_{1.1} & 0 & \dots & 0 \\ p_{2.1} & p_{2.2} & \dots & 0 \\ ⋮ & ⋮ & ⋮ \\ p_{N, 1} & p_{N . 2} & \dots & p_{N, N} \end{matrix}] = [\begin{matrix} C & 0 \\ R & T \end{matrix}]$

From Equation (18), the sum of the probabilities of each row in P ′ is equal to 1. According to the conclusion of formula (20) in the lemma, the following formula can be obtained: $\begin{matrix} p_{2.1} > 0, R = {(p_{2.1}, p_{3.1}, p_{4.1}, \dots, p_{N, 1})}^{T}, \\ T = [\begin{matrix} p_{2.2} & \dots & 0 \\ ⋮ & ⋮ \\ p_{N, 2} & \dots & p_{N, N} \end{matrix}] \\ \neq 0, C = (p_{1.1}) = (1) \neq 0 \end{matrix}$

It can be seen from the above that P ′ is an N-order reducible random matrix, that is, the Markov state transition matrix, which satisfies the conditions of Theorem 1, and there is the following formula: $P^{' \infty} = lim_{k \to \infty} [\begin{matrix} C^{k} & \dots & 0 \\ \sum_{i = 1}^{k - 1} T^{i} {RC}^{k - i} & \dots & T^{k} \end{matrix}] = [\begin{matrix} C^{\infty} & \dots & 0 \\ R^{\infty} & \dots & 0 \end{matrix}]$

Since C ^′∞ = C = ( 1 ), T ^∞ = 0 , there must be R ^∞ = (1, 1, …, 1) ^T. This is because the sum of the probabilities of any row in P^′ is always 1, and there is formula:

$P^{' \infty} = [\begin{matrix} 1 & 0 & \dots & 0 \\ 1 & 0 & \dots & 0 \\ ⋮ & ⋮ & { \end{matrix} ⋮ 1 0 \dots 0]$ , and it is a stable random array.

The above formula shows that when k⟶∞, _pi,1 = 1, i = 1,2,...,N, that is, regardless of the initial state of each villager, it can finally transfer to the global optimal state 1 with probability 1. So there is the formula: $lim_{t \to \infty} p {f (X_{i}^{t}) \to f (X^{*})} = 1, i = 1, 2, \dots, N$

Therefore, PIDO algorithm has global convergence.

4 PIDO algorithm performance analysis

The hardware environment is used in the EBDO algorithm test experiment: The CPU is Intel Core™ I5, M520@ 2.40 GHz, the memory is 8 GB, and the operating system is Windows 10. The optimization termination condition has a maximum criterion and a norm criterion. In high-dimensional optimization problems, the norm criterion is used. The greater the adaptability, the better the solution. The method of judging whether the approximate global optimal solution has been obtained is the termination condition of the random algorithm. Within the range of the maximum number of iterations, one of the following conditions can be selected as the termination condition: (1) The maximum fitness value and the average fitness value change little and tend to be stable; (2) The distance between adjacent generation populations is less than the acceptable value. The PlatEMO test platform (https://github.com/BIMK/PlatEMO) is used for testing [30].

4.1 Analysis of the influence of the number of villagers on PIDO algorithm performance

The following takes the famous Bump optimization problem and the Michalewicz optimization problem as examples to analyze the performance of PIDO algorithm under different villagers’ number settings. These two optimization problems are very difficult to solve.

Bump optimization problem: ${\begin{matrix} min E_{0} (X) = \frac{| \sum_{i = 1}^{n} \cos^{4} (x_{i}) - 2 \prod_{i = 1}^{n} \cos^{2} (x_{i}) |}{\sqrt{\sum_{i = 1}^{n} {ix}_{i}^{2}}} \\ s . t . \prod_{i = 1}^{n} x_{i} ⩾ 0.75 \\ \sum_{i = 1}^{n} x_{i} ⩽ 7.5 n \\ 0 < x_{i} ⩽ 10 \end{matrix}$

Michalewicz optimization problem: $min E_{1} (X) = - \sum_{i = 1}^{n} sin (x_{i}) {[\sin (\frac{{ix}_{i}^{2}}{π})]}^{2 m}, (m = 15)$

The PIDO algorithm is used to solve the Bump optimization problem and the Michalewicz optimization problem. The calculation results are shown in Table 4. When calculating, G = 10⁷, ɛ= 10^–7, Q₀ = 0.9, L = 3, T = 3.

Table 4
Comparisons of calculation results under different parameter settings

n N Error Bump problem Michalewicz problem

objective function value Average calculation time/s objective function value Average calculation time/s

30 100 0.010 –0.805 921 803 –29.628 760 111

200 0.010 –0.811 803 3 099 –29.628 872 203

300 0.010 –0.808 967 2 145 –29.628 872 378

400 0.010 –0.813 265 5 885 –29.628 872 337

60 100 0.010 –0.828 497 1 121 –59.612 750 173

200 0.010 –0.825 695 3 064 –59.621 135 749

300 0.010 –0.824 975 3 167 –59.621 135 639

400 0.010 –0.826 308 4 353 –59.621 135 952

100 100 0.001 –0.835 829 2 002 –99.611 630 551

200 0.001 –0.835 526 2 593 –99.616 749 1 054

300 0.001 –0.833 598 5 294 –99.618 139 1 705

300 50 0.001 –0.843 502 2 623 –298.601 193 2 517

100 0.001 –0.832 157 4 826 –299.306 641 3 203

200 0.001 –0.848 444 9 993 –299.512 019 7 473

500 50 0.001 –0.825 474 4 200 –495.287 649 7 346

100 0.001 –0.834 220 7 739 –497.759 218 7 096

200 0.001 –0.823 085 13 096 –498.673 640 19 243

800 50 0.001 –0.725 698 8 090 –786.152 839 4 124

0.005 –0.783 280 6 328 –782.634 723 6 924

100 0.001 –0.730 634 14 779 –786.056 194 9 495

0.005 –0.796 418 14 352 –788.491 170 9 111

200 0.001 –0.722 817 38 722 –785.077 567 14 949

0.005 –0.817 508 28 829 –791.492 633 18 546

1000 50 0.001 –0.736 951 10 188 –983.373 683 10 143

0.005 –0.745 096 8 972 –988.308 603 5 805

100 0.001 –0.726 373 18 324 –976.184 758 12 885

0.005 –0.761 237 20 259 –977.959 383 14 192

200 0.001 –0.713 034 31 774 –976.270 965 18 256

0.005 –0.802 655 37 739 –983.158 225 22 457

n	N	Error	Bump problem	Michalewicz problem
30	100	0.010	–0.805 921	803	–29.628 760	111
	200	0.010	–0.811 803	3 099	–29.628 872	203
	300	0.010	–0.808 967	2 145	–29.628 872	378
	400	0.010	–0.813 265	5 885	–29.628 872	337
60	100	0.010	–0.828 497	1 121	–59.612 750	173
	200	0.010	–0.825 695	3 064	–59.621 135	749
	300	0.010	–0.824 975	3 167	–59.621 135	639
	400	0.010	–0.826 308	4 353	–59.621 135	952
100	100	0.001	–0.835 829	2 002	–99.611 630	551
	200	0.001	–0.835 526	2 593	–99.616 749	1 054
	300	0.001	–0.833 598	5 294	–99.618 139	1 705
300	50	0.001	–0.843 502	2 623	–298.601 193	2 517
	100	0.001	–0.832 157	4 826	–299.306 641	3 203
	200	0.001	–0.848 444	9 993	–299.512 019	7 473
500	50	0.001	–0.825 474	4 200	–495.287 649	7 346
	100	0.001	–0.834 220	7 739	–497.759 218	7 096
	200	0.001	–0.823 085	13 096	–498.673 640	19 243
800	50	0.001	–0.725 698	8 090	–786.152 839	4 124
		0.005	–0.783 280	6 328	–782.634 723	6 924
	100	0.001	–0.730 634	14 779	–786.056 194	9 495
		0.005	–0.796 418	14 352	–788.491 170	9 111
	200	0.001	–0.722 817	38 722	–785.077 567	14 949
		0.005	–0.817 508	28 829	–791.492 633	18 546
1000	50	0.001	–0.736 951	10 188	–983.373 683	10 143
		0.005	–0.745 096	8 972	–988.308 603	5 805
	100	0.001	–0.726 373	18 324	–976.184 758	12 885
		0.005	–0.761 237	20 259	–977.959 383	14 192
	200	0.001	–0.713 034	31 774	–976.270 965	18 256
		0.005	–0.802 655	37 739	–983.158 225	22 457

It can be seen from Table 4:

When N increases, the CPU time consumed greatly increases.

For a given n, if N increases, the CPU consumption time also increases greatly, but the accuracy of the average optimal objective function value increases at the same time.

For a given n and N, if E0 increases, the accuracy of the average optimal target function value also increases, but the CPU consumed time may increase or decrease.

Therefore, if n > 500, N = 100∼200 is sufficient; if n < 500, N = 200 is sufficient.

4.2 Comparison with other swarm intelligence algorithms

In this paper, six benchmark function optimization problems provided by the internationally universal intelligent optimization algorithm test package CEC2019 [31] are used to test the performance of the PIDO algorithm, as shown in Table 5. The test includes 28 well-designed benchmark function optimization problems.

Table 5
6 benchmark function optimization problems

Benchmark function number Benchmark function name Theoretical global optimal solution Global optimal objective function value

F14 Schwefel’s Function unknown –3 721.476 4

F17 Lunacek bi-Ratrigin Function O 300.000 0

F19 Rotated Expanded Griewank’s plus O 500.000 0

F21 Rosenbrock’s Function O 900.000 0

F25 Composition Function 1 unknown 1 200.000 0

F28 Composition Function 2 unknown –4 756.706 8

Search range [–100, 100]ⁿ

Benchmark function number	Benchmark function name	Theoretical global optimal solution	Global optimal objective function value
F14	Schwefel’s Function	unknown	–3 721.476 4
F17	Lunacek bi-Ratrigin Function	O	300.000 0
F19	Rotated Expanded Griewank’s plus	O	500.000 0
F21	Rosenbrock’s Function	O	900.000 0
F25	Composition Function 1	unknown	1 200.000 0
F28	Composition Function 2	unknown	–4 756.706 8
Search range [–100, 100]ⁿ

In Table 5, O is the theoretical global optimal solution that can be set arbitrarily. The theoretical global optimal solutions of F14, F25, and F28 have not yet been found. The global optimal objective function value in Table 5 is the best objective function value found by the algorithm so far. In this paper, PIDO algorithm is used to solve the 6 function optimization problems in Table 5. The parameters are N = 200, n = 50, G = 10⁷, ɛ=10^–7, Q₀ = 0.9, E₀ = 0.01, L = 3, T = 3. The seven intelligent optimization algorithms compared with PIDO algorithm are all selected from the algorithms published in recent international journals. These algorithms are shown in Table 6, namely RCGA (real-coded genetic algorithm), DASA (differential ant-stigmergy algorithm), NP -PSO (non-parametric particle swarm optimization), MpBBO (metropolis biogeography-based optimization), MDE-LiGO (modified differential evolution with locality induced genetic operators), SLADE (symmetric Latin-based adaptive differential evolution), ABC (artificial bee colony) algorithm). The termination conditions for these seven algorithms are: evolutionary algebra G = 10⁷ or optimal solution error ɛ= 10^–7.

Table 6

Parameters setting of compared intelligent optimization algorithms

Intelligent optimization algorithm	parameter settings
RCGA^[3]	k_max = 100, pr = 0.5, φ ₀ = 0.5,λ= 0.5, C_I = 200
DASA^[4]	b = 10; m = 38;ρ= 0.24; s + = 0.025, s_– = 0.015;ɛ= 1.0E-16
NP-PSO^[6]	N = 200
MpBBO^[9]	k = 100,elit = 4,m_max = 0.01,R_m = round(0.5k,k), T = (LCR,ECR,ICR,ILCR,LCR),E = I = 1
MDE-LiGO^[10]	N = 200; p = 10, SU = 5, c = 0.1,μ _Cr = 0.5,μ _F = 0.5,μ _Br = 0.5
SLADE^[11]	N = 200;θ _CR = 0.5, a = 0.1,μ _F = 0.5
ABC^[12]	Bees = 200, Trys = 200n, q = 0.15

The meaning of each parameter in Table 6 can be found in related literature. Aiming at the 6 optimization problems listed in Table 5, PIDO algorithm and the 7 algorithms listed in Table 6 are used to solve them, and each algorithm is solved 51 times independently. Table 7 shows the average, standard deviation, median, minimum, maximum, and fitness evaluation times of the optimal objective function values obtained by the eight algorithms.

Table 7

Optimal solutions obtained by 8 compared algorithms

Benchmark function	Evaluation index	PIDO	RCGA	DASA	NP-PSO	MpBBO	MDE-LiGO	SLADE	ABC
F14	Average	3.72E+03	–3.72E+03	–3.72E+03	–3.72E+03	–8.37E+02	–3.72E+03	–3.72E+03	–3.72E+03
	Median	–3.72E+03	–3.72E+03	–3.72E+03	–3.72E+03	–8.37E+02	–3.72E+03	–3.72E+03	–3.72E+03
	Standard deviation	1.07E–05	6.90E–06	3.06E–06	3.17E–06	3.19E–05	3.05E–06	3.32E–06	2.83E–06
	Minimum	–3.72E+03	–3.72E+03	–3.72E+03	–3.72E+03	–8.37E+02	–3.72E+03	–3.72E+03	–3.72E+03
	Max	–3.72E+03	–3.72E+03	–3.72E+03	–3.72E+03	–8.37E+02	–3.72E+03	–3.72E+03	–3.72E+03
	Evaluation times	8.79E+03	1.22E+04	1.66E+06	1.14E+05	8.76E+05	1.12E+05	2.04E+05	1.68E+05
	Rank 1	1	2	4	7	8	6	5	3
	Rank 2	1	2	4	7	8	6	5	3
F17	Average	3.00E+02	3.00E+02	3.00E+02	3.00E+02	3.02E+02	3.00E+02	3.00E+02	3.00E+02
	Median	3.00E+02	3.00E+02	3.00E+02	3.00E+02	3.02E+02	3.00E+02	3.00E+02	3.00E+02
	Standard deviation	0.00E+00	1.79E–03	2.78E–03	2.07E–03	3.03E–03	2.99E–03	3.08E–03	0.00E+00
	Minimum	3.00E+02	3.00E+02	3.00E+02	3.00E+02	3.02E+02	3.00E+02	3.00E+02	3.00E+02
	Max	3.00E+02	3.00E+02	3.00E+02	3.00E+02	3.02E+02	3.00E+02	3.00E+02	3.00E+02
	Evaluation times	1.00E+03	1.16E+04	8.49E+05	2.45E+05	7.10E+05	3.37E+05	9.55E+05	2.98E+02
	Rank 1	1	5	4	3	8	7	6	1
	Rank 2	2	5	4	3	8	7	6	1
F19	Average	5.12E+02	9.18E+02	5.53E+02	1.13E+04	5.38E+02	5.02E+02	5.02E+02	5.28E+02
	Median	5.12E+02	9.18E+02	5.53E+02	1.13E+04	5.38E+02	5.02E+02	5.02E+02	5.28E+02
	Standard deviation	2.93E–03	3.38E–03	3.15E–03	2.85E+02	2.84E–03	2.96E–03	3.27E–03	3.14E–03
	Minimum	5.12E+02	9.18E+02	5.53E+02	1.06E+04	5.38E+02	5.02E+02	5.02E+02	5.28E+02
	Max	5.12E+02	9.18E+02	5.53E+02	1.18E+04	5.38E+02	5.02E+02	5.02E+02	5.28E+02
	Evaluation times	6.80E+07	8.26E+07	1.12E+07	2.13E+07	1.43E+07	1.42E+06	1.01E+07	6.70E+06
	Rank 1	3	7	6	8	5	2	1	4
	Rank 2	3	7	6	8	5	2	1	4
F21	Average	9.00E+02	6.14E+03	9.14E+02	9.07E+02	9.21E+02	9.00E+02	9.00E+02	9.04E+02
	Median	9.00E+02	6.14E+03	9.14E+02	9.07E+02	9.21E+02	9.00E+02	9.00E+02	9.04E+02
	Standard deviation	0.00E+00	2.72E–02	3.42E–03	2.60E–03	3.18E–03	0.00E+00	0.00E+00	2.98E–03
	Minimum	9.00E+02	6.14E+03	9.14E+02	9.07E+02	9.21E+02	9.00E+02	9.00E+02	9.04E+02
	Max	9.00E+02	6.14E+03	9.14E+02	9.07E+02	9.21E+02	9.00E+02	9.00E+02	9.04E+02
	Evaluation times	7.44E+05	9.80E+06	4.27E+06	7.63E+06	9.43E+05	3.80E+05	4.79E+05	1.03E+07
	Rank 1	1	8	6	5	7	1	1	4
	Rank 2	3	8	6	5	7	1	2	4
F25	Average	1.25E+03	1.27E+03	1.28E+03	1.26E+03	1.29E+03	1.26E+03	1.28E+03	1.27E+03
	Median	1.25E+03	1.27E+03	1.28E+03	1.26E+03	1.29E+03	1.26E+03	1.28E+03	1.27E+03
	Standard deviation	2.89E–02	3.28E–02	2.67E–02	3.14E–02	3.11E–02	2.66E–02	3.30E–02	3.32E–02
	Minimum	1.25E+03	1.27E+03	1.28E+03	1.26E+03	1.29E+03	1.26E+03	1.28E+03	1.27E+03
	Max	1.25E+03	1.27E+03	1.28E+03	1.26E+03	1.29E+03	1.26E+03	1.28E+03	1.27E+03
	Evaluation times	2.05E+06	2.47E+05	2.92E+03	3.34E+05	3.62E+05	4.41E+05	2.44E+05	7.93E+05
	Rank 1	1	5	7	3	8	2	6	4
	Rank 2	1	5	7	3	8	2	6	4
F28	Average	–4.75E+03	7.84E+03	3.97E+03	–2.49E+03	4.02E+03	3.06E+03	3.03E+03	4.12E+03
	Median	–4.75E+03	7.84E+03	3.97E+03	–2.49E+03	4.02E+03	3.06E+03	3.03E+03	4.12E+03
	Standard deviation	2.87E–02	3.21E–02	3.25E–02	3.06E–05	2.61E–02	3.01E–02	3.42E–02	3.13E–02
	Minimum	–4.75E+03	7.84E+03	3.97E+03	–2.48E+03	4.02E+03	3.06E+03	3.03E+03	4.12E+03
	Max	–4.75E+03	7.84E+03	3.97E+03	–2.48E+03	4.02E+03	3.06E+03	3.03E+03	4.12E+03
	Evaluation times	3.98E+06	3.15E+06	1.44E+06	2.72E+06	1.27E+06	2.02E+06	1.54E+06	9.17E+05
	Rank 1	1	8	5	2	6	4	3	7
	Rank 2	1	8	5	2	6	4	3	7
Overall ranking 1		8	35	32	28	42	22	22	23
Overall ranking 2		11	35	32	28	42	22	23	23
Final overall ranking 1		1	7	6	5	8	2	2	4
Final overall ranking 2		1	7	6	5	8	2	3	3

Note: The color word indicates that the global theoretical optimal solution has been reached.

Overall ranking 1 is the ranking of these algorithms based on the average optimal objective function value, and overall ranking 2 is the ranking of these algorithms based on the average optimal objective function value and the number of fitness evaluations. Final overall ranking 1 and final overall ranking 2 are based on rankings based on overall ranking 1 and overall ranking 2, respectively.

The non-parametric Friedman test [32, 33] is a non-parametric test based on the results obtained from the PIDO algorithm and the results obtained from the seven compared algorithms. The results of the Friedman test are shown in Table 8. Significance = 1 means that PIDO algorithm performance is 99% statistically different from the compared algorithm, and significance = 0 means that there is no significant difference. In Table 8, the number of cases with significance = 1 and the number of cases with significance = 0 respectively represent the number of the PIDO algorithm and the seven compared algorithms that are significantly different and not significantly different from the number of benchmark function optimization problems.

Table 8

Results comparison of Friedman test (α= 0.01)

Benchmark optimization problem	PIDO algorithm vs.
	RCGA	DASA	NP-PSO	MpBBO	MDE-LiGO	SLADE	ABC
F14	0.001	0	0	0	0	0	0
Significance	1	1	1	1	1	1	1
F17	0	0	0	0	0	0	1
Significance	1	1	1	1	1	1	0
F19	0	0	0	0	0	0	0
Significance	1	1	1	1	1	1	1
F21	0	0	0	0	1	1	0
Significance	1	1	1	1	0	0	1
F25	0	0	0	0	0	0	0
Significance	1	1	1	1	1	1	1
F28	0	0	0	0	0	0	0
Significance	1	1	1	1	1	1	1
Cases number of	6	6	6	6	5	5	5
significance = 1
Cases number of	0	0	0	0	1	1	1
significance = 0

It can be seen from Table 7 that PIDO, RCGA, DASA, NP-PSO, MpBBO, MDE-LiGO, SLADE, ABC are sorted according to the accuracy of the average optimal objective function value and the results are as follows: PIDO > SLADE = MDE-LiGO > ABC > NP-PSO > DASA > RCGA > MpBBO

The results of sorting according to the average optimal objective function value and the number of fitness evaluations are as follows: PIDO > MDE-LiGO > SLADE > ABC > NP-PSO > DASA > RCGA > MpBBO

It can be seen from Table 8 that the total number of significant cases for the PIDO algorithm to solve the 6 benchmark function optimization problems is 39, which is significantly larger than the total number of insignificant cases of 3, indicating that the performance of PIDO algorithm is significantly better than that of the seven compared algorithms.

Figure 3 (a)∼Fig. 3 (f) illustrate the sample convergence curve of PIDO, RCGA, DASA, NPPSO, MpBBO, MDE-LiGO, SLADE, ABC algorithm to solve the 6 benchmark function optimization problems in Table 5, where the level and the vertical axis adopts a logarithmic scale. It can be seen from Fig. 3 and Table 7 that when PIDO algorithm solves the 6 benchmark function optimization problems, it can find the best global optimal solution. From a comprehensive point of view, the comprehensive performance of PIDO algorithm is better than that of the seven compared algorithms, indicating that it has high solution accuracy and fast calculation speed.

Fig. 3

Convergence curves of samples.

4.3 Multi-objective engineering problem test

Three real-world engineering design optimization problems ie, a 3-objective car impact (CSI) problem [34], a 5-objective water resource management (WRM) problem [35] and a 10-objective general avaition aircraft (GAA) problem [36] are solved.

Inverted generational distance (IGD) [37] is used as the metrics for benchmarking. The reference sets for IGD calculation are obtained from the PlatEMO framework. For the engineering design problems, the reference sets have been taken from [38]. Performance of the proposed algorithm is compared with cK-RVEA1 (the best performing approach reported in [39]) on constrained engineering design problems. The statistics of the IGD metric across 25 independent runs are presented in Table 9. Here, ↑, ↓ and ≈ are indicates whether the proposed algorithm is statistically better, worse and equivalent to cK-RVEA1. From Table 9, it can be observed that the proposed method comprehensively outperforms cK-RVEA1 in the problem instances.

Table 9
Test results for IGD metric obtained by the proposed method and cK-RVEA1 for constrained engineering problems

Prooblem M proposed method cK-RVEA1

Min Mean Max Min Mean Max

CSI 3 0.2734 0.2922 0.3233 ↑ 0.3337 0.4724 0.6424

WRM 5 24973.8 28969.3 39621.3 ↑ 47806.5 83425.4 123370.7

GRAA 10 25.496 49,591 121.504 ↑ 53.920 113.828 264.735

Prooblem	M	proposed method	cK-RVEA1
CSI	3	0.2734	0.2922	0.3233	↑	0.3337	0.4724	0.6424
WRM	5	24973.8	28969.3	39621.3	↑	47806.5	83425.4	123370.7
GRAA	10	25.496	49,591	121.504	↑	53.920	113.828	264.735

4.4 Analysis of local optimization ability and global optimization ability

In this paper, the method in [40] is used to analyze the balance between the local and global optimization capabilities of PIDO algorithm. The penetration behavior is used to illustrate the global optimization capability of the PIDO algorithm, and the expansion behavior is used to illustrate the local optimization capability of the PIDO algorithm.

Figure 4 shows the coordination relationship between penetration and expansion behaviors in the time intervals [0,694] and [0,966] when the PIDO algorithm solves the benchmark optimization problem F3 [31]. It can be seen from Fig. 4 that the penetration behavior is dominant from 0 s to 120 s, and the expansion behavior is dominant after 120 s.

Fig. 4

Coordination of penetration and expansion behavior.

Figure 5 illustrates the coordination between the expansion and sticky behaviors in the time intervals [0,694] and [0,966] when the PIDO algorithm solves F3. It can be seen from Fig. 5 that the expansion behavior is always dominant, and the sticky behavior is not dominant at all, but the dominance of the expansion behavior decreases with the passage of time; on the contrary, the number of occurrences of the sticky behavior increases with time.

Fig. 5

Coordination of expansion and stickiness behavior.

5 Conclusion

In order to solve some optimization problems with multi-local optimal solutions, a plague infectious disease optimization (PIDO) algorithm is proposed by using a dynamic model of plague epidemic with impulse vaccination and time delay. The approach includes constraint handling which is often overlooked in contemporary algorithms. The performance of the approach is benchmarked extensively on a set of unconstrained and constrained problems with regular and irregular PFs (Pareto fronts). A statistical comparison with the existing techniques highlights the efficacy and potential of the approach. Pareto frontier analysis, the pareto frontier is the solution closest to the true value in the feasible region of the pareto. The solution requires: (1) The distribution is as uniform as possible; (2) The more within a certain range, the better. There are two ways to solve the feasible solution of pareto: (1) According to the importance of the constraint function, it is weighted to synthesize a single function (sometimes with a penalty factor) to solve it. (2) Using some common evolutionary algorithms (EA) to solve the problem, the current common There are GA, PSO, etc. The algorithm idea here has Pareto optimal characteristics. It refers to an ideal state of resource allocation. Assuming that an inherent group of people and allocable resources change from one state of distribution to another, without making anyone worse off, at least a person becomes better, this is Pareto improvement or Pareto optimization. Improving economic efficiency means reducing waste. If no one in the economy can make his own situation better without making others worse off, then this state has reached the optimal allocation of resources. The efficiency defined in this way is called the Pareto optimal efficiency. In other words, if a person can improve his own situation without harming the interests of others, he has achieved a Pareto improvement in resource allocation.

There are the following characteristics in PIDO:

PIDO algorithm includes 9 operators such as S_S, S_E, E_E, E_I, E_R, I_I, I_R, R_R, R_S, etc. These operators can significantly increase the search capability of the algorithm.

The random method is used to select the parameters of the plague epidemic model in PIDO algorithm and the relevant parameters in each operator, which can not only make the model better express the actual situation, but also greatly reduce the number of input parameters.

The operators S_S, S_E, E_E, E_I, E_R, I_I, I_R, R_R, R_S in the PIDO algorithm are constructed by using the plague epidemic model, and they do not need to be related to the problem to be solved, so the PIDO algorithm is universal.

In PIDO algorithm, the operators S_S, E_E, I_I, and R_R can enable villagers with high HHI index to transmit strong feature information to villagers with low HHI index, so that the latter can develop in a good direction. Operators S_E, E_I, E_R, I_R and R_S not only enable villagers in different states to exchange information, but also enable villagers to obtain the processed characteristic information of other villagers, reducing the probability of villagers falling into local optimal traps, thereby increasing PIDO Algorithm’s local exploitation ability. Pulse vaccination has the characteristic of making villagers jump out of the local optimal trap. Therefore, PIDO algorithm can realize information exchange between villagers from multiple angles, which is meaningful to improve the exploration ability of PIDO algorithm.

When the PIDO algorithm is used for evolutionary calculation, since only a small part of the villagers’ features are attacked by the plague virus, each time the villagers exchange information, only a small part of the features participate in the calculation, so when solving the optimization of higher dimensionality, the convergence speed of PIDO algorithm can be significantly improved, and PIDO algorithm is suitable for solving optimization problems with higher dimensions.

The evolution process of PIDO algorithm has Markov characteristics and “not bad at every step” characteristics, so as to ensure that the algorithm has global convergence.

The future improvement directions of the PIDO algorithm are as follows:

The plague infectious disease dynamic model is used to optimize the relevant parameters of PIDO algorithm to make these parameter settings more reasonable.

In-depth study of the dynamic characteristics of operators such as S_S, S_E, E_E, E_I, E_R, I_I, I_R, R_R, R_S, etc.

In-depth study of the dynamic characteristics of the villagers in the process of PIDO algorithm solving.

Footnotes

Acknowledgments

This work was supported by the Scientific Research Project (NO.20B337) of Hunan Provincial Education Department, China & Hunan Natural Science Foundation (Grant No: 2019JJ40154), Hunan Province, China.

References

Lin

S.J.

, Dong

, Chen

M.Z.

, et al., Summary of new group intelligent optimization algorithms, Applications 54(12) (2018), 1–9.

Wolpert

D.H.

and Macready

W.G.

, No free lunch theorems for optimization, IEEE Transactions on Evolutionary Computation 1(1) (1997), 67–82.

Chuang

Y.C.

, Chen

C.T.

and wang

C.H.

, A simple and efficient real-coded genetic algorithm for constrained optimization, Applied Soft Computing 38 (2016), 87–105.

Korošec

, Šilc

and Filipic

, The differential ant-stigmergy algorithm, Information Sciences 192 (2012), 82–97.

Liu

Z.Q.

, You

X.M.

and Liu

, Ant colony algorithm for heuristic dynamic pheromone update strategy, Computer Engineering and Applications 54(20) (2018), 20–27.

Beheshti

and Shamsuddin

S.M.

, Non-parametric particle swarm optimization for global optimization, Applied Soft Computing 28 (2015), 345–359.

Yuan

, Ge

H.W.

and Jiang

D.Y.

, Partical swarm optimization algorithm with adaptive filter based on health degree, Journal of Frontiers of Computer Science and Technology 12(2) (2018), 332–340.

Aburomman

A.A.

and Reaz

M.B.I.

, A novel SVM-kNN-PSO ensemble method for intrusion detection system, Applied Soft Computing 38(C) (2016), 360–372.

Al-Roomi

A.R.

and El-Hawary

M.E.M.E.

, Metropolis biogeographybased optimization, Information Sciences 360 (2016), 73–95.

10.

Mukherjee

, Debchoudhury

and Das

, Modified differential evolution with locality induced genetic operators for dynamic optimization, European Journal of Operational Research 253(2) (2016), 337–355.

11.

Zhao

Z.W.

, Yang

J.M.

, Hu

Z.Y.

, et al., A differential evolution algorithm with self-adaptive strategy and control parameters based on symmetric Latin hypercube design for unconstrained optimization problems, European Journal of Operational Research 250(1) (2016), 30–45.

12.

Mernik

, Liu

S.H.

, Karaboga

, et al., On clarifying misconceptions when comparing variants of the artificial bee colony algorithm by offering a new implementation, Information Sciences 291 (2015), 115–127.

13.

Khan

and Das

, Multiuser Detection for MIMO-OFDM system in Underwater Communica-tion Using a Hybrid Bionic Binary Spotted Hyena Optimizer, Journal of Bionic Engineering 18(2) (2021), 462–472. doi: 10.1007/s42235-021-0018-y

14.

Lei

X.S.

, Zhang

X.Y.

and Hao

Y.K.

, An Intelligent Ellipsoid Calibration Method Based on the Grey Wolf Algorithm for Magnetic Compass, Journal of Bionic Engineering 18(2) (2021), 453–461. doi: 10.1007/s42235-021-0033-z

15.

Azqandi

M.S.

, Delavar

and Arjmand

, An enhanced time evolutionary optimization for solving engineering design problems, Engineering with Computers 36(2) (2020), 763–781. doi: 10.1007/s003 66-019-00729-w

16.

Mal

and Bhattacharjee

, Optimization of Chaboche kinematic hardening parameters for 20MnMoNi55 reactor pressure vessel steel by sequenced genetic algorithms maintaining the hierarchy of dependence, Engineering Optimization 53(2) (2021), 335–347. doi: 10.1080/030521 5X.2020.1726340.

17.

Habib

, Singh

H.K.

, et al., A Multiple Surrogate Assisted Decomposition Based Evolutionary Algorithm for Expensive Multi/Many-Objective Optimization, IEEE Transactions on Evolutionary Computation 23(06) (2019), 1000–1014. doi: 10.1109/TEVC.2019.2899030

18.

Chen

L.S.

, Meng

X.Z.

and Jiao

J.J.

, Biodynamics. Beijing: Science Press, (2009), China.

19.

Guo

Z.Y.

, Xiao

, Li

H.T.

, et al., Study on construction of plague transmission dynamics model, Progress in Modern Biomedicine 12(20) (2012), 3838–3841.

20.

Grassly

N.C.

and Fraser

, Mathematical models of infectious disease transmission, Nature Reviews Microbiology 42(2) (2008), 3181–3195.

21.

Monecke

, Monecke

and Monecke

, Modelling the black death: a historical case study and implications for the epidemiology of bubonic plague, International Journal of Medical Microbiology 45(5) (2009), 2473–2486.

22.

Bjorn

P.Z.

and Hartmut

, The history of the plague and the research on the causative agent Yersinia pestis, International Journal of Hygiene and Environmental Health 39(7) (2004), 2135–2145.

23.

Yan

, Liu

G.C.

, Lan

X.Y.

, et al., Effects of meteorological factors on population density of meriones unguiculatus in plague foci of Hebei Province, Chinese Journal of Hygienic Insecticides & Equipments 23(2) (2017), 159–161.

24.

, Li

G.C.

, Si

X.Y.

, et al., Nonlinear effects of climate driven plague in meriones unguiculatus natural foci in Inner Mongolia, Chinese Journal of Vector Biology and Control 27(4) (2016), 321–325.

25.

Liao

L.F.

, He

Z.M.

, Hou

Y.Y.

, et al., Biological characteristics and application on cricetus migratorius, Laboratory Animal Science 33(3) (2016), 46–51.

26.

Y.J.

and Yang

R.F.

, Animal model of plague: classics and questions, Acta Laboratorium Animalis Scientia Sinica 19(2) (2011), 161–171.

27.

Zhou

X.L.

, Liu

Z.C.

, Cong

X.B.

, et al., Progress in prediction and early warning of plague, Chinese Journal of Control of Endemic Diseases 27(5) (2012), 336–338.

28.

Liu

Z.C.

, Zhou

X.L.

, Zhang

B.Y.

, et al., Establishment of predictive model and warning index of animal plague, Chinese Journal of Control of Endemic Diseases 30(1) (2015), 1–3.

29.

Iosifescu

, Finite Markov processes and their applications, NewYork: Dover Publications, (1980).

30.

Tian

, Cheng

, Zhang

X.Y.

and Jin

Y.C.

, PlatEMO: A MATLAB platform for evolutionary multi-objective optimization, IEEE Computational Intelligence Magazine 12(4) (2017), 73–87.

31.

Liang

J.J.

, Qu

B.Y.

, et al., Problem Definitions and Evaluation Criteria for the CEC 2019 Special Session on Multimodal Multiobjective Optimization, China: Zhengzhou University (2019). doi:10.13140/RG.2.2.33423. 64164. http://www5.zzu.edu.cn/cilab/info/1005/1300.htm.

32.

Črepinšek

, Liu

S.H.

and Mernik

, Replication and comon of computational experiments in applied evolutionary computing: common pitfalls and guidelines to avoid them, Applied Soft Computing 19 (2014), 161–170 paris.

33.

Derrac

, García

, Molina

, et al., A practical tutorial on the use of non-parametric statistical tests as a methodology for comparing evolutionary and swarm intelligence algorithms, Swarm and Evolutionary Computation 1 (2011), 3–18.

34.

Jain

and Deb

, An evolutionary many-objective optimization algorithm using reference-point based non-dominated sorting approach part II: Handling constraints and extending to an adaptive approach, IEEE Trans Evol Comput, 18(4) (2014), 602–622.

35.

Ray

, Tai

and Seow

K.C.

, An evolutionary algorithm for multi-objective optimization, Eng Optim, 33(3) (2001), 399–424.

36.

Simpson

T.W.

, Chen

, Allen

J.K.

and Mistree

, Conceptual design of a family of products through the use of the robust concept exploration method, Proc. AIAA/USAF/NASA/ISSMO Symp. Multidiscipl. Anal. Optim., (1996), 1535–1545.

37.

Bosman

P.A.N.

and Thierens

, The balance between proximity and diversity in multiobjective evolutionary algorithms, IEEE Trans. Evol. Comput. 7(2) (2003), 174–188.

38.

Bhattacharjee

K.S.

, Singh

H.K.

and Ray

, A novel decomposition-based evolutionary algorithm for engineering design optimization, J Mech Design 139(4) (2017), 1–13.

39.

Chugh

, Sindhya

, Miettinen

, Hakanen

and Jin

, On constraint handling in surrogate-assisted evolutionary many-objective optimization, Proc. Parallel Problem Solving Nat. (2016), 214–224.

40.

Ding

and Lu

, A microbial kinetic optimization approach, Journal of Intelligent & Fuzzy Systems 40(3) (2021), 5697–5713. doi:10.3233/JIFS-201828