[Python] How to solve power flow equations using GEKKO in Python?

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 == BusType.REF:
m.fix(Vm, val=ppc["bus"][i, idx_bus.VM])
m.fix(theta, val=ppc["bus"][i, idx_bus.VA])
continue

if i in gen_bus_indices and bus_types == BusType.PQ:
gen_idx = np.argwhere(gen_bus_indices == i).item()
m.fix(Pg, val=ppc["gen"][gen_idx, idx_gen.PG])
m.fix(Qg, val=ppc["gen"][gen_idx, idx_gen.QG])
continue

m.fix(Pg, val=0)
m.fix(Qg, 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 + sum([Vm * Vm[k] * (Gbus[i, k] * m.cos(theta - theta[k]) + Bbus[i, k] * m.sin(theta - theta[k])) for k in range(nb)]) for i in range(nb)
]
)

# reactive power conservation
m.Equations(
[
0 == -Q + sum([Vm * Vm[k] * (Gbus[i, k] * m.sin(theta - theta[k]) - Bbus[i, k] * m.cos(theta - 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

Continue reading...