Using pivot perturbation for ill-conditioned state estimation

[TonyXiang8787]

Closes #853

Mathematical approach

This PR attempts to fix the problem by using pivot perturbation. This is an alternative method than partial/full pivoting method in the LU factorization. In this way, the structure of the sparse matrix is not changed. The methods are described in [Li99] and [Schenk04].

This PR applies the methods with following aspects:

1. The perturbation threshold is set to $\epsilon = 10^{-13}$, same as used in MKL PARDISO. This applied for both perturbed pivot and iterative refinement.
2. The infinite-norm of the matrix is calculated block-wise. First the infinite norm of each block is calculated, then the norm of the whole matrix by treating each block as a single value. When calculating the whole matrix, the diagonal blocks are ignored. In Power System State Estimation, we sometimes have huge diagonal values. We do not want those values to influence the norm of the matrix.
3. The backward error is calculated by capping the denominators to a minimum to avoid scenarios where both $|x|$ and $|b|$ are zero.

$$ D_{\max} = \max_i {(|A||x| + |b|)_i} $$

$$berr = \max_{ i } {\frac{|r|_i}{\max ( (|A||x| + |b|)_i, \epsilon_{berr} D_{\max}) }}$$ $$\epsilon_{berr} = 10^{-4}$$

5. The iterative refinement is capped at 5 times. The convergence criteria is $berr < \epsilon$.

Implementation

We adjust LUSparseSolver to allow optionally pivot perturbation. To do this, the original matrix has to be cached to do the iterative refinement. Optional member data are defined in the class to store the cache. This perturbation is an opt-in feature which should only be enabled if we suspect ill-conditioned scenario. Currently, it is only enabled in State Estimation (both NR and IL) when the following two conditions are met:

1. The system is sufficiently observable. We need this to make sure that we do not create fake observability by perturbation.
2. The system could be possibly ill-conditioned. That means the following conditions are met:
   1. The bus does not have (zero) injection measurement.
   2. All upstream branches, i.e. branches between this bus and the buses with higher bus numbers, do not have branch measurement.
   3. The bus does not have a phasor voltage measurement.

Check list

• Implement the core
• Unit test
• Extensive validation test
• Documentation
• Run Nightly

[figueroa1395]

Since this is both a big and a very important PR, this is a friendly reminder to run the Nightly build manually once it is "ready" to be merged.

[mgovers]

For future reference purposes only

• important changes to mathematics after the initially proposed PR:
  • back error calculation: Using pivot perturbation for ill-conditioned state estimation #862 (comment)
• Other discussions worth of reference:
  • absence of early-out found in literature Using pivot perturbation for ill-conditioned state estimation #862 (comment)
  • element-wise absolute value of matrix: Using pivot perturbation for ill-conditioned state estimation #862 (comment)
  • block-wise off-diagonal infinite norm: Using pivot perturbation for ill-conditioned state estimation #862 (comment)
  • different implementation of iteration: Using pivot perturbation for ill-conditioned state estimation #862 (comment)

[TonyXiang8787]

Nightly build in https://github.com/PowerGridModel/power-grid-model/actions/runs/12924710668

[mgovers]

Running the benchmark yields comparatively similar results. I found no significant regressions. In fact, the new implementation often appears to give slightly better results, in particular for ClangCL, but it is most most likely within run differences on my machine and therefore not significant.

output.zip

[mgovers]

Profiling results

using the Visual Studio 2022 Profiler on msvc-release on power_grid_model_benchmark_cpp

Very similar results between main and this branch up to the IterativePFLUSolver (Root in the profiles below), and even into the SparseLUSolver. After that, the implementation changes a lot and therefore, the implementation is also different.

Note: the threaded version has almost identical behavior, as expected; see included .diagsession files

main

this branch

Diagsessions

msvc-release-profile-report.zip
s/assets/61aeb8e5-e172-4d3c-bc42-b8d03dd72ff8)

[TonyXiang8787]

Full run in https://github.com/PowerGridModel/power-grid-model/actions/runs/12945646002

[sonarqubecloud]

Quality Gate passed

Issues
3 New issues
0 Accepted issues

Measures
0 Security Hotspots
99.6% Coverage on New Code
0.4% Duplication on New Code

See analysis details on SonarQube Cloud

[mgovers]

I hereby approve the current set of changes