CN104795811A - Power system interval state estimation method - Google Patents

Abstract

The invention discloses a power system interval state estimation method. According to the power system interval state estimation method, interval analysis is applied to state estimation with the consideration of parameter uncertainty, the measured quantity and parameters of a power grid can be reasonably utilized to estimate the distribution interval of the state quantity, and therefore more sufficient electrical information is provided for dispatching operating personnel. The power system section state estimation method comprises the steps that firstly, system analysis is conducted on the influence of parameter uncertainty on state analysis, and a state estimation model with the consideration of parameter uncertainty is provided; secondly, a variable substitution technology is provided, nonlinear state estimation is converted into two-step linearization, and conservatism in the interval estimation process is eliminated through an affine technology; finally, an IEEE14 node system is tested. The result shows that the power system interval state estimation method has high estimation accuracy.

Description

A kind of electric power system state of section method of estimation

Technical field

The present invention relates to a kind of electric power system state of section method of estimation, belong to power system monitoring, uncertainty analysis and control technology field.

Background technology

In the research calculated electric system simulation, choosing of method and model is no doubt important, and the accuracy that model parameter obtains is also very large on the impact of simulation result.By installing dynamic monitor in a large number in systems in which, can contribute to obtaining system element parameter relatively accurately.But nonetheless, the uncertainty of simulation parameters or objective reality, it is approximate that reason has be simulation mathematical model itself at two: one to be to entity, and model parameter is only the half quantification index describing entity attribute; Two is, when obtaining the concrete numerical value of model parameter, because the existence of error in measurement and the finite accuracy of computer represent, can not obtain the exact value of parameter.

Uncertain problem in state estimation is one of difficult point of state estimation always.At present, consider that probabilistic algorithm mainly contains 3 kinds: 1, Fuzzy Analysis Method in field of power, the method is a kind of mathematical tool for research boundary fuzzy problem; 2, probability analysis method, namely utilizes the random information of the mode treatment variable of probability; 3, Interval Analytical Method, utilizes that intervl mathematics and Novel Interval Methods process place are prolonged clearly, the indefinite information of intension, as long as namely know its approximate range and do not know the information of its determined value.Because interval analysis does not need artificial hypothesis, avoid the impact of subjective factor on result of calculation of people to a certain extent, thus be widely used on process uncertain problem.

Interval algorithm is due to relativity problem and the impact of wrapping up effect, and estimated result is too conservative sometimes, thus loses practical value.Affine arithmetic, by recording the dependence between each Uncertainty, makes full use of these extraneous informations, can obtain the result that obtain a bit more accurate than conventional Interval Computation.Document proposes to utilize affine technology to overcome correlation in interval arithmetic, simultaneously in order to improve computational efficiency, adopting optimum theory to solve Interval Power Flow problem, greatly reducing algorithm complex.

In mathematical form, state estimation is a Nonlinear Optimization Problem containing constraints, and the optimization problem of interval nonlinear equation group is still one of difficult point.Can consider using uncertain variables as inequality constraints process, and meet for target function with maximum constrained, the method is only a kind of approximate model, and error is larger.Another kind of feasible method is that to utilize phasor measurement unit (pharos measurement unit, PMU) to measure function be that the feature of linear model carries out interval estimation, and obtain the bound of quantity of state, computational speed is fast.But, because PMU measure configuration is still limited, also cannot realize the whole network and cover installation, simultaneously in order to make full use of existing traditional data collection and supervisory control system (supervisory controland data acquisition, SCADA) measurement information, reduce the information loss being measured by traditional SCADA and be converted into PMU measurement participant status and estimate to cause, the fast algorithm of a large amount of scholar to state estimation is studied, wherein, the bilinearity technology that Antonio proposes, namely by arranging intermediate variable, nonlinear state Eq is made to be converted into the inearized model of 2 steps, its estimated result and WLS estimate basically identical, greatly can reduce amount of calculation.

Interval analysis and affine arithmetic

After Moore in 1966 proposes interval analysis, cause the concern of numerous scholars, its objective is the bound obtaining the uncertain result of calculation caused.In order to distinguish point value number, the present invention is for given interval number meet for all set, that is:

$\tilde{x}=[\underset{\‾}{x},\stackrel{\‾}{x}]=\{x\∈R|\underset{\‾}{x}\≤x\≤\stackrel{\‾}{x}\}$

In formula, be respectively lower bound and the upper bound of interval number, and meet r is set of real numbers.Distinguishingly, when time, be point value real number.

The operation method of interval number is as follows:

Addition: $\tilde{x}+\tilde{y}=[\underset{\‾}{x}+\underset{\‾}{y},\stackrel{\‾}{x}+\stackrel{\‾}{y}]$

Subtraction: $\tilde{x}-\tilde{y}=[\underset{\‾}{x}-\stackrel{\‾}{y},\stackrel{\‾}{x}-\underset{\‾}{y}]$

Multiplication: $\tilde{x}\×\tilde{y}=[\min (\underset{\‾}{x}\underset{\‾}{y},\underset{\‾}{x}\stackrel{\‾}{y},\stackrel{\‾}{x}\underset{\‾}{y},\stackrel{\‾}{x}\stackrel{\‾}{y}),\max (\underset{\‾}{x}\underset{\‾}{y},\underset{\‾}{x}\stackrel{\‾}{y},\stackrel{\‾}{x}\underset{\‾}{y},\stackrel{\‾}{x}\stackrel{\‾}{y})]$

Division: $\tilde{x}\÷\tilde{y}=\tilde{x}\×[\min (\frac{1}{\underset{\‾}{y}},\frac{1}{\stackrel{\‾}{y}}),\max (\frac{1}{\underset{\‾}{y}},\frac{1}{\stackrel{\‾}{y}})]$ (wherein, $0\∉\tilde{y}$ )

Interval friendship: $\tilde{x}\∩\tilde{y}=[\max \left(\underset{\‾}{x,}\underset{\‾}{y}\right),\min (\stackrel{\‾}{x},\stackrel{\‾}{y})]$

Interval also: $\tilde{x}\∪\tilde{y}=[\min \left(\underset{\‾}{x,}\underset{\‾}{y}\right),\max (\stackrel{\‾}{x},\stackrel{\‾}{y})]$

Interval comparison: $\tilde{x}<\tilde{y}\&DoubleLeftRightArrow;\stackrel{\&OverBar;}{x}<\underset{\&OverBar;}{y}$

Interval comprises: $\tilde{x}\&SubsetEqual;\tilde{y}\&DoubleLeftRightArrow;\underset{\&OverBar;}{x}\&GreaterEqual;\underset{\&OverBar;}{y}$ And $\stackrel{\&OverBar;}{x}\&le;\stackrel{\&OverBar;}{y}$

In order to overcome deficiency too conservative in interval analysis, Comba and Stolfi proposes affine arithmetic.In affine arithmetic, an interval variable affine form represents, processes interval variable by affine technology:

$\tilde{x}=x_0+x_1{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_1+\&CenterDot;\&CenterDot;\&CenterDot;x_n{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_n=x_0+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_i$

In formula, represent noise unit, x _idetermine the size of noise unit.If same noise unit appear in different interval variables, then mean that their uncertainty has certain contact and interdependency.

Affine number and interval number can be changed mutually, assuming that some interval numbers are get simultaneously $x_o=(\underset{\&OverBar;}{x}+\stackrel{\&OverBar;}{x})/2=\mathrm{mid}\left(\tilde{x}\right)$ For mid point, $x_1=(\stackrel{\&OverBar;}{x}-\underset{\&OverBar;}{x})/2=\mathrm{rad}\left(\tilde{x}\right)$ For interval width, then affine form be:

$\tilde{x}=x_0+x_1{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_1$

In formula, ${\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_1\&Element;[-1,1].$

Otherwise, for affine number definition then the interval number form of its correspondence is:

$\tilde{x}=[x_0-r_{\mathrm{ad}},x_0+r_{\mathrm{ad}}]$

Given affine number then affine number with algorithm as follows.

Addition and subtraction: $\tilde{x}\&PlusMinus;\tilde{y}=(x_0+y_0)+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}(x_i\&PlusMinus;y_i)\stackrel{\&OverBar;}{{\mathrm{\&epsiv;}}_i}$

Multiplication: $\tilde{x}\&times;\tilde{y}=(x_0+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_i)\&times;(y_0+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{y}_i{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_i)=$

$x_0{y}_0+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}(x_0{y}_i+y_0{x}_i){\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_i+\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_i\right)\&times;\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{y}_i{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_i\right)$

From multiplication formula, 2 affine numbers are multiplied the quadratic polynomial that can produce about noise source, and its simplified calculation method is as follows:

$\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_i\right)\&times;\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{y}_i{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_i\right)=\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\right|x_i\left|\right)\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\right|y_i\left|\right)$

Similar with interval number, the division of affine number only need be converted into the multiplication of its inverse, that is:

Summary of the invention:

Goal of the invention: the present invention proposes a kind of electric power system state of section method of estimation based on affine arithmetic and bilinearity technology, has possessed better estimated accuracy.

Technical scheme: the present invention proposes a kind of electric power system state of section method of estimation, comprises the following steps,

The parameter of input electric power system, topology, measurement information, and the uncertain region determining branch parameters, measurement amount;

Set up one-phase linear state estimation model, and solve one-phase state of section estimation problem with affine arithmetic, obtain the interval value of one-phase quantity of state;

Intermediate variable nonlinear transformation, obtains the interval value of intermediate variable quantity of state;

Set up two-stage linear state estimation model, and solve two-stage state of section estimation problem with affine arithmetic, obtain the interval value of two-stage quantity of state;

The interval estimation value of each node voltage amplitude, phase angle in output network.

Preferably, described one-phase linear state estimation comprises:

For every bar branch road of connection bus i and bus j, be defined as follows variable:

K _ij＝V _iV _jcosθ _ij

L _ij＝V _iV _jsinθ _ij

In formula: V _i, V _jbe respectively the voltage magnitude of bus i, j, θ _ij=θ _i-θ _j, θ _i, θ _jbe respectively the voltage phase angle of bus i, j.

For the every bar bus in system, definition voltage magnitude square is new variable:

U _i＝V _i ²

Assuming that system comprises N bar bus, b bar branch road, then one-phase linear state estimation is introduced N+2b and is tieed up quantity of state y:

y＝{U _i,K _ij,L _ij}

Thus m dimension measurement vectorial z and quantity of state y can be expressed as linear relationship:

z＝By+e _z

Preferably, described intermediate variable nonlinear transformation comprises:

The nonlinear transformation of intermediate variable, for waiting dimension conversion, is defined as follows N and ties up variable α _i, b ties up variable α _ij, θ _ij:

α _i＝lnU _i＝2lnV _i

α _ij＝ln(K _ij+L _ij)＝α _i+α _j

${\mathrm{\&theta;}}_{\mathrm{ij}}=\arctan \left(\frac{L_{\mathrm{ij}}}{K_{\mathrm{ij}}}\right)={\mathrm{\&theta;}}_i-{\mathrm{\&theta;}}_j$

Make u={ α _i, α _ij, θ _ij, then N+2b ties up variable u and y in non-linear relation

Preferably, described two-stage linear state estimation comprises:

Definition 2N-1 ties up quantity of state x=[α θ] ^t(being fixed as 0 with reference to the phase angle of bus), then two-stage quantity of state x and intermediate variable u is following linear relationship:

u＝Cx+e _u

$C=\left[\begin{array}{cc}I& 0\\ \left|A^T\right|& 0\\ 0& A_r^T\end{array}\right]$

In formula: I is unit battle array, A is node incidence matrices, and Ar does not comprise the node incidence matrices with reference to bus.

Beneficial effect: interval analysis affine arithmetic and bilinearity technology are applied in the state estimation considering parameter uncertainty by the present invention, first, on the impact of state estimation, network analysis is carried out to parameter uncertainty, proposes the state estimation model considering parameter uncertainty.Propose substitution of variable technology and nonlinear state Eq is converted into two step linearisations, and utilize affine technology to overcome conservative in interval estimation.The present invention comparatively reasonably utilizes the distributed area of power grid measurement amount and the parameter Estimation amount of doing well, thus provides management and running personnel more sufficient electric information, and has good estimated accuracy.

Accompanying drawing illustrates:

Fig. 1: the flow chart of a kind of electric power system state of section of the present invention method of estimation;

Fig. 2: the disaggregation of Interval linear equation group and shell.

Embodiment:

Be described in detail below in conjunction with the techniqueflow of accompanying drawing to invention:

Mathematically, Power system state estimation is the nonlinear optimal problem containing constraints.Consider parameter uncertainty, because state of section is estimated to be difficult to direct solution, therefore, the present invention adopts the method for two step linearisation state estimation, that is: the thought of substitution of variable is utilized, one group of intermediate variable is set, thus electric power system nonlinear state Eq is converted into the linear state estimation problem of two steps, and comprise the nonlinear transformation of a step variable between two linear state estimations.

One-phase linear state estimation

For every bar branch road of connection bus i and bus j, be defined as follows variable:

K _ij＝V _iV _jcosθ _ij

L _ij＝V _iV _jsinθ _ij

In formula: V _i, V _jbe respectively the voltage magnitude of bus i, j, θ _ij=θ _i-θ _j, θ _i, θ _jbe respectively the voltage phase angle of bus i, j.

For the every bar bus in system, definition voltage magnitude square is new variable:

U _i＝V _i ²

Assuming that system comprises N bar bus, b bar branch road, then one-phase linear state estimation is introduced N+2b and is tieed up quantity of state y:

y＝{U _i,K _ij,L _ij}

The measurement amount that then SCADA system provides becomes following linear relationship with quantity of state y:

Branch power measures:

$P_{\mathrm{ij}}^m=(g_{\mathrm{si}}+g_{\mathrm{ij}})U_i-g_{\mathrm{ij}}{L}_{\mathrm{ij}}-b_{\mathrm{ij}}{K}_{\mathrm{ij}}+e_{P_{\mathrm{ij}}}$

$Q_{\mathrm{ij}}^m=-(b_{\mathrm{si}}+b_{\mathrm{ij}})U_i+b_{\mathrm{ij}}{L}_{\mathrm{ij}}-g_{\mathrm{ij}}{K}_{\mathrm{ij}}+e_{Q_{\mathrm{ij}}}$

In formula: g _ij, b _ijbe respectively conductance, the susceptance of branch road π type equivalent electric circuit, g _si, b _sibe respectively bus i side conductance, susceptance over the ground.

Node injects and measures:

$P_i^m=\underset{j\&Element;i}{\mathrm{\&Sigma;}}{P}_{\mathrm{ij}}+e_{P_i}$

$Q_i^m=\underset{j\&Element;i}{\mathrm{\&Sigma;}}{Q}_{\mathrm{ij}}+e_{Q_i}$

Voltage magnitude measures:

${\left({V_i}^2\right)}^m=U_i+e_{U_i}$

In formula: e is error in measurement vector, and supposition e Normal Distribution.

Thus m dimension measurement vectorial z and quantity of state y can be expressed as linear relationship:

z＝By+e _z

When considering parameter uncertainty, after utilizing affine technology to calculate above-mentioned one-phase linear state estimation function, the interval estimation value of one-phase quantity of state can be obtained.

Intermediate variable nonlinear transformation

The nonlinear transformation of intermediate variable, for waiting dimension conversion, is defined as follows N and ties up variable α _i, b ties up variable α _ij, θ _ij:

α _i＝lnU _i＝2lnV _i

α _ij＝ln(K _ij+L _ij)＝α _i+α _j

${\mathrm{\&theta;}}_{\mathrm{ij}}=\arctan \left(\frac{L_{\mathrm{ij}}}{K_{\mathrm{ij}}}\right)={\mathrm{\&theta;}}_i-{\mathrm{\&theta;}}_j$

Make u={ α _i, α _ij, θ _ij, then N+2b ties up variable u and y in non-linear relation

Two-stage linear state estimation

Definition 2N-1 ties up quantity of state x=[α θ] ^t(being fixed as 0 with reference to the phase angle of bus), then two-stage quantity of state x and intermediate variable u is following linear relationship:

u＝Cx+e _u

$C=\left[\begin{array}{cc}I& 0\\ \left|A^T\right|& 0\\ 0& A_r^T\end{array}\right]$

In formula: I is unit battle array, A is node incidence matrices, and Ar does not comprise the node incidence matrices with reference to bus.When considering parameter uncertainty, after utilizing affine technology to calculate two-stage linear state estimation function, the interval estimation value of two-stage state variable can be obtained.

When measuring Noise, the object of state estimation is filtering, to estimate system running state more accurately.Bilinearity State Estimation Theory linear filtering in two steps, in measurement noise Normal Distribution, under the condition not containing measurement rough error, the precision of bilinearity state estimation is suitable with traditional WLS, and advantage is to substantially increase computational efficiency.

Affine technology is utilized to estimate POWER SYSTEM STATE when considering parameter uncertainty

Under the condition of given network configuration, metric data, the measurement equation of electric power system is:

z _i＝h _i(x)+r _ii＝1,2,...,m

In formula, x is system state variables, comprises node voltage amplitude and phase angle; h _ix () is No. i-th measurement function measured; z _ibe No. i-th measuring value measured; r _ibe No. i-th measurement residuals.

Usually, Interval linear equation group can form as follows:

[A][x]＝[b]

In formula, [A] ∈ I ^{n × n}(R) be the interval real coefficient matrix on n × n rank, [x] ∈ I ^{n × 1}(R) be the interval solutions vector on rank, n × 1, [b] ∈ I ^{n × 1}(R) be the interval constant vector on rank, n × 1.

Solve for above-mentioned Interval linear equation group, it is interval that its mathematical meaning is to find all solutions met the following conditions:

Σ([A],[b])＝{x|Ax＝b,A∈[A],b∈[b]}

The feature of Interval linear equation group solution is described with a simple mathematical example.

Solve Interval linear equation group, the distribution of its disaggregation S and shell (hull) are as shown in annex Fig. 2.

As far as possible the target that solves of Interval linear equation group obtains the smallest interval vector comprising disaggregation S, also known as the shell of Interval linear equation group for this reason, and the region that namely in annex Fig. 2, rectangle surrounds:

$\left[x\right]=\left(\begin{array}{c}[-\mathrm{6,6}]\\ [-\mathrm{4,4}]\end{array}\right)$

Theoretically, overlap with S if shell can be tried to achieve, be then optimal solution, but due to the conservative in interval computation process, be often greater than S.In addition, owing to having certain " correlation " in parameter uncertainty in state of section estimation, and these " correlations " are considered as same uncertain variables in pure interval computation process, thus cause interval extension.Such as, for interval function:

f([x])＝[x]-2[x]x∈[1,2]

The computational process of pure interval computation is [1,2]-[2,4]=[-3,0], but, in fact f ([x])=-[x]=[-2,-1], this is because in pure interval computation process, do not consider that two uncertain variables are actually same variable, thus cause interval extension.After utilizing affine technology to calculate this function, can correctly be separated.

Sample calculation analysis

In order to verify the validity of the inventive method, carry out programming realization based on Matlab platform, and launch test on the PC compatible of 3.2GHz dominant frequency, 4G internal memory.First, with IEEE 14 node for research object is analyzed.The uncertainty interval of each metric data and system parameters is arranged by table 1.Compare with interior some optimal load flow algorithm (interior pointmethod optimal power flow, IPMOPF), Monte Carlo method (emulating for about 20000 times), checking the present invention put forward the accuracy of algorithm.

Table 2 and table 3 give when the interval of parameter uncertainty is 5% standard deviation, and the state of section estimated result of distinct methods to IEEE 14 node system compares, and wherein, node voltage amplitude is perunit value.For No. 14 buses, analyze the interval computation performance of distinct methods: Monte Carlo method is by carrying out a large amount of simulation calculation to uncertain variables, the interval that can obtain close to actual value distributes, its result is [1.033 67,1.037 35] ∠ [-0.289 16,-0.272 00], node voltage amplitude interval width is 0.003 68; When adopting IPMOPF to calculate, result is [1.030 33,1.04168] ∠ [-0.292 81 ,-0.266 00], and voltage magnitude width is 0.011 35; Utilize institute of the present invention extracting method, result is [1.031 78,1.04 0 22] ∠ [-0.289 32 ,-0.269 43], and the interval width of voltage magnitude is 0.00844, and comparing IPMOPM can reduce 25.64%.Thus compared to IPMOPM, the inventive method interval estimation precision is higher.

In addition, time of Monte Carlo simulation is 525.2s, IPMOPM computing time is 0.25s, and the inventive method computing time is 0.11s, and thus the computational efficiency of the inventive method is higher.

Table 1 measures the uncertain standard deviation with electrical network parameter

The voltage magnitude interval estimation results contrast (5% standard deviation) of table 2 IEEE 14 node system

The voltage phase angle interval estimation results contrast (5% standard deviation) of table 3 IEEE 14 node system

Claims (4)

1. an electric power system state of section method of estimation, is characterized in that, comprises the following steps,

The parameter of input electric power system, topology, measurement information, and the uncertain region determining branch parameters, measurement amount;

Set up one-phase linear state estimation model, and solve one-phase state of section estimation problem with affine arithmetic, obtain the interval value of one-phase quantity of state;

Intermediate variable nonlinear transformation, obtains the interval value of intermediate variable quantity of state;

Set up two-stage linear state estimation model, and solve two-stage state of section estimation problem with affine arithmetic, obtain the interval value of two-stage quantity of state;

The interval estimation value of each node voltage amplitude, phase angle in output network.

2. electric power system state of section method of estimation according to claim 1, is characterized in that, described one-phase linear state estimation comprises:

For every bar branch road of connection bus i and bus j, be defined as follows variable:

K _ij＝V _iV _jcosθ _ij

L _ij＝V _iV _jsinθ _ij

In formula: V _i, V _jbe respectively the voltage magnitude of bus i, j, θ _ij=θ _i-θ _j, θ _i, θ _jbe respectively the voltage phase angle of bus i, j.

For the every bar bus in system, definition voltage magnitude square is new variable:

U _i＝V _i ²

Assuming that system comprises N bar bus, b bar branch road, then one-phase linear state estimation is introduced N+2b and is tieed up quantity of state y:

y＝{U _i,K _ij,L _ij}

Thus m dimension measurement vectorial z and quantity of state y can be expressed as linear relationship:

z＝By+e _z。

3. electric power system state of section method of estimation according to claim 1, is characterized in that, described intermediate variable nonlinear transformation comprises:

The nonlinear transformation of intermediate variable, for waiting dimension conversion, is defined as follows N and ties up variable α _i, b ties up variable α _ij, θ _ij:

α _i＝lnU _i＝2lnV _i

α _ij＝ln(K _ij+L _ij)＝α _i+α _j

Make u={ α _i, α _ij, θ _ij, then N+2b ties up variable u and y in non-linear relation

4. electric power system state of section method of estimation according to claim 1, is characterized in that, described two-stage linear state estimation comprises:

Definition 2N-1 ties up quantity of state x=[α θ] ^t(being fixed as 0 with reference to the phase angle of bus), then two-stage quantity of state x and intermediate variable u is following linear relationship:

u＝Cx+e _u

In formula: I is unit battle array, A is node incidence matrices, and Ar does not comprise the node incidence matrices with reference to bus.