Repository: CURENT/andes (language: python)

Issue:
Fixed some minor function docstring errors in andes/core/model/model.py
Fixed docstring errors in model.py file


Discussion / hints:

Fix patch:
diff --git a/andes/core/common.py b/andes/core/common.py
index 7f543f8a4..ef676d329 100644
--- a/andes/core/common.py
+++ b/andes/core/common.py
@@ -209,18 +209,18 @@ def add(self, dct=None, **kwargs):
 
         Existing configs will NOT be overwritten.
         """
-        def warn_upper(s):
+        def warn_upper_case(s):
             if any(x.isupper() for x in s):
                 logger.warning("Config fields must be in lower case, found %s", s)
 
         if dct is not None:
             for s in dct.keys():
-                warn_upper(s)
+                warn_upper_case(s)
 
             self._add(**dct)
 
         for s in kwargs.keys():
-            warn_upper(s)
+            warn_upper_case(s)
 
         self._add(**kwargs)
 
@@ -245,21 +245,37 @@ def add_extra(self, dest, dct=None, **kwargs):
             self.__dict__[dest][key] = value
 
     def _add(self, **kwargs):
+        """
+        Internal function for adding new config keys.
+
+        This function does not perform input data consistency check.
+        """
+
         for key, val in kwargs.items():
             # skip existing entries that are already loaded (from config files)
             if key in self.__dict__:
                 continue
 
-            if isinstance(val, str):
+            self._set(key, val)
+
+    def _set(self, key, val):
+        """
+        Set a pair of key and value to the config dict.
+
+        This function does not perform consistency check on input data and will
+        not warn of non-existent fields.
+        """
+
+        if isinstance(val, str):
+            try:
+                val = int(val)
+            except ValueError:
                 try:
-                    val = int(val)
+                    val = float(val)
                 except ValueError:
-                    try:
-                        val = float(val)
-                    except ValueError:
-                        pass
+                    pass
 
-            self.__dict__[key] = val
+        self.__dict__[key] = val
 
     def as_dict(self, refresh=False):
         """
@@ -343,3 +359,18 @@ def check(self):
     @property
     def tex_names(self):
         return self._tex
+
+    def update(self, dct: dict = None, **kwargs):
+        """
+        Update existing configuration fields.
+
+        Values are checked after the update.
+        """
+
+        if dct is not None:
+            kwargs.update(dct)
+
+        for key, val in kwargs.items():
+            self._set(key, val)
+
+        self.check()
diff --git a/andes/core/model/model.py b/andes/core/model/model.py
index 97024c694..eb99cb53b 100644
--- a/andes/core/model/model.py
+++ b/andes/core/model/model.py
@@ -45,7 +45,7 @@ class Model:
 
 
     After subclassing `ModelData`, subclass `Model`` to complete a DAE model.
-    Subclasses of `Model` defines DAE variables, services, and other types of parameters,
+    Subclasses of `Model` define DAE variables, services, and other types of parameters,
     in the constructor ``__init__``.
 
     Attributes
@@ -56,7 +56,7 @@ class Model:
 
     Examples
     --------
-    Take the static PQ as an example, the subclass of `Model`, `PQ`, should looks like ::
+    Take the static PQ as an example, the subclass of `Model`, `PQ`, should look like ::
 
         class PQ(PQData, Model):
             def __init__(self, system, config):
@@ -77,7 +77,7 @@ def __init__(self, system, config):
     neither used in power flow nor time-domain simulation. **A very common pitfall is forgetting to set the flag**.
 
     Next, the group name can be provided. A group is a collection of models with common parameters and variables.
-    Devices idx of all models in the same group must be unique. To provide a group name, use ::
+    Devices' idx of all models in the same group must be unique. To provide a group name, use ::
 
         self.group = 'StaticLoad'
 
@@ -87,15 +87,15 @@ def __init__(self, system, config):
     If not provided with a group class name, the model will be placed in the `Undefined` group.
 
     Next, additional configuration flags can be added.
-    Configuration flags for models are load-time variables specifying the behavior of a model.
-    It can be exported to an `andes.rc` file and automatically loaded when creating the `System`.
+    Configuration flags for models are load-time variables, specifying the behavior of a model.
+    They can be exported to an `andes.rc` file and automatically loaded when creating the `System`.
     Configuration flags can be used in equation strings, as long as they are numerical values.
     To add config flags, use ::
 
         self.config.add(OrderedDict((('pq2z', 1), )))
 
     It is recommended to use `OrderedDict` instead of `dict`, although the syntax is verbose.
-    Note that booleans should be provided as integers (1, or 0), since `True` or `False` is interpreted as
+    Note that booleans should be provided as integers (1 or 0), since `True` or `False` is interpreted as
     a string when loaded from the `rc` file and will cause an error.
 
     Next, it's time for variables and equations! The `PQ` class does not have internal variables itself.
@@ -117,17 +117,17 @@ def __init__(self, system, config):
             self.v.e_str = "u * q"
 
     where the `e_str` attribute is the equation string attribute. `u` is the connectivity status.
-    Any parameter, config, service or variables can be used in equation strings.
+    Any parameter, config, service or variable can be used in equation strings.
 
     Three additional scalars can be used in equations:
-    - ``dae_t`` for the current simulation time can be used if the model has flag `tds`.
+    - ``dae_t`` for the current simulation time (can be used if the model has flag `tds`).
     - ``sys_f`` for system frequency (from ``system.config.freq``).
     - ``sys_mva`` for system base mva (from ``system.config.mva``).
 
     The above example is overly simplified. Our `PQ` model wants a feature to switch itself to
     a constant impedance if the voltage is out of the range `(vmin, vmax)`.
     To implement this, we need to introduce a discrete component called `Limiter`, which yields three arrays
-    of binary flags, `zi`, `zl`, and `zu` indicating in range, below lower limit, and above upper limit,
+    of binary flags, `zi`, `zl`, and `zu` indicating in-range, below lower-limit, and above upper-limit,
     respectively.
 
     First, create an attribute `vcmp` as a `Limiter` instance ::
@@ -149,7 +149,7 @@ def __init__(self, system, config):
     Note that `PQ.a.e_str` can use the three variables from `vcmp` even before defining `PQ.vcmp`, as long as
     `PQ.vcmp` is defined, because `vcmp_zi` is just a string literal in `e_str`.
 
-    The two equations above implements a piecewise power injection equation. It selects the original power demand
+    The two equations above implement a piece-wise power injection equation. It selects the original power demand
     if within range, and uses the calculated power when out of range.
 
     Finally, to let ANDES pick up the model, the model name needs to be added to `models/__init__.py`.
@@ -408,7 +408,7 @@ def idx2uid(self, idx):
 
     def _one_idx2uid(self, idx):
         """
-        Helper function for checking if an idx exist and
+        Helper function for checking if an idx exists and
         converting it to uid.
         """
 
@@ -526,8 +526,8 @@ def get_inputs(self, refresh=False):
         ----------
         refresh : bool
             Refresh the values in the dictionary.
-            This is only used when the memory address of arrays changed.
-            After initialization, all array assignments are inplace.
+            This is only used when the memory addresses of arrays change.
+            After initialization, all array assignments are in place.
             To avoid overhead, refresh should not be used after initialization.
 
         Returns
@@ -898,7 +898,7 @@ def f_update(self):
         Notes
         -----
         In-place equations: added to the corresponding DAE array.
-        Non-inplace equations: in-place set to internal array to
+        Non-in-place equations: in-place set to internal array to
         overwrite old values (and avoid clearing).
         """
         if callable(self.calls.f):
@@ -1153,7 +1153,7 @@ def set_in_use(self):
         Models with no back references will have internal variable addresses assigned but external addresses
         being empty.
 
-        For internal equations that has external variables, the row indices will be non-zeros, while the col
+        For internal equations that have external variables, the row indices will be non-zeros, while the col
         indices will be empty, which causes an error when updating Jacobians.
 
         Setting `self.in_use` to False when `len(back_ref_instance.v) == 0` avoids this error. See COI.
@@ -1322,7 +1322,7 @@ def numba_jitify(self, parallel=False, cache=True, nopython=False):
         Warnings
         --------
         This feature is experimental and does not guarantee a speed up.
-        In fact, the program will likely end up slower due to compilation.
+        In fact, the program will likely end up being slower due to compilation.
         """
         if self.system.config.numba != 1:
             return
diff --git a/andes/main.py b/andes/main.py
index d370c3dfb..3f81a6e26 100644
--- a/andes/main.py
+++ b/andes/main.py
@@ -231,7 +231,6 @@ def remove_output(recursive=False):
                     file.endswith('_out.lst') or \
                     file.endswith('_out.npy') or \
                     file.endswith('_out.npz') or \
-                    file.endswith('_out.csv') or \
                     file.endswith('_prof.prof') or \
                     file.endswith('_prof.txt'):
                 found = True
diff --git a/andes/models/__init__.py b/andes/models/__init__.py
index c7a9bded2..f08d2d99e 100644
--- a/andes/models/__init__.py
+++ b/andes/models/__init__.py
@@ -11,6 +11,7 @@
 
 file_classes = list([
     ('info', ['Summary']),
+    ('misc', ['Output']),
     ('timer', ['Toggler', 'Fault', 'Alter']),
     ('timeseries', ['TimeSeries']),
     ('bus', ['Bus']),
diff --git a/andes/models/group.py b/andes/models/group.py
index 4cefb87ea..5cf0afe2e 100644
--- a/andes/models/group.py
+++ b/andes/models/group.py
@@ -772,3 +772,13 @@ class Information(GroupBase):
     def __init__(self):
         GroupBase.__init__(self)
         self.common_params = []
+
+
+class OutputSelect(GroupBase):
+    """
+    Group for selecting outputs.
+    """
+
+    def __init__(self):
+        super().__init__()
+        self.common_params = []
diff --git a/andes/models/misc/__init__.py b/andes/models/misc/__init__.py
new file mode 100644
index 000000000..9fcf6d1f5
--- /dev/null
+++ b/andes/models/misc/__init__.py
@@ -0,0 +1,5 @@
+"""
+Package for miscellaneous models.
+"""
+
+from andes.models.misc.output import Output  # NOQA
diff --git a/andes/models/misc/output.py b/andes/models/misc/output.py
new file mode 100644
index 000000000..23082220a
--- /dev/null
+++ b/andes/models/misc/output.py
@@ -0,0 +1,61 @@
+"""
+Module for specifying output variables as part of the data file.
+"""
+
+import numpy as np
+from andes.core.model import ModelData, Model
+from andes.core.param import DataParam
+
+
+class OutputData(ModelData):
+    """
+    Data for outputs.
+    """
+
+    def __init__(self):
+        ModelData.__init__(self, three_params=False)
+
+        self.model = DataParam(info='Name of the model', mandatory=True)
+        self.varname = DataParam(info='Variable name', )
+        self.dev = DataParam(info='Device name', )
+
+
+class Output(OutputData, Model):
+    """
+    Model for specifying output models and/or variables.
+    """
+
+    def __init__(self, system, config):
+        OutputData.__init__(self)
+        Model.__init__(self, system, config)
+
+        self.group = 'OutputSelect'
+        self.xidx = []
+        self.yidx = []
+
+    def in1d(self, addr, v_code):
+        """
+        Helper function for finding boolean flags to indicate intersections.
+
+        Parameters
+        ----------
+        idx : array-like
+            indices to find
+        v_code : str
+            variable code in 'x' and 'y'
+        """
+
+        if v_code == 'x':
+            return np.in1d(self.xidx, addr)
+        if v_code == 'y':
+            return np.in1d(self.yidx, addr)
+
+        raise NotImplementedError("v_code <%s> not recognized" % v_code)
+
+    def to_output_addr(self, addr, v_code):
+        """
+        Convert DAE-based variable address to relative output addresses.
+        """
+
+        bool_intersect = self.in1d(addr, v_code)
+        return np.where(bool_intersect)
diff --git a/andes/models/synchronous/genrou.py b/andes/models/synchronous/genrou.py
index 9368d8350..dda57383b 100644
--- a/andes/models/synchronous/genrou.py
+++ b/andes/models/synchronous/genrou.py
@@ -26,7 +26,7 @@ def __init__(self):
                            tex_name='x_q',
                            z=True,
                            )
-        self.xd2 = NumParam(default=0.204, info='d-axis sub-transient reactance',
+        self.xd2 = NumParam(default=0.3, info='d-axis sub-transient reactance',
                             tex_name=r"x''_d", z=True)
 
         self.xq1 = NumParam(default=0.5, info='q-axis transient reactance',
@@ -152,7 +152,7 @@ def __init__(self):
 
         # begin variables and equations
         self.psi2q = Algeb(tex_name=r"\psi_{aq}", info='q-axis air gap flux',
-                           v_str='psi2q0',
+                           v_str='u * psi2q0',
                            e_str='gq1*e1d + (1-gq1)*e2q - psi2q',
                            )
 
@@ -196,7 +196,7 @@ def __init__(self):
 
         self.e1d = State(info='d-axis transient voltage',
                          tex_name=r"e'_d",
-                         v_str='e1d0',
+                         v_str='u * e1d0',
                          e_str='-XaqI1q',
                          t_const=self.Tq10,
                          )
@@ -210,7 +210,7 @@ def __init__(self):
 
         self.e2q = State(info='q-axis sub-transient voltage',
                          tex_name=r"e''_q",
-                         v_str='e2q0',
+                         v_str='u * e2q0',
                          e_str='(-e2q + e1d + (xq1 - xl) * Iq)',
                          t_const=self.Tq20,
                          )
@@ -223,6 +223,13 @@ def __init__(self):
 class GENROU(GENROUData, GENBase, GENROUModel, Flux0):
     """
     Round rotor generator with quadratic saturation.
+
+    Notes
+    -----
+    Parameters:
+
+    - ``xd2`` and ``xq2`` must be equal to pass initialization.
+
     """
 
     def __init__(self, system, config):
diff --git a/andes/plot.py b/andes/plot.py
index 464bec3dd..0c5c9ddae 100644
--- a/andes/plot.py
+++ b/andes/plot.py
@@ -10,10 +10,85 @@
 import numpy as np
 
 from andes.core.var import BaseVar
-from andes.shared import plt, set_font, set_latex
+from andes.shared import find_executable, mpl, plt
+from andes.utils.paths import get_dot_andes_path  # NOQA
 
 logger = logging.getLogger(__name__)
-DPI = 80
+DPI = None
+
+
+def set_latex():
+    """
+    Enables LaTeX for matplotlib based on the `with_latex` option and `dvipng` availability.
+
+    Returns
+    -------
+    bool
+        True for LaTeX on, False for off
+    """
+
+    if find_executable('dvipng'):
+        mpl.rc('text', usetex=True)
+
+        no_warn_file = os.path.join(get_dot_andes_path(), '.no_warn_latex')
+        if not os.path.isfile(no_warn_file):
+            print('Using LaTeX for rendering. If an error occurs:')
+            print('a) If you are using `andes plot`, disable with option "-d",')
+            print('b) If you are using `plot()`, set "latex=False".')
+
+            try:
+                with open(os.path.join(get_dot_andes_path(), '.no_warn_latex'), 'w') as f:
+                    f.write('0')
+            except OSError:
+                pass
+
+        return True
+
+    return False
+
+
+def set_font(family='serif', size=12, style='normal', weight='normal'):
+    """
+    Sets the font for matplotlib.
+
+    Parameters
+    ----------
+    family : str
+        Font family.
+    size : int
+        Font size.
+    style : str
+        Font style.
+    weight : str
+        Font weight.
+    """
+
+    mpl.rc('font', family=family, size=size, style=style, weight=weight)
+
+
+def set_style(style='default'):
+    """
+    Set matplotlib style.
+
+    Parameters
+    ----------
+    style : str
+        `default`, `ieee` (require `scienceplots`), or other available styles
+        (see `matplotlib.pyplot.style.available`).
+    """
+
+    if style is None:
+        style = 'default'
+
+    if style == 'ieee':
+        try:
+            plt.style.use(['science', 'ieee'])
+        except (NameError, OSError):
+            logger.error("Please install `scienceplots` with `pip` to use the 'ieee' style.")
+
+    else:
+        plt.style.use(style)
+        set_font()
 
 
 class TDSData:
@@ -55,6 +130,12 @@ def __init__(self, full_name=None, mode='file', dae=None, path=None):
         self._process_names()
 
     def _process_names(self):
+        """
+        Helper function to process the file names.
+
+        Use "Untitled" if name is not set for case file.
+        """
+
         if self.full_name is None:
             logger.info("Input file name not detected. Using `Untitled`.")
             self.full_name = 'Untitled'
@@ -73,21 +154,32 @@ def _process_names(self):
         self._csv_file = os.path.join(self._path, self.file_name + '.csv')
 
     def load_dae(self):
-        """Load from DAE time series"""
+        """
+        Load from DAE time series.
+        """
+
         dae = self.dae
+        system = self.dae.system
+
         self.t = dae.ts.t
-        self.nvars = dae.n + dae.m + dae.o + 1
 
+        if system.Output.n == 0:
+            self.nvars = dae.n + dae.m + dae.o + 1
+        else:
+            self.nvars = len(system.Output.xidx) + len(system.Output.yidx) + dae.o + 1
+
+        self._uname = ['Time [s]'] + dae.x_name_output + dae.y_name_output + dae.z_name
+        self._fname = ['Time [s]'] + dae.x_tex_name_output + dae.y_tex_name_output + dae.z_tex_name
         self._idx = list(range(self.nvars))
-        self._uname = ['Time [s]'] + dae.x_name + dae.y_name + dae.z_name
-        self._fname = ['Time [s]'] + dae.x_tex_name + dae.y_tex_name + dae.z_tex_name
-        self._data = dae.ts.txyz
 
         if dae.system.files.lst is not None:
             self.full_name = dae.system.files.lst
         else:
             self.full_name = dae.system.files.case
 
+        if len(self.t) > 0:
+            self._data = dae.ts.txyz
+
     def load_lst(self):
         """
         Load the lst file into internal data structures `_idx`, `_fname`, `_uname`, and counts the number of
@@ -269,21 +361,44 @@ def _process_yidx(self, yidx, a):
 
         Indexing by ``a`` is considered.
         """
+
+        dae = self.dae
+        system = self.dae.system
+
         if isinstance(yidx, BaseVar):
             yidx = [yidx]
 
         if isinstance(yidx, list) and isinstance(yidx[0], BaseVar):
             all_yidx = np.array([], dtype=int)
+
             for item in yidx:
                 if item.n == 0:
-                    logger.warning("Variable <%s> contains no values, ignored.", item.name)
+                    logger.info("Parent model <%s> of variable <%s> contains no device, ignored.",
+                                item.owner.class_name, item.name)
                     continue
+
+                if system.Output.n > 0:
+                    output_addr = system.Output.to_output_addr(item.a, item.v_code)[0]
+                    if len(output_addr) == 0:
+                        logger.info("<%s.%s> contains no saved data, skipped.", item.owner.class_name, item.name)
+                        continue
+
+                    if len(output_addr) != len(item.a):
+                        logger.info("<%s.%s> is partially stored as set in <Output>. Showing all saved data.",
+                                    item.owner.class_name, item.name)
+
+                    nx = len(system.Output.xidx)
+                else:
+                    output_addr = item.a
+                    nx = dae.n
+
+                # states are offset by 1 for Time. Algebs are offset by 1 + nx
                 if item.v_code == 'y':
-                    offs = self.dae.n + 1
+                    offset = nx + 1
                 else:
-                    offs = 1
+                    offset = 1
 
-                new_yidx = item.a + offs
+                new_yidx = output_addr + offset
 
                 if a is not None:
                     new_yidx = np.take(new_yidx, a)
@@ -291,6 +406,9 @@ def _process_yidx(self, yidx, a):
 
             yidx = all_yidx
 
+        elif isinstance(yidx, int):
+            yidx = [yidx]
+
         # a list of integers will remain unchanged
 
         return yidx
@@ -508,7 +626,7 @@ def plot_data(self, xdata, ydata, *, xheader=None, yheader=None, xlabel=None, yl
                   hline1=None, hline2=None, vline1=None, hline=None, vline=None,
                   vline2=None, set_xlim=True, set_ylim=True, autoscale=False, figsize=None,
                   legend_bbox=None, legend_loc=None, legend_ncol=1,
-                  mask=True, color=None,
+                  mask=True, color=None, style='default',
                   **kwargs):
         """
         Plot lines for the supplied data and options.
@@ -549,7 +667,8 @@ def plot_data(self, xdata, ydata, *, xheader=None, yheader=None, xlabel=None, yl
         >>>                      ylabel='Ipcmd [pu]')
 
         """
-        set_font(family='serif', size=font_size)
+
+        set_style(style)
 
         if not isinstance(ydata, np.ndarray):
             raise TypeError("ydata must be a numpy array. Retrieve with get_values().")
diff --git a/andes/routines/eig.py b/andes/routines/eig.py
index 63d67e85b..6c1a26d27 100644
--- a/andes/routines/eig.py
+++ b/andes/routines/eig.py
@@ -3,24 +3,21 @@
 """
 
 import logging
-import scipy.io
-import numpy as np
-
 from math import ceil, pi
+
+import numpy as np
+import scipy.io
 from scipy.linalg import solve
 
 from andes.io.txt import dump_data
-from andes.utils.misc import elapsed
+from andes.plot import set_latex, set_style
 from andes.routines.base import BaseRoutine
+from andes.shared import div, matrix, plt, sparse, spdiag, spmatrix
+from andes.utils.misc import elapsed
 from andes.variables.report import report_info
 
-from andes.plot import set_latex
-
-from andes.shared import matrix, set_font, spmatrix, sparse, plt
-from andes.shared import div, spdiag
-
-
 logger = logging.getLogger(__name__)
+DPI = None
 
 
 class EIG(BaseRoutine):
@@ -334,8 +331,8 @@ def stats(self):
 
     def plot(self, mu=None, fig=None, ax=None,
              left=-6, right=0.5, ymin=-8, ymax=8, damping=0.05,
-             line_width=0.5, dpi=100, figsize=None, base_color='black',
-             show=True, latex=True):
+             line_width=0.5, dpi=DPI, figsize=None, base_color='black',
+             show=True, latex=True, style='default'):
         """
         Plot utility for eigenvalues in the S domain.
 
@@ -360,7 +357,7 @@ def plot(self, mu=None, fig=None, ax=None,
         line_width : float, optional
             default line width, by default 0.5
         dpi : int, optional
-            figure dpi, by default 100
+            figure dpi
         figsize : [type], optional
             default figure size, by default None
         base_color : str, optional
@@ -379,7 +376,7 @@ def plot(self, mu=None, fig=None, ax=None,
 
         """
 
-        set_font(family='Times New Roman', size=12)
+        set_style(style)
 
         if mu is None:
             mu = self.mu
diff --git a/andes/routines/tds.py b/andes/routines/tds.py
index 6d467acaf..6afe875dd 100644
--- a/andes/routines/tds.py
+++ b/andes/routines/tds.py
@@ -49,10 +49,14 @@ def __init__(self, system=None, config=None):
                                      ('qrt', 0),
                                      ('kqrt', 1.0),
                                      ('store_z', 0),
-                                     ('store_f', 0.0),
-                                     ('store_h', 0.0),
-                                     ('store_i', 0.0),
-                                     ('no_tqdm', 0.0),
+                                     ('store_f', 0),
+                                     ('store_h', 0),
+                                     ('store_i', 0),
+                                     ('limit_store', 0),
+                                     ('max_store', 900),
+                                     ('save_every', 1),
+                                     ('save_mode', 'auto'),
+                                     ('no_tqdm', 0),
                                      )))
         self.config.add_extra("_help",
                               method='DAE solution method',
@@ -77,6 +81,10 @@ def __init__(self, system=None, config=None):
                               store_f='store RHS of diff. equations',
                               store_h='store RHS of external diff. equations',
                               store_i='store RHS of external algeb. equations',
+                              limit_store='limit in-memory timeseries storage',
+                              max_store='maximum steps of data stored in memory before offloading',
+                              save_every='save results for one step every "save_every" steps',
+                              save_mode='automatically or manually save output data when done',
                               no_tqdm='disable tqdm progressbar and outputs',
                               )
         self.config.add_extra("_alt",
@@ -101,6 +109,10 @@ def __init__(self, system=None, config=None):
                               store_f=(0, 1),
                               store_h=(0, 1),
                               store_i=(0, 1),
+                              limit_store=(0, 1),
+                              max_store='positive integer',
+                              save_every='integer',
+                              save_mode=('auto', 'manual'),
                               no_tqdm=(0, 1),
                               )
 
@@ -186,6 +198,7 @@ def init(self):
         #   to the new array after extending `dae.y`.
         system.set_address(models=system.exist.pflow_tds)
         system.set_dae_names(models=system.exist.tds)
+        system.set_output_subidx(models=system.exist.pflow_tds)
 
         system.dae.clear_ts()
         system.store_sparse_pattern(models=system.exist.pflow_tds)
@@ -344,6 +357,10 @@ def run(self, no_summary=False, **kwargs):
         self.qrt_start = time.time()
         self.headroom = 0.0
 
+        # write variable list file at the beginning
+        if not system.files.no_output:
+            system.dae.write_lst(self.system.files.lst)
+
         t0, _ = elapsed()
 
         while (system.dae.t - self.h < self.config.tf) and (not self.busted):
@@ -362,17 +379,36 @@ def run(self, no_summary=False, **kwargs):
                 self.call_stats.append((system.dae.t.tolist(), self.niter, step_status))
 
             if step_status:
-                dae.store()
+                if config.save_every != 0:
+                    if config.save_every == 1:
+                        dae.store()
+                    else:
+                        if dae.kcount % config.save_every == 0:
+                            dae.store()
+
+                # offload if exceeds `max_store`
+                if self.config.limit_store and len(dae.ts._ys) >= self.config.max_store:
+
+                    # write to file if enabled
+                    if not system.files.no_output:
+                        self.save_output()
+                        logger.info("Offload data from memory to file for t=%.2f - %.2f sec",
+                                    dae.ts.t[0], dae.ts.t[-1])
+
+                    # clear storage in memory anyway
+                    dae.ts.reset()
 
                 self.streaming_step()
+
                 if self.check_criteria() is False:
-                    self.err_msg = 'Violated stability criteria.'
+                    self.err_msg = 'Violated stability criteria. To turn off, set [TDS].criteria = 0.'
                     self.busted = True
 
                 # check if the next step is critical time
                 self.do_switch()
                 self.calc_h()
                 dae.t += self.h
+                dae.kcount += 1
 
                 # show progress in percentage
                 perc = max(min((dae.t - config.t0) / (config.tf - config.t0) * 100, 100), 0)
@@ -420,15 +456,24 @@ def run(self, no_summary=False, **kwargs):
 
         t1, s1 = elapsed(t0)
         self.exec_time = t1 - t0
-        logger.info('Simulation completed in %s.', s1)
+        logger.info('Simulation to t=%.2f sec completed in %s.', config.tf, s1)
 
         if config.qrt:
             logger.debug('QRT headroom time: %.4g s.', self.headroom)
 
-        # need to unpack data in case of resumed simulations.
-        system.dae.ts.unpack()
-        if not system.files.no_output:
+        # in case of resumed simulations,
+        # manually unpack data to update arrays in `dae.ts`
+        # disable warning in case data has just been dumped
+        system.dae.ts.unpack(warn_empty=False)
+
+        if (not system.files.no_output) and (config.save_mode == 'auto'):
+            t0, _ = elapsed()
             self.save_output()
+            _, s1 = elapsed(t0)
+
+            np_file = self.system.files.npz
+            logger.info('Outputs to "%s" and "%s".', self.system.files.lst, np_file)
+            logger.info('Outputs written in %s.', s1)
 
         # end data streaming
         if system.config.dime_enabled:
@@ -668,21 +713,13 @@ def save_output(self, npz=True):
             True if files are written. False otherwise.
         """
 
-        t0, _ = elapsed()
-
-        self.system.dae.write_lst(self.system.files.lst)
-
         if npz is True:
-            np_file = self.system.files.npz
             self.system.dae.write_npz(self.system.files.npz)
         else:
-            np_file = self.system.files.npy
             self.system.dae.write_npy(self.system.files.npy)
 
-        _, s1 = elapsed(t0)
+        self.system.dae.ts.idx_ptr = len(self.system.dae.ts.t)
 
-        logger.info('Outputs to "%s" and "%s".', self.system.files.lst, np_file)
-        logger.info('Outputs written in %s.', s1)
         return True
 
     def do_switch(self):
diff --git a/andes/shared.py b/andes/shared.py
index 4d5a8d00a..8ae53eb06 100644
--- a/andes/shared.py
+++ b/andes/shared.py
@@ -11,15 +11,13 @@
 #     For example, NumPy must not be imported with LazyImport.
 
 import math
-import os
-import psutil
-
 
 import coloredlogs  # NOQA
 import numpy as np  # NOQA
+import psutil
 from kvxopt import printing  # NOQA
 from kvxopt import (div, klu, matrix, mul, sparse, spdiag, spmatrix,  # NOQA
-                    umfpack,)
+                    umfpack)
 from kvxopt.lapack import gesv  # NOQA
 
 from andes.utils.lazyimport import LazyImport
@@ -27,7 +25,6 @@
 printing.options['dformat'] = '%.1f'
 printing.options['width'] = -1
 
-from andes.utils.paths import get_dot_andes_path  # NOQA
 
 # --- SYSTEM INFO ---
 NCPUS_PHYSICAL = psutil.cpu_count(logical=False)
@@ -81,52 +78,3 @@
 
 # --- Shared functions ---
 find_executable = LazyImport('from distutils.spawn import find_executable')
-
-
-def set_latex():
-    """
-    Enables LaTeX for matplotlib based on the `with_latex` option and `dvipng` availability.
-
-    Returns
-    -------
-    bool
-        True for LaTeX on, False for off
-    """
-
-    if find_executable('dvipng'):
-        mpl.rc('text', usetex=True)
-
-        no_warn_file = os.path.join(get_dot_andes_path(), '.no_warn_latex')
-        if not os.path.isfile(no_warn_file):
-            print('Using LaTeX for rendering. If an error occurs:')
-            print('a) If you are using `andes plot`, disable with option "-d",')
-            print('b) If you are using `plot()`, set "latex=False".')
-
-            try:
-                with open(os.path.join(get_dot_andes_path(), '.no_warn_latex'), 'w') as f:
-                    f.write('0')
-            except OSError:
-                pass
-
-        return True
-
-    return False
-
-
-def set_font(family='serif', size=12, style='normal', weight='normal'):
-    """
-    Sets the font for matplotlib.
-
-    Parameters
-    ----------
-    family : str
-        Font family.
-    size : int
-        Font size.
-    style : str
-        Font style.
-    weight : str
-        Font weight.
-    """
-
-    mpl.rc('font', family=family, size=size, style=style, weight=weight)
diff --git a/andes/system.py b/andes/system.py
index 98808a1f7..df5e10213 100644
--- a/andes/system.py
+++ b/andes/system.py
@@ -2148,6 +2148,64 @@ def fix_address(self):
 
         return True
 
+    def set_output_subidx(self, models):
+        """
+        Process :py:class:`andes.models.misc.Output` data and store the
+        sub-indices into ``dae.xy``.
+
+        Parameters
+        ----------
+        models : OrderedDict
+            Models currently in use for the routine
+        """
+
+        export_vars = dict(x=list(), y=list())  # indices of export x and y
+
+        for model, var, dev in zip(self.Output.model.v,
+                                   self.Output.varname.v,
+                                   self.Output.dev.v):
+
+            # check validity of model name
+            if model not in models:
+                logger.info("Output model <%s> invalid or contains no device. Skipped.",
+                            model)
+                continue
+            mdl_instance = models[model]
+            mdl_all_vars = mdl_instance.cache.all_vars
+
+            # check validity of var name
+            if var is not None and (var not in mdl_all_vars):
+                logger.info("Output model <%s> contains no variable <%s>. Skipped.",
+                            model, var)
+                continue
+
+            # check validity of dev idx
+            if (dev is not None) and (dev not in mdl_instance.idx.v):
+                logger.info("Output model <%s> contains no device <%s>. Skipped.",
+                            model, dev)
+                continue
+
+            # TODO: dev-based indexing is not fully supported
+            # for multi-index variables, such as those in COI.
+
+            if var is None:
+                for item in mdl_all_vars.values():
+                    if dev is None:
+                        export_vars[item.v_code].extend(item.a)
+                    else:
+                        uid = mdl_instance.idx2uid(dev)
+                        export_vars[item.v_code].append(item.a[uid])
+            else:  # with variable name
+                item = mdl_all_vars[var]
+                if dev is None:
+                    export_vars[item.v_code].extend(item.a)
+                else:  # with exact index
+                    uid = mdl_instance.idx2uid(dev)
+                    export_vars[item.v_code].append(item.a[uid])
+
+        self.Output.xidx = sorted(np.unique(export_vars['x']))
+        self.Output.yidx = sorted(np.unique(export_vars['y']))
+
 
 def load_pycode_from_path(pycode_path):
     """
diff --git a/andes/variables/dae.py b/andes/variables/dae.py
index 14757e2dd..60ad3724b 100644
--- a/andes/variables/dae.py
+++ b/andes/variables/dae.py
@@ -33,16 +33,27 @@ def __init__(self, dae=None):
         self._hs = OrderedDict()
         self._is = OrderedDict()
 
-    def unpack_np(self):
+        self.idx_ptr = 0  # index pointer to the beginning of data that should be written
+
+    def unpack_np(self, attr, warn_empty=True):
         """
         Unpack dict data into numpy arrays.
         """
+
         n_steps = len(self._ys)
 
-        self.t = np.array(list(self._ys.keys()))
+        if attr is None or 't' in attr:
+            self.t = np.array(list(self._ys.keys()))
 
         def _dict2array(src, dest):
-            nx = len(self.__dict__[src][0]) if len(self.__dict__[src]) else 0
+            """
+            Helper function to convert data stord in a dict to an array.
+            """
+            if len(self.__dict__[src]):
+                nx = len(self.__dict__[src][self.t[0]])
+            else:
+                nx = 0
+
             self.__dict__[dest] = np.zeros((n_steps, nx))
 
             if len(self.__dict__[src]) > 0:
@@ -53,51 +64,72 @@ def _dict2array(src, dest):
                  ('_fs', 'f'), ('_hs', 'h'), ('_is', 'i'))
 
         for a, b in pairs:
-            _dict2array(a, b)
+            if attr is None or b in attr:
+                _dict2array(a, b)
 
-        self.xy = np.hstack((self.x, self.y))
-        self.txy = np.hstack((self.t.reshape((-1, 1)), self.xy))
-        self.txyz = np.hstack((self.t.reshape((-1, 1)), self.xy, self.z))
+        if attr is None or attr == 'xy':
+            self.xy = np.hstack((self.x, self.y))
+        if attr is None or attr == 'txy':
+            self.txy = np.hstack((self.t.reshape((-1, 1)), self.x, self.y))
+        if attr is None or attr == 'txyz':
+            self.txyz = np.hstack((self.t.reshape((-1, 1)), self.x, self.y, self.z))
 
         if n_steps == 0:
-            logger.warning("TimeSeries does not contain any time stamp.")
+            if warn_empty:
+                logger.warning("TimeSeries does not contain any time stamp.")
             return False
+
         return True
 
-    def unpack(self, df=False):
+    def unpack(self, df=False, attr=None, warn_empty=True):
         """
         Unpack dict-stored data into arrays and/or dataframes.
 
         Parameters
         ----------
         df : bool
-            True to construct DataFrames `self.df` and `self.df_z` (time-consuming).
+            True to construct DataFrames `self.df` and `self.df_z`
+            (time-consuming).
+        attr : str, optional
+            Attribute name to unpack. If None, unpack all.
 
         Returns
         -------
         True when done.
         """
 
-        self.unpack_np()
+        self.unpack_np(attr=attr, warn_empty=warn_empty)
         if df is True:
-            self.unpack_df()
+            self.unpack_df(attr=attr)
 
         return True
 
-    def unpack_df(self):
+    def unpack_df(self, attr):
         """
         Construct pandas dataframes.
         """
 
-        self.df_x = pd.DataFrame.from_dict(self._xs, orient='index',
-                                           columns=self.dae.x_name)
-        self.df_y = pd.DataFrame.from_dict(self._ys, orient='index',
-                                           columns=self.dae.y_name)
-        self.df_z = pd.DataFrame.from_dict(self._zs, orient='index',
-                                           columns=self.dae.z_name)
+        uxname = self.dae.x_name_output
+        uyname = self.dae.y_name_output
+        uzname = self.dae.z_name
+
+        if attr is None or 'x' in attr:
+            self.df_x = pd.DataFrame.from_dict(self._xs, orient='index',
+                                               columns=uxname)
+
+        if attr is None or 'y' in attr:
+            self.df_y = pd.DataFrame.from_dict(self._ys, orient='index',
+                                               columns=uyname)
+
+        if attr is None or 'z' in attr:
+            self.df_z = pd.DataFrame.from_dict(self._zs, orient='index',
+                                               columns=uzname)
 
-        self.df_xy = pd.concat((self.df_x, self.df_y), axis=1)
-        self.df_xyz = pd.concat((self.df_xy, self.df_z), axis=1)
+        if attr is None or attr == 'df_xy':
+            self.df_xy = pd.concat((self.df_x, self.df_y), axis=1)
+
+        if attr is None or attr == 'df_xyz':
+            self.df_xyz = pd.concat((self.df_x, self.df_y, self.df_z), axis=1)
 
         return True
 
@@ -160,6 +192,7 @@ def _access_array(self, array_name, indices=None):
 
         The function checks for empty arrays and shows warnings.
         """
+
         if np.count_nonzero(self.__dict__[array_name]) == 0:
             logger.error("TimeSeries matrix <%s> contains no element. Check if `[TDS] store_%s = 1`",
                          array_name, array_name)
@@ -173,18 +206,40 @@ def _access_array(self, array_name, indices=None):
 
     @property
     def df(self):
+        """
+        Short-hand for the xy dataframe.
+        """
+
         return self.df_xy
 
     def __getattr__(self, attr):
         if attr in super().__getattribute__('_public'):
             df = True if attr.startswith("df") else False
-            self.unpack(df=df)
+            self.unpack(df=df, attr=attr)
 
         return super().__getattribute__(attr)
 
     def __getstate__(self):
         return self.__dict__
 
+    def reset(self):
+        """
+        Reset the internal storage and erase all data.
+        """
+        self._xs = OrderedDict()
+        self._ys = OrderedDict()
+        self._zs = OrderedDict()
+        self._fs = OrderedDict()
+        self._hs = OrderedDict()
+        self._is = OrderedDict()
+
+        self.unpack_np(attr=None, warn_empty=False)
+        self.unpack_df(attr=None)
+
+        self.idx_ptr = 0
+
+        logger.debug("TimeSeries storage is cleared.")
+
 
 class DAE:
     r"""
@@ -243,6 +298,7 @@ def __init__(self, system):
         self.system = system
         self.t = np.array(0.0, dtype=float)
         self.ts = DAETimeSeries(self)
+        self.kcount = 0  # time step count
 
         self._array_and_counter = {
             'f': 'n',  # differential equation RHS
@@ -252,7 +308,7 @@ def __init__(self, system):
             'z': 'o',  # limiter flags
             'h': 'p',  # RHS of external states
             'i': 'q',  # RHS of external algebraic variables
-            }
+        }
 
         self.m, self.n, self.o, self.p, self.q = 0, 0, 0, 0, 0
 
@@ -277,7 +333,31 @@ def __init__(self, system):
 
         self.tpl = dict()  # sparsity templates with constants
 
+        self._write_append = False  # True if data should be appended when writing to output
+        self._lst_written = False
+
     def request_address(self, array_name: str, ndevice, nvar, collate=False):
+        """
+        Interface for requesting addresses for a model.
+
+        Parameters
+        ----------
+        array_name : str
+            array name in 'x' and 'y'
+        ndevice : int
+            number of devices
+        nvar : int
+            number of variables
+        collate : bool, optional
+            False if the same variable for different devices are contiguous.
+            True if variables for the same devices should collate. Note: setting
+            ``collate`` to True will degrade the performance.
+
+        Returns
+        -------
+        list
+            A list of arrays for each variable.
+        """
 
         out = []
         counter_name = self._array_and_counter[array_name]
@@ -307,13 +387,16 @@ def clear_arrays(self):
         """
         Reset equation and variable arrays to empty.
         """
+
         self.clear_fg()
         self.clear_xy()
         self.clear_z()
 
     def clear_fg(self):
-        """Resets equation arrays to empty
         """
+        Resets equation arrays to empty.
+        """
+
         self.f[:] = 0
         self.g[:] = 0
 
@@ -321,6 +404,7 @@ def clear_xy(self):
         """
         Reset variable arrays to empty.
         """
+
         self.x[:] = 0
         self.y[:] = 0
 
@@ -328,12 +412,14 @@ def clear_z(self):
         """
         Reset status arrays to empty
         """
+
         self.z[:] = 0
 
     def clear_ijv(self):
         """
         Clear stored triplets.
         """
+
         self.triplets.clear_ijv()
 
     def restore_sparse(self, names=None):
@@ -347,6 +433,7 @@ def restore_sparse(self, names=None):
         names : None or list
             List of Jacobian names to restore sparsity pattern
         """
+
         if names is None:
             names = jac_names
         elif isinstance(names, str):
@@ -362,6 +449,7 @@ def reset(self):
         """
         Reset array sizes to zero and clear all arrays.
         """
+
         self.set_t(0.0)
         self.m = 0
         self.n = 0
@@ -467,41 +555,47 @@ def _compare_pattern(self, name):
 
     def store(self):
         """
-        Store values and equations to in internal TimeSeries storage.
+        Store values for the current time step to the TimeSeries storage. Values
+        include variables, equation RHS and discrete states.
         """
+
+        system = self.system
+        tds = self.system.TDS
         ts = self.ts
         t = self.t.tolist()
-        tds = self.system.TDS
 
-        z_vals = self.system.get_z(self.system.exist.pflow_tds) if tds.config.store_z else None
-        f_vals = self.f if tds.config.store_f else None
-        h_vals = self.h if tds.config.store_h else None
-        i_vals = self.i if tds.config.store_i else None
-
-        ts._xs[t] = np.array(self.x)
-        ts._ys[t] = np.array(self.y)
+        if system.Output.n > 0:
+            # select variables based on `Output`
+            ts._xs[t] = self.x[system.Output.xidx]
+            ts._ys[t] = self.y[system.Output.yidx]
+        else:
+            ts._xs[t] = np.array(self.x)
+            ts._ys[t] = np.array(self.y)
 
-        if z_vals is not None:
+        if tds.config.store_z:
+            z_vals = system.get_z(system.exist.pflow_tds)
             ts._zs[t] = np.array(z_vals)
-        if f_vals is not None:
-            ts._fs[t] = np.array(f_vals)
-        if h_vals is not None:
-            ts._hs[t] = np.array(h_vals)
-        if i_vals is not None:
-            ts._is[t] = np.array(i_vals)
+
+        if tds.config.store_f:
+            ts._fs[t] = np.array(self.f)
+        if tds.config.store_h:
+            ts._hs[t] = np.array(self.h)
+        if tds.config.store_i:
+            ts._is[t] = np.array(self.i)
 
     def resize_arrays(self):
         """
         Resize arrays to the new sizes `m` and `n`, and `o`.
 
         If ``m > len(self.y)`` or ``n > len(self.x``, arrays will be extended.
-        Otherwise, new empty arrays will be sliced, starting from 0 to the given size.
+        Otherwise, new empty arrays will be sliced, starting from 0 to the given
+        size.
 
         Warnings
         --------
         This function should not be called directly. Instead, it is called in
-        ``System.set_address`` which re-points variables used in power flow
-        to the new array for dynamic analyses.
+        ``System.set_address`` which re-points variables used in power flow to
+        the new array for dynamic analyses.
         """
         self.x = self._extend_or_slice(self.x, self.n)
         self.y = self._extend_or_slice(self.y, self.m)
@@ -549,44 +643,131 @@ def alloc_or_extend_names(self):
 
     @property
     def xy(self):
-        """Return a concatenated array of [x, y]."""
+        """
+        Return a concatenated array of [x, y].
+        """
+
         return np.hstack((self.x, self.y))
 
     @property
     def xyz(self):
-        """Return a concatenated array of [x, y]."""
+        """
+        Return a concatenated array of [x, y].
+        """
+
         return np.hstack((self.x, self.y, self.z))
 
     @property
     def fg(self):
-        """Return a concatenated array of [f, g]."""
+        """
+        Return a concatenated array of [f, g].
+        """
+
         return np.hstack((self.f, self.g))
 
+    @property
+    def x_name_output(self):
+        """
+        Return a list of state var names selected by Output.
+        """
+        if self.system.Output.n == 0:
+            return self.x_name
+        else:
+            return [self.x_name[i] for i in self.system.Output.xidx]
+
+    @property
+    def y_name_output(self):
+        """
+        Return a list of algeb var names selected by Output.
+        """
+
+        if self.system.Output.n == 0:
+            return self.y_name
+        else:
+            return [self.y_name[i] for i in self.system.Output.yidx]
+
+    @property
+    def x_tex_name_output(self):
+        """
+        Return a list of state var LaTeX names selected by Output.
+        """
+
+        if self.system.Output.n == 0:
+            return self.x_tex_name
+        else:
+            return [self.x_tex_name[i] for i in self.system.Output.xidx]
+
+    @property
+    def y_tex_name_output(self):
+        """
+        Return a list of algeb var LaTeX names selected by Output.
+        """
+
+        if self.system.Output.n == 0:
+            return self.y_tex_name
+        else:
+            return [self.y_tex_name[i] for i in self.system.Output.yidx]
+
     @property
     def xy_name(self):
-        """Return a concatenated list of all variable names without format."""
+        """
+        Return a concatenated list of all variable names without format.
+        """
+
         return self.x_name + self.y_name
 
     @property
     def xyz_name(self):
-        """Return a concatenated list of all variable names without format."""
+        """
+        Return a concatenated list of all variable names without format.
+        """
+
         return self.x_name + self.y_name + self.z_name
 
     @property
     def xy_tex_name(self):
-        """Return a concatenated list of all variable names in LaTeX format."""
+        """
+        Return a concatenated list of all variable names in LaTeX format.
+        """
+
         return self.x_tex_name + self.y_tex_name
 
     @property
     def xyz_tex_name(self):
-        """Return a concatenated list of all variable names in LaTeX format."""
+        """
+        Return a concatenated list of all variable names in LaTeX format.
+        """
+
         return self.x_tex_name + self.y_tex_name + self.z_tex_name
 
     def get_name(self, arr):
+        """
+        Helper function for geting the list of variable names based on the
+        array name.
+
+        Parameters
+        ----------
+        arr : str
+            Array name in 'f', 'g', 'x', 'y', 'z'.
+        """
+
         mapping = {'f': 'x', 'g': 'y', 'x': 'x', 'y': 'y', 'z': 'z'}
         return self.__dict__[mapping[arr] + '_name']
 
     def print_array(self, name, values=None, tol=None):
+        """
+        Debug helper to print array values and names.
+
+        Parameters
+        ----------
+        name : str
+            array name in 'f', 'g', 'x', 'y'
+        values : array-like, optional
+            substitute array values to use
+        tol : float, optional
+            tolerance value to use. Values below `tol` will not be displayed
+        """
+
         if values is None:
             values = self.__dict__[name]
 
@@ -618,18 +799,28 @@ def write_lst(self, lst_path):
             succeed flag
         """
 
+        if self._lst_written is True:
+            return
+
+        system = self.system
+
         out = ''
         template = '{:>6g}, {:>25s}, {:>35s}\n'
 
         # header line
         out += template.format(0, 'Time [s]', 'Time [s]')
 
-        # output variable indices
-        idx = list(range(self.m + self.n + self.o))
+        if system.Output.n == 0:
+            # output variable indices
+            idx = list(range(self.m + self.n + self.o))
 
-        # variable names concatenated
-        uname = self.xyz_name
-        fname = self.xyz_tex_name
+            # variable names concatenated
+            uname = self.xyz_name
+            fname = self.xyz_tex_name
+        else:
+            idx = list(range(len(system.Output.xidx) + len(system.Output.yidx) + self.o))
+            uname = self.x_name_output + self.y_name_output + self.z_name
+            fname = self.x_tex_name_output + self.y_tex_name_output + self.z_tex_name
 
         for e, i in enumerate(idx):
             # `idx` in the lst file is always consecutive
@@ -637,18 +828,63 @@ def write_lst(self, lst_path):
 
         with open(lst_path, 'w') as f:
             f.write(out)
+
+        self._lst_written = True
+
         return True
 
     def write_npy(self, file_path):
         """
         Write TDS data into NumPy uncompressed format.
         """
+
         txyz_data = self.ts.txyz
         np.save(file_path, txyz_data)
 
     def write_npz(self, file_path):
         """
         Write TDS data into NumPy compressed format.
+
+        The function supports writing out all values at once or writing them out
+        incrementally.
         """
-        txyz_data = self.ts.txyz
-        np.savez_compressed(file_path, data=txyz_data)
+
+        tds = self.system.TDS
+        ts = self.ts
+
+        if not tds.config.limit_store:
+            # write the whole TimeSeries in one step
+            txyz_data = self.ts.txyz
+            np.savez_compressed(file_path, data=txyz_data)
+
+        else:
+            # create a new npz file and write for the first time
+            if self._write_append is False:
+                txyz_data = self.ts.txyz[ts.idx_ptr:, :]
+                np.savez_compressed(file_path, data=txyz_data)
+                self._write_append = True
+                ts.idx_ptr = len(self.ts.t)
+
+            # write and append to an existing npz file
+            else:
+                self.ts.unpack()
+                txyz_data = self.ts.txyz[ts.idx_ptr:, :]
+
+                # skip if no new data
+                if len(txyz_data) == 0:
+                    logger.debug("No new data to write to file. Skipped.")
+                    return
+
+                data = np.load(file_path)['data']
+
+                # in most cases, append new data to the existing
+                if len(data) > 0:
+                    data = np.vstack((data, txyz_data))
+                    logger.debug("Appended new data to output file.")
+
+                # in case the previous step stopped at tf=0
+                else:
+                    data = txyz_data
+
+                np.savez_compressed(file_path, data=data)
+                ts.idx_ptr = len(self.ts.t)
diff --git a/docs/source/release-notes.rst b/docs/source/release-notes.rst
index f47b84696..6839a7942 100644
--- a/docs/source/release-notes.rst
+++ b/docs/source/release-notes.rst
@@ -6,6 +6,53 @@ Release notes
 
 The APIs before v3.0.0 are in beta and may change without prior notice.
 
+v1.7 Notes
+==========
+
+v1.7.0 (2022-05-22)
+-------------------
+Allow incrementally offloading simulation data from memory to the output file:
+
+- ``[TDS].limit_store`` is a boolean value to enable the limit for in-memory
+  time-series storage. If set to ``1``, data will be offloaded to the ``npz``
+  file every ``[TDS].max_store`` steps. Offloaded data will then be erased from
+  memory.
+- If you need to interact with the time-series data in memory, you need to keep
+  ``[TDS].limit_store`` to ``0``.
+
+Allow specifying models, variables, and/or devices to output:
+
+- See :ref:`Output`. The `model` field is mandatory. Leaving `varname` or `dev`
+  blank indicates the selection of all applicable elements. For example,
+  specifying `model` and `varname` without `dev` means that the variable for all
+  devices will be exported.
+- Plot tool works with in-memory time-series data specified by Output.
+
+Simulation output control:
+
+- Allow controlling the save frequency for output data in ``[TDS].save_every``.
+  The default value is ``1``, which means that every step will be saved. Setting
+  it to ``4``, for example, will save data every four steps. This setting
+  applies to the in-memory storage and the output data file.
+- Setting ``save_every = 0`` will immediately discard all data after each
+  simulation step.
+- Added the option ``[TDS].save_mode`` to change the automatic simulation data
+  dumping to manual. Accepted values are ``auto`` and ``manual``. This option
+  shall only be adjusted to ``manual`` when one is manually stepping the
+  simulation and wants to avoid writing to the output file when the simulation
+  reaches ``TDS.config.tf``. One will need to call ``TDS.save_output()`` when
+  the full simulation concludes to avoid losing unsaved data.
+
+Other changes:
+
+- Fix the initialization of offline synchronous generators.
+- Allow styles to be set for plots using the argument ``style``. To generate
+  figures for IEEE publications, use ``style=ieee`` (require package
+  ``scienceplots``).
+- Moved the writing of the ``lst`` file to the first step of simulation.
+- ``andes misc -C`` will not remove ``_out.csv`` file as it is considered data
+  for post-processing just like exported figures.
+
 v1.6 Notes
 ==========
 


Test patch:
diff --git a/tests/test_output.py b/tests/test_output.py
new file mode 100644
index 000000000..e887e36ed
--- /dev/null
+++ b/tests/test_output.py
@@ -0,0 +1,57 @@
+"""
+Test output selection.
+"""
+
+import numpy as np
+import unittest
+
+import andes
+
+
+class TestOutput(unittest.TestCase):
+    """
+    Class for Output test.
+    """
+
+    def test_output_xyname(self):
+        """
+        Test x_name and y_name for Output
+        """
+
+        ss = andes.load(andes.get_case("5bus/pjm5bus.json"),
+                        no_output=True,
+                        setup=False,
+                        default_config=True)
+
+        ss.add("Output", {"model": "GENCLS", "varname": "omega"})
+        ss.add("Output", {"model": "GENCLS", "varname": "delta", "dev": 2})
+        ss.add("Output", {"model": "Bus"})
+
+        ss.setup()
+
+        ss.PFlow.run()
+        ss.TDS.config.tf = 0.1
+        ss.TDS.run()
+        ss.TDS.load_plotter()
+
+        nt = len(ss.dae.ts.t)
+        nx = 5
+        ny = 10
+
+        self.assertEqual(len(ss.dae.x_name_output), nx)
+        self.assertEqual(len(ss.dae.y_name_output), ny)
+
+        self.assertEqual(ss.dae.ts.x.shape[1], nx)
+        self.assertEqual(ss.dae.ts.y.shape[1], ny)
+
+        np.testing.assert_array_equal(ss.Output.xidx, [1, 4, 5, 6, 7])
+        np.testing.assert_array_equal(ss.Output.yidx, [0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
+
+        # test loaded data by plot
+        self.assertEqual(ss.TDS.plt._data.shape, (nt, nx + ny + 1))
+        np.testing.assert_array_equal(ss.TDS.plt._process_yidx(ss.GENCLS.omega, a=None),
+                                      [2, 3, 4, 5])
+        np.testing.assert_array_equal(ss.TDS.plt._process_yidx(ss.GENCLS.delta, a=None),
+                                      [1])
+        np.testing.assert_array_equal(ss.TDS.plt._process_yidx(ss.TG2.pout, a=None),
+                                      [])