Question Details

No question body available.

Tags

python mathematical-optimization gekko pandapower

Answers (2)

Accepted Answer Available
Accepted Answer
September 19, 2025 Score: 1 Rep: 14,543 Quality: High Completeness: 80%

There are two problems:

  1. Angles: MATPOWER stores VA in degrees but GEKKO sin/cos is in radians. The theta[i] is fixed to a degree value, so the trig terms may be incorrect.
  2. Per-unit scaling: 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 idxbus, idxgen
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 getbustypes(ppc: dict) -> ArrayLike:
    return ppc["bus"][:, idxbus.BUSTYPE]
def replacepvwithpq(ppc: dict) -> dict:
    bustypes = getbustypes(ppc)
    ppc["bus"][bustypes == BusType.PV, idxbus.BUSTYPE] = 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"][:, idxbus.PD] /= base
    ppc["bus"][:, idxbus.QD] /= base
    ppc["gen"][:, idxgen.PG] /= base
    ppc["gen"][:, idxgen.QG] /= base
    ppc = replacepvwithpq(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=-2np.pi,ub=2np.pi)
    Pg = m.Array(m.Var, nb)
    Qg = m.Array(m.Var, nb)
    genbusindices = ppc["gen"][:, idxgen.GENBUS]
    bustypes = getbustypes(ppc)
    # fix variables that are actually constant
    for i in range(nb):
        if bustypes[i] == BusType.REF:
            m.fix(Vm[i], val=ppc["bus"][i, idxbus.VM])
            m.fix(theta[i], val=np.deg2rad(ppc["bus"][i, idxbus.VA]))
            continue
        if i in genbusindices and bustypes[i] == BusType.PQ:
            genidx = np.argwhere(genbusindices == i).item()
            m.fix(Pg[i], val=ppc["gen"][genidx, idxgen.PG])
            m.fix(Qg[i], val=ppc["gen"][genidx, idxgen.QG])
            continue
        m.fix(Pg[i], val=0)
        m.fix(Qg[i], val=0)
    # add parameters
    Pd = ppc["bus"][:, idxbus.PD]
    Qd = ppc["bus"][:, idxbus.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
 ---------------------------------------------------
September 20, 2025 Score: 2 Rep: 53 Quality: Low Completeness: 60%

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 idxbus, idxgen
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)
    genbusindices = ppc["gen"][:, idxgen.GENBUS].astype(int)
    bustypes = ppc["bus"][:, idxbus.BUSTYPE].astype(int)
    # fix variables that are actually constant
    for i in range(nb):
        if bustypes[i] == BusType.REF:
            m.fix(Vm[i], val=ppc["bus"][i, idxbus.VM])
            m.fix(theta[i], val=np.deg2rad(ppc["bus"][i, idxbus.VA]))
        if bustypes[i] == BusType.PQ:
            if i in genbusindices:
                m.fix(Pg[i], val=sum(ppc["gen"][genbusindices == i, idxgen.PG]) / ppc["baseMVA"])
                m.fix(Qg[i], val=sum(ppc["gen"][genbusindices == i, idxgen.QG]) / ppc["baseMVA"])
            else:
                m.fix(Pg[i], val=0)
                m.fix(Qg[i], val=0)
        if bustypes[i] == BusType.PV:
            m.fix(Vm[i], val=ppc["bus"][i, idxbus.VM])
            if i in genbusindices:
                m.fix(Pg[i], val=sum(ppc["gen"][genbusindices == i, idxgen.PG]) / ppc["baseMVA"])
            else:
                m.fix(Pg[i], val=0)
    # add parameters
    Pd = ppc["bus"][:, idxbus.PD] / ppc["baseMVA"]
    Qd = ppc["bus"][:, idxbus.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
    thetagekko = np.rad2deg(np.array([theta[i].value[0] for i in range(nb)]))
    Vmgekko = np.array([Vm[i].value[0] for i in range(nb)])
    # extract generator values (account for multiple generators )
    genbusindices = list(ppc["gen"][:, idxgen.GENBUS].astype(int))
    ng = ppc["gen"].shape[0]
    Pggekko = np.zeros((ng,))
    Qggekko = np.zeros((ng,))
    occurrenceflags = {key: False for key in np.unique(genbusindices)}
    for i in range(ng):
        busidx = genbusindices[i]
        count = genbusindices.count(busidx)
        Pggekko[i] = Pg[busidx].value[0]  ppc["baseMVA"]  int(not occurrenceflags[busidx])
        Qggekko[i] = Qg[busidx].value[0] * ppc["baseMVA"] / count
        if count > 1:
            occurrenceflags[busidx] = True
    # compare with pypower solution
    ppopt = ppoption(OUTALL=0, VERBOSE=0)
    solvedppc, success = runpf(ppc, ppopt)
    if success == 0:
        raise ValueError("PYPOWER power flow did not converge.")
    assert all(np.isclose(thetagekko, solvedppc["bus"][:, idxbus.VA])), "Voltage angles in GEKKO don't match PYPOWER"
    assert all(np.isclose(Vmgekko, solvedppc["bus"][:, idxbus.VM])), "Voltage magnitudes in GEKKO don't match PYPOWER"
    assert all(np.isclose(Pggekko, solvedppc["gen"][:, idxgen.PG])), "Generator active powers in GEKKO don't match PYPOWER"
    assert all(np.isclose(Qggekko, solvedppc["gen"][:, idxgen.QG])), "Generator reactive powers in GEKKO don't match PYPOWER"
if name == "main":
    main()