1. A multi-objective optimization method for the thickness of an automobile composite bumper anti-collision beam is characterized by comprising the following steps:

the method comprises the following steps of firstly, measuring the thicknesses of layers of an upper plate, a lower plate, a front plate, a rear plate and a rib plate of a plurality of anti-collision beams respectively to obtain the maximum value and the minimum value so as to determine a value range, and determining the layer sequence corresponding to the thickness of each layer;

generating an initial sample data set, and determining an objective function f (x) and a constraint function g (x) for optimizing the thickness of the anti-collision beam layer;

the sample in which the initial sample data is concentrated comprises the layer thicknesses of an upper plate, a lower plate, a front plate, a rear plate and a rib plate, the mass of the anti-collision beam corresponding to the layer thicknesses, the collision force of the anti-collision beam, the absorption energy of the anti-collision beam and the intrusion amount of the anti-collision beam;

f(x)＝(m,F)，g(x)＝(E,D)；

wherein m represents the mass of the minimum anti-collision beam, F represents the collision force of the minimum anti-collision beam, E represents the intrusion amount of 60mm or less, and D represents the energy absorption of 400J or more;

thirdly, setting initial parameters of a kernel function of the Gaussian random process according to the target function and the constraint function, and training to obtain a GPR model;

obtaining a Pareto front edge solution set through a GPR model, screening out a new sample in the Pareto front edge solution set, and adding the new sample into the initial sample data set;

and circulating the second step to the fourth step until the iteration times are met to obtain an optimized Pareto front edge solution set, and obtaining the optimized layer thicknesses of the upper plate, the lower plate, the front plate, the rear plate and the rib plate of the anti-collision beam according to the optimized Pareto front edge solution set.

2. The multi-objective optimization method for the thickness of the automotive composite bumper anti-collision beam according to claim 1, characterized in that in the step one, the value ranges of the thicknesses of the layers of the upper plate, the lower plate, the front plate, the rear plate and the rib plate of the anti-collision beam are determined as follows:

in the formula, x₁Thickness of the layer, x, of the rear panel of the impact beam₂Thickness of the layer, x, of the front panel of the impact beam₃Thickness of the layer of the rib representing the impact beam, x₄Thickness of the layer, x, of the lower plate of the impact beam₅The ply thickness of the upper panel of the impact beam is shown.

3. The method for multi-objective optimization of automotive composite bumper impact beam thickness of claim 2, wherein x is₄＝x₅。

4. The method for multi-objective optimization of thickness of an automotive composite bumper impact beam according to claim 3, wherein in step one, the sequence of plies for each ply thickness is:

when the thickness of the layer is 3mm, the layer laying sequence is as follows: [ (0/90/+ -45)₂(±45)₂]_s；

When the thickness of the layer is 2.5mm, the layer laying sequence is as follows: [ (0/90/+ -45)₂(±45)]_s；

When the thickness of the layer is 2mm, the layer laying sequence is as follows: [ (0/90/+ -45)]_2s；

When the thickness of the layer is 1.5mm, the layer laying sequence is as follows: [ 0/90/+ -45]_s；

When the thickness of the paving layer is 1mm, the paving sequence is as follows: [ 0/90/+ -. 45]_s。

5. The method for multi-objective optimization of automotive composite bumper impact beam thickness according to claim 3 or 4, wherein in step four, a new sample is screened out in the Pareto frontier solution set, comprising the steps of:

step 1, calculating an Euler distance-based multi-target expectation improvement criterion EIM of each individual in a Pareto frontier solution set_eWherein:

in the formula, S is Pareto front edge solution set, subscript represents target number, superscript represents number of solutions; EIM_e(x) For an expected improvement criterion based on euler distances, EI represents an expected improvement matrix, each line of EI represents an improvement value for all objective functions at the jth leading edge position, each line of EI represents an improvement value for the jth objective function at all leading edge positions; f. of_i ^j(x) For the ith objective function value of the j th individual in the Pareto leading edge solution set,s-f_i(x) A predicted value and an error function value corresponding to a GPR model of the ith target function; phi () is a probability density function and a cumulative distribution function of the standard normal distribution;

step 2, EIM selection_eThe largest individual was taken as the new sample.

6. The method for multi-objective optimization of automotive composite bumper impact beam thickness according to claim 5, characterized in that in said second step, initial samples are generated by means of an optimal Latin-ultrasonic test.

7. The method for multi-objective optimization of automotive composite bumper impact beam thickness according to claim 6, characterized in that the number of initial samples generated by the optimal Latin hyper-square test is 50.

Background

At present, multi-objective optimization is mostly based on an optimization method of a proxy model, the proxy model with sufficient precision needs to be established before optimization, and under normal conditions, the precision of the proxy model is only up to more than 0.9, and the optimization result is accurate and credible. However, for high-dimensional highly nonlinear mathematical models, such as automobile collision and the like, hundreds or even thousands of sets of design of experiments (DOEs) are required to be used for simulation and emulation. This cost is not affordable. In addition, most of the traditional proxy models do not directly reflect the uncertainty of the model, so that the optimization result generates larger deviation, and when the optimization is performed for higher dimensionality, the optimization cannot be performed smoothly. In addition, for one-step optimization design, modeling is disposable, no approximate model lifting process exists, and if a model has a large approximate error, optimization based on the model is inevitably inaccurate. At present, no optimization method for realizing multi-objective optimization by combining sequence optimization and Gaussian process regression model exists.

Disclosure of Invention

The invention designs and develops a multi-objective optimization method for the thickness of an automobile composite bumper anti-collision beam, which adopts a third-generation non-dominated sequencing sequence optimization method of Gaussian process regression, only needs to obtain a small number of sampling points, updates the sampling points in gradual optimization, adds the sampling points to perform model reconstruction and performs multi-objective optimization to obtain a relatively accurate Pareto front; therefore, the thickness optimization result of the carbon fiber composite anti-collision beam meeting the requirement of collision performance can be quickly and accurately obtained.

The technical scheme provided by the invention is as follows:

a multi-objective optimization method for the thickness of an automobile composite bumper anti-collision beam comprises the following steps:

the method comprises the following steps of firstly, measuring the thicknesses of layers of an upper plate, a lower plate, a front plate, a rear plate and a rib plate of a plurality of anti-collision beams respectively to obtain a maximum value and a minimum value so as to determine a value range (taking the maximum value of a measurement result as an upper limit of the value range and the minimum value as a lower limit of the value range), and determining the layer sequence corresponding to the thickness of each layer;

generating an initial sample data set, and determining an objective function f (x) and a constraint function g (x) for optimizing the thickness of the anti-collision beam layer;

the sample in which the initial sample data is concentrated comprises the layer thicknesses of an upper plate, a lower plate, a front plate, a rear plate and a rib plate, the mass of the anti-collision beam corresponding to the layer thicknesses, the collision force of the anti-collision beam, the absorption energy of the anti-collision beam and the intrusion amount of the anti-collision beam;

f(x)＝(m,F)，g(x)＝(E,D)；

wherein m represents the mass of the minimum anti-collision beam, F represents the collision force of the minimum anti-collision beam, E represents the intrusion amount of 60mm or less, and D represents the energy absorption of 400J or more;

thirdly, setting initial parameters of a kernel function of the Gaussian random process according to the target function and the constraint function, and training to obtain a GPR model;

obtaining a Pareto front edge solution set through a GPR model, screening out a new sample in the Pareto front edge solution set, and adding the new sample into the initial sample data set;

and circulating the second step to the fourth step until the iteration times are met to obtain an optimized Pareto front edge solution set, and obtaining the optimized layer thicknesses of the upper plate, the lower plate, the front plate, the rear plate and the rib plate of the anti-collision beam according to the optimized Pareto front edge solution set.

Preferably, in the first step, the thickness value ranges of the layers of the upper plate, the lower plate, the front plate, the rear plate and the rib plate of the anti-collision beam are determined as follows:

in the formula, x₁Thickness of the layer, x, of the rear panel of the impact beam₂Thickness of the layer, x, of the front panel of the impact beam₃Thickness of the layer of the rib representing the impact beam, x₄Thickness of the layer, x, of the lower plate of the impact beam₅The ply thickness of the upper panel of the impact beam is shown.

Preferably, x is characterized in that₄＝x₅。

Preferably, in the step one, the sequence of the plies corresponding to the thickness of each ply is as follows:

when the thickness of the layer is 3mm, the layer laying sequence is as follows: [ (0/90/+ -45)₂(±45)₂]_s；

When the thickness of the layer is 2.5mm, the layer laying sequence is as follows: [ (0/90/+ -45)₂(±45)]_s；

When the thickness of the layer is 2mm, the layer laying sequence is as follows: [ (0/90/+ -45)]_2s；

When the thickness of the layer is 1.5mm, the layer laying sequence is as follows: [ 0/90/+ -45]_s；

When the thickness of the paving layer is 1mm, the paving sequence is as follows: [ 0/90/+ -. 45]_s。

Preferably, in the fourth step, a new sample is screened out from the Pareto frontier solution set, which includes the following steps:

step 1, calculating an Euler distance-based multi-target expectation improvement criterion EIM of each individual in a Pareto frontier solution set_eWherein:

in the formula, S is a Pareto front edge solution set, subscripts represent the number of targets, and superscripts represent the number of solutions; EIM_e(x) For an expected improvement criterion based on euler distances, EI represents an expected improvement matrix, each line of EI represents an improvement value for all objective functions at the jth leading edge position, each line of EI represents an improvement value for the jth objective function at all leading edge positions; f. of_i ^j(x) For the ith objective function value of the j th individual in the Pareto leading edge solution set,s-f_i(x) A predicted value and an error function value corresponding to a GPR model of the ith target function; phi () is a probability density function and a cumulative distribution function of the standard normal distribution;

step 2, EIM selection_eThe largest individual was taken as the new sample.

Preferably, in the second step, an initial sample is generated by an optimal latin hyper square test.

Preferably, the number of initial samples generated by the optimal latin hyper square test is 50.

The invention has the beneficial effects that:

the multi-objective optimization method for the thickness of the bumper anti-collision beam made of the automobile composite material adopts a third-generation non-dominated sorting sequence optimization method of Gaussian process regression, only a small number of sampling points are needed to be obtained, the sampling points are updated in gradual optimization, the sampling points are added to perform model reconstruction, and multi-objective optimization is performed, so that a relatively accurate Pareto front edge can be obtained; in the sequence optimization process, the error information of the model is fully considered, and the overall improvement of the optimization result is realized in the process of modeling for multiple times; meanwhile, the simulation time and cost are shortened, and the accuracy is better; therefore, the thickness optimization result of the carbon fiber composite anti-collision beam meeting the requirement of collision performance can be quickly and accurately obtained.

Drawings

Fig. 1 is a schematic structural view of an automotive composite bumper impact beam according to the present invention.

FIG. 2 is a schematic diagram of calculating the distance between a reference line and a normalized target in solution space according to the present invention.

Fig. 3 is a schematic diagram of an optimized Pareto front surface according to the present invention.

Detailed Description

The present invention is further described in detail below with reference to the attached drawings so that those skilled in the art can implement the invention by referring to the description text.

The invention provides a multi-objective optimization method for the thickness of an anti-collision beam of an automobile composite bumper, which comprises the following specific processes:

as shown in FIG. 1, a design variable x is set₁Represents the thickness, x, of the impact beam rear panel 110₂Represents the thickness, x, of the impact beam front panel 120₃Represents the thickness, x, of the rib 130 of the impact beam₄Represents the thickness, x, of the lower bumper beam 140₅Representing the thickness of the upper plate 150 of the anti-collision beam, and defining the thickness of the upper plate and the thickness of the lower plate to be the same; the design variable adopts a discrete value taking mode, and the value taking interval is 0.5.

Therefore, the range of the design variables of the thickness of the paving layer is as follows:

the sequence of the layers at each thickness is shown in table 1:

TABLE 1 sequence of plies corresponding to the thickness of each ply

Thickness of the mat	Sequence of layering
		3mm	[(0/90/±45)2(±45)2]s
2.5mm	[(0/90/±45)2(±45)]s
		2mm	[(0/90/±45)]2s
1.5mm	[0/90/±45/±45]s
		1mm	[0/90/±45]s

Defining an optimized mathematical model as follows:

Minimize f_m(x),f_F(x)

Subjectto g_E(x)≥400J

g_D≤60.0mm

1.0≤x₁≤2.5

1.0≤x₂≤3

1.0≤x₃≤2.5

1.0≤x₄≤2。

after the mathematical model is constructed, performing multi-objective optimization design according to a target function and a constraint function, and optimizing a third-generation non-dominated sorting genetic sequence by adopting Gaussian process regression, wherein the method mainly comprises the following steps:

step 1: performing an optimal Latin hyper-Square test (OLHS) on the complex system in an initial design space to generate initial sample points, and obtaining an objective function f (x) and a constraint function g (x) at the initial sample points, as shown in Table 2;

wherein f (x) is ═ f₁(x),f₂(x),...f_j(x)...,f_n(x)]；

g(x)＝[g₁(x),g₂(x),...g_j(x)...,g_t(x)]；

In the formula, n and t are the numbers of the objective function and the constraint function respectively, and x is a design variable.

TABLE 2 initial sample Table

According to the layering structure division, the number of design variables of the anti-collision beam is 5, because the upper laminated plate and the lower laminated plate adopt the same layering mode, the design variables can be reduced to 4, the variable design range is given, and according to the single-layer thickness of the laminated plate, 50 groups of optimal Latin hypercube sampling design is carried out.

Obtaining the objective function f (x) and the constraint function g (x) after finite element solution:

step 2: setting initial parameters of a kernel function of a Gaussian random process according to a target function and a constraint function; the model training method is maximum likelihood estimation, and the posterior density function of the hyperparameter theta is deduced according to Bayes as follows:p (Y | X) is independent of the parameter θ, so that the posterior distribution of the distribution with respect to the parameter for a given observation is equivalent to a constant and does not affect the mean and variance constants of the posterior distribution, when the maximum posterior estimate of the hyperparameter isMaximize P (Y | X, theta) P (theta), when P (theta) is tested as uniformIn distribution, the maximum a posteriori estimate of the hyperparameters is to find θ, the maximum likelihood estimate, which maximizes P (Y | X, θ). Taking the logarithm of the boundary likelihood as a likelihood function ((maximizing) log likelihood margin LML), and solving the partial derivative of the hyper-parameter:

and carrying out regression model hyper-parameter tuning and model training according to the initial sample, wherein one key setting is the selection of a high-dimensional mapping kernel function, and a single kernel method cannot obtain a more accurate regression model for data mixed with various numerical noises. So the initialization prior mean is set to 0; meanfunc [ ]; the kernel function is a composite kernel, and specifically comprises:

square Exponential (SE) function:

marten temperature covariance (Matern) function:

wherein the content of the first and second substances,for the smoothing coefficient, l is the length scale parameter, Γ (x) is the euler second integral, defined in the real domain as:

as a Bessel function of the second kind, the following areStandard solution function y (x) of the partial equation:

rational Quadratic (RQ) function:(l is a length scale parameter, and alpha is a proportional mixing coefficient);

the variance kernel selects rational quadratic and square exponential kernels set to

Initial value of sigma is 2, l₁Initial value of 0.25, initial value of l2 of 2, alpha of 10^-3。

And optimizing the hyper-parameters by adopting a conjugate gradient method;

in the formula, gp is a regression function, hyp and Hyper are kernel function Hyper-parameters, including a mean function meanfunc and a covariance function covfunc. Minimize is the conjugate gradient optimization function, Sample _ x is the original design point, here the design variable [ x₁,x₂,x₃,x₄]. Sample _ F is the objective function of the Sample point, i.e. [ m, F ] in the data set]It is desirable to minimize the mass of the impact beam and minimize the impact force of the impact beam. Sample _ g is a constraint function; as [ D, E ] in the data set](ii) a The restraint invasion amount is not more than 60mm, and the absorption energy is not less than 400J.

And step 3: generating reference point information and initializing individual information; and performing simulated binary intersection and polynomial variation to generate a new individual.

Firstly, generating an initial reference point, wherein the reference point generation mode is as follows: wherein M is the dimension of the objective function, and H is the number of segments;

1) generatingAn initial reference point set with dimension M-1; here the dimension is 2.

X＝(1/H)*nchoosek(0:H+M-2,M-1)；

Wherein nchoosek is a sampling function, the first parameter represents the total number of samples, and the second parameter represents the number of samples;

2) for theA is changed as follows_ij＝x_ij- (j-1)/H, obtaining a transformation set A;

3) with transformation set a, the reference points are generated according to the following equation:

generally, since the number of segments H and the dimension M increase to make the reference point set too large, an inner and outer layer mechanism is adopted: h1 is the number of boundary segments, H2 is the number of inner layer segments;

h1 and H2 are required to satisfy the requirementsN is the total number of the initially generated individuals, and the number of the inner layers is less than or equal to the number of the boundary segments; the boundary layer generates the reference point S1 according to the above 3 processes, the inner layer first generates the reference point S (3 processes), the mapping coefficient is taken to be 0.5, S2 is generated by the following formula, and S1, S2 are combined;

for the 2 target, H1 ═ 149, H2 ═ 0. The population number is set to be initialized to 150, the iteration number is set to be 200,

obtaining the predicted values and the corresponding variances of the target function and the constraint function by using a GPR model, and constructing fitness information;

based on a kernel method, Bayesian linear regression is mapped into nonlinear regression of a high-dimensional space, and a function f (x) is regarded as a Gaussian random process which follows Gaussian distribution and is represented by the following formula:

P(f_*|X,y,X_*,σ²)～N(σ^-2φ(X_*)^TA^-1φy,φ(X_*)^TA^-1φ(X_*))

μ_*＝k(X_*,X)(K+σ²I)^-1y

σ_*＝k(X_*,X_*)-k(X_*,X)(K+σ²I)^-1k(X,X_*)

acquiring the hyperparameters in step 2, estimating target functions, constraint function information and model errors at unknown individual positions by combining distribution information, and constructing a fitness function for non-dominated comparison by utilizing a penalty mechanism;

in the formula, a₁₁,a₁₂,…,a_n1,…a_ntFor the penalty coefficient, PopDec is initialized individual position information, sf (x), sg (x) is an error function of the predicted values of the objective function and the constraint function, and F is a fitness function; carrying out population individual target function, constraint function prediction value and variance estimation according to the steps: PopDec is individual position information generated through genetic manipulation, and has the same dimension as a variable.

[f1,sf1]＝gp(hyp1,'infGaussLik',[],covfunc2,@likGauss,Sample_x,Sample_f1,PopDec)；

[f2,sf2]＝gp(hyp2,'infGaussLik',[],covfunc2,@likGauss,Sample_x,Sample_f2,PopDec)；

[g1,sg1]＝gp(hyp3,'infGaussLik',[],covfunc2,@likGauss,Sample_x,Sample_g1,PopDec)；

[g2,sg2]＝gp(hyp4,'infGaussLik',[],covfunc2,@likGauss,Sample_x,Sample_g2,PopDec)；

f12＝[f1 f2]；

f12_interval＝[sf1 sf2]；

g12＝[-g1+400g2-60]；

g12_interval＝[sg1 sg2]；

z＝[f12 g12]；

f_penalty＝100*max(0,-g1+400)+100*max(0,g2-60)；

y＝f12+f_penalty；

y is a constructed fitness function, and Sample is a design variable in original training data, and Sample _ f1, Sample _ f2, Sample _ g1 and Sample _ g2 correspond to mass, impact force, energy absorption and intrusion response. f. of₁Is a predicted value of the mass model of the anti-collision beam, g₁Predicted value f of energy absorption model of anti-collision beam₂Predicted value of collision force model of the anti-collision beam, g₂And predicting the predicted value of the energy absorption model of the anti-collision beam. Sf, Sg are the target function and the constraint function predicted value variance, respectively.

The generation of individual positions is performed according to genetic manipulation:

ILower＝repmat(Lower,N,1)

IUpper＝repmat(Upper,N,1)

Indiv_Pos＝unifrnd(ILower,IUpper)

wherein Upper and Lower limits are used as variables,

lower＝[1 1 1 1]；upper＝[2.5 3 2.5 2]；

repmat is a repeated constructor, and N is 150; the unifrnd is a random number function, and generates an NxN random number matrix between Lower and Upper to ensure the randomness and diversity of the sample; then, according to the genetic operation process, firstly determining the size of the mating pool, and screening according to the mating pool, wherein the specific implementation method comprises the following steps:

randomly generating individual serial numbers (repeatable, serial numbers between 1 and N) with the row number of 2 and the column number of N, wherein the value is always less than 0 according to the function value of the particle violating the constraint, and the smaller the value is, the higher the satisfaction degree of the constraint is, randomly generating two groups of individuals firstly, and optimizing the two groups of individuals according to the violation value:

Tool＝randi(N,M,N)；

sorting the default values in ascending order to obtain an index corresponding to the original position:

[～,index1]＝Sort(Popcon，[],1)；

in this example, Popcon ═ 0; then, the indexes are sorted in an ascending order to obtain index2, and at this time, index2 includes the size sorting sequence number corresponding to the individual at the original position:

[～,index2]＝sortrows(index1)

according to index2 information, the size sorting of individuals in a mating pool can be known, and the smaller the sequence number is, the higher the possibility of meeting the constraint is, so that two individuals are compared, and the sorting sequence number corresponding to the smaller individual is selected; performing the operation of simulating binary intersection and polynomial variation on the filial generation individuals; simulating binary crossing, simulating single-point crossing in binary coding of genetic algorithm:

the main following is that the mean values of parent and offspring before and after intersection are equal:

define the Spread Factor (Spread Factor. beta.), get c₁And c₂:

Wherein the propagation factor β can be determined by the following formula:

wherein o is a distribution coefficient; also, for polynomial variation, the following formula is followed:

indiv_k'＝indiv_k+δ(Upper_k-Lower_k)

δ₁＝(indiv_k-Lower_k)/(Upper_k-Lower_k)

δ₁＝(indiv_k-Upper_k)/(Upper_k-Lower_k)

in this example, r is a random operator, o, η ═ 20, and if the random operator is 0.5 or less, the propagation factor and the coefficient of variation are calculated according to the above equation. Otherwise, the propagation factor and the coefficient of variation are calculated according to a second expression.

And 4, step 4: for the individuals generated in step 3 (total number 2 × N), according to fitness evaluation, dividing the non-dominant grade to define a reserve set S_t(N(S_t) > N), definition of F_kA non-dominated solution set with a level k; to reserve solution S_tNormalization is carried out, a high-efficiency non-dominant sorting method ENS (efficient non-dominant sort) is adopted to reduce the time and space complexity in calculation when the non-dominant grades are divided, MaxFNo is defined as the current highest dominant grade (the dominant grade is higher when the number is smaller), and the initial value is set to be 0; FrontNo is the leading edge set and is initially set to INF, and the main steps involved are:

1) sorting in a descending order according to the first dimension of fitness of all individuals to obtain sorting information;

2) defining a large-cycle end condition as that the total number of individuals which are not INF in FrontNo exceeds N, starting to add one to MaxFNo each cycle, and recording the maximum non-dominant grade when the cutoff condition is met; for the current individual, performing non-dominant comparison with the individual in the leading edge solution set before, wherein the comparison is performed from the second dimension, and if the individual is dominant, performing non-dominant comparison on the next individual; if the individual is not dominated, FrontNo corresponding to the individual position is the current MaxFNo;

3) and (5) finishing the step 2, acquiring all individuals with non-dominant grades smaller than MaxFNo according to the sorting information and MaxFNo, and defining as P_t+1The individual with a non-dominant rank equal to MaxFNo is defined as F_m,S_tIs P_t+1And F_mA union of (1);

for the acquired S_tThe normalization operation was performed as follows:

4) ideal point position acquisition, for minimum problem, Z_min＝min(f₁,f₂,…f_n)；

5) The coordinate system conversion, with the ideal point as the origin, needs to obtain the extreme points corresponding to each target axis and the intercept at each target axis for the multi-target problem after the function target f' is changed to f-Zmin, and the obtaining mechanism is:

XA＝1

dir is a direction vector with the size of nxn, the main diagonal is 1, and the rest positions are 1 e-6; x is a set of n extreme points of size nxn, A is a set of target intercepts of size nx1, 1 is a column vector of all 1's and size nx 1.

If the rank of matrix X is less than n, then the n limit points cannot form an n-dimensional hyperplane. Even if a hyperplane can be built, it is possible that no intercept or some intercept a may be obtained in some directions_iNot meet the requirements ofIn all of the above cases, for each i e {1,2,. eta., n },is set as S_tAt target f_iThe maximum value of (c). For S_tNormalizing each target function:

and 5: as shown in fig. 2, taking three dimensions as an example, the distance between the reference line (the connecting line between the reference point and the ideal point) and the normalized target in the solution space is calculated, and the distance matrix D is generated according to the following formula:

wherein q represents S_tThe number of the middle individuals, r, corresponds to the number of the reference points; s_tEach individual considers the distance from each reference point, and selects the reference point with the minimum distance for attachment;

step 6: and under the condition of considering the front MaxFNo-1 front edge grade, selecting fewer reference points to be used as the basis for selecting the last front edge population individual according to the selected number of the reference points, and further realizing selection. Selecting an individual according to the niche preserving operator, S obtained in step 5_tAll individuals correspond to reference point information and define rho_jIs P_t+1And the total number of individuals attached to the reference point j is judged as follows:

1) reference point selection, J_min＝{j:argminρ_jIf J_minThere are multiple, random choices

2) Carrying out F_m(final leading edge solution set) selection,there are two cases:

firstly, if F_mNo individuals found in the above-mentioned group were found to be present inCorrespond to thatScreening a new reference point from the positions of the residual reference points according to the step 1 when the current iteration step is not considered; if F_mThere are individuals j and associations, as well as two cases: the first case is that the reference point value for the minimum number of individuals is 0, meaning that for the previous non-dominant front where no individuals are present correspond to it, in order to maintain the diversity of the individuals, from F_mSelecting individuals to join in P (if there are multiple individuals, calculating the distance between each individual and the reference line corresponding to the reference point according to step 5, and selecting the individual with the minimum distance)_t+1，ρ_jThe value of (a) increases with the addition of the individual; in the second case, the reference point value corresponding to the minimum number of individuals is not 0, indicating that there already exists an individual corresponding to the reference point, from F_mA plurality of individuals (if a plurality of individuals exist) attached to the reference point are randomly selected and added; rho_jThe value of (a) increases with the addition of the individual;

the above process is repeated until P_t+1The number of individuals in (1) reaches the required number of individuals N; according to the set iteration times, the individual continuously evolves and evolves to obtain an approximate Pareto approximate leading edge solution;

and 7: taking the approximate Pareto leading edge solution set as a real Pareto leading edge solution set of a real function, and obtaining a multi-target expectation improvement criterion of each individual based on the Euler distance:

the current optimal solution is an approximate Pareto front composed of multiple multidimensional points together:

in the formula, S is the condition of approximate Pareto solution, subscript represents the number of targets, and superscript represents the number of approximate solutions; EIM_e(x) For an expected improvement criterion based on euler distances, EI represents an expected improvement matrix, each line of EI represents an improvement value for all objective functions at the jth leading edge position, each line of EI represents an improvement value for the jth objective function at all leading edge positions; f. of_i ^j(x) To approximate the i-th objective function value of the j-th individual in the Pareto frontier solution set,s-f_i(x) A predicted value and an error function value corresponding to the Gaussian process regression model of the ith target function; phi () is a probability density function and a cumulative distribution function of the standard normal distribution;

selection of EIM_eAdding the individual with the maximum value as a new sample point into the original sample point to generate a new sample point set; and (3) taking the new sample point set as an initial sample point, repeating the steps 1-7 until the iteration times are met, and outputting a final sample point set.

Setting an optimization result according to the steps 1-7 as shown in fig. 3, wherein fig. 3 is a graph for optimizing a Pareto frontage, and the optimization is carried out according to a Toposis criterion (C) based on an entropy method_iCloser to 1, better individual) selects one target solution in frontier Pareto:

entropy method (determining weights): if the front pareto solution contains m individuals, each individual is an n-dimensional objective function. The following formula:

the above formula is normalized according to the column:

calculating respective entropy weights of the n objective functions:

toposis method: judging each element of PY, if the PY belongs to an extremely large index (the larger the index is better), not processing the PY, and if the PY belongs to an extremely small index, normalizing the index:

obtaining an equidirectional processing matrix:

the optimal scheme is the maximum value of each column after the homodromous treatmentThe worst case is the minimum value of each column

Calculating the closeness of each individual to the optimal scheme and the worst scheme:

[x₁,x₂,x₃,x₄,x₅]＝[1,1.5262,1.4610,1.1295,1.1295]

[x₁,x₂,x₃,x₄,x₅]＝[1.0,1.5,1.5,1.0,1.0]

[E,F,D,M]＝[435.5802,10.6882,59.9988,0.8668]

the optimization result meets the target and the constraint condition.

While embodiments of the invention have been described above, it is not limited to the applications set forth in the description and the embodiments, which are fully applicable in various fields of endeavor to which the invention pertains, and further modifications may readily be made by those skilled in the art, it being understood that the invention is not limited to the details shown and described herein without departing from the general concept defined by the appended claims and their equivalents.