Merge pull request #79 from CURENT/misc

Several fix

ams/models/group.py

@@ -200,7 +200,7 @@ class StaticGen(GroupBase):
     def __init__(self):
         super().__init__()
         self.common_params.extend(('Sn', 'Vn', 'p0', 'q0', 'ra', 'xs', 'subidx',
-                                   'bus', 'pmax', 'pmin', 'pg0', 'ctrl'))
+                                   'bus', 'pmax', 'pmin', 'pg0', 'ctrl', 'R10'))
         self.common_vars.extend(('p', 'q'))
 
 
@@ -217,7 +217,7 @@ class StaticLoad(GroupBase):
 
     def __init__(self):
         super().__init__()
-        self.common_params.extend(('p0', 'q0'))
+        self.common_params.extend(('p0', 'q0', 'zone'))
 
 
 class StaticShunt(GroupBase):

ams/models/line.py

@@ -1,6 +1,6 @@
 from andes.models.line.line import LineData
 from andes.core.param import NumParam
-from andes.shared import deg2rad
+from andes.shared import deg2rad, np, spmatrix
 
 from ams.core.model import Model
 
@@ -12,6 +12,9 @@ class Line(LineData, Model):
     The model is also used for two-winding transformer. Transformers can set the
     tap ratio in ``tap`` and/or phase shift angle ``phi``.
 
+    Note that the bus admittance matrix is built on fly and is not stored in the
+    object.
+
     Notes
     -----
     There is a known issue that adding Algeb ``ud`` will cause Line.algebs run into
@@ -39,3 +42,169 @@ def __init__(self, system=None, config=None) -> None:
         self.rate_a.default = 999.0
         self.rate_b.default = 999.0
         self.rate_c.default = 999.0
+
+        # NOTE: following parameters are prepared for building matrices
+        # they are initialized here but populated in ``System.setup()``.
+        self.a1a = None
+        self.a2a = None
+
+    # NOTE: following code are minly copied from `andes.models.line.Line`
+    # and they are not fully verified
+    # potential issues:
+    # `build_Bp` contains 'fdxb', which is not included in the input parameters,
+    # and the results are the negative of `Bbus` from `makeBdc` in PYPOWER
+    # `build_Bpp` ignores the line resistance for all three methods
+    # `build_Bdc` results are the negative of `Bbus` from `makeBdc` in PYPOWER
+    # `build_y` results have inignorable differences at diagonal elements with `makeYbus` in PYPOWER
+
+    def build_y(self):
+        """
+        Build bus admittance matrix. Copied from ``andes.models.line.line.Line``.
+
+        Returns
+        -------
+        Y : spmatrix
+            Bus admittance matrix.
+        """
+
+        nb = self.system.Bus.n
+
+        y1 = self.u.v * (self.g1.v + self.b1.v * 1j)
+        y2 = self.u.v * (self.g2.v + self.b2.v * 1j)
+        y12 = self.u.v / (self.r.v + self.x.v * 1j)
+        m = self.tap.v * np.exp(1j * self.phi.v)
+        m2 = self.tap.v**2
+        mconj = np.conj(m)
+
+        # build self and mutual admittances into Y
+        Y = spmatrix((y12 + y1 / m2), self.a1a, self.a1a, (nb, nb), 'z')
+        Y -= spmatrix(y12 / mconj, self.a1a, self.a2a, (nb, nb), 'z')
+        Y -= spmatrix(y12 / m, self.a2a, self.a1a, (nb, nb), 'z')
+        Y += spmatrix(y12 + y2, self.a2a, self.a2a, (nb, nb), 'z')
+
+        return Y
+
+    def build_Bp(self, method='fdpf'):
+        """
+        Function for building B' matrix.
+
+        Parameters
+        ----------
+        method : str
+            Method for building B' matrix. Choose from 'fdpf', 'fdbx', 'dcpf'.
+
+        Returns
+        -------
+        Bp : spmatrix
+            B' matrix.
+        """
+        nb = self.system.Bus.n
+
+        if method not in ("fdpf", "fdbx", "dcpf"):
+            raise ValueError(f"Invalid method {method}; choose from 'fdpf', 'fdbx', 'dcpf'")
+
+        # Build B prime matrix -- FDPF
+        # `y1`` neglects line charging shunt, and g1 is usually 0 in HV lines
+        # `y2`` neglects line charging shunt, and g2 is usually 0 in HV lines
+        y1 = self.u.v * self.g1.v
+        y2 = self.u.v * self.g2.v
+
+        # `m` neglected tap ratio
+        m = np.exp(self.phi.v * 1j)
+        mconj = np.conj(m)
+        m2 = np.ones(self.n)
+
+        if method in ('fdxb', 'dcpf'):
+            # neglect line resistance in Bp in XB method
+            y12 = self.u.v / (self.x.v * 1j)
+        else:
+            y12 = self.u.v / (self.r.v + self.x.v * 1j)
+
+        Bdc = spmatrix((y12 + y1) / m2, self.a1a, self.a1a, (nb, nb), 'z')
+        Bdc -= spmatrix(y12 / mconj, self.a1a, self.a2a, (nb, nb), 'z')
+        Bdc -= spmatrix(y12 / m, self.a2a, self.a1a, (nb, nb), 'z')
+        Bdc += spmatrix(y12 + y2, self.a2a, self.a2a, (nb, nb), 'z')
+        Bdc = Bdc.imag()
+
+        for item in range(nb):
+            if abs(Bdc[item, item]) == 0:
+                Bdc[item, item] = 1e-6 + 0j
+
+        return Bdc
+
+    def build_Bpp(self, method='fdpf'):
+        """
+        Function for building B'' matrix.
+
+        Parameters
+        ----------
+        method : str
+            Method for building B'' matrix. Choose from 'fdpf', 'fdbx', 'dcpf'.
+
+        Returns
+        -------
+        Bpp : spmatrix
+            B'' matrix.
+        """
+
+        nb = self.system.Bus.n
+
+        if method not in ("fdpf", "fdbx", "dcpf"):
+            raise ValueError(f"Invalid method {method}; choose from 'fdpf', 'fdbx', 'dcpf'")
+
+        # Build B double prime matrix
+        # y1 neglected line charging shunt, and g1 is usually 0 in HV lines
+        # y2 neglected line charging shunt, and g2 is usually 0 in HV lines
+        # m neglected phase shifter
+        y1 = self.u.v * (self.g1.v + self.b1.v * 1j)
+        y2 = self.u.v * (self.g2.v + self.b2.v * 1j)
+
+        m = self.tap.v
+        m2 = abs(m)**2
+
+        if method in ('fdbx', 'fdpf', 'dcpf'):
+            # neglect line resistance in Bpp in BX method
+            y12 = self.u.v / (self.x.v * 1j)
+        else:
+            y12 = self.u.v / (self.r.v + self.x.v * 1j)
+
+        Bpp = spmatrix((y12 + y1) / m2, self.a1a, self.a1a, (nb, nb), 'z')
+        Bpp -= spmatrix(y12 / np.conj(m), self.a1a, self.a2a, (nb, nb), 'z')
+        Bpp -= spmatrix(y12 / m, self.a2a, self.a1a, (nb, nb), 'z')
+        Bpp += spmatrix(y12 + y2, self.a2a, self.a2a, (nb, nb), 'z')
+        Bpp = Bpp.imag()
+
+        for item in range(nb):
+            if abs(Bpp[item, item]) == 0:
+                Bpp[item, item] = 1e-6 + 0j
+
+        return Bpp
+
+    def build_Bdc(self):
+        """
+        The MATPOWER-flavor Bdc matrix for DC power flow.
+
+        The method neglects line charging and line resistance. It retains tap ratio.
+
+        Returns
+        -------
+        Bdc : spmatrix
+            Bdc matrix.
+        """
+
+        nb = self.system.Bus.n
+
+        y12 = self.u.v / (self.x.v * 1j)
+        y12 = y12 / self.tap.v
+
+        Bdc = spmatrix(y12, self.a1a, self.a1a, (nb, nb), 'z')
+        Bdc -= spmatrix(y12, self.a1a, self.a2a, (nb, nb), 'z')
+        Bdc -= spmatrix(y12, self.a2a, self.a1a, (nb, nb), 'z')
+        Bdc += spmatrix(y12, self.a2a, self.a2a, (nb, nb), 'z')
+        Bdc = Bdc.imag()
+
+        for item in range(nb):
+            if abs(Bdc[item, item]) == 0:
+                Bdc[item, item] = 1e-6
+
+        return Bdc

ams/opt/omodel.py

@@ -6,7 +6,6 @@
 from typing import Any, Optional, Union
 from collections import OrderedDict
 import re
-import ast
 
 import numpy as np
 import scipy.sparse as spr
@@ -573,8 +572,7 @@ def v2(self):
 
         try:
             logger.debug(f"Value code: {code}")
-            out = ast.literal_eval(code)
-            return out
+            return eval(code)
         except Exception as e:
             logger.error(f"Error in calculating constr <{self.name}>.")
             logger.error(f"Original error: {e}")
@@ -661,8 +659,7 @@ def v2(self):
 
         try:
             logger.debug(f"Value code: {code}")
-            out = ast.literal_eval(code)
-            return out
+            return eval(code)
         except Exception as e:
             logger.error(f"Error in calculating obj <{self.name}>.")
             logger.error(f"Original error: {e}")
@@ -746,9 +743,9 @@ def __init__(self, routine):
         self.constrs = OrderedDict()
         self.obj = None
         self.initialized = False
-        self._parsed = False
+        self.parsed = False
 
-    def _parse(self, no_code=True):
+    def parse(self, no_code=True):
         """
         Parse the optimization model from the symbolic description.
 
@@ -815,8 +812,8 @@ def _parse(self, no_code=True):
             else:
                 logger.warning(f"{rtn.class_name} has no objective function!")
                 _, s = elapsed(t0)
-                self._parsed = False
-                return self._parsed
+                self.parsed = False
+                return self.parsed
         _, s = elapsed(t0)
         logger.debug(f"Parse Objective in {s}")
 
@@ -827,8 +824,8 @@ def _parse(self, no_code=True):
         _, s = elapsed(t0)
         logger.debug(f"Parse Expressions in {s}")
 
-        self._parsed = True
-        return self._parsed
+        self.parsed = True
+        return self.parsed
 
     def init(self, no_code=True):
         """
@@ -850,8 +847,8 @@ def init(self, no_code=True):
         """
         t_setup, _ = elapsed()
 
-        if not self._parsed:
-            self._parse(no_code=no_code)
+        if not self.parsed:
+            self.parse(no_code=no_code)
 
         if self.rtn.type == 'PF':
             _, s_setup = elapsed(t_setup)

ams/routines/routine.py

@@ -102,15 +102,19 @@ def get(self, src: str, idx, attr: str = 'v',
         if idx_all is None:
             raise ValueError(f"<{self.class_name}> item <{src}> has no idx.")
 
+        is_format = False  # whether the idx is formatted as a list
+        idx_u = None
         if isinstance(idx, (str, int)):
-            idx = [idx]
-
-        if isinstance(idx, np.ndarray):
-            idx = idx.tolist()
-
-        loc = [idx_all.index(idxe) if idxe in idx_all else None for idxe in idx]
+            idx_u = [idx]
+            is_format = True
+        elif isinstance(idx, np.ndarray):
+            idx_u = idx.tolist()
+        elif isinstance(idx, list):
+            idx_u = idx.copy()
+
+        loc = [idx_all.index(idxe) if idxe in idx_all else None for idxe in idx_u]
         if None in loc:
-            idx_none = [idxe for idxe in idx if idxe not in idx_all]
+            idx_none = [idxe for idxe in idx_u if idxe not in idx_all]
             msg = f"Var <{self.class_name}.{src}> does not contain value with idx={idx_none}"
             raise ValueError(msg)
         out = getattr(item, attr)[loc]
@@ -132,7 +136,8 @@ def get(self, src: str, idx, attr: str = 'v',
             out = out[:, loc_h]
             if out.shape[1] == 1:
                 out = out[:, 0]
-        return out
+
+        return out[0] if is_format else out
 
     def set(self, src: str, idx, attr: str = "v", value=0.0):
         """
@@ -286,13 +291,6 @@ def init(self, force=False, no_code=True, **kwargs):
         logger.info(msg)
         return self.initialized
 
-    def prepare(self):
-        """
-        Prepare the routine.
-        """
-        logger.debug("Generating code for %s", self.class_name)
-        self.syms.generate_symbols()
-
     def solve(self, **kwargs):
         """
         Solve the routine optimization model.
@@ -516,7 +514,7 @@ def update(self, params=None, mat_make=True,):
             if no system matrices are changed.
         """
         t0, _ = elapsed()
-        re_setup = False
+        re_init = False
         # sanitize input
         sparams = []
         if params is None:
@@ -532,12 +530,13 @@ def update(self, params=None, mat_make=True,):
                 param.update()
         for param in sparams:
             if param.optz is None:  # means no_parse=True
-                re_setup = True
+                re_init = True
                 break
         if mat_make:
             self.system.mats.make()
-        if re_setup:
+        if re_init:
             logger.warning(f"<{self.class_name}> reinit OModel due to non-parametric change.")
+            self.om.parsed = False
             _ = self.om.init(no_code=True)
         results = self.om.update(params=sparams)
         t0, s0 = elapsed(t0)
@@ -663,7 +662,7 @@ def _post_add_check(self):
         # --- reset optimization model status ---
         self.om.initialized = False
         # --- reset OModel parser status ---
-        self.om._parsed = False
+        self.om.parsed = False
 
     def addRParam(self,
                   name: str,

ams/system.py

@@ -429,6 +429,11 @@ def setup(self):
                 algeb.a = np.arange(a0, a0 + algeb.owner.n)
                 a0 += algeb.owner.n
 
+        # NOTE: this is a temporary solution for building Y matrix
+        # consider refator this part if any other similar cases occur in the future
+        self.Line.a1a = self.Bus.get(src='a', attr='a', idx=self.Line.bus1.v)
+        self.Line.a2a = self.Bus.get(src='a', attr='a', idx=self.Line.bus2.v)
+
         _, s = elapsed(t0)
         logger.info('System set up in %s.', s)
 

docs/source/release-notes.rst

@@ -15,6 +15,10 @@ v0.9.6 (2024-xx-xx)
 - Refactor DCPF, PFlow, and ACOPF
 - Add a loss factor in ``RTED.dc2ac()``
 - Add ``DCOPF.dc2ac()``
+- Fix OModel parse status to ensure no_parsed params can be updated
+- Fix and rerun ex2
+- Format ``Routine.get()`` return type to be consistent with input idx type
+- Remove unused ``Routine.prepare()``
 
 v0.9.5 (2024-03-25)
 -------------------