FAST AND FLEXIBLE HOLOMORPHIC EMBEDDING-BASED METHOD FOR ASSESSING POWER SYSTEM LOAD MARGINS

Abstract

A fast and flexible holomorphic embedding-based method for assessing power system load margins includes the following steps: S1, acquiring required decrypted power system data from a partner indirectly; S2, establishing a system of continuation power flow equations; S3, solving the continuation power flow equation by FFHE; and S4, designing and planning a scheduling policy for the power system based on the solved load margin. Compared with prediction through a linear function, the present method is considerably precise and efficient by utilizing rational approximants obtained by expanding arc-length series without repeatedly applying a local solver inefficiently and multifariously for correction. An efficient solution is developed to solve this type of nonlinear equations efficiently. Compared with existing methods, the present method is significantly improved in computational efficiency, computational accuracy, solvable system scale and the like.

Claims

1. A fast and flexible holomorphic embedding-based method for assessing power system load margins, comprising the following steps: step S1, acquiring required decrypted power system data from a partner indirectly; step S2, establishing a system of continuation power flow equations; step S3, solving the system of continuation power flow equations by FFHE; and step S4, designing and planning a scheduling policy for the power system based on the solved load margin.

2. The fast and flexible holomorphic embedding-based method for assessing power system load margins according to claim 1, wherein the establishing a system of continuation power flow equations in step S2 specifically comprises: a power flow problem is described as the following algebraic equations:
S.sub.i=V.sub.iΣ.sub.kY*.sub.i,kV*.sub.k,  (1) where S.sub.i=P.sub.i+jQ.sub.i represents a power injection at a bus, P.sub.i, Q.sub.i represent active and reactive powers, respectively, at bus i, Y.sub.i,k=G.sub.i,k+jB.sub.i,k is an admittance between corresponding buses, G.sub.i,k, B.sub.i,k represent a conductance and a susceptance, respectively, between corresponding buses, Y*.sub.i,k, V*.sub.k represent conjugates of Y.sub.i,k, V.sub.k, respectively, and V.sub.i represents a bus voltage; equation (1) can be converted to the rectangular form as:
P.sub.i=V.sub.i.sup.RΣ.sub.k{G.sub.i,kV.sub.k.sup.R−B.sub.i,kV.sub.k.sup.I}+V.sub.i.sup.IΣ.sub.k{B.sub.i,kV.sub.k.sup.R+G.sub.i,kV.sub.k.sup.I};
Q.sub.i=−V.sub.i.sup.RΣ.sub.k{B.sub.i,kV.sub.k.sup.R+G.sub.i,kV.sub.k.sup.I}+V.sub.l.sup.iΣ.sub.k{G.sub.i,kV.sub.k.sup.R−B.sub.i,kV.sub.k.sup.I};  (2) wherein V.sub.k.sup.R is a real part of the bus voltage, while V.sub.k.sup.I is an imaginary part of the bus voltage; and by equation (2), a continuation power flow equation can be constructed as
V.sub.i.sup.RΣ.sub.k{G.sub.i,kV.sub.k.sup.R−B.sub.i,kV.sub.k.sup.I}+V.sub.i.sup.IΣ.sub.k{B.sub.i,kV.sub.k.sup.R+G.sub.i,kV.sub.k.sup.I}−P.sub.i−λ.Math.(P.sub.targ,i−P.sub.base,i)=0;  (3)
−V.sub.i.sup.RΣ.sub.k{B.sub.i,kV.sub.k.sup.R+G.sub.i,kV.sub.k.sup.I}+V.sub.i.sup.IΣ.sub.k{G.sub.i,kV.sub.k.sup.R−B.sub.i,kV.sub.k.sup.I}−Q.sub.i−λ.Math.(Q.sub.targ,i−Q.sub.base,i)=0; wherein λ is a real number, P.sub.targ,i, P.sub.base,i represent a target value and a base value of the active power at a load bus, respectively, and Q.sub.targ,i, Q.sub.base,i represent a target value and a base value of the reactive power at the load bus, respectively.

3. The fast and flexible holomorphic embedding-based method for assessing power system load margins according to claim 1, wherein step S3 specifically comprises: step S301, initializing and setting an order of series expansion q.sub.max, a threshold of the acceptable mismatch e between two sides of the equation, and a given entry of at least one part in the initial point X.sub.0=(V.sup.R(0), V.sup.I(0), λ(0)); step S302, introducing a new parameter s, such as an arc-length parameter, and obtaining an embedding system; $\begin{array}{cc}{.Math.}_k\left\{{\left(\frac{{\mathrm{dV}}_k^R}{\mathrm{ds}}\right)}^2\left(s\right)+{\left(\frac{{\mathrm{dV}}_k^I}{\mathrm{ds}}\right)}^2\left(s\right)\right\}+{\left(\frac{d\lambda }{\mathrm{ds}}\left(s\right)\right)}^2=1;& \left(4\right)\end{array}$ $V_i^R\left(s\right){.Math.}_k\left\{G_{i,k}{V}_k^R\left(s\right)-B_{i,k}{V}_k^I\left(s\right)\right\}+V_i^I\left(s\right){.Math.}_k\left\{B_{i,k}{V}_k^R\left(s\right)+G_{i,k}^I{V}_k^I\left(s\right)\right\}-P_i-\lambda \left(s\right).Math.\left(P_{\mathrm{targ},i}-P_{\mathrm{base},i}\right)=0;$ $-V_i^R\left(s\right){.Math.}_k\left\{B_{i,k}{V}_k^R\left(s\right)+G_{i,k}{V}_k^I\left(s\right)\right\}+V_i^I\left(s\right){.Math.}_k\left\{G_{i,k}{V}_k^R\left(s\right)-B_{i,k}{V}_k^I\left(s\right)\right\}-Q_i-\lambda \left(s\right).Math.\left(Q_{\mathrm{targ},i}-Q_{\mathrm{base},i}\right)=0;$ step S303, representing λ(s), V.sup.R(s), V.sup.I(s) as series expansions in s;
λ(s)=Σ.sub.q≥0a.sub.0,qs.sup.q; V.sub.k.sup.R(s)=Σ.sub.q≥0a.sub.k,qs.sup.q, V.sub.k.sup.I(s)=Σ.sub.q≥0b.sub.k,qs.sup.q, k≥1; step S304, inserting the series expansion into the embedded system (4), and obtaining a system of equations by taking series coefficients in the series expansions as the unknowns:
Σ.sub.k≠ref{(Σ.sub.q≥0(1+q).Math.a.sub.k,q+1s.sup.q).sup.2+(Σ.sub.q≥0(1+q).Math.b.sub.k,q+1s.sup.q).sup.2}+(Σ.sub.q≥0(1+q).Math.a.sub.0,q+1s.sup.q)=1;
{Σ.sub.q≥0a.sub.i,qs.sup.q}.Math.Σ.sub.k{G.sub.i,kΣ.sub.q≥0a.sub.k,qs.sup.q−B.sub.i,kΣ.sub.q≥0b.sub.k,qs.sup.q}+{Σ.sub.q≥0b.sub.i,qs.sup.q}.Math.Σ.sub.k{B.sub.i,kΣ.sub.q≥0a.sub.k,qs.sup.q+G.sub.i,kΣ.sub.q≥0b.sub.k,qs.sup.q}−P.sub.i−{Σ.sub.q≥0a.sub.0,qs.sup.q}(P.sub.targ,i−P.sub.base,i)=0;
{(Σ.sub.q≥0a.sub.i,qs.sup.q}.Math.Σ.sub.k{B.sub.i,kΣ.sub.q≥0a.sub.k,qs.sup.q+G.sub.i,kΣ.sub.q≥0b.sub.k,qs.sup.q}+{Σ.sub.q≥0b.sub.i,qs.sup.q}.Math.Σ.sub.k{G.sub.i,kΣ.sub.q≥0a.sub.k,qs.sup.q−B.sub.i,kΣ.sub.q≥0b.sub.k,qs.sup.q}−Q.sub.i−{Σ.sub.q≥0a.sub.0,qs.sup.q}.Math.(Q.sub.targ,i−Q.sub.base,i)=0. step S305, comparing the coefficients of the terms with a same order of s, wherein when q=0,
a.sub.k,0=V.sub.k.sup.R(0), b.sub.k,0=V.sub.k.sup.I(0), 1≤k≤n; a.sub.0,0=λ(0)=λ.sub.0; a.sub.i,0.Math.Σ.sub.k{G.sub.i,ka.sub.k,0−B.sub.i,kb.sub.k,0}+b.sub.i,0.Math.Σ.sub.k{B.sub.i,ka.sub.k,0+G.sub.i,kb.sub.k,0}−P.sub.i−a.sub.0,0.Math.(P.sub.targ,i−P.sub.base,i)=0;
−a.sub.i,0.Math.Σ.sub.k{B.sub.i,ka.sub.k,0+G.sub.i,kb.sub.k,0}+b.sub.i,0.Math.Σ.sub.k{G.sub.i,ka.sub.k,0−B.sub.i,kb.sub.k,0}−Q.sub.i−a.sub.0,0.Math.(Q.sub.targ,i−Q.sub.base,i)=0. selecting at least one of a.sub.0,0, a.sub.k,0, b.sub.k,0 as a known initial value, to solve all a.sub.0,0, a.sub.k,0, b.sub.k,0; step S306, comparing the coefficients of the terms with a same order of s, wherein when q=1,
Σ.sub.k{a.sub.k,1.sup.2+b.sub.k,1.sup.2}+(a.sub.0,1).sup.2=1  (5)
a.sub.i,0.Math.Σ.sub.k{G.sub.i,ka.sub.k,1−B.sub.i,kb.sub.k,1}+a.sub.i,1.Math.Σ.sub.k(G.sub.i,ka.sub.k,0−B.sub.i,kb.sub.k,0)+b.sub.i,0.Math.Σ.sub.k{B.sub.i,k,a.sub.k,1+G.sub.i,kb.sub.k,1}+b.sub.i,1.Math.Σ.sub.k{B.sub.i,ka.sub.k,0+G.sub.i,kb.sub.k,0}−a.sub.0,1.Math.(P.sub.targ,i−P.sub.base,i)=0;  (6)
−a.sub.i,0.Math.Σ.sub.k(B.sub.i,ka.sub.k,1+G.sub.i,kb.sub.k,1)−a.sub.i,1.Math.Σ.sub.k{B.sub.i,ka.sub.k,0+G.sub.i,kb.sub.k,0}+b.sub.i,0.Math.Σ.sub.k{G.sub.i,ka.sub.k,1−B.sub.i,kb.sub.k,1}+b.sub.i,1.Math.Σ.sub.k(G.sub.i,ka.sub.k,0−B.sub.i,kb.sub.k,0)−a.sub.0,1.Math.(Q.sub.targ,i−Q.sub.base,i)=0;  (7) wherein a.sub.0,0, a.sub.k,0, b.sub.k,0 herein are known; in order to obtain solutions of the above-mentioned nonlinear equations (5)-(7), a novel solving method based on equivalent reduction of an equation and a quadratic formula of a quadratic equation can be utilized; step S307, deriving rational approximants using the series expansions of λ(s), V.sup.R(s), V.sup.I(s); step S308, substituting a value of the rational approximant at s=s.sub.0 into the equation (3), and making a comparison to figure out whether the mismatch between a left side and a right side of the equation is smaller than a preset acceptable threshold e; if yes, expanding so till the mismatch is not smaller than the preset acceptable threshold; if no, shrinking s.sub.0 till the mismatch is smaller than the preset acceptable threshold, i.e., returning s.sub.0 as great as possible with the mismatch smaller than the threshold e; step S309, taking (V.sup.R(s.sub.0), V.sup.I(s.sub.0), λ(s.sub.0)) as a new initial point X.sub.0; and step S310, repeating step S302 to step S309 till pinpointing an SNB point.

Description

BRIEF DESCRIPTION OF DRAWINGS

[0040] FIG. 1 is a processing flowchart of FFHE for solving a continuation power flow equation of the present invention.

[0041] FIG. 2 is a schematic diagram showing that a scheduling center provides power system load margins of the present invention.

DETAILED DESCRIPTION OF EMBODIMENTS

[0042] As shown in FIG. 1, the present invention provides a fast and flexible holomorphic embedding-based method for assessing power system load margins, mainly including the following steps:

[0043] Step S1, as the original power system data are confidential, the required decrypted power system data are acquired from a partner indirectly.

[0044] Step S2, a system of continuation power flow equations is established.

[0045] For brevity, a PQ bus is only discussed here, while a PV bus can be processed similarly.

[0046] A power flow problem is commonly described as the following algebraic equations:

S.sub.i=V.sub.iΣ.sub.kY*.sub.i,kV*.sub.k,  (1)

[0047] where S.sub.i=P.sub.i+jQ.sub.i represents a power injection at a bus, P.sub.i, Q.sub.i represent active and reactive powers, respectively, at bus i, Y.sub.i,k=G.sub.i,k+jB.sub.i,k is an admittance between corresponding buses, G.sub.i,k, B.sub.i,k represent a conductance and a susceptance, respectively, between corresponding buses, Y*.sub.i,k, V*.sub.k represent conjugates of Y.sub.i,k, V.sub.k, respectively, and V.sub.i represents a bus voltage; equation (1) can be converted to the rectangular form as:

P.sub.i=V.sub.i.sup.RΣ.sub.k{G.sub.i,kV.sub.k.sup.R−B.sub.i,kV.sub.k.sup.I}+V.sub.i.sup.IΣ.sub.k{B.sub.i,kV.sub.k.sup.R+G.sub.i,kV.sub.k.sup.I};

Q.sub.i=V.sub.i.sup.RΣ.sub.k{B.sub.i,kV.sub.k.sup.R+G.sub.i,kV.sub.k.sup.I}+V.sub.i.sup.IΣ.sub.k{G.sub.i,kV.sub.k.sup.R−B.sub.i,kV.sub.k.sup.I};  (2)

[0048] where V.sub.k.sup.R is a real part of the bus voltage, while V.sub.k.sup.I is an imaginary part of the bus voltage; and

[0049] by equation (2), a continuation power flow equation can be constructed as

V.sub.i.sup.RΣ.sub.k{G.sub.i,kV.sub.k.sup.R−B.sub.i,kV.sub.k.sup.I}+V.sub.i.sup.IΣ.sub.k{B.sub.i,kV.sub.k.sup.R+G.sub.i,kV.sub.k.sup.I}−P.sub.i−λ.Math.(P.sub.targ,i−P.sub.base,i)=0;

−V.sub.i.sup.RΣ.sub.k{B.sub.i,kV.sub.k.sup.R+G.sub.i,kV.sub.k.sup.I}+V.sub.l.sup.iΣ.sub.k{G.sub.i,kV.sub.k.sup.R−B.sub.i,kV.sub.k.sup.I}−Q.sub.i−λ.Math.(Q.sub.targ,i−Q.sub.base,i)=0;  (3)

[0050] where λ is a real number, P.sub.targ,i, P.sub.base,i represent a target value and a base value of the active power at a load bus, respectively, and Q.sub.targ,i, Q.sub.base,i represent a target value and a base value of the reactive power at the load bus, respectively.

[0051] Step S3, a continuation power flow equation is solved by FFHE.

[0052] For the convenience of description, equation (3) is denoted by f(V.sup.R, V.sup.I, λ)=0, where V.sup.R=(V.sub.1.sup.R, V.sub.2.sup.R, . . . , V.sub.k.sup.R, . . . ).sup.T, V.sup.I=(V.sub.1.sup.I, V.sub.2.sup.I, . . . , V.sub.k.sup.I, . . . ).sup.T. In order to solve the equation by FFHE, parameter s is introduced, and an arc-length parametrization equation is selected simultaneously therewith.

[00002]$$\begin{gathered} \begin{array}{cc}{.Math.}_k\left\{{\left(\frac{{\mathrm{dV}}_k^R}{\mathrm{ds}}\right)}^2\left(s\right)+{\left(\frac{{\mathrm{dV}}_k^I}{\mathrm{ds}}\right)}^2\left(s\right)\right\}+{\left(\frac{d\lambda }{\mathrm{ds}}\left(s\right)\right)}^2=1; \\ f\left(V^R\left(s\right),V^I\left(s\right),\lambda \left(s\right)\right)=0.& \left(4\right)\end{array} \end{gathered}$$

[0053] While the series expansions of λ(s), V.sup.R(s), V.sup.I(s):

λ(s)=Σ.sub.q≥0a.sub.0,qs.sup.q; V.sub.k.sup.R(s)=Σ.sub.q≥0a.sub.k,qs.sup.q, V.sub.k.sup.I(s)=Σ.sub.q≥0b.sub.k,qs.sup.q, k≥1

[0054] are inserted into equation (4), and a system of equations with the unknowns a.sub.0,q, a.sub.k,q, b.sub.k,q (k≥1) can be determined by comparing the coefficients of the terms with a same order of s.

[0055] After the above steps are completed, the series expansions of λ(s), V.sup.R(s), V.sup.I(s) can be obtained. Rational approximants are derived according to the existing series expansions, so as to expand a convergence domain. A value of the rational approximant at s=s.sub.0 is substituted into the equation (3), and whether a mismatch between a left side and a right side of the equation is smaller than a preset acceptable threshold is compared. s.sub.0 as great as possible is found out, so that the mismatch is smaller than the threshold. By taking V.sup.R(s.sub.0), V.sup.I(s.sub.0), λ(s.sub.0) as a new initial point, the above steps are repeated till the SNB point is pinpointed.

[0056] Step S4, a scheduling policy of the power system are designed and planned based on the solved load margin.

[0057] As it needs to reserve a certain load margin to cope with change of the load during the operation of a power distribution network, it needs to determine a constraint condition of the scheduling model based on the load margin and design an optimum scheduling policy. When the load margin is as low as a certain degree, a part of load is removed immediately, and at the time, a scheduler shall draw a clear distinction between the primary and the secondary to guarantee the power supply of main areas. Areas with probable voltage collapse can be configured with enough automatic load shedding devices to prevent accidents.

[0058] During the specific implementation, the above-mentioned step S3 includes:

[0059] Step S301, an order of series expansion q.sub.max, a threshold of the acceptable mismatch e between two sides of the equation, and a given entry of at least one part in the initial point X.sub.0=(V.sup.R(0), V.sup.I(0), λ(0)) are initialized and set.

[0060] Step S302, a new parameter s is introduced and an embedding system is obtained:

[00003]${.Math.}_k\left\{{\left(\frac{{\mathrm{dV}}_k^R}{\mathrm{ds}}\right)}^2\left(s\right)+{\left(\frac{{\mathrm{dV}}_k^I}{\mathrm{ds}}\right)}^2\left(s\right)\right\}+{\left(\frac{d\lambda }{\mathrm{ds}}\left(s\right)\right)}^2=1;$$V_i^R\left(s\right){.Math.}_k\left\{G_{i,k}{V}_k^R\left(s\right)-B_{i,k}{V}_k^I\left(s\right)\right\}+V_i^I\left(s\right){.Math.}_k\left\{B_{i,k}{V}_k^R\left(s\right)+G_{i,k}^I{V}_k^I\left(s\right)\right\}-P_i-\lambda \left(s\right).Math.\left(P_{\mathrm{targ},i}-P_{\mathrm{base},i}\right)=0;$$-V_i^R\left(s\right){.Math.}_k\left\{B_{i,k}{V}_k^R\left(s\right)+G_{i,k}{V}_k^I\left(s\right)\right\}+V_i^I\left(s\right){.Math.}_k\left\{G_{i,k}{V}_k^R\left(s\right)-B_{i,k}{V}_k^I\left(s\right)\right\}-Q_i-\lambda \left(s\right).Math.\left(Q_{\mathrm{targ},i}-Q_{\mathrm{base},i}\right)=0;$

[0061] Step S303, λ(s), V.sup.R(s), V.sup.I(s) are represented as series expansions in s:

λ(s)=Σ.sub.q≥0a.sub.0,qs.sup.q; V.sub.k.sup.R(s)=Σ.sub.q≥0a.sub.k,qs.sup.q, V.sub.k.sup.I(s)=Σ.sub.q≥0b.sub.k,qs.sup.q, k≥1

[0062] Step S304, the series expansion is inserted into the embedding system, and a system of equations is obtained by taking series coefficients in the series expansions as the unknowns:

Σ.sub.k≠ref{(Σ.sub.q≥0(1+q).Math.a.sub.k,q+1s.sup.q).sup.2+(Σ.sub.q≥0(1+q).Math.b.sub.k,q+1s.sup.q).sup.2}+(Σ.sub.q≥0(1+q).Math.a.sub.0,q+1s.sup.q).sup.2=1;

{Σ.sub.q≥0a.sub.i,qs.sup.q}.Math.Σ.sub.k{G.sub.i,kΣ.sub.q≥0a.sub.k,qs.sup.q−B.sub.i,kΣ.sub.q≥0b.sub.k,qs.sup.q}+{Σ.sub.q≥0b.sub.i,qs.sup.q}.Math.Σ.sub.k{B.sub.i,kΣ.sub.q≥0a.sub.k,qs.sup.q+G.sub.i,kΣ.sub.q≥0b.sub.k,qs.sup.q}−P.sub.i−{Σ.sub.q≥0a.sub.0,qs.sup.q}.Math.(P.sub.targ,i−P.sub.base,i)=0;

−{Σ.sub.q≥0a.sub.i,qs.sup.q}.Math.Σ.sub.k{B.sub.i,kΣ.sub.q≥0a.sub.k,qs.sup.q+G.sub.i,kΣ.sub.q≥0b.sub.k,qs.sup.q}+{Σ.sub.q≥0b.sub.i,qs.sup.q}.Math.Σ.sub.k{G.sub.i,kΣ.sub.q≥0a.sub.k,qs.sup.q−B.sub.i,kΣ.sub.q≥0b.sub.k,qs.sup.q}−Q.sub.i−{Σ.sub.q≥0a.sub.0,qs.sup.q}.Math.(Q.sub.targ,i−Q.sub.base,i)=0;

[0063] Step S305, the coefficients of the terms with a same order of s are compared, wherein when q=0,

a.sub.k,0=V.sub.k.sup.R(0), b.sub.k,0=V.sub.k.sup.I(0), 1≤k≤n; a.sub.0,0=λ(0)=λ.sub.0;

a.sub.i,0.Math.Σ.sub.k{G.sub.i,ka.sub.k,0−B.sub.i,kb.sub.k,0}+b.sub.i,0.Math.Σ.sub.k{B.sub.i,ka.sub.k,0+G.sub.i,kb.sub.k,0}−P.sub.i−a.sub.0,0.Math.(P.sub.targ,i−P.sub.base,i)=0;

−a.sub.i,0.Math.Σ.sub.k{B.sub.i,ka.sub.k,0+G.sub.i,kb.sub.k,0}+b.sub.i,0.Math.Σ.sub.k{G.sub.i,ka.sub.k,0−B.sub.i,kb.sub.k,0}−Q.sub.i−a.sub.0,0.Math.(Q.sub.targ,i−Q.sub.base,i)=0;

[0064] at least one of a.sub.0,0, a.sub.k,0, b.sub.k,0 is selected as a known initial value, to solve all a.sub.0,0, a.sub.k,0, b.sub.k,0.

[0065] Step S306, the coefficients of the terms with a same order of s are compared, wherein when q=1,

Σ.sub.k{a.sub.k,1.sup.2+b.sub.k,1.sup.2}+(a.sub.0,1).sup.2=1;  (5)

a.sub.i,0.Math.Σ.sub.k{G.sub.i,ka.sub.k,1−B.sub.i,kb.sub.k,1}+a.sub.i,1.Math.Σ.sub.k{G.sub.i,ka.sub.k,0−B.sub.i,kb.sub.k,0}+b.sub.i,0.Math.Σ.sub.k{B.sub.i,ka.sub.k,1+G.sub.i,kb.sub.k,1}+b.sub.i,1.Math.Σ.sub.k{B.sub.i,ka.sub.k,0+G.sub.i,kb.sub.k,0}−a.sub.0,1.Math.(P.sub.targ,i−P.sub.base,i)=0;  (6)

−a.sub.i,0.Math.Σ.sub.k{B.sub.i,ka.sub.k,1+G.sub.i,kb.sub.k,1}−a.sub.i,1.Math.Σ.sub.k{B.sub.i,ka.sub.k,0+G.sub.i,kb.sub.k,0}+b.sub.i,0.Math.Σ.sub.k{G.sub.i,ka.sub.k,1−B.sub.i,kb.sub.k,1}+b.sub.i,1.Math.Σ.sub.k{G.sub.i,ka.sub.k,0−B.sub.i,kb.sub.k,0}−a.sub.0,1.Math.(Q.sub.targ,i−Q.sub.base,i)=0;  (7)

[0066] where a.sub.0,0, a.sub.k,0, b.sub.k,0 herein are known.

[0067] For the convenience of description, the to-be-solved equations (5), (6) and (7) are rewritten into

Σ.sub.i=1.sup.nx.sub.i.sup.2=1, Ax=B.  (8)

[0068] where x=(x.sub.1, . . . , x.sub.n).sup.T is an unknown vector, A is a matrix with (n−1) rows and n columns, and B is a column vector. At the time, the equation set (8) is solved to obtain a.sub.0,1, a.sub.k,1, b.sub.k,1.

[0069] How to solve the equation (8) is introduced below:

[0070] First, A is decomposed into A=[A.sub.1, A.sub.2], marked as {tilde over (x)}.sub.1={tilde over (x)}.sub.1, {tilde over (x)}.sub.2=(x.sub.2, . . . , {tilde over (x)}.sub.n).sup.T.Math.{tilde over (x)}.sub.2=A.sub.2.sup.−1(B−A.sub.1{tilde over (x)}.sub.1) can be deduced from A.sub.1{tilde over (x)}.sub.1+A.sub.2{tilde over (x)}.sub.2=B, and a quadratic equation with one unknown regarding {tilde over (x)}.sub.1 can be obtained in combination with the equation (8)

{tilde over (x)}.sub.1.sup.2+(B−A.sub.1{tilde over (x)}.sub.1).sup.T(A.sub.2.sup.−1).sup.T(A.sub.2.sup.−1)(B−A.sub.1{tilde over (x)}.sub.1)=1,

[0071] can be rewritten as

[00004]$c_0{\tilde{x}}_1^2+c_1{\tilde{x}}_1+c_2=0,$$\mathrm{where}$$c_0=1+{A_1^T\left(A_2^{-1}\right)}^T\left(A_2^{-1}\right)A_1,$$c_1=-2.Math.{A_1^T\left(A_2^{-1}\right)}^T\left(A_2^{-1}\right)B,$$c_2={B^T\left(A_2^{-1}\right)}^T\left(A_2^{-1}\right)B-1.$$\mathrm{then}$${\tilde{x}}_1=\frac{1}{2c_0}.Math.\left(-c_1\pm \sqrt{c_1^2-4c_0{c}_2}\right),$$\mathrm{thus}$$x=\left({\tilde{x}}_1;A_2^{-1}\left(B-A_1{\tilde{x}}_1\right)\right).$

[0072] Step S307, as the system of equations met by a.sub.0,q, a.sub.k,q, b.sub.k,q, (q.sub.max≥q≥2) is a linear equation set, it can be solved fast with a classic method.

[0073] Step S308, Rational approximants are derived using the series expansions of λ(s), V.sup.R(s), V.sup.I(s).

[0074] Step S309, a value of the rational approximant at s=s.sub.0 is substituted into the equation (3), and a mismatch between a left side and a right side of the equation is compared to a preset acceptable threshold e to figure out whether the mismatch is smaller than the preset acceptable threshold e. If yes, s.sub.0 is expanded till the mismatch is not smaller than the preset acceptable threshold; if no, s.sub.0 is shrunk till the mismatch is smaller than the preset acceptable threshold. That is, s.sub.0 as great as possible and smaller than the threshold e is found.

[0075] Step S310, (V.sup.R(s.sub.0), V.sup.I(s.sub.0), λ(s.sub.0)) is taken as a new initial point X.sub.0.

[0076] Step S311, S302 to S310 are repeated till an SNB point is pinpointed.

[0077] Detailed description on the preferred specific embodiments of the present invention is made above. It is to be understood that those of ordinary skill in the art can make various modifications and variations in accordance with concept of the present invention without creative efforts. Therefore, technical solutions capable of being obtained by technicians in the technical field through logical analysis, inference or limited experiments in accordance with concept of the present invention based on the prior art shall fall into the scope of protection determined by claims.