95 views (last 30 days)

John D'Errico
on 18 Mar 2017

Is that the right way to do it? Very possibly there is no good way. If you are getting warning messages, that generally means your matrix is numerically singular. Finding the inverse of a numerically singular matrix will not be well posed, no matter what computation you use.

So the very first thing you need to do is test the condition number of the matrix. If it is truly very large and sparse, then condest may be the best tool, to give at least an estimate of the approximate condition number. Of course, a lot of people think their matrices are large and sparse, when they are neither truly large or truly sparse. So cond may suffice for you, to tell you if the matrix is singular. Again, if your matrix is singular, then you are wasting your time to compute the diagonal of the inverse, since the inverse matrix will be numerical garbage.

A better solution may depend on how the matrix was created, using a little mathematics. But that is something we are not able to know, since you have told us nothing of value.

John D'Errico
on 20 Mar 2017

As is often the case, people think they have large sparse matrices, when they don't. :)

No matter how sparse it is, a 22x22 matrix is not large. Not even worth using sparse storage to store it. The point is, just make it a full matrix. Things get easier then.

YB = full(Ybus_pu);

>> rank(YB)

ans =

21

>> cond(YB)

ans =

3.22109878812164e+17

The matrix is singular. Many people don't understand what that means. The condition number is roughly 3e17. What that means is if you try to solve a linear system of equations, OR compute the inverse matrix, the system will amplify any noise in your problem by roughly a factor of 3e17.

Is there noise in your problem? YES, there is! The noise comes from how those numbers are stored. They have random junk in the least significant bits of the numbers. That last bit will be corrupted, even if the numbers themselves were computed with no "error". And ANY floating point computations end up corrupting those least significant bits. Even just the process of solving for the inverse.

So accept that there is junk in your matrix entries down in the least significant bits, that is on the order of eps*YB(i,j). In double precision, eps is:

eps

ans =

2.22044604925031e-16

But remember that the condition number of your matrix, thus the extent of any amplification of the noise, is 3e17.

What does this tell you? It says that the elements of the inverse are complete junk. They will be completely corrupted by the noise in those least significant bits of the matrix.

I know. you don't believe me. Lets do a little test. First, compute the diagonal elements of the inverse matrix directly.

diag(inv(YB))

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.161271e-18.

ans =

-6991140573317.74 + 2649900922810.41i

-6991140537171.36 + 2649900963812.15i

-6991140537171.34 + 2649900963812.15i

-6991140537171.26 + 2649900963812.17i

-6991140537171.32 + 2649900963812.16i

-6991140537171.28 + 2649900963812.18i

-6991140537171.24 + 2649900963812.2i

-6991140537171.21 + 2649900963812.21i

-6991140537171.22 + 2649900963812.21i

-6991140537171.18 + 2649900963812.22i

-6991140537171.3 + 2649900963812.17i

-6991140537171.29 + 2649900963812.18i

-6991140537171.21 + 2649900963812.2i

-6991140537171.27 + 2649900963812.19i

-6991140537171.25 + 2649900963812.2i

-6991140537171.24 + 2649900963812.22i

-6991140537171.21 + 2649900963812.23i

-6991140537171.24 + 2649900963812.2i

-6991140537171.23 + 2649900963812.23i

-6991140537171.21 + 2649900963812.24i

-6991140537171.18 + 2649900963812.25i

-6991140573317.74 + 2649900922810.41i

Now, perturb the matrix elements by a TINY amount, on the order of eps.

diag(inv(YB + randn(size(YB)).*eps(YB)))

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.469158e-18.

ans =

1341274048446.89 + 46993008610.608i

1341274044976.97 + 46993074953.964i

1341274044976.98 + 46993074953.9822i

1341274044977.06 + 46993074953.9981i

1341274044976.99 + 46993074954.0002i

1341274044977.03 + 46993074954.0155i

1341274044977.06 + 46993074954.0322i

1341274044977.09 + 46993074954.0497i

1341274044977.09 + 46993074954.0424i

1341274044977.12 + 46993074954.0529i

1341274044977.01 + 46993074954.0161i

1341274044977.02 + 46993074954.0311i

1341274044977.09 + 46993074954.0471i

1341274044977.03 + 46993074954.0459i

1341274044977.04 + 46993074954.0592i

1341274044977.04 + 46993074954.0724i

1341274044977.08 + 46993074954.0829i

1341274044977.06 + 46993074954.057i

1341274044977.06 + 46993074954.0869i

1341274044977.08 + 46993074954.0947i

1341274044977.1 + 46993074954.1026i

1341274048446.89 + 46993008610.608i

The tiny permutations in those elements results in crap that was as large as the elements of the original inverse.

Essentially, if you think of this as a signal to noise thing, there is NO signal remaining in the elements of that inverse.

It does not matter how much you want to compute the elements of the inverse matrix when it is singular. The numbers you will produce are COMPLETELY MEANINGLESS. There is no information content remaining.

Ok, some might now say, but a 22x22 matrix is small. Just use the symbolic toolbox. The condition number is still 3e17. And the noise in your elements is of the same magnitude, because they are created in double precision. That means you will still see amplification of that noise by roughly the condition number.

Sorry, but you can't succeed via that route.

Ok, suppose you go back to the original matrix, and created it in full symbolic form. So never go through double precision. Now all the entries of the matrix are symbolic, and have no corruption in the least significant bits. Can we possibly now survive?

The question is why is your matrix singular. If I look at the singular values, of this thing, I see what is one effectively zero singular value. (Compare it to the largest singular value. It is relatively near eps.)

svd(YB)

ans =

5235.34878600224

224.456140741978

208.213133465402

183.772675608031

144.79532279094

122.428521425073

117.857092301615

100.073512943469

71.6504865943657

66.4367745517612

49.3938713012852

43.2493963956952

35.2213230481699

21.088537030153

18.9310735629174

16.497370799831

14.506749928581

8.95375281077995

4.90625096086107

2.01074591372002

9.2676852076552e-06

1.62533009087101e-14

So I have no idea how that matrix was generated. It may well be that even if you built it in symbolic form, it would still be singular!

Walter Roberson
on 18 Mar 2017

Walter Roberson
on 20 Mar 2017

If you have the symbolic toolbox, then you can proceed symbolically:

dinv = diag( inv( sym(Ybus_pu) ) );

double(dinv)

The values are mostly close to -8327187525072.06 + 2366252476427.26i with the "ones" and the decimals varying -- the first 12 places are pretty constant for most of the entries.

Walter Roberson
on 20 Mar 2017

If you take

Y = sym(Ybus_pu);

Y1 = Y;

Y1(1,1) = Y1(1,1) + 8.11130830789689e-14;

dinv1 = diag(inv(Y1));

r41 = real(dinv1(4));

Y2 = Y;

Y2(1,1) = Y2(1,1) + 1.41747416292681e-13;

dinv2 = diag(inv(Y2));

r42 = real(dinv2(4));

then r41 will be about -15814803937051 and r42 will be about 15828133351471 . This indicates that a change of 1E-14 to 1E-13 can change the sign of the result completely.

But

eps(real(Ybus_pu))

is 4.54747350886464e-13 . which is about 4 to 8 times larger than those shifts.

This tells us that the answers you get out through the process are essentially numeric garbage, completely different with a variation in values in the input smaller than MATLAB double precision can represent.

Your situation is hopeless unless you can generate those bus values to higher precision such as by using the Symbolic Toolbox when you create them.

Opportunities for recent engineering grads.

Apply TodayFind the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
## 0 Comments

Sign in to comment.