CN104252571B - WLAV robust state estimation methods based on many prediction correction interior points - Google Patents

Abstract

Based on the WLAV robust state estimation methods of many prediction correction interior points, it is related to a kind of Operation of Electric Systems and control method.State estimation uses weighted least-squares state estimation, and estimated accuracy is difficult to further improve.The present invention comprises the following steps：1）Obtain parameters of electric power system；2）Obtain detection data；3）Initialization；4）Apply for each measurement Hessian matrix memory headroom and solve；5）Calculate duality gap C_Gap=α^Tl+β^TU, judges whether to meet C_Gap＜ ε or K ＜ K_max；6）Prediction step：Discontinuous Factors μ=0 is set, is predicted according to formula, obtains affine direction；7）Correction step；8）Judge number of corrections counter t ＜ 4;9）It is rightEquation solution, obtains Δ λ_co, and correct；10）Calculate new iteration step lengthIfMore than former step-length, by formula Δ λ_new=Δ λ_af+ωΔλ_coIt is updated, t=t+1, and performs step 8)；Otherwise, iteration calculator K=K+1, and perform step 5).The technical program reduces iterations, improves processing speed；Further improve convergence of algorithm characteristic.

Description

WLAV robust state estimation methods based on many predictor-corrector interior point methods

Technical field

The present invention relates to a kind of Operation of Electric Systems and control method.

Background technology

With the implementation of transferring electricity from the west to the east, the rapid implementation of electricity market, extra-high voltage, remote, alternating current-direct current mixing transmission of electricity skill Art quickly growing in China's power network, electric power system dispatching center automatization level is also required to step up.State estimation is made For the core and foundation stone of modern energy management systems (EMS), pass through the measurement number transmitted to data acquisition monitoring (SCADA) system According to being handled in real time, electric power system data precision can be increased substantially.It is traditional based on least square method and by least square The derivative quick decoupling method state estimation procedure of method has practical operating experiences for many years at the scene.Yet with least square method It is the optimal estimation to Gaussian distributed sample, and the larger bad data of error is there may be in actual SCADA, these is not Good data disobey Gaussian Profile.Therefore, Legacy Status is difficult to further improve using weighted least-squares precision of state estimation.

For how to suppress the problem of bad data influences in state estimation, domestic and foreign scholars propose a large amount of solution party Method, mainly including following 2 aspects：One is to seek new bad data recognition method, and it is rejected from effective measurement system； Two be to propose new Robust filter device.When there is the bad data of bad leverage measurement or multiple strong correlations in SCADA, it is difficult to Pick out and completely, and Robust filter can effectively suppress the influence of bad data without raw data detection.Wherein, WLAV State estimation is due to ensure that multiple measurement residuals are 0, so as to effectively abandon bad data using correct measuring value, in recent years To be paid close attention to by many person scholars.WLAV is introduced Power system state estimation by Irving and Owen earliest, and they turn WLAV Linear programming is turned to be solved.Kotiuga thinks that the essence of WLAV estimations is the interpolation for measuring collection with Vidyasagar, therefore With certain exclusion bad data ability.Weiization proposes to be based on primal dual interior point method (primal-dual interior Point method, PDIPM) WLAV state estimations, this method has preferable numerical stability, but it is inclined to there is iterations Many shortcomings.

The content of the invention

The technical problem to be solved in the present invention and propose technical assignment be prior art to be improved with being improved, WLAV robust state estimation methods based on many predictor-corrector interior point methods are provided, bad data shadow in holddown estimation is reached Ring and improve the purpose of processing speed.Therefore, the present invention takes following technical scheme.

WLAV robust state estimation methods based on many predictor-corrector interior point methods, it is characterised in that comprise the following steps：

1）Parameters of electric power system is obtained, including：Bus numbering, the WLAV robust states based on many predictor-corrector interior point methods Method of estimation, compensating electric capacity, the branch road number of transmission line of electricity, headend node and endpoint node numbering, series resistance, series reactance, Shunt conductance, shunt susceptance, transformer voltage ratio and impedance；

2）Obtain detection data, including voltage magnitude, generator active power, generator reactive power, load wattful power Rate, reactive load power, circuit head end active power, circuit head end reactive power, line end active power and circuit end Hold reactive power；

3）Initialization, including：Initial value (rectangular coordinate system) is set to quantity of state, Lagrange multiplier and penalty factor are set Initial value, node order optimization, formation bus admittance matrix, recovery iteration calculator K=1, set maximum iteration K_maxIf, Put convergence precision and require ε；

4）Apply for each measurement Hessian matrix memory headroom and solve；

5）Calculate duality gap C_Gap=α^Tl+β^TU, judges whether to meet C_Gap＜ ε or K ＜ K_max, if so, then output calculates knot Really, circulation is exited；If it is not, then performing step 6)；

6）Prediction step：Discontinuous Factors μ=0 is set, is predicted according to below equation, obtains affine direction：

In formula：X is quantity of state, includes node voltage amplitude and phase angle；L, u are slack variable, i.e., former variable；Y, α, β are glug Bright day multiplier, i.e. dual variable；Δ x Δs x, Δ y, Δ w, Δ l, Δ u, Δ α, Δ β are respectively x, y, w, l, u, α, β correction; ▽_xH (x) is h (x) Jacobian matrix,For h (x) Hessian matrix, μ is Discontinuous Factors, L=diag (l₁,…,l_m)、U =diag (u₁,…,u_m), A=diag (a₁,…,a_m), B=diag (β₁,…,β_m), e=[1 ..., 1]^T,Utilize damped Newton method pairSolve, by Formula solves Δ λ, and former, dual variable is modified：λ^(k+1)=λ^(k)+ α Δs λ, α are iteration step length, and its size is defined below：

Calculate the complementary gap in affine direction：

C_gap ^af=(α+α^*Δα_af)^T(l+α^*Δl_af)+(β+α^*Δβ_af)^T(u+α^*Δu_af)

Dynamic estimation Center Parameter：

μ=min { (C_gap ^af/C_gap)³,0.1}C_gap/2m

7）Correction step, sets number of corrections counter t=1；

8）Number of corrections counter t ＜ 4 are judged, if so, then performing step 9);If it is not, then performing step 5);

9）It is rightEquation is solved, and obtains Δ λ_co, and use dynamic select orientation Shared weight, i.e. Δ λ in total Newton direction_new=Δ λ_af+ωΔλ_co；Weight is scanned for using 2 terrace works, I.e.：1st stage, in [α_pα_d, 1] and middle progress linear search, and allow to find different optimal weights in former, dual spaces；True Behind fixed efficient best initial weights search subinterval, the 2nd stage found optimal weights in the subinterval, it is similar with the stage 1, Linear search is carried out in subinterval, is eventually found optimalWith

10）Calculate new iteration step lengthIfMore than former step-length, by formula Δ λ_new=Δ λ_af+ ωΔλ_coIt is updated, number of corrections counter t=t+1, and performs step 8)；Otherwise, iteration calculator K=K+1, and perform Step 5).

The technical program retains the order of information of nonlinear terms, and increases iteration step length using the mode repeatedly corrected, subtracts Few iterations, improves processing speed.Meanwhile, processing is weighted to orientation, is searched out in total Newton direction Optimal proportion, so as to ensure that iteration point is drawn close to centrode, further improves convergence of algorithm characteristic.

As further improving and supplementing to above-mentioned technical proposal, present invention additionally comprises following additional technical feature.

t_max=4。

Duality gap C_GapComputing formula is α^Tl+β^Tu。

When the 6th step is predicted and walked, first to L (x, y, l, u, α, β)=c^T(l+u)-y^T[z-h(x)+l-u]-α^Tl-β^Tu

Derivation is carried out according to KKT conditions, is obtained：

Beneficial effect：The technical program can effectively suppress the influence of bad data, effectively be got rid of using correct measuring value Bad data is abandoned, with preferable numerical stability.And retain the order of information of nonlinear terms, and utilize the mode repeatedly corrected Increase iteration step length, further reduce iterations, improve processing speed.Meanwhile, processing is weighted to orientation, is found The optimal proportion gone out in total Newton direction, so as to ensure that iteration point is drawn close to centrode, further improves convergence property.

Brief description of the drawings

Fig. 1 is the inventive method flow chart.

Fig. 2（a）It is the circuit Π shape equivalent circuit diagrams that the present invention is used.

Fig. 2（b）It is the transformer Π shape equivalent circuit diagrams that the present invention is used.

Fig. 3 is the IEEE-14 standard example systems that the present invention is applied.

Fig. 4 is that the iteration step length of the different interior point methods of IEEE-57 node systems in example two compares.

Fig. 5 is that the duality gap convergence process of the different interior point methods of IEEE-57 node systems in example two compares.

Fig. 6 is ω and α in the search of certain iteration phase 1 of IEEE-57 node systems in example two_p,α_dRelation.

Embodiment

Technical scheme is described in further detail below in conjunction with Figure of description.

As shown in figure 1, the present invention comprises the following steps：

1）Parameters of electric power system is obtained, including：Bus numbering, the WLAV robust states based on many predictor-corrector interior point methods Method of estimation, compensating electric capacity, the branch road number of transmission line of electricity, headend node and endpoint node numbering, series resistance, series reactance, Shunt conductance, shunt susceptance, transformer voltage ratio and impedance；

2）Obtain detection data, including voltage magnitude, generator active power, generator reactive power, load wattful power Rate, reactive load power, circuit head end active power, circuit head end reactive power, line end active power and circuit end Hold reactive power；

3）Initialization, including：Initial value (rectangular coordinate system) is set to quantity of state, Lagrange multiplier and penalty factor are set Initial value, node order optimization, formation bus admittance matrix, recovery iteration calculator K=1, set maximum iteration K_maxIf, Put maximum correction number of times t_max, set convergence precision to require ε；

4）Apply for each measurement Hessian matrix memory headroom and solve；

5）Calculate duality gap C_Gap, judge whether to meet C_Gap＜ ε, if so, then exporting result of calculation, exit circulation；If It is no, then perform step 6)；

6）Prediction step：Discontinuous Factors μ=0 is set, is predicted according to below equation, obtains affine direction：

In formula：X is quantity of state, includes node voltage amplitude and phase angle；L, u are slack variable, i.e., former variable；Y, α, β are glug Bright day multiplier, i.e. dual variable；Δ x Δs x, Δ y, Δ w, Δ l, Δ u, Δ α, Δ β are respectively x, y, w, l, u, α, β correction; ▽_xH (x) is h (x) Jacobian matrix,Hessian matrix, μ is Discontinuous Factors, L=diag (l₁,…,l_m)、U =diag (u₁,…,u_m), A=diag (a₁,…,a_m), B=diag (β₁,…,β_m), e=[1 ..., 1]^T,Profit With damped Newton method pairSolve, byFormula solves Δ λ, and right Former, dual variable is modified：λ^(k+1)=λ^(k)+ α Δs λ, α are iteration step length, and its size is defined below：

Calculate the complementary gap in affine direction：

C_gap ^af=(α+α^*Δα_af)^T(l+α^*Δl_af)+(β+α^*Δβ_af)^T(u+α^*Δu_af)

Dynamic estimation Center Parameter：

μ=min { (C_gap ^af/C_gap)³,0.1}C_gap/2m

7）Correction step, sets number of corrections counter t=1；

8）Judge number of corrections counter t<t_tmax, if so, then performing step 9);If it is not, then performing step 5);

9）It is rightEquation is solved, and obtains Δ λ_co, and use dynamic select orientation Shared weight, i.e. Δ λ in total Newton direction_new=Δ λ_af+ωΔλ_co；Weight is scanned for using 2 terrace works, I.e.：1st stage, in [α_pα_d, 1] and middle progress linear search, and allow to find different optimal weights in former, dual spaces；True Behind fixed efficient best initial weights search subinterval, the 2nd stage found optimal weights in the subinterval, it is similar with the stage 1, Linear search is carried out in subinterval, is eventually found optimalWith

10）Calculate new iteration step lengthIfMore than former step-length, by formula Δ λ_new=Δ λ_af+ω Δλ_coIt is updated, number of corrections counter t=t+1, and performs step 8)；Otherwise, iteration calculator K=K+1, and perform step It is rapid 5).

Under conditions of given network connection, branch parameters and measurement system, non-linear measurement equation is represented by：

Z=h (x)+r

In formula, z is measuring value vector, and the overwhelming majority is obtained by remote measurement, also has sub-fraction manually to set； H (x) is the measurement function set up by the basic circuit law such as kirchhoff；X is system state variables；R is measurement with chance error Difference.If system has n node, using node complex voltage real part and imaginary part as state variable, balance nodes are not involved in iteration.

Such as Fig. 2（a）, Fig. 2（b）It is shown.In Power system state estimation, the type of measurement configuration is than conventional Load Flow Many, the injecting power for not only including each node measures P_i、Q_i, the power measurement P of branch road can also be included_ij、Q_ij、P_ji、Q_jiWith And the voltage magnitude of node measures V_i, measurement equation is shown below：

Node i voltage：

Node i injecting power：

Top power on circuit i-j：

End power on circuit i-j：

Top power on transformer lines i-j：

End power on transformer lines i-j：

By taking the top power on node injecting power and circuit i-j as an example, Jacobian matrix H (x) correspondences of vector are measured Element be：

Node i injecting power：

(j=1,2 ..., i-1, i+1 ..., n)

(j=1,2 ..., i-1, i+1 ..., n)

Top power on circuit i-j

The Hessian matrix of each measurement equation is obtained, the step only need to be calculated and stored before iteration, nothing in iterative process It need to be calculated again.By taking the top power on circuit i-j as an example, the corresponding elements of Hessian matrix He (x) for measuring vector are：

In formula, e_i、e_jRespectively node i, node j complex voltage real parts；G_ij、B_ijFor the element of admittance matrix；g、b、y_cFor Parameter in circuit Π shape models；K is the non-standard no-load voltage ratio of transformer；b_TFor the susceptance of transformer standard side.Circuit and transformation The Π shape equivalent circuits of device are as shown in Figure 2.

In order to verify effectiveness of the invention, it is compared with different conditions method of estimation.Wherein WLS state estimations are referred to as Method 1, WLS+BD is referred to as method 2, and predictor-corrector interior point method (predictor-collector PDIPM, PCPDIPM) is referred to as Method 3, many predictor-corrector interior point methods (multiple PCPDIPM, MPCPDIPM) are referred to as method 4, linear primal dual interior point method Referred to as method 5.More to meet actual physical meaning, each method state variable estimated result is provided by polar form, converted Method is：

Two embodiments of the present invention are described below：

The present invention divides 3 kinds of situations to be calculated using the standard example of the IEEE-14 nodes shown in Fig. 3.1) system is not containing Good data, and normal measure takes actual value, for analyzing identification capability of the distinct methods to bad data；2) system is not without Good data, and all measurements carry random error, for analyzing the distinct methods under the normal measurement for meeting normal distribution State estimation performance；3) system contain bad data, and measurement is with random error, i.e. situation 1), synthesis 2), this is also Most common a kind of pattern in practice.Table 1 gives different methods of estimation in situation 1) under bad data debate knowledge result, table 2 Distinct methods state variable estimated result is given to compare.

The distinct methods situation 1 of table 1) under state estimation result

The distinct methods 3 of table 2）State variable error in the case of kind

From table 1,2, method 2 does not identify bad data 5, No. 7, and all measurement evaluated errors of method 3 and 4 Respectively less than 1%, bad data has been picked out well, meets engine request.Although method 5 also can pick out bad data, Evaluated error is significantly greater than method 3 and 4.Experiment shows above, and (method 4) of the invention compares linear primal dual interior point method, larger Solving precision is improved to degree, with stronger robustness.

Example two：

In order to verify the computational efficiency of this method, 6 systems to 14 to 3012 nodes are tested, wherein between antithesis Gap convergence precision is set to 10^-6.Table 3 is given under different system, utilizes the iterations of 3 kinds of different interior point methods, Yi Jiyu PDIPM compares reduced iterations percentage.Table 4 gives the calculating time of correlation method and the reduction compared with PDIPM Calculate percentage of time.

The distinct methods iterations of table 3 compares

Note：Wp-2383 in table, 3012 be Polish network system；Boldface is best result in table, similarly hereinafter

The distinct methods CPU of the table 4 calculating times compare

From table 3,4, method 3 is compared with PDIPM, and iterations can at least reduce 13.64%, and the calculating time can subtract Few 13.84%, and method 4 is on the basis of method 3, by repeatedly correcting, iterations can reduce 20.83%, and the time of calculating can More than 16.28% is reduced, with faster convergence rate.

Further, Fig. 4 gives the iteration step length of the different Interior-point methods of IEEE-57 node systems, and Fig. 5 is given between antithesis The convergence process of gap.As shown in Figure 4, PDIPM step-lengths in 3 times iterative process are larger, are gradually reduced afterwards, method 3 due to Correction step is added, is increased compared to PDIPM step-lengths, but at iteration initial stage, iteration step length diminishes on the contrary, and this is due to school Positive direction is not pointing at optimal direction, and method 4 uses many predictions of weighting, and finds out optimal correction weights guarantor using 2 terrace works Demonstrate,prove orientation and point to optimal direction, it is only necessary to which iteration can restrain for 11 times, and iteration step length remains higher value.As shown in Figure 5, Restrained after PDIPM iteration 16 times, the duality gap convergence rate of method 3 is very fast, it is only necessary to can restrain for 13 times, compared to PDIPM iteration Number of times reduces 18.75%, and uses after method 4, due in each iterative process, being repeatedly corrected, so as to obtain most Excellent iteration step length, duality gap speed is significantly improved, and only iteration can restrain for 11 times, and 31.25% is reduced compared with PDIPM.

Fig. 6 gives certain iteration ω in interval [α_pα_d, 1] find best initial weights search procedure, wherein α_p=0.5524, α_d=0.7606.In order to obtain the optimal step size of former, dual variable, former, dual variable is allowed to take different optimal power herein Value.Fig. 6 is the 1st stage Search process, finds out α_pOptimal subinterval [0.47818,0.59414], α_dOptimal subinterval be [0.82606,0.94202].Optimal step size and weight are mostly near parabolic relation, therefore during 1 stage Search, set-point Number does not require too many, as long as ensure that obtained subinterval is global optimum.Search is set in the 2nd stage subinterval Count as 4, actually often error is larger, and the present embodiment is set to 10, to ensure Search Results closer to true optimal Weights.This final optimizing result is α_p=0.88857, α_d=0.80510, obtained necessarily compared to former iteration step length (i.e. ω=1) The increase of degree.The result of the test of example two shows that the 2 stage linear search methods that this method is used find the effective of optimal weights Property.

The WLAV robust state estimation methods based on many predictor-corrector interior point methods shown in figure 1 above -6 are of the invention Specific embodiment, has embodied substantive distinguishing features of the present invention and progress, can be according to actual use needs, in opening for the present invention Under showing, carry out the equivalent modifications in terms of shape, structure to it, this programme protection domain row.

Claims (4)

1. a kind of WLAV robust state estimation methods based on many predictor-corrector interior point methods, it is characterised in that comprise the following steps：

1) parameters of electric power system is obtained, including：Bus numbering, a kind of WLAV robust states based on many predictor-corrector interior point methods Method of estimation, compensating electric capacity, the branch road number of transmission line of electricity, headend node and endpoint node numbering, series resistance, series reactance, Shunt conductance, shunt susceptance, transformer voltage ratio and impedance；

2) detection data are obtained, including it is voltage magnitude, generator active power, generator reactive power, load active power, negative Lotus reactive power, circuit head end active power, circuit head end reactive power, line end active power and line end are idle Power；

3) initialize, including：Initial value (rectangular coordinate system) is set to quantity of state, Lagrange multiplier and penalty factor are set just Value, node order optimization, formation bus admittance matrix, recovery iteration calculator K=1, set maximum iteration K_max, set Maximum correction number of times t_max, set convergence precision to require ε；

4) apply for each measurement Hessian matrix memory headroom and solve；

5) duality gap C is calculated_Gap, judge whether to meet C_Gap＜ ε, if so, then exporting result of calculation, exit circulation；If it is not, then Perform step 6)；

6) prediction step：Discontinuous Factors μ=0 is set, is predicted according to below equation, obtains affine direction：

$\left[\begin{array}{cccccc}A& L& 0& 0& 0& 0\\ 0& -I& 0& 0& 0& -I\\ 0& 0& B& U& 0& 0\\ 0& 0& 0& -I& 0& I\\ 0& 0& 0& 0& {\&dtri;}_x^{2}h\left(x\right)y& {\&dtri;}_x^{T}h\left(x\right)\\ 0& 0& 0& 0& {\&dtri;}_{x}h\left(x\right)& H\end{array}\right]\left[\begin{array}{c}\mathrm{\&Delta;}\mathrm{\&alpha;}\\ \mathrm{\&Delta;}l\\ \mathrm{\&Delta;}\mathrm{\&beta;}\\ \mathrm{\&Delta;}u\\ \mathrm{\&Delta;}x\\ \mathrm{\&Delta;}y\end{array}\right]=\left[\begin{array}{c}-ALe\\ -c+y+\mathrm{\&alpha;}\\ -BUe\\ -c-y+\mathrm{\&beta;}\\ -{\&dtri;}_{x}h\left(x\right)y\\ -L_y+M\end{array}\right]+\left[\begin{array}{c}\mathrm{\&mu;}e\\ 0\\ \mathrm{\&mu;}e\\ 0\\ 0\\ 0\end{array}\right]+\left[\begin{array}{c}-\mathrm{\&Delta;}A\mathrm{\&Delta;}Le\\ 0\\ -\mathrm{\&Delta;}B\mathrm{\&Delta;}Ue\\ 0\\ 0\\ 0\end{array}\right]$

In formula：X is quantity of state, includes node voltage amplitude and phase angle；L, u are slack variable, i.e., former variable；Y, α, β are that glug is bright Day multiplier, i.e. dual variable；Δ x, Δ y, Δ w, Δ l, Δ u, Δ α, Δ β are respectively x, y, w, l, u, α, β correction； For h (x) Jacobian matrix,For h (x) Hessian matrix, μ is Discontinuous Factors, L=diag (l₁,…,l_m), U=diag (u₁,…,u_m), A=diag (α₁,…,α_m), B=diag (β₁,…,β_m), e=[1 ..., 1]^T, Utilize damped Newton method pairSolve, byFormula solves Δ λ, and to it is former, Dual variable is modified：λ^(k+1)=λ^(k)+ α Δs λ, α are iteration step length, and its size is defined below：

${\mathrm{\&alpha;}}_p=0.9995min\{\underset{i}{min}(\frac{-l_i}{{\mathrm{\&Delta;l}}_i},{\mathrm{\&Delta;l}}_i<0;\frac{-u_i}{{\mathrm{\&Delta;u}}_i},{\mathrm{\&Delta;u}}_i<0),1\}$

${\mathrm{\&alpha;}}_d=0.9995min\{\underset{i}{min}(\frac{-{\mathrm{\&alpha;}}_i}{{\mathrm{\&Delta;\&alpha;}}_i},{\mathrm{\&Delta;\&alpha;}}_i<0;\frac{-{\mathrm{\&beta;}}_i}{{\mathrm{\&Delta;\&beta;}}_i},{\mathrm{\&Delta;\&beta;}}_i<0),1\}$

Calculate the complementary gap in affine direction：

C_gap ^af=(α+α^*Δα_af)^T(l+α^*Δl_af)+(β+α^*Δβ_af)^T(u+α^*Δu_af)

Dynamic estimation Center Parameter：

μ=min { (C_gap ^af/C_gap)³,0.1}C_gap/2m

7) correction step, sets number of corrections counter t=1；

8) number of corrections counter t is judged<t_tmax, if so, then performing step 9)；If it is not, then performing step 5)；

9) it is rightEquation is solved, and obtains Δ λ_co, and using dynamic select orientation total Newton direction in shared weight, i.e. Δ λ_new=Δ λ_af+ωΔλ_co；Weight is scanned for using 2 terrace works, i.e.,：1st Stage, in [α_pα_d, 1] and middle progress linear search, and allow to find different optimal weights in former, dual spaces；It is determined that efficiently Best initial weights search subinterval after, the 2nd stage found optimal weights in the subintervalIt is similar with the stage 1, in sub-district Between carry out linear search, eventually find optimalWith

10) new iteration step length is calculatedIfMore than former step-length, by formula Δ λ_new=Δ λ_af+ωΔλ_co It is updated, number of corrections counter t=t+1, and performs step 8)；Otherwise, iteration calculator K=K+1, and perform step 5)。

2. a kind of WLAV robust state estimation methods based on many predictor-corrector interior point methods according to claim 1, it is special Levy and be:t_max=4.

3. a kind of WLAV robust state estimation methods based on many predictor-corrector interior point methods according to claim 1, it is special Levy and be:Duality gap C_GapComputing formula is α^Tl+β^Tu。

4. a kind of WLAV robust state estimation methods based on many predictor-corrector interior point methods according to claim 1, it is special Levy and be:When the 6th step is predicted and walked, first to L (x, y, l, u, α, β)=c^T(l+u)-y^T[z-h(x)+l-u]-α^Tl-β^Tu

Derivation is carried out according to KKT conditions, is obtained：

$\left\{\begin{array}{c}{L}_{\mathrm{\&alpha;}}=-l\&DoubleRightArrow;L_{\mathrm{\&alpha;}}^{\mathrm{\&mu;}}=ALe-\mathrm{\&mu;}e=0\\ L_l=c-y-\mathrm{\&alpha;}=0\\ L_{\mathrm{\&beta;}}=-u\&DoubleRightArrow;L_{\mathrm{\&beta;}}^{\mathrm{\&mu;}}=BUe-\mathrm{\&mu;}e=0\\ L_u=c+y-\mathrm{\&beta;}=0\\ L_x={\&dtri;}_{x}h\left(x\right)y=0\\ L_y=-z+h\left(x\right)-l+u=0\end{array}\right.$

In formula, z is measuring value vector.