[FEATURE] Numerically stable source injection in power flow.

power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp

@@ -190,6 +190,7 @@ template <symmetry_tag sym_type> struct MathModelParam {
     std::vector<BranchCalcParam<sym>> branch_param;
     ComplexTensorVector<sym> shunt_param;
     ComplexTensorVector<sym> source_param;
+    std::vector<std::pair<DoubleComplex, DoubleComplex>> source_param_y0_y1;
 };
 
 struct MathModelParamIncrement {

power_grid_model_c/power_grid_model/include/power_grid_model/common/three_phase_tensor.hpp

@@ -284,6 +284,16 @@ inline std::pair<DoubleComplex, DoubleComplex> inv_sym_param(DoubleComplex const
     return {(s + m) * det_1, -m * det_1};
 }
 
+// inverse of symmetric tensor
+template <symmetry_tag sym> inline ComplexTensor<sym> inv_sym_tensor(ComplexTensor<sym> const& x) {
+    if constexpr (is_symmetric_v<sym>) {
+        return 1.0 / x;
+    } else {
+        auto [s, m] = inv_sym_param(x(0), x(1));
+        return ComplexTensor<sym>(s, m);
+    }
+}
+
 // is nan
 template <class Derived> inline bool is_nan(Eigen::ArrayBase<Derived> const& x) { return x.isNaN().all(); }
 inline bool is_nan(std::floating_point auto x) { return std::isnan(x); }

power_grid_model_c/power_grid_model/include/power_grid_model/component/source.hpp

@@ -42,7 +42,9 @@ class Source : public Appliance {
         y0_ref_ = y1_ref_ / z01_ratio;
     }
 
-    // getter for calculation param, y_ref
+    // getter for y0 y1 ref
+    std::pair<DoubleComplex, DoubleComplex> get_y0_y1() const { return {y0_ref_, y1_ref_}; }
+
     template <symmetry_tag sym> ComplexTensor<sym> math_param() const {
         // internal element_admittance
         if constexpr (is_symmetric_v<sym>) {

power_grid_model_c/power_grid_model/include/power_grid_model/main_model_impl.hpp

@@ -874,6 +874,7 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
             math_param[i].branch_param.resize(state_.math_topology[i]->n_branch());
             math_param[i].shunt_param.resize(state_.math_topology[i]->n_shunt());
             math_param[i].source_param.resize(state_.math_topology[i]->n_source());
+            math_param[i].source_param_y0_y1.resize(state_.math_topology[i]->n_source());
         }
         // loop all branch
         for (Idx i = 0; i != static_cast<Idx>(state_.comp_topo->branch_node_idx.size()); ++i) {
@@ -917,6 +918,8 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
             // assign parameters
             math_param[math_idx.group].source_param[math_idx.pos] =
                 state_.components.template get_item_by_seq<Source>(i).template math_param<sym>();
+            math_param[math_idx.group].source_param_y0_y1[math_idx.pos] =
+                state_.components.template get_item_by_seq<Source>(i).get_y0_y1();
         }
         return math_param;
     }

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/common_solver_functions.hpp

@@ -8,7 +8,9 @@
 #include "y_bus.hpp"
 
 #include "../calculation_parameters.hpp"
+#include "../common/common.hpp"
 #include "../common/exception.hpp"
+#include <numeric>
 
 namespace power_grid_model::math_solver::detail {
 
@@ -57,22 +59,94 @@ template <symmetry_tag sym> inline void copy_y_bus(YBus<sym> const& y_bus, Compl
     });
 }
 
-template <symmetry_tag sym>
-inline void calculate_source_result(IdxRange const& sources, Idx bus_number, YBus<sym> const& y_bus,
-                                    PowerFlowInput<sym> const& input, SolverOutput<sym>& output) {
+// Calculates current and power injection of source i for multiple symmetric sources at a node.
+// The current injection of source i to the bus in phase space is:
+//     i_inj_i = (y_ref_i * z_ref_t) [ (u_ref_i * y_ref_t - i_ref_t) + i_inj_t]
+// Phase space is used because y_ref is a DoubleComplex in the symmetric case.
+// The equation is organized in such a way to avoid the numerical instability induced in the substraction if y_ref_i
+// is large.
+inline void calculate_multiple_source_result(IdxRange const& sources, YBus<symmetric_t> const& y_bus,
+                                             PowerFlowInput<symmetric_t> const& input,
+                                             ComplexValue<symmetric_t> const& i_inj_t,
+                                             SolverOutput<symmetric_t>& output, Idx const& bus_number) {
+    ComplexVector const y_ref = y_bus.math_model_param().source_param;
+    DoubleComplex const y_ref_t =
+        std::accumulate(sources.begin(), sources.end(), DoubleComplex{},
+                        [&](DoubleComplex sum, Idx const source) { return y_ref[source] + sum; });
+    DoubleComplex const i_ref_t =
+        std::accumulate(sources.begin(), sources.end(), DoubleComplex{}, [&](DoubleComplex sum, Idx const source) {
+            return input.source[source] * y_ref[source] + sum;
+        });
+    for (Idx const source : sources) {
+        DoubleComplex const y_ref_i_over_y_ref_t = y_ref[source] / y_ref_t;
+        DoubleComplex const i_inj_i_lhs = y_ref_i_over_y_ref_t * (input.source[source] * y_ref_t - i_ref_t);
+        output.source[source].i = i_inj_i_lhs + (y_ref_i_over_y_ref_t * i_inj_t);
+        output.source[source].s = output.u[bus_number] * conj(output.source[source].i);
+    }
+}
+
+// Calculates current and power injection for source i for multiple asymmetric sources at a node.
+// The current injection of source i to the bus in 012 domain is:
+// i_inj_0 = y_ref_i_0 / y_ref_t_0 * i_inj_t_0
+// i_inj_1 = y_ref_i_1 * [(u_ref_i_1 - i_ref_t_1 / y_ref_t_1) + i_inj_t_1 / y_ref_t_1]
+// i_inj_2 = y_ref_i_2 / y_ref_t_2 * i_inj_t_2
+// Note: u_ref_i_0 = u_ref_i_2 = 0
+// 012 domain is used to solve the numerical instability as for the symmetric case.
+inline void calculate_multiple_source_result(IdxRange const& sources, YBus<asymmetric_t> const& y_bus,
+                                             PowerFlowInput<asymmetric_t> const& input,
+                                             ComplexValue<asymmetric_t> const& i_inj_t,
+                                             SolverOutput<asymmetric_t>& output, Idx const& bus_number) {
+    std::vector<std::pair<power_grid_model::DoubleComplex, power_grid_model::DoubleComplex>> const y0_y1 =
+        y_bus.math_model_param().source_param_y0_y1;
+    ComplexValue<asymmetric_t> const y_ref_t_012 =
+        std::accumulate(sources.begin(), sources.end(), ComplexValue<asymmetric_t>{},
+                        [&](ComplexValue<asymmetric_t> sum, Idx const source) {
+                            sum(0) += y0_y1[source].first;
+                            sum(1) += y0_y1[source].second;
+                            sum(2) = sum(1); // y1 = y2
+                            return sum;
+                        });
+    DoubleComplex const i_ref_1_t =
+        std::accumulate(sources.begin(), sources.end(), DoubleComplex{}, [&](DoubleComplex sum, Idx const source) {
+            return input.source[source] * y0_y1[source].second + sum;
+        });
+    ComplexValue<asymmetric_t> const i_inj_t_012 = dot(get_sym_matrix_inv(), i_inj_t);
+
     for (Idx const source : sources) {
-        ComplexValue<sym> const u_ref{input.source[source]};
-        ComplexTensor<sym> const y_ref = y_bus.math_model_param().source_param[source];
-        output.source[source].i = dot(y_ref, u_ref - output.u[bus_number]);
+        DoubleComplex const y_ref_i_1_over_y_ref_t_1 = y0_y1[source].second / y_ref_t_012[1];
+        DoubleComplex const i_inj_i_1_lhs =
+            y_ref_i_1_over_y_ref_t_1 * (input.source[source] * y_ref_t_012[1] - i_ref_1_t);
+        DoubleComplex const i_inj_i_1 = i_inj_i_1_lhs + (y_ref_i_1_over_y_ref_t_1 * i_inj_t_012(1));
+
+        DoubleComplex const i_inj_i_0 = (y0_y1[source].first / y_ref_t_012[0]) * i_inj_t_012(0);
+        DoubleComplex const i_inj_i_2 = (y0_y1[source].second / y_ref_t_012[2]) * i_inj_t_012(2);
+        ComplexValue<asymmetric_t> const i_inj_012{i_inj_i_0, i_inj_i_1, i_inj_i_2};
+        output.source[source].i = dot(get_sym_matrix(), i_inj_012);
         output.source[source].s = output.u[bus_number] * conj(output.source[source].i);
     }
 }
 
+// Current implementation avoids numerical instability. For details refer to:
+// https://github.com/PowerGridModel/power-grid-model-workshop/blob/experiment/source-calculation/source_calculation/source_calculation.ipynb
+template <symmetry_tag sym>
+inline void calculate_source_result(IdxRange const& sources, Idx const& bus_number, YBus<sym> const& y_bus,
+                                    PowerFlowInput<sym> const& input, SolverOutput<sym>& output,
+                                    ComplexValue<sym> const& i_load_gen_bus) {
+    ComplexValue<sym> const i_inj_t = conj(output.bus_injection[bus_number] / output.u[bus_number]) - i_load_gen_bus;
+    if (sources.size() == 1) {
+        output.source[*sources.begin()].i = i_inj_t;
+        output.source[*sources.begin()].s = output.u[bus_number] * conj(output.source[*sources.begin()].i);
+    } else if (!sources.empty()) {
+        calculate_multiple_source_result(sources, y_bus, input, i_inj_t, output, bus_number);
+    }
+}
+
 template <symmetry_tag sym, class LoadGenFunc>
     requires std::invocable<std::remove_cvref_t<LoadGenFunc>, Idx> &&
              std::same_as<std::invoke_result_t<LoadGenFunc, Idx>, LoadGenType>
 inline void calculate_load_gen_result(IdxRange const& load_gens, Idx bus_number, PowerFlowInput<sym> const& input,
-                                      SolverOutput<sym>& output, LoadGenFunc&& load_gen_func) {
+                                      SolverOutput<sym>& output, LoadGenFunc&& load_gen_func,
+                                      ComplexValue<sym>& i_load_gen_bus) {
     for (Idx const load_gen : load_gens) {
         switch (LoadGenType const type = load_gen_func(load_gen); type) {
             using enum LoadGenType;
@@ -93,6 +167,7 @@ inline void calculate_load_gen_result(IdxRange const& load_gens, Idx bus_number,
             throw MissingCaseForEnumError("Power injection", type);
         }
         output.load_gen[load_gen].i = conj(output.load_gen[load_gen].s / output.u[bus_number]);
+        i_load_gen_bus += output.load_gen[load_gen].i;
     }
 }
 
@@ -114,11 +189,12 @@ inline void calculate_pf_result(YBus<sym> const& y_bus, PowerFlowInput<sym> cons
     output.load_gen.resize(load_gens_per_bus.element_size());
     output.bus_injection.resize(sources_per_bus.size());
 
+    output.bus_injection = y_bus.calculate_injection(output.u);
     for (auto const& [bus_number, sources, load_gens] : enumerated_zip_sequence(sources_per_bus, load_gens_per_bus)) {
-        calculate_source_result<sym>(sources, bus_number, y_bus, input, output);
-        calculate_load_gen_result<sym>(load_gens, bus_number, input, output, load_gen_func);
+        ComplexValue<sym> i_load_gen_bus{};
+        calculate_load_gen_result<sym>(load_gens, bus_number, input, output, load_gen_func, i_load_gen_bus);
+        calculate_source_result<sym>(sources, bus_number, y_bus, input, output, i_load_gen_bus);
     }
-    output.bus_injection = y_bus.calculate_injection(output.u);
 }
 
 template <symmetry_tag sym>

tests/data/power_flow/source-big-sk/input.json

@@ -39,15 +39,17 @@
         "status": 0,
         "u_ref": 1.0,
         "sk": 2e50,
-        "rx_ratio": 0.1
+        "rx_ratio": 0.1,
+        "z01_ratio": 1.0
       },
       {
         "id": 7,
         "node": 1,
         "status": 0,
         "u_ref": 1.0,
         "sk": 3e50,
-        "rx_ratio": 0.1
+        "rx_ratio": 0.1,
+        "z01_ratio": 1.0
       }
     ],
     "sym_load": [

tests/data/power_flow/source-big-sk/params.json

@@ -1,8 +1,5 @@
 {
   "calculation_method": ["newton_raphson", "iterative_current"],
   "rtol": 1e-6,
-  "atol": 1e-6,
-  "fail": {
-    "reason": "Source injection calculation #458"
-  }
+  "atol": 1e-6
 }