CN101034808A - Distributed computing method of the features of the power system - Google Patents

Abstract

The invention relates to the distributed computing method of the electrical power system characteristic value, belongs to the electrical power system distributional simulation technology area, its characteristic lies in, it contains the interconnection electrical network syncopate method taped with boundary district and distributed computing method of the electrical power system characteristic value. The distributed computing method of electrical power system characteristic value has includes the distributional solve generator equivalent conductance, the distributional solve boundary equation of compatibility, the distributional solve amount of dispersion of district mechanical torque, the distributional solve system characteristic value. In the computational process, only need the little exchange boundary pitch point condition quantity and so on data of the various districts and the boundary district, is suitable for the electrical power system distributed environment, has a better usability.

Description

The distributed computing method of electric power system characteristic value

Technical field

The distributed computing method of electric power system characteristic value is to belong to electric power system distributed simulation technology field.

Background technology

In extensive interconnected electric power system, waving relatively of rotor usually taken place between generator or the electric power generator group, and cause lasting vibration when lacking damping.The frequency range of vibration is generally between 0.2 ~ 2.5Hz, so be called low-frequency oscillation, or electromechanical oscillations.Low-frequency oscillation in recent years also happens occasionally in China, has had a strong impact on power delivery and safety and stability between electrical network.In this case, the research to low-frequency oscillation problem receives much concern.

The operation and the management of China Da Qu interconnected network has the characteristics of " differentiated control, hierarchical control, distribution process ", under this environment, utilize the electric power system distributed computing technology, calculation task by each control centre of composition decomposition, finish the incorporate low-frequency oscillation analytical work of the whole network, the pattern of vibrating between can better analyzed area.But the achievement in research of the distribution analysis method of relevant low-frequency oscillation is also fewer up to now.Some scholars are the research that purpose has been carried out a series of parallel algorithms at the small interference stability analysis field to improve computational speed, yet because the environment difference that parallel computation and wide-area distribution type calculate is very big, parallel algorithm can not be applied directly in distributed calculating of electric power system and the analysis.

This patent is based on the distributed nature value-based algorithm of self-excitation method, and this method can be calculated the relevant characteristic value of low frequency oscillation mode by each subregion, thereby provides important references for the low-frequency oscillation problem of studying extensive interconnected electric power system.This method only needs low volume datas such as each subregion and border subregion exchange boundary node state amount in computational process, be applicable to the distributed environment of electric power system, has better practicability.

Eigenvalue calculation method based on self-excitation method

Self-excitation method is the Eigenvalue Analysis method that is used to analyze the dynamo-electric pattern of electric power system small interference stability.Its basic thought is from the frequency-domain analysis theory of linear system.In electric power system, for generator wherein, the differential equation of describing its dynamic characteristic is:

$\left\{\begin{array}{c}\mathrm{\Δ}\stackrel{\·}{\mathrm{\ω}}=\frac{1}{2H}(\mathrm{\Δ}{T}_m-{\mathrm{\ΔT}}_e)\\ \mathrm{\Δ}\stackrel{\·}{\mathrm{\δ}}={\mathrm{\ω}}_0\mathrm{\Δ\ω}\end{array}\right.---\left(1\right)$

Wherein ω is a rotor velocity, ω ₀Be system's reference angle speed, δ is a rotor angle,

Be the derivative of angular speed deviation,

The derivative of angular deviation (this paper hereinafter Δ is meant the deviation of variable, adds some points and all represent the derivative of this variable in the variable top).H is the inertia constant of generator, T _mBe machine torque, T _eBe electromagnetic torque, Δ T _mBe the deviation of machine torque, Δ T _eBe the deviation of electromagnetic torque, and Δ T _e=K _SΔ δ+K _DΔ ω, K _SBe synchronous torque coefficient, K _DBe damping coefficient, K _SComprised of the influence of whole network to this generator.

Can try to achieve with Δ T with Laplace transform according to formula (1) _mBe input, Δ ω for the ssystem transfer function of output is

$G\left(s\right)=\frac{\mathrm{\Δ\ω}\left(s\right)}{{\mathrm{\ΔT}}_m\left(s\right)}=\frac{1}{2\mathrm{Hs}+K_D\left(s\right)+\frac{K_S\left(s\right)}{s}}---\left(2\right)$

S in the formula is the transformation factor of Laplace transform, and when the transfer function denominator was zero, the s value of obtaining was the characteristic value of corresponding system in frequency-domain analysis.The basic principle of self-excitation method is calculated relevant with the dynamo-electric pattern of this generator characteristic value by the limit of calculating formula (2) exactly.This limit can be passed through the Newton solution by iterative method, and iteration form is

$s_{k+1}=s_k-{\left[\frac{\mathrm{\Δ}{T}_m\left(s\right)}{4H_e}\right]}_{s=s_k}---\left(3\right)$

H wherein _eFor " equivalent inertia coeffeicent ", be defined as

$H_e=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{H}_i{\left|{\mathrm{\Δ\ω}}_i\right|}^2---\left(4\right)$

H in the formula _iWith Δ ω _iBe respectively the inertia constant of i platform generator and the variable quantity of rotor velocity, n is the number of generator in the system.When $\left|{\mathrm{\&Delta;T}}_m\left(s\right)\right|<{\mathrm{\&epsiv;}}_{{\mathrm{\&Delta;T}}_m}$ (ε wherein _{Δ Tm}Certain threshold value of position, desirable 10 ^-6), iteration convergence, the s of this moment _kThe limit of being asked, the i.e. characteristic value of system exactly.Just can carry out the low-frequency oscillation analysis after trying to achieve the system features value.

What traditional self-excitation method adopted is the centralized calculation mode, this mode is had relatively high expectations for computing capability, computational speed is relatively slow, and centralized self-excitation method each step in iterative process all needs by finding the solution the value that system-wide network equation solves the disturbance quantity Δ Tm of machine torque, the operation and the management of China's electric power system has the characteristics of " differentiated control; hierarchical control; distribution process ", each subdispatch center has the parameter and the dynamic data of electrical network in the compass of competency, and the higher level control centre has the parameter and the dynamic data of interregional interconnection network.The data that can share between the control centre at the same level are limited, so the difficulty of integral data is very big.Above reason has restricted the development of traditional self-excitation method.Under distributed environment, need improve traditional self-excitation method to adapt to the characteristics of electric power system actual motion and management.The following chapters and sections of this paper will be on the basis of the interconnected network cutting method of band edge circle subregion, basic principle from self-excitation method, the distributed computing method of distributed eigenvalue calculation method and left and right sides characteristic vector is proposed, to adapt to the actual demand that the distributed low-frequency oscillation of electric power system is analyzed.

Summary of the invention

The present invention is based on traditional self-excitation method and carry out the general step of eigenvalue calculation,, proposed the computational methods of distributed characteristic value in conjunction with the characteristics of China's electric power system hierarchical management.This method only needs low volume datas such as each subregion and border subregion exchange boundary node state amount, can the distributed characteristic value of finding the solution.Practical Calculation to the IEEE39 node system shows that this method has drawn the result consistent with conventional method, and number of communications is also less, is applicable to the low-frequency oscillation analysis of electric power system distributed environment, has better practicability.

Generally, calculate, can adopt to comprise the responsible coordination calculation services of calculating and the distributed computing environment (DCE) of the service of the subregion eigenvalue calculation in a plurality of control centre coordinated in order to realize extensive interconnected network distributed nature value.Coordinate to communicate by network between calculation services and the subregion calculation services, the exchange data necessary is finished integrated eigenvalue calculation.According to the actual conditions of China's electric power system, we consider the electric power system distributed environment as accompanying drawing 1.In this environment, comprise up and down two-level scheduler center, the electric power data network by wide area between higher level control centre and each subdispatch center is connected.Details and the real-time status that network in the zone is arranged grasped at the subdispatch center; The higher level control centre grasp have interregional interconnection and directly under the factory station details and real-time status.Each subdispatch center and higher level control centre all need to dispose the analytical calculation system that has realized algorithm of the present invention.The invention is characterized in that it has following implementation step (idiographic flow is as shown in Figure 2):

Step 1 initialization (cutting of interconnected network)

The cutting method of interconnected network is the basis of realizing Distributed Calculation.The Distributed Calculation of wide area interconnected network requires the topological connection relation from electric power system, according to the operation of system's reality and the situation of management system is carried out cutting by the higher level control centre.This paper adopt document (Chen Ying. distributed computing method research [D] in the power grid. Beijing: Tsing-Hua University, the interconnected network cutting method that 2006.04.) proposes based on the band edge circle subregion of power-balance condition.In the system after the cutting, each subdispatch center has the parameter and the dynamic data (hereinafter also the subdispatch center being called subregion) of electrical network in the compass of competency, and the higher level control centre has the parameter and the dynamic data (hereinafter also the higher level control centre is called and coordinates side or border subregion) of interregional interconnection network.

With electric power system shown in Figure 3 is example, and system comprises two subregion A ₁, A ₂, they are connected by interconnection l.B ₁Represent the boundary node of subregion 1, the set of remembering them is A ₁ ^BThe network that other nodes in the subregion except that the node of border are formed is called internal node, and the set of remembering them is A ₁ ^InThe control centre of subregion 1 grasps A ₁ ^BAnd A ₁ ^InDetailed data.Subregion 2 is similar with it.During cutting, with the dummy node at interconnection l and two ends

(the corresponding B of difference ₁, B ₂) regard a subregion separately as, be called " border subregion ".

After using above-mentioned cutting method, the boundary node of i subregion satisfies following relation

$\left\{\begin{array}{c}{u}_x^{B_i}+{\mathrm{ju}}_y^{B_i}=u_x^{{\tilde{B}}_i}+{\mathrm{ju}}_y^{{\tilde{B}}_i}\\ (i_x^{B_i}+{\mathrm{ji}}_y^{B_i})+(i_x^{{\tilde{B}}_i}+{\mathrm{ji}}_y^{{\tilde{B}}_i})=0\end{array}\right.---\left(5\right)$

U wherein _x, u _y, i _x, i _yBe respectively the component of voltage and current on x, y reference axis under the synchronous coordinate system, i=1,2.

The generator of the selected research of step 2 is obtained the characteristic value of its frequency of oscillation between 0.2 ~ 2.5Hz

This characteristic value refers to the characteristic value of electric power system lienarized equation state matrix, the hereinafter unified characteristic value that is called for short.The key of self-excitation method is disturbance quantity (the Δ T that calculates machine torque in each step iterative process _m) value.After using " the interconnected network cutting method of band edge circle subregion " division network, use distributed self-excitation method to calculate little interference characteristic value.Fig. 4 is the flow chart that calculates the electrical network characteristic value with distributed method.

The disturbance quantity Δ T of step 201 Distributed Calculation machinery torque _m

In each step iterative process, calculate the disturbance quantity Δ T of machine torque _mValue be the key (formula (3)) of self-excitation method, this paper calculates Δ T according to detailed process shown in Figure 5 _m

First step initialization partition information

The computation sequence in this step as shown in Figure 6.Suppose that above-mentioned flow process is by subregion 1 initiation.Simultaneously, according to the situation of electric power system actual motion and management, we can think that the self-excitation machine of selection was also in this subregion when subregion 1 was initiated the calculating of self-excitation method.Subregion 1 is selected to be transmitted to each subregion by the border subregion after the s initial value.In addition, because electric power system is a time-varying system, so when subregion 1 initiates to calculate, also need the passing time reference coordinate.After each subregion was received this time coordinate, then unified network admittance matrix and the linearisation state matrix of choosing the own subregion in this moment participated in calculating.Initialized flow process as shown in Figure 6, the dotted arrow among the figure is represented the flow direction of by stages data.

The second step subregion 1 need be according to the characteristic value s of current iteration _kCalculate the equivalent current deviation amount Δ I of k platform generator _Ek(s)

Δ I _e(s) the equivalent electric current correction of expression, k is a generator numbering in the subregion 1, the method that adopts Prabha Kundur to mention at " PowerSystem Stability and Control " book is calculated equivalent electric current correction

${\mathrm{\&Delta;I}}_{e_k}\left(s_k\right)=c_1+\frac{c_2{\mathrm{\&omega;}}_0}{s_k}+C_r{(s_{k}I-A_{\mathrm{rr}})}^{-1}(a_{r1}+\frac{a_{r2}{\mathrm{\&omega;}}_0}{s_k})$

(variable of representing with bold signs in this paper formula is a vector, and what represent with ordinary symbol then is scalar)

A in the formula ₁₁, a ₁₂, a _1r, ω ₀, a _R1, a _R2, A _RrBe state variable in the formula (6) hereinafter $\left[\begin{array}{c}{\mathrm{\&Delta;\&omega;}}_k\\ {\mathrm{\&Delta;\&delta;}}_k\\ {\mathrm{\&Delta;x}}_k\end{array}\right]$ Element, the b of coefficient matrix ₁, B _rBe the coefficient of voltage variety Δ v, these coefficients are datum under the situation that network parameter is determined.

The 3rd each subregion of step also needs to calculate respectively according to the characteristic value sk of current iteration the equivalent admittance of each generator in this subregion

The equivalent admittance battle array of k platform generator (self-excitation machine) is:

$Y_{e_k}\left(s\right)=Y_{D_k}-C_r{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}{B}_r$

I platform generator (i ≠ k, i=1,2 ... h, h are total platform number of system's generator) equivalent admittance battle array be:

$Y_{e_i}\left(s\right)=-C_i{(\mathrm{sI}-A_i)}^{-1}{B}_i+Y_{D_i}$

I is and A in the formula _RrThe unit matrix that dimension equates, A _i, C _iBe Δ x in the formula (11) _iCoefficient matrix, B _i, Y _DiCoefficient matrix for Δ v in the formula (11).

The distributed border coordination equation of finding the solution of the 4th step

As step composition decomposition calculating busbar voltage vector, the border bus current vector of Fig. 5, marked the process of this part iteration among the figure with dash area according to the front.Wherein the computing formula of correction is:

For subregion i:

${\mathrm{\&Delta;v}}_i^{\mathrm{In}}={\left(Y_i^{\mathrm{InIn}}\right)}^{-1}({\mathrm{\&Delta;I}}_e\left(s\right){|}_i^{\mathrm{In}}-Y_i^{\mathrm{InB}}\&CenterDot;{\mathrm{\&Delta;v}}_i^B)$

${\mathrm{\&Delta;I}}_e\left(s\right){|}_i^B=Y_i^{\mathrm{BIn}}{\mathrm{\&Delta;v}}_i^{\mathrm{In}}+Y_i^{\mathrm{BB}}{\mathrm{\&Delta;v}}_i^B$

For coordinating side:

ΔI _e(s)| ^B-Y ^BΔv ^B＝δ ^B

Δ v in the formula _iThe correction value of representing the voltage vector of i subregion, subscript In represents internal node, subscript B represents boundary node (hereinafter unified this identification principle of using), Δ I _e(s) | _i ^InThe equivalent electric current correction value of representing i all internal nodes of subregion, its element is combined by the electric current correction of each node in i the subregion, Δ I _e(s) | _i ^BThe equivalent electric current correction value of representing i all boundary nodes of subregion, Δ I _e(s) | ^BThe equivalent electric current correction value of all boundary nodes of side is coordinated in expression.δ ^BFor calculating Δ I _e(s) | ^B-Y ^BΔ v ^BDifference, this value is used for judging whether boundary voltage, electric current restrain.The admittance battle array of the admittance battle array of subregion i (being not counted in border subregion circuit) is:

M is a subregion interior nodes number, Y _{I, j}Be node i, j (i, j=1,2 ... m) transadmittance between, the row, column sequence number corresponding node sequence number of this admittance battle array rearranges this admittance battle array according to internal node, boundary node, can obtain:

$Y_i=\left[\begin{array}{cc}{Y}^{\mathrm{InIn}}& Y^{\mathrm{InB}}\\ Y^{\mathrm{BIn}}& Y^{\mathrm{BB}}\end{array}\right]$

Y _i ^InIn, Y _i ^BB, Y _i ^InBAnd Y _i ^BInThe matrix of representing matrix Y corresponding internal node, boundary node in subregion i and incidence matrices are between the two directly got Y respectively _iMiddle respective element promptly.Y ^BThe matrix of representing matrix Y in the subregion of border.

If ‖ δ ^B‖＞ε _δ, then calculate and do not restrain voltage, the current value of the node that need revise the boundary.Makeover process adopts JFNG function (document " A Jacobian-Free Newton-GMRES (m) Method with Adaptive Preconditioner andItsApplication for Power Flow Calculations; IEEE TRANSACTIONS ON POWER SYSTEMS; VOL.21; NO.3; AUGUST 2006 " discloses this function), imports all partition boundaries voltage variety Δ v ^B, obtain all partition boundaries voltage variety Δ v ^BCorrection, can be expressed as: Δ v ^B=Δ v ^B+ Δ (Δ v ^B), Δ (Δ v wherein ^B)=JFNG (Δ v ^B), recomputated for the 4th step then.

If ‖ δ ^B‖＜ε _δ, then calculate convergence.Program changes following the 5th step over to.

It is original that above formula belongs to this paper, encloses its derivation below.

The complete inearized model equation of generator can be designated as

$\left[\begin{array}{c}\mathrm{\&Delta;}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_k\\ \mathrm{\&Delta;}{\stackrel{\&CenterDot;}{\mathrm{\&delta;}}}_k\\ \mathrm{\&Delta;}{\stackrel{\&CenterDot;}{x}}_k\end{array}\right]\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{1r}\\ {\mathrm{\&omega;}}_0& 0& 0\\ a_{r1}& a_{r2}& A_{\mathrm{rr}}\end{array}\right]\left[\begin{array}{c}{\mathrm{\&Delta;\&omega;}}_k\\ {\mathrm{\&Delta;\&delta;}}_k\\ {\mathrm{\&Delta;x}}_k\end{array}\right]+\left[\begin{array}{c}{b}_1\\ 0\\ B_r\end{array}\right]\mathrm{\&Delta;v}+\left[\begin{array}{c}\frac{1}{2H_k}\\ 0\\ 0\end{array}\right]{\mathrm{\&Delta;T}}_m---\left(6\right)$

${\mathrm{\&Delta;i}}_k=\left[c_1{c}_2{C}_r\right]\left[\begin{array}{c}{\mathrm{\&Delta;\&omega;}}_k\\ {\mathrm{\&Delta;\&delta;}}_k\\ {\mathrm{\&Delta;x}}_k\end{array}\right]-Y_{D_k}\mathrm{\&Delta;v}---\left(7\right)$

Δ ω _k, Δ δ _kDefinition with above, Δ i _kBe generator machine end injection current correction, subscript k represents the label of this generator in system.Δ v is a network bus voltage vector correction, Δ x _kFor in the model equation except Δ ω _kWith Δ δ _kOutside the vector corrected amount (depend on used model, Δ xk is a vector of determining when model is determined) of all state variables.c ₁, c ₂, C _rFormula (7) state variable $\left[\begin{array}{c}{\mathrm{\&Delta;\&omega;}}_k\\ {\mathrm{\&Delta;\&delta;}}_k\\ {\mathrm{\&Delta;x}}_k\end{array}\right]$ Coefficient, Y _DkCoefficient for voltage correction amount v, these coefficients determine it is that (detail formula is referring to PRABHA KUNDUR.Power System Stability and Control.McGraw-HillCompanies for known quantity at network parameter, Inc., 1994. the 12nd chapter formula 12.8 joints), following formula is carried out Laplace transform and cancellation Δ δ _kWith Δ x _k, can obtain

${\mathrm{\&Delta;T}}_m=2H[s-a_{11}-\frac{a_{12}{\mathrm{\&omega;}}_0}{s}-a_{1r}{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}(a_{r1}+\frac{a_{r2}{\mathrm{\&omega;}}_0}{s})]{\mathrm{\&Delta;\&omega;}}_k-2H[b_1+a_{1r}{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}{B}_r]\mathrm{\&Delta;v}---\left(8\right)$

${\mathrm{\&Delta;i}}_k=\mathrm{\&Delta;}{I}_{e_k}\left(s\right)-Y_{e_k}\left(s\right)\mathrm{\&Delta;v}---\left(9\right)$

Wherein

${\mathrm{\&Delta;I}}_{e_k}\left(s\right)=c_1+\frac{c_2{\mathrm{\&omega;}}_0}{s}{C}_r{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}(a_{r1}+\frac{a_{r2}{\mathrm{\&omega;}}_0}{s})---\left(10\right)$

$Y_{e_k}\left(s\right)=Y_{D_k}-C_r{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}{B}_r$

The definition of each amount is the same in the formula.

For the generator of being studied, can establish Δ ω _k=1.0p.u. (p.u. represents perunit value).So, in formula (8), all the local variable of generator k except that Δ v, so composition decomposition is found the solution Δ v, just become realization Distributed Calculation Δ T _mKey.

For other generators on the network except generator k, its linearizing state equation also has and the similar form of formula (8), the Δ T in the formula _m=0.If brief note is following form

$\mathrm{\&Delta;}{\stackrel{\&CenterDot;}{x}}_i=A_i{\mathrm{\&Delta;x}}_i+B_i\mathrm{\&Delta;v}$

$\mathrm{\&Delta;}{i}_i=C_i{\mathrm{\&Delta;x}}_i+{Y_D}_i\mathrm{\&Delta;v}---\left(11\right)$

Wherein subscript i represents the label of this generator in system, i=1,2 ... h, i ≠ k, h are the total platform number of system's generator, Δ x _iBe the variable quantity of equipment state variable, Δ i _iBe the variable quantity of equipment to the injection current of network, Δ v is the variable quantity of network bus voltage vector, A _i, C _iBe Δ x _iCoefficient matrix, B _i, Y _DiCoefficient matrix for Δ v in the formula (11).Can obtain according to Laplace transform erased condition variable so

${\mathrm{\&Delta;i}}_i=[C_i{(\mathrm{sI}-A_i)}^{-1}{B}_i-Y_{D_i}]\mathrm{\&Delta;v}=-Y_{e_i}\left(s\right)\mathrm{\&Delta;v}---\left(12\right)$

Restriction relation below the voltage and current vector of the whole network satisfies:

ΔI＝Y _NΔv (13)

Wherein Δ I is the variable quantity of injection current vector, and Δ v is the variable quantity of busbar voltage vector, Y in the formula _NBe the total system network admittance matrix, form can be expressed as:

H is a total system node number, Y _{I, j}Be node i, j (i, j=1,2 ... h) transadmittance between.

Wushu (9) and formula (12) merge with formula (13), can get equivalent current:

ΔI _e(s)＝(Y _N+Y _e(s))Δv (14)

Δ I wherein _e(s) have only the position of generator k correspondence to have nonzero element Δ I in _Ek(s), Y _e(s) be Y _Ei(s) diagonal matrix, i=1,2 ..., h, when i=k, Y _Ei(s) determine by formula (9), otherwise determine by formula (12).

Following formula can be write the form of an accepted way of doing sth (15) difference functions so that iterative

g(Δv)＝ΔI _e(s)-(Y _N+Y _e(s))Δv (15)

Note Y=Y _N+ Y _e(s), according to the interconnected network cutting method of the band edge circle subregion of step 1, if system is divided into border subregion and n subregion, so for subregion i=1 ..., n has

$\left[\begin{array}{c}{g}_i^{\mathrm{In}}\left(\mathrm{\&Delta;v}\right)\\ g_i^B\left(\mathrm{\&Delta;v}\right)\end{array}\right]\left[\begin{array}{c}{\mathrm{\&Delta;I}}_{e_i}^{\mathrm{In}}\left(s\right)\\ {\mathrm{\&Delta;I}}_{e_i}^B\left(s\right)\end{array}\right]-\left[\begin{array}{cc}{Y}_i^{\mathrm{InIn}}& Y_i^{\mathrm{InB}}\\ Y_i^{\mathrm{BIn}}& Y_i^{\mathrm{BB}}\end{array}\right]\left[\begin{array}{c}{\mathrm{\&Delta;v}}_i^{\mathrm{In}}\\ {\mathrm{\&Delta;v}}_i^B\end{array}\right]=0---\left(6\right)$

To contain the amount of In be the variable of subregion internal node to subscript in the formula, and the amount that subscript contains B is the variable of partition boundaries node.For the border subregion, have

$g^{\tilde{B}}\left(\mathrm{\&Delta;v}\right)={\mathrm{\&Delta;I}}_e^{\tilde{B}}\left(s\right)-Y^B{\mathrm{\&Delta;v}}^{\tilde{B}}=0$

Subscript in the formula

Amount be physical quantity in the subregion of border, so The equivalent current correction value of expression border subregion,

The voltage correction value of expression border subregion.

Obtain easily from formula (16) and (17) ${\mathrm{\&Delta;v}}_i^{\mathrm{In}}={\left(Y_i^{\mathrm{InIn}}\right)}^{-1}({\mathrm{\&Delta;I}}_{e_i}^{\mathrm{In}}\left(s\right)-Y_i^{\mathrm{InB}}\&CenterDot;{\mathrm{\&Delta;v}}_i^B)$ And ${\mathrm{\&Delta;I}}_e^B\left(s\right)-Y^B{\mathrm{\&Delta;v}}^B=0$ (in the iterative process ${\mathrm{\&Delta;I}}_e^B\left(s\right)-Y^B{\mathrm{\&Delta;v}}^B={\mathrm{\&delta;}}^B,$ δ B is the difference of computational process).

The 5th step subregion 1 calculates Δ T _m

Subregion 1 calculates Δ T _m

${\mathrm{\&Delta;T}}_m=2H[s-a_{11}-\frac{a_{12}{\mathrm{\&omega;}}_0}{s}-a_{1r}{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}(a_{r1}+\frac{a_{r2}{\mathrm{\&omega;}}_0}{s})]-2H[b_1+a_{1r}{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}{B}_r]\mathrm{\&Delta;v}$

The explanation of each variable sees above in the formula.

Step 202 subregion 1 is judged Δ T _mConvergence

If | Δ T _m|＜10 ^-6, calculate convergence, then the internal layer iteration finishes, and the s that solve this moment is characteristic value to be asked, and eigenvalue calculation finishes;

If | Δ T _m|＞10 ^-6, calculate and do not restrain,, algorithm branches at step 203.

Step 203 is calculated the equivalent inertia constant H of the whole network _e

Computational process as shown in Figure 7, each subregion ${H_e}_i=\underset{j=1}{\overset{m_i}{\mathrm{\&Sigma;}}}{H}_j{\left|{\mathrm{\&Delta;\&omega;}}_j\right|}^2,$ The equivalent inertia constant of the whole network $H_e=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{H_e}_i$

H wherein _eBe equivalent inertia coeffeicent, H _EiBe the equivalent inertia coeffeicent of i subregion, H _jWith Δ ω _jBe respectively the inertia constant of j platform generator and the variable quantity of rotor velocity, n is the number of subregion.When $\left|{\mathrm{\&Delta;T}}_m\left(s\right)\right|<{{\mathrm{\&epsiv;}}_{\mathrm{\&Delta;T}}}_m$ (wherein ${{\mathrm{\&epsiv;}}_{\mathrm{\&Delta;T}}}_m>0$ For the convergence criterion of setting, can get ${{\mathrm{\&epsiv;}}_{\mathrm{\&Delta;T}}}_m={10}^{-6}$ ) time, iteration convergence, the s of this moment _kThe limit of being asked, the i.e. characteristic value of system exactly.

Step 204 subregion 1 computation of characteristic values s

Then Ci Shi characteristic value can by $s_{k+1}=s_k-{\left[\frac{\mathrm{\&Delta;}{T}_m\left(s\right)}{4H_e}\right]}_{s=s_k}$ Calculate.Program forwards step (201) to.

The definition of each amount is the same in the formula.

This section utilizes the IEEE39 node system to test above-mentioned distributed algorithm.The IEEE39 node system has comprised 39 buses, and 46 circuits, test macro are divided into borderline region and 3 zones, and the boundary node number that each zone comprises is as shown in table 1.

Table 1 test macro Area Node number

Table 1 Test systems and its partitions

Systematic name	Zone 1	Zone 2	The zone	Borderline region
					IEEE39	15	12	12	8

Said system is carried out distributed small interference stability eigenvalue calculation.Choosing parameter is: ε _{Δ Tm}=1e-6, ε _δ=1e-6.The generator of being studied is numbered 39, belongs to zone 1.Initial value s=5.0j.The self-excitation method convergence process contrast that the process of iteration convergence and conventional serial are calculated is as shown in table 2.

Table 2 iteration convergence process relatively

Table 2 Iteration comparison between serialized and distributed methods

Iterative step n	Traditional algorithm sn	Distributed algorithm
				External iteration sn	Border coordination equation iteration number of communications
		1 2 3 4 5 6	-0.0897+6.0744i -0.1182+4.6237i -0.1052+4.3468i -0.0983+4.3302i -0.0977+4.3299i -0.0977+4.3299i			-0.0897+6.0744i -0.1182+4.6237i -0.1052+4.3468i -0.0983+4.3302i -0.0977+4.3299i -0.0977+4.3299i	6 6 5 5 5 5

In distributed algorithm, the border coordination equation iteration of internal layer has been communicated by letter 6 times in the 1st step external iteration, subsequently, coordinating side calculates by using last preconditioning matrix, the iteration number of communications of the inner boundary equation of comptability reduces to some extent, only needs 5 communication when external iteration is near convergence at last.The data volume of being transmitted in the communication process is relevant with the scale of borderline region, and in this test macro, borderline region node number is 8, so only needs to transmit the data of 8 nodes in the communication process.

Under the electric power system distributed environment, network communication delay is the topmost performance bottleneck of Distributed Calculation.The distributed characteristic value algorithm that this paper proposed, owing in iterative process, need repeatedly swap data of each zone and borderline region, may be on the performance not as good as traditional serial algorithm.But, the test example in also as can be seen, the algorithm that this paper proposed is being only required under the condition of transmitting low volume datas such as boundary node state amount, still has convergence preferably, number of communications is also less, can satisfy practical power systems is finished calculating in a dispatching cycle requirement generally speaking.Because algorithm does not need the detailed network data in each regional exchange area in computational process, thereby has practical meaning under the management system of China's electric power system layering and zoning, can be used as the basis of online small interference stability parser.The characteristic value that characteristic value that the distributed algorithm convergence obtains and conventional serial algorithm obtain is consistent, and has proved the correctness of distributed algorithm.

Select generator that bus #31 connected as research object again, calculated characteristics vector sum correlation factor.Right characteristic vector can be calculated in finding the solution the process of characteristic value, does not need extra communication process.In this example, the border coordination equation has been communicated by letter 6 times in the left eigenvector computational process, has transmitted the quantity of state data of 8 boundary nodes; In the correlation factor computational process each the zone with the borderline region coordinating communication 2 times, transmitted 10 correlation factor results.Because under the electric power system distributed environment, network communication delay is the topmost performance bottleneck of Distributed Calculation, so above-mentioned communication process is the main factor that influences the computational efficiency of distributed algorithm.The signal intelligence of adding up from above as can be seen, the algorithm that this paper proposed is only required low volume datas such as transmitting the boundary node vector, only need a spot of communication can finish calculating, can satisfy practical power systems is finished calculating in a dispatching cycle requirement generally speaking.Because algorithm does not need the detailed network data in each regional exchange area in computational process, thereby has practical meaning under the management system of China's electric power system layering and zoning.

Description of drawings:

The electric power system distributed environment schematic diagram of Fig. 1 emulation.A higher level control centre of expression and three subdispatch centers among the figure link to each other by electric power networks between them.

Fig. 2 uses self-excitation method and carries out the flow chart that low-frequency oscillation is analyzed.The work of this paper is the distributed electric power system characteristic value of finding the solution, and this work is the core procedure of low-frequency oscillation distribution analysis method.

The partition method of Fig. 3 band edge circle subregion.The zone of both sides is two subregions among the figure, and middle zone is to coordinate side.For the system of reality, subregion is one often more than two and coordinate side all the time.

The overall flow of Fig. 4 distributed method computation of characteristic values.

Fig. 5 calculates Δ T _mFlow process.This figure is actual to be to calculate Δ T among Fig. 4 _mSpecific implementation.

Before Fig. 6 computation of characteristic values, the initialized flow process of each subregion and coordination side.

Fig. 7 calculates H _eFlow process.

Figure 86 node 3 partition systems are used for showing concrete execution mode.

Embodiment:

Example system is described:

With one 6 node 3 partition systems is example explanation calculation process.System comprises 2 generators, 1 load, 6 buses, 3 circuits and 3 transformer branch roads, as shown in Figure 8.System parameters (PSAT form):

The information of node:

Bus.con＝[...

1 1.0 1 0 1 1；

2 1.0 1 0 2 1；

3 1.0 1 0 3 1；

4 1.0 1 0 1 1；

5 1.0 1 0 2 1；

6 1.0 1 0 3 1；]；

Interconnection information:

Line.con＝[...

6 5 100 1.0 50 0 0 0.032 0.161

0.306 0 0 0 0 0；

5 4 100 1.0 50 0 0 0.01 0.085

0.176 0 0 0 0 0；

6 4 100 1.0 50 0 0 0.017 0.092

0.158 0 0 0 0 0；

2 5 100 1.0 50 0 1.025 0 0.0625

0 0 0 0 0 0；

3 6 100 1.0 50 0 1.025 0 0.0586

0 0 0 0 0 0；

1 4 100 1.0 50 0 1.025 0 0.0576

0 0 0 0 0 0]；

Slack bus information:

SW.con＝[...

1 100 1.0 1.04 0 99 -99 1.1 0.9

0.8 1]；

The Pv nodal information:

PV.con＝[...

2 100 1.0 1.63 1.025 99 -99 1.1 0.9

1；

3 100 1.0 0.85 1.025 99 -99 1.1 0.9

1]；

The Pq nodal information:

PQ.con＝[...

6 100 1.0 0.9 0.3 1.2 0.8 0；

5 100 1.0 1 0.35 1.2 0.8 0；

4 100 1.0 1.25 0.5 1.2 0.8 0]；

Synchronous motor information:

Syn.con＝[...

2 100 1.0 50 2 0.0521 0 0.8958 0.1198

0 6 0 0.8645 0.1969 0 0.535 0 12.8

0 0 0 1 1 0.002；

3 100 1.0 50 2 0.0742 0 1.3125 0.1813

0 5.89 0 1.2578 0.25 0 0.6 0 6.02

0 0 0 1 1 0.002；

]；

％D

Syn.con(1，19)＝10.0；

Syn.con(2，19)＝10.0；

Step 1:

System is divided into 3 subregions, is respectively subregion 1, subregion 2 and subregion 3.Subregion comprises generator 1, bus 2 and bus 5; Subregion 2 comprises generator 2, bus 3 and bus 6; Subregion 3 comprises load, bus 1 and bus 4.Boundary node with each subregion: bus 4, bus 5 and bus 6 are as the bus of border subregion.

In the calculating below, selected generator 2 belongs to subregion 2 as the self-excitation machine.

S=10.0i. is chosen in primary iteration

Step 201:

First step initialization partition information:

Subregion 1 admittance battle array:

Y _InIn：

017.361

-17.3610

Y _InB：

0-17.361

17.3610

Y _BIn：

0-17.361

17.3610

Y _BB：

017.361

-17.3610

Subregion 2 admittance battle arrays:

Y _InIn：

10.651+0.228i7.419-0.456i

-15.905+0.114i-10.651-0.228i

Y _InB：

0-16

160

Y _BIn：

0-16

160

Y _BB：

016

-160

Subregion 3 admittance battle arrays:

Y _InIn：

2.352 20.019

-20.422-2.352

Y _InB：

0-17.065

17.0650

Y _BIn：

0-17.065

17.0650

Y _BB：

017.065

-17.0650

The border subregion:

3.307379 21.94778 -1.36519 -11.6041 -1.94219 -10.5107

-21.9478 3.307379 11.6041 -1.36519 10.51068 -1.94219

-1.36519 -11.6041 2.552792 17.33823 -1.1876 -5.97513

11.6041 -1.36519 -17.3382 2.552792 5.975135 -1.1876

-1.94219 -10.5107 -1.1876 -5.97513 3.129796 16.25382

10.51068 -1.94219 5.975135 -1.1876 -16.2538 3.129796

Second step is according to the characteristic value s of current iteration _kCalculate equivalent electric current Δ I _Ek(s)

When s=10.0i, ${{\mathrm{\&Delta;I}}_e}_k\left(s\right)=\left[\begin{array}{c}-185.45i\\ -34.03i\end{array}\right]$

The 3rd each subregion of step also needs the characteristic value s according to current iteration _kCalculate the Y of each generator in this subregion respectively _Ej(s)

${Y_e}_1\left(s\right)=\left[\begin{array}{cc}11.0558-0.3401i& -9.3908+0.6808i\\ 0.2975-0.1699i& -11.0558+0.3401i\end{array}\right]$

${Y_e}_2\left(s\right)=\left[\begin{array}{cc}2.3517& 2.9539\\ -3.3568& -2.3517\end{array}\right]$

The distributed border coordination equation of finding the solution of the 4th step

Δ v=behind the iteration convergence

0.3231 -0.7833i

-0.2991 -8.4029i

0.3367 -7.1604i

-0.2974 +6.0561i

0.2933 +5.1534i

-0.2834 -22.5695i

0.3231 -0.7833i

-0.2991 -8.4029i

0.3296 -2.8327i

-0.3001 -2.4659i

0.3120 +1.0625i

-0.2920 -14.8989i

The 5th step was calculated Δ T _m

ΔT _m＝11.9908+8.2126i

Step 202 subregion 1 is judged Δ T _mConvergence

Current | Δ T _m|＞10 ^-6, do not restrain

Step 203 is calculated the inertia constant H of the whole network _e

Subregion 1: ${H_e}_1=10.0$

Subregion 2: ${H_e}_2=11.234$

Last total H _e=21.234

Step 204 subregion 1 computation of characteristic values s

s＝-0.9959+9.3179i

Above data are results of first step iteration, owing to calculate not convergence as yet, need to carry out repeatedly iteration below up to convergence, refuse record as space is limited.When iteration convergence, the final characteristic value that can provide, this characteristic value are that the referable analyst uses.

Claims (1)

1. The computational methods of 1, electric power system distributed nature value is characterized in that, this method contains following steps successively:

   The interconnected interconnected network of step (1) Da Qu based on the interconnected network cutting method of the band edge circle subregion of power-balance condition by following mode is the interconnected network cutting: a plurality of subdispatch centers that respectively call subregion (with two be example, called after A ₁And A ₂), and one be called the higher level control centre of coordinating side or border subregion, the parameter of electrical network and dynamic data in the subregion compass of competency, border subregion are administered the parameter and the dynamic data of the network that each interregional interconnection 1 formed;

   B ₁, B ₂Represent subregion A respectively ₁And A ₂Boundary node, A ₁ ^B, A ₂ ^BRepresent the set of boundary node separately respectively, other node in the subregion except that the node of border is called internal node, and corresponding internal node set is designated as A ₁ ^In, A ₂ ^InInterconnection 1 and two ends thereof are corresponded respectively to each partition boundaries Node B ₁, B ₂Dummy node

   

   

   Be considered as a subregion separately, be called the border subregion;

   The boundary node of described each subregion satisfies following relation:

   $\left\{\begin{array}{c}{u}_x^{B_i}+j{u}_y^{B_i}=u_x^{{\tilde{B}}_i}+{\mathrm{ju}}_y^{{\tilde{B}}_i}\\ (i_x^{B_i}+{\mathrm{ji}}_y^{B_i})+(i_x^{{\tilde{B}}_i}+{\mathrm{ji}}_y^{{\tilde{B}}_i})=0\end{array}\right.,$

   U wherein _x, u _y, i _x, i _yBe respectively the component of voltage and current on x, y reference axis under the synchronous coordinate system, i=1,2;

   Step (2) subregion A ₁Calculating according to the following steps: the generator of selected research, the frequency of oscillation of obtaining this generator be between 0.2 ~ 2.5Hz the time, the characteristic value s of the state matrix of electric power system lienarized equation;

   The deviation delta T of step (2.1) calculating machine torque _m,, its step is as follows:

   Each partition information of step (2.1.1) initialization:

   Step (2.1.1.1) subregion A ₁Inscribe oneself network admittance battle array and linearisation state matrix when determining t, obtain the initial estimate s of characteristic value s ₀

   Step (2.1.1.2) subregion A ₁Value s ₀Be transmitted to each subregion with time coordinate t by the border subregion;

   Each subregion of step (2.1.1.3) is after receiving this time coordinate, and unified network admittance battle array and the linearisation state matrix of choosing the own subregion of this moment t participates in following calculating:

   Step (2.1.2) subregion A ₁Characteristic value s with current iteration _kCalculate the equivalent current deviation amount Δ I of k platform generator in this subregion _Ek(s), step is as follows:

   Step (2.1.2.1) is designated as the complete inearized model equation of k platform generator:

   $\left[\begin{array}{c}\mathrm{\&Delta;}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_k\\ \mathrm{\&Delta;}{\stackrel{\&CenterDot;}{\mathrm{\&delta;}}}_k\\ \mathrm{\&Delta;}{\stackrel{\&CenterDot;}{x}}_k\end{array}\right]=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{1r}\\ {\mathrm{\&omega;}}_0& 0& 0\\ a_{r1}& a_{r2}& A_{\mathrm{rr}}\end{array}\right]\left[\begin{array}{c}\mathrm{\&Delta;}{\mathrm{\&omega;}}_k\\ \mathrm{\&Delta;}{\mathrm{\&delta;}}_k\\ \mathrm{\&Delta;}{x}_k\end{array}\right]+\left[\begin{array}{c}{b}_1\\ 0\\ B_r\end{array}\right]\mathrm{\&Delta;v}+\left[\begin{array}{c}\frac{1}{2H_k}\\ 0\\ 0\end{array}\right]\mathrm{\&Delta;}{T}_m$

   $\mathrm{\&Delta;}{i}_k=\left[\begin{array}{ccc}{c}_1& c_2& C_r\end{array}\right]\left[\begin{array}{c}{\mathrm{\&Delta;\&omega;}}_k\\ \mathrm{\&Delta;}{\mathrm{\&delta;}}_k\\ \mathrm{\&Delta;}{x}_k\end{array}\right]-Y_{D_k}\mathrm{\&Delta;v}$

   Wherein: Δ i _kBe generator machine end injection current vector departure, subscript k represents the label of this generator in system,

   Δ v is a network bus voltage vector departure,

   Δ x _kFor in the model equation except Δ ω _kWith Δ δ _kOutside the vectorial departure of all state variables, be set point,

   And the coefficient a of Δ v ₁₁, a ₁₂, a _1r, ω ₀, a _R1, a _R2, A _Rr, c ₁, c ₂, C _rAnd Y _Dk, b ₁, B _r

   When the interconnected network network equation is determined known quantity,

   H _kBe the inertia constant of k platform generator, be known quantity,

   ω ₀For generator reference angle speed, be set point,

   ω _kBe the rotor velocity of generator k, δ _kBe rotor angle,

   Δ ω _k, Δ δ _kBe respectively angular speed deviation and rotor angle deviation;

   Step (2.1.2.2) is calculated as follows Δ I _Ek(s _k):

   $\mathrm{\&Delta;}{I}_{e_k}\left(s_k\right)=c_1+\frac{c_2{\mathrm{\&omega;}}_0}{s_k}+C_r{(s_{k}I-A_{\mathrm{rr}})}^{-1}(a_{r1}+\frac{a_{r2}{\mathrm{\&omega;}}_0}{s_k}),$

   Wherein, I is and A _RrThe unit matrix that dimension equates;

   Step (2.1.2.3) comprises A ₁At each interior subregion according to current iteration characteristic value s _kCalculate the equivalent admittance of each generator in each subregion:

   Equivalent admittance as the k platform generator of self-excitation machine: $Y_{e_k}\left(s\right)=Y_{D_k}-C_r{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}{B}_r,$

   The equivalent admittance of i platform generator: $Y_{e_i}\left(s\right)=-C_i{(\mathrm{sI}-A_i)}^{-1}{B}_i+Y_{D_i},$

   I ≠ k, i=1,2 ... h, h are total platform number of system's generator;

   Step (2.1.3) comprises subregion A ₁Find the solution border convergence equation at each interior subregion, its step is as follows:

   Step (2.1.3.1) border subregion is set Value, pass through then $\mathrm{\&Delta;}{v}^{{\tilde{B}}_i}=\mathrm{\&Delta;}{v}^{B_i}$ Send to each subregion;

   Each subregion of step (2.1.3.2) calculates following each amount successively, obtains δ ^B:

   I subregion inside, the voltage deviation Δ v of node _i ^In:

   $\mathrm{\&Delta;}{v}_i^{\mathrm{In}}={\left(Y_i^{\mathrm{InIn}}\right)}^{-1}({\mathrm{\&Delta;I}}_e\left(s\right){|}_i^{\mathrm{In}}-Y_i^{\mathrm{InB}}\&CenterDot;\mathrm{\&Delta;}{v}_i^B),$

   The equivalent current deviation delta I of i partition boundaries node _e(s) | _i ^B:

   $\mathrm{\&Delta;}{I}_e\left(s\right){|}_i^B=Y_i^{\mathrm{BIn}}\mathrm{\&Delta;}{v}_i^{\mathrm{In}}+Y_i^{\mathrm{BB}}\mathrm{\&Delta;}{v}_i^B,$

   The equivalent current deviation delta I of the whole boundary nodes of each subregion _e(s) | ^B:

   ΔI _e(s)| ^B-Y ^BΔ _v ^B＝δ ^B，

   δ wherein ^BBeing used for judging whether boundary voltage, electric current restrain, is a border coordination parameter.Δ I _e(s) | ^BThe vector that obtains for the boundary node equivalent current that makes up each subregion;

   The admittance battle array of the admittance battle array of subregion i (being not counted in border subregion circuit) is:

   

   M is a subregion interior nodes number, Y _{I, j}Be node i, j (i, j=1,2 ... m) transadmittance between, the row, column sequence number corresponding node sequence number of this admittance battle array rearranges this admittance battle array according to internal node, boundary node, can obtain:

   $Y_i=\left[\begin{array}{cc}{Y}^{\mathrm{InIn}}& Y^{\mathrm{InB}}\\ Y^{\mathrm{BIn}}& Y^{\mathrm{BB}}\end{array}\right]$

   Y _i ^InIn, Y _i ^BB, Y _i ^InBAnd Y _i ^BInRepresent the matrix of corresponding internal node, boundary node among the subregion i and incidence matrices between the two respectively, directly get Y _iMiddle respective element promptly;

   Step (2.1.3.3) is judged ‖ δ ^BWhether ‖ is less than ε _δ, ε _δ=10 ^-6:

   If ‖ δ ^B‖＞ε _δ, execution in step (2.1.3.4) then,

   If ‖ δ ^B‖＜ε _δ, execution in step (2.1.4) then;

   Step (2.1.3.4) order

   Make Δ v ^B=Δ v ^B+ Δ (Δ v ^B), Δ (Δ v wherein ^B)=JFNG (Δ v ^B) (JFNG is disclosed function), repeated execution of steps (2.1.3.1)～step (2.1.3.3) is up to ‖ δ again ^B‖＜ε _δTill;

   Step (2.1.4) subregion A ₁Calculate the Δ T of k platform generator _m

   $\mathrm{\&Delta;}{T}_m=2H[s-a_{11}-\frac{a_{12}{\mathrm{\&omega;}}_0}{s}-a_{1r}{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}(a_{r1}+\frac{a_{r2}{\mathrm{\&omega;}}_0}{s})]-2H[b_1+a_{1r}{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}{B}_r]\mathrm{\&Delta;v}$ Step (2.2)

   Judge Δ T _mConvergence:

   If | Δ T _m|＜10 ^-6, calculate convergence, the s that solve this moment _kBe characteristic value to be asked,

   If | Δ T _m|＞10 ^-6, then do not restrain, change step (2.3);

   Step (2.3) is calculated as follows the equivalent inertia constant H of the whole network _e:

   $H_e=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{H}_{e_i},$

   H _EiBe the equivalent inertia constant of each subregion, $H_{e_i}=\underset{j=1}{\overset{m_i}{\mathrm{\&Sigma;}}}{H}_j{\left|\mathrm{\&Delta;}{\mathrm{\&omega;}}_j\right|}^2,$

   H _jWith Δ ω _jBe respectively the inertia constant of j platform generator and the variable quantity of rotor velocity, n is the number of subregion, and m is the total generating board number of system;

   Step (2.4) subregion A ₁Calculate the characteristic value that iteration is used:

   $s_{k+1}=s_k-{\left[\frac{\mathrm{\&Delta;}{T}_m\left(s\right)}{4H_e}\right]}_{s=s_k}$

   Recomputate step (2.1)～step (2.2) with this characteristic value.

Images