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
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
---------------------------------------------------