An aircraft system manufacturing process optimization method based on improved SPEA2 algorithm

Abstract

The application discloses an aircraft system manufacturing process optimization method based on an improved SPEA2 algorithm, and relates to the field of aviation part manufacturing. First, a manufacturing process mathematical model is constructed; second, the mathematical model is optimized by using the improved SPEA2 algorithm; the improved SPEA2 algorithm is a neighborhood search method for increasing solution region density on the basis of the SPEA2 algorithm, and can organically combine global search with local search; finally, the aircraft system manufacturing process with the shortest manufacturing period and the lowest resource consumption is obtained through decoding. The application maps the aircraft system manufacturing process corresponding to the mathematical model into the SPEA2 algorithm through encoding and decoding, realizes optimization of the aircraft system manufacturing process, saves aircraft system manufacturing time, and reduces manufacturing cost; the improved SPEA2 algorithm effectively ensures the diversity of the population, and can accelerate the convergence speed of the SPEA2 algorithm.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of aircraft component manufacturing, and particularly relates to an aircraft system manufacturing process optimization method based on an improved SPEA2 algorithm. BACKGROUND

[0002] At present, it is still difficult to plan the large number of parallel and staggered manufacturing processes objectively existing in the manufacturing process of the aircraft in a more scientific and reasonable manner, and the aircraft system manufacturing process often has typical problems such as long cycle and low efficiency. On the one hand, due to the strict requirements of aircraft integration verification, the sequence of manufacturing of each component has a complex coupling relationship, which is manifested as complex timing constraints between manufacturing processes. On the other hand, even for the same aircraft component manufacturing, different process planning will lead to differences in types and quantities of manufacturing resources. At the same time, the use of various types of manufacturing resources is often subject to various limitations of the actual manufacturing process. Therefore, before the aircraft is formally put into manufacturing, it is necessary to deeply analyze the logical relationship between the manufacturing of each component and to clearly define the actual constraints on the consumption of manufacturing resources under different process planning, which is the key to solving the complex problems of coupling relationship and actual constraints in the current aircraft system manufacturing process.

[0003] Among various algorithms, SPEA2 (Strength Pareto Evolutionary Algorithm 2) is a classic multi-objective evolutionary algorithm. The core of SPEA2 lies in its unique fitness function, which combines the degree to which an individual is dominated and its ability to dominate other individuals to evaluate its fitness. In addition, the algorithm uses a fixed-size external archive to save the current best non-dominated solution, and continuously generates new solutions through selection, crossover and mutation operations. SPEA2 has the advantages of strong search ability, balanced convergence and diversity, easy implementation and tuning, etc., and has good effect when facing low-complexity parallel process optimization problems,

[0004] However, when dealing with complex aircraft system manufacturing process optimization problems, there are problems such as slow convergence speed of the algorithm, poor diversity of population distribution, and easy to fall into local optimum. The improved SPEA2 algorithm can effectively make up for the shortcomings of the traditional SPEA2 algorithm in dealing with aircraft system manufacturing process optimization problems by introducing a neighborhood search method based on solution region density, so as to adapt to the actual needs of aircraft system manufacturing process optimization. SUMMARY

[0005] In order to solve the above technical problems, the purpose of the present application is to provide an aircraft system manufacturing process optimization method based on improved SPEA2 algorithm, first, a manufacturing process mathematical model is constructed, second, the improved SPEA2 algorithm is used to optimize the mathematical model, the improved SPEA2 algorithm is a neighborhood search method for increasing solution region density based on the SPEA2 algorithm, which can combine global search and local search organically, finally, the shortest manufacturing cycle and the lowest resource consumption of the aircraft system manufacturing process are obtained through decoding. Figure 1

[0006] The technical solution process of the present application is implemented according to the following steps:

[0007] Step one, constructing a mathematical model of the aircraft system manufacturing process:

[0008] In order to optimize the manufacturing cycle and resource consumption of the aircraft manufacturing process, when establishing the mathematical model, the time sequence constraint, the resource constraint and the interaction and influence of different manufacturing methods in the same manufacturing process on the aircraft manufacturing process need to be considered.

[0009] S1.1, constructing a mathematical model of time sequence constraint:

[0010] Firstly, the time sequence constraint relationship is sorted out:

[0011] In order to clearly describe various constraint conditions and accurately express the optimization objective of the shortest manufacturing cycle and the lowest resource consumption, the aircraft manufacturing process includes J manufacturing processes, wherein, process i and process j are any two manufacturing processes in the J manufacturing processes; the start time of process i is S i , the end time is F i , the duration of process i is d i ; the start time of process j is S j , the end time is F j , the duration of process j is d j , the weight between process i and process j is l ij , and the weight between process j and process i is l ji ;

[0012] ​The generalized precedence relations of the timing constraints of the aircraft system manufacturing process are divided into four types: Start-Start (S-S), Finish-Start (F-S), Start-Finish (S-F) and Finish-Finish (F-F), and the four types of process relations are accompanied by minimum time lags and maximum time lags.

[0013] As Figure 2 shown:

[0014] In Figure 2 (a) S-S relations, the directed connection is drawn from the left side of process i and introduced from the left side of process j, and the number pair on the directed connection indicates that process j should start within to this time after process i starts, where are the minimum time lag and the maximum time lag under S-S relations, respectively;

[0015] In Figure 2 (b) F-S relations, the directed connection is drawn from the right side of process i and introduced from the left side of process j, and the number pair on the directed connection indicates that process j should start within to this time after process i ends, where are the minimum time lag and the maximum time lag under F-S relations, respectively;

[0016] In Figure 2 (c) S-F relations, the directed connection is drawn from the left side of process i and introduced from the right side of process j, and the number pair on the directed connection indicates that process j should end within to this time after process i starts, where are the minimum time lag and the maximum time lag under S-F relations, respectively;

[0017] In Figure 2 (d) F-F relations, the directed connection is drawn from the right side of process i and introduced from the right side of process j, and the number pair on the directed connection indicates that process j should end within to this time after process i ends, where are the minimum time lag and the maximum time lag under F-F relations, respectively;

[0018] Then, the timing constraint model is determined: the end time of the process is represented by the start time plus the manufacturing period, and the generalized precedence relations between process i and process j are represented and uniformly transformed into formulas (1) to (8).

[0019] For SS relationship, the start time of process j is within the start time of process i. back, Before, therefore need to meet The relationship between process i and process j is transformed into the following formula:

[0020]

[0021] For the FS relationship, the start time of process j is within the start time of process i. after, Before, therefore need to meet The relationship between process i and process j is transformed into the following formula:

[0022]

[0023] For SF relationship, the start time of process j is within the start time of process i. after, Before, therefore need to meet The relationship between process i and process j is transformed into the following formula:

[0024]

[0025] For FF relationship, the start time of process j is within the start time of process i. after, Before, therefore need to meet The relationship between process i and process j is transformed into the following formula:

[0026]

[0027] S1.2, Construct a mathematical model of manufacturing method constraints:

[0028] Among the J processes in the aircraft manufacturing process, any process j has M j The selection method of the manufacturing method of the aircraft manufacturing process is shown in formula (9);

[0029]

[0030] Among them, x jm represents the manufacturing method selected by process j, x jm When the value is 1, it means that process j selects manufacturing method m; jm When the value is 0, it means that process j does not select manufacturing method m;

[0031] Because M jEach manufacturing process can only choose one manufacturing method, and cannot adopt multiple manufacturing methods at the same time, so x jm The constraints that need to be met are shown in formula (10):

[0032]

[0033] S1.3, Constructing a mathematical model of resource constraints:

[0034] For updateable resources: for The set of processes being executed at any given moment. The completion of the aircraft system manufacturing process requires consumption of three renewable resources. When process j is executed under manufacturing mode m, the number of processes per day is The demand for this resource is o jmk , No. The daily supply of renewable resources is Since the use of renewable resources at any time cannot exceed the maximum supply, the mathematical model of the aircraft system manufacturing process must meet the constraints shown in formula (11):

[0035]

[0036] For non-renewable resources: The aircraft system manufacturing process needs to consume W non-renewable resources to complete. When the jth process is executed under the manufacturing mode m, the jth process consumes W non-renewable resources every day. w The demand for non-renewable resources is o jmw , No. w The supply of non-renewable resources is O w Since the use of non-renewable resources cannot exceed the maximum supply during the entire process of aircraft system manufacturing execution, the constraint conditions shown in formula (12) must be met:

[0037]

[0038] When the objective function f1 conforms to the following formula (13), the optimized aircraft manufacturing cycle is the shortest;

[0039] f1=min max(F j )(j=1,2,…,J) (13)

[0040] When the aircraft manufacturing cycle is the shortest, and the objective function f2 conforms to the following formula (14), the resource consumption in the aircraft manufacturing process is the lowest;

[0041]

[0042] Among them, the unit price of the kth renewable resource is C k , the unit price of the lth non-renewable resource is Cl ;

[0043] S1.4. Finally, while comprehensively considering the timing constraints, resource constraints, and manufacturing method constraints, a mathematical model for optimizing the aircraft system manufacturing process with the shortest manufacturing cycle and the lowest resource consumption as the optimization goals is established as follows:

[0044] f1=min max(F j )(j=1,2,…,J)

[0045]

[0046] Step 2: Optimize the mathematical model using the improved SPEA2 algorithm:

[0047] The improved SPEA2 algorithm adds a neighborhood search method with solution area density to the original SPEA2 algorithm, organically combining global search with local search.

[0048] S2.1, encoding and decoding of the SPEA2 algorithm;

[0049] S2.2, population initialization: Generate an initial population P0 of size N in a random manner, create an archive set Q0 of size η, and set the algorithm parameters, including the current number of iterations t, the maximum number of iterations T, and the crossover probability Λ c , mutation probability Λ m ;

[0050] S2.3, calculate individual fitness function;

[0051] S2.4, environmental selection;

[0052] S2.5, determine whether to terminate the program.

[0053] Furthermore, the encoding and decoding process of the SPEA2 algorithm is as follows:

[0054] Chromosome coding: Using a three-dimensional coding method based on priority, the three-dimensional coding structure of the chromosome consists of three dimensions: priority, pattern, and time lag coefficient of the aircraft system manufacturing process; Figure 3 As shown in the figure, the first dimension is the priority of process j, and its value range is between [0,1]. The larger the value, the higher the priority. The second dimension is the manufacturing method of process j, and its value range is [0,M j ); the third dimension is the time lag coefficient of process j, which is a real number in the range [0,1] and represents the ratio of the interval between process j and the previous process to the maximum and minimum time lag difference; the priority, manufacturing method and time lag coefficient of process J in aircraft system manufacturing are obtained through chromosome three-dimensional coding;

[0055] chromosome decoding: V a = (v1, v2,... v a ) represents a compatible chain table with a items process, where 1≤a≤J, H a = (h1, h2,... h a ) represents a sequence list of manufacturing methods corresponding to the compatible chain table, Z a = (z1, z2,... z a ) represents a sequence list of time delay coefficients corresponding to the compatible chain table;

[0056] S2.1.1, traverse all unsorted processes, when all immediate previous processes of process a+1 have been located in the compatible chain table V a , then process a+1 is placed into the qualified process set U;

[0057] S2.1.2, traverse the qualified process set U, select the process a+1 with the highest priority in the qualified process set U and place it into the compatible chain table, update the compatible chain table and sequentially empty the qualified process set U;

[0058] S2.1.3, repeat the above S2.1.1 and S2.1.2 until the compatible chain table contains J items processes;

[0059] S2.1.4, after obtaining the complete compatible chain table, H a and Z a are directly obtained according to the decoding of the chromosome;

[0060] S2.1.5, according to the three tables V a , H a and Z a , a serial scheduling method is used to generate an aircraft system manufacturing scheme, and the generated aircraft system manufacturing scheme includes the start time, end time and manufacturing method of each process;

[0061] Further, the individual fitness function calculation method is as follows:

[0062] P t represents the population of iteration t times, Q t represents the archive set of iteration t times, the individual α dominates the individual β, and the strength of the number of individuals dominated by the individual α is

[0063]

[0064] In the formula, |. | represents the size of the set; represents that α dominates β;

[0065] The individual alpha rough fitness is R(alpha), which considers not only the information of the non-dominated set, but also the information of the dominated individual alpha. The solution formula of R(alpha) is as follows:

[0066]

[0067] The smaller the R(alpha) value is, the more the number of the dominated individuals of the individual alpha is. When R(alpha)=0, the individual alpha is a non-dominated individual. When multiple individuals have the same rough fitness, in order to reflect the difference between different individuals, the crowdedness D(alpha) between the individual alpha and other individuals with the same rough fitness is defined, and the calculation method is as follows:

[0068] Firstly, the k-neighbor method is used to evaluate the density between individuals, and the distance z t (α) between the individual alpha and other individuals in the population P t and the archive set Q j is calculated, and is arranged in ascending order as:

[0069] z1(α)≤z2(α)≤…≤z n (α) (18)

[0070] Wherein, n=N+η;

[0071] Then, the kth in the above n distances is selected as [] represents the rounding operation;

[0072] The density D(alpha) of the individual alpha is:

[0073]

[0074] In the formula, 0<D(alpha)<1, k≤n;

[0075] Finally, the exact fitness Ψ(alpha) of the individual alpha is defined as:

[0076] Ψ(i)=R(alpha)+D(alpha) (20)

[0077] In summary, by defining the strength function S(alpha) and the density function D(alpha) of all individuals in the population P t and the archive set Q t , the fitness ψ(alpha) of each individual is finally obtained;

[0078] Further, the environment selection method is as follows:

[0079] In the iteration process of the above step S2.3, the appropriate individual is selected into the archive set Q tThe non-dominated individuals are preferentially selected into the archive set, and when the number of non-dominated individuals is less than the size of the archive set, the dominated individuals with the minimum fitness are selected into the archive set in sequence; the environmental selection method ensures the diversity of the population and avoids the result from falling into local optimum;

[0080] In the environmental selection process, the individuals with the fitness value smaller than 1 are preferentially stored in the archive set, i.e.

[0081] Q t+1 ={α|α∈P t +Q t ∩Ψ(α)<1} (21)

[0082] When the number of individuals in Q t+1 is less than the size of the archive set, i.e. Q t+1 <η, then (η-|Q t |) individuals with the minimum fitness are selected from the previous generation P t and Q t+1 into Q t+1 ; when Q t+1 >η, the number of individuals in Q t+1 is reduced by using the pruning process, and the pruning process is as follows:

[0083] Firstly, the distance between any two individuals in the archive set is calculated to generate a distance table;

[0084] Then, the two points with the minimum distance ε and γ are selected from the distance table;

[0085] Then, the second nearest neighbors of ε and γ are selected, i.e. the points with the second minimum distance to ε and γ, and the point with the smaller distance to the second nearest neighbors is deleted from the two points ε and γ;

[0086] Further, the judgment method of whether to terminate the program is as follows:

[0087] If t≥T, all non-dominated individuals in Q t+1 are taken as the return result;

[0088] Otherwise, if t

[0089] S2.5.1, the archive set Q t is updated as follows:

[0090] Firstly, the function value of the objective function is normalized:

[0091] In order to avoid the numerical difference between the objective function f1 and the objective function f2 from affecting the optimization result, the following formula is used for processing:

[0092]

[0093] Where φ represents one of the objective functions f1 and f2, f t,φ is the function value to be normalized in the objective function φ of the tth generation in the archive set, and f t,φ The maximum and minimum values ​​in f' t,φ is the normalized objective function φ;

[0094] Secondly, the normalized objective function is divided into solution areas and assigned to the target areas:

[0095] Solution region division: The numerical range of the objective function φ is [0,1], the coordinates of the extreme point in the solution space are (1,1), and a curve is constructed with the extreme point as the center and the objective function range as the radius. Rays are drawn from the extreme point to the equal division points of the curve to divide the solution space into r solution regions with equal angles; for example Figure 4 As shown, the solution area B v (1≤v≤r) corresponding angle θ v (1≤v≤r) are equal;

[0096] Solution region belonging determination: For individual α in the archive set, the solution region to which individual α belongs is determined by calculating the angle between individual α and the boundary of the solution space. The specific method is:

[0097]

[0098] Among them, B v Indicates the solution region to which it belongs, f′ t,1 represents the normalized t-generation objective function f1, f′ t,2 represents the normalized t-generation objective function f2, Indicates rounding down;

[0099] Then, the solution region is classified: the sparse solution region and the crowded solution region in the archive set are determined by the number of solutions in the solution region; the upper bound of the number of solutions in the sparse solution region is defined as BN min , the number of solutions in the crowded solution region is lower bounded by BN max , the number of individuals is greater than 0 and less than BN min The solution region is defined as the sparse solution region, and the solution in the sparse solution region is defined as a sparse solution. Similarly, the number of individuals is greater than BN max The solution region with a value less than η is defined as the crowded solution region, and the solutions in the crowded solution region are defined as crowded solutions;

[0100] Then, the neighborhood solution is generated:

[0101] In the sparse solution region: For each sparse solution μ represents the number of decision variables, is the ξth decision variable of the sparse solution , 1≤ξ≤μ; the neighborhood interval is divided into ρ parts, and the neighborhood search strategy is used to generate neighborhood solutions {y1, y2,..., y σ} (1≤σ≤ρ) around each sparse solution, and the generation method is as follows:

[0102]

[0103] , where τ is the neighborhood radius, y σ,ξ is the ξth decision variable of the neighborhood solution y σ ;

[0104] In the crowded solution region: for each crowded solution , randomly select three solutions Υ=(Υ1,Υ2,...,Υ μ ), e=(e1,e2,...,e μ ) and Φ=(Φ1,Φ2,...,Φ μ ) in the crowded solution region, and perturb the neighborhood of the crowded solution x, and obtain two neighborhood solutions

[0105]

[0106] , where is the Θth decision variable of the crowded solution , is the Θth decision variable of the crowded solution , is the Θth decision variable of the crowded solution , Υ Θ , e Θ , Φ Θ are the Θth decision variables of the crowded solutions Υ, e, Φ respectively, and Δ is the perturbation factor;

[0107] Finally, the solution region processing:

[0108] For each sparse solution region, combine the sparse solutions and neighborhood solutions in the solution region to form a new population P', and select the solution with the lowest rough fitness value in the population in turn, and if the rough fitness is the same, select the solution with the smallest density D(α) according to the k-neighbor method in turn, until the number of solutions reaches the number boundary BN min of the sparse solution region; for each crowded region, combine the crowded solutions and neighborhood solutions in the solution region to form a new population P'', and select the solution with the lowest rough fitness R(α) in the population in turn, and if the rough fitness is the same, select the solution with the smallest density D(α) according to the k-neighbor method in turn, until the number of solutions reaches the number boundary BNmax ;

[0109] To ensure that the size of the archive set remains unchanged, if the size of the archive set is less than η after the solution area processing, individuals with high fitness are selected from the solutions deleted in the solution area processing results and added to the archive set until the archive set size reaches η. If the size of the archive set is greater than η after the solution area processing, the archive set is reduced through the same pruning process as in S2.4 until the archive set size reaches η.

[0110] At this point, the archive set is updated;

[0111] S2.5.2, Individual Selection:

[0112] Use the binary tournament method to extract the updated archive set Q t Select individuals from the mating pool G t , from the archive set Q t Two individuals are randomly selected, and the one with higher fitness among the two individuals enters the mating pool G t , the individual with lower fitness among the two individuals is put back into the archive set Q t , re-execute the selection operation of the binary tournament method until the mating pool G t The size of the population reaches N;

[0113] S2.5.3, for mating pool G t The individuals in the crossover and mutation:

[0114] Crossover operation: According to the crossover probability Λ c , for the mating pool G t The individuals in are cross-operated with each other; Figure 5 As shown, for any two individuals g1 and g2, two integers q1, q2∈[1,J] are randomly generated; the genes q1 and q2 in the chromosome of individual g1 are kept in their original positions, and the genes with unchanged positions are copied to the offspring individual c1. Except for the genes with unchanged positions in the chromosome of individual g1, the other genes are replaced by genes with corresponding positions in g2; the genes q1 and q2 in the chromosome of parent individual g2 are kept in their original positions and copied to the offspring individual c2. Except for the genes with unchanged positions in the chromosome of individual g2, the other remaining genes are replaced by genes with corresponding positions in the parent g1.

[0115] Mutation operation: according to the mutation probability Λ m , for the mating pool G t The individual g3 in the mutated state is mutated to obtain the offspring individual c3; Figure 6 As shown, a mutation probability Λ is set to determine whether each gene on the chromosome mutates g , for each gene in the parent individual g3, randomly generate a random number ο∈[0,1], if ο<Λg For the priority of the gene corresponding to the random number, a random real number in the interval [0, 1] is taken instead of the original gene, for the manufacturing mode of the gene corresponding to the random number, a random integer in the interval [0, M J ) is taken instead of the original gene, for the time delay coefficient of the gene corresponding to the random number, a random real number in the interval [0, 1] is taken instead of the original gene;

[0116] S2.5.4, after completing the crossover and mutation, the next generation population P t+1 , let the iteration number be t+1, go to S2.3, continue the loop until t≥T, then take all non-dominated individuals in Q t+1 as the return result.

[0117] Advantages of the present application:

[0118] 1. The manufacturing process optimization method of the present application involves multiple restrictions of time sequence constraints, resource constraints and manufacturing mode selection, realizes the dual objectives of cycle minimization and resource consumption optimization of the aircraft system manufacturing process;

[0119] 2. The present application comprehensively considers the actual demand and complex situation in the aircraft manufacturing process, so that the improved SPEA2 algorithm can more comprehensively solve the practical problems in the optimization of the aircraft system manufacturing process;

[0120] 3. The improved SPEA algorithm adopts a neighborhood search strategy based on solution region density, realizes the organic combination of global search and local search, effectively ensures the diversity of the population, can make the final solution set excellent in diversity and uniformity, and due to the enhanced search effect, can accelerate the convergence speed of the SPEA algorithm. BRIEF DESCRIPTION OF DRAWINGS

[0121] Figure 1 is the flowchart of the present application;

[0122] Figure 2 is the directed connection diagram of the four generalized priority relations and the corresponding flow relation;

[0123] Figure 2 , (a) is the S-S relation, (b) is the F-S relation, (c) is the S-F relation, and (d) is the F-F relation;

[0124] Figure 3 is a schematic diagram of chromosome coding;

[0125] Figure 4 is a schematic diagram of the division of the solution space;

[0126] Figure 5 is a schematic diagram of chromosome crossover;

[0127] Figure 6 is a schematic diagram of chromosome variation;

[0128] Figure 7 is a time sequence constraint diagram of component manufacturing process;

[0129] Figure 8 is an improved SPEA2 optimization result diagram;

[0130] Figure 9 is a HUDS Gantt chart. DETAILED DESCRIPTION

[0131] The technical solutions of the present application will be specifically described below with reference to the drawings.

[0132] The present application provides an aircraft system manufacturing process optimization method based on improved SPEA2, which plans the start and end time of each aircraft process to shorten the construction period and reduce resource consumption. A neighborhood search method based on solution region density is used to re-search regions with too small or too large solution region density. The manufacturing process of the airborne head-up display system (HUDS) is taken as an example for optimization. Figure 7 The manufacturing implementation process of each component realizes the constraint relationship in time sequence. As can be seen from the figure, HUDS is divided into three subsystems, namely data processing subsystem (DPS), control subsystem (CS), and display subsystem (DS). The numbered boxes in the figure represent a process, and the numbers in the box represent the process number. The arrows between the process block diagrams represent the generalized priority relationship between processes. The number pair on the arrow represents the minimum and maximum time lag. For example, for two processes numbered 1 and 2 in the figure, the arrow is drawn from the left side of box 1 and drawn into the left side of box 2, indicating that the two processes have an S-S relationship. The number pair (3, 5) on the arrow indicates that the minimum time lag is 3 days and the maximum time lag is 5 days. For example, for two processes numbered 5 and 9 in the figure, the arrow is drawn from the right side of box 5 and drawn into the left side of box 9, indicating that the two processes have an F-S relationship. The number pair (2, 4) on the arrow indicates that the minimum time lag is 2 days and the maximum time lag is 4 days. From the figure, it can be seen that the manufacturing process of each component is divided into a plurality of processes, and the processes are connected by arrows to form a process network. The arrows between the processes represent the generalized priority relationship between the processes. The number pair on the arrow represents the minimum and maximum time lag between the two processes. The time lag is the time interval between the start time of the process on the left side of the arrow and the end time of the process on the right side of the arrow. The minimum time lag is the minimum time interval between the start time of the process on the left side of the arrow and the end time of the process on the right side of the arrow. The maximum time lag is the maximum time interval between the start time of the process on the left side of the arrow and the end time of the process on the right side of the arrow. The time lag, the minimum time lag, and the maximum time lag are all positive numbers. Figure 7It can also be seen that the number of timing constraints within each subsystem is relatively large, and the number of timing constraints between subsystems is relatively small. In addition, Table 1 shows the consumption of different types of resources by each component of the airborne head-up display system (HUDS) under different manufacturing methods. As can be seen from Table 1, each process has two manufacturing methods, and in general, the first manufacturing method requires more time but less resources, and the second manufacturing method requires less time but more resources. In the case of time and resource constraints, it is necessary to reasonably select the manufacturing method for each process to achieve the completion of the manufacturing of the HUDS in a shorter time with less resource consumption. In addition, as can be seen from Table 1, the number of renewable resources is 3, including technical personnel, display equipment, and manufacturing equipment; the number of non-renewable resources is 3, including peripheral elements, core elements, and hardware materials. Therefore, the manufacturing processes of the 30 components are planned to achieve the shortest manufacturing cycle and the lowest resource consumption.

[0133] Table 1 Different manufacturing methods of HUDS component-level manufacturing processes

[0134]

[0135]

[0136]

[0137] The following are the specific implementation steps of the present application:

[0138] Step 1, build a mathematical model of the aircraft system manufacturing process:

[0139] S1.1, build a mathematical model of timing constraints:

[0140] Define that the aircraft manufacturing process includes J manufacturing processes, J = 30, wherein process i and process j are any two manufacturing processes in the J manufacturing processes; the start time of process i is S i , the end time is F i , the duration of process i is d i ; the start time of process j is S j , the end time is F j , the duration of process j is d j , the weight between process i and process j is l ij , and the weight between process j and process i is l jiThe generalized precedence relations of the aircraft system manufacturing process timing constraints are divided into four types: Start-Start (SS), Finish-Start (FS), Start-Finish (SF), and Finish-Finish (FF). These four types of process relations are also accompanied by minimum and maximum time delays.

[0141] The end time of the process is expressed as the start time plus the manufacturing cycle, and the generalized priority relationship between process i and process j is expressed and uniformly converted as follows:

[0142] S i +l ij ≤S j ,

[0143]

[0144] S1.2, Construct a mathematical model of manufacturing method constraints:

[0145] Among the J processes in the aircraft manufacturing process, any process j has M j Manufacturing method, M j =2; The selection of manufacturing methods in the aircraft manufacturing process is as follows:

[0146]

[0147] Among them, x jm represents the manufacturing method selected by process j, x jm When the value is 1, it means that process j selects manufacturing method m; jm When the value is 0, it means that process j does not select manufacturing method m;

[0148] Because M j Each manufacturing process can only choose one manufacturing method, and cannot adopt multiple manufacturing methods at the same time, so x jm The constraints that need to be met are as shown in the formula:

[0149]

[0150] S1.3, Constructing a mathematical model of resource constraints:

[0151] For updateable resources: for The set of processes being executed at the moment, the aircraft system manufacturing process needs to consume Ξ kinds of renewable resources to complete, Ξ = 3, when process j is executed under manufacturing mode m, the number of renewable resources per day is The demand for this resource is o jmk , No. The daily supply of renewable resources is Since the use of renewable resources at any time cannot exceed the maximum supply, the mathematical model of the aircraft system manufacturing process must meet the following constraints:

[0152]

[0153] For non-renewable resources: The aircraft system manufacturing process needs to consume W non-renewable resources, W = 3, and the jth process is executed under the manufacturing mode m. w The demand for non-renewable resources is o jmw , No. w The supply of non-renewable resources is O w Since the use of non-renewable resources cannot exceed the maximum supply during the entire process of aircraft system manufacturing execution, the following constraints must be met:

[0154]

[0155] When the objective function f1 conforms to the following formula, the optimized aircraft manufacturing cycle is the shortest;

[0156] f1=minmax(F j )(j=1,2,…,J)

[0157] When the aircraft manufacturing cycle is the shortest, and the objective function f2 meets the following formula, the resource consumption in the aircraft manufacturing process is the lowest;

[0158]

[0159] Among them, the unit price of the kth renewable resource is C k , the unit price of the lth non-renewable resource is C l ;

[0160] Therefore, the mathematical model for aircraft system manufacturing process optimization is as follows:

[0161] f1=minmax(F j )(j=1,2,…,J)

[0162]

[0163] Step 2: Optimize the mathematical model using the improved SPEA2 algorithm:

[0164] The improved SPEA2 algorithm adds a neighborhood search method with solution area density to the original SPEA2 algorithm, organically combining global search with local search.

[0165] S2.1, encoding and decoding of SPEA2 algorithm:

[0166] The three-dimensional coding mode of the chromosome code adopts priority, and the three-dimensional coding structure of the chromosome is priority (Priority), manufacturing mode (Pattern) and time lag coefficient (Timelag) of the manufacturing process of the aircraft system; the first dimension is the priority of the process j, and the value range is between [0, 1], and the larger the value is, the higher the priority is; the second dimension is the manufacturing mode of the process j, and the value is an integer in [0, M j ); the third dimension is the time lag coefficient of the process j, and the value is a real number in [0, 1], which represents the interval time of the process j and the previous process of the process j accounts for the proportion of the maximum and minimum time lag difference; the priority, manufacturing mode and time lag coefficient of the J processes of the aircraft system manufacturing are obtained through the three-dimensional coding of the chromosome;

[0167] Chromosome decoding: V a =(v1,v2,...v a ) represents a compatible chain table with a processes, wherein 1≤a≤J, H a =(h1,h2,...h a ) represents a sequence list of manufacturing modes corresponding to the compatible chain table, Z a =(z1,z2,...z a ) represents a sequence list of time lag coefficients corresponding to the compatible chain table;

[0168] S2.1.1, traverse all unsorted processes, when all the immediately preceding processes of the process a+1 have been located in the compatible chain table V a , then the process a+1 is placed in the qualified process set U;

[0169] S2.1.2, traverse the qualified process set U, select the process a+1 with the highest priority in the qualified process set U, place it in the compatible chain table, and update the compatible chain table, and sequentially empty the qualified process set U;

[0170] S2.1.3, repeat the above S2.1.1 and S2.1.2 until the compatible chain table contains J processes;

[0171] S2.1.4, after obtaining the complete compatible chain table, H a and Z a are directly obtained according to the decoding of the chromosome;

[0172] S2.1.5, according to the three tables V a , H a and Z a , a aircraft system manufacturing scheme is generated by a serial scheduling method, and the generated aircraft system manufacturing scheme includes the start time, end time and manufacturing mode of each process;

[0173] S2.2, population initialization:

[0174] An initial population P0of size N=100 is generated randomly, and an archive set Q0of size η=50 is created. Set the current iteration number t=0, the maximum iteration number T=200, the crossover probability Λ c =0.9, and the mutation probability Λ m =0.1.

[0175] S2.3, calculation of individual fitness function:

[0176] P t represents the population at iteration t, Q t represents the archive set at iteration t, the individual dominated by individual α is β, and the strength of the number of individuals dominated by individual α is

[0177]

[0178] In the formula, |. | represents the size of the set; represents that α dominates β;

[0179] The rough fitness of individual α is R(α), and the solution formula of R(α) is as follows:

[0180]

[0181] When multiple individuals have the same rough fitness, in order to reflect the difference between different individuals, the crowding degree D(α) between individual α and other individuals with the same rough fitness is defined, and the calculation method is as follows:

[0182] First, the k-neighborhood method is used to evaluate the density between individuals, and the distance z j (α) between individual α and other individuals in the population P t and the archive set Q t is calculated and arranged in ascending order as:

[0183] z1(α)≤z2(α)≤…≤z n (α)

[0184] Wherein, n=N+η;

[0185] Then, the kth in the above n distances is selected as [] represents the rounding operation;

[0186] The density D(α) of individual α is:

[0187]

[0188] In the formula, 0<D(α)<1, k≤n.

[0189] Finally, the exact fitness Ψ(α) of individual α is defined as:

[0190] Ψ(i)=R(α)+D(α)

[0191] In summary, by t and archive set Q t All individuals in the network define the intensity function S(α) and the density function D(α), and finally obtain the fitness of each individual ψ(α);

[0192] S2.4, Environment Selection:

[0193] In the iterative process of the previous step S2.3, suitable individuals are selected through environmental selection to enter the archive set Q t , non-dominated individuals are preferentially selected to enter the archive set. When the number of non-dominated individuals is smaller than the size of the archive set, the dominated individuals with the smallest fitness are selected in turn to enter the archive set. This method of environmental selection ensures the diversity of the population and avoids the result from falling into the local optimum.

[0194] In the process of environmental selection, individuals with an exact fitness value Ψ(α) smaller than 1 are first stored in the archive set, i.e.

[0195] Q t+1 ={α|α∈P t +Q t ∩Ψ(α)<1}

[0196] When Q t+1 The number of individuals in the archive set is less than the size η of the archive set, that is, Q t+1 <η, then in the previous generation P t and Q t Select (η-|Q t+1 |) excellent individuals with small fitness values ​​enter Q t+1 In; when Q t+1 When >η, the pruning process is used to reduce Q t+1 The number of individuals in the pruning process is as follows:

[0197] First, the distance between any two objects in the archive is calculated to generate a distance table. Then, the two points ε and γ with the smallest distance are selected from the distance table. Then, the second nearest neighbor of ε and γ is selected, that is, the point with the second smallest distance to ε and γ. The distances of the second nearest neighbor are compared, and the point with the smallest distance to the second nearest neighbor is deleted.

[0198] S2.5, determine whether to terminate the program:

[0199] If t≥T, then Q t+1all non-dominated individuals as the return result; otherwise, if t < T, continue with the following steps:

[0200] S2.5.1, updating the archive set Q t

[0201] First, the function value of the objective function is normalized by the following formula:

[0202]

[0203] where φ represents one of the objective functions f1 and f2, f t,φ is the function value to be normalized in the objective function φ of the tth generation in the archive set, and are the maximum and minimum values in f t,φ , respectively, and f' t,φ is the normalized objective function φ;

[0204] Second, the normalized objective function is subjected to solution region division and attribution determination:

[0205] Solution region division: the value range of the objective function φ is [0, 1], and the extreme point coordinates of the solution space are (1, 1). A curve is constructed with the extreme point as the center and the objective function range value as the radius. From the extreme point, a ray is drawn to the equally divided point of the curve, and the solution space is equally divided into r = 10 solution regions;

[0206] Solution region attribution determination: for an individual α in the archive set, the angle between the individual α and the boundary of the solution space is calculated to determine the solution region to which the individual α belongs. The specific method is as follows:

[0207]

[0208] where B v represents the solution region to which the individual α belongs, f' t,1 represents the normalized tth generation objective function f1, and f' t,2 represents the normalized tth generation objective function f2, represents the floor function;

[0209] Then, solution region classification: the number of solutions in the solution region is used to determine the sparse solution region and the crowded solution region in the archive set. The upper bound BN min of the number of solutions in the sparse solution region is defined as 4, and the lower bound BN max of the number of solutions in the crowded solution region is defined as 8. The solution region with a number of individuals greater than 0 and less than BN min is defined as the sparse solution region, and the solutions in the sparse solution region are defined as sparse solutions. Similarly, the solution region with a number of individuals greater than BN max ​The solution region smaller than η is defined as the crowded solution region, and the solution in the crowded solution region is defined as the crowded solution;

[0210] Then, the neighborhood solution is generated:

[0211] In the sparse solution region: for each sparse solution μ represents the number of decision variables, which is 3, respectively, the priority of the process, the execution mode and the time delay coefficient, is the ξth decision variable of the sparse solution , 1≤ξ≤μ; the neighborhood interval is divided into ρ parts, ρ takes the value of 6, and the neighborhood search strategy is used to generate neighborhood solutions {y1, y2,..., y σ}(1≤σ≤ρ) around each sparse solution, and the generation method is as follows:

[0212]

[0213] Where τ is the neighborhood radius, where the neighborhood radius of the priority dimension takes the value of 0.1, the neighborhood radius of the execution mode dimension takes the value of 1, and the neighborhood radius of the time delay coefficient dimension takes the value of 0.1, y σ,ξ is the ξth decision variable of the neighborhood solution y σ ;

[0214] In the crowded solution region: for each crowded solution In the crowded solution region, three solutions Υ=(Υ1,Υ2,...,Υ μ ) are randomly selected, e=(e1,e2,...,e μ ) and Φ=(Φ1,Φ2,...,Φ μ ) are used to perturb the neighborhood of the crowded solution x, and two neighborhood solutions

[0215]

[0216] Where, is the Θth decision variable of the crowded solution , is the Θth decision variable of the crowded solution , is the Θth decision variable of the crowded solution Υ Θ , e Θ , Φ Θ are the Θth decision variables of the crowded solutions Υ, e, Φ, respectively, and Δ is the perturbation factor, which takes the value of 0.03;

[0217] Finally, the solution region is processed:

[0218] For each sparse solution region, the sparse solutions in the solution region and the neighborhood solutions are combined to form a new population P', and the solution with the lowest rough fitness value is selected from the population in turn, and if the rough fitness is the same, the solution with the minimum density D(a) is selected according to the k-neighbor method, until the number of solutions reaches the number boundary BN of the sparse solution region min ; for each crowded solution region, the crowded solutions in the solution region and the neighborhood solutions are combined to form a new population P", and the solution with the lowest rough fitness R(a) is selected from the population in turn, and if the rough fitness is the same, the solution with the minimum density D(a) is selected according to the k-neighbor method, until the number of solutions reaches the number boundary BN of the crowded solution region max ;

[0219] To ensure that the size of the archive set is unchanged, if the size of the archive set after the solution region processing is less than η, the individuals with high fitness are selected from the deleted solutions in the solution region processing result in turn to enter the archive set, until the size of the archive set reaches η, and if the size of the archive set after the solution region processing is greater than η, the archive set is reduced through the same pruning process as in S2.4 until the size of the archive set reaches η;

[0220] At this point, the update of the archive set is completed;

[0221] S2.5.2, selection:

[0222] The binary tournament method is used to select individuals from the updated archive set Q t into the mating pool G t Two individuals are randomly selected from the archive set Q t , and the individual with high fitness enters the mating pool G t , and the individual with low fitness is returned to the archive set Q t , and the selection operation of the binary tournament method is re-executed until the size of the mating pool G t reaches the population size N;

[0223] S2.5.3, the individuals in the mating pool G t are crossed and mutated:

[0224] Cross operation: according to the cross probability Λ c , the individuals in the mating pool G tThe individuals in the mating pool G are crossed two by two; for any two individuals g1 and g2, two integers q1, q2 ∈ [1, J] are randomly generated; the genes at positions q1 and q2 in the chromosome of individual g1 are kept in place and copied to the child individual c1, and the genes at other positions in the chromosome of individual g1 are replaced by the genes at the corresponding positions in individual g2; the genes at positions q1 and q2 in the chromosome of parent individual g2 are kept in place and copied to the child individual c2, and the remaining genes in the chromosome of individual g2 are replaced by the genes at the corresponding positions in parent individual g1;

[0225] Mutation operation: according to mutation probability Λ m The individuals g3 in the mating pool G t are subjected to mutation operation to obtain child individuals c3; a mutation probability Λ is set for determining whether each gene on the chromosome is mutated g For each gene in parent individual g3, a random number ο ∈ [0, 1] is randomly generated, and if ο < Λ g For the priority of the gene corresponding to the random number, a random real number in the interval [0, 1] is used to replace the original gene, for the manufacturing method of the gene corresponding to the random number, a random integer in the interval [0, M J ) is used to replace the original gene, and for the time delay coefficient of the gene corresponding to the random number, a random real number in the interval [0, 1] is used to replace the original gene.

[0226] S2.5.4, after crossing and mutation, the next generation population P is formed t+1 Let the iteration number be t+1, go to S2.3, and continue the loop until t≥T, then all non-dominated individuals in Q t+1 are taken as the return result.

[0227] The improved SPEA2 proposed in this paper is applied to the optimization and solution of this example, and the obtained set of Pareto solutions is shown in FIG. 1 Figure 8 Each point in the figure represents a Pareto optimal solution, the horizontal coordinate represents the manufacturing time of the component level of the airborne head-up display system, and the vertical coordinate represents the resource consumption of the system. For each Pareto solution, its process order and manufacturing method can be determined through model calculation. From the distribution of the solutions, it can be seen that in the actual optimization process, the improved SPEA2 proposed in this paper can obtain a more superior solution set than the original SPEA2 algorithm, with more solutions, more uniform distribution of solutions, better diversity and convergence, proving that the improved SPEA2 can obtain a better solution than the original SPEA2, and proving that the proposed method can be effectively applied to the optimization of aircraft system manufacturing processes.

[0228] In addition, as shown in FIG. 2 Figure 9The HUDS (Head-up Display System, HUDS) Gantt chart shows the feasible detailed scheduling scheme corresponding to the solution of the manufacturing cycle of 97 days and the resource consumption of 66861 yuan in the Pareto diagram, which is represented in the form of a Gantt chart. The horizontal coordinate in the figure represents the execution days of the scheme, and the vertical coordinate represents the manufacturing of 30 components. Each rectangle represents the manufacturing process of a component being executed, where orange yellow represents the manufacturing process using manufacturing mode 1, and blue represents the manufacturing process using manufacturing mode 2. As can be seen from the figure, this aspect can reasonably arrange the start time and end time of the 30 processes under the conditions of meeting the timing constraints and resource constraints, and can reasonably select the two manufacturing modes to achieve the purpose of shortening the manufacturing cycle and reducing the resource consumption.

Claims

1. An aircraft system manufacturing process optimization method based on improved SPEA2 algorithm, characterized in that, The manufacturing process optimization method comprises the following steps: Step one, constructing a mathematical model of aircraft system manufacturing process; Step two, optimizing the mathematical model by using improved SPEA2 algorithm; The optimization method of the mathematical model is as follows: S2.1, coding and decoding of SPEA2 algorithm; S2.2, population initialization: generate an initial population P0 of size N in a random manner, create an archive set Q0 of size η at the same time, set algorithm parameters, the algorithm parameters include the current iteration number t, the maximum iteration number T, the crossover probability Λ c , the mutation probability Λ m ; S2.3, calculating individual fitness function; S2.4, environment selection; S2.5, judging whether to terminate the program; In S2.1, the coding and decoding process of the SPEA2 algorithm is as follows: Chromosome coding: the three-dimensional coding mode of priority is adopted, and the three-dimensional coding structure of the chromosome is the priority, manufacturing mode and time lag coefficient of the aircraft system manufacturing process. The first dimension is the priority of the process j, and the value range is between [0, 1]. The larger the value is, the higher the priority is. The second dimension is the manufacturing mode of the process j, and the value is an integer in [0, M j The third dimension is the time lag coefficient of the process j, and the value is a real number in [0, 1], which represents the proportion of the interval time of the process j and the previous process j to the maximum and minimum time lag difference. The priority, manufacturing mode and time lag coefficient of the J processes of the aircraft system manufacturing are obtained through the three-dimensional coding of the chromosome. Chromosome decoding: V a = (v1, v2,... v a ) denotes a compatible chain list with a items procedure, where 1≤a≤J, H a = (h1, h2,... h a ) denotes a sequential list of manufacturing ways corresponding to the compatible chain list, Z a = (z1, z2,... z a ) denotes a sequential list of time lag coefficients corresponding to the compatible chain list; S2.1.1, traverse all unsorted procedures, when all immediate predecessor procedures of procedure a+1 have all been located in the compatible list V a then procedure a+1 is placed in the eligible procedure set U; S2.1.2, traverse the qualified process set U, select the process a+1 with the highest priority in the qualified process set U and place it into the compatible chain table, update the compatible chain table and empty the qualified process set U in turn; S2.1.3, repeat the above S2.1.1 and S2.1.2 until the compatible chain table contains J processes; S2.1.4, after obtaining the complete list of compatible chains, H is obtained directly from the decoding of the chromosome a and Z a ; S2.1.5, according to three tables V a , H a and Z a , generating an aircraft system manufacturing plan by a serial scheduling method, the generated aircraft system manufacturing plan including a start time, an end time and a manufacturing manner of each flow; Step three, optimizing the aircraft system manufacturing process according to the mathematical model.

2. The method of claim 1, wherein the method is characterized by, In step one, the construction steps of the mathematical model are as follows: S1.1, constructing a mathematical model of time sequence constraint: Firstly, the time constraint relationship is combed: the aircraft manufacturing process includes J processes, wherein process i and process j are any two manufacturing processes in the J processes; the start time of process i is S i , the end time is F i , the duration of process i is d i ; the start time of process j is S j , the end time is F j , the duration of process j is d j , the weight between process i and process j is l ij , and the weight between process j and process i is l ji ; The generalized priority relationship of the time sequence constraint of the aircraft system manufacturing process includes four types of start-start, end-start, start-end and end-end, and the four types of process relationship are accompanied by minimum time lag and maximum time lag; In S-S relations, the directed connection leads out from the left side of process i and into the left side of process j, the pair of numbers on the directed connection denotes that process j is to start after process i has started to within this time, where are the minimum and maximum time lags, respectively, under S-S relations In the F-S relationship, the directed connection is drawn from the right side of process i and into the left side of process j, with the pair of numbers on the directed connection denotes that process j is to start after process i has finished to within this time, where are the minimum and maximum time lags, respectively, under the F-S relationship; In S-F relations, the directed connection is drawn from the left side of process i and into the right side of process j, the pair of numbers on the directed connection denotes that process j is to be started after process i has finished to within this time, where are the minimum and maximum time lags, respectively, under S-F relations In F-F relations, the directed connection leads out from the right side of process i and into the right side of process j, the pair of numbers on the directed connection denotes that process j is to be started after process i has finished to within this time, where are the minimum and maximum time lags, respectively, under F-F relations; Then, the time sequence constraint model is determined: the end time of the process is represented by the start time and the manufacturing period, and the generalized priority relationship between process i and process j is represented and uniformly transformed into formulas (1) to (8) as shown: For S-S relationship, the start time of process j is before the start time of process i, so the following condition must be met The relationship between process i and process j is converted into the following formula:​​ For F-S relationship, the start time of process j is before the start time of process i After, Before, so it needs to meet The relationship between process i and process j is converted to the following formula: For S-F relationship, the start time of process j is before the start time of process i After, Before, so it needs to satisfy The relationship between process i and process j is converted to the following formula: For F-F relationship, the start time of process j is after the start time of process i later, earlier, so it needs to satisfy The relationship between process i and process j is converted to the following formula: S1.2, constructing a mathematical model of manufacturing mode constraint: In the J-item process of the aircraft manufacturing process, any item process j has M j manufacturing methods; the selection method of the manufacturing methods of the aircraft manufacturing process is shown in formula (9); wherein x jm represents the manufacturing mode selected by flow j, x jm takes the value of 1 when flow j selects manufacturing mode m; x jm takes the value of 0 when flow j does not select manufacturing mode m; Because M j Each manufacturing process can only choose one manufacturing method, and multiple manufacturing methods cannot be adopted at the same time, so x jm The constraint condition needs to be met, as shown in equation (10): S1.3, constructing a mathematical model of resource constraint: For renewable resources: For the set of processes being executed at time t, the airplane system manufacturing process needs to consume Ξ kinds of renewable resources, the daily demand of the i th resource for process j executed in mode m is o jmk , and the daily supply of the i th renewable resource is Since the usage of any renewable resource at any time cannot exceed the maximum supply, the airplane system manufacturing process mathematical model needs to satisfy the constraint condition as shown in equation (11): For non-renewable resources: W kinds of non-renewable resources are consumed in the process of completing the aircraft system manufacturing, and the daily demand of the wth kind of non-renewable resource is o jmw when the jth process is performed in the manufacturing mode m w Since the use of non-renewable resources in the entire process of aircraft system manufacturing cannot exceed the maximum supply, the constraint condition shown in equation (12) must be met: When the target function f1 meets the following formula (13), the optimized aircraft manufacturing period is the shortest; f1= min max(F j )(j = 1, 2,..., J) (13) When the aircraft manufacturing period is the shortest, when the target function f2 meets the following formula (14), the resource consumption in the aircraft manufacturing process is the lowest; wherein the unit price of the kth renewable resource is C k , and the unit price of the lth non-renewable resource is C l ; S1.4, finally, while considering the time sequence constraint, resource constraint and manufacturing mode constraint, the optimized target of the aircraft system manufacturing process optimization mathematical model is to minimize the manufacturing period and the resource consumption, and the mathematical model is as follows: f1 = minmax(F j )(j = 1, 2,..., J) 3. The method of claim 1, wherein the method is characterized by: In S2.3, the individual fitness function calculation method is as follows: P t represents the population iterated t times, Q t represents the archive set iterated t times, individual a dominates individual b, and the strength of the domination of the number of individuals dominated by individual a is where |. | denotes the size of a set; denotes that a dominates β; The rough fitness of individual a is R(a), which not only considers the information of the non-dominated set, but also considers the information of the dominated individual a, and the solution formula of R(a) is as follows: The smaller the value of R(a) represents the more individual a dominates, when R(a)=0, individual a is a non-dominated individual; when multiple individuals have the same rough fitness, in order to reflect the difference between different individuals, the crowdedness D(a) between individual a and other individuals with the same rough fitness is defined, and the calculation method is as follows: First, the density between individuals is assessed using the k- nearest neighbor method, calculating the distance z t between individual a and other individuals in population P t and archive set Q j ( a ), in ascending order as: z1(a) < z2(a) <... < z n (α) (18) Wherein, n=N+η; Then, the kth of the n distances is selected as [] denotes a rounding operation; The density D(a) of individual a is: In the formula, 0<D(a)<1, k≤n; Finally, the exact fitness Ψ(a) of individual a is defined as: Ψ(i)=R(a)+D(a) (20) In summary, by defining the strength function S(α) and the density function D(α) for all individuals in the population P t and the archive set Q t the fitness ψ(α) of each individual is finally obtained.

4. The method of claim 1, wherein the method is characterized by, In S2.4, the environment selection method is as follows: Selecting suitable individuals into the archive set Q by environmental selection t , and the non-dominated individuals are selected into the archive set Q preferentially, and when the number of the non-dominated individuals is less than the size of the archive set, the dominated individuals with the minimum fitness are selected into the archive set Q in turn. In the environment selection process, the individual with a value of exact fitness Ψ(a) less than 1 is preferentially stored in the archive set, that is, Q t+1 = {a | a e P t + Q t ∩ Ψ(a) < 1} (21) When Q t+1 The number of individuals in the archive set is less than the size η of the archive set, that is, Q t+1 <η, then in the previous generation P t and Q t Select (η-|Q t+1 |) excellent individuals with small fitness values ​​enter Q t+1 In; when Q t+1 When >η, the pruning process is used to reduce Q t+1 The number of individuals in the pruning process is as follows: Firstly, the distance between any two individuals in the archive set is calculated to generate a distance table; Then, select the two points ε, γ with the minimum distance from the distance table; Then, select the second nearest neighbors of ε, γ, i.e. the two points with the second minimum distance to ε, γ, compare the distances of the second nearest neighbors, and delete the point with the smaller distance to the second nearest neighbors from the two points ε, γ.

5. The method of claim 1, wherein the method is characterized by: In S2.5, the termination procedure is determined as follows: If t ≥ T, then Q t+1 all non-dominated individuals in the population as the return result; Conversely, if t < T, continue with the following steps: S2.5.

1. Updating the archive set Q t Update: First, normalize the function values of the objective functions: In order to avoid the numerical difference between the objective function f1 and the objective function f2 being too large to affect the optimization result, the following formula is used for processing: Where φ represents one of the objective functions f1 and f2, f t,φ is the function value to be normalized in the objective function φ of the tth generation in the archive set, and f t,φ The maximum and minimum values ​​in f′ t,φ is the normalized objective function φ; Secondly, the normalized objective function is divided into solution regions and attributed: Solution region division: the numerical range of the objective function φ is [0, 1], the extreme point coordinates of the solution space are (1, 1), a curve is constructed with the extreme point as the center and the range value of the objective function as the radius, and the solution space is equally divided into r solution regions from the extreme point to the equal point of the curve; solution region B v The corresponding angle θ v Equal, where 1≤v≤r; Solution region attribution: for an individual α in the archive set, the angle between the individual α and the boundary of the solution space is calculated to determine the solution region to which the individual α belongs, and the specific method is as follows: wherein B v represents the assigned solution region, f′ t,1 represents the normalized t-th generation objective function f1, f′ t,2 represents the normalized t-th generation objective function f2, represents the floor function; Then, the solution region is classified: the sparse solution region and the crowded solution region in the archive set are determined by the number of solutions in the solution region; the upper bound of the number of solutions in the sparse solution region is defined as BN min , the number of solutions in the crowded solution region is lower bounded by BN max , the number of individuals is greater than 0 and less than BN min The solution region is defined as the sparse solution region, and the solution in the sparse solution region is defined as a sparse solution. Similarly, the number of individuals is greater than BN max The solution region with a value less than η is defined as the crowded solution region, and the solutions in the crowded solution region are defined as crowded solutions; Then, the neighborhood solution is generated: In sparse solution region: for each sparse solution μ denotes the number of decision variables, is the ξth decision variable of the sparse solution 1≤ξ≤μ; the neighborhood interval is divided into ρ equal parts, and the neighborhood search strategy is used to generate neighborhood solutions {y1, y2,..., y σ}(1≤σ≤ρ) around each sparse solution, in the following way: where τ is the neighborhood radius, y σ,ξ is the ξth decision variable of the neighborhood solution y σ . In the crowded solution region: for each crowded solution In the crowded solution region, randomly select three solutions Y = (Y1, Y2,..., Y μ ), e = (e1, e2,..., e μ ) and Φ = (Φ1, Φ2,..., Φ μ ) perturb the neighborhood of the crowded solution x, and get two neighborhood solutions where is the Θth decision variable of the congestion solution is the Θth decision variable of the congestion solution is the Θth decision variable of the congestion solution Θ , e Θ , Φ Θ are the Θth decision variables of the congestion solutions Y, e, Φ, respectively, and Δ is a perturbation factor;​​​ Finally, the solution region is processed: For each sparse solution region, the sparse solutions in the solution region and the neighborhood solutions are combined to form a new population P', and the solution with the lowest rough fitness value is selected in the population in turn, and if the rough fitness is the same, the solution with the minimum density D(α) is selected in turn according to the k-neighbor method, until the number of solutions reaches the number boundary BN of the sparse solution region min For each crowded solution region, the crowded solutions in the solution region and the neighborhood solutions are combined to form a new population P", and the solution with the lowest rough fitness R(α) is selected in the population in turn, and if the rough fitness is the same, the solution with the minimum density D(α) is selected in turn according to the k-neighbor method, until the number of solutions reaches the number boundary BN of the crowded solution region max ​ In order to ensure that the size of the archive set does not change, if the size of the archive set after the solution region processing is less than η, select the individuals with high fitness from the deleted solutions in the solution region processing result into the archive set in turn until the size of the archive set reaches η, if the size of the archive set after the solution region processing is greater than η, reduce the archive set through the same pruning process as in S2.4 until the size of the archive set reaches η; At this point, the update of the archive set is completed; S2.5.2, individual selection: selecting individuals from the updated archive set Q t into a mating pool G t by randomly selecting two individuals from the archive set Q t , the individual with the higher fitness of the two individuals into the mating pool G t , the individual with the lower fitness of the two individuals back into the archive set Q t , re-executing the selection operation of the binary tournament method until the size of the mating pool G t reaches the population size N; S2.5.

3. Crossing and mutating individuals in the mating pool G t S2.5.3.

1. Crossing two individuals S2.5.4, after crossover and mutation, the next generation population P is formed t+1 , let the iteration number be t+1, go to S2.3, continue the loop until t≥T, then all non-dominated individuals in Q t+1 are returned as the result.

6. The method of claim 5, wherein the method is characterized by, In S2.5.3, the specific method of the crossover operation is as follows: According to the cross probability Λ c , the individuals in the mating pool G t are crossed two by two; for any two individuals g1 and g2, two integers q1, q2 ∈ [1, J] are randomly generated; the genes at positions q1 and q2 in the chromosome of individual g1 are kept in place and are copied to the child individual c1, and the genes at positions other than the positions kept in place in the chromosome of individual g1 are taken from the corresponding positions of individual g2; the genes at positions q1 and q2 in the chromosome of parent individual g2 are kept in place and are copied to the child individual c2, and the remaining genes at positions other than the positions kept in place in the chromosome of individual g2 are taken from the corresponding positions of parent individual g1.

7. The method of claim 5, wherein the method is characterized by: In S2.5.3, the specific method of the mutation operation is as follows: According to the mutation probability Λ m , the individual g3 in the mating pool G t is subjected to a mutation operation to obtain a child individual c3; a mutation probability Λ is set for determining whether each gene on the chromosome is mutated g , for each gene in the parent individual g3, a random number ο∈[0, 1] is randomly generated, and if ο<Λ g , the priority of the gene corresponding to the random number is replaced by a random real number in the interval [0, 1], the manufacturing method of the gene corresponding to the random number is replaced by a random integer in the interval [0, M J ), and the time delay coefficient of the gene corresponding to the random number is replaced by a random real number in the interval [0, 1].

Images