A closed-loop control method for helicopter vibration based on cone constraint

Abstract

The present invention provides a closed-loop control method for helicopter vibration based on cone constraint, belonging to the field of vibration control. First, an aeroelastic coupling model of the rotor is established according to the structural characteristics of the helicopter rotor, and at the same time, dissimilar physical quantities between the blades are set, and the blade cone mutual difference value and the hub 1 / rev vibration harmonic are calculated. Secondly, a quadratic performance index function of the control problem is established with the blade collective pitch change amount and the rotor 1 / rev vibration harmonic component as the control input and output. Thirdly, according to the relationship between the cone and the blade collective pitch, and considering the limitation of the collective pitch change amount at the same time, a constraint equation of the control input is established. Fourthly, the transfer matrix between the control input and output is estimated to obtain the optimal control problem, which is transformed into a quadratic programming problem. Finally, the optimal control quantity is obtained by solving. The present invention is not only applicable to numerical simulation verification, but also can be applied to real-time control experiments, improving the identification efficiency, having strong operability and feasibility, and being 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] Since the active vibration control technology has high control potential, it has been extensively studied. Among them, Chinese invention patent [CN 116424552 B] proposes an active torsional blade vibration control method and system. This invention changes the aerodynamic layout of the blade through the force-electric coupling effect of the actively twisted blade to improve the aerodynamic load, thereby reducing the hub vibration, and can calculate the optimal torsional change amount according to the vibration level of the helicopter. However, the actuator installed in this method is integrated into the blade, that is, it is necessary to embed active fiber composite materials or thick fiber piezoelectric materials in the blade skin. This not only requires remanufacturing the blade, but also the embedded actuator will affect the material properties of the blade, and even requires redesigning the blade structure, which will bring a great deal of extra work. In addition, this control method does not consider the rotor cone deviation, and the performance of the controller is not fully exerted.

[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 values. Therefore, the present invention proposes an active control method that can fully exert the vibration control effect. This method can consider the rotor cone deviation value and does not make any changes to the blade. Only a slight adjustment of the blade pitch angle is required to maximize the performance of the active vibration control algorithm. Summary of the Invention

[0006] To solve the above technical problems, the present invention proposes a closed-loop control method for helicopter vibration based on cone constraint. According to the independent blade control (IBC) idea, the change amount of the collective pitch angle of the blade is used as the control input, the control target is the 1 / rev harmonic component of the vibration at the rotor hub center, 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 real experiments. Therefore, when performing numerical simulation in the present invention, it is necessary to first establish a rotor aeroelastic coupling model with blade dissimilarity to calculate the rotor cone deviation value and the 1 / rev harmonic component of the vibration velocity at the hub center in the control method. In the experiment, a cone sensor and a vibration sensor are directly used to measure the two data of the test rotor system. In addition, the vibration measured by the vibration sensor in the experiment is measured by vibration velocity. Therefore, in order to maintain the unity of the control target, it is also necessary to convert it to vibration velocity in the numerical simulation stage. Then, after obtaining the control target data, a quadratic performance index function for this vibration problem is established. Using the change amount of the collective pitch angle of the rotor blade as the control input, the 1 / rev harmonic component of the vibration velocity at the rotor hub center as the control target, and the rotor cone deviation value as the control constraint, a constrained optimization problem is constructed. The least squares identification algorithm is used to estimate the transfer matrix, and the quadratic programming optimization algorithm is used to solve for the optimal control quantity.

[0007] To achieve the above object, the technical solution adopted by the present invention is as follows:

[0008] A closed-loop control method for helicopter vibration based on cone constraint. The method includes the following steps:

[0009] Step 1: Set dissimilar physical quantities between the blades, establish a rotor aeroelastic coupling model, and calculate the 1 / rev harmonic component of the vibration velocity at the rotor hub center 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 two data items are directly measured by a cone sensor and a vibration sensor to replace the numerical calculation of this step. Specifically:

[0010] First, set dissimilar physical quantities between the blades. Wear, deformation, aging and other loss problems occurring 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 flap stiffness , the lead-lag stiffness , the torsional stiffness , the linear density dm, the lift coefficient , the drag coefficient , and the pitching moment coefficient . These physical quantities will all change to varying degrees after the blades are worn, etc., and will further affect the dynamic response of the rotor, resulting in an increase in the 1 / rev harmonic component of the vibration velocity at the rotor hub center.

[0011] Then, based on the multi-body dynamics theory, considering the dissimilar physical quantities of the blades, establish a rotor aeroelastic coupling model and solve the rotor cone deviation value. Considering the constraint equation , a set of differential-algebraic equations (DAEs) is obtained, and the linear term coefficients of the generalized acceleration and the generalized coordinates are extracted to obtain the overall rotor aeroelastic coupling model, which is arranged in the standard form of differential-algebraic equations as follows:

[0012] (1)

[0013] In the formula, are respectively the mass matrix, stiffness matrix, generalized nonlinear term and rotor aerodynamic force term of the rotor aeroelastic coupling model; represents the collective pitch angle of the rotor and its change amount, which directly affects the aerodynamic force term , where is the control input; are respectively the generalized coordinates and generalized acceleration; the superscript T represents the matrix transpose symbol; represents the constraint equation of the rotor aeroelastic coupling model; represents the transpose of the Jacobian matrix of the constraint equation, that is, ; is the Lagrange multiplier.

[0014] By using the generalized-α method to solve the rotor aeroelastic coupling model in Equation (1), the dynamic response of each blade can be obtained. Furthermore, the tip height value of each blade can be calculated. Assume that 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; represents the tip height values of the 1st, 2nd, …, th blades; represents the tip height value of the e-th blade.

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

[0018] (3)

[0019] In the formula, is the azimuth angle of the m-th blade; is the number of blades; the superscript H represents the hub; the superscript m represents the m-th blade; the subscripts x, y, z represent the forces or torques in the x, y, z coordinate axis directions; respectively represent the blade root shear forces of the m-th blade in the x, y, z directions; respectively represent the blade root torques of the m-th blade in the x, y, z directions; respectively represent the vibration shear forces at the hub center in the x, y, z directions; respectively represent the vibration torques at the hub center in the x, y, z directions.

[0020] Finally, convert the vibration load at the hub center to the vibration velocity at the hub center, and use Fourier transform to extract the 1 / rev harmonic component of the vibration velocity. In engineering, vibration is generally measured using vibration sensors, and the vibration magnitudes are characterized by the vibration velocities in the x, y, z directions at the hub center. Therefore, combining the relationship between the vibration load and the vibration velocity, convert the hub vibration load to the vibration velocity, as follows:

[0021] (4)

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

[0023] For the calculated vibration velocity perform frequency domain analysis using Fourier transform to extract the 1 / rev harmonic component, that is: The 1 / rev harmonic component of includes the sine component ; The 1 / rev harmonic component of includes the sine component ; The 1 / rev harmonic component of includes the sine component . Thus, the control objective is obtained:

[0024] (5)

[0025] In the formula, represents the 1 / rev harmonic component of the vibration velocity, which is also the control objective; the superscript T represents the transpose symbol.

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

[0027] Based on the optimal control theory, establish a performance index function for the control problem. Generally, it is constructed as a quadratic function, which has characteristics such as clear physical meaning and convenient mathematical processing. Therefore, set the performance index function as:

[0028] (6)

[0029] In the formula: represents the value of the performance index function; is the control objective, that is, the 1 / rev harmonic component of the rotor vibration velocity in formula (5); the subscript n represents the nth closed-loop control period; is the control input in the nth closed-loop control period, that is, the change in the total pitch angle in formula (1); is the change in the control input from the (n - 1)th closed-loop control period to the nth closed-loop control period; are respectively the control objective , the control input , and the control input change corresponding weight matrices.

[0030] In addition, it is necessary to establish the relationship between the 1 / rev harmonic component of the vibration velocity at the hub center and the change in the collective pitch angle within two adjacent control cycles which is also the control objective and the control input Specifically: Specifically:

[0031] (7)

[0032] In the formula, is the transfer matrix between the control objective and the control input ; represents the control objective and control input in the (n - 1)-th control cycle.

[0033] Step 3: According to the relationship between the rotor cone deviation value and the change in the collective pitch angle, and considering the limitation of the change in the collective pitch angle, establish the constraint equation of the control input. Specifically:

[0034] Establish the relationship between the rotor cone deviation value and the change in the collective pitch angle as follows:

[0035] (8)

[0036] In the formula: represents the rotor cone deviation value in the n-th closed-loop control cycle, that is, Equation (2); is the transfer matrix between the rotor cone deviation value and the control input ; represents the rotor cone deviation value when the control input is .

[0037] The rotor cone deviation value shall not exceed the set value , so there is:

[0038] (9)

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

[0040] Substitute Equation (8) into Equation (9) to obtain a set of inequality constraints on the control input , which are sorted out as:

[0041] (10)

[0042] In the formula, the inequality coefficient matrix is ; the right-hand side term of the inequality is , where is the transpose of the column vector composed of

[0043] In addition, the control input is the change amount of the total rotor pitch angle. When the rotor is working normally, the change of 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] where is the maximum allowable change amount of the total pitch angle of a single blade.

[0046] Step 4: Use the recursive least squares identification algorithm to calculate the transfer matrices between the control objective and the control input, and between 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 matrices in formulas (7) and (8) and . It is applicable to online identification, can update the identification parameters in real time in combination with the closed-loop control algorithm, and this method has simple derivation and is easy to implement.

[0048] First, determine the parameters to be identified, that is, the transfer matrices in formulas (7) and (8) and . The calculation processes of the two are exactly the same. Taking the transfer matrix in formula (7) as an example, the identification process is described. Assume that the matrix has the dimension of i×j, and connect its columns end to end to obtain a column vector with (s = i·j) parameters, and obtain the parameter to be identified, denoted as , then the discrete observation equation for the (n + 1)-th closed-loop control period is written as:

[0049] (12)

[0050] where is the observation vector for the (n + 1)-th closed-loop control period; is the observation matrix for the (n + 1)-th closed-loop control period; is the measurement noise for the (n + 1)-th closed-loop control period; is the parameter to be identified for the (n + 1)-th closed-loop control period.

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

[0052] (13)

[0053] where denotes the total observation vector for the first n closed-loop control cycles; denotes the total observation matrix for the first n closed-loop control cycles; denotes the total measurement noise for the first n closed-loop control cycles

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

[0055] (14)

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

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

[0058] (15)

[0059] where is the observation matrix for the nth closed-loop control cycle.

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

[0061] (16)

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

[0063] Equation (16) uses the identification result of the previous control cycle, and the observation matrix only uses the data of the current time step, with a smaller calculation dimension, which 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] After obtaining the transfer matrices and combining equations (6), (7), (10), (11), the control problem for the nth closed-loop control cycle is transformed into a standard quadratic programming problem and organized as:

[0065] (17)

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

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

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

[0069] Then, it is necessary to determine whether the optimal control input obtained in the current closed-loop control period makes the vibration speed 1 / rev harmonic meet the error conditions. Therefore, substitute the current control input into Equation (1) to calculate the rotor cone deviation value , and then use the Fourier transform to calculate the 1 / rev harmonic amplitude of the vibration speed, and determine whether the following formula holds:

[0070] (18)

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

[0072] If the optimal control quantity in the current closed-loop control period 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, take the optimal control quantity in the current closed-loop control period as the control input initial value of the (n + 1)-th closed-loop control period, substitute it into Equation (1), recalculate the data of the control objective and control constraints, form the performance index function in Step 2, then calculate the constraint equation in Step 3, execute the recursive least squares identification algorithm in Step 4, and finally use the active set method in Step 5 to solve the standard quadratic programming problem to obtain the optimal control input of the (n + 1)-th closed-loop control period Substitute it into Equation (1) and then determine again whether Equation (18) holds. Repeat this process in a loop.

[0073] The beneficial effects of the present invention compared with the prior art are as follows:

[0074] (1) A closed-loop control method for rotor vibration that is applicable to both numerical simulation and real tests proposed by the present invention performs closed-loop control on the rotor cone deviation and the 1 / rev harmonic of the vibration velocity at the rotor hub center. Compared with the traditional manual adjustment method, it can identify the vibration state of the rotor, actively reduce the amplitude of the 1 / rev harmonic of the vibration velocity at the hub center, and at the same time constrain the rotor cone deviation value. It is not only applicable to numerical simulation verification but also can 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. For the traditional batch least squares method, the identification efficiency decreases as the observed data increases, and it is difficult to be applied to scenarios with high real-time requirements. This method only uses the current observed data and the identification result of the previous step to maintain the update, which also reduces the computational amount, improves the identification efficiency, and can be applied to real-time control tests.

[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 is the calculation flow chart of the present invention.

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

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

[0080] Figure 4 is the comparison of the bar charts of the vibration velocity amplitudes at different frequencies in the x direction of the hub center in the non-control and closed-loop control states in Embodiment 1 of the present invention.

[0081] Figure 5 is the comparison of the bar charts of the vibration velocity amplitudes at different frequencies in the y direction of the hub center in the non-control and closed-loop control states in Embodiment 1 of the present invention.

[0082] Figure 6 is the comparison of the bar charts of the vibration velocity amplitudes at different frequencies in the z direction of the hub center in the non-control and closed-loop control states in Embodiment 1 of the present invention.

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

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

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

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

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

[0088] Figure 12 This is the change curve of the maximum value of the rotor cone deviation of the rotor vibration closed-loop control test in Embodiment 2 of the present invention before, during, and after control. Detailed implementation mode

[0089] To make the purpose, technical solutions, and advantages of the implementation of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, rather than all of them. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative work fall within the protection scope of the present invention.

[0090] Figure 1 The calculation flowchart of the present invention is given. Figure 2 This is the flowchart of the closed-loop control algorithm of the present invention, and the following embodiment verification is carried out accordingly.

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

[0092] Step 1: Set the dissimilar physical quantities between the blades, establish aeroelastic coupling model of the rotor, and calculate the rotor cone deviation value and the 1 / rev harmonic component of the hub center vibration velocity. 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. Physically, this is manifested as a random change in this physical quantity. In this embodiment, the rand random function in matlab is used to generate the coefficients of random changes for the following physical quantities, for the flapping stiffness and the lag stiffness and the torsional stiffness and other three stiffness parameters to produce changes within ±5%, for the linear density dm to produce changes within ±1.2%, and for the lift coefficient , the drag coefficient , and the pitching moment coefficient and other three aerodynamic parameters to produce changes within ±8%. The other parameters of the BO105 helicopter are shown in Table 1. Figure 3 is the numerical simulation diagram of the BO105 helicopter rotor established in this embodiment.

[0094] Table 1 BO105 helicopter rotor parameters

[0095]

[0096] Then, based on the multi-body dynamics theory, considering the dissimilar physical quantities of the blades, a rotor aeroelastic coupling model as shown in Equation (1) is established and solved using the generalized-α method to obtain the rotor cone deviation value as shown in Equation (2).

[0097] Furthermore, solve the vibration load at the hub center as shown in Equation (3).

[0098] Finally, use Equation (4) to convert the vibration load at the hub center into the vibration velocity at the hub center. The transformation matrix used in Equation (4) in this embodiment is:

[0099] (19)

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

[0101] For the calculated vibration velocity use the fast Fourier transform function FFT in the matlab software to perform Fourier transform and extract the 1 / rev harmonic component, from which the control target data as shown in Equation (5) can be obtained.

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

[0103] Based on the optimal control theory, a performance index function as shown in Equation (6) is set, where the weight matrices are taken as follows respectively:

[0104] (20)

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

[0106] In addition, it is necessary to establish the relationship between the 1 / rev harmonic component of the vibration velocity of the hub center and the change in the collective pitch angle in two adjacent control cycles, as shown in Equation (7).

[0107] Step 3: According to the relationship between the rotor cone deviation value and the change in the collective pitch angle, and considering the limitation of the change in the collective pitch angle, establish the constraint equation of the control input. Specifically:

[0108] Establish the relationship between the rotor cone deviation value and the change in the collective pitch angle, as shown in Equation (8); the rotor cone deviation value shall not exceed the set value , where in this embodiment is taken, and the cone deviation constraint inequality is arranged into the form of Equation (10).

[0109] In addition, the control input is the change in the collective pitch angle of the rotor and cannot change too much. Therefore, there is an inequality constraint as shown in Equation (11), where the maximum allowable change in the collective pitch angle of a single blade is taken as 0.175°.

[0110] Step 4: Use the recursive least squares identification algorithm to calculate the transfer matrices between the control objective and the control input, and between 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] Adopt the recursive least squares identification method, that is, Equation (16), to calculate the transfer matrices and in Equations (7) and (8). Then, substitute them into Equation (17) to obtain the standard quadratic programming problem of the control problem in the nth closed-loop control cycle.

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

[0113] First, the active set method is used to solve the standard quadratic programming problem to calculate the optimal control quantity that satisfies the constraints in the nth closed-loop control period. . Then, substitute the current control input into Equation (1) to obtain the rotor cone deviation value . Next, use the Fourier transform to calculate the 1 / rev harmonic amplitude of the vibration velocity . Then, determine whether the optimal control input obtained in the current closed-loop control period makes the 1 / rev harmonic of the vibration velocity satisfy the error condition of Equation (18), where the error limit in Equation (18) is taken as follows:

[0114] (21)

[0115] If the optimal control quantity in the current closed-loop control period can make Equation (18) hold, keep the output of the current optimal control quantity and monitor in real time whether Equation (18) holds; if Equation (18) does not hold, use the optimal control quantity in the current closed-loop control period as the initial value of the control input in the (n + 1)th closed-loop control period, substitute it into Equation (1), recalculate the data of the control objective and control constraints to form the performance index function in Step 2, then calculate the constraint equation in Step 3, execute the recursive least squares identification algorithm in Step 4, and finally use the active set method in Step 5 to solve the standard quadratic programming problem to obtain the optimal control input in the (n + 1)th closed-loop control period, and determine whether Equation (18) holds, and so on in a loop.

[0116] Finally, execute the above steps to make the 1 / rev amplitudes of the vibration velocities in the three directions of the hub center and the rotor cone deviation value satisfy the error conditions, and obtain the result graph as Figures 4 to 7 shown.

[0117] As Figures 4 to 6 shown, it is the comparison of the amplitudes of the first four-order frequencies of the vibration harmonics in the three directions before and after the closed-loop control. First, it can be clearly seen before the control that after considering the differences in the dissimilar physical quantities of the blades in Step 1, the 1 / rev frequency amplitude increases significantly, and the amplitudes of the 2 / rev and 3 / rev harmonics of non-integer multiples N / rev also appear. After applying the closed-loop control algorithm proposed in this paper, the 1 / rev harmonic amplitudes of the vibration velocities in the x, y, and z directions can be reduced by more than 80% significantly, which shows that the control strategy proposed in the present invention can correctly adjust the vibration imbalance characteristics generated by dissimilar blades.

[0118] As Figure 7It shows the variation trend of the maximum value of the rotor cone deviation during the closed-loop control process. From the curve before control, it can be seen that the dissimilarity of the blades results in a relatively large deviation value. During the control process, with the execution of online identification and closed-loop control, the identification results become gradually accurate, and the cone deviation value decreases rapidly. Eventually, it reaches near the set maximum value, which is about 54% lower than that before control.

[0119] Embodiment 2: Experimental verification of the closed-loop control of helicopter rotor vibration based on cone constraint

[0120] Similarly, perform the above-mentioned Steps 1 to 5 for experimental verification of the closed-loop control method of helicopter vibration based on cone constraint. Among them, there are two differences from Embodiment 1: The first difference is in Step 1, where the rotor aeroelastic coupling model derived from the 1 / rev harmonic of the hub vibration velocity and the rotor cone deviation value for calculating is replaced by the vibration sensor and cone sensor in the test system, that is, directly use the experimental data collected by the sensor for subsequent online identification and closed-loop control calculations; The second difference is that in order to combine the algorithm with the test hardware, the online identification algorithm and closed-loop control algorithm in Steps 2 to 5 need to be adapted into a real-time operable simulation software using the Simulink simulation platform and embedded into the test hardware for experimental verification, and other steps remain unchanged. As Figure 8 shown, it gives the software and hardware combined architecture of the test system. Thus, the experiment is carried out to obtain the result graph as Figures 9 to 12 shown.

[0121] Figures 9 to 11 It gives the change curves of the 1 / rev amplitude of the vibration velocity in the x, y, and z directions before, during, and after the closed-loop control. During the control process, the 1 / rev amplitudes in the three directions show a significant decrease. Eventually, after the control stabilizes, the control effect can reach more than 60% in all three directions, and the reduction amplitude in the y direction reaches 72%; Figure 12 It is the change curve of the maximum value of the rotor cone deviation before, during, and after the control. Although in the first half of the time after the control input stabilizes, the rotor may be affected by accidental gusts and there are large fluctuations, but the mean value in the second half also decreases significantly and gradually reaches near the constraint limit value, and the overall mean value also reaches a control effect of 61%. Therefore, the experimental results also prove the effectiveness of the present invention.

[0122] The above-mentioned embodiments only represent the implementation manners of the present invention, but should not be construed as limiting the scope of the patent of the present invention. It should be noted that for those skilled in the art, without departing from the concept of the present invention, several deformations and improvements can still be made, and these all belong to the protection scope of the present invention.

Claims

1. A closed-loop control method for helicopter vibration based on cone constraint, characterized in that It includes the following steps: Step 1: Set dissimilar physical quantities between the blades, establish a rotor aeroelastic coupling model, and thereby calculate the 1 / rev harmonic component of the vibration velocity at the rotor hub center and the rotor cone deviation value; Step 2: Take the change amount of the blade collective pitch angle as the control input and the 1 / rev harmonic component of the rotor vibration velocity as the control target, and establish a quadratic performance index function for the control problem; Step 3: Based on the relationship between the rotor cone deviation value and the change amount of the collective pitch angle, and considering the limitation of the change amount of the collective pitch angle at the same time, establish a constraint equation for the control input; Step 4: Use the recursive least squares identification algorithm to calculate the transfer matrices between the control target and the control input, and between the rotor cone deviation value and the control input, obtain a complete optimal control problem, and transform it into a standard quadratic programming problem; Step 5: Use the active set method to solve the optimal control problem, and judge whether the optimal control input obtained in the current closed-loop control period makes the 1 / rev harmonic of the vibration velocity and the rotor cone deviation value meet the error conditions; The specific content of Step 2 is as follows: Step 2.1, Based on the optimal control theory, establish a performance index function for the control problem, and set the performance index function as: where: J represents the value of the performance index function; z n is the control target, i.e., the 1 / rev harmonic component of the rotor vibration velocity in Equation (5); the subscript n represents the nth closed-loop control period; x n is the control input in the nth closed-loop control period, i.e., the change in the collective pitch angle in Equation (1); Δx n = x n - x n-1 is the change in the control input from the (n - 1)th closed-loop control period to the nth closed-loop control period; W z , W x , W Δx are the weight matrices corresponding to the control target z n , the control input x n , and the change in the control input Δx n respectively; Step 2.2, establish the relationship between the 1 / rev harmonic component z of the hub center vibration velocity and the change amount x of the collective pitch angle within two adjacent control cycles, that is, the relationship between the control target z and the control input x, specifically as follows: n and n the relationship between n and n the control input x, specifically as follows: z n = z n-1 + T z ·(x n - x n-1 ) In Equation (7), T z is the transfer matrix between the control target z n and the control input x n ; z n-1 , x n-1 represent the control target and the control input in the (n - 1)-th control period.

2. The helicopter vibration closed-loop control method based on cone constraint according to claim 1, characterized in that The specific content of Step 1 is as follows: Step 1.1, set dissimilar physical quantities between the blades, including the flap stiffness EI y , the lag stiffness EI z , the torsion stiffness GJ, the linear density dm, the lift coefficient c l0 , the drag coefficient c d0 , the pitch moment coefficient c m0 , these physical quantities will affect the dynamic response of the rotor, resulting in an increase in the 1 / rev harmonic component of the vibration velocity at the hub center; Step 1.2, Based on the multi-body dynamics theory, establish a rotor aeroelastic coupling model based on the dissimilar physical quantities of the blades, and solve the rotor cone deviation value; Considering the constraint equation Φ = 0, a set of differential-algebraic equations (DAEs) is obtained, and the linear term coefficients of the generalized acceleration and the generalized coordinate q are extracted to obtain the overall rotor aeroelastic coupling model, which is arranged in the standard form of differential-algebraic equations as follows: wherein, M, K, Q non , F a are respectively the mass matrix, stiffness matrix, generalized nonlinear term, and rotor aerodynamic force term of the rotor aeroelastic coupling model; θ0, x represent the collective pitch angle of the rotor and its variation, which directly affect the aerodynamic force term F a , where x is the control input; q, are respectively the generalized coordinates and generalized accelerations; the superscript T represents the matrix transpose symbol; Φ = 0 represents the constraint equation 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 q of each blade, and then the blade tip height value β of each blade is calculated. Assuming that the number of rotor blades is N b , arbitrarily select the e-th blade as the reference blade, and calculate the rotor cone deviation value, that is, the difference between the blade tip height of other blades and the reference blade, as follows: In the formula, Δβ represents the rotor cone deviation value vector; represents the tip height values of the 1st, 2nd, …, N b blade tips; β e represents the tip height value of the e-th blade; Step 1.3, Solve the vibration load at the rotor hub center, which is composed of the blade root loads superimposed, specifically as follows: In the formula, is the azimuth angle of the m-th blade; N b is the number of blades; the superscript H represents the hub; the superscript m represents the m-th blade; the subscripts x, y, z represent the forces or torques in the x, y, z coordinate axis directions; respectively represent the root shear forces of the m-th blade in the x, y, z directions; respectively represent the root torques of the m-th blade in the x, y, z directions; respectively represent the vibration shear forces of the hub center in the x, y, z directions; respectively represent the vibration torques of the hub center in the x, y, z directions; Step 1.4, Convert the vibration load at the rotor hub center into the vibration velocity at the rotor hub center, and use Fourier transform to extract the 1 / rev harmonic component of the vibration velocity; In engineering, vibration is measured using vibration sensors, and the vibration velocities in the x, y, and z directions at the rotor hub center are used to characterize the magnitude of the vibration; Therefore, combining the relationship between the vibration load and the vibration velocity, convert the hub vibration load into the vibration velocity, as follows: where a x , a y , a z represent the vibration velocities in the x, y, and z directions of the hub center, respectively; the superscript T represents the transpose symbol; T a is a 3×6 dimensional transformation matrix; For the calculated vibration velocity a x , a y , a z Perform frequency domain analysis using Fourier transform to extract the 1 / rev harmonic component, that is: a x The 1 / rev harmonic component of includes the sine component a x1s and the cosine component a x1c ; a y The 1 / rev harmonic component of includes the sine component a y1s and the cosine component a y1c ; a z The 1 / rev harmonic component of includes the sine component a z1s and the cosine component a z1c ; Thus, the control objective is obtained: z = [a x1s , a x1c , a y1s , a y1c , a z1s , a z1c T In Equation (5), z represents the 1 / rev harmonic component of the vibration velocity, which is also the control target; the superscript T represents the transpose symbol.​ 3. A closed-loop control method for helicopter vibration based on cone constraint according to claim 1, characterized in that, The specific content of Step 3 is as follows: Establish the relationship between the rotor cone deviation value and the change amount of the collective pitch angle, as follows: Δβ n = T β ·x n + Δβ0 (8) where: Δβ n represents the rotor cone deviation value in the nth closed-loop control cycle, i.e., Equation (2); T β is the transfer matrix between the rotor cone deviation value and the control input x n ; Δβ0 represents the rotor cone deviation value when the control input x = 0; Rotor cone deviation value Δβ n It shall not exceed the set value β0, so there is: |Δβ n |≤β0 (9) In the formula, β0 is the maximum value of the set rotor cone deviation value; Substitute Equation (8) into Equation (9) to obtain a set of inequality constraints on the control input x n as follows: A·x n ≤b (10) In the formula, the inequality coefficient matrix The right-hand side term of the inequality is where is the transpose of the column vector composed of β0.

4. A closed-loop control method for helicopter vibration based on cone constraint according to claim 3, characterized in that In the said step 3, the control input x n is the change amount of the total rotor pitch angle. When the rotor is operating normally, the change of the control input cannot affect the normal rotation of the rotor, and there are the following restrictions: |x n | ≤ x0 (11) In the formula, x0 is the maximum allowable change amount of the collective pitch angle of a single blade.

5. A closed-loop control method for helicopter vibration based on cone constraint according to claim 4, characterized in that, The specific content of Step 4 is as follows: Using the recursive least squares identification method, calculate the transfer matrices T in equations (7) and (8). z and T β ; First, determine the parameters to be identified, i.e., the transfer matrices T in Equations (7) and (8). z and T β ; the calculation processes of both are exactly the same. The transfer matrix T in Equation (7) is used for illustration. z Suppose the dimension of matrix T z is i×j. Connect its columns end to end to obtain a column vector with (s = i·j) parameters, which is the parameter to be identified, denoted as α n . Then, the discrete observation equation for the (n + 1)-th closed-loop control period is written as: y n+1 = g n+1 ·α n+1 + v n+1 (12) where y n+1 is the observation vector of the (n + 1)-th closed-loop control period; g n+1 is the observation matrix of the (n + 1)-th closed-loop control period; v n+1 is the measurement noise of the (n + 1)-th closed-loop control period; α n+1 is the parameter to be identified in the (n + 1)-th closed-loop control period; Then the total observation equation for the first n closed-loop control periods is: Y n = G n ·α n + V n (13) where Y n represents the total observation vector in the previous n closed-loop control cycles; G n represents the total observation matrix in the previous n closed-loop control cycles; V n represents the total measurement noise in the previous n closed-loop control cycles According to the least squares estimation criterion, considering the total observation vector for the first n + 1 closed-loop control periods, the parameter to be identified is written as: wherein, is the parameter to be identified estimated by the identification algorithm; G n+1 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 transpose of the matrix; In formula (14), which is called the information matrix and expands to: where g n is the observation matrix for the nth closed-loop control period; According to the matrix inversion formula, the parameter to be identified for the (n + 1)-th closed-loop control period is written as: In the formula, is the parameter to be identified in the (n + 1)-th closed-loop control period, that is, the latest identification result; is the parameter to be identified in the n-th closed-loop control period, that is, the identification result of the previous control period; K n is the gain matrix in the n-th closed-loop control period; After obtaining the transfer matrices T z and T β and combining equations (6), (7), (10), and (11), the control problem in the nth closed-loop control period is transformed into a standard quadratic programming problem, which is organized as follows: where f(x n ) is the function value of the quadratic performance index function; Q is the quadratic coefficient matrix of x n ; c is the linear coefficient vector of x n .

6. The helicopter vibration closed-loop control method based on cone constraint according to claim 5, characterized in that, The specific content of Step 5 is as follows: Step 5.1, First, use the active 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 period makes the vibration velocity 1 / rev harmonic satisfy the error condition. Therefore, substitute the current control input x n into Equation (1) to calculate the rotor cone deviation value Δβ, and then use Fourier transform to calculate the 1 / rev harmonic amplitude a x1a , a y1a , a z1a , and determine whether the following equation holds: Wherein, Δβ is the rotor cone deviation vector; β0 is the maximum value of the set rotor cone deviation; a x1a , a y1a , a z1a are respectively the 1 / rev harmonic amplitudes of the vibration velocities in the x, y, and z directions of the hub center; a x0 , a y0 , a z0 are respectively the maximum values of the 1 / rev harmonic amplitudes of the vibration velocities in the x, y, and z directions of the set hub center; If the optimal control quantity in the current closed-loop control period can make Equation (18) hold, then the current optimal control quantity output is maintained, and whether Equation (18) holds is monitored in real time; if Equation (18) does not hold, the optimal control quantity x of the current closed-loop control period n is used as the control input x of the (n + 1)-th closed-loop control period n+1 initial value, substitute it into Equation (1), recalculate the data of the control objective and control constraints, form the performance index function in Step 2, then calculate the constraint equation in Step 3, execute the recursive least squares identification algorithm in Step 4, and finally use the active set method in Step 5 to solve the standard quadratic programming problem to obtain the optimal control input x of the (n + 1)-th closed-loop control period n+1 , substitute it into Equation (1) and judge again whether Equation (18) holds, and so on in a loop.

Images