Mechanical vibration micro-grid electromechanical coupling system stability analysis method

Abstract

The invention relates to the technical field of power grid stability analysis, in particular to a mechanical vibration micro-grid electromechanical coupling system stability analysis method, which comprises the following steps of: detecting a tooth surface contact radius, a relative speed and a Young modulus to generate a contact stiffness initial value, correcting nonlinear stiffness by interval difference, and calculating the stability of a micro-grid electromechanical coupling system. The method comprises the following steps: calculating a coupling response vector by combining rotor displacement and bearing clearance, inputting a Lyapunov exponent method to solve an exponent vector, constructing a model length ratio to detect derivative characteristics, extracting a multidimensional state vector in a positive value period and the coupling response vector to generate a compensation action, and constructing an FPGA simulation path injection disturbance execution error feedback back-transmission fuzzy control algorithm. According to the method, a contact stiffness dynamic correction and nonlinear difference method is introduced, six-degree-of-freedom displacement and bearing clearance parameters are combined, an electromechanical coupling model is established, a stability boundary is recognized through Lyapunov regression, a critical state is detected through derivative features, and multi-dimensional state vector fuzzy control compensation is fused. And the response and the immunity are optimized by combining FPGA simulation and a disturbance mechanism.

Description

Technical Field

[0001] The present invention relates to the technical field of power grid stability analysis, and in particular to a stability analysis method for a mechanically vibrating microgrid electromechanical coupling system. Background Art

[0002] The field of grid stability analysis involves modeling and evaluating the dynamic characteristics of voltage, current, frequency, and power in power systems to ensure that the system can maintain its operating state or quickly recover to a stable operating condition under various disturbances. Its core issues include transient stability analysis of power systems, small disturbance stability assessment, electromagnetic transient simulation, frequency response capability determination, and coordinated control research under multi-source grid-connected conditions. The technical field generally covers the hierarchical response relationship from large power grids to regional power grids and microgrids. In particular, it systematically studies the modeling methods and stability assessment methods for the dynamic feedback between energy conversion devices and physical coupling mechanisms in the context of high renewable energy penetration. Traditional microgrid electromechanical coupling system stability analysis methods address the dynamic response problems caused by the coupling characteristics of the mechanical structure and electrical system in distributed generation systems. These methods typically use a set of mathematical equations based on equivalent circuit parameters to linearly approximate the motor-side state variables. Eigenvalue analysis is performed under steady-state conditions by mapping the input torque, output current, and shaft speed. Frequency sweeps are then used to identify the stability boundary of the system at a given operating point. This method relies on the analytical construction of fixed-frequency gain maps and the characteristic spectrum distribution of the state derivative matrix.

[0003] Traditional methods are based on linear approximation models of equivalent circuit parameters, which make it difficult to accurately characterize the coupling relationship between the dynamic deformation of mechanical structures and electrical parameters. When the shaft system displacement suddenly changes or the bearing clearance changes, the fixed frequency domain gain diagram cannot reflect the nonlinear evolution law of the contact stiffness, resulting in the deviation of the eigenvalue analysis results from the actual working conditions. The frequency scanning method under steady state is adopted, and there is a lack of real-time tracking mechanism for the dynamic characteristics of the impedance response. When the system is in a critical stable state, the lack of first-order derivative zero point detection leads to a lag in the stability boundary judgment. The state derivative matrix analysis that relies on fixed operating points cannot adapt to transient parameter fluctuations under multi-source disturbances, resulting in a decrease in the matching degree between the control strategy and the real-time working conditions, affecting the dynamic adjustment accuracy of the electromechanical coupling system. Summary of the Invention

[0004] The purpose of the present invention is to solve the shortcomings of the prior art and to propose a stability analysis method for a microgrid electromechanical coupling system with mechanical vibration.

[0005] To achieve the above objectives, the present invention adopts the following technical solution: a method for analyzing the stability of a microgrid electromechanical coupling system with mechanical vibration, comprising the following steps:

[0006] S1: Detect the tooth surface contact radius, relative speed and Young's modulus, generate the initial value of contact stiffness through product operation, use the interval combined difference method of forward difference, center difference and backward difference to correct the nonlinear change of stiffness, combine the rotor six-degree-of-freedom displacement data and bearing clearance, and use the finite difference method to output the contact stiffness correction value and coupling response vector;

[0007] S2: Input the coupling response vector into the Lyapunov index method to perform trajectory sequence regression, calculate the Lyapunov index vector, construct the Z1 / Z2 mode length ratio based on the impedance response, and perform first-order derivative zero point detection and second-order derivative sign determination;

[0008] S3: According to the positive period of the Lyapunov index vector, the vibration acceleration, electromagnetic torque and winding temperature are collected, combined into a multi-dimensional state vector, and the multi-dimensional state vector and the coupling response vector are input into the fuzzy control algorithm to generate a compensation action vector;

[0009] S4: Construct an FPGA simulation path based on the compensation action vector, perform a 10μs step state synchronization update, inject geometric dimension deviation parameters to generate an error disturbance version, output an error feedback table through a speed / load combination test and return the fuzzy control algorithm.

[0010] As a further solution of the present invention, the interval combined difference method calculates the stiffness change gradient through forward difference, verifies the gradient direction through central difference, and corrects the nonlinear error through backward difference. The synergistic effect of the three reduces the stiffness correction error by 12%-15%.

[0011] The Z1 / Z2 module length ratio is a dynamic coupling factor between the real module value and the imaginary module value of the impedance;

[0012] The positive time period determination threshold is that the exponential component is greater than 0.5 for three consecutive sampling periods.

[0013] As a further solution of the present invention, the step S1 is specifically as follows:

[0014] S101: Detect the tooth surface contact radius, relative speed, and Young's modulus, establish a geometric relationship model between the contact radius and relative speed, construct a Young's modulus compensation function, substitute the three sets of parameters into the stiffness calculation formula, perform multiplication operation, and generate the initial value of the contact stiffness;

[0015] S102: Based on the initial value of the contact stiffness, a differential iterative calculation model is established, the stiffness change gradient is calculated by forward difference, the gradient direction is verified by central difference, and the nonlinear error is corrected by backward difference to obtain a stiffness correction coefficient;

[0016] S103: Call the rotor six-degree-of-freedom displacement data and bearing clearance parameters, construct a displacement-stiffness coupling equation, substitute the stiffness correction coefficient into the finite difference method calculation framework, solve the stiffness correction term through the Jacobian matrix, perform matrix iterative operation in combination with the displacement response coupling term, and synchronously output the contact stiffness correction value and the coupling response vector.

[0017] As a further solution of the present invention, the stiffness calculation formula is:

[0018] Where E represents the Young's modulus of the gear material in GPa, R c Represents the instantaneous curvature radius of the tooth surface contact point, unit: mm, v r Represents the tangential relative speed of the gear meshing point, unit m / s, denominator is the coupling correction term between contact geometry and kinematic state.

[0019] As a further solution of the present invention, the step S2 is specifically as follows:

[0020] S201: Based on the coupled response vector, the trajectory sequence is regressed in the time domain using the Lyapunov exponent method. A sliding window mechanism is introduced when calculating the exponential components of multiple dimensions. The exponential decay curve is fitted using the least squares method to generate a Lyapunov exponent vector.

[0021] S202: Call the Lyapunov exponent vector, extract the real and imaginary moduli of the impedance response, and construct the Z1 / Z2 modulus ratio using the formula:

[0022]

[0023] Obtaining dynamic impedance ratios through calculation and generating a dynamic impedance ratio sequence;

[0024] Among them, R represents the dynamic impedance ratio, Z1 represents the real part modulus of the bearing support structure impedance, unit N·s / m, Z2 represents the imaginary part modulus of the rotor system impedance, unit N·s / m, δ represents the dynamic coupling factor of the bearing clearance, θ represents the impedance phase difference, unit radian, ε represents the ambient temperature interference coefficient, range 0.1-0.3, the denominator is the joint correction term for environmental disturbance and dynamic coupling;

[0025] S203: Based on the dynamic impedance ratio sequence, a five-point central difference method is used to calculate the first-order derivative, the zero point position is detected by sign change, and a cubic spline interpolation is applied to obtain the sign distribution of the second-order derivative to generate a guide change feature identification group.

[0026] As a further solution of the present invention, the sliding window mechanism sets the window length to 2.5 times the system base frequency period.

[0027] As a further solution of the present invention, the step S3 is specifically as follows:

[0028] S301: establishing a time synchronization acquisition window for vibration acceleration, electromagnetic torque, and winding temperature according to the positive value period of the Lyapunov exponent vector, processing the original signal using a sliding average filter, performing a standard deviation normalization operation on the three sets of parameters, and combining them into a multidimensional state vector;

[0029] S302: Calling the coupling response vector, constructing a three-dimensional input domain of acceleration-torque-temperature, designing a triangular membership function to cover multiple parameter variation ranges, establishing a nine-rule fuzzy rule base based on an IF-THEN structure, performing fuzzy matching between the multidimensional state vector and the coupling response vector, and generating a fuzzy decision matrix;

[0030] S303: Based on the fuzzy decision matrix, a centroid method is used to perform defuzzification operations, a weighted average of the membership of multiple control variables is calculated, a compensation action quantization model is established, and a compensation action vector is generated through a multiplication and accumulation operation of the membership and the domain value.

[0031] As a further solution of the present invention, the triggering mechanism of the time synchronization acquisition window is that any two dimensional components in the Lyapunov exponent vector simultaneously exceed 0.5;

[0032] The vertex position of the triangle membership function is dynamically adjusted according to the bearing type. The vertex offset of rolling bearings is set to ±15%, and that of sliding bearings is set to ±8%.

[0033] As a further solution of the present invention, the step S4 is specifically as follows:

[0034] S401: Based on the compensation action vector, construct an FPGA hardware description language model, design a state synchronization update mechanism, set timing constraints, and generate an FPGA simulation path;

[0035] S402: calling the FPGA simulation path, establishing a mapping relationship between geometric dimension deviation parameters and the simulation path, generating a deviation parameter combination using a Latin hypercube sampling method, injecting the parameter combination into the simulation path for disturbance propagation calculation, and outputting an error disturbance version;

[0036] S403: Execute a speed / load combination test, input the error disturbance version into the test platform, collect torque fluctuation and vibration spectrum data, construct an error index evaluation matrix, calculate the comprehensive error value of multiple test points through weighted summation, generate the error feedback table and send it back to the fuzzy control algorithm.

[0037] As a further solution of the present invention, the timing constraint condition sets the clock period to 10ns and the setup time margin to ≥2ns;

[0038] The Latin hypercube sampling method adopts a dimensional orthogonalization strategy to improve the sampling efficiency by 40%-60% compared with the Monte Carlo method;

[0039] The weight distribution of the weighted summation is as follows: the torque fluctuation weight is 0.6, the vibration spectrum weight is 0.4, and the comprehensive error value determination threshold is 0.75.

[0040] Compared with the prior art, the advantages and positive effects of the present invention are:

[0041] In the present invention, by introducing the dynamic correction mechanism of contact stiffness and the nonlinear interval difference method, the transient stiffness change characteristics of the mechanical coupling system are effectively captured. By combining the six-degree-of-freedom displacement data and the bearing clearance parameters, a more accurate electromechanical coupling dynamic model is established. The Lyapunov index vector trajectory regression analysis is adopted to realize the real-time identification of the multi-dimensional stability boundary. The derivative feature detection based on the impedance response modulus ratio is enhanced to enhance the prediction ability of the critical state of the system. The vibration acceleration, electromagnetic torque and temperature parameters are integrated to construct a multi-dimensional state vector. The dynamic compensation strategy is generated by the fuzzy control algorithm. The FPGA high-speed simulation path and the error disturbance injection mechanism are combined to form a closed-loop parameter optimization system, which significantly improves the dynamic response accuracy and anti-interference ability of the system under nonlinear working conditions, and simultaneously realizes the stability control and adaptive adjustment of the electromechanical coupling system under complex load fluctuations. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] Figure 1 It is a schematic diagram of the main steps of the present invention;

[0043] Figure 2 is a flow chart of the steps of S1 of the present invention;

[0044] Figure 3 This is a flow chart of the steps of S2 of the present invention;

[0045] Figure 4 This is a flow chart of the steps of S3 of the present invention;

[0046] Figure 5 This is a flow chart of the steps of S4 of the present invention. DETAILED DESCRIPTION

[0047] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention.

[0048] In the description of the present invention, it should be understood that the terms "length," "width," "up," "down," "front," "back," "left," "right," "vertical," "horizontal," "top," "bottom," "inside," "outside," and the like, indicating positions or relationships, are based on the positions or relationships shown in the accompanying drawings and are intended only to facilitate the description of the present invention and simplify the description. They do not indicate or imply that the devices or elements referred to must have a specific orientation, be constructed, or operate in a specific orientation. Therefore, they should not be construed as limiting the present invention. Furthermore, in the description of the present invention, "plurality" means two or more, unless otherwise expressly and specifically defined.

[0049] Example 1

[0050] See also Figure 1 The present invention provides a technical solution: a method for analyzing the stability of a microgrid electromechanical coupling system with mechanical vibration, comprising the following steps:

[0051] S1: Detect the tooth surface contact radius, relative speed and Young's modulus, generate the initial value of contact stiffness through product operation, use the interval combined difference method of forward difference, center difference and backward difference to correct the nonlinear change of stiffness, combine the rotor six-degree-of-freedom displacement data and bearing clearance, and use the finite difference method to output the contact stiffness correction value and coupling response vector;

[0052] S2: Input the coupled response vector into the Lyapunov index method to perform trajectory sequence regression, calculate the Lyapunov index vector, construct the Z1 / Z2 mode length ratio based on the impedance response, and perform first-order derivative zero point detection and second-order derivative sign determination;

[0053] S3: According to the positive period of the Lyapunov index vector, the vibration acceleration, electromagnetic torque and winding temperature are collected and combined into a multi-dimensional state vector, which is then input into the fuzzy control algorithm together with the coupling response vector to generate a compensation action vector;

[0054] S4: Build an FPGA simulation path based on the compensation action vector, perform a 10μs step-size state synchronization update, inject geometric dimension deviation parameters to generate an error perturbation version, output the error feedback table through the speed / load combination test and feed it back to the fuzzy control algorithm.

[0055] The interval combined difference method calculates the stiffness change gradient through forward difference, verifies the gradient direction through central difference, and corrects the nonlinear error through backward difference. The synergistic effect of the three reduces the stiffness correction error by 12%-15%.

[0056] The Z1 / Z2 mode length ratio is the dynamic coupling factor between the real and imaginary impedance modes.

[0057] The threshold for determining the positive period is that the exponential component is greater than 0.5 for three consecutive sampling periods.

[0058] See also Figure 2 , the steps of S1 are as follows:

[0059] S101: Detect the tooth surface contact radius, relative speed, and Young's modulus, establish a geometric relationship model between the contact radius and relative speed, construct a Young's modulus compensation function, substitute the three sets of parameters into the stiffness calculation formula, perform multiplication operation, and generate the initial value of the contact stiffness;

[0060] After detecting the tooth surface contact radius, relative speed and Young's modulus, the laser displacement sensor array deployed inside the gearbox is aligned with the gear meshing line. The three-dimensional coordinates of the meshing point are continuously collected at a sampling frequency of 10kHz and substituted into the curvature calculation model of Hertz contact theory to obtain the instantaneous curvature radius R of the tooth surface contact point at time t = 0.1s. c At the same time, a dual-channel laser Doppler vibrometer is used to synchronously measure the tangential linear velocity of the driving wheel and the driven wheel at the meshing point, and the driving wheel speed v1 = 22.5 m / s and the driven wheel speed v2 = 10.5 m / s are obtained respectively. The tangential relative velocity v is obtained by performing a difference operation on the two. r =|v1-v2|=12m / s. In addition, from the design database of the gearbox of this type of wind turbine generator (model HZD-2.5MW), the Young's modulus E of the gear material (grade 42CrMo steel) is retrieved as 210GPa. The geometric relationship model established is to put R c With v r The instantaneous measurement value is mapped to a two-dimensional state plane in real time, and the Young's modulus compensation function is constructed based on the internal oil temperature of the gearbox fed back by the temperature sensor installed in the bearing seat. When the oil temperature deviates from the standard operating temperature by 60°C and increases by 10°C, a compensation coefficient of -0.2% is applied to the Young's modulus E. Assuming that the current oil temperature is 60°C and the compensation coefficient is 0, these three sets of parameters are then substituted into the stiffness calculation formula Perform product operation. In this formula, K0 is the initial value of contact stiffness, in N / m, and E represents the Young's modulus of the gear material. The actual measured value is 210 GPa, which is converted to 2.1×10 11 Pa, R c Represents the instantaneous curvature radius of the tooth surface contact point, its measured value is 20mm, and the unit is uniformly converted to 0.02m, v r Represents the tangential relative speed of the gear meshing point, and its measured value is 12m / s. The denominator The operation logic is to square the curvature radius representing the static geometric characteristics and the relative velocity representing the dynamic kinematic characteristics and take the square root of the sum to obtain a comprehensive geometric-kinematic state coupling correction factor. This factor performs a nonlinear correction on the product term representing the linear relationship in the numerator. Substituting the obtained and unified numerical values ​​into the formula, the operation of the numerator is (2.1×10 11 Pa)×(0.02m)×(12m / s)=5.04×10 10 N·m / s, the calculation of the denominator is Finally, by dividing the numerator and denominator, the initial value of contact stiffness K0≈4.2×10 9 N / m.

[0061] S102: Based on the initial value of the contact stiffness, a differential iterative calculation model is established. The stiffness change gradient is calculated through forward difference, the gradient direction is verified using central difference, and the nonlinear error is corrected in combination with backward difference to obtain the stiffness correction coefficient.

[0062] Based on the obtained initial value of contact stiffness K0 = 4.2 × 10 9 N / m, a differential iterative calculation model is established. The model runs in a calculation domain with a preset time step of 100. The time step Δt is set to 0.001s according to the system sampling frequency. First, the stiffness change gradient is calculated by forward difference. Specifically, the stiffness value K(t) = 4.2×10 9 N / m, and the predicted stiffness value K(t+Δt)=4.203×10 9 N / m, perform gradient calculation: Secondly, the central difference is used to verify the gradient direction, and the stiffness history value K(t-Δt) = 4.198×10 9 N / m, perform central difference gradient calculation: Since the forward differential gradient is 3×10 9 N / (m·s) and central differential gradient 2.5×10 9 The signs of N / (m·s) are all positive, indicating that the stiffness is monotonically increasing in this range and the gradient direction is valid. Subsequently, the nonlinear error is corrected by combining the backward difference, and the backward difference gradient is calculated as:

[0063] The three gradient values ​​are weighted averaged, where the weight coefficient is set based on a large amount of experimental data analysis. For the working condition where the system changes relatively slowly, the central difference weight is the highest. The setting of this coefficient is based on the comparison of the root mean square error between the correction results and the actual values ​​under different weight combinations in 100 sets of simulation tests, and the set of weights with the smallest error is selected, namely the forward weight w f =0.3, center weight w c =0.4, backward weight w b = 0.3, the calculated comprehensive gradient is 0.3×(3×10 9 )+0.4×(2.5×10 9 )+0.3×(2×10 9 )=2.5×10 9 N / (m·s), and multiplying the comprehensive gradient by the time step Δt yields the stiffness correction ΔK = 2.5 × 10 9 × 0.001 = 2.5 × 10 6 N / m, the ratio of this correction value to the initial value of contact stiffness is the stiffness correction coefficient, and its value is The final stiffness correction coefficient is 0.000595.

[0064] S103: Call the rotor six-degree-of-freedom displacement data and bearing clearance parameters to construct the displacement-stiffness coupling equation, substitute the stiffness correction coefficient into the finite difference method calculation framework, solve the stiffness correction term through the Jacobian matrix, perform matrix iterative operation combined with the displacement response coupling term, and synchronously output the contact stiffness correction value and coupling response vector;

[0065] First, the rotor six-degree-of-freedom displacement data and bearing clearance parameters are called. The six-degree-of-freedom displacement data is acquired in real time by three pairs of orthogonal eddy current sensors evenly distributed on the bearing seats at both ends of the rotor. In the current sampling period, the displacement vector collected and integrated is X=

[0066] 0.01mm,0.015mm,0.005mm,0.001rad,0.0012rad,0.0008rad] T , these six components correspond to the translational displacements in the x, y, and z directions and the rotational angular displacements around these three axes, respectively. The bearing clearance parameters are retrieved from the database of the equipment digital twin model according to the bearing model used by the rotor (NN3020K). Its radial clearance is 15 μm. The constructed displacement-stiffness coupling equation is a 7×7 matrix equation. The variables of this equation include the displacements of six degrees of freedom and a stiffness term to be calculated. The stiffness correction coefficient of 0.000595 obtained in step S102 is substituted into the finite difference method calculation framework. During the specific execution, this coefficient is combined with the initial value of the contact stiffness 4.2×10 9 Multiply N / m to obtain the stiffness correction term ΔK = 2.5 × 106 N / m, and use it as the initial perturbation value of the seventh element on the diagonal of the Jacobian matrix (corresponding to the stiffness variable). The other elements of the Jacobian matrix are composed of the partial derivatives of the force on each degree of freedom with respect to the displacement of each degree of freedom in the system dynamics equation, for example The Jacobian matrix is ​​solved by the Gauss-Seidel iterative method, and the convergence criterion is set as the norm of the difference between the two iterative results is less than 10 -6 After 5 iterations, the matrix solution converges and the complete stiffness correction matrix is ​​obtained. The matrix iteration operation is combined with the displacement response coupling term. The displacement response coupling term here refers to the product of the six-degree-of-freedom displacement vector X and the corrected stiffness matrix. This operation synchronously outputs the contact stiffness correction value and the coupling response vector. The contact stiffness correction value is the sum of the initial stiffness value and the correction term accumulated in all iterative steps, that is, 4.2×10 9 +2.5×10 6 +…(correction value of subsequent iterations) = 4.20251×10 9 N / m, the coupled response vector is a force and torque vector F containing six components, and its calculated value is F = [42025N, 63037N, 21012N, 42.0Nm, 50.4Nm, 33.6Nm] T .

[0067] The stiffness calculation formula is

[0068] Where E represents the Young's modulus of the gear material in GPa, R c Represents the instantaneous curvature radius of the tooth surface contact point, unit: mm, v r Represents the tangential relative speed of the gear meshing point, unit m / s, denominator is the coupling correction term between contact geometry and kinematic state.

[0069] See also Figure 3 , the steps of S2 are as follows:

[0070] S201: Based on the coupled response vector, the Lyapunov exponent method is used to perform time-domain regression on the trajectory sequence. A sliding window mechanism is introduced when calculating the exponential components of multiple dimensions. The exponential decay curve is fitted using the least squares method to generate the Lyapunov exponent vector.

[0071] The coupled response vector F based on the S103 output is [42025, 63037, 21012, 42.0, 50.4, 33.6] TThe multidimensional trajectory sequence consisting of a total of 4096 sampling points is collected continuously and subsequently. The Lyapunov index method is used to perform time domain regression analysis on the attractor trajectory reconstructed in the phase space of the multidimensional time series. A sliding window mechanism is introduced when calculating the multi-dimensional exponential components. The window length of this mechanism is set according to the system fundamental frequency. By performing fast Fourier transform (FFT) analysis on the electromagnetic torque signal, the system fundamental frequency is identified as the gear meshing frequency of 400Hz. Based on this, the window time length T is calculated. w = 2.5 / 400Hz = 0.00625s, the corresponding number of sampling points is 0.00625s × 10kHz = 62.5, rounded up to 63 data points, set the window to slide over a sequence of 4096 points with a step size of 10 data points each time, and apply the Wolf algorithm to calculate the maximum Lyapunov exponent for each data segment within the window. When executing it, select an initial point in the reconstructed phase space and find a point with a minimum distance to the initial point (for example, a Euclidean distance of 10). -6 ) and track the two points in the window time T w The evolution of the window is calculated, the logarithm of their separation distance is calculated, and the logarithm of the separation distance is divided by the evolution time. The separation rate calculated for all windows is averaged. At the same time, the least square method is used to calculate the separation rate of each response component (such as F x ,F y )’s trajectory divergence trend is fitted as A·e λt The exponential curve is an exponential curve, where the fitted exponent λ is the Lyapunov exponent component of this dimension. After calculating and fitting the six dimensional components separately, the Lyapunov exponent vector is finally generated. For example, at a certain moment when the system tends to be unstable, the calculated vector is L = [0.62, 0.58, -0.15, -0.21, -0.85, -1.50].

[0072] S202: Call the Lyapunov exponent vector to extract the real and imaginary moduli of the impedance response and construct the Z1 / Z2 modulus ratio using the formula:

[0073]

[0074] Obtaining dynamic impedance ratios through calculation and generating a dynamic impedance ratio sequence;

[0075] Among them, R represents the dynamic impedance ratio, Z1 represents the real part modulus of the bearing support structure impedance, unit N·s / m, Z2 represents the imaginary part modulus of the rotor system impedance, unit N·s / m, δ represents the dynamic coupling factor of the bearing clearance, θ represents the impedance phase difference, unit radian, ε represents the ambient temperature interference coefficient, range 0.1-0.3, the denominator is the joint correction term for environmental disturbance and dynamic coupling;

[0076] The Lyapunov index vector L generated by S201 is called, and the real part modulus of the bearing support structure impedance and the imaginary part modulus of the rotor system impedance synchronized with the timestamp of the vector are extracted from the electromechanical coupling system impedance online monitoring module. The dynamic excitation force and response speed of the support structure are measured in real time by the force sensor and acceleration sensor installed on the bearing seat. The real part modulus of the impedance |Z1| is calculated to be 8.5×10 4 N·s / m. By analyzing the cross-correlation spectrum between the rotor vibration displacement sensor signal and the motor drive current signal, the phase information and amplitude information are extracted, and the impedance imaginary modulus |Z2| of the rotor system is calculated to be 5.2×10 4 N·s / m, then construct the Z1 / Z2 modulus ratio and substitute it into the formula The dynamic impedance ratio is obtained by calculation. The benefit of this formula is that by introducing a joint correction term including the dynamic coupling factor δ of the bearing clearance, the impedance phase difference θ and the ambient temperature interference coefficient ε, the calculated dynamic impedance ratio R not only reflects the inherent impedance matching characteristics of the system, but also integrates the dynamic influence of nonlinear operating conditions and external environmental changes. In the formula, R represents the dimensionless dynamic impedance ratio, |Z1| is the real part modulus of the bearing support structure impedance, |Z2| is the imaginary part modulus of the rotor system impedance, and δ is the dynamic coupling factor of the bearing clearance, which is calculated based on the real-time radial force F of the bearing. r and bearing dynamic load rating C r The quantification formula is: Measure F under current working conditions r =80kN, refer to the manual to get C r =150kN, δ≈0.8 is calculated, θ is the impedance phase difference, obtained by comparing the phase angles of Z1 and Z2 (15° and 75° respectively), |θ|=|15°-75°|=60°, and converted to radians, that is ε is the ambient temperature interference coefficient, and its setting reference is: at the standard ambient temperature of 25℃, ε takes the base value of 0.1, and the ε value increases by 0.05 for every 5℃ deviation of the ambient temperature from the standard. This rule is verified by experiments. The system is placed in a temperature control box and operated at 15℃, 25℃, 35℃, and 45℃ respectively. The random fluctuation component of the system response is recorded, and its normalized variance is calculated. It is found that the variance is approximately linearly related to the temperature deviation value, and the proportional coefficient is about 0.01 / ℃. Based on this, if the current ambient temperature sensor reading is 35℃, then Substitute the above parameter assignments into the formula for calculation: This calculation is performed continuously at a frequency of 100 Hz, generating a dynamic impedance ratio series R(t). The result 1.658 indicates the dynamic impedance ratio state of the system at the current moment. A larger value generally means that the real damping has a more significant effect than the imaginary inertia / elastic force.

[0077] S203: Based on the dynamic impedance ratio sequence, the first-order derivative is calculated using the five-point central difference method, the zero point position is detected by the sign change, and the second-order derivative sign distribution is obtained by applying cubic spline interpolation to generate a guide change feature identification group;

[0078] Based on the dynamic impedance ratio sequence R(t) generated by S202, for example, a segment containing 10 consecutive data points is intercepted: [1.650, 1.658, 1.663, 1.665, 1.664, 1.660, 1.655, 1.651, 1.648, 1.646]. The first-order derivative of each data point in the sequence is calculated using the five-point central difference method. The fourth point in the sequence, R4=

[0079] For example, 1.665, the calculation of its first-order derivative uses the values ​​of the two points before and after it, according to the formula Where h is the sampling time interval of 0.01s, substitute the value for calculation: By performing this calculation for each point in the sequence, a first-order derivative sequence is obtained, and the zero point position is determined by detecting whether the sign of the derivative sequence value changes from positive to negative or from negative to positive. For example, if the calculation finds that R′(t4)>0 and R′(t5)<0, then the zero point is determined to be between the two sampling points t4 and t5. Subsequently, cubic spline interpolation is applied to obtain the sign distribution of the second-order derivative. Specifically, a piecewise cubic polynomial is constructed at four points near the above zero point (for example, t3, t4, t5, t6). The polynomial has the characteristic of continuous second-order derivatives at the connection points, so that it passes through these data points smoothly. Then, the second-order derivative of the interpolation polynomial at the zero point is calculated and its sign is determined. If the calculated value of the second-order derivative is negative, it indicates that the zero point corresponds to a local maximum point. If it is positive, it is a local minimum point. The position, type (maximum / minimum), and corresponding first-order derivative sign change and second-order derivative sign of each identified zero point are combined into a feature identifier, for example [t=0.045s, ZeroCrossing, Max, R′from(+)to(-), R″<0]. Finally, all identified feature identifiers are summarized to generate a derivative change feature identifier group.

[0080] The sliding window mechanism sets the window length to 2.5 times the system base frequency period.

[0081] See also Figure 4 , the specific steps of S3 are:

[0082] S301: Based on the positive period of the Lyapunov exponent vector, a time synchronization acquisition window for vibration acceleration, electromagnetic torque, and winding temperature is established. The original signal is processed using a sliding average filter, and the three sets of parameters are normalized by standard deviation and combined into a multidimensional state vector.

[0083] According to the positive component in the Lyapunov exponent vector L = [0.62, 0.58, -0.15, -0.21, -0.85, -1.50], a positive period judgment is performed. The threshold of this judgment is that any two-dimensional components in the exponent vector are greater than 0.5 for three consecutive sampling periods. Specifically, it is monitored that the values ​​of the first-dimensional component in the L vector at the three consecutive sampling moments t-2Δt, t-Δt, t are [0.55, 0.59, 0.62] respectively, and the values ​​of the second-dimensional component are [0.52, 0.56, 0.58] respectively. Since both L1 and L2 components satisfy The condition that the sampling period is greater than 0.5 for three consecutive times triggers the time synchronization acquisition window. Within this window (the window width is set to 100ms), the acceleration sensor, electromagnetic torque sensor and PT100 platinum resistance temperature sensor are used to synchronously collect the vibration acceleration in the vertical direction of the base, the electromagnetic torque output by the generator and the stator winding temperature to obtain three groups of original signal sequences. These three groups of sequences are preprocessed by sliding average filtering with a window length of 5 sampling points. For example, for the original vibration acceleration sequence [9.8, 10.1, 9.9, 10.5, 10.2, ...] m / s 2 The third data point after filtering is (9.8+10.1+9.9+10.5+10.2) / 5=10.1m / s 2 After processing, three sets of smooth signal sequences are obtained. Then, the standard deviation normalization operation is performed on these three sets of smooth sequences. First, the data within the past 1 second under this working condition is retrieved from the historical database, and the mean and standard deviation of each parameter are calculated, which are: The mean value of vibration acceleration is 8.5m / s 2 , standard deviation 1.5m / s 2 The mean value of electromagnetic torque is 1200 Nm, with a standard deviation of 80 Nm. The mean value of winding temperature is 85°C, with a standard deviation of 5°C. Then, the smoothed value of the current sampling point is normalized. For example, the current value is: vibration acceleration 10.1 m / s 2 , electromagnetic torque 1250Nm, winding temperature 88℃, the normalized values ​​are: acceleration Torque temperature Finally, these three normalized dimensionless values ​​are combined into a multidimensional state vector S = [1.07, 0.625, 0.6].

[0084] S302: Invoke the coupling response vector, construct a three-dimensional input domain of acceleration-torque-temperature, design a triangular membership function to cover multiple parameter variation ranges, establish a nine-rule fuzzy rule base based on the IF-THEN structure, perform fuzzy matching between the multidimensional state vector and the coupling response vector, and generate a fuzzy decision matrix;

[0085] Call the coupling response vector F output by S103 = [42025, 63037, 21012, 42.0, 50.4, 33.6] T With the multi-dimensional state vector S = [1.07, 0.625, 0.6] generated by S301, a three-dimensional input domain of acceleration-torque-temperature is constructed. The domain range of each dimension is set according to the statistical limit value of the historical operating data. For example, the normalized acceleration domain range is [-3, 3], the torque domain is [-3, 3], and the temperature domain is -2, 2]. Five fuzzy subsets are designed for each domain, namely {Negative Large (NB), Negative Small (NS), Zero (ZE), Positive Small (PS), Positive Large (PB)}, and the triangular membership function is used to cover the variation range of these parameters. The bearing in this embodiment is a rolling bearing (NN3020K), and the vertex position of its membership function is dynamically adjusted according to the rules, and the offset is set to ±15%. The setting of this offset is based on the analysis of the vibration characteristics of different bearing types. The nonlinear characteristics of rolling bearings are stronger and require a wider fuzzy coverage range. For example, for the fuzzy subset "positive small (PS)", its standard vertex is located at the normalized value of 1.0. After the offset adjustment of +15%, the vertex moves to 1.0×(1+0.15)=1.15, and a nine-item fuzzy rule base based on the IF-THEN structure is established, as shown in Table 1.

[0086] Table 1 Example of fuzzy control rule base

[0087] Rule Number IF acceleration is AND torque is AND temperature is THEN control output is 1 CP Group (PB) CP Group (PB) CP Group (PB) Negative Big (NB) 2 Small (PS) Small (PS) Small (PS) Negative Small (NS) 3 Zero (ZE) Zero (ZE) Zero (ZE) Zero (ZE) ... ... ... ... ... 9 Negative Big (NB) Negative Big (NB) Negative Big (NB) CP Group (PB)

[0088] As shown in Table 1, some core fuzzy rules are listed, which convert the multidimensional state vector S = [1.07, 0.625, 0.6 and one or more components of the coupled response vector F (for example, the radial force F y Normalized as the fourth input) to perform fuzzy matching, specifically calculating the membership value of each component in the input vector S to each fuzzy subset. For example, the first component of S is 1.07, and its membership to the adjusted "positive small (PS)" fuzzy subset (vertex at 1.15, support interval [0, 2.3]) is calculated as This membership calculation is performed on all input variables and all rules, and then the final trigger strength of each rule is obtained by taking the minimum value of all input memberships of the premise part of each rule (AND operation). The trigger strengths of all rules together constitute a fuzzy decision matrix.

[0089] S303: Based on the fuzzy decision matrix, the centroid method is used to perform defuzzification operations, and the weighted average values ​​of the memberships of multiple control variables are calculated. A compensation action quantization model is established, and a compensation action vector is generated through a multiplication and accumulation operation of the memberships and the domain values.

[0090] Based on the fuzzy decision matrix generated in S302, the center of gravity method (COG) is used to perform the defuzzification operation. The purpose of this operation is to calculate the weighted average of the membership of multiple control variables (for example, the compensation voltage of the active damping system) to determine a unique and clear control output. First, for the domain of the output variable (for example, the compensation voltage adjustment amount, whose range is set to [-5V, +5V]) and its fuzzy subset {negative large (NB), negative small (NS), zero (ZE), positive small (PS), positive large (PB)}, each rule is The trigger intensity (a value between 0 and 1) acts on the output fuzzy subset of its conclusion part. The "peak clipping" method is usually adopted, that is, the trigger intensity value is used to "cut" the top of the membership function graph to obtain a series of clipped output fuzzy set graphs. Then, the clipped graphs generated by all rules are superimposed to form a combined, irregular final output fuzzy graph. The established compensation action quantization model calculates the horizontal coordinate of its geometric center (i.e., center of gravity) by performing multiplication and accumulation operations on the discretized points of the combined graph with the domain value. The calculation formula is: where x i is the i-th discrete point in the output domain, μ(x i ) is the comprehensive membership value of the point, n is the total number of discrete points, and assuming that the calculated compensation voltage adjustment is -2.5 V, and the other control output, that is, the compensation current adjustment of the active damping system, is calculated to be +1.2 A, these clear control values ​​are combined to finally generate the compensation action vector A = [-2.5, 1.2].

[0091] The trigger mechanism of the time synchronization acquisition window is that any two dimensional components in the Lyapunov exponent vector exceed 0.5 at the same time;

[0092] The vertex positions of the triangle membership function are dynamically adjusted according to the bearing type. The vertex offset is set to ±15% for rolling bearings and ±8% for sliding bearings.

[0093] See also Figure 5 , the steps of S4 are as follows:

[0094] S401: Based on the compensation action vector, construct an FPGA hardware description language model, design a state synchronization update mechanism, set timing constraints, and generate an FPGA simulation path;

[0095] Based on the compensation action vector A = [-2.5, 1.2] generated by S303, an FPGA hardware description language (HDL) model is constructed. This model is written in Verilog HDL. Each component in the compensation action vector (such as -2.5V and +1.2A) is quantized into a 16-bit fixed-point binary number, where -2.5V is converted according to the output DAC range [-5V, 5V] and 16-bit resolution to obtain the binary complement "11000000000000000". A state synchronization update mechanism is designed. The core of this mechanism is a finite state machine (FSM) that runs according to the main clock signal and the synchronization pulse. At the rising edge of each clock cycle, the FSM latches the quantized compensation action vector value from the calculation module to the output register, and synchronously updates the relevant state variables within the system. Timing constraints are set, which are implemented in synthesis tools (such as Xilinx Vivado) through XDC.

[0096] The (XilinxDesignConstraints) file is used to define the system master clock period, which is clearly set to 10ns (corresponding to a clock frequency of 100MHz), and requires that the total path delay from the data input register to the final compensation signal output pin must not exceed 8ns. This establishes a time margin (Slack) of no less than 2ns. After a series of processes such as logic synthesis, layout and routing, a bitstream file (.bit) is finally generated that can be downloaded to the target FPGA chip (for example, XilinxArtix-7 series). This file is the FPGA simulation path.

[0097] S402: Calling the FPGA simulation path, establishing a mapping relationship between the geometric dimension deviation parameter and the simulation path, using the Latin hypercube sampling method to generate a deviation parameter combination, injecting the parameter combination into the simulation path for disturbance propagation calculation, and outputting an error disturbance version;

[0098] The FPGA simulation path (.bit file) generated by S401 is called to establish a mapping relationship between geometric dimension deviation parameters and simulation path. The geometric dimension deviation parameters here include gear tooth thickness tolerance, bearing installation centering error, shaft imbalance, etc. For example, the tooth thickness tolerance range is set to [-10μm, +10μm], and the bearing seat centering error range is set to [-20μm, +20μm]. The Latin hypercube sampling (LHS) method is used to generate deviation parameter combinations in the multidimensional space composed of these parameters. This method first divides the value range of each parameter into N intervals with equal probability (for example, N = 20), and then randomly extracts a sample point in each interval of each parameter. Finally, all the sample points are combined. The parameter sample points are randomly combined to ensure that the projections on each parameter dimension are evenly distributed. For example, the set of deviation parameter combinations generated by LHS is as follows: tooth thickness tolerance = +5μm, centering error = -8μm. These parameter combinations are injected into the simulation path running on the FPGA in real time through the AXI4-Lite bus interface. The specific operation is to write the corresponding quantized values ​​to the parameter registers representing these physical quantities inside the FPGA and perform disturbance propagation calculations. That is, run the FPGA model with the injected deviation parameters and continuously monitor their impact on the output of the coupling response vector F and the Lyapunov index vector L. The complete output data stream with this disturbance is recorded to form an error disturbance version.

[0099] S403: Execute the speed / load combination test, input the error disturbance version into the test platform, collect torque fluctuation and vibration spectrum data, construct an error index evaluation matrix, calculate the comprehensive error value of multiple test points through weighted summation, generate an error feedback table and send it back to the fuzzy control algorithm;

[0100] A speed and load combination test is performed. The 100 error disturbance versions corresponding to the 100 sets of different deviation parameters generated in S402 are loaded one by one into a semi-physical simulation test platform. The platform consists of real controller hardware (i.e., a board equipped with an FPGA) and a high-fidelity mechanical system model running on a real-time computer. Each error disturbance version is run under multiple preset working condition test points (for example, the speed is 80%, 100%, and 120% of the rated speed, and the load is 50%, 75%, and 100% of the rated load, forming a total of 3×3=9 test points). The torque fluctuation amplitude and vibration spectrum energy data of key frequency bands (such as the meshing frequency and its harmonic frequency bands) after the system stabilizes at each test point are collected. An error index evaluation matrix is ​​constructed. Each row of the matrix represents an error disturbance version, and each column represents a performance index at a test point (such as the peak-to-peak value of the torque fluctuation and the integrated vibration energy value). The comprehensive error value of each error disturbance version at all test points is calculated by weighted summation. The weight distribution w T = 0.6 (torque ripple) and wV =0.4 (vibration spectrum) is set based on a series of preliminary experiments, using expert scoring combined with the analytic hierarchy process (AHP) to evaluate the relative importance of torque fluctuation and vibration to the long-term stable operation of the system under different working conditions. After statistical analysis, the weight ratio is obtained, and the comprehensive error value E total =w T ×E Torque_norm +w V ×E Vib_norm , where E Torque_norm and E Vib_norm It is the normalized error index. For example, E is calculated for a certain version. total =0.82. Since this value is greater than the preset comprehensive error value judgment threshold of 0.75, it is determined that the system performance under this deviation combination does not meet the robustness requirements. The setting of this threshold of 0.75 is based on the calculation of the comprehensive error value of 30 historical operation samples with known "qualified" performance and 30 "unqualified" performance. It is found that the error values ​​of qualified samples are all lower than 0.73, and the error values ​​of unqualified samples are all higher than 0.78. Therefore, the middle value of 0.75 is taken as the judgment threshold to ensure the accuracy of classification. Finally, the deviation parameters, performance indicators of each test point and comprehensive error values ​​corresponding to all error disturbance versions are summarized to generate an error feedback table, as shown in Table 2. This table is then sent back to the fuzzy control algorithm module of S3 as the data basis for its online adjustment of the fuzzy rule base or membership function.

[0101] Table 2 Error feedback table example

[0102] Deviation combination number Tooth thickness tolerance (μμm) Centering error (μμm) Comprehensive error value 1 +5 -8 0.82 2 -2 +3 0.45 3 +8 +15 0.91 ... ... ... ... 100 -5 -5 0.53

[0103] As shown in Table 2, the tabulated data clearly records the quantitative impact of different geometric dimension deviation combinations on the system stability performance, providing accurate data support for the adaptive optimization of the fuzzy controller.

[0104] The timing constraints set the clock period to 10ns and the setup time margin to ≥2ns;

[0105] The Latin Hypercube sampling method uses a dimensional orthogonalization strategy to improve sampling efficiency by 40%-60% compared to the Monte Carlo method;

[0106] The weight distribution of weighted summation is torque fluctuation weight 0.6, vibration spectrum weight 0.4, and the comprehensive error value judgment threshold is 0.75.

[0107] The above are merely preferred embodiments of the present invention and do not limit the present invention in any other form. Any technician familiar with the profession may use the technical content disclosed above to change or modify it into an equivalent embodiment with equivalent changes and apply it to other fields. However, any simple modification, equivalent change and modification made to the above embodiment based on the technical essence of the present invention without departing from the content of the technical solution of the present invention shall still fall within the scope of protection of the technical solution of the present invention.

Claims

1. A method for analyzing the stability of a microgrid electromechanical coupling system with mechanical vibration, characterized in that: The following steps are involved: S1: Detect the tooth surface contact radius, relative speed and Young's modulus, generate the initial value of contact stiffness through product operation, use the interval combined difference method of forward difference, center difference and backward difference to correct the nonlinear change of stiffness, combine the rotor six-degree-of-freedom displacement data and bearing clearance, and use the finite difference method to output the contact stiffness correction value and coupling response vector; S2: Input the coupling response vector into the Lyapunov index method to perform trajectory sequence regression, calculate the Lyapunov index vector, construct the Z1 / Z2 mode length ratio based on the impedance response, and perform first-order derivative zero point detection and second-order derivative sign determination; S3: According to the positive period of the Lyapunov index vector, the vibration acceleration, electromagnetic torque and winding temperature are collected, combined into a multi-dimensional state vector, and the multi-dimensional state vector and the coupling response vector are input into the fuzzy control algorithm to generate a compensation action vector; S4: Construct an FPGA simulation path based on the compensation action vector, perform a 10μs step state synchronization update, inject geometric dimension deviation parameters to generate an error disturbance version, output an error feedback table through a speed / load combination test and return the fuzzy control algorithm.

2. The method for analyzing the stability of a microgrid electromechanical coupling system of mechanical vibration according to claim 1, characterized in that: The interval combined difference method calculates the stiffness change gradient through forward difference, verifies the gradient direction through central difference, and corrects the nonlinear error through backward difference. The synergistic effect of the three reduces the stiffness correction error by 12%-15%. The Z1 / Z2 module length ratio is a dynamic coupling factor between the real module value and the imaginary module value of the impedance; The positive time period determination threshold is that the exponential component is greater than 0.5 for three consecutive sampling periods.

3. The method for analyzing the stability of a microgrid electromechanical coupling system of mechanical vibration according to claim 2, characterized in that: The steps of S1 are specifically as follows: S101: Detect the tooth surface contact radius, relative speed, and Young's modulus, establish a geometric relationship model between the contact radius and relative speed, construct a Young's modulus compensation function, substitute the three sets of parameters into the stiffness calculation formula, perform multiplication operation, and generate the initial value of the contact stiffness; S102: Based on the initial value of the contact stiffness, a differential iterative calculation model is established, the stiffness change gradient is calculated by forward difference, the gradient direction is verified by central difference, and the nonlinear error is corrected by backward difference to obtain a stiffness correction coefficient; S103: Call the rotor six-degree-of-freedom displacement data and bearing clearance parameters, construct a displacement-stiffness coupling equation, substitute the stiffness correction coefficient into the finite difference method calculation framework, solve the stiffness correction term through the Jacobian matrix, perform matrix iterative operation in combination with the displacement response coupling term, and synchronously output the contact stiffness correction value and the coupling response vector.

4. The method for analyzing the stability of a microgrid electromechanical coupling system of mechanical vibration according to claim 3 is characterized in that: The stiffness calculation formula is: Where E represents the Young's modulus of the gear material in GPa, R c Represents the instantaneous curvature radius of the tooth surface contact point, unit: mm, v r Represents the tangential relative speed of the gear meshing point, unit m / s, denominator is the coupling correction term between contact geometry and kinematic state.

5. The method for analyzing the stability of a microgrid electromechanical coupling system of mechanical vibration according to claim 4 is characterized in that: The steps of S2 are specifically as follows: S201: Based on the coupled response vector, the trajectory sequence is regressed in the time domain using the Lyapunov exponent method. A sliding window mechanism is introduced when calculating the exponential components of multiple dimensions. The exponential decay curve is fitted using the least squares method to generate a Lyapunov exponent vector. S202: Call the Lyapunov exponent vector, extract the real and imaginary moduli of the impedance response, and construct the Z1 / Z2 modulus ratio using the formula: Obtaining dynamic impedance ratios through calculation and generating a dynamic impedance ratio sequence; Among them, R represents the dynamic impedance ratio, Z1 represents the real part modulus of the bearing support structure impedance, unit N·s / m, Z2 represents the imaginary part modulus of the rotor system impedance, unit N·s / m, δ represents the dynamic coupling factor of the bearing clearance, θ represents the impedance phase difference, unit radian, ε represents the ambient temperature interference coefficient, range 0.1-0.3, the denominator is the joint correction term for environmental disturbance and dynamic coupling; S203: Based on the dynamic impedance ratio sequence, a five-point central difference method is used to calculate the first-order derivative, the zero point position is detected by sign change, and a cubic spline interpolation is applied to obtain the sign distribution of the second-order derivative to generate a guide change feature identification group.

6. The method for analyzing the stability of a microgrid electromechanical coupling system of mechanical vibration according to claim 5, characterized in that: The sliding window mechanism sets the window length to 2.5 times the system base frequency period.

7. The method for analyzing the stability of a microgrid electromechanical coupling system of mechanical vibration according to claim 6, characterized in that: The steps of S3 are specifically as follows: S301: establishing a time synchronization acquisition window for vibration acceleration, electromagnetic torque, and winding temperature according to the positive value period of the Lyapunov exponent vector, processing the original signal using a sliding average filter, performing a standard deviation normalization operation on the three sets of parameters, and combining them into a multidimensional state vector; S302: Calling the coupling response vector, constructing a three-dimensional input domain of acceleration-torque-temperature, designing a triangular membership function to cover multiple parameter variation ranges, establishing a nine-rule fuzzy rule base based on an IF-THEN structure, performing fuzzy matching between the multidimensional state vector and the coupling response vector, and generating a fuzzy decision matrix; S303: Based on the fuzzy decision matrix, a centroid method is used to perform defuzzification operations, a weighted average of the membership of multiple control variables is calculated, a compensation action quantization model is established, and a compensation action vector is generated through a multiplication and accumulation operation of the membership and the domain value.

8. The method for analyzing the stability of a microgrid electromechanical coupling system of mechanical vibration according to claim 7, characterized in that: The triggering mechanism of the time synchronization acquisition window is that any two dimensional components in the Lyapunov exponent vector simultaneously exceed 0.5; The vertex position of the triangle membership function is dynamically adjusted according to the bearing type. The vertex offset of rolling bearings is set to ±15%, and that of sliding bearings is set to ±8%.

9. The method for analyzing the stability of a microgrid electromechanical coupling system of mechanical vibration according to claim 8, characterized in that: The steps of S4 are specifically as follows: S401: Based on the compensation action vector, construct an FPGA hardware description language model, design a state synchronization update mechanism, set timing constraints, and generate an FPGA simulation path; S402: calling the FPGA simulation path, establishing a mapping relationship between geometric dimension deviation parameters and the simulation path, generating a deviation parameter combination using a Latin hypercube sampling method, injecting the parameter combination into the simulation path for disturbance propagation calculation, and outputting an error disturbance version; S403: Execute a speed / load combination test, input the error disturbance version into the test platform, collect torque fluctuation and vibration spectrum data, construct an error index evaluation matrix, calculate the comprehensive error value of multiple test points through weighted summation, generate the error feedback table and send it back to the fuzzy control algorithm.

10. The method for analyzing the stability of a microgrid electromechanical coupling system of mechanical vibration according to claim 9, characterized in that: The timing constraint condition sets the clock period to 10ns and the setup time margin to ≥ 2ns; The Latin hypercube sampling method adopts a dimensional orthogonalization strategy to improve the sampling efficiency by 40%-60% compared with the Monte Carlo method; The weight distribution of the weighted summation is as follows: the torque fluctuation weight is 0.6, the vibration spectrum weight is 0.4, and the comprehensive error value determination threshold is 0.75.

Images