Fix problem in state estimation

[rdzman]

It seems that there is a bug in computing the S matrix that affects the compulation of the derivative matrix H. As, such a new function has been developed to compute S. Moreover, the code is hardwired to solve the problem 6.7 in the book "Computational Methods for Electric Power Systems" by Mariesa Crow, 1st edition. In particular, it only accepts the measurements in that problem. So, the code is expanded to accommodate any sets of measurements. A test file "test_se_Ex_6_17.m" is provided that solves example 6.17 in the above 2nd edit of the book. The file shows the difference in computing H between the original, which is edited to allow more measurements, and the modified code. The modified code produced results that coincide with the solution of example 6.17 in the book.

[rdzman]

Moving MATPOWER/matpower-extras#6 here, since this is the master repo for the state estimation code.

It's just pulled in as a subrepo over there.

[rdzman]

Thanks again, @aldalahmeh, for this contribution. I finally have a chance to begin reviewing it for inclusion in the next release. But rather than merging this pull request directly, which simply adds some new files, I would like to merge your fixes into the existing files. However, before I do that, I'd like to get a better understanding of exactly what issues are being addressed. Keep in mind that the state estimation code was contributed, so I'm not as familiar with it as I am with the MATPOWER code I wrote myself.

If I understand correctly, from your comments and going through the code changes, it seems there are five primary changes in doSEmod() relative to doSE().

1. Addition of real and reactive load as measurements.
   • Just wondering, in the previous version, was the load assumed to be a known quantity, or what?
2. Removal of voltage magnitude from state to be estimated at PV buses.
3. The measurement variances can be specified individually for each measurement. They need not be identical for each type of measurement.
4. The Rinv matrix is now the inverse of the former Rinv matrix.
   • Was this a bug? Or did you just invert the meaning of the sigma inputs?
5. A new method of computing the estimates for the power flows on branches.
   • Can you explain, how the new version is conceptually different from the old? The old simply computes the branch injections from the vector of bus voltages. Can you help me understand what the new version is computing, and why it is correct and the old incorrect?

Finally, what is the effect of these changes/fixes on the other examples, e.g. test_se_14bus.m, test_se_14bus_err.m? And do you have other working examples that can be included? E.g. I believe you mentioned a wiki with some examples in a previous PR, but I couldn't find it.

Once again, @aldalahmeh, thanks for the contribution.

[aldalahmeh]

Dear @rdzman ,

I hope that my contribution is useful for the community. I reply to your points below in order:

1. In state estimation, the load is assumed to be measured with some error. So, the true load values are taken from the power flow solution then noise is added. This is captured in the z structure in line 63 in doSEmod.m file.
2. As mentioned earlier, voltages magnitudes are known at PV buses.
3. True.
4. I just looked into it, and turns out that it was a bug, since in least squares estimation (LSE) provided a diagonal covariance matrix R = diag(1/sigma^_1, 1/sigma^_2, ... , 1/sigma^_N ) then the matrix used in the LSE is R^{-1} = diag(sigma^_1, sigma^_2, ... , sigma^_N ). So, a simpler fix (and less confusing one) is to change line 125 to R_inv = diag(sigma_square);. Would you please do that?
5. I could not pinpoint the error in computing the S matrix and hence the power flow. So I had to write the code from scratch in a way I understood it then, I just coded equations (3.90) and (3.91) in the book mentioned earlier, which are the real and reactive power flows respectively. I did it in a complex and vectorized way in line 8 in the cmptSmat function.

The final outputs of test_se_14bus.m and test_se_14bus_err.m are changed. I have a working example on a IEEE 14 bus test case, which can be added. On a different note, I found that there a certain pieces of information missing that may make things easier to understand regarding the naming convention of from-to & to-from bus injections. I suggest to add a file explaining such concepts. What do you think?

[rdzman]

Thanks, @aldalahmeh, for your response. Here is some more discussion on each of the points, in order.

1. Since we can only estimate net power injections at a bus, not separate load and generation quantities, it seems we can really only make use of a load measurement if we assume any generation at the same bus is known (or non-existent, i.e. equal to zero). And conversely for generation measurements, we can only use them if we assume the load at the bus is known. It appears that your version still includes, as in the original code, the assumption that when estimating generation, the load at the bus is a known quantity, pulled from the PD and QD columns of bus. It also assumes there is no generation at buses that have load measurements. If we are going to include measurements for both load and generation, we need to think through how to handle cases that have both load and generation at a bus, or multiple generators, for that matter. My initial thought is that the correct thing to do is to discard any measurements at such a bus unless we have measurements for all injections (loads and gens) at the bus, in which case they should be combined and treated internally as a single measurement of the net bus injection.

2. It seems to me that completely eliminating voltage magnitudes at PV buses may impose an unnecessary restriction on the application. Is it true in real-world state estimation usage that voltage magnitudes are assumed to be known? Could we not arrive at the same solution by keeping them as state variables and simply specifying the PV bus voltage setpoints as measurements with a very small sigma?

3. 👍🏼

4. I looked into this a bit more and, from the references I can find, if sigma is the standard deviation of the errors associated with the measurments, then sigma.^2 is the variance and the corresponding variance-covariance matrix, by definition, would be R = diag(sigma.^2). In the weighted LSE we use the inverse as the matrix of weights, so R_inv = diag(1./sigma.^2) as in the original code. It seems to me that if you are saying the variance-covariance matrix has 1/sigma on the diagonal, then you are redefining sigma to be something other than the standard deviation (e.g. the inverse of the SD). Conclusion: I believe the original code was correct, but feel free to convince me if you still believe otherwise.

5. I'm still trying to understand this one. What exactly does the S matrix (Sij in cmptSmat()) represent? That is, what is entry (i, j)? Maybe you could e-mail me a scan/photo of the relevant pages from the book you reference. (I don't have a copy.) Conceptually, the original code is very straightforward. It is computing the branch power flows from the voltages using code that has been thoroughly validated, and the Hessian terms are consistent with the computed flows. If your code is giving a different value for the flows, then it is computing something slightly different (which presumably needs a different Hessian). The only thing I can think of, given that it is computed entirely from the voltages and only Ybus (as opposed to Yf and Yt) is that it must include bus shunt elements somehow in the computation of branch flows. I don't understand how that could be correct.

Thanks again, and hopefully we can sort these things out and finalize any updates to the SE code.

[rdzman]

On a different note, I found that there a certain pieces of information missing that may make things easier to understand regarding the naming convention of from-to & to-from bus injections. I suggest to add a file explaining such concepts. What do you think?

Certainly. Please share what you think is missing that could be helpful.

[aldalahmeh]

@rdzman In response to the interesting points you have raised:

1. Your observation is true. The original code either assumes a generation measurement or load measurement (both of which are noisy). However, if you look at it from a measurement point of view, the power flow at the specific bus is measured. So having multiple generators and/or loads can be accommodated if the net power is consider. I think this can be done, do suggest adding it to this version?

2. From what I read in literature, the PV bus has a known voltage. However, in theory you can assume the voltage to be uncertain to some degree, as you suggested. However, we need a separate equations for the voltage magnitude and phase (can we have that?). Otherwise, the system will be rank deficient and hence unobservable leading to the failure of the Newton-Raphson algorithm.

3. Great.

4. I believe the definition of the R matrix is different as per the mention book. See the attached snapshot from it below.
   

5. I found the bug in the original code by studying the mentioned example in the book (please ignore my earlier reference to eqs (3.90) and (3.91) it was not in place). See the snapshot of the example below
   

Consider the nonlinear functions h_2(x), h_3(x), h_4(x) and h_5(x), where the measurements (z) are power measurements being power flows (with some abuse of defition) e.g. P_13 & Q_21 or power demands (P_3 & Q_2). Note in the case of h_2(x) and h_3(x), there exit terms involving the square bus voltages, which cannot be computed using original code
Sfe = V(f) .* conj(Yf * V); Ste = V(t) .* conj(Yt * V);

Conceptually, h_2(x), h_3(x), h_4(x) and h_5(x) are power flows. Let us take P_13 for example. It is in fact the power flow from bus 1 to bus 3 (see snapshot below) minus the power demand at the bus it self (hence the squared bus voltage term).

So I computed the complex power flow matrix Sij in the function cmpt_Smat. I coded this matrix in the one-liner below long time ago and it took me a while to decrypt it!
Sij = V * V'.*conj(Ybus) - (V .* conj(V) * ones(1,nb)).*conj(Ybus);

What it does is computes the complex power flow as described above, where the first term is the power flow between the buses and the second term is the bus power (demand).

[rdzman]

1. The complexity here is with not with the math, but with handling each of the special cases appropriately and making sure it's clear to the user what that is. The cases as I see it are ...
   (a) load measurement only, no generator at the same bus - straightforward, net injection corresponds only to the load
   (b) generator measurement only, single generator with no load at same bus - straightforward, net injection corresponds only to the generator
   (c) load measurement, with gen at same bus with no measurement - do we discard this, or make an assumption about the generation to get the relationship between estimated net injection and estimated load?
   (d) gen measurement, with another gen or load at same bus with no measurement - similar to (c)
   (e) multiple injections (load, gens) at a bus, each with a measurement - estimates can be created for each measured quantity based on single estimated net injection, using standard deviations of individual measurements to distribute the error across the multiple injections
   So, I think we have two options. The first is to simply document the current assumptions and go with your modified code. The second is to try to do something more reasonable with each of the special cases above, and document that. I'm open to either, but don't have the time to do it myself.
2. I think the original code handles this just fine without rank deficiency. In theory, I think that fixing the voltage magnitude at a PV bus (like your code does) should be equivalent to using the original code and providing a measurement for it with a very small standard deviation. My only concern would be whether that would cause numerical issues. I would love input here from someone with more real-world state estimation experience.
4. I could be wrong, but given my understanding of LSE from other sources, I believe my previous post is correct and the book has a typo. I think it should say that the "inverse of the covariance matrix of the measurements" is ... and (6.157) should say "R^-1 = ...".
5. I confess I didn't actually follow all the math and explanations here, but I think the "bug" is not in the code but in the data. The diagram appears to represent the lines as simple series impedance values (no line charging capacitance considered). Apparently, the modified way of computing these flows ignores the line charging capacitance. You can confirm this by setting mpc.branch(:, BR_B) = 0 and running the original code. This gives the same result as your modified code. Note also that it does it in 3 iterations instead of 5, which I attribute to the fact that the Hessian is correct for the function being computed. In any case, I still believe the original code is correct here.

[aldalahmeh]

1. The complexity here is with not with the math, but with handling each of the special cases appropriately and making sure it's clear to the user what that is. The cases as I see it are ...
   (a) load measurement only, no generator at the same bus - straightforward, net injection corresponds only to the load
   (b) generator measurement only, single generator with no load at same bus - straightforward, net injection corresponds only to the generator
   (c) load measurement, with gen at same bus with no measurement - do we discard this, or make an assumption about the generation to get the relationship between estimated net injection and estimated load?

Case (c) is implemented, as you noted earlier, in doSEmod in line 147, where the generated power and demand power are accounted for.

(d) gen measurement, with another gen or load at same bus with no measurement - similar to (c)
(e) multiple injections (load, gens) at a bus, each with a measurement - estimates can be created for each measured quantity based on single estimated net injection, using standard deviations of individual measurements to distribute the error across the multiple injections

I believe the issue here is the way the mpc is coded, which actually encodes the power network including the generators and loads, it does not allow having multiple generators nor loads.

So, I think we have two options. The first is to simply document the current assumptions and go with your modified code. The second is to try to do something more reasonable with each of the special cases above, and document that. I'm open to either, but don't have the time to do it myself.

I believe is difficult to change the way mpc is coded, so I suggest to code a layer over the mpc to accommodate (d) and (e) cases above, which adds the generators (loads) into a single one and then encoded in the mpc structure.

2. I think the original code handles this just fine without rank deficiency. In theory, I think that fixing the voltage magnitude at a PV bus (like your code does) should be equivalent to using the original code and providing a measurement for it with a very small standard deviation. My only concern would be whether that would cause numerical issues. I would love input here from someone with more real-world state estimation experience.

I think you are referring to adding the voltage and magnitude of the same bus as a measurement, in which the rank will be affected. In contrast to adding them as states. However, as you states, input from a practitioner is always informative.

4. I could be wrong, but given my understanding of LSE from other sources, I believe my previous post is correct and the book has a typo. I think it should say that the "inverse of the covariance matrix of the measurements" is ... and (6.157) should say "R^-1 = ...".

I agree with you that the book's notation is rather confusing. In LSE literature a weighing matrix is used and usually denoted as W. In our case, let W = diag(1./sigma.^2), which is actually the inverse of the covariance matrix. And as you correctly noted, the book denoted R= diag(1./sigma.^2), which usually denotes the covariance matrix causing the confusion here.

So I suggest changing the matrix R used in the code to W.

5. I confess I didn't actually follow all the math and explanations here, but I think the "bug" is not in the code but in the data. The diagram appears to represent the lines as simple series impedance values (no line charging capacitance considered). Apparently, the modified way of computing these flows ignores the line charging capacitance. You can confirm this by setting mpc.branch(:, BR_B) = 0 and running the original code. This gives the same result as your modified code. Note also that it does it in 3 iterations instead of 5, which I attribute to the fact that the Hessian is correct for the function being computed. In any case, I still believe the original code is correct here.

The example included is actually an extension to a previous example, which is example 3.9 with screenshot included below. I encoded the system in the mpc structure below. I example discusses power flow using Newton-Raphson method. I implemented the example using the standard runpf function and everything checks up. The lines originally have charging capacitance and there is no reason not to accommodate it, as the original code did. Unless it is stated that the original code only addresses such special case.

https://www.dropbox.com/scl/fi/8cxr2a3o1rpzhk18y0tq9/Power_flow_case3bus_Ex3_9.m?rlkey=8vyfhyk9stjhszkzvzbwvu5se&dl=0

[rdzman]

Case (c) is implemented, as you noted earlier, in doSEmod in line 147, where the generated power and demand power are accounted for.

Unless I'm mistaken, line 147 in doSEmod relates to generation measurements, where the load is assumed to be known. So it's not (c), but rather a specific case of (d) with a specific way to handle it (assume load is known).

I believe the issue here is the way the mpc is coded, which actually encodes the power network including the generators and loads, it does not allow having multiple generators nor loads.

I'm afraid I don't follow. The mpc does allow for a single load and multiple generators at each bus.

So, I think we have two options. The first is to simply document the current assumptions and go with your modified code. The second is to try to do something more reasonable with each of the special cases above, and document that. I'm open to either, but don't have the time to do it myself.

I believe is difficult to change the way mpc is coded, so I suggest to code a layer over the mpc to accommodate (d) and (e) cases above, which adds the generators (loads) into a single one and then encoded in the mpc structure.

I'm sorry, but again I don't follow what you mean about "the way mpc is coded" or your proposed solution.

4. I could be wrong, but given my understanding of LSE from other sources, I believe my previous post is correct and the book has a typo. I think it should say that the "inverse of the covariance matrix of the measurements" is ... and (6.157) should say "R^-1 = ...".

I agree with you that the book's notation is rather confusing. In LSE literature a weighing matrix is used and usually denoted as W. In our case, let W = diag(1./sigma.^2), which is actually the inverse of the covariance matrix. And as you correctly noted, the book denoted R= diag(1./sigma.^2), which usually denotes the covariance matrix causing the confusion here.

So I suggest changing the matrix R used in the code to W.

Let's just stick with the original code, where Rinv, the inverse of the variance/covariance matrix, was computed and used correctly as the LSE weight matrix.

The lines originally have charging capacitance and there is no reason not to accommodate it, as the original code did. Unless it is stated that the original code only addresses such special case.

Just to confirm, you agree that the original code is correct in this regard?

[rdzman]

To summarize where I currently stand on the 5 primary changes proposed by this PR.

1. Addition of load measurements – This is not consistent with the current assumption that loads are known quantities (used in computing generation measurement errors). More design work is required to allow for consistent handling of the different cases of multiple injections (generators, load) at a bus.
2. Removal of voltage magnitude from PV bus state – I'm not convinced this is necessary / desirable. Would like input from someone with more state estimation experience.
3. Add option to specify individual measurement variances – Seems like a nice feature addition.
4. Invert weight matrix – Original code was correct. Let's not change anything.
5. Different method for computing branch flow estimates – Original code was correct. Let's not change anything.

Since the proposed changes are pretty much independent of one another, my recommendation moving forward is to separate each into its own pull request consisting of changes to the original files (as opposed to additions of modified versions). For example, #3 could be accepted and merged in quickly. #1 and #2 will need more work, but they could be tracked and discussed independently.

[aldalahmeh]

Case (c) is implemented, as you noted earlier, in doSEmod in line 147, where the generated power and demand power are accounted for.

Unless I'm mistaken, line 147 in doSEmod relates to generation measurements, where the load is assumed to be known. So it's not (c), but rather a specific case of (d) with a specific way to handle it (assume load is known).

I believe the issue here is the way the mpc is coded, which actually encodes the power network including the generators and loads, it does not allow having multiple generators nor loads.

I'm afraid I don't follow. The mpc does allow for a single load and multiple generators at each bus.

So, I think we have two options. The first is to simply document the current assumptions and go with your modified code. The second is to try to do something more reasonable with each of the special cases above, and document that. I'm open to either, but don't have the time to do it myself.

I believe is difficult to change the way mpc is coded, so I suggest to code a layer over the mpc to accommodate (d) and (e) cases above, which adds the generators (loads) into a single one and then encoded in the mpc structure.

I'm sorry, but again I don't follow what you mean about "the way mpc is coded" or your proposed solution.

I think we branched in the discussion here. Going back to beginning, the original code is hard coded for a subset of the states and the proposed code extended it to the general case. Let us consider this case and then look at other extensions in the future.

4. I could be wrong, but given my understanding of LSE from other sources, I believe my previous post is correct and the book has a typo. I think it should say that the "inverse of the covariance matrix of the measurements" is ... and (6.157) should say "R^-1 = ...".

I agree with you that the book's notation is rather confusing. In LSE literature a weighing matrix is used and usually denoted as W. In our case, let W = diag(1./sigma.^2), which is actually the inverse of the covariance matrix. And as you correctly noted, the book denoted R= diag(1./sigma.^2), which usually denotes the covariance matrix causing the confusion here.
So I suggest changing the matrix R used in the code to W.

Let's just stick with the original code, where Rinv, the inverse of the variance/covariance matrix, was computed and used correctly as the LSE weight matrix.

Very well.

The lines originally have charging capacitance and there is no reason not to accommodate it, as the original code did. Unless it is stated that the original code only addresses such special case.

Just to confirm, you agree that the original code is correct in this regard?

I'm afraid that I do not agree, if the code was correct, then it would produce the correct solution as in the book's Example 6.17. And it is not a bug in the data since the way the data is coded works perfectly with Example 3.9.

[aldalahmeh]

As the 5 points you mentioned, you cannot add variance in #3 without adding the actual measurement in point #1.

[rdzman]

I'm sorry. I'm not trying to be difficult. I just need to understand clearly what each change actually does. Let's focus on one point at a time.

First, about item 3, you said ...

As the 5 points you mentioned, you cannot add variance in #3 without adding the actual measurement in point 1.

What I meant by point 3 was that, in the original code, each element of the sigma input is a scalar, defining the standard deviation of all measurements of that type (one standard deviation per measurement type). In your version, the elements of the sigma input can be vectors, with a different value for each individual measurement. It is simply implementing a new option for specifying the standard deviations of the existing measurements.

Is there something that prevents this change from being done without making any changes to the set of measurements?

[aldalahmeh]

I appreciate your time and effort in revising this PR. I agree, let us go through all the points.

In response to your comment:

What I meant by point 3 was that, in the original code, each element of the sigma input is a scalar, defining the standard deviation of all measurements of that type (one standard deviation per measurement type). In your version, the elements of the sigma input can be vectors, with a different value for each individual measurement. It is simply implementing a new option for specifying the standard deviations of the existing measurements.

Is there something that prevents this change from being done without making any changes to the set of measurements?

No, there nothing that prevents that. The original code used the same variance for all the PD (or any kind of measurements) whereas in the modified code each measurement (e.g. different PDs at different buses) can have a different variance.

If I may take us to the next point, which is I would like to note that in the original code the demand real and reactive (PD & QD) where not included, whereas they are in the modified code. I believe it is a reasonable generalization.

[rdzman]

Ok, I think we agree on point 3. Would you be willing to create a new PR with just those changes? As I said, that's very straightforward and could be merged in quickly.

On the next point you raise (point 1) about adding load measurements, let's take the discussion to #2.

[rdzman]

And let's take the discussion of point 5 to #3.