Helicopter vibration closed-loop control method based on cone constraint

Abstract

The invention provides a helicopter vibration closed-loop control method based on cone constraint, and belongs to the field of vibration control. The method comprises the following steps: firstly, establishing an aeroelastic coupling model of a rotor according to structural characteristics of the helicopter rotor, setting dissimilar physical quantities among blades, and calculating a blade cone mutual difference value and a hub 1 / rev vibration harmonic wave; secondly, establishing a quadratic performance index function of a control problem by taking a blade total pitch variation and a rotor wing 1 / rev vibration harmonic component as control input and output; thirdly, according to the relation between the cone and the blade total pitch, the limitation of the total pitch variation is considered at the same time, and a constraint equation of control input is established; and 4, estimating a transfer matrix between control input and output to obtain an optimal control problem, and converting the optimal control problem into a quadratic programming problem. And finally, solving to obtain the optimal control quantity. The method is not only suitable for numerical simulation verification, but also can be applied to a real-time control test, improves the identification efficiency, has very high operability and feasibility, and is convenient for practical application.

Description

Technical Field

[0001] The invention belongs to the field of vibration control and relates to a helicopter vibration closed-loop control method based on cone constraints. Background Art

[0002] When the helicopter is in flight, the rotor is in a high-speed rotating, highly unstable aerodynamic environment. At the same time, there is also an aeroelastic coupling between the flexible blades. The rotor will produce violent vibrations and noise, which converge and overlap at the center of the hub and are transmitted to the helicopter body and cabin through the hub shaft. When the rotor blades are exactly the same, the hub is equivalent to a filter, and only high-order harmonics of integer multiples of the rotation frequency, that is, kN / rev frequency (k is a positive integer, N is the number of blades), can be transmitted to the fuselage. At this time, the high-order harmonics are the main vibrations; however, when the blades are not exactly the same, all harmonics including non-kN / rev frequencies will be transmitted to the fuselage. At this time, the 1 / rev vibration harmonic, that is, the harmonic equal to the rotor rotation frequency, accounts for a larger proportion, and the cone formed by the rotation trajectory of each blade no longer overlaps. This phenomenon is called rotor cone deviation, which will significantly affect the sensitivity of airborne equipment and cause crew fatigue. Reducing 1 / rev vibration and reducing rotor cone deviation is a long-term and necessary task. Due to errors in the manufacturing and installation processes, or unavoidable problems such as blade wear and deformation during later use, each blade is not always exactly the same. The traditional method is generally to manually adjust the blade collective pitch angle, counterweight or trailing edge flap slightly, and this process requires adjustment and test flight, and repeated execution to ensure that the vibration is reduced to the ideal range. This process is time-consuming, labor-intensive and costly. Reducing the rotor cone deviation is also an issue that deserves attention. Eric Bechhoefer and others believe that it is related to the vertical load and may also excite high-order harmonics. Reducing the cone deviation will also optimize the 2 / rev component of the hub vibration load.

[0003] In response to the above-mentioned helicopter vibration problem, the Chinese invention patent [CN 106081078 A] discloses a helicopter rotor blade vibration reduction device. The invention uses a special vibration absorber device, which is installed inside the blade near the blade tip, and absorbs the blade flapping vibration through the movement of the mass block, thereby achieving the purpose of reducing the vibration of the helicopter. However, the invention belongs to a passive vibration reduction device. The vibration absorber is installed in the blade, which limits the volume and weight of the mass block therein. If it is too large, it will affect the normal operation of the blade, and if it is too small, it will not achieve the purpose of vibration reduction. Therefore, it is difficult for this vibration absorber device to further improve the vibration reduction performance. Therefore, as people's demand for helicopter comfort becomes higher and higher, passive vibration reduction technology has been difficult to achieve the ideal vibration level.

[0004] Active vibration control technology has been widely studied due to its high control potential. Among them, the Chinese invention patent [CN 116424552 B] proposes an active torsion blade vibration control method and system. The invention changes the aerodynamic layout of the blade by actively twisting the blade under the action of electromechanical coupling to improve the aerodynamic load and reduce the hub vibration. It can also calculate the optimal torsion change according to the vibration level of the helicopter. However, the driver installed in this method is integrated into the blade, that is, it is necessary to embed active fiber composite materials or coarse fiber piezoelectric materials in the blade skin. This not only requires the blade to be remade, but the embedded driver will also affect the material properties of the blade, and even requires the blade structure to be redesigned, which will bring a lot of additional workload. In addition, this control method does not take into account the rotor cone deviation, and the controller performance is not fully utilized.

[0005] In view of the above difficulties, it is necessary to develop an active vibration control algorithm to reduce helicopter vibration and rotor cone deviation. Therefore, the present invention proposes an active control method that can give full play to the vibration control effect, and the method can take into account the rotor cone deviation value, and does not make any changes to the blades, only slightly adjusting the blade pitch angle, so as to maximize the performance of the active vibration control algorithm. Summary of the invention

[0006] In order to solve the above technical problems, the present invention proposes a closed-loop control method for helicopter vibration based on cone constraints. According to the independent blade control (IBC) concept, the total pitch angle change of the blade is used as the control input, the control target is the 1 / rev harmonic component of the vibration at the center of the rotor hub, and the control constraint is the rotor cone deviation value. The implementation of the control method is divided into two scenarios. The first is to verify the feasibility and effectiveness of the control method at the numerical simulation level, and the second is to apply the fully verified control method to the real test. Therefore, the present invention needs to first establish a rotor aeroelastic coupling model with blade dissimilarity during numerical simulation, so as to calculate the rotor cone deviation value and the 1 / rev harmonic component of the hub center vibration velocity in the control method, and directly use the cone sensor and the vibration sensor to measure the two data of the test rotor system in the test. In addition, the vibration measured by the vibration sensor in the test is measured in terms of vibration velocity. Therefore, in order to maintain the unity of the control target, it is also necessary to convert it into vibration velocity in the numerical simulation stage. Then, after obtaining the control target data, a quadratic performance index function of the vibration problem is established. The total pitch angle change of the rotor blades is taken as the control input, the vibration velocity 1 / rev harmonic component of the rotor hub center is taken as the control target, and the rotor cone deviation value is taken as the control constraint. A constrained optimization problem is constructed, and the transfer matrix is ​​estimated using the least squares identification algorithm. The quadratic programming optimization algorithm is used to solve and obtain the optimal control quantity.

[0007] In order to achieve the above object, the technical solution adopted by the present invention is:

[0008] A helicopter vibration closed-loop control method based on cone constraints includes the following steps:

[0009] Step 1: Set the dissimilar physical quantities between the blades, establish the rotor aeroelastic coupling model, and calculate the harmonic component of the rotor hub center vibration velocity 1 / rev and the rotor cone deviation value. This step is applied to the numerical simulation scenario of the control method. When applied to the test, the cone sensor and the vibration sensor are directly used to measure the two data to replace the numerical calculation of this step. Specifically:

[0010] First, set the dissimilar physical quantities between the blades. Wear, deformation, aging and other loss problems that occur during the use of the blades will affect the material properties and aerodynamic characteristics of the blades themselves. Therefore, these problems are reflected in specific physical quantities, namely, the flapping stiffness. , shimmy stiffness , torsional stiffness , linear density dm, lift coefficient , drag coefficient , pitch moment coefficient These physical quantities will change to varying degrees after the blades are worn, which will affect the dynamic response of the rotor and cause an increase in the harmonic component of the vibration velocity 1 / rev at the center of the hub.

[0011] Then, based on the multi-body dynamics theory and considering the dissimilar physical quantities of the blades, the rotor aeroelastic coupling model is established and the rotor cone deviation value is solved. , we obtain a set of differential algebraic equations (DAEs) to extract the generalized acceleration and generalized coordinates The linear term coefficient of is obtained, and the overall aeroelastic coupling model of the rotor is obtained, which is organized into the standard differential algebraic equation form as follows:

[0012] (1)

[0013] In the formula, They are the mass matrix, stiffness matrix, generalized nonlinear term and rotor aerodynamic term of the rotor aeroelastic coupling model; Indicates the rotor's total pitch angle and its variation, which directly affects the aerodynamic terms ,in It is the control input; are generalized coordinates and generalized acceleration respectively; the superscript T represents the matrix transpose symbol; Represents the constraint equations of the rotor aeroelastic coupling model; represents the transpose of the Jacobian matrix of the constraint equation, that is ; is the Lagrange multiplier.

[0014] The generalized-α method is used to solve the rotor aeroelastic coupling model (1), and the dynamic response of each blade can be obtained: , and then the tip height of each blade can be calculated , assuming the number of rotor blades is , arbitrarily select the e-th blade as the reference blade, and calculate the rotor cone deviation value, that is, the difference between the tip height of other blades and the reference blade, as follows:

[0015] (2)

[0016] In the formula, Represents the rotor cone deviation value vector; Indicates the 1st, 2nd, ..., The tip height of the propeller blade; Indicates the tip height value of the e-th blade.

[0017] Furthermore, the vibration load at the hub center is solved. The vibration load at the hub center of the rotor is the main source of helicopter vibration. It is the vibration load generated by each blade during high-speed rotation, which is transmitted from the blade root to the hub center and superimposed. The hub center vibration load is calculated from the superposition of the blade root load, specifically:

[0018] (3)

[0019] In the formula, is the azimuth angle of the mth blade; is the number of blades; the superscript H indicates the hub; the superscript m indicates the mth blade; the subscripts x, y, and z indicate the force or moment in the direction of the x, y, and z coordinate axes; They represent the root shear force of the mth blade in the x, y, and z directions respectively; They represent the root moments of the mth blade in the x, y, and z directions respectively; Respectively represent the vibration shear force of the hub center in the x, y, and z directions; They represent the vibration moments of the hub center in the x, y, and z directions respectively.

[0020] Finally, the hub center vibration load is converted to the hub center vibration velocity, and the 1 / rev harmonic component of the vibration velocity is extracted using Fourier transform. In engineering, vibration is generally measured using a vibration sensor, and the vibration velocity in the three directions of x, y, and z at the hub center is used to characterize the magnitude of the vibration. Therefore, combined with the relationship between vibration load and vibration velocity, the hub vibration load is converted to vibration velocity as follows:

[0021] (4)

[0022] In the formula, They represent the vibration speeds in the three directions of x, y, and z at the hub center respectively; the superscript T represents the transposition symbol; is a 3×6 dimensional transformation matrix.

[0023] The calculated vibration velocity Use Fourier transform to do frequency domain analysis and extract the 1 / rev harmonic component, that is: The 1 / rev harmonic components include sinusoidal components and cosine components ; The 1 / rev harmonic components include sinusoidal components and cosine components ; The 1 / rev harmonic components include sinusoidal components and cosine components The control target is obtained as follows:

[0024] (5)

[0025] In the formula, It represents the 1 / rev harmonic component of the vibration velocity, which is also the control target. The superscript T represents the transposition sign.

[0026] Step 2: Taking the blade total pitch angle change as the control input and the rotor vibration speed 1 / rev harmonic component as the control target, establish the quadratic performance index function of the control problem. Specifically:

[0027] Based on the optimal control theory, the performance index function of the control problem is established. It is generally constructed as a quadratic function, which has the characteristics of clear physical meaning and easy mathematical processing. Therefore, the performance index function is set as:

[0028] (6)

[0029] Where: Represents the performance indicator function value; is the control target, i.e., the harmonic component of the rotor vibration speed 1 / rev in equation (5); the subscript n represents the nth closed-loop control cycle; is the control input of the nth closed-loop control cycle, i.e., the total pitch angle change in equation (1); is the change in the control input from the n-1th closed-loop control cycle to the nth closed-loop control cycle; The control targets are , control input , Control input change The corresponding weight matrix.

[0030] In addition, it is necessary to establish the 1 / rev harmonic component of the hub center vibration speed within two adjacent control cycles and the total pitch angle change The relationship between and control input The relationship between them is specifically:

[0031] (7)

[0032] In the formula, To control the target and control input The transfer matrix between them; Represents the control target and control input of the n-1th control cycle.

[0033] Step 3: Based on the relationship between the rotor cone deviation value and the total pitch angle change, and considering the limit of the total pitch angle change, establish the constraint equation of the control input. Specifically:

[0034] The relationship between the rotor cone deviation value and the total pitch angle change is established as follows:

[0035] (8)

[0036] Where: represents the rotor cone deviation value of the nth closed-loop control cycle, that is, formula (2); is the rotor cone deviation value and control input The transfer matrix between them; Indicates that when the control input The rotor cone deviation value at .

[0037] Rotor cone deviation Cannot exceed the set value , then:

[0038] (9)

[0039] In the formula, It is the maximum value of the set rotor cone deviation.

[0040] Substituting equation (8) into equation (9), we can get a set of control input The inequality constraints of can be rearranged as:

[0041] (10)

[0042] In the formula, the inequality coefficient matrix ; The right-hand side of the inequality is ,in for The transpose of the column vector.

[0043] In addition, the control input is the change in the rotor collective pitch angle. When the rotor is working normally, the change in the control input cannot affect the normal rotation of the rotor and cannot change too much. Therefore, there are the following restrictions:

[0044] (11)

[0045] In the formula, The maximum allowable change in the total pitch angle of a single blade.

[0046] Step 4: Use the recursive least squares identification algorithm to calculate the transfer matrix between the control target and the control input, the rotor cone deviation value and the control input, obtain the complete optimal control problem, and transform it into a standard quadratic programming problem. Specifically:

[0047] The present invention adopts the recursive least squares identification method to calculate the transfer matrix in equation (7) and equation (8): and It is suitable for online identification and can update the identification parameters in real time when combined with the closed-loop control algorithm. The method is simple to derive and easy to implement.

[0048] First, determine the parameters to be identified, that is, the transfer matrix in equations (7) and (8) and The calculation process of the two is exactly the same, with the transfer matrix in equation (7) As an example, the identification process is explained. Assume that the matrix The dimension is i×j. Connect them head to tail and get a column vector with (s=i·j) parameters. The parameters to be identified are recorded as , then the discrete observation equation of the n+1th closed-loop control cycle is written as:

[0049] (12)

[0050] in, is the observation vector of the n+1th closed-loop control cycle; is the observation matrix of the n+1th closed-loop control cycle; is the measurement noise of the n+1th closed-loop control cycle; is the parameter to be identified in the n+1th closed-loop control cycle.

[0051] Similarly, the total observation equation for the first n closed-loop control cycles is:

[0052] (13)

[0053] In the formula, Represents the total observation vector of the first n closed-loop control cycles; Represents the total observation matrix of the first n closed-loop control cycles; Represents the total measurement noise of the first n closed-loop control cycles

[0054] According to the least squares estimation criterion, considering the total observation vector of the first n+1 closed-loop control cycles, the parameters to be identified are written as:

[0055] (14)

[0056] In the formula, is the parameter to be identified estimated by the identification algorithm; Represents the total observation matrix of the first n+1 closed-loop control cycles; the superscript -1 represents the inverse of the matrix; the superscript T represents the matrix transpose.

[0057] In formula (14), It is called the information matrix and can be expanded as:

[0058] (15)

[0059] In the formula, is the observation matrix of the nth closed-loop control cycle.

[0060] According to the matrix inversion formula, the parameters to be identified in the n+1th closed-loop control cycle are written as:

[0061] (16)

[0062] In the formula, is the parameter to be identified in the n+1th closed-loop control cycle, i.e., the latest identification result; is the parameter to be identified in the nth closed-loop control cycle, i.e., the identification result of the previous control cycle; is the gain matrix of the nth closed-loop control cycle.

[0063] Formula (16) uses the identification result of the previous control cycle, and the observation matrix only uses the data of the current time step. The calculation dimension is smaller and it can be applied to the parameter estimation of time advancement. As the closed-loop control process progresses, the identification result will become more and more accurate.

[0064] In the transfer matrix and Finally, by combining equations (6), (7), (10), and (11), the control problem of the nth closed-loop control cycle is transformed into a standard quadratic programming problem, which can be summarized as follows:

[0065] (17)

[0066] In the formula, is the function value of the quadratic performance index function; for The quadratic coefficient matrix of ; for The linear coefficient vector of .

[0067] Step 5: Use the effective set method to solve the optimal control problem and determine whether the optimal control input obtained in the current closed-loop control cycle makes the vibration speed 1 / rev harmonic and the rotor cone deviation value meet the error conditions. Specifically:

[0068] First, the effective set method is used to solve the standard quadratic programming problem and calculate the optimal control quantity that satisfies the constraints. This method can accurately handle constraints and is suitable for small and medium-sized problems. The effective set method transforms the quadratic programming problem with inequality constraints into a series of equality constraint sub-problems by dynamically identifying and adjusting constraints (effective sets).

[0069] Then, it is necessary to determine whether the optimal control input obtained in the current closed-loop control cycle makes the vibration speed 1 / rev harmonic meet the error condition, so the current control input Substitute into equation (1) to calculate the rotor cone deviation value: , and then use Fourier transform to calculate the 1 / rev harmonic amplitude of the vibration velocity , and determine whether the following formula is true:

[0070] (18)

[0071] In the formula, is the rotor cone deviation vector; is the maximum value of the rotor cone deviation; are the 1 / rev harmonic amplitudes of the vibration velocity in the three directions x, y, and z of the hub center; They are the maximum values ​​of the 1 / rev harmonic amplitude of the vibration velocity in the three directions of x, y, and z of the hub center.

[0072] If the optimal control quantity of the current closed-loop control cycle can make equation (18) hold, then keep the current optimal control quantity output and monitor in real time whether equation (18) holds; if equation (18) does not hold, change the optimal control quantity of the current closed-loop control cycle to As the control input of the n+1th closed-loop control cycle The initial value is substituted into formula (1), and the control target and control constraint data are recalculated to form the performance index function of step 2. Then the constraint equation of step 3 is calculated, and the recursive least squares identification algorithm of step 4 is executed. Finally, step 5 uses the effective set method to solve the standard quadratic programming problem and obtain the optimal control input of the n+1th closed-loop control cycle. , substitute it into formula (1) and judge again whether formula (18) is true, and repeat the cycle.

[0073] Compared with the prior art, the present invention has the following beneficial effects:

[0074] (1) The present invention proposes a closed-loop control method for rotor vibration that is suitable for both numerical simulation and real tests. The closed-loop control is performed on the rotor cone deviation and the 1 / rev harmonic of the rotor hub center vibration velocity. Compared with the traditional manual adjustment method, the method can actively reduce the amplitude of the 1 / rev harmonic of the vibration velocity of the hub center by identifying the vibration state of the rotor, and at the same time constrain the rotor cone deviation value. It is not only suitable for numerical simulation verification, but can also be conveniently applied to real tests.

[0075] (2) The present invention provides an online identification algorithm based on the recursive least squares method. This method can update the parameters to be identified in real time according to the observed data. The traditional batch least squares method will have a lower identification efficiency as the amount of observed data increases, making it difficult to apply to scenarios with high real-time requirements. This method only uses the current observed data and the identification results of the previous step to keep up to date, which also reduces the amount of calculation and improves the identification efficiency. It can be applied to real-time control experiments.

[0076] In summary, the method described in the present invention has strong operability and feasibility and is convenient for practical application. BRIEF DESCRIPTION OF THE DRAWINGS

[0077] Figure 1 It is a calculation flow chart of the present invention.

[0078] Figure 2 This is a flow chart of the rotor vibration closed-loop control algorithm of the present invention.

[0079] Figure 3 This is a simulation diagram of the BO105 rotor in an embodiment of the present invention.

[0080] Figure 4 The figure is a comparison of the different frequency amplitude bar graphs of the vibration velocity in the x direction of the hub center in Example 1 of the present invention under the uncontrolled and closed-loop controlled states.

[0081] Figure 5 The figure is a comparison of the different frequency amplitude bar graphs of the vibration velocity in the y direction of the hub center in Example 1 of the present invention under the uncontrolled and closed-loop controlled states.

[0082] Figure 6 The figure is a comparison of the different frequency amplitude bar graphs of the vibration velocity in the z direction of the hub center in Example 1 of the present invention under the uncontrolled and closed-loop controlled states.

[0083] Figure 7This is a curve diagram of the maximum value of the rotor cone deviation in Example 1 of the present invention before, during and after closed-loop control.

[0084] Figure 8 This is a diagram of the rotor vibration closed-loop control test system in Example 2 of the present invention.

[0085] Fig. 9 This is a curve showing the variation of the harmonic amplitude of the vibration velocity 1 / rev of the hub center in the x-direction before, during and after control in the rotor vibration closed-loop control test in Example 2 of the present invention.

[0086] Fig.10 This is a curve showing the variation of the harmonic amplitude of the vibration velocity 1 / rev of the hub center in the y direction before, during and after control in the rotor vibration closed-loop control test in Example 2 of the present invention.

[0087] Fig.11 This is a curve showing the variation of the harmonic amplitude of the vibration velocity 1 / rev of the hub center in the z direction before, during and after control in the rotor vibration closed-loop control test in Example 2 of the present invention.

[0088] Fig.12 This is a curve showing the variation of the maximum value of the rotor cone deviation before, during and after control in the rotor vibration closed-loop control test in Example 2 of the present invention. DETAILED DESCRIPTION

[0089] In order to make the purpose, technical scheme and advantages of the implementation of the present invention clearer, the technical scheme in the embodiment of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiment of the present invention. Obviously, the described embodiment is a part of the embodiment of the present invention, not all 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.

[0090] Figure 1 The calculation flow chart of the present invention is given. Figure 2 This is a flow chart of the closed-loop control algorithm of the present invention, and the following example verification is carried out based on it.

[0091] Example 1: Numerical simulation verification of BO105 helicopter rotor vibration closed-loop control based on cone constraints

[0092] Step 1: Set the dissimilar physical quantities between the blades and establish the rotor aeroelastic coupling model, thereby calculating the rotor cone deviation value and the harmonic component of the hub center vibration velocity 1 / rev. Specifically:

[0093] First, set the dissimilar physical quantities between the blades. Wear, deformation, aging and other loss problems that occur during the use of the blades are all random phenomena, which are reflected in the physical quantity as random changes to the physical quantity. This embodiment uses the rand random function in MATLAB to generate random change coefficients for the following physical quantities, and for the flapping stiffness , shimmy stiffness , torsional stiffness The three stiffness parameters are changed within ±5%, the linear density dm is changed within ±1.2%, and the lift coefficient , drag coefficient , pitch moment coefficient The other three aerodynamic parameters vary within ±8%. Other parameters of BO105 helicopter are shown in Table 1. Figure 3 Numerical simulation diagram of the BO105 helicopter rotor established for this embodiment.

[0094] Table 1 BO105 helicopter rotor parameters

[0095]

[0096] Then, based on the multi-body dynamics theory and considering the dissimilar physical quantities of the blades, the rotor aeroelastic coupling model is established as shown in formula (1). The generalized -α method is used to solve it, and the rotor cone deviation value as shown in formula (2) can be obtained.

[0097] Furthermore, the hub center vibration load is solved as shown in equation (3).

[0098] Finally, the hub center vibration load is converted into the hub center vibration velocity using equation (4). The conversion matrix in equation (4) used in this embodiment is: for:

[0099] (19)

[0100] In the formula, the elements in the last column of the matrix are all 0 in this embodiment because the last column corresponds to the hub center torque in formula (4). The influence of this load on the vibration velocity of the hub center is very small and can generally be ignored.

[0101] The calculated vibration velocity The fast Fourier transform function FFT in MATLAB software is used to perform Fourier transform and extract the 1 / rev harmonic component, thereby obtaining the control target data of formula (5).

[0102] Step 2: Taking the blade total pitch angle change as the control input and the rotor vibration speed 1 / rev harmonic component as the control target, establish the quadratic performance index function of the control problem. Specifically:

[0103] Based on the optimal control theory, the performance index function is set as shown in formula (6), where the weight matrix Take them as:

[0104] (20)

[0105] In the formula, are 6-dimensional and 3-dimensional unit matrices respectively. Multiply the identity matrix by It is to increase the weight value of the control target to improve the vibration control effect.

[0106] In addition, it is necessary to establish the 1 / rev harmonic component of the hub center vibration speed within two adjacent control cycles and the total angle change The relationship between is shown in formula (7).

[0107] Step 3: Based on the relationship between the rotor cone deviation value and the total pitch angle change, and considering the limit of the total pitch angle change, establish the constraint equation of the control input. Specifically:

[0108] The relationship between the rotor cone deviation value and the total pitch angle change is established as shown in formula (8); Cannot exceed the set value , in which in this embodiment , and the cone deviation constraint inequality is organized into the form of equation (10).

[0109] In addition, the control input is the change in the rotor collective pitch angle, which cannot change too much. Therefore, there is an inequality constraint as shown in equation (11), where the maximum change allowed in the collective pitch angle of a single blade in this embodiment is Take 0.175°.

[0110] Step 4: Use the recursive least squares identification algorithm to calculate the transfer matrix between the control target and the control input, the rotor cone deviation value and the control input, obtain the complete optimal control problem, and transform it into a standard quadratic programming problem. Specifically:

[0111] The recursive least squares identification method, i.e., equation (16), is used to calculate the transfer matrix in equations (7) and (8): and Then, substituting into equation (17), we obtain the standard quadratic programming problem for the control problem of the nth closed-loop control cycle.

[0112] Step 5: Use the effective set method to solve the optimal control problem and determine whether the optimal control input obtained in the current closed-loop control cycle makes the vibration speed 1 / rev harmonic and the rotor cone deviation value meet the error conditions. Specifically:

[0113] First, the standard quadratic programming problem is solved using the effective set method to calculate the optimal control quantity that satisfies the constraints in the nth closed-loop control cycle. Then, the current control input Substituting into equation (1), the rotor cone deviation value is , and then use Fourier transform to calculate the 1 / rev harmonic amplitude of the vibration velocity , determine the optimal control input obtained in the current closed-loop control cycle Whether the vibration speed 1 / rev harmonic satisfies the error condition, that is, equation (18), where the error limit in equation (18) is as follows:

[0114] (twenty one)

[0115] If the optimal control quantity of the current closed-loop control cycle can make equation (18) hold, then keep the current optimal control quantity output and monitor in real time whether equation (18) holds; if equation (18) does not hold, change the optimal control quantity of the current closed-loop control cycle to As the control input of the n+1th closed-loop control cycle The initial value is substituted into formula (1), and the control target and control constraint data are recalculated to form the performance index function of step 2. Then the constraint equation of step 3 is calculated, and the recursive least squares identification algorithm of step 4 is executed. Finally, step 5 uses the effective set method to solve the standard quadratic programming problem and obtain the optimal control input of the n+1th closed-loop control cycle. , and judge whether formula (18) is established, and repeat this cycle.

[0116] Finally, the above steps are performed so that the 1 / rev amplitude of the vibration velocity in the three directions of the hub center and the rotor cone deviation value meet the error conditions, and the following is obtained: Figure 4~Figure 7 The result graph is shown.

[0117] like Figure 4~Figure 6 As shown, the amplitude comparison of the first four frequencies of the vibration harmonics in the three directions before and after closed-loop control is shown. First, it can be clearly seen before control that after considering the difference in the dissimilar physical quantities of the blades in step 1, the 1 / rev frequency amplitude is significantly increased, and non-integer multiples of N / rev 2 / rev and 3 / rev harmonic amplitudes are also generated. After applying the closed-loop control algorithm proposed in this article, the 1 / rev harmonic amplitude of the vibration speed in the three directions of x, y, and z can be greatly reduced by more than 80%, which shows that the control strategy proposed in this invention can correctly adjust the vibration imbalance characteristics generated by dissimilar blades.

[0118] like Figure 7The changing trend of the maximum value of the rotor cone deviation during the closed-loop control process is given. From the curve before control, it can be seen that the blade dissimilarity produces a large deviation value. During the control process, with the execution of online identification and closed-loop control, the identification results gradually become more accurate, and the cone deviation value also decreases rapidly, eventually reaching near the set maximum value, which is about 54% lower than before control.

[0119] Example 2: Experimental verification of helicopter rotor vibration closed-loop control based on cone constraints

[0120] Similarly, the above steps 1 to 5 are performed to experimentally verify the helicopter vibration closed-loop control method based on cone constraints. There are two differences from Example 1: The first difference is that in step 1, the rotor aeroelastic coupling model used to calculate the 1 / rev harmonic of the hub vibration speed and the rotor cone deviation value is replaced by the vibration sensor and cone sensor in the test system, that is, the experimental data collected by the sensor is directly used for subsequent online identification and closed-loop control calculations; the second difference is that in order to realize the combination of the algorithm and the test hardware, the online identification algorithm and closed-loop control algorithm in steps 2 to 5 need to be adapted into real-time simulation software using the Simulink simulation platform, and embedded in the test hardware for experimental verification, while other steps remain unchanged. Figure 8 The hardware and software architecture of the experimental system is shown in Figure 1. Figure 9~Figure 12 Result graph.

[0121] Figure 9~Figure 11 The variation curves of the vibration velocity 1 / rev amplitude in the x, y, and z directions before, during, and after closed-loop control are given. The 1 / rev amplitude in the three directions dropped significantly during the control process. After the control was finally stabilized, the control effect in the three directions was more than 60%, and the y direction achieved a 72% reduction. Fig.12 The curve of the maximum value of the rotor cone deviation before, during and after the control is shown. Although the rotor may be affected by a similar accidental gust of wind in the first half of the period after the control input is stable, and there is a large fluctuation, the average value in the second half is also greatly reduced, gradually reaching the constraint limit value, and the overall average value also reaches 61% of the control effect. Therefore, the test results also prove the effectiveness of the present invention.

[0122] The above-described embodiments merely express the implementation methods of the present invention, but they cannot be understood as limiting the scope of the patent of the present invention. It should be pointed out that for those skilled in the art, several modifications and improvements can be made without departing from the concept of the present invention, which all belong to the protection scope of the present invention.

Claims

1. A helicopter vibration closed-loop control method based on cone constraints, characterized in that: The following steps are involved: Step 1: Set the dissimilar physical quantities between the blades and establish the rotor aeroelastic coupling model, thereby calculating the harmonic component of the rotor hub center vibration velocity 1 / rev and the rotor cone deviation value; Step 2: Taking the blade total pitch angle change as the control input and the rotor vibration speed 1 / rev harmonic component as the control target, establish the quadratic performance index function of the control problem; Step 3: Based on the relationship between the rotor cone deviation value and the total pitch angle change, and considering the limit of the total pitch angle change, establish the constraint equation of the control input; Step 4: Use the recursive least squares identification algorithm to calculate the transfer matrix between the control target and the control input, the rotor cone deviation value and the control input, obtain the complete optimal control problem, and transform it into a standard quadratic programming problem; Step 5: Use the effective set method to solve the optimal control problem and determine whether the optimal control input obtained in the current closed-loop control cycle makes the vibration speed 1 / rev harmonic and the rotor cone deviation value meet the error conditions.

2. The helicopter vibration closed-loop control method based on cone constraint according to claim 1, characterized in that: The step 1 is specifically as follows: Step 1.1, set the dissimilar physical quantities between the blades, including the flapping stiffness , shimmy stiffness , torsional stiffness , linear density dm, lift coefficient , drag coefficient , pitch moment coefficient , these physical quantities will affect the dynamic response of the rotor, resulting in an increase in the harmonic component of the vibration velocity 1 / rev at the hub center; Step 1.2, based on the multi-body dynamics theory and the blade dissimilar physical quantities, a rotor aeroelastic coupling model is established, and the rotor cone deviation value is solved; Consider the constraint equation , we obtain a set of differential algebraic equations (DAEs) and extract the generalized acceleration and generalized coordinates The linear term coefficient of is obtained, and the overall aeroelastic coupling model of the rotor is obtained, which is organized into the standard differential algebraic equation form as follows: (1) In the formula, They are the mass matrix, stiffness matrix, generalized nonlinear term and rotor aerodynamic term of the rotor aeroelastic coupling model; Indicates the rotor's total pitch angle and its variation, which directly affects the aerodynamic terms ,in It is the control input; are the generalized coordinates and generalized acceleration respectively; The superscript T represents the matrix transpose symbol; Represents the constraint equations of the rotor aeroelastic coupling model; represents the transpose of the Jacobian matrix of the constraint equation, that is ; is the Lagrange multiplier; The generalized-α method is used to solve the rotor aeroelastic coupling model shown in formula (1) to obtain the dynamic response of each blade: , and then calculate the tip height of each blade , assuming the number of rotor blades is , arbitrarily select the e-th blade as the reference blade, and calculate the rotor cone deviation value, that is, the difference between the tip height of other blades and the reference blade, as follows: (2) In the formula, Represents the rotor cone deviation value vector; Indicates the 1st, 2nd, ..., The tip height of the propeller blade; represents the tip height value of the e-th blade; Step 1.3, solve the hub center vibration load, which is the superposition of the propeller root load, as follows: (3) In the formula, is the azimuth angle of the mth blade; is the number of blades; the superscript H indicates the hub; the superscript m indicates the mth blade; the subscripts x, y, and z indicate the force or moment in the direction of the x, y, and z coordinate axes; They represent the root shear force of the mth blade in the x, y, and z directions respectively; They represent the root moments of the mth blade in the x, y, and z directions respectively; Respectively represent the vibration shear force of the hub center in the x, y, and z directions; Respectively represent the vibration moments of the hub center in the x, y, and z directions; Step 1.4, convert the hub center vibration load into the hub center vibration velocity, and use Fourier transform to extract the 1 / rev harmonic component of the vibration velocity; in engineering, vibration is measured using a vibration sensor, and the vibration velocity in the three directions of x, y, and z at the hub center is used to characterize the magnitude of the vibration; therefore, combined with the relationship between vibration load and vibration velocity, the hub vibration load is converted into vibration velocity as follows: (4) In the formula, They represent the vibration speeds in the three directions of x, y, and z at the hub center respectively; the superscript T represents the transposition symbol; is a transformation matrix of 3×6 dimensions; The calculated vibration velocity Use Fourier transform to do frequency domain analysis and extract the 1 / rev harmonic component, that is: The 1 / rev harmonic components include sinusoidal components and cosine components ; The 1 / rev harmonic components include sinusoidal components and cosine components ; The 1 / rev harmonic components include sinusoidal components and cosine components ; The control target is obtained as follows: (5) In the formula, It represents the 1 / rev harmonic component of the vibration velocity, which is also the control target. The superscript T represents the transposition sign.

3. The helicopter vibration closed-loop control method based on cone constraint according to claim 2, characterized in that: The step 2 is specifically as follows: Step 2.1, based on the optimal control theory, establish the performance index function of the control problem, and set the performance index function as: (6) Where: Represents the performance indicator function value; is the control target, i.e., the harmonic component of the rotor vibration speed 1 / rev in equation (5); the subscript n represents the nth closed-loop control cycle; is the control input of the nth closed-loop control cycle, i.e., the total pitch angle change in equation (1); is the change in the control input from the n-1th closed-loop control cycle to the nth closed-loop control cycle; Control objectives , control input , Control input change The corresponding weight matrix; Step 2.2: Establish the harmonic component of the hub center vibration speed 1 / rev in two adjacent control cycles and the total pitch angle change The relationship between and control input The relationship between them is specifically: (7) In the formula, To control the target and control input The transfer matrix between them; Represents the control target and control input of the n-1th control cycle.

4. The helicopter vibration closed-loop control method based on cone constraint according to claim 3 is characterized in that: The step 3 is specifically as follows: The relationship between the rotor cone deviation value and the total pitch angle change is established as follows: (8) Where: represents the rotor cone deviation value of the nth closed-loop control cycle, that is, formula (2); is the rotor cone deviation value and control input The transfer matrix between them; Indicates that when the control input Rotor cone deviation value at ; Rotor cone deviation Cannot exceed the set value , then: (9) In the formula, is the maximum value of the set rotor cone deviation value; Substituting equation (8) into equation (9), we get a set of control input The inequality constraints are: (10) In the formula, the inequality coefficient matrix ; The right-hand side of the inequality is ,in for The transpose of the column vector.

5. The helicopter vibration closed-loop control method based on cone constraint according to claim 4, characterized in that: In step 3, the control input is the change in the rotor collective pitch angle. When the rotor is working normally, the change in control input cannot affect the normal rotation of the rotor. There are the following restrictions: (11) In the formula, The maximum allowable change in the total pitch angle of a single blade.

6. The helicopter vibration closed-loop control method based on cone constraint according to claim 5, characterized in that: The step 4 is specifically as follows: Using the recursive least squares identification method, the transfer matrix in equations (7) and (8) is calculated: and ; First, determine the parameters to be identified, that is, the transfer matrix in equations (7) and (8) and The calculation process of the two is exactly the same, with the transfer matrix in formula (7) To illustrate; assume the matrix The dimension is i×j. Connect them head to tail and get a column vector with (s=i·j) parameters. The parameters to be identified are recorded as , then the discrete observation equation of the n+1th closed-loop control cycle is written as: (12) in, is the observation vector of the n+1th closed-loop control cycle; is the observation matrix of the n+1th closed-loop control cycle; is the measurement noise of the n+1th closed-loop control cycle; is the parameter to be identified in the n+1th closed-loop control cycle; Then the total observation equation for the first n closed-loop control cycles is: (13) In the formula, Represents the total observation vector of the first n closed-loop control cycles; Represents the total observation matrix of the first n closed-loop control cycles; Represents the total measurement noise of the first n closed-loop control cycles According to the least squares estimation criterion, considering the total observation vector of the first n+1 closed-loop control cycles, the parameters to be identified are written as: (14) In the formula, is the parameter to be identified estimated by the identification algorithm; It represents the total observation matrix of the first n+1 closed-loop control cycles; the superscript -1 represents the inverse of the matrix; the superscript T represents the matrix transpose; In formula (14), It is called the information matrix and is expanded as: (15) In the formula, is the observation matrix of the nth closed-loop control cycle; According to the matrix inversion formula, the parameters to be identified in the n+1th closed-loop control cycle are written as: (16) In the formula, is the parameter to be identified in the n+1th closed-loop control cycle, i.e., the latest identification result; is the parameter to be identified in the nth closed-loop control cycle, i.e., the identification result of the previous control cycle; is the gain matrix of the nth closed-loop control cycle; In the transfer matrix and Finally, by combining equations (6), (7), (10), and (11), the control problem of the nth closed-loop control cycle is transformed into a standard quadratic programming problem, which can be summarized as follows: (17) In the formula, is the function value of the quadratic performance index function; for The quadratic coefficient matrix of ; for The linear coefficient vector of .

7. A helicopter vibration closed-loop control method based on cone constraint according to claim 6, characterized in that: The step 5 is specifically as follows: Step 5.1, first, use the effective set method to solve the standard quadratic programming problem and calculate the optimal control quantity that satisfies the constraints; Step 5.2, then, it is necessary to determine whether the optimal control input obtained in the current closed-loop control cycle makes the vibration speed 1 / rev harmonic meet the error condition, so the current control input Substitute into equation (1) to calculate the rotor cone deviation value: , and then use Fourier transform to calculate the 1 / rev harmonic amplitude of the vibration velocity , and determine whether the following formula is true: (18) In the formula, is the rotor cone deviation vector; is the maximum value of the rotor cone deviation; are the 1 / rev harmonic amplitudes of the vibration velocity in the three directions x, y, and z of the hub center; are the maximum values ​​of the 1 / rev harmonic amplitude of the vibration velocity in the three directions x, y, and z of the hub center; If the optimal control quantity of the current closed-loop control cycle can make equation (18) hold, then keep the current optimal control quantity output and monitor in real time whether equation (18) holds; if equation (18) does not hold, change the optimal control quantity of the current closed-loop control cycle to As the control input of the n+1th closed-loop control cycle The initial value is substituted into formula (1), and the control target and control constraint data are recalculated to form the performance index function of step 2. Then the constraint equation of step 3 is calculated, and the recursive least squares identification algorithm of step 4 is executed. Finally, step 5 uses the effective set method to solve the standard quadratic programming problem and obtain the optimal control input of the n+1th closed-loop control cycle. , substitute it into formula (1) and judge again whether formula (18) is true, and repeat the cycle.

Images