adding option to use SOC overestimators from paper

egret/models/ac_relaxations.py

@@ -2,19 +2,24 @@
 import egret.model_library.transmission.branch as libbranch
 from egret.data.data_utils import map_items, zip_items
 from collections import OrderedDict
-try:
-    import coramin
-    coramin_available = True
-except ImportError:
-    coramin_available = False
+import coramin
+from coramin.relaxations.custom_block import declare_custom_block
+from coramin.relaxations.relaxations_base import BaseRelaxationData, ComponentWeakRef
+from pyomo.core.expr.numvalue import NumericConstant, is_constant
+import pyomo.environ as pe
+import math
+from math import sqrt
+import logging
+from pyomo.common.collections.orderedset import OrderedSet
 
 
-def _relaxation_helper(model, md, include_soc, use_linear_relaxation):
-    if not coramin_available:
-        raise ImportError('Cannot create relaxation unless coramin is available.')
+logger = logging.getLogger(__name__)
+
+
+def _relaxation_helper(model, md, include_soc, use_linear_relaxation, use_fbbt=True):
     coramin.relaxations.relax(model,
                               in_place=True,
-                              use_fbbt=True,
+                              use_fbbt=use_fbbt,
                               fbbt_options={'deactivate_satisfied_constraints': True,
                                             'max_iter': 2})
     if not use_linear_relaxation:
@@ -31,10 +36,195 @@ def _relaxation_helper(model, md, include_soc, use_linear_relaxation):
                                    use_outer_approximation=use_linear_relaxation)
 
 
-def create_soc_relaxation(model_data, use_linear_relaxation=True, include_feasibility_slack=False):
+def _check_bounds(v, lb, ub):
+    if lb is None:
+        lb = -math.inf
+    if ub is None:
+        ub = math.inf
+    if ub < lb:
+        raise ValueError('Variable lb is greater than ub')
+    if not v.is_fixed() and (ub - lb) <= 1e-8:
+        logger.warning('Variable is not fixed, but (ub - lb) is small; var: {0}; lb: {1}; ub: {2}'.format(str(v), lb, ub))
+
+
+def _get_bounds(v):
+    if v.is_fixed():
+        lb = ub = v.value
+    else:
+        lb = v.lb
+        ub = v.ub
+    _check_bounds(v, lb, ub)
+    return lb, ub
+
+
+@declare_custom_block(name='SOCEdgeCuts')
+class SOCEdgeCutsData(BaseRelaxationData):
+    """
+    Relaxation for
+
+      0 >= vmsq_1 * vmsq_2 - c**2 - s**2
+
+    based on
+
+    Kocuk, B., Dey, S.S. & Sun, X.A. Matrix minor reformulation
+    and SOCP-based spatial branch-and-cut method for the AC
+    optimal power flow problem. Math. Prog. Comp. 10, 557–596
+    (2018). https://doi.org/10.1007/s12532-018-0150-9
+    """
+    def __init__(self, component):
+        BaseRelaxationData.__init__(self, component)
+        self._cref = ComponentWeakRef(None)
+        self._sref = ComponentWeakRef(None)
+        self._vmsq_1_ref = ComponentWeakRef(None)
+        self._vmsq_2_ref = ComponentWeakRef(None)
+        self._aux_var = NumericConstant(0)
+
+    @property
+    def _c(self):
+        return self._cref.get_component()
+
+    @property
+    def _s(self):
+        return self._sref.get_component()
+
+    @property
+    def _vmsq_1(self):
+        return self._vmsq_1_ref.get_component()
+
+    @property
+    def _vmsq_2(self):
+        return self._vmsq_2_ref.get_component()
+
+    def get_rhs_vars(self):
+        return [self._c, self._s, self._vmsq_1, self._vmsq_2]
+
+    def get_rhs_expr(self):
+        return self._vmsq_1 * self._vmsq_2 - self._c**2 - self._s**2
+
+    def vars_with_bounds_in_relaxation(self):
+        return self.get_rhs_vars()
+
+    def set_input(self, c, s, vmsq_1, vmsq_2, persistent_solvers=None):
+        self._set_input(relaxation_side=coramin.utils.RelaxationSide.UNDER,
+                        persistent_solvers=persistent_solvers,
+                        use_linear_relaxation=True)
+        self._cref.set_component(c)
+        self._sref.set_component(s)
+        self._vmsq_1_ref.set_component(vmsq_1)
+        self._vmsq_2_ref.set_component(vmsq_2)
+
+    def build(self, c, s, vmsq_1, vmsq_2, persistent_solvers=None):
+        self.set_input(c=c, s=s, vmsq_1=vmsq_1, vmsq_2=vmsq_2, persistent_solvers=persistent_solvers)
+        self.rebuild()
+
+    def _build_relaxation(self):
+        clb, cub = _get_bounds(self._c)
+        slb, sub = _get_bounds(self._s)
+        vmsq_1_lb, vmsq_1_ub = _get_bounds(self._vmsq_1)
+        vmsq_2_lb, vmsq_2_ub = _get_bounds(self._vmsq_2)
+
+        if None in {clb, cub, slb, sub, vmsq_1_lb, vmsq_1_ub, vmsq_2_lb, vmsq_2_ub}:
+            return None
+
+        if vmsq_1_lb == vmsq_1_ub and vmsq_2_lb == vmsq_2_ub:
+            if clb == cub or slb == sub:
+                rhs = [vmsq_1_lb * vmsq_2_lb]
+            else:
+                rhs = [math.sqrt(vmsq_1_lb * vmsq_2_lb)]
+        elif vmsq_1_lb == vmsq_1_ub:
+            if clb == cub or slb == sub:
+                rhs = [vmsq_1_lb * self._vmsq_2]
+            else:
+                m = (math.sqrt(vmsq_2_ub) - math.sqrt(vmsq_2_lb)) / (vmsq_2_ub - vmsq_2_lb)
+                b = math.sqrt(vmsq_2_ub) - m * vmsq_2_ub
+                rhs = [math.sqrt(vmsq_1_lb) * (m * self._vmsq_2 + b)]
+        elif vmsq_2_lb == vmsq_2_ub:
+            if clb == cub or slb == sub:
+                rhs = [vmsq_2_lb * self._vmsq_1]
+            else:
+                m = (math.sqrt(vmsq_1_ub) - math.sqrt(vmsq_1_lb)) / (vmsq_1_ub - vmsq_1_lb)
+                b = math.sqrt(vmsq_1_ub) - m * vmsq_1_ub
+                rhs = [math.sqrt(vmsq_2_lb) * (m * self._vmsq_1 + b)]
+        else:
+            if clb == cub or slb == sub:
+                rhs = [vmsq_1_lb * self._vmsq_2 + self._vmsq_1 * vmsq_2_lb - vmsq_1_lb * vmsq_2_lb,
+                       vmsq_1_ub * self._vmsq_2 + self._vmsq_1 * vmsq_2_ub - vmsq_1_ub * vmsq_2_ub]
+            else:
+                cm = [sqrt(vmsq_1_lb) / (sqrt(vmsq_2_lb) + sqrt(vmsq_2_ub)),
+                      sqrt(vmsq_1_ub) / (sqrt(vmsq_2_lb) + sqrt(vmsq_2_ub))]
+                bm = [sqrt(vmsq_2_lb) / (sqrt(vmsq_1_lb) + sqrt(vmsq_1_ub)),
+                      sqrt(vmsq_2_ub) / (sqrt(vmsq_1_lb) + sqrt(vmsq_1_ub))]
+                am = [sqrt(vmsq_1_lb * vmsq_2_lb) - bm[0] * vmsq_1_lb - cm[0] * vmsq_2_lb,
+                      sqrt(vmsq_1_ub * vmsq_2_ub) - bm[1] * vmsq_1_ub - cm[1] * vmsq_2_ub]
+                rhs = list()
+                for i in [0, 1]:
+                    rhs.append(am[i] + bm[i] * self._vmsq_1 + cm[i] * self._vmsq_2)
+
+        if clb == cub and slb == sub:
+            lhs = [clb**2 + slb**2]
+        elif clb == cub:
+            m = (sub**2 - slb**2) / (sub - slb)
+            b = sub**2 - m * sub
+            lhs = [clb**2 + m * self._s + b]
+        elif slb == sub:
+            m = (cub**2 - clb**2) / (cub - clb)
+            b = cub**2 - m * cub
+            lhs = [m * self._c + b + slb**2]
+        else:
+            if sqrt(cub ** 2 + sub ** 2) + sqrt(clb ** 2 + slb ** 2) - sqrt(cub ** 2 + slb ** 2) - sqrt(
+                    clb ** 2 + sub ** 2) <= 0:
+                cn = [(sqrt(clb ** 2 + sub ** 2) - sqrt(clb ** 2 + slb ** 2)) / (sub - slb),
+                      (sqrt(cub ** 2 + sub ** 2) - sqrt(cub ** 2 + slb ** 2)) / (sub - slb)]
+                bn = [(sqrt(cub ** 2 + slb ** 2) - sqrt(clb ** 2 + slb ** 2)) / (cub - clb),
+                      (sqrt(cub ** 2 + sub ** 2) - sqrt(clb ** 2 + sub ** 2)) / (cub - clb)]
+                an = [sqrt(clb ** 2 + slb ** 2) - bn[0] * clb - cn[0] * slb,
+                      sqrt(cub ** 2 + sub ** 2) - bn[1] * cub - cn[1] * sub]
+            else:
+                cn = [(sqrt(clb ** 2 + sub ** 2) - sqrt(clb ** 2 + slb ** 2)) / (sub - slb),
+                      (sqrt(cub ** 2 + sub ** 2) - sqrt(cub ** 2 + slb ** 2)) / (sub - slb)]
+                bn = [(sqrt(cub ** 2 + sub ** 2) - sqrt(clb ** 2 + sub ** 2)) / (cub - clb),
+                      (sqrt(cub ** 2 + slb ** 2) - sqrt(clb ** 2 + slb ** 2)) / (cub - clb)]
+                an = [sqrt(clb ** 2 + slb ** 2) - bn[0] * clb - cn[0] * slb,
+                      sqrt(cub ** 2 + sub ** 2) - bn[1] * cub - cn[1] * sub]
+            lhs = list()
+            for i in [0, 1]:
+                lhs.append(an[i] + bn[i] * self._c + cn[i] * self._s)
+
+        self.underestimators = pe.ConstraintList()
+        for _lhs in lhs:
+            for _rhs in rhs:
+                if is_constant(_lhs) and is_constant(_rhs):
+                    continue
+                self.underestimators.add(_lhs >= _rhs)
+
+    def is_rhs_convex(self):
+        return False
+
+    def is_rhs_concave(self):
+        return False
+
+    @property
+    def use_linear_relaxation(self):
+        return self._use_linear_relaxation
+
+    @use_linear_relaxation.setter
+    def use_linear_relaxation(self, val):
+        if val is not True:
+            raise ValueError('The SOCEdgeCuts class can only produce linear relaxations')
+        self._use_linear_relaxation = True
+
+
+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)
     if use_linear_relaxation:
-        _relaxation_helper(model=model, md=md, include_soc=True, use_linear_relaxation=use_linear_relaxation)
+        _relaxation_helper(model=model,
+                           md=md,
+                           include_soc=True,
+                           use_linear_relaxation=use_linear_relaxation,
+                           use_fbbt=use_fbbt)
     else:
         branch_attrs = md.attributes(element_type='branch')
         bus_pairs = zip_items(branch_attrs['from_bus'], branch_attrs['to_bus'])
@@ -44,21 +234,59 @@ def create_soc_relaxation(model_data, use_linear_relaxation=True, include_feasib
     return model, md
 
 
-def create_atan_relaxation(model_data, use_linear_relaxation=True, include_feasibility_slack=False):
+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)
     del model.ineq_soc
     del model._con_ineq_soc
-    _relaxation_helper(model=model, md=md, include_soc=True, use_linear_relaxation=use_linear_relaxation)
+    if use_soc_edge_cuts:
+        del model.ineq_soc_ub
+        del model._con_ineq_soc_ub
+    _relaxation_helper(model=model,
+                       md=md,
+                       include_soc=True,
+                       use_linear_relaxation=use_linear_relaxation,
+                       use_fbbt=use_fbbt)
+    if use_soc_edge_cuts:
+        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())
+        model.ineq_soc_ub_set = pe.Set(initialize=list(unique_bus_pairs))
+        model.ineq_soc_ub = SOCEdgeCuts(model.ineq_soc_ub_set)
+        for from_bus, to_bus in model.ineq_soc_ub_set:
+            model.ineq_soc_ub[from_bus, to_bus].build(c=model.c[from_bus, to_bus],
+                                                      s=model.s[from_bus, to_bus],
+                                                      vmsq_1=model.vmsq[from_bus],
+                                                      vmsq_2=model.vmsq[to_bus])
     return model, md
 
 
-def create_polar_acopf_relaxation(model_data, include_soc=True, use_linear_relaxation=True, include_feasibility_slack=False):
+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)
-    _relaxation_helper(model=model, md=md, include_soc=include_soc, use_linear_relaxation=use_linear_relaxation)
+    _relaxation_helper(model=model,
+                       md=md,
+                       include_soc=include_soc,
+                       use_linear_relaxation=use_linear_relaxation,
+                       use_fbbt=use_fbbt)
     return model, md
 
 
-def create_rectangular_acopf_relaxation(model_data, include_soc=True, use_linear_relaxation=True, include_feasibility_slack=False):
+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)
-    _relaxation_helper(model=model, md=md, include_soc=include_soc, use_linear_relaxation=use_linear_relaxation)
+    _relaxation_helper(model=model,
+                       md=md,
+                       include_soc=include_soc,
+                       use_linear_relaxation=use_linear_relaxation,
+                       use_fbbt=use_fbbt)
     return model, md

egret/models/acopf.py

@@ -21,7 +21,7 @@
 import egret.model_library.transmission.gen as libgen
 from egret.model_library.defn import FlowType, CoordinateType
 from egret.data.data_utils import map_items, zip_items
-from math import pi, radians, degrees
+from math import pi, radians, degrees, cos, sin
 from collections import OrderedDict
 from egret.data.networkx_utils import get_networkx_graph
 import networkx
@@ -106,8 +106,52 @@ def _create_base_power_ac_model(model_data, include_feasibility_slack=False):
                             initialize={k: v**2 for k, v in bus_attrs['vm'].items()},
                             bounds=zip_items({k: v**2 for k, v in bus_attrs['v_min'].items()},
                                              {k: v**2 for k, v in bus_attrs['v_max'].items()}))
-    libbranch.declare_var_c(model=model, index_set=unique_bus_pairs, initialize=1)
-    libbranch.declare_var_s(model=model, index_set=unique_bus_pairs, initialize=0)
+
+    branch_w_index = {(v['from_bus'], v['to_bus']): v for v in branches.values()}
+
+    def _c_bounds_rule(m, from_bus, to_bus):
+        bdat = branch_w_index[(from_bus, to_bus)]
+        if bdat['angle_diff_max'] < 0:
+            the_mid = bdat['angle_diff_max']
+        elif bdat['angle_diff_min'] > 0:
+            the_mid = bdat['angle_diff_min']
+        else:
+            the_mid = 0
+        lb = (bus_attrs['v_min'][from_bus]
+              * bus_attrs['v_min'][to_bus]
+              * cos(max(abs(bdat['angle_diff_min']),
+                        abs(bdat['angle_diff_max']))*pi/180))
+        ub = (bus_attrs['v_max'][from_bus]
+              *  bus_attrs['v_max'][to_bus]
+              *  cos(the_mid*pi/180))
+        return lb, ub
+    libbranch.declare_var_c(model=model,
+                            index_set=unique_bus_pairs,
+                            initialize=1,
+                            bounds=_c_bounds_rule)
+
+    def _s_bounds_rule(m, from_bus, to_bus):
+        bdat = branch_w_index[(from_bus, to_bus)]
+
+        if bdat['angle_diff_min'] >= 0:
+            lb = (bus_attrs['v_min'][from_bus]
+                  * bus_attrs['v_min'][to_bus]
+                  * sin(bdat['angle_diff_min']*pi/180))
+        else:
+            lb = (bus_attrs['v_max'][from_bus]
+                  * bus_attrs['v_max'][to_bus]
+                  * sin(bdat['angle_diff_min']*pi/180))
+        if bdat['angle_diff_max'] >= 0:
+            ub = (bus_attrs['v_max'][from_bus]
+                  * bus_attrs['v_max'][to_bus]
+                  * sin(bdat['angle_diff_max']*pi/180))
+        else:
+            ub = (bus_attrs['v_min'][from_bus]
+                  * bus_attrs['v_min'][to_bus]
+                  * sin(bdat['angle_diff_max']*pi/180))
+        return lb, ub
+    libbranch.declare_var_s(model=model, index_set=unique_bus_pairs, initialize=0,
+                            bounds=_s_bounds_rule)
 
     ### include the feasibility slack for the bus balances
     p_rhs_kwargs = {}
@@ -277,10 +321,15 @@ def create_atan_acopf_model(model_data, include_feasibility_slack=False):
         else:
             cycle_basis_bus_pairs.add((b2, b1))
 
+    branches = dict(md.elements(element_type='branch'))
+    branch_w_index = {(v['from_bus'], v['to_bus']): v for v in branches.values()}
+    def _dva_bounds_rule(m, from_bus, to_bus):
+        bdat = branch_w_index[(from_bus, to_bus)]
+        return max(-pi/2, bdat['angle_diff_min'] * pi / 180), min(pi/2, bdat['angle_diff_max'] * pi / 180)
     libbranch.declare_var_dva(model=model,
                               index_set=list(cycle_basis_bus_pairs),
                               initialize=0,
-                              bounds=(-pi/2, pi/2))
+                              bounds=_dva_bounds_rule)
     libbranch.declare_eq_dva_arctan(model=model, index_set=list(cycle_basis_bus_pairs))
     libbranch.declare_eq_dva_cycle_sum(model=model, cycle_basis=cycle_basis, valid_bus_pairs=cycle_basis_bus_pairs)
     libbranch.declare_ineq_soc(model=model, index_set=list(unique_bus_pairs), use_outer_approximation=False)