Merge pull request #72 from CURENT/misc

Quick fix on DCOPF series

ams/cases/5bus/pjm5bus_uced.json

@@ -337,7 +337,7 @@
       "trans": 0.0,
       "tap": 1.0,
       "phi": 0.0,
-      "rate_a": 0.0,
+      "rate_a": 2.0,
       "rate_b": 0.0,
       "rate_c": 0.0,
       "owner": null,
@@ -367,7 +367,7 @@
       "trans": 0.0,
       "tap": 1.0,
       "phi": 0.0,
-      "rate_a": 0.0,
+      "rate_a": 2.0,
       "rate_b": 0.0,
       "rate_c": 0.0,
       "owner": null,
@@ -397,7 +397,7 @@
       "trans": 0.0,
       "tap": 1.0,
       "phi": 0.0,
-      "rate_a": 0.0,
+      "rate_a": 2.0,
       "rate_b": 0.0,
       "rate_c": 0.0,
       "owner": null,
@@ -427,7 +427,7 @@
       "trans": 0.0,
       "tap": 1.0,
       "phi": 0.0,
-      "rate_a": 0.0,
+      "rate_a": 2.0,
       "rate_b": 0.0,
       "rate_c": 0.0,
       "owner": null,

ams/io/matpower.py

@@ -374,6 +374,8 @@ def system2mpc(system) -> dict:
         branch[:, 8] = system.Line.tap.v
         branch[:, 9] = system.Line.phi.v * rad2deg
         branch[:, 10] = system.Line.u.v
+        branch[:, 11] = system.Line.amin.v * rad2deg
+        branch[:, 12] = system.Line.amax.v * rad2deg
 
     # --- GCost ---
     # NOTE: adjust GCost sequence to match the generator sequence

ams/models/__init__.py

@@ -8,7 +8,7 @@
 ams_file_classes = list([
     ('info', ['Summary']),
     ('bus', ['Bus']),
-    ('static', ['PQ', 'PV', 'Slack']),
+    ('static', ['PQ', 'Slack', 'PV']),
     ('shunt', ['Shunt']),
     ('line', ['Line']),
     ('distributed', ['PVD1', 'ESD1']),

ams/routines/dcopf.py

@@ -85,11 +85,14 @@ def __init__(self, system, config):
                              name='rate_a', tex_name=r'R_{ATEA}',
                              unit='MVA', model='Line',)
         # --- line angle difference ---
-        self.amax = NumOp(u=self.rate_a, fun=np.ones_like,
-                          rfun=np.dot, rargs=dict(b=np.pi),
-                          name='amax', tex_name=r'\theta_{bus, max}',
-                          info='max line angle difference',
-                          no_parse=True,)
+        self.amax = RParam(model='Line', src='amax',
+                           name='amax', tex_name=r'\theta_{bus, max}',
+                           info='max line angle difference',
+                           no_parse=True,)
+        self.amin = RParam(model='Line', src='amin',
+                           name='amin', tex_name=r'\theta_{bus, min}',
+                           info='min line angle difference',
+                           no_parse=True,)
         # --- shunt ---
         self.gsh = RParam(info='shunt conductance',
                           name='gsh', tex_name=r'g_{sh}',
@@ -169,7 +172,7 @@ def __init__(self, system, config):
                                 e_str='Bf@aBus + Pfinj - rate_a',)
         self.alflb = Constraint(info='line angle difference lower bound',
                                 name='alflb', is_eq=False,
-                                e_str='-CftT@aBus - amax',)
+                                e_str='-CftT@aBus + amin',)
         self.alfub = Constraint(info='line angle difference upper bound',
                                 name='alfub', is_eq=False,
                                 e_str='CftT@aBus - amax',)

ams/routines/ed.py

@@ -150,6 +150,8 @@ def __init__(self, system, config):
 
         self.dud.expand_dims = 1
         self.ddd.expand_dims = 1
+        self.amin.expand_dims = 1
+        self.amax.expand_dims = 1
 
         # --- Data Section ---
         self.ugt = NumOp(u=self.ug, fun=np.transpose,
@@ -182,7 +184,7 @@ def __init__(self, system, config):
         self.plf.info = '2D Line flow'
         self.plflb.e_str = '-Bf@aBus - Pfinj@tlv - rate_a@tlv'
         self.plfub.e_str = 'Bf@aBus + Pfinj@tlv - rate_a@tlv'
-        self.alflb.e_str = '-CftT@aBus - amax@tlv'
+        self.alflb.e_str = '-CftT@aBus + amin@tlv'
         self.alfub.e_str = 'CftT@aBus - amax@tlv'
 
         self.plfc.e_str = 'Bf@aBus + Pfinj@tlv'

ams/routines/uc.py

@@ -140,6 +140,8 @@ def __init__(self, system, config):
         self.sd.model = 'UCTSlot'
 
         self.ug.expand_dims = 1
+        self.amin.expand_dims = 1
+        self.amax.expand_dims = 1
 
         # --- Model Section ---
         # --- gen ---

ams/system.py

@@ -402,16 +402,18 @@ def setup(self):
         # --- model parameters range check ---
         # TODO: there might be other parameters check?
         # Line rate
-        adjusted_rate = []
+        adjust_rate = []
         default_rate = 999
         for rate in [self.Line.rate_a, self.Line.rate_b, self.Line.rate_c]:
-            rate.v[rate.v == 0] = default_rate
-            adjusted_rate.append(rate.name)
-        adjusted_rate = ', '.join(adjusted_rate)
-        msg = f"Zero line rates detacted in {adjusted_rate}, "
-        msg += f"adjusted to {default_rate}."
-        msg += "\nIf expect a line outage, please set 'u' to 0."
-        logger.info(msg)
+            if np.all(rate.v == 0):
+                adjust_rate.append(rate.name)
+                rate.v[rate.v == 0] = default_rate
+        if adjust_rate:
+            adjusted_rate = ', '.join(adjust_rate)
+            msg = f"Zero line rates detacted in {adjusted_rate}, "
+            msg += f"adjusted to {default_rate}."
+            msg += "\nIf expect a line outage, please set 'u' to 0."
+            logger.info(msg)
         # === no device addition or removal after this point ===
         self.calc_pu_coeff()   # calculate parameters in system per units
 

docs/source/release-notes.rst

@@ -13,6 +13,11 @@ v0.9.5 (2024-xx-xx)
 -------------------
 
 - Add more plots in demo_AGC
+- Improve line rating adjustment
+- Adjust static import sequence in `models.__init__.py`
+- Adjust pjm5bus case line rate_a
+- Fix formulation of constraint line angle diff
+- Add amin amax conversion to MATPOWER converter
 
 v0.9.4 (2024-03-16)
 -------------------

tests/test_addressing.py

@@ -27,14 +27,14 @@ def test_ieee14_address(self):
 
         # external variable indices
         np.testing.assert_array_equal(ss.PV.ud.a,
-                                      np.array([28, 29, 30, 31]))
+                                      np.array([31, 32, 33, 34]))
         np.testing.assert_array_equal(ss.PV.p.a,
-                                      np.array([32, 33, 34, 35]))
+                                      np.array([35, 36, 37, 38]))
         np.testing.assert_array_equal(ss.PV.q.a,
-                                      np.array([36, 37, 38, 39]))
+                                      np.array([39, 40, 41, 42]))
         np.testing.assert_array_equal(ss.Slack.ud.a,
-                                      np.array([40]))
+                                      np.array([28]))
         np.testing.assert_array_equal(ss.Slack.p.a,
-                                      np.array([41]))
+                                      np.array([29]))
         np.testing.assert_array_equal(ss.Slack.q.a,
-                                      np.array([42]))
+                                      np.array([30]))

tests/test_known_good.py

@@ -1,23 +1,52 @@
 import logging
 import unittest
 
+import numpy as np
+
 import ams
 
 logger = logging.getLogger(__name__)
 
 
 class TestKnownResults(unittest.TestCase):
-    # NOTE: DCOPF objective values are from MATPOWER v8.0b1
-    sets = (
-        (ams.get_case('matpower/case14.m'), 7642.59177699),
-        (ams.get_case('matpower/case39.m'), 41263.9486962204),
-        (ams.get_case('matpower/case118.m'), 125947.8814179),
-    )
-
-    def test_known_results(self):
-        for case, obj in self.sets:
-            sp = ams.load(case, setup=True, no_output=True)
-            sp.DCOPF.run(solver='ECOS')
-            msg = f'Case: {case}; Result: {sp.DCOPF.obj.v}; Expected: {obj}'
-            logger.info(msg)
-            self.assertAlmostEqual(sp.DCOPF.obj.v, obj, places=2)
+    # NOTE: Baseline comes from MATPOWER v8.0b1
+    # NOTE: when using ECOS in this case, the accuracy is ralatively low
+
+    def test_known_case14(self):
+        case = ams.get_case('matpower/case14.m')
+        sp = ams.load(case, setup=True, no_output=True)
+        sp.DCOPF.run(solver='ECOS')
+        self.assertAlmostEqual(sp.DCOPF.obj.v, 7642.59177715644, places=2)
+
+        aBus_mp = np.array([
+            -0.0, -0.088451789, -0.22695322, -0.18549695,
+            -0.15942929, -0.25995018, -0.24348912, -0.24348912,
+            -0.27468284, -0.2795552, -0.27334388, -0.27941265,
+            -0.28242718, -0.30074116,])
+        aBus_mp -= aBus_mp[0]
+
+        aBus = sp.DCOPF.aBus.v - sp.DCOPF.aBus.v[0]
+        self.assertTrue(np.allclose(aBus, aBus_mp, atol=1e-4))
+
+        pg_mp = np.array([
+            2.2096769, 0.38032305, 0.0, 0.0, 0.0])
+        pg = sp.DCOPF.pg.v
+        self.assertTrue(np.allclose(pg, pg_mp, atol=1e-4))
+
+    def test_known_case39(self):
+        case = ams.get_case('matpower/case39.m')
+        sp = ams.load(case, setup=True, no_output=True)
+        sp.DCOPF.run(solver='ECOS')
+        self.assertAlmostEqual(sp.DCOPF.obj.v, 41263.94078592735, places=2)
+
+        pi = sp.DCOPF.pi.v / sp.config.mva
+        self.assertTrue(np.all(np.isclose(pi, 13.51692, atol=1e-2)))
+
+    def test_known_case118(self):
+        case = ams.get_case('matpower/case118.m')
+        sp = ams.load(case, setup=True, no_output=True)
+        sp.DCOPF.run(solver='ECOS')
+        self.assertAlmostEqual(sp.DCOPF.obj.v, 125947.8814181522, places=2)
+
+        pi = sp.DCOPF.pi.v / sp.config.mva
+        self.assertTrue(np.all(np.isclose(pi, 39.38136794580036, atol=1e-2)))