A physically feasible broadband impedance equivalent modeling method for strong electromagnetic equipment

Abstract

The present invention discloses a physically implementable strong electromagnetic equipment wide-band impedance equivalent modeling method. First, according to the vector matching method, the internal impedance of the strong electromagnetic equipment is expressed in the form of a rational function; the obtained rational function is subjected to a network port passivation process by using the residue perturbation method to obtain a circuit network function; the obtained circuit network function is converted into a corresponding equivalent circuit by using the impedance equivalent relationship between the rational function and the parameters R, L, and C to form the circuit; for non-positive parameters in the equivalent circuit, the non-positive parameters are corrected to positive values ​​by using a particle swarm algorithm to obtain a final equivalent circuit so as to perform physical equivalence with actual R, L, and C components. The method can establish a two-port equivalent circuit composed entirely of passive components through the measured impedance data of loads such as motors, can be directly used for equivalent physical simulation of the impedance of strong electromagnetic equipment, and can lay a foundation for experimental verification.

Description

Technical Field

[0001] The invention relates to an improved impedance equivalent simulation technology, and more specifically to a physically implementable strong electromagnetic equipment wide-band impedance equivalent modeling method, belonging to the technical field of electromagnetic compatibility prediction. Background Art

[0002] With the rapid development and progress of new materials, new power electronic devices, and new control technologies, ships have gradually adopted integrated power systems to replace traditional power systems, providing a basis for the use of high-power and high-energy equipment, and developing towards the second-generation integrated power system. In addition, the power density is further improved, the voltage level is further improved, the volume and weight are further reduced, and the degree of intelligence is further improved. Control systems are more susceptible to electromagnetic interference generated by themselves and surrounding strong electromagnetic equipment related to power generation, transmission and transformation, which brings more serious electromagnetic interference risks. The electromagnetic compatibility of equipment, control systems and operating software will also directly affect the realization of control objectives and the safe and stable operation of the system. Therefore, in the electromagnetic compatibility design stage of the platform, it is necessary to simulate and evaluate the electromagnetic compatibility performance of the power system in the platform by performing electromagnetic compatibility equivalent modeling of various strong electromagnetic equipment, and propose targeted control measures for weak links. At the same time, it is also necessary to perform physical equivalent simulation of strong electromagnetic equipment as the load of the power grid, so as to realize the semi-physical electromagnetic compatibility physical simulation of the power system and each strong electromagnetic equipment in the laboratory environment, test the electromagnetic interference state at key locations on the power grid, and verify the electromagnetic interference suppression effect of control measures (such as high-power EMI filters). However, due to the complexity of the technical system, composition, and connection method of the strong electromagnetic equipment on the ship platform, it is difficult to obtain the internal specific information of the power generation, transmission, and transformation system, especially complex equipment such as diesel generator sets and permanent magnet propulsion motors. In order to solve this problem, some scholars have adopted the vector fitting method to achieve modeling and simulation of the electromagnetic compatibility characteristics of the motor through experimental test data.

[0003] Vector Fitting is a stable and effective rational function fitting method proposed by Norwegian scholar B. Gustavsen in 1999. It avoids the ill-conditioned and unbalanced weighting problems that occur in the rational approximation process, and has become a popular macromodel extraction method and a tool for frequency domain linear system identification. It is essentially an improvement on the Sanathan-Koerner iterative method based on rational partial fractions. This method can calculate broadband network functions from frequency domain responses, time domain responses, or a combination of the two, and achieves its robustness through iterative least squares estimation and rational partial fraction basis functions. So far, the vector matching method has been widely used in the fields of electromagnetic transient simulation of power systems, complete signal simulation in microwave systems, and has achieved good results.

[0004] However, when applying the vector matching method to model strong electromagnetic equipment, it can only correspond to each frequency point, and the problem of port passivity of the equivalent circuit network model cannot be solved, which can easily lead to model divergence. Moreover, the equivalent circuit network model established often has the phenomenon of "non-positive value" of components, which cannot be used to build physically realizable strong electromagnetic equipment simulation loads, and is difficult to use for electromagnetic compatibility semi-physical tests of platform power grids. Building a physically realizable circuit network model is of great significance for the physical simulation of the integrated power system on ships, and can be used to build a low-power conducted electromagnetic interference model for verification. Summary of the invention

[0005] In view of the above-mentioned deficiencies in the prior art, the purpose of the present invention is to provide a physically implementable method for wide-band impedance equivalent modeling of strong electromagnetic equipment. This method can establish a two-port equivalent circuit composed entirely of passive components through the measured impedance data of loads such as motors, which can be directly used for equivalent physical simulation of the impedance of strong electromagnetic equipment and lay the foundation for experimental verification.

[0006] The technical solution of the present invention is achieved in this way:

[0007] A physically implementable method for equivalent modeling of broadband impedance of strong electromagnetic equipment is carried out in the following steps:

[0008] 1) According to the vector matching method, the internal impedance of the strong electromagnetic device to be impedance equivalent is expressed in the form of a rational function;

[0009] 2) For the rational function obtained in step 1), the network port is subjected to passivation processing by using the residue perturbation method to obtain the circuit network function;

[0010] 3) The circuit network function obtained by passivation is converted into the corresponding equivalent circuit using the impedance equivalent relationship between the rational function and the parameters R, L, and C.

[0011] 4) For the non-positive parameters of R, L, and C in the equivalent circuit, the particle swarm algorithm is used to correct the non-positive parameters in the equivalent circuit to positive values, that is, to obtain the final equivalent circuit so that the actual R, L, and C components can be used for physical equivalence.

[0012] Among them, in step 1), the rational function representing the internal impedance of the strong electromagnetic device is:

[0013]

[0014] In the formula, the residue r n and the extreme point p nis a real number or a conjugate complex number pair; d and e are real numbers; Z(s) is the internal impedance of a strong electromagnetic device that changes with frequency; s indicates that the function is a complex frequency domain function, n = 1, 2, ..., N, N is the order of the polynomial (1);

[0015] Based on the wide frequency band test data Z(s k ), s k (k=1,2,…,p) is the test frequency, p is the number of test frequency; determine the residue r in the rational function by iteration n , extreme point p n , real number d, real number e are as follows:

[0016] 1.1) Remember is the initial pole, and the auxiliary function σ(s) is introduced, as shown in formula (2);

[0017] Multiply the auxiliary function σ(s) by the rational function Z(s) of the internal impedance of the strong electromagnetic device, and assume that during the iteration, σZ(s) and σ(s) have the same poles.

[0018]

[0019]

[0020] In the formula, is the residue of the auxiliary function σ(s);

[0021] 1.2) Using equations (2) and (3), establish the parameter r n , The linear equation of d and e is:

[0022]

[0023] 1.3) The test data Z(s) k ) and the corresponding frequency s k Substituting into equation (4), we get the linear equation, and the matrix form of the linear equation is as follows:

[0024] Ax=b (5)

[0025] Among them, the coefficient matrix A P×(2N+2) It is expressed as follows:

[0026]

[0027] The vector x to be solved and the known vector b are as follows:

[0028]

[0029] 1.4) Solve the linear equation matrix (5) using the least squares method to obtain r n , d and e;

[0030] 1.5) Rewrite the above function σZ(s) and auxiliary function σ(s) into zero-pole form, that is:

[0031]

[0032] In the formula, z n and are the zeros of σZ(s) and σ(s) respectively;

[0033] 1.6) According to formulas (9)-(10), the rational function Z(s) of the internal impedance of the strong electromagnetic device is updated to obtain:

[0034]

[0035] From (11), we can see that the pole of the function Z(s) is the zero of the auxiliary function σ(s), and the zero of the auxiliary function σ(s) is the initial pole of the next iteration Z(s);

[0036] 1.7) According to equation (2) and equation (10), equation (12) is established, that is:

[0037]

[0038] Writing the right side of equation (12) in the form of a determinant yields:

[0039]

[0040] 1.8) Change the determinant H on the right side of equation (13) into the following form:

[0041]

[0042] In the formula, P is the An N-order diagonal matrix, U is a 1×N-dimensional column vector with 1 elements, is the residue containing σ(s) 1×N order row vector; In formula (14), the diagonal matrix P, the row vector They are as follows:

[0043]

[0044] 1.9) Since the characteristic equation of matrix H is |sE-H|=0, E is the unit determinant, the matrix H is replaced by equation (14) combined with equations (15) and (16), and then becomes |PR|=0; further using the matrix R formed by the residues of equation (1), as shown in equation (17); according to equation (13) and eig|PR|=0, the zero point of the auxiliary function σ(s) is obtained That is, get the initial value of the next iteration

[0045]

[0046] 1.10) Repeat steps 1.4) to 1.9) to obtain r n , d, e, and Determine whether the data error of the function Z(s) obtained in two consecutive iterations is less than the threshold. If so, end the iteration. Otherwise, define the latest pole as the initial pole. Return to step 1.4) and recalculate until the required r is obtained. n 、p n , d and e.

[0047] Among them, the steps of step 2) passive processing are as follows:

[0048] 2.1) Establish the following constraint equation for network passivity:

[0049]

[0050] That is, all eigenvalues ​​of the real part of the impedance matrix of the model are positive for all frequencies; the following constraint equation is established:

[0051]

[0052] ΔX represents the correction deviation added to X as the quantity to be determined, X is the unified representation of the required correction quantity in equations (19)-(22), and ΔX is an unknown quantity that needs to be solved by establishing a quadratic programming equation;

[0053] 2.2) According to constraints (19)-(22), establish the quadratic programming equation for passive optimization:

[0054]

[0055] -RΔx≤λ(24)

[0056] In the formula, Δx is the increment of parameter x, and x is r n 、p n , d, e; λ=eig(Re(Z)), R=MRe(Z),

[0057] 2.3) Convert the rational function formula (1) into the form of matrix multiplication as shown in formula (25);

[0058] Z(s)=C(sI-A) -1 B+D+sE(25)

[0059] Among them, Z(s) represents the form of converting Z(s) into a matrix, I is the N-order unit matrix, C=[r1 r2 … r N ], D = [d], E = [e], the matrix A is expressed as follows:

[0060]

[0061] 2.4) Perform positive definite processing on matrix D and matrix E, that is:

[0062]

[0063] Wherein, the elements in the matrix D and the matrix E include real numbers d and e respectively; T D is the normalized orthogonal matrix of D; D is the eigenvalue matrix of D; T E is the normalized orthogonal matrix of E; E is the eigenvalue matrix of E;

[0064] Λ D and Λ E Replace the negative eigenvalues ​​in with zero, update the matrix D and the matrix E, and get the matrix D mod and the matrix E mod ;

[0065] 2.5) In equation (25), the complex matrices A and C are simplified to the Jordan standard form matrix A mod and the matrix C mod ;

[0066] 2.6) Determine whether the eigenvalue of the real part of Z(s) at each frequency point without correction is greater than 0. If so, it is a passive system; otherwise, collect the frequency point information s that does not meet the conditions w , and perform the following steps;

[0067] 2.7) The frequency information s collected in step 2.6) is w , arrange them in ascending order, and take the median between every two frequency points to establish a new frequency list s mod , and check the passivity at these frequencies. If the real part of Z(s) is modThe first frequency eigenvalue of is less than 0, then the starting frequency of passive optimization is the starting frequency of the vector matching method. If the real part of Z(s) is in s mod If the last frequency eigenvalue of is less than 0, the termination frequency of passive optimization is the termination frequency of vector matching method. Otherwise, the passive optimization interval should at least satisfy s w The interval satisfied;

[0068] 2.8) Use the quadratic programming algorithm to calculate the least squares solution that satisfies the constraints (23) and (24) until the correction parameters that satisfy the constraints are found, including the residue r n,mod , first-order term e mod , constant term d mod Otherwise, continue to repeat steps 2.4)-2.8); the rational function of the internal impedance after passivation is

[0069]

[0070] Among them, r n,mod =r n +Δr n , d mod =d+Δd,e mod =e+Δe, these correction values ​​are solved by quadratic programming of equations (19)-(22).

[0071] Step 3) converting into a corresponding equivalent circuit as follows;

[0072] 3.1) For the constant term d and the linear term e of a rational function, the equivalent circuit in impedance form is a series circuit of R and L, and the values ​​of the parameters R and L are R = d and L = e respectively;

[0073] 3.2) When the poles and residues of a rational function are real numbers, an RC parallel circuit is used to represent the equivalent, where:

[0074]

[0075] 3.3) When the poles and residues of a rational function are conjugate complex numbers, the corresponding equivalent circuit is a second-order circuit, which includes two resistors, one inductor and one capacitor, wherein one resistor and one inductor are connected in series and then connected in parallel with the capacitor and another resistor. The values ​​of each circuit element are calculated according to formula (30);

[0076]

[0077] In the formula

[0078]

[0079] 3.4) Decompose equation (1) into three fractions (3.1), 3.2) and 3.3), each fraction is equivalent to a corresponding subcircuit, and finally all subcircuits are connected in series to obtain the equivalent circuit corresponding to the network function.

[0080] Step 4) The method of using the particle swarm algorithm to correct the non-positive parameters in the equivalent circuit to positive values ​​is:

[0081] 4.1) Set the objective function to

[0082]

[0083] Among them, ω i is the passive frequency of the component, M is the number of frequencies involved in the optimization, and Z mod is the rational function obtained by passivation, that is, the rational function of internal impedance (28), Z fit is the optimized resistance impedance function;

[0084] 4.2) Set the inertia weight to

[0085]

[0086] Among them, ω max and ω min are the maximum and minimum values ​​of the inertia weight respectively;

[0087] 4.3) Set the learning factor of the asynchronous change as shown in formula (34),

[0088]

[0089] Among them, c 1s and c 2s are the initial values ​​of the individual and population learning factors c1 and c2, respectively, c 1e and c 2e are the termination values ​​of the learning factors c1 and c2 respectively;

[0090] 4.4) PSO is initialized as a group of random particles, and then finds the optimal solution through iteration. In each iteration, the particles update themselves by tracking two extreme values; the first extreme value is the optimal solution found by the particle itself, which is called the individual extreme value pbest, and the other extreme value is the optimal solution found by the entire population, which is the global extreme value gbest; when these two extreme values ​​are found, the particles update their speed and new position according to the following formula:

[0091]

[0092] x j =x j +v j(36)

[0093] Where j = 1, 2, ..., N max , N max is the total number of particles in this group, v j is the speed of each particle, rand() is a random number between (0,1), x j is the position of each particle, c1 and c2 are the learning factors of individuals and populations, ω PSO is the inertia factor;

[0094] 4.5) According to formulas (35) and (36), the particle speed and position are adjusted until the optimal position searched by the particle swarm is met and the predetermined minimum adaptation threshold is met or the end condition of the maximum number of iterations is reached; if the end condition is not met, the iteration is continued.

[0095] The present invention firstly performs vector matching through the measured impedance characteristics, fits its rational function expression in the complex frequency domain, then performs passive optimization of the equivalent circuit network, corrects the poles, residues and other information in the rational function expression, ensures that the circuit is a strict "passive" network, and then converts the network function into a corresponding equivalent circuit. Finally, the particle swarm algorithm is used to find suitable values ​​for the "non-positive" elements in the circuit for replacement, and finally forms a strong electromagnetic device impedance equivalent circuit network modeling and physical simulation method based on the vector matching method and the particle swarm algorithm, which can be directly used for the equivalent physical simulation of the impedance of the strong electromagnetic device, can be used for the semi-physical verification test of the electromagnetic compatibility of the power system in the scenes such as ships, and for the modeling model prediction of the electromagnetic equipment whose specific circuit parameters are difficult to obtain, so as to provide a prediction model and data support for the prediction of the conducted interference of the strong electromagnetic equipment and the power system and the targeted electromagnetic interference suppression.

[0096] Compared with the prior art, the present invention has the following beneficial effects:

[0097] 1. The present invention uses a vector matching method to express the measured impedance data using a rational polynomial composed of poles and residues. Then, the circuit elements can be mapped according to the numerical states of the poles and residues (divided into real numbers, complex numbers, constant terms, and linear terms). Finally, the measured impedance can be equivalent to a two-port circuit network.

[0098] 2. The present invention utilizes the residue perturbation method and the particle swarm algorithm to realize "port passivity" and "component passivity" step by step. The equivalent circuit network finally established is composed entirely of positive R, L, and C components, and can be physically equivalent with actual R, L, and C components.

[0099] 3. The wide-band impedance equivalent modeling method of strong electromagnetic equipment proposed in the present invention can directly use the impedance test data of strong electromagnetic equipment such as motors to build its equivalent circuit model without understanding the internal topology and circuit specific information of strong electromagnetic equipment such as motors.

[0100] 4. The modeling ideas and methods proposed in the present invention can also be applied to wide-band circuit models of other types of strong electromagnetic equipment such as AC motors and converters, and have good versatility. BRIEF DESCRIPTION OF THE DRAWINGS

[0101] Figure 1 It is the equivalent circuit diagram of the constant term and the first-order term of the present invention.

[0102] Figure 2 This is the equivalent circuit diagram of the present invention when the poles and residues are real numbers.

[0103] Figure 3 This is an equivalent circuit diagram when the poles and residues of the present invention are conjugate complex number pairs.

[0104] Figure 4 It is the equivalent circuit diagram of the impedance parameters of the present invention.

[0105] Figure 5 Schematic diagram of the amplitude-frequency characteristic fitting result of an embodiment of the present invention.

[0106] Figure 6 Schematic diagram of the phase-frequency characteristic fitting result of an embodiment of the present invention. DETAILED DESCRIPTION

[0107] The present invention is described in detail below with reference to the accompanying drawings and specific embodiments.

[0108] The present invention firstly uses the impedance characteristic number of the strong electromagnetic equipment obtained by measurement using a measuring instrument, and then uses the vector matching method to perform vector matching, fits its rational function expression in the complex frequency domain, and then performs passive optimization of the equivalent circuit network, and performs "passive" correction on the poles, residues and other information in the rational function expression to ensure that the circuit is a strict "passive" network, and then converts the network function into a corresponding equivalent circuit network, and finally uses the particle swarm algorithm to optimize the "non-positive" elements in the circuit into suitable "positive" elements for replacement, and finally forms a passive circuit fitting method using the vector matching method and the particle swarm algorithm, which can be directly used for equivalent physical simulation of the impedance of strong electromagnetic equipment, can be used for semi-physical verification tests of electromagnetic compatibility of power systems in scenes such as ships, and can be used for modeling and model prediction of other electromagnetic equipment whose specific circuit parameters are difficult to obtain, so as to provide prediction models and data support for conducting interference prediction of strong electromagnetic equipment and power systems and targeted electromagnetic interference suppression.

[0109] Based on this method, it can be used to establish a physically realizable circuit model. The detailed steps of its implementation method are as follows:

[0110] (1) According to the vector matching method, the internal impedance of a strong electromagnetic device (such as a permanent magnet propulsion motor) is expressed in the form of a rational function.

[0111] According to the vector matching method, the internal impedance model of a strong electromagnetic device can be expressed by the following rational function:

[0112]

[0113] In the formula, the residue r n and the extreme point p n is a real number or a conjugate complex number pair; d and e are real numbers; Z(s) is the internal impedance of a strong electromagnetic device that changes with frequency; s indicates that the function is a complex frequency domain function, n = 1, 2, ..., N, N is the order of the polynomial (1).

[0114] Then, based on the wide-band test data Z(s k ), s k (k=1,2,…,p) is the test frequency, and p is the number of test frequencies. Determine the residue r in the internal impedance model of the strong electromagnetic device by iteration. n , extreme point p n , real number d, real number e include:

[0115] 1.1) Remember is the initial pole, which can be randomly selected and is a known quantity; at the same time, in order to solve the coefficients of Z(s), an auxiliary function σ(s) is introduced.

[0116] Multiply the auxiliary function σ(s) by the internal impedance function Z(s) of the strong electromagnetic device, and assume that during the iteration, σZ(s) and σ(s) have the same poles.

[0117]

[0118] In the formula, is the residue of the auxiliary function σ(s).

[0119] 1.2) Using equations (2) and (3), establish the parameter r n , The linear equation of d and e is:

[0120]

[0121] 1.3) The test data Z(s) k ) and the corresponding frequency s k Substituting into equation (4), we get the linear equation, and the matrix form of the linear equation is as follows:

[0122] Ax=b (5)

[0123] Among them, the coefficient matrix A P×(2N+2) It is expressed as follows:

[0124]

[0125] The vector x to be solved and the known vector b are as follows:

[0126]

[0127] b=[Z(s1) Z(s2) … Z(s k ) … Z(s p )] T (8)

[0128] 1.4) Solve the linear equation matrix (5) using the least squares method to obtain r n , d and e. It should be noted that the above parameters are at the initial extremes Assuming that n , d and e are not accurate, so it is necessary to and r n , d and e are iteratively corrected.

[0129] 1.5) Further rewrite the above function σZ(s) and auxiliary function σ(s) into zero-pole form, that is:

[0130]

[0131] In the formula, z n and are the zeros of σZ(s) and σ(s) respectively;

[0132] 1.6) According to formulas (9)-(10), the output function Z(s) of the internal impedance model of the strong electromagnetic device is updated to obtain:

[0133]

[0134] From (11), we can see that the pole of the function Z(s) is the zero of the auxiliary function σ(s), and the zero of the auxiliary function σ(s) is the initial pole of the next iteration Z(s).

[0135] 1.7) According to equation (2) and equation (10), equation (12) is established, that is:

[0136]

[0137] Writing the right side of equation (12) in the form of a determinant yields:

[0138]

[0139] 1.8) Change the determinant H on the right side of (13) into the following form:

[0140]

[0141] In the formula, P is the An N-order diagonal matrix, U is a 1×N-dimensional column vector with 1 elements, is the residue containing σ(s) 1×N order row vector; in (14), the diagonal matrix P, the row vector They are as follows:

[0142]

[0143] 1.9) Since the characteristic equation of matrix H is |sE-H|=0, and E is the unit determinant, replace the matrix H with (14) combined with (15) and (16), and it becomes |PR|=0. Further, the matrix R formed by the residues of formula (1) is shown in formula (17). According to formula (13) and eig|PR|=0, the zeros of the auxiliary function σ(s) can be obtained. The initial value of the next iteration can be obtained

[0144]

[0145] 1.10) Repeat steps 1.4) to 1.9) to obtain r n , d, e, and Determine whether the data error of the function Z(s) obtained in two consecutive iterations is less than the threshold (which can be set by yourself, the typical value is 10^-6). If so, end the iteration; otherwise, define the latest pole as the initial pole. Return to step 1.4) and recalculate until the required r is obtained. n 、p n , d and e. This allows Z(s) to be written in the form of a rational function (1).

[0146] (2) For the rational function obtained in the previous step, the residue perturbation method is used to perform “passivation” processing on the network port.

[0147] Since the vector matching method is a mathematical fit of the measurement results of the port broadband characteristics, in some cases the model will not meet the passivity requirements, causing the simulation of the entire system model to diverge. If a circuit network is made a passive network, its impedance matrix Z(s) must meet the following two conditions:

[0148] ①Z(s) should be an analytical function;

[0149] ②Re(Z(s))>0.

[0150] Based on the above conditions, the steps for passive optimization of the internal impedance model Z(s) of strong electromagnetic equipment are as follows:

[0151] 2.1) Establish the constraint equations for network passivity. A passive model should satisfy:

[0152]

[0153] That is, all eigenvalues ​​of the real part of the impedance matrix of the model are positive for all frequencies. Most passive optimization methods rely on a post-processing step of the model, assuming that only small corrections are required, thereby establishing the following constraint equations:

[0154]

[0155] eig{Re(Z+ΔZ)}>0 (20)

[0156] eig(d+Δd)>0 (21)

[0157] eig(e+Δe)>0 (22)

[0158] ΔX represents the correction deviation added to X as the quantity to be determined, X is the unified representation of the required correction quantity in (19)-(22), ΔX is an unknown quantity that needs to be solved by establishing a quadratic programming equation. Constraint (19) is to ensure that the impedance correction is very small, constraint (20) is to ensure that the circuit can meet the passive requirements after correction, and constraints (21) and (22) are to ensure that matrices d and e are positive values.

[0159] 2.2) According to constraints (19)-(22), establish the quadratic programming equation for passive optimization:

[0160]

[0161] -RΔx≤λ(24)

[0162] In the formula, Δx is the increment of parameter x, and x can be r n 、p n , d, e. Among them, λ=eig(Re(Z)), R=MRe(Z),

[0163] 2.3) Convert the above pole-residue form of internal impedance function (1) into the form of matrix multiplication as shown in (25);

[0164] Z(s)=C(sI-A) -1 B+D+sE(25)

[0165] Among them, Z(s) represents the form of converting Z(s) into a matrix, I is the N-order unit matrix, C=[r1 r2 … r N ], D = [d], E = [e], the matrix A is expressed as follows:

[0166]

[0167] 2.4) Perform positive definite processing on matrix D and matrix E, that is:

[0168]

[0169] Wherein, the elements in the matrix D and the matrix E include real numbers d and e respectively; T D is the normalized orthogonal matrix of D; D is the eigenvalue matrix of D; T E is the normalized orthogonal matrix of E; E is the eigenvalue matrix of E.

[0170] Λ D and Λ E Replace the negative eigenvalues ​​in with zero, update the matrix D and the matrix E, and get the matrix D mod and the matrix E mod ; 2.5) In formula (25), considering that matrices A and C are both complex matrices, the eigenvalue calculation is relatively complicated. In order to reduce the amount of calculation, this patent uses the Jordan decomposition theorem to simplify the complex matrices A and C into the Jordan standard form matrix A mod and the matrix C mod (The Jordan canonical form is a quasi-diagonal matrix consisting of several Jordan blocks with the main diagonal as eigenvalues, the lower (or upper) secondary diagonal all 1, and the rest all 0) arranged diagonally.)

[0171] 2.6) First, determine whether the eigenvalue of the real part of Z(s) at each frequency point without correction is greater than 0. If the requirement is met, it means that this system is a passive system. Otherwise, collect the frequency point information s that do not meet the conditions. w , and perform the following steps;

[0172] 2.7) The frequency information s collected in step 2.6) is w, arrange them in ascending order, and take the median between every two frequency points to establish a new frequency list s mod , and check the passivity at these frequencies. If the real part of Z(s) is mod The first frequency eigenvalue of is less than 0, then the starting frequency of passive optimization should be the starting frequency of the vector matching method. If the real part of Z(s) is in s mod If the last frequency eigenvalue of is less than 0, the end frequency of passive optimization should be the end frequency of vector matching method. Otherwise, the passive optimization interval should at least satisfy s w The interval satisfied;

[0173] 2.8) Use the quadratic programming algorithm to calculate the least squares solution that satisfies the constraints (23) and (24) until the correction parameters that satisfy the constraints are found, including the residue r n,mod , first-order term e mod , constant term d mod Otherwise, steps 2.4)-2.8 will be repeated. The internal impedance model after passive optimization correction should be:

[0174]

[0175] Among them, r n,mod =r n +Δr n , d mod =d+Δd,e mod =e+Δe, these correction quantities need to be solved by solving the quadratic programming of (19)-(22).

[0176] (3) Based on the “passivation of circuit network functions”, the circuit network model of the internal impedance of strong electromagnetic equipment can be given by using the impedance equivalent relationship between the internal impedance model Z(s) and the parameters R, L, and C.

[0177] 1) Constant term d and linear term e.

[0178] For the constant term and the first-order term of equation (28), the equivalent circuit in impedance form is as follows: Figure 1 As shown, the values ​​of parameters R and L are R=d and L=e respectively.

[0179] 2) The poles and residues are real numbers.

[0180] When the poles and residues of equation (28) are real numbers, an RC parallel circuit can be used to represent the equivalent circuit. The topological structure of the RC circuit is as follows: Figure 2 As shown. Among them:

[0181]

[0182] 3) The poles and residues are conjugate complex numbers.

[0183] When the poles and residues of equation (28) appear in the form of conjugate complex pairs, the corresponding equivalent circuit is a second-order circuit containing two energy storage elements. The topological structure of the equivalent circuit is as follows: Figure 3 As shown, the figure contains two resistor elements, one inductor element and one capacitor element. Figure 3 The values ​​of the circuit components are calculated according to formula (30).

[0184]

[0185] In the formula

[0186]

[0187] Therefore, if the circuit network function in equation (1) is regarded as an impedance parameter, it can be decomposed into three fractions. Each fraction is equivalent to a corresponding subcircuit. Finally, all subcircuits are connected in series to obtain the equivalent circuit corresponding to the network function, as shown in Figure 4 shown.

[0188] (4) The non-positive parameters of R, L, and C in the equivalent circuit are processed by "component passivity" using the particle swarm algorithm, and the non-positive parameters in the equivalent circuit are corrected to positive values ​​so that physical equivalence can be achieved using actual R, L, and C components.

[0189] Since there are some non-positive components in the fitting circuit obtained in the previous step, some methods are needed to select appropriate positive values ​​for these components, so that actual equivalent circuits can be built through capacitors, inductors and resistors in the later stage to achieve physical equivalence of wide-band impedance. Here, the particle swarm algorithm is selected to operate and construct a circuit network function containing these parameters to be optimized.

[0190] Particle swarm optimization (PSO) is a random search algorithm based on group collaboration developed by simulating the foraging behavior of bird flocks. In PSO, each solution to the optimization problem is called a "particle". All particles have a fitness value determined by the function being optimized, and each particle also has a speed that determines the direction and distance they fly. Then the particles follow the current optimal particle to search in the solution space.

[0191] 4.1) Set the objective function to

[0192]

[0193] Among them, ω i is the passive frequency of the component, M is the number of frequencies involved in the optimization, and Z mod is the rational function obtained by “network function passivation” in the first two steps, Z fitis the optimized resistance impedance function;

[0194] By randomly setting the initial values ​​of these non-positive components within the restricted range of R, L, and C components (i.e., particles), and using the particle swarm algorithm to find the minimum value of the objective function, the optimal solution that satisfies the "component passivity" parameter can be obtained.

[0195] 4.2) In order to balance the global exploration ability and local search ability, the inertia weight is set to

[0196]

[0197] Among them, ω max and ω min They are the maximum and minimum values ​​of the inertia weight, which are selected as 0.9 and 0.4 respectively in this patent.

[0198] 4.3) Set the asynchronously changing learning factor as shown in formula (34). In the initial stage of optimization iteration, the social learning ability of particles is weak and the self-learning ability is strong, which is convenient for fast search and enhances the global search ability. In the later stage of optimization iteration, the self-learning ability of particles is weak and the social learning ability is strong, which is conducive to local fine search and converges to the global optimal solution with higher accuracy.

[0199]

[0200] Among them, c 1s and c 2s are the initial values ​​of the learning factors c1 and c2, c 1e and c 2e are the termination values ​​of the learning factors c1 and c2 respectively. This patent selects c 1s =1,c 1e =2,c 2s =1.495,c 2e = 2. c1 and c2 are the learning factors of individuals and populations.

[0201] 4.4) PSO is initialized as a group of random particles (random solutions), and then finds the optimal solution through iteration. In each iteration, the particles update themselves by tracking two "extreme values". The first is the optimal solution found by the particle itself, which is called the individual extreme value pbest. The other extreme value is the optimal solution found by the entire population, which is the global extreme value gbest. When these two optimal values ​​are found, the particles update their speed and new position according to the following formula:

[0202]

[0203] x j =x j +v j(36)

[0204] Where j = 1, 2, ..., N max , N max is the total number of particles in this group, v j is the speed of each particle, rand() is a random number between (0,1), x j is the position of each particle, c1 and c2 are the learning factors of individuals and populations, ω PSO is the inertia weight.

[0205] 4.5) According to formulas (35) and (36), adjust the particle speed and position until the optimal position searched by the particle swarm is met and the predetermined minimum adaptation threshold is met or the end condition of the maximum number of iterations is reached; if the end condition is not met, continue the iteration.

[0206] A specific example is given below to help further illustrate the present invention.

[0207] First, an impedance analyzer is used to measure the impedance of a three-phase permanent magnet synchronous motor. The test frequency band is 100kHz to 110MHz. Since the test frequency band is wide, in order to ensure the fitting accuracy, this patent divides it into 100kHz-30MHz and 30Mz-110MHz for fitting. The modeling and simulation are carried out in sequence through the steps mentioned in this patent. The parameters of each component are shown in Table 1 and Table 2. The fitting results of the amplitude-frequency characteristics and phase-frequency characteristics of the equivalent circuit are shown in Figure 5 and Figure 6 As shown in the figure, it can be found that the equivalent circuit established by the vector matching method and the particle swarm algorithm has good amplitude and phase fitting results at each frequency point and has high accuracy.

[0208] Table 1 Circuit component parameters in the range of 1100kHz-30MHz

[0209] R(Ω) L(H) <![CDATA[R0(Ω)]]> <![CDATA[C0(F)]]> <![CDATA[R1(Ω)]]> <![CDATA[R2(Ω)]]> <![CDATA[C1(F)]]> <![CDATA[L1(H)]]> 7.7081 1.2313e-07 212102 3.9444e-10 302.544 1814.13 6.1132e-10 3.8952e-04 78826.1 0.9877 0.20719 45.4627 2.8106e-09 2.2177e-08

[0210] Table 2 Circuit component parameters in the range of 30Mz-110MHz

[0211] R(Ω) L(H) R1(Ω) R2(Ω) C1(F) L1(H) 2.8616 1.4337e-07 3.34157 992.6077 2.7505e-10 3.0614 1.35177e-05 43559.7367 0.4688 0.52968 0.82877 64764.485 1.1434e-09 1.05e-09

[0212] In the semi-physical test of electromagnetic compatibility for integrated power systems, strong electromagnetic equipment such as high-power motors must be replaced as loads. To this end, it is necessary to model and physically simulate strong electromagnetic equipment such as motors. The present invention proposes an equivalent physical simulation method for the wide-band impedance of strong electromagnetic equipment based on circuit networks. It only needs to know the impedance characteristics of strong electromagnetic equipment such as motors to equate loads such as motors to a "black box" circuit network, and solves the common "passivity" in circuit fitting and the negative value of components in the circuit network. The problem ensures the stable simulation of the interconnected network, and can be directly used for the equivalent physical simulation of the impedance of strong electromagnetic equipment. It can be used for semi-physical verification tests of electromagnetic compatibility of power systems in scenes such as ships, as well as modeling and model prediction of other electromagnetic equipment that is difficult to obtain specific circuit parameters, providing prediction models and data support for the prediction of conducted interference of strong electromagnetic equipment and power systems and targeted electromagnetic interference suppression.

[0213] Finally, it should be noted that the above embodiments are only used to illustrate the technical solution of the present invention rather than to limit it. Although the applicant has described the present invention in detail with reference to the preferred embodiments, ordinary technicians in the field should understand that modifications or equivalent substitutions of the technical solution of the present invention without departing from the purpose and scope of the technical solution should be included in the scope of the claims of the present invention.

Claims

1. A physically implementable method for equivalent modeling of broadband impedance of strong electromagnetic equipment, characterized by: performing the following steps: 1) According to the vector matching method, the internal impedance of the strong electromagnetic device to be impedance equivalent is expressed in the form of a rational function; 2) For the rational function obtained in step 1), the network port is subjected to passivation processing by using the residue perturbation method to obtain the circuit network function; 3) The circuit network function obtained by passivation is converted into the corresponding equivalent circuit using the impedance equivalent relationship between the rational function and the parameters R, L, and C. 4) For the non-positive parameters of R, L, and C in the equivalent circuit, the particle swarm algorithm is used to correct the non-positive parameters in the equivalent circuit to positive values, that is, to obtain the final equivalent circuit so that the actual R, L, and C components can be used for physical equivalence; In step 1), the rational function representing the internal impedance of the strong electromagnetic device is: In the formula, the residue r n and the extreme point p n is a real number or a conjugate complex number pair; d and e are real numbers; Z(s) is the internal impedance of a strong electromagnetic device that changes with frequency; s indicates that the function is a complex frequency domain function, n = 1, 2, ..., N, N is the order of the polynomial (1); Based on the wide frequency band test data Z(s k ), s k (k=1,2,…,p) is the test frequency, p is the number of test frequency; determine the residue r in the rational function by iteration n , extreme point p n , real number d, real number e are as follows: 1.1) Remember is the initial pole, and the auxiliary function σ(s) is introduced, as shown in formula (2); Multiply the auxiliary function σ(s) by the rational function Z(s) of the internal impedance of the strong electromagnetic device, and assume that during the iteration, σZ(s) and σ(s) have the same poles. In the formula, is the residue of the auxiliary function σ(s); 1.2) Using equations (2) and (3), establish the parameter r n , The linear equation of d and e is: 1.3) The test data Z(s) k ) and the corresponding frequency s kSubstituting into equation (4), we get the linear equation, and the matrix form of the linear equation is as follows: Ax=b (5) Among them, the coefficient matrix A P×(2N+2) It is expressed as follows: The vector x to be solved and the known vector b are as follows: b=[Z(s1) Z(s2)…Z(s k )…Z(s p )] T (8) 1.4) Solve the linear equation matrix (5) using the least squares method to obtain r n , d and e; 1.5) Rewrite the above function σZ(s) and auxiliary function σ(s) into zero-pole form, that is: In the formula, z n and are the zeros of σZ(s) and σ(s) respectively; 1.6) According to formulas (9)-(10), the rational function Z(s) of the internal impedance of the strong electromagnetic device is updated to obtain: From (11), we can see that the pole of the function Z(s) is the zero of the auxiliary function σ(s), and the zero of the auxiliary function σ(s) is the initial pole of the next iteration Z(s); 1.7) According to equation (2) and equation (10), equation (12) is established, that is: Writing the right side of equation (12) in the form of a determinant yields: 1.8) Change the determinant H on the right side of equation (13) into the following form: In the formula, P is the An N-order diagonal matrix, U is a 1×N-dimensional column vector with elements of 1, is the residue containing σ(s) 1×N order row vector; In formula (14), the diagonal matrix P, the row vector They are as follows: 1.9) Since the characteristic equation of matrix H is |sE-H|=0, E is the unit determinant, the matrix H is replaced by equation (14) combined with equations (15) and (16), and then becomes |PR|=0; further using the matrix R formed by the residues of equation (1), as shown in equation (17); according to equation (13) and eig|PR|=0, the zero point of the auxiliary function σ(s) is obtained That is, get the initial value of the next iteration 1.10) Repeat steps 1.4) to 1.9) to obtain r n , d, e, and Determine whether the data error of the function Z(s) obtained in two consecutive iterations is less than the threshold. If so, end the iteration. Otherwise, define the latest pole as the initial pole. Return to step 1.4) and recalculate until the required r is obtained. n 、p n , d and e; Step 2) The steps of the passivation process are as follows: 2.1) Establish the following constraint equation for network passivity: That is, all eigenvalues ​​of the real part of the impedance matrix of the model are positive for all frequencies; the following constraint equation is established: eig{Re(Z+ΔZ)}>0 (20) eig(d+Δd)>0 (21) eig(e+Δe)>0 (22) ΔX represents the correction deviation added to X as the quantity to be determined, X is the unified representation of the required correction quantity in equations (19)-(22), and ΔX is an unknown quantity that needs to be solved by establishing a quadratic programming equation; 2.2) According to constraints (19)-(22), establish the quadratic programming equation for passive optimization: -RΔx≤λ (24) In the formula, Δx is the increment of parameter x, and x is r n 、p n , d, e; λ=eig(Re(Z)), R=MRe(Z), 2.3) Convert the rational function formula (1) into the form of matrix multiplication as shown in formula (25); Z(s)=C(sI-A) -1 B+D+sE (25) Among them, Z(s) represents the form of converting Z(s) into a matrix, I is the N-order unit matrix, C=[r1r2…r N ], D = [d], E = [e], the matrix A is expressed as follows: 2.4) Perform positive definite processing on matrix D and matrix E, that is: Wherein, the elements in the matrix D and the matrix E include real numbers d and e respectively; T D is the normalized orthogonal matrix of D; D is the eigenvalue matrix of D; T E is the normalized orthogonal matrix of E; E is the eigenvalue matrix of E; Λ D and Λ E Replace the negative eigenvalues ​​in with zero, update the matrix D and the matrix E, and get the matrix D mod and the matrix E mod ; 2.5) In equation (25), the complex matrices A and C are simplified to the Jordan standard form matrix A mod and the matrix C mod ; 2.6) Determine whether the eigenvalue of the real part of Z(s) at each frequency point without correction is greater than 0. If so, it is a passive system; otherwise, collect the frequency point information s that does not meet the conditions w , and perform the following steps; 2.7) The frequency information s collected in step 2.6) is w , arrange them in ascending order, and take the median between every two frequency points to establish a new frequency list s mod , and check the passivity at these frequencies. If the real part of Z(s) is mod The first frequency eigenvalue of is less than 0, then the starting frequency of passive optimization is the starting frequency of the vector matching method. If the real part of Z(s) is in s mod If the last frequency eigenvalue of is less than 0, the termination frequency of passive optimization is the termination frequency of vector matching method. Otherwise, the passive optimization interval should at least satisfy s w The interval satisfied; 2.8) Use the quadratic programming algorithm to calculate the least squares solution that satisfies the constraints (23) and (24) until the correction parameters that satisfy the constraints are found, including the residue r n,mod , first-order term e mod , constant term d mod Otherwise, continue to repeat steps 2.4)-2.8); the rational function of the internal impedance after passivation is Among them, r n,mod =r n +Δr n , d mod =d+Δd,e mod =e+Δe, these correction values ​​are solved by quadratic programming of equations (19)-(22).

2. According to a physically implementable method for modeling strong electromagnetic equipment wide-band impedance equivalents according to claim 1, it is characterized in that: Step 3) converting into a corresponding equivalent circuit as follows; 3.1) For the constant term d and the linear term e of a rational function, the equivalent circuit in impedance form is a series circuit of R and L, and the values ​​of the parameters R and L are R = d and L = e respectively; 3.2) When the poles and residues of a rational function are real numbers, an RC parallel circuit is used to represent the equivalent, where: 3.3) When the poles and residues of a rational function are conjugate complex numbers, the corresponding equivalent circuit is a second-order circuit, which includes two resistors, one inductor and one capacitor, wherein one resistor and one inductor are connected in series and then connected in parallel with the capacitor and another resistor. The values ​​of each circuit element are calculated according to formula (30); In the formula 3.4) Decompose equation (1) into three fractions (3.1), 3.2) and 3.3), each fraction is equivalent to a corresponding subcircuit, and finally all subcircuits are connected in series to obtain the equivalent circuit corresponding to the network function.

3. According to a physically implementable method for strong electromagnetic equipment wide-band impedance equivalent modeling as described in claim 2, it is characterized in that: Step 4) The method of using the particle swarm algorithm to correct the non-positive parameters in the equivalent circuit to positive values ​​is: 4.1) Set the objective function to Among them, ω i is the passive frequency of the component, M is the number of frequencies involved in the optimization, and Z mod is the rational function obtained by passivation, that is, the rational function of internal impedance (28), Z fit is the optimized resistance impedance function; 4.2) Set the inertia weight to Among them, ω max and ω min are the maximum and minimum values ​​of the inertia weight respectively; 4.3) Set the learning factor of the asynchronous change as shown in formula (34), Among them, c 1s and c 2s are the initial values ​​of the individual and population learning factors c1 and c2, respectively, c 1e and c 2e are the termination values ​​of the learning factors c1 and c2 respectively; 4.4) PSO is initialized as a group of random particles, and then finds the optimal solution through iteration. In each iteration, the particles update themselves by tracking two extreme values; the first extreme value is the optimal solution found by the particle itself, which is called the individual extreme value pbest, and the other extreme value is the optimal solution found by the entire population, which is the global extreme value gbest; when these two extreme values ​​are found, the particles update their speed and new position according to the following formula: x j =x j +v j (36) Where j = 1, 2, ..., N max , N max is the total number of particles in this group, v j is the speed of each particle, rand() is a random number between (0,1), x j is the position of each particle, c1 and c2 are the learning factors of individuals and populations, ω PSO is the inertia factor; 4.5) According to formulas (35) and (36), the particle speed and position are adjusted until the optimal position searched by the particle swarm is met and the predetermined minimum adaptation threshold is met or the end condition of the maximum number of iterations is reached; if the end condition is not met, the iteration is continued.

Images