Merge pull request #81 from CURENT/develop

Prep for patch release v0.9.6

ams/interop/andes.py

@@ -56,8 +56,10 @@
              'pq': 'PQ', }
 
 
-def to_andes(system, setup=False, addfile=None,
-             **kwargs):
+def to_andes(system, addfile=None,
+             setup=False, no_output=False,
+             default_config=True,
+             verify=True, tol=1e-3):
     """
     Convert the AMS system to an ANDES system.
 
@@ -80,12 +82,19 @@ def to_andes(system, setup=False, addfile=None,
     ----------
     system : System
         The AMS system to be converted to ANDES format.
-    setup : bool, optional
-        Whether to call `setup()` after the conversion. Default is True.
     addfile : str, optional
         The additional file to be converted to ANDES dynamic mdoels.
-    **kwargs : dict
-        Keyword arguments to be passed to `andes.system.System`.
+    setup : bool, optional
+        Whether to call `setup()` after the conversion. Default is True.
+    no_output : bool, optional
+        To ANDES system.
+    default_config : bool, optional
+        To ANDES system.
+    verify : bool
+        If True, the converted ANDES system will be verified with the source
+        AMS system using AC power flow.
+    tol : float
+        The tolerance of error.
 
     Returns
     -------
@@ -97,9 +106,8 @@ def to_andes(system, setup=False, addfile=None,
     >>> import ams
     >>> import andes
     >>> sp = ams.load(ams.get_case('ieee14/ieee14_uced.xlsx'), setup=True)
-    >>> sa = sp.to_andes(setup=False,
-    ...                  addfile=andes.get_case('ieee14/ieee14_full.xlsx'),
-    ...                  overwrite=True, no_output=True)
+    >>> sa = sp.to_andes(addfile=andes.get_case('ieee14/ieee14_full.xlsx'),
+    ...                  setup=False, overwrite=True, no_output=True)
 
     Notes
     -----
@@ -109,13 +117,15 @@ def to_andes(system, setup=False, addfile=None,
     """
     t0, _ = elapsed()
 
-    adsys = andes_System()
+    adsys = andes_System(no_output=no_output,
+                         default_config=default_config)
     # FIXME: is there a systematic way to do this? Other config might be needed
     adsys.config.freq = system.config.freq
     adsys.config.mva = system.config.mva
 
     for mdl_name, mdl_cols in pflow_dict.items():
         mdl = getattr(system, mdl_name)
+        mdl.cache.refresh("df_in")  # refresh cache
         for row in mdl.cache.df_in[mdl_cols].to_dict(orient='records'):
             adsys.add(mdl_name, row)
 
@@ -139,8 +149,14 @@ def to_andes(system, setup=False, addfile=None,
     # finalize
     system.dyn = Dynamic(amsys=system, adsys=adsys)
     system.dyn.link_andes(adsys=adsys)
+
     if setup:
         adsys.setup()
+    elif verify:
+        logger.warning('PFlow verification is skipped due to no setup.')
+        return adsys
+    if verify:
+        verify_pf(amsys=system, adsys=adsys, tol=tol)
     return adsys
 
 
@@ -531,7 +547,7 @@ def send(self, adsys=None, routine=None):
             logger.info(f'Send <{rtn.class_name}> results to ANDES <{hex(id(sa))}>...')
 
         # NOTE: if DC type, check if results are converted
-        if (rtn.type != 'ACED') and (not rtn.is_ac):
+        if (rtn.type != 'ACED') and (not rtn.converted):
             logger.error(f'<{rtn.class_name}> AC conversion failed or not done yet!')
 
         # --- Mapping ---
@@ -553,16 +569,26 @@ def send(self, adsys=None, routine=None):
             idx_ads = var_ams.get_idx()  # use AMS idx as target ANDES idx
 
             # --- special scenarios ---
+            # 0. send PV bus voltage to StaticGen.v0 if not PFlow yet and AC converted
+            cond_vpv = (mname_ads == 'Bus') and (pname_ads == 'v0')
+            if cond_vpv and (not self.is_tds) and (rtn.converted):
+                # --- StaticGen ---
+                stg_idx = sp.StaticGen.get_idx()
+                bus_stg = sp.StaticGen.get(src='bus', attr='v', idx=stg_idx)
+                vBus = rtn.get(src='vBus', attr='v', idx=bus_stg)
+                sa.StaticGen.set(value=vBus, idx=stg_idx, src='v0', attr='v')
+                logger.info(f'*Send <{vname_ams}> to StaticGen.v0')
+
             # 1. gen online status; in TDS running, setting u is invalid
             cond_ads_stg_u = (mname_ads in ['StaticGen', 'PV', 'Sclak']) and (pname_ads == 'u')
             if cond_ads_stg_u and (self.is_tds):
-                logger.info(f'Skip sending {vname_ams} to StaticGen.u during TDS')
+                logger.info(f'*Skip sending {vname_ams} to StaticGen.u during TDS')
                 continue
 
             # 2. Bus voltage
             cond_ads_bus_v0 = (mname_ads == 'Bus') and (pname_ads == 'v0')
             if cond_ads_bus_v0 and (self.is_tds):
-                logger.info(f'Skip sending {vname_ams} t0 Bus.v0 during TDS')
+                logger.info(f'*Skip sending {vname_ams} t0 Bus.v0 during TDS')
                 continue
 
             # 3. gen power reference; in TDS running, pg should go to TurbineGov
@@ -602,7 +628,7 @@ def send(self, adsys=None, routine=None):
                     var_dest = 'TurbineGov.pref0'
                 if len(dg_ams) > 0:
                     var_dest += ' and DG.pref0'
-                logger.warning(f'Send <{vname_ams}> to {var_dest}')
+                logger.warning(f'*Send <{vname_ams}> to {var_dest}')
                 continue
 
             # --- other scenarios ---
@@ -743,7 +769,7 @@ def receive(self, adsys=None, routine=None, no_update=False):
         if no_update and (len(pname_to_update) > 0):
             logger.info(f'Please update <{rtn.class_name}> parameters: {pname_to_update}')
         elif len(pname_to_update) > 0:
-            rtn.update(params=pname_to_update, mat_make=False)
+            rtn.update(params=pname_to_update, build_mats=False)
         return True
 
 
@@ -933,3 +959,43 @@ def make_link_table(adsys):
             ]
     out = ssa_key[cols].sort_values(by='stg_idx', ascending=False).reset_index(drop=True)
     return out
+
+
+def verify_pf(amsys, adsys, tol=1e-3):
+    """
+    Verify the power flow results between AMS and ANDES.
+    Note that this function will run PFlow in both systems.
+
+    Parameters
+    ----------
+    sp : ams.System
+        The AMS system.
+    sa : andes.System
+        The ANDES system.
+
+    Returns
+    -------
+    bool
+        True if the power flow results are consistent; False otherwise.
+    """
+    amsys.PFlow.run()
+    if adsys.is_setup:
+        adsys.PFlow.run()
+    else:
+        logger.info('ANDES system is not setup, quit verification.')
+        return False
+    v_check = np.allclose(amsys.Bus.v.v, adsys.Bus.v.v, atol=tol)
+    a_check = np.allclose(amsys.Bus.a.v, adsys.Bus.a.v, atol=tol)
+    check = v_check and a_check
+
+    v_diff_max = np.max(np.abs(amsys.Bus.v.v - adsys.Bus.v.v))
+    a_diff_max = np.max(np.abs(amsys.Bus.a.v - adsys.Bus.a.v))
+    diff_msg = f'Voltage diff max: {v_diff_max}, Angle diff max: {a_diff_max}'
+    logger.debug(diff_msg)
+    if check:
+        logger.info('Power flow results are consistent.')
+    else:
+        msg = 'Power flow results are inconsistent!'
+        logger.warning(msg)
+        logger.warning(diff_msg)
+    return check

ams/models/group.py

@@ -118,7 +118,7 @@ class RenGen(GroupBase):
 
     def __init__(self):
         super().__init__()
-        self.common_params.extend(('bus', 'gen', 'Sn'))
+        self.common_params.extend(('bus', 'gen', 'Sn', 'q0'))
         self.common_vars.extend(('Pe', 'Qe'))
 
 
@@ -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',))
+        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