New energy system frequency domain modeling and equivalent circuit order reduction analysis method

Abstract

The invention discloses a new energy system frequency domain modeling and equivalent circuit order reduction analysis method, and relates to the technical field of new energy power system modeling and stability analysis. The method is used for solving the problems that in an existing frequency domain modeling means, the automation degree of data acquisition is low, fitting precision and stability are poor, equivalent circuit universality is poor, and order reduction analysis is not systematic. The method comprises the following steps: acquiring system frequency domain response data based on disturbance injection; constructing a system rational function model by adopting an improved vector fitting algorithm; and then model order reduction is realized through state space modeling and an SVD-based balanced truncation method, and response consistency and stability analysis is carried out in a frequency domain and a time domain. The method can be widely applied to modeling analysis of new energy equipment such as wind power equipment, voltage source converters and high-voltage direct-current power transmission equipment, and has the advantages of being high in automation degree, high in modeling precision, good in visualization, clear in structure, suitable for simulation and control design and the like.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of new energy power system modeling and stability analysis, in particular to a new energy system frequency domain modeling and equivalent circuit order reduction analysis method. BACKGROUND

[0002] With the further improvement of ecological environment, social policy and other aspects of new energy consumption requirements, the power grid development presents the characteristics of high proportion of new energy and high proportion of power electronics. The power system gradually evolves from the traditional "source-grid-load" structure to the AC-DC hybrid system containing a large number of power electronic devices. Among them, various types of new energy power generation equipment based on power electronic converters represented by wind power and photovoltaic, flexible DC transmission devices and new types of power electronic loads represented by electric vehicles are widely connected to the power grid and microgrid, forming a new type of power system.

[0003] However, with the increasing proportion of new energy generation, multi-machine, machine-grid interaction has frequently occurred wide frequency oscillation phenomenon. There are also related reports of wide frequency oscillation phenomenon from 1Hz to 1000Hz on the microgrid level. Therefore, establishing an accurate, analyzable and scalable new modeling and analysis tool has become a research hotspot and key requirement in the current power system engineering technology field.

[0004] In view of the new type of power system instability problem, many studies have carried out related work from various aspects. The current mainstream modeling method is mainly time domain simulation, which depends on the fine modeling of the controller structure and parameters, and has the problems of high computational complexity and poor model universality. As a more systematic and visualized means, frequency domain modeling can provide a clear characterization of the small signal stability, resonance characteristics and interconnection impedance characteristics of the system, and is gradually becoming an important means of new energy system modeling and analysis.

[0005] There are still the following problems in the existing frequency domain modeling technology: the automaticity of frequency domain data extraction is low, the traditional method often relies on manual configuration of disturbance signal and obtains system response through frequency-by-frequency simulation, lacks integrated control mechanism for disturbance injection, simulation control, data acquisition and post-processing, and cannot meet the batch modeling demand in engineering environment; the fitting algorithm lacks stability and precision control mechanism, high-dimensional frequency domain data often has complex multi-resonance and multi-pole characteristics, and traditional fitting methods (such as empirical fitting based on Bode diagram or single function fitting) often have overfitting, underfitting or numerical instability problems when dealing with multi-variable and wide-band response; the equivalent circuit model construction lacks universality and SPICE compatibility, the existing rational function-circuit model conversion method relies on manual matching of base function structure, and cannot automatically map the fitting result to a standardized, clear structure and visual circuit topology, which seriously limits the reusability of the model in electromagnetic simulation and control design; there is a lack of systematic order reduction and stability analysis mechanism, high-order rational function can improve the fitting accuracy, but there is a significant computational burden in real-time control or hardware-in-the-loop (HIL) simulation. The existing order reduction method cannot be effectively integrated with the modeling process, and there is also a lack of systematic analysis framework for the stability and fidelity (such as pole retention, frequency response consistency, etc.) of the reduced model; the overall process lacks modularity, automation and expansion capability, and there is currently no complete process tool chain that can automatically process the frequency domain disturbance injection, response data acquisition, rational function fitting, circuit modeling, model order reduction and verification in a closed loop.

[0006] In summary, in view of the increasingly complex modeling requirements in new energy power systems, a new frequency domain modeling and analysis method with full-process automatic processing capability, high-dimensional system compatibility and model order reduction and stability maintenance mechanism is urgently needed to improve model quality, accelerate engineering deployment and provide a solid foundation for system stability research and controller design. SUMMARY

[0007] The present application aims to solve the problems of low automaticity of frequency domain data extraction, lack of stability and precision control mechanism in fitting algorithm, lack of universality and SPICE compatibility in equivalent circuit model construction, lack of systematic order reduction and stability analysis mechanism, and lack of modularity, automation and expansion capability in the overall process of the existing frequency domain modeling technology, and proposes a new energy system frequency domain modeling and equivalent circuit order reduction analysis method. The present application can realize full-process automatic operation from disturbance injection, simulation control, spectrum analysis, function fitting, circuit equivalence, system order reduction to stability analysis, and is suitable for new energy power systems composed of typical power electronic devices such as wind power, voltage source converter (VSC) and high voltage direct current transmission (HVDC).

[0008] The object of the application can be achieved by the following technical solutions: a new energy system frequency domain modeling and equivalent circuit order reduction analysis method, comprising:

[0009] S1: a frequency sweeping method based on measured data or m file and simulink model interaction, automatically injecting disturbance signals at multiple frequency points to obtain frequency domain response data of the system;

[0010] S2: performing rational function approximation fitting on the frequency domain response data extracted in step S1 based on an improved vector fitting method, and simplifying high-order and non-analytic frequency domain response data into a low-order rational function model;

[0011] S3: converting the rational function model obtained in step S2 into a SPICE compatible equivalent circuit model composed of first-order and second-order equivalent circuits as basic components; then converting the rational function model into a state space model, using the balanced truncation order reduction method based on SVD to reduce the order of the state space model of the equivalent circuit, and simultaneously converting the reduced order model into an equivalent circuit;

[0012] S4: performing system stability and precision trade-off analysis based on the rational function models before and after model order reduction obtained in steps S2 and S3 and the equivalent circuits.

[0013] Compared with the prior art, the application has the following advantages:

[0014] 1. The application realizes full-process automatic modeling, realizes full-automatic process from disturbance injection, frequency domain response data extraction, vector fitting modeling, circuit equivalent conversion to state space order reduction and stability evaluation, greatly reduces manual operation, and improves modeling efficiency and engineering adaptability.

[0015] 2. The equivalent circuit has strong compatibility, the equivalent circuit construction method is based on standard RLC and controlled source components, can be converted into a circuit structure compatible with the SPICE simulation environment, has good visualization and reusability, and is convenient for subsequent electromagnetic simulation and control system debugging.

[0016] 3. The order reduction of the application is controllable and the stability is guaranteed, the model order reduction is realized by the state space balanced truncation method, the main dynamic characteristics are retained, the theoretical error upper bound is provided, and the stability and fidelity of the reduced order model are ensured through pole tracking, frequency domain and time domain consistency verification.

[0017] 4. The application has wide application range and is suitable for various new energy systems including wind power, voltage source converter (VSC), high voltage direct current transmission system (HVDC), etc., and provides efficient and reliable technical support for system modeling, simulation verification and controller design in actual engineering. BRIEF DESCRIPTION OF DRAWINGS

[0018] For the convenience of those skilled in the art to understand, the present application is further described below in conjunction with the drawings.

[0019] Figure 1 The structure diagram of embedding the simulink model for the system to be measured of the present application Figure 1 ;

[0020] Figure 2 The structure diagram of embedding the simulink model for the system to be measured of the present application Figure 2 ;

[0021] Figure 3 The internal structure diagram of the measurement data export module for the system to be measured of the present application

[0022] Figure 4 The internal structure diagram of the disturbance injection module of the present application

[0023] Figure 5 The flow chart of the whole process of the present application

[0024] Figure 6 The schematic diagram of the swept frequency waveform of the present application

[0025] Figure 7 The schematic diagram of the vector fitting effect of the present application

[0026] Figure 8 The structure diagram of three types of equivalent sub-circuit of the present application

[0027] Figure 9 The schematic diagram of the automatic identification and modeling of the equivalent circuit of the present application. DETAILED DESCRIPTION

[0028] The technical solutions of the present application will be described below in conjunction with the embodiments, obviously, the described embodiments are only some of the embodiments of the present application, not all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor belong to the scope of protection of the present application.

[0029] Please refer to Figures 1-3 The new energy system frequency domain modeling and equivalent circuit order reduction analysis method shown in the drawing, the equivalent power grid structure based on Matlab / Simulink is embedded in the measured power electronic system, the structure includes: the measured power electronic system 1 (including wind power system, inverter, SVG, etc.), three-phase voltage and current measurement module 2, disturbance injection module 3, equivalent infinite grid 4 and measurement data export module 5.

[0030] Please refer to Figure 4As shown, the internal part of the disturbance injection module 3 includes a sine wave generation module 31 and a controlled voltage source 32.

[0031] Please refer to Figures 5-9 As shown, a new energy system frequency domain modeling and equivalent circuit reduction analysis method, the whole process includes:

[0032] S1: Based on the measured data or m file and simulink model interaction of frequency sweep method, automatic disturbance signal injection at multiple frequency points, obtain the frequency domain response data of the system.

[0033] S2: Based on the improved vector fitting method for step S1 to extract the frequency domain response data of the rational function approximation fitting, the high order, non-analytic frequency domain response data is simplified to the low order (by pole p k , the residue c k And other parameters defined) rational function model.

[0034] S3: The system rational function obtained by fitting in step S2 is converted into a SPICE compatible equivalent circuit model composed of first and second order equivalent circuits (including controlled sources) as basic components. Then the rational function model is converted into a state space model, and the state space model of the equivalent circuit is reduced in model using the balanced truncation reduction method based on SVD method, and the reduced model is converted into an equivalent circuit.

[0035] S4: Based on the rational function model and equivalent circuit obtained in steps S2 and S3 before and after model reduction, the system stability and accuracy trade-off analysis is carried out.

[0036] In steps S1 and S2, the frequency domain response data of the power electronic system to be measured (such as wind power system, voltage source converter VSC, VSG, etc.) is obtained based on Python and MATLAB, and vector fitting is carried out. The frequency domain response data acquisition and fitting steps include:

[0037] S11: The frequency domain data of the system is obtained based on the engineering actual measurement data or the frequency sweep acquisition of the m file and the simulink model interaction;

[0038] S12: The frequency sweep model is based on embedding the power electronic system to be measured (such as wind power system, voltage source converter VSC, etc.) into equivalent grid environment, constructing complete electromagnetic transient simulation model on Simulink platform, and setting three-phase controlled voltage source interface which can inject disturbance;

[0039] S13: Control the simulink model through the corresponding m file, superimpose positive and negative sequence three-phase balanced disturbance on the system steady state operating point twice, the disturbance is injected in the form of sine wave, and set the amplitude and phase sequence of the disturbance to be injected, wherein the amplitude of the disturbance signal is small to ensure that the system operates in the small signal linear region;

[0040] S14: Set the frequency range [ω min , ω max ] and the frequency distribution interval Δω, the m file automatically controls the simulation to operate at multiple frequency points f k , and extracts the steady state response data of the system. The injection mode of each frequency point or multiple frequency points in parallel can be selected, and after the simulation is completed, the oscilloscope exports the dq axis voltage and current waveform data under the positive and negative sequence disturbance to the Matlab workspace for subsequent construction of the small signal impedance matrix of the system;

[0041] S15: Since the simulation obtains a time domain discrete sampling signal, in order to efficiently extract the amplitude and phase characteristics at the disturbance frequency f k , the present application adopts the fast Fourier transform (FFT, Fast Fourier Transform) method for frequency spectrum analysis. Taking the current signal as an example, the sampling signal is:

[0042] i[n] = i(nΔt), n = 0, 1,..., N-1

[0043] Where Δt is the sampling interval, N is the total number of sampling points, and the corresponding sampling frequency is

[0044] The FFT transform formula is:

[0045]

[0046] Where: I[k] represents the complex spectrum value of the kth frequency point (corresponding to the frequency );

[0047] The amplitude is:

[0048]

[0049] The phase is:

[0050]

[0051] In the frequency sweep simulation, only the I[k*] and V[k*] corresponding to the frequency point k* closest to the disturbance frequency in the two sets of FFT outputs (positive sequence disturbance, negative sequence disturbance) are extracted, and the steady state voltage and current on the complex plane at the frequency can be obtained.

[0052] In the dq-axis small-signal analysis, the impedance matrix Z dq is a 2x2 matrix, containing four components: Z dd (d-axis self-impedance), Z dq (d-axis to q-axis mutual impedance), Z qd (q-axis to d-axis mutual impedance), and Z qq (q-axis self-impedance). Thus the dq-axis small-signal impedance matrix is:

[0053]

[0054] The linear relationship in the system can be considered as:

[0055] V = -Z dq ·I

[0056] The present application adopts injecting two groups of linearly independent voltage disturbances for frequency scanning, and respectively injects positive sequence voltage disturbance and negative sequence voltage disturbance into the system, the frequency of the positive sequence disturbance is 2xpi x h 注入 + 2xpi x f0, the frequency of the negative sequence disturbance is abs(2xpi x h 注入 - 2xpi x f0), when the disturbance frequency is greater than the fundamental frequency, the disturbance phase sequence is ABC, when the disturbance frequency is less than the fundamental frequency, the disturbance phase sequence is ACB, and then the experimental data are collected to carry out calculation:

[0057] For the primary disturbance test, we have:

[0058]

[0059] or

[0060]

[0061] The dq-axis voltage and current signals after the two disturbances are collected, and are represented by a 2x2 matrix:

[0062] The voltage matrix V dq is:

[0063]

[0064] Wherein:

[0065] x is the sweep point.

[0066] The current matrix I dq is:

[0067]

[0068] Wherein:

[0069] x is the sweep point.

[0070] Voltage matrix V dq and current matrix I dq The data of each sweep point is analyzed by the Fast Fourier Transform (FFT) method.

[0071] Calculate the inverse matrix of the current matrix

[0072]

[0073] Calculate Finally, the complete impedance matrix can be constructed:

[0074]

[0075] S21: Extract the system impedance data information obtained in S15 as a csv file format, and use the VectorFitting program based on Python to process and construct the rational function model in the following form:

[0076] where s = σ + jω.

[0077] To achieve the desired fitting effect, the model function should match the given frequency response at each frequency sample (ω n ) of the model function:

[0078] F s = jω k ) = 0 F sampled (ω k )

[0079] S22: Select the initial poles according to the distribution range of the frequency samples:

[0080] a n = -α + jβ, a n+1 = -α - jβ

[0081] where β is linearly distributed in the frequency range; α = β / 100 is very small in real part, which ensures that the poles are located in the left half of the complex plane, but does not cause numerical instability. When the initial poles are complex numbers, in order to ensure that the final fitting result is a real function, the zero points obtained by fitting need to be conjugate pairs.

[0082] 1. There is a pair of conjugate poles:

[0083] ​

[0084] 2. Construct a new real matrix system, preserving the conjugate relationship by the following block matrix form:

[0085]

[0086] 3. Construct a new system matrix:

[0087]

[0088] 4. Solve the eigenvalues of H, the result is a conjugate pair.

[0089] S23: Fit the frequency response data f(s k ) or f(jω k ) using a rational function approximation model:

[0090]

[0091] Introduce an auxiliary function σ(s) to linearize the problem of solving the nonlinear rational function to be fitted, and specify:

[0092]

[0093]

[0094] where is a set of known poles selected initially; is the remainder corresponding to , which is unknown, so each term linearly depends on

[0095] Then construct an approximate relationship, which can be obtained:

[0096]

[0097] S24: For the approximate relationship obtained in S23, there are K frequency points (i.e. s1, s2, …, s K ), substitute each frequency point into the relationship. For the i-th frequency point s i , we write:

[0098]

[0099] The above formula can also be rewritten as:

[0100]

[0101] Establish this relationship for multiple frequency points to form an overdetermined linear equation system (the number of observation points > the number of unknowns), and construct a matrix form linear system:

[0102] Ax = b

[0103] After simplification, we can obtain a linear expression for all unknown parameters x:

[0104]

[0105] There are a total of 2N+2 unknowns.

[0106] For each frequency point s i Construct a row, and the row vector is represented as:

[0107]

[0108] Each A i It is a single-row (2N+2) dimensional vector.

[0109] The right-hand vector b is:

[0110] b i =f(s) i )

[0111] Solving for x using the least squares method yields all the residual coefficients c. n , With constant terms d and h.

[0112] S25: Based on the σ(s)f(s) and σ(s) obtained in S24, the approximate relationship can be converted into a fitted form:

[0113] (σf) fit (s)≈σ fit (s)f(s)

[0114] Human-made settings:

[0115] The fitting result of σf (σf) fit (s) has zeros (z1, z2, ..., z). N+1 );

[0116] σ fit (s) has a zero point

[0117] They all use the same starting pole. As the denominator (i.e., the initial poles remain unchanged).

[0118] Therefore, the fitting function (σf) can be... fit (s), σ fit The expression for (s) can be rewritten as a fractional product:

[0119]

[0120] Divide both sides of the product by σ.fit (s), so that a new fitting expression is obtained:

[0121]

[0122] The specific expression is:

[0123]

[0124] The poles of f(s) are changed from a n to The zeros are changed from c n to z n , and h is a proportional coefficient. Sometimes the new poles will fall in the right half of the complex plane (unstable region), in which case the pole reflection is used, only the real part is taken as negative, moved to the left half plane, and the system remains stable, that is, if is replaced by

[0125] After updating the poles, we will get the new poles obtained in the last round as Then repeat the above fitting process until the fitting error is small enough, that is, the poles change is no longer significant.

[0126] Using equivalent circuit conversion and SVD-based balanced truncation reduction method, the fitted system is simplified for analysis, including:

[0127] S31: Based on the rational function of the frequency response data fitted out in S25, it is identified and divided into three categories: constant term and proportional term, real pole and residue, and complex conjugate pole and residue. As shown in Figure 8 (a) is a real pole and residue RL series circuit, (b) is a constant term and proportional term RC parallel circuit, and (c) is a complex conjugate pole and residue series RLC circuit with VCCS.

[0128] S32: Convert the three types of parts in S31 into equivalent sub-circuits respectively, all of which are represented by equivalent admittance Y ( s ).

[0129] 1. The constant term and proportional term (d e ) are constructed by a parallel RC circuit.

[0130] The target response of the constant term and proportional term is:

[0131] H target ( s ) = d i + s e i

[0132] The admittance of the parallel RC circuit is:

[0133]

[0134] Let and C = e i , so that the admittance matches the target response.

[0135] 2. Real pole and residue term Constructed from a series RL circuit.

[0136] The target response of the real pole and residue term is:

[0137]

[0138] The admittance of the series RL circuit is:

[0139]

[0140] Let and , so that the admittance matches the target response.

[0141] 3. Complex conjugate pole and residue term Constructed from a series RLC circuit with controlled source.

[0142] In vector fitting, complex poles ( p k = p' k + jp" k ) and residues ( c k = c' k + jc" k ) always come in complex conjugate pairs and . Therefore, the target response of the equivalent circuit of a complex conjugate pair is:

[0143]

[0144] The admittance of the series RLC circuit with controlled source is:

[0145] A series RLC circuit combined with a parallel current source controlled by the voltage across the capacitor provides an admittance Y RCL,I ( s ), which can be matched to the target response:

[0146]

[0147] Let and b = -2(c' k p' k + c" k p" )k ),in To match the admittance with the target response.

[0148] S33: Based on the classification of rational functions and the target corresponding matching formula in step S32, the equivalent circuit is automatically modeled and the circuit element parameters are calculated using the CircuitDrawer class written in Python, generating a SPICE-compatible equivalent circuit model.

[0149] S34: Convert the rational function of the frequency response data fitted in step S25 into a state-space equation model. Specifically, convert each term... Transforming it into a first-order dynamic subsystem, we construct a state-space model of the following form:

[0150] y(t)=Cx(t)+Du(t)

[0151] Wherein, the state matrix A∈R n×n From the pole p k Construct it as a diagonal matrix, with the input matrix B∈R n×m Generally, it is formed by stacking all-one column vectors or by setting the output matrix C∈R based on the actual residual normalization of the system. p×n From the residual matrix c k The output direct-pass term D is taken as d. For the case containing the linear term se, it is treated as a derivative channel, corresponding to adding an additional first-order integral subsystem to the state space. Its output item is ex e This expands the dimension of the entire state space.

[0152] S35: After completing the state-space model construction, apply Gramians-based equilibrium truncation to the system. First, calculate the system's controllability using GramianW. c and observability GramianW o Their definitions are as follows:

[0153]

[0154] This step requires the system to be stable, i.e., all eigenvalues ​​Re(λ(A)) < 0. The resulting Gramsians are then subjected to Cholesky decomposition, as follows:

[0155]

[0156] Then the intermediate matrix Applying the Singular Value Decomposition (SVD) operation, we obtain:

[0157]

[0158] where ∑ = diag(σ1,...,σ n ) is the Hankel singular value diagonal matrix of the system, representing the importance of each state to the overall input and output response of the system. Further, construct the balanced transformation matrix:

[0159] T = R c V∑ -1 / 2 ,

[0160] This applies a similarity transformation to the system, converting it to the "balanced coordinate system", where the system matrix is transformed to:

[0161]

[0162] At this point, the system states are sorted in the new coordinate system according to the size of the Hankel singular value, with optimal truncation properties. If the order of the system needs to be reduced from the original n to r, only the state components corresponding to the first r principal singular value directions can be retained, and the corresponding reduced-order system can be represented as:

[0163]

[0164] The system error introduced after reduction has a theoretical bound control. In the sense, the error between the reduced-order system and the original system satisfies the following upper bound:

[0165]

[0166] This error bound has adjustability, and the truncation order r can be flexibly determined according to the singular value spectrum variation law, effectively reducing the system complexity and improving the simulation efficiency and controller design feasibility while preserving the key dynamic behavior.

[0167] Based on the rational function models and equivalent circuits obtained before and after the model reduction of S2 and S3, the system stability and accuracy trade-off analysis is carried out, including:

[0168] S41: Based on the rational function models and equivalent circuits obtained before and after the model reduction of S2 and S3, the system stability and accuracy trade-off analysis is carried out.

[0169] S42: Pole position comparison and analysis:

[0170] Let the state space model of the original system be:

[0171] y(t) = Cx(t) + Du(t)

[0172] The system transfer function is:

[0173] G(s) = C(sI-A)-1 B+D

[0174] After reducing it to an r-order system through S35, the reduced model is obtained:

[0175] yr(t)=C r x r (t)+D r u(t)

[0176] Calculate the pole set {λ i} of the two systems again. Analyze whether it meets:

[0177]

[0178] Ensure that the reduced system is still stable and the main poles have not shifted significantly.

[0179] S43: Frequency domain response consistency verification

[0180] Draw the Bode diagram and Nyquist curve of the system before and after reduction, and analyze the consistency of the key frequency band (such as the resonance frequency f r , cutoff frequency f c ).

[0181] For transfer function G(s), its frequency response is:

[0182] G(jω)=C(jωI-A) -1 B+D

[0183] Pay attention to the following error indicators:

[0184] |G(jω)-G r (jω)|≤ε , ω∈[ω min ,ω max ]

[0185] Where ε is the frequency response tolerance error. The goal is to maintain the error within ±1dB in the target frequency band.

[0186] S44: Time domain response simulation verification:

[0187] Apply a unit step signal u(t) = 1(t) to the system, and analyze the output response y(t), y r (t) of the original model and the reduced model.

[0188] Compare the following indicators: rise time t r ; overshoot M p ; steady-state value deviation:

[0189] Δy ∞= | y(∞) - y r (∞) | ; adjustment time t s ;

[0190] Desired to keep the error:

[0191]

[0192] Where δ is the upper bound of the allowable error, T is the simulation period.

[0193] S45: stability and precision trade-off analysis:

[0194] In the process of order reduction, the states corresponding to the first r dominant Hankel singular values are retained, and the order is reduced from n to r. We use the error bound formula of balanced truncation method to evaluate the accuracy of order reduction:

[0195]

[0196] Where σ i is the Hankel singular value of the system. Through the inequality, the maximum error upper bound under different order reduction depths is quantified, and a trade-off relationship between system accuracy and simplification degree is established.

[0197] The preferred embodiments of the application disclosed above are only used to help illustrate the application. The preferred embodiments do not describe all the details, nor limit the application to the specific embodiments. Obviously, many modifications and variations can be made according to the content of the present application. The present application is selected and described in detail, in order to better explain the principles and practical application of the application, so that those skilled in the art can well understand and utilize the application. The application is limited only by the claims and their full scope and equivalents.

Claims

1. A method for frequency domain modeling and reduced order analysis of equivalent circuit of new energy system, characterized in that, Comprise: S1: based on the measured data or m file and simulink model interaction of frequency sweep method, automatic disturbance signal injection on multiple frequency points, obtain the frequency domain response data of system; S2: based on the improved vector fitting method for step S1 extracted frequency domain response data for rational function approximation fitting, the high order, non-analytic frequency domain response data is simplified to low order rational function model; S3: the rational function model obtained by fitting in step S2 is converted into a SPICE compatible equivalent circuit model composed of first and second order equivalent circuits as basic components; the rational function model is converted into a state space model, and the state space model of the equivalent circuit is reduced in model using the balanced truncation reduction method based on SVD method, and the reduced model is converted into an equivalent circuit; S4: based on the rational function model and equivalent circuit obtained before and after model reduction in step S2 and step S3, the system stability and precision trade-off analysis is carried out.

2. The method of claim 1, wherein, The frequency domain response data of system obtained in step S1 is based on engineering actual measurement data or m file and simulink model interaction of sweep acquisition; the sweep basis is to embed the power electronic system to be tested into the equivalent grid environment, build a complete electromagnetic transient simulation model on the Simulink platform, and then connect the sweep program to the grid point; The process of sweep in step S1 is: through the corresponding m file, the simulink model is controlled, the positive and negative sequence three-phase balanced disturbance is superimposed on the system steady state operating point twice, the disturbance is injected in the form of sine wave, and the amplitude and phase sequence of the disturbance to be injected are set, wherein the amplitude of the disturbance signal is small to ensure that the system operates in the small signal linear region; the frequency range and frequency distribution interval are set, the m file automatically controls the simulation to run on multiple frequency points, and the steady state response data of the system is extracted; After the simulation is completed, the oscilloscope exports the dq axis voltage and current waveform data of the measurement point under positive and negative sequence disturbance to the Matlab workspace, and the fast Fourier transform method is used to analyze the frequency spectrum of the time domain discrete sampling signal, and the small signal impedance matrix corresponding to each frequency point is solved by the response of two groups of linearly independent disturbances.

3. The method of claim 2, wherein, In the sweep simulation, the complex spectrum value corresponding to the frequency point closest to the disturbance frequency in the fast Fourier transform formula output is extracted, and the dq axis steady state response current or voltage at the frequency is obtained on the complex plane; the four components of the small signal impedance matrix are obtained by solving the nonlinear equation set, and the complete impedance matrix is finally constructed.

4. The method of claim 3, wherein, The construction process of rational function model in step S2 is: S21: extract the system impedance data information as a csv file format, use the VectorFitting program to process, and construct a rational function model, which matches the given frequency response at each frequency sample; S22: Select the initial pole according to the distribution range of the frequency sample; when the initial pole is a complex number, make the fitted zero point appear in pairs of conjugate; set a pair of conjugate poles, construct a new real matrix system, preserve the conjugate relationship through block matrix form, construct a new system matrix, and solve the eigenvalue of H, the result is a conjugate pair; S23: Fit the frequency response data using the rational function approximation model, introduce an auxiliary function to linearize the problem of solving the nonlinear rational function to be fitted; then construct an approximate relationship; S24: Set K frequency points for the approximate relationship obtained in step S23, and substitute each frequency point into the relationship; establish the relationship for multiple frequency points to form an over-determined linear equation system, and construct a matrix form linear system; After sorting, a linear expression about all unknown parameters x is obtained, which has 2N+2 unknowns; For each frequency point, a row is constructed, and the row vector is represented; the least squares method is used to solve the unknown parameters x, that is, all residual coefficients and constant terms are obtained; S25: Convert the approximate relationship to the fitting form based on the auxiliary function, rewrite the expression of the fitting function, and obtain a new fitting expression in the form of a fraction product, so that the new pole may fall in the right half of the complex plane, in which case the pole reflection is used, only the real part is taken to the left half plane, and the system stability is maintained; After updating the pole, the new pole obtained in the last round is used as the starting pole in the next round, and then the above fitting process is repeated until the fitting error is small enough, that is, the pole change is no longer significant.

5. The method of claim 4, wherein, Using equivalent circuit conversion and SVD-based balanced truncation reduction method, the fitted system is simplified and analyzed, including: S31: Based on the rational function of the frequency response data fitted in step S25, it is identified and divided into constant term and proportional term, real pole and residue, complex conjugate pole and residue three categories; S32: Convert the three parts in step S31 into equivalent sub-circuits respectively, all in the form of equivalent admittance, the constant term and the proportional term are constructed by a parallel RC circuit; the pole and the residue term are constructed by a series RLC circuit; the complex conjugate pole and the residue term are constructed by a series RLC circuit with a controlled source; in vector fitting, complex poles and residues always appear in complex conjugate pairs; the admittance of a series RLC circuit with a controlled source; a series RLC circuit combined with a parallel current source controlled by the voltage on the capacitor, which provides an admittance matching the target response; S33: Based on the classification of the rational function in step S32 and the target response matching formula, the CircuitDrawer class based on Python is used to automatically model the equivalent circuit and calculate the circuit element parameters, and a SPICE compatible equivalent circuit model is generated; S34: Convert the rational function of the frequency response data fitted in step S25 into a state space equation model; S35: After the completion of the state space model construction, the system is subjected to balanced truncation processing based on Gramians; first, the controllability Gramian W c and the observability Gramian W o are calculated, all eigenvalues are less than zero, then the resulting Gramians are subjected to Cholesky decomposition, and then singular value decomposition operation is applied to the intermediate matrix; the balanced transformation matrix is constructed to subject the system to similar transformation and convert it to the balanced coordinate system, at this time the system state is sorted in the new coordinate system according to the Hankel singular value size; if the order of the system needs to be reduced from the original n order to r order, only the state components in the direction of the corresponding first r dominant singular values are retained, the system error introduced after order reduction has a theoretical limit control, and the error between the reduced order system and the original system satisfies the following upper bound, the error limit has adjustability, and the truncation order r is flexibly determined according to the singular value spectrum variation law.

6. The method of claim 5, wherein, Based on the rational function model and the equivalent circuit of the model before and after reduction obtained in steps S2 and S3, the system stability and precision trade-off analysis is performed, including: Pole position comparison analysis: The state space model and transfer function of the original system are set, and the reduced order model is obtained after the order reduction to r order system through step S35. Then, the pole sets of the two systems are calculated; Frequency domain response consistency verification: The Bode diagram and Nyquist curve of the system before and after order reduction are drawn respectively, and the consistency of the key frequency band is analyzed, including the resonance frequency and the cutoff frequency; Time domain response simulation verification: The unit step signal is applied to the system, and the output response of the original model and the reduced order model is analyzed. The rise time, overshoot, steady state value deviation, adjustment time and expected holding error indicators are compared; Stability and precision trade-off analysis: In the process of order reduction, the states corresponding to the first r main Hankel singular values are retained, the order is reduced from n to r, the error bound formula of balanced truncation method is used to evaluate the order reduction accuracy, the maximum error upper limit under different order reduction depth is quantified, and the trade-off relationship between system accuracy and simplification degree is established.

Images