generalized DCOPF formulation for standard IEEE test cases

Author: abodh

DCOPF.py

@@ -0,0 +1,223 @@
+# import libraries
+import networkx
+from pyomo.environ import *
+import numpy as np
+import scipy.io as readmat
+import pdb
+import math
+import os
+
+# get the current working directory
+dir_file = os.getcwd()
+
+# testcase to use for the optimization
+test_case = '300'
+
+# degree to radian conversion factor
+deg_to_rad = math.pi / 180
+
+# extract the scenarios and their probabilities from external file
+scenario_file = dir_file + r"/scenario_files/failure_scenarios.mat"
+probability_file = dir_file + r"/scenario_files/scenario_probabilities_49.csv"
+
+# create a concrete pyomo model
+model = ConcreteModel()
+
+# load matpower test case in mat format
+matpower_mat_file = readmat.loadmat(dir_file + '/test_cases/case' + test_case + '.mat',
+                                    struct_as_record=False,
+                                    squeeze_me=True)
+# ensure that all the saved mat file are saved under workspace var name 'matpower_testcase'
+test_case = matpower_mat_file['matpower_testcase']
+
+# load the bus, generation, and line data
+model.bus = test_case.bus
+model.line = test_case.branch
+model.gen = test_case.gen
+model.gen_cost = test_case.gencost
+
+# initialize the parameters
+model.nEdges = len(model.line)  # total number of edges
+model.nNodes = len(model.bus)  # total number of nodes
+model.nGen = len(model.gen)  # total number of generators
+
+# gen dataset also has specific MVA base; use that instead when updating this version
+model.Pbase = 100  # MVA base 100 MVA to convert the system to pu system
+
+# let us assume that infinite capacity is equal to 100 GW
+Inf_transfer_Pmax = 10e6
+Pmax_line = 10e6/model.Pbase
+
+# # Pmax, Gmax, and max load for bounds set
+# if max(model.line[:, 7]) == 0:
+#     Pmax_line = Inf_transfer_Pmax / model.Pbase
+# else:
+#     Pmax_line = max(model.line[:, 7]) / model.Pbase
+
+Gmax = max(model.gen[:, 8]) / model.Pbase
+max_load = max(model.bus[:, 2]) / model.Pbase
+
+# variable ranges
+model.x_ij = range(0, model.nEdges)  # edges variable range
+model.b_i = range(0, model.nNodes)  # buses variable range
+
+###############################################################################################################
+####################################### Variables #############################################################
+###############################################################################################################
+
+# declaring pyomo variables
+
+# first we declare steady state variables i.e. power flow needs to be maintained while in the steady state
+# since these should be the same for all scenarios, they are first stage variables
+# they have ss at the end for representation
+
+# although bounds are mentioned here to maintain the standard, they will be redefined as per gen bus
+model.bus_gen_ss = Var(model.b_i, bounds=(0, Gmax), within=Reals, initialize=0)  # bus generation variable
+model.Pij_ss = Var(model.x_ij, bounds=(-Pmax_line, Pmax_line), within=Reals,
+                   initialize=0)  # active power flowing through each lines
+model.theta_ss = Var(model.b_i, bounds=(-2 * math.pi, 2 * math.pi), within=Reals, initialize=0)  # angle of each bus
+
+# model.xij = Var(model.x_ij, bounds=(0,1), within=Binary, initialize=0)  # connectivity of each line
+# model.load_shed = Var(model.b_i, bounds=(0, max_load), within=Reals)  # real active power shed at each bus
+
+
+###############################################################################################################
+####################################### Constraints ###########################################################
+###############################################################################################################
+
+# pyomo constraints
+# creates a list of constraints as placeholders
+#################### bus power balance constraints ############################
+model.power_balance = ConstraintList()
+
+# bus data col 3: active power demand, col 5: shunt conductance
+for bus_num_idx in range(model.nNodes):
+    bus_num = model.bus[bus_num_idx, 0]
+
+    # identify the list of generators connected to each bus
+    gens = np.where(model.gen[:, 0] == bus_num)[0].tolist()
+    to_bus_list = np.where(model.line[:, 1] == bus_num)[0].tolist()
+    from_bus_list = np.where(model.line[:, 0] == bus_num)[0].tolist()
+
+    model.power_balance.add(sum(model.bus_gen_ss[np.where(model.bus[:, 0] == model.gen[gen_num, 0])[0][0]]
+                                for gen_num in gens) +
+                            sum(model.Pij_ss[to_bus] for to_bus in to_bus_list) -
+                            sum(model.Pij_ss[from_bus] for from_bus in from_bus_list) ==
+                            model.bus[bus_num_idx, 2] / model.Pbase + model.bus[bus_num_idx, 4] / model.Pbase)
+
+################## generator power limit constraint ###########################
+model.gen_limit = ConstraintList()
+# generator should generate power between its min and max active power limit
+# col 9: PMAX and col 10: PMIN (Note: in Python number starts from 0)
+for gen_num in range(model.nGen):
+    model.gen_limit.add(model.bus_gen_ss[np.where(model.bus[:, 0] == model.gen[gen_num, 0])[0][0]] <=
+                        model.gen[gen_num, 8] / model.Pbase)
+    model.gen_limit.add(model.bus_gen_ss[np.where(model.bus[:, 0] == model.gen[gen_num, 0])[0][0]] >=
+                        model.gen[gen_num, 9] / model.Pbase)
+
+# make sure non-generating bus do not generate anything
+for bus_num_idx in range(model.nNodes):
+    bus_num = model.bus[bus_num_idx, 0]
+
+    if not np.any(np.equal(model.gen[:, 0], bus_num)):
+        model.gen_limit.add(model.bus_gen_ss[bus_num_idx] == 0)
+
+####################### active power flow constraint on each line ################
+model.power_flow = ConstraintList()
+'''
+Note: in Python number starts from 0
+
+linedata:
+col 4: reactance (X)
+col 9: transformer tap ratio
+col 10: transformer phase shift (in degrees)
+
+busdata:
+col 9: voltage angle (in degrees) -> this is a variable here so no need to use as parameter
+'''
+
+for line_num in range(model.nEdges):
+
+    # MATPOWER keeps 0 for transmission lines without transformer
+    # here we need to ensure tap ratio for transmission line is 1
+    if model.line[line_num, 8] == 0:
+        model.line[line_num, 8] = 1
+
+    reciprocal_term = 1 / (model.line[line_num, 3] * model.line[line_num, 8])
+    model.power_flow.add(model.Pij_ss[line_num] ==
+                         reciprocal_term * (model.theta_ss[np.where(model.bus[:, 0] == model.line[line_num, 0])[0][0]] -
+                                            model.theta_ss[np.where(model.bus[:, 0] == model.line[line_num, 1])[0][0]] -
+                                            (model.line[line_num, 9] * deg_to_rad)))
+
+################### thermal limit (MVA_limits) ############################
+# since the flow can be bi-directional, limits range from neg to positive value
+# col 6: max MVA limit of the line (0 means unlimited capacity)
+# this constraint tightens the -inf, inf bound set during variable initialization
+for line_num in range(model.nEdges):
+    if model.line[line_num, 5] == 0:
+        model.line[line_num, 5] = Pmax_line
+    model.power_flow.add(model.Pij_ss[line_num] <= model.line[line_num, 5] / model.Pbase)
+    model.power_flow.add(model.Pij_ss[line_num] >= - model.line[line_num, 5] / model.Pbase)
+
+################### angle difference between two buses on each line ################
+model.angle_limit = ConstraintList()
+# from bus and to bus reference is obtained via line
+# col 12: min angle difference (degree), col 13: max angle difference (degree)
+for angle_num in range(model.nEdges):
+    model.angle_limit.add((model.theta_ss[np.where(model.bus[:, 0] == model.line[angle_num, 0])[0][0]] -
+                           model.theta_ss[np.where(model.bus[:, 0] == model.line[angle_num, 1])[0][0]])
+                          <= model.line[angle_num, 12] * deg_to_rad)
+    model.angle_limit.add((model.theta_ss[np.where(model.bus[:, 0] == model.line[angle_num, 0])[0][0]] -
+                           model.theta_ss[np.where(model.bus[:, 0] == model.line[angle_num, 1])[0][0]])
+                          >= model.line[angle_num, 11] * deg_to_rad)
+
+# the angle can be anywhere from -2pi to 2pi hence we need to maintain 0 angle at reference (slack) bus
+
+# identifying slack bus
+slack_bus = np.where(model.bus[:, 1] == 3)[0][0]
+
+# ensure the angle at reference bus is 0
+model.angle_limit.add(model.theta_ss[slack_bus] == 0)
+
+
+# ###################### max load shedding  ############################################
+# for load_bus_num in range(model.nNodes):
+#     model.c.add(model.load_shed[load_bus_num] >= 0)
+#     model.c.add(model.load_shed[load_bus_num] <= model.bus[load_bus_num, 2])
+
+
+############################# Overall Objective ###########################
+
+def overall_objective(model):
+    expr = 0
+    expr = sum(model.gen_cost[gen_num, 4] * (model.bus_gen_ss[np.where(model.bus[:, 0] == model.gen[gen_num, 0])[0][0]]
+                                             * model.Pbase) ** 2 +
+               model.gen_cost[gen_num, 5] * (model.bus_gen_ss[np.where(model.bus[:, 0] == model.gen[gen_num, 0])[0][0]]
+                                             * model.Pbase) +
+               model.gen_cost[gen_num, 6] for gen_num in range(model.nGen))
+    # expr = sum(model.load_shed[bus_num] for bus_num in range(model.nNodes))
+    return expr
+
+
+model.min_gen_cost = Objective(rule=overall_objective, sense=minimize)
+
+####################### Solve the economic dispatch problem ################
+
+# create lp file for debugging the model
+model.write('check.lp', io_options={'symbolic_solver_labels': True})
+solver = SolverFactory('gurobi')
+results = solver.solve(model, tee=True)
+
+# Observing the results
+
+for k in range(0, model.nGen):
+    print("G[%d] = %f" % (model.gen[k, 0], value(model.bus_gen_ss[np.where(model.bus[:, 0] ==
+                                                                           model.gen[k, 0])[0][0]])))
+
+sum_P = 0;
+for k in range(0, model.nEdges):
+    print("Pij[%d] = %f" % (k + 1, value(model.Pij_ss[k])))
+    sum_P = sum_P + value(model.Pij_ss[k])
+
+for k in range(0, model.nNodes):
+    print("theta[%d] = %f" % (k + 1, value(model.theta_ss[k]) * 1 / deg_to_rad))

power_system_test_cases/case118.mat

@@ -0,0 +1,14 @@
+1	3	0.0	0.0	0	0	1	1.06	0.00	0	1	1.06	0.94
+2	2	21.7	12.7	0	0	1	1.05	-4.98	0	1	1.06	0.94
+3	2	94.2	19.0	0	0	1	1.01	-12.72	0	1	1.06	0.94
+4	1	47.8	-3.9	0	0	1	1.02	-10.33	0	1	1.06	0.94
+5	1	7.6	1.6	0	0	1	1.02	-8.78	0	1	1.06	0.94
+6	2	11.2	7.5	0	0	1	1.07	-14.22	0	1	1.06	0.94
+7	1	0.0	0.0	0	0	1	1.06	-13.37	0	1	1.06	0.94
+8	2	0.0	0.0	0	0	1	1.09	-13.36	0	1	1.06	0.94
+9	1	29.5	16.6	0	19	1	1.06	-14.94	0	1	1.06	0.94
+10	1	9.0	5.8	0	0	1	1.05	-15.10	0	1	1.06	0.94
+11	1	3.5	1.8	0	0	1	1.06	-14.79	0	1	1.06	0.94
+12	1	6.1	1.6	0	0	1	1.06	-15.07	0	1	1.06	0.94
+13	1	13.5	5.8	0	0	1	1.05	-15.16	0	1	1.06	0.94
+14	1	14.9	5.0	0	0	1	1.04	-16.04	0	1	1.06	0.94

power_system_test_cases/case14.mat

@@ -0,0 +1,5 @@
+2	0	0	3	0.0430292599000000	20	0
+2	0	0	3	0.250000000000000	20	0
+2	0	0	3	0.0100000000000000	40	0
+2	0	0	3	0.0100000000000000	40	0
+2	0	0	3	0.0100000000000000	40	0

power_system_test_cases/case300.mat

@@ -0,0 +1,5 @@
+1	232	-16.9	10	0	1	100	1	332.4	0	0	0	0	0	0	0	0	0	0	0	0
+2	40	42.4	50	-40	1	100	1	140.0	0	0	0	0	0	0	0	0	0	0	0	0
+3	0	23.4	40	0	1	100	1	100.0	0	0	0	0	0	0	0	0	0	0	0	0
+6	0	12.2	24	-6	1	100	1	100.0	0	0	0	0	0	0	0	0	0	0	0	0
+8	0	17.4	24	-6	1	100	1	100.0	0	0	0	0	0	0	0	0	0	0	0	0

power_system_test_cases/case39.mat

@@ -0,0 +1,20 @@
+1	2	0.01938	0.05917	0.0528	0	0	0	1	0	1	-360	360
+1	5	0.05403	0.22304	0.0492	0	0	0	1	0	1	-360	360
+2	3	0.04699	0.19797	0.0438	0	0	0	1	0	1	-360	360
+2	4	0.05811	0.17632	0.0340	0	0	0	1	0	1	-360	360
+2	5	0.05695	0.17388	0.0346	0	0	0	1	0	1	-360	360
+3	4	0.06701	0.17103	0.0128	0	0	0	1	0	1	-360	360
+4	5	0.01335	0.04211	0.0000	0	0	0	1	0	1	-360	360
+4	7	0.00000	0.20912	0.0000	0	0	0	0.978	0	1	-360	360
+4	9	0.00000	0.55618	0.0000	0	0	0	0.969	0	1	-360	360
+5	6	0.00000	0.25202	0.0000	0	0	0	0.932	0	1	-360	360
+6	11	0.09498	0.19890	0.0000	0	0	0	1	0	1	-360	360
+6	12	0.12291	0.25581	0.0000	0	0	0	1	0	1	-360	360
+6	13	0.06615	0.13027	0.0000	0	0	0	1	0	1	-360	360
+7	8	0.00000	0.17615	0.0000	0	0	0	1	0	1	-360	360
+7	9	0.00000	0.11001	0.0000	0	0	0	1	0	1	-360	360
+9	10	0.03181	0.08450	0.0000	0	0	0	1	0	1	-360	360
+9	14	0.12711	0.27038	0.0000	0	0	0	1	0	1	-360	360
+10	11	0.08205	0.19207	0.0000	0	0	0	1	0	1	-360	360
+12	13	0.22092	0.19988	0.0000	0	0	0	1	0	1	-360	360
+13	14	0.17093	0.34802	0.0000	0	0	0	1	0	1	-360	360