Complex Frequency Domain Admittance Matrix Modeling Method for Industrial Parks Based on Vector Fitting Technology

Abstract

The present invention provides a modeling method for the complex frequency domain admittance matrix of an industrial park based on the vector fitting technology. The method includes the following steps: S101. Measure the real-time amplitude-frequency and phase-frequency characteristic curves of each new energy unit under the grid-connected operation state of the industrial park by using the method of disturbing the grid connection point; S102. Obtain the number of peaks of the new energy unit, determine the order of the complex frequency domain impedance of the new energy unit, and further determine the order of vector fitting; S103. Use the vector fitting technology to segmentally identify the complex frequency domain impedance model of each new energy unit to obtain a segmented accurate parameter representation of the complex frequency domain impedance model; S104. Based on the complex frequency domain impedance models of each new energy unit, combine the line topology information and line parameters in the industrial park to establish the complex frequency domain admittance matrix of the industrial park. The present invention effectively solves the problems of "gray box" and "black box" of the converter caused by the confidentiality of manufacturers, and obtains a specific and accurate complex frequency domain admittance model of the converter to analyze the grid connection stability problem of the converter.

Description

Technical Field

[0001] The present invention relates to the technical field of vector fitting, and particularly to a method for modeling the complex frequency domain admittance matrix of an industrial park based on vector fitting technology. Background Art

[0002] The complex frequency domain admittance matrix is often used to analyze the resonance stability problem of power systems. By solving the roots where the determinant of the complex frequency domain admittance matrix is zero, the stability of the system can be judged according to the positive or negative real part of the characteristic roots. The parameters of traditional equipment such as wires and transformers are easy to know, but for the equipment of new energy power generation systems, due to the confidentiality of the internal control structure and parameters of the equipment by suppliers, it is difficult to know the specific parameters, which brings difficulties to the establishment of the complex frequency domain admittance matrix.

[0003] At present, the methods for establishing the complex frequency domain admittance matrix include two strategies: directly establishing based on the topological structure and deriving based on the state space equation. The complex frequency domain admittance matrix is constructed according to the construction theory of the complex frequency domain admittance matrix when the system topology and various parameters are clear in the direct establishment; the state space equation derivation method also requires clear system parameters, constructs the state space equation in the time domain and then transforms it into the complex frequency domain admittance matrix. However, at present, both strategies need to be applied when the parameters and topological structure of each device in the system are clear, and are no longer applicable when the internal control structure and parameters of new energy devices are not clear. Summary of the Invention

[0004] In view of this, the purpose of the present invention is to provide a method for modeling the complex frequency domain admittance matrix of an industrial park based on vector fitting technology to solve or at least partially solve the above problems existing in the prior art.

[0005] To achieve the above purpose, the present invention provides a method for modeling the complex frequency domain admittance matrix of an industrial park based on vector fitting technology, and the method includes the following steps:

[0006] S101. Measure the real-time amplitude-frequency and phase-frequency characteristic curves of each new energy unit under the grid-connected operation state of the industrial park by using the method of disturbing the grid connection point, and obtain the frequency points of the impedance amplitude and phase of each new energy unit;

[0007] S102. Based on the frequency points of the impedance amplitude and phase of each new energy unit, obtain the number of peaks of each new energy unit, determine the order of the complex frequency domain impedance of the new energy unit according to the number of peaks of each new energy unit, and further determine the order of vector fitting;

[0008] S103. Based on the order of vector fitting, use vector fitting technology to segmentally identify the complex frequency domain impedance model of each new energy unit, and obtain the segmented accurate parameters of the complex frequency domain impedance model;

[0009] S104. Based on the complex frequency domain impedance models of each new energy unit, combined with the line topology information and line parameters in the industrial park, establish the complex frequency domain admittance matrix of the industrial park.

[0010] Further, the step S101 specifically includes the following steps:

[0011] S11. Collect the three-phase voltage and current operation data of each new energy unit, use a feedforward neural network to fit the three-phase voltage and current operation data of each new energy unit, and establish a digital twin model;

[0012] S12. Apply small-signal voltage disturbances with different disturbance frequencies at the grid connection point of each new energy unit where is the frequency, is the th frequency point;

[0013] S13. Calculate the difference between the digital twin model and the grid-connected current after the small-signal voltage disturbance of each new energy unit to obtain the small-signal current response, and perform Fourier analysis on it to obtain the current response of each new energy unit under voltage disturbance in the complex frequency domain, which is expressed as follows:

[0014]

[0015] where is the frequency domain signal, is the discrete Fourier analysis, is the function for performing Fourier analysis, that is, the small-signal current response, is the sampling point sequence number, is the number of sampling points, is the natural logarithm, is the imaginary unit, is a natural number;

[0016] S14. Based on the Fourier analysis results, calculate the ratio of the small-signal current response at the disturbance frequency to the small-signal voltage disturbance intensity to obtain the corresponding complex frequency domain impedance ;

[0017] S15. Modify the disturbance frequency to obtain the amplitude-frequency and phase-frequency characteristic curves of the impedance Z of the new energy unit;

[0018] S16. When each new energy unit is operating, inject small-signal voltage disturbances in real-time and cyclically, obtain the complex frequency domain impedance , and judge whether to update the impedance Z of the new energy unit according to the data of the complex frequency domain impedance .

[0019] Further, step S102 includes: plotting Bode plots for each frequency point obtained in step S101, and determining the order of fitting to be selected according to the number of peaks in the Bode plots. .

[0020] Further, step S103 specifically includes the following steps:

[0021] S31: Segment the transfer function of the complex frequency-domain impedance of each new energy unit.

[0022] S32: Define the form of the transfer function of the complex frequency-domain impedance of the new energy unit in the frequency band, and perform vector fitting based on the form of the transfer function of the complex frequency-domain impedance of the new energy unit, which is expressed as follows:

[0023]

[0024] where, is the transfer function of the complex frequency-domain impedance of the new energy unit, is the order of each segment in the piecewise fitting, is the th pole, , , , are the parameters to be obtained through vector fitting respectively, is the independent variable;

[0025] S33: Based on the vector fitting results, integrate the transfer functions of the complex frequency-domain impedance of the new energy unit in frequency bands to obtain the piecewise transfer function of the complex frequency-domain impedance of the new energy unit as , which is expressed as follows:

[0026]

[0027] where, , are the transfer functions of the first segment and the th segment after vector fitting respectively, is the transfer function of the complex frequency-domain impedance of the new energy unit, is the frequency, and are the complex frequency-domain impedances the th and the th peak frequencies corresponding to on the Bode plot;

[0028] S34: Without segmentation, repeat step S32 to obtain the entire-segment transfer function of the complex frequency-domain impedance of the new energy unit as , the piecewise transfer function and the whole transfer function are weighted to obtain the final transfer function of the complex frequency domain impedance of each new energy unit as , which is expressed as follows:

[0029]

[0030] Wherein, and are the weight coefficients respectively.

[0031] Furthermore, vector fitting is performed based on the form of the transfer function of the complex frequency domain impedance of the new energy unit, which specifically includes the following steps:

[0032] S3.1. Uniformly select N initial poles from the fitting frequency band, and define an unknown function , which satisfies that after multiplying it with the fitting transfer function, it has the same poles as , and is expressed as follows:

[0033]

[0034] Wherein, is the fitting transfer function, is 's conjugate value;

[0035] S3.2. Substitute the frequency points in step S101 to solve the vector fitting linear equation system formula, which is expressed as follows:

[0036]

[0037]

[0038]

[0039]

[0040] Wherein, is the coefficient matrix in the vector fitting algorithm, is the number of frequency sampling points, is the unknown vector of the vector fitting equation, is the right side term of the vector fitting algorithm, is the complex frequency variable, is the frequency response value of the objective function, is the th residue of the pole, is the th conjugate value of the residue of the pole, is the transpose operation on the matrix;

[0041] S3.3. Based on steps S3.1 and S3.2, the fitted transfer function is obtained and expressed as follows:

[0042]

[0043] Wherein, is the fitted transfer function, and are respectively and zeros of The zero of is the pole of the fitted transfer function and can be obtained by solving the eigenvalue ;

[0044] S3.4. Repeat the calculation process of steps S3.2 - S3.3 to gradually obtain a more accurate vector fitting result.

[0045] Furthermore, the complex frequency domain admittance matrix of the industrial park includes:

[0046] Establish multiple nodes according to the topological structure of the industrial park. The diagonal element of each node is the algebraic sum of the complex frequency domain admittances of all components connected to this node;

[0047] The off - diagonal element represents the negative value of the complex frequency domain admittance between two nodes;

[0048] Each element of the complex frequency domain admittance matrix is symmetric about the main diagonal, that is = .

[0049] Compared with the prior art, the beneficial effects of the present invention are:

[0050] The present invention proposes a method for modeling the complex frequency domain admittance matrix of an industrial park based on vector fitting technology. By measuring the real - time amplitude - frequency and phase - frequency characteristic curves of each new energy unit under the grid - connected operation state of the industrial park, determining the order of the complex frequency domain impedance of the new energy unit, and then determining the order of vector fitting, using vector fitting technology to segmentally identify the complex frequency domain impedance model of each new energy unit, establishing the complex frequency domain admittance matrix of the industrial park, effectively solving the problems of "gray box" and "black box" of the converter caused by the confidentiality of manufacturers, obtaining a specific and accurate complex frequency domain admittance model of the converter to analyze the grid - connection stability problem of the converter. Brief Description of the Drawings

[0051] To more clearly illustrate the technical solutions in the embodiments of the present invention, the following briefly introduces the drawings required for the description of the embodiments. Obviously, the drawings in the following description are only the preferred embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can be obtained based on these drawings.

[0052] Figure 1 Schematic diagram of the process of the complex frequency domain admittance matrix modeling method provided by the embodiment of the present invention;

[0053] Figure 2 Topological diagram of a four-node industrial park provided by the embodiment of the present invention;

[0054] Figure 3 Schematic diagram of the energy storage side control structure of the photovoltaic power generation equipment provided by the embodiment of the present invention;

[0055] Figure 4 Bode diagram of the vector fitting result of the photovoltaic power generation equipment provided by the embodiment of the present invention. Detailed implementation manners

[0056] The following describes the principles and features of the present invention with reference to the accompanying drawings. The listed embodiments are only used to explain the present invention and are not used to limit the scope of the present invention.

[0057] Refer to Figure 1 , this embodiment provides a method for modeling the complex frequency domain admittance matrix of an industrial park based on vector fitting technology. The method includes the following steps:

[0058] S101. Measure the real-time amplitude-frequency and phase-frequency characteristic curves of each new energy unit under the grid-connected operation state of the industrial park by using the method of disturbing the grid connection point, and obtain the frequency points of the impedance amplitude and phase of each new energy unit, which specifically includes the following steps:

[0059] S11. Collect the three-phase voltage and current operation data of each new energy unit, and use a feedforward neural network to fit the three-phase voltage and current operation data of each new energy unit to establish a digital twin model.

[0060] S12. Apply a small-signal voltage disturbance with a disturbance frequency of at the grid connection point to each new energy unit , where ranges from 1 to 4000, corresponding to being 1 to 4000 Hz.

[0061] S13. Calculate the difference between the digital twin model and the grid-connected current of each new energy unit after the small-signal voltage disturbance to obtain the small-signal current response , and perform Fourier analysis on it to obtain the current responses of each new energy unit under voltage disturbance in the complex frequency domain, which are expressed as follows:

[0062] (1)

[0063] Among them, is the frequency-domain signal, is the discrete Fourier transform, is the function for performing Fourier analysis, that is, the small-signal current response , is the sampling point number, is the number of sampling points, is the natural logarithm, is the imaginary unit, is the natural number.

[0064] S14. Based on the Fourier analysis results, calculate the ratio of the small-signal current response to the small-signal voltage disturbance intensity at the disturbance frequency to obtain the corresponding complex frequency domain impedance .

[0065] S15. Modify the disturbance frequency to obtain the amplitude-frequency and phase-frequency characteristic curves of the impedance Z of the new energy unit in the range of 1 - 4000 Hz.

[0066] S16. When each new energy unit is operating, continuously inject a small-signal voltage disturbance with a frequency disturbance of 1 - 4000 Hz in real-time to obtain the complex frequency domain impedance , and determine whether to update the impedance Z of the new energy unit according to the data of the complex frequency domain impedance ; taking the disturbance frequency as an example, for 1 < <4000, if , then update to ; otherwise, keep unchanged and continue to inject disturbances at the next frequency point, where is the frequency.

[0067] S102. Based on the frequency points of the impedance amplitude and phase of each new energy unit, obtain the number of peaks of each new energy unit, determine the order of the complex frequency domain impedance of the new energy unit according to the number of peaks of each new energy unit, and then determine the order of vector fitting, specifically including:

[0068] Plot the Bode diagram for each frequency point obtained in step S101. According to the impedance part in the Bode diagram, observe its peaks and determine the order of fitting to be selected according to the number of peaks in the Bode diagram , where the The frequency at which the peak is located is .

[0069] S103. Based on the order of vector fitting, use the vector fitting technology to segmentally identify the complex frequency domain impedance model of each new energy unit, and obtain the accurate segmented parameters of the complex frequency domain impedance model, which specifically includes the following steps:

[0070] S31. Segment the transfer function of the complex frequency domain impedance of each new energy unit. If the fitted order = 1, no segmentation is performed; if > 1, the first segment is , the th segment is , the th segment is , where is the frequency corresponding to the th peak of the Bode plot of the complex frequency domain impedance curve, is the vector fitting order.

[0071] S32. Define the transfer function form of the complex frequency domain impedance of the new energy unit in the frequency band, and perform vector fitting based on the transfer function form of the complex frequency domain impedance of the new energy unit, which is expressed as follows:

[0072] (2)

[0073] Among them, is the transfer function of the complex frequency domain impedance of the new energy unit, is the order of each segment in the segmented fitting (its value is generally selected as 1), is the th pole, , , , are the parameters to be obtained through vector fitting respectively, is the independent variable.

[0074] Performing vector fitting based on the transfer function form of the complex frequency domain impedance of the new energy unit specifically includes the following steps:

[0075] S3.1. Uniformly select N initial poles from the fitting frequency band, and define an unknown function that satisfies that after multiplying it with the fitting transfer function, it has the same poles as , which is expressed as follows:

[0076] (3)

[0077] Among them, To fit the transfer function, is the conjugate value of; Multiply the second line formula on the right side of formula (3) by as shown below:

[0078] (4)

[0079] where, is the frequency domain sampling point, , a, b, are the unknown parameters to be solved respectively.

[0080] S3.2. Substitute the frequency points in step S101 to solve the vector fitting linear equation system formula, as shown below:

[0081] (5)

[0082]

[0083]

[0084]

[0085] where, is the coefficient matrix in the vector fitting algorithm, is the number of frequency sampling points, is the unknown vector of the vector fitting equation, is the right side term of the vector fitting algorithm, is the complex frequency variable, is the frequency response value of the objective function, is the th residue of the pole, is the th conjugate value of the residue of the pole, is the transpose operation on the matrix.

[0086] S3.3. Based on steps S3.1 and S3.2, after solving , a, b, parameters, obtain the expression of the fitted transfer function according to formula (3), as shown below:

[0087] (6)

[0088] where, and are respectively and zeros, and formula (6) can be expressed as:

[0089] (7)

[0090] Among them, the zero point of is the pole of the fitted transfer function and can be obtained by solving the eigenvalues, expressed as follows: which is shown as:

[0091] (8)

[0092] Among them, is an instruction in Matlab (used to obtain the eigenvalues and eigenvectors of matrix ), is a diagonal matrix whose diagonal elements are the initial poles ; is a column vector with all elements being 1; is the residue and is a row vector.

[0093] S3.4. Repeat the calculation process of formulas (5) - (8) in steps S3.2 - S3.3 to gradually obtain a more accurate vector fitting result.

[0094] S33. Based on the vector fitting result, integrate the transfer functions of the complex frequency domain impedances of the new energy units in several frequency bands to obtain the piecewise transfer function of the complex frequency domain impedance of the new energy units as which is shown as:

[0095]

[0096] Among them, , are the transfer functions of the first and the th segments respectively after vector fitting, is the transfer function of the complex frequency domain impedance of the new energy unit, is the frequency, and are the complex frequency domain impedance the th and the th peak frequencies corresponding to the Bode plot.

[0097] S34. In the case of no segmentation, repeat step S32 to obtain the overall transfer function of the complex frequency domain impedance of the new energy unit as . Weight the piecewise transfer function and the overall transfer function to obtain the final transfer function of the complex frequency domain impedance of each new energy unit as , which is expressed as follows:

[0098]

[0099] Among them, and are weight coefficients respectively, and generally take the value of 0.5.

[0100] S104. Based on the complex frequency domain impedance models of each new energy unit, combined with the line topology information and line parameters in the industrial park, establish the complex frequency domain admittance matrix of the industrial park, which specifically includes the following:

[0101] Establish multiple nodes according to the topological structure of the industrial park;

[0102] For the diagonal elements , the diagonal element of each node is the algebraic sum of the complex frequency domain admittances of all components connected to this node;

[0103] For the non - diagonal elements , the non - diagonal element represents the negative value of the complex frequency domain admittance between two nodes;

[0104] Symmetry. Each element of the complex frequency domain admittance matrix is symmetric about the main diagonal, that is = ;

[0105] Based on the above characteristics of the complex frequency domain admittance matrix of the industrial park, the four - node complex frequency domain admittance matrix is expressed as follows:

[0106]

[0107] If a certain node is connected in parallel with a certain grounding branch , the complex frequency domain admittance matrix can be transformed into:

[0108]

[0109] The topological diagram of the four - node industrial park is as shown in Figure 2 , including photovoltaic power generation equipment, wind power generation equipment, hydrogen production equipment, infinite power source and connecting wires; the wind power generation capacity is 500kW, the photovoltaic power generation capacity is 10kW, the hydrogen production equipment capacity is 10kW, and the specific parameters of the wires and generators are shown in Table 1 below:

[0110] Table 1

[0111]

[0112] The energy storage - side control structure of the said photovoltaic power generation equipment is as shown in Figure 3 . In the figure, is the DC - side voltage, , , constitute a filter, is a passive damping, is an infinite power source, is the inverter output current, , is the given output current, is the system angular frequency, usually taken as ; is the PWM modulation signal, , are the current sampling delay and voltage sampling delay respectively, is the frequency domain expression of the delay, is the sampling delay, expressed as follows:

[0113]

[0114] is the phase-locked loop transfer function, expressed as follows:

[0115]

[0116] wherein, and are the proportional control coefficient and integral control coefficient of the phase-locked loop respectively; the controller voltage loop first samples the grid voltage, and after a sampling delay, the result is sent to the phase-locked loop, and the phase angle information output by the phase-locked loop is provided to the current loop;

[0117] The current loop is the control link of the system output. The sampled current is sent to the control system after sampling delay, and the current information is converted to the coordinate system through Clark transformation; at the same time, the control signal is also converted to the coordinate system through Clark transformation, compared with the sampled signal difference to obtain an error, the error is corrected by the PR controller, and then through the inverse Clark transformation, and finally used as the PWM modulator signal to control the PWM output.

[0118] As Figure 4 shown, the vector fitting results of the photovoltaic power generation equipment are shown. The disturbance injection frequency points are selected according to 1:5:1000, and the fitting order is set to 8th order. The scatter points in the figure are the actual measurement data, and the solid line is the vector fitting result. It can be seen that the vector fitting result is consistent with the actual result.

[0119] Based on the complex frequency domain impedance models of each new energy unit, this embodiment combines the line topology information and line parameters in the four-node industrial park, and according to the four-node complex frequency domain admittance matrix, establishes the complex frequency domain admittance matrix of the four-node industrial park, which specifically includes the following steps:

[0120] S41. Collect and analyze the complex frequency domain impedance information of the photovoltaic power generation equipment, wind power generation equipment, and hydrogen production equipment at the corresponding disturbance frequency points respectively;

[0121] S42. Draw a Bode plot based on the complex frequency domain impedance information, and determine the order of fitting to be selected according to the number of peaks in the Bode plot;

[0122] S43. Perform vector fitting on the four-node industrial park according to the order of fitting to be selected to obtain the impedance result;

[0123] The impedance results of various equipment in the four-node industrial park can be expressed as:

[0124] Photovoltaic power generation equipment:

[0125] (1820063171227027 / 35184372088832 - 3247379264791407i / 281474976710656) / (s + (2919633481258255 / 140737488355328 - 2249213456603405i / 140737488355328)) +(1820063171227027 / 35184372088832 + 3247379264791407i / 281474976710656) / (s +(2919633481258255 / 140737488355328 + 2249213456603405i / 140737488355328)) +1968152767294939 / (70368744177664*(s - 6329008279353853 / 1125899906842624)) +4311938400802685 / (8796093022208*(s - 8524990996596231 / 35184372088832)) -(2642684564939107 / 274877906944 + 5340351738696959i / 1099511627776) / (s +(6620793037775075 / 8796093022208 - 8308519720416377i / 17592186044416)) -(2642684564939107 / 274877906944 - 5340351738696959i / 1099511627776) / (s +(6620793037775075 / 8796093022208 + 8308519720416377i / 17592186044416)) +(8362356439062255 / 549755813888 - 3159139800790623i / 549755813888) / (s +(2082175259461917 / 4398046511104 - 6421300226697625i / 1099511627776)) +(8362356439062255 / 549755813888 + 3159139800790623i / 549755813888) / (s +(2082175259461917 / 4398046511104 + 6421300226697625i / 1099511627776))+4721327007789329 / 18014398509481984

[0126] Hydrogen production equipment:

[0127] 2803409786950655 / (8796093022208*(s - 5919311367852555 / 70368744177664)) - (2662951035825651 / 274877906944 - 8274656406909023i / 1099511627776) / (s + (6390597773090793 / 8796093022208 + 3437403568279981i / 8796093022208)) - (2662951035825651 / 274877906944 + 8274656406909023i / 1099511627776) / (s + (6390597773090793 / 8796093022208 - 3437403568279981i / 8796093022208)) + 5764726097722303 / (140737488355328*(s + 4315188373007271 / 562949953421312)) + (8381884063123789 / 549755813888 - 6457246162570513i / 1099511627776) / (s + (1043763906397199 / 2199023255552 - 1606099629524209i / 274877906944)) + (8381884063123789 / 549755813888 + 6457246162570513i / 1099511627776) / (s + (1043763906397199 / 2199023255552 + 1606099629524209i / 274877906944)) + 3187544514323967 / 18014398509481984

[0128] Wind power generation equipment:

[0129] 4224399476078235 / (17592186044416*(s + 8204319024818089 / 140737488355328)) + (7117966373798719 / 35184372088832 - 5325168461125201i / 140737488355328) / (s + (2443071200439595 / 17592186044416 - 3023421061775005i / 8796093022208)) + (7117966373798719 / 35184372088832 + 5325168461125201i / 140737488355328) / (s + (2443071200439595 / 17592186044416 + 3023421061775005i / 8796093022208)) + 4503261624015589 / (281474976710656*(s - 5956685792391061 / 4503599627370496)) + 2613785018373549 / (4398046511104*(s - 6021643760586325 / 8796093022208)) + (8300149847756637 / 549755813888 - 6260927296545137i / 1099511627776) / (s + (2071349039847839 / 4398046511104 - 802539422664075i / 137438953472)) + (8300149847756637 / 549755813888 + 6260927296545137i / 1099511627776) / (s + (2071349039847839 / 4398046511104 + 802539422664075i / 137438953472)) + 5730719476778187 / (35184372088832*(s - 5406623362818137 / 70368744177664)) - (2763423522391657 / 274877906944 + 2536174117970495i / 549755813888) / (s + (433264647378109 / 549755813888 -(243797366785953i / 549755813888) - (2763423522391657 / 274877906944 - 2536174117970495i / 549755813888) / (s + (433264647378109 / 549755813888 + 243797366785953i / 549755813888)) + 7786431676834377 / 18014398509481984

[0130] S44. Based on step S43, first establish the complex frequency domain admittance matrix without new energy units to obtain the four-node complex frequency domain admittance matrix. Each row of the four-node complex frequency domain admittance matrix is expressed as follows:

[0131] The first row:

[0132] [1 / ((17*s) / (25000*pi) + 1 / 100) + 1 / ((23*s) / (25000*pi) + 17 / 1000) + 1 / ((161*s) / (100000*pi) + 4 / 125) + 1 / ((3617*s) / (1000000*pi) + 319 / 5000) + s / (625*pi),-1 / ((17*s) / (25000*pi)+1 / 100),-1 / ((23*s) / (25000*pi)+17 / 1000),-1 / ((161*s) / (100000*pi) + 4 / 125)]

[0133] The second row:

[0134] [-1 / ((17*s) / (25000*pi) + 1 / 100), 1 / ((17*s) / (25000*pi) + 1 / 100) + (11*s) / (25000*pi),0,0]

[0135] The third row:

[0136] [-1 / ((23*s) / (25000*pi)+17 / 1000),0,1 / ((23*s) / (25000*pi)+17 / 1000)+ (79*s) / (200000*pi),0]

[0137] The fourth row:

[0138] [-1 / ((161*s) / (100000*pi)+4 / 125),0,0, 1 / ((161*s) / (100000*pi) + 4 / 125)+ (153*s) / (200000*pi)]

[0139] After establishing the complex frequency domain admittance matrix after adding new energy units, the second diagonal element corresponds to the photovoltaic power generation equipment, the third corresponds to the hydrogen production equipment, and the fourth corresponds to the wind power generation equipment, which are specifically expressed as follows:

[0140] The first row:

[0141] [1 / ((17*s) / (25000*pi) + 1 / 100) + 1 / ((23*s) / (25000*pi) + 17 / 1000) + 1 / ((161*s) / (100000*pi) + 4 / 125) + 1 / ((3617*s) / (1000000*pi) + 319 / 5000) + s / (625*pi),-1 / ((17*s) / (25000*pi)+1 / 100),-1 / ((23*s) / (25000*pi)+17 / 1000),-1 / ((161*s) / (100000*pi) + 4 / 125)]

[0142] The second row:

[0143] [-1 / ((17*s) / (25000*pi)+1 / 100),1 / ((17*s) / (25000*pi)+1 / 100)+(11*s) / (25000*pi)+(1820063171227027 / 35184372088832 - 3247379264791407i / 281474976710656) / (s+(2919633481258255 / 140737488355328 - 2249213456603405i / 140737488355328))+(1820063171227027 / 35184372088832 + 3247379264791407i / 281474976710656) / (s+(2919633481258255 / 140737488355328 + 2249213456603405i / 140737488355328))+1968152767294939 / (70368744177664*(s - 6329008279353853 / 1125899906842624))+4311938400802685 / (8796093022208*(s - 8524990996596231 / 35184372088832))-(2642684564939107 / 274877906944 + 5340351738696959i / 1099511627776) / (s+(6620793037775075 / 8796093022208 - 8308519720416377i / 17592186044416))-(2642684564939107 / 274877906944 - 5340351738696959i / 1099511627776) / (s+(6620793037775075 / 8796093022208 + 8308519720416377i / 17592186044416))+(8362356439062255 / 549755813888 - 3159139800790623i / 549755813888) / (s+(2082175259461917 / 4398046511104 - 6421300226697625i / 1099511627776))+(8362356439062255 / 549755813888 +3159139800790623i / 549755813888) / (s + (2082175259461917 / 4398046511104 + 6421300226697625i / 1099511627776)) + 4721327007789329 / 18014398509481984,0,0]

[0144] The third line:

[0145] [-1 / ((23*s) / (25000*pi)+17 / 1000),0, 1 / ((23*s) / (25000*pi)+17 / 1000)+(79*s) / (200000*pi)+2803409786950655 / (8796093022208*(s - 5919311367852555 / 70368744177664)) - (2662951035825651 / 274877906944 - 8274656406909023i / 1099511627776) / (s + (6390597773090793 / 8796093022208 + 3437403568279981i / 8796093022208)) - (2662951035825651 / 274877906944 + 8274656406909023i / 1099511627776) / (s + (6390597773090793 / 8796093022208 - 3437403568279981i / 8796093022208)) + 5764726097722303 / (140737488355328*(s + 4315188373007271 / 562949953421312)) + (8381884063123789 / 549755813888 - 6457246162570513i / 1099511627776) / (s + (1043763906397199 / 2199023255552 - 1606099629524209i / 274877906944)) + (8381884063123789 / 549755813888 + 6457246162570513i / 1099511627776) / (s + (1043763906397199 / 2199023255552 + 1606099629524209i / 274877906944)) + 3187544514323967 / 18014398509481984,0]

[0146] Fourth line:

[0147] [-1 / ((161*s) / (100000*pi)+4 / 125),0,0, 1 / ((161*s) / (100000*pi) + 4 / 125)+ (153*s) / (200000*pi)+ 4224399476078235 / (17592186044416*(s +8204319024818089 / 140737488355328)) + (7117966373798719 / 35184372088832 -5325168461125201i / 140737488355328) / (s + (2443071200439595 / 17592186044416 -3023421061775005i / 8796093022208)) + (7117966373798719 / 35184372088832 +5325168461125201i / 140737488355328) / (s + (2443071200439595 / 17592186044416 +3023421061775005i / 8796093022208)) + 4503261624015589 / (281474976710656*(s -5956685792391061 / 4503599627370496)) + 2613785018373549 / (4398046511104*(s -6021643760586325 / 8796093022208)) + (8300149847756637 / 549755813888 -6260927296545137i / 1099511627776) / (s + (2071349039847839 / 4398046511104 -802539422664075i / 137438953472)) + (8300149847756637 / 549755813888 +6260927296545137i / 1099511627776) / (s + (2071349039847839 / 4398046511104 +802539422664075i / 137438953472)) + 5730719476778187 / (35184372088832*(s -5406623362818137 / 70368744177664)) - (2763423522391657 / 274877906944((+2536174117970495i / 549755813888) / (s + (433264647378109 / 549755813888 - 243797366785953i / 549755813888))) - ((2763423522391657 / 274877906944 - 2536174117970495i / 549755813888) / (s + (433264647378109 / 549755813888 + 243797366785953i / 549755813888))) + 7786431676834377 / 18014398509481984

[0148] The above are only the preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principles of the present invention shall be included within the protection scope of the present invention.

Claims

1. The complex frequency domain admittance matrix modeling method of industrial parks based on vector fitting technology is characterized by: The method comprises the following steps: S101, using the grid-connected point disturbance method to measure the real-time amplitude-frequency and phase-frequency characteristic curves of each new energy unit in the grid-connected operation state of the industrial park, and obtaining the impedance amplitude and phase frequency points of each new energy unit, specifically including the following steps: S11. Collect the three-phase voltage and current operation data of each new energy unit, use a feedforward neural network to fit the three-phase voltage and current operation data of each new energy unit, and establish a digital twin model; S12, apply disturbance frequencies f at different frequency points to each new energy unit at the grid connection point p The small signal voltage disturbance is , where f is the frequency and p is the pth frequency point; S13. Calculate the difference between the digital twin model and the grid-connected current of each renewable energy unit after the small signal voltage disturbance to obtain the small signal current response, and perform Fourier analysis on it to obtain the current response of each renewable energy unit under the voltage disturbance in the complex frequency domain, which is expressed as follows: Wherein, X(k) is the frequency domain signal, DFT is discrete Fourier analysis, x(m) is the function for Fourier analysis, i.e., the small signal current response, m is the sampling point number, M is the number of sampling points, e is the natural logarithm, j is the imaginary unit, and k is a natural number; S14. Calculate the disturbance frequency f based on the Fourier analysis results p The ratio of the small signal current response and the small signal voltage disturbance intensity at , and the corresponding complex frequency domain impedance Z p ; S15. Modify the disturbance frequency f p , and obtain the amplitude-frequency and phase-frequency characteristic curves of the impedance Z of the new energy unit; S16. When each new energy unit is running, a small signal voltage disturbance is injected in real time to obtain the complex frequency domain impedance Z p ′, and according to the complex frequency domain impedance Z p ' data to determine whether to update the impedance Z of the new energy unit; S102, based on the frequency points of the impedance amplitude and phase of each new energy unit, obtaining the peak number of each new energy unit, determining the order of the complex frequency domain impedance of the new energy unit according to the peak number of each new energy unit, and then determining the order of vector fitting; S103, based on the order of vector fitting, using vector fitting technology to segmentally identify the complex frequency domain impedance model of each new energy unit, and obtain segmented accurate parameters of the complex frequency domain impedance model, specifically including the following steps: S31, segmenting the transfer function of the complex frequency domain impedance of each new energy unit; S32, defining the transfer function form of the complex frequency domain impedance of the new energy generator set in the i-th frequency band, and performing vector fitting based on the transfer function form of the complex frequency domain impedance of the new energy generator set, which is expressed as follows: Among them, H i (s) is the transfer function of the complex frequency domain impedance of the new energy unit, N is the order of each segment in the segmented fitting, n is the nth pole, a, b, c, d n They are the parameters to be determined through vector fitting, and s is the independent variable; The vector fitting is performed based on the transfer function form of the complex frequency domain impedance of the new energy unit, and specifically includes the following steps: S3.

1. Evenly select N initial poles from the fitting frequency band Define an unknown function H i (s)′, which satisfies the multiplication of the fitting transfer function with H i (s) have the same poles, which are expressed as follows: in, To fit the transfer function, is the conjugate value of c; S3.2, bring in the frequency points in step S101, and solve the vector fitting linear equations formula, which is expressed as follows: B k =f(s k ) Among them, A k is the coefficient matrix in the vector fitting algorithm, k is the number of frequency sampling points, x is the unknown vector of the vector fitting equation, B is the right-hand side term of the vector fitting algorithm, and s k is a complex frequency variable, f(s k ) is the frequency response value of the objective function, c N is the residual of the Nth pole, is the conjugate value of the residual of the Nth pole, and T is the transposition operation of the matrix; S3.

3. Based on steps S3.1 and S3.2, the fitting transfer function is obtained, which is expressed as follows: in, To fit the transfer function, {z n }and They are and H i The zero point of (S)′, H i Zero point of (s)′ To fit the poles of the transfer function, we can obtain S3.4, repeat the calculation process of steps S3.2-S3.3 to gradually obtain more accurate vector fitting results; S33. Based on the vector fitting results, the transfer functions of the complex frequency domain impedance of the new energy generator set in N frequency bands are integrated to obtain the piecewise transfer function of the complex frequency domain impedance of the new energy generator set as H par (s), which is represented as follows: Among them, H1(s), H N (s) are the transfer functions of the first and Nth segments after vector fitting, respectively. i (s) is the transfer function of the complex frequency domain impedance of the new energy unit, f is the frequency, f i and f N is the complex frequency domain impedance Z p The frequencies corresponding to the i-th and N-th peaks on the Bode plot; S34, without segmentation, repeat step S32 to obtain the entire segment transfer function of the complex frequency domain impedance of the new energy unit as H all (s), weight the piecewise transfer function and the whole-segment transfer function, and obtain the final transfer function of the complex frequency domain impedance of each new energy unit as H(s), which is expressed as follows: H(s)=αH all (s)+βH par (s) Among them, α and β are weight coefficients respectively; S104, based on the complex frequency domain impedance model of each new energy unit, combined with the line topology information and line parameters in the industrial park, establish a complex frequency domain admittance matrix of the industrial park, specifically including the following steps: S41, respectively collecting and analyzing the complex frequency domain impedance information of the photovoltaic power generation equipment, the wind power generation equipment, and the hydrogen production equipment at the corresponding disturbance frequency points; S42, drawing a Bode plot according to the complex frequency domain impedance information, and determining the order of fitting to be selected according to the number of peaks of the Bode plot; S43, performing vector fitting on the industrial park according to the order of fitting to be selected to obtain impedance results; S44. Based on step S43, first establish the complex frequency domain admittance matrix without the new energy unit and the representation of each row of the complex frequency domain admittance matrix, and then establish the complex frequency domain admittance matrix after adding the new energy unit. The second row of diagonal elements is photovoltaic power generation equipment, the third row is hydrogen production equipment, and the fourth row is wind power generation equipment.

2. The industrial park complex frequency domain admittance matrix modeling method based on vector fitting technology according to claim 1 is characterized in that: The step S102 includes: drawing a Bode diagram for each frequency point obtained in step S101, and determining the order N to be selected for fitting according to the number of peaks in the Bode diagram.

3. The industrial park complex frequency domain admittance matrix modeling method based on vector fitting technology according to claim 1 is characterized in that: The complex frequency domain admittance matrix of the industrial park includes: According to the topological structure of the industrial park, multiple nodes are established, and the diagonal element Y of each node ii is the algebraic sum of the complex frequency domain admittances of all components connected to the node; The off-diagonal elements represent the negative values ​​of the complex frequency domain admittance between two nodes; The elements of the complex frequency domain admittance matrix are symmetric about the main diagonal, that is, Y ji =Y ij .

Images