Menu
I want to solve the AC power flow equations (wiki) using GEKKO in Python. This is a nonlinear system of equations describing conservation of energy on a power grid.
I'm struggling to get GEKKO to solve them. I get @error: Solution Not Found. I have checked that this particular instance does converge using a library specialized for power flow equations. The reason I'm trying to solve the system using GEKKO is to do a proof of concept and later expand the problem into a mixed integer nonlinear program.
What am I doing wrong? Is there a way to get more insight into why the solver isn't converging?
Here's the code (the main function is where the problem is defined):
import copy
from enum import IntEnum
import hydra
import numpy as np
from gekko import GEKKO
from numpy._typing import ArrayLike
from omegaconf import DictConfig
from pypower import idx_bus, idx_gen
import pandapower.networks as pn
from pandapower import runpp
from pypower.makeYbus import makeYbus
from src import CONFIGS_PATH
class BusType(IntEnum):
PQ = 1
PV = 2
REF = 3
ISOLATED = 4
def get_bus_types(ppc: dict) -> ArrayLike:
return ppc["bus"][:, idx_bus.BUS_TYPE]
def replace_pv_with_pq(ppc: dict) -> dict:
bus_types = get_bus_types(ppc)
ppc["bus"][bus_types == BusType.PV, idx_bus.BUS_TYPE] = BusType.PQ
return ppc
def main():
# load power grid data
net = pn.case14()
runpp(net)
ppc = copy.deepcopy(net._ppc)
ppc = replace_pv_with_pq(ppc)
# define variables
nb = ppc["bus"].shape[0]
m = GEKKO(remote=False)
Vm = m.Array(m.Var, nb, lb=0, value=1)
theta = m.Array(m.Var, nb)
Pg = m.Array(m.Var, nb)
Qg = m.Array(m.Var, nb)
gen_bus_indices = ppc["gen"][:, idx_gen.GEN_BUS]
bus_types = get_bus_types(ppc)
# fix variables that are actually constant
for i in range(nb):
if bus_types[i] == BusType.REF:
m.fix(Vm[i], val=ppc["bus"][i, idx_bus.VM])
m.fix(theta[i], val=ppc["bus"][i, idx_bus.VA])
continue
if i in gen_bus_indices and bus_types[i] == BusType.PQ:
gen_idx = np.argwhere(gen_bus_indices == i).item()
m.fix(Pg[i], val=ppc["gen"][gen_idx, idx_gen.PG])
m.fix(Qg[i], val=ppc["gen"][gen_idx, idx_gen.QG])
continue
m.fix(Pg[i], val=0)
m.fix(Qg[i], val=0)
# add parameters
Pd = ppc["bus"][:, idx_bus.PD]
Qd = ppc["bus"][:, idx_bus.QD]
P = Pg - Pd
Q = Qg - Qd
Ybus, _, _ = makeYbus(ppc["baseMVA"], ppc["bus"], ppc["branch"])
Gbus = Ybus.real.toarray()
Bbus = Ybus.imag.toarray()
# active power conservation
m.Equations(
[
0 == -P[i] + sum([Vm[i] * Vm[k] * (Gbus[i, k] * m.cos(theta[i] - theta[k]) + Bbus[i, k] * m.sin(theta[i] - theta[k])) for k in range(nb)]) for i in range(nb)
]
)
# reactive power conservation
m.Equations(
[
0 == -Q[i] + sum([Vm[i] * Vm[k] * (Gbus[i, k] * m.sin(theta[i] - theta[k]) - Bbus[i, k] * m.cos(theta[i] - theta[k])) for k in range(nb)]) for i in range(nb)
]
)
m.solve(disp=True)
if __name__ == "__main__":
main()
and I get the following output:
----------------------------------------------------------------
APMonitor, Version 1.0.3
APMonitor Optimization Suite
----------------------------------------------------------------
--------- APM Model Size ------------
Each time step contains
Objects : 0
Constants : 0
Variables : 56
Intermediates: 0
Connections : 56
Equations : 28
Residuals : 28
Number of state variables: 28
Number of total equations: - 28
Number of slack variables: - 0
---------------------------------------
Degrees of freedom : 0
solver 3 not supported
using default solver: APOPT
----------------------------------------------
Steady State Optimization with APOPT Solver
----------------------------------------------
Iter Objective Convergence
0 8.58297E-18 1.17750E+01
1 1.57354E-16 1.17750E+01
2 4.53362E-17 1.17236E+01
...(this goes on for awhile)
Maximum iterations
---------------------------------------------------
Solver : IPOPT (v3.12)
Solution time : 0.336599999980535 sec
Objective : 0.000000000000000E+000
Unsuccessful with error code 0
---------------------------------------------------
Creating file: infeasibilities.txt
@error: Solution Not Found
209
asked May 25 '26 15:05
Thank you to the answer above! In case anyone needs it, here's a full working example of solving the PF equations in GEKKO that includes all bus types (REF/Slack, PQ, and PV):
import copy
from enum import IntEnum
import numpy as np
import pandapower.networks as pn
from gekko import GEKKO
from pandapower import runpp
from pypower import idx_bus, idx_gen
from pypower.makeYbus import makeYbus
from pypower.ppoption import ppoption
from pypower.runpf import runpf
class BusType(IntEnum):
PQ = 1
PV = 2
REF = 3
ISOLATED = 4
def main():
# load power grid data
net = pn.case14()
# net = pn.case30()
# net = pn.case57()
# net = pn.case118()
# net = pn.case300() # @error: Max Equation Length
runpp(net)
if not net.converged:
raise ValueError("Pandapower power flow did not converge.")
ppc = copy.deepcopy(net._ppc)
# define variables
m = GEKKO(remote=False)
# define variables
nb = ppc["bus"].shape[0]
Vm = m.Array(m.Var, nb, lb=0, value=1)
theta = m.Array(m.Var, nb, lb=-np.pi, ub=np.pi)
Pg = m.Array(m.Var, nb)
Qg = m.Array(m.Var, nb)
gen_bus_indices = ppc["gen"][:, idx_gen.GEN_BUS].astype(int)
bus_types = ppc["bus"][:, idx_bus.BUS_TYPE].astype(int)
# fix variables that are actually constant
for i in range(nb):
if bus_types[i] == BusType.REF:
m.fix(Vm[i], val=ppc["bus"][i, idx_bus.VM])
m.fix(theta[i], val=np.deg2rad(ppc["bus"][i, idx_bus.VA]))
if bus_types[i] == BusType.PQ:
if i in gen_bus_indices:
m.fix(Pg[i], val=sum(ppc["gen"][gen_bus_indices == i, idx_gen.PG]) / ppc["baseMVA"])
m.fix(Qg[i], val=sum(ppc["gen"][gen_bus_indices == i, idx_gen.QG]) / ppc["baseMVA"])
else:
m.fix(Pg[i], val=0)
m.fix(Qg[i], val=0)
if bus_types[i] == BusType.PV:
m.fix(Vm[i], val=ppc["bus"][i, idx_bus.VM])
if i in gen_bus_indices:
m.fix(Pg[i], val=sum(ppc["gen"][gen_bus_indices == i, idx_gen.PG]) / ppc["baseMVA"])
else:
m.fix(Pg[i], val=0)
# add parameters
Pd = ppc["bus"][:, idx_bus.PD] / ppc["baseMVA"]
Qd = ppc["bus"][:, idx_bus.QD] / ppc["baseMVA"]
P = Pg - Pd
Q = Qg - Qd
Ybus, _, _ = makeYbus(ppc["baseMVA"], ppc["bus"], ppc["branch"])
Gbus = Ybus.real.toarray()
Bbus = Ybus.imag.toarray()
# active power conservation
m.Equations(
[
0 == -P[i] + sum([Vm[i] * Vm[k] * (Gbus[i, k] * m.cos(theta[i] - theta[k]) + Bbus[i, k] * m.sin(theta[i] - theta[k])) for k in range(nb)]) for i in range(nb)
]
)
# reactive power conservation
m.Equations(
[
0 == -Q[i] + sum([Vm[i] * Vm[k] * (Gbus[i, k] * m.sin(theta[i] - theta[k]) - Bbus[i, k] * m.cos(theta[i] - theta[k])) for k in range(nb)]) for i in range(nb)
]
)
m.options.SOLVER = 1
m.options.RTOL = 1e-8
m.solve(disp=True)
# extract values
theta_gekko = np.rad2deg(np.array([theta[i].value[0] for i in range(nb)]))
Vm_gekko = np.array([Vm[i].value[0] for i in range(nb)])
# extract generator values (account for multiple generators )
gen_bus_indices = list(ppc["gen"][:, idx_gen.GEN_BUS].astype(int))
ng = ppc["gen"].shape[0]
Pg_gekko = np.zeros((ng,))
Qg_gekko = np.zeros((ng,))
occurrence_flags = {key: False for key in np.unique(gen_bus_indices)}
for i in range(ng):
bus_idx = gen_bus_indices[i]
count = gen_bus_indices.count(bus_idx)
Pg_gekko[i] = Pg[bus_idx].value[0] * ppc["baseMVA"] * int(not occurrence_flags[bus_idx])
Qg_gekko[i] = Qg[bus_idx].value[0] * ppc["baseMVA"] / count
if count > 1:
occurrence_flags[bus_idx] = True
# compare with pypower solution
ppopt = ppoption(OUT_ALL=0, VERBOSE=0)
solved_ppc, success = runpf(ppc, ppopt)
if success == 0:
raise ValueError("PYPOWER power flow did not converge.")
assert all(np.isclose(theta_gekko, solved_ppc["bus"][:, idx_bus.VA])), "Voltage angles in GEKKO don't match PYPOWER"
assert all(np.isclose(Vm_gekko, solved_ppc["bus"][:, idx_bus.VM])), "Voltage magnitudes in GEKKO don't match PYPOWER"
assert all(np.isclose(Pg_gekko, solved_ppc["gen"][:, idx_gen.PG])), "Generator active powers in GEKKO don't match PYPOWER"
assert all(np.isclose(Qg_gekko, solved_ppc["gen"][:, idx_gen.QG])), "Generator reactive powers in GEKKO don't match PYPOWER"
if __name__ == "__main__":
main()
There are two problems:
theta[i] is fixed to a degree value, so the trig terms may be incorrect.Pg,Qg,Pd,Qd in the MATPOWER/Pandapower ppc are in MW/MVAr, while the Y-bus may be per-unit. The power-balance equations must use per-unit injections (divide by baseMVA).Below are changes that solve successfully. I've removed a few packages from the import list that weren't required.
import copy
from enum import IntEnum
import numpy as np
from gekko import GEKKO
from numpy._typing import ArrayLike
from pypower import idx_bus, idx_gen
import pandapower.networks as pn
from pandapower import runpp
from pypower.makeYbus import makeYbus
class BusType(IntEnum):
PQ = 1
PV = 2
REF = 3
ISOLATED = 4
def get_bus_types(ppc: dict) -> ArrayLike:
return ppc["bus"][:, idx_bus.BUS_TYPE]
def replace_pv_with_pq(ppc: dict) -> dict:
bus_types = get_bus_types(ppc)
ppc["bus"][bus_types == BusType.PV, idx_bus.BUS_TYPE] = BusType.PQ
return ppc
def main():
# load power grid data
net = pn.case14()
runpp(net)
ppc = copy.deepcopy(net._ppc)
# per-unit scaling for powers
base = ppc["baseMVA"]
ppc["bus"][:, idx_bus.PD] /= base
ppc["bus"][:, idx_bus.QD] /= base
ppc["gen"][:, idx_gen.PG] /= base
ppc["gen"][:, idx_gen.QG] /= base
ppc = replace_pv_with_pq(ppc)
# define variables
nb = ppc["bus"].shape[0]
m = GEKKO(remote=False)
Vm = m.Array(m.Var, nb, lb=0, value=1)
theta = m.Array(m.Var, nb, lb=-2*np.pi,ub=2*np.pi)
Pg = m.Array(m.Var, nb)
Qg = m.Array(m.Var, nb)
gen_bus_indices = ppc["gen"][:, idx_gen.GEN_BUS]
bus_types = get_bus_types(ppc)
# fix variables that are actually constant
for i in range(nb):
if bus_types[i] == BusType.REF:
m.fix(Vm[i], val=ppc["bus"][i, idx_bus.VM])
m.fix(theta[i], val=np.deg2rad(ppc["bus"][i, idx_bus.VA]))
continue
if i in gen_bus_indices and bus_types[i] == BusType.PQ:
gen_idx = np.argwhere(gen_bus_indices == i).item()
m.fix(Pg[i], val=ppc["gen"][gen_idx, idx_gen.PG])
m.fix(Qg[i], val=ppc["gen"][gen_idx, idx_gen.QG])
continue
m.fix(Pg[i], val=0)
m.fix(Qg[i], val=0)
# add parameters
Pd = ppc["bus"][:, idx_bus.PD]
Qd = ppc["bus"][:, idx_bus.QD]
P = Pg - Pd
Q = Qg - Qd
Ybus, _, _ = makeYbus(ppc["baseMVA"], ppc["bus"], ppc["branch"])
Gbus = Ybus.real.toarray()
Bbus = Ybus.imag.toarray()
# active power conservation
m.Equations(
[
0 == -P[i] + sum([Vm[i] * Vm[k] * (Gbus[i, k] * m.cos(theta[i] - theta[k]) + Bbus[i, k] * m.sin(theta[i] - theta[k])) for k in range(nb)]) for i in range(nb)
]
)
# reactive power conservation
m.Equations(
[
0 == -Q[i] + sum([Vm[i] * Vm[k] * (Gbus[i, k] * m.sin(theta[i] - theta[k]) - Bbus[i, k] * m.cos(theta[i] - theta[k])) for k in range(nb)]) for i in range(nb)
]
)
m.options.SOLVER = 1
m.solve(disp=True)
if __name__ == "__main__":
main()
Here is the output from the script:
----------------------------------------------------------------
APMonitor, Version 1.0.3
APMonitor Optimization Suite
----------------------------------------------------------------
--------- APM Model Size ------------
Each time step contains
Objects : 0
Constants : 0
Variables : 56
Intermediates: 0
Connections : 56
Equations : 28
Residuals : 28
Number of state variables: 28
Number of total equations: - 28
Number of slack variables: - 0
---------------------------------------
Degrees of freedom : 0
----------------------------------------------
Steady State Optimization with APOPT Solver
----------------------------------------------
Iter Objective Convergence
0 1.46423E-20 1.17750E-01
1 4.96153E-22 1.10222E-02
2 9.25933E-25 3.73276E-04
3 9.25933E-25 3.73276E-04
Successful solution
---------------------------------------------------
Solver : APOPT (v1.0)
Solution time : 0.02789999999999998 sec
Objective : 0.
Successful solution
---------------------------------------------------
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With
