CN105846429A - Power flow optimization method for electric power system based on subarea division and class extension variable relaxation - Google Patents

Abstract

The invention discloses a power flow optimization method for an electric power system based on subarea division and class extension variable relaxation. The method comprises the following steps of decomposing the electric power system to obtain a plurality of subsystems which are mutually independent; coordinating the subsystems to obtain a decomposition coordination model of the electric power system based on subarea division; adopting an auxiliary problem principle to carry out distributed parallel optimal power flow computation on the decomposition coordination model of the electric power system; and adopting a class extension variable interior point method to solve the subsystems respectively to obtain an optimal solution of power flow of the electric power system. The method conforms to the characteristics of scattered data and structural distribution of the large-scale interconnected electric power systems; the dimensionality of the subsystems is effectively reduced; the optimal computation of each subarea only needs data and boundary node data of the area and does not need a detailed mathematic model outside the area, so that the bottleneck of data transmission is avoided and the optimization speed is accelerated; through the class extension variable interior point method adopted by a subsystem optimization problem, a solution space of the optimization problem is expanded and the convergence is improved.

Description

Based on the electric system tide optimization method that sub area division and class expansion variable are lax

Technical field

The present invention relates to Optimal Technology of Power Systems field, especially a kind of lax based on sub area division and class expansion variable Electric system tide optimization method.

Background technology

The basic ideas of composition decomposition parallel optimization algorithm based on sub area division are to manage to have large scale optimization problem Effect is decomposed into a series of subproblem, and carries out distributed solving and coordinating and optimizing, be applied to the earliest daily trading planning optimization and point During cloth optimal load flow calculates.Solve and be based primarily upon 3 kinds of considerations: the independence 1. retaining each region is asked to carry out total optimization The calculating of topic；The most in a distributed fashion or parallel mode solves optimal power flow problems, improve and calculate speed；3. subdispatcher's energy is made Enough dispatch electrical network in phase, and without knowing the service data in other regions.

In recent years, composition decomposition parallel optimization algorithm based on sub area division has had certain research, parallel except using Outside the intelligent algorithm of processing mode, its implementation the most substantially has three classes: first kind method is based on optimization problem Karush Kuhn Tucker (KKT) condition obtains the system update equation of the whole network, then uses diagonal angle edged model, approximation newton This system update equation is decoupled by direction equal matrix theory.But the coupling that such method computational efficiency is limited by between subnet Degree, the number of Region Decomposition, and there is the cooperation layer burden problem such as heavily, it is more suitable for the scene with ground Distributed Parallel Computing. Equations of The Second Kind method, by external network is carried out equivalence, carries out independently solving of each partition zone optimizing, and introduces outer layer coordination link Realize the optimization of the whole network.Such method need when carrying out local computing the operational factor of adjoining subnetworks to carry out external network equivalent, Have impact on the independence that subnet optimization calculates, and the calculating of this equivalent parameters and update relatively cumbersome.3rd class method is based on drawing Ge Lang Relaxation Theory, being relaxed to object function by the coupling constraint between subnet and use certain decomposition algorithm to realize it can Divide property.As a example by decomposition algorithm based on Auxiliary Problem Principle (Auxiliary Problem Principle, APP), these classics Former saddle-point problem, by structure auxiliary function, is divided into the extreme-value problem that can independently solve, efficiently solves concentration by method The problem that formula parallel optimization algorithm brings.

Nonlinear interior-point method has preferable convergence and calculates speed faster, has become as and solves electric power system optimization The powerful of problem, but its solution space unsatisfactory.Class expansion variable interior point method is by the inequality in same kind Constraint equation increases identical class expansion variable, and forces class expansion variable to tend to 0 realization with penalty factor in object function The reduction of constraints, expands the solution space of optimization problem, improves convergence.

Summary of the invention

The present invention is directed to the deficiencies in the prior art, propose a kind of power train lax based on sub area division and class expansion variable System tide optimization method, meets the dispersion of Large-Scale Interconnected electric power system data, the feature of structure distribution；Significantly reduce subsystem Dimension so that the optimization problem scale of each subsystem is less, improves the calculating speed of network optimum trend；Each subregion Optimize to calculate and have only to data and the boundary node data of one's respective area, it is not necessary to mathematical model detailed outside region, it is to avoid number According to the bottleneck of transmission, accelerate optimal speed；The class expansion variable interior point method that subsystem optimization problem uses expands optimization and asks The solution space of topic, improves convergence.

In order to realize foregoing invention purpose, the present invention provides techniques below scheme:

A kind of electric system tide optimization method lax based on sub area division and class expansion variable, described method include as Lower step:

Step 1. decomposes power system, obtains multiple separate subsystem；

Step 2. coordinates each described subsystem, obtains power system decomposition-coordination model based on sub area division；

Step 3. uses Auxiliary Problem Principle that described power system decomposition-coordination model carries out distributed parallel optimum tide Stream calculation；

Step 4. uses class expansion variable interior point method to solve each described subsystem respectively, obtains the trend of described power system Optimal solution.

Further, before described step 1, including:

Determine the optimal load flow mathematical model of described power system:

$f\left(x\right)=\left\{\begin{array}{c}\underset{x}{\min }f\left(x\right)\\ \begin{array}{cc}s.t.& h\left(x\right)=0\\ & \underset{\‾}{g}\≤g\left(x\right)\≤\stackrel{\‾}{g}\end{array}\end{array}\right.---\left(1\right)$

In formula (1): f (x) is object function, x is the column vector that independent variable is constituted；H (x) and g (x) is respectively equality constraint With inequality constraints condition；gWithIt is respectively lower limit and the higher limit of inequality constraints condition g (x).

Further, described step 1 includes:

1-1., according to the geographical distribution of described power system, determines boundary node at interconnection；

1-2. judges whether described boundary node is active node；

The most then enter step 1-3；

If it is not, then enter step 1-4；

1-3. increases virtual passive bus on described interconnection, using passive bus with it as new boundary node, Enter step 1-4；

Described power system is decomposed at each described boundary node by 1-4., obtains multiple separate described subsystem System.

Further, described step 2 includes:

2-1. using described boundary node parameter as coordination variable；

Coordinated by the described coordination variable of exchange between each described subsystem of 2-2., obtain based on sub area division Power system decomposition-coordination model.

Further, described step 2-2 includes:

A. coordinated by the described coordination variable of exchange between each described subsystem, obtain institute based on sub area division State

Object function minf (x) of power system decomposition-coordination model:

Minf (x)=f₁(x₁)+f₂(x₂)….+f_n(x_n) (2)

In formula (2), x₁、x₂...x_nIt is respectively described subsystems；f₁(x₁)、f₂(x₂)...f_n(x_n) it is respectively each institute State the object function of subsystem；

B. the boundary constraint of described power system decomposition-coordination model is determined.

Further, described step 3 includes:

3-1., according to the kernel function in described Auxiliary Problem Principle, constructs the auxiliary of described power system decomposition-coordination model Function；

The primal problem that described power system decomposition-coordination model solves is converted into the saddle solving described auxiliary problem (AP) by 3-2. Point problem；

Described auxiliary problem (AP) is iterated by 3-3. according to the threshold value arranged, until loop ends, obtains described power train The distributed parallel optimal load flow result of calculation of system decomposition-coordination model.

Further, described step 4 includes:

4-1. increases identical class expansion variable in the generic inequality constraints equation in each described subsystem, and Described class expansion variable is respectively provided with by upper and lower bound；

4-2. increases the expansion variable quadratic term of band penalty factor in the object function of each described subsystem, expands inequality The solution space of constraint；

Logarithm barrier function inequality constraints as equality constraint, is processed equation about with method of Lagrange multipliers by 4-3. Bundle, obtains the Lagrangian after introducing class expansion variable；

4-4. derivation KKT condition, obtains the update equation group of the Lagrangian after described introducing class expansion variable；

4-5. solves described update equation, obtains the trend optimal solution of described power system.

Further, described step 4-5 includes:

(1) described equation group is initialized: putting iterations is 0, given independent variable, relaxation factor, Lagrange multiplier, class The initial value of expansion variable；

(2) calculating target function and the Jacobian matrix of constraint and extra large gloomy matrix, forms linear system；

(3) solve described linear system, obtain the correction of each iteration；

(4) the affine step-length of former variable and the affine step-length of dual variable are calculated；

(5) each variable and class expansion variable are updated；

(6) according to the convergence set, it is judged that whether optimal load flow restrains；

The most then stop calculating, obtain the trend optimal solution of described power system；

If it is not, then return step (2).

The present invention proposes a kind of power system Auxiliary Problem Principle (APP) and class expansion variable interior point method combined also Row optimal load flow algorithm, utilizes APP to set up decomposition-coordination model based on sub area division, utilizes class expansion variable interior point method Process the optimization problem of each subsystem, thus improve Algorithm Convergence while reducing calculating dimension, shortening the calculating time. The simulation analysis of IEEE 118 meshed network demonstrates reasonability and the effectiveness of this algorithm.It has the beneficial effect that

(1) method that the present invention proposes meets the dispersion of Large-Scale Interconnected electric power system data, the feature of structure distribution；

(2) the method significantly reduces the dimension of subsystem so that the optimization problem scale of each subsystem is less, improves The calculating speed of network optimum trend；

(3) optimizing of each subregion calculates data and the boundary node data having only to one's respective area, it is not necessary to outside region in detail Thin mathematical model, it is to avoid the bottleneck of data transmission, accelerates optimal speed；

(4) the class expansion variable interior point method that subsystem optimization problem uses expands the solution space of optimization problem, improves Convergence.

Accompanying drawing explanation

Fig. 1 is the schematic flow sheet of the trend parallel optimization method of the power system based on sub area division of the present invention；

Fig. 2 be the present invention as a example by two subregions in power system decompose before structural representation；

Fig. 3 be the present invention as a example by two subregions in power system decompose after structural representation；

Fig. 4 is the calculation flow chart of the concrete application examples of the present invention.

Detailed description of the invention

Describing the present invention below in conjunction with the accompanying drawings, the description of this part is only exemplary and explanatory, should not Protection scope of the present invention is had any restriction effect.

As it is shown in figure 1, the present invention provides a kind of electric power system tide lax based on sub area division and class expansion variable excellent Change method, comprises the steps:

Step 1. decomposes power system, obtains multiple separate subsystem；

Step 2. coordinates subsystems, obtains power system decomposition-coordination model based on sub area division；

Step 3. uses Auxiliary Problem Principle that power system decomposition-coordination model carries out distributed parallel optimal load flow meter Calculate；

Step 4. uses class expansion variable interior point method to solve each subsystem respectively, obtains the trend optimal solution of power system.

Wherein, before step 1, including:

Determine the optimal load flow mathematical model of power system:

$f\left(x\right)=\left\{\begin{array}{c}\underset{x}{\min }f\left(x\right)\\ \begin{array}{cc}s.t.& h\left(x\right)=0\\ & \underset{\‾}{g}\≤g\left(x\right)\≤\stackrel{\‾}{g}\end{array}\end{array}\right.---\left(1\right)$

In formula (1): f (x) is object function, x is the column vector that independent variable is constituted；H (x) and g (x) is respectively equality constraint With inequality constraints condition；gWithIt is respectively lower limit and the higher limit of inequality constraints condition g (x).

Wherein, step 1 includes:

1-1., according to the geographical distribution of power system, determines boundary node at interconnection；

1-2. judges whether boundary node is active node；

The most then enter step 1-3；

If it is not, then enter step 1-4；

1-3. increases virtual passive bus on interconnection, using passive bus with it as new boundary node, enters Step 1-4；

Power system is decomposed at each boundary node by 1-4., obtains multiple separate subsystem.

Wherein, step 2 includes:

2-1. using boundary node parameter as coordination variable；

Coordinated by exchange coordination variable between 2-2. subsystems, obtain power system based on sub area division Decomposition-coordination model.

Wherein, step 2-2 includes:

A. coordinated by exchange coordination variable between subsystems, obtain power system based on sub area division and divide Object function minf (x) of solution Coordination Model:

Minf (x)=f₁(x₁)+f₂(x₂)….+f_n(x_n) (2)

In formula (2), x1, x2...xn are respectively subsystems；F1 (x1), f 2 (x2) ... f n (xn) is respectively each son The object function of system；

B. the boundary constraint of power system decomposition-coordination model is determined.

Wherein, step 3 includes:

3-1. is according to the kernel function in Auxiliary Problem Principle, the auxiliary function of structure power system decomposition-coordination model；

The primal problem that power system decomposition-coordination model solves is converted into the saddle-point problem solving auxiliary problem (AP) by 3-2.；

Auxiliary problem (AP) is iterated by 3-3. according to the threshold value arranged, until loop ends, obtains power system and decomposes association The distributed parallel optimal load flow result of calculation of mode transfer type.

Wherein, step 4 includes:

4-1. increases identical class expansion variable in the generic inequality constraints equation in each subsystem, and class expands Exhibition variable is respectively provided with by upper and lower bound；

4-2. increases the expansion variable quadratic term of band penalty factor in the object function of each subsystem, expands inequality constraints Solution space；

Logarithm barrier function inequality constraints as equality constraint, is processed equation about with method of Lagrange multipliers by 4-3. Bundle, obtains the Lagrangian after introducing class expansion variable；

4-4. derivation KKT condition, obtains the update equation group of the Lagrangian after introducing class expansion variable；

4-5. solves update equation, obtains the trend optimal solution of power system.

Wherein, step 4-5 includes:

(1) initialized equations group: putting iterations is 0, given independent variable, relaxation factor, Lagrange multiplier, class extension The initial value of variable；

(2) calculating target function and the Jacobian matrix of constraint and extra large gloomy matrix, forms linear system；

(3) solve linear system, obtain the correction of each iteration；

(4) the affine step-length of former variable and the affine step-length of dual variable are calculated；

(5) each variable and class expansion variable are updated；

(6) according to the convergence set, it is judged that whether optimal load flow restrains；

The most then stop calculating, obtain the trend optimal solution of power system；

If it is not, then return step (2).

The present invention introduces subregion and the parallel optimization model thereof of bulk power grid as a example by two subregions, as follows:

1, optimal load flow mathematical model

Optimal load flow calculates mathematically can be attributed to constrained nonlinear programming problem, and its mathematical model is generally:

$\left\{\begin{array}{c}\underset{x}{\min }f\left(x\right)\\ \begin{array}{cc}s.t.& h\left(x\right)=0\\ & \underset{\‾}{g}\≤g\left(x\right)\≤\stackrel{\‾}{g}\end{array}\end{array}\right.---\left(1\right)$

In formula: f (x) is object function, the generating expense that the object function of optimal load flow is usually system total (sends meritorious Idle cost) minimum, or the circuit active power loss minimum that system is total；X is the column vector that independent variable is constituted, and represents generating Machine is meritorious, idle exerts oneself and control variable or the state variable such as node voltage amplitude angle；H (x) and g (x) is expressed as Equality constraint and inequality constraints condition.

2, sub area division

The complexity that optimal load flow calculates is closely related with electrical network scale, for reducing the difficulty of problem solving, Ke Yikao Consider the actual geographic according to electrical network to be distributed, at some interconnection, whole electrical network is decomposed into multiple relatively independent region, point Do not calculate.The present invention uses composition decomposition method, is i.e. become by one big system Direct Resolution by " duplication " boundary node Multiple subsystems, subsystems is relatively independent, with the parameter of boundary node as coordination variable, the solution of each subregion only and oneself Built-in variable relevant with Boundary Variables, only rely on exchange boundary node data between subregion and coordinate.Below with two subregions As a example by introduce subregion and the parallel optimization process thereof of bulk power grid:

As in figure 2 it is shown, the two of electrical network S subsystem X₁, X₂By border X_bBeing connected, the branch road being connected is referred to as interconnection and props up Road.As it is shown on figure 3, " duplication " boundary node X after decomposition_b, carry out network separation.Use X_bi(i=1,2) each subsystem is represented " replicate " boundary node, C_iRepresent each subregion after separating, represent the set of all subregions, i.e. E={C with E_i(i=1, 2)。

Each boundary node can be regarded as a virtual electromotor node, x_biIt it is the i-th subregion virtual synchronous generator joint The parametric variable of point.Virtual synchronous generator compensate for because the trend causing the Direct Resolution of electrical network is uneven, also makes simultaneously Obtain C_iDefine a relatively independent subsystem.

System after decomposition to be made is equivalent with original system, should specify x_b1And x_b2It is being electrically same point, is having identical Electric parameters, i.e.WithEqual, i.e. x_b1=x_b2.Wherein,WithRepresent the virtual synchronous generator injecting power of boundary node；WithRepresent voltage magnitude and phase Angle is when carrying out subregion, if the boundary node of actual electric network is active node, first should increase virtual nothing on interconnection branch road Source node, replaces original boundary node, carries out subregion the most again.Meanwhile, when the constraint in subregion is relevant with district's exogenousd variables, The node at this variable place can be regarded as boundary node, make subregion constraint and in this subregion and boundary node change Measure relevant, with district outside variable unrelated.

So, former object function f (x) can be write as the object function sum of each subregion, is designated as f₁(x₁)+f₂(x₂), x₁With x₂Represent the inner parameter variable in 2 regions respectively.The equation of the whole network and inequality constraints are decomposed into the constraint in each region.For Boundary node, owing to should have identical electric parameters, it should increase an equality constraint on boundary node between subsystem.

After using above-mentioned decomposition-coordination approach, obtain 2 subregion equivalence optimal load flow decomposition-coordination models as follows:

Object function:

Minf (x)=f₁(x₁)+f₂(x₂) (2)

Boundary constraint:

h₁(x₁)=0, h₂(x₂)=0

${\underset{\‾}{g}}_1\≤g_1\left(x_1\right)\≤{\stackrel{\‾}{g}}_1,{\underset{\‾}{g}}_2\≤g_2\left(x_2\right)\≤{\stackrel{\‾}{g}}_2---\left(3\right)$

θ (x)=x_b1-x_b2=0

In formula, x=[x₁,x_b1,x₂,x_b2]^T, f_i(x_i) it is subregion C_iObject function, θ (x) is global restriction.

For the power system of multi partition, the main thought that its decomposition-coordination model is set up may be used to 2 partition systems Explanation.

3, parallel optimal load flow computational methods based on Auxiliary Problem Principle

Do not consider the constraints of subsystem internal, only consider the boundary point constraint produced due to subregion, now according to formula (2), formula (3), the Lagrangian of optimization problem can be written as:

L (x, λ)=f₁(x₁)+f₂(x₂)+<λ, θ (x)>(4)

In order to improve convergence, use Augmented Lagrange method, on former Lagrangian, i.e. increase by one Individual quadratic term:

$L(x,\mathrm{\&lambda;})=f_1\left(x_1\right)+f_2\left(x_2\right)+<\mathrm{\&lambda;},\mathrm{\&theta;}\left(x\right)>+\frac{c}{2}<\mathrm{\&theta;}\left(x\right),\mathrm{\&theta;}\left(x\right)>---\left(5\right)$

In formula, λ is Lagrange multiplier vector, and c is penalty factor, is a constant.

It can be seen that introduce quadratic term can't affect final result of calculation, when last iteration convergence, between subsystem Boundary node electric parameters tends to consistent, and now the value of quadratic term tends to 0.Owing to θ (x) is the letter being correlated with 2 area variables Count, and the quadratic term itself added is nondecomposable, thus destroy the separability between two subsystems so that optimize meter Calculation cannot directly complete in 2 subregions independently of one another.In this case, introduce Auxiliary Problem Principle (APP) to solve This problem.

3.1 Auxiliary Problem Principles (APP)

The thought of APP is, by structure kernel function, will solve the saddle-point problem of L (x, λ) and change into and solve series of iterations The auxiliary saddle-point problem of form, and the auxiliary problem (AP) of iteration is decomposable problem each time.

As a example by minimizing, aided algorithm principle is described.Assume E₁X () can be micro-, and Lipschitz is continuous, function E₂ X () not necessarily can be micro-, for primal problem minE₁(x)+E₂X () constructs an auxiliary problem (AP): minG (x)+ε E₂X (), if may certify that There is x^*Make

G′(x^*)=ε E '₁(x^*) (6)

Set up, then x^*Also the solution of primal problem is become, therefore, structure auxiliary function G (x, v), wherein, Ψ (x) is kernel function:

G (x, v)=Ψ (x)+< ε E '₁(v)-Ψ′(v),x> (7)

Then primal problem can be converted into the saddle-point problem solving following auxiliary problem (AP), it may be assumed that

x^k+1=arg min Ψ (x)+< ε E '₁(x^k)-Ψ′(x^k),x>+εE₂(x) (8)

In formula, argmin represents makes auxiliary problem (AP) Ψ (x)+< ε E '₁(x^k)-Ψ′(x^k),x>+εE₂When () obtains minima x The value of x.When | | x^k-x^k+1| | stop less than when specifying threshold value, otherwise continue iteration.

3.2 parallel optimal load flow computational methods based on APP

For primal problem L (x, λ), order

L (x, λ)=f (x)+J (x, λ) (9)

F (x)=f₁(x₁)+f₂(x₂) (10)

$J(x,\mathrm{\&lambda;})=<\mathrm{\&lambda;},\mathrm{\&theta;}\left(x\right)>+\frac{c}{2}<\mathrm{\&theta;}\left(x\right),\mathrm{\&theta;}\left(x\right)>---\left(11\right)$

Wherein, J (x, λ) can be micro-, and f (x) not necessarily can be micro-.Structure auxiliary function is:

$G^{\left(x^k,{\mathrm{\&lambda;}}^k\right)}\left(x,\mathrm{\&lambda;}\right)=\mathrm{\&Psi;}\left(x,\mathrm{\&lambda;}\right)+<{\mathrm{\&epsiv;J}}^{\&prime;}\left(x^k,{\mathrm{\&lambda;}}^k\right)-{\mathrm{\&Psi;}}^{\&prime;}\left(x^k,{\mathrm{\&lambda;}}^k\right),\left(\begin{array}{c}x\\ \mathrm{\&lambda;}\end{array}\right)>---\left(12\right)$

In formula, ε ＞ 0；Ψ (x, λ) is the kernel function needing structure, is configured to:

$\mathrm{\&Psi;}(x,\mathrm{\&lambda;})=\mathrm{\&epsiv;}\frac{b}{2}<x,x>-\frac{1}{2}<\mathrm{\&lambda;},\mathrm{\&lambda;}>---\left(13\right)$

ε takes ε=ε during minimizing and maximum problem respectively₁=1 and ε=ε₂；B is positive number.

Two-layer algorithm model according to Auxiliary Problem Principle, in kth time iteration, the saddle-point problem solving L (x, λ) is suitable In alternately solving following minimum and maximum problem:

Solution minimum problem:

$x^{k+1}=\arg min\&lsqb;f\left(x\right)+\frac{b}{2}<x,x>-<{\mathrm{bx}}^k,x>+<{\mathrm{\&lambda;}}^k+c\mathrm{\&theta;}\left(x^k\right),\mathrm{\&theta;}\left(x\right)>\&rsqb;---\left(14\right)$

Solution maximum problem:

λ^k+1=λ^k+ε₂θ(x^k+1) (15)

It is obvious that the crossed product that will not occur again in formula (14) between different variable, can be further divided.For 2 points Sound zone system, substitutes into formula (14), formula (15) by f (x), θ (x) and x, makes β=1+b, arrange:

$\begin{array}{c}{x}^{k+1}=\arg min\&lsqb;\frac{\mathrm{\&beta;}}{2}{x_1}^T{x}_1+\frac{\mathrm{\&beta;}}{2}{x_2}^T{x}_2+\frac{\mathrm{\&beta;}}{2}{x_{b1}}^T{x}_{b1}+\frac{\mathrm{\&beta;}}{2}{x_{b2}}^T{x}_{b2}-\mathrm{\&beta;}{\left(x_1^k\right)}^T{x}_1-\mathrm{\&beta;}{\left(x_2^k\right)}^T{x}_2-\\ \mathrm{\&beta;}{\left(x_{b1}^k\right)}^T{x}_{b1}-\mathrm{\&beta;}{\left(x_{b2}^k\right)}^T{x}_{b2}+<{\mathrm{\&lambda;}}^k+c\left(x_{b1}^k-x_{b2}^k\right),\left(x_{b1}-x_{b2}\right)>+\\ f_1\left(x_1\right)+f_2\left(x_2\right)\&rsqb;\end{array}---\left(16\right)$

${\mathrm{\&lambda;}}^{k+1}={\mathrm{\&lambda;}}^k+{\mathrm{\&epsiv;}}_2(x_{b1}^k-x_{b2}^k)---\left(17\right)$

Being not difficult to find out, formula (16) can be solved formula (18) by each Paralleled and realize:

$\begin{array}{c}\left(x_i^{k+1},x_{bi}^{k+1}\right)=\arg \underset{\left(x_i,x_{bi}\right)\&Element;C_i}{\min }\{f_i\left(x_i\right)+\frac{\mathrm{\&beta;}}{2}{x_i}^T{x}_i+\frac{\mathrm{\&beta;}}{2}{x_{bi}}^T{x}_{bi}-\mathrm{\&beta;}{\left(x_i^k\right)}^T{x}_i-\mathrm{\&beta;}{\left(x_{bi}^k\right)}^T{x}_{bi}+\\ {\&lsqb;{\mathrm{\&lambda;}}^k+c\left(x_{b1}^k-x_{b2}^k\right)\&rsqb;}^T{\mathrm{qx}}_{bi}\}\end{array}---\left(18\right)$

In formula, i=1,2；During i=1, q=1；During i=2, q=-1.

4, the class expansion variable interior point method of subsystem optimization problem solves

Above the optimization problem of a macroreticular is decomposed into the optimization problem of 2 subsystems, and for subsystems Define respective object function and constraints, with subsystem C₁As a example by, its sub-optimization problem and being constrained to:

$\left\{\begin{array}{c}min{f}_1(x_1,x_{b1})=\{f_1\left(x_1\right)+\frac{\mathrm{\&beta;}}{2}{x_1}^T{x}_1+\frac{\mathrm{\&beta;}}{2}{x_{b1}}^T{x}_{b1}-\mathrm{\&beta;}{\left(x_1^k\right)}^T{x}_1-\mathrm{\&beta;}{\left(x_{b1}^k\right)}^T{x}_{b1}+\\ {\&lsqb;{\mathrm{\&lambda;}}^k+c(x_{b1}^k-x_{b2}^k)\&rsqb;}^T{x}_{b1}\}\\ \begin{array}{cc}s.t.& h_1(x_1,x_{b1})=0\\ & {\underset{\&OverBar;}{g}}_1\&le;g_1(x_1,x_{b1})\&le;{\stackrel{\&OverBar;}{g}}_1\end{array}\end{array}\right.---\left(19\right)$

In order to improve convergence, class expansion variable interior point method is used to solve the optimization problem of subsystem here.Class Expansion variable interior point method passes through to increase identical class expansion variable in the inequality constraints equation of kind, and at target letter Number forces the quadratic sum of class expansion variable to tend to 0 with penalty factor, thus expands the solution space of inequality constraints.When former problem Time feasible, this algorithm can converge to the optimal solution of former problem；When former problem is without solving, can automatically arrive in bigger feasible zone Optimizing, quickly obtains approximate solution, can readily obtain the measure of adjustment from result of calculation.

The introducing of 4.1 class expansion variables

Inequality constraints category in formula (19) is increased expansion variable, it is contemplated that same class constraint is from upper and lower bound pine The demand relaxed is generally independent of one another, and class expansion variable need to be respectively provided with by upper and lower bound.For make the solution of optimization problem as far as possible to The limit value of former setting is drawn close, and increases the expansion variable quadratic term of band penalty factor in object function.Then formula (19) becomes:

$\left\{\begin{array}{c}\min f_1\left(x_1,x_{b1}\right)+\underset{k=1}{\overset{n}{\&Sigma;}}{\mathrm{\&xi;}}_k\left({\mathrm{\&alpha;}}_k^2+{\mathrm{\&beta;}}_k^2\right)\\ \begin{array}{cc}s.t& h_1\left(x_1,x_{b1}\right)=0\\ & {\underset{\&OverBar;}{g}}_1-T\mathrm{\&alpha;}\&le;g_1\left(x_1,x_{b1}\right)\&le;{\stackrel{\&OverBar;}{g}}_1+T\mathrm{\&beta;}\end{array}\end{array}\right.---\left(20\right)$

In formula (20), n is the kind of the inequality constraints arranging expansion variable, and α, β are class expansion variable column vector, respectively Representing lower expansion variable and upper expansion variable, dimension is n；α_k、β_kRepresent the kth element of α, β respectively；ξ_kExpand for kth class The penalty factor of exhibition variable quadratic term, T is inequality constraints and the relational matrix of class expansion variable, ties up matrix for r × n, and form is such as Shown in lower:

$T=\left[\begin{array}{cccc}1& 0& \mathrm{...}& 0\\ 1& 0& \mathrm{...}& 0\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ 0& 0& \mathrm{...}& 1\\ 0& 0& \mathrm{...}& 0\end{array}\right]$

If formula (20) adds in i-th inequality constraints jth class expansion variable, then T_i,jIt is 1, its of T the i-th row Its element is 0.

4.2 subsystem optimization problems solve scheme

According to the core concept of interior point method, it is equality constraint with logarithm barrier function inequality constraints, with Lagrange Multiplier method processes equality constraint, and the Lagrangian after introducing class expansion variable corresponding with formula (20) is:

$\begin{array}{c}{L}_{1g}=f_1\left(x_1,x_{b1}\right)+\underset{k=1}{\overset{n}{\&Sigma;}}{\mathrm{\&xi;}}_k\left({\mathrm{\&alpha;}}_k^2+{\mathrm{\&beta;}}_k^2\right)-y^T{h}_1\left(x_1,x_{b1}\right)-z^T\left(g_1\left(x_1,x_{b1}\right)-l-{\underset{\&OverBar;}{g}}_1+Ta\right)-\mathrm{\&mu;}\underset{j=1}{\overset{r}{\&Sigma;}}\ln l_j-\\ w^T\left(g_1\left(x_1,x_{b1}\right)+u-{\stackrel{\&OverBar;}{g}}_1-T\mathrm{\&beta;}\right)-\mathrm{\&mu;}\underset{j=1}{\overset{r}{\&Sigma;}}\ln u_j\end{array}---\left(21\right)$

Derivation KKT condition, can obtain update equation group is:

$\left[\begin{array}{cccc}{H}^{\&prime;}& {\&dtri;}_x{h}_1\left(x_1,x_{b1}\right)& F1& F2\\ {\&dtri;}_x^T{h}_1\left(x_1,x_{b1}\right)& 0& 0& 0\\ -F{1}^T& 0& F3& 0\\ -F{2}^T& 0& 0& F4\end{array}\right]\left[\begin{array}{c}\mathrm{\&Delta;}x\\ \mathrm{\&Delta;}y\\ \mathrm{\&Delta;}\mathrm{\&alpha;}\\ \mathrm{\&Delta;}\mathrm{\&beta;}\end{array}\right]=\left[\begin{array}{c}{L_x}^{\&prime;}\\ -L_y\\ -{L_{\mathrm{\&alpha;}}}^{\&prime;}\\ -{L_{\mathrm{\&beta;}}}^{\&prime;}\end{array}\right]---\left(22\right)$

$\left[\begin{array}{cc}L& Z\\ I& 0\end{array}\right]\left[\begin{array}{c}\mathrm{\&Delta;}z\\ \mathrm{\&Delta;}l\end{array}\right]=\left[\begin{array}{c}-L_1\\ -{L_z}^{\&prime;}-{\mathrm{ZL}}^{-1}\left({\&dtri;}_x^T{g}_1\right(x_1,x_{b1})\mathrm{\&Delta;}x+T\mathrm{\&Delta;}\mathrm{\&alpha;})\end{array}\right]---\left(23\right)$

$\left[\begin{array}{cc}U& W\\ I& 0\end{array}\right]\left[\begin{array}{c}\mathrm{\&Delta;}w\\ \mathrm{\&Delta;}u\end{array}\right]=\left[\begin{array}{c}-L_u\\ {L_w}^{\&prime;}+{\mathrm{WU}}^{-1}\left({\&dtri;}_x^T{g}_1\right(x_1,x_{b1})\mathrm{\&Delta;}x-T\mathrm{\&Delta;}\mathrm{\&beta;})\end{array}\right]---\left(24\right)$

In equation group (22)～(24), have:

$\begin{array}{c}{H}^{\&prime;}=-\&lsqb;{\&dtri;}_x^2{f}_1(x_1,x_{b1})-{\&dtri;}_x^2{h}_1\left(x_1,x_{b1}\right)y-{\&dtri;}_x^2{g}_1\left(x_1,x_{b1}\right)\left(z+w\right)\&rsqb;+\\ {\&dtri;}_x{g}_1\left(x_1,x_{b1}\right)\left({\mathrm{WU}}^{-1}-{\mathrm{ZL}}^{-1}\right){\&dtri;}_x^T{g}_1\left(x_1,x_{b1}\right)\end{array}---\left(25\right)$

F1=-_xg₁(x₁,x_b1)ZL^-1T (26)

F2=-_xg₁(x₁,x_b1)WU^-1T (27)

F3=2 ξ+T^TZL^-1T (28)

F4=2 ξ-T^TWU^-1T (29)

L_x'=L_x+▽_xg₁(x₁,x_b1)(L_Z'-L_W') (30)

L_α'=L_α+T^TL_z' (31)

L_β'=L_β+T^TL_w' (32)

L_z'=L^-1(ZL_z+L_l) (33)

L_w'=U^-1(WL_w-L_u) (34)

Solving equations step is as follows:

(1) initialize: putting iterations is 0, given independent variable, relaxation factor, Lagrange multiplier, class expansion variable Initial value；

(2) calculating target function and the Jacobian matrix of constraint and extra large gloomy matrix, form the line shown in formula (22)～(24) Sexual system；

(3) solve formula (22)～(24), obtain the correction of each iteration；

(4) affine step-length λ of former variable is calculated_afpStep-length λ affine with dual variable_afd:

$\left\{\begin{array}{c}{\mathrm{\&lambda;}}_{afp}=\min (0.9995\underset{i}{\min }(\frac{-l_i}{{\mathrm{\&Delta;l}}_{afi}},\frac{-u_i}{{\mathrm{\&Delta;u}}_{afi}}),1)\\ {\mathrm{\&Delta;l}}_{afi}<0,{\mathrm{\&Delta;u}}_{afi}<0\\ {\mathrm{\&lambda;}}_{afd}=\min (0.9995\underset{i}{\min }(\frac{-z_i}{{\mathrm{\&Delta;z}}_{afi}},\frac{-w_i}{{\mathrm{\&Delta;w}}_{afi}}),1)\\ {\mathrm{\&Delta;z}}_{afi}<0,{\mathrm{\&Delta;w}}_{afi}>0\end{array}\right.---\left(35\right)$

Update Center Parameter σ and the predictive value μ of barrier parameter_af:

$\left\{\begin{array}{c}\mathrm{\&sigma;}=min\left({\left(\frac{{\mathrm{\&rho;}}_{af}}{\mathrm{\&rho;}}\right)}^3,0.1\right)\\ {\mathrm{\&mu;}}_{af}=\mathrm{\&sigma;}\frac{\mathrm{\&rho;}}{2r}\end{array}\right.---\left(36\right)$

In formula (36), duality gap ρ and affine duality gap ρ_afFor:

ρ=l^Tz-u^Tw (37)

ρ_af=(l+ λ_afpΔl_af)^T(z+λ_afdΔz_af)-(u+λ_afpΔu_af)^T(w+λ_afdΔw_af) (38)

(5) each variable and class expansion variable are updated as the following formula:

$\left\{\begin{array}{c}x=x+{\mathrm{\&lambda;}}_p\mathrm{\&Delta;}x\\ l=l+{\mathrm{\&lambda;}}_p\mathrm{\&Delta;}l\\ u=u+{\mathrm{\&lambda;}}_p\mathrm{\&Delta;}u\\ \mathrm{\&alpha;}=\mathrm{\&alpha;}+{\mathrm{\&lambda;}}_p\mathrm{\&Delta;}\mathrm{\&alpha;}\\ \mathrm{\&beta;}=\mathrm{\&beta;}+{\mathrm{\&lambda;}}_p\mathrm{\&Delta;}\mathrm{\&beta;}\\ y=y+{\mathrm{\&lambda;}}_d\mathrm{\&Delta;}y\\ z=z+{\mathrm{\&lambda;}}_d\mathrm{\&Delta;}z\\ w=w+{\mathrm{\&lambda;}}_d\mathrm{\&Delta;}w\end{array}\right.---\left(39\right)$

(6) according to the convergence set, it is judged that whether optimal load flow restrains.If convergence, stop calculating, otherwise turn step Rapid 2.

As shown in Figure 4, the concrete application examples 5 calculating side as a example by IEEE 118 bus test system, to present invention discussion Method is verified and analyzes.Wherein, system carrying out two subregions, node 23,68,69 as separation.Test system basic Data message is as shown in table 1.

Table 1IEEE 118 bus test system master data

With the minimum object function of system losses, for the calculating effect of comparison algorithm, use following three kinds of modes to survey Test system carries out optimal load flow calculating.

(1) the non-subregion of electrical network, uses nonlinear preprocessing to calculate.

(2) electrical network two subregion, uses nonlinear preprocessing to calculate.

(3) electrical network two subregion, uses class expansion variable interior point method to calculate.

Wherein, using border difference θ (x) ＜ 0.03 and iterations>50 as global optimum convergence criterions, with ρ< 1D-5 and | | h (x) | | the convergence criterion that ∞ < 2D-4 optimizes as subsystem.Optimum results is as shown in table 2.

Table 2 optimum results contrasts

As can be seen from Table 2, on the calculating time, it is decreased obviously after have employed partitioning algorithm.In method 1, electrical network During non-subregion, although iterations is less, but iteration is time-consumingly too much every time, therefore total time is still far above partitioning algorithm.In method 3 In, do not save total calculating time only with partitioning algorithm, and due to the introducing of class expansion variable, expand optimization problem Solution space so that the generating expense phase ratio method 2 after optimization has declined, and total iterations is without notable change.

A kind of electric system tide optimization method lax based on sub area division and class expansion variable of the present invention, its useful effect Fruit is:

(1) method that the present invention proposes meets the dispersion of Large-Scale Interconnected electric power system data, the feature of structure distribution；

(2) the method significantly reduces the dimension of subsystem so that the optimization problem scale of each subsystem is less, improves The calculating speed of network optimum trend；

(3) optimizing of each subregion calculates data and the boundary node data having only to one's respective area, it is not necessary to outside region in detail Thin mathematical model, it is to avoid the bottleneck of data transmission, accelerates optimal speed；

(4) the class expansion variable interior point method that subsystem optimization problem uses expands the solution space of optimization problem, improves Convergence.

Finally should be noted that: above example is only in order to illustrate that technical scheme is not intended to limit, to the greatest extent The present invention has been described in detail by pipe with reference to above-described embodiment, and those of ordinary skill in the field are it is understood that still The detailed description of the invention of the present invention can be modified or equivalent, and any without departing from spirit and scope of the invention Amendment or equivalent, it all should be contained in the middle of scope of the presently claimed invention.

Claims (8)

1. the electric system tide optimization method relaxed based on sub area division and class expansion variable, it is characterised in that described Method comprises the steps:

Step 1. decomposes power system, obtains multiple separate subsystem；

Step 2. coordinates each described subsystem, obtains power system decomposition-coordination model based on sub area division；

Step 3. uses Auxiliary Problem Principle that described power system decomposition-coordination model is carried out distributed parallel optimal load flow meter Calculate；

Step 4. uses class expansion variable interior point method to solve each described subsystem respectively, and the trend obtaining described power system is optimum Solve.

2. the method for claim 1, it is characterised in that before described step 1, determines the optimum tide of described power system Stream mathematical model:

$f\left(x\right)=\left\{\begin{array}{cc}\underset{x}{\min }& f\left(x\right)\\ s.t.& h\left(x\right)=0\\ & \underset{\&OverBar;}{g}\&le;g\left(x\right)\&le;\stackrel{\&OverBar;}{g}\end{array}\right.---\left(1\right)$

In formula (1): f (x) is object function, x is the column vector that independent variable is constituted；H (x) and g (x) are respectively equality constraint and not Equality constraint；gWithIt is respectively lower limit and the higher limit of inequality constraints condition g (x).

3. method as claimed in claim 2, it is characterised in that described step 1 includes:

1-1., according to the geographical distribution of described power system, determines boundary node at interconnection；

1-2. judges whether described boundary node is active node；

The most then enter step 1-3；

If it is not, then enter step 1-4；

1-3. increases virtual passive bus on described interconnection, using passive bus with it as new boundary node, enters Step 1-4；

Described power system is decomposed at each described boundary node by 1-4., obtains multiple separate described subsystem.

4. method as claimed in claim 3, it is characterised in that described step 2 includes:

2-1. using described boundary node parameter as coordination variable；

Coordinated by the described coordination variable of exchange between each described subsystem of 2-2., obtain electric power based on sub area division System decomposition Coordination Model.

5. method as claimed in claim 4, it is characterised in that described step 2-2 includes:

A. coordinated by the described coordination variable of exchange between each described subsystem, obtain described electricity based on sub area division Object function minf (x) of Force system decomposition-coordination model:

Minf (x)=f₁(x₁)+f₂(x₂)….+f_n(x_n) (2)

In formula (2), x1, x2...xn are respectively described subsystems；f₁(x₁)、f₂(x₂)...f_n(x_n) it is respectively each described son The object function of system；

B. the boundary constraint of described power system decomposition-coordination model is determined.

6. method as claimed in claim 5, it is characterised in that described step 3 includes:

3-1., according to the kernel function in described Auxiliary Problem Principle, constructs the auxiliary letter of described power system decomposition-coordination model Number；

The primal problem that described power system decomposition-coordination model solves is converted into and solves the saddle point of described auxiliary problem (AP) and ask by 3-2. Topic；

Described auxiliary problem (AP) is iterated by 3-3. according to the threshold value arranged, until loop ends, obtains described power system and divides Solve the distributed parallel optimal load flow result of calculation of Coordination Model.

7. method as claimed in claim 6, it is characterised in that described step 4 includes:

4-1. increases identical class expansion variable in the generic inequality constraints equation in each described subsystem, and described Class expansion variable is respectively provided with by upper and lower bound；

4-2. increases the expansion variable quadratic term of band penalty factor in the object function of each described subsystem, expands inequality constraints Solution space；

Logarithm barrier function inequality constraints as equality constraint, is processed equality constraint with method of Lagrange multipliers by 4-3., Obtain the Lagrangian after introducing class expansion variable；

4-4. derivation KKT condition, obtains the update equation group of the Lagrangian after described introducing class expansion variable；

4-5. solves described update equation, obtains the trend optimal solution of described power system.

8. method as claimed in claim 7, it is characterised in that described step 4-5 includes:

(1) described equation group is initialized: putting iterations is 0, given independent variable, relaxation factor, Lagrange multiplier, class extension The initial value of variable；

(2) calculating target function and the Jacobian matrix of constraint and extra large gloomy matrix, forms linear system；

(3) solve described linear system, obtain the correction of each iteration；

(4) the affine step-length of former variable and the affine step-length of dual variable are calculated；

(5) each variable and class expansion variable are updated；

(6) according to the convergence set, it is judged that whether optimal load flow restrains；

The most then stop calculating, obtain the trend optimal solution of described power system；

If it is not, then return step (2).