Motor multistage multi-target joint optimization design method based on MotorCAD and ANSYS

Abstract

The invention discloses a motor multi-stage multi-objective joint optimization design method based on MotorCAD and ANSYS, and relates to the technical field of motor multi-objective optimization. According to the method, on the basis that MotorCAD quickly calculates the motor performance, multi-target joint optimization design is carried out after an optimal prognosis meta-model is generated based on ANSYS optiSLang sensitivity analysis layering, and the advantages of the MotorCAD road algorithm for quickly analyzing the motor performance and the advantages of ANSYS optiSLang integrated simulation process building platform and multidisciplinary design optimization are fused; according to the method, a Python customized script is utilized to perform parallel computing and automatic calling and analysis of a data result on MotorCAD, the computing cost is reduced, a self-programming algorithm is adopted to perform motor optimization design, and the scheme is compared with a scheme obtained in a software optimization algorithm library, so that the optimization accuracy is improved.

Description

Technical Field

[0001] The invention belongs to the technical field of motor multi-objective optimization, and in particular relates to a motor multi-level multi-objective joint optimization design method based on MotorCAD and ANSYS. Background Art

[0002] Global optimization design methods require repeated calls to the calculation of the objective function. When such methods are used for the optimization design of motors, finite element analysis is usually required to obtain the motor's operating performance. A large number of finite element calculations increase the cost of optimization design, especially for optimization problems based on multi-physics field coupling problems or considering robustness and reliability indicators, which are difficult to achieve.

[0003] The motor design software MotorCAD integrates the magnetic circuit method, thermal circuit method, thermal network method, finite element analysis method, and intelligent optimization algorithm, which can realize the electromagnetic, thermal, and magnetothermal coupling design of the motor. The performance analysis is fast, but the effective optimization accuracy is not high. ANSYS optiSLang is a simulation process integration and design optimization tool, including multidisciplinary optimization, parameter sensitivity analysis, reliability evaluation, etc. The optimization process is fast and convenient, breaking the software barriers, and can realize the integration and optimization of the simulation process.

[0004] Regarding the optimal design of motors, the shortcomings of current performance analysis and optimization design methods are as follows:

[0005] (1) The accurate finite element analysis process is complicated and the calculation cost is too high;

[0006] (2) Each software runs independently, which has many human factors and is prone to errors. Summary of the invention

[0007] In order to overcome the problems existing in the existing motor optimization design, the present invention provides a motor multi-level multi-objective joint optimization design method based on MotorCAD and ANSYS, which aims to achieve multi-level multi-objective rapid comprehensive optimization design of the motor.

[0008] The steps of the multi-level and multi-objective joint optimization design method of the motor based on MotorCAD and ANSYS are as follows:

[0009] Step 1: Analyze the topology and working scenarios of the motor, and determine several design variables, several objective functions and constraints;

[0010] Step 2: According to the motor stator and rotor structural parameters, winding settings, and material selection, use MotorCAD to perform parametric modeling on the motor to obtain a motor parametric model;

[0011] Step 3: Set the solution conditions of the motor parametric model in MotorCAD, and solve the motor's electric field, magnetic field, and temperature field under multidisciplinary coupling to obtain the electromagnetic characteristics, thermal characteristics, and mechanical characteristics of the motor;

[0012] Step 4: Customize the Python script to create the interface between ANSYSoptiSLang and MotorCAD based on the determined design variables, objective function and constraints;

[0013] Step 5: Based on the interface between ANSYS optiSLang and MotorCAD, ANSYS optiSLang uses Python scripts to call the electromagnetic characteristics, thermal characteristics, and mechanical properties of the motor obtained by MotorCAD, and extracts several determined design variables, several objective functions, and constraints;

[0014] Step 6: According to the sensitivity of each design variable relative to several objective functions, the sensitivity coefficient of each design variable is calculated respectively, and a strong sensitivity layer and a weak sensitivity layer are established respectively according to the strength of the sensitivity coefficient;

[0015] S com (x a )=w1|S1(x a )|+w2|S2(x a )+…+w n |S m (x a )| (1)

[0016] w1+w2+…+w m =1 (2)

[0017] Among them, x a is the ath design variable, S com (x a ) is the sensitivity coefficient of the a-th design variable, w m is the weight coefficient of the mth objective function, m is the number of objective functions, S m (x a ) Design variable x a Sensitivity factor to the mth objective function;

[0018] Step 7: Sampling the strong sensitivity layer and the weak sensitivity layer to obtain sample data respectively, and establishing the best prognostic meta-models of the strong sensitivity layer and the weak sensitivity layer respectively based on the sample data of the strong sensitivity layer and the weak sensitivity layer;

[0019] Step 8: Use the prognostic coefficient to evaluate the prediction quality of the best prognostic meta-model in the strong sensitivity layer and the weak sensitivity layer respectively;

[0020] Specifically, for the strong sensitive layer, if the prognosis coefficient is greater than or equal to the set threshold, it is considered that the accuracy meets the requirement, and step 9 is executed. If it is less than the set threshold, it is considered that the accuracy does not meet the requirement, then increase the sample points and return to step 7. For the weak sensitive layer, if the prognosis coefficient is greater than or equal to the set threshold, it is considered that the accuracy meets the requirement, and go to step 9. If it is less than the set threshold, increase the sample points and return to step 7.

[0021] Step 9: Output the best prognostic meta-model that meets the prediction quality requirements;

[0022] Step 10: In the ANSYS optiSLang algorithm library, select the evolutionary algorithm or particle swarm algorithm to optimize several objective functions of the motor based on the best prognosis metamodel to obtain the Pareto frontier corresponding to each objective function;

[0023] Step 11: using a wave search algorithm, optimizing several objective functions of the motor based on the best prognosis metamodel, and obtaining the Pareto front corresponding to each objective function;

[0024] Step 11.1: Set the number of simulated electromagnetic wave particles to n, the number of design variables to d, and initialize the positions of simulated electromagnetic wave particles;

[0025] Specifically: generate s random numbers k1, k2, ...k f ,…,k s , where k f is the fth random number, s≤n, and the simulated electromagnetic wave particle position is initialized according to the method shown in equations (3) and (4):

[0026]

[0027] W i =lb+x i *(ub-lb) (4)

[0028] Among them, x f is a random number k f The value after homogenization, x i * is a value randomly selected from the uniformed values ​​of random numbers, ub is the upper bound of the search space, lb is the lower bound of the search space, and W i is the position of the ith simulated electromagnetic wave particle;

[0029] Step 11.2: Use matrix W to represent the positions of all simulated electromagnetic wave particles;

[0030]

[0031] Among them, W i =[W i1 ,Wi2 ,…W id ] is the position of the ith simulated electromagnetic wave particle, each simulated electromagnetic wave particle position represents a set of solutions for the design variables in the objective function, W ij represents the solution of the jth design variable in the position of the i-th simulated electromagnetic wave particle;

[0032] Step 11.3: Set the current number of iterations t = 0 and the maximum number of iterations to t max ;

[0033] Step 11.4: Calculate the fitness of each simulated electromagnetic wave particle at present;

[0034] The fitness of the simulated electromagnetic wave particle population is expressed as a vector, which includes the fitness of each simulated electromagnetic wave particle:

[0035] F=[f([W 11 ,W 12 ,…,W 1d ]); f([W 21 ,W 22 ,…,W 2d ]); f([W n1 ,W n2 ,…W nd ])] (6)

[0036] Among them, F is the fitness of the simulated electromagnetic wave particle population, f([W i1 ,W i2 ,…W id ]) is the fitness of the ith simulated electromagnetic wave particle;

[0037] The fitness of the simulated electromagnetic wave particle is calculated as follows: when the objective function is to find the maximum value, the fitness of the objective function is the value of the objective function when the design variable is solved in the position of the simulated electromagnetic wave particle; when the objective function is to find the minimum value, the fitness of the objective function is the reciprocal of the value of the objective function when the design variable is solved in the position of the simulated electromagnetic wave particle; finally, the fitness of each objective function is weighted and summed to obtain the fitness of the simulated electromagnetic wave particle;

[0038] Step 11.5: Sort all simulated electromagnetic wave particle positions according to the size of fitness, and take the simulated electromagnetic wave particle position with the largest current fitness as the current best position;

[0039] Step 11.6: Update the current positions of all simulated electromagnetic wave particles within the search space using equations (7) and (8), calculate the fitness of the updated simulated electromagnetic wave particles and sort them according to the fitness, and update the current best position and the global best position;

[0040] The method for updating the current best position and the global best position is: updating the simulated electromagnetic wave particle position with the largest current fitness as the current best position, and taking the current best position and the last best position with the larger fitness value as the global best position;

[0041]

[0042] in, is the calculated new position of the i-th simulated electromagnetic wave particle, W min is a vector composed of the minimum values ​​of each dimension of the matrix W, r1 is a random number between 0 and 1, W max is the vector of the maximum values ​​of each dimension of matrix W, W i ' is the updated position of the ith simulated electromagnetic wave particle, is the calculated new position of the i-th simulated electromagnetic wave particle The fitness of the corresponding simulated electromagnetic wave particle, f mean is the average fitness of all simulated electromagnetic wave particles;

[0043] Step 11.7: simulate the electromagnetic wave to diffuse outward, rearrange the simulated electromagnetic wave particle position according to the degree of proximity to the global optimal position, and update the current positions of all simulated electromagnetic wave particles according to equations (9) and (10), calculate the fitness of the updated simulated electromagnetic wave particles and sort them according to the fitness, and update the current optimal position and the global optimal position; the degree of proximity is the difference in fitness;

[0044]

[0045]

[0046] Where W i " is the position of the ith simulated electromagnetic wave particle updated after the simulated electromagnetic wave diffuses outward, is the calculated new position of the i-th simulated electromagnetic wave particle, W best is the global optimal position, W li For the matrix W, best The degree of proximity of the rearranged matrix, m i is a column vector whose elements obey normal distribution and are arranged in order, σ is the waveform size control coefficient, f max To simulate the maximum fitness of the electromagnetic wave particle position, is the calculated new position of the i-th simulated electromagnetic wave particle The fitness of the corresponding simulated electromagnetic wave particle;

[0047] Step 11.8: simulate the reflected electromagnetic wave. When the position of the simulated electromagnetic wave particle whose fitness is lower than the set threshold value encounters an obstacle and reflects to the global optimal position, the position of the simulated electromagnetic wave particle whose fitness is lower than the set threshold value is updated according to equations (11) and (12). The positions of the remaining simulated electromagnetic wave particles continue to simulate the outward diffusion of electromagnetic waves according to the method in step 11.7. The fitness of the updated simulated electromagnetic wave particles is calculated and sorted according to the fitness value. The current optimal position and the global optimal position are updated.

[0048]

[0049] in, is the calculated new position of the ith simulated electromagnetic wave particle, β is the reflection intensity coefficient, r2 is a random number between 0 and 1, and n w2 is the number of simulated electromagnetic wave particles that simulate the reflected electromagnetic wave, W fi is the matrix after W is rearranged in order of fitness value from small to large, W″′ i is the updated position of the ith simulated electromagnetic wave particle, is the calculated new position of the i-th simulated electromagnetic wave particle The fitness of the corresponding simulated electromagnetic wave particle;

[0050] Step 11.9: simulate the reception of electromagnetic waves, update the current positions of all simulated electromagnetic wave particles according to equations (13) and (14), calculate the fitness of the updated simulated electromagnetic wave particle positions and sort them according to the fitness, and update the current best position and the global best position;

[0051]

[0052] In the formula, is the calculated new position of the ith simulated electromagnetic wave particle, δ is the receiving coefficient, r3, r4, r5, r6 are random numbers between 0 and 1, is a random number that follows a normal distribution, n w3 is the number of simulated electromagnetic wave particles for simulating the received electromagnetic wave, λ is the correction factor, is the historical optimal position obtained by convolution, W i "" is the updated position of the ith simulated electromagnetic wave particle, is the calculated new position of the i-th simulated electromagnetic wave particle The corresponding fitness of the simulated electromagnetic wave particle, is the updated position W of the i-th simulated electromagnetic wave particle i "'The corresponding fitness of the simulated electromagnetic wave particle is l1=l2=...=l t =1 / t;

[0053] Step 11.10: Update the current positions of all simulated electromagnetic wave particles according to equations (16) and (17) through step size experiments combined with the fitted gradient descent method to obtain the final global optimal position of the current iteration round;

[0054] g i =(f(W i ””+ε)-f(W i ””-ε)) / 2ε (16)

[0055] W i ””'=W i ””-α·g i (17)

[0056] In the formula, g i is the gradient, ε is 10 -6 , W i ””’ is the updated position of the ith simulated electromagnetic wave particle, α is the step length, which is determined by step length experiments;

[0057] The step length test is: when the number of iterations is 0, α is the set initial value, and for the t-th iteration, the simulated electromagnetic wave particle position W at the t-1th iteration is calculated. i ””' and simulated electromagnetic wave particle position W i "" fitness, if the simulated electromagnetic wave particle position W i The fitness of ""' is greater than the simulated electromagnetic wave particle position W i ""'s fitness, then α=α t-1 / c, otherwise α=α t-1 *c, c is the scaling factor, α t-1 is the step size at the t-1th iteration;

[0058] Step 11.11: Determine whether the current number of iterations reaches the maximum number of iterations. If so, output the final global optimal position of all iterations to obtain the Pareto frontier. Otherwise, return to step 11.6;

[0059] Step 12: Find and select a set of optimal solutions from the Pareto fronts obtained in Step 10 and Step 11 for ANSYS finite element verification, and take the set of optimal solutions with the best verification effect as the final optimal solution;

[0060] The specific operation of selecting the optimal solution is as follows: set a reference point, evenly select points on the Pareto front, and the point closest to the reference point is taken as the optimal solution.

[0061] The present invention uses MotorCAD to analyze the multi-disciplinary performance of the motor and combines ANSYS optiSLang sensitivity analysis to perform multi-objective multi-level optimization. The intelligent optimization method has the following advantages:

[0062] (1) It combines the advantages of MotorCAD algorithm for rapid analysis of motor performance with the advantages of ANSYS optiSLang integrated simulation process building platform and multidisciplinary design optimization.

[0063] (2) Use Python customized scripts to perform parallel calculations on MotorCAD, automatically call and analyze data results, and reduce computing costs.

[0064] (3) A self-programming algorithm is used to optimize the motor design and compared with the solutions obtained in the software optimization algorithm library to increase the accuracy of the optimization. BRIEF DESCRIPTION OF THE DRAWINGS

[0065] Figure 1 Flow chart of the conventional motor optimization method according to an embodiment of the present invention;

[0066] Figure 2 Embodiment of the present invention Flow chart of the method of the present invention. DETAILED DESCRIPTION

[0067] like Figure 1 The figure shows a flow chart of a traditional motor optimization method. The motor multi-level and multi-objective joint optimization design method based on MotorCAD and ANSYS in the present invention is based on the rapid calculation of motor performance by MotorCAD, and then the multi-objective joint optimization design is performed after the best prognosis meta-model is generated hierarchically based on ANSYS optiSLang sensitivity analysis.

[0068] The following Figure 2 The following is a multi-objective optimization design of a permanent magnet synchronous motor with three objectives: output torque, waveform distortion rate, and efficiency. The multi-level and multi-objective joint optimization design method of the motor based on MotorCAD and ANSYS is described in detail:

[0069] Step 1: Analyze the topology and working scenarios of the motor, and determine several design variables, several objective functions and constraints;

[0070] In this embodiment, the design variable is determined as the ratio of the inner and outer diameters of the stator R st , core length l, magnetic steel thickness H, tooth width W th 、Rotor pole arc coefficient α p , the objective function is the output torque T out , torque ripple T ri , efficiency η; where the constraints stipulate that the output torque is greater than a threshold, the torque ripple is less than a threshold, and the efficiency is greater than a threshold;

[0071] The output torque T out , torque ripple T riThe calculation formula of the efficiency η is as follows:

[0072] T out =p[ψ f i q +(L d -L q )i d i q ] (18)

[0073]

[0074] Where p is the pole pair number, ψ f is the rotor permanent magnet flux, L d is the d-axis inductance, L q is the q-axis inductance, i d is the d-axis current, i q is the q-axis current; T max , T min and T av Respectively represent the maximum, minimum and average value of the output torque; P out , P Fe and P Cu They represent output loss, armature iron loss and armature copper loss respectively;

[0075] Step 2: According to the motor stator and rotor structural parameters, winding settings, and material selection, use MotorCAD to perform parametric modeling on the motor to obtain a motor parametric model;

[0076] In this embodiment, the motor is parametrically modeled using MotorCAD according to the size and structure parameters of the motor stator and rotor, the size of the motor, the type of bearings, the definition of the winding and the material selection;

[0077] Step 3: Set the solution conditions of the motor parametric model in MotorCAD, and solve the motor's electric field, magnetic field, and temperature field under multidisciplinary coupling to obtain the electromagnetic characteristics, thermal characteristics, and mechanical characteristics of the motor;

[0078] In this implementation, the solution conditions include the excitation source, the driving simulation conditions, and the field coupling conditions;

[0079] Step 4: Customize the Python script to create the interface between ANSYSoptiSLang and MotorCAD based on the determined design variables, objective function and constraints;

[0080] Step 5: Based on the interface between ANSYS optiSLang and MotorCAD, ANSYS optiSLang uses Python scripts to call the electromagnetic characteristics, thermal characteristics, and mechanical properties of the motor obtained by MotorCAD, and extracts several determined design variables, several objective functions, and constraints;

[0081] Step 6: According to the sensitivity of each design variable relative to several objective functions, the sensitivity coefficient of each design variable is calculated respectively, and a strong sensitivity layer and a weak sensitivity layer are established respectively according to the strength of the sensitivity coefficient;

[0082] In step 6, in order to accurately analyze the impact of each design variable on the objective function, a suitable hierarchical optimization model is established, and the sensitivity analysis method is used to analyze the impact of each design variable on the objective function. The sensitivity coefficient S is introduced com (x a ), the weight coefficient method is used to comprehensively evaluate each optimization variable, which can be described as:

[0083] S com (x a )=w1|S out (x a )|+w2|S ri (x a |)+w3|S η (x a )| (21)

[0084] w1+w2+w3=1 (22)

[0085] Among them, x a is the ath design variable, S com (x a ) is the sensitivity coefficient of the ath design variable, S out (x a ), S ri (x a ), S η (x a ) are the design variables x a Output torque T for the objective function out , torque ripple T ri , efficiency η, w1, w2, w3 are the output torque T out , torque ripple T ri , the weight coefficient of efficiency η;

[0086] In this embodiment, a design variable with a sensitivity coefficient greater than 0.2 is a strong sensitive layer, and a design variable with a sensitivity coefficient less than 0.2 is a weak sensitive layer.

[0087] Step 7: Sample data are obtained by sampling the strong sensitivity layer and the weak sensitivity layer respectively, and the optimal prognosis metamodel (Metamodel of Optimal Prognosis, MOP) of the strong sensitivity layer and the weak sensitivity layer are established respectively based on the sample data of the strong sensitivity layer and the weak sensitivity layer;

[0088] Establish the best prognosis metamodel for parameter reduction, fully explore the potential feasible solutions in the design space, find the most suitable predictive best prognosis metamodel, which can be used to quickly evaluate the responses corresponding to different parameters, and generate competing best prognosis metamodels according to different strategies after exploring each response in the sample space based on adaptive sampling.

[0089] Step 8: Use the prognostic coefficient to evaluate the prediction quality of the best prognostic meta-model in the strong sensitivity layer and the weak sensitivity layer respectively;

[0090] Specifically, for the strong sensitive layer, if the prognosis coefficient is greater than or equal to the set threshold, it is considered that the accuracy meets the requirement, and step 9 is executed. If it is less than the set threshold, it is considered that the accuracy does not meet the requirement, then increase the sample points and return to step 7. For the weak sensitive layer, if the prognosis coefficient is greater than or equal to the set threshold, it is considered that the accuracy meets the requirement, and go to step 9. If it is less than the set threshold, increase the sample points and return to step 7.

[0091] The coefficient of prognosis (CoP) is a measure of whether the best prognostic metamodel can effectively predict unknown data points. The closer the value is to 1, the more accurate the model prediction is;

[0092]

[0093] Among them, RMSE pred is the root mean square error, which is equivalent to the squared error of the prediction error. h is the actual value, is the estimated value, h is the number of the data point, N is the total number of data points, cv is the ratio of the standard deviation of the data point to the mean, y is a random variable, τ 2 y is the variance of the random variable y;

[0094] In this embodiment, if the prognosis coefficient of the strong sensitive layer is greater than 0.9, step 9 is executed to perform the next step of optimization algorithm optimization; if the prognosis coefficient of the weak sensitive layer is greater than 0.8, step 9 is executed to perform the next step of optimization algorithm calculation; if the accuracy does not meet the requirements, return to step 7 to increase sample points within the entire parameter range and then continue to step 8;

[0095] Step 9: Output the best prognostic meta-model that meets the prediction quality requirements;

[0096] Step 10: In the ANSYS optiSLang algorithm library, select the evolutionary algorithm or particle swarm algorithm to optimize several objective functions of the motor based on the best prognosis metamodel to obtain the Pareto frontier corresponding to each objective function;

[0097] Step 11: using a wave search algorithm, optimizing several objective functions of the motor based on the best prognosis metamodel, and obtaining the Pareto front corresponding to each objective function;

[0098] Step 11.1: Set the number of simulated electromagnetic wave particles to n, the number of design variables to d, and initialize the positions of simulated electromagnetic wave particles;

[0099] Specifically: generate s random numbers k1, k2, ...k f ,…,k s , where k f is the fth random number, s≤n, and the simulated electromagnetic wave particle position is initialized according to the method shown in equations (25) and (26):

[0100]

[0101] Among them, x f is a random number k f The value after homogenization, x i * is a value randomly selected from the uniformed values ​​of random numbers, ub is the upper bound of the search space, lb is the lower bound of the search space, and W i is the position of the ith simulated electromagnetic wave particle;

[0102] Step 11.2: Use matrix W to represent the positions of all simulated electromagnetic wave particles;

[0103]

[0104] Among them, W i =[W i1 ,W i2 ,…W id ] is the position of the ith simulated electromagnetic wave particle, each simulated electromagnetic wave particle position represents a set of solutions for the design variables in the objective function, W ij represents the solution of the jth design variable in the position of the i-th simulated electromagnetic wave particle;

[0105] Step 11.3: Set the current number of iterations t = 0 and the maximum number of iterations to t max ;

[0106] Step 11.4: Calculate the fitness of each simulated electromagnetic wave particle at present;

[0107] The fitness of the simulated electromagnetic wave particle population is expressed as a vector, which includes the fitness of each simulated electromagnetic wave particle:

[0108] F=[f([W 11 ,W 12 ,…,W 1d ]); f([W 21 ,W 22 ,…,W 2d ]); f([W n1 ,W n2 ,…W nd ])] (28)

[0109] Among them, F is the fitness of the simulated electromagnetic wave particle population, [f([W i1 ,W i2 ,…W id ]) is the fitness of the ith simulated electromagnetic wave particle;

[0110] The fitness of the simulated electromagnetic wave particle is calculated as follows: when the objective function is to find the maximum value, the fitness of the objective function is the value of the objective function when the design variable is solved in the position of the simulated electromagnetic wave particle; when the objective function is to find the minimum value, the fitness of the objective function is the reciprocal of the value of the objective function when the design variable is solved in the position of the simulated electromagnetic wave particle; finally, the fitness of each objective function is weighted and summed to obtain the fitness of the simulated electromagnetic wave particle;

[0111] Step 11.5: Sort all simulated electromagnetic wave particle positions according to the size of fitness, and take the simulated electromagnetic wave particle position with the largest current fitness as the current best position;

[0112] Step 11.6: Update the current positions of all simulated electromagnetic wave particles within the search space using equations (29) and (30), calculate the fitness of the updated simulated electromagnetic wave particles and sort them according to the fitness, and update the current best position and the global best position;

[0113] The method for updating the current best position and the global best position is: updating the simulated electromagnetic wave particle position with the largest current fitness as the current best position, and taking the current best position and the last best position with the larger fitness value as the global best position;

[0114] The advantage of formula (29) is that the search range can be gradually narrowed, and the newly generated points are generated within the narrowed range, which improves the search efficiency. Formula (30) is an improved greedy mechanism that controls the position of the population near the global optimal position;

[0115]

[0116] in, is the calculated new position of the i-th simulated electromagnetic wave particle, W min is a vector composed of the minimum values ​​of each dimension of the matrix W, r1 is a random number between 0 and 1, W max is the vector of the maximum values ​​of each dimension of matrix W, W i ' is the updated position of the ith simulated electromagnetic wave particle, is the calculated new position of the i-th simulated electromagnetic wave particle The fitness of the corresponding simulated electromagnetic wave particle, f mean is the average fitness of all simulated electromagnetic wave particles;

[0117] Step 11.7: simulate the electromagnetic wave to diffuse outward, rearrange the simulated electromagnetic wave particle positions according to the degree of proximity between them and the global optimal position, and update all current simulated electromagnetic wave particle positions according to equations (31) and (32), calculate the fitness of the updated simulated electromagnetic wave particle positions and sort them according to the fitness, and update the current optimal position and the global optimal position; the degree of proximity is the difference in fitness;

[0118] The function of formula (31) is to simulate the outward diffusion of electromagnetic waves, reduce the possibility of falling into the local optimum, and improve the search efficiency. Formula (32) is an improved greedy mechanism, which makes the particle position no worse than the current particle position when the particle position fluctuates outward;

[0119]

[0120] Where W i " is the position of the ith simulated electromagnetic wave particle updated after the simulated electromagnetic wave diffuses outward, is the calculated new position of the i-th simulated electromagnetic wave particle, W best is the global optimal position, W li For the matrix W, best The degree of proximity of the rearranged matrix, m i is a column vector whose elements obey normal distribution and are arranged in order, σ is the waveform size control coefficient, f max To simulate the maximum fitness of the electromagnetic wave particle position, is the calculated new position of the i-th simulated electromagnetic wave particle The fitness of the corresponding simulated electromagnetic wave particle;

[0121] Step 11.8: simulate the reflected electromagnetic wave. When the position of the simulated electromagnetic wave particle whose fitness is lower than the set threshold value encounters an obstacle and reflects to the global optimal position, the position of the simulated electromagnetic wave particle whose fitness is lower than the set threshold value is updated according to equations (33) and (34). The positions of the remaining simulated electromagnetic wave particles continue to simulate the outward diffusion of electromagnetic waves according to the method in step 11.7. The fitness of the updated simulated electromagnetic wave particles is calculated and sorted according to the fitness value. The current optimal position and the global optimal position are updated.

[0122] Formula (33) simulates the situation that some particles (particles with lower fitness values) encounter obstacles and move towards W. best Reflection, while the remaining particles (particles with higher fitness values) move away from W best Continue to spread outward;

[0123]

[0124] in, is the calculated new position of the ith simulated electromagnetic wave particle, β is the reflection intensity coefficient, r2 is a random number between 0 and 1, and n w2 is the number of simulated electromagnetic wave particles that simulate the reflected electromagnetic wave, W fi is the matrix after W is rearranged in order of fitness value from small to large, W″′ i is the updated position of the ith simulated electromagnetic wave particle, is the calculated new position of the i-th simulated electromagnetic wave particle The fitness of the corresponding simulated electromagnetic wave particle;

[0125] Step 11.9: simulate the reception of electromagnetic waves, update the current positions of all simulated electromagnetic wave particles according to equations (35) and (36), calculate the fitness of the updated simulated electromagnetic wave particle positions and sort them according to the fitness, and update the current best position and the global best position;

[0126] Receiving electromagnetic waves; It simulates that the radar usually receives electromagnetic waves normally, but sometimes it is interfered and needs to be corrected and processed, as shown in formula (35); Its function is to make the particle group search in the current optimal direction;

[0127]

[0128] In the formula, is the calculated new position of the ith simulated electromagnetic wave particle, δ is the receiving coefficient, r3, r4, r5, r6 are random numbers between 0 and 1, is a random number that follows a normal distribution, n w3 is the number of simulated electromagnetic wave particles for simulating the received electromagnetic wave, λ is the correction factor, is the historical optimal position obtained by convolution, W i"" is the updated position of the ith simulated electromagnetic wave particle, is the calculated new position of the i-th simulated electromagnetic wave particle The corresponding fitness of the simulated electromagnetic wave particle, is the updated position W of the i-th simulated electromagnetic wave particle i "'The corresponding fitness of the simulated electromagnetic wave particle is l1=l2=...=l t =1 / t;

[0129] Step 11.10: According to equations (38) and (39), the positions of all currently simulated electromagnetic wave particles are updated by step size experiments combined with the fitted gradient descent method to obtain the final global optimal position of the current iteration round;

[0130] g i =(f(W i ””+ε)-f(W i ””-ε)) / 2ε (38)

[0131] W i ””'=W i ””-α·g i (39)

[0132] In the formula, g i is the gradient, ε is 10 -6 , W i ””’ is the updated position of the ith simulated electromagnetic wave particle, α is the step length, which is determined by step length experiments;

[0133] The step length test is: when the number of iterations is 0, α is the set initial value, and for the t-th iteration, the simulated electromagnetic wave particle position W at the t-1th iteration is calculated. i ””' and simulated electromagnetic wave particle position W i "" fitness, if the simulated electromagnetic wave particle position W i The fitness of ""' is greater than the simulated electromagnetic wave particle position W i ""'s fitness, then α=α t-1 / c, otherwise α=α t-1 *c, c is the scaling factor, α t-1 is the step size at the t-1th iteration;

[0134] Step 11.11: Determine whether the current number of iterations reaches the maximum number of iterations. If so, output the final global optimal position of all iterations to obtain the Pareto frontier. Otherwise, return to step 11.6;

[0135] Step 12: Find and select a set of optimal solutions from the Pareto fronts obtained in Step 10 and Step 11 for ANSYS finite element verification, and take the set of optimal solutions with the best verification effect as the final optimal solution;

[0136] The specific operation of selecting the optimal solution is as follows: set a reference point, evenly select points on the Pareto front, and the point closest to the reference point is taken as the optimal solution.

[0137] The method of the present invention can effectively solve the complex multi-objective optimization design problem of the motor and quickly find the global optimal design solution.

Claims

1. A motor multi-level multi-objective joint optimization design method based on MotorCAD and ANSYS, characterized in that: Here are the steps: Step 1: Analyze the topology and working scenarios of the motor, and determine several design variables, several objective functions and constraints; Step 2: According to the motor stator and rotor structural parameters, winding settings, and material selection, use MotorCAD to perform parametric modeling on the motor to obtain a motor parametric model; Step 3: Set the solution conditions of the motor parametric model in MotorCAD, and solve the motor's electric field, magnetic field, and temperature field under multidisciplinary coupling to obtain the electromagnetic characteristics, thermal characteristics, and mechanical characteristics of the motor; Step 4: Customize the Python script to create the interface between ANSYSoptiSLang and MotorCAD based on the determined design variables, objective function and constraints; Step 5: Based on the interface between ANSYS optiSLang and MotorCAD, ANSYS optiSLang uses Python scripts to call the electromagnetic characteristics, thermal characteristics, and mechanical properties of the motor obtained by MotorCAD, and extracts several determined design variables, several objective functions, and constraints; Step 6: According to the sensitivity of each design variable relative to several objective functions, the sensitivity coefficient of each design variable is calculated respectively, and a strong sensitivity layer and a weak sensitivity layer are established respectively according to the strength of the sensitivity coefficient; Step 7: Sampling the strong sensitivity layer and the weak sensitivity layer to obtain sample data respectively, and establishing the best prognostic meta-models of the strong sensitivity layer and the weak sensitivity layer respectively based on the sample data of the strong sensitivity layer and the weak sensitivity layer; Step 8: Use the prognostic coefficient to evaluate the prediction quality of the best prognostic meta-model in the strong sensitivity layer and the weak sensitivity layer respectively; Step 9: Output the best prognostic meta-model that meets the prediction quality requirements; Step 10: In the ANSYS optiSLang algorithm library, select the evolutionary algorithm or particle swarm algorithm to optimize several objective functions of the motor based on the best prognosis metamodel to obtain the Pareto frontier corresponding to each objective function; Step 11: using a wave search algorithm, optimizing several objective functions of the motor based on the best prognosis metamodel, and obtaining the Pareto front corresponding to each objective function; Step 12: Find and select a set of optimal solutions from the Pareto fronts obtained in steps 10 and 11 for ANSYS finite element verification, and take the set of optimal solutions with the best verification effect as the final optimal solution.

2. The motor multi-level multi-objective joint optimization design method based on MotorCAD and ANSYS according to claim 1, characterized in that: The sensitivity coefficient of each design variable described in step 6 is calculated as: S com (x a )=w1|S1(x a )|+w2|S2(x a )|+…+w n |S m (x a )|| (1) w1+w2+…+w m =1 (2) Among them, x a is the ath design variable, S com (x a ) is the sensitivity coefficient of the a-th design variable, w m is the weight coefficient of the mth objective function, m is the number of objective functions, S m (x a ) Design variable x a Sensitivity factor to the mth objective function.

3. The motor multi-level multi-objective joint optimization design method based on MotorCAD and ANSYS according to claim 1, characterized in that: The step 8 is specifically as follows: for the strong sensitive layer, if the prognosis coefficient is greater than or equal to the set threshold, it is considered that the accuracy meets the requirement, and step 9 is executed; if it is less than the set threshold, it is considered that the accuracy does not meet the requirement, then sample points are added and the process returns to step 7; for the weak sensitive layer, if the prognosis coefficient is greater than or equal to the set threshold, it is considered that the accuracy meets the requirement, and the process returns to step 9; if it is less than the set threshold, then sample points are added and the process returns to step 7.

4. The motor multi-level multi-objective joint optimization design method based on MotorCAD and ANSYS according to claim 1, characterized in that: The step 11 specifically includes: Step 11.1: Set the number of simulated electromagnetic wave particles to n, the number of design variables to d, and initialize the positions of simulated electromagnetic wave particles; Specifically: generate s random numbers k1, k2, ...k f ,…,k s , where k f is the fth random number, s≤n, and the simulated electromagnetic wave particle position is initialized according to the method shown in equations (3) and (4): W i =lb+x i * (ub-lb) (4) Among them, x f is a random number k f The value after homogenization, x i * is a value randomly selected from the uniformed values ​​of random numbers, ub is the upper bound of the search space, lb is the lower bound of the search space, and W i is the position of the i-th simulated electromagnetic wave particle; Step 11.2: Use matrix W to represent the positions of all simulated electromagnetic wave particles; Among them, W i =[W i1 ,W i2 ,…W id ] is the position of the ith simulated electromagnetic wave particle, each simulated electromagnetic wave particle position represents a set of solutions for the design variables in the objective function, W ij represents the solution of the jth design variable in the position of the i-th simulated electromagnetic wave particle; Step 11.3: Set the current number of iterations t = 0 and the maximum number of iterations to t max ; Step 11.4: Calculate the fitness of each simulated electromagnetic wave particle at present; The fitness of the simulated electromagnetic wave particle population is expressed as a vector, which includes the fitness of each simulated electromagnetic wave particle: F=[f([W 11 ,W 12 ,…,W 1d ]);f([W 21 ,W 22 ,…,W 2d ]);f([W n1 ,W n2 ,…W nd ])] (6) Among them, F is the fitness of the simulated electromagnetic wave particle population, f([W i1 ,W i2 ,…W id ]) is the fitness of the ith simulated electromagnetic wave particle; The fitness of the simulated electromagnetic wave particle is calculated as follows: when the objective function is to find the maximum value, the fitness of the objective function is the value of the objective function when the design variable is solved in the position of the simulated electromagnetic wave particle; when the objective function is to find the minimum value, the fitness of the objective function is the reciprocal of the value of the objective function when the design variable is solved in the position of the simulated electromagnetic wave particle; finally, the fitness of each objective function is weighted and summed to obtain the fitness of the simulated electromagnetic wave particle; Step 11.5: Sort all simulated electromagnetic wave particle positions according to the size of fitness, and take the simulated electromagnetic wave particle position with the largest current fitness as the current best position; Step 11.6: Update the current positions of all simulated electromagnetic wave particles within the search space using equations (7) and (8), calculate the fitness of the updated simulated electromagnetic wave particles and sort them according to the fitness, and update the current best position and the global best position; The method for updating the current best position and the global best position is: updating the simulated electromagnetic wave particle position with the largest current fitness as the current best position, and taking the current best position and the last best position with the larger fitness value as the global best position; in, is the calculated new position of the i-th simulated electromagnetic wave particle, W min is a vector composed of the minimum values ​​of each dimension of the matrix W, r1 is a random number between 0 and 1, W max is the vector of the maximum values ​​of each dimension of matrix W, W i ' is the updated position of the ith simulated electromagnetic wave particle, is the calculated new position of the i-th simulated electromagnetic wave particle The corresponding fitness of the simulated electromagnetic wave particle, f mean is the average fitness of all simulated electromagnetic wave particles; Step 11.7: simulate the electromagnetic wave to diffuse outward, rearrange the simulated electromagnetic wave particle position according to the degree of proximity to the global optimal position, and update the current positions of all simulated electromagnetic wave particles according to equations (9) and (10), calculate the fitness of the updated simulated electromagnetic wave particles and sort them according to the fitness, and update the current optimal position and the global optimal position; the degree of proximity is the difference in fitness; Where W i " is the position of the ith simulated electromagnetic wave particle updated after the simulated electromagnetic wave diffuses outward, is the calculated new position of the i-th simulated electromagnetic wave particle, W best is the global optimal position, W li For the matrix W, best The degree of proximity of the rearranged matrix, m i is a column vector whose elements obey normal distribution and are arranged in order, σ is the waveform size control coefficient, f max To simulate the maximum fitness of the electromagnetic wave particle position, is the calculated new position of the i-th simulated electromagnetic wave particle The fitness of the corresponding simulated electromagnetic wave particle; Step 11.8: simulate the reflected electromagnetic wave. When the position of the simulated electromagnetic wave particle whose fitness is lower than the set threshold value encounters an obstacle and reflects to the global optimal position, the position of the simulated electromagnetic wave particle whose fitness is lower than the set threshold value is updated according to equations (11) and (12). The positions of the remaining simulated electromagnetic wave particles continue to simulate the outward diffusion of electromagnetic waves according to the method in step 11.

7. The fitness of the updated simulated electromagnetic wave particles is calculated and sorted according to the fitness value. The current optimal position and the global optimal position are updated. in, is the calculated new position of the ith simulated electromagnetic wave particle, β is the reflection intensity coefficient, r2 is a random number between 0 and 1, and n w2 is the number of simulated electromagnetic wave particles that simulate the reflected electromagnetic wave, W fi is the matrix after W is rearranged in order of fitness value from small to large, W' i " is the updated position of the ith simulated electromagnetic wave particle, is the calculated new position of the i-th simulated electromagnetic wave particle The fitness of the corresponding simulated electromagnetic wave particle; Step 11.9: simulate the reception of electromagnetic waves, update the current positions of all simulated electromagnetic wave particles according to equations (13) and (14), calculate the fitness of the updated simulated electromagnetic wave particle positions and sort them according to the fitness, and update the current best position and the global best position; In the formula, is the calculated new position of the ith simulated electromagnetic wave particle, δ is the receiving coefficient, r3, r4, r5, r6 are random numbers between 0 and 1, ζ is a random number that obeys the normal distribution, n w3 is the number of simulated electromagnetic wave particles for simulating the received electromagnetic wave, λ is the correction factor, is the historical optimal position obtained by convolution, W i "" is the updated position of the ith simulated electromagnetic wave particle, is the calculated new position of the i-th simulated electromagnetic wave particle The corresponding fitness of the simulated electromagnetic wave particle, is the updated position W of the i-th simulated electromagnetic wave particle i "'The corresponding fitness of the simulated electromagnetic wave particle is l1=l2=...=l t =1 / t; Step 11.10: Update the current positions of all simulated electromagnetic wave particles according to equations (16) and (17) through step size experiments combined with the fitted gradient descent method to obtain the final global optimal position of the current iteration round; g i =(f(W i ””+ε)-f(W i ””-ε)) / 2ε (16) W i ””'=W i ””-α·g i (17) In the formula, g i is the gradient, ε is 10 -6 , W i ””’ is the updated position of the ith simulated electromagnetic wave particle, α is the step length, which is determined by step length experiments; The step length test is: when the number of iterations is 0, α is the set initial value, and for the t-th iteration, the simulated electromagnetic wave particle position W at the t-1th iteration is calculated. i ””' and simulated electromagnetic wave particle position W i "" fitness, if the simulated electromagnetic wave particle position W i The fitness of ""' is greater than the simulated electromagnetic wave particle position W i ""'s fitness, then α=α t-1 / c, otherwise α=α t-1 *c, c is the scaling factor, α t-1 is the step size at the t-1th iteration; Step 11.11: Determine whether the current number of iterations has reached the maximum number of iterations. If so, output the final global optimal position of all iterations to obtain the Pareto frontier. Otherwise, return to step 11.

6.

5. The motor multi-level multi-objective joint optimization design method based on MotorCAD and ANSYS according to claim 1, characterized in that: The specific operation of selecting the optimal solution in step 12 is as follows: set a reference point, evenly select points on the Pareto front, and the point closest to the reference point is taken as the optimal solution.

Images