Large-scale wind power blade fatigue test load optimization method based on hybrid optimization algorithm

Abstract

The invention discloses a large-scale wind power blade fatigue test load optimization method based on a hybrid optimization algorithm, and relates to the technical field of wind power generation equipment testing. The method comprises the following steps: step 1, determining test parameters and physical characteristic parameters according to the structural characteristics of the wind power blade, and defining optimization variables and an optimization variable range; step 2, establishing a blade region segmentation weighting objective function, setting a differentiation weight and a threshold value for each region, and constructing an objective function of a differentiation punishment mechanism; 3, according to the objective function in the step 2, a Bayesian global optimization algorithm is adopted, a Gaussian process proxy model of the objective function is constructed for optimization, and a global optimal solution of an optimization variable is output; and 4, taking the Bayesian optimization result in the step 3 as a starting point, performing local accurate optimization by applying an interior point method, and outputting a final local optimization result of the optimization variable.

Description

Technical Field

[0001] The invention relates to the technical field of wind power generation equipment testing, and in particular to a large wind power blade fatigue test load optimization method based on a hybrid optimization algorithm. Background Art

[0002] With the rapid development of wind power generation technology, fatigue testing of large wind turbine blades has become a key link to ensure the safe and reliable operation of wind power equipment. The traditional wind turbine blade fatigue test method mainly relies on manual experience to adjust the exciter parameters, which is not only time-consuming and labor-intensive, but also with the continuous increase in blade length, wind turbine blade fatigue test load optimization has become a new topic.

[0003] In the existing technology, a single optimization algorithm has many limitations when solving the optimization problem of wind turbine blade fatigue testing, and it is difficult to balance test accuracy and test efficiency. On the one hand, traditional gradient-based optimization methods (such as the steepest descent method, conjugate gradient method, etc.) are prone to fall into local optimality and it is difficult to find the global optimal solution in a complex parameter space; on the other hand, although heuristic algorithms (such as genetic algorithms, particle swarm algorithms, etc.) have certain global search capabilities, they have slow convergence speed, low efficiency, and limited accuracy, and it is difficult to meet the needs of high-precision testing. Therefore, it is urgent to develop a new method that can efficiently and accurately optimize the load distribution of wind turbine blade fatigue testing to improve the test accuracy and efficiency of large wind turbine blades. Summary of the invention

[0004] The embodiment of the present invention provides a large wind turbine blade fatigue test load optimization method based on a hybrid optimization algorithm to solve the problem that it is difficult to balance test accuracy and test efficiency in a traditional large wind turbine blade fatigue test load optimization algorithm.

[0005] In order to solve the above technical problems, an embodiment of the present invention provides a large wind turbine blade fatigue test load optimization method based on a hybrid optimization algorithm, the method comprising: Step 1: According to the structural characteristics of the wind turbine blade, determine the test parameters and physical characteristic parameters, and define the optimization variables and the range of optimization variables. The test parameters include: wind turbine blade length, test position point distribution and target bending moment distribution, the physical characteristic parameters include: blade damping ratio, and the optimization variables include: exciter position, exciter mass and excitation force amplitude.

[0006] Step 2: Establish a segmented weighted objective function for the leaf area, set differentiated weights and thresholds for each area, and construct an objective function for the differentiated penalty mechanism. The leaf is divided into five areas: the root area, the root-middle transition area, the middle area, the middle-tip transition area, and the tip area. The objective function is the accumulation of the penalty mechanism for each area.

[0007] Step 3: Based on the objective function of step 2, a Bayesian global optimization algorithm is used to construct a Gaussian process proxy model of the objective function for optimization, and the global optimal solution of the optimization variables is output. The global optimal solution includes: the global optimal solution of the exciter position, the global optimal solution of the exciter mass, and the global optimal solution of the excitation force amplitude.

[0008] Step 4: Using the Bayesian optimization result of step 3 as the starting point, apply the interior point method to perform local precise optimization and output the final local optimization results of the optimization variables. The final local optimization results of the optimization variables include: the local optimal solution of the exciter position, the local optimal solution of the exciter mass and the local optimal solution of the excitation force amplitude.

[0009] Optionally, step 1 includes step 1.1 and step 1.2.

[0010] Step 1.1: Set the test parameters and physical characteristic parameters according to the structural characteristics of the wind turbine blade.

[0011] Step 1.2: Define the exciter position, exciter mass and excitation force amplitude as optimization variables, and define the variation range of the optimized exciter position, exciter mass and excitation force amplitude.

[0012] Among them, the range of the exciter position change is defined according to the blade length: According to the loading method and loading position, the range of the vibration exciter mass change is defined as , unit: kilogram (kg), according to the loading method, the loading position of the exciter, and the mass of the exciter, the amplitude variation range of the exciter is defined as , unit: Newton (N).

[0013] Optionally, step 2 includes steps 2.1 to 2.4: Step 2.1, divide the blade test area into: root area, root-middle transition area, middle area, middle-tip transition area, tip area.

[0014] Among them, the proportion of the root area to the entire leaf position is: , the root-middle transition area accounts for the proportion of the entire leaf position: , the proportion of the middle area to the entire blade position is: , the proportion of the middle-tip transition area to the entire blade position is: , the tip area accounts for the proportion of the entire blade position: ;in, .

[0015] Step 2.2: Calculate the relative error between the test load and the target load in each area in step 2.1.

[0016] Among them, for any test point selected on the blade , and its relative error is defined as: ; in, Indicates leaf The relative error of the test points is Indicates any test point selected on the blade when the stress ratio is -1 The test load amplitude is Indicates the test point corresponding to the test load amplitude when the stress ratio is -1 given in the test report The target load amplitude.

[0017] Step 2.3: Based on the relative error of each region calculated in step 2.2, differentiated weights and thresholds are assigned to each region, and the relative error penalty formula for each region is constructed as follows: ; in, Indicates the area identifier, which includes: root area root, root-middle transition area root-mid, middle area mid, middle-tip transition area mid-tip, tip area tip, Represents the index set of test points in the corresponding area, Indicates leaf The relative error of the test points is The weight coefficient representing the relative error penalty for the corresponding area exceeding the threshold, represents the tolerance threshold of the corresponding area, Represents the mechanism by which the corresponding region penalizes relative errors that exceed a threshold.

[0018] Step 2.4: After the relative error penalty mechanism of each region is constructed in step 2.3, the first objective function of the total relative error penalty mechanism of all regions of the blade is obtained, which is expressed as: ; in, Represents the first objective function of the total relative error penalty mechanism for all regions, represents the relative error penalty mechanism for root areas exceeding the threshold, represents the relative error penalty mechanism for exceeding the threshold in the root-middle transition region, Represents the relative error penalty mechanism for the central area exceeding the threshold, Relative error penalty mechanism for the mid-to-apex transition region exceeding the threshold, Relative error penalty mechanism for peak areas exceeding the threshold.

[0019] Optionally, step 2.3 may further include step 2.3a.

[0020] Step 2.3a: Special treatment is given to the negative relative error of the root in step 2.3. Only the negative relative error value of the root area is retained and penalized by an exponential function: ; ; in, Represents a set of negative relative error values ​​in the root area, Represents the test point index set in the root area, Represents the weight coefficient of the negative relative error penalty in the root area, Indicates leaf The relative error of the test points is represents the exponential coefficient, Represents the negative relative error penalty mechanism in the root region.

[0021] After executing step 2.3a, use the root area total relative error penalty mechanism Update the relative error penalty mechanism for root area exceeding the threshold in step 2.4 above .

[0022] Among them, the total relative error penalty mechanism for the root area is expressed as: ; Similarly, the first objective function in step 2.4 is updated to the second objective function; ; in, represents the total relative error penalty mechanism in the root area, represents the relative error penalty mechanism for root areas exceeding the threshold, represents the negative relative error penalty mechanism in the root area, Represents the first objective function of the total relative error penalty mechanism for all regions, The second objective function representing the total relative error penalty mechanism for all regions.

[0023] Optionally, step 2.3 may further include step 2.3b.

[0024] Step 2.3b: Introduce an exponential penalty exceeding the threshold of 0.15 for the central region with the largest relative error in step 2.3: ; ; in, Represents the set of relative error values ​​in the middle area that exceed the threshold of 0.15. Indicates leaf The relative error of the test points is Represents the test point index set in the middle area, represents the weight coefficient of the relative error penalty exceeding the threshold of 0.15 in the middle area, represents the exponential coefficient, Relative error penalty mechanism for exceeding the threshold of 0.15 in the middle region.

[0025] After executing step 2.3b, use the total relative error penalty mechanism in the middle area Update the relative error penalty mechanism for the root middle area exceeding the threshold in step 2.4 above .

[0026] Among them, the total relative error penalty mechanism for the central area is expressed as: ; Similarly, the first objective function in step 2.4 is updated to the third objective function; ; in, represents the total relative error penalty mechanism in the middle area, Represents the relative error penalty mechanism for the central area exceeding the threshold, Relative error penalty mechanism for the central region exceeding the threshold of 0.15. Represents the first objective function of the total relative error penalty mechanism for all regions, The third objective function representing the total relative error penalty mechanism for all regions.

[0027] Optionally, after step 2.3 and before step 2.4, step 2.5 is also included.

[0028] Step 2.5: Calculate the load gradients of adjacent test points in all regions, and perform secondary penalties on all gradients to smoothly control the load distribution gradients. The formula is as follows: ; ; Similarly, the first objective function in step 2.4 is updated to the fourth objective function; ; in, , represents the relative error gradient between adjacent test points, represents the total number of test points, Indicates leaf The relative error of the test points is represents the weight coefficient for smoothing the load distribution gradient in all regions, represents the penalty mechanism for smooth control of the load distribution gradient in all regions, Represents the first objective function of the total relative error penalty mechanism for all regions, The fourth objective function representing the total relative error penalty mechanism for all regions.

[0029] Optionally, after step 2.3, step 2.3a, step 2.3b, step 2.5 and step 2.4 are performed in sequence, and the first objective function in step 2.4 is updated to the fifth objective function.

[0030] Step 2.3a: Special treatment is given to the negative relative error of the root in step 2.3. Only the negative relative error value of the root area is retained and penalized by an exponential function: ; ; in, Represents a set of negative relative error values ​​in the root area, Represents the test point index set in the root area, Represents the weight coefficient of the negative relative error penalty in the root area, Indicates leaf The relative error of the test points is represents the exponential coefficient, Represents the negative relative error penalty mechanism in the root region.

[0031] Step 2.3b: Introduce an exponential penalty exceeding the threshold of 0.15 for the central region with the largest relative error in step 2.3: ; ; in, Represents the set of relative error values ​​in the middle area that exceed the threshold of 0.15. Indicates leaf The relative error of the test points is Represents the test point index set in the middle area, represents the weight coefficient of the relative error penalty exceeding the threshold of 0.15 in the middle area, represents the exponential coefficient, Relative error penalty mechanism for exceeding the threshold of 0.15 in the middle region.

[0032] Step 2.5: Calculate the load gradients of adjacent test points in all regions, and perform quadratic penalties on all gradients to smoothly control the load distribution gradients.

[0033] The formula is as follows: ; ; in, , represents the relative error gradient between adjacent test points, represents the total number of test points, Indicates leaf The relative error of the test points is represents the weight coefficient for smoothing the load distribution gradient in all regions, Represents the penalty mechanism for smoothing the gradient of the load distribution in all regions.

[0034] ; in, Represents the first objective function of the total relative error penalty mechanism for all regions, represents the negative relative error penalty mechanism in the root area, Relative error penalty mechanism for the central region exceeding the threshold of 0.15. represents the penalty mechanism for smooth control of the load distribution gradient in all regions, The fifth objective function representing the total relative error penalty mechanism for all regions.

[0035] Optionally, after step 2 and before step 3, step 2.6 is also included.

[0036] Step 2.6: Perform variance control on the loads in all regions, and update the fifth objective function to the sixth objective function.

[0037] The variance control penalty formula is: ; ; in, Indicates the area identification, which includes: root area root, root-middle transition area root-mid, middle area mid, middle-tip transition area mid-tip, tip area tip, Represents the index set of test points in the corresponding area, Indicates the number of elements in a collection. represents the mean of the relative error in the corresponding area, Indicates leaf The relative error of the test points is represents the variance of the relative error in the corresponding area, represents the weight coefficient of the variance in the corresponding region, represents the penalty mechanism for variance control within all regions.

[0038] The fifth objective function in step 2.4 is then updated to the sixth objective function.

[0039] ; in, Represents the first objective function of the total relative error penalty mechanism for all regions, represents the negative relative error penalty mechanism in the root area, Relative error penalty mechanism for the central region exceeding the threshold of 0.15. represents the penalty mechanism for smooth control of the load distribution gradient in all regions, The fifth objective function of the relative error penalty mechanism for all regions is represented by, The sixth objective function representing the total relative error penalty mechanism for all regions.

[0040] Optionally, step 3 includes the following steps 3.1 to 3.3.

[0041] Step 3.1: Set the Bayesian optimization parameters, enable parallel computing, and select the acquisition function you want to improve.

[0042] Among them, the Bayesian optimization parameters include: the maximum number of objective function evaluations is , the maximum running time is hours, the exploration-exploitation balance ratio is , the size of the Gaussian process activity set is ; The expected improvement as a collection function is expressed as: ; in, represents the acquisition function, represents the expected function, represents the currently observed optimal function value, Represents the optimization variable The corresponding objective function value.

[0043] Step 3.2: Use the Gaussian process model as a proxy model for Bayesian optimization to solve the mean of the multivariate normal distribution and the covariance matrix of the multivariate normal distribution.

[0044] Among them, let the objective function is a realization of a Gaussian process, for any finite set of points , function value The mean is , the covariance matrix is The multivariate normal distribution of : ,in, represents the mean of the multivariate normal distribution, Represents the covariance matrix of the multivariate normal distribution.

[0045] Step 3.3: According to the objective function in step 2, initialize the Gaussian process model, execute the Bayesian optimization process, and output the global optimal solution of Bayesian optimization. .

[0046] Among them, the global optimal solution Includes: Vibrator position The global optimal solution and the quality of the shaker The global optimal solution and exciting force amplitude of The global optimal solution of .

[0047] Optionally, step 4 includes steps 4.1 to 4.3.

[0048] Step 4.1: Set the interior point optimization parameters and enable parallel computing.

[0049] Among them, the optimization parameters of the interior point method include: the maximum number of function evaluations is , the maximum number of iterations is , the objective function gradient tolerance is , the variable change step tolerance is , and enable parallel computing.

[0050] Step 4.2: Global optimal solution obtained based on Bayesian optimization in step 3.3 , using the strategy of reducing the search space, the search range of the interior point method is limited to of( ) neighborhood to find the local optimal solution.

[0051] in, is the proportional coefficient, and its value range is , the lower and upper bound expressions of the optimization range are: ; in, and are the lower and upper bounds of the optimization range of the optimization variables of the Bayesian global optimization set in step 1, and Respectively represent the range-limiting parameters, and They respectively represent the lower and upper bounds of the optimization range of the optimization variable of the interior point method local optimization.

[0052] Step 4.3: Use Optimization Toolbox to perform interior point iterative optimization and output the final local optimization results of the optimization variables. , that is, the final local optimization results of the optimization variables include: the local optimal solution of the exciter position, the local optimal solution of the exciter mass and the local optimal solution of the exciting force amplitude.

[0053] Among them, step 4.3 is implemented by steps 4.3a to 4.3c.

[0054] Step 4.3a: By introducing a logarithmic barrier function, the constrained optimization problem is transformed into a series of sub-problems with no constraints or only equality constraints, and the fatigue test parameters of the wind turbine blade are optimized.

[0055] ; ; in, Output variable when the comprehensive objective function takes the minimum value The value of represents the optimization variable vector, Corresponding to the position of the exciter , Vibrator quality and the excitation force amplitude ; is the comprehensive objective function defined, and are the lower and upper bounds of the optimization variable, respectively.

[0056] Step 4.3b: By reducing the obstacle parameter to approach the optimal solution, the interior point method transforms the original optimization problem into: ; in, Output variable when the objective function containing logarithmic obstacle term reaches the minimum value The value of represents the objective function including logarithmic obstacle term, is the comprehensive objective function defined, is the barrier parameter, , gradually decreases as the algorithm iterates. When the obstacle problem is solved Converge to the optimal solution of the original problem, Indicates optimization variables, and Respectively The lower and upper bounds of the optimization variables.

[0057] Step 4.3c: Solve step 4.3b and output the final local optimization result of the optimization variable. .

[0058] Beneficial effects of the present invention: 1. This application combines Bayesian global optimization with interior point method local precise optimization to form a two-stage hybrid optimization strategy, which can not only effectively explore the global parameter space, but also perform high-precision local optimization in local areas, significantly improving the optimization efficiency and result accuracy.

[0059] 2. This application innovatively proposes a regional segmented weighted objective function, sets differentiated weights, thresholds and penalty mechanisms for different functional areas of wind turbine blades, strengthens the processing of negative relative errors in the root area and ultra-high relative errors in the middle area, and effectively improves the problem of traditional methods that it is difficult to balance the load matching accuracy of different areas of the blade, thereby further improving the test accuracy.

[0060] 3. This application ensures the smooth transition of the test load between different regions of the blade and the consistency within the region through load distribution gradient smoothing control and regional variance constraint mechanism, avoids sudden changes in load distribution during the test, reduces the risk of unreasonable stress distribution, and further improves the reliability of the test.

[0061] 4. The Gaussian process proxy model used in this application can find the global optimal solution within a limited number of function evaluations, greatly reducing computational overhead, accelerating the optimization process, and also improving testing efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] Figure 1 One of the flow charts of the optimization method provided in this application; Figure 2 The second flowchart of the optimization method provided for this application; Figure 3 The third flowchart of the optimization method provided in this application; Figure 4 Flow chart 4 of the optimization method provided in this application; Figure 5 Flow chart 5 of the optimization method provided in this application; Figure 6 Flow chart six of the optimization method provided for this application; Figure 7 Flow chart seven of the optimization method provided for this application; Figure 8 These are the comparison diagrams of load optimization results and error result diagrams for the wind turbine blade test of this application, where (a) is the bending moment amplitude comparison diagram, and (b) is the bending moment amplitude comparison error diagram. DETAILED DESCRIPTION

[0063] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0064] Example 1 like Figure 1 As shown, this embodiment discloses a large wind turbine blade fatigue test load optimization method based on a hybrid optimization algorithm, and the method includes the following steps 1 to 4.

[0065] Step 1: According to the structural characteristics of the wind turbine blade, determine the test parameters and physical characteristic parameters, and define the optimization variables and the optimization variable range.

[0066] Among them, the test parameters include: wind turbine blade length, test position point distribution and target bending moment distribution, the physical characteristic parameters include: blade damping ratio, and the optimization variables include: exciter position, exciter mass and excitation force amplitude.

[0067] Optionally, step 1 specifically includes the following steps 1.1 and 1.2.

[0068] Step 1.1: Set the test parameters and physical characteristic parameters according to the structural characteristics of the wind turbine blade.

[0069] Step 1.2: Define the exciter position, exciter mass and excitation force amplitude as optimization variables, and define the variation range of the optimized exciter position, exciter mass and excitation force amplitude.

[0070] Among them, the range of the exciter position change is defined according to the blade length: , unit: meter (m), according to the loading method and loading position, the range of the vibration exciter mass change is defined as , unit: kilogram (kg), according to the loading method, the loading position of the exciter, and the mass of the exciter, the amplitude variation range of the exciter is defined as , unit: Newton (N).

[0071] It should be noted that step 1 may also include the following step 1.3, which is used by the designer to parallelly calculate and run the position, and may be set before step 1.1, before step 1.2, or after step 1.2, and this application does not limit this.

[0072] Step 1.3: Initialize the parallel computing pool and set it to run locally.

[0073] In addition, in addition to being set to run locally, the above step 1 can also be set to run in a distributed network, run on a server, etc., which is not limited in this application.

[0074] Step 2: Establish a segmented weighted objective function for the leaf area, set differentiated weights and thresholds for each area, and construct an objective function for the differentiated penalty mechanism.

[0075] Among them, the blade is divided into five areas, namely: root area, root-middle transition area, middle area, middle-tip transition area, and tip area. The objective function is the accumulation of the penalty mechanism of each area.

[0076] Optional, combined Figure 1 ,like Figure 2 As shown, step 2 includes steps 2.1 to 2.4.

[0077] Step 2.1, divide the blade test area into: root area, root-middle transition area, middle area, middle-tip transition area, tip area.

[0078] Among them, the proportion of the root area to the entire leaf position is: , the root-middle transition area accounts for the proportion of the entire leaf position: , the proportion of the middle area to the entire blade position is: , the proportion of the middle-tip transition area to the entire blade position is: , the tip area accounts for the proportion of the entire blade position: ;in, .

[0079] For example, when , , , The blade test area can be divided into: root area (0~20%), root-middle transition area (20~35%), middle area (35~55%), middle-tip transition area (55~70%), tip area (70~100%). to Determined according to the geometric dimensions of the specific blade.

[0080] Step 2.2: Calculate the relative error between the test load and the target load in each area in step 2.1.

[0081] Among them, for any test point selected on the blade , and its relative error is defined as: ; in, Indicates leaf The relative error of the test points is Indicates any test point selected on the blade when the stress ratio is -1 The test load amplitude is Indicates the test point corresponding to the test load amplitude when the stress ratio is -1 given in the test report The target load amplitude.

[0082] Step 2.3: Based on the relative error of each region calculated in step 2.2, differentiated weights and thresholds are assigned to each region, and a relative error penalty formula for each region is constructed.

[0083] The relative error penalty formula for each region is expressed as: ; in, Indicates the area identifier, which includes: root area root, root-middle transition area root-mid, middle area mid, middle-tip transition area mid-tip, tip area tip, Represents the index set of test points in the corresponding area, Indicates leaf The relative error of the test points is The weight coefficient representing the relative error penalty for the corresponding area exceeding the threshold, represents the tolerance threshold of the corresponding area, Indicates the mechanism by which the corresponding region penalizes the relative error exceeding the threshold. The contexts are the same region identifier and have a one-to-one correspondence.

[0084] Step 2.4: After the relative error penalty mechanism for each region is constructed in step 2.3, the first objective function of the total relative error penalty mechanism for all regions of the blade is obtained.

[0085] The first objective function is expressed as: ; in, Represents the first objective function of the total relative error penalty mechanism for all regions, represents the relative error penalty mechanism for root areas exceeding the threshold, represents the relative error penalty mechanism for exceeding the threshold in the root-middle transition region, Represents the relative error penalty mechanism for the central area exceeding the threshold, Relative error penalty mechanism for the mid-to-apex transition region exceeding the threshold, Relative error penalty mechanism for peak areas exceeding the threshold.

[0086] Optionally, step 2.3 may also include step 2.3a, which further optimizes the root relative error penalty mechanism based on step 2.3 (i.e., considering the impact of negative relative error on the root relative error penalty mechanism), thereby further improving the rationality and accuracy of the relative error penalty mechanism. In practice, it can be selectively executed according to actual needs.

[0087] Step 2.3a: Special treatment is given to the negative relative error of the root in step 2.3. Only the negative relative error value of the root area is retained and penalized by an exponential function: ; ; in, Represents a set of negative relative error values ​​in the root area, Represents the test point index set in the root area, Represents the weight coefficient of the negative relative error penalty in the root area, Indicates leaf The relative error of the test points is represents the exponential coefficient, Represents the negative relative error penalty mechanism in the root region.

[0088] After executing step 2.3a, use the root area total relative error penalty mechanism Update the relative error penalty mechanism for root area exceeding the threshold in step 2.4 above .

[0089] Among them, the total relative error penalty mechanism for the root area is expressed as: ; Similarly, the first objective function in step 2.4 is updated to the second objective function; ; in, represents the total relative error penalty mechanism in the root area, represents the relative error penalty mechanism for root areas exceeding the threshold, represents the negative relative error penalty mechanism in the root area, Represents the first objective function of the total relative error penalty mechanism for all regions, The second objective function representing the total relative error penalty mechanism for all regions.

[0090] Optionally, step 2.3 may also include step 2.3b: step 2.3b is a further optimization of the relative error penalty mechanism for the central region based on step 2.3 (i.e., introducing an exponential penalty exceeding a threshold of 0.15 for the central region with the largest relative error in step 2.3), thereby further improving the rationality and accuracy of the relative error penalty mechanism. In practice, it can be selectively executed according to actual needs.

[0091] Step 2.3b: Introduce an exponential penalty exceeding the threshold of 0.15 for the central region with the largest relative error in step 2.3: ; ; in, Represents the set of relative error values ​​in the middle area that exceed the threshold of 0.15. Indicates leaf The relative error of the test points is Represents the test point index set in the middle area, represents the weight coefficient of the relative error penalty exceeding the threshold of 0.15 in the middle area, represents the exponential coefficient, Relative error penalty mechanism for exceeding the threshold of 0.15 in the middle region.

[0092] After executing step 2.3b, use the total relative error penalty mechanism in the middle area Update the relative error penalty mechanism for the root middle area exceeding the threshold in step 2.4 above .

[0093] Among them, the total relative error penalty mechanism for the central area is expressed as: ; Similarly, the first objective function in step 2.4 is updated to the third objective function; ; in, represents the total relative error penalty mechanism in the middle area, Represents the relative error penalty mechanism for the central area exceeding the threshold, Relative error penalty mechanism for the central region exceeding the threshold of 0.15. Represents the first objective function of the total relative error penalty mechanism for all regions, The third objective function representing the total relative error penalty mechanism for all regions.

[0094] Optional, combined Figure 2 ,like Figure 3As shown, after step 2.3 and before step 2.4, step 2.5 is also included.

[0095] Step 2.5: Calculate the load gradients of adjacent test points in all regions, and perform quadratic penalties on all gradients to smoothly control the load distribution gradients.

[0096] The control formula is as follows: ; ; Similarly, the first objective function in step 2.4 is updated to the fourth objective function; ; in, , represents the relative error gradient between adjacent test points, represents the total number of test points, Indicates leaf The relative error of the test points is represents the weight coefficient for smoothing the load distribution gradient in all regions, represents the penalty mechanism for smooth control of the load distribution gradient in all regions, Represents the first objective function of the total relative error penalty mechanism for all regions, The fourth objective function representing the total relative error penalty mechanism for all regions.

[0097] Optional, combined Figure 1 , Figure 2 ,like Figure 4 As shown, after step 2.3, step 2.3a, step 2.3b, step 2.5 and step 2.4 are executed in sequence (that is, to further improve the accuracy of the penalty mechanism, the optimization steps 2.3a and 2.3b of the two penalty mechanisms are both executed), and the first objective function in step 2.4 is updated to the fifth objective function.

[0098] Step 2.3a: Special treatment is given to the negative relative error of the root in step 2.3, and only the negative relative error value of the root area is retained, and penalized by an exponential function.

[0099] The penalty formula is as follows: ; ; in, Represents a set of negative relative error values ​​in the root area, Represents the test point index set in the root area, Represents the weight coefficient of the negative relative error penalty in the root area, Indicates leaf The relative error of the test points is represents the exponential coefficient, Represents the negative relative error penalty mechanism in the root region.

[0100] Step 2.3b: Introduce an exponential penalty exceeding the threshold of 0.15 for the central region with the largest relative error in step 2.3.

[0101] The penalty formula is as follows: ; ; in, Represents the set of relative error values ​​in the middle area that exceed the threshold of 0.15. Indicates leaf The relative error of the test points is Represents the test point index set in the middle area, represents the weight coefficient of the relative error penalty exceeding the threshold of 0.15 in the middle area, represents the exponential coefficient, Relative error penalty mechanism for exceeding the threshold of 0.15 in the middle region.

[0102] Step 2.5: Calculate the load gradients of adjacent test points in all regions, and perform quadratic penalties on all gradients to smoothly control the load distribution gradients.

[0103] The formula is as follows: ; ; in, , represents the relative error gradient between adjacent test points, represents the total number of test points, Indicates leaf The relative error of the test points is represents the weight coefficient for smoothing the load distribution gradient in all regions, Represents the penalty mechanism for smoothing the gradient of the load distribution in all regions.

[0104] ; in, Represents the first objective function of the total relative error penalty mechanism for all regions, represents the negative relative error penalty mechanism in the root area, Relative error penalty mechanism for the central region exceeding the threshold of 0.15. represents the penalty mechanism for smooth control of the load distribution gradient in all regions, The fifth objective function representing the total relative error penalty mechanism for all regions.

[0105] Optional, combined Figure 4 ,like Figure 5 As shown, after step 2 and before step 3, step 2.6 is also included.

[0106] Step 2.6: Perform variance control on the loads in all regions, and update the fifth objective function to the sixth objective function.

[0107] The variance control penalty formula is: ; ; in, Indicates the area identification, which includes: root area root, root-middle transition area root-mid, middle area mid, middle-tip transition area mid-tip, tip area tip, Represents the index set of test points in the corresponding area, Indicates the number of elements in a collection. represents the mean of the relative error in the corresponding area, Indicates leaf The relative error of the test points is represents the variance of the relative error in the corresponding area, represents the weight coefficient of the variance in the corresponding region, represents the penalty mechanism for variance control within all regions.

[0108] The fifth objective function in step 2.4 is then updated to the sixth objective function.

[0109] ; in, Represents the first objective function of the total relative error penalty mechanism for all regions, represents the negative relative error penalty mechanism in the root area, Relative error penalty mechanism for the central region exceeding the threshold of 0.15. represents the penalty mechanism for smooth control of the load distribution gradient in all regions, The fifth objective function of the relative error penalty mechanism for all regions is represented by, The sixth objective function representing the total relative error penalty mechanism for all regions.

[0110] Step 3: Based on the objective function of step 2, a Bayesian global optimization algorithm is used to construct a Gaussian process proxy model of the objective function for optimization, and the global optimal solution of the optimization variables is output.

[0111] Among them, the global optimal solution of the optimization variables includes: the global optimal solution of the exciter position, the global optimal solution of the exciter mass and the global optimal solution of the exciting force amplitude.

[0112] It should be noted that the objective function in step 3 may be any one of the first objective function to the sixth objective function in step 2.

[0113] Optional, combined Figure 1 ,like Figure 6 As shown, step 3 includes the following steps 3.1 to 3.3.

[0114] Step 3.1: Set the Bayesian optimization parameters, enable parallel computing, and select the acquisition function you want to improve.

[0115] Among them, the Bayesian optimization parameters include: the maximum number of objective function evaluations is , the maximum running time is hours, the exploration-exploitation balance ratio is , the size of the Gaussian process activity set is The expected improvement as a collection function is expressed as: ; in, represents the acquisition function, represents the expected function, represents the currently observed optimal function value, Represents the optimization variable The corresponding objective function value may specifically be any one of the first objective function to the sixth objective function.

[0116] It should be noted that the expected improvement acquisition function helps select the point that is most likely to improve the current optimal result as the next evaluation point by quantifying the expected improvement of the objective function value at a candidate point compared to the current optimal value.

[0117] Step 3.2: Use the Gaussian process model as a proxy model for Bayesian optimization to solve the mean of the multivariate normal distribution and the covariance matrix of the multivariate normal distribution.

[0118] Among them, let the objective function is a realization of a Gaussian process, for any finite set of points , function value The mean is , the covariance matrix is The multivariate normal distribution of : ; in, represents the mean of the multivariate normal distribution, Represents the covariance matrix of the multivariate normal distribution.

[0119] Step 3.3: According to the objective function in step 2, initialize the Gaussian process model, execute the Bayesian optimization process, and output the global optimal solution of Bayesian optimization. , which is the global optimal solution Includes: Vibrator position The global optimal solution and the quality of the shaker The global optimal solution and exciting force amplitude of The global optimal solution of .

[0120] Among them, Statistics and Machine Learning Toolbox is called to perform Bayesian iterative optimization. Each iteration requires: updating the Gaussian process model, maximizing the acquisition function to determine the next evaluation point, evaluating the objective function, and updating the optimal solution.

[0121] Step 4: Taking the Bayesian optimization result of step 3 as the starting point, apply the interior point method to perform local precise optimization and output the final local optimization result of the optimization variable.

[0122] The final local optimization results of the optimization variables include: the local optimal solution of the exciter position, the local optimal solution of the exciter mass and the local optimal solution of the exciting force amplitude.

[0123] Optional, combined Figure 1 ,like Figure 7 As shown, the above step 4 includes steps 4.1 to 4.3.

[0124] Step 4.1: Set the interior point optimization parameters and enable parallel computing.

[0125] Among them, the optimization parameters of the interior point method include: the maximum number of function evaluations is , the maximum number of iterations is , the objective function gradient tolerance is , the variable change step tolerance is , and enable parallel computing.

[0126] Step 4.2: Global optimal solution obtained based on Bayesian optimization in step 3.3 , using the strategy of reducing the search space, the search range of the interior point method is limited to of( ) neighborhood to find the local optimal solution.

[0127] in, is the proportional coefficient, and its value range is , the lower and upper bound expressions of the optimization range are: ; in, and are the lower and upper bounds of the optimization range of the optimization variables of the Bayesian global optimization set in step 1, and Respectively represent the range-limiting parameters, and They respectively represent the lower and upper bounds of the optimization range of the optimization variable of the interior point method local optimization.

[0128] Step 4.3: Use Optimization Toolbox to perform interior point iterative optimization and output the final local optimization results of the optimization variables. .

[0129] That is, the final local optimization results of the optimization variables include: the local optimal solution of the exciter position, the local optimal solution of the exciter mass and the local optimal solution of the exciting force amplitude.

[0130] Among them, step 4.3 is implemented by steps 4.3a to 4.3c.

[0131] Step 4.3a: By introducing a logarithmic barrier function, the constrained optimization problem is transformed into a series of sub-problems with no constraints or only equality constraints, and the fatigue test parameters of the wind turbine blade are optimized.

[0132] ; ; in, Output variable when the comprehensive objective function takes the minimum value The value of represents the optimization variable vector, Corresponding to the position of the exciter , Vibrator quality and the excitation force amplitude ; is the comprehensive objective function defined, and are the lower and upper bounds of the optimization variable, respectively.

[0133] Step 4.3b: By reducing the obstacle parameter to approach the optimal solution, the interior point method transforms the original optimization problem into: ; in, Output variable when the objective function containing logarithmic obstacle term reaches the minimum value The value of represents the objective function including logarithmic obstacle term, is the comprehensive objective function defined, is the barrier parameter, , gradually decreases as the algorithm iterates. When the obstacle problem is solved Converge to the optimal solution of the original problem, Indicates optimization variables, and Respectively The lower and upper bounds of the optimization variables.

[0134] Step 4.3c: Solve step 4.3b and output the final local optimization result of the optimization variable. .

[0135] It can be understood that, on the one hand, the present application combines Bayesian global optimization with the local precise optimization of the interior point method to form a two-stage hybrid optimization strategy, which can not only effectively explore the global parameter space, but also perform high-precision local optimization in the local area, significantly improving the optimization efficiency and result accuracy. On the other hand, the present application innovatively proposes a regional segmented weighted objective function, sets differentiated weights, thresholds and penalty mechanisms for different functional areas of wind turbine blades, and especially strengthens the processing of negative relative errors in the root area and ultra-high relative errors in the middle area, effectively improving the problem that it is difficult to take into account the load matching accuracy of different areas of the blade in the traditional method, thereby further improving the test accuracy. In addition, the present application ensures the smooth transition of the test load between the various areas of the blade and the consistency within the area through the load distribution gradient smoothing control and regional variance constraint mechanism, avoids sudden changes in load distribution during the test process, reduces the risk of unreasonable stress distribution, and further improves the reliability of the test. On the other hand, the Gaussian process proxy model used in the present application can find the global optimal solution within a limited number of function evaluations, greatly reduces the computational overhead, accelerates the optimization process, and also improves the test efficiency. The present application is particularly suitable for processing computationally intensive wind turbine blade fatigue test optimization problems.

[0136] Example 2 In order to more clearly understand the technical solution of the present invention, an exemplary description is given below taking a 90-meter wind turbine blade as an example.

[0137] This embodiment provides a large wind turbine blade fatigue test load optimization method based on a hybrid optimization algorithm, which specifically includes: Step 1: According to the structural characteristics of the wind turbine blade, determine the test parameters and physical characteristic parameters, and define the optimization variables and the optimization variable range.

[0138] Step 1.1. In this embodiment, a 90-meter wind turbine blade is selected as the research object. When the blade is subjected to the fatigue test in the flapping direction, it is cut at 93% of the blade span, and the tip of the blade is cut off. The length of the wind turbine blade subjected to the full-size structural fatigue test is , 64 test positions were selected on the blade to optimize the matching between the test load and the target load, and the damping ratio was set to a fixed value of 0.0175.

[0139] Step 1.2, define the exciter position, exciter mass and excitation force amplitude as test parameters, and define the range of variation of the optimized exciter position, exciter mass and excitation force amplitude: Define the position of the exciter according to the blade length The range of variation is: [30m, 80m].

[0140] Define the vibration exciter mass according to the rope traction resonance loading method and loading position of the wind turbine blade fatigue test The range of variation is: [1000kg, 6000kg].

[0141] According to the wind turbine blade fatigue test rope traction resonance loading method, the exciter loading position, the quality of the exciter, the amplitude of the exciter is defined as The range of variation is: [8000N, 30000N].

[0142] Step 2: Establish a segmented weighted objective function for the leaf area, set differentiated weights and thresholds for each area, and construct an objective function for the differentiated penalty mechanism.

[0143] Step 2.1, the blade test area can be divided into: root area (0~20%), root-middle transition area (20~35%), middle area (35~55%), middle-tip transition area (55~70%), tip area (70~100%).

[0144] Step 2.2: Calculate the relative error between the test load and the target load in each area in step 2.1.

[0145] For any test point selected on the blade , and its relative error is defined as: ; in, Indicates any test point selected on the blade when the stress ratio is -1 The test load amplitude is Indicates the test point corresponding to the test load amplitude when the stress ratio is -1 given in the test report The target load amplitude.

[0146] Step 2.3: Based on the relative error of each region calculated in step 2.2, differentiated weights and thresholds are assigned to each region, and the relative error penalty formula for each region is constructed as follows: Root zone: ; Root-Mid Transition Zone: ; Central region: ; Mid-tip transition area: ; Tip area: ; in, Represents the index set of test points in the corresponding area, Represents the region identifier, which includes: root region root, root-middle transition region root-mid, middle region mid, middle-tip transition region mid-tip, tip region tip. The coefficient is the weight coefficient of the relative error penalty for the corresponding region exceeding the threshold. The threshold is the tolerance threshold of the corresponding region. Indicates leaf The relative error of each test point. The contexts are the same region identifier and have a one-to-one correspondence.

[0147] Step 2.3a: Special treatment is given to the negative relative error of the root in step 2.3. Only the negative relative error value of the root area is retained and penalized by an exponential function: ; ; in, Represents a set of negative relative error values ​​in the root area, Represents the test point index set in the root area, Indicates leaf The relative error of each test point is the weight coefficient of the negative relative error penalty in the root area, and 100 is the exponential coefficient.

[0148] Step 2.3b: Introduce an exponential penalty exceeding the threshold of 0.15 for the central region with the largest relative error in step 2.3: ; ; in, Represents the set of relative error values ​​in the middle area that exceed the threshold of 0.15. Indicates leaf The relative error of the test points is represents the test point index set in the middle area, the coefficient 50 is the weight coefficient of the relative error penalty exceeding the threshold value 0.15 in the middle area, and 10 is the exponential coefficient.

[0149] Step 2.5: Calculate the load gradients of adjacent test points in all regions, and perform quadratic penalties on all gradients to smoothly control the load distribution gradients.

[0150] Calculate the load gradients of adjacent test points in all regions and perform a quadratic penalty on all gradients as follows: ; ; in, ;letter Represents the set of relative error values ​​in the middle area that exceed the threshold of 0.15. Indicates leaf The relative error of the test points is represents the test point index set in the middle area, the coefficient is the weight coefficient of the relative error penalty exceeding the threshold of 0.15 in the middle area, and 10 is the exponential coefficient.

[0151] Step 2.6: Perform variance control on the loadings in all regions.

[0152] After the load is subjected to the super-relative error penalty constraint and the load distribution gradient is smoothly controlled, in order to ensure the uniformity and coordination of the load distribution, the variance of the load in all regions needs to be controlled. The variance control penalty formula is: ; ; in, Indicates the area identification (such as root area root, root-middle transition area root-mid, middle area mid, middle-tip transition area mid-tip, tip area tip), Represents the index set of test points in the corresponding area, Indicates the number of elements in a collection. represents the mean of the relative error in the corresponding area, Indicates leaf The relative error of the test points is Represents the variance of the relative error in the corresponding area, and the coefficient is the weight coefficient of the variance in the corresponding area. The contexts are the same region identifier and have a one-to-one correspondence.

[0153] In addition, in actual working conditions, the following steps 2.7 and 2.8 may also be included.

[0154] Step 2.7: Constrain the position of the vibrator to ensure that the position of the vibrator is within the range of [30m, 80m] defined in step 1. The constraint formula is: ; in, represents the shaker position constraint function, represents the indicator function, Indicates the position of the exciter in the span direction of the blade.

[0155] Step 2.8, combining steps 2.2 to 2.7, finally obtaining the objective function (i.e., the updated objective function 6 of the embodiment, where step 2.7 introduces This is to ensure that the position of the vibrator is within the range of [30m, 80m] defined in step 1, and to avoid optimizing outside this range, which can be ignored) is: ; It should be noted that Example 2 is implemented by splitting step 2.6 in Example 1 into steps 2.6 and 2.8 in Example 2. Step 2.7 is added here to ensure that the position of the exciter is within the range of [30m, 80m] defined in step 1. In actual working conditions, step 2.7 may be added or not. If it is not added, if an error occurs in the solution, return to step 1 and re-optimize within the defined interval.

[0156] Step 3: Based on the objective function of step 2, a Bayesian global optimization algorithm is used to construct a Gaussian process proxy model of the objective function for optimization, and the global optimal solution of the optimization variables is output.

[0157] It should be noted that the objective function of this embodiment 2 adopts the objective function 6 in the above step 2.8.

[0158] Specifically, the Bayesian global optimization process is as follows: Step 3.1. Set the Bayesian optimization parameters: the maximum number of objective function evaluations is 300, the maximum running time is 20 hours, the exploration-exploitation balance ratio is 0.5, the Gaussian process active set size is 350, parallel computing is enabled, and the expected improvement (EI) acquisition function is selected.

[0159] The expected improvement acquisition function is used to balance exploration and exploitation and calculate the expected improvement of the potential evaluation point relative to the current optimal value. It is defined as: ; in, represents the currently observed optimal function value, Represents the optimization variable The corresponding objective function value.

[0160] Step 3.2: Use the Gaussian process as a proxy model for Bayesian optimization to solve the mean of the multivariate normal distribution and the covariance matrix of the multivariate normal distribution.

[0161] Specifically, its theoretical basis is to set the objective function is a realization of a Gaussian process, that is, the function values ​​at any finite number of points jointly obey the multivariate normal distribution. Formally, for any finite set of points , function value The mean is , the covariance matrix is The multivariate normal distribution of : ; in, represents the mean of the multivariate normal distribution, Represents the covariance matrix of the multivariate normal distribution.

[0162] Step 3.3: According to the objective function in step 2, initialize the Gaussian process model, execute the Bayesian optimization process, and output the global optimal solution of Bayesian optimization. .

[0163] Specifically, the Bayesian optimization process is as follows: First, construct the objective function evaluation interface; second, initialize the Gaussian process model; then, call Statistics and Machine Learning Toolbox to perform Bayesian iterative optimization, where each iteration requires: updating the Gaussian process model, maximizing the acquisition function to determine the next evaluation point, evaluating the objective function, and updating the optimal solution; finally, output the optimal solution of Bayesian optimization .

[0164] Step 4: Taking the Bayesian optimization result of step 3 as the starting point, apply the interior point method to perform local precise optimization and output the final local optimization result of the optimization variable.

[0165] It should be noted that although Bayesian optimization can effectively explore the global parameter space, its sampling nature and the approximate characteristics of the proxy model limit the optimization accuracy. In order to further improve the optimization accuracy of wind turbine blade fatigue test parameters, the interior point method is introduced for local precise optimization. As an efficient algorithm for solving nonlinear constrained optimization problems, the interior point method has the characteristics of fast convergence speed and good numerical stability, and is suitable for precise search in the local area provided by Bayesian optimization.

[0166] Specifically, the implementation of step 4 can be achieved through the following steps 4.1 to 4.5.

[0167] Step 4.1: Set the interior point optimization parameters and enable parallel computing.

[0168] Specifically, set the interior point optimization parameters: the maximum number of function evaluations is 350, the maximum number of iterations is 200, and the objective function gradient tolerance is , the variable change step tolerance is , and enable parallel computing.

[0169] Step 4.2: Global optimal solution based on Bayesian optimization , using the strategy of reducing the search space, the search range of the interior point method is limited to of( ) neighborhood to find the local optimal solution.

[0170] Specifically, the lower and upper bound expressions of the optimization range are: ; in, and are the lower and upper bounds of the optimization range of the optimization variables of the Bayesian global optimization set in step 1, 0.95 and 1.05 represent the range-limiting parameters, and They respectively represent the lower and upper bounds of the optimization range of the optimization variable of the interior point method local optimization.

[0171] Step 4.3: Perform iterative optimization using the interior point method and output the final local optimization results of the optimization variables. .

[0172] Specifically, use Optimization Toolbox to perform interior point iterative optimization. The optimization steps are as follows: First, the constrained optimization problem is transformed into a series of sub-problems with no constraints or only equality constraints by introducing a logarithmic barrier function. For the optimization of wind turbine blade fatigue test parameters, the standard form of constrained optimization problem is considered: ; ; in, Output variable when the comprehensive objective function takes the minimum value The value of represents the optimization variable vector, Corresponding to the position of the exciter , Vibrator quality and the excitation force amplitude ; is the defined comprehensive objective function; and are the lower and upper bounds of the optimization variable, respectively.

[0173] Then, by reducing the obstacle parameter to approach the optimal solution, the interior point method transforms the original problem into: ; in, Output variable when the objective function containing logarithmic obstacle term reaches the minimum value The value of represents the objective function including logarithmic obstacle term, is the comprehensive objective function defined, is the barrier parameter, which decreases gradually as the algorithm iterates. When the obstacle problem is solved Converge to the optimal solution of the original problem, Indicates optimization variables, and Respectively The lower and upper bounds of the optimization variables.

[0174] Finally, output the final optimization results .

[0175] It should be noted that, usually, after the above steps 1 to 4 are completed, the optimization method of the present application is completed. However, in order to further illustrate the superiority of the present application, the following step 5 is added to verify the result and compare the optimization result. The details are as follows:

[0176] Step 5: Conduct experimental verification on the optimization results.

[0177] According to the parameter settings of steps 1 to 4, and combined with the above optimization process, the final optimized parameters are: exciter position 47m, exciter mass 6000kg, excitation force amplitude 24043N. The final optimization results are as follows: Figure 8 As shown (where Figure 8 (a) is a comparison diagram between the optimized bending moment amplitude and the target bending moment amplitude; Figure 8 (b) is the relative error diagram between the optimized bending moment amplitude and the target bending moment amplitude). Figure 8 From (a), it can be seen that the optimized bending moment distribution (solid line part) and the target bending moment distribution (dot part) show a good matching relationship. The optimized bending moment amplitude is always slightly higher than the target bending moment amplitude, which meets the safety requirements of the test standard and verifies the effectiveness of the regional segmented weighting strategy. The optimized bending moment amplitude distribution curve is smooth and has no mutation points, which confirms the effectiveness of gradient control. Figure 8As can be seen from (b), the relative errors of the 64 test points selected along the span of the blade are all controlled within 12%. The relative error results are better than the existing technical level, and the relative error consistency in each region is good, verifying the effect of regional variance control. In addition, although the maximum number of objective function evaluations for Bayesian global optimization is set to 300, Bayesian optimization can find the global optimal solution of the final output in the first 160 evaluations, while the interior point method can achieve high-precision convergence in only 20 iterations. This shows that the use of the Gaussian process surrogate model effectively reduces the consumption of computing resources and accelerates the optimization convergence process, while the use of the interior point method can achieve high-precision local optimization with fewer iterations. The use of Bayesian global optimization and interior point method local optimization at the same time avoids the problem of falling into the local optimum, further reflecting the superiority of using a hybrid optimization algorithm to optimize the load distribution of full-scale structural fatigue tests of large wind turbine blades.

[0178] It can be understood that, on the one hand, the present application combines Bayesian global optimization with the local precise optimization of the interior point method to form a two-stage hybrid optimization strategy, which can not only effectively explore the global parameter space, but also perform high-precision local optimization in the local area, significantly improving the optimization efficiency and result accuracy. On the other hand, the present application innovatively proposes a regional segmented weighted multi-objective function, sets differentiated weights, thresholds and penalty mechanisms for different functional areas of wind turbine blades, and especially strengthens the processing of negative relative errors in the root area and ultra-high relative errors in the middle area, effectively improving the problem that it is difficult to take into account the load accuracy of different areas of the blade in the traditional method, thereby further improving the test accuracy. In addition, the present application ensures the smooth transition of the test load between the various areas of the blade and the consistency within the area through the load distribution gradient smoothing control and regional variance constraint mechanism, avoids sudden changes in load distribution during the test process, reduces the risk of unreasonable stress distribution, and further improves the reliability of the test. On the other hand, the Gaussian process proxy model used in the present application can find the global optimal solution within a limited number of function evaluations, greatly reduces the computational overhead, accelerates the optimization process, and also improves the test efficiency. The present application is particularly suitable for processing computationally intensive wind turbine blade fatigue test optimization problems.

[0179] This application uses specific examples to illustrate the principles and implementation methods of the present invention. The description of the above embodiments is only used to help understand the core idea of ​​the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A large wind turbine blade fatigue test load optimization method based on a hybrid optimization algorithm, characterized in that: The method comprises: Step 1: According to the structural characteristics of the wind turbine blade, determine the test parameters and physical characteristic parameters, and define the optimization variables and the optimization variable range; Among them, the test parameters include: wind turbine blade length, test position point distribution and target bending moment distribution, the physical characteristic parameters include: blade damping ratio, and the optimization variables include: exciter position, exciter mass and excitation force amplitude; Step 2: Establish a segmented weighted objective function for the leaf area, set the differentiated weight and threshold for each area, and construct the objective function of the differentiated penalty mechanism; The blade is divided into five regions: root region, root-middle transition region, middle region, middle-tip transition region, and tip region; the objective function is the accumulation of penalty mechanisms for each region; Step 3: Based on the objective function of step 2, a Bayesian global optimization algorithm is used to construct a Gaussian process proxy model of the objective function for optimization, and the global optimal solution of the optimization variables is output; Among them, the global optimal solution includes: the global optimal solution of the exciter position, the global optimal solution of the exciter mass and the global optimal solution of the exciting force amplitude; Step 4: Using the Bayesian optimization result of step 3 as the starting point, apply the interior point method to perform local precise optimization and output the final local optimization result of the optimization variable; The final local optimization results of the optimization variables include: the local optimal solution of the exciter position, the local optimal solution of the exciter mass and the local optimal solution of the exciting force amplitude.

2. The large wind turbine blade fatigue test load optimization method based on hybrid optimization algorithm according to claim 1 is characterized in that: Step 1 includes step 1.1 and step 1.2; Step 1.1, according to the structural characteristics of the wind turbine blade, set the test parameters and physical characteristic parameters; Step 1.2, define the exciter position, exciter mass and excitation force amplitude as optimization variables, and define the variation range of the optimized exciter position, exciter mass and excitation force amplitude; Among them, the range of the exciter position change is defined according to the blade length: According to the loading method and loading position, the range of the vibration exciter mass change is defined as According to the loading method, the loading position of the exciter, and the mass of the exciter, the amplitude variation range of the exciter is defined as .

3. The large wind turbine blade fatigue test load optimization method based on hybrid optimization algorithm according to claim 1 is characterized in that: Step 2 includes: Step 2.1, dividing the blade test area into: root area, root-middle transition area, middle area, middle-tip transition area, tip area; Among them, the proportion of the root area to the entire leaf position is: , the root-middle transition area accounts for the proportion of the entire leaf position: , the proportion of the middle area to the entire blade position is: , the proportion of the middle-tip transition area to the entire blade position is: , the tip area accounts for the proportion of the entire blade position: ;in, ; Step 2.2, calculate the relative error between the test load and the target load in each area in step 2.1; Among them, for any test point selected on the blade , and its relative error is defined as: ； in, Indicates leaf The relative error of the test points is Indicates any test point selected on the blade when the stress ratio is -1 The test load amplitude is Indicates the test point corresponding to the test load amplitude when the stress ratio is -1 given in the test report The target load amplitude; Step 2.3: Based on the relative error of each region calculated in step 2.2, assign differentiation weights and thresholds to each region, and construct the error penalty formula for each region as follows: ； in, Indicates the area identifier, which includes: root area root, root-middle transition area root-mid, middle area mid, middle-tip transition area mid-tip, tip area tip; Represents the index set of test points in the corresponding area, Indicates leaf The relative error of the test points is represents the weight coefficient of the error penalty for the corresponding area exceeding the threshold, represents the tolerance threshold of the corresponding area, The mechanism that represents the corresponding region penalizes errors exceeding the threshold; Step 2.4: After constructing the error penalty mechanism for each region in step 2.3, the first objective function of the total error penalty mechanism for all regions of the blade is obtained, which is expressed as: ； in, represents the first objective function of the total error penalty mechanism for all regions, represents the error penalty mechanism for root area exceeding the threshold, represents the error penalty mechanism for exceeding the threshold in the root-middle transition region, represents the error penalty mechanism for the central area exceeding the threshold, represents the error penalty mechanism for the mid-to-apex transition region exceeding the threshold, Indicates the error penalty mechanism for the peak area exceeding the threshold.

4. The large wind turbine blade fatigue test load optimization method based on hybrid optimization algorithm according to claim 3 is characterized in that: In step 2.3, further comprising step 2.3a; Step 2.3a: Special treatment is given to the root negative error in step 2.

3. Only the negative error value in the root area is retained and penalized by an exponential function: ； ； in, Represents a set of negative error values ​​in the root area, Represents the test point index set in the root area, represents the weight coefficient of the negative error penalty in the root area, Indicates leaf The relative error of the test points is represents the exponential coefficient, Represents the negative error penalty mechanism in the root area; After executing step 2.3a, use the total error penalty mechanism in the root area Update the error penalty mechanism for root area exceeding the threshold in step 2.4 ; Among them, the total error penalty mechanism for the root area is expressed as: ; Similarly, the first objective function in step 2.4 is updated to the second objective function; ； in, represents the total error penalty mechanism in the root region, represents the error penalty mechanism for root area exceeding the threshold, Represents the negative error penalty mechanism in the root area; Represents the first objective function of the total error penalty mechanism for all regions; The second objective function representing the total error penalty mechanism for all regions.

5. The large wind turbine blade fatigue test load optimization method based on hybrid optimization algorithm according to claim 3 is characterized in that: In step 2.3, further comprising step 2.3b; Step 2.3b: introduce an exponential penalty exceeding the threshold of 0.15 for the central region with the largest error in step 2.3; ； ； in, represents the set of error values ​​in the middle area that exceed the threshold of 0.

15. Indicates leaf The relative error of the test points is Represents the test point index set in the middle area, represents the weight coefficient of the error penalty exceeding the threshold of 0.15 in the middle area, represents the exponential coefficient, It represents the error penalty mechanism for the error exceeding the threshold of 0.15 in the middle area; After executing step 2.3b, use the total error penalty mechanism in the middle area Update the error penalty mechanism for the central area exceeding the threshold in step 2.4 above ; Among them, the total error penalty mechanism for the central area is expressed as: ; Similarly, the first objective function in step 2.4 is updated to the third objective function; ； in, represents the total error penalty mechanism in the middle region, represents the error penalty mechanism for the central area exceeding the threshold, represents the error penalty mechanism for the error exceeding the threshold of 0.15 in the middle area, represents the first objective function of the total error penalty mechanism for all regions, The third objective function representing the total error penalty mechanism for all regions.

6. The large wind turbine blade fatigue test load optimization method based on hybrid optimization algorithm according to claim 3 is characterized in that: After step 2.3 and before step 2.4, step 2.5 is further included; Step 2.5, calculate the load gradients of adjacent test points in all regions, and perform quadratic penalties on all gradients to smoothly control the load distribution gradients; The formula is as follows: ； ； Similarly, the first objective function in step 2.4 is updated to the fourth objective function; ； in, , represents the error gradient between adjacent test points, represents the total number of test points, Indicates leaf The relative error of the test points is represents the weight coefficient for smooth control of the load distribution gradient in all regions, represents the penalty mechanism for smooth control of the load distribution gradient in all regions; Represents the first objective function of the total error penalty mechanism for all regions; The fourth objective function representing the total error penalty mechanism for all regions.

7. The large wind turbine blade fatigue test load optimization method based on hybrid optimization algorithm according to claim 3 is characterized in that: After step 2.3, step 2.3a, step 2.3b, step 2.5 and step 2.4 are sequentially performed, and the first objective function in step 2.4 is updated to the fifth objective function; Step 2.3a: Specially handle the root negative error in step 2.3, retain only the negative error value in the root area, and penalize it through an exponential function; ； ； in, Represents a set of negative error values ​​in the root area, Represents the test point index set in the root area, represents the weight coefficient of the negative error penalty in the root area, Indicates leaf The relative error of the test points is represents the exponential coefficient, Represents the negative error penalty mechanism in the root area; Step 2.3b: introduce an exponential penalty exceeding the threshold of 0.15 for the central region with the largest error in step 2.3; ； ； in, represents the set of error values ​​in the middle area that exceed the threshold of 0.

15. Indicates leaf The relative error of the test points is Represents the test point index set in the middle area, represents the weight coefficient of the error penalty exceeding the threshold of 0.15 in the middle area, represents the exponential coefficient, It represents the error penalty mechanism for the error exceeding the threshold of 0.15 in the middle area; Step 2.5, calculate the load gradients of adjacent test points in all regions, and perform quadratic penalties on all gradients to smoothly control the load distribution gradients; ； ； in, , represents the error gradient between adjacent test points, represents the total number of test points, Indicates leaf The relative error of the test points is represents the weight coefficient for smooth control of the load distribution gradient in all regions, represents the penalty mechanism for smooth control of the load distribution gradient in all regions; ； in, represents the first objective function of the total error penalty mechanism for all regions, represents the negative error penalty mechanism in the root area, represents the error penalty mechanism for the error exceeding the threshold of 0.15 in the middle area, represents the penalty mechanism for smooth control of the load distribution gradient in all regions, The fifth objective function representing the total error penalty mechanism for all regions.

8. The large wind turbine blade fatigue test load optimization method based on hybrid optimization algorithm according to claim 7 is characterized in that: After step 2 and before step 3, step 2.6 is also included; Step 2.6, perform variance control on the loads in all regions, and update the fifth objective function to the sixth objective function; Among them, the variance control penalty formula is: ； ； in, Indicates the area identification, which includes: root area root, root-middle transition area root-mid, middle area mid, middle-tip transition area mid-tip, tip area tip, Represents the index set of test points in the corresponding area, Indicates the number of elements in a collection. represents the mean value of the error in the corresponding area, Indicates leaf The relative error of the test points is represents the variance of the error in the corresponding region, represents the weight coefficient of the variance in the corresponding region, represents the penalty mechanism for variance control in all regions; Similarly, the fifth objective function in step 2.4 is updated to the sixth objective function; ； in, represents the first objective function of the total error penalty mechanism for all regions, represents the negative error penalty mechanism in the root area, represents the error penalty mechanism for the error exceeding the threshold of 0.15 in the middle area, represents the penalty mechanism for smooth control of the load distribution gradient in all regions, The fifth objective function of the total error penalty mechanism for all regions is represented by, The sixth objective function representing the total error penalty mechanism for all regions.

9. The large wind turbine blade fatigue test load optimization method based on hybrid optimization algorithm according to claim 1 is characterized in that: Step 3 includes: Step 3.1, set the Bayesian optimization parameters, enable parallel computing, and select the acquisition function you want to improve; Among them, the Bayesian optimization parameters include: the maximum number of objective function evaluations is , the maximum running time is hours, the exploration-exploitation balance ratio is , the size of the Gaussian process activity set is ; The expected improvement as a collection function is expressed as: ； in, represents the acquisition function, represents the expected function, represents the currently observed optimal function value, Represents the optimization variable The corresponding objective function value; Step 3.2, using the Gaussian process model as a proxy model for Bayesian optimization to solve the mean of the multivariate normal distribution and the covariance matrix of the multivariate normal distribution; Among them, let the objective function is a realization of a Gaussian process, for any finite set of points , function value The mean is , the covariance matrix is The multivariate normal distribution of : ; in, represents the mean of the multivariate normal distribution, Represents the covariance matrix of the multivariate normal distribution; Step 3.3: According to the objective function in step 2, initialize the Gaussian process model, execute the Bayesian optimization process, and output the global optimal solution of Bayesian optimization. ; Among them, the global optimal solution Includes: Vibrator position The global optimal solution and the quality of the shaker The global optimal solution and exciting force amplitude of The global optimal solution of .

10. The large wind turbine blade fatigue test load optimization method based on hybrid optimization algorithm according to claim 9, characterized in that: Step 4 includes; Step 4.1, set the interior point optimization parameters and enable parallel computing; Among them, the optimization parameters of the interior point method include: the maximum number of function evaluations is , the maximum number of iterations is , the objective function gradient tolerance is , the variable change step tolerance is , and enable parallel computing; Step 4.2: Global optimal solution obtained based on Bayesian optimization in step 3.3 , using the strategy of reducing the search space, the search range of the interior point method is limited to of( ) Find the local optimal solution in the neighborhood; in, is the proportional coefficient, and its value range is , the lower and upper bound expressions of the optimization range are: ； in, and are the lower and upper bounds of the optimization range of the optimization variables of the Bayesian global optimization set in step 1, and Respectively represent the range-limiting parameters, and They represent the lower and upper bounds of the optimization range of the optimization variable of the interior point method local optimization respectively; Step 4.3: Perform iterative optimization using the interior point method and output the final local optimization results of the optimization variables. ; The final local optimization results of the optimization variables include: the local optimal solution of the exciter position, the local optimal solution of the exciter mass and the local optimal solution of the excitation force amplitude; Wherein, step 4.3 is implemented by steps 4.3a to 4.3c; Step 4.3a, by introducing a logarithmic barrier function, the constrained optimization problem is transformed into a series of sub-problems with no constraints or only equality constraints, and the fatigue test parameters of the wind turbine blade are optimized; ； ； in, Output variable when the comprehensive objective function takes the minimum value The value of represents the optimization variable vector, Corresponding to the position of the exciter , Vibrator quality and the excitation force amplitude ; is the defined comprehensive objective function; and are the lower and upper bounds of the optimization variable, respectively; Step 4.3b: By reducing the obstacle parameter to approach the optimal solution, the interior point method transforms the original optimization problem into: ； in, Output variable when the objective function containing logarithmic obstacle term reaches the minimum value The value of represents the objective function including logarithmic obstacle term, is the comprehensive objective function defined, is the barrier parameter, , gradually decreases as the algorithm iterates. When the obstacle problem is solved Converge to the optimal solution of the original problem, Indicates optimization variables, and Respectively lower and upper bounds of the optimization variables; Step 4.3c: Solve step 4.3b and output the final local optimization result of the optimization variable .

Images