Method and system for model identification of a magnetic levitation rotor system based on a frequency response model

Abstract

The application discloses a kind of based on frequency response model's magnetic suspension rotor system model identification method and system, comprising: using harmonic signal excitation magnetic suspension rotor system, real-time acquisition magnetic suspension rotor system in each input-output channel under the time-domain signal data sample generated by excitation;The time-domain data sample is converted into frequency domain data, and the frequency response data set for model identification is formed;Unbiased estimation is carried out to the frequency response data set using optimization algorithm, and empirical transfer function estimation model is established;With normalized root mean square error evaluation empirical transfer function estimation model precision, and the prediction output of this model is compared with the time-domain signal data sample, verify the accuracy and reliability of model.The application is excited by accurate harmonic and real-time data acquisition, combined with data transformation technology, effectively convert time-domain signal into frequency domain data, significantly improve the efficiency and accuracy of performance analysis and control design of magnetic suspension rotor system.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of magnetic bearing identification methods, and particularly relates to a magnetic suspension rotor system model identification method and system based on a frequency response model. BACKGROUND

[0002] Magnetic suspension bearings are a new type of high-performance bearing that uses controllable electromagnetic fields to stably suspend a rotor at a given position, and have the outstanding advantages of non-contact and actively controllable stiffness, and are an ideal supporting method for high-precision, long-life high-speed rotor systems. Accurate knowledge of the system model is the basis for using advanced control methods to improve the control accuracy, stability, and reliability of magnetic bearing systems. For example, applying advanced control methods such as robust control and adaptive control to magnetic bearing systems can improve the control accuracy and stability of the system; fault diagnosis and fault-tolerant control of the magnetic bearing system can greatly improve the reliability of the system, and the use of these methods is premised on knowing the model of the magnetic bearing system, and identification is an important method for obtaining the system model.

[0003] Existing identification methods include a magnetic force parameter identification method based on a Gaussian modulation function method, a composite identification method, a frequency domain identification method, and a system modal parameter identification method based on the combination of random subspace and empirical mode decomposition, etc., which are important research results in the field of magnetic suspension rotor system model identification in recent years. These methods each have their own characteristics, such as the Gaussian modulation function method which avoids direct processing of differential signals and simultaneously avoids the integral initial value problem; the composite identification method which corrects the reduced-order mathematical model obtained by finite element analysis through system identification of the magnetic suspension flexible rotor from dynamic characteristic analysis data by the variable LEVY method; the frequency domain identification method which obtains the frequency characteristics of the system by model identification of the actual system, and performs model fitting on the identification data; and the method combining random subspace and empirical mode decomposition which can effectively identify the first and second order natural frequencies of the rotor system under operating conditions.

[0004] However, although these methods improve the accuracy and efficiency of model identification to some extent, there are still some shortcomings. First, many methods still face challenges when dealing with complex systems. For example, there is a large difference between the theoretical model of the magnetic suspension bearing-rotor system and the actual characteristics, making it difficult for the controller designed based on the theoretical model to obtain good control performance; second, some methods may be affected by experimental condition limitations in actual application. For example, in the process of unbalance response method for magnetic suspension rotor system parameter identification, high-frequency noise filtering methods need to be considered; in addition, although some methods can improve the accuracy of model identification, they have high computational complexity, which may affect the speed and efficiency of model identification, for example, although the joint simulation model of the flexible magnetic suspension rotor system established by using virtual prototyping technology can more accurately reflect the operating state of the rotor, the modeling and simulation process is relatively complex. SUMMARY

[0005] In view of the above defects or improvement needs of the prior art, the purpose of the present application is to provide a magnetic levitation rotor system model identification method and system, which effectively converts time domain signals into frequency domain data by precise harmonic excitation and real-time data acquisition, combined with data transformation technology, and uses an optimization algorithm to make unbiased estimation of the frequency response of the system, establish an accurate empirical transfer function estimation model, and through normalized root mean square error evaluation, ensure that the identified model has high prediction accuracy and reliability, significantly improving the performance analysis and control design efficiency and accuracy of the magnetic levitation rotor system.

[0006] In the first aspect, the present application provides a magnetic levitation rotor system model identification method based on a frequency response model, characterized in that it comprises:

[0007] S100, exciting the magnetic levitation rotor system with a harmonic signal, and collecting real-time time domain signal data samples under each input and output channel generated by the excitation in the magnetic levitation rotor system;

[0008] S200, converting the time domain data samples into frequency domain data to obtain a frequency response data set for model identification, and establishing a frequency response model based on the frequency response data set;

[0009] S300, using an optimization algorithm to make unbiased estimation of the frequency response model, and establishing an empirical transfer function estimation model;

[0010] S400, evaluating the accuracy of the empirical transfer function estimation model with normalized root mean square error, and comparing the prediction output of the model with the time domain signal data samples to verify the accuracy and reliability of the model.

[0011] Further, in step S100, the sine displacement harmonic instruction The formula is:

[0012] Where t=cnt·T s ,

[0013] In the formula, is the harmonic instruction, t h and T s are preset values, respectively representing the amplitude of the excitation signal, the duration of single excitation, and the control step size; and cnt are variables, the former is the frequency of the sinusoidal signal at the nth excitation, and the latter is the count value; the variable superscript k = 1, 2, 3, 4 respectively represents the relevant signal or transfer function of the X, Y degrees of freedom of the DE end radial magnetic bearing A and the X, Y degrees of freedom of the NDE end radial magnetic bearing B, and is marked as DE-X, DE-Y, NDE-X and NDE-Y respectively.

[0014] Further, the harmonic signal excitation in step S100 needs to be performed multiple times, and is combined with the magnetic bearing control host computer to complete the acquisition of the signal, specifically including:

[0015] S101, initialize the excitation channel, set k = 1, and the signal excitation will be applied to the control loop in the DE-X direction;

[0016] S102, initialize the excitation frequency, set n = 1, and the host computer sends the initial sinusoidal harmonic parameters to the k-channel control loop, and the controller generates sinusoidal harmonic instructions according to the sinusoidal displacement harmonic instruction formula;

[0017] S103, the four radial degrees of freedom of the magnetic suspension rotor system produce harmonic response under the action of the sinusoidal harmonic excitation, and the controller stores the control current instructions of the excitation channel and the displacement voltage feedback signals on the four feedback channels

[0018] S104, after the response ends, the controller sends the stored time domain data to the host computer end for re-saving, and the host computer updates and sends the excitation frequency to the controller;

[0019] S105, the system enters the next harmonic response stage, and repeats steps S102-S104 until the excitation frequency reaches the upper limit value 1000Hz.

[0020] Further, step S200 includes:

[0021] S201, cooperate with the magnetic bearing control host computer to analyze the time domain data, and convert the time domain signal to a frequency domain signal;

[0022] S202, set k = k + 1, enter the frequency response test of the next excitation channel, and complete the data acquisition of the 4x4 frequency response data matrix between the DE-X, NDE-X, DE-Y and NDE-Y four input and output channels.

[0023] Further, the frequency response model in step S200 is:

[0024]

[0025] In the formula, H frdrepresenting a frequency response model, ωn represents the angular frequency of the nth frequency test point; N f N represents the total number of frequency test points; and represents the subsystem frequency response model of the ith input of the system to the jth output.

[0026] Further, the step S300 comprises:

[0027] S301, determining the structure of the empirical transfer function estimation model, simplifying the complexity of model identification;

[0028] S302, constructing a covariance-based weighted cost function for indicating the iteration direction of the parameter set;

[0029] S303, selecting an initial parameter set and an initial damping coefficient, setting an upper limit of iteration times and a convergence threshold;

[0030] S304, calculating the Jacobian matrix of the weighted cost function;

[0031] S305, updating the initial parameter set according to the Jacobian matrix, wherein the damping coefficient is adaptively adjusted with the iteration process to optimize the convergence speed;

[0032] S306, iteratively updating the parameters, comparing the difference of the cost function values of adjacent iterations according to the updated parameter set to perform the convergence test, when the difference is less than the preset threshold, considering that the model has converged, terminating the algorithm, and outputting the model identification result.

[0033] Further, the structure of the empirical transfer function estimation model needs to be established in the step S301, which can be specifically written as:

[0034]

[0035] In the formula, G represents the empirical transfer function estimation model; θ ij represents the parameter variable to be estimated in the sub-transfer function

[0036] Further, in the model identification problem, a common structure form is to set it as a standard rational polynomial form of the transfer function, that is:

[0037]

[0038] In the formula, s η , s η-1 ,..., s represents different powers of s, which is a variable in the complex frequency domain; η is the order of the denominator polynomial; α η-1 ​​..., a1, a0 are coefficients of denominator polynomial of transfer function; β ξ ..., β1, β0 are coefficients of numerator polynomial of transfer function. ξ-1 ..., β1, β0 are coefficients of numerator polynomial of transfer function.

[0039] Further, step S301 is to simplify the complexity of model identification, and the experience transfer function is converted into each order modal combination transfer function form, that is:

[0040]

[0041] In the formula, K gin is the gain coefficient of the subsystem model; p α , z β are the first-order pole and zero point of the system respectively; n α , n δ indicate the modal order considered in the identification process; n β , n γ are the number of first-order and second-order links of the numerator respectively; ζ Mγ , ω Mγ , ζ Mδ , ω Mδ represent the damping ratio and natural frequency of the second-order links of the numerator and denominator respectively.

[0042] Further, after the power amplifier is equivalent to a first-order low-pass filter link and the displacement sensor is equivalent to a fixed gain, the experience transfer function can be finally expressed as:

[0043]

[0044] Among the parameter set θ to be identified, there are 24 parameters in total, and the specific expression is:

[0045]

[0046] In the formula, K ij indicates the gain coefficient of the subsystem model; indicates the pole, which can reflect the bandwidth of the power amplifier; [p 1x , p 2x , p 1y , p 2y ] reflect the ratio of displacement stiffness coefficient to mass inertia of each radial degree of freedom magnetic bearing; indicates the zero point, which is introduced by mechanical coupling effect.

[0047] Further, step S302 specifically includes:

[0048] Define the unit error matrix:

[0049]

[0050] Further, a cost function V(θ) is constructed to obtain the overall difference between the two models in the whole frequency test section, including:

[0051]

[0052] wherein, is the output error weight matrix at frequency point ; tr(·) represents the trace operator of a matrix; is the frequency weight coefficient, used to adjust the influence degree of the overall deviation of the model at different frequency points in the cost function calculation; N' f is the number of calculation frequency points.

[0053] Further, step S304 includes:

[0054] The Jacobian matrix J(θ) contains the partial derivatives of V(θ) with respect to each parameter in the parameter set θ, and the elements J (q) in the matrix can be expressed as:

[0055]

[0056] Further, the parameter set update formula in step S305 is:

[0057] θ (k+1) = θ (k) - [J(θ (k) ) T J(θ (k) )+λ·diag(J(θ (k) ) T J(θ (k) )] -1 · (θ (k) ) T E(θ (k) )

[0058] wherein, θ (k) and θ (k+1) are the parameter sets at the current iteration step and the next iteration step, respectively; λ is a dynamically adjusted damping factor.

[0059] Further, the calculation formula of the normalized root mean square error in step S400 is:

[0060]

[0061] wherein, represents the normalized root mean square error.

[0062] In a second aspect, the application further provides a magnetic suspension rotor system model identification system based on a frequency response model, which is implemented by using any step of the magnetic suspension rotor system model identification method based on a frequency response model.

[0063] A first module is configured to collect time domain signal data samples under each input and output channel generated by excitation in the magnetic suspension rotor system in real time.

[0064] A second module is configured to convert the time domain data samples obtained by the first module into frequency domain data to obtain a frequency response data set for model identification, and establish a frequency response model based on the frequency response data set.

[0065] A third module is configured to use an optimization algorithm to perform unbiased estimation on the frequency response model obtained by the second module, and establish an empirical transfer function estimation model.

[0066] A fourth module is configured to evaluate the precision of the empirical transfer function estimation model obtained by the third module in terms of normalized root mean square error, and compare the predicted output of the empirical transfer function with the time domain signal data samples to verify the accuracy and reliability of the model.

[0067] The application has the following beneficial effects:

[0068] 1. The application provides a magnetic suspension rotor system model identification method, which effectively converts time domain signals into frequency domain data by using precise harmonic excitation and real-time data collection, combining data transformation techniques, and using an optimization algorithm to perform unbiased estimation on the frequency response of the system to establish an accurate empirical transfer function estimation model, and ensuring that the identified model has high prediction accuracy and reliability by evaluating the normalized root mean square error, thereby significantly improving the performance analysis and control design efficiency and accuracy of the magnetic suspension rotor system.

[0069] 2. The application provides a magnetic suspension rotor system model identification method, which improves the accuracy and robustness of model identification by designing a covariance-based weighted cost function to adapt to the correlation and noise level between different channels.

[0070] 3. The application provides a magnetic suspension rotor system model identification method, which realizes fast and stable convergence of model parameters by adaptively adjusting the damping coefficient and the fine control of the iteration process.

[0071] Additional aspects and advantages of the application will be in part apparent and in part pointed out hereinafter. BRIEF DESCRIPTION OF DRAWINGS

[0072] The above and / or additional aspects and advantages of the present application will become apparent and more readily appreciated from the following description, taken in conjunction with the following drawings of which:

[0073] Figure 1 is a flow chart of a magnetic suspension rotor system model identification method in an embodiment of the present application;

[0074] Figure 2 is a principle block diagram of a magnetic bearing system model identification in an embodiment of the present application;

[0075] Figure 3 is a flow chart of a frequency response data set collection of a magnetic suspension rotor system in an embodiment of the present application;

[0076] Figure 4 is a frequency response test result of a 75kW motor magnetic suspension rotor system in an embodiment of the present application;

[0077] Figure 5 is a comparison of amplitude frequency response curves of an identified model, a theoretical model and a frequency response test result in an embodiment of the present application. DETAILED DESCRIPTION

[0078] The present application will be further described by way of illustration with reference to the accompanying drawings and embodiments. It is to be understood that the embodiments described herein are merely illustrative of the present application and should not be construed as limiting the present application. It is also to be understood that only the parts and aspects of the drawings that are necessary for illustrating the present application are shown.

[0079] It is to be understood that the singular forms "a," "an," and "the" include plural referents unless the context clearly dictates otherwise. It is further understood that the terms "comprising," "including," "containing," or "having" and the like, when used herein, mean "including but not limited to," and are not intended to (and do not) exclude other moieties, additives, components, integers or steps.

[0080] It is to be understood that the singular forms "a," "an," and "the" include plural referents unless the context clearly dictates otherwise. It is further understood that the terms "comprising," "including," "containing," or "having" and the like, when used herein, mean "including but not limited to," and are not intended to (and do not) exclude other moieties, additives, components, integers or steps.

[0081] The present application provides a magnetic suspension rotor system model identification method, which effectively converts time domain signals into frequency domain data by precise harmonic excitation and real-time data acquisition, combined with data transformation technology, and utilizes an optimization algorithm to make unbiased estimation of the frequency response of the system, establish an accurate empirical transfer function estimation model, and ensure that the identified model has high prediction accuracy and reliability through normalized root mean square error evaluation, thereby significantly improving the performance analysis and control design efficiency and accuracy of the magnetic suspension rotor system.

[0082] Embodiment 1

[0083] As shown in Figure 1 The present application provides a magnetic suspension rotor system model identification method and system based on frequency response model, which comprises:

[0084] S100, a harmonic signal is used to excite the magnetic suspension rotor system, and time domain signal data samples under each input and output channel of the magnetic suspension rotor system are collected in real time;

[0085] S200, the time domain data samples are converted into frequency domain data to form a frequency response data set for model identification;

[0086] S300, an optimization algorithm is used to make unbiased estimation of the frequency response data set to establish an empirical transfer function estimation model;

[0087] S400, the accuracy of the empirical transfer function estimation model is evaluated by normalized root mean square error, and the prediction output of the model is compared with the time domain signal data samples to verify the accuracy and reliability of the model.

[0088] The step S100 specifically comprises:

[0089] The frequency response data of the magnetic suspension rotor system is the original form of its dynamic characteristics and the data basis for model identification. In this embodiment, the system is excited by a harmonic signal to generate a harmonic response, and specific time domain signal data samples under each input and output channel of the MIMO system are collected in real time. Finally, the collected time domain data samples are converted into frequency domain data to obtain the frequency response data set required for model identification.

[0090] Furthermore, three types of harmonic excitation are commonly used in model identification: stepped sine signals, linear frequency sweep signals, and superimposed sine signals. Although frequency response testing under stepped sine signal excitation is time-consuming, its frequency and amplitude design flexibility is large, the collected data has a high signal-to-noise ratio, the data post-processing process is simple, and the hardware storage capacity requirements are lower. Therefore, this embodiment uses a stepped sine signal as the excitation signal for the frequency response test and combines it with the developed magnetic bearing control host computer to complete the collection of frequency response data sets for each input and output of the magnetic levitation rotor system.

[0091] Furthermore, due to the open-loop instability of the system itself, the experiment needs to be carried out under the condition of position closed loop, that is, the rotor is statically suspended. Figure 2 As shown in the experiment, the position instruction of the controller r Superimposed sinusoidal displacement harmonic command Apply a step-by-step sinusoidal excitation signal and collect the control current command output by the position controller when the harmonic response occurs in real time And the displacement feedback voltage signal output by the displacement sensor Among them, harmonic instruction include:

[0092]

[0093] Where, t h and T s are preset values, representing the amplitude of the excitation signal, the duration of a single excitation, and the control step length respectively; and cnt are variables, the former is the frequency of the sinusoidal signal at the nth excitation, and the latter is the count value.

[0094] Furthermore, in the frequency response test, in order to avoid problems such as power amplifier saturation, excessive displacement vibration, and modal excitation, which would cause serious nonlinear distortion in the system and invalidate the identification results, it is necessary to reasonably configure the values ​​of the above parameters. After several experimental experiments, the excitation signal parameters for the magnetic levitation rotor of the 75kW motor were finally determined, as shown in Table 4.1:

[0095] Table 4.1

[0096]

[0097] Furthermore, if Figure 3 As shown, to obtain the complete frequency response model H frd , multiple frequency response test experiments need to be carried out, and the real-time acquisition of the input and output channels of the magnetic levitation rotor MIMO system during the frequency sweep process is required. and time domain data samples, and finally the data is converted in frequency domain and calculated. The embodiment cooperates with the frequency response test module of the magnetic bearing control host computer to automatically realize harmonic excitation generation, data acquisition and post-processing, shortens the test and data processing time, and the specific operation process includes:

[0098] S101, excitation channel initialization, setting k=1, signal excitation is applied to the control loop in the DE-X direction;

[0099] S102, excitation frequency initialization, setting n=1, the host computer sends the initial sinusoidal harmonic parameters to the k-channel control loop, and the controller generates sinusoidal harmonic instructions according to the formula 4.1;

[0100] S103, the 4 radial degrees of freedom of the magnetic suspension rotor system generate harmonic response under the action of sinusoidal harmonic excitation, and the controller stores the control current instructions of the excitation channel and the displacement voltage feedback signals on the 4 feedback channels

[0101] S104, after the response ends, the controller sends the stored time domain data to the host computer end to save again, and the host computer updates and sends the excitation frequency to the controller, the system enters the next harmonic response stage, and the step S102 is repeated until the excitation frequency reaches the upper limit value 1000Hz;

[0102] Further, the step S200 includes:

[0103] S201, cooperating with the magnetic bearing control host computer, analyzing the time domain data, and converting the time domain signal into a frequency domain signal;

[0104] S202, setting k=k+1, entering the frequency response test of the next excitation channel, completing the collection of the 4x4 frequency response data matrix among DE-X, NDE-X, DE-Y and NDE-Y four input and output channels and establishing a data set.

[0105] Further, the frequency response model of the MIMO magnetic suspension rotor system (only considering the radial degrees of freedom) obtained is the frequency domain input-output characteristic from the control current instruction vector output by the position controller to the displacement feedback voltage signal vector output by the displacement sensor, that is:

[0106]

[0107] In the formula, is the angular frequency of the nth frequency test point; N f is the total number of sweep points, which is 220 in the frequency response experiment in the embodiment; and indicates the i-th input to the jth output The subsystem frequency response model at any frequency point ω needs to be obtained by Discrete Fourier Transform (DFT) of the input / output time domain signals, the calculation process is shown in equation (4.3).

[0108]

[0109] In the formula, N m is the total number of sampling points of the time domain signal, that is, N m = (t h - 0.1) / T s ; T s is the signal sampling time interval; m is the sampling count value; kk represents the frequency index in the frequency domain after Discrete Fourier Transform (DFT), which is used to identify the index value of a certain frequency point. Note that, in order to avoid the influence of transient process, the data in the initial 0.1s period of the sampling data is discarded.

[0110] Further, as shown in Figure 4 , according to the obtained frequency response model matrix H frd (ω), the frequency response characteristics (amplitude frequency and phase frequency) of the 75kW motor MIMO magnetic suspension rotor system are drawn, in the figure, the four frequency response characteristic curves of the diagonal elements represent the input-output response of the system in a single control degree of freedom, which is actually approximately equivalent to the frequency response characteristics of the SISO system in the DE-X, NDE-X, DE-Y and NDE-Y degrees of freedom after the magnetic bearing system is forced to decouple; the frequency response characteristics of the non-diagonal elements reflect the coupling effect between different degrees of freedom:

[0111] Among them, the sub-diagonal and super-diagonal elements reflect the mechanical coupling effect between the DE and NDE ends in the X (or Y) degree of freedom, and their amplitude frequency characteristic amplitudes are generally lower than -25dB, indicating that the mechanical coupling effect is weak, which proves the rationality of the dispersion PID control of the rotor from the side; the remaining non-diagonal elements reflect the mutual coupling effect between the X and Y degrees of freedom. In theory, when the rotor is statically suspended, they should be completely decoupled, but in fact, due to the influence of factors such as machining and assembly accuracy, magnetic circuit coupling, etc., the experimental results show that there is still a weak coupling effect (amplitude below -40dB).

[0112] It is worth noting that the phase of each frequency response characteristic curve of the diagonal elements at the frequency point of 825Hz appears a significant change and a weak damping resonance phenomenon, which reveals that the first order bending modal frequency of this magnetic suspension rotor system is 825Hz, which is far lower than the highest rotation frequency of the motor, so the rotor mainly presents the characteristics of a rigid rotor during normal operation.

[0113] The specific steps of step S300 include:

[0114] wherein, in obtaining the frequency response model H frd After that, the Levenberg-Marquardt optimization algorithm is used herein to make unbiased estimation on it to obtain the empirical transfer function estimation model Note that the L-M method is an estimation method of least square estimation of regression parameters in nonlinear regression, which combines the advantages of gradient descent method and Newton method, and can effectively balance the advantages of global search and local fine search, so it is widely used in parameter estimation tasks of complex models in system identification.

[0115] Further, the specific steps of the model identification (i.e. parameter estimation) include:

[0116] S301, determine the structure of the empirical transfer function estimation model, simplify the complexity of model identification, including:

[0117] Further, the dimension of the ETFE parameter model to be estimated should be consistent with the frequency response model H frd , specifically represented by the 4x4 transfer function matrix as shown below.

[0118]

[0119] wherein, represents the empirical transfer function estimation model; θ ij represents the parameter variables to be estimated in the sub-transfer function .

[0120] Further, in the model identification problem, the common construction form is to set it as the standard rational polynomial form of the transfer function, that is:

[0121]

[0122] wherein, β ξ , β ξ-1 ,..., β0 are the coefficients of the numerator polynomial of the transfer function; s η , s η-1 ,..., s represent different powers of s, which is a variable in the complex frequency domain; η is the order of the denominator polynomial; a η-1 ,..., a1, η0 are the coefficients of the denominator polynomial of the transfer function.

[0123] The problem of such structure is that there is a potential conflict between the complexity of model parameters and the full characterization of dynamic characteristics: the higher the preset system order (ξ and η), the higher the model accuracy, but the number of parameters to be estimated (ξ0,..., ξ n ) and (η0,..., η n-1 ) will be more, and the parameter value range is difficult to limit, making the model identification difficult; the identification of low-order models is relatively easy, but the identified model often cannot fully characterize the dynamic characteristics of the MIMO system. From the perspective of structural dynamics, the standard rational polynomial form can be converted into a modal combination transfer function form for each order, that is:

[0124]

[0125] In the formula, K gain is the gain coefficient; p α , z β are the parameters of the system zero point; n α , n δ represent the modal order considered in the identification process; n β , n γ are the number of first-order and second-order elements in the numerator, respectively; ζ Mγ , ω Mγ , ζ Mδ , ω Mδ represent the damping ratio and natural frequency of the first-order and second-order elements in the numerator and denominator, respectively.

[0126] The multiplication part in the formula is composed of first-order and second-order elements, the former corresponds to the rigid modal of the rotor and the equivalent model of the electrical elements such as displacement sensor and switching power amplifier, and the latter corresponds to the flexible modal model of the rotor. Compared with the rational polynomial form, the modal combination transfer function form has the following two advantages:

[0127] 1. The model complexity is limited by the number of modes n α +n δ . Designers can directly determine the order of the ETFE model according to the number of rotor modes expected to be obtained through identification;

[0128] 2. The parameter value range does not change with the system order, but is only related to the inherent dynamic characteristics of the system, which is easier to limit. For example, the zero-pole parameters [z β , p a ] of the first-order element, which are related to the bearing stiffness coefficient, mass, current loop bandwidth and other physical quantities; the damping ratio and undamped natural frequency [ζ M(δ，γ) , ω M(δ，γ) ] of the second-order element essentially directly reflect the damping ratio (underdamped form, typical value distributed between 0-0.1) and modal frequency information of the flexible modal.

[0129] Further, since the control problem of rigid magnetic suspension rotor is mainly discussed in this embodiment, the second order link in the above formula is ignored, and the power amplifier is equivalent to a first order low pass filter link and the displacement sensor is equivalent to a fixed gain, then the empirical transfer function estimation model can be expressed as:

[0130]

[0131] Among the parameter set θ to be identified, there are 24 parameters, and the expression is:

[0132]

[0133] In the formula, K ij is the gain coefficient of the subsystem model; the pole reflects the power amplifier characteristics (i.e. the bandwidth of the power amplifier); [p 1x , p 2x , p 1y , p 2y ] reflect the ratio of displacement stiffness coefficient to mass inertia of each radial degree of freedom magnetic bearing, that is zero point is introduced by mechanical coupling effect, and its value is also associated with , and for the magnetic suspension rotor system with weak mechanical coupling effect, should be similar to the size of poles p 1x，y , p 2x，y . In fact, the bandwidth of the power amplifier is generally hundreds of hertz, and the displacement stiffness coefficient of the magnetic bearing of the 75kW motor used in this embodiment is generally in the order of n×10 5 N / m, and the rotor weight is about tens of kilograms. In addition, the value range of K ij can be estimated according to formula (4.9), wherein the value range of each parameter in formula (4.8) is listed in table 4.2.

[0134] In fact, the bandwidth of the power amplifier is generally hundreds of hertz, and the displacement stiffness coefficient of the magnetic bearing of the 75kW motor used in this embodiment is generally in the order of n×10 5 N / m, and the rotor weight is about tens of kilograms.

[0135]

[0136] Table 4.2

[0137]

[0138]

[0139] Step S302 specifically includes:

[0140] Since the L-M method requires the cost function to indicate the iteration direction of the parameter set θ, a covariance-based weighted cost function form is constructed here:

[0141]

[0142] where, represents the response error of the ETFE model at the specific parameter set θ, and the frequency response model H frd at the frequency test point . Note that it is a complex matrix, so it reflects the difference between the amplitude and phase of the two response models. Further, a (scalar) error cost function V(θ) is constructed to obtain the overall difference between the two models in the entire frequency test section, and the specific expression includes:

[0143]

[0144] where, is the output error weight matrix at the frequency point , which realizes the importance allocation of each subsystem in the MIMO system in the identification process by weighting each element in the error matrix ; tr(·) represents the trace operator of the matrix, which obtains the sum of the modulus of each element in the weighted error matrix (at the frequency point ), reflecting the overall deviation of the model; is the frequency weight coefficient, which is used to adjust the influence of the overall deviation of the model at different frequency points in the calculation of the cost function; N′ f is the number of calculated frequency points.

[0145] Further, it can be seen from Figure 5 that after 650 Hz, the frequency response characteristics of the system have been affected by the resonance of the first-order bending mode. Since this paper only identifies the rigid magnetic suspension rotor model, only the model deviation in the range of 5-650 Hz is considered in the calculation of V(θ). For an adaptive weight matrix is introduced here:

[0146]

[0147] Further, this weighting method can automatically adjust the sensitivity to each output error in the model estimation process, which is considered to provide the best linear unbiased (minimum variance) estimation in statistics. In essence, it calculates the inverse of the covariance matrix of the unit error, and by giving higher weights to those error items with smaller variances, it reduces the overfitting of high-noise data points in the model estimation process. On the other hand, it can be observed that Figure 5It can be seen that the gain of the amplitude-frequency characteristic is high (around 0 dB) in the low frequency band, while the gain value is very low (lower than -30 dB, at the noise level) in the high frequency band. In order to balance the contributions of errors in different frequency bands to V(θ), the frequency weight coefficient of the frequency band needs to be appropriately increased. In this identification problem, The value of θ obeys the conditional expression shown in the following table:

[0148]

[0149] In summary, the model identification problem of the rigid magnetic suspension rotor can be explained as follows: within a limited parameter range (as shown in Table 4.2), find an optimal parameter set that minimizes the following formula:

[0150]

[0151] wherein

[0152]

[0153] S303, select an initial parameter set and an initial damping coefficient, and set an upper limit of iteration times and an error convergence threshold.

[0154] S304, calculate the Jacobian matrix of the weighted cost function, specifically including:

[0155] The Jacobian matrix J(θ) contains the partial derivative of the cost function V(θ) with respect to each parameter θ ij The calculation formula of the element J (q) in the matrix can be expressed as:

[0156]

[0157] S305, update the initial parameter set according to the Jacobian matrix, specifically including:

[0158] During the iteration process of the L-M algorithm, the update rule of the parameter set θ follows the formula:

[0159] θ (k+1) =θ (k) -[J(θ (k) ) T J(θ (k) )+λ·diag(J(θ (k) ) T J(θ (k) ))] -1 ·J(θ(k)) T E(θ (k) ) (4.16)

[0160] Here, λ is a dynamically adjusted damping factor that helps control the stability and convergence rate of the optimization process. Its value generally changes adaptively during the iteration process. Here, a larger λ value is selected in the initial stage and gradually reduced as the iteration progresses to accelerate the convergence of the algorithm.

[0161] S306: When the number of iterations reaches the set value, it is considered that the convergence condition has been met and the algorithm terminates. The parameter set of the last generation is recorded as Substitute it into formula (4.7) and output the model identification result of the LM algorithm:

[0162] Finally, after 100 iterations, an optimal model matrix was obtained: The transfer functions corresponding to each element in the matrix are listed in equations (4.17) to (4.24). And theoretical model and frequency response test results H frd The amplitude-frequency characteristic curve of Figure 5 shown.

[0163] Further analysis of the four amplitude-frequency characteristics of the diagonal elements reveals that within the 0-650Hz identification frequency band, the identification model and the frequency response test results are almost identical. Although the (rigid) identification model cannot follow the actual amplitude-frequency curve after 650Hz due to the influence of flexible modes, the overall gain roll-off trend remains consistent. However, there is a certain gap between the theoretical model and the frequency response test results. Its amplitude-frequency curve consistently lies above the test results, and this amplitude deviation becomes more pronounced as the frequency increases.

[0164] Furthermore, regarding the amplitude-frequency curves of the superdiagonal and subdiagonal elements, both the identification model and the theoretical model are unable to achieve a good fit for the amplitude-frequency characteristics of the frequency response test results. This is because the mechanical coupling effect of the magnetically suspended rotor of the 75kW motor is weak, so the signal-to-noise ratio of the frequency response test data is low, and there are random fluctuations in the amplitude within a small neighborhood frequency range. The transfer function of the rigid rotor model is of a lower order (5th order), so it cannot fit the actual test curve. However, compared with the theoretical model, the identification model still has a higher degree of fit, and its amplitude-frequency curve is in the middle of the frequency response test results; while the amplitude-frequency curves of the two superdiagonal elements of the theoretical model show a clear amplitude "concave" phenomenon in the mid-frequency range (around 200Hz). This is due to the mismatch of the theoretical parameters, which causes imaginary zeros to appear in the numerator of the transfer function.

[0165]

[0166] S400, the fitting quality between the identification model, the theoretical model and the frequency response test model is quantitatively evaluated by a normalized root mean square error; the calculation method of the NRMSE error is referred to formula (4.25), which can provide a standardized error measurement for data of different scales or orders, and is therefore widely used in performance evaluation of the identification model;

[0167]

[0168] Further, the error calculation results are listed in Table 4.3 and Table 4.4. The NRMSE error of the diagonal element of the identification model is close to 0, and the error of the theoretical model is about 6 times that of the identification model.

[0169] Further, this result shows that the error of the diagonal element of the identification model is significantly smaller than the fluctuation level of the test data, and has better model accuracy and reliability. On the other hand, due to the limitation of the response amplitude and the model order, the NRMSE of the superdiagonal element and the subdiagonal element of the identification model is relatively large, and the error distribution is between 1.3 and 1.6, but compared with the theoretical model, the error performance is still better.

[0170] Table 4.3

[0171]

[0172] Table 4.4

[0173]

[0174] Further, the application also provides a magnetic suspension rotor system model identification system based on a frequency response model, which is realized by using any step of the above-mentioned magnetic suspension rotor system model identification method based on a frequency response model, and includes:

[0175] A first module is configured to collect time domain signal data samples under each input and output channel generated by excitation in the magnetic suspension rotor system in real time;

[0176] A second module is configured to convert the time domain data samples obtained by the first module into frequency domain data, obtain a frequency response data set for model identification, and establish a frequency response model based on the frequency response data set;

[0177] A third module is configured to use an optimization algorithm to perform unbiased estimation on the frequency response model obtained by the second module, and establish an empirical transfer function estimation model;

[0178] And a fourth module is configured to evaluate the precision of the empirical transfer function estimation model obtained by the third module by using a normalized root mean square error, compare the predicted output of the empirical transfer function with the time domain signal data samples, and verify the accuracy and reliability of the model.

[0179] It should be understood that although the steps in the flowcharts of the accompanying drawings are shown in a sequential order following the arrows, the steps are not necessarily executed in the order indicated by the arrows. Unless explicitly stated otherwise herein, the execution of the steps is not strictly limited to the order indicated by the arrows, and can be executed in other orders. Moreover, at least some of the steps in the flowcharts of the accompanying drawings can include multiple sub-steps or multiple stages, which are not necessarily executed at the same time, but can be executed at different times, and the execution order is not necessarily sequential, but can be round-robin or alternating with at least some of the other steps or sub-steps or stages of other steps.

[0180] The above only describes some embodiments of the present application, and it should be pointed out that for those skilled in the art, without departing from the principles of the present application, a number of improvements and refinements can be made, which should also be considered as the protection scope of the present application.

Claims

1. A magnetic levitation rotor system model identification method based on a frequency response model, characterized in that: include: S100, using a harmonic signal to excite the magnetic levitation rotor system, and collecting time domain signal data samples of each input and output channel generated by the excitation in the magnetic levitation rotor system in real time; S200, converting the time domain signal data samples into frequency domain data to obtain a frequency response data set for model identification, and establishing a frequency response model based on the frequency response data set; S300, using an optimization algorithm to perform unbiased estimation on the frequency response model and establish an empirical transfer function estimation model; S400 , evaluating the accuracy of the empirical transfer function estimation model using normalized root mean square error, and comparing the predicted output of the model with the time domain signal data sample to verify the accuracy and reliability of the model.

2. The method for identifying a magnetic levitation rotor system model based on a frequency response model according to claim 1, characterized in that: In step S100, a step sine signal is selected as the excitation signal, and the sinusoidal displacement harmonic instruction of the step sine signal is The formula is: where t = cnt·T s , Where, is the harmonic instruction, t h and T s are preset values, representing the amplitude of the excitation signal, the duration of a single excitation, and the control step length respectively; and cnt are variables, the former is the sinusoidal signal frequency at the nth excitation, and the latter is the count value; the variable superscript k=1 represents the relevant signal or transmission function of the X degree of freedom of the radial magnetic bearing A at the DE end of the magnetic levitation rotor, marked with DE-X, the variable superscript k=2 represents the relevant signal or transmission function of the Y degree of freedom of the radial magnetic bearing A at the DE end of the magnetic levitation rotor, marked with DE-Y, the variable superscript k=3 represents the relevant signal or transmission function of the X degree of freedom of the radial magnetic bearing B at the NDE end of the magnetic levitation rotor, marked with NDE-X, the variable superscript k=4 represents the relevant signal or transmission function of the Y degree of freedom of the radial magnetic bearing B at the NDE end of the magnetic levitation rotor, marked with NDE-Y, and they are marked with DE-X, DE-Y, NDE-X and NDE-Y respectively.

3. The method for identifying a magnetic levitation rotor system model based on a frequency response model according to claim 2, characterized in that: The harmonic signal excitation in step S100 needs to be performed multiple times, and the magnetic bearing control host computer is used to complete the signal acquisition, which specifically includes: S101, excitation channel initialization, set k = 1, signal excitation will be applied to the control loop in the DE-X direction; S102, initializing the excitation frequency, setting n=1, the host computer sends the initial sinusoidal harmonic parameters to the k-channel control loop, and the controller generates sinusoidal harmonic instructions beat by beat according to the sinusoidal displacement harmonic instruction formula; S103: The four radial degrees of freedom of the magnetic levitation rotor system generate harmonic responses under the action of sinusoidal harmonic excitation, and the controller stores the control current instructions of the excitation channel in real time. and displacement voltage feedback signals on 4 feedback channels S104. After the response is completed, the controller sends the stored time domain signal data to the host computer for re-save, and the host computer updates and sends the excitation frequency To the controller; S105: The system enters the next harmonic response stage and repeats steps S102-S104 until the excitation frequency reaches the upper limit of 1000 Hz.

4. The method for identifying a magnetic levitation rotor system model based on a frequency response model according to claim 3, wherein: Step S200 includes: S201, cooperating with the magnetic bearing control host computer to analyze the time domain signal data and convert the time domain signal into a frequency domain signal; S202, set k = k + 1, enter the frequency response test of the next excitation channel, and complete the data acquisition of the 4×4 frequency response data matrix between the four input and output channels DE-X, NDE-X, DE-Y and NDE-Y.

5. The method for identifying a magnetic levitation rotor system model based on a frequency response model according to claim 4, characterized in that: The frequency response model in step S200 is: Where H frd represents the frequency response model, is the angular frequency of the nth frequency test point; N f is the total number of sweep points; It means the system's i-th input To the jth output Subsystem frequency response model.

6. The method for identifying a magnetic levitation rotor system model based on a frequency response model according to claim 5, characterized in that: Step S300 includes: S301, determining the structure of the empirical transfer function estimation model to simplify the complexity of model identification; S302, constructing a covariance-based weighted cost function to indicate the iteration direction of the parameter set; S303, selecting an initial parameter set and an initial damping coefficient, and setting an upper limit on the number of iterations and a convergence threshold; S304, calculating the Jacobian matrix of the weighted cost function; S305, updating the initial parameter set according to the Jacobian matrix, wherein the damping coefficient is adaptively adjusted as the iteration progresses to optimize the convergence speed; S306, iteratively update the parameters, and compare the difference of the cost function values ​​of adjacent iterations according to the updated parameter set to perform a convergence test. When the difference is less than a preset threshold, the model is considered to have converged, the algorithm is terminated, and the model identification result is output.

7. The method for identifying a magnetic levitation rotor system model based on a frequency response model according to claim 6, characterized in that: The structure of the empirical transfer function estimation model in step S301 is specifically written as: Where, represents the empirical transfer function estimation model; θ ij Represents the sub-transfer function The parameter variables to be estimated in ; Furthermore, in the model identification problem, The construction form is to set it to the standard rational polynomial form of the transfer function, that is: Where s η , s η-1 ,...,s represents different powers of S, which is a variable in the complex frequency domain; η is the order of the denominator polynomial; a η-1 , ..., a1, η0 are the coefficients of the denominator polynomial of the transfer function; β ξ , β ξ-1 , ..., β0 are the coefficients of the numerator polynomial of the transfer function.

8. The method for identifying a magnetic levitation rotor system model based on a frequency response model according to claim 7, characterized in that: In step S301, in order to simplify the complexity of model identification, the empirical transfer function Converted into the combined transfer function form of each order mode, that is: Where K gain is the gain coefficient of the subsystem model; p α , z β are the first-order poles and zeros of the system respectively; n α , n δ Indicates the modal order considered in the identification process; n β , n γ are the number of first-order and second-order links of the molecule, respectively; ζ Mγ represents the damping ratio of the second-order link of the molecule, ω Mγ represents the natural frequency of the second-order link of the molecule, ζ Mδ Represents the damping ratio of the denominator second-order link, ω Mδ Represents the natural frequency of the denominator second-order link.

9. The method for identifying a magnetic levitation rotor system model based on a frequency response model according to claim 8, characterized in that: After the power amplifier is equivalent to a first-order low-pass filter link and the displacement sensor is equivalent to a fixed gain, the empirical transfer function The final expression is: Among them, there are 24 parameters in the parameter set θ to be identified, and the specific expression is: Where K ij Represents the gain coefficient of the subsystem model; Indicates the pole, reflecting the bandwidth of the power amplifier; [p 1x , p 2x , p 1y , p 2y ] reflects the ratio of the displacement stiffness coefficient to the mass inertia of the magnetic bearing in each radial degree of freedom; Represents the zero point, introduced by the mechanical coupling effect.

10. The method for identifying a magnetic levitation rotor system model based on a frequency response model according to claim 6, characterized in that: Step S302 specifically includes: Define the element error matrix: Furthermore, an error cost function V(θ) is constructed to obtain the overall difference between the two models in the entire frequency test section, including: Where, At the frequency point The output error weight matrix of ; tr(·) represents the trace operator of the matrix; is the frequency weight coefficient, which is used to adjust the influence of the overall deviation of the model at different frequency points in the cost function calculation; N′ f is the number of frequency points to be calculated.

11. The method for identifying a magnetic levitation rotor system model based on a frequency response model according to claim 10, characterized in that: Step S304 includes: The Jacobian matrix J(θ) contains the partial derivatives of V(θ) with respect to each parameter in the parameter set θ. The elements J in the matrix are (q) The calculation formula is expressed as: Where E is the unit error matrix.

12. The method for identifying a magnetic levitation rotor system model based on a frequency response model according to claim 11, characterized in that: The parameter set update formula in step S305 is: i (k+1) =θ (k) -[J(θ (k) ) T J(θ (k) )+λ·diag(J(θ (k) ) T J(θ (k) ))] -1 ·J(θ (k) ) T E(θ (k) ) Where θ (k) and θ (k+1) are the parameter sets for the current iteration step and the next iteration step respectively; λ is the dynamically adjusted damping factor.

13. A method for identifying a magnetic levitation rotor system model based on a frequency response model according to any one of claims 10 to 12, characterized in that: The calculation formula of the normalized root mean square error in step S400 is: Where, stands for normalized root mean square error.

14. A magnetic levitation rotor system model identification system based on a frequency response model, characterized in that: The method for identifying a magnetic levitation rotor system model based on a frequency response model according to any one of claims 1 to 13 is implemented, comprising: The first module is used to collect time domain signal data samples of each input and output channel generated by the excitation in the magnetic levitation rotor system in real time; A second module is configured to convert the time domain signal data samples obtained by the first module into frequency domain data, obtain a frequency response data set for model identification, and establish a frequency response model based on the frequency response data set; The third module is used to use an optimization algorithm to perform unbiased estimation on the frequency response model obtained in the second module and establish an empirical transfer function estimation model; And a fourth module is used to evaluate the accuracy of the empirical transfer function estimation model obtained by the third module using the normalized root mean square error, and compare the predicted output of the empirical transfer function with the time domain signal data sample to verify the accuracy and reliability of the model.

Images