Add flat start option for Newton-Raphson power flow

docs/algorithms/pf-algorithms.md

@@ -127,6 +127,31 @@ For each iteration the following steps are executed:
 - Using LU decomposition, solve $J(i) \Delta x(i)  =  \Delta y(i)$ for $\Delta x(i)$
 - Compute $x(i+1)$ from $\Delta x(i) =  x(i+1) - x(i)$
 
+### Initialization
+
+The start $x(0)$ is set by the `calculation_initialization` option
+({py:class}`CalculationInitialization <power_grid_model.enum.CalculationInitialization>`):
+
+- `default`: the default initialization, which is `linear` for Newton-Raphson.
+- `linear`: a linear voltage guess.
+  Every load and generator is replaced by the constant admittance $-\overline{S}$ at 1 p.u. (the specified reactive
+  power of a regulated generator is left out), and the resulting linear network is solved once.
+- `flat`: every node at 1 p.u. with the reference angle of the sources (`u_ref_angle`, the angle of their average
+  voltage if there are several) plus its topological phase shift, a node with a source at the source's reference
+  voltage.
+- `average_source`: every node at the average reference voltage of all sources, with its topological phase shift.
+  This is the start of the [iterative current](#iterative-current-power-flow) power flow.
+
+In all cases, voltage regulated nodes then start at their reference magnitude `u_ref`, keeping the angle of the
+start.
+
+The linear guess is a good start where load and generation are close to the source.
+In a meshed transmission grid with much of the generation far from the source, it replaces that generation by
+negative conductances, and the guess can lie outside the region where Newton-Raphson converges.
+The calculation then diverges, or stops at a singular Jacobian, although the power flow has a solution that a flat
+start reaches.
+A flat start is also the common start of other power flow tools, which makes results and iteration counts comparable.
+
 ### PV nodes and reactive-power limits
 
 In Newton-Raphson power flow, an active `voltage_regulator` replaces the reactive-power equation of the regulated node

docs/user_manual/calculations.md

@@ -340,6 +340,9 @@ find the range of loading conditions that are relevant for your use case and the
 only use linear methods within this range for the specific grid configuration.
 
 Non convergence of newton raphson is a good signal of unpractical or unfeasible systems.
+In a meshed grid with much of the generation far from the source, it can also come from the start: try
+`calculation_initialization=CalculationInitialization.flat` (see
+[Initialization](../algorithms/pf-algorithms.md#initialization)) before concluding that the system has no solution.
 This signal can be ignored when using linear methods.
 Similarly, having atleast some results from linear methods can aid in finding data errors or the reason
 for non convergence of newton raphson method.

power_grid_model_c/power_grid_model/include/power_grid_model/calculation.hpp

@@ -80,9 +80,14 @@ template <typename T, typename sym> struct Calculator;
 
 template <symmetry_tag sym> struct Calculator<power_flow_t, sym> {
     template <typename State>
-    static auto preparer(State const& state, ComponentToMathCoupling& /*comp_coup*/,
-                         MainModelOptions const& /*options*/) {
-        return [&state](Idx n_math_solvers) { return main_core::prepare_power_flow_input<sym>(state, n_math_solvers); };
+    static auto preparer(State const& state, ComponentToMathCoupling& /*comp_coup*/, MainModelOptions const& options) {
+        return [&state, initialization = options.calculation_initialization](Idx n_math_solvers) {
+            auto input = main_core::prepare_power_flow_input<sym>(state, n_math_solvers);
+            for (auto& math_input : input) {
+                math_input.initialization = initialization;
+            }
+            return input;
+        };
     }
     static auto solver(CalculationMethod calculation_method, MainModelOptions const& options, bool cache_run,
                        Logger& logger) {

power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp

@@ -279,6 +279,7 @@ template <symmetry_tag sym_type> struct PowerFlowInput {
     ComplexValueVector<sym> s_injection; // Specified injection power of each load_gen
     std::vector<VoltageRegulatorCalcParam<sym>> voltage_regulator;
     IntSVector load_gen_status;
+    CalculationInitialization initialization{CalculationInitialization::default_initialization};
 };
 
 template <symmetry_tag sym_type> struct StateEstimationInput {

power_grid_model_c/power_grid_model/include/power_grid_model/common/enum.hpp

@@ -109,6 +109,15 @@ enum class FaultPhase : IntS {
 
 enum class ShortCircuitVoltageScaling : IntS { minimum = 0, maximum = 1 };
 
+// Initialization of the calculation, only applicable if supported by the calculation type and method.
+// For the Newton-Raphson power flow:
+// - linear (default): a linear guess with loads and generators as constant admittances
+// - flat: every bus at 1 p.u. with the source angle and its topological phase shift, source buses at their reference
+//   voltage
+// - average_source: every bus at the average reference voltage of all sources with its topological phase shift, as
+//   used by the iterative current power flow
+enum class CalculationInitialization : IntS { default_initialization = 0, flat = 1, linear = 2, average_source = 3 };
+
 enum class CType : IntS { c_int32 = 0, c_int8 = 1, c_double = 2, c_double3 = 3 };
 
 enum class SerializationFormat : IntS { json = 0, msgpack = 1 };

power_grid_model_c/power_grid_model/include/power_grid_model/main_model_fwd.hpp

@@ -32,6 +32,7 @@ struct MainModelOptions {
     Idx threading{sequential};
 
     ShortCircuitVoltageScaling short_circuit_voltage_scaling{ShortCircuitVoltageScaling::maximum};
+    CalculationInitialization calculation_initialization{CalculationInitialization::default_initialization};
 };
 
 } // namespace power_grid_model

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_current_pf_solver.hpp

@@ -94,7 +94,7 @@ class IterativeCurrentPFSolver : public IterativePFSolver<sym_type, IterativeCur
     // Add source admittance to Y bus and set variable for prepared y bus to true
     void initialize_derived_solver(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input,
                                    SolverOutput<sym>& output) {
-        make_flat_start(input, output.u);
+        this->make_average_source_start(input, output.u);
 
         auto const& sources_per_bus = this->sources_per_bus_.get();
         IdxVector const& bus_entry = y_bus.lu_diag();
@@ -203,26 +203,6 @@ class IterativeCurrentPFSolver : public IterativePFSolver<sym_type, IterativeCur
                                       ComplexValue<sym>{input.source[source_number]});
         }
     }
-
-    void make_flat_start(PowerFlowInput<sym> const& input, ComplexValueVector<sym>& output_u) {
-        std::vector<double> const& phase_shift = this->phase_shift_.get();
-        // average u_ref of all sources
-        DoubleComplex const u_ref = [&]() {
-            DoubleComplex sum_u_ref = 0.0;
-            for (auto const& [bus, sources] : enumerated_zip_sequence(this->sources_per_bus_.get())) {
-                for (Idx const source : sources) {
-                    sum_u_ref += input.source[source] * std::exp(1.0i * -phase_shift[bus]); // offset phase shift
-                }
-            }
-            return sum_u_ref / static_cast<double>(input.source.size());
-        }();
-
-        // assign u_ref as flat start
-        for (Idx i = 0; i != this->n_bus_; ++i) {
-            // consider phase shift
-            output_u[i] = ComplexValue<sym>{u_ref * std::exp(1.0i * phase_shift[i])};
-        }
-    }
 };
 
 } // namespace iterative_current_pf

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_pf_solver.hpp

@@ -93,6 +93,29 @@ template <symmetry_tag sym, typename DerivedSolver> class IterativePFSolver {
                                     [this](Idx i) { return (load_gen_type_.get())[i]; });
     }
 
+    // Average reference voltage of all sources, with the topological phase shift of the source buses offset
+    DoubleComplex average_source_voltage(PowerFlowInput<sym> const& input) const {
+        std::vector<double> const& phase_shift = phase_shift_.get();
+        DoubleComplex sum_u_ref = 0.0;
+        for (auto const& [bus, sources] : enumerated_zip_sequence(sources_per_bus_.get())) {
+            for (Idx const source : sources) {
+                sum_u_ref += input.source[source] * std::exp(1.0i * -phase_shift[bus]); // offset phase shift
+            }
+        }
+        return sum_u_ref / static_cast<double>(input.source.size());
+    }
+
+    // Initialize every bus to the average reference voltage of all sources, with the topological phase shift of each
+    // bus accounted for
+    void make_average_source_start(PowerFlowInput<sym> const& input, ComplexValueVector<sym>& output_u) const {
+        std::vector<double> const& phase_shift = phase_shift_.get();
+        DoubleComplex const u_ref = average_source_voltage(input);
+        for (Idx i = 0; i != n_bus_; ++i) {
+            // consider phase shift
+            output_u[i] = ComplexValue<sym>{u_ref * std::exp(1.0i * phase_shift[i])};
+        }
+    }
+
   private:
     Idx n_bus_;
     std::reference_wrapper<DoubleVector const> phase_shift_;

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_pf_solver.hpp

@@ -246,12 +246,8 @@ class NewtonRaphsonPFSolver : public IterativePFSolver<sym_type, NewtonRaphsonPF
           voltage_regulators_per_load_gen_{std::ref(topo.voltage_regulators_per_load_gen)},
           clamped_regulators_per_load_gen_(topo.load_gen_type.size(), RealValue<sym>{nan}) {}
 
-    // Initialize the unknown variable in polar form using a linear voltage guess solved in the real domain.
-    // This implementation reuses the class-level Jacobian matrix (data_jac_) and RHS vector (del_x_pq_)
-    // to eliminate temporary memory allocations.
-    // The complex system (G + jB)(Ur + jUi) = Ir + jIi is mapped to the real system:
-    // [[G, -B], [B, G]] [Ur, Ui]^T = [Ir, Ii]^T
-    // This mapping ensures that G aligns with the N/M sub-blocks as in the NR Jacobian.
+    // Initialize the unknown variable in polar form, from a linear voltage guess or a flat start
+    // (PowerFlowInput::initialization). Either way, PV buses then start at their reference magnitude.
     void initialize_derived_solver(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input,
                                    SolverOutput<sym>& output) {
         // Reset reused buffers to zero
@@ -265,6 +261,54 @@ class NewtonRaphsonPFSolver : public IterativePFSolver<sym_type, NewtonRaphsonPF
         const bool has_usable_limits = set_bus_types_and_q_limits(input);
         limit_check_countdown_ = has_usable_limits ? limit_check_at_iteration : no_limit_check;
 
+        switch (input.initialization) {
+            using enum CalculationInitialization;
+        case flat:
+            make_flat_start(input, output.u);
+            break;
+        case average_source:
+            this->make_average_source_start(input, output.u);
+            break;
+        default:
+            make_linear_start(y_bus, input, output.u);
+            break;
+        }
+
+        set_reference_voltage_for_pv_buses(output.u);
+
+        // get magnitude and angle of start voltage
+        for (Idx i = 0; i != this->n_bus_; ++i) {
+            x_[i].v() = cabs(output.u[i]);
+            x_[i].theta() = arg(output.u[i]);
+        }
+    }
+
+    // Flat start: every bus at 1 p.u. with the reference angle of the sources (the angle of their average voltage) plus
+    // its topological phase shift, a bus with a source at that source's reference voltage (the mean if there are
+    // several). PV buses are set to their reference magnitude afterwards.
+    void make_flat_start(PowerFlowInput<sym> const& input, ComplexValueVector<sym>& u) const {
+        std::vector<double> const& phase_shift = this->phase_shift_.get();
+        double const source_angle = arg(this->average_source_voltage(input));
+        for (auto const& [bus, sources] : enumerated_zip_sequence(this->sources_per_bus_.get())) {
+            if (sources.empty()) {
+                u[bus] = ComplexValue<sym>{std::exp(1.0i * (source_angle + phase_shift[bus]))};
+                continue;
+            }
+            DoubleComplex u_ref{};
+            for (Idx const source : sources) {
+                u_ref += input.source[source];
+            }
+            u[bus] = ComplexValue<sym>{u_ref / static_cast<double>(std::ranges::size(sources))};
+        }
+    }
+
+    // Linear voltage guess solved in the real domain, with every load_gen as a constant admittance.
+    // This implementation reuses the class-level Jacobian matrix (data_jac_) and RHS vector (del_x_pq_)
+    // to eliminate temporary memory allocations.
+    // The complex system (G + jB)(Ur + jUi) = Ir + jIi is mapped to the real system:
+    // [[G, -B], [B, G]] [Ur, Ui]^T = [Ir, Ii]^T
+    // This mapping ensures that G aligns with the N/M sub-blocks as in the NR Jacobian.
+    void make_linear_start(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input, ComplexValueVector<sym>& u) {
         // Map network admittance to real-domain system
         IdxVector const& map_lu_y_bus = y_bus.map_lu_y_bus();
         ComplexTensorVector<sym> const& ydata = y_bus.admittance();
@@ -290,16 +334,7 @@ class NewtonRaphsonPFSolver : public IterativePFSolver<sym_type, NewtonRaphsonPF
         sparse_solver_.prefactorize_and_solve(data_jac_, perm_, del_x_pq_, del_x_pq_);
 
         for (Idx const i : std::views::iota(Idx{}, this->n_bus_)) {
-            output.u[i] =
-                ComplexValue<sym>{RealValue<sym>{del_x_pq_[i].u_real()}, RealValue<sym>{del_x_pq_[i].u_imag()}};
-        }
-
-        set_reference_voltage_for_pv_buses(output.u);
-
-        // get magnitude and angle of start voltage
-        for (Idx i = 0; i != this->n_bus_; ++i) {
-            x_[i].v() = cabs(output.u[i]);
-            x_[i].theta() = arg(output.u[i]);
+            u[i] = ComplexValue<sym>{RealValue<sym>{del_x_pq_[i].u_real()}, RealValue<sym>{del_x_pq_[i].u_imag()}};
         }
     }
 

power_grid_model_c/power_grid_model_c/include/power_grid_model_c/basics.h

@@ -229,6 +229,22 @@ enum PGM_ShortCircuitVoltageScaling {
     PGM_short_circuit_voltage_scaling_maximum = 1, /**< voltage scaling for maximum short circuit currents */
 };
 
+/**
+ * @brief Enumeration of the calculation initializations.
+ *
+ */
+enum PGM_CalculationInitialization {
+    /** default initialization of the selected calculation method */
+    PGM_calculation_initialization_default = 0,
+    /** flat start: every bus at 1 p.u. with the source angle and its phase shift, source buses at their reference
+     * voltage, PV buses at their reference magnitude */
+    PGM_calculation_initialization_flat = 1,
+    /** start from a linear voltage guess (loads and generators as constant admittances) */
+    PGM_calculation_initialization_linear = 2,
+    /** every bus at the average reference voltage of all sources with its phase shift */
+    PGM_calculation_initialization_average_source = 3,
+};
+
 /**
  * @brief Enumeration of tap changing strategies.
  *

power_grid_model_c/power_grid_model_c/include/power_grid_model_c/options.h

@@ -29,6 +29,7 @@ extern "C" {
  *   - max_iter: 20
  *   - threading: -1
  *   - short_circuit_voltage_scaling: PGM_short_circuit_voltage_scaling_maximum
+ *   - calculation_initialization: PGM_calculation_initialization_default
  *   - experimental_features: PGM_experimental_features_disabled
  *
  * @param handle
@@ -113,6 +114,18 @@ PGM_API void PGM_set_threading(PGM_Handle* handle, PGM_Options* opt, PGM_Idx thr
 PGM_API void PGM_set_short_circuit_voltage_scaling(PGM_Handle* handle, PGM_Options* opt,
                                                    PGM_Idx short_circuit_voltage_scaling) PGM_NOEXCEPT;
 
+/**
+ * @brief Specify how the calculation is initialized
+ *
+ * Only applicable if supported for the calculation type and method.
+ *
+ * @param handle
+ * @param opt pointer to option instance
+ * @param calculation_initialization See #PGM_CalculationInitialization
+ */
+PGM_API void PGM_set_calculation_initialization(PGM_Handle* handle, PGM_Options* opt,
+                                                PGM_Idx calculation_initialization) PGM_NOEXCEPT;
+
 /**
  * @brief Specify the tap changing strategy for power flow calculations
  *

power_grid_model_c/power_grid_model_c/src/model.cpp

@@ -191,6 +191,10 @@ constexpr auto get_short_circuit_voltage_scaling(PGM_Options const& opt) {
     return safe_enum<ShortCircuitVoltageScaling>(opt.short_circuit_voltage_scaling);
 }
 
+constexpr auto get_calculation_initialization(PGM_Options const& opt) {
+    return safe_enum<CalculationInitialization>(opt.calculation_initialization);
+}
+
 constexpr auto extract_calculation_options(PGM_Options const& opt) {
     return MainModel::Options{.calculation_type = get_calculation_type(opt),
                               .calculation_symmetry = get_calculation_symmetry(opt),
@@ -200,7 +204,8 @@ constexpr auto extract_calculation_options(PGM_Options const& opt) {
                               .err_tol = opt.err_tol,
                               .max_iter = opt.max_iter,
                               .threading = opt.threading,
-                              .short_circuit_voltage_scaling = get_short_circuit_voltage_scaling(opt)};
+                              .short_circuit_voltage_scaling = get_short_circuit_voltage_scaling(opt),
+                              .calculation_initialization = get_calculation_initialization(opt)};
 }
 
 class BadCalculationRequest : public PowerGridError {

power_grid_model_c/power_grid_model_c/src/options.cpp

@@ -50,6 +50,12 @@ void PGM_set_short_circuit_voltage_scaling(PGM_Handle* handle, PGM_Options* opt,
         safe_ptr_get(opt).short_circuit_voltage_scaling = short_circuit_voltage_scaling;
     });
 }
+void PGM_set_calculation_initialization(PGM_Handle* handle, PGM_Options* opt,
+                                        PGM_Idx calculation_initialization) noexcept {
+    call_with_catch(handle, [opt, calculation_initialization] {
+        safe_ptr_get(opt).calculation_initialization = calculation_initialization;
+    });
+}
 void PGM_set_tap_changing_strategy(PGM_Handle* handle, PGM_Options* opt, PGM_Idx tap_changing_strategy) noexcept {
     call_with_catch(handle,
                     [opt, tap_changing_strategy] { safe_ptr_get(opt).tap_changing_strategy = tap_changing_strategy; });

power_grid_model_c/power_grid_model_c/src/options.hpp

@@ -22,6 +22,7 @@ struct PGM_Options {
     Idx max_iter{20};
     Idx threading{-1};
     Idx short_circuit_voltage_scaling{PGM_short_circuit_voltage_scaling_maximum};
+    Idx calculation_initialization{PGM_calculation_initialization_default};
     Idx tap_changing_strategy{PGM_tap_changing_strategy_disabled};
     Idx experimental_features{PGM_experimental_features_disabled};
 };

power_grid_model_c/power_grid_model_cpp/include/power_grid_model_cpp/options.hpp

@@ -35,6 +35,10 @@ class Options {
         handle_.call_with(PGM_set_short_circuit_voltage_scaling, get(), short_circuit_voltage_scaling);
     }
 
+    void set_calculation_initialization(Idx calculation_initialization) {
+        handle_.call_with(PGM_set_calculation_initialization, get(), calculation_initialization);
+    }
+
     void set_tap_changing_strategy(Idx tap_changing_strategy) {
         handle_.call_with(PGM_set_tap_changing_strategy, get(), tap_changing_strategy);
     }

src/power_grid_model/__init__.py

@@ -17,6 +17,7 @@
     AngleMeasurementType,
     Branch3Side,
     BranchSide,
+    CalculationInitialization,
     CalculationMethod,
     CalculationType,
     ComponentAttributeFilterOptions,
@@ -37,6 +38,7 @@
     "AttributeType",
     "Branch3Side",
     "BranchSide",
+    "CalculationInitialization",
     "CalculationMethod",
     "CalculationType",
     "ComponentAttributeFilterOptions",

src/power_grid_model/_core/enum.py

@@ -118,6 +118,32 @@ class ShortCircuitVoltageScaling(IntEnum):
     maximum = 1
 
 
+class CalculationInitialization(IntEnum):
+    """The way the calculation is initialized.
+
+    Only applicable if supported for the selected calculation type and method.
+    """
+
+    default = 0
+    """
+    The default initialization of the selected calculation method (linear for Newton-Raphson).
+    """
+    flat = 1
+    """
+    Flat start: every node at 1 p.u. with the source reference angle and its phase shift, source nodes at their
+    reference voltage and voltage regulated nodes at their reference magnitude.
+    """
+    linear = 2
+    """
+    Start from a linear voltage guess, with loads and generators as constant admittances.
+    """
+    average_source = 3
+    """
+    Every node at the average reference voltage of all sources, with its phase shift. This is the start of the
+    iterative current power flow.
+    """
+
+
 class _ExperimentalFeatures(IntEnum):
     """Experimental features"""