diff --git a/andes/cases/GBnetwork/README.md b/andes/cases/GBnetwork/README.md index ac10babbd..c24e2574f 100644 --- a/andes/cases/GBnetwork/README.md +++ b/andes/cases/GBnetwork/README.md @@ -1,4 +1,4 @@ GBnetwork power flow case downloaded from -https://www.maths.ed.ac.uk/optenergy/NetworkData/fullGB/ + Dynamic data is randomly generated and does not represent the actual system. diff --git a/andes/cases/ieee14/README.md b/andes/cases/ieee14/README.md index ff1fba78b..859425af2 100644 --- a/andes/cases/ieee14/README.md +++ b/andes/cases/ieee14/README.md @@ -2,7 +2,7 @@ IEEE 14-Bus System Power flow data is created by J. Conto. -Data is available at https://drive.google.com/drive/folders/0B7uS9L2Woq_7YzYzcGhXT2VQYXc +Data is available at Dynamic data is created by H. Cui for ANDES. diff --git a/andes/cases/ieee14/ieee14_pvd1.xlsx b/andes/cases/ieee14/ieee14_pvd1.xlsx new file mode 100644 index 000000000..93e834b5e Binary files /dev/null and b/andes/cases/ieee14/ieee14_pvd1.xlsx differ diff --git a/andes/cases/ieee14/pert.py b/andes/cases/ieee14/pert.py new file mode 100644 index 000000000..151dd266c --- /dev/null +++ b/andes/cases/ieee14/pert.py @@ -0,0 +1,14 @@ +""" +An empty pert file. +""" + + +def pert(t, system): + """ + Perturbation function called at each step. + + The function is named "pert" and takes two arguments: + ``t`` for simulation time, and ``system`` for the system object. + + """ + pass diff --git a/andes/cases/kundur/kundur_vsc.xlsx b/andes/cases/kundur/kundur_vsc.xlsx new file mode 100644 index 000000000..bd4a351a9 Binary files /dev/null and b/andes/cases/kundur/kundur_vsc.xlsx differ diff --git a/andes/cases/nordic44/README.md b/andes/cases/nordic44/README.md index 97b1b1f2b..68d43081c 100644 --- a/andes/cases/nordic44/README.md +++ b/andes/cases/nordic44/README.md @@ -1,5 +1,5 @@ Nordic 44-Bus System -Source: https://github.com/ALSETLab/Nordic44-Nordpool/tree/master/nordic44/models +Source: Most dynamic models in this system have not been supported. diff --git a/andes/cases/wecc/wecc_full.xlsx b/andes/cases/wecc/wecc_full.xlsx new file mode 100644 index 000000000..afeade470 Binary files /dev/null and b/andes/cases/wecc/wecc_full.xlsx differ diff --git a/andes/cases/wecc/wecc_gencls.xlsx b/andes/cases/wecc/wecc_gencls.xlsx new file mode 100644 index 000000000..4a09240b7 Binary files /dev/null and b/andes/cases/wecc/wecc_gencls.xlsx differ diff --git a/andes/core/block.py b/andes/core/block.py index b2b5cf101..1fe436029 100644 --- a/andes/core/block.py +++ b/andes/core/block.py @@ -1467,6 +1467,8 @@ class GainLimiter(Block): │ │ _____/ └─────┘ lower + TODO: Add an extra gain block "R" for y. + Parameters ---------- u : str, BaseVar @@ -1474,7 +1476,7 @@ class GainLimiter(Block): """ - def __init__(self, u, K, upper, lower, no_upper=False, no_lower=False, + def __init__(self, u, K, lower, upper, no_lower=False, no_upper=False, name=None, tex_name=None, info=None): Block.__init__(self, name=name, tex_name=tex_name, info=info) self.u = dummify(u) @@ -1518,6 +1520,66 @@ def define(self): self.y.e_str += f' - {self.name}_y' +class LimiterGain(Block): + """ + Limiter followed by a gain. + + Exports the limited output `y`, unlimited output `x`, and HardLimiter `lim`. :: + + upper ┌─────┐ + /¯¯¯¯¯ │ │ + u -> / -> │ K │ -> y + _____/ │ │ + lower └─────┘ + + The intermediate variable before the gain is not saved. + + Parameters + ---------- + u : str, BaseVar + Input variable, or an equation string for constructing an anonymous variable + + """ + + def __init__(self, u, K, lower, upper, no_lower=False, no_upper=False, + name=None, tex_name=None, info=None): + Block.__init__(self, name=name, tex_name=tex_name, info=info) + self.u = u + self.K = dummify(K) + self.upper = dummify(upper) + self.lower = dummify(lower) + + if (no_upper and no_lower) is True: + raise ValueError("no_upper or no_lower cannot both be True") + + self.no_lower = no_lower + self.no_upper = no_upper + + self.lim = HardLimiter(u=self.u, lower=self.lower, upper=self.upper, + no_upper=no_upper, no_lower=no_lower, + tex_name='lim') + + self.y = Algeb(info='Gain output after limiter', tex_name='y', discrete=self.lim) + + self.vars = {'lim': self.lim, 'y': self.y} + + def define(self): + """ + TODO: write docstring + """ + self.y.e_str = f'{self.K.name} * {self.u.name} * {self.name}_lim_zi' + self.y.v_str = f'{self.K.name} * {self.u.name} * {self.name}_lim_zi' + + if not self.no_upper: + self.y.e_str += f' + {self.K.name} * {self.name}_lim_zu*{self.upper.name}' + self.y.v_str += f' + {self.K.name} * {self.name}_lim_zu*{self.upper.name}' + if not self.no_lower: + self.y.e_str += f' + {self.K.name} * {self.name}_lim_zl*{self.lower.name}' + self.y.v_str += f' + {self.K.name} * {self.name}_lim_zl*{self.lower.name}' + + self.y.e_str += f' - {self.name}_y' + + class Piecewise(Block): """ Piecewise block. Outputs an algebraic variable `y`. diff --git a/andes/core/discrete.py b/andes/core/discrete.py index 86d8cf592..33c0dbca6 100644 --- a/andes/core/discrete.py +++ b/andes/core/discrete.py @@ -203,8 +203,10 @@ class Limiter(Discrete): True to only use the upper limit no_upper : bool True to only use the lower limit - equal: bool + equal : bool True to include equal signs in comparison (>= or <=). + no_warn : bool + Disable initial limit warnings zu : 0 or 1 Default value for `zu` if not enabled zl : 0 or 1 @@ -224,7 +226,8 @@ class Limiter(Discrete): """ def __init__(self, u, lower, upper, enable=True, name=None, tex_name=None, info=None, - no_upper=False, no_lower=False, equal=True, zu=0.0, zl=0.0, zi=1.0): + no_upper=False, no_lower=False, equal=True, no_warn=False, + zu=0.0, zl=0.0, zi=1.0): Discrete.__init__(self, name=name, tex_name=tex_name, info=info) self.u = u self.lower = dummify(lower) @@ -233,6 +236,7 @@ def __init__(self, u, lower, upper, enable=True, name=None, tex_name=None, info= self.no_upper = no_upper self.no_lower = no_lower self.equal = equal + self.no_warn = no_warn self.zu = np.array([zu]) self.zl = np.array([zl]) diff --git a/andes/core/model.py b/andes/core/model.py index 147ebd6d7..f07de476f 100644 --- a/andes/core/model.py +++ b/andes/core/model.py @@ -14,7 +14,7 @@ import os import logging from collections import OrderedDict, defaultdict -from typing import Iterable, Sized +from typing import Iterable, Sized, Callable, Union from andes.core.common import ModelFlags, JacTriplet, Config from andes.core.discrete import Discrete @@ -504,7 +504,8 @@ def __init__(self, system, config): 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. - An addition variable `dae_t` for the current simulation time can be used if the model has flag `tds`. + The addition variable `dae_t` for the current simulation time can be used if the model has flag `tds`. + The additional variable `sys_f` is for system frequency (from ``system.config.freq``). 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)`. @@ -572,7 +573,9 @@ def __init__(self, system=None, config=None): self.services_ext = OrderedDict() # external services (to be retrieved) self.services_ops = OrderedDict() # operational services (for special usages) - self.tex_names = OrderedDict((('dae_t', 't_{dae}'),)) + self.tex_names = OrderedDict((('dae_t', 't_{dae}'), + ('sys_f', 'f_{sys}'), + )) # Model behavior flags self.flags = ModelFlags() @@ -877,8 +880,9 @@ def refresh_inputs(self): for key, val in self.config.as_dict(refresh=True).items(): self._input[key] = np.array(val) - # update`dae_t` + # update`dae_t` and `sys_f` self._input['dae_t'] = self.system.dae.t + self._input['sys_f'] = self.system.config.freq def refresh_inputs_arg(self): """ @@ -1621,22 +1625,27 @@ def numba_jitify(self, parallel=False, cache=False): if self.system.config.numba != 1: return - try: - import numba - except ImportError: - return - if self.flags.jited is True: return - self.calls.f = numba.jit(self.calls.f, parallel=parallel, cache=cache) - self.calls.g = numba.jit(self.calls.g, parallel=parallel, cache=cache) + self.calls.f = self._jitify_func_only(self.calls.f, parallel=parallel, cache=cache) + self.calls.g = self._jitify_func_only(self.calls.g, parallel=parallel, cache=cache) for jname in self.calls.j: - self.calls.j[jname] = numba.jit(self.calls.j[jname], parallel=parallel, cache=cache) + self.calls.j[jname] = self._jitify_func_only(self.calls.j[jname], + parallel=parallel, cache=cache) self.flags.jited = True + def _jitify_func_only(self, func: Union[Callable, None], parallel=False, cache=False): + try: + import numba + except ImportError: + return + + if func is not None: + return numba.jit(func, parallel=parallel, cache=cache) + class SymProcessor: """ @@ -1673,7 +1682,7 @@ def __init__(self, parent): self.parent = parent # symbols that are input to lambda functions - # including parameters, variables, services, configs and "dae_t" + # including parameters, variables, services, configs, and scalars (dae_t, sys_f) self.inputs_dict = OrderedDict() self.vars_dict = OrderedDict() self.iters_dict = OrderedDict() @@ -1798,6 +1807,7 @@ def generate_symbols(self): self.tex_names[Symbol(var)] = Symbol(self.parent.__dict__[var].tex_name) self.inputs_dict['dae_t'] = Symbol('dae_t') + self.inputs_dict['sys_f'] = Symbol('sys_f') # build ``non_vars_dict`` by removing ``vars_dict`` keys from a copy of ``inputs`` self.non_vars_dict = OrderedDict(self.inputs_dict) diff --git a/andes/core/service.py b/andes/core/service.py index 254011601..026426cda 100644 --- a/andes/core/service.py +++ b/andes/core/service.py @@ -173,6 +173,9 @@ def __init__(self, self.u = dummify(u) def check(self, **kwargs): + """ + Check status and set event flags. + """ if not np.all(self.v == self.u.v): self.owner.system.TDS.custom_event = True logger.debug(f"Event flag set at t={self.owner.system.dae.t:.6f} sec.") @@ -255,6 +258,9 @@ def __init__(self, self.n_ext = 0 # number of extended events def assign_memory(self, n): + """ + Assign memory for internal data. + """ VarService.assign_memory(self, n) self.t_final = np.zeros_like(self.v) self.v_event = np.zeros_like(self.v) @@ -265,6 +271,11 @@ def assign_memory(self, n): self.t_ext.v = np.ones_like(self.u.v) * self.t_ext.v def check(self, **kwargs): + """ + Check if an extended event is in place. + + Supplied as a ``v_numeric`` to ``VarService``. + """ dae_t = self.owner.system.dae.t if dae_t == 0.0: diff --git a/andes/io/__init__.py b/andes/io/__init__.py index 683c00b71..1a8c86cda 100644 --- a/andes/io/__init__.py +++ b/andes/io/__init__.py @@ -73,7 +73,7 @@ def guess(system): if testlines(fid): true_format = item files.input_format = true_format - logger.debug(f'Input format guessed as {true_format}.') + logger.debug('Input format guessed as %s.', true_format) break if not true_format: @@ -85,7 +85,7 @@ def guess(system): for key, val in input_formats.items(): if add_ext[1:] in val: files.add_format = key - logger.debug(f'Addfile format guessed as {key}.') + logger.debug('Addfile format guessed as %s.', key) break return true_format @@ -106,32 +106,32 @@ def parse(system): # exit when no input format is given if not system.files.input_format: if not guess(system): - logger.error('Input format is not specified and cannot be inferred.') + logger.error('Input format unknown for file "%s".', system.files.case) return False # try parsing the base case file - logger.info(f'Parsing input file "{system.files.case}"') + logger.info('Parsing input file "%s"...', system.files.case) input_format = system.files.input_format parser = importlib.import_module('.' + input_format, __name__) if not parser.read(system, system.files.case): - logger.error(f'Error parsing case file {system.files.fullname} with {input_format} format parser.') + logger.error('Error parsing file "%s" with <%s> parser.', system.files.fullname, input_format) return False _, s = elapsed(t) - logger.info(f'Input file parsed in {s}.') + logger.info('Input file parsed in %s.', s) # Try parsing the addfile t, _ = elapsed() if system.files.addfile: - logger.info(f'Parsing additional file "{system.files.addfile}"') + logger.info('Parsing additional file "%s"...', system.files.addfile) add_format = system.files.add_format add_parser = importlib.import_module('.' + add_format, __name__) if not add_parser.read_add(system, system.files.addfile): - logger.error(f'Error parsing addfile {system.files.addfile} with {input_format} parser.') + logger.error('Error parsing addfile "%s" with %s parser.', system.files.addfile, input_format) return False _, s = elapsed(t) - logger.info(f'Addfile parsed in {s}.') + logger.info('Addfile parsed in %s.', s) return True @@ -175,7 +175,7 @@ def dump(system, output_format, full_path=None, overwrite=False, **kwargs): ret = writer.write(system, system.files.dump, overwrite=overwrite, **kwargs) _, s = elapsed(t) if ret: - logger.info(f'Format conversion completed in {s}.') + logger.info('Format conversion completed in %s.', s) return True else: logger.error('Format conversion failed.') diff --git a/andes/io/psse.py b/andes/io/psse.py index 6181e6628..d32abb215 100644 --- a/andes/io/psse.py +++ b/andes/io/psse.py @@ -153,7 +153,7 @@ def _read_dyr_dict(file): # concatenate multi-line device data input_concat_dict = defaultdict(list) multi_line = list() - for i, line in enumerate(input_list): + for line in input_list: if line == '': continue if '/' not in line: diff --git a/andes/models/__init__.py b/andes/models/__init__.py index 6becf2429..d5b86546e 100644 --- a/andes/models/__init__.py +++ b/andes/models/__init__.py @@ -1,4 +1,4 @@ -from collections import OrderedDict # NOQA +from collections import OrderedDict # Notes: @@ -12,7 +12,7 @@ ('pv', ['PV', 'Slack']), ('shunt', ['Shunt']), ('line', ['Line']), - ('area', ['Area', 'ACE']), + ('area', ['Area', 'ACE', 'ACEc']), ('synchronous', ['GENCLS', 'GENROU']), ('governor', ['TG2', 'TGOV1', 'TGOV1DB', 'IEEEG1']), ('exciter', ['EXDC2', 'IEEEX1', 'ESDC2A', 'EXST1', 'ESST3A', 'SEXS']), @@ -22,6 +22,7 @@ ('coi', ['COI', ]), ('dcbase', ['Node', 'Ground', 'R', 'L', 'C', 'RCp', 'RCs', 'RLs', 'RLCs', 'RLCp']), ('vsc', ['VSCShunt']), - ('renewable', ['REGCA1', 'REECA1', 'REPCA1', 'WTDTA1', 'WTDS', 'WTARA1', 'WTPTA1', 'WTTQA1', 'PVD1']), + ('renewable', ['REGCA1', 'REECA1', 'REPCA1', 'WTDTA1', 'WTDS', 'WTARA1', 'WTPTA1', 'WTTQA1']), + ('distributed', ['PVD1']), ('experimental', ['PI2', 'TestDB1', 'TestPI', 'TestLagAWFreeze']), ]) diff --git a/andes/models/area.py b/andes/models/area.py index 524c9f9ec..bd6d05d6c 100644 --- a/andes/models/area.py +++ b/andes/models/area.py @@ -65,14 +65,17 @@ def __init__(self): self.bias = NumParam(default=1.0, info='bias parameter', tex_name=r'\beta', unit='MW/0.1Hz', power=True) - self.busf = IdxParam(info='Optional BusFreq idx', model='BusFreq', + self.busf = IdxParam(info='Optional BusFreq device idx', model='BusFreq', default=None) -class ACE(ACEData, Model): +class ACEc(ACEData, Model): """ Area Control Error model. + Continuous frequency sampling. + System base frequency from ``system.config.freq`` is used. + Note: area idx is automatically retrieved from `bus`. """ @@ -84,13 +87,10 @@ def __init__(self, system, config): self.group = 'Calculation' self.config.add(OrderedDict([('freq_model', 'BusFreq'), - ('interval', 4.0), - ('offset', 0.0), ])) - self.config.add_extra('_help', {'freq_model': 'default freq. measurement model', - 'interval': 'sampling time interval', - 'offset': 'sampling time offset'}) - + self.config.add_extra('_help', + {'freq_model': 'default freq. measurement model', + }) self.config.add_extra('_alt', {'freq_model': ('BusFreq',)}) self.area = ExtParam(model='Bus', src='area', indexer=self.bus, export=False) @@ -98,9 +98,42 @@ def __init__(self, system, config): self.busf.model = self.config.freq_model self.busfreq = DeviceFinder(self.busf, link=self.bus, idx_name='bus') - self.f = ExtAlgeb(model='FreqMeasurement', src='f', indexer=self.busfreq, - export=False, info='Bus frequency', + self.f = ExtAlgeb(model='FreqMeasurement', + src='f', + indexer=self.busfreq, + export=False, + info='Bus frequency', + unit='p.u. (Hz)' ) + self.ace = Algeb(info='area control error', + unit='p.u. (MW)', + tex_name='ace', + e_str='10 * bias * sys_f * (f - 1) - ace', + ) + + +class ACE(ACEc): + """ + Area Control Error model. + + Discrete frequency sampling. + System base frequency from ``system.config.freq`` is used. + + Frequency sampling period (in seconds) can be specified in + ``ACE.config.interval``. The sampling start time (in seconds) + can be specified in ``ACE.config.offset``. + + Note: area idx is automatically retrieved from `bus`. + """ + + def __init__(self, system, config): + ACEc.__init__(self, system, config) + + self.config.add(OrderedDict([('interval', 4.0), + ('offset', 0.0), + ])) + self.config.add_extra('_help', {'interval': 'sampling time interval', + 'offset': 'sampling time offset'}) self.fs = Sampling(self.f, interval=self.config.interval, @@ -109,7 +142,4 @@ def __init__(self, system, config): info='Sampled freq.', ) - self.ace = Algeb(info='area control error', unit='MW (p.u.)', - tex_name='ace', - e_str='10 * bias * (fs_v - 1) - ace', - ) + self.ace.e_str = '10 * bias * sys_f * (fs_v - 1) - ace' diff --git a/andes/models/dcbase.py b/andes/models/dcbase.py index 6e2402df6..c432daf28 100644 --- a/andes/models/dcbase.py +++ b/andes/models/dcbase.py @@ -272,7 +272,8 @@ def __init__(self, system, config): info='Inductance current', unit='p.u.', v_str='0', - e_str='-u * (v1 - v2) / L', + e_str='-u * (v1 - v2)', + t_const=self.L, ) self.v1.e_str = '-IL' self.v2.e_str = '+IL' @@ -297,7 +298,8 @@ def __init__(self, system, config): info='Capacitor current', unit='p.u.', v_str='0', - e_str='-u * Idc / C', + e_str='-u * Idc', + t_const=self.C ) self.Idc = Algeb(tex_name='I_{dc}', info='Current from node 2 to 1', @@ -334,8 +336,9 @@ def __init__(self, system, config): self.IL = State(tex_name='I_L', info='Inductance current', unit='p.u.', - e_str='u * (v1 - v2 - R * IL) / L', + e_str='u * (v1 - v2 - R * IL)', v_str='(v1 - v2) / R', + t_const=self.L, ) self.Idc = Algeb(tex_name='I_{dc}', info='Current from node 2 to 1', @@ -370,8 +373,9 @@ def __init__(self, system, config): self.vC = State(tex_name='v_C', info='Capacitor current', unit='p.u.', - e_str='-u * (Idc - vC/R) / C', + e_str='-u * (Idc - vC/R)', v_str='v1 - v2', + t_const=self.C, ) self.Idc = Algeb(tex_name='I_{dc}', info='Current from node 2 to 1', @@ -416,13 +420,15 @@ def __init__(self, system, config): info='Inductance current', unit='p.u.', v_str='0', - e_str='u * vC / L', + e_str='u * vC', + t_const=self.L, ) self.vC = State(tex_name='v_C', info='Capacitor current', unit='p.u.', - e_str='-u * (Idc - vC/R - IL) / C', + e_str='-u * (Idc - vC/R - IL)', v_str='v1 - v2', + t_const=self.C, ) self.Idc = Algeb(tex_name='I_{dc}', info='Current from node 2 to 1', @@ -459,8 +465,9 @@ def __init__(self, system, config): self.vC = State(tex_name='v_C', info='Capacitor current', unit='p.u.', - e_str='-u * Idc / C', + e_str='-u * Idc', v_str='v1 - v2', + t_const=self.C, ) self.Idc = Algeb(tex_name='I_{dc}', info='Current from node 2 to 1', @@ -504,14 +511,16 @@ def __init__(self, system, config): self.IL = State(tex_name='I_L', info='Inductance current', unit='p.u.', - e_str='u * (v1 - v2 - R * IL - vC) / L', + e_str='u * (v1 - v2 - R * IL - vC)', v_str='0', + t_const=self.L, ) self.vC = State(tex_name='v_C', info='Capacitor current', unit='p.u.', - e_str='u * IL / C', + e_str='u * IL', v_str='v1 - v2', + t_const=self.C, ) self.Idc = Algeb(tex_name='I_{dc}', info='Current from node 2 to 1', diff --git a/andes/models/distributed.py b/andes/models/distributed.py new file mode 100644 index 000000000..b7e38d672 --- /dev/null +++ b/andes/models/distributed.py @@ -0,0 +1,428 @@ +""" +Distributed energy resource models. +""" + +from andes.core.model import Model, ModelData +from andes.core.param import NumParam, IdxParam +from andes.core.block import Lag, DeadBand1, LimiterGain +from andes.core.var import ExtAlgeb, Algeb + +from andes.core.service import ConstService, ExtService, VarService +from andes.core.service import DataSelect, DeviceFinder +from andes.core.discrete import Switcher, Limiter + + +class PVD1Data(ModelData): + """ + Data for distributed PV. + """ + def __init__(self): + ModelData.__init__(self) + + self.bus = IdxParam(model='Bus', + info="interface bus id", + mandatory=True, + ) + + self.gen = IdxParam(info="static generator index", + mandatory=True, + ) + + self.Sn = NumParam(default=100.0, tex_name='S_n', + info='Model MVA base', + unit='MVA', + ) + + self.fn = NumParam(default=60.0, tex_name='f_n', + info='nominal frequency', + unit='Hz', + ) + + self.busf = IdxParam(info='Optional BusFreq idx', + model='BusFreq', + default=None, + ) + + self.xc = NumParam(default=0.0, tex_name='x_c', + info='coupling reactance', + unit='p.u.', + z=True, + ) + + self.pqflag = NumParam(info='P/Q priority for I limit; 0-Q priority, 1-P priority', + mandatory=True, + unit='bool', + ) + + # --- parameters found from ESIG.energy --- + self.igreg = IdxParam(model='Bus', + info='Remote bus idx for droop response, None for local', + ) + + self.qmx = NumParam(default=0.33, tex_name='q_{mx}', + info='Max. reactive power command', + power=True, + ) + + self.qmn = NumParam(default=-0.33, tex_name='q_{mn}', + info='Min. reactive power command', + power=True, + ) + + self.v0 = NumParam(default=0.8, tex_name='v_0', + info='Lower limit of deadband for Vdroop response', + unit='pu', non_zero=True, + ) + self.v1 = NumParam(default=1.1, tex_name='v_1', + info='Upper limit of deadband for Vdroop response', + unit='pu', non_zero=True, + ) + + self.dqdv = NumParam(default=-1.0, tex_name='dq/dv', + info='Q-V droop characteristics (negative)', + power=True, non_zero=True + ) + + self.fdbd = NumParam(default=-0.01, tex_name='f_{dbd}', + info='frequency deviation deadband', + ) + + self.ddn = NumParam(default=0.0, tex_name='D_{dn}', + info='Gain after f deadband', + ) + + self.ialim = NumParam(default=1.3, tex_name='I_{alim}', + info='Apparent power limit', + current=True, + ) + + self.vt0 = NumParam(default=0.88, tex_name='V_{t0}', + info='Voltage tripping response curve point 0', + ) + + self.vt1 = NumParam(default=0.90, tex_name='V_{t1}', + info='Voltage tripping response curve point 1', + ) + + self.vt2 = NumParam(default=1.1, tex_name='V_{t2}', + info='Voltage tripping response curve point 2', + ) + + self.vt3 = NumParam(default=1.2, tex_name='V_{t3}', + info='Voltage tripping response curve point 3', + ) + + self.vrflag = NumParam(default=0.0, tex_name='z_{VR}', + info='Voltage tripping is latching (0) or partially self-resetting (0-1)', + ) + + self.ft0 = NumParam(default=59.5, tex_name='f_{t0}', + info='Frequency tripping response curve point 0', + ) + + self.ft1 = NumParam(default=59.7, tex_name='f_{t1}', + info='Frequency tripping response curve point 1', + ) + + self.ft2 = NumParam(default=60.3, tex_name='f_{t2}', + info='Frequency tripping response curve point 2', + ) + + self.ft3 = NumParam(default=60.5, tex_name='f_{t3}', + info='Frequency tripping response curve point 3', + ) + + self.frflag = NumParam(default=0.0, tex_name='z_{FR}', + info='Frequency tripping is latching (0) or partially self-resetting (0-1)', + ) + + self.tip = NumParam(default=0.02, tex_name='T_{ip}', + info='Inverter active current lag time constant', + unit='s', + ) + + self.tiq = NumParam(default=0.02, tex_name='T_{iq}', + info='Inverter reactive current lag time constant', + unit='s', + ) + + self.gammap = NumParam(default=1.0, tex_name=r'\gamma_p', + info='Ratio of P from PVD1 w.r.t to that from PV generator', + vrange='(0, 1]', + ) + + self.gammaq = NumParam(default=1.0, tex_name=r'\gamma_q', + info='Ratio of Q from PVD1 w.r.t to that from PV generator', + vrange='(0, 1]', + ) + + +class PVD1Model(Model): + """ + Model implementation of PVD1. + """ + def __init__(self, system, config): + Model.__init__(self, system, config) + self.flags.tds = True + self.group = 'DG' + + self.SWPQ = Switcher(u=self.pqflag, options=(0, 1), tex_name='SW_{PQ}', cache=True) + + self.buss = DataSelect(self.igreg, self.bus, + info='selected bus (bus or igreg)', + ) + + self.busfreq = DeviceFinder(self.busf, link=self.buss, idx_name='bus') + + # --- initial values from power flow --- + # v : bus voltage magnitude + # a : bus voltage angle + # p0s : active power from connected static PV generator + # q0s : reactive power from connected static PV generator + + self.v = ExtAlgeb(model='Bus', src='v', indexer=self.buss, tex_name='V', + info='bus (or igreg) terminal voltage', + unit='p.u.', + e_str='-Iqout_y * v * u', + ) + + self.a = ExtAlgeb(model='Bus', src='a', indexer=self.buss, tex_name=r'\theta', + info='bus (or igreg) phase angle', + unit='rad.', + e_str='-Ipout_y * v * u', + ) + + self.p0s = ExtService(model='StaticGen', + src='p', + indexer=self.gen, + tex_name='P_{0s}', + info='Initial P from static gen', + ) + self.q0s = ExtService(model='StaticGen', + src='q', + indexer=self.gen, + tex_name='Q_{0s}', + info='Initial Q from static gen', + ) + # --- calculate the initial P and Q for this distributed device --- + self.p0 = ConstService(v_str='gammap * p0s', tex_name='P_0', + info='Initial P for the PVD1 device', + ) + self.q0 = ConstService(v_str='gammaq * q0s', tex_name='Q_0', + info='Initial Q for the PVD1 device', + ) + + # frequency measurement variable `f` + self.f = ExtAlgeb(model='FreqMeasurement', src='f', indexer=self.busfreq, export=False, + info='Bus frequency', unit='p.u.', + ) + + self.fHz = Algeb(info='frequency in Hz', + v_str='fn * f', e_str='fn * f - fHz', + unit='Hz', + tex_name='f_{Hz}', + ) + + # --- frequency branch --- + self.FL1 = Limiter(u=self.fHz, lower=self.ft0, upper=self.ft1, + info='Under frequency comparer', no_warn=True, + ) + self.FL2 = Limiter(u=self.fHz, lower=self.ft2, upper=self.ft3, + info='Over frequency comparer', no_warn=True, + ) + + self.Kft01 = ConstService(v_str='1/(ft1 - ft0)', tex_name='K_{ft01}') + + self.Ffl = Algeb(info='Coeff. for under frequency', + v_str='FL1_zi * Kft01 * (fHz - ft0) + FL1_zu', + e_str='FL1_zi * Kft01 * (fHz - ft0) + FL1_zu - Ffl', + tex_name='F_{fl}', + discrete=self.FL1, + ) + + self.Kft23 = ConstService(v_str='1/(ft3 - ft2)', tex_name='K_{ft23}') + + self.Ffh = Algeb(info='Coeff. for over frequency', + v_str='FL2_zl + FL2_zi * (1 + Kft23 * (ft2 - fHz))', + e_str='FL2_zl + FL2_zi * (1 + Kft23 * (ft2 - fHz)) - Ffh', + tex_name='F_{fh}', + discrete=self.FL2, + ) + + self.Fdev = Algeb(info='Frequency deviation', + v_str='fn - fHz', e_str='fn - fHz - Fdev', + unit='Hz', tex_name='f_{dev}', + ) + + self.Dfdbd = ConstService(v_str='fdbd * ddn', info='Deadband lower limit after gain', + tex_name='D_{fdbd}', + ) + self.DB = DeadBand1(u=self.Fdev, center=0.0, lower=self.Dfdbd, upper=0.0, + info='frequency deviation deadband with gain', + ) # outputs `Pdrp` + self.DB.db.no_warn = True + + # --- Voltage flags --- + self.VL1 = Limiter(u=self.v, lower=self.vt0, upper=self.vt1, + info='Under voltage comparer', no_warn=True, + ) + self.VL2 = Limiter(u=self.v, lower=self.vt2, upper=self.vt3, + info='Over voltage comparer', no_warn=True, + ) + + self.Kvt01 = ConstService(v_str='1/(vt1 - vt0)', tex_name='K_{vt01}') + + self.Fvl = Algeb(info='Coeff. for under voltage', + v_str='VL1_zi * Kvt01 * (v - vt0) + VL1_zu', + e_str='VL1_zi * Kvt01 * (v - vt0) + VL1_zu - Fvl', + tex_name='F_{vl}', + discrete=self.VL1, + ) + + self.Kvt23 = ConstService(v_str='1/(vt3 - vt2)', tex_name='K_{vt23}') + + self.Fvh = Algeb(info='Coeff. for over voltage', + v_str='VL2_zl + VL2_zi * (1 + Kvt23 * (vt2 - v))', + e_str='VL2_zl + VL2_zi * (1 + Kvt23 * (vt2 - v)) - Fvh', + tex_name='F_{vh}', + discrete=self.VL2, + ) + # --- sensed voltage with lower limit of 0.01 --- + + self.VLo = Limiter(u=self.v, lower=0.01, upper=999, no_upper=True, + info='Voltage lower limit (0.01) flag', + ) + + self.vp = Algeb(tex_name='V_p', + info='Sensed positive voltage', + v_str='v * VLo_zi + 0.01 * VLo_zl', + e_str='v * VLo_zi + 0.01 * VLo_zl - vp', + ) + + self.Pext = Algeb(tex_name='P_{ext}', + info='External power signal', + v_str='0', + e_str='0 - Pext' + ) + + self.Psum = Algeb(tex_name='P_{tot}', + info='Sum of P signals', + v_str='Pext + p0 + DB_y', + e_str='Pext + p0 + DB_y - Psum', + ) # `p0` is the initial `Pref`, and `DB_y` is `Pdrp` (f droop) + + self.Vcomp = VarService(v_str='abs(v*exp(1j*a) + (1j * xc) * (Ipout_y + 1j * Iqout_y))', + info='Voltage before Xc compensation', + tex_name='V_{comp}' + ) + + self.Vqu = ConstService(v_str='v1 - (q0 - qmn) / dqdv', + info='Upper voltage bound => qmx', + tex_name='V_{qu}', + ) + + self.Vql = ConstService(v_str='v0 + (qmx - q0) / dqdv', + info='Lower voltage bound => qmn', + tex_name='V_{ql}', + ) + + self.VQ1 = Limiter(u=self.Vcomp, lower=self.Vql, upper=self.v0, + info='Under voltage comparer for Q droop', + no_warn=True, + ) + + self.VQ2 = Limiter(u=self.Vcomp, lower=self.v1, upper=self.Vqu, + info='Over voltage comparer for Q droop', + no_warn=True, + ) + + Qsum = 'VQ1_zl * qmx + VQ2_zu * qmn + ' \ + 'VQ1_zi * (qmx + dqdv *(Vqu - Vcomp)) + ' \ + 'VQ2_zi * (q0 + dqdv * (v1 - Vcomp)) + ' \ + 'q0' + + self.Qsum = Algeb(info='Total Q (droop + initial)', + v_str=Qsum, + e_str=f'{Qsum} - Qsum', + tex_name='Q_{sum}', + discrete=(self.VQ1, self.VQ2), + ) + + self.Ipul = Algeb(info='Ipcmd before hard limit', + v_str='Psum / vp', + e_str='Psum / vp - Ipul', + tex_name='I_{p,ul}', + ) + + self.Iqul = Algeb(info='Iqcmd before hard limit', + v_str='Qsum / vp', + e_str='Qsum / vp - Iqul', + tex_name='I_{q,ul}', + ) + + # --- Ipmax, Iqmax and Iqmin --- + Ipmaxsq = "(Piecewise((0, Le(ialim**2 - Iqcmd_y**2, 0)), ((ialim**2 - Iqcmd_y ** 2), True)))" + Ipmaxsq0 = "(Piecewise((0, Le(ialim**2 - (q0 / v)**2, 0)), ((ialim**2 - (q0 / v) ** 2), True)))" + self.Ipmaxsq = VarService(v_str=Ipmaxsq, tex_name='I_{pmax}^2') + self.Ipmaxsq0 = ConstService(v_str=Ipmaxsq0, tex_name='I_{pmax0}^2') + + self.Ipmax = Algeb(v_str='SWPQ_s1 * ialim + SWPQ_s0 * sqrt(Ipmaxsq0)', + e_str='SWPQ_s1 * ialim + SWPQ_s0 * sqrt(Ipmaxsq) - Ipmax', + tex_name='I_{pmax}', + ) + + Iqmaxsq = "(Piecewise((0, Le(ialim**2 - Ipcmd_y**2, 0)), ((ialim**2 - Ipcmd_y ** 2), True)))" + Iqmaxsq0 = "(Piecewise((0, Le(ialim**2 - (p0 / v)**2, 0)), ((ialim**2 - (p0 / v) ** 2), True)))" + self.Iqmaxsq = VarService(v_str=Iqmaxsq, tex_name='I_{qmax}^2') + self.Iqmaxsq0 = ConstService(v_str=Iqmaxsq0, tex_name='I_{qmax0}^2') + + self.Iqmax = Algeb(v_str='SWPQ_s0 * ialim + SWPQ_s1 * sqrt(Iqmaxsq0)', + e_str='SWPQ_s0 * ialim + SWPQ_s1 * sqrt(Iqmaxsq) - Iqmax', + tex_name='I_{qmax}', + ) + + self.Iqmin = VarService(v_str='-Iqmax', tex_name='I_{qmin}') + + # --- Ipcmd, Iqcmd --- + + self.Ipcmd = LimiterGain(u=self.Ipul, K='Fvl * Fvh * Ffl * Ffh', + lower=0.0, upper=self.Ipmax, + info='Ip with limiter and coeff.', + tex_name='I^{pcmd}', + ) + + self.Iqcmd = LimiterGain(u=self.Iqul, K='Fvl * Fvh * Ffl * Ffh', + lower=self.Iqmin, upper=self.Iqmax, + info='Iq with limiter and coeff.', + tex_name='I^{qcmd}', + ) + + self.Ipout = Lag(u=self.Ipcmd_y, T=self.tip, K=1.0, + info='Output Ip filter', + ) + + self.Iqout = Lag(u=self.Iqcmd_y, T=self.tiq, K=1.0, + info='Output Iq filter', + ) + + def v_numeric(self, **kwargs): + """ + Disable the corresponding `StaticGen`s. + """ + self.system.groups['StaticGen'].set(src='u', idx=self.gen.v, attr='v', value=0) + + +class PVD1(PVD1Data, PVD1Model): + """ + WECC Distributed PV model. + + Power rating specified in `Sn`. + Frequency and voltage recovery latching has not been implemented. + + Reference: + [1] ESIG, WECC Distributed and Small PV Plants Generic Model (PVD1), [Online], + Available: https://www.esig.energy/wiki-main-page/wecc-distributed-and-small-pv-plants-generic-model-pvd1/ + """ + def __init__(self, system, config): + PVD1Data.__init__(self) + PVD1Model.__init__(self, system, config) diff --git a/andes/models/renewable.py b/andes/models/renewable.py index d37a796a5..b073f74fc 100644 --- a/andes/models/renewable.py +++ b/andes/models/renewable.py @@ -2102,173 +2102,3 @@ class WTTQA1(WTTQA1Data, WTTQA1Model): def __init__(self, config, system): WTTQA1Data.__init__(self) WTTQA1Model.__init__(self, system, config) - - -class PVD1Data(ModelData): - """ - Data for distributed PV. - """ - def __init__(self): - ModelData.__init__(self) - - self.bus = IdxParam(model='Bus', - info="interface bus id", - mandatory=True, - ) - - self.gen = IdxParam(info="static generator index", - mandatory=True, - ) - - self.Sn = NumParam(default=100.0, tex_name='S_n', - info='Model MVA base', - unit='MVA', - ) - - self.xc = NumParam(default=0.0, tex_name='x_c', - info='coupling reactance', - unit='p.u.', - z=True, - ) - - self.igreg = IdxParam(model='Bus', - info='Remote bus idx for droop response, None for local', - ) - - self.qmx = NumParam(default=0.33, tex_name='q_{mx}', - info='Max. reactive power command', - power=True, - ) - - self.qmn = NumParam(default=-0.33, tex_name='q_{mn}', - info='Min. reactive power command', - power=True, - ) - - self.v0 = NumParam(default=0.0, tex_name='v_0', - info='Lower limit of deadband for Vdroop response', - unit='pu', - ) - self.v1 = NumParam(default=0.0, tex_name='v_1', - info='Upper limit of deadband for Vdroop response', - unit='pu', - ) - - self.dqdv = NumParam(default=0.0, tex_name='dq/dv', - info='Q-V droop characteristics', - power=True, - ) - - self.fdbd = NumParam(default=-0.01, tex_name='f_{dbd}', - info='frequency deviation deadband', - ) - - self.ddn = NumParam(default=0.0, tex_name='D_{dn}', - info='Gain after f deadband', - ) - - self.ialim = NumParam(default=1.3, tex_name='I_{alim}', - info='Apparent power limit', - current=True, - ) - - self.vt0 = NumParam(default=0.88, tex_name='V_{t0}', - info='Voltage tripping response curve point 0', - ) - - self.vt1 = NumParam(default=0.90, tex_name='V_{t1}', - info='Voltage tripping response curve point 1', - ) - - self.vt2 = NumParam(default=1.1, tex_name='V_{t2}', - info='Voltage tripping response curve point 2', - ) - - self.vt3 = NumParam(default=1.2, tex_name='V_{t3}', - info='Voltage tripping response curve point 3', - ) - - self.vrflag = NumParam(default=0.0, tex_name='z_{VR}', - info='Voltage tripping is latching (0) or partially self-resetting (0-1)', - ) - - self.ft0 = NumParam(default=59.5, tex_name='f_{t0}', - info='Frequency tripping response curve point 0', - ) - - self.ft1 = NumParam(default=59.7, tex_name='f_{t1}', - info='Frequency tripping response curve point 1', - ) - - self.ft2 = NumParam(default=60.3, tex_name='f_{t2}', - info='Frequency tripping response curve point 2', - ) - - self.ft3 = NumParam(default=60.5, tex_name='f_{t3}', - info='Frequency tripping response curve point 3', - ) - - self.frflag = NumParam(default=0.0, tex_name='z_{FR}', - info='Frequency tripping is latching (0) or partially self-resetting (0-1)', - ) - - self.tip = NumParam(default=0.02, tex_name='T_{ip}', - info='Inverter active current lag time constant', - unit='s', - ) - - self.tiq = NumParam(default=0.02, tex_name='T_{iq}', - info='Inverter reactive current lag time constant', - unit='s', - ) - - -class PVD1Model(Model): - """ - Model implementation of PVD1. - """ - def __init__(self, system, config): - Model.__init__(self, system, config) - - self.flags.tds = True - self.group = 'DG' - - self.buss = DataSelect(self.igreg, self.bus, - info='selected bus (bus or igreg)', - ) - - # initial voltages from selected bus - self.v = ExtAlgeb(model='Bus', src='v', indexer=self.buss, tex_name='V', - info='bus (or igreg) terminal voltage', - unit='p.u.', - ) - - self.a = ExtAlgeb(model='Bus', src='a', indexer=self.buss, tex_name=r'\theta', - info='bus (or igreg) phase angle', - unit='rad.', - ) - - self.p0 = ExtService(model='StaticGen', - src='p', - indexer=self.gen, - tex_name='P_0', - ) - self.q0 = ExtService(model='StaticGen', - src='q', - indexer=self.gen, - tex_name='Q_0', - ) - - -class PVD1(PVD1Data, PVD1Model): - """ - Distributed PV model. (TODO: work in progress) - - Power rating specified in `Sn`. - - Reference: ESIG, WECC Distributed and Small PV Plants Generic Model (PVD1), [Online], - Available: https://www.esig.energy/wiki-main-page/wecc-distributed-and-small-pv-plants-generic-model-pvd1/ - """ - def __init__(self, system, config): - PVD1Data.__init__(self) - PVD1Model.__init__(self, system, config) diff --git a/andes/plot.py b/andes/plot.py index dc0e16a12..33f143890 100644 --- a/andes/plot.py +++ b/andes/plot.py @@ -416,7 +416,7 @@ def guess_event_time(self): def bqplot_data(self, xdata, ydata, *, xheader=None, yheader=None, xlabel=None, ylabel=None, left=None, right=None, ymin=None, ymax=None, legend=True, grid=False, fig=None, latex=True, dpi=150, line_width=1.0, greyscale=False, savefig=None, save_format=None, - show=True, title=None, + title=None, **kwargs): """ Plot with ``bqplot``. Experimental and incomplete. diff --git a/andes/routines/base.py b/andes/routines/base.py index 15a682aa6..935087edb 100644 --- a/andes/routines/base.py +++ b/andes/routines/base.py @@ -43,8 +43,8 @@ def init(self): def run(self, **kwargs): raise NotImplementedError - def summary(self): + def summary(self, **kwargs): raise NotImplementedError - def report(self): + def report(self, **kwargs): raise NotImplementedError diff --git a/andes/routines/eig.py b/andes/routines/eig.py index bb78384d1..558b79a15 100644 --- a/andes/routines/eig.py +++ b/andes/routines/eig.py @@ -1,13 +1,17 @@ -import logging -import scipy.io +""" +Module for eigenvalue analysis. +""" +import logging from math import ceil, pi -from cvxopt import mul, div, spdiag -from cvxopt.lapack import gesv + +import scipy.io +from andes.shared import mul, div, spdiag, gesv from andes.io.txt import dump_data from andes.utils.misc import elapsed from andes.routines.base import BaseRoutine +from andes.variables.report import report_info from andes.shared import np, matrix, spmatrix, plt, mpl from andes.plot import set_latex @@ -142,10 +146,10 @@ def run(self, **kwargs): if system.PFlow.converged is False: logger.warning('Power flow not solved. Eig analysis will not continue.') return succeed - else: - if system.TDS.initialized is False: - system.TDS.init() - system.TDS._itm_step() + + if system.TDS.initialized is False: + system.TDS.init() + system.TDS.itm_step() if system.dae.n == 0: logger.error('No dynamic model. Eig analysis will not continue.') @@ -153,7 +157,7 @@ def run(self, **kwargs): else: if sum(system.dae.Tf != 0) != len(system.dae.Tf): self.singular_idx = np.argwhere(np.equal(system.dae.Tf, 0.0)).ravel().astype(int) - logger.info(f"System contains {len(self.singular_idx)} zero time constants. ") + logger.info("System contains %d zero time constants. ", len(self.singular_idx)) logger.debug([system.dae.x_name[i] for i in self.singular_idx]) self.x_name = np.array(system.dae.x_name) @@ -183,6 +187,9 @@ def run(self, **kwargs): 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=150, show=True, latex=True): + """ + Plot utility for eigenvalues in the S domain. + """ mpl.rc('font', family='Times New Roman', size=12) if mu is None: @@ -206,11 +213,10 @@ def plot(self, mu=None, fig=None, ax=None, left=-6, right=0.5, ymin=-8, ymax=8, if len(p_mu_real) > 0: logger.warning( - 'System is not stable due to {} positive eigenvalues.'.format( - len(p_mu_real))) + 'System is not stable due to %d positive eigenvalues.', len(p_mu_real)) else: logger.info( - 'System is small-signal stable in the initial neighbourhood.') + 'System is small-signal stable in the initial neighborhood.') set_latex(latex) @@ -257,33 +263,25 @@ def export_state_matrix(self): 'x_tex_name': np.array(system.dae.x_tex_name, dtype=np.object), } scipy.io.savemat(system.files.mat, mdict=out) - logger.info(f'State matrix saved to "{system.files.mat}".') + logger.info('State matrix saved to "%s"', system.files.mat) return True - def report(self, x_name=None): + def report(self, x_name=None, **kwargs): """ - Save eigenvalue analysis reports + Save eigenvalue analysis reports. Returns ------- None """ - from andes.variables.report import report_info - - system = self.system - mu = self.mu - part_fact = self.part_fact if x_name is None: x_name = self.x_name - text = [] - header = [] - rowname = [] - data = [] + text, header, rowname, data = list(), list(), list(), list() - neig = len(mu) - mu_real = mu.real() - mu_imag = mu.imag() + neig = len(self.mu) + mu_real = self.mu.real() + mu_imag = self.mu.imag() n_positive = sum(1 for x in mu_real if x > 0) n_zeros = sum(1 for x in mu_real if x == 0) n_negative = sum(1 for x in mu_real if x < 0) @@ -302,7 +300,7 @@ def report(self, x_name=None): freq = [0] * neig ufreq = [0] * neig damping = [0] * neig - for idx, item in enumerate(mu): + for idx, item in enumerate(self.mu): if item.imag == 0: freq[idx] = 0 ufreq[idx] = 0 @@ -315,7 +313,7 @@ def report(self, x_name=None): # obtain most associated variables var_assoc = [] for prow in range(neig): - temp_row = part_fact[prow, :] + temp_row = self.part_fact[prow, :] name_idx = list(temp_row).index(max(temp_row)) var_assoc.append(x_name[name_idx]) @@ -323,7 +321,7 @@ def report(self, x_name=None): for prow in range(neig): temp_row = [] for pcol in range(neig): - temp_row.append(round(part_fact[prow, pcol], 5)) + temp_row.append(round(self.part_fact[prow, pcol], 5)) pf.append(temp_row) # opening info section @@ -372,5 +370,5 @@ def report(self, x_name=None): rowname.append(x_name) data.append(pf[start:end]) - dump_data(text, header, rowname, data, system.files.eig) - logger.info(f'Report saved to "{system.files.eig}".') + dump_data(text, header, rowname, data, self.system.files.eig) + logger.info('Report saved to "%s".', self.system.files.eig) diff --git a/andes/routines/pflow.py b/andes/routines/pflow.py index 064256ce8..d18ba2076 100644 --- a/andes/routines/pflow.py +++ b/andes/routines/pflow.py @@ -1,3 +1,8 @@ +""" +Module for power flow calculation. +""" + +import logging from collections import OrderedDict from andes.utils.misc import elapsed @@ -5,7 +10,6 @@ from andes.variables.report import Report from andes.shared import np, matrix, sparse, newton_krylov, IP_ADD -import logging logger = logging.getLogger(__name__) @@ -64,6 +68,7 @@ def init(self): self.y_sol = None self.system.init(self.models, routine='pflow') + logger.info('Power flow initialized.') # force precompile if numba is on - improves timing accuracy if self.system.config.numba: @@ -71,7 +76,6 @@ def init(self): self.system.g_update(self.models) self.system.j_update(models=self.models) - logger.info('Power flow initialized.') return self.system.dae.xy def nr_step(self): @@ -126,9 +130,15 @@ def summary(self): Output a summary for the PFlow routine. """ ipadd_status = 'CVXOPT normal (ipadd not available)' + + # extract package name, `cvxopt` or `kvxopt` + sp_module = sparse.__module__ + if '.' in sp_module: + sp_module = sp_module.split('.')[0] + if IP_ADD: if self.system.config.ipadd: - ipadd_status = 'Fast in-place' + ipadd_status = f'Fast in-place ({sp_module})' else: ipadd_status = 'CVXOPT normal (ipadd disabled in config)' @@ -162,18 +172,18 @@ def run(self, **kwargs): self.niter = 0 while True: mis = self.nr_step() - logger.info(f'{self.niter}: |F(x)| = {mis:<10g}') + logger.info('%d: |F(x)| = %.10g', self.niter, mis) if mis < self.config.tol: self.converged = True break - elif self.niter > self.config.max_iter: + if self.niter > self.config.max_iter: break - elif np.isnan(mis).any(): + if np.isnan(mis).any(): logger.error('NaN found in solution. Convergence not likely') self.niter = self.config.max_iter + 1 break - elif mis > 1e4 * self.mis[0]: + if mis > 1e4 * self.mis[0]: logger.error('Mismatch increased too fast. Convergence not likely.') break self.niter += 1 @@ -185,12 +195,12 @@ def run(self, **kwargs): max_idx = np.argmax(np.abs(system.dae.xy)) name = system.dae.xy_name[max_idx] logger.error('Mismatch is not correctable possibly due to large load-generation imbalance.') - logger.error(f'Largest mismatch on equation associated with <{name}>') + logger.error('Largest mismatch on equation associated with <%s>', name) else: - logger.error(f'Power flow failed after {self.niter + 1} iterations for {system.files.case}.') + logger.error('Power flow failed after %d iterations for "%s".', self.niter + 1, system.files.case) else: - logger.info(f'Converged in {self.niter+1} iterations in {s1}.') + logger.info('Converged in %d iterations in %s.', self.niter + 1, s1) # make a copy of power flow solutions self.x_sol = system.dae.x.copy() diff --git a/andes/routines/tds.py b/andes/routines/tds.py index cbb2f2f63..59635bcdf 100644 --- a/andes/routines/tds.py +++ b/andes/routines/tds.py @@ -1,7 +1,12 @@ +""" +ANDES module for time-domain simulation. +""" + import sys import os import time import importlib +import logging from collections import OrderedDict from andes.routines.base import BaseRoutine @@ -10,7 +15,6 @@ from andes.shared import tqdm, np, pd from andes.shared import matrix, sparse, spdiag -import logging logger = logging.getLogger(__name__) @@ -102,6 +106,13 @@ def __init__(self, system=None, config=None): self.qrt_start = None self.headroom = 0.0 + # internal storage for iterations + self.x0 = None + self.y0 = None + self.f0 = None + self.Ac = None + self.inc = None + def init(self): """ Initialize the status, storage and values for TDS. @@ -167,12 +178,17 @@ def init(self): # if `dae.n == 1`, `calc_h_first` depends on new `dae.gy` self.calc_h() + # allocate for internal variables + self.x0 = np.zeros_like(system.dae.x) + self.y0 = np.zeros_like(system.dae.y) + self.f0 = np.zeros_like(system.dae.f) + _, s1 = elapsed(t0) if self.initialized is True: - logger.info(f"Initialization for dynamics was successful in {s1}.") + logger.info("Initialization for dynamics was successful in %s.", s1) else: - logger.error(f"Initialization for dynamics failed in {s1}.") + logger.error("Initialization for dynamics failed in %s.", s1) if system.dae.n == 0: tqdm.write('No dynamic component loaded.') @@ -270,7 +286,7 @@ def run(self, no_pbar=False, no_summary=False, **kwargs): step_status = False # call the stepping method of the integration method (or data replay) if self.data_csv is None: - step_status = self._itm_step() # compute for the current step + step_status = self.itm_step() # compute for the current step else: step_status = self._csv_step() @@ -300,7 +316,7 @@ def run(self, no_pbar=False, no_summary=False, **kwargs): # if the ending time has passed if time.time() - rt_end > 0: - logger.debug(f'Simulation over-run at t={dae.t:4.4g} s.') + logger.debug('Simulation over-run at t=%4.4g s.', dae.t) else: self.headroom += (rt_end - time.time()) @@ -320,7 +336,7 @@ def run(self, no_pbar=False, no_summary=False, **kwargs): if self.busted: logger.error(self.err_msg) - logger.error(f"Simulation terminated at t={system.dae.t:.4f}.") + logger.error("Simulation terminated at t=%.4f s.", system.dae.t) system.exit_code += 1 elif system.dae.t == self.config.tf: succeed = True # success flag @@ -329,10 +345,10 @@ def run(self, no_pbar=False, no_summary=False, **kwargs): system.exit_code += 1 _, s1 = elapsed(t0) - logger.info(f'Simulation completed in {s1}.') + logger.info('Simulation completed in %s.', s1) if config.qrt: - logger.debug(f'QRT headroom time: {self.headroom} s.') + logger.debug('QRT headroom time: %.4g s.', self.headroom) system.TDS.save_output() @@ -346,7 +362,7 @@ def run(self, no_pbar=False, no_summary=False, **kwargs): return succeed - def _itm_step(self): + def itm_step(self): """ Integrate with Implicit Trapezoidal Method (ITM) to the current time. @@ -366,12 +382,12 @@ def _itm_step(self): self.niter = 0 self.converged = False - self.x0 = np.array(dae.x) - self.y0 = np.array(dae.y) - self.f0 = np.array(dae.f) + self.x0[:] = dae.x + self.y0[:] = dae.y + self.f0[:] = dae.f while True: - self._fg_update(models=system.exist.pflow_tds) + self.fg_update(models=system.exist.pflow_tds) # lazy Jacobian update @@ -636,7 +652,7 @@ def test_init(self): Update f and g to see if initialization is successful. """ system = self.system - self._fg_update(system.exist.pflow_tds) + self.fg_update(system.exist.pflow_tds) system.j_update(models=system.exist.pflow_tds) # warn if variables are initialized at limits @@ -648,26 +664,27 @@ def test_init(self): if np.max(np.abs(system.dae.fg)) < self.config.tol: logger.debug('Initialization tests passed.') return True - else: - fail_idx = np.where(abs(system.dae.fg) >= self.config.tol) - fail_names = [system.dae.xy_name[int(i)] for i in np.ravel(fail_idx)] + # otherwise, show suspect initialization error + fail_idx = np.where(abs(system.dae.fg) >= self.config.tol) + fail_names = [system.dae.xy_name[int(i)] for i in np.ravel(fail_idx)] - title = 'Suspect initialization issue! Simulation may crash!' - err_data = {'Name': fail_names, - 'Var. Value': system.dae.xy[fail_idx], - 'Eqn. Mismatch': system.dae.fg[fail_idx], - } - tab = Tab(title=title, - header=err_data.keys(), - data=list(map(list, zip(*err_data.values())))) + title = 'Suspect initialization issue! Simulation may crash!' + err_data = {'Name': fail_names, + 'Var. Value': system.dae.xy[fail_idx], + 'Eqn. Mismatch': system.dae.fg[fail_idx], + } + tab = Tab(title=title, + header=err_data.keys(), + data=list(map(list, zip(*err_data.values()))), + ) - logger.error(tab.draw()) + logger.error(tab.draw()) - if system.options.get('verbose') == 1: - breakpoint() - system.exit_code += 1 - return False + if system.options.get('verbose') == 1: + breakpoint() + system.exit_code += 1 + return False def save_output(self, npz=True): """ @@ -685,19 +702,18 @@ def save_output(self, npz=True): """ if self.system.files.no_output: return False - else: - t0, _ = elapsed() - self.system.dae.write_lst(self.system.files.lst) + t0, _ = elapsed() - if npz is True: - self.system.dae.write_npz(self.system.files.npz) - else: - self.system.dae.write_npy(self.system.files.npy) + self.system.dae.write_lst(self.system.files.lst) + if npz is True: + self.system.dae.write_npz(self.system.files.npz) + else: + self.system.dae.write_npy(self.system.files.npy) - _, s1 = elapsed(t0) - logger.info(f'TDS outputs saved in {s1}.') - return True + _, s1 = elapsed(t0) + logger.info('TDS outputs saved in %s.', s1) + return True def do_switch(self): """ @@ -740,7 +756,7 @@ def do_switch(self): return ret - def _fg_update(self, models): + def fg_update(self, models): """ Update `f` and `g` equations. """ @@ -773,7 +789,7 @@ def _fg_wrapper(self, xy): system.dae.y[:] = xy[system.dae.n:] system.vars_to_models() - self._fg_update(system.exist.pflow_tds) + self.fg_update(system.exist.pflow_tds) return system.dae.fg @@ -782,19 +798,23 @@ def _load_pert(self): Load perturbation files to ``self.callpert``. """ system = self.system - if system.files.pert: - if not os.path.isfile(system.files.pert): - logger.warning(f'Pert file not found at "{system.files.pert}".') - return False - - sys.path.append(system.files.case_path) - _, full_name = os.path.split(system.files.pert) - name, ext = os.path.splitext(full_name) - - module = importlib.import_module(name) - self.callpert = getattr(module, 'pert') - logger.info(f'Perturbation file "{system.files.pert}" loaded.') - return True + if not system.files.pert: + return False + + if not os.path.isfile(system.files.pert): + logger.warning('Pert file not found at "%s".', system.files.pert) + return False + + pert_path, full_name = os.path.split(system.files.pert) + logger.debug('Pert file "%s" located at path %s', full_name, pert_path) + + sys.path.append(pert_path) + name, _ = os.path.splitext(full_name) + + module = importlib.import_module(name) + self.callpert = getattr(module, 'pert') + logger.info('Perturbation file "%s" loaded.', system.files.pert) + return True def _load_csv(self, csv_file): """ @@ -811,7 +831,7 @@ def _load_csv(self, csv_file): data = df.to_numpy() if data.ndim != 2: - raise ValueError("Data from CSV is not 2-dimentional (time versus variable)") + raise ValueError("Data from CSV is not 2-dimensional (time versus variable)") if data.shape[0] < 2: logger.warning("CSV data does not contain more than one time step.") @@ -831,7 +851,8 @@ def _debug_g(self, y_idx): Index of the equation into the `g` array. Diff. eqns. are not counted in. """ y_idx = y_idx.tolist() - logger.debug(f'Max. algebraic mismatch associated with {self.system.dae.y_name[y_idx]} [y_idx={y_idx}]') + logger.debug('Max. algebraic mismatch associated with <%s> [y_idx=%d]', + self.system.dae.y_name[y_idx], y_idx) assoc_vars = self.system.dae.gy[y_idx, :] vars_idx = np.where(np.ravel(matrix(assoc_vars)))[0] @@ -839,9 +860,7 @@ def _debug_g(self, y_idx): logger.debug(f'{"y_index":<10} {"Variable":<20} {"Derivative":<20}') for v in vars_idx: v = v.tolist() - logger.debug(f'{v:<10} {self.system.dae.y_name[v]:<20} {assoc_vars[v]:<20g}') - - pass + logger.debug('%10d %20s %20g', v, self.system.dae.y_name[v], assoc_vars[v]) def _debug_ac(self, xy_idx): """ @@ -860,8 +879,8 @@ def _debug_ac(self, xy_idx): eqns_idx = np.where(np.ravel(matrix(assoc_eqns)))[0] vars_idx = np.where(np.ravel(matrix(assoc_vars)))[0] - logger.debug(f'Max. correction is for variable {self.system.dae.xy_name[xy_idx]} [{xy_idx}]') - logger.debug(f'Associated equation value is {self.system.dae.fg[xy_idx]:<20g}') + logger.debug('Max. correction is for variable %s [%d]', self.system.dae.xy_name[xy_idx], xy_idx) + logger.debug('Associated equation value is %20g', self.system.dae.fg[xy_idx]) logger.debug('') logger.debug(f'{"xy_index":<10} {"Equation":<20} {"Derivative":<20} {"Eq. Mismatch":<20}') diff --git a/andes/shared.py b/andes/shared.py index 22d392874..45548500e 100644 --- a/andes/shared.py +++ b/andes/shared.py @@ -1,38 +1,61 @@ """ Shared constants and delayed imports. """ - -# # Known issues of LazyImport :: # # 1) High overhead when called hundreds of thousands times. # For example, NumPy must not be imported with LazyImport. -# from andes.utils.lazyimport import LazyImport import math -import coloredlogs # NOQA -import numpy as np # NOQA -from tqdm import tqdm # NOQA +import coloredlogs # NOQA +import numpy as np # NOQA +from numpy import ndarray # NOQA +from tqdm import tqdm # NOQA -import cvxopt # NOQA -from cvxopt import umfpack # NOQA -from cvxopt import spmatrix, matrix, sparse, spdiag # NOQA -from numpy import ndarray # NOQA - -from andes.utils.texttable import Texttable # NOQA -from andes.utils.paths import get_dot_andes_path # NOQA +# Library preference: +# KVXOPT + ipadd > CVXOPT + ipadd > KXVOPT > CVXOPT (+ KLU or not) try: - from cvxoptklu import klu + import kvxopt + from kvxopt import umfpack # test if shared libs can be found + from kvxopt import spmatrix as kspmatrix + KIP_ADD = False + if hasattr(kspmatrix, 'ipadd'): + KIP_ADD = True except ImportError: - klu = None + kvxopt = None + kspmatrix = None + KIP_ADD = False + -if hasattr(spmatrix, 'ipadd'): - IP_ADD = True +from cvxopt import spmatrix as cspmatrix +if hasattr(cspmatrix, 'ipadd'): + CIP_ADD = True else: - IP_ADD = False + CIP_ADD = False + + +if kvxopt is None or (KIP_ADD is False and CIP_ADD is True): + from cvxopt import umfpack # NOQA + from cvxopt import spmatrix, matrix, sparse, spdiag # NOQA + from cvxopt import mul, div # NOQA + from cvxopt.lapack import gesv # NOQA + try: + from cvxoptklu import klu # NOQA + except ImportError: + klu = None + IP_ADD = CIP_ADD +else: + from kvxopt import umfpack, klu # NOQA + from kvxopt import spmatrix, matrix, sparse, spdiag # NOQA + from kvxopt import mul, div # NOQA + from kvxopt.lapack import gesv # NOQA + IP_ADD = KIP_ADD + +from andes.utils.texttable import Texttable # NOQA +from andes.utils.paths import get_dot_andes_path # NOQA # --- constants --- @@ -49,12 +72,14 @@ pd = LazyImport('import pandas') cupy = LazyImport('import cupy') -plt = LazyImport('from matplotlib import pyplot') mpl = LazyImport('import matplotlib') -Pool = LazyImport('from pathos.multiprocessing import Pool') -Process = LazyImport('from multiprocess import Process') unittest = LazyImport('import unittest') yaml = LazyImport('import yaml') + +plt = LazyImport('from matplotlib import pyplot') +Pool = LazyImport('from pathos.multiprocessing import Pool') +Process = LazyImport('from multiprocess import Process') + newton_krylov = LazyImport('from scipy.optimize import newton_krylov') fsolve = LazyImport('from scipy.optimize import fsolve') solve_ivp = LazyImport('from scipy.integrate import solve_ivp') diff --git a/andes/system.py b/andes/system.py index d0f84766d..5ed458880 100644 --- a/andes/system.py +++ b/andes/system.py @@ -18,6 +18,7 @@ import os import sys import inspect +import dill from collections import OrderedDict from typing import List, Dict, Tuple, Union, Optional @@ -308,7 +309,6 @@ def _prepare_mp(self, quick=False): Function is not working. Serialization failed for `conj`. """ from andes.shared import Pool - import dill dill.settings['recurse'] = True # consistency check for group parameters and variables @@ -813,11 +813,10 @@ def j_update(self, models: OrderedDict, info=None): f'j_size={j_size}') raise e - msg = f"Jacobian updated at t={self.dae.t}" if info: - msg += f' due to {info}' - - logger.debug(msg) + logger.debug("Jacobian updated at t=%.6f due to %s.", self.dae.t, info) + else: + logger.debug("Jacobian updated at t=%.6f.", self.dae.t) def store_sparse_pattern(self, models: OrderedDict): """ @@ -989,9 +988,9 @@ def find_models(self, flag: Optional[Union[str, Tuple]], skip_zero: bool = True) def dill(self): """ - Serialize generated numerical functions in `System.calls` with package `dill`. + Serialize generated numerical functions in ``System.calls`` with package ``dill``. - The serialized file will be stored to ``~/andes/calls.pkl``, where `~` is the home directory path. + The serialized file will be stored to ``~/.andes/calls.pkl``, where `~` is the home directory path. Notes ----- @@ -999,7 +998,6 @@ def dill(self): """ logger.debug("Dumping calls to calls.pkl with dill") - import dill dill.settings['recurse'] = True pkl_path = get_pkl_path() @@ -1008,7 +1006,9 @@ def dill(self): @staticmethod def _load_pkl(): - import dill + """ + Helper function to open and load dill-pickled functions. + """ dill.settings['recurse'] = True pkl_path = get_pkl_path() @@ -1027,7 +1027,7 @@ def _load_pkl(): def undill(self): """ - Deserialize the function calls from ``~/andes.calls.pkl`` with dill. + Deserialize the function calls from ``~/.andes/calls.pkl`` with ``dill``. If no change is made to models, future calls to ``prepare()`` can be replaced with ``undill()`` for acceleration. @@ -1055,11 +1055,13 @@ def undill(self): # try to replace equations and jacobian calls with saved code if pycode is not None and self.config.use_pycode: for model in self.models.values(): - model.calls.f = pycode.__dict__[model.class_name].__dict__.get("f_update") - model.calls.g = pycode.__dict__[model.class_name].__dict__.get("g_update") + pycode_model = pycode.__dict__[model.class_name] + + model.calls.f = pycode_model.__dict__.get("f_update") + model.calls.g = pycode_model.__dict__.get("g_update") for jname in model.calls.j: - model.calls.j[jname] = pycode.__dict__[model.class_name].__dict__[f'{jname}_update'] + model.calls.j[jname] = pycode_model.__dict__.get(f'{jname}_update') logger.info("Using generated Python code for equations and Jacobians.") else: logger.debug("Using undilled lambda functions.") @@ -1222,7 +1224,6 @@ def import_models(self): ``system.models['Bus']`` points the same instance as ``system.Bus``. """ - # non-JIT models for file, cls_list in file_classes.items(): for model_name in cls_list: the_module = importlib.import_module('andes.models.' + file) @@ -1235,13 +1236,6 @@ def import_models(self): group_name = self.__dict__[model_name].group self.__dict__[group_name].add_model(model_name, self.__dict__[model_name]) - # *** JIT import code *** - # import JIT models - # for file, pair in jits.items(): - # for cls, name in pair.items(): - # self.__dict__[name] = JIT(self, file, cls, name) - # *********************** - def import_routines(self): """ Import routines as defined in ``routines/__init__.py``. @@ -1365,6 +1359,10 @@ def _store_calls(self): self.calls[name] = mdl.calls def _list2array(self): + """ + Helper function to call models' ``list2array`` method, which + usually performs memory preallocation. + """ self.call_models('list2array', self.models) def set_config(self, config=None): diff --git a/andes/utils/func.py b/andes/utils/func.py index 173d39dde..f2f36a94b 100644 --- a/andes/utils/func.py +++ b/andes/utils/func.py @@ -10,8 +10,8 @@ def list_flatten(input_list): """ if len(input_list) > 0 and isinstance(input_list[0], (list, np.ndarray)): return functools.reduce(operator.iconcat, input_list, []) - else: - return input_list + + return input_list def interp_n2(t, x, y): diff --git a/docs/source/install.rst b/docs/source/install.rst index d7686f5a6..a4e5cbb17 100644 --- a/docs/source/install.rst +++ b/docs/source/install.rst @@ -165,46 +165,37 @@ Change directory to the ANDES source code folder that contains ``setup.py`` and Performance Packages (Advanced) =============================== -The following two forks of ``cvxopt``, ``cvxoptklu``, ``cvxopt`` with ``spmatrix.ipadd`` +The following two forks of ``cvxopt``: ``kvxopt`` and ``cvxopt`` with ``spmatrix.ipadd`` are optional but can significantly boost the performance of ANDES. -**Installation requires a C compiler**, ``openblas`` and ``SuiteSparse`` libraries. .. note:: Performance packages can be safely skipped and will not affect the functionality of ANDES. -.. warning:: - - We have not tried to compile either package on Windows. - Refer to the CVXOPT installation instructions for Windows at - http://cvxopt.org/install/index.html#windows - -cxvoptklu ---------- -``cvxoptklu`` is a fork of the CVXOPT with KLU by Uriel Sandoval (@sanurielf). -In addition to UMFPACK, ``cvxoptklu`` interfaces ``cvxopt`` to KLU, which is +KVXOPT +------ +``KVXOPT`` is a fork of the CVXOPT with KLU by Uriel Sandoval (@sanurielf). +In addition to UMFPACK, ``KVXOPT`` interfaces ``cvxopt`` to KLU, which is roughly 20% faster than UMFPACK for circuit simulation based on our testing. -To install ``cvxoptklu``, on Debian GNU/Linux, one can do +To install ``KVXOPT`` run .. code:: bash - sudo apt install libopenblas-dev libsuitesparse-dev - pip install cvxoptklu + python -m pip install kvxopt -On macOS, one can install with homebrew using +CVXOPT with ipadd +----------------- -.. code:: bash +**Installation requires a C compiler**, ``openblas`` and ``SuiteSparse`` libraries. - brew install openblas suitesparse - pip install cvxoptklu +.. warning:: -To install from source code, use the repository at -https://github.com/cuihantao/cvxoptklu. + We have not tried to compile CVXOPT on Windows. + Refer to the CVXOPT installation instructions for Windows at + http://cvxopt.org/install/index.html#windows -CVXOPT with ipadd ------------------ To install our fork of ``cvxopt`` with ``spmatrix.ipadd``, one need to clone the repository and compile from source. diff --git a/docs/source/modeling.rst b/docs/source/modeling.rst index 6a058bef5..3ac7ad67b 100644 --- a/docs/source/modeling.rst +++ b/docs/source/modeling.rst @@ -48,7 +48,7 @@ It takes several seconds up to a minute to finish the generation. .. warning:: When models are modified (such as adding new models or changing equation strings), code generation needs to be executed again for consistency. It can be more conveniently triggered from command line with - ``andes prepare -qi``. + ``andes prepare -i``. .. autofunction:: andes.system.System.prepare :noindex: diff --git a/docs/source/release-notes.rst b/docs/source/release-notes.rst index 463ef25ff..6037b5354 100644 --- a/docs/source/release-notes.rst +++ b/docs/source/release-notes.rst @@ -8,6 +8,23 @@ The APIs before v3.0.0 are in beta and may change without prior notice. v1.2 Notes ---------- +v1.2.2 (2020-11-01) +``````````````````` +New Models: + +- ``PVD1`` model, WECC distributed PV model. + Supports multiple PVD1 devices on the same bus. +- Added ``ACEc`` model, ACE calculation with continuous freq. + +Changes and fixes: + +- Renamed `TDS._itm_step` to `TDS.itm_step` as a public API. +- Allow variable `sys_f` (system frequency) in equation strings. +- Fixed ACE equation. + measurement. +- Support ``kvxopt`` as a drop-in replacement for ``cvxopt`` + to bring KLU to Windows (and other platforms). +- Added ``kvxopt`` as a dependency for PyPI installation. v1.2.1 (2020-10-11) ``````````````````` diff --git a/docs/source/tutorial.rst b/docs/source/tutorial.rst index 85a658fd3..b195711b2 100644 --- a/docs/source/tutorial.rst +++ b/docs/source/tutorial.rst @@ -94,14 +94,14 @@ Generated code are stored in the folder ``.andes/calls.pkl`` in your home direct In addition, ``andes selftest`` implicitly calls the code generation. If you are using ANDES as a package in the user mode, you won't need to call it again. -For developers, ``andes prepare`` needs to be called immediately following any model equation -modification. Otherwise, simulation results will not reflect the new equations and will likely lead to an error. - Option ``-q`` or ``--quick`` (enabled by default) can be used to speed up the code generation. It skips the generation of :math:`\LaTeX`-formatted equations, which are only used in documentation and the interactive mode. -Option ``-i`` or ``--incremental`` can be used to further speed up the code generation during model development. +For developers, ``andes prepare`` needs to be called immediately following any model equation +modification. Otherwise, simulation results will not reflect the new equations and will likely lead to an error. +Option ``-i`` or ``--incremental``, instead of ``-q``, can be used to further speed up the code generation +during model development. ``andes prepare -i`` only generates code for models with modified equations. andes run @@ -247,7 +247,7 @@ For example, to convert ``case5.m`` into the ``xlsx`` format, run andes run case5.m --convert xlsx -The output messages will look like +The output messages will look like :: Parsing input file CASE5 Power flow data for modified 5 bus, 5 gen case based on PJM 5-bus system diff --git a/examples/5. profiling.ipynb b/examples/5. profiling.ipynb index eec9ec03c..0322067b5 100644 --- a/examples/5. profiling.ipynb +++ b/examples/5. profiling.ipynb @@ -101,7 +101,7 @@ "\n", " ncalls tottime percall cumtime percall filename:lineno(function)\n", " 1 0.007 0.007 1.597 1.597 /Users/hcui7/repos/andes/andes/routines/tds.py:160(run)\n", - " 602 0.104 0.000 1.533 0.003 /Users/hcui7/repos/andes/andes/routines/tds.py:251(_itm_step)\n", + " 602 0.104 0.000 1.533 0.003 /Users/hcui7/repos/andes/andes/routines/tds.py:251(itm_step)\n", " 2318 0.010 0.000 1.141 0.000 /Users/hcui7/repos/andes/andes/routines/tds.py:590(_fg_update)\n", " 11649 0.092 0.000 1.085 0.000 /Users/hcui7/repos/andes/andes/system.py:1014(call_models)\n", " 2323 0.002 0.000 0.596 0.000 /Users/hcui7/repos/andes/andes/system.py:700(g_update)\n", @@ -375,7 +375,7 @@ } ], "source": [ - "%lprun -f ss.TDS._itm_step ss.TDS.run()" + "%lprun -f ss.TDS.itm_step ss.TDS.run()" ] }, { diff --git a/tests/test_paths.py b/tests/test_paths.py index ee0b67a8c..a4d64659d 100644 --- a/tests/test_paths.py +++ b/tests/test_paths.py @@ -8,16 +8,37 @@ class TestPaths(unittest.TestCase): def setUp(self) -> None: self.kundur = 'kundur/' self.matpower = 'matpower/' + self.ieee14 = andes.get_case("ieee14/ieee14.raw") def test_tree(self): list_cases(self.kundur, no_print=True) list_cases(self.matpower, no_print=True) def test_addfile_path(self): - ieee14 = andes.get_case("ieee14/ieee14.raw") - path, case = os.path.split(ieee14) - andes.load('ieee14.raw', addfile='ieee14.dyr', input_path=path, default_config=True) + path, case = os.path.split(self.ieee14) + andes.load('ieee14.raw', addfile='ieee14.dyr', + input_path=path, default_config=True, + ) - andes.run('ieee14.raw', addfile='ieee14.dyr', input_path=path, + andes.run('ieee14.raw', addfile='ieee14.dyr', + input_path=path, no_output=True, default_config=True, ) + + def test_pert_file(self): + """Test path of pert file""" + path, case = os.path.split(self.ieee14) + + # --- with pert file --- + ss = andes.run('ieee14.raw', pert='pert.py', + input_path=path, no_output=True, default_config=True, + ) + ss.TDS.init() + self.assertIsNotNone(ss.TDS.callpert) + + # --- without pert file --- + ss = andes.run('ieee14.raw', + input_path=path, no_output=True, default_config=True, + ) + ss.TDS.init() + self.assertIsNone(ss.TDS.callpert)