Merge pull request #201 from cuihantao/develop

Prepare for v1.5.0 release.

TODO.md

@@ -1,12 +1,15 @@
 # v1.5
 
-# Milestones
+## Milestones
 
 * Separate array-enabled parameters from NumParam into
   `BaseConvParam`, which implements `iconvert` and `oconvert`,
   and then `ListParam`
 * Improve the robustness of accessing the fields of `pycode` by properly
   handling `KeyError`s.
+* Numba code generation for all models. Add a function for triggering
+  all the JIT compilation before simulation.
+* Unit Model Test framework.
 
 # Later Versions
 
@@ -25,11 +28,13 @@
 * Eigenvalue analysis report sort options: by damping, frequency, or eigenvalue
 * Root loci plots
 
-# v1.4
+# Previous
+
+## v1.4
 
 * [X] Disallow linking ExtAlgeb to State and ExtState to Algeb (check in prepare).
 
-# v1.3
+## v1.3
 
 ## Milestones
 

andes/core/block.py

@@ -211,8 +211,9 @@ def define(self):
 
     def export(self):
         """
-        Method for exporting instances defined in this class in a dictionary. This method calls the ``define``
-        method first and returns ``self.vars``.
+        Method for exporting instances defined in this class in a dictionary.
+
+        This method calls the ``define`` method first and returns ``self.vars``.
 
         Returns
         -------

andes/core/common.py

@@ -148,8 +148,7 @@ def append_ijv(self, j_full_name, ii, jj, vv):
         vv : array-like
             Value indices
         """
-        if len(ii) == 0 and len(jj) == 0:
-            return
+
         self.ijac[j_full_name].append(ii)
         self.jjac[j_full_name].append(jj)
         self.vjac[j_full_name].append(vv)

andes/core/model.py

@@ -622,7 +622,9 @@ def __init__(self, system=None, config=None):
 
         self.tex_names = OrderedDict((('dae_t', 't_{dae}'),
                                       ('sys_f', 'f_{sys}'),
-                                      ('sys_mva', 'S_{b,sys}')
+                                      ('sys_mva', 'S_{b,sys}'),
+                                      ('__ones', 'O_{nes}'),
+                                      ('__zeros', 'Z_{eros}'),
                                       ))
 
         # Model behavior flags
@@ -792,6 +794,17 @@ def _one_idx2uid(self, idx):
 
         return self.uid[idx]
 
+    def set_backref(self, name, from_idx, to_idx):
+        """
+        Helper function for setting idx-es to ``BackRef``.
+        """
+
+        if name not in self.services_ref:
+            return
+
+        uid = self.idx2uid(to_idx)
+        self.services_ref[name].v[uid].append(from_idx)
+
     def get(self, src: str, idx, attr: str = 'v', allow_none=False, default=0.0):
         """
         Get the value of an attribute of a model property.
@@ -955,6 +968,12 @@ def refresh_inputs(self):
         self._input['sys_f'] = self.system.config.freq
         self._input['sys_mva'] = self.system.config.mva
 
+        # two vectors with the length of the number of devices.
+        #   Useful in the choices of `PieceWise`, which need to be vectors
+        #   for numba to compile.
+        self._input['__ones'] = np.ones(self.n)
+        self._input['__zeros'] = np.zeros(self.n)
+
     def refresh_inputs_arg(self):
         """
         Refresh inputs for each function with individual argument list.
@@ -1263,16 +1282,10 @@ def j_update(self):
             for idx, fun in enumerate(self.calls.vjac[jname]):
                 try:
                     self.triplets.vjac[jname][idx][:] = ret[idx]
-                except ValueError as e:
+                except (ValueError, IndexError, FloatingPointError) as e:
                     row_name, col_name = self._jac_eq_var_name(jname, idx)
-                    logger.error('%s shape error: j_idx=%s, d%s / d%s',
-                                 jname, idx, row_name, col_name)
-
-                    raise e
-                except FloatingPointError as e:
-                    row_name, col_name = self._jac_eq_var_name(jname, idx)
-                    logger.error('%s eval error: j_idx=%s, d%s / d%s',
-                                 jname, idx, row_name, col_name)
+                    logger.error('%s: error calculating or storing Jacobian <%s>: j_idx=%s, d%s / d%s',
+                                 self.class_name, jname, idx, row_name, col_name)
 
                     raise e
 
@@ -1673,12 +1686,20 @@ def init(self, routine):
             kwargs = self.get_inputs(refresh=True)
 
             for idx, name in enumerate(self.calls.init_seq):
-                # single variable
+                # single variable - do assignment
                 if isinstance(name, str):
                     instance = self.__dict__[name]
                     _eval_discrete(instance)
                     if instance.v_str is not None:
-                        instance.v[:] = self.calls.ia[name](*self.ia_args[name])
+                        if not instance.v_str_add:
+                            # assignment is for most variable initialization
+                            instance.v[:] = self.calls.ia[name](*self.ia_args[name])
+                        else:
+                            # in-place add initial values.
+                            # Voltage compensators can set part of the `v` of exciters.
+                            # Exciters will then set the bus voltage part.
+                            instance.v[:] += self.calls.ia[name](*self.ia_args[name])
+
                     # single variable iterative solution
                     if name in self.calls.ii:
                         self.solve_iter(name, kwargs)

andes/core/symprocessor.py

@@ -5,8 +5,9 @@
 import os
 import logging
 import pprint
-
+import sympy
 import numpy as np
+
 from collections import OrderedDict, defaultdict
 
 from sympy import Symbol, Matrix
@@ -16,6 +17,10 @@
 from andes.core.npfunc import safe_div
 from andes.shared import dilled_vars
 from andes.utils.paths import get_pycode_path
+from andes.utils.sympy import FixPiecewise
+
+
+sympy.Piecewise = FixPiecewise
 
 logger = logging.getLogger(__name__)
 
@@ -129,10 +134,14 @@ def generate_symbols(self):
             if var in self.parent.__dict__ and self.parent.__dict__[var].tex_name is not None:
                 self.tex_names[Symbol(var)] = Symbol(self.parent.__dict__[var].tex_name)
 
+        # additional variables by conventions
         self.inputs_dict['dae_t'] = Symbol('dae_t')
         self.inputs_dict['sys_f'] = Symbol('sys_f')
         self.inputs_dict['sys_mva'] = Symbol('sys_mva')
+        self.inputs_dict['__ones'] = Symbol('__ones')
+        self.inputs_dict['__zeros'] = Symbol('__zeros')
 
+        # custom functions
         self.lambdify_func[0]['Indicator'] = lambda x: x
         self.lambdify_func[0]['imag'] = np.imag
         self.lambdify_func[0]['real'] = np.real
@@ -221,7 +230,9 @@ def generate_services(self):
             s_args[name] = [str(i) for i in expr.free_symbols]
             s_calls[name] = lambdify(s_args[name], s_syms[name], modules=self.lambdify_func)
 
+        # TODO: below triggers DeprecationWarning with SymPy 1.9
         self.s_matrix = Matrix(list(s_syms.values()))
+
         self.calls.s = s_calls
         self.calls.s_args = s_args
 
@@ -371,6 +382,12 @@ def generate_pycode(self, pycode_path, yapf_pycode):
         Create output source code file for generated code.
 
         Generated code are stored at ``~/.andes/pycode``.
+
+        Notes
+        -----
+        In the current implementation, each model saves a ``.py`` file.
+        In systems with slow disk access (such as networked file systems),
+        this function can be the bottleneck.
         """
 
         pycode_path = get_pycode_path(pycode_path, mkdir=True)
@@ -379,6 +396,9 @@ def generate_pycode(self, pycode_path, yapf_pycode):
         header = \
             """from collections import OrderedDict  # NOQA
 
+import numpy
+
+
 from numpy import nan, pi, sin, cos, tan, sqrt, exp, select         # NOQA
 from numpy import greater_equal, less_equal, greater, less, equal   # NOQA
 from numpy import logical_and, logical_or, logical_not              # NOQA

andes/core/var.py

@@ -50,6 +50,14 @@ class BaseVar:
         local-storage of the corresponding equation value
     e_str : str
         the string/symbolic representation of the equation
+    v_str : str
+        explicit initialization equation
+    v_str_add : bool
+        True if the value of `v_str` will be added to the variable.
+        Useful when other models access this variable and set part
+        of the initial value
+    v_iter : str
+        implicit iterative equation in the form of 0 = v_iter
     """
 
     def __init__(self,
@@ -63,6 +71,7 @@ def __init__(self,
                  discrete: Optional[Discrete] = None,
                  v_setter: Optional[bool] = False,
                  e_setter: Optional[bool] = False,
+                 v_str_add: Optional[bool] = False,
                  addressable: Optional[bool] = True,
                  export: Optional[bool] = True,
                  diag_eps: Optional[float] = 0.0,
@@ -84,6 +93,7 @@ def __init__(self,
         self.discrete = discrete
         self.v_setter = v_setter        # True if this variable sets the variable value
         self.e_setter = e_setter        # True if this var sets the equation value
+        self.v_str_add = v_str_add
 
         self.addressable = addressable  # True if this var needs to be assigned an address FIXME: not in use
         self.export = export            # True if this var's value needs to exported
@@ -411,10 +421,16 @@ def set_arrays(self, dae, inplace=True, alloc=True):
         when ``e_str`` exists..
         """
 
-        if self.e_str is None:
+        if self.e_str is None or (self.n == 0):
             return
 
-        slice_idx = slice(self.r[0], self.r[-1] + 1)
+        try:
+            slice_idx = slice(self.r[0], self.r[-1] + 1)
+        except IndexError as e:
+            print(self.owner.class_name)
+            print(self.name)
+            raise e
+
         if isinstance(self, ExtState):
             self.e = dae.h[slice_idx]
         elif isinstance(self, ExtAlgeb):
@@ -451,6 +467,7 @@ def link_external(self, ext_model):
         self.parent = ext_model
 
         if isinstance(ext_model, GroupBase):
+            # determine the number of elements based on `indexer.v`
             if self.indexer.n > 0 and isinstance(self.indexer.v[0], (list, np.ndarray)):
                 self._n = [len(i) for i in self.indexer.v]  # number of elements in each sublist
                 self._idx = np.concatenate([np.array(i) for i in self.indexer.v])

andes/models/__init__.py

@@ -23,10 +23,12 @@
     ('area', ['Area', 'ACE', 'ACEc']),
     ('dynload', ['ZIP', 'FLoad']),
     ('synchronous', ['GENCLS', 'GENROU']),
-    ('governor', ['TG2', 'TGOV1', 'TGOV1N', 'TGOV1DB', 'IEEEG1', 'IEESGO']),
+    ('governor', ['TG2', 'TGOV1', 'TGOV1DB', 'TGOV1N', 'TGOV1NDB',
+                  'IEEEG1', 'IEESGO']),
+    ('vcomp', ['IEEEVC']),
     ('exciter', ['EXDC2', 'IEEEX1', 'ESDC2A', 'EXST1', 'ESST3A', 'SEXS',
                  'IEEET1', 'EXAC1', 'EXAC4', 'ESST4B', 'AC8B', 'IEEET3',
-                 'ESAC1A']),
+                 'ESAC1A', 'ESST1A']),
     ('pss', ['IEEEST', 'ST2CUT']),
     ('motor', ['Motor3', 'Motor5']),
     ('measurement', ['BusFreq', 'BusROCOF', 'PMU']),

andes/models/exciter/__init__.py

@@ -11,5 +11,6 @@
 from andes.models.exciter.ac8b import AC8B      # NOQA
 from andes.models.exciter.ieeet3 import IEEET3      # NOQA
 from andes.models.exciter.esac1a import ESAC1A      # NOQA
+from andes.models.exciter.esst1a import ESST1A      # NOQA
 
-from andes.models.exciter.saturation import ExcQuadSat, ExcExpSat  # NOQA
+from andes.models.exciter.saturation import ExcQuadSat, ExcExpSat  # NOQA

andes/models/exciter/ac8b.py

@@ -46,33 +46,31 @@ def __init__(self):
                            )
 
         self.VPMAX = NumParam(info='PID maximum limit',
-                              tex_name='V_{PMAX}',
+                              tex_name=r'V_{PMAX}',
                               default=999,
                               unit='p.u.')
         self.VPMIN = NumParam(info='PID minimum limit',
-                              tex_name='V_{PMIN}',
+                              tex_name=r'V_{PMIN}',
                               default=-999,
                               unit='p.u.')
 
-        self.VRMAX = NumParam(info='Maximum excitation limit',
-                              tex_name='V_{RMAX}',
+        self.VRMAX = NumParam(info='Maximum regulator limit',
+                              tex_name=r'V_{RMAX}',
                               default=7.3,
                               unit='p.u.',
                               vrange=(1, 10))
-        self.VRMIN = NumParam(info='Minimum excitation limit',
-                              tex_name='V_{RMIN}',
+        self.VRMIN = NumParam(info='Minimum regulator limit',
+                              tex_name=r'V_{RMIN}',
                               default=1,
                               unit='p.u.',
                               vrange=(-1, 1.5))
 
-        # TODO: check default value for VFEMAX
-        self.VFEMAX = NumParam(info='Maximum VFE',
+        self.VFEMAX = NumParam(info='Exciter field current limit',
                                tex_name=r'V_{FEMAX}',
                                default=999,
                                unit='p.u.')
 
-        # TODO: check default value for VEMIN
-        self.VEMIN = NumParam(info='Minimum excitation output',
+        self.VEMIN = NumParam(info='Minimum exciter voltage output',
                               tex_name=r'V_{EMIN}',
                               default=-999,
                               unit='p.u.')
@@ -99,7 +97,7 @@ def __init__(self):
                            unit='p.u.',
                            )
         self.SE1 = NumParam(info='Value at first saturation point',
-                            tex_name='S_{E1}',
+                            tex_name=r'S_{E1}',
                             default=0.,
                             unit='p.u.',
                             )
@@ -109,7 +107,7 @@ def __init__(self):
                            unit='p.u.',
                            )
         self.SE2 = NumParam(info='Value at second saturation point',
-                            tex_name='S_{E2}',
+                            tex_name=r'S_{E2}',
                             default=1.,
                             unit='p.u.',
                             )
@@ -169,6 +167,8 @@ def __init__(self, system, config):
 
         ExcVsum.__init__(self)
 
+        self.vref.v_str = 'v'
+
         self.vi = Algeb(info='Total input voltages',
                         tex_name='V_i',
                         unit='p.u.',