option to keep variables for out of service elements

.github/workflows/egret.yml

@@ -29,14 +29,14 @@ jobs:
             matrix:
                 os: [ubuntu-latest]
                 python-version: [3.7, 3.8, 3.9]
-                pyomo-version: [6.1.2]
+                pyomo-version: [6.4.0]
                 include:
                     - os: macos-latest
                       python-version: 3.7
-                      pyomo-version: 6.1.2
+                      pyomo-version: 6.4.0
                     - os: windows-latest
                       python-version: 3.7
-                      pyomo-version: 6.1.2
+                      pyomo-version: 6.4.0
                     - os: ubuntu-latest
                       python-version: 3.7
                       pyomo-version: main
@@ -86,13 +86,13 @@ jobs:
                   #  - ipopt needs: libopenblas-dev gfortran liblapack-dev
                   sudo apt-get install libopenblas-dev libgfortran4 libgomp1 gfortran liblapack-dev
                   curl --max-time 150 --retry 8 \
-                    -L $URL/idaes-solvers-ubuntu1804-64.tar.gz \
+                    -L $URL/idaes-solvers-ubuntu1804-x86_64.tar.gz \
                     > $IPOPT_TAR
                 else
                   # windows
                    curl --max-time 150 --retry 8 \
-                    -L $URL/idaes-solvers-windows-64.tar.gz \
-                    $URL/idaes-lib-windows-64.tar.gz > $IPOPT_TAR
+                    -L $URL/idaes-solvers-windows-x86_64.tar.gz \
+                    $URL/idaes-lib-windows-x86_64.tar.gz > $IPOPT_TAR
                 fi
                 cd $IPOPT_DIR
                 tar -xzi < $IPOPT_TAR

README.md

@@ -25,7 +25,7 @@ EGRET is available under the BSD License (see [LICENSE.txt](https://github.com/g
 ### Requirements
 
 * Python 3.7 or later
-* Pyomo version 6.1.2 or later
+* Pyomo version 6.4.0 or later
 * pytest
 * Optimization solvers for Pyomo - specific requirements depends on the models being solved. EGRET is tested with Gurobi or CPLEX for MIP-based problems (e.g., unit commitment) and Ipopt (with HSL linear solvers) for NLP problems.
 

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,88 @@ def _s_bounds_rule(m, from_bus, to_bus):
 
     model.obj = pe.Objective(expr=obj_expr)
 
+    out_of_service_gens_set = set(out_of_service_gens)
+
+    for gen_name, e in model.pg_operating_cost.items():
+        if gen_name in out_of_service_gens_set:
+            e.expr = 0
+
+    if len(pw_pg_cost_gens) > 0:
+        if pw_cost_model == 'delta':
+            for gen_name, ndx in model.delta_pg_set:
+                if gen_name in out_of_service_gens_set:
+                    model.delta_pg[gen_name, ndx].set_value(0, skip_validation=True)
+                    model.delta_pg[gen_name, ndx].fix()
+                    model.pg_delta_pg_con[gen_name].deactivate()
+        else:
+            for gen_name, ndx in model.pg_piecewise_cost_set:
+                if gen_name in out_of_service_gens_set:
+                    model.pg_cost[gen_name].set_value(0, skip_validation=True)
+                    model.pg_cost[gen_name].fix()
+                    model.pg_piecewise_cost_cons[gen_name, ndx].deactivate()
+
+    if q_costs is not None:
+        for gen_name, e in model.qg_operating_cost.items():
+            if gen_name in out_of_service_gens_set:
+                e.expr = 0
+        if len(pw_qg_cost_gens) > 0:
+            if pw_cost_model == 'delta':
+                for gen_name, ndx in model.delta_qg_set:
+                    if gen_name in out_of_service_gens_set:
+                        model.delta_qg[gen_name, ndx].set_value(0, skip_validation=True)
+                        model.delta_qg[gen_name, ndx].fix()
+                        model.qg_delta_qg_con[gen_name].deactivate()
+            else:
+                for gen_name, ndx in model.qg_piecewise_cost_set:
+                    if gen_name in out_of_service_gens_set:
+                        model.qg_cost[gen_name].set_value(0, skip_validation=True)
+                        model.qg_cost[gen_name].fix()
+                        model.qg_piecewise_cost_cons[gen_name, ndx].deactivate()
+
+    for gen_name in out_of_service_gens:
+        model.pg[gen_name].set_value(0, skip_validation=True)
+        model.qg[gen_name].set_value(0, skip_validation=True)
+        model.pg[gen_name].fix()
+        model.qg[gen_name].fix()
+        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].set_value(0, skip_validation=True)
+        model.pt[branch_name].set_value(0, skip_validation=True)
+        model.qf[branch_name].set_value(0, skip_validation=True)
+        model.qt[branch_name].set_value(0, skip_validation=True)
+        model.pf[branch_name].fix()
+        model.pt[branch_name].fix()
+        model.qf[branch_name].fix()
+        model.qt[branch_name].fix()
+        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 +417,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 +460,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,43 @@ def create_btheta_dcopf_model(model_data, include_angle_diff_limits=False, inclu
 
     model.obj = pe.Objective(expr=obj_expr)
 
+    out_of_service_gens_set = set(out_of_service_gens)
+
+    for gen_name, e in model.pg_operating_cost.items():
+        if gen_name in out_of_service_gens_set:
+            e.expr = 0
+
+    if len(pw_pg_cost_gens) > 0:
+        if pw_cost_model == 'delta':
+            for gen_name, ndx in model.delta_pg_set:
+                if gen_name in out_of_service_gens_set:
+                    model.delta_pg[gen_name, ndx].set_value(0, skip_validation=True)
+                    model.delta_pg[gen_name, ndx].fix()
+                    model.pg_delta_pg_con[gen_name].deactivate()
+        else:
+            for gen_name, ndx in model.pg_piecewise_cost_set:
+                if gen_name in out_of_service_gens_set:
+                    model.pg_cost[gen_name].set_value(0, skip_validation=True)
+                    model.pg_cost[gen_name].fix()
+                    model.pg_piecewise_cost_cons[gen_name, ndx].deactivate()
+
+    for gen_name in out_of_service_gens:
+        model.pg[gen_name].set_value(0, skip_validation=True)
+        model.pg[gen_name].fix()
+        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].set_value(0, skip_validation=True)
+        model.pf[branch_name].fix()
+        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'):