option to keep variables for out of service elements

egret/model_library/transmission/tx_utils.py

@@ -17,11 +17,40 @@
 from math import radians
 import logging
 import copy
+from collections import OrderedDict
+from egret.data.data_utils import map_items, zip_items
 
 
 logger = logging.getLogger(__name__)
 
 
+def get_unique_bus_pairs(md):
+    branch_attrs = md.attributes(element_type='branch')
+    bus_pairs = zip_items(branch_attrs['from_bus'], branch_attrs['to_bus'])
+    unique_bus_pairs = list(OrderedDict((val, None) for idx, val in bus_pairs.items()))
+    return unique_bus_pairs
+
+
+def _get_out_of_service_gens(md):
+    out_of_service_gens = list()
+    for gen_name, gen_dict in md.elements(element_type='generator'):
+        if not gen_dict['in_service']:
+            out_of_service_gens.append(gen_name)
+            gen_dict['in_service'] = True
+
+    return out_of_service_gens
+
+
+def _get_out_of_service_branches(md):
+    out_of_service_branches = list()
+    for branch_name, branch_dict in md.elements(element_type='branch'):
+        if not branch_dict['in_service']:
+            out_of_service_branches.append(branch_name)
+            branch_dict['in_service'] = True
+
+    return out_of_service_branches
+
+
 def dicts_of_vr_vj(buses):
     """
     Create dictionaries of vr and vj values from the bus vm and va values

egret/models/ac_relaxations.py

@@ -217,8 +217,10 @@ def use_linear_relaxation(self, val):
 def create_soc_relaxation(model_data,
                           use_linear_relaxation=True,
                           include_feasibility_slack=False,
-                          use_fbbt=True):
-    model, md = _create_base_power_ac_model(model_data, include_feasibility_slack=include_feasibility_slack)
+                          use_fbbt=True,
+                          keep_vars_for_out_of_service_elements=False):
+    model, md = _create_base_power_ac_model(model_data, include_feasibility_slack=include_feasibility_slack,
+                                            keep_vars_for_out_of_service_elements=keep_vars_for_out_of_service_elements)
     if use_linear_relaxation:
         _relaxation_helper(model=model,
                            md=md,
@@ -238,8 +240,10 @@ def create_atan_relaxation(model_data,
                            use_linear_relaxation=True,
                            include_feasibility_slack=False,
                            use_soc_edge_cuts=False,
-                           use_fbbt=True):
-    model, md = create_atan_acopf_model(model_data=model_data, include_feasibility_slack=include_feasibility_slack)
+                           use_fbbt=True,
+                           keep_vars_for_out_of_service_elements=False):
+    model, md = create_atan_acopf_model(model_data=model_data, include_feasibility_slack=include_feasibility_slack,
+                                        keep_vars_for_out_of_service_elements=keep_vars_for_out_of_service_elements)
     del model.ineq_soc
     del model._con_ineq_soc
     if use_soc_edge_cuts:
@@ -268,8 +272,10 @@ def create_polar_acopf_relaxation(model_data,
                                   include_soc=True,
                                   use_linear_relaxation=True,
                                   include_feasibility_slack=False,
-                                  use_fbbt=True):
-    model, md = create_psv_acopf_model(model_data, include_feasibility_slack=include_feasibility_slack)
+                                  use_fbbt=True,
+                                  keep_vars_for_out_of_service_elements=False):
+    model, md = create_psv_acopf_model(model_data, include_feasibility_slack=include_feasibility_slack,
+                                       keep_vars_for_out_of_service_elements=keep_vars_for_out_of_service_elements)
     _relaxation_helper(model=model,
                        md=md,
                        include_soc=include_soc,
@@ -282,8 +288,10 @@ def create_rectangular_acopf_relaxation(model_data,
                                         include_soc=True,
                                         use_linear_relaxation=True,
                                         include_feasibility_slack=False,
-                                        use_fbbt=True):
-    model, md = create_rsv_acopf_model(model_data, include_feasibility_slack=include_feasibility_slack)
+                                        use_fbbt=True,
+                                        keep_vars_for_out_of_service_elements=False):
+    model, md = create_rsv_acopf_model(model_data, include_feasibility_slack=include_feasibility_slack,
+                                       keep_vars_for_out_of_service_elements=keep_vars_for_out_of_service_elements)
     _relaxation_helper(model=model,
                        md=md,
                        include_soc=include_soc,

egret/models/acopf.py

@@ -68,7 +68,14 @@ def _validate_and_extract_slack_penalties(model_data):
     return model_data.data['system']['load_mismatch_cost'], model_data.data['system']['q_load_mismatch_cost']
 
 
-def _create_base_power_ac_model(model_data, include_feasibility_slack=False, pw_cost_model='delta'):
+def _create_base_power_ac_model(model_data, include_feasibility_slack=False, pw_cost_model='delta', keep_vars_for_out_of_service_elements=False):
+    if keep_vars_for_out_of_service_elements:
+        out_of_service_gens = tx_utils._get_out_of_service_gens(model_data)
+        out_of_service_branches = tx_utils._get_out_of_service_branches(model_data)
+    else:
+        out_of_service_gens = list()
+        out_of_service_branches = list()
+
     md = model_data.clone_in_service()
     tx_utils.scale_ModelData_to_pu(md, inplace=True)
 
@@ -86,8 +93,7 @@ def _create_base_power_ac_model(model_data, include_feasibility_slack=False, pw_
         tx_utils.inlet_outlet_branches_by_bus(branches, buses)
     gens_by_bus = tx_utils.gens_by_bus(buses, gens)
 
-    bus_pairs = zip_items(branch_attrs['from_bus'], branch_attrs['to_bus'])
-    unique_bus_pairs = list(OrderedDict((val, None) for idx, val in bus_pairs.items()))
+    unique_bus_pairs = tx_utils.get_unique_bus_pairs(md)
 
     model = pe.ConcreteModel()
 
@@ -287,16 +293,44 @@ def _s_bounds_rule(m, from_bus, to_bus):
 
     model.obj = pe.Objective(expr=obj_expr)
 
+    for gen_name in out_of_service_gens:
+        model.pg[gen_name].fix(0)
+        model.qg[gen_name].fix(0)
+        model_data.data['elements']['generator'][gen_name]['in_service'] = False
+        md.data['elements']['generator'][gen_name]['in_service'] = False
+    for branch_name in out_of_service_branches:
+        model.pf[branch_name].fix(0)
+        model.pt[branch_name].fix(0)
+        model.qf[branch_name].fix(0)
+        model.qt[branch_name].fix(0)
+        model.eq_pf_branch[branch_name].deactivate()
+        model.eq_pt_branch[branch_name].deactivate()
+        model.eq_qf_branch[branch_name].deactivate()
+        model.eq_qt_branch[branch_name].deactivate()
+        model.ineq_sf_branch_thermal_limit[branch_name].deactivate()
+        model.ineq_st_branch_thermal_limit[branch_name].deactivate()
+        model.ineq_angle_diff_branch_lb[branch_name].deactivate()
+        model.ineq_angle_diff_branch_ub[branch_name].deactivate()
+        model_data.data['elements']['branch'][branch_name]['in_service'] = False
+        md.data['elements']['branch'][branch_name]['in_service'] = False
+
+    unique_bus_pairs_set = set(unique_bus_pairs)
+    ubp = set(tx_utils.get_unique_bus_pairs(md))
+    ubp_diff = unique_bus_pairs_set - ubp
+    for from_bus, to_bus in ubp_diff:
+        model.c[from_bus, to_bus].setlb(None)
+        model.c[from_bus, to_bus].setub(None)
+        model.s[from_bus, to_bus].setlb(None)
+        model.s[from_bus, to_bus].setub(None)
+
     return model, md
 
 
-def create_atan_acopf_model(model_data, include_feasibility_slack=False, pw_cost_model='delta'):
+def create_atan_acopf_model(model_data, include_feasibility_slack=False, pw_cost_model='delta', keep_vars_for_out_of_service_elements=False):
     model, md = _create_base_power_ac_model(model_data, include_feasibility_slack=include_feasibility_slack,
-                                            pw_cost_model=pw_cost_model)
+                                            pw_cost_model=pw_cost_model, keep_vars_for_out_of_service_elements=keep_vars_for_out_of_service_elements)
 
-    branch_attrs = md.attributes(element_type='branch')
-    bus_pairs = zip_items(branch_attrs['from_bus'], branch_attrs['to_bus'])
-    unique_bus_pairs = OrderedSet(val for val in bus_pairs.values())
+    unique_bus_pairs = tx_utils.get_unique_bus_pairs(md)
     for fb, tb in unique_bus_pairs:
         assert (tb, fb) not in unique_bus_pairs
 
@@ -339,13 +373,11 @@ def _dva_bounds_rule(m, from_bus, to_bus):
     return model, md
 
 
-def create_psv_acopf_model(model_data, include_feasibility_slack=False, pw_cost_model='delta'):
+def create_psv_acopf_model(model_data, include_feasibility_slack=False, pw_cost_model='delta', keep_vars_for_out_of_service_elements=False):
     model, md = _create_base_power_ac_model(model_data, include_feasibility_slack=include_feasibility_slack,
-                                            pw_cost_model=pw_cost_model)
+                                            pw_cost_model=pw_cost_model, keep_vars_for_out_of_service_elements=keep_vars_for_out_of_service_elements)
     bus_attrs = md.attributes(element_type='bus')
-    branch_attrs = md.attributes(element_type='branch')
-    bus_pairs = zip_items(branch_attrs['from_bus'], branch_attrs['to_bus'])
-    unique_bus_pairs = list(OrderedDict((val, None) for idx, val in bus_pairs.items()).keys())
+    unique_bus_pairs = tx_utils.get_unique_bus_pairs(md)
 
     # declare the polar voltages
     libbranch.declare_var_dva(model=model,
@@ -384,13 +416,11 @@ def create_psv_acopf_model(model_data, include_feasibility_slack=False, pw_cost_
     return model, md
 
 
-def create_rsv_acopf_model(model_data, include_feasibility_slack=False, pw_cost_model='delta'):
+def create_rsv_acopf_model(model_data, include_feasibility_slack=False, pw_cost_model='delta', keep_vars_for_out_of_service_elements=False):
     model, md = _create_base_power_ac_model(model_data, include_feasibility_slack=include_feasibility_slack,
-                                            pw_cost_model=pw_cost_model)
+                                            pw_cost_model=pw_cost_model, keep_vars_for_out_of_service_elements=keep_vars_for_out_of_service_elements)
     bus_attrs = md.attributes(element_type='bus')
-    branch_attrs = md.attributes(element_type='branch')
-    bus_pairs = zip_items(branch_attrs['from_bus'], branch_attrs['to_bus'])
-    unique_bus_pairs = list(OrderedDict((val, None) for idx, val in bus_pairs.items()).keys())
+    unique_bus_pairs = tx_utils.get_unique_bus_pairs(md)
 
     # declare the rectangular voltages
     neg_v_max = map_items(op.neg, bus_attrs['v_max'])

egret/models/dcopf.py

@@ -49,7 +49,15 @@ def _include_feasibility_slack(model, bus_names, bus_p_loads, gens_by_bus, gen_a
                     for bus_name in bus_names)
     return p_rhs_kwargs, penalty_expr
 
-def create_btheta_dcopf_model(model_data, include_angle_diff_limits=False, include_feasibility_slack=False, pw_cost_model='delta'):
+def create_btheta_dcopf_model(model_data, include_angle_diff_limits=False, include_feasibility_slack=False, pw_cost_model='delta',
+                              keep_vars_for_out_of_service_elements=False):
+    if keep_vars_for_out_of_service_elements:
+        out_of_service_gens = tx_utils._get_out_of_service_gens(model_data)
+        out_of_service_branches = tx_utils._get_out_of_service_branches(model_data)
+    else:
+        out_of_service_gens = list()
+        out_of_service_branches = list()
+
     md = model_data.clone_in_service()
     tx_utils.scale_ModelData_to_pu(md, inplace = True)
 
@@ -209,6 +217,21 @@ def create_btheta_dcopf_model(model_data, include_angle_diff_limits=False, inclu
 
     model.obj = pe.Objective(expr=obj_expr)
 
+    for gen_name in out_of_service_gens:
+        model.pg[gen_name].fix(0)
+        model_data.data['elements']['generator'][gen_name]['in_service'] = False
+        md.data['elements']['generator'][gen_name]['in_service'] = False
+    for branch_name in out_of_service_branches:
+        model.pf[branch_name].fix(0)
+        model.eq_pf_branch[branch_name].deactivate()
+        model.ineq_pf_branch_thermal_lb[branch_name].deactivate()
+        model.ineq_pf_branch_thermal_ub[branch_name].deactivate()
+        if include_angle_diff_limits:
+            model.ineq_angle_diff_branch_lb[branch_name].deactivate()
+            model.ineq_angle_diff_branch_ub[branch_name].deactivate()
+        model_data.data['elements']['branch'][branch_name]['in_service'] = False
+        md.data['elements']['branch'][branch_name]['in_service'] = False        
+
     return model, md
 
 def create_ptdf_dcopf_model(model_data, include_feasibility_slack=False, base_point=BasePointType.FLATSTART, ptdf_options=None, pw_cost_model='delta'):

egret/models/tests/test_acopf.py

@@ -14,6 +14,7 @@
 import math
 import unittest
 from pyomo.opt import SolverFactory, TerminationCondition
+import pyomo.environ as pe
 from egret.models.acopf import *
 from egret.data.model_data import ModelData
 from parameterized import parameterized
@@ -84,6 +85,29 @@ def test_acopf_model(self, test_case, soln_case, p_and_v_soln_case, include_kwar
         self.assertTrue(comparison)
         _test_p_and_v(self, p_and_v_soln_case, md)
 
+    def test_keep_vars(self):
+        fname = os.path.join(current_dir, 'transmission_test_instances/pglib-opf-master/pglib_opf_case5_pjm.m')
+        md = ModelData.read(fname)
+        md.data["elements"]["generator"]["1"]["in_service"] = False
+        md.data["elements"]["branch"]["2"]["in_service"] = False
+
+        m1, _ = create_psv_acopf_model(md, keep_vars_for_out_of_service_elements=False)
+        m2, _ = create_psv_acopf_model(md, keep_vars_for_out_of_service_elements=True)
+
+        opt = SolverFactory('ipopt')
+        res1 = opt.solve(m1)
+        res2 = opt.solve(m2)
+
+        self.assertEqual(res1.solver.termination_condition, TerminationCondition.optimal)
+        self.assertEqual(res2.solver.termination_condition, TerminationCondition.optimal)
+
+        obj1 = pe.value(m1.obj)
+        obj2 = pe.value(m2.obj)
+
+        self.assertAlmostEqual(obj1, obj2)
+        self.assertTrue(m2.pg["1"].fixed)
+
+
 class TestArctanACOPF(unittest.TestCase):
     show_output = True
 
@@ -106,6 +130,28 @@ def test_acopf_model(self, test_case, soln_case):
         self.assertTrue(res.solver.termination_condition == TerminationCondition.optimal)
         self.assertAlmostEqual(pe.value(model.obj)/md_soln.data['system']['total_cost'], 1, 4)
 
+    def test_keep_vars(self):
+        fname = os.path.join(current_dir, 'transmission_test_instances/pglib-opf-master/pglib_opf_case5_pjm.m')
+        md = ModelData.read(fname)
+        md.data["elements"]["generator"]["1"]["in_service"] = False
+        md.data["elements"]["branch"]["2"]["in_service"] = False
+
+        m1, _ = create_atan_acopf_model(md, keep_vars_for_out_of_service_elements=False)
+        m2, _ = create_atan_acopf_model(md, keep_vars_for_out_of_service_elements=True)
+
+        opt = SolverFactory('ipopt')
+        res1 = opt.solve(m1)
+        res2 = opt.solve(m2)
+
+        self.assertEqual(res1.solver.termination_condition, TerminationCondition.optimal)
+        self.assertEqual(res2.solver.termination_condition, TerminationCondition.optimal)
+
+        obj1 = pe.value(m1.obj)
+        obj2 = pe.value(m2.obj)
+
+        self.assertAlmostEqual(obj1, obj2)
+        self.assertTrue(m2.pg["1"].fixed)
+
 class TestRSVACOPF(unittest.TestCase):
     show_output = True
 
@@ -134,6 +180,29 @@ def test_acopf_model(self, test_case, soln_case, p_and_v_soln_case, include_kwar
         self.assertTrue(comparison)
         _test_p_and_v(self, p_and_v_soln_case, md)
 
+    def test_keep_vars(self):
+        fname = os.path.join(current_dir, 'transmission_test_instances/pglib-opf-master/pglib_opf_case5_pjm.m')
+        md = ModelData.read(fname)
+        md.data["elements"]["generator"]["1"]["in_service"] = False
+        md.data["elements"]["branch"]["2"]["in_service"] = False
+
+        m1, _ = create_rsv_acopf_model(md, keep_vars_for_out_of_service_elements=False)
+        m2, _ = create_rsv_acopf_model(md, keep_vars_for_out_of_service_elements=True)
+
+        opt = SolverFactory('ipopt')
+        res1 = opt.solve(m1)
+        res2 = opt.solve(m2)
+
+        self.assertEqual(res1.solver.termination_condition, TerminationCondition.optimal)
+        self.assertEqual(res2.solver.termination_condition, TerminationCondition.optimal)
+
+        obj1 = pe.value(m1.obj)
+        obj2 = pe.value(m2.obj)
+
+        self.assertAlmostEqual(obj1, obj2)
+        self.assertTrue(m2.pg["1"].fixed)
+
+
 class TestRIVACOPF(unittest.TestCase):
     show_output = True
 

egret/models/tests/test_dcopf.py

@@ -73,6 +73,28 @@ def test_ptdf_dcopf_model(self, test_case, soln_case, include_kwargs=False):
         comparison = math.isclose(md.data['system']['total_cost'], md_soln.data['system']['total_cost'], rel_tol=1e-6)
         self.assertTrue(comparison)
 
+    def test_keep_vars(self):
+        fname = os.path.join(current_dir, 'transmission_test_instances/pglib-opf-master/pglib_opf_case5_pjm.m')
+        md = ModelData.read(fname)
+        md.data["elements"]["generator"]["1"]["in_service"] = False
+        md.data["elements"]["branch"]["2"]["in_service"] = False
+
+        m1, _ = create_btheta_dcopf_model(md, keep_vars_for_out_of_service_elements=False)
+        m2, _ = create_btheta_dcopf_model(md, keep_vars_for_out_of_service_elements=True)
+
+        opt = SolverFactory('ipopt')
+        res1 = opt.solve(m1)
+        res2 = opt.solve(m2)
+
+        self.assertEqual(res1.solver.termination_condition, TerminationCondition.optimal)
+        self.assertEqual(res2.solver.termination_condition, TerminationCondition.optimal)
+
+        obj1 = pe.value(m1.obj)
+        obj2 = pe.value(m2.obj)
+
+        self.assertAlmostEqual(obj1, obj2)
+        self.assertTrue(m2.pg["1"].fixed)
+
 
 def poly_cost_to_pw_cost(md: ModelData, num_points=10):
     gen_attrs = md.attributes(element_type='generator')