Complex magnetic field excitation-oriented electrical steel dynamic magnetic field intensity calculation method

Abstract

The invention discloses a complex magnetic field excitation-oriented electrical steel dynamic magnetic field intensity calculation method, which comprises the following steps of: firstly, performing infinite interpolation expansion on the basis of a small number of actually measured symmetric hysteresis loops to reduce the experimental data volume required by static Preisach hysteresis model identification; an improved dynamic modeling method is provided, model parameters related to the magnetization rotation point and the magnetic bias quantity are introduced, and the method has higher adaptability compared with a traditional model under complex excitation. An improved particle swarm optimization algorithm is adopted to optimize parameters in the dynamic Preisach hysteresis model, a population eradication and loop initialization mechanism is introduced into the improved particle swarm optimization algorithm, the optimization precision is obviously improved compared with that of a classical particle swarm algorithm, and more efficient model parameter identification and more accurate dynamic hysteresis loop simulation are achieved; based on electrical steel hysteresis loop test data under various excitation conditions, it is more comprehensively verified that the provided modeling method and parameter optimization algorithm have better comprehensive performance.

Description

Technical Field

[0001] This invention relates to the fields of magnetic materials and electromagnetic calculation technology, specifically to a method for calculating the dynamic magnetic field strength of electrical steel under complex magnetic field excitation. Background Technology

[0002] As core functional materials for key equipment such as motor drive systems, high-frequency transformers, and electromagnetic sensors, the accurate modeling and parameter identification of the dynamic hysteresis characteristics of magnetic materials have always been a key research direction in the field of electrical engineering. With the rapid development of technologies such as new energy power generation, electric vehicles, and high-speed rail transit, electrical equipment is evolving towards higher frequencies and higher power densities. Complex magnetic field excitation scenarios (such as high-frequency pulse width modulation, multiple harmonic superposition, and DC bias) pose unprecedented challenges to the modeling of magnetic materials.

[0003] The dynamic hysteresis modeling method proposed in the literature (Zhao X, Yang L, Wang R, et al. A dynamic model for prediction of hysteresis and loss property of soft material under multiple harmonic magnetizations[J].IEEE Transactions on Industry Applications,2022,58(6):7163-7170.) treats the abnormal loss coefficient as a function dependent on the degree of material saturation, which improves the accuracy of iron loss prediction under dynamic magnetization conditions to a certain extent. The "Dynamic Hysteresis Modeling Method for Iron Core Based on Permeability-Capacitance Analogy and Analytical Preisach Model" (Three Gorges University.202311265657.6[P].2023-12-26.) uses the classical field separation method to characterize the dynamic permeability of the iron core and introduces the fractional derivative of RL to correct the eddy current loss expression. The proposed method can be used for electromagnetic coupling calculation of high-frequency transformer cores. However, existing methods have failed to fully consider the nonlinear coupling of eddy current losses, abnormal losses and hysteresis effects, resulting in large errors in the magnetic field prediction of the model under complex excitation conditions. This directly causes problems such as distortion in the calculation of motor iron losses and deviation in the prediction of transformer temperature rise, which seriously restricts the design accuracy and operational reliability of high power density electrical equipment.

[0004] At the parameter identification level, existing methods usually require a large number of repetitive experiments to obtain hysteresis loop data in a wide frequency domain and multiple excitation modes. The experimental cycle is long and the cost is high, which makes it difficult to meet the design requirements of rapid iteration of power electronic equipment. In addition, most traditional optimization algorithms are prone to getting trapped in local optima during the model parameter identification process, which makes it difficult to meet the global simulation accuracy of the model under complex excitation conditions. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a method for calculating the dynamic magnetic field strength of electrical steel under complex magnetic field excitation. This invention addresses the problem of identifying dynamic Preisach hysteresis model parameters based on limited experimental data by focusing on two aspects: dynamic hysteresis modeling methods and parameter optimization algorithms. This reduces the reliance of model identification on experimental data, improves the modeling mechanism to enhance applicability, and refines traditional optimization algorithms to improve optimization accuracy. Thus, it achieves accurate simulation and prediction of dynamic hysteresis loops and iron losses under complex excitation conditions, enabling accurate simulation of the dynamic hysteresis loops of electrical steel.

[0006] To address the above problems, this invention provides a method for calculating the dynamic magnetic field strength of electrical steel under complex magnetic field excitation, characterized by the following steps:

[0007] Step 1: Under quasi-static symmetrical alternating excitation test conditions, collect no more than 10 quasi-static measured hysteresis loop data from a single sample of the target electrical steel. Then, using an infinite interpolation expansion method, expand the number of hysteresis loops based on the collected quasi-static measured hysteresis loops, obtaining a total of no less than 29 hysteresis loop data. Based on this hysteresis loop data, derive the functional expression of the Everett function and establish a static Preisach hysteresis model based on the Everett function.

[0008] Step 2: Design a dynamic Preisach hysteresis model based on the static Preisach hysteresis model;

[0009] The output value of the dynamic Preisach hysteresis model is calculated as follows:

[0010]

[0011] γ=γ L (B L (t))×γ G (B G (t)) (12)

[0012] B G (t)=δB * (t) / B m (13)

[0013] B L (t)=(B * (t)-B rev (t)) / (B dc +δ(|B m |-|B dc |)-B rev(t)) (14)

[0014] Among them, B * (t) represents the input sequence of this dynamic hysteresis model, H hy (B * (t) represents the output value of the static Preisach hysteresis model established based on the method described in step 1; σ represents the conductivity of the target electrical steel monolithic sample, d represents the thickness of the target electrical steel monolithic sample, and C0 represents the parameter to be identified; B L (t) and B G (t) represents the normalized values ​​of local and global magnetic flux density, respectively, and serves as two intermediate quantities; B rev (t) represents the magnetic flux density waveform B. * (t) is the extreme value before time t and closest to time t; B dc For B * The magnetic flux density bias value of (t), B m For B * The amplitude of (t), δ=sign(dB) * (t) / dt), sign() is the sign function, which is used to determine the sign of a number. It returns 1 for positive numbers, 0 for zero, and -1 for negative numbers;

[0015] Step 3: Construct the function γ in equation (12) using a quadratic polynomial. L () and γ G (), such that a0=γ L (0), a2=γ L (1) and b0 = γ G (-1),b1=γ G (0), b2=γ G (1); The function γ can be derived. L () and γ G The expression for () is:

[0016] γ L (x)=(2a0-4a1+2a2)x 2 +(-3a0+4a1-a2)x+a0 (15)

[0017]

[0018] Among them, a0, a1, a2, b0, b1, and b2 are 6 model parameters to be identified, and together with C0 in equation (11), there are a total of 7 parameters to be identified.

[0019] Step 4: Optimize the 7 parameters to be identified in Step 3 using an improved particle swarm optimization algorithm;

[0020] The improved particle swarm optimization algorithm combines the traditional particle swarm optimization (PSO) algorithm with population extinction and cyclic initialization mechanisms. The specific process is as follows:

[0021] Set the population size N. n The maximum number of iterations N in a single round of the PSO algorithm it And the maximum number of population extinctions N ex The optimization interval for C0 is set to [0, 0.5], and the optimization intervals for a0, a1, a2, b0, b1, and b2 are all set to [0.5, 1.5]. A random assignment method within the optimization interval is used to perform N-value calculations on the seven parameters to be identified (a0, a1, a2, b0, b1, b2, and C0). n The initialization yields N n The initial population of particles; based on the values ​​of a0, a1, a2, b0, b1, b2 of each particle, γ is obtained by combining equations (15) and (16). L (), γ G The function formula of (), then γ L (), γ G Substituting the function formula of () and the corresponding value of C0 into formula (11), and using Nl measured hysteresis loops under dynamic alternating excitation conditions as the identification group, the corresponding magnetic flux density sampling sequence B is extracted. * Substituting (t) into (11) as input, we obtain the magnetic field strength sequence H(t) under the corresponding excitation condition, and apply this sequence to the calculation of the objective function;

[0022] Set the objective function as follows:

[0023]

[0024] Where Nl is the number of measured hysteresis loops in the identification group, Nl is no greater than 6, l is the hysteresis loop symbol in the identification group; Np is the number of discrete points in a single cycle on a hysteresis loop. and These represent the magnetic field strength at the p-th discrete time point in the l-th measured and model-simulated hysteresis loop, respectively. The input sequence B is calculated using formula (11). * (t) represents the magnetic flux density sequence corresponding to the l-th measured hysteresis loop in the identification group; |·| indicates the sign of the absolute value.

[0025] According to formula (17), calculate the objective function Y for each particle in the initial population. FitThe position of the particle with the smallest objective function value is the optimal individual position and the global historical optimal individual position during the first iteration update, used to guide the first iteration update of the population. After the first iteration update is completed, the objective function value of each particle in the population that has completed the update is calculated, and the particle with the smallest objective function value is selected as the optimal individual in the population that has completed the first update. The objective function value of the optimal individual in the population that has completed the first update is compared with the objective function value of the optimal individual in the initial population, and the individual with the smallest objective function value is selected as the global historical optimal individual. After each iteration update is completed, the optimal individual in the current population is obtained by calculating the objective function value, and this optimal individual is compared with the global historical optimal individual in the previous iteration update process, and the individual with the smallest objective function value is selected as the current global historical optimal individual. The position of the current optimal individual in the population and the position of the current global historical optimal individual are used to guide the next iteration update of the population. The specific process is as follows:

[0026] Let the position coordinate vector of the k-th particle in the population be represented by X. k =(X k1 ,X k2 ,…,X k7 ), X k Each element represents the values ​​of a0, a1, a2, b0, b1, b2, and C0 of the particle in the optimization space; during the optimization process, the position update method of the k-th particle is shown in equation (18);

[0027] X k (T+1)=ω×X k (T)+c1×r×X P (T)+c2×r×X G (T) (18)

[0028] Where T represents the number of iterations completed so far, ω, c1, and c2 are the weight factor, individual learning rate, and group learning rate, respectively, and are set to 0.5, 1.5, and 1.5; r is a random number in the interval [0,1], and X... k (T) represents the individual position of particle k at the time of the T-th iteration update, X P (T), X G (T) represents the best individual position in the population and the global historical best individual position at the time of the Tth iteration update, respectively;

[0029] The position of each particle in the current population is iteratively updated according to formula (18). The values ​​of parameters a0, a1, a2, b0, b1, b2, and C0 corresponding to the particle are determined according to the particle position. Then, the objective function value of each updated particle is calculated. The particle with the smallest objective function value is taken as the optimal individual in the population when the (T+1)th iteration update is completed. At the same time, the position X of the particle is set. P The objective function value corresponding to (T+1) and the global historical best individual position X G The objective function values ​​corresponding to (T) are compared, and the position of the particle with the smallest objective function value is taken as the globally historical best individual position X when the (T+1)th iteration update is completed. G (T+1); then increment the iteration count by 1, repeat the process of iterating and updating each particle in the current population according to formula (18), and repeat the process of obtaining the best individual position in the population and the global historical best individual position when the current iteration update is completed, and repeat this process until the iteration update count reaches the maximum iteration count N. it This means completing one round of PSO algorithm computation.

[0030] When one round of the PSO algorithm operation is completed, after completing N... it In the next updated population, N was selected using a roulette wheel method. n Ten surviving particles directly become part of the initial population for the next round; this process is called population extinction. The remaining N... n -N n The initial 10 individuals are randomly generated within the optimization interval and are called new individuals. During initialization, it is ensured that the normalized Euclidean distance from all new individuals to any surviving individual is greater than 0.05 to guarantee the diversity of the initial population. This process is called cyclic initialization. After the particle initialization is completed, the current round of PSO algorithm calculation continues. When the current round of PSO algorithm calculation is completed, the above population extinction and cyclic initialization process continues, and the next round of PSO algorithm calculation is performed, repeating continuously. This process continues until the maximum number of population extinctions, N, is reached. ex At that time, the final loop initialization is performed, and the final round of PSO algorithm operation is performed. After the final round of PSO algorithm operation is completed, the global historical best individual position of the entire population is the optimal solution, that is, the optimal values ​​of the 7 parameters to be identified are obtained.

[0031] Step 5: Substitute the optimal values ​​of the 7 parameters to be identified obtained in Step 4 into Formula (11) to obtain the optimal dynamic hysteresis model; obtain the magnetic flux density waveform of the dynamic alternating excitation condition to be applied to the target electrical steel single sample in Step 1, and then input the magnetic flux density waveform into the above-mentioned optimal dynamic hysteresis model to obtain the magnetic field strength waveform sequence of the target electrical steel single sample under the dynamic alternating excitation condition.

[0032] Compared with existing technologies, the advantages of this invention are as follows: This invention provides a method for calculating the dynamic magnetic field strength of electrical steel under complex magnetic field excitation. Firstly, it uses an infinite interpolation expansion method based on a small number of measured symmetrical hysteresis loops, reducing the amount of experimental data required for identifying the static Preisach hysteresis model. Secondly, the proposed improved dynamic modeling method introduces the magnetization rotation point and bias magnetization, exhibiting stronger adaptability under complex excitations compared to traditional models. Thirdly, an improved particle swarm optimization algorithm is used to optimize the parameters in the dynamic Preisach hysteresis model. This improved algorithm introduces population extinction and cyclic initialization mechanisms, significantly improving the optimization accuracy compared to the classic particle swarm optimization algorithm, achieving more efficient model parameter identification and more accurate dynamic hysteresis loop simulation. Finally, based on test data of electrical steel hysteresis loops under various excitation conditions, the proposed modeling method and parameter optimization algorithm are more comprehensively verified to have superior overall performance. Attached Figure Description

[0033] The exemplary embodiments of the present invention can be more fully understood by referring to the accompanying drawings.

[0034] Figure 1 This is a schematic diagram of the operation process of an embodiment of the method for calculating the dynamic magnetic field strength of electrical steel for complex magnetic field excitation according to the present invention.

[0035] Figure 2 This is a comparison of the simulation results of a static Preisach hysteresis model before and after employing the hysteresis loop extension method in one embodiment; wherein, Figure 2 Figure (1) shows the simulation results of the static Preisach hysteresis model before the hysteresis loop extension method. Figure 2 Figure (2) shows the simulation effect of the static Preisach hysteresis model after adopting the hysteresis loop expansion method.

[0036] Figure 3 This is a comparison chart of the simulation effects of a conventional embodiment and the improved dynamic Preisach hysteresis model of the present invention. Figure 3 Figure (1) shows the simulation results using the traditional dynamic Preisach hysteresis model. Figure 3 Figure (2) shows the simulation effect of the improved dynamic Preisach hysteresis model of the present invention.

[0037] Figure 4 This is a schematic diagram comparing the convergence of a conventional particle swarm optimization algorithm (PSO) with the improved particle swarm optimization algorithm (EPSO) of this invention in one embodiment.

[0038] Figure 5 This is a comparison chart of the prediction results of the optimal dynamic hysteresis model obtained by using the traditional Particle Swarm Optimization (PSO) algorithm and the improved Particle Swarm Optimization (EPSO) algorithm of this invention, according to one embodiment; wherein, Figure 5 Figure (1) shows the prediction results of the optimal dynamic hysteresis model obtained by the traditional particle swarm optimization algorithm (PSO). Figure 5 (2) is the prediction effect diagram of the optimal dynamic hysteresis model obtained by using the improved particle swarm optimization algorithm (EPSO) of this invention.

[0039] Figure 6 This is a binarized Boolean plot of the goodness-of-fit of the dynamic hysteresis loop under different harmonic excitation conditions using an optimal dynamic hysteresis model of an embodiment of the present invention for calculating the dynamic magnetic field strength of electrical steel under complex magnetic field excitation. Detailed Implementation

[0040] This invention provides a method for calculating the dynamic magnetic field strength of electrical steel under complex magnetic field excitation, the method specifically including the following steps:

[0041] Step 1: Under quasi-static symmetrical alternating excitation test conditions, collect 8 to 10 quasi-static measured hysteresis loop data points from a single sample of the target electrical steel. Then, using an infinite interpolation expansion method, expand the number of hysteresis loops based on the collected quasi-static measured hysteresis loops, obtaining a total of no less than 29 hysteresis loop data points. Based on this hysteresis loop data, derive the Everett function expression and establish a static Preisach hysteresis model based on the Everett function.

[0042] Interpolation augmentation of hysteresis loop data is achieved through the basic magnetization curve and weighting coefficients:

[0043]

[0044] Wherein, the superscript i represents the sequence number of the measured hysteresis loop data, and f() is the basic magnetization curve function. and represents the flux density amplitude of the i-th and (i+1)-th measured hysteresis loops, respectively. and represents the magnitude of the magnetic field strength of the i-th and i+1-th measured hysteresis loops, respectively. and These represent the amplitudes of the magnetic flux density and magnetic field strength of the expanded hysteresis loop, respectively. f() is obtained by solving for all the vertices of the measured hysteresis loops using piecewise cubic Hermite interpolation, and is essentially a single-valued mapping from B to H. and These represent the weighting coefficients of the magnetic flux density (B) waveform and the magnetic field strength (H) waveform of the extended hysteresis loop, respectively.

[0045]

[0046] Where t represents a discrete sampling time point. and Let B(t) be the B-waveform sequence and H-waveform sequence of the extended hysteresis loop. i B(t) i+1 Let H(t) represent the B-wavelength sequences of the i-th and (i+1)-th measured hysteresis loops, respectively. i H(t) i+1 These represent the H-waveform sequences of the i-th and i+1-th measured hysteresis loops, respectively.

[0047] Based on the measured and extended hysteresis loops, the Everett function of the static Preisach hysteresis model is constructed as follows:

[0048]

[0049] Here, the bivariate function E(α,β) represents the Everett function. α () represents the measured or extended hysteresis loop descending branch function with magnetic flux density amplitude α, and is a single-valued mapping from B to H. α can take the magnetic flux density amplitude of all measured and extended hysteresis loops. For a specific α, β can take any value satisfying equation (6):

[0050] -α≤β≤α (6)

[0051] After obtaining the scatter set (α, β, E(α, β)), the symmetry of the Everett function along the line α = -β is used to obtain a supplementary scatter set (-β, -α, E(α, β)). Linear interpolation is then performed using the two scatter sets to obtain the Everett function expression E(α, β) with α and β as variables.

[0052] The quasi-static measured hysteresis loop is a symmetrical hysteresis loop, therefore H α() represents the descent branch function of the symmetric hysteresis loop with a magnetic flux density amplitude of α. As can be seen from equation (5), a greater number of extended loops means a denser number of Everett function data points. The interpolation extension method can supplement the amount of data when there is only a small amount of measured data, thus ensuring the accuracy of the static Preisach hysteresis model.

[0053] For a given magnetic flux density input sequence B(t), the output value of the static Preisach hysteresis model based on the Everett function is calculated by (7):

[0054]

[0055] In the formula, B(t) is the magnetic flux density input sequence of the static Preisach hysteresis model, and H hy (B(t)) represents the output value of the model at time t; N REV Let E1 be the total number of inversion points of the Everett function at time t, and let E1 = E(α1, β1); H hy (B(t-Δt)) is the output value of the model at time t-Δt, where Δt is the discrete interval of the time sampling points, determined by the experimental testing system of the hysteresis loop, and is a constant; E ∑+ and E ∑- Two intermediate quantities for calculation are defined by (8) and (9) respectively:

[0056]

[0057] In the formula, (α) n ,β n ) represents the coordinates of the nth inversion point of the Everett function.

[0058] At time t, α n ,β n The method for determining it is as follows: Assume that at time t and before, the B(t) sequence has a total of N occurrences. br Take a number of extreme values ​​and arrange them in chronological order to form an extreme value sequence. Add an element B(0) before the first element of the extreme value sequence, which is the first value in the B(t) sequence. If the value of this element is less than zero, it is recorded as a minimum value; otherwise, it is recorded as a maximum value.

[0059] Check the Kth sequence br Is element K in the Kth order? br -1 and the Kth br Within the interval defined by -2 element values, if not, delete the Kth element in the extreme value sequence. br -1 and the Kth br -2 elements are used to form a new extreme value sequence. Then, the Kth element in the new sequence is... brRepeat the same operation on each element until the sequence satisfies the Kth element. br The element is always in the Kth position. br -1 and the Kth br The sequence continues until the value of the first element is within the range defined by -2 element values, or until the number of elements in the sequence is less than 3; then, it checks whether the absolute value of the first element of the sequence is the maximum absolute value among all elements in the sequence. If not, the element is deleted; during this process, K br It can be any value greater than 2 and not greater than the total length of the latest extreme value sequence; this process is called erasure.

[0060] After erasure, the number of elements in the extreme value sequence is N. REV For the first element, if it is a maximum value, then denote its value as α1 and its negative as β1; if it is a minimum value, then denote its value as β1 and its negative as α1; for the nth element, 1 <n≤N REV If it is a maximum value, then the value of this element is denoted as α. n Let the value of the (n-1)th element be β. n If it is a local minimum, then the value of this element is β. n Let the value of the (n-1)th element be α. n Thus, a total of N can be obtained. REV Everett function inverts point coordinates (α) n ,β n ).

[0061] Step 2: Design a dynamic Preisach hysteresis model based on the static Preisach hysteresis model.

[0062] Traditional dynamic modeling methods treat the total magnetic field strength as a superposition of three components: hysteresis, eddy current, and anomaly.

[0063]

[0064] In the formula, σ is the conductivity of the target electrical steel single-piece sample, d is the thickness of the target electrical steel single-piece sample, C0 is the parameter to be identified, and H... hy is the hysteresis component of the magnetic field strength.

[0065] The dynamic model designed in this invention additionally considers the influence of magnetization rotation point and bias, making it applicable to complex excitation conditions. The output value calculation method of this dynamic Preisach hysteresis model is as follows:

[0066]

[0067] γ=γ L (B L (t))×γ G(B G (t)) (12)

[0068] B G (t)=δB * (t) / B m (13)

[0069] B L (t)=(B * (t)-B rev (t)) / (B dc +δ(|B m |-|B dc |)-B rev (t)) (14)

[0070] Among them, B * (t) represents the input sequence of this dynamic hysteresis model, H hy (B * (t) represents the output value of the static Preisach hysteresis model established based on the method described in step 1. σ is the conductivity of the target electrical steel monolithic sample, d is the thickness of the target electrical steel monolithic sample, and C0 is the parameter to be identified; B L (t) and B G (t) represents the normalized values ​​of local and global magnetic flux density, respectively, and serves as two intermediate quantities. B rev (t) represents the magnetic flux density waveform B. * (t) is the extreme value of t before time t and closest to time t. B dc For B * The magnetic flux density bias value of (t), B m For B * The amplitude of (t), δ=sign(dB) * (t) / dt).

[0071] As can be seen from equations (11), (12), (13) and (14), the dynamic method of the present invention utilizes the information of magnetization rotation point and bias magnetic quantity, and considers the coupling relationship between hysteresis nonlinearity and eddy current effect. Therefore, it is more adaptable and more accurate than the traditional model.

[0072] Step 3: Construct the function γ in equation (12) using a quadratic polynomial. L () and γ G (), such that a0=γ L (0), a2=γ L (1) and b0 = γ G (-1),b1=γ G (0), b2=γ G (1). The function γ can be derived. L() and γ G The expression for () is:

[0073] γ L (x)=(2a0-4a1+2a2)x 2 +(-3a0+4a1-a2)x+a0 (15)

[0074]

[0075] Among them, a0, a1, a2, b0, b1, b2 are 6 model parameters to be identified, and together with C0 in equation (11), there are a total of 7 parameters to be identified.

[0076] Step 4: Optimize the 7 parameters to be identified in Step 3 using an improved Particle Swarm Optimization (EPSO) algorithm.

[0077] The improved particle swarm optimization algorithm introduces population extinction and cyclic initialization mechanisms, which can fully utilize initialization to increase population diversity and improve the probability of escaping local optima. This improved particle swarm optimization algorithm combines the traditional particle swarm optimization algorithm with population extinction and cyclic initialization mechanisms. The specific process is as follows:

[0078] Set the population size N. n The maximum number of iterations N in a single round of the PSO algorithm (the traditional particle swarm optimization algorithm) it And the maximum number of population extinctions N ex The optimization interval for C0 is set to [0, 0.5], and the optimization intervals for a0, a1, a2, b0, b1, and b2 are all set to [0.5, 1.5]. A random assignment method within the optimization interval is used to perform N-value calculations on the seven parameters to be identified (a0, a1, a2, b0, b1, b2, and C0). n The initialization yields N n The initial population of particles. Based on the values ​​of a0, a1, a2, b0, b1, b2 of each particle, γ is obtained by combining equations (15) and (16). L (), γ G The function formula of (), then γ L (), γ G Substituting the function formula of () and the corresponding value of C0 into formula (11), and using 4-5 measured hysteresis loops under dynamic alternating excitation conditions as the identification group, the corresponding magnetic flux density sampling sequence B is extracted. * Substituting (t) into (11) as input, we obtain the magnetic field strength sequence H(t) under the corresponding excitation conditions, and apply this sequence to the calculation of the objective function.

[0079] Set the objective function as follows:

[0080]

[0081] Where Nl is the number of measured hysteresis loops in the identification group, and l is the symbol of the hysteresis loop in the identification group. Np is the number of discrete points in a single cycle on a hysteresis loop. and These represent the magnetic field strength at the p-th discrete time point in the l-th measured and model-simulated hysteresis loop, respectively. The input sequence B is calculated using formula (11). * (t) represents the magnetic flux density sequence corresponding to the l-th measured hysteresis loop in the identification group; |·| indicates the absolute value sign. In this example, Nl is 4, and the excitation conditions of the corresponding measured hysteresis loop are shown in Table 2.

[0082] According to formula (17), calculate the objective function Y for each particle in the initial population. Fit The position of the particle with the smallest objective function value is the optimal individual position and the global historical optimal individual position during the first iteration update, used to guide the first iteration update of the population. After the first iteration update is completed, the objective function value of each particle in the population that has completed the update is calculated, and the particle with the smallest objective function value is selected as the optimal individual in the population that has completed the first update. The objective function value of the optimal individual in the population that has completed the first update is compared with the objective function value of the optimal individual in the initial population, and the individual with the smallest objective function value is selected as the global historical optimal individual. After each iteration update is completed, the optimal individual in the current population is obtained by calculating the objective function value, and this optimal individual is compared with the global historical optimal individual in the previous iteration update process, and the individual with the smallest objective function value is selected as the current global historical optimal individual. The position of the current optimal individual in the population and the position of the current global historical optimal individual are used to guide the next iteration update of the population. The specific process is as follows:

[0083] Let the position coordinate vector of the k-th particle in the population be represented by X. k =(X k1 ,X k2 ,…,X k7 ), X k Each element represents the values ​​of a0, a1, a2, b0, b1, b2, and C0 of the particle in the optimization space; during the optimization process, the position update method of the k-th particle is shown in equation (18);

[0084] X k (T+1)=ω×X k (T)+c1×r×X P (T)+c2×r×X G (T) (18)

[0085] Where T represents the number of iterations completed so far, ω, c1, and c2 are the weight factor, individual learning rate, and group learning rate, respectively, and are set to 0.5, 1.5, and 1.5; r is a random number in the interval [0,1], and X... k (T) represents the individual position of particle k at the time of the T-th iteration update, X P (T), X G (T) represents the best individual position in the population and the global historical best individual position at the time of the Tth iteration update, respectively;

[0086] The position of each particle in the current population is iteratively updated according to formula (18). The values ​​of parameters a0, a1, a2, b0, b1, b2, and C0 corresponding to the particle are determined according to the particle position. Then, the objective function value of each updated particle is calculated. The particle with the smallest objective function value is taken as the optimal individual in the population when the (T+1)th iteration update is completed. At the same time, the position X of the particle is set. P The objective function value corresponding to (T+1) and the global historical best individual position X G The objective function values ​​corresponding to (T) are compared, and the position of the particle with the smallest objective function value is taken as the globally historical best individual position X when the (T+1)th iteration update is completed. G (T+1); then increment the iteration count by 1, repeat the process of iterating and updating each particle in the current population according to formula (18), and repeat the process of obtaining the best individual position in the population and the global historical best individual position when the current iteration update is completed, and repeat this process until the iteration update count reaches the maximum iteration count N. it This means completing one round of PSO algorithm computation.

[0087] When one round of the PSO algorithm operation is completed, after completing N... it In the next updated population, N was selected using a roulette wheel method. n Ten surviving particles directly become part of the initial population for the next round; this process is called extinction. The remaining N... n -N n The initial 10 individuals are randomly generated within the optimization interval and are called new individuals. However, during initialization, it must be ensured that the normalized Euclidean distance from all new individuals to any surviving individual is greater than 0.05 to guarantee the diversity of the initial population. This process is called cyclic initialization. After particle initialization is complete, the current round of PSO algorithm calculation continues. When the current round of PSO algorithm calculation is complete, the above population extinction and cyclic initialization process continues, and the next round of PSO algorithm calculation is performed, repeating continuously. The maximum number of population extinctions, N, is reached. exAt this point, the final loop initialization is performed, followed by the final round of PSO algorithm calculation. After the final round of PSO algorithm calculation is completed, the global historical best individual position of the entire population is the optimal solution, which yields the optimal values ​​of the 7 parameters to be identified.

[0088] Step 5: Substitute the optimal values ​​of the 7 parameters to be identified obtained in Step 4 into Formula (11) to obtain the optimal dynamic hysteresis model; obtain the magnetic flux density waveform of the dynamic alternating excitation condition to be applied to the target electrical steel single sample in Step 1, and then input the magnetic flux density waveform into the above-mentioned optimal dynamic hysteresis model to obtain the magnetic field strength waveform sequence of the target electrical steel single sample under the dynamic alternating excitation condition.

[0089] The magnetic flux density waveform of the proposed dynamic alternating excitation condition is a sinusoidal waveform or / and a waveform with harmonics or / and a bias, with a fundamental frequency in the range of 0.01Hz-3kHz.

[0090] The embodiments of the present invention will be illustrated below with examples.

[0091] This invention verifies the hysteresis loop interpolation extension method described in step 1. To verify the improvement of the accuracy of the static Preisach hysteresis model by loop extension, quasi-static tests were conducted on non-oriented silicon steel sheet 20CS1500. Seventeen hysteresis loops with magnetic flux density amplitudes ranging from 0.1T to 1.7T in 0.1T steps were measured. Only eight of these loops were selected for static Preisach hysteresis model identification, and the rest were used as a verification group. First, without extension, only the selected eight loops were used to complete the Everett function identification. Three loops were selected from the verification group for model accuracy verification, and the results are as follows. Figure 2 As shown in (1) in the example. Secondly, this example expands the selected 8 loop interpolations to 29 and uses them all for Everett function identification, where... In the interval The top section is divided into three equal parts. Three identical hysteresis loops are selected from the validation group for model accuracy validation, and the results are as follows: Figure 2 As shown in (2) of the diagram. Figure 2 It is evident that the interpolation extension of the hysteresis loop significantly improves the identification accuracy of the static model, and this improvement is particularly noticeable when experimental data is limited.

[0092] The parameter identification results of the dynamic Preisach hysteresis model described in steps 2, 3, and 4 are verified. To verify the accurate simulation of the actual hysteresis loop by the proposed dynamic Preisach hysteresis model, and to improve the efficiency and accuracy of the particle swarm optimization (EPSO) algorithm in identifying model parameters, firstly, the population size N of the EPSO algorithm is... n and the number of iterations Nit All are set to 30, and the number of extinctions N is also set to 30. ex For 4 times, ω, c1, and c2 are 0.5, 1.5, and 1.5 respectively. Next, the traditional dynamic model (10) is embedded into EPSO to identify the only undetermined parameter C0. The optimization interval is set to [0, 0.5], and the identification result is C0 = 0.2291. Then, the optimization intervals of a0, a1, a2, b0, b1, and b2 are uniformly set to [0.5, 1.5], and other algorithm settings remain unchanged. The parameters of the dynamic Preisach hysteresis model (11) described in this invention are identified, and the results are listed in Table 1. In this example, EPSO uses the objective function in formula (17) with 4 measured hysteresis loops, i.e., Nl = 4. The corresponding excitation conditions are listed in Table 2, and the number of discrete points Np for each loop is 300. To verify the improvement in modeling accuracy, the simulation results of the two models before and after the improvement are displayed in Table 2. Figure 3 The simulated excitation condition's magnetic flux density waveform is: a 400Hz sine wave superimposed with a 30% content and a 120° phase shift of the 7th harmonic, with a magnetic flux density amplitude of 1.3T. According to... Figure 3 It is evident that the dynamic model described in this invention has a significant advantage in hysteresis loop fitting accuracy. This is mainly because the modeling method described in step 2 considers the coupling effect between hysteresis nonlinearity and eddy current field and introduces relevant model parameters.

[0093] Table 1. Parameter identification results of the dynamic Preisach hysteresis model

[0094]

[0095] Table 2 shows the excitation conditions corresponding to the four measured hysteresis loops in the identification group.

[0096]

[0097] Furthermore, the performance improvement of EPSO over PSO needs to be verified separately. To ensure that the total number of iterations for both algorithms is consistent, in this example, the number of iterations N for PSO is... it Set to 150, other parameters are the same as EPSO, extinction count N ex No settings are required. Similarly, the measured hysteresis loops shown in Table 2 are used to identify the parameters of the dynamic Preisach hysteresis model (11). Figure 4Comparing the iterative curves of the objective function values ​​of the historical best individuals in the parameter optimization process of PSO and EPSO, it can be seen that the optimal solution of EPSO is already better than that of PSO in the early stage of convergence, and the best individual is kept updated before 110 generations, while traditional PSO converges after about 30 iterations. This shows that EPSO has better global search ability and convergence accuracy. The reason is that EPSO's population extinction mechanism preserves relatively excellent individuals for the algorithm, while the cyclic initialization process enhances the differences between particles. The combination of these two improvements improves the algorithm's ability to escape local optima while ensuring the quality of the globally optimal individuals. The optimization results of PSO are shown in Table 1, showing the difference between the results of the two algorithms. To make the verification results more convincing, the two sets of model parameters in Table 1 were used for dynamic hysteresis loop simulation. The simulation excitation conditions were: 800Hz SPWM wave excitation and bias magnetic flux density of 0.3T. The simulation results corresponding to the two sets of parameters were compared with the experimental results. Figure 5 Among them, the model parameters corresponding to (1) are identified by PSO, and the model parameters corresponding to (2) are identified by EPSO. It is easy to see that the model parameters obtained by EPSO achieve higher hysteresis loop prediction accuracy.

[0098] The simulation results obtained in step 5 were visualized and verified. Taking the visualization verification under harmonic excitation conditions as an example, based on the systematic magnetic property test of non-oriented silicon steel sheet 20CS1500, a total of 540 sets of experimental data under different harmonic configurations were stored. In order to achieve global visualization verification, each harmonic configuration was assigned to a pixel block and stacked in sequence to form a grid. Using the parameters obtained from EPSO in Table 1, the dynamic hysteresis loop under each set of excitation conditions was simulated using the improved dynamic Preisach hysteresis model. Based on the simulation and experimental results, the determination coefficient R was calculated by equation (19). 2 value.

[0099]

[0100] Where Np is the number of discrete points in a single cycle on a hysteresis loop. H mea (p) and H mod (p) represents the magnetic field strength at the p-th discrete time point in the measured and model-simulated hysteresis loops, respectively. R is the average magnetic field strength at all discrete points on the measured loop. Generally, R 2 The closer the value is to 1, the better the fit.

[0101] Set color separation threshold The value is 0.985, for all R values. 2 Value higher than The pixels are filled with white, and the rest are filled with black, generating... Figure 6 The binary Boolean image shown. Figure 6 The graph is a 15×36 binary Boolean image. The horizontal axis is divided into 15 equal pixel blocks, representing 15 magnetic flux density amplitudes from 0.2T to 1.6T in 0.1T increments. The vertical axis is divided into 36 equal pixel blocks, corresponding to any combination of 3 fundamental frequencies, 3 harmonic orders, 2 harmonic contents, and 2 harmonic phase shifts (a total of 36 combinations). The fundamental frequencies are 200Hz, 400Hz, and 800Hz; the harmonic orders are 3, 5, and 7; the harmonic contents are 30% and 50%; and the harmonic phase shifts are 0° and 120°. Figure 6 The order of these excitation conditions is marked. For example, the pixel block located in the 7th column from left to right and the 9th row from top to bottom represents the excitation conditions as follows: a 200Hz sine wave superimposed with a 30% content and a 120° phase shift of the 3rd harmonic, and a magnetic flux density amplitude of 0.8T.

[0102] There are 16 black blocks in the image. If these are considered as cases of poor model accuracy, they only account for 2.96% of all 540 sets of validation data. The model parameters were identified using only 4 sets of hysteresis data from Table 2. In this example, the redundancy of a large amount of experimental data is to verify the comprehensiveness of the model simulation accuracy. In practical applications, if there is also redundancy in the experimental data, a binary Boolean image can be generated, and the measured data under the excitation conditions corresponding to the black blocks can be selected to supplement the calculation of the EPSO objective function and the model parameters can be optimized again. In this example, based on the position of the black blocks, the measured data corresponding to the excitation conditions of fundamental frequency 200Hz, magnetic flux density amplitude 1.4T, and superimposed 5th harmonic (content 30%, phase shift 120°) are added to the parameter optimization process in step 4. That is, 5 sets of measured hysteresis are used to calculate the objective function in equation (17). After the second optimization, the number of black blocks in the model parameters obtained by EPSO is reduced to 14. However, the model's global accuracy will not increase indefinitely with the increase of identification data. The parameter optimization and model calculation results achieved by five sets of measured laps are already close to optimal. In engineering applications, whether to perform secondary optimization can be decided based on specific experimental conditions.

[0103] Using the method described above, the goodness-of-fit of the model under numerous excitation conditions can be integrated into a single image, by adjusting the color separation threshold. It can quickly evaluate the accuracy of the model itself and perform secondary optimization; by comparing the area of ​​color blocks in the binary Boolean images generated by different models, it can also intuitively judge the performance of different models.

[0104] This invention proposes a method for infinite interpolation expansion based on a small number of measured symmetrical hysteresis loops, reducing the amount of experimental data required for identifying static Preisach hysteresis models. The proposed improved dynamic modeling method introduces the magnetization rotation point and bias, exhibiting stronger adaptability under complex excitations compared to traditional models. The improved particle swarm optimization algorithm incorporates population extinction and cyclic initialization mechanisms, significantly improving optimization accuracy compared to the classic PSO algorithm. By embedding the dynamic Preisach hysteresis model into the improved PSO algorithm, more efficient model parameter identification and more accurate dynamic hysteresis loop simulation are achieved. Based on hysteresis loop test data of electrical steel under various excitation conditions, the overall performance of the proposed modeling method and parameter optimization algorithm is more comprehensively verified and enhanced. In practical applications, this method can provide accurate magnetization model support for high-precision magnetic field calculations of power electronic equipment cores.

[0105] Any aspects not covered in this invention are applicable to existing technologies.

Claims

1. A method for calculating the dynamic magnetic field strength of electrical steel under complex magnetic field excitation, characterized in that, The method specifically includes the following steps: Step 1: Under quasi-static symmetrical alternating excitation test conditions, collect no more than 10 quasi-static measured hysteresis loop data from a single sample of the target electrical steel. Then, using an infinite interpolation expansion method, expand the number of hysteresis loops based on the collected quasi-static measured hysteresis loops, obtaining a total of no less than 29 hysteresis loop data. Based on this hysteresis loop data, derive the functional expression of the Everett function and establish a static Preisach hysteresis model based on the Everett function. Step 2: Design a dynamic Preisach hysteresis model based on the static Preisach hysteresis model; The output value of the dynamic Preisach hysteresis model is calculated as follows: c = c L (B L (t))×γ G (B G (t)) (12) B G (t)=δB * (t) / B m (13) B L (t)=(B * (t)-B rev (t)) / (B dc +δ(|B m |-|B dc |)-B rev (t)) (14) Among them, B * (t) represents the input sequence of this dynamic hysteresis model, H hy (B * (t) represents the output value of the static Preisach hysteresis model established based on the method described in step 1; σ represents the conductivity of the target electrical steel monolithic sample, d represents the thickness of the target electrical steel monolithic sample, and C0 represents the parameter to be identified; B L (t) and B G (t) represents the normalized values ​​of local and global magnetic flux density, respectively, and serves as two intermediate quantities; B rev (t) represents the magnetic flux density waveform B. * (t) is the extreme value before time t and closest to time t; B dc For B * The magnetic flux density bias value of (t), B m For B * The amplitude of (t), δ=sign(dB) * (t) / dt), sign() is the sign function, which is used to determine the sign of a number. It returns 1 for positive numbers, 0 for zero, and -1 for negative numbers; Step 3: Construct the function γ in equation (12) using a quadratic polynomial. L () and γ G (), such that a0=γ L (0), a2=γ L (1) and b0 = γ G (-1),b1=γ G (0), b2=γ G (1); The function γ can be derived. L () and γ G The expression for () is: γ L (x)=(2a0-4a1+2a2)x 2 +(-3a0+4a1-a2)x+a0 (15) Among them, a0, a1, a2, b0, b1, and b2 are 6 model parameters to be identified, and together with C0 in equation (11), there are a total of 7 parameters to be identified. Step 4: Optimize the 7 parameters to be identified in Step 3 using an improved particle swarm optimization algorithm; The improved particle swarm optimization algorithm combines the traditional particle swarm optimization (PSO) algorithm with population extinction and cyclic initialization mechanisms. The specific process is as follows: Set the population size N. n The maximum number of iterations N in a single round of the PSO algorithm it And the maximum number of population extinctions N ex The optimization interval for C0 is set to [0, 0.5], and the optimization intervals for a0, a1, a2, b0, b1, and b2 are all set to [0.5, 1.5]. A random assignment method within the optimization interval is used to perform N-value calculations on the seven parameters to be identified (a0, a1, a2, b0, b1, b2, and C0). n The initialization yields N n The initial population of particles; based on the values ​​of a0, a1, a2, b0, b1, b2 of each particle, γ is obtained by combining equations (15) and (16). L (), γ G The function formula of (), then γ L (), γ G Substituting the function formula of () and the corresponding value of C0 into formula (11), and using Nl measured hysteresis loops under dynamic alternating excitation conditions as the identification group, the corresponding magnetic flux density sampling sequence B is extracted. * Substituting (t) into (11) as input, we obtain the magnetic field strength sequence H(t) under the corresponding excitation condition, and apply this sequence to the calculation of the objective function; Set the objective function as follows: Where Nl is the number of measured hysteresis loops in the identification group, Nl is no greater than 6, l is the hysteresis loop symbol in the identification group; Np is the number of discrete points in a single cycle on a hysteresis loop. and These represent the magnetic field strength at the p-th discrete time point in the l-th measured and model-simulated hysteresis loop, respectively. The input sequence B is calculated using formula (11). * (t) represents the magnetic flux density sequence corresponding to the l-th measured hysteresis loop in the identification group; |·| indicates the sign of the absolute value. According to formula (17), calculate the objective function Y for each particle in the initial population. Fit The position of the particle with the smallest objective function value is the optimal individual position and the global historical optimal individual position during the first iteration update, used to guide the first iteration update of the population. After the first iteration update is completed, the objective function value of each particle in the population that has completed the update is calculated, and the particle with the smallest objective function value is selected as the optimal individual in the population that has completed the first update. The objective function value of the optimal individual in the population that has completed the first update is compared with the objective function value of the optimal individual in the initial population, and the individual with the smallest objective function value is selected as the global historical optimal individual. After each iteration update is completed, the optimal individual in the current population is obtained by calculating the objective function value, and this optimal individual is compared with the global historical optimal individual in the previous iteration update process, and the individual with the smallest objective function value is selected as the current global historical optimal individual. The position of the current optimal individual in the population and the position of the current global historical optimal individual are used to guide the next iteration update of the population. The specific process is as follows: Let the position coordinate vector of the k-th particle in the population be represented by X. k =(X k1 ,X k2 ,…,X k7 ), X k Each element represents the values ​​of a0, a1, a2, b0, b1, b2, and C0 of the particle in the optimization space; during the optimization process, the position update method of the k-th particle is shown in equation (18); X k (T+1)=ω×X k (T)+c1×r×X P (T)+c2×r×X G (T) (18) Where T represents the number of iterations completed so far, ω, c1, and c2 are the weight factor, individual learning rate, and group learning rate, respectively, and are set to 0.5, 1.5, and 1.5; r is a random number in the interval [0,1], and X... k (T) represents the individual position of particle k at the time of the T-th iteration update, X P (T), X G (T) represents the best individual position in the population and the global historical best individual position at the time of the Tth iteration update, respectively; The position of each particle in the current population is iteratively updated according to formula (18). The values ​​of parameters a0, a1, a2, b0, b1, b2, and C0 corresponding to the particle are determined according to the particle position. Then, the objective function value of each updated particle is calculated. The particle with the smallest objective function value is taken as the optimal individual in the population when the (T+1)th iteration update is completed. At the same time, the position X of the particle is set. P The objective function value corresponding to (T+1) and the global historical best individual position X G The objective function values ​​corresponding to (T) are compared, and the position of the particle with the smallest objective function value is taken as the globally historical best individual position X when the (T+1)th iteration update is completed. G (T+1); then increment the iteration count by 1, repeat the process of iterating and updating each particle in the current population according to formula (18), and repeat the process of obtaining the best individual position in the population and the global historical best individual position when the current iteration update is completed, and repeat this process until the iteration update count reaches the maximum iteration count N. it This means completing one round of PSO algorithm computation. When one round of the PSO algorithm operation is completed, after completing N... it In the next updated population, N was selected using a roulette wheel method. n Ten surviving particles directly become part of the initial population for the next round; this process is called population extinction. The remaining N... n -N n The initial 10 individuals are randomly generated within the optimization interval and are called new individuals. During initialization, it is ensured that the normalized Euclidean distance from all new individuals to any surviving individual is greater than 0.05 to guarantee the diversity of the initial population. This process is called cyclic initialization. After the particle initialization is completed, the current round of PSO algorithm calculation continues. When the current round of PSO algorithm calculation is completed, the above population extinction and cyclic initialization process continues, and the next round of PSO algorithm calculation is performed, repeating continuously. This process continues until the maximum number of population extinctions, N, is reached. ex At that time, the final loop initialization is performed, and the final round of PSO algorithm operation is performed. After the final round of PSO algorithm operation is completed, the global historical best individual position of the entire population is the optimal solution, that is, the optimal values ​​of the 7 parameters to be identified are obtained. Step 5: Substitute the optimal values ​​of the 7 parameters to be identified obtained in Step 4 into Formula (11) to obtain the optimal dynamic hysteresis model; obtain the magnetic flux density waveform of the dynamic alternating excitation condition to be applied to the target electrical steel single sample in Step 1, and then input the magnetic flux density waveform into the above-mentioned optimal dynamic hysteresis model to obtain the magnetic field strength waveform sequence of the target electrical steel single sample under the dynamic alternating excitation condition.

2. The method for calculating the dynamic magnetic field strength of electrical steel under complex magnetic field excitation as described in claim 1, characterized in that, In step 1, 8 to 10 quasi-static measured hysteresis loop data are collected from a single sample of the target electrical steel.

3. The method for calculating the dynamic magnetic field strength of electrical steel under complex magnetic field excitation as described in claim 1, characterized in that, In step 1, the interpolation augmentation of the hysteresis loop data is achieved through the basic magnetization curve and weighting coefficients: Wherein, the superscript i represents the sequence number of the measured hysteresis loop data, and f() is the basic magnetization curve function; and These represent the magnetic flux density amplitudes of the i-th and (i+1)-th measured hysteresis loops, respectively. and These represent the magnetic field strength amplitudes of the i-th and (i+1)-th measured hysteresis loops, respectively. and These represent the amplitudes of the magnetic flux density and magnetic field strength of the expanded hysteresis loop, respectively. f() is obtained by solving for all the vertices of the measured hysteresis loops using piecewise cubic Hermite interpolation, and is essentially a single-valued mapping from B to H. and These represent the weighting coefficients of the magnetic flux density (B) waveform and the magnetic field strength (H) waveform of the extended hysteresis loop, respectively. Where t represents a discrete sampling time point; and Let B(t) be the B-waveform sequence and H-waveform sequence of the extended hysteresis loop. i B(t) i+1 Let H(t) represent the B-wavelength sequences of the i-th and (i+1)-th measured hysteresis loops, respectively. i H(t) i+1 These represent the H-waveform sequences of the i-th and (i+1)-th measured hysteresis loops, respectively.

4. The method for calculating the dynamic magnetic field strength of electrical steel under complex magnetic field excitation as described in claim 1, characterized in that, In step 1, based on the measured and extended hysteresis loops, the Everett function of the static Preisach hysteresis model is constructed as follows: Wherein, the bivariate function E(α,β) represents the Everett function; H α () represents the measured or extended hysteresis loop descending branch function with magnetic flux density amplitude α, and is a single-valued mapping from B to H; α can take the magnetic flux density amplitude of all measured and extended hysteresis loops; for a specific α, β can take any value satisfying equation (6): -α≤β≤α (6) After obtaining the scatter set (α,β,E(α,β)), the Everett function is used to obtain a supplementary scatter set (-β,-α,E(α,β)) by utilizing the symmetry of the Everett function along the line α=-β; linear interpolation is then performed using the two scatter sets to obtain the Everett function expression E(α,β) with α and β as variables.

5. The method for calculating the dynamic magnetic field strength of electrical steel under complex magnetic field excitation as described in claim 1, characterized in that, In step 1, for a given magnetic flux density input sequence B(t), the output value of the static Preisach hysteresis model based on the Everett function is calculated by (7): In the formula, B(t) is the magnetic flux density input sequence of the static Preisach hysteresis model, and H hy (B(t)) represents the output value of the model at time t; N REV Let E1 be the total number of inversion points of the Everett function at time t, and let E1 = E(α1, β1); H hy (B(t-Δt)) is the output value of the model at time t-Δt, where Δt is the discrete interval of the time sampling points, determined by the experimental testing system of the hysteresis loop, and is a constant; E ∑+ and E ∑- Two intermediate quantities for calculation are defined by (8) and (9) respectively: In the formula, (α) n ,β n ) represents the coordinates of the nth inversion point of the Everett function.

6. The method for calculating the dynamic magnetic field strength of electrical steel under complex magnetic field excitation as described in claim 5, characterized in that, At time t, α n ,β n The method for determining it is as follows: Assume that at time t and before, the B(t) sequence has a total of N occurrences. br There are several extreme values. Arrange these extreme values ​​in chronological order to form an extreme value sequence. Add an element B(0) before the first element of the extreme value sequence. That is, add the first value in the B(t) sequence. If the value of this element is less than zero, it is recorded as a minimum value. Otherwise, it is recorded as a maximum value. Check the Kth sequence br Is element K in the Kth order? br -1 and the Kth br Within the interval defined by -2 element values, if not, delete the Kth element in the extreme value sequence. br -1 and the Kth br -2 elements are used to form a new extreme value sequence; then, the Kth element in the new sequence is... br Repeat the same operation on each element until the sequence satisfies the Kth element. br The element is always in the Kth position. br -1 and the Kth br The sequence continues until the value of the first element is within the range defined by -2 element values, or until the number of elements in the sequence is less than 3; then, it checks whether the absolute value of the first element of the sequence is the maximum absolute value among all elements in the sequence. If not, the element is deleted; during this process, K br It can be any value greater than 2 and not greater than the total length of the latest extreme value sequence; This process is called erasure; After erasure, the number of elements in the extreme value sequence is N. REV For the first element, if it is a maximum value, then denote the value of this element as α1 and the negative of this element as β1. If it is a local minimum, then denote the value of this element as β1, and denote the negative of this element as α1; for the nth element, 1 <n≤N REV If it is a maximum value, then the value of this element is denoted as α. n Let the value of the (n-1)th element be β. n If it is a local minimum, then the value of this element is β. n Let the value of the (n-1)th element be α. n Thus, a total of N can be obtained. REV Everett function inverts the coordinates of the point (α) n ,β n ).

7. The method for calculating the dynamic magnetic field strength of electrical steel under complex magnetic field excitation as described in claim 1, characterized in that, In step 4, Nl is 4.

8. The method for calculating the dynamic magnetic field strength of electrical steel under complex magnetic field excitation as described in claim 7, characterized in that, The excitation conditions for the four measured hysteresis loops of the identification group are shown in Table 2. Table 2 shows the excitation conditions corresponding to the four measured hysteresis loops in the identification group.

9. The method for calculating the dynamic magnetic field strength of electrical steel under complex magnetic field excitation as described in claim 1, characterized in that, In step 5, the magnetic flux density waveform of the dynamic alternating excitation condition to be applied is a sinusoidal waveform or / and a waveform with harmonics or / and a bias, with a fundamental frequency in the range of 0.01Hz-3kHz.

Images