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

Abstract

The invention discloses a WLAV robust state estimation method based on a multi-prediction-correction interior point method, and relates to a power system operation and control method. State estimation adopts weighted least square state estimation, and it is difficult to further improve the estimation accuracy. The present invention includes the following steps: 1) Acquiring power system parameters; 2) Obtaining detection data; 3) Initializing; 4) Applying for each quantity to measure Hessian matrix memory space and solving; 5) Calculating the dual gap C _Gap = α ^T l + β ^T u, judging whether C _Gap <ε or K<K _max is satisfied; 6) Prediction step: set the disturbance factor μ = 0, predict according to the formula, and find the affine direction; 7) Correction step; 8) Judgment correction times counter t<4; 9) for Solve the equation, get Δλ _co , and correct it; 10) Calculate the new iteration step size like If it is greater than the original step size, update it according to the formula Δλ _new =Δλ _af +ωΔλ _co , t=t+1, and execute step 8); otherwise, iterate the calculator K=K+1, and execute step 5). The technical proposal reduces the number of iterations, improves the processing speed, and further improves the convergence characteristic of the algorithm.

Description

WLAV Robust State Estimation Method Based on Multi-Prediction-Correction Interior Point Method

technical field

The invention relates to a power system operation and control method.

Background technique

With the implementation of west-to-east power transmission and the rapid implementation of the power market, UHV, long-distance, and AC-DC hybrid transmission technologies are developing rapidly in my country's power grids, and the automation level of power system dispatching centers also needs to be gradually improved. As the core and cornerstone of modern energy management system (EMS), state estimation can greatly improve the accuracy of power system data through real-time processing of measurement data transmitted by SCADA system. The traditional state estimation program based on the least square method and the fast decoupling method derived from the least square method has many years of actual operation experience in the field. However, since the least squares method is the best estimate for samples that obey the Gaussian distribution, there may be bad data with large errors in the actual SCADA, and these bad data do not obey the Gaussian distribution. Therefore, it is difficult to further improve the accuracy of traditional state estimation using weighted least squares.

Aiming at the problem of how to suppress the influence of bad data in state estimation, scholars at home and abroad have proposed a large number of solutions, mainly including the following two aspects: one is to seek new bad data identification methods and eliminate them from the effective measurement system; A new robust estimator is proposed. When there are bad leverage measurements or multiple strongly correlated bad data in SCADA, it is difficult to fully identify them, while robust estimation does not require bad data detection and can effectively suppress the influence of bad data. Among them, WLAV state estimation has attracted the attention of many scholars in recent years because it can ensure that multiple measurement residuals are 0, so that the correct measurement values can be effectively used to discard bad data. Irving and Owen first introduced WLAV into power system state estimation, and they transformed WLAV into linear programming for solution. Kotiuga and Vidyasagar believe that the essence of WLAV estimation is the interpolation of measurement sets, so it has a certain ability to exclude bad data. Wei Hua proposed a WLAV state estimation based on the primal-dual interior point method (PDIPM). This method has good numerical stability, but has the disadvantage of too many iterations.

Contents of the invention

The technical problem to be solved and the technical task proposed by the present invention are to improve and improve the existing technical solutions, and provide a WLAV robust state estimation method based on the multi-prediction-correction interior point method, so as to suppress the influence of bad data in the state estimation, and The purpose of improving processing speed. For this reason, the present invention takes the following technical solutions.

The WLAV robust state estimation method based on multi-prediction-correction interior point method is characterized in that comprising the following steps:

1) Obtain power system parameters, including: bus number, WLAV tolerance state estimation method based on multi-prediction-correction interior point method, compensation capacitor, branch number of transmission line, head-end node and end node number, series resistance, series connection Reactance, parallel conductance, parallel susceptance, transformer ratio and impedance;

2) Obtain detection data, including voltage amplitude, generator active power, generator reactive power, load active power, load reactive power, line head end active power, line head end reactive power, line end active power and line Terminal reactive power;

3) Initialization, including: setting the initial value of the state quantity (rectangular coordinate system), setting the initial value of the Lagrange multiplier and the penalty factor, optimizing the node order, forming the node admittance matrix, and restoring the iterative calculator K=1, Set the maximum number of iterations K _max , and set the convergence accuracy requirement ε;

4) Apply each quantity to measure the memory space of the Hessian matrix and solve it;

5) Calculate the dual gap C _Gap = α ^T l + β ^T u, judge whether C _Gap <ε or K < K _max , if yes, output the calculation result, and exit the loop; if not, execute step 6);

6) Prediction step: set the disturbance factor μ=0, predict according to the following formula, and find the affine direction:

In the formula: x is the state quantity, including node voltage amplitude and phase angle; l and u are slack variables, namely original variables; y, α, β are Lagrangian multipliers, namely dual variables; ΔxΔx, Δy, Δw , Δl, Δu, Δα, Δβ are the corrections of x, y, w, l, u, α, β respectively; ▽ _x h(x) is the Jacobian matrix of h(x), is the Hessian matrix of h(x), μ is the disturbance factor, L=diag(l ₁ ,…,l _m ), U=diag(u ₁ ,…,u _m ), A=diag(a ₁ ,…, a _m ), B=diag(β ₁ ,…,β _m ), e=[1,…,1] ^T , Using damped Newton's method for solve, by Δλ is solved by the formula, and the original and dual variables are corrected: λ ^(k+1) = λ ^(k) + αΔλ, α is the iteration step size, and its size is determined as follows:

Compute the complementary gap in the affine direction:

C _gap ^af ＝(α+α ^* Δα _af ) ^T (l+α ^* Δl _af )+(β+α ^* Δβ _af ) ^T (u+α ^* Δu _af )

Dynamic estimation center parameters:

μ＝min{(C _gap ^af /C _gap ) ³ ,0.1}C _gap /2m

7) Correction step, set the correction times counter t=1;

8) Judging that the number of correction times counter t<4, if yes, execute step 9); if not, execute step 5);

9) yes Solve the equation to get Δλ _co , and use the dynamic selection of the weight of the correction direction in the total Newton direction, that is, Δλ _new = Δλ _af + ωΔλ _co ; use the two-stage method to search for the weight, namely: the first stage, Perform a linear search in [α _p α _d ,1], and allow different optimal weights to be found in the original and dual spaces; after determining the efficient optimal weight search sub-interval, the second stage finds the optimal Excellent weight , similar to stage 1, conduct a linear search in the subinterval, and finally find the optimal with

10) Calculate the new iteration step size like If it is greater than the original step length, update it according to the formula Δλ _new =Δλ _af +ωΔλ _co , correct the number of times counter t=t+1, and execute step 8); otherwise, iterate the calculator K=K+1, and execute step 5).

The technical solution retains the high-order information of the nonlinear item, and uses multiple corrections to increase the iteration step size, reduce the number of iterations, and improve the processing speed. At the same time, the correction direction is weighted to find the optimal proportion in the total Newton direction, so as to ensure that the iteration point is close to the center trajectory, and further improve the convergence characteristics of the algorithm.

As a further improvement and supplement to the above technical solutions, the present invention also includes the following additional technical features.

t _max =4.

The formula for calculating the duality gap C _Gap is α ^T l+β ^T u.

In the 6th step of prediction, firstly, L(x,y,l,u,α,β)=c ^T (l+u)-y ^T [zh(x)+lu]-α ^T l-β ^T u

Derivation according to the KKT conditions, we get:

Beneficial effects: the technical solution can effectively suppress the influence of bad data, effectively use correct measurement values to discard bad data, and has better numerical stability. And retain the high-order information of nonlinear items, and use multiple corrections to increase the iteration step size, further reduce the number of iterations, and improve the processing speed. At the same time, the correction direction is weighted, and the optimal proportion in the total Newton direction is found, so as to ensure that the iteration point is close to the center trajectory, and further improve the convergence characteristics.

Description of drawings

Fig. 1 is a flow chart of the method of the present invention.

Fig. 2 (a) is the circuit Π-shaped equivalent circuit diagram adopted by the present invention.

Fig. 2 (b) is the transformer Π-shaped equivalent circuit diagram adopted in the present invention.

Fig. 3 is the IEEE-14 standard calculation example system applied by the present invention.

Fig. 4 is a comparison of the iteration step size of different interior point methods for the IEEE-57 node system in Calculation Example 2.

Fig. 5 is a comparison of the dual gap convergence process of different interior point methods for the IEEE-57 node system in Calculation Example 2.

Figure 6 shows the relationship between ω and α _p , α _d in an iterative stage 1 search of the IEEE-57 node system in Calculation Example 2.

detailed description

The technical solution of the present invention will be further described in detail below in conjunction with the accompanying drawings.

As shown in Figure 1, the present invention comprises the following steps:

1) Obtain power system parameters, including: bus number, WLAV tolerance state estimation method based on multi-prediction-correction interior point method, compensation capacitor, branch number of transmission line, head-end node and end node number, series resistance, series connection Reactance, parallel conductance, parallel susceptance, transformer ratio and impedance;

2) Obtain detection data, including voltage amplitude, generator active power, generator reactive power, load active power, load reactive power, line head end active power, line head end reactive power, line end active power and line Terminal reactive power;

3) Initialization, including: setting the initial value of the state quantity (rectangular coordinate system), setting the initial value of the Lagrange multiplier and the penalty factor, optimizing the node order, forming the node admittance matrix, and restoring the iterative calculator K=1, Set the maximum number of iterations K _max , set the maximum number of corrections t _max , and set the convergence accuracy requirement ε;

4) Apply each quantity to measure the memory space of the Hessian matrix and solve it;

5) Calculate the dual gap C _Gap , judge whether C _Gap <ε is satisfied, if yes, output the calculation result, and exit the loop; if not, execute step 6);

6) Prediction step: set the disturbance factor μ=0, predict according to the following formula, and find the affine direction:

In the formula: x is the state quantity, including node voltage amplitude and phase angle; l and u are slack variables, namely original variables; y, α, β are Lagrangian multipliers, namely dual variables; ΔxΔx, Δy, Δw , Δl, Δu, Δα, Δβ are the corrections of x, y, w, l, u, α, β respectively; ▽ _x h(x) is the Jacobian matrix of h(x), The Hessian matrix, μ is the disturbance factor, L=diag(l ₁ ,…,l _m ), U=diag(u ₁ ,…,u _m ), A=diag(a ₁ ,…,a _m ), B = diag(β ₁ ,...,β _m ), e=[1,...,1] ^T , Using damped Newton's method for solve, by Δλ is solved by the formula, and the original and dual variables are corrected: λ ^(k+1) = λ ^(k) + αΔλ, α is the iteration step size, and its size is determined as follows:

Compute the complementary gap in the affine direction:

C _gap ^af ＝(α+α ^* Δα _af ) ^T (l+α ^* Δl _af )+(β+α ^* Δβ _af ) ^T (u+α ^* Δu _af )

Dynamic estimation center parameters:

μ＝min{(C _gap ^af /C _gap ) ³ ,0.1}C _gap /2m

7) Correction step, set the correction times counter t=1;

8) Judging the number of corrections counter t<t _tmax , if yes, execute step 9); if not, execute step 5);

9) yes Solve the equation to get Δλ _co , and use the dynamic selection of the weight of the correction direction in the total Newton direction, that is, Δλ _new = Δλ _af + ωΔλ _co ; use the two-stage method to search for the weight, namely: the first stage, Perform a linear search in [α _p α _d ,1], and allow different optimal weights to be found in the original and dual spaces; after determining the efficient optimal weight search sub-interval, the second stage finds the optimal Excellent weight , similar to stage 1, conduct a linear search in the subinterval, and finally find the optimal with

10) Calculate the new iteration step size like If it is greater than the original step length, update it according to the formula Δλ _new =Δλ _af +ωΔλ _co , correct the number of times counter t=t+1, and execute step 8); otherwise, iterate the calculator K=K+1, and execute step 5).

Under the condition of given network connection, branch parameters and measurement system, the nonlinear measurement equation can be expressed as:

z=h(x)+r

In the formula, z is the measurement value vector, most of which are obtained through telemetry, and a small part is manually set; h(x) is the measurement function established by basic circuit laws such as Kirchhoff; x is the system state variable; r is the measurement random error. Assuming that the system has n nodes, the real and imaginary parts of the node complex voltage are used as state variables, and the balance nodes do not participate in the iteration.

As shown in Figure 2(a) and Figure 2(b). In power system state estimation, there are more types of measurement configuration than conventional power flow, including not only the injected power measurements P _i , Q _i of each node, but also the branch power measurements P _ij , Q _ij , P _ji , Q _ji and the node voltage amplitude measurement V _i , the measurement equation is as follows:

Node i voltage:

Node i injected power:

Start-end power on line i-j:

Terminal power on line i-j:

Power at the beginning of the transformer line i-j:

Terminal power on transformer line i-j:

Taking the node injection power and the head end power on the line i-j as an example, the elements corresponding to the Jacobian matrix H(x) of the measurement vector are:

Node i injected power:

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

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

Head end power on line i-j

To find the Hessian matrix of each measurement equation, this step only needs to be calculated and stored before the iteration, and no further calculation is required during the iteration. Taking the initial power on the line i-j as an example, the elements corresponding to the Hessian matrix He(x) of the measurement vector are:

In the formula, e _i and e _j are the real part of the complex voltage of node i and node j respectively; G _ij and B _ij are the elements of the admittance matrix; g, b, y _c are the parameters in the line Π-shaped model; K is the transformer The non-standard transformation ratio; b _T is the susceptance of the standard side of the transformer. The Π-shaped equivalent circuit of the line and transformer is shown in Figure 2.

In order to verify the effectiveness of the present invention, it is compared with different state estimation methods. Among them, WLS state estimation is called method 1, WLS+BD is called method 2, predictor-collector PDIPM (PCPDIPM) is called method 3, and multiple prediction-correction interior point method (multiple PCPDIPM, MPCPDIPM) is called is method 4, and the linear primal-dual interior point method is called method 5. In order to be more in line with the actual physical meaning, the state variable estimation results of each method are given in the form of polar coordinates, and the transformation method is:

Introduce two embodiment of the present invention below:

The present invention adopts the standard calculation example of the IEEE-14 node shown in FIG. 3 and performs calculations in three cases. 1) The system contains bad data, and the normal measurements take the real value, which is used to analyze the ability of different methods to identify bad data; 2) The system does not contain bad data, and all measurements have random errors, which are used to analyze the The state estimation performance of different methods under the normal measurement that conforms to the normal distribution; 3) The system contains bad data, and the measurement also has random errors, that is, the synthesis of cases 1) and 2), which is also the most common one in practice mode. Table 1 shows the bad data identification results of different estimation methods in case 1), and Table 2 shows the comparison of state variable estimation results with different methods.

Table 1 State estimation results under different methods 1)

Table 2 State variable error under different methods 3)

It can be seen from Tables 1 and 2 that method 2 did not identify bad data No. 5 and No. 7, while methods 3 and 4 all measured and estimated errors were less than 1%, and bad data were well identified, meeting engineering requirements. Although method 5 can also identify bad data, the estimation error is significantly larger than methods 3 and 4. The above tests show that, compared with the linear primal-dual interior point method, the present invention (method 4) greatly improves the solution accuracy and has stronger tolerance to errors.

Calculation example two:

In order to verify the computational efficiency of this method, 6 systems from 14 to 3012 nodes were tested, and the dual gap convergence accuracy was set to 10 ^-6 . Table 3 shows the number of iterations using three different interior point methods under different systems, and the percentage of iterations reduced compared with PDIPM. Table 4 gives the calculation time of the corresponding method and the percentage of calculation time reduction compared with PDIPM.

Table 3 Comparison of iteration times of different methods

Note: Wp-2383, 3012 in the table is the Polish power grid system; the part in bold is the best result in the table, the same below

Table 4 CPU calculation time comparison of different methods

It can be seen from Tables 3 and 4 that compared with PDIPM, method 3 can reduce the number of iterations by at least 13.64%, and the calculation time can be reduced by 13.84%. On the basis of method 3, the number of iterations can be reduced by 20.83% through multiple corrections. , the calculation time can be reduced by more than 16.28%, and it has a faster convergence speed.

Further, Figure 4 shows the iteration step size of different interior point methods for IEEE-57 node systems, and Figure 5 shows the convergence process of the dual gap. It can be seen from Figure 4 that the step size of PDIPM is relatively large in the first three iterations, and then gradually decreases. Method 3 increases the step size compared with PDIPM due to the addition of correction steps, but in the initial stage of iteration, the iteration step size changes instead. Small, this is because the correction direction does not point to the optimal direction, and method 4 uses weighted multi-prediction, and uses the 2-stage method to find the optimal correction weight to ensure that the correction direction points to the optimal direction, only needing 11 iterations Convergence, the iteration step size always maintains a large value. It can be seen from Figure 5 that PDIPM converges after 16 iterations. Method 3 has a faster convergence speed of the dual gap, and it only needs 13 iterations to converge, which is 18.75% less than the number of PDIPM iterations. During the process, the correction is performed several times, so as to obtain the optimal iteration step size, the dual gap speed is significantly improved, and it can converge after only 11 iterations, which is 31.25% lower than that of PDIPM.

Figure 6 shows the search process for a certain iteration ω to find the optimal weight in the interval [α _p α _d ,1], where α _p =0.5524, α _d =0.7606. In order to obtain the optimal step size of the original and dual variables, this paper allows the original and dual variables to take different optimal weights. Figure 6 shows the search process in the first stage, find out the optimal subinterval of α _p [0.47818, 0.59414], and the optimal subinterval of α _d is [0.82606, 0.94202]. Most of the optimal step size and weight have an approximate parabolic relationship, so in the first stage of search, the number of set points does not require too much, as long as the obtained sub-interval can be guaranteed to be globally optimal. Setting the number of search points to 4 in the sub-interval of the second stage often results in large errors. In this embodiment, it is set to 10 to ensure that the search results are closer to the real optimal weight. The final optimization results of this time are α _p =0.88857, α _d =0.80510, which have been increased to a certain extent compared with the original iteration step size (ie ω=1). The experimental results of the second example show the validity of the two-stage linear search method used in this method to find the optimal weight.

The WLAV robust state estimation method based on the multi-prediction-correction interior point method shown in the above Figures 1-6 is a specific embodiment of the present invention, which has reflected the substantive features and progress of the present invention, and can be used according to actual needs. Under the enlightenment of the present invention, equivalent modifications to it in terms of shape, structure, etc., all fall within the scope of protection of this proposal.

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.