Skip to content

Support voltage-regulator Q limits and PV→PQ switching #1452

Open
scud-soptim wants to merge 6 commits into
PowerGridModel:mainfrom
scud-soptim:feature/q_limit_handling_1236
Open

Support voltage-regulator Q limits and PV→PQ switching #1452
scud-soptim wants to merge 6 commits into
PowerGridModel:mainfrom
scud-soptim:feature/q_limit_handling_1236

Conversation

@scud-soptim

Copy link
Copy Markdown
Contributor

Implements #1236

Summary

This PR adds reactive-power limit handling for voltage-regulated nodes in Newton-Raphson power-flow calculations.

An active voltage_regulator now models the regulated node as a PV node by enforcing the configured voltage magnitude u_ref. If the reactive power required to maintain this voltage violates q_min or q_max, the regulated object is clamped to the violated limit and the effective node type is switched from PV to PQ.

The PR also exposes the effective node type in steady-state node output and documents the voltage-regulator behavior, including PV→PQ switching and reactive-power allocation when multiple regulators are connected to the same node.

Main changes

  • Add effective bus_type to steady-state node output.

    • 0: PQ
    • 1: PV
    • 2: Source/Slack
  • Extend voltage-regulator steady-state output with a limit indicator.

    • -1: lower reactive-power limit q_min reached
    • 0: no limit violation
    • 1: upper reactive-power limit q_max reached
  • Implement PV-node handling in Newton-Raphson power flow for active voltage regulators.

  • Implement PV→PQ switching when the required reactive power violates configured Q limits.

  • Clamp the regulated object to the violated Q limit after switching to PQ.

  • Propagate final bus type and Q-limit status to steady-state output.

  • Add and update tests for voltage-regulator Q-limit behavior, PV→PQ switching, batch updates, and output metadata.

  • Update user documentation for voltage-regulator behavior and power-flow algorithm details.

Reactive-power allocation with multiple voltage regulators

When multiple active voltage regulators are connected to the same node, they jointly regulate the same node voltage. The required reactive power is initially distributed equally over the active regulated objects. If one regulator reaches its individual q_min or q_max, it is clamped at that limit and removed from the remaining allocation. The unallocated reactive power is redistributed over the remaining regulators until all required reactive power has been allocated or all available regulators have reached a limit.

For asymmetric calculations, limit checks use the total three-phase reactive power. The phase distribution follows the available phase proportions where possible; if no usable phase proportion is available, the value is distributed equally over the three phases.

Documentation updates

The documentation now describes:

  • voltage-regulator PV-node behavior;
  • PV→PQ switching on reactive-power limit violation;
  • the meaning of node.bus_type;
  • the meaning of voltage_regulator.limit_violated;
  • reactive-power allocation for multiple voltage regulators at the same node;
  • the fact that voltage-controlled buses are supported through active voltage_regulator components in Newton-Raphson power flow.

Testing

Added and updated tests cover:

  • node output metadata for bus_type;
  • voltage-regulator output metadata for limit_violated;
  • Newton-Raphson PV-node handling;
  • PV→PQ switching on lower and upper Q-limit violations;
  • multiple voltage regulators at the same node;
  • batch calculation scenarios with voltage-regulator updates;
  • reference power-flow test data for the new PV→PQ switching case.

Notes for reviewers

  • The implementation is limited to Newton-Raphson power flow.
  • This PR does not change the public input model beyond the voltage-regulator Q-limit behavior already represented by q_min and q_max.
  • The reactive-power value is reported on the regulated load/generator output. The voltage_regulator output only reports the limit status.
  • If this PR fully resolves issue [FEATURE] Q limit handling for PV nodes #1236, the PR description can be changed from Related to #1236 to Closes #1236.

frie-soptim and others added 5 commits June 25, 2026 09:53
- add bus_type to node output
- implement feature in newton_raphson solver
- implmentet iterative equal distribution of Q at the end
  - try distributing equally
  - if individual regulator limit is violated, then save unallocated
    amount and distribute again in next iteration

Signed-off-by: Eduard Fried <eduard.fried@soptim.de>
Signed-off-by: Eduard Fried <eduard.fried@soptim.de>
Signed-off-by: Eduard Fried <eduard.fried@soptim.de>
Signed-off-by: SCUD-SOPTIM <udo.schmitz@soptim.de>
Signed-off-by: Eduard Fried <eduard.fried@soptim.de>
@figueroa1395 figueroa1395 added documentation Improvements or additions to documentation feature New feature or request improvement Improvement on internal implementation labels Jun 26, 2026
@figueroa1395

Copy link
Copy Markdown
Member

Hello @scud-soptim and @frie-soptim,

Thank you once more for your contribution!

Since this is a large and very involved PR, I'll give some initial comments:

  • We will start now the reviewing progress, but it will take us some time. Please be patient.
  • We'll address any question/comment presented during the review process itself.
  • As usual, feel free to ping us at any point in time if there is any urgency or you need help.

Additionally, we notice that this PR only includes the Q-limit handling and not yet the voltage-setpoint controller as described in #1236. We agree with this decision. Let's keep the scope of this PR strictly to the Q-limit handling and leave the voltage-setpoint controller for a follow up discussion. That said, I'll add a few additional remarks about this in the issue itself.

@figueroa1395 figueroa1395 left a comment

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.

@scud-soptim , @frie-soptim As always, great job!

Some additional comments:

  • CI is broken, but it's just minor stuff.
  • I see merge conflicts due to a recently merged documentation PR. Let us know if you want us to help resolving those.
  • There's some TODOs in power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/common_solver_functions.hpp, I agree that those should be implemented.
  • The logic looks good, no comment yet about that part.
  • The rest of the review will follow tomorrow :)

- Voltage controlled bus: a bus with known $P$ and $U$.
Note: this bus is not supported by power-grid-model yet.
- Voltage-controlled bus: a bus with known active power $P$ and voltage magnitude $U$. In power-grid-model this is
modeled by an active `voltage_regulator` on a regulated load or generator and is supported by the Newton-Raphson

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.

Maybe add a hyperlink for voltage_regulator here.

Comment on lines +34 to +38
},
{
"data_type": "IntS",
"names": ["bus_type"],
"description": "effective bus type after calculation"

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.

I'm not sure about this new output attribute. As you describe in the documentation below, doesn't the voltage_regulator.limit_violated output attribute already give enough information to determine if its bus is in PQ or PV mode? Or if multiple regulators are attached to the same node, just "loop" over such attributes?

The reason I'm hesitant is because this would now be an output attribute present for all calculation types and methods, but it's only useful for Newton Raphson Power Flow.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

We added the bus_type only here, but it would have been of interest in the previous issue (#185) as well. Because the bus type PV declares that the voltage is regulated, i.e. fixed to a set point. Yes, you can look if there are voltage regulators connected to loads/gens, which are connected to a certain node. But you would also have to check the status of these components as well. Basically, you would have to re-do what is done internally anyway. And even if you do that, with this PR the bus the can be switched back to PQ and the voltage won't be fixed anymore. So you have to also look at the limit_violated properties.

So, in order to find out, why the voltage at a node is at a particular value, you have to execute this complex logic to "re-calculate" the bus type yourself. This might be easy to look up for a small test network you created yourself, but not so easy for a larger unknown or unfamiliar network, You can implement these steps in Python, so you don't have to look up manually, but then you have to know what to look for. Right now we know this, but will we still remember in six months.

Later, maybe a new type for Soft PV nodes will be added and this would have to be signaled to the user somehow as well.

So, the question is, is a small piece of explicit information so costly, that it must be omitted from the output?
Yes, it is only needed for the newton_raphson method, but is there maybe a way to prevent or filter out the additional output for other calculation methods?

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.

You made a very good point. Let me get back to you on this, I'll think about it.

@figueroa1395 figueroa1395 Jul 23, 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.

I've brought this up with the team and even though this would indeed be very helpful to you, having an additional output value "polluting" all output when it's only relevant for NRPF, with voltage regulators, and with q-limit handling involved is something we don't want to do. That means increasing the memory footprint for everyone when it isn't needed most of the time. Additionally, even though this isn't a problem with columnar data format as one can exclude the unnecessary/unneeded output attribute types, this is always a problem for row data format. In short, we don't want bus_type as an additional output attribute type; keeping it in the core only is fine if useful.

That said, we have two proposals:

  • What we can do is create a utility function in power-grid-model-ds which determines runtime which buses are PV or PQ, which can be potentially extended in the future if more types are needed. This would provide a working solution for you right now. Would this be satisfactory?
  • We still understand that having such a feature available at the C-API or at the Python API directly via power-grid-model is desirable and we are exploring options on how to allow this via some sort of opt-in functionality (potentially available only in columnar data - but yet to be explored). Moreover, such a feature may be relevant not only for this but other "new" attributes as well. If this would be interesting to you, I'd ask you to please create a separate feature request and we can brainstorm together there on a long term better solution for this.

```

```{warning}
`bus_type` output is available only for the [Newton-Raphson power flow](./calculations.md#newton-raphson-power-flow) method; for other methods, `bus_type` is set to `0` (PQ) by default.

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.

This is what I meant in https://github.com/PowerGridModel/power-grid-model/pull/1452/changes#r3504924406. How important is this attribute for you? Rather than bus_type only being available for Newton-Raphson Power Flow, the case is that it's only relevant for that method, but it's "available", albeit useless, for the rest.

Comment on lines +1263 to +1264
Voltage regulation is supported only by the [Newton-Raphson power flow](./calculations.md#newton-raphson-power-flow)
method.

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.

Can you leave this as a warning or note please?

return node.template get_null_output<sym>();
}
auto bus_type = BusType::pq;
if (solver_output[math_id.group].bus.size() == solver_output[math_id.group].u.size()) {

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.

Isn't this always true?

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

It is only true for the newton raphson solver, because the bus array was resized only there.

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.

This might not be relevant after #1452 (comment)

}

// apply distributed Q to load_gens
for (auto const& [idx, load_gen] : enumerate(load_gens)) {

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.

I think regulating_gens already contains the relevant load_gens. That would also remove the need for the status check below.

Comment on lines +264 to +267
for (auto const& [idx, load_gen] : enumerate(load_gens)) {
if (!loadgen_to_regulator.contains(load_gen)) {
continue;
}

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.

You can also pass as argument to this function the list of relevant ones.

Comment on lines +344 to +349
if (std::abs(base_total) > numerical_tolerance) {
double const scale = q_scalar / base_total;
return RealValue<asymmetric_t>{base_distribution(0) * scale, base_distribution(1) * scale,
base_distribution(2) * scale};
}
return RealValue<asymmetric_t>{q_scalar / 3.0, q_scalar / 3.0, q_scalar / 3.0};

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.

This logic looks very similar to what's done in allocate_q_bus_limit_violated when distributing. With some minor refactoring it can probably be made into a helper function.


// auto num_regulating_gens = static_cast<double>(std::ranges::distance(regulating_gens)); // <- fails in windows
// build
double num_regulating_gens = 0.0;

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.

I just realized, why is this a double?

// When a regulator hits its limit, it is removed from the active set and the unallocated Q is
// redistributed among the remaining regulators. Unallocated Q remaining after all regulators hit
// their limits should not occur, as then the bus would have been switched to a PQ bus.
while (std::abs(total_q(q_remaining)) > numerical_tolerance && num_regulating_gens > 0.0) {

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.

Could this get into an infinite loop?

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

It shouldn't, because q_remaining is initialized with 0 in the loop and is only changed when limits were violated. At the same time, num_regulating_gens is decremented when a limit was violated, so that the value should reach 0 eventually.
However, if num_regulating_gens was calculated wrongly, then there maybe still be a small unallocated Q. but no more regulators to put in on, which could result in an infinite loop.
For such cases, I will add a check that q_remaining from iteration i should not be equal to q_remainingfrom iteration i+1.

@figueroa1395 figueroa1395 Jul 23, 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.

I agree with the check to be added.

Also, small additional question, remind me why num_regulating_gens is a double? shouldn't it be an int?

@figueroa1395 figueroa1395 left a comment

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.

Some additional comments

Comment on lines +77 to +82
// finalize
{
// Timer const sub_timer{log, LogEvent::calculate_math_result}; // TODO(mgovers): need new event ?!?
derived_solver.finalize_derived_result(input, output);
}

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.

Is the TODO a question to us?

Also, maybe doing something along the lines of

template <typename DerivedSolver>
void finalize(DerivedSolver& derived_solver, PowerFlowInput<sym> const& input, SolverOutput<sym> const& output) {
    if constexpr (requires { derived_solver.finalize_derived_result(input, output); }) {
        derived_solver.finalize_derived_result(input, output);
    }
}

makes sense so we don't have to include no-op finalize_derived_result in other solvers other than the Newton-Raphson one.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

The question was about whether a new LogEvent should be added for the "finalize result" step. But if i don't add a generic new step and instead do it like this:

        // calculate math result
        {
            Timer const sub_timer{log, LogEvent::calculate_math_result};
            if constexpr (requires { derived_solver.finalize_result(input, output); }) {
                derived_solver.finalize_result(input, output);
            }
            calculate_result(y_bus, input, output);
        }

then the timings for my function will be included in the calculate_math_result LogEvent.

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.

I understand now. Let me get back to you on this one.

@figueroa1395 figueroa1395 Jul 23, 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.

Since this is only relevant to NRPF and it's part of step of redistributing the reactive power, it's fine to keep it part of calculate_math_result as you proposed. I'd consider it part of that LogEvent. That said, if you'd prefer a separate LogEvent because you'd like to bench mark this step yourself, that's also fine. I leave it up to you.

On the other hand, I'd also like you to take a look my proposal in #1452 (comment) so we don't have to include a no-op step in the other solvers.


void parameters_changed(bool changed) { parameters_changed_ = parameters_changed_ || changed; }

void finalize_derived_result(PowerFlowInput<sym> const& /*input*/, SolverOutput<sym>& /*output*/) {

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.

// initialize with BusType::pq, set calculated value in newton_raphson solver
BusType bus_type{BusType::pq};
// 0: no violation, -1: q_min violated, 1: q_max violated
IntS q_limit_violated{0};

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.

Let's make this use a named variable or enum please.

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.

Very nice cleanup :)

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.

Thanks for the useful description of what's going on.

Additionally, would it make sense that for the node with two regulators, one violates q_max and the other violates q_min but they compensate each other? If this makes sense, can you add such a test? The same for the situation in which both violate opposite boundaries but they aren't able to compensate and the node becomes PQ instead.

@figueroa1395 figueroa1395 left a comment

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.

To be continued. I'm almost done now :)

Comment on lines +411 to +412
// TODO(frie-soptim): remove experimental feature, but add option to globally disable q-limit handling even if
// limits are set?!

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.

  • Let's remove the experimental feature flag as the last step indeed.
  • How do you propose to add this switch? Would leaving the q-limits un-initialized or set as nan enough? So when that's the case the regulator behaves as if it had an "infinite" pool of q to regulate? Or is this switch important for you?

If so, perhaps a bool to turn on-off the q-limits at the voltage regulator component-level is enough?

Comment on lines +332 to +334
// todo(frie-soptim): use "max_iter == 1" to indicate that this is the base case for a batch run,
// and thus q-limit handling does not apply here
// => maybe find a better solution

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.

Indeed using max_iter isn't ideal. This would be problematic in the hypothetical case of running with 1 iteration only and q-limit-handling enabled. It's also not very readable.

Instead, I have the following suggestion:
In

the corresponding solver is invoked and at that point Options is available, which also contains cache_run. With this in mind, you can modify the signature of run_power_flow in the math_solver at
SolverOutput<sym> run_power_flow(PowerFlowInput<sym> const& input, double err_tol, Idx max_iter, Logger& log,
so that it now additionally takes cache_run as an argument and propagate this "flag" downstream to the important code block. Don't forget to update the other solver's signatures and base class if necessary, but only NRPF should "consume" cache_run.

This way we have minimal changes and have a readable reason of why to ignore q-limit-handling (it's the first cache run). What do you think? Other suggestions are welcome.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

I tried to propagate a separate cache_run value to to NR solver, and it works. Maybe some function signatures could be left unchanged if cache_run were moved into the MainModelOptions object, but I assume this is probably not wanted.

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.

MainModelOptions is not wanted indeed, that's a proxy for the actual Options for the user to set, and mixing things up is not a good idea. If the proposed way works and you agree, let's stick to that then.

Comment on lines +384 to +385
// TODO(frie-soptim): sets max_iter to 1, which must handled in the newton-raphson solver for batch mode. Is
// there a better solution?

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.

See my other comment.

auto& bus_control = bus_control_[bus_idx];
if (!sources.empty()) {
bus_types_[bus_idx] = BusType::slack;
bus_control.type = BusType::slack;

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.

Reinforcing on earlier comments. slack type is only here for convenience but can be easily deduced as well starting from your comment. In addition, unless we would also make such changes in other solvers, I don't think we should add BusType to output types.

In addition, we had added such a check at construction time already, right? As in PGM would throw if a regulator is added where there is a source, right?

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

I extended SolverOutput with

template <symmetry_tag sym_type> struct SolverOutput {
    ....
    std::vector<ComplexValue<sym>> u;
    std::vector<ComplexValue<sym>> bus_injection;
    std::vector<BusSolverOutput> bus; // <-- NEW
    ....
};

struct BusSolverOutput {
    BusType bus_type{BusType::pq};
    LimitViolation q_limit_violated{LimitViolation::none};
};

because I needed at least the q_limit_violated value on the bus to control the Q distribution in common_solver_functions::calculate_voltage_regulator_result.

The problem is that common_solver_functions::calculate_pf_result is called from the class IterativePFSolver and not from NewtonRaphsonPFSolver, and in the the base class I don't have immediate access to the NR's internal array.

Alternatively, I would have to update the generator/load output.S values in the NT solver instead of in the calculate_pf_result function. And I would have to make sure that calculate_pf_result does not overwrite the calculated Q with the specified values. And, at least to me, it feels wrong to do that in the solver.

Therefore, I extended the SolverOuput as the option with the presumably least impact. But if you have a better proposal how to pass on NR values to the common_solver_functions::calculate_voltage_regulator_result function, then I would do that instead.

@frie-soptim frie-soptim Jul 20, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Also, we added a check that there cannot be an active Source AND and active VoltageRegulator in the same node in calculation_preparation::check_sources_and_voltage_regulators. If one of them is disabled, then it is OK.

@figueroa1395 figueroa1395 Jul 21, 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.

Let me see if I can come up with something else. I'll get back to you

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.

I'm a bit confused with the answer. I see you don't use output.bus_type on common_solver_functions::calculate_voltage_regulator_result, or did I miss it? How is that related with the question above? I think you have already access to q_limit_violated there, no?

On the other hand, maybe check by adding a small test, if not there yet, what happens when there is a regulator on a source. This should hard crash and also be added to the validator (I believe it is, but let's confirm).

std::vector<BusType> bus_types_;
struct QLimitState {
bool has_q_limits{false}; // check before using limits
IntS limit_violated{0}; // -1: q_min, +1: q_max

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.

Let's add an enum for this instead.

for (Idx const load_gen : load_gens) {
for (Idx const regulator : regulators_per_load_gen.get_element_range(load_gen)) {
if (input.load_gen_status[load_gen] != 0 && input.voltage_regulator[regulator].status != 0) {
specified_regulating_q += specified_q(load_gen, bus, input);

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.

I don't quite remember if having an asym and a sym load on the same node is allowed nor if this would handle it correctly. Could you make sure it's unit tested to make sure?

// initialize bus state, in case solver instance is reused in batching
std::ranges::fill(bus_control_, BusControlState{});

q_limits_need_initial_check_ = set_u_ref_and_bus_types_and_q_limits(input, output.u);

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.

Why the variable name change from has_usable_q_limits? Later in prepare_matrix_and_rhs it becomes confusing, does that need to have q-limits available or do they need to also be in some sort of state that requires intiial check?

Comment on lines +636 to +638
// = (Q_specified_regulating + Q_specified_nonregulating - del_x_pq_[bus].q())
// - Q_specified_nonregulating
// = Q_specified_regulating - del_x_pq_[bus].q()

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.

These couple of lines made me a bit confused. With
Since Q_specified_all = Q_specified_regulating + Q_specified_nonregulating
and the above equation we can reach the desired result. Could you remove them?

const bool buses_switched = !q_limits_need_initial_check_ && enforce_q_limits(input);
if (buses_switched) {
// rebuild jacobian and del_pq after PV -> PQ switching
build_jacobian_and_rhs(y_bus, input, u);

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.

I believe only add_loads and add_sources would be "affected" or important in this additional re-build. Would it be possible to do that instead of a full rebuild? So a partial one instead?

Or does that make the logic too complex and not worth the effort given the small gains? Please think about this.

On the other hand, I see that a first full build is needed because you use del_x_pq_[bus].q() to determine q_limits. Would it be possible to consider using a "base reference" from the previous run or an approximated one from the initial "guess" start? This is thinking on how not to do two builds.

Comment on lines +583 to +588
return 1;
}
if (!is_nan(q_min) && q_total < q_min - numerical_tolerance) {
return -1;
}
return 0;

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.

Let's enum this.

@scud-soptim

Copy link
Copy Markdown
Contributor Author

Hi, we'll answer your questions when Frie returns from vacation.

@figueroa1395 figueroa1395 left a comment

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.

Double checking, but is there any additional validation needed on the python side? Say anything regarding the q-limits for instance? What would happen if q_min is set to 1 and q_max to -1 by mistake? I imagine this should be catched by both the validator and a hard error thrown?

Can you add some unit testing handling sym and asym regulated loads on a node? one violating limits, both violating limits, and both not violating limits?

I think I should be done with the overall review. If I have anything else to comment on, I'll mention it once Frie is back so this doesn't grow even more.

// The specified Q values from regulating generators were already added to del_x_pq_[bus].q()
// in add_loads(). To extract the ACTUAL Q that regulators are providing (including their
// contribution to voltage control), we must "subtract back out" this mismatch term.
// This works because: actual = specified - mismatch.

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.

Does this mean Q_calc_network = Q_specified_all - del_x_pq_[bus].q() as indicated above, or is this line referring to something different and I am miss-understanding things?

Comment on lines +674 to +675
clamped_regulators_per_load_gen_[load_gen] = RealValue<asymmetric_t>{
q_limit_scalar / 3.0, q_limit_scalar / 3.0, q_limit_scalar / 3.0};

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.

Is this fallback physically meaningful enough? Could it happen that the input q is zero but the regulators fall out of the limit anyway? Could you test this?

Comment on lines +681 to +683
// TODO(frie-soptim): don't forget to reset regulators if the bus is switched back to PV in the future
// clamped_regulators_per_load_gen_[load_gen] = RealValue<sym>{nan};

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.

This would only be relevant in case we enable switching back from PQ to PV after the first PV to PQ transition has already happened because the limits were violated, right?

So for now, this isn't relevant, correct?

Comment on lines +692 to +697
static bool is_value_nan(RealValue<sym> const& value) {
if constexpr (is_symmetric_v<sym>) {
return is_nan(value);
} else {
// all three phases are set together, need to check only one phase for NaN
return is_nan(value(0));

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.

Why is this needed? Wouldn't the overload of is_nan present in common.hpp already work for this?

// load_gen hit a Q-limit, use the clamped value
del_x_pq_[bus_number].p() += real(input.s_injection[load_number]);
del_x_pq_[bus_number].q() += clamped_q;
// TODO(frie-soptim): prevent regulation of const_i/const_y load_gens in validation

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.

What do you mean?

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

The thought was that Q for const_i and const_y is dependent on voltage. But the limit is just a fix MVar value. So when this limit is initially violated, then the bus is switched to PQ and the voltage now becomes flexible. Will the adjusted Q (with the new calculated Voltage) still violated the original limit? Maybe yes, maybe no.

Also, it is unclear if the limit itself is voltage dependent or not. It is probably just a fixed MVar value independent of voltage. But, if I would do bus type switching "manually", by disabling the regulator and setting the q_specified value on the load_gen to the limit value. The bus will then be a PQ bus and the limit/q_specified will now be treated as voltage dependent, which contradicts the previous assumption that the limit is not voltage dependent.

So, the behavior of voltage regulation with limit consideration for const_y/const_i load_gens is quite unclear, and it probably does not reflect reality anyway. If a generator regulates voltage by freely moving Q in some limit range, then how can this Q be voltage dependent via an additional equation?

Therefore the idea was to exclude such cases via validation or exceptions. Now, it is currently possible to connect voltage regulators with const_i/const_y load_gens, so maybe this can be excused by saying the user knows what he is doing if he calculates such a network, instead of simply disallowing it. The validation could then prohibit the usage of limits for such cases.

@figueroa1395 figueroa1395 Jul 23, 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.

We want to look at the math with a bit more detail on how should we handle this case. Naively, we see no issue and we think const_y and const_i load_gens should just allot the excess Q if the q-limit is violated in PV mode and increase the V respectively, but again, this is naive.

In short, we will look at this math and come back to you later. We'll ping you so you don't miss it once we have an answer..

C.C. @nitbharambe

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.

We just went over this again and we have a conclusion. Let's throw an exception when other load_gen type other than const_p is used with voltage regulators. That is, we only support it for PQ-PV nodes.

This can be thrown at construction time since load_gen.type isn't updatable. There should also be a check included to the validator and this unit tested.

This way we don't have to worry about this logic at all. Do you agree?

Comment on lines +104 to +105
CHECK(imag(output.load_gen[0].s) > -1.0);
CHECK(imag(output.load_gen[0].s) < 1.0);

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.

I imaging checking for the exact value is not practical? But double checking why a range instead? Also, why not check for q directly?

TEST_CASE_TEMPLATE_INVOKE(test_math_solver_pf_id, NewtonRaphsonPFSolver<symmetric_t>);
TEST_CASE_TEMPLATE_INVOKE(test_math_solver_pf_id, NewtonRaphsonPFSolver<asymmetric_t>);

TEST_CASE("Newton-Raphson PV reactive-power limits switch one-way to PQ") {

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.

I think the name is a bit missleading since in the first subcase it never switches to PQ since no limits are violated.

}
}

TEST_CASE("Newton-Raphson PV reactive-power limit switches are deterministic") {

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.

Why the "deterministic" in the name? Curious.

CHECK(output.bus[3].bus_type == BusType::pv); // Two regulators
}

SUBCASE("Single regulator controls voltage") {

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.

Can you explicitly check the other remaining two don't?

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.

Incredibly extended testing! Kudos for the awesome work!

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

....we answer next week to your comments....

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.

@scud-soptim No problem and no rush, I just wanted to finish the overall review before moving on. We won't be making any more comments/reviews until you have had the time to read and address the current ones.

Signed-off-by: Martijn Govers <martijn.govers@alliander.com>
@figueroa1395

figueroa1395 commented Jul 14, 2026

Copy link
Copy Markdown
Member

Hello @scud-soptim and @frie-soptim,

Once small remark for when you return. We have enabled mandatory cryptographic signing for the PGM-org. Could you please amend your commits and either ssh or gpg sign them? This is additional to DCO signing which was already in place.

This was announced during the last community meeting, where @scud-soptim was present.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

documentation Improvements or additions to documentation feature New feature or request improvement Improvement on internal implementation

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants