Sparse matrix solver: Trigger pivot perturbation based on each possibly ill-conditioned one

[figueroa1395]

This PR aims to improve the pivot perturbation scheme currently used, based on a finer grained criteria.

Before, we would judge if the whole matrix was ill-conditioned (prior to performing the LU decomposition), whereas with this PR we aim to judge the ill-condition on a per pivot basis. Furthermore, because of this we had to use the off diagonal infinite norm of the matrix to judge the pivot's size in order to avoid introducing fake observability, but now, we would be able to use the full infinite norm without having that issue.

In addition, with this PR we aim to resolve the edge cases reported in #1118 (to be solved only once the full observability check is in place, see #1085), as well as make pivot perturbation more robust for radial grids. The scheme used is described in more detail in #1118 (comment).

[figueroa1395]

Converting back to draft as the validation cases created in #1118 are still unresolved if I force in the perturbation... Several things happen:

1. The ill-conditioned check is not strong enough, as the validation cases in Sparse matrix error: validation cases #1118 have enough sensors. In order for this to work after meshed observability is implemented, we would first need to come up with a stronger ill-condition.
2. Even if I manually mark the pivots as ill-conditioned, they aren't considered small enough according to the current threshold (related to the matrix norm), so the pivot perturbation is not triggered still. This raises questions and is tied to the following point.
3. Even if the pivot is not considered small, we are still triggering SparseMatrixError after Gaussian elimination. This other threshold is based on the maximum pivot present in the matrix after the transformations are performed. The problem here is that we end up indirectly comparing two distinct things: the norm of the original matrix (used to determine if a pivot perturbation is to be triggered) with the max pivot of the matrix after several transformations have been done.
   This is problematic because even though the norm of the matrix (including the diagonal) should always be bigger than the largest pivot for the original matrix, we are now comparing it with the max pivot of a transformed matrix for which the pivots can become a lot larger than the original norm. This is the edge case we are currently running into. So we end up doing a comparison that isn't fair.

To summarize, this PR still has value, as it improves how we do perturbation (more fine grained and the whole matrix norm is now taken into account without generating fake observability), but it doesn't resolve (at least yet) the edge cases we are running into. I don't have a concrete proposal on how to go forward, so thoughts are welcome, but I will try some experimentation with how we determine the max pivot (maybe calculating it from the original matrix makes more sense, but it still doesn't seem satisfactory) and comment back here if I find any relevant results.

[figueroa1395]

c1fa814 and a9f0158 is a quick and dirty implementation of an idea I had to address this issue. This idea is actually independent of pivot perturbation - hence it is independent of the meshed-ness of the grid and the observability check.

In short, I replaced the max_pivot (previously calculated after gaussian elimination has taken place), with one that is calculated based on the original matrix. Now the validation tests cases seem to pass locally (but other test cases based on big sk are now broken). Leaving here to collect thoughts in case it is worth to keep investigating.

Edit: This doesn't work for ill-conditioned cases generated by the opposite compared to the ones in the validation tests, i.e. for cases in which the source short circuit voltage is extremely high (very large sk). This idea is scratched.

[figueroa1395]

The final idea from this exploration is applied in 5d379a9. In short, the threshold to raise SparseMatrixError when the matrix is ill-conditioned has been split into two: relative and absolute thresholds.

• Relative threshold: It replaces the original, epsilon * max_pivot where max_pivot is given by the max pivot in the dense block after gaussian elimination. The new threshold is given by epsilon * matrix_norm where matrix_norm is calculated based in the whole original matrix (before any factorization or gaussian decomposition takes place). I will give a somewhat heuristic explanation for why I think this strategy makes sense.

The condition number $\kappa(\mathbf{A})$ of a matrix $\mathbf{A}$ is defined as

$$\kappa(\mathbf{A}) = \|\mathbf{A}\| \times \|\mathbf{A^{-1}}\| .$$

Considering we are interested in the LU decomposition of the matrix, we can re-write the above as

$$\kappa(\mathbf{A}) = \|\mathbf{A}\| \times \|\mathbf{U^{-1}} \times \mathbf{L^{-1}}\| \leq \|\mathbf{A}\| \times \|\mathbf{U^{-1}} \| \times \| \mathbf{L^{-1}}\| \approx \|\mathbf{A}\| \times \|\mathbf{U^{-1}} \| ,$$

as $\mathbf{L^{-1}}$ has 1's in the diagonal and is usually well behaved. Considering now

$$\|\mathbf{U^{-1}}\| \sim 1 / \text{min pivot} ,$$

as that is the element that affects the norm the most, we can approximate the condition number by

$$\kappa(\mathbf{A}) \lesssim \|\mathbf{A}\| \times 1 / \text{min pivot} ,$$

and for the matrix to not be ill-conditioned, it needs to satisfy

$$\kappa(\mathbf{A}) < 1 / \epsilon ,$$

which translates into the introduced relative threshold check

$$\text{min pivot} < \|\mathbf{A}\| \times \epsilon .$$

This works with the norm of the original matrix two ways. The condition number is defined over the original matrix, as this is the original problem we are trying to solve, and we homogenize the perturbation threshold with this one. This also takes care of both extreme ill-conditioned cases, when sk is very large or when we have links.

• Absolute threshold: This one just checks against machine epsilon. This is by no means needed, but I still thought it is worth it. Furthermore, all tests pass with this one removed.

With the introduction of the above, the introduced ill condition validation cases pass, while the old ones do as well.

Note 1: This still makes use of the (not so strong) "topological" check of possibly ill conditioned pivots during the observability check, as the full norm of the matrix (including its diagonal) is to be used. This shouldn't be problem once meshed observability is introduced, as it can be made into a standalone routine.

Note 2: I will cleanup the PR only after input is received, as this is very much experimental from my side.

[TonyXiang8787]

@figueroa1395 if you use relative threshold to raise SparseMatrixError, does it mean the matrix norm has to be always calculated? Even for power flow equations where you do not foresee ill-conditioned pivots at all?

[figueroa1395]

@figueroa1395 if you use relative threshold to raise SparseMatrixError, does it mean the matrix norm has to be always calculated? Even for power flow equations where you do not foresee ill-conditioned pivots at all?

@TonyXiang8787 The norm is only calculated with initialize_pivot_perturbation(), which is only called when there are possibly ill conditioned pivots as determined (on a per pivot basis) at the observability level, hence the norm is only relevant for state estimation.

In short, no. Power flow calculations should not be affected.

Note: I just manually checked, and the power flow validation cases never trigger it.

[sonarqubecloud]

Quality Gate failed

Failed conditions
6.2% Duplication on New Code (required ≤ 3%)

See analysis details on SonarQube Cloud

[Jerry-Jinfeng-Guo]

do-not-merge added due to the experimental nature of the PR.