# Copyright (c) 1996-2016, Power Systems Engineering Research Center (PSERC) by Ray Zimmerman, PSERC Cornell
# Modifications Copyright (c) 2026, Liangyu Zhang
# SPDX-License-Identifier: BSD-3-Clause
import numpy as np
from ..corex import MatpowerConfig
from .make_vcorr import make_vcorr
from .make_zpv import make_zpv
from .mpoption import get_zip_weights
[docs]
def calc_v_y_sum(Vslack, nb, nl, f, Zb, Ybf, Ybt, Yd, Sd, pv, Pg, Vg, mpopt: MatpowerConfig):
"""Solve radial power flow by the admittance summation method.
Mirrors MATPOWER's ``calc_v_y_sum`` radial helper. It performs iterative
backward/forward sweeps using equivalent admittances and current
injections, optionally applies PV-bus voltage correction, and returns
voltages, flows, slack power, iteration count, and convergence status.
Parameters
----------
Vslack : complex
Slack-bus voltage.
nb : int
Number of buses.
nl : int
Number of branches.
f : array_like
One-based parent-bus indices for each branch.
Zb : array_like
Branch series impedances.
Ybf : array_like
Branch shunt admittances at the from end.
Ybt : array_like
Branch shunt admittances at the to end.
Yd : array_like
Bus shunt admittances.
Sd : array_like
Complex bus power demands.
pv : array_like
One-based PV bus indices.
Pg : array_like
Real generation at PV buses.
Vg : array_like
Target voltage magnitudes at PV buses.
mpopt : dict
MATPOWER options struct.
Returns
-------
tuple
``(V, Qpv, Sf, St, Sslack, iterations, success)``.
"""
tol = mpopt.pf.tol
iter_max = mpopt.pf.radial.max_it
vcorr = mpopt.pf.radial.vcorr == 1.0
nb = int(np.asarray(nb).reshape(-1)[0])
nl = int(np.asarray(nl).reshape(-1)[0])
f = np.asarray(f, dtype=int).reshape(-1)
Zb = np.asarray(Zb, dtype=complex).reshape(-1)
Ybf = np.asarray(Ybf, dtype=complex).reshape(-1)
Ybt = np.asarray(Ybt, dtype=complex).reshape(-1)
Yd = np.asarray(Yd, dtype=complex).reshape(-1)
Sd = np.asarray(Sd, dtype=complex).reshape(-1).copy()
pv = np.asarray(pv, dtype=int).reshape(-1)
Pg = np.asarray(Pg, dtype=float).reshape(-1)
Vg = np.asarray(Vg, dtype=float).reshape(-1)
f = f - 1
pv = pv - 1
Sd[pv] = Sd[pv] - Pg
V = Vslack * np.ones(nb, dtype=complex)
Vold = V.copy()
iter_count = 0
success = 0.0
pw, qw = get_zip_weights(mpopt)
Sdz = np.real(Sd) * pw[2] + 1j * np.imag(Sd) * qw[2]
Sdi = np.real(Sd) * pw[1] + 1j * np.imag(Sd) * qw[1]
Sdp = np.real(Sd) * pw[0] + 1j * np.imag(Sd) * qw[0]
f = np.r_[0, f]
Zb = np.r_[0.0 + 0.0j, Zb]
nl = nl + 1
if pv.size:
Zpv = make_zpv(pv + 1, nb, nl, f + 1, Zb, Yd)
Bpv = np.linalg.inv(np.imag(Zpv))
else:
Bpv = np.empty((0, 0))
npv = pv.size
Qpv = np.zeros(npv)
Ye = np.conj(Sdz) + Yd
D = np.zeros(nl, dtype=complex)
for k in range(nl - 1, 0, -1):
D[k] = 1.0 / (1.0 + Zb[k] * Ye[k])
i = f[k]
Ye[i] = Ye[i] + D[k] * Ye[k]
Je = np.zeros(nl, dtype=complex)
while success == 0.0 and iter_count < iter_max:
iter_count += 1
Vm = np.abs(V)
S = Sdp + Sdi * Vm
Je = np.conj(S / V)
for k in range(nl - 1, 0, -1):
i = f[k]
Je[i] = Je[i] + D[k] * Je[k]
for k in range(1, nl):
i = f[k]
V[k] = D[k] * (V[i] - Zb[k] * Je[k])
DU = np.abs(V - Vold)
DU[np.isnan(DU)] = np.inf
if np.max(DU) > tol:
Vold = V.copy()
if pv.size:
DE = (Vg / np.abs(V[pv]) - 1.0) * np.real(V[pv])
DD = Bpv @ DE
if vcorr:
DC = DD * np.imag(V[pv]) / np.real(V[pv])
V = V + np.asarray(make_vcorr(DC + 1j * DD, pv + 1, nb, nl, f + 1, Zb)).reshape(-1)
DQ = DD * np.abs(V[pv]) ** 2 / np.real(V[pv])
Qpv = Qpv + DQ
Sdp[pv] = Sdp[pv] - 1j * DQ
else:
success = 1.0
Sslack = V[0] * np.conj(Je[0]) + np.conj(Ye[0]) * np.abs(V[0]) ** 2
f = f[1:]
Zb = Zb[1:]
I = (V[f] - V[1:]) / Zb
Sf = V[f] * np.conj(I) + np.conj(Ybf) * np.abs(V[f]) ** 2
St = V[1:] * np.conj(I) - np.conj(Ybt) * np.abs(V[1:]) ** 2
return (
V.reshape(-1, 1),
Qpv.reshape(-1, 1),
Sf.reshape(-1, 1),
St.reshape(-1, 1),
Sslack,
float(iter_count),
success,
)