CN112069668A - Electromagnetic transient rapid simulation method based on differential quadrature method - Google Patents

Electromagnetic transient rapid simulation method based on differential quadrature method Download PDF

Info

Publication number
CN112069668A
CN112069668A CN202010871833.0A CN202010871833A CN112069668A CN 112069668 A CN112069668 A CN 112069668A CN 202010871833 A CN202010871833 A CN 202010871833A CN 112069668 A CN112069668 A CN 112069668A
Authority
CN
China
Prior art keywords
matrix
formula
differential
value
electromagnetic transient
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202010871833.0A
Other languages
Chinese (zh)
Other versions
CN112069668B (en
Inventor
张静
叶婧
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Three Gorges University CTGU
Original Assignee
China Three Gorges University CTGU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China Three Gorges University CTGU filed Critical China Three Gorges University CTGU
Priority to CN202010871833.0A priority Critical patent/CN112069668B/en
Publication of CN112069668A publication Critical patent/CN112069668A/en
Application granted granted Critical
Publication of CN112069668B publication Critical patent/CN112069668B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/12—Simultaneous equations, e.g. systems of linear equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13—Differential equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/04—Constraint-based CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
    • Y02E60/00—Enabling technologies; Technologies with a potential or indirect contribution to GHG emissions mitigation

Abstract

The electromagnetic transient rapid simulation method based on the differential quadrature method applies the differential quadrature method to electromagnetic transient simulation, and provides an electromagnetic transient rapid calculation method based on the differential quadrature method according to the V transformation characteristic of a weighting coefficient matrix in the differential quadrature method. The invention avoids direct inversion calculation of the matrix which is subjected to dimension increase due to the introduction of inner points, blocks the matrix according to the characteristics of the dimension increase matrix, and converts the inversion calculation of a large-scale matrix into simple prior generation operation and the inversion calculation of a small-scale matrix, thereby improving the electromagnetic transient calculation efficiency based on a differential quadrature method. The calculation example shows that the algorithm can obtain the effective acceleration ratio for the nonlinear power system.

Description

Electromagnetic transient rapid simulation method based on differential quadrature method
Technical Field
The invention relates to the field of electromagnetic transient simulation calculation, in particular to an electromagnetic transient rapid simulation method based on a differential quadrature method.
Background
With the development of modern industrial technology and the introduction of distributed energy, the power grid structure becomes huge and complex, and the power grid can be subjected to the impact of operation overvoltage and lightning overvoltage at any time and can also be subjected to the influence caused by various types of faults. In order to protect expensive electrical equipment in a power grid from being damaged under a steady state or fault condition, the research on electromagnetic transient simulation calculation is particularly important. In electromagnetic transient numerical calculation, the calculation efficiency is low due to the introduction of a differential equation representing details; the precision and the speed cannot be considered comprehensively. In recent years, people have developed a plurality of researches on four aspects of modeling technologies of a parallel computing architecture, a novel computing device simulation platform, a hybrid computing architecture simulation mode and a fast solving mode and have achieved fruitful results. It can be seen that the large-step computation and the parallel computation play an important role in improving the electromagnetic transient computation efficiency.
Many documents develop researches on the problem of large-step simulation, wherein large-step electromagnetic transient simulation based on time scale change is proposed, namely, a time domain signal is converted into a complex signal through Hilbert change, and the complex signal is subjected to rotation conversion to obtain a slowly-changing signal, so that large-step simulation is realized. On the other hand, the introduction of dynamic phasors brings more possibilities for electromagnetic transient large step simulation. The essence of large-step simulation realized by the dynamic phasor modeling method is that the down-conversion treatment is performed, so that the fast change is changed into the slow change, but the dynamic phasor also has a lot of defects, the modeling of the dynamic phasor on the nonlinear element is still very complex, and when various non-negligible frequency components appear in the nonlinear system, huge calculation amount is increased. In general, studies on dynamic phasors have yet to be increased.
Meanwhile, parallel computation is also active in electromagnetic transient numerical computation as a means for improving simulation speed. And discretizing the differential equation by using an implicit trapezoidal method and an implicit Euler method, wherein the implicit Euler method is used for discretizing at a sudden change point to avoid virtual numerical value oscillation, and the implicit trapezoidal method is used for discretizing at other moments. Based on the point of taking different step lengths, the discrete algebraic equation successfully realizes matrix diagonalization, and all time points are decoupled with each other to obtain the observable acceleration ratio.
In recent years, more and more numerical methods are used in electromagnetic transient simulation calculations. The classical implicit trapezoidal method is widely applied to electromagnetic transient calculation due to the A stability and the second-order calculation precision. However, the implicit trapezoidal method does not have L stability, i.e., the oscillation characteristic under the condition of sudden change cannot be effectively suppressed, and the simulation step size is small, so that more multi-level and multi-order methods are applied to electromagnetic transient calculation. The differential quadrature method is simple in principle and easy to implement. Has A stability and can effectively inhibit oscillation. The calculation result is s-level s-order precision, and compared with an implicit trapezoidal method, a larger simulation step length can be adopted. However, the introduction of interior points follows the s-level differential integration method, and since the number of the interior points needs to be introduced in each step, the dimension of the coefficient matrix is increased to s times of the original dimension, so that a great deal of time is spent on matrix inversion, and the significance of large-step calculation is lost.
Disclosure of Invention
The method aims at the problem of dimension disaster caused by introduction of interior points in electromagnetic transient simulation calculation by an s-level differential quadrature method. The invention provides an electromagnetic transient rapid simulation method based on a differential quadrature method, which is characterized in that on the basis of the differential quadrature method, the characteristic of V transformation is combined, block solution is carried out according to the characteristic of a matrix after dimension increase, and the speed of matrix inversion is improved; the advantage of large-step-length simulation adopted by the differential quadrature method is fully exerted, and the electromagnetic transient calculation speed is improved.
The technical scheme adopted by the invention is as follows:
the electromagnetic transient rapid simulation method based on the differential quadrature method comprises the following steps:
step 1: the electromagnetic transient process is a transformation process of voltage and current caused by the change of a magnetic field and an electric field, and is represented by a set of differential equations with initial values as follows:
Figure BDA0002651329910000021
wherein: x is formed by Rm×1X is the set of state variables in the system, g (t) is related to the time variable, x0Then is t0Time (t ═ 0)Initial value of state variable x.
The differential quadrature method is as follows:
if the function f (x) is smoothly derivable over an interval, f (x) is at grid point ciThe derivative over i e (1, s) can be represented by a linear weighted sum of the function values at all grid points over the interval, i.e.:
Figure BDA0002651329910000022
order:
gij∈G (3)
wherein: f. of(1)(ci) I e (1, s) represents the derivative value of the function f (x) at any grid point over the interval, gijWeighting coefficients for differential integration, f (c)j) J e (1, s) represents the function value of the function f (x) at any grid point in the interval, s being the number of grid points in the interval.
The mathematical expression of the differential integration method is derived as follows:
Figure BDA0002651329910000031
in the formula, tnIs a starting time, tn+1For the end time, x is the set of state variables in the system, and x ∈ Rm×1，xnIs tnValue of the time-of-day state variable x, xn+1Is tn+1The value of the state variable to be solved at the moment,
Figure BDA0002651329910000032
indicates an arbitrary grid point c on the intervaliThe value of the state variable, h is the simulation step length, aijFor the weighting coefficients, A ≡ a is definedij]And A ═ G-1And s is the number of interior points.
Discretization (1) by differential integration and written in matrix form:
Figure BDA0002651329910000033
step 2: using A ═ G-1Equally transforming equation (5) in step 1 may result:
Figure BDA0002651329910000034
in formula (6), G is a weighting coefficient matrix of differential quadrature method, ImIs m-dimensional unit matrix, e is s-dimensional unit column vector, xnIs tnThe value of the state variable x at the moment, h is the simulation step length.
In the formula (6), the mathematical expression is explained as follows:
Figure BDA0002651329910000035
(1) when formula (1) is described as a non-linear problem, the definition:
Figure BDA0002651329910000036
in formula (7): s is the number of the inner points,
Figure BDA0002651329910000037
indicates the grid point c on the intervalkValue of state variable
Solving a matrix equation (6) after dispersion by using a Newton iteration method, and taking an average value vector of s inner points as an initial value of the Jacobian matrix to obtain:
Figure BDA0002651329910000041
in formula (8):
Figure BDA0002651329910000042
in the formula: g is differential quadrature weightingCoefficient matrix, ImIs an m-dimensional identity matrix, h is a simulation step length, IsIs an s-dimensional unit matrix, xi is the iteration number,
Figure BDA0002651329910000043
denotes the ξ -th average of the interior points in the iteration,
Figure BDA0002651329910000044
denotes the ξ +1 th iteration at grid point ckThe value of the upper state variable x.
In formula (9):
Figure BDA0002651329910000045
and step 3: the V transform form of the differential quadrature method is:
Figure BDA0002651329910000046
in the formula (10), G is a differential multiplication coefficient matrix, V is a Van der Mond matrix, and V is-1Being the inverse of the van der Mond matrix
Order:
Figure BDA0002651329910000047
in the formula (11), V-1Is the inverse of the van der Mond matrix, ImIs an m-dimensional identity matrix, and
Figure BDA0002651329910000048
applying the formulas (9) and (10) to the formula (8) in the step 2 to obtain:
Figure BDA0002651329910000049
in the formula:
Figure BDA00026513299100000410
in the formula (13), ImIs an m-dimensional identity matrix, IsIs an s-dimensional unit matrix and is,
Figure BDA0002651329910000051
xi is the number of iterations, step 4:
Figure BDA0002651329910000052
unfolding and blocking formula (12) obtained in step 3 to obtain:
Figure BDA0002651329910000053
Figure BDA0002651329910000054
Figure BDA0002651329910000055
in the formula (15), ImIs an m-dimensional unit matrix and is a matrix,
Figure BDA0002651329910000056
ξ is the number of iterations.
Equation (12) can be finally written as:
Figure BDA0002651329910000057
observed to obtain, H2Is a lower triangular matrix, so that y can be transformed by simple prior generation operation2Is expressed as y1Without the need to perform complex trigonometric decomposition, i.e.:
y2＝φ(y1)； (19)
bringing formula (19) into H3y1+H4y2＝r2Finally, the formula (18) is converted to contain y alone1The algebraic system of equations, i.e.:
Figure BDA0002651329910000058
Q∈Rm×m (21)
q is the calculated coefficient matrix.
Using LU decomposition to solve y1Through y2＝φ(y1) To give out y2Finally, Δ x is solved back. And calculating the next step length by judging whether the current step length is smaller than the convergence precision or not, if the current step length is smaller than the convergence precision, calculating the next step length, and if the current step length is larger than the convergence precision, performing the next iteration by taking the average value of the interior points as an initial value until the current step length is converged into the precision range.
And 5: when the problem described by equation (1) is a linear time-varying problem, that is:
f(x)＝Utx (22)
Utfor a matrix of coefficients relating only to time variables, we define:
Figure BDA0002651329910000061
in the formula (23), s is the number of inner points, U (t)n+cih) I ∈ (1, s) is tn+ciValue of the time-h coefficient matrix, ciIs the grid point, h is the simulation step length.
Order:
Figure BDA0002651329910000062
in the formula (24), V-1Is the inverse of the van der Mond matrix, ImIs an m-dimensional unit matrix and is a matrix,
Figure BDA0002651329910000063
the linear time-varying differential equation is dispersed by using a differential quadrature method and is introduced into U to obtain:
Figure BDA0002651329910000064
wherein:
Figure BDA0002651329910000065
Figure BDA0002651329910000066
Figure BDA0002651329910000067
in the formula (27), V-1Is the inverse of the van der Mond matrix, ImIs an m-dimensional identity matrix.
In the formula (28), G is a weighting coefficient matrix of differential integration method, ImIs m-dimensional unit matrix, e is s-dimensional unit column vector, xnIs tnThe value of the state variable x at time, h is the simulation step length, g (t) ═ gT(tn+c1h) … gT(tn+csh)]T。
Matrix partitioning processing is performed on equation (25) to obtain:
Figure BDA0002651329910000068
Figure BDA0002651329910000071
Figure BDA0002651329910000072
equation (24) is finally written as:
Figure BDA0002651329910000073
H2is a lower triangular matrix, so that y can be transformed by simple prior generation operation2Is expressed as y1The linear expression of (a), namely:
y2＝φ(y1)； (33)
bringing formula (33) into H3y1+H4y2＝r2Finally, the formula (32) is converted to contain y alone1The algebraic system of equations, i.e.:
Figure BDA0002651329910000074
Q∈Rm×m (35)
using LU decomposition to solve y1Through y2＝φ(y1) To give out y2Finally, the solution is reversely solved by the formula (24)
Figure BDA0002651329910000075
The calculation of the next step is started.
The invention applies a differential quadrature method to electromagnetic transient simulation, and provides an electromagnetic transient fast calculation method based on the differential quadrature method according to the V transformation characteristic of a weighting coefficient matrix in the differential quadrature method. According to the invention, direct inverse calculation of the matrix which is subjected to dimension increase due to the introduction of the inner points is avoided, the matrix is partitioned according to the characteristics of the dimension increase matrix, and the inverse calculation of the large-scale matrix is converted into simple prior generation operation and the inverse calculation of the small-scale matrix, so that the electromagnetic transient calculation efficiency based on a differential quadrature method is improved. The calculation example shows that the algorithm can obtain the effective acceleration ratio for the nonlinear power system.
Drawings
Fig. 1 is a diagram of a high-voltage transmission line system with nonlinear loads.
FIG. 2 is a pi model equivalent circuit diagram.
Fig. 3 is a diagram of the voltage at the head end and the tail end of the transmission line.
Fig. 4 is a precision error map.
Detailed Description
The electromagnetic transient rapid simulation method based on the differential quadrature method comprises the following steps:
step 1: the electromagnetic transient process is a transformation process of voltage and current caused by the change of a magnetic field and an electric field, and is represented by a set of differential equations with initial values as follows:
Figure BDA0002651329910000081
wherein: x is formed by Rm×1X is the set of state variables in the system, g (t) is related to the time variable, x0Then is t0(t ═ 0) the initial value of the state variable x at time.
The differential quadrature method is as follows:
if the function f (x) is smoothly derivable over an interval, f (x) is at grid point ciThe derivative over i e (1, s) can be represented by a linear weighted sum of the function values at all grid points over the interval, i.e.:
Figure BDA0002651329910000082
order:
gij∈G (3)
wherein: f. of(1)(ci) I e (1, s) represents the derivative value of the function f (x) at any grid point over the interval, gijWeighting coefficients for differential integration, f (c)j) J e (1, s) represents the function value of the function f (x) at any grid point in the interval, s being the number of grid points in the interval.
The mathematical expression of the differential integration method is derived as follows:
Figure BDA0002651329910000083
in the formula, tnIs a starting time, tn+1For the end time, x is the set of state variables in the system, and x ∈ Rm×1，xnIs tnValue of the time-of-day state variable x, xn+1Is tn+1The value of the state variable to be solved at the moment,
Figure BDA0002651329910000084
indicates an arbitrary grid point c on the intervaliThe value of the state variable, h is the simulation step length, aijFor the weighting coefficients, A ≡ a is definedij]And A ═ G-1And s is the number of interior points.
Discretization (1) by differential integration and written in matrix form:
Figure BDA0002651329910000085
step 2:
using A ═ G-1Equally transforming equation (5) in step 1 may result:
Figure BDA0002651329910000091
in formula (6), G is a weighting coefficient matrix of differential quadrature method, ImIs m-dimensional unit matrix, e is s-dimensional unit column vector, xnIs tnThe value of the state variable x at the moment, h is the simulation step length.
In the formula (6), the mathematical expression is explained as follows:
Figure BDA0002651329910000092
wherein:
Figure BDA0002651329910000093
indicates an arbitrary grid point c on the intervaliForm of pointThe value of the state variable is set to,
Figure BDA0002651329910000094
representing the function f (x) at any grid point in the interval, gT(tn+cjh) J e (1, s) represents the function value at the grid point that is related to only the time variable.
(1) When formula (1) is described as a non-linear problem, the definition:
Figure BDA0002651329910000095
in formula (7): s is the number of the inner points,
Figure BDA0002651329910000096
indicates the grid point c on the intervalkValue of state variable
Solving a matrix equation (6) after dispersion by using a Newton iteration method, and taking an average value vector of s inner points as an initial value of the Jacobian matrix to obtain:
Figure BDA0002651329910000097
in formula (8):
Figure BDA0002651329910000098
in the formula: g is a weighting coefficient matrix of the differential quadrature method, ImIs an m-dimensional identity matrix, h is a simulation step length, IsIs an s-dimensional unit matrix, xi is the iteration number,
Figure BDA0002651329910000101
denotes the ξ -th average of the interior points in the iteration,
Figure BDA0002651329910000102
denotes the ξ +1 th iteration at grid point ckThe value of the upper state variable x.
In formula (9):
Figure BDA0002651329910000103
wherein:
Figure BDA0002651329910000104
indicates an arbitrary grid point c on the intervaliThe value of the state variable is set to,
Figure BDA0002651329910000105
representing the function f (x) at any grid point in the interval, gT(tn+cjh) J e (1, s) represents the function value at the grid point that is related to only the time variable.
And step 3:
the V transform form of the differential quadrature method is:
Figure BDA0002651329910000106
in the formula (10), G is a differential multiplication coefficient matrix, V is a Van der Mond matrix, and V is-1Being the inverse of the van der Mond matrix
Order:
Figure BDA0002651329910000107
in the formula (11), V-1Is the inverse of the van der Mond matrix, ImIs an m-dimensional identity matrix, and
Figure BDA0002651329910000108
applying the formulas (9) and (10) to the formula (8) in the step 2 to obtain:
Figure BDA0002651329910000109
in the formula:
Figure BDA00026513299100001010
in the formula (13), ImIs an m-dimensional identity matrix, IsIs an s-dimensional unit matrix and is,
Figure BDA00026513299100001011
xi is the number of iterations,
and 4, step 4:
Figure BDA0002651329910000111
in equation (14), a represents an element in the matrix, which is related to the selected mesh type.
Unfolding and blocking formula (12) obtained in step 3 to obtain:
Figure BDA0002651329910000112
Figure BDA0002651329910000113
Figure BDA0002651329910000114
in the formula (15), ImIs an m-dimensional unit matrix and is a matrix,
Figure BDA0002651329910000115
ξ is the number of iterations.
Equation (12) can be finally written as:
Figure BDA0002651329910000116
observed to obtain, H2Is a lower triangular matrix, so that y can be transformed by simple prior generation operation2Is expressed as y1Without the need to perform complex trigonometric decomposition, i.e.:
y2＝φ(y1)； (19)
bringing formula (19) into H3y1+H4y2＝r2Finally, the formula (18) is converted to contain y alone1The algebraic system of equations, i.e.:
Figure BDA0002651329910000117
Q∈Rm×m (21)
q is the calculated coefficient matrix.
Using LU decomposition to solve y1Through y2＝φ(y1) To give out y2Finally, Δ x is solved back. And calculating the next step length by judging whether the current step length is smaller than the convergence precision or not, if the current step length is smaller than the convergence precision, calculating the next step length, and if the current step length is larger than the convergence precision, performing the next iteration by taking the average value of the interior points as an initial value until the current step length is converged into the precision range.
And 5:
when the problem described by equation (1) is a linear time-varying problem, that is:
f(x)＝Utx (22)
Utfor a matrix of coefficients relating only to time variables, we define:
Figure BDA0002651329910000121
in the formula (23), s is the number of inner points, U (t)n+cih) I ∈ (1, s) is tn+ciValue of the time-h coefficient matrix, ciIs the grid point, h is the simulation step length.
Order:
Figure BDA0002651329910000122
in the formula (24), the reaction mixture is,V-1is the inverse of the van der Mond matrix, ImIs an m-dimensional unit matrix and is a matrix,
Figure BDA0002651329910000123
the linear time-varying differential equation is dispersed by using a differential quadrature method and is introduced into U to obtain:
Figure BDA0002651329910000124
wherein:
Figure BDA0002651329910000125
Figure BDA0002651329910000126
Figure BDA0002651329910000127
in the formula (27), V-1Is the inverse of the van der Mond matrix, ImIs an m-dimensional identity matrix.
In the formula (28), G is a weighting coefficient matrix of differential integration method, ImIs m-dimensional unit matrix, e is s-dimensional unit column vector, xnIs tnThe value of the state variable x at time, h is the simulation step length, g (t) ═ gT(tn+c1h) … gT(tn+csh)]T。
Matrix partitioning processing is performed on equation (25) to obtain:
Figure BDA0002651329910000131
Figure BDA0002651329910000132
Figure BDA0002651329910000133
equation (24) is finally written as:
Figure BDA0002651329910000134
H2is a lower triangular matrix, so that y can be transformed by simple prior generation operation2Is expressed as y1The linear expression of (a), namely:
y2＝φ(y1)； (33)
bringing formula (33) into H3y1+H4y2＝r2Finally, the formula (32) is converted to contain y alone1The algebraic system of equations, i.e.:
Figure BDA0002651329910000135
Q∈Rm×m (35)
using LU decomposition to solve y1Through y2＝φ(y1) To give out y2Finally, the solution is reversely solved by the formula (24)
Figure BDA0002651329910000139
The calculation of the next step is started.
In step 1, the G matrix expression in formula (2) is:
Figure BDA0002651329910000136
wherein:
Figure BDA0002651329910000137
Figure BDA0002651329910000138
wherein V is a Van der Monde matrix, s is the number of interval grid points, ckAnd k ∈ (1, s) is an interval grid point.
For the selection of the grid points, there are Legendre grids, Chebyshev grids, uniform grids and the like, and the function variables in the formula (5) are as follows:
Figure BDA0002651329910000141
in the step 2, substituting a ═ G-1Equal transformation is achieved. This has the advantage of converting formula (4) to formula (5) in a simpler form.
In step 3, an iterative format in which the average phasor of s interior points is used as an initial value is adopted, so that the following formula holds:
Figure BDA0002651329910000142
in the formula (39), V is a van der Mond matrix, IsIs an s-dimensional identity matrix, ImIs an m-dimensional identity matrix, and
Figure BDA0002651329910000143
ξ is the number of iterations.
And finally in the form of formula (12) in step 3.
In the step 3, the function is to utilize the characteristic of V transformation of differential quadrature method to realize the special structure
Figure BDA0002651329910000149
And the separation is carried out, so that the matrix of inversion is more sparse, and the calculation amount is reduced. And aim at
Figure BDA00026513299100001410
And (5) a matrix characteristic, namely partitioning the matrix.
In the step 4, the process of the step,
Figure BDA0002651329910000144
Figure BDA0002651329910000145
H3＝[a1,sIm]； (42)
Figure BDA0002651329910000146
Figure BDA0002651329910000147
Figure BDA0002651329910000148
Figure BDA0002651329910000151
Figure BDA0002651329910000152
in the formula ImIs an m-dimensional unit matrix and is a matrix,
Figure BDA0002651329910000153
ξ is the number of iterations.
In the step 4, the function is to perform special block processing on the formula (12) according to the characteristics of the formula (12), so that the triangular decomposition of the s × m dimensional matrix needs to be converted into the triangular decomposition of the m dimensional matrix in the serial process, the calculation speed of matrix inversion is increased, and the electromagnetic transient calculation speed based on the differential integration method is further increased.
In the step 5:
Figure BDA0002651329910000154
Figure BDA0002651329910000155
H3＝[a1,sIm]； (50)
H4＝[0 … -U] (51)
Figure BDA0002651329910000156
Figure BDA0002651329910000157
Figure BDA0002651329910000158
Figure BDA0002651329910000159
in the formula ImIs an m-dimensional unit matrix and is a matrix,
Figure BDA00026513299100001510
taking the simulation of the high-voltage transmission line containing the nonlinear load as an example, the calculation method is verified, and the method comprises the following steps:
fig. 1 is a diagram of a high-voltage transmission line system with nonlinear loads. Head end of power transmission line L, power supply and resistor RsInductance LsSerially connected, the tail end of the power transmission line is connected with a non-linear load consisting of a fixed resistor RLAnd a non-linear inductance LLIs represented by a combination of LLMiddle magnetic linkage phi and iLThe mathematical relationship of (1) is as follows: phi is atanhbiL. The whole length of the transmission line is 100 km. The resistance, inductance and capacitance to ground of the circuit unit length are respectively R0，L0，C0. The above parameters are given in the form of table 1.
TABLE 1 parameter table
Figure BDA0002651329910000161
FIG. 2 is a pi-type equivalent circuit diagram. The electromagnetic transient process of each section of the transmission line is represented by the following formula:
Figure BDA0002651329910000162
in the formula (56), L0，R0And C0Respectively an inductor, a resistor and a capacitor of a unit length of the line. u. ofk(t)、uk+1(t) is the voltage of node k and the voltage of node k +1 in the equivalent model of the transmission line shown in FIG. 2, ik-1(t)、ik(t) the current injected into node k and the current flowing out of node k to the transmission line, respectively, for the transmission line.
The total simulation time is 0.004 seconds, and the switch-on is carried out at the moment when t is 0. The constraint equation of the head end and the tail end of the transmission line is as follows:
Figure BDA0002651329910000163
in the formula (57), LSAnd RsFor connecting inductances and resistances in series to the transmission line, i0Current injected into the transmission line for the power supply, v1The voltage of the head end of the power transmission line L is shown as e (t), and the voltage of the power supply is shown as e (t); v. ofN+1For the voltage of the L-end node of the transmission line, iN+1For current flowing from the end node to the load, RLIs the equivalent resistance of the load, iLIs the current flowing through the non-linear part of the load; l isLThe inductance of the nonlinear part of the load is a constant coefficient of a and b, and specific values are given in table 1.
The transmission line equation and the constraint equations at the head end and the tail end are combined, and finally, the transmission line electromagnetic transient mathematical model containing the nonlinear load is shown as the following formula:
Figure BDA0002651329910000171
x is a system state variable, UxIs a coefficient matrix containing state variables. g (t) is a variable in the system that is only time dependent.
And (3) respectively dispersing the above formula by using a classical implicit trapezoidal method and a serial differential quadrature method, and solving the dispersed nonlinear equation set by using Newton-Raphson iteration. The calculation method of the invention carries out discrete solution on the above formula. Simulating the step length by an implicit trapezoidal method for 1 mu s, and subsequently calling the step length as a small step length; the serial differential integration method and the calculation method simulate the step length of 10 mu s, and are subsequently called as large step length. Set convergence accuracy to 10-4The three algorithms only need to iterate twice within each step length to complete convergence.
Fig. 3 is a diagram showing voltages at the head end and the tail end of the power transmission line obtained by simulation using the calculation method of the present invention, and fig. 4 is a diagram showing an error curve of simulation results of the calculation method (large step size) and the implicit trapezoidal method (small step size) of the present invention. From fig. 3 and fig. 4, it can be seen that even though the step size of the differential integration method is 10 times that of the implicit trapezoidal method, the accuracy of the simulation error between the step size of the implicit trapezoidal method and the step size of the implicit trapezoidal method is very small, which indicates that the differential integration method can ensure the simulation accuracy on the premise of adopting a large step size. It should be noted that since the serial differential integration method and the calculation method of the present invention are consistent in the selection of the discrete method and the step size, the accuracy error of the two methods is small, and the accuracy is considered to be the same, and no error curve is given.
The ratio of the calculation time of the calculation method of the invention to the calculation time of other methods is defined as an acceleration ratio, and the specific definition is as follows: k1The acceleration ratio (the step length is consistent) of the calculation method and the traditional serial differential integration method is shown; k2Represents the acceleration ratio of the algorithm and the small step implicit trapezoidal method; k2The acceleration ratio of the traditional serial differential integration method and the small-step implicit trapezoidal method is shown. The large step size and the small step size are specifically defined in the embodiment.
The acceleration ratio is given in tabular form:
TABLE 2 acceleration ratio
Figure BDA0002651329910000172
As can be seen from Table 2, in the conventional serial manner, the triangular decomposition of the s × m dimensional matrix is performed in each iterative solution process, and even if large step calculation is adopted, the triangular decomposition of the dimension-increasing matrix causes low calculation efficiency, and the K reflected in the table above is3The value, i.e. its acceleration ratio, is not significant. For the calculation method, on the basis of adopting large step length calculation, the triangular decomposition of the s multiplied by m dimensional matrix is converted into the prior generation operation of the (s-1) multiplied by m dimensional matrix and the triangular decomposition of the m dimensional matrix in each iteration process, and compared with the traditional serial differential integration method, the acceleration advantage is accumulated along with the iteration of one time and another, and finally, a considerable acceleration ratio is formed and is reflected as K in a table1，K2The value is obtained. Therefore, the algorithm of the invention overcomes the problem of increased calculated amount caused by introducing inner points in the serial algorithm of the differential quadrature method, and improves the electromagnetic transient calculation speed based on the differential quadrature method.
The comparison result of the precision error diagram and the acceleration ratio shows that the differential integration method can still keep better calculation precision under large-step simulation. The acceleration ratio shows that the calculation method can effectively improve the calculation efficiency based on the differential integration method, so that the advantage of large-step simulation is fully exerted.

Claims (7)

1. The electromagnetic transient rapid simulation method based on the differential quadrature method is characterized by comprising the following steps:
step 1: the electromagnetic transient process is a transformation process of voltage and current caused by the change of a magnetic field and an electric field, and is represented by a set of differential equations with initial values as follows:
Figure FDA0002651329900000011
wherein: x is formed by Rm×1X is the set of state variables in the system, g (t) is related to the time variable, x0Then is t0(t ═ 0) the initial value of the state variable x at time;
the differential quadrature method is as follows:
if the function f (x) is smoothly derivable over an interval, f (x) is at grid point ciThe derivative over i e (1, s) can be represented by a linear weighted sum of the function values at all grid points over the interval, i.e.:
Figure FDA0002651329900000012
order: gij∈G (3)
Wherein: f. of(1)(ci) I e (1, s) represents the derivative value of the function f (x) at any grid point over the interval, gijWeighting coefficients for differential integration, f (c)j) J e (1, s) represents a function value of the function f (x) at any grid point in the interval, and s is the grid point number in the interval;
the mathematical expression of the differential integration method is derived as follows:
Figure FDA0002651329900000013
in the formula, tnIs a starting time, tn+1For the end time, x is the set of state variables in the system, and x ∈ Rm×1，xnIs tnValue of the time-of-day state variable x, xn+1Is tn+1The value of the state variable to be solved at the moment,
Figure FDA0002651329900000014
indicates an arbitrary grid point c on the intervaliThe value of the state variable, h is the simulation step length, aijFor the weighting coefficients, A ≡ a is definedij]And A ═ G-1And s is the number of inner points;
discretization (1) by differential integration and written in matrix form:
Figure FDA0002651329900000021
step 2: using A ═ G-1Equally transforming equation (5) in step 1 may result:
Figure FDA0002651329900000022
in formula (6), G is a weighting coefficient matrix of differential quadrature method, ImIs m-dimensional unit matrix, e is s-dimensional unit column vector, xnIs tnThe value of a state variable x at the moment, and h is a simulation step length;
in the formula (6), the mathematical expression is explained as follows:
Figure FDA0002651329900000023
(1) when formula (1) is described as a non-linear problem, the definition:
Figure FDA0002651329900000024
in formula (7): s is the number of the inner points,
Figure FDA0002651329900000025
indicates the grid point c on the intervalkValue of state variable
Solving a matrix equation (6) after dispersion by using a Newton iteration method, and taking an average value vector of s inner points as an initial value of the Jacobian matrix to obtain:
Figure FDA0002651329900000026
in formula (8):
Figure FDA0002651329900000027
in the formula: g is a weighting coefficient matrix of the differential quadrature method, ImIs an m-dimensional identity matrix, h is a simulation step length, IsIs an s-dimensional unit matrix, xi is the iteration number,
Figure FDA0002651329900000028
denotes the ξ -th average of the interior points in the iteration,
Figure FDA0002651329900000029
denotes the ξ +1 th iteration at grid point ckThe value of the upper state variable x;
in formula (9):
Figure FDA0002651329900000031
and step 3: the V transform form of the differential quadrature method is:
Figure FDA0002651329900000032
in the formula (10), G is a differential multiplication coefficient matrix, V is a Van der Mond matrix, and V is-1Being the inverse of the van der Mond matrix
Order:
Figure FDA0002651329900000033
in the formula (11), V-1Is the inverse of the van der Mond matrix, ImIs an m-dimensional identity matrix, and
Figure FDA0002651329900000034
applying the formulas (9) and (10) to the formula (8) in the step 2 to obtain:
Figure FDA0002651329900000035
in the formula:
Figure FDA0002651329900000036
in the formula (13), ImIs an m-dimensional identity matrix, IsIs an s-dimensional unit matrix and is,
Figure FDA0002651329900000037
xi is the number of iterations,
and 4, step 4:
Figure FDA0002651329900000038
unfolding and blocking formula (12) obtained in step 3 to obtain:
Figure FDA0002651329900000041
Figure FDA0002651329900000042
Figure FDA0002651329900000043
in the formula (15), ImIs an m-dimensional unit matrix and is a matrix,
Figure FDA0002651329900000044
xi is the iteration number;
equation (12) can be finally written as:
Figure FDA0002651329900000045
observed to obtain, H2Is a lower triangular matrix, so that y can be transformed by simple prior generation operation2Is expressed as y1Without the need to perform complex trigonometric decomposition, i.e.:
y2＝φ(y1)； (19)
bringing formula (19) into H3y1+H4y2＝r2Finally, the formula (18) is converted to contain y alone1The algebraic system of equations, i.e.:
Figure FDA0002651329900000046
Q∈Rm×m (21)
q is a coefficient matrix obtained by calculation;
using LU decomposition to solve y1Through y2＝φ(y1) To give out y2Finally, solving the triangle x in an inverse manner; calculating the next step length by judging whether the current step length is smaller than the convergence precision or not, if the current step length is smaller than the convergence precision, calculating the next step length, and if the current step length is larger than the convergence precision, performing the next iteration by taking the average value of the interior points as an initial value until the current step length is converged into the precision range;
and 5: when the problem described by equation (1) is a linear time-varying problem, that is:
f(x)＝Utx (22)
Utfor a matrix of coefficients relating only to time variables, we define:
Figure FDA0002651329900000051
in the formula (23), s is the number of inner points, U (t)n+cih) I ∈ (1, s) is tn+ciMoment of coefficient at h timeValue of the matrix, ciIs a grid point, h is a simulation step length;
order:
Figure FDA0002651329900000052
in the formula (24), V-1Is the inverse of the van der Mond matrix, ImIs an m-dimensional unit matrix and is a matrix,
Figure FDA0002651329900000053
the linear time-varying differential equation is dispersed by using a differential quadrature method and is introduced into U to obtain:
Figure FDA0002651329900000054
wherein:
Figure FDA0002651329900000055
Figure FDA0002651329900000056
Figure FDA0002651329900000057
in the formula (27), V-1Is the inverse of the van der Mond matrix, ImIs an m-dimensional unit matrix;
in the formula (28), G is a weighting coefficient matrix of differential integration method, ImIs m-dimensional unit matrix, e is s-dimensional unit column vector, xnIs tnThe value of the state variable x at time, h is the simulation step length, g (t) ═ gT(tn+c1h) … gT(tn+csh)]T；
Matrix partitioning processing is performed on equation (25) to obtain:
Figure FDA0002651329900000058
Figure FDA0002651329900000059
Figure FDA00026513299000000510
equation (24) is finally written as:
Figure FDA00026513299000000511
H2is a lower triangular matrix, so that y can be transformed by simple prior generation operation2Is expressed as y1The linear expression of (a), namely:
y2＝φ(y1)； (33)
bringing formula (33) into H3y1+H4y2＝r2Finally, the formula (32) is converted to contain y alone1The algebraic system of equations, i.e.:
Figure FDA0002651329900000061
Q∈Rm×m (35)
using LU decomposition to solve y1Through y2＝φ(y1) To give out y2Finally, the solution is reversely solved by the formula (24)
Figure FDA0002651329900000062
The calculation of the next step is started.
2. The method for electromagnetic transient rapid simulation based on differential integration according to claim 1, wherein:
in step 1, the G matrix expression in formula (2) is:
Figure FDA0002651329900000063
wherein:
Figure FDA0002651329900000064
αs＝[α1,α2,…,αs]T＝(1/s)V-1cs；
Figure FDA0002651329900000065
wherein V is a Van der Monde matrix, s is the number of interval grid points, ckK is an interval grid point (1, s);
for the selection of the grid points, there are Legendre grids, Chebyshev grids, uniform grids and the like, and the function variables in the formula (5) are as follows:
Figure FDA0002651329900000066
3. the method for electromagnetic transient rapid simulation based on differential integration according to claim 1, wherein:
in the step 2, substituting a ═ G-1Realizing equal transformation; this has the advantage of converting formula (4) to formula (5) in a simpler form.
4. The method for electromagnetic transient rapid simulation based on differential integration according to claim 1, wherein:
in step 3, an iterative format in which the average phasor of s interior points is used as an initial value is adopted, so that the following formula holds:
Figure FDA0002651329900000071
in the formula (39), V is a van der Mond matrix, IsIs an s-dimensional identity matrix, ImIs an m-dimensional identity matrix, and
Figure FDA0002651329900000072
xi is the iteration number;
and finally in the form of formula (12) in step 3.
5. The method for electromagnetic transient rapid simulation based on differential integration according to claim 1, wherein:
in the step 4, the process of the step,
Figure FDA0002651329900000073
Figure FDA0002651329900000074
H3＝[a1,sIm]； (42)
Figure FDA0002651329900000075
Figure FDA0002651329900000076
Figure FDA0002651329900000077
Figure FDA0002651329900000078
Figure FDA0002651329900000079
in the formula ImIs an m-dimensional unit matrix and is a matrix,
Figure FDA00026513299000000710
ξ is the number of iterations.
6. The method for electromagnetic transient rapid simulation based on differential integration according to claim 1, wherein:
in the step 5:
Figure FDA0002651329900000081
Figure FDA0002651329900000082
H3＝[a1,sIm]； (50)
H4＝[0 … -U] (51)
Figure FDA0002651329900000083
Figure FDA0002651329900000084
Figure FDA0002651329900000085
Figure FDA0002651329900000086
in the formula ImIs an m-dimensional unit matrix and is a matrix,
Figure FDA0002651329900000087
7. the differential quadrature method is applied to electromagnetic transient simulation calculation.
CN202010871833.0A 2020-08-26 2020-08-26 Matrix calculation method based on differential product method and V transformation in electromagnetic transient simulation Active CN112069668B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010871833.0A CN112069668B (en) 2020-08-26 2020-08-26 Matrix calculation method based on differential product method and V transformation in electromagnetic transient simulation

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010871833.0A CN112069668B (en) 2020-08-26 2020-08-26 Matrix calculation method based on differential product method and V transformation in electromagnetic transient simulation

Publications (2)

Publication Number Publication Date
CN112069668A true CN112069668A (en) 2020-12-11
CN112069668B CN112069668B (en) 2023-06-30

Family

ID=73660079

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010871833.0A Active CN112069668B (en) 2020-08-26 2020-08-26 Matrix calculation method based on differential product method and V transformation in electromagnetic transient simulation

Country Status (1)

Country Link
CN (1) CN112069668B (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116720338A (en) * 2023-05-30 2023-09-08 杭州盛星能源技术有限公司 Electromagnetic transient parallel iteration real-time simulation compensation method and device

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP0884690A2 (en) * 1997-06-11 1998-12-16 Siemens Aktiengesellschaft Computer aided simulation method for determining the electromagnetic field of an object
US20100303205A1 (en) * 2009-05-26 2010-12-02 Agile Planet Inc System and Method for Radiation Therapy Imaging and Treatment Workflow Scheduling and Optimization
NL2005687C2 (en) * 2010-11-12 2012-05-15 Tu Delft Method for determining a spring constant for a deformable scanning probe microscope element, and scanning probe microscope and calibration device arranged for determining a spring constant for a probe element.
CN103646152A (en) * 2013-12-23 2014-03-19 南方电网科学研究院有限责任公司 Electromagnetic transient simulation method of electric system based on matrix exponential
CN105404610A (en) * 2015-10-23 2016-03-16 三峡大学 Electromagnetic transient calculation method based on two-stage three-order single diagonally implicit Runge-Kutta method
CN105468864A (en) * 2015-12-14 2016-04-06 三峡大学 Electromagnetic transient numerical computation method of high-voltage power transmission line based on increment dimension precise integration
US20170250540A1 (en) * 2009-09-15 2017-08-31 Rajiv Kumar Varma Multivariable modulator controller for power generation facility
CN107679287A (en) * 2017-09-11 2018-02-09 三峡大学 Electro-magnetic transient numerical computation method based on the implicit Taylor series Method of the rank of 3 step 4
CN109472756A (en) * 2018-11-15 2019-03-15 昆明理工大学 Image de-noising method based on shearing wave conversion and with directionality local Wiener filtering

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP0884690A2 (en) * 1997-06-11 1998-12-16 Siemens Aktiengesellschaft Computer aided simulation method for determining the electromagnetic field of an object
US20100303205A1 (en) * 2009-05-26 2010-12-02 Agile Planet Inc System and Method for Radiation Therapy Imaging and Treatment Workflow Scheduling and Optimization
US20170250540A1 (en) * 2009-09-15 2017-08-31 Rajiv Kumar Varma Multivariable modulator controller for power generation facility
NL2005687C2 (en) * 2010-11-12 2012-05-15 Tu Delft Method for determining a spring constant for a deformable scanning probe microscope element, and scanning probe microscope and calibration device arranged for determining a spring constant for a probe element.
CN103646152A (en) * 2013-12-23 2014-03-19 南方电网科学研究院有限责任公司 Electromagnetic transient simulation method of electric system based on matrix exponential
CN105404610A (en) * 2015-10-23 2016-03-16 三峡大学 Electromagnetic transient calculation method based on two-stage three-order single diagonally implicit Runge-Kutta method
CN105468864A (en) * 2015-12-14 2016-04-06 三峡大学 Electromagnetic transient numerical computation method of high-voltage power transmission line based on increment dimension precise integration
CN107679287A (en) * 2017-09-11 2018-02-09 三峡大学 Electro-magnetic transient numerical computation method based on the implicit Taylor series Method of the rank of 3 step 4
CN109472756A (en) * 2018-11-15 2019-03-15 昆明理工大学 Image de-noising method based on shearing wave conversion and with directionality local Wiener filtering

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
JALILI-MARANDI V等: "Interfacing techniques for transient stability and electromagnetic transient programs", 《IEEE TRANSACTIONS ON POWER DELIVERY》 *
汪芳宗等: "基于矩阵对角化的电磁暂态时间并行计算方法", 《电网技术》 *
王成山 等: "电力系统电磁暂态仿真算法研究进展", 《电力系统自动化》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116720338A (en) * 2023-05-30 2023-09-08 杭州盛星能源技术有限公司 Electromagnetic transient parallel iteration real-time simulation compensation method and device
CN116720338B (en) * 2023-05-30 2024-02-02 杭州盛星能源技术有限公司 Electromagnetic transient parallel iteration real-time simulation compensation method and device

Also Published As

Publication number Publication date
CN112069668B (en) 2023-06-30

Similar Documents

Publication Publication Date Title
Gustavsen et al. A robust approach for system identification in the frequency domain
Chiprout et al. Asymptotic waveform evaluation
Nelms et al. Using a personal computer to teach power system transients
CN114091363B (en) Quantum algorithm-based computational fluid dynamics simulation method, device and equipment
Polyuga et al. Effort-and flow-constraint reduction methods for structure preserving model reduction of port-Hamiltonian systems
Galiana et al. On the application of a pre-conditioned conjugate gradient algorithm to power network analysis
CN113191105A (en) Electrical simulation method based on distributed parallel operation method
CN112069668A (en) Electromagnetic transient rapid simulation method based on differential quadrature method
Zhang et al. Time‐weighted MR for switched LPV systems via balanced realisation
Farhan et al. Fast simulation of microwave circuits with nonlinear terminations using high-order stable methods
Ramirez et al. Reduced-sample numerical Laplace transform for transient and steady-state simulations: Application to networks involving power electronic converters
Zhang et al. Fast nonlinear model order reduction via associated transforms of high-order Volterra transfer functions
CN106934123A (en) A kind of circuit transient response computational methods based on recursive convolution
CN107944682B (en) Load flow calculation admittance matrix calculation method based on Matlab matrix operation
Lang et al. The harmonic balance method
Zhang et al. Implicit integration factor method for the nonlinear Dirac equation
Bechtold et al. Model order reduction: an advanced, efficient and automated computational tool for microsystems
Kolka et al. Program for multi-domain symbolic analysis
Bhawal et al. New results and techniques for computation of stored energy in lossless/all-pass systems
Cingoz et al. Reduced order, high-fidelity modeling of energy storage units in vehicular power systems
CN107665184B (en) Load flow calculation admittance matrix calculation method based on incidence matrix operation
Liu et al. Krylov subspace method for passive reduced‐order modelling of electrical networks
Wang et al. A Novel Electro-magnetic Transient Analysis Method Based on the Projection Operator for Power System
Bizzarri et al. Reduction of harmonic balance equations through Galerkin's method
Zhang et al. A Time-Domain Explicit Integration Algorithm for Fast Overvoltage Computation of High-Voltage Transmission Line

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant