<jupyter_start><jupyter_text>PF benchmark  Manual solution<jupyter_code>from pypower.api import ext2int, bustypes, makeBdc, makeSbus, dcpf, int2ext, case9

from numpy import r_, c_, ix_, zeros, pi, ones, exp, argmax, union1d
from numpy import flatnonzero as find

from pypower.idx_bus import PD, QD, VM, VA, GS, BUS_TYPE, PV, PQ, REF
from pypower.idx_brch import PF, PT, QF, QT
from pypower.idx_gen import PG, QG, VG, QMAX, QMIN, GEN_BUS, GEN_STATUS

# Read PyPower case
ppc = case9()

## add zero columns to branch for flows if needed 
if ppc["branch"].shape[1] < QT: 
    ppc["branch"] = c_[ppc["branch"], 
                    zeros((ppc["branch"].shape[0], 
                            QT - ppc["branch"].shape[1] + 1))] 

## convert to internal indexing
ppc = ext2int(ppc)
baseMVA, bus, gen, branch = \
    ppc["baseMVA"], ppc["bus"], ppc["gen"], ppc["branch"]

## get bus index lists of each type of bus
ref, pv, pq = bustypes(bus, gen)

## generator info
on = find(gen[:, GEN_STATUS] > 0)      ## which generators are on?
gbus = gen[on, GEN_BUS].astype(int)    ## what buses are they at?

##-----  run the power flow  -----
## initial state
Va0 = bus[:, VA] * (pi / 180)

## build B matrices and phase shift injections
B, Bf, _, _ = makeBdc(baseMVA, bus, branch)

## compute complex bus power injections [generation - load]
Pbus = makeSbus(baseMVA, bus, gen).real

## "run" the power flow
pvpq = np.matrix(r_[pv, pq]) 

# initialize result vector 
# Va = np.copy(Va0)
Va = np.zeros_like(Va0)

B_dense = B.todense()
B_bar = B_dense[pvpq.T, pvpq]
P_bar = Pbus[pvpq]

## update angles for non-reference buses
Va[pvpq] = np.dot(np.linalg.inv(B_bar), P_bar.T).reshape(1,8)
 
# update data matrices with solution
branch[:, [QF, QT]] = zeros((branch.shape[0], 2))
branch[:, PF] = (Bf * Va) * baseMVA
branch[:, PT] = -branch[:, PF]
bus[:, VM] = ones(bus.shape[0])
bus[:, VA] = Va * (180 / pi)
## update Pg for slack generator (1st gen at ref bus)
## (note: other gens at ref bus are accounted for in Pbus)
##      Pg = Pinj + Pload
##      newPg = oldPg + newPinj - oldPinj
refgen = find(gbus == ref)
gen[on[refgen[0]], PG] = gen[on[refgen[0]], PG] + (B[ref, :] * Va - Pbus[ref])* baseMVA


##-----  output results  -----
## convert back to original bus numbering & print results
ppc["bus"], ppc["gen"], ppc["branch"] = bus, gen, branch
results = int2ext(ppc)
r_bus = results["bus"]
r_branch = results["branch"]
r_gen = results["gen"]<jupyter_output><empty_output><jupyter_text>PyPower solution<jupyter_code>ppc = case9()
pp_dcpf, success = pp.rundcpf(ppc)
pp_bus = pp_dcpf["bus"]
pp_branch = pp_dcpf["branch"]
pp_gen = pp_dcpf["gen"]<jupyter_output>PYPOWER Version 5.1.4, 27-June-2018 -- DC Power Flow

Converged in 0.00 seconds
================================================================================
|     System Summary                                                           |
================================================================================

How many?                How much?              P (MW)            Q (MVAr)
---------------------    -------------------  -------------  -----------------
Buses              9     Total Gen Capacity     820.0           0.0 to 0.0
Generators         3     On-line Capacity       820.0           0.0 to 0.0
Committed Gens     3     Generation (actual)    314.0               0.0
Loads              3     Load                   315.0               0.0
  Fixed            3       Fixed                315.0               0.0
  Dispatchable     0       Dispatchable           0.0 of 0.0        0.0
Shunts             0     Shunt (inj)              0.0               0.0
Branches    [...]<jupyter_text>Comparison<jupyter_code>print("Bus voltage magitude:")
print("Manual: ", r_bus[:,VM])
print("Pypower: ", pp_bus[:,VM])

print("Bus voltage angle:")
print("Manual: ", r_bus[:,VA])
print("Pypower: ", pp_bus[:,VA])

print("Real power generation:")
print("Manual: ", r_gen[:,PG])
print("PyPower: ", pp_gen[:,PG])

print("Real power flow on branch:")
print("Manual: ", r_branch[:,PF])
print("PyPower: ", pp_branch[:,PF])<jupyter_output>Bus voltage magitude:
Manual:  [1. 1. 1. 1. 1. 1. 1. 1. 1.]
Pypower:  [1. 1. 1. 1. 1. 1. 1. 1. 1.]
Bus voltage angle:
Manual:  [ 0.          9.79601886  5.06056005 -2.21115872 -3.73809125  2.20665727
  0.82244106  3.95901132 -4.06340049]
Pypower:  [ 0.          9.79601886  5.06056005 -2.21115872 -3.73809125  2.20665727
  0.82244106  3.95901132 -4.06340049]
Real power generation:
Manual:  [ 66 163  85]
PyPower:  [ 66 163  85]
Real power flow on branch:
Manual:  [  67.          28.9673913  -61.0326087   85.          23.9673913
  -76.0326087 -163.          86.9673913  -38.0326087]
PyPower:  [  67.          28.9673913  -61.0326087   85.          23.9673913
  -76.0326087 -163.          86.9673913  -38.0326087]<jupyter_text>OPF benchmark  Manual solution<jupyter_code>import numpy as np
from numpy import flatnonzero as find
from pypower.api import case9, ext2int, bustypes, makeBdc, rundcopf, ppoption, case39
from pypower.idx_bus import BUS_TYPE, REF, VA, PD, LAM_P, LAM_Q, MU_VMAX, MU_VMIN
from pypower.idx_gen import PG, MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PMAX, PMIN, GEN_BUS
from pypower.idx_brch import PF, PT, QF, QT, RATE_A, MU_SF, MU_ST
from pypower.idx_cost import COST

# Run PyPower case to get initial condition
ppc = case9()
ppopt = ppoption(VERBOSE=0)
pp_dcpf = rundcopf(ppc, ppopt)
pp_bus = pp_dcpf["bus"]
pp_branch = pp_dcpf["branch"]
pp_gen = pp_dcpf["gen"]
## convert to internal indexing
ppc = ext2int(ppc)
baseMVA, bus, gen, branch, gencost = \
    ppc["baseMVA"], ppc["bus"], ppc["gen"], ppc["branch"], ppc["gencost"]

## get bus index lists of each type of bus
ref, pv, pq = bustypes(bus, gen)
pvpq = np.matrix(np.r_[pv, pq])

## generator info
gbus = gen[:, GEN_BUS].astype(int)    ## what buses are they at?
refgen = find(gbus == ref)

## build B matrices and phase shift injections
B, Bf, _, _ = makeBdc(baseMVA, bus, branch)
B = B.todense()
Bf = Bf.todense()

# Problem dimensions
NG = gen.shape[0] # Number of generators
NB = bus.shape[0] # Number of buses
NBr = branch.shape[0] # Number of lines
NL = 3 # Number of loads

CG = np.zeros((NB,NG))
CG[gbus,range(NG)] = 1

# Generator capacity limit p.u.
Pmax = gen[:, PMAX]/baseMVA
Pmin = gen[:, PMIN]/baseMVA

# Line flow limit p.u.
Lmax = branch[:, RATE_A]/baseMVA
Lmin = - branch[:, RATE_A]/baseMVA

# Quadratic cost coefficients
# Convert to p.u.
Cg2 = gencost[:, COST]*baseMVA**2
Cg1 = gencost[:, COST+1]*baseMVA
Cg0 = gencost[:, COST+2]

# Generator and demand set points
# Power in p.u., Va in rad
Pg0 = pp_gen[:, PG]/baseMVA
Pd0 = pp_bus[:, PD]/baseMVA
Pinj0 = (np.matmul(CG, Pg0) - Pd0)
Va0 = pp_bus[:, VA]/180*np.pi
from pyomo.environ import *

model = ConcreteModel(name="OPF")

# Define sets
model.buses = Set(initialize=range(NB))
model.lines = Set(initialize=range(NBr))
model.gens = Set(initialize=range(NG))
model.loads = Set(initialize=range(NL))

# Define parameters
def Pd_init(model, i):
    return Pd0[i]
model.Pd = Param(model.buses, initialize=Pd_init)

def B_init(model, i, j):
    return B[i,j]
model.B = Param(model.buses, model.buses, initialize=B_init)

def Bf_init(model, i, j):
    return Bf[i,j]
model.Bf = Param(model.lines, model.lines, initialize=Bf_init)

def Pmax_init(model, i):
    return Pmax[i]
model.Pmax = Param(model.gens, initialize=Pmax_init)

def Pmin_init(model, i):
    return Pmin[i]
model.Pmin = Param(model.gens, initialize=Pmin_init)

def Lmax_init(model, i):
    return Lmax[i]
model.Lmax = Param(model.lines, initialize=Lmax_init)

def Lmin_init(model, i):
    return Lmin[i]
model.Lmin = Param(model.lines, initialize=Lmin_init)

def Cg0_init(model, i):
    return Cg0[i]
model.Cg0 = Param(model.gens, initialize=Cg0_init)

def Cg1_init(model, i):
    return Cg1[i]
model.Cg1 = Param(model.gens, initialize=Cg1_init)

def Cg2_init(model, i):
    return Cg2[i]
model.Cg2 = Param(model.gens, initialize=Cg2_init)

def CG_init(model, i, j):
    return CG[i,j]
model.CG = Param(model.buses, model.gens, initialize=CG_init)
# Define variables
def Pg_init(model, i):
    return Pg0[i]
model.Pg = Var(range(NG), within=PositiveReals, initialize=Pg_init, bounds=(0,3))

def Pinj_init(model, i):
    return Pinj0[i]
model.Pinj = Var(range(NB), within=Reals, initialize=Pinj_init, bounds=(-2,2))

def Va_init(model, i):
    return Va0[i]
model.Va = Var(range(NB), within=Reals, initialize=Va_init, bounds=(-1,1))

model.dual = Suffix(direction=Suffix.IMPORT)
# Define objectives
def Gen_cost(model):
    cost0 = sum(model.Cg0[i] for i in model.gens)
    cost1 = sum(model.Cg1[i] * model.Pg[i] for i in model.gens)
    cost2 = sum(model.Cg2[i] * model.Pg[i]**2 for i in model.gens)
    return (cost0 + cost1 + cost2)
model.obj = Objective(rule=Gen_cost, sense=minimize)
# Define Constraints
def Pmax_rule(model, i):
    return model.Pg[i] <= model.Pmax[i]
model.c_Pmax = Constraint(model.gens, rule=Pmax_rule)

def Pmin_rule(model, i):
    return model.Pg[i] >= model.Pmin[i]
model.c_Pmin = Constraint(model.gens, rule=Pmin_rule)

def Lmax_rule(model, i):
    return sum(model.Bf[i,j]*model.Va[j] for j in model.buses) <= model.Lmax[i]
model.c_Lmax = Constraint(model.lines, rule=Lmax_rule)

def Lmin_rule(model, i):
    return sum(model.Bf[i,j]*model.Va[j] for j in model.buses) >= model.Lmin[i]
model.c_Lmin = Constraint(model.lines, rule=Lmin_rule)

def Pinj_rule(model, i):
    return model.Pinj[i] == (sum(model.CG[i,j]*model.Pg[j] for j in model.gens) - model.Pd[i])
model.c_Pinj = Constraint(model.buses, rule=Pinj_rule)

def PF_rule(model, i):
    return model.Pinj[i] == sum(model.B[i,j]*model.Va[j] for j in model.buses)
model.c_PF = Constraint(model.buses, rule=PF_rule)
# Solve the problem
solver = SolverFactory("mindtpy")
solver.solve(model)
model.display()
# Print dual variables (shadow prices)
model.dual.pprint()
# Calculate average price
price = 0
for i in range(NB):
    price = price + model.dual[model.c_PF[i]]
price = price/NB/baseMVA
print("Average price = %.3f" %price)
# Generation and cost
r_pg = np.zeros((NG, 1))
r_cost = np.zeros((NG, 1))
r_mc = np.zeros((NG, 1))

for i in range(NG):
    r_pg[i] = model.Pg[i]()
    r_cost[i] = Cg0[i] + Cg1[i]*r_pg[i] + Cg2[i]*r_pg[i]**2
    r_mc[i] = Cg1[i]/baseMVA + Cg2[i]*r_pg[i]/baseMVA

print("Pg = \n", r_pg)
print("Cost = \n", r_cost)
print("MC = \n", r_mc)
print(model.obj())
print(np.sum(r_cost))<jupyter_output>5216.026607747273
5216.026607747273<jupyter_text>PyPower solution<jupyter_code>import pypower.api as pp
ppc = case9()
pp_dcopf = pp.rundcopf(ppc)
pp_bus = pp_dcopf["bus"]
pp_branch = pp_dcopf["branch"]
pp_gen = pp_dcopf["gen"]
pp_pg = pp_gen[:, PG]/baseMVA
pp_va = pp_bus[:, VA]/180*np.pi
pp_pd = pp_bus[:, PD]/baseMVA
pp_pinj = (np.matmul(CG, pp_pg) - pp_pd)
# Power injection at buses
pp_pbus = np.matmul(B, pp_va)*baseMVA
print(pp_pbus)
pp_pinj = (np.matmul(CG, pp_pg) - pp_pd)*baseMVA
print(pp_pinj)
pp_cost = np.sum(Cg0) + np.sum(Cg1 * pp_pg) + np.sum(Cg2 * pp_pg**2)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17063203777697059, 'fraction_numerical': 0.045677442789683984, 'mean_word_length': 3.5773067331670823, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '45821', 'license': 'MIT', 'n_tokens_mistral': 4680, 'n_tokens_neox': 4383, 'n_words': 1060, 'path': 'pypower_test_case9.ipynb', 'repo_name': 'boyuan276/stochasticOPF', 'size': '11012'}