CN108054753B - Direct-drive wind power plant cluster division method considering low-voltage ride through characteristics - Google Patents

Direct-drive wind power plant cluster division method considering low-voltage ride through characteristics Download PDF

Info

Publication number
CN108054753B
CN108054753B CN201711394587.9A CN201711394587A CN108054753B CN 108054753 B CN108054753 B CN 108054753B CN 201711394587 A CN201711394587 A CN 201711394587A CN 108054753 B CN108054753 B CN 108054753B
Authority
CN
China
Prior art keywords
wind turbine
direct
value
turbine generator
voltage
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.)
Expired - Fee Related
Application number
CN201711394587.9A
Other languages
Chinese (zh)
Other versions
CN108054753A (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.)
Hefei University of Technology
Original Assignee
Hefei University of Technology
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 Hefei University of Technology filed Critical Hefei University of Technology
Priority to CN201711394587.9A priority Critical patent/CN108054753B/en
Publication of CN108054753A publication Critical patent/CN108054753A/en
Application granted granted Critical
Publication of CN108054753B publication Critical patent/CN108054753B/en
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • H—ELECTRICITY
    • H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
    • H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
    • H02J3/00—Circuit arrangements for ac mains or ac distribution networks
    • H—ELECTRICITY
    • H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
    • H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
    • H02J3/00—Circuit arrangements for ac mains or ac distribution networks
    • H02J3/24—Arrangements for preventing or reducing oscillations of power in networks
    • H02J3/386—
    • H—ELECTRICITY
    • H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
    • H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
    • H02J2203/00—Indexing scheme relating to details of circuit arrangements for AC mains or AC distribution networks
    • H02J2203/20—Simulating, e g planning, reliability check, modelling or computer assisted design [CAD]
    • 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
    • Y02E10/00—Energy generation through renewable energy sources
    • Y02E10/70—Wind energy
    • Y02E10/76—Power conversion electric or electronic aspects

Landscapes

  • Engineering & Computer Science (AREA)
  • Power Engineering (AREA)
  • Control Of Eletrric Generators (AREA)
  • Wind Motors (AREA)

Abstract

The invention discloses a direct-drive wind power plant cluster division method considering low-voltage ride through characteristics, which comprises the following steps of: 1. establishing a detailed model of a direct-drive wind power plant; 2. collecting input wind energy of each wind turbine generator set; 3. collecting terminal voltage values of each wind turbine generator set at the moment of failure; 4. setting an active current reference value of a grid-side converter when the voltage of a power grid drops; 5. judging the conduction condition of each wind turbine generator set unloading circuit; 6. and the direct-drive wind power plant cluster division is completed by applying an immune random K value and an improved K mean value clustering algorithm of a sensitive clustering center. The method can accurately represent the conduction condition of each unloading circuit of the wind power generation set in the direct-drive wind power plant, and improves the dividing quality of the wind power plant groups, so that the fitting precision of the equivalent model of the direct-drive wind power plant to the output external characteristics of the wind power plant is effectively improved.

Description

Direct-drive wind power plant cluster division method considering low-voltage ride through characteristics
Technical Field
The invention relates to the technical field of equivalent modeling of direct-drive wind power plants, in particular to a method for dividing a cluster of a direct-drive wind power plant considering low-voltage ride through characteristics.
Background
The direct-drive permanent magnet wind turbine generator system saves a gear box with high failure rate, has the advantages of small mechanical loss, high operation efficiency, low maintenance cost, wide speed regulation range and the like, is connected with a power grid through a full-power converter, is not directly coupled with the power grid, has stronger low-voltage ride through capability compared with a double-fed wind turbine generator system, and increasingly becomes the mainstream wind turbine type of the current wind power plant. With the rapid development of the wind power generation technology, the scale of the grid-connected wind power plant is gradually enlarged, and if each wind power generation unit adopts a detailed model, various problems of too long simulation calculation time, too large occupied memory, difficult convergence of tide and the like occur. Therefore, it is necessary to establish a proper equivalent model of the wind power plant.
At present, equivalent methods of direct-drive wind power plants mainly comprise a single-machine characterization method and a multi-machine characterization method. The single-machine characterization method is suitable for the wind power plant with smaller scale, more centralized distribution of the wind turbine generators on the geographical position and smaller difference between the electrical distance of the wind turbine generators and the input wind speed. As wind farms are continually scaled up, the single-machine characterization method usually produces large errors. The multi-machine characterization method adds a coherent cluster division link, and the main idea is to divide the cluster by adopting a clustering algorithm based on the cluster division principle that the wind turbine generators have similar operating points. At present, for a direct-drive permanent magnet wind turbine, common grouping indexes mainly include input wind speed, actually measured active data of the wind turbine, and comprehensive grouping indexes composed of wind speed, active power, voltage, current and the like. The kmeans clustering algorithm is the most commonly adopted clustering algorithm for cluster division, and has the advantages of simplicity, easiness in implementation, high running speed and capability of processing a large-scale data set.
However, when a cluster is divided by using the existing multi-machine characterization method, the voltage drop difference at the wind turbine generator end is mostly not fully considered, the influence of the conduction condition of the unloading circuit on the transient characteristic difference of the direct-drive permanent magnet wind turbine generator is ignored, and the output external characteristic of the wind power plant cannot be accurately fitted. In addition, the traditional kmeans clustering algorithm has the advantages that the k value of the clustering number needs to be given in advance; the method is seriously dependent on the selection of an initial clustering center, so that a clustering result is easy to fall into a local optimal solution; easily influenced by isolated point data and the like. The defects of the traditional kmeans clustering algorithm greatly reduce the clustering quality, and cause inaccuracy and instability of cluster division results.
Disclosure of Invention
The invention provides a direct-drive wind power plant cluster division method considering low-voltage ride through characteristics in order to avoid the defects of the prior art, so that the conduction condition of each wind power plant unloading circuit in a wind power plant can be accurately represented, the cluster division quality of the direct-drive wind power plant is improved, and the fitting precision of an equivalent model of the direct-drive wind power plant on the output external characteristics of the wind power plant can be effectively improved.
The invention adopts the following technical scheme for solving the technical problems:
the invention relates to a direct-drive wind power plant cluster division method considering low-voltage ride through characteristics, which is characterized by comprising the following steps of:
step 1, establishing a detailed model of a direct-drive wind power plant;
the direct-drive wind power plant is composed of n direct-drive permanent magnet wind power units with the same type number, and the detailed model comprises the following components: the system comprises a single machine model, a current collection circuit model, a generator-end transformer model and a main transformer model of each wind turbine generator set in a wind power plant; wherein, the stand-alone model of wind turbine generator system includes: the system comprises a wind turbine model, a permanent magnet synchronous generator model, a full-power converter and a control system model thereof, a variable pitch control system model, an unloading circuit and a control system model thereof;
step 2, collecting input wind energy of each wind turbine generator set;
collecting the real-time wind speeds of n direct-drive permanent magnet wind turbine generators of the direct-drive wind power plant, and recording the input wind energy of each wind turbine generator as { P }1w,P2w,···,Pjw,···,Pnw}，PjwRepresenting the input wind energy of the jth wind power generation set;
step 3, collecting terminal voltage values of each wind turbine generator at the moment of failure;
set at t0At the moment, a three-phase short-circuit fault is arranged at an outlet of the direct-drive wind power plant, and t is acquired0The terminal voltage value of each wind turbine generator set at any moment is marked as { U1(t0),U2(t0),···,Uj(t0),···,Un(t0)}，Uj(t0) Represents t0The voltage value of the generator end of the jth wind turbine generator set at the moment is taken as a per unit value; at t1The three-phase short-circuit fault is removed at any moment;
step 4, setting an active current reference value of the grid-side converter when the voltage of the power grid drops;
step 4.1, obtaining the instruction value i of the grid-side converter direct-current voltage PI regulator in the steady state by using the formula (1)dref1：
Figure BDA0001518163920000021
In the formula (1), KdpAnd KdiProportional coefficients and integral coefficients of a direct current bus voltage control outer ring in the grid-side converter are respectively set; u shapedcTaking a per unit value for the direct-current bus voltage of the wind turbine generator in a steady state;
Figure BDA0001518163920000022
taking a per unit value as a direct current bus voltage reference value of the wind turbine generator;
step 4.2, during the three-phase short circuit fault of the direct-drive wind power plant, when the voltage of the grid-connected point of the wind turbine generator is lower than 20% of the nominal voltage, the wind turbine generator is cut out of the power grid;
when the grid-connected point voltage of the wind turbine generator is within 20% -90% of the nominal voltage, calculating a reactive current reference value i of the grid-side converter by using a formula (2)qref2(t)：
iqref2(t)＝1.5(0.9-Ug(t))IN(0.2pu≤Ug(t)≤0.9pu) (2)
In the formula (2), Ug(t) is a per unit value of the grid-connected point voltage of the wind turbine generator; i isNThe voltage is a per unit value of rated current of the grid-side converter;
step 4.3, obtaining an active current reference value limiting value i by using the formula (3)dref2(t)：
Figure BDA0001518163920000031
In the formula (3), imaxThe maximum current value allowed by the grid-side converter;
step 4.4, selecting the instruction value i of the direct current voltage PI regulator in the steady statedref1And an active current reference value limit value idref2(t) the smaller value is used as the active current reference value i of the grid-side converter when the voltage of the power grid dropsdref(t)；
Step 5, judging the conduction condition of each wind turbine generator unloading circuit;
step 5.1, at the t1At any moment, the direct current bus voltage U of the jth wind turbine generator set network side power converter is obtained by using the formula (4)jdc(t1)：
Figure BDA0001518163920000032
In the formula (4), CjIs the direct current bus capacitance value of the jth wind power generator set，SBThe reference capacity of the grid-side converter;
step 5.2, compare t1DC bus voltage U of jth wind turbine generator set at momentjdc(t1) And unloading circuit action threshold Udc_inSize of (1), if Ujdc(t1) Greater than Udc_inIf so, judging that the jth wind turbine generator set unloading circuit is conducted; otherwise, judging that the unloading circuit of the jth wind turbine generator set is not conducted;
step 6, dividing the direct-drive wind power plant cluster by applying an immune random K value and a sensitive clustering center improved K-means clustering algorithm;
step 6.1, dividing all the wind turbine generators which are conducted by the unloading circuits in the n direct-drive permanent magnet wind turbine generators into a cluster;
and 6.2, dividing m wind turbine generators with the rest of non-conducted unloading circuits into K wind turbine groups by taking the terminal voltage values of the wind turbine generators at the fault moment as grouping indexes and applying an immune random K value and a sensitive clustering center to improve a K-means clustering algorithm.
The method for dividing the direct-drive wind power plant cluster is also characterized in that in the step 6.2, an immune random K value and sensitive clustering center improved K-means clustering algorithm is carried out according to the following steps:
step 6.2.1, the terminal voltage values { U } of the m wind turbine generators with the residual unloading circuits not conducted1(t0),U2(t0),···,Up(t0),···,Uq(t0),···,Um(t0) As a sample data set, it is denoted as S ═ x1,x2,···,xp,···,xq,···,xm}; wherein, Up(t0) And Uq(t0) Respectively representing the terminal voltage values of the p-th and q-th wind turbine generators, taking the per unit value, and converting the U value into the voltage valuep(t0) And Uq(t0) Respectively denoted as data object xpAnd xq，p,q＝1,2,···,m,p≠q；
Step 6.2.2, calculate data object xpAnd xqIs (x) ofp,xq)；
Step 6.2.3, obtaining the average distance MeanDist between any two data objects in the sample data set S by using the formula (5):
Figure BDA0001518163920000041
in the formula (5), the reaction mixture is,
Figure BDA0001518163920000042
the number of combinations of 2 data objects is selected from n data objects;
step 6.2.4, define: with data object xpThe area centered at the average distance MeanDist is called data object xpOf the data object x, the data object xpIs called data object xpA density parameter based on the distance MeanDist;
obtaining data object x using equation (6)pDensity parameter (x) ofpMeanDist) to obtain density parameters for m data objects:
Figure BDA0001518163920000043
in the formula (6), u (MeanDist-dist (x)p,xq) Represents a function and has:
Figure BDA0001518163920000044
step 6.2.5, adding the density parameters of the first M data objects with the largest density parameter in the density parameters of the M data objects into the candidate point set D, where D is { density (x)p,MeanDist),p＝1,2,···,M}；
Step 6.2.6, defining a loop variable r, and initializing r to be 1; defining the clustering number as k, and initializing k to be 1; defining the maximum similarity between the initial average classes as AMSr-1And initializing the AMSr-1(ii) a Defining the cluster center set as Ar-1And initialized to an empty set; with the saidThe alternative point set D is used as the r-1 st alternative point set Dr-1；
Step 6.2.7, from the r-1 st candidate point set Dr-1Selecting the mean value of two data objects with the maximum density parameter as the r-th clustering center crPut into clustering center set Ar-1In the method, the r-th clustering center set A is obtainedr(ii) a Simultaneously, the density parameters corresponding to the two selected data objects are selected from the r-1 st candidate point set Dr-1Deleting to obtain the r-th candidate point set Dr；
Step 6.2.8, from the r-th candidate point set DrThe (r) th cluster center set ArThe data object with the farthest cluster center distance in (1) is used as the (r + 1) th cluster center cr+1Put into the r-th clustering center set ArIn the method, the r +1 th clustering center set A is obtainedr+1(ii) a At the same time, cluster the r +1 th cluster center cr+1Corresponding density parameter is selected from the r-th candidate point set DrDeleting to obtain the r +1 th alternative point set Dr+1；
Step 6.2.9, assigning k +1 to k;
step 6.2.10, calculating the rest data objects in the sample data set S and the r +1 th clustering center set A respectivelyr+1The Euclidean distance of each clustering center is determined, and each data object is distributed to the class where the clustering center with the nearest Euclidean distance is located, so that k classes are obtained;
step 6.2.11, obtaining the data object x of any ith class in the k classes by using the formula (8)pClustering center c to class iiMean value s of the distance betweeni：
Figure BDA0001518163920000051
In the formula (8), PiThe total number of the data objects in the ith class; 1,2, ·, k;
step 6.2.12, obtaining the clustering center c of the ith class by using the formula (9)iCluster center c with class jjA distance d betweeni,j：
di,j＝dist(ci,cj) (9)
In formula (9), i, j ≠ 1,2, ·, k, i ≠ j;
step 6.2.13, obtaining the maximum similarity AMS between the r-th average class by using the formula (10)r：
Figure BDA0001518163920000052
Step 6.2.14, determine AMSr＜AMSr-1If true, go to step 6.2.15; otherwise, go to step 6.2.19;
step 6.2.15, obtaining the updated cluster center c using equation (11)i′：
Figure BDA0001518163920000053
In the formula (11), PiIs the total number of data objects in class i, xpIs any data object of the ith class, i ═ 1,2, ·, k;
step 6.2.16, calculating the r +1 th candidate point set Dr+1The sum of Euclidean distances between the data object corresponding to any one density parameter and all updated clustering centers is obtained, so that a set of the sums of the Euclidean distances between the data object corresponding to all density parameters and all updated clustering centers is obtained, and the data object corresponding to the maximum value is selected from the set of the sums of the Euclidean distances to serve as the (r + 2) th clustering center cr+2And put into the r +1 th clustering center set Ar+1To obtain the r +2 th clustering center set Ar+2；
Step 6.2.17, clustering the r +2 th clustering center cr+2The corresponding density parameter is selected from the r +1 th candidate point set Dr+1Deleting to obtain the r +2 th alternative point set Dr+2；
6.2.18, assigning r +1 to r, and turning to step 6.2.9;
step 6.2.19, use AMSr-1Corresponding k clustering centers are used as kmeans clustersAnd performing a kmeans clustering algorithm on the sample data set S by using an initial clustering center of the class algorithm to obtain k machine groups.
Compared with the prior art, the invention has the beneficial effects that:
1. according to the direct-drive wind power plant cluster division method, the unloading circuit conduction condition of each wind power plant in the direct-drive wind power plant is used as a clustering basis, the voltage value of the generator end of each wind power plant in the direct-drive wind power plant at the fault moment is used as a clustering index, and the direct-drive wind power plant cluster division is completed by applying an immune random K value and an improved K mean value clustering algorithm of a sensitive clustering center, so that the output external characteristic of the direct-drive wind power plant is effectively fitted, the actual running state of the direct-drive wind power plant is accurately reflected, and the accuracy and the effectiveness of the direct-drive wind power plant cluster division result are.
2. According to the method, the reactive current reference value of the grid-side converter is adjusted according to the drop depth of the grid-connected point voltage of the wind turbine generator when the short-circuit fault occurs on the grid side, so that the active current reference value of the grid-side converter is set, the switching process of the control strategy of the direct-drive permanent magnet wind turbine generator in the low-voltage ride-through mode is reflected, and the requirements that the wind turbine generator does not break off the grid and continuously operates and provides reactive support for the grid when the short-circuit fault occurs on the grid side are met.
3. The improved K-means clustering algorithm of the immune random K value and the sensitive clustering center is adopted, and the optimal clustering number K is automatically determined by comparing the value of the maximum similarity index AMS between average classes, so that the problem that the K value of the traditional kmeans algorithm needs to be given in advance is solved, and the clustering quality is improved.
4. The improved K mean value clustering algorithm of the immune random K value and the sensitive clustering center can dynamically adjust the current clustering center and dynamically add the next clustering center, and the finally obtained K clustering centers are the optimal K clustering centers obtained through K-1 times of dynamic distribution and are generated according to the characteristics of data objects strictly, so that the volatility of clustering results is effectively avoided, the problem of the sensitive initial clustering centers of the traditional kmeans algorithm is solved, the global searching capability is improved, and the accuracy and the stability of clustering are improved to a greater extent.
5. The improved K-means clustering algorithm of the immune random K value and the sensitive clustering center can independently divide the isolated point data into a cluster, reduce the influence of the isolated point data on cluster division results, avoid the clustering center from moving away from a data dense area to tend to the isolated point data, and improve the clustering quality.
Drawings
FIG. 1 is a flow chart of a direct-drive wind farm cluster partitioning method proposed by the present invention;
FIG. 2 is a topological structure diagram of a direct-drive wind power plant in the invention;
FIG. 3 is a topological structure diagram of a direct-drive permanent magnet wind turbine generator set in the invention;
fig. 4 is a block diagram of a grid-side converter control strategy in the present invention;
FIG. 5 is a flow chart of an improved K-means clustering algorithm for immune random K values and sensitive clustering centers in the invention.
Detailed Description
Embodiments of the present invention will be described in detail below with reference to the accompanying drawings.
In the embodiment, the direct-drive wind power plant cluster division method considering the low-voltage ride-through characteristic considers the influence of the unloading circuit conduction condition on the transient response characteristic of the direct-drive permanent-magnet wind power plant when a fault occurs on the power grid side, judges the unloading circuit conduction condition of each wind power plant according to the change of the direct-current bus voltage, uses the voltage value at the generator end of each wind power plant in the direct-drive wind power plant at the fault time as a cluster division index, and uses the immune random K value and the improved K-means clustering algorithm of the sensitive clustering center to complete the cluster division of the direct-drive wind power plant, so that the cluster division accuracy and stability are improved, the multi-machine equivalent modeling method of the wind power plant is optimized, and a foundation is laid for the electromagnetic.
The specific flow chart of the direct-drive wind power plant cluster dividing method provided by the invention is shown in FIG. 1, and mainly comprises the following steps:
step 1, establishing a detailed model of a direct-drive wind power plant;
the direct-drive wind power plant consists of n direct-drive permanent magnet wind power units with the same type number, and the detailed model comprises the following components: the system comprises a single machine model, a current collection circuit model, a generator-end transformer model and a main transformer model of each wind turbine generator set in a wind power plant; the topological structure diagram of the direct-drive wind power plant is shown in fig. 2, each direct-drive permanent magnet wind turbine generator is connected to a medium-voltage bus of the wind power plant through an overhead line after being boosted by a generator-end transformer of the direct-drive permanent magnet wind turbine generator, and then is connected to a power grid through a double-circuit overhead line from an outlet of the direct-drive wind power plant after being boosted secondarily by a main transformer of the wind power;
the topological structure diagram of the direct-drive permanent magnet wind turbine generator is shown in fig. 3, and the single-machine model of the wind turbine generator comprises: the system comprises a wind turbine model, a permanent magnet synchronous generator model, a full-power converter and a control system model thereof, a variable pitch control system model, an unloading circuit and a control system model thereof;
step 2, collecting input wind energy of each wind turbine generator set;
collecting the real-time wind speeds of n direct-drive permanent magnet wind turbine generators of the direct-drive wind power plant every 5 minutes, randomly selecting the real-time wind speeds of 50 sets of each wind turbine generator within continuous 10 hours, calculating the average value, and recording the input wind energy of each wind turbine generator as { P1w,P2w,···,Pjw,···,Pnw}，PjwRepresenting the input wind energy of the jth wind power plant.
Step 3, collecting terminal voltage values of each wind turbine generator at the moment of failure;
set at t0At the moment, a three-phase short-circuit fault is set at an outlet of the direct-drive wind power plant, and t is acquired0The terminal voltage value of each wind turbine generator set at any moment is marked as { U1(t0),U2(t0),···,Uj(t0),···,Un(t0)}；Uj(t0) Represents t0The voltage value of the generator end of the jth wind turbine generator set at the moment is taken as a per unit value; at t1The three-phase short-circuit fault is cut off at any moment;
the three-phase short-circuit fault point is arranged at the outlet of the direct-drive wind power plant, so that the influence of line impedance in the wind power plant on the conduction condition of the unloading circuit of the wind turbine generator can be fully reflected.
Step 4, setting an active current reference value of the grid-side converter when the voltage of the power grid drops;
a block diagram of the grid-side converter control strategy is shown in fig. 4. According to the control strategy shown in fig. 4, both the active current reference value and the reactive current reference value of the grid-side converter in a steady state adopt the instruction value of a voltage outer loop PI regulator; when the voltage of the grid-connected point of the wind turbine generator drops, the reactive current reference value of the grid-side converter is set according to the voltage drop degree of the grid-connected point of the wind turbine generator, so that the active current reference value of the grid-side converter is set.
Step 4.1, when the grid voltage is normal, the grid-side converter executes an active priority maximum power tracking control mode, and the instruction value i of the grid-side converter direct-current voltage PI regulator in a steady state is obtained by using the formula (1)dref1：
Figure BDA0001518163920000081
In the formula (1), KdpAnd KdiProportional coefficients and integral coefficients of a direct current bus voltage control outer ring in the grid-side converter are respectively set; u shapedcTaking a per unit value for the direct-current bus voltage of the wind turbine generator in a steady state;
Figure BDA0001518163920000082
taking a per unit value as a direct current bus voltage reference value of the wind turbine generator;
step 4.2, during the three-phase short circuit fault of the direct-drive wind power plant, when the voltage of the grid-connected point of the wind turbine generator is lower than 20% of the nominal voltage, the wind turbine generator is cut out of the power grid;
when the grid-connected point voltage of the wind turbine generator is within a range of 20% -90% of the nominal voltage, the grid-side converter needs to be switched to a static reactive compensation operation mode with reactive priority because the active priority control mode of the grid-side converter is difficult to fully exert the reactive support capability of the direct-drive permanent magnet wind turbine generator on the power grid. At the moment, the grid-side converter does not need to adopt a PI control mode, the reactive current reference value of the grid-side converter is adjusted according to the drop depth of the grid-connected point voltage of the wind turbine generator, and the grid is calculated by using the formula (2)Reactive current reference value i of side converterqref2(t)：
iqref2(t)＝1.5(0.9-Ug(t))IN(0.2pu≤Ug(t)≤0.9pu) (2)
In the formula (2), Ug(t) is a per unit value of the grid-connected point voltage of the wind turbine generator; i isNThe voltage is a per unit value of rated current of the grid-side converter;
step 4.3, because the static reactive compensation operation mode takes the reactive current as the main control object, the active current reference value needs to be limited, and the active current reference value limiting value i is obtained by using the formula (3)dref2(t)：
Figure BDA0001518163920000091
In the formula (3), imaxThe maximum current value allowed by the grid-side converter;
step 4.4, selecting the instruction value i of the direct current voltage PI regulator in the steady statedref1And an active current reference value limit value idref2(t) the smaller value is used as the active current reference value i of the grid-side converter when the voltage of the power grid dropsdref(t)；
When the instruction value i of the DC voltage PI regulatordref1Less than the active current reference value limit value idref2(t) when the active power reference value i of the grid-side converter is obtained, the outer ring of the direct-current side voltage of the grid-side converter can still adjust the direct-current side voltagedref(t) is still set to the command value i of the DC voltage PI regulatordref1(ii) a When the instruction value i of the DC voltage PI regulatordref1Greater than the active current reference value limit value idref2(t), the outer ring of the direct current side voltage cannot effectively keep the direct current side voltage stable, at the moment, a direct current side unloading circuit needs to be put into the outer ring of the direct current side voltage, redundant energy accumulated on the direct current side is consumed, the direct current side voltage is kept in a safe range, and an active power reference value i of the grid-side converterdref(t) setting to an active current reference value limit value idref2(t)。
Step 5, judging the conduction condition of each wind turbine generator unloading circuit;
step 5.1, during the three-phase short circuit fault period of the direct-drive wind power plant, the input wind energy of each wind power generation set is kept unchanged, and at t1At any moment, the wind energy P input by the wind turbine generator set at the failure foreground can be usedjwThe active power sent into the power grid by the wind turbine generator system grid side converter is obtained by using the formula (4) to obtain the direct current bus voltage U of the jth wind turbine generator system grid side power converterjdc(t1)：
Figure BDA0001518163920000092
In the formula (4), CjIs the DC bus capacitance value S of the jth wind turbine generator setBThe reference capacity of the grid-side converter;
step 5.2, compare t1DC bus voltage U of jth wind turbine generator set at momentjdc(t1) And unloading circuit action threshold Udc_inSize of (1), if Ujdc(t1) Greater than Udc_inIf so, judging that the jth wind turbine generator set unloading circuit is conducted; otherwise, judging that the unloading circuit of the jth wind turbine generator set is not conducted;
step 6, dividing the direct-drive wind power plant cluster by applying an immune random K value and a sensitive clustering center improved K-means clustering algorithm;
step 6.1, dividing all the wind turbine generators which are conducted by the unloading circuits in the n direct-drive permanent magnet wind turbine generators into a cluster;
and 6.2, during the three-phase short circuit fault of the direct-drive wind power plant, due to the fault isolation effect of the full-power converter, the terminal voltage value of each wind turbine generator is basically kept unchanged, the terminal voltage value of each wind turbine generator at the fault moment is used as a clustering index, an immune random K value and a sensitive clustering center are applied to improve a K mean value clustering algorithm, and m wind turbine generators which are not conducted by the rest unloading circuits are divided into K turbine groups.
FIG. 5 is a flow chart of an improved K-means clustering algorithm for immune random K-values and sensitive clustering centers, which mainly comprises the following steps:
step 6.2.1, the rest of m wind turbine generators with the unloading circuits not conductedTerminal voltage value { U1(t0),U2(t0),···,Up(t0),···,Uq(t0),···,Um(t0) As a sample data set, it is denoted as S ═ x1,x2,···,xp,···,xq,···,xm}; wherein, Up(t0) And Uq(t0) Respectively representing the terminal voltage values of the p-th and q-th wind turbine generators, taking the per unit value, and converting the U value into the voltage valuep(t0) And Uq(t0) Respectively denoted as data object xpAnd xq，p,q＝1,2,···,m,p≠q；
Step 6.2.2, calculate data object xpAnd xqIs (x) ofp,xq)；
The similarity between the data objects is measured by Euclidean distance, and the smaller the Euclidean distance between the data objects is, the more similar the data objects are;
step 6.2.3, obtaining the average distance MeanDist between any two data objects in the sample data set S by using the formula (5):
Figure BDA0001518163920000101
in the formula (5), the reaction mixture is,
Figure BDA0001518163920000102
the number of combinations of 2 data objects is selected from n data objects;
step 6.2.4, define: with data object xpThe area centered at the average distance MeanDist is called data object xpOf the data object x, the data object xpIs called data object xpA density parameter based on the distance MeanDist;
obtaining data object x using equation (6)pDensity parameter (x) ofpMeanDist) to obtain density parameters for m data objects:
Figure BDA0001518163920000103
in the formula (6), u (MeanDist-dist (x)p,xq) Represents a function and has:
Figure BDA0001518163920000104
data object xpThe greater the number of data objects in the neighborhood, i.e. data object xpDensity parameter (x) ofpMeanDist) larger, the data object x is illustrated as beingpThe better the clustering effect as the clustering center.
Step 6.2.5, adding the density parameters of the first M data objects with the largest density parameter in the density parameters of the M data objects into the candidate point set D, where D is { density (x)pMeanDist), p ═ 1,2, ·, M }; typically, M takes the value of M/2;
a cluster center is selected from M candidate data objects with high density parameters, so that the cluster center is guaranteed to be a data object in a relatively dense area in the same class, and high similarity in the class is guaranteed.
Step 6.2.6, defining a loop variable r, and initializing r to be 1; defining the clustering number as k, and initializing k to be 1; defining the maximum similarity between the initial average classes as AMSr-1And initializing the AMSr-1(ii) a Defining the cluster center set as Ar-1And initialized to an empty set; taking the candidate point set D as the r-1 st candidate point set Dr-1；
Step 6.2.7, from the r-1 st candidate point set Dr-1Selecting the mean value of two data objects with the maximum density parameter as the r-th clustering center crPut into clustering center set Ar-1In the method, the r-th clustering center set A is obtainedr(ii) a Simultaneously, the density parameters corresponding to the two selected data objects are selected from the r-1 st candidate point set Dr-1Deleting to obtain the r-th candidate point set Dr；
Step 6.2.8, from the r-th candidate point set DrMiddle selection and r-th clustering centerSet ArThe data object with the farthest cluster center distance in (1) is used as the (r + 1) th cluster center cr+1Put into the r-th clustering center set ArIn the method, the r +1 th clustering center set A is obtainedr+1(ii) a At the same time, cluster the r +1 th cluster center cr+1Corresponding density parameter is selected from the r-th candidate point set DrDeleting to obtain the r +1 th alternative point set Dr+1；
Step 6.2.9, assigning k +1 to k;
step 6.2.10, calculating the rest data objects in the sample data set S and the r +1 th clustering center set A respectivelyr+1The Euclidean distance of each clustering center is determined, and each data object is distributed to the class where the clustering center with the nearest Euclidean distance is located, so that k classes are obtained;
step 6.2.11, obtaining the data object x of any ith class in the k classes by using the formula (8)pClustering center c to class iiMean value s of the distance betweeni：
Figure BDA0001518163920000111
In the formula (8), PiThe total number of the data objects in the ith class; 1,2, ·, k;
step 6.2.12, obtaining the clustering center c of the ith class by using the formula (9)iCluster center c with class jjA distance d betweeni,j：
di,j＝dist(ci,cj) (9)
In formula (9), i, j is 1,2, L, k, i ≠ j;
step 6.2.13, obtaining the maximum similarity AMS between the r-th average class by using the formula (10)r：
Figure BDA0001518163920000121
The average inter-class maximum similarity AMS represents the average of the maximum similarities among the classes, when the AMS obtains the minimum value, the clustering effect is the best at the moment, and the k value is the best clustering number at the moment.
Step 6.2.14, determine AMSr＜AMSr-1If true, go to step 6.2.15; otherwise, go to step 6.2.19;
step 6.2.15, obtaining updated cluster center c 'by using formula (11)'i：
Figure BDA0001518163920000122
In the formula (11), PiIs the total number of data objects in class i, xpIs any data object of the ith class, i ═ 1,2, ·, k;
step 6.2.16, calculating the r +1 th candidate point set Dr+1The sum of Euclidean distances between the data object corresponding to any one density parameter and all updated clustering centers is obtained, so that a set of the sums of the Euclidean distances between the data object corresponding to all density parameters and all updated clustering centers is obtained, and the data object corresponding to the maximum value is selected from the set of the sums of the Euclidean distances to serve as the (r + 2) th clustering center cr+2And put into the r +1 th clustering center set Ar+1To obtain the r +2 th clustering center set Ar+2；
Selecting a cluster center c 'corresponding to all updates from M candidate data objects of the high-density parameter'iThe data object with the largest sum of Euclidean distances is used as a new clustering center to be added into the current clustering center set, so that the clustering centers of different classes are mutually exclusive as much as possible, and the low similarity among the different classes is ensured.
Step 6.2.17, clustering the r +2 th clustering center cr+2The corresponding density parameter is selected from the r +1 th candidate point set Dr+1Deleting to obtain the r +2 th alternative point set Dr+2；
6.2.18, assigning r +1 to r, and turning to step 6.2.9;
step 6.2.19, use AMSr-1The corresponding k clustering centers are used as initial clustering centers of a kmeans clustering algorithm, and the samples are subjected to the clusteringAnd (5) carrying out a kmeans clustering algorithm on the data set S to obtain k machine groups.

Claims (1)

1. A direct-drive wind power plant cluster division method considering low-voltage ride through characteristics is characterized by comprising the following steps:
step 1, establishing a detailed model of a direct-drive wind power plant;
the direct-drive wind power plant is composed of n direct-drive permanent magnet wind power units with the same type number, and the detailed model comprises the following components: the system comprises a single machine model, a current collection circuit model, a generator-end transformer model and a main transformer model of each wind turbine generator set in a wind power plant; wherein, the stand-alone model of wind turbine generator system includes: the system comprises a wind turbine model, a permanent magnet synchronous generator model, a full-power converter and a control system model thereof, a variable pitch control system model, an unloading circuit and a control system model thereof;
step 2, collecting input wind energy of each wind turbine generator set;
collecting the real-time wind speeds of n direct-drive permanent magnet wind turbine generators of the direct-drive wind power plant, and recording the input wind energy of each wind turbine generator as { P }1w,P2w,···,Pjw,···,Pnw}，PjwRepresenting the input wind energy of the jth wind power generation set;
step 3, collecting terminal voltage values of each wind turbine generator at the moment of failure;
set at t0At the moment, a three-phase short-circuit fault is arranged at an outlet of the direct-drive wind power plant, and t is acquired0The terminal voltage value of each wind turbine generator set at any moment is marked as { U1(t0),U2(t0),···,Uj(t0),···,Un(t0)}，Uj(t0) Represents t0The voltage value of the generator end of the jth wind turbine generator set at the moment is taken as a per unit value; at t1The three-phase short-circuit fault is removed at any moment;
step 4, setting an active current reference value of the grid-side converter when the voltage of the power grid drops;
step 4.1, obtaining the instruction value i of the grid-side converter direct-current voltage PI regulator in the steady state by using the formula (1)dref1：
Figure FDA0002409030720000011
In the formula (1), KdpAnd KdiProportional coefficients and integral coefficients of a direct current bus voltage control outer ring in the grid-side converter are respectively set; u shapedcTaking a per unit value for the direct-current bus voltage of the wind turbine generator in a steady state;
Figure FDA0002409030720000012
taking a per unit value as a direct current bus voltage reference value of the wind turbine generator;
step 4.2, during the three-phase short circuit fault of the direct-drive wind power plant, when the voltage of the grid-connected point of the wind turbine generator is lower than 20% of the nominal voltage, the wind turbine generator is cut out of the power grid;
when the grid-connected point voltage of the wind turbine generator is within 20% -90% of the nominal voltage, calculating a reactive current reference value i of the grid-side converter by using a formula (2)qref2(t)：
iqref2(t)＝1.5(0.9-Ug(t))IN，0.2pu≤Ug(t)≤0.9pu (2)
In the formula (2), Ug(t) is a per unit value of the grid-connected point voltage of the wind turbine generator; i isNThe voltage is a per unit value of rated current of the grid-side converter;
step 4.3, obtaining an active current reference value limiting value i by using the formula (3)dref2(t)：
Figure FDA0002409030720000021
In the formula (3), imaxThe maximum current value allowed by the grid-side converter;
step 4.4, selecting the instruction value i of the direct current voltage PI regulator in the steady statedref1And an active current reference value limit value idref2(t) the smaller value is used as the active current reference value i of the grid-side converter when the voltage of the power grid dropsdref(t)；
Step 5, judging the conduction condition of each wind turbine generator unloading circuit;
step 5.1, at the t1At any moment, the direct current bus voltage U of the jth wind turbine generator set network side power converter is obtained by using the formula (4)jdc(t1)：
Figure FDA0002409030720000022
In the formula (4), CjIs the DC bus capacitance value S of the jth wind turbine generator setBThe reference capacity of the grid-side converter;
step 5.2, compare t1DC bus voltage U of jth wind turbine generator set at momentjdc(t1) And unloading circuit action threshold Udc_inSize of (1), if Ujdc(t1) Greater than Udc_inIf so, judging that the jth wind turbine generator set unloading circuit is conducted; otherwise, judging that the unloading circuit of the jth wind turbine generator set is not conducted;
step 6, dividing the direct-drive wind power plant cluster by applying an immune random K value and a sensitive clustering center improved K-means clustering algorithm;
step 6.1, dividing all the wind turbine generators which are conducted by the unloading circuits in the n direct-drive permanent magnet wind turbine generators into a cluster;
step 6.2, dividing m wind turbine generators with the rest of the unloading circuits not conducted into K groups by taking the terminal voltage values of the wind turbine generators at the fault moment as grouping indexes and applying an immune random K value and a sensitive clustering center to improve a K-means clustering algorithm;
step 6.2.1, the terminal voltage values { U } of the m wind turbine generators with the residual unloading circuits not conducted1(t0),U2(t0),···,Up(t0),···,Uq(t0),···,Um(t0) As a sample data set, it is denoted as S ═ x1,x2,···,xp,···,xq,···,xm}; wherein, Up(t0) And Uq(t0) Respectively representing the terminal voltage values of the p-th and q-th wind turbine generators, taking the per unit value, and converting the U value into the voltage valuep(t0) And Uq(t0) Respectively denoted as data object xpAnd xq，p,q＝1,2,···,m,p≠q；
Step 6.2.2, calculate data object xpAnd xqIs (x) ofp,xq)；
Step 6.2.3, obtaining the average distance MeanDist between any two data objects in the sample data set S by using the formula (5):
Figure FDA0002409030720000031
in the formula (5), the reaction mixture is,
Figure FDA0002409030720000032
the number of combinations of 2 data objects is selected from n data objects;
step 6.2.4, define: with data object xpThe area centered at the average distance MeanDist is called data object xpOf the data object x, the data object xpIs called data object xpA density parameter based on the distance MeanDist;
obtaining data object x using equation (6)pDensity parameter (x) ofpMeanDist) to obtain density parameters for m data objects:
Figure FDA0002409030720000033
in the formula (6), u (MeanDist-dist (x)p,xq) Represents a function and has:
Figure FDA0002409030720000034
step 6.2.5, adding the density parameters of the first M data objects with the largest density parameter in the density parameters of the M data objects into the candidate point set D, where D is { density (x)p,MeanDist),p＝1,2,···,M}；
Step 6.2.6, defining a loop variable r, and initializing r to be 1; defining the clustering number as k, and initializing k to be 1; defining the maximum similarity between the initial average classes as AMSr-1And initializing the AMSr-1(ii) a Defining the cluster center set as Ar-1And initialized to an empty set; taking the candidate point set D as the r-1 st candidate point set Dr-1；
Step 6.2.7, from the r-1 st candidate point set Dr-1Selecting the mean value of two data objects with the maximum density parameter as the r-th clustering center crPut into clustering center set Ar-1In the method, the r-th clustering center set A is obtainedr(ii) a Simultaneously, the density parameters corresponding to the two selected data objects are selected from the r-1 st candidate point set Dr-1Deleting to obtain the r-th candidate point set Dr；
Step 6.2.8, from the r-th candidate point set DrThe (r) th cluster center set ArThe data object with the farthest cluster center distance in (1) is used as the (r + 1) th cluster center cr+1Put into the r-th clustering center set ArIn the method, the r +1 th clustering center set A is obtainedr+1(ii) a At the same time, cluster the r +1 th cluster center cr+1Corresponding density parameter is selected from the r-th candidate point set DrDeleting to obtain the r +1 th alternative point set Dr+1；
Step 6.2.9, assigning k +1 to k;
step 6.2.10, calculating the rest data objects in the sample data set S and the r +1 th clustering center set A respectivelyr+1The Euclidean distance of each clustering center is determined, and each data object is distributed to the class where the clustering center with the nearest Euclidean distance is located, so that k classes are obtained;
step 6.2.11, obtaining the data object x of any ith class in the k classes by using the formula (8)pClustering center c to class iiMean value s of the distance betweeni：
Figure FDA0002409030720000041
In the formula (8), PiThe total number of the data objects in the ith class; 1,2, ·, k;
step 6.2.12, obtaining the clustering center c of the ith class by using the formula (9)iCluster center c with class jjA distance d betweeni,j：
di,j＝dist(ci,cj) (9)
In formula (9), i, j ≠ 1,2, ·, k, i ≠ j;
step 6.2.13, obtaining the maximum similarity AMS between the r-th average class by using the formula (10)r：
Figure FDA0002409030720000042
Step 6.2.14, determine AMSr＜AMSr-1If true, go to step 6.2.15; otherwise, go to step 6.2.19;
step 6.2.15, obtaining updated cluster center c 'by using formula (11)'i：
Figure FDA0002409030720000043
In the formula (11), PiIs the total number of data objects in class i, xpIs any data object of the ith class, i ═ 1,2, ·, k;
step 6.2.16, calculating the r +1 th candidate point set Dr+1The sum of Euclidean distances between the data object corresponding to any one density parameter and all updated clustering centers is obtained, so that a set of the sums of the Euclidean distances between the data object corresponding to all density parameters and all updated clustering centers is obtained, and the data object corresponding to the maximum value is selected from the set of the sums of the Euclidean distances to serve as the (r + 2) th clustering center cr+2And put into the r +1 th clustering center set Ar+1To obtain the r +2 th clustering center set Ar+2；
Step 6.2.17, clustering the r +2 th clustering center cr+2Is correspondingly provided withFrom the r +1 th order candidate point set Dr+1Deleting to obtain the r +2 th alternative point set Dr+2；
6.2.18, assigning r +1 to r, and turning to step 6.2.9;
step 6.2.19, use AMSr-1And performing the K mean value clustering algorithm on the sample data set S to obtain K machine groups.
CN201711394587.9A 2017-12-21 2017-12-21 Direct-drive wind power plant cluster division method considering low-voltage ride through characteristics Expired - Fee Related CN108054753B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711394587.9A CN108054753B (en) 2017-12-21 2017-12-21 Direct-drive wind power plant cluster division method considering low-voltage ride through characteristics

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711394587.9A CN108054753B (en) 2017-12-21 2017-12-21 Direct-drive wind power plant cluster division method considering low-voltage ride through characteristics

Publications (2)

Publication Number Publication Date
CN108054753A CN108054753A (en) 2018-05-18
CN108054753B true CN108054753B (en) 2020-05-22

Family

ID=62131163

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711394587.9A Expired - Fee Related CN108054753B (en) 2017-12-21 2017-12-21 Direct-drive wind power plant cluster division method considering low-voltage ride through characteristics

Country Status (1)

Country Link
CN (1) CN108054753B (en)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110704995B (en) * 2019-11-28 2020-05-01 电子科技大学中山学院 Cable layout method and computer storage medium for multiple types of fans of multi-substation
CN111367254B (en) * 2020-02-26 2021-05-07 哈尔滨工业大学 Photovoltaic power station analytic single machine equivalence method, system and equipment
CN113346483B (en) * 2021-05-20 2022-05-27 华中科技大学 Low-voltage ride-through operation control method and system of power electronic transformer

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105243230A (en) * 2015-11-02 2016-01-13 合肥工业大学 Equivalent modeling method for direct-drive wind farm of considering off-load circuit

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105243230A (en) * 2015-11-02 2016-01-13 合肥工业大学 Equivalent modeling method for direct-drive wind farm of considering off-load circuit

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于免疫离群数据和敏感初始中心的K-means算法的风电场机群划分;潘险险等;《中国电机工程学报》;20161020;第36卷(第20期);第5461-5468页 *

Also Published As

Publication number Publication date
CN108054753A (en) 2018-05-18

Similar Documents

Publication Publication Date Title
Zou et al. A fuzzy clustering algorithm-based dynamic equivalent modeling method for wind farm with DFIG
CN111384725A (en) Short-circuit current calculation method and system for new energy power station
CN108054753B (en) Direct-drive wind power plant cluster division method considering low-voltage ride through characteristics
CN110970905B (en) Multi-power-supply reactive power coordination method and system for optimizing voltage control capability of wind power plant
CN114726009B (en) Wind power plant group reactive power hierarchical optimization control method and system considering power prediction
CN108879769B (en) Combined planning method for offshore wind farm grid connection and onshore power grid extension
CN111027179B (en) Equivalent modeling method for double-fed wind power plant considering auxiliary frequency modulation service
CN114123344A (en) Power system inertia evaluation method and device based on adaptive recursive least squares
CN109412203A (en) A kind of more feed-in electric system grid connection capacity Limit Calculation Methods of wind-powered electricity generation
Shao et al. Aggregation modeling algorithm of PMSM wind farm considering submarine cable power loss
CN114513004A (en) New energy station equivalence method based on improved k-means algorithm and application
CN107591843B (en) Double-fed wind field reactive power output optimization method in system recovery process
CN106410862A (en) Wind power plant single machine equivalent method based on active recovery slope correction
CN110571794B (en) Transient model equivalent calculation method suitable for doubly-fed wind power plant
CN113346484A (en) Power distribution network elasticity improving method and system considering transient uncertainty
CN114865649B (en) Wind-solar-storage integrated station reactive power regulation method and device and electronic equipment
CN116488267A (en) Modeling-based wind farm reactive capacity limit simulation calculation method and device
Asadollah et al. Decentralized reactive power and voltage control of wind farms with type-4 generators
CN116341202A (en) Marine wind farm submarine cable transmission limit calculation method and device
CN115483715A (en) Virtual synchronous generator self-adaptive control method and system for centralized photovoltaic power station
Ma et al. Online clustering modeling of photovoltaic power plant with LVRT control function
CN114825425A (en) New energy acceptance capacity assessment method and device for voltage drop induced frequency safety
Zhu et al. A study of dynamic equivalence method for multiple wind farms in urban power grids
Cao et al. Equivalence method for wind farm based on clustering of output power time series data
Ding et al. Research on simplified modeling of large-scale wind farms based on equivalent transfer function and aggregate equivalent

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
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20200522