Skip to content

Fix short-circuit source impedance scaling #1547

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
shin4141 wants to merge 1 commit into
base: main
Choose a base branch
from
+3,390 −3,330
Open

Fix short-circuit source impedance scaling #1547

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 20 additions & 3 deletions ...id_model_c/power_grid_model/include/power_grid_model/math_solver/short_circuit_solver.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,13 +92,31 @@ template <symmetry_tag sym> class ShortCircuitSolver {
auto& diagonal_element = mat_data_[diagonal_position];
auto& u_bus = output.u_bus[bus_number];

detail::add_sources<sym>(sources, bus_number, y_bus, input.source, diagonal_element, u_bus);
add_sources(sources, y_bus, input, diagonal_element, u_bus);

add_faults(faults, bus_number, y_bus, input, diagonal_element, u_bus, infinite_admittance_fault_counter,
fault_type, phase_1, phase_2);
}
}

// Source admittance is calculated from sk for a specific u_ref in PGM's power-flow solvers. In short-circuit
// calculations, sk is defined at rated voltage, so the source-equivalent impedance scales with c and the
// corresponding admittance scales with 1 / c = 1 / cabs(input.source[source_number]).
static ComplexTensor<sym> source_admittance(YBus<sym> const& y_bus, ShortCircuitInput const& input,
Idx source_number) {
double const scaling = 1.0 / cabs(input.source[source_number]);
return scaling * y_bus.math_model_param().source_param[source_number].template y_ref<sym>();
}
Comment on lines +105 to +109

@nitbharambe nitbharambe Aug 28, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nitpick: Could you please move the above comment of calculation_input_preparation regarding the reason to this place?

I would like a explaining docstring here along the lines:
"source admittance is calculated on sk for a specific u_ref for powerflow solvers within PGM.
However for generic short circuit calculation, sk is defined at rated voltage, so its source-equivalent impedance scales with c and the corresponding admittance scales with 1 / c = 1 / cabs(input.source[source_number])"

@shin4141 shin4141 Sep 10, 2026

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done. I moved the explanation to ShortCircuitSolver::source_admittance() and documented
why the short-circuit source admittance is y_ref / c, with
c = cabs(input.source[source_number]).


static void add_sources(IdxRange const& sources, YBus<sym> const& y_bus, ShortCircuitInput const& input,
ComplexTensor<sym>& diagonal_element, ComplexValue<sym>& u_bus) {
for (Idx const source_number : sources) {
ComplexTensor<sym> const y_source = source_admittance(y_bus, input, source_number);
diagonal_element += y_source;
u_bus += dot(y_source, ComplexValue<sym>{input.source[source_number]});
}
}

void add_faults(IdxRange const& faults, Idx bus_number, YBus<sym> const& y_bus, ShortCircuitInput const& input,
ComplexTensor<sym>& diagonal_element, ComplexValue<sym>& u_bus,
IdxVector& infinite_admittance_fault_counter, FaultType const& fault_type, IntS phase_1,
Expand Down Expand Up @@ -307,8 +325,7 @@ template <symmetry_tag sym> class ShortCircuitSolver {
ComplexValue<sym> i_source_bus{}; // total source current in to the bus
ComplexValue<sym> i_source_inject{}; // total raw source current as a Norton equivalent
for (Idx const source_number : sources) {
ComplexTensor<sym> const y_source =
y_bus.math_model_param().source_param[source_number].template y_ref<sym>();
ComplexTensor<sym> const y_source = source_admittance(y_bus, input, source_number);
ComplexValue<sym> const i_source_inject_single =
dot(y_source, ComplexValue<sym>{input.source[source_number]});
output.source[source_number].i = i_source_inject_single - dot(y_source, output.u_bus[bus_number]);
Expand Down
93 changes: 68 additions & 25 deletions tests/cpp_unit_tests/math_solver/test_math_solver_sc.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -158,6 +158,7 @@ TEST_CASE("Short circuit solver") {
double const vref = 1.1;
DoubleComplex const yref{10.0 - 50.0i};
DoubleComplex const zref{1.0 / yref};
DoubleComplex const zref_scaled{vref * zref};
// line
DoubleComplex const y0{1.0 - 2.0i};
DoubleComplex const y0_0{0.5 + 0.5i};
Expand Down Expand Up @@ -185,7 +186,7 @@ TEST_CASE("Short circuit solver") {
YBus<asymmetric_t> const y_bus_asym{topo_sc, param_sc_asym};
ShortCircuitSolver<asymmetric_t> solver{y_bus_asym, topo_sc};
auto sc_input = create_sc_test_input(three_phase, FaultPhase::abc, y_fault, vref, fault_buses);
auto sc_output_ref = create_sc_test_output<asymmetric_t>(three_phase, z_fault, z0, z0_0, vref, zref);
auto sc_output_ref = create_sc_test_output<asymmetric_t>(three_phase, z_fault, z0, z0_0, vref, zref_scaled);
auto output = solver.run_short_circuit(y_bus_asym, sc_input);
assert_sc_output<asymmetric_t>(output, sc_output_ref);

Expand All @@ -198,7 +199,8 @@ TEST_CASE("Short circuit solver") {
YBus<asymmetric_t> const y_bus_asym{topo_sc, param_sc_asym};
ShortCircuitSolver<asymmetric_t> solver{y_bus_asym, topo_sc};
auto sc_input = create_sc_test_input(three_phase, FaultPhase::abc, y_fault_solid, vref, fault_buses);
auto sc_output_ref = create_sc_test_output<asymmetric_t>(three_phase, z_fault_solid, z0, z0_0, vref, zref);
auto sc_output_ref =
create_sc_test_output<asymmetric_t>(three_phase, z_fault_solid, z0, z0_0, vref, zref_scaled);
auto output = solver.run_short_circuit(y_bus_asym, sc_input);
assert_sc_output<asymmetric_t>(output, sc_output_ref);
}
Expand All @@ -207,7 +209,7 @@ TEST_CASE("Short circuit solver") {
YBus<symmetric_t> const y_bus_sym{topo_sc, param_sc_sym};
ShortCircuitSolver<symmetric_t> solver{y_bus_sym, topo_sc};
auto sc_input = create_sc_test_input(three_phase, FaultPhase::abc, y_fault, vref, fault_buses);
auto sc_output_ref = create_sc_test_output<symmetric_t>(three_phase, z_fault, z0, z0_0, vref, zref);
auto sc_output_ref = create_sc_test_output<symmetric_t>(three_phase, z_fault, z0, z0_0, vref, zref_scaled);
auto output = solver.run_short_circuit(y_bus_sym, sc_input);
assert_sc_output<symmetric_t>(output, sc_output_ref);

Expand All @@ -220,7 +222,8 @@ TEST_CASE("Short circuit solver") {
YBus<symmetric_t> const y_bus_sym{topo_sc, param_sc_sym};
ShortCircuitSolver<symmetric_t> solver{y_bus_sym, topo_sc};
auto sc_input = create_sc_test_input(three_phase, FaultPhase::abc, y_fault_solid, vref, fault_buses);
auto sc_output_ref = create_sc_test_output<symmetric_t>(three_phase, z_fault_solid, z0, z0_0, vref, zref);
auto sc_output_ref =
create_sc_test_output<symmetric_t>(three_phase, z_fault_solid, z0, z0_0, vref, zref_scaled);
auto output = solver.run_short_circuit(y_bus_sym, sc_input);
assert_sc_output<symmetric_t>(output, sc_output_ref);
}
Expand All @@ -229,7 +232,8 @@ TEST_CASE("Short circuit solver") {
YBus<asymmetric_t> const y_bus_asym{topo_sc, param_sc_asym};
ShortCircuitSolver<asymmetric_t> solver{y_bus_asym, topo_sc};
auto sc_input = create_sc_test_input(single_phase_to_ground, FaultPhase::a, y_fault, vref, fault_buses);
auto sc_output_ref = create_sc_test_output<asymmetric_t>(single_phase_to_ground, z_fault, z0, z0_0, vref, zref);
auto sc_output_ref =
create_sc_test_output<asymmetric_t>(single_phase_to_ground, z_fault, z0, z0_0, vref, zref_scaled);
auto output = solver.run_short_circuit(y_bus_asym, sc_input);
assert_sc_output<asymmetric_t>(output, sc_output_ref);

Expand All @@ -243,7 +247,7 @@ TEST_CASE("Short circuit solver") {
ShortCircuitSolver<asymmetric_t> solver{y_bus_asym, topo_sc};
auto sc_input = create_sc_test_input(single_phase_to_ground, FaultPhase::a, y_fault_solid, vref, fault_buses);
auto sc_output_ref =
create_sc_test_output<asymmetric_t>(single_phase_to_ground, z_fault_solid, z0, z0_0, vref, zref);
create_sc_test_output<asymmetric_t>(single_phase_to_ground, z_fault_solid, z0, z0_0, vref, zref_scaled);
auto output = solver.run_short_circuit(y_bus_asym, sc_input);
assert_sc_output<asymmetric_t>(output, sc_output_ref);
}
Expand All @@ -252,7 +256,7 @@ TEST_CASE("Short circuit solver") {
YBus<asymmetric_t> const y_bus_asym{topo_sc, param_sc_asym};
ShortCircuitSolver<asymmetric_t> solver{y_bus_asym, topo_sc};
auto sc_input = create_sc_test_input(two_phase, FaultPhase::bc, y_fault, vref, fault_buses);
auto sc_output_ref = create_sc_test_output<asymmetric_t>(two_phase, z_fault, z0, z0_0, vref, zref);
auto sc_output_ref = create_sc_test_output<asymmetric_t>(two_phase, z_fault, z0, z0_0, vref, zref_scaled);
auto output = solver.run_short_circuit(y_bus_asym, sc_input);
assert_sc_output<asymmetric_t>(output, sc_output_ref);

Expand All @@ -264,7 +268,7 @@ TEST_CASE("Short circuit solver") {
YBus<asymmetric_t> const y_bus_asym{topo_sc, param_sc_asym};
ShortCircuitSolver<asymmetric_t> solver{y_bus_asym, topo_sc};
auto sc_input = create_sc_test_input(two_phase, FaultPhase::bc, y_fault_solid, vref, fault_buses);
auto sc_output_ref = create_sc_test_output<asymmetric_t>(two_phase, z_fault_solid, z0, z0_0, vref, zref);
auto sc_output_ref = create_sc_test_output<asymmetric_t>(two_phase, z_fault_solid, z0, z0_0, vref, zref_scaled);
auto output = solver.run_short_circuit(y_bus_asym, sc_input);
assert_sc_output<asymmetric_t>(output, sc_output_ref);
}
Expand All @@ -273,7 +277,8 @@ TEST_CASE("Short circuit solver") {
YBus<asymmetric_t> const y_bus_asym{topo_sc, param_sc_asym};
ShortCircuitSolver<asymmetric_t> solver{y_bus_asym, topo_sc};
auto sc_input = create_sc_test_input(two_phase_to_ground, FaultPhase::bc, y_fault, vref, fault_buses);
auto sc_output_ref = create_sc_test_output<asymmetric_t>(two_phase_to_ground, z_fault, z0, z0_0, vref, zref);
auto sc_output_ref =
create_sc_test_output<asymmetric_t>(two_phase_to_ground, z_fault, z0, z0_0, vref, zref_scaled);
auto output = solver.run_short_circuit(y_bus_asym, sc_input);
assert_sc_output<asymmetric_t>(output, sc_output_ref);

Expand All @@ -287,7 +292,7 @@ TEST_CASE("Short circuit solver") {
ShortCircuitSolver<asymmetric_t> solver{y_bus_asym, topo_sc};
auto sc_input = create_sc_test_input(two_phase_to_ground, FaultPhase::bc, y_fault_solid, vref, fault_buses);
auto sc_output_ref =
create_sc_test_output<asymmetric_t>(two_phase_to_ground, z_fault_solid, z0, z0_0, vref, zref);
create_sc_test_output<asymmetric_t>(two_phase_to_ground, z_fault_solid, z0, z0_0, vref, zref_scaled);
auto output = solver.run_short_circuit(y_bus_asym, sc_input);
assert_sc_output<asymmetric_t>(output, sc_output_ref);
}
Expand Down Expand Up @@ -330,25 +335,25 @@ TEST_CASE("Short circuit solver") {
ShortCircuitSolver<asymmetric_t> solver{y_bus_asym, topo_comp};
ShortCircuitSolver<symmetric_t> sym_solver{y_bus_sym, topo_comp};

DoubleComplex const if_comp = vref / (zref + z_fault);
DoubleComplex const uf_comp = vref - if_comp * zref;
DoubleComplex const if_comp_solid = vref / (zref + z_fault_solid);
DoubleComplex const uf_comp_solid = vref - if_comp_solid * zref;
DoubleComplex const if_comp = vref / (zref_scaled + z_fault);
DoubleComplex const uf_comp = vref - if_comp * zref_scaled;
DoubleComplex const if_comp_solid = vref / (zref_scaled + z_fault_solid);
DoubleComplex const uf_comp_solid = vref - if_comp_solid * zref_scaled;

DoubleComplex const if_b_comp = (vref * (a * a - a)) / (2.0 * zref + z_fault);
DoubleComplex const uf_b_comp = vref * a * a - if_b_comp * zref;
DoubleComplex const uf_c_comp = vref * a + if_b_comp * zref;
DoubleComplex const if_b_comp = (vref * (a * a - a)) / (2.0 * zref_scaled + z_fault);
DoubleComplex const uf_b_comp = vref * a * a - if_b_comp * zref_scaled;
DoubleComplex const uf_c_comp = vref * a + if_b_comp * zref_scaled;

DoubleComplex const if_b_comp_solid = (vref * (a * a - a)) / (2.0 * zref + z_fault_solid);
DoubleComplex const uf_b_comp_solid = vref * a * a - if_b_comp_solid * zref;
DoubleComplex const uf_c_comp_solid = vref * a + if_b_comp_solid * zref;
DoubleComplex const if_b_comp_solid = (vref * (a * a - a)) / (2.0 * zref_scaled + z_fault_solid);
DoubleComplex const uf_b_comp_solid = vref * a * a - if_b_comp_solid * zref_scaled;
DoubleComplex const uf_c_comp_solid = vref * a + if_b_comp_solid * zref_scaled;

DoubleComplex const uf_b_2phg = (vref * (a * a + a)) * z_fault / (zref + 2.0 * z_fault);
DoubleComplex const if_b_2phg = (vref * a * a - uf_b_2phg) / zref;
DoubleComplex const if_c_2phg = (vref * a - uf_b_2phg) / zref;
DoubleComplex const uf_b_2phg = (vref * (a * a + a)) * z_fault / (zref_scaled + 2.0 * z_fault);
DoubleComplex const if_b_2phg = (vref * a * a - uf_b_2phg) / zref_scaled;
DoubleComplex const if_c_2phg = (vref * a - uf_b_2phg) / zref_scaled;
DoubleComplex const uf_b_2phg_solid = 0.0 + 0.0i;
DoubleComplex const if_b_2phg_solid = vref * a * a / zref;
DoubleComplex const if_c_2phg_solid = vref * a / zref;
DoubleComplex const if_b_2phg_solid = vref * a * a / zref_scaled;
DoubleComplex const if_c_2phg_solid = vref * a / zref_scaled;

SUBCASE("Source on 3ph sym fault") {
ShortCircuitSolverOutput<symmetric_t> sc_output_ref;
Expand Down Expand Up @@ -462,4 +467,42 @@ TEST_CASE("Short circuit solver") {
}
}

TEST_CASE("Short circuit source admittance follows the voltage factor") {
MathModelTopology topology;
topology.slack_bus = 0;
topology.phase_shift = {0.0};
topology.sources_per_bus = {from_sparse, {0, 1}};
topology.shunts_per_bus = {from_sparse, {0, 0}};
topology.load_gens_per_bus = {from_sparse, {0, 0}};

DoubleComplex const source_admittance{10.0 - 50.0i};
MathModelParam<symmetric_t> parameters;
parameters.source_param = {SourceCalcParam{.y1 = source_admittance, .y0 = source_admittance}};

YBus<symmetric_t> const y_bus{topology, std::move(parameters)};
ShortCircuitSolver<symmetric_t> solver{y_bus, topology};

DoubleComplex const solid_fault_admittance{std::numeric_limits<double>::infinity(),
std::numeric_limits<double>::infinity()};
auto const run_with_voltage_factor = [&](double voltage_factor) {
ShortCircuitInput input;
input.source = {{voltage_factor, 0.0}};
input.fault_buses = {from_sparse, {0, 1}};
input.faults = {
{.y_fault = solid_fault_admittance, .fault_type = FaultType::three_phase, .fault_phase = FaultPhase::abc}};
return solver.run_short_circuit(y_bus, input);
};

auto const low_voltage_minimum_output = run_with_voltage_factor(0.95);
auto const high_voltage_minimum_output = run_with_voltage_factor(1.0);
auto const maximum_output = run_with_voltage_factor(1.1);

CHECK(cabs(low_voltage_minimum_output.fault[0].i_fault - source_admittance) < numerical_tolerance);
CHECK(cabs(high_voltage_minimum_output.fault[0].i_fault - source_admittance) < numerical_tolerance);
CHECK(cabs(maximum_output.fault[0].i_fault - source_admittance) < numerical_tolerance);
CHECK(cabs(maximum_output.source[0].i - source_admittance) < numerical_tolerance);
CHECK(cabs(maximum_output.fault[0].i_fault - low_voltage_minimum_output.fault[0].i_fault) < numerical_tolerance);
CHECK(cabs(maximum_output.fault[0].i_fault - high_voltage_minimum_output.fault[0].i_fault) < numerical_tolerance);
}

} // namespace power_grid_model::math_solver
4 changes: 2 additions & 2 deletions tests/data/short_circuit/branch_source_single_phase_to_ground_c_maximum/sc_output.json
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,11 +29,11 @@
"id": 5,
"energized": 1,
"i_f": [
63508.52961085883410,
57735.026918962576,
0.0,
0.0
]
}
]
}
}
}
6 changes: 3 additions & 3 deletions tests/data/short_circuit/branch_source_single_phase_to_ground_c_maximum/sc_output_batch.json
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@
"id": 5,
"energized": 1,
"i_f": [
63508.52961085883410,
57735.026918962576,
0.0,
0.0
]
Expand All @@ -43,7 +43,7 @@
"id": 1,
"energized": 1,
"u_pu": [
1.083914433075477,
1.0817685322751267,
1.1,
1.1
]
Expand All @@ -63,7 +63,7 @@
"id": 5,
"energized": 1,
"i_f": [
6257.9828971464731,
6245.59353309911,
0.0,
0.0
]
Expand Down
Loading