Question about how runpf.m computes the Jacobian matrix without calling newtonpf.m

[lmrklz]

The runpf.m function runs the power flow normally without any errors. However, it seems that it does not call the default newtonpf.m function, which is supposed to compute the Jacobian matrix in the Newton-Raphson method.

I even tested by deleting newtonpf.m from the directory, and runpf still ran without any issues, which confused me.

According to the documentation, the default power flow algorithm used by runpf is newtonpf, so I expected the Jacobian matrix calculation to depend directly on this function.

My questions are:

Why does runpf keep running even without the newtonpf.m file?

How is runpf computing the system Jacobian matrix in this case?

[rdzman]

The reason is the newtonpf.m is part of the legacy MATPOWER implementation that was superseded by MP-Core in MATPOWER 8. If you turn off MP-Core with have_feature('mp_core', 0), then runpf will call the legacy newtonpf to solve the power flow.

In the new approach, the non-linear constraints are added to the math model (an MP-Opt-Model object) in e.g. mp.math_model_pf_acps.add_node_balance_constraints(), using an implementation found in mp.mm_shared_pfcpf_acps.node_balance_equations(). The solve() method of the MP-Opt-Model object then evaluates and assembles the various Jacobian terms into the full matrix to be used in the Newton method.

It is much more flexible than the legacy approach, allowing for all kinds of customization, but I realize it's not as simple and transparent as the legacy version was.

[lmrklz]

O motivo é que isso newtonpf.mfaz parte da implementação legada do MATPOWER, que foi substituída pelo MP-Core no MATPOWER 8. Se você desativar o MP-Core com have_feature('mp_core', 0), então runpfchamará a implementação legada newtonpfpara resolver o fluxo de energia.

Na nova abordagem, as restrições não lineares são adicionadas ao modelo matemático (um objeto MP-Opt-Modelmp.math_model_pf_acps.add_node_balance_constraints() ), por exemplo , usando uma implementação encontrada em mp.mm_shared_pfcpf_acps.node_balance_equations(). O solve()método do objeto MP-Opt-Model então avalia e monta os vários termos do Jacobiano na matriz completa a ser usada no método de Newton.

É muito mais flexível do que a abordagem anterior, permitindo todos os tipos de personalização, mas reconheço que não é tão simples e transparente quanto a versão anterior.

Thank you for your response.
Could you please clarify why larger systems — such as those with around 70,000 buses — are not converging in version 8.1?

[rdzman]

If you are referring to case_ACTIVSg70k, I'm not sure why the MIPS solver fails, but fmincon, IPOPT, Knitro are all able to solve it fine.

[lmrklz]

Se você está se referindo a case_ACTIVSg70k, não tenho certeza do porquê o solucionador MIPS falha, mas fmincon, IPOPT e Knitro conseguem resolvê-lo sem problemas.

When I try to solve it using newtonpf it fails.

[rdzman]

Here's what I get ...

>> mpopt = mpoption('verbose', 2, 'out.all', 0);
>> runpf('case_ACTIVSg70k', mpopt)

MATPOWER Version 8.1, 12-Jul-2025
Power Flow -- AC-polar-power formulation

 it    max residual        max ∆x
----  --------------  --------------
  0      1.301e+02           -    
  1      1.384e+01       2.054e+00
  2      9.483e+00       1.865e+00
  3      8.115e-01       1.283e-01
  4      1.090e-02       2.230e-02
  5      2.641e-06       4.059e-04
  6      2.051e-11       1.079e-07
Newton's method converged in 6 iterations.
PF successful

What are you seeing?

[lmrklz]

"Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.482208e-18.

In mplinsolve (line 243)
In newtonpf (line 115)
In runpf (line 339)"

[rdzman]

That is just a warning. Does it also not converge? Could you share the full output of the two commands above?

[lmrklz]

Isso é apenas um aviso. O problema também não ocorre? Você poderia compartilhar a saída completa dos dois comandos acima?

I am trying to simulate the classical 10-bus system presented by Kundur (1994) using runpf with the newtonpf solver, but the power flow does not converge within 10 iterations, regardless of the parameter changes I apply. However, I have already simulated this same system in other software tools, where it converged without any issues.
Therefore, I would like to ask if you could provide some guidance or suggestions regarding possible reasons for this behavior. I am sharing the code used to model the system below, in case there is any mistake or inconsistency that I may have overlooked.

function mpc = case10kundur
%CASE10KUNDUR Power flow data for custom 10-bus Kundur system (Level 1)
% System fully renumbered to buses 1..10.
% Includes 5 parallel lines between buses 5 and 6.
%
% MATPOWER Case Format : Version 2

%%----- Power Flow Data -----%%
mpc.version = '2';

%% System MVA base
mpc.baseMVA = 100;

%%----- Bus Data -----%%
% bus_i type Pd Qd Gs Bs area Vm Va baseKV zone Vmax Vmin
mpc.bus = [
1 3 0 0 0 0 1 1.00 0 345 1 1.1 0.7; % slack
2 2 0 0 0 0 1 1.00 0 345 1 1.1 0.7; % PV
3 2 0 0 0 0 1 1.00 0 345 1 1.1 0.7; % PV
4 1 0 0 0 0 1 1.00 0 345 1 1.1 0.7;
5 1 0 0 0 0 1 1.00 0 345 1 1.1 0.7;
6 1 0 0 0 7.63 1 1.00 0 345 1 1.1 0.7; % shunt 763 MVAr
7 1 3271 1015 0 6.00 1 1.00 0 345 1 1.1 0.7; % load + shunt 600 MVAr
8 1 0 0 0 17.10 1 1.00 0 345 1 1.1 0.7; % shunt 1710 MVAr
9 1 0 0 0 0 1 1.00 0 345 1 1.1 0.7;
10 1 3384 971 0 0 1 1.00 0 345 1 1.1 0.7; % load
];

%%----- Generator Data -----%%
% bus Pg Qg Qmax Qmin Vg mBase status Pmax Pmin
mpc.gen = [
1 3981 0 9999 -9999 0.9800 100 1 9000 0; % G1
2 1736 0 9999 -9999 0.9646 100 1 9000 0; % G2
3 1154 0 9999 -9999 1.0400 100 1 9000 0; % G3
];

%%----- Branch Data (lines + transformers) -----%%
% fbus tbus r x b rateA rateB rateC ratio angle status angmin angmax

mpc.branch = [

%%==================== LINES ====================%%

%% Line 4–5
4 5 0.0000 0.0040 0 250 250 250 0 0 1 -360 360;

%% 5 parallel lines between 5–6
5 6 0.0015 0.0288 0.2346 250 250 250 0 0 1 -360 360;
5 6 0.0015 0.0288 0.2346 250 250 250 0 0 1 -360 360;
5 6 0.0015 0.0288 0.2346 250 250 250 0 0 1 -360 360;
5 6 0.0015 0.0288 0.2346 250 250 250 0 0 1 -360 360;
5 6 0.0015 0.0288 0.2346 250 250 250 0 0 1 -360 360;
%% Line 8–9
8 9 0.0010 0.0030 0 250 250 250 0 0 1 -360 360;

%%================ TRANSFORMERS ================%%

%% T1: 1–4
1 4 0 0.0020 0 250 250 250 0.8857 0 1 -360 360;

%% T2: 2–5
2 5 0 0.0045 0 250 250 250 0.8857 0 1 -360 360;

%% T3: 3–6
3 6 0 0.0125 0 250 250 250 0.9024 0 1 -360 360;

%% T4: 6–7
6 7 0 0.0030 0 250 250 250 1.0664 0 1 -360 360;

%% T5: 6–8
6 8 0 0.0026 0 250 250 250 1.0800 0 1 -360 360;

%% T6: 9–10 (tap variable – level 1)
9 10 0 0.0010 0 250 250 250 0.9750 0 1 -360 360;
];