diff --git a/README.md b/README.md index bb31ccd3..bfafabf7 100644 --- a/README.md +++ b/README.md @@ -28,7 +28,10 @@ Python Software for Power System Scheduling Modeling and Co-Simulation with Dyna # Why AMS -With the built-in interface with dynamic simulation engine, ANDES, AMS enables Dynamics Interfaced Stability Constrained Production Cost and Market Operation Modeling. +With the built-in interface with ANDES, AMS enables **Dynamics Incorporated Stability-Constrained Scheduling**. + +AMS is a **Modeling Framework** that provides a descriptive way to formulate scheduling problems. +The optimization problems are then handled by **CVXPY** and solved with third-party solvers. AMS produces credible scheduling results and competitive performance. The following results show the comparison of DCOPF between AMS and other tools. diff --git a/ams/__init__.py b/ams/__init__.py index ae078bea..379d7b01 100644 --- a/ams/__init__.py +++ b/ams/__init__.py @@ -1,6 +1,7 @@ from . import _version __version__ = _version.get_versions()['version'] +from ams import opt # NOQA from ams import benchmarks # NOQA from ams.main import config_logger, load, run # NOQA @@ -10,4 +11,4 @@ __author__ = 'Jining Wang' -__all__ = ['System', 'get_case', 'System'] +__all__ = ['System', 'get_case'] diff --git a/ams/core/documenter.py b/ams/core/documenter.py index 727c2d32..21b85501 100644 --- a/ams/core/documenter.py +++ b/ams/core/documenter.py @@ -223,9 +223,10 @@ def get(self, max_width=78, export='plain'): # add tables self.parent.syms.generate_symbols() out += self._obj_doc(max_width=max_width, export=export) - out += self._constr_doc(max_width=max_width, export=export) out += self._expr_doc(max_width=max_width, export=export) + out += self._constr_doc(max_width=max_width, export=export) out += self._var_doc(max_width=max_width, export=export) + out += self._exprc_doc(max_width=max_width, export=export) out += self._service_doc(max_width=max_width, export=export) out += self._param_doc(max_width=max_width, export=export) out += self.config.doc(max_width=max_width, export=export) @@ -316,17 +317,26 @@ def _constr_doc(self, max_width=78, export='plain'): rest_dict=rest_dict) def _expr_doc(self, max_width=78, export='plain'): - # expression documentation + # Expression documentation if len(self.parent.exprs) == 0: return '' # prepare temporary lists - names, var_names, info = list(), list(), list() + names, info = list(), list() + units, sources, units_rest = list(), list(), list() for p in self.parent.exprs.values(): names.append(p.name) - var_names.append(p.var) info.append(p.info if p.info else '') + units.append(p.unit if p.unit else '') + units_rest.append(f'*{p.unit}*' if p.unit else '') + + slist = [] + if p.owner is not None and p.src is not None: + slist.append(f'{p.owner.class_name}.{p.src}') + elif p.owner is not None and p.src is None: + slist.append(f'{p.owner.class_name}') + sources.append(','.join(slist)) # expressions based on output format expressions = [] @@ -342,13 +352,67 @@ def _expr_doc(self, max_width=78, export='plain'): expressions = math_wrap(expressions, export=export) plain_dict = OrderedDict([('Name', names), - ('Variable', var_names), ('Description', info), + ('Unit', units), + ]) + rest_dict = OrderedDict([('Name', names), + ('Description', info), + ('Expression', expressions), + ('Unit', units_rest), + ('Source', sources), + ]) + + # convert to rows and export as table + return make_doc_table(title=title, + max_width=max_width, + export=export, + plain_dict=plain_dict, + rest_dict=rest_dict) + + def _exprc_doc(self, max_width=78, export='plain'): + # ExpressionCalc documentation + if len(self.parent.exprcs) == 0: + return '' + + # prepare temporary lists + names, info = list(), list() + units, sources, units_rest = list(), list(), list() + + for p in self.parent.exprcs.values(): + names.append(p.name) + info.append(p.info if p.info else '') + units.append(p.unit if p.unit else '') + units_rest.append(f'*{p.unit}*' if p.unit else '') + + slist = [] + if p.owner is not None and p.src is not None: + slist.append(f'{p.owner.class_name}.{p.src}') + elif p.owner is not None and p.src is None: + slist.append(f'{p.owner.class_name}') + sources.append(','.join(slist)) + + # expressions based on output format + expressions = [] + if export == 'rest': + for p in self.parent.exprcs.values(): + expr = _tex_pre(self, p, self.parent.syms.tex_map) + logger.debug(f'{p.name} math: {expr}') + expressions.append(expr) + + title = 'ExpressionCalcs\n----------------------------------' + else: + title = 'ExpressionCalcs' + expressions = math_wrap(expressions, export=export) + + plain_dict = OrderedDict([('Name', names), + ('Description', info), + ('Unit', units), ]) rest_dict = OrderedDict([('Name', names), - ('Variable', var_names), ('Description', info), ('Expression', expressions), + ('Unit', units_rest), + ('Source', sources), ]) # convert to rows and export as table @@ -423,6 +487,8 @@ def _var_doc(self, max_width=78, export='plain'): slist = [] if p.owner is not None and p.src is not None: slist.append(f'{p.owner.class_name}.{p.src}') + elif p.owner is not None and p.src is None: + slist.append(f'{p.owner.class_name}') sources.append(','.join(slist)) # symbols based on output format @@ -490,6 +556,8 @@ def _param_doc(self, max_width=78, export='plain'): slist = [] if p.owner is not None and p.src is not None: slist.append(f'{p.owner.class_name}.{p.src}') + elif p.owner is not None and p.src is None: + slist.append(f'{p.owner.class_name}') sources.append(','.join(slist)) # symbols based on output format diff --git a/ams/core/matprocessor.py b/ams/core/matprocessor.py index 151dac3d..bf0e4ac1 100644 --- a/ams/core/matprocessor.py +++ b/ams/core/matprocessor.py @@ -13,7 +13,7 @@ from andes.shared import tqdm, tqdm_nb from andes.utils.misc import elapsed, is_notebook -from ams.opt.omodel import Param +from ams.opt import Param from ams.shared import pd, sps logger = logging.getLogger(__name__) @@ -140,7 +140,9 @@ def class_name(self): class MatProcessor: """ - Class for matrix processing in AMS system. + Class for matrices processing in AMS system. + The connectivity matrices `Cft`, `Cg`, `Cl`, and `Csh` ***have taken*** the + devices connectivity into account. The MParams' row names and col names are assigned in `System.setup()`. """ diff --git a/ams/core/param.py b/ams/core/param.py index 50c0cc4c..dbefba49 100644 --- a/ams/core/param.py +++ b/ams/core/param.py @@ -10,7 +10,7 @@ import numpy as np from scipy.sparse import issparse -from ams.opt.omodel import Param +from ams.opt import Param logger = logging.getLogger(__name__) diff --git a/ams/core/service.py b/ams/core/service.py index d94ce5e8..2520cce0 100644 --- a/ams/core/service.py +++ b/ams/core/service.py @@ -10,7 +10,7 @@ from andes.core.service import BaseService -from ams.opt.omodel import Param +from ams.opt import Param logger = logging.getLogger(__name__) @@ -33,8 +33,6 @@ class RBaseService(BaseService, Param): Description. vtype : Type, optional Variable type. - model : str, optional - Model name. no_parse: bool, optional True to skip parsing the service. sparse: bool, optional @@ -106,8 +104,8 @@ class ValueService(RBaseService): Description. vtype : Type, optional Variable type. - model : str, optional - Model name. + no_parse: bool, optional + True to skip parsing the service. sparse: bool, optional True to return output as scipy csr_matrix. """ @@ -155,8 +153,8 @@ class ROperationService(RBaseService): Description. vtype : Type, optional Variable type. - model : str, optional - Model name. + no_parse: bool, optional + True to skip parsing the service. sparse: bool, optional True to return output as scipy csr_matrix. """ @@ -194,6 +192,8 @@ class LoadScale(ROperationService): Unit. info : str, optional Description. + no_parse: bool, optional + True to skip parsing the service. sparse: bool, optional True to return output as scipy csr_matrix. """ @@ -236,7 +236,7 @@ class NumOp(ROperationService): Note that the scalar output is converted to a 1D array. - The optional kwargs are passed to the input function. + The `rargs` are passed to the input function. Parameters ---------- @@ -397,8 +397,6 @@ class NumOpDual(NumOp): Description. vtype : Type, optional Variable type. - model : str, optional - Model name. rfun : Callable, optional Function to apply to the output of ``fun``. rargs : dict, optional @@ -407,6 +405,8 @@ class NumOpDual(NumOp): Expand the dimensions of the output array along a specified axis. array_out : bool, optional Whether to force the output to be an array. + no_parse: bool, optional + True to skip parsing the service. sparse: bool, optional True to return output as scipy csr_matrix. """ @@ -467,6 +467,10 @@ class MinDur(NumOpDual): Unit. info : str, optional Description. + vtype : Type, optional + Variable type. + no_parse: bool, optional + True to skip parsing the service. sparse: bool, optional True to return output as scipy csr_matrix. """ @@ -537,8 +541,12 @@ class NumHstack(NumOp): Description. vtype : Type, optional Variable type. - model : str, optional - Model name. + rfun : Callable, optional + Function to apply to the output of ``fun``. + rargs : dict, optional + Keyword arguments to pass to ``rfun``. + no_parse: bool, optional + True to skip parsing the service. sparse: bool, optional True to return output as scipy csr_matrix. """ @@ -624,8 +632,12 @@ class ZonalSum(NumOp): Description. vtype : Type Variable type. - model : str - Model name. + rfun : Callable, optional + Function to apply to the output of ``fun``. + rargs : dict, optional + Keyword arguments to pass to ``rfun``. + no_parse: bool, optional + True to skip parsing the service. sparse: bool, optional True to return output as scipy csr_matrix. """ @@ -676,7 +688,7 @@ class VarSelect(NumOp): A numerical matrix to select a subset of a 2D variable, ``u.v[:, idx]``. - For example, if nned to select Energy Storage output + For example, if need to select Energy Storage output power from StaticGen `pg`, following definition can be used: ```python class RTED: @@ -709,6 +721,8 @@ class RTED: Keyword arguments to pass to ``rfun``. array_out : bool, optional Whether to force the output to be an array. + no_parse: bool, optional + True to skip parsing the service. sparse: bool, optional True to return output as scipy csr_matrix. """ @@ -801,8 +815,12 @@ class VarReduction(NumOp): A description of the operation. vtype : Type, optional The variable type. - model : str, optional - The model name associated with the operation. + rfun : Callable, optional + Function to apply to the output of ``fun``. + rargs : dict, optional + Keyword arguments to pass to ``rfun``. + no_parse: bool, optional + True to skip parsing the service. sparse: bool, optional True to return output as scipy csr_matrix. """ @@ -818,13 +836,11 @@ def __init__(self, rfun: Callable = None, rargs: dict = None, no_parse: bool = False, - sparse: bool = False, - **kwargs): + sparse: bool = False,): super().__init__(name=name, tex_name=tex_name, unit=unit, info=info, vtype=vtype, u=u, fun=None, rfun=rfun, rargs=rargs, - no_parse=no_parse, sparse=sparse, - **kwargs) + no_parse=no_parse, sparse=sparse,) self.fun = fun @property @@ -859,8 +875,12 @@ class RampSub(NumOp): Description. vtype : Type Variable type. - model : str - Model name. + rfun : Callable, optional + Function to apply to the output of ``fun``. + rargs : dict, optional + Keyword arguments to pass to ``rfun``. + no_parse: bool, optional + True to skip parsing the service. sparse: bool, optional True to return output as scipy csr_matrix. """ diff --git a/ams/core/symprocessor.py b/ams/core/symprocessor.py index c1bc280e..365f7a5e 100644 --- a/ams/core/symprocessor.py +++ b/ams/core/symprocessor.py @@ -89,6 +89,9 @@ def __init__(self, parent): (r'(== 0|<= 0)$', ''), # remove the comparison operator (r'cp\.(Minimize|Maximize)', r'float'), # remove cp.Minimize/Maximize (r'\bcp.\b', 'np.'), + (r'\bexp\b', 'np.exp'), + (r'\blog\b', 'np.log'), + (r'\bconj\b', 'np.conj'), ]) self.status = { @@ -124,9 +127,7 @@ def generate_symbols(self, force_generate=False): # Vars for vname, var in self.parent.vars.items(): - tmp = sp.symbols(f'{var.name}') - # tmp = sp.symbols(var.name) - self.inputs_dict[vname] = tmp + self.inputs_dict[vname] = sp.symbols(f'{vname}') self.sub_map[rf"\b{vname}\b"] = f"self.om.{vname}" self.tex_map[rf"\b{vname}\b"] = rf'{var.tex_name}' self.val_map[rf"\b{vname}\b"] = f"rtn.{vname}.v" @@ -163,6 +164,13 @@ def generate_symbols(self, force_generate=False): if not service.no_parse: self.val_map[rf"\b{sname}\b"] = f"rtn.{sname}.v" + # Expressions + for ename, expr in self.parent.exprs.items(): + self.inputs_dict[ename] = sp.symbols(f'{ename}') + self.sub_map[rf"\b{ename}\b"] = f"self.om.{ename}" + self.val_map[rf"\b{ename}\b"] = f"rtn.{ename}.v" + self.tex_map[rf"\b{ename}\b"] = f'{expr.tex_name}' + # Constraints # NOTE: constraints are included in sub_map for ExpressionCalc # thus, they don't have the suffix `.v` diff --git a/ams/interface.py b/ams/interface.py index fe93dd32..d019d752 100644 --- a/ams/interface.py +++ b/ams/interface.py @@ -17,6 +17,7 @@ # Models used in ANDES PFlow +# FIXME: add DC models, e.g. Node pflow_dict = OrderedDict([ ('Bus', create_entry('Vn', 'vmax', 'vmin', 'v0', 'a0', 'xcoord', 'ycoord', 'area', 'zone', @@ -49,10 +50,66 @@ 'pq': 'PQ', } +def sync_adsys(amsys, adsys): + """ + Sync parameters value of PFlow models between AMS and ANDES systems. + + Parameters + ---------- + amsys : AMS.system.System + The AMS system. + adsys : ANDES.system.System + The ANDES system. + + Returns + ------- + bool + True if successful. + """ + for mname, params in pflow_dict.items(): + ad_mdl = adsys.__dict__[mname] + am_mdl = amsys.__dict__[mname] + idx = am_mdl.idx.v + for param in params: + if param in ['idx', 'name']: + continue + # NOTE: when setting list values to DataParam, sometimes run into error + try: + ad_mdl.set(src=param, attr='v', idx=idx, + value=am_mdl.get(src=param, attr='v', idx=idx)) + except Exception: + logger.debug(f"Skip updating {mname}.{param}") + continue + return True + + +def _to_andes_pflow(system, no_output=False, default_config=True, **kwargs): + """ + Helper function to convert the AMS system to an ANDES system with only + power flow models. + """ + + adsys = andes_System(no_outpu=no_output, default_config=default_config, **kwargs) + # FIXME: is there a systematic way to do this? Other config might be needed + adsys.config.freq = system.config.freq + adsys.config.mva = system.config.mva + + for mdl_name, mdl_cols in pflow_dict.items(): + mdl = getattr(system, mdl_name) + mdl.cache.refresh("df_in") # refresh cache + for row in mdl.cache.df_in[mdl_cols].to_dict(orient='records'): + adsys.add(mdl_name, row) + + sync_adsys(amsys=system, adsys=adsys) + + return adsys + + def to_andes(system, addfile=None, setup=False, no_output=False, default_config=True, - verify=True, tol=1e-3): + verify=True, tol=1e-3, + **kwargs): """ Convert the AMS system to an ANDES system. @@ -110,17 +167,8 @@ def to_andes(system, addfile=None, """ t0, _ = elapsed() - adsys = andes_System(no_output=no_output, - default_config=default_config) - # FIXME: is there a systematic way to do this? Other config might be needed - adsys.config.freq = system.config.freq - adsys.config.mva = system.config.mva - - for mdl_name, mdl_cols in pflow_dict.items(): - mdl = getattr(system, mdl_name) - mdl.cache.refresh("df_in") # refresh cache - for row in mdl.cache.df_in[mdl_cols].to_dict(orient='records'): - adsys.add(mdl_name, row) + # --- convert power flow models --- + adsys = _to_andes_pflow(system, no_output=no_output, default_config=default_config, **kwargs) _, s = elapsed(t0) diff --git a/ams/main.py b/ams/main.py index 3df5861e..93d987be 100644 --- a/ams/main.py +++ b/ams/main.py @@ -18,6 +18,7 @@ from ._version import get_versions from andes.main import _find_cases +from andes.main import config_logger as ad_config_logger from andes.shared import Pool, Process, coloredlogs, unittest, NCPUS_PHYSICAL from andes.utils.misc import elapsed, is_interactive @@ -71,6 +72,7 @@ def config_logger(stream_level=logging.INFO, *, Original author: Hantao Cui License: GPL3 """ + ad_config_logger(stream_level) lg = logging.getLogger('ams') lg.setLevel(logging.DEBUG) diff --git a/ams/opt/__init__.py b/ams/opt/__init__.py index d442bbcb..83b5268f 100644 --- a/ams/opt/__init__.py +++ b/ams/opt/__init__.py @@ -2,4 +2,11 @@ Module for optimization modeling. """ -from ams.opt.omodel import Var, Constraint, Objective, OModel # NOQA +from ams.opt.optbase import OptzBase, ensure_symbols, ensure_mats_and_parsed # NOQA +from ams.opt.var import Var # NOQA +from ams.opt.exprcalc import ExpressionCalc # NOQA +from ams.opt.param import Param # NOQA +from ams.opt.constraint import Constraint # NOQA +from ams.opt.objective import Objective # NOQA +from ams.opt.expression import Expression # NOQA +from ams.opt.omodel import OModelBase, OModel # NOQA diff --git a/ams/opt/constraint.py b/ams/opt/constraint.py new file mode 100644 index 00000000..8ab5a508 --- /dev/null +++ b/ams/opt/constraint.py @@ -0,0 +1,172 @@ +""" +Module for optimization Constraint. +""" +import logging + +from typing import Optional +import re + +import numpy as np + +import cvxpy as cp + +from ams.utils import pretty_long_message +from ams.shared import _prefix, _max_length + +from ams.opt import OptzBase, ensure_symbols, ensure_mats_and_parsed + +logger = logging.getLogger(__name__) + + +class Constraint(OptzBase): + """ + Base class for constraints. + + This class is used as a template for defining constraints. Each + instance of this class represents a single constraint. + + Parameters + ---------- + name : str, optional + A user-defined name for the constraint. + e_str : str, optional + A mathematical expression representing the constraint. + info : str, optional + Additional informational text about the constraint. + is_eq : str, optional + Flag indicating if the constraint is an equality constraint. False indicates + an inequality constraint in the form of `<= 0`. + + Attributes + ---------- + is_disabled : bool + Flag indicating if the constraint is disabled, False by default. + rtn : ams.routines.Routine + The owner routine instance. + is_disabled : bool, optional + Flag indicating if the constraint is disabled, False by default. + dual : float, optional + The dual value of the constraint. + code : str, optional + The code string for the constraint + """ + + def __init__(self, + name: Optional[str] = None, + e_str: Optional[str] = None, + info: Optional[str] = None, + is_eq: Optional[bool] = False, + ): + OptzBase.__init__(self, name=name, info=info) + self.e_str = e_str + self.is_eq = is_eq + self.is_disabled = False + self.dual = None + self.code = None + + @ensure_symbols + def parse(self): + """ + Parse the constraint. + """ + # parse the expression str + sub_map = self.om.rtn.syms.sub_map + code_constr = self.e_str + for pattern, replacement in sub_map.items(): + try: + code_constr = re.sub(pattern, replacement, code_constr) + except TypeError as e: + raise TypeError(f"Error in parsing constr <{self.name}>.\n{e}") + # parse the constraint type + code_constr += " == 0" if self.is_eq else " <= 0" + # store the parsed expression str code + self.code = code_constr + msg = f" - Constr <{self.name}>: {self.e_str}" + logger.debug(pretty_long_message(msg, _prefix, max_length=_max_length)) + return True + + def _evaluate_expression(self, code, local_vars=None): + """ + Helper method to evaluate the expression code. + + Parameters + ---------- + code : str + The code string representing the expression. + + Returns + ------- + cp.Expression + The evaluated cvxpy expression. + """ + return eval(code, {}, local_vars) + + @ensure_mats_and_parsed + def evaluate(self): + """ + Evaluate the constraint. + """ + msg = f" - Constr <{self.name}>: {self.code}" + logger.debug(pretty_long_message(msg, _prefix, max_length=_max_length)) + try: + local_vars = {'self': self, 'cp': cp, 'sub_map': self.om.rtn.syms.val_map} + self.optz = self._evaluate_expression(self.code, local_vars=local_vars) + except Exception as e: + raise Exception(f"Error in evaluating Constraint <{self.name}>.\n{e}") + + def __repr__(self): + enabled = 'OFF' if self.is_disabled else 'ON' + out = f"{self.class_name}: {self.name} [{enabled}]" + return out + + @property + def e(self): + """ + Return the calculated constraint LHS value. + Note that `v` should be used primarily as it is obtained + from the solver directly. + + `e` is for debugging purpose. For a successfully solved problem, + `e` should equal to `v`. However, when a problem is infeasible + or unbounded, `e` can be used to check the constraint LHS value. + """ + if self.code is None: + logger.info(f"Constraint <{self.name}> is not parsed yet.") + return None + + val_map = self.om.rtn.syms.val_map + code = self.code + for pattern, replacement in val_map.items(): + try: + code = re.sub(pattern, replacement, code) + except TypeError as e: + raise TypeError(e) + + try: + logger.debug(pretty_long_message(f"Value code: {code}", + _prefix, max_length=_max_length)) + local_vars = {'self': self, 'np': np, 'cp': cp, 'val_map': val_map} + return self._evaluate_expression(code, local_vars) + except Exception as e: + logger.error(f"Error in calculating constr <{self.name}>.\n{e}") + return None + + @property + def v(self): + """ + Return the CVXPY constraint LHS value. + """ + if self.optz is None: + return None + if self.optz._expr.value is None: + try: + shape = self._expr.shape + return np.zeros(shape) + except AttributeError: + return None + else: + return self.optz._expr.value + + @v.setter + def v(self, value): + raise AttributeError("Cannot set the value of the constraint.") diff --git a/ams/opt/exprcalc.py b/ams/opt/exprcalc.py new file mode 100644 index 00000000..430d8fef --- /dev/null +++ b/ams/opt/exprcalc.py @@ -0,0 +1,139 @@ +""" +Module for optimization ExpressionCalc. +""" +import logging + +from typing import Optional +import re + +import numpy as np + +import cvxpy as cp + +from ams.utils import pretty_long_message +from ams.shared import _prefix, _max_length + +from ams.opt import OptzBase, ensure_symbols, ensure_mats_and_parsed + +logger = logging.getLogger(__name__) + + +class ExpressionCalc(OptzBase): + """ + Class for calculating expressions. + + Note that `ExpressionCalc` is not a CVXPY expression, and it should not be used + in the optimization model. + It is used to calculate expression values **after** successful optimization. + """ + + def __init__(self, + name: Optional[str] = None, + info: Optional[str] = None, + unit: Optional[str] = None, + e_str: Optional[str] = None, + model: Optional[str] = None, + src: Optional[str] = None, + ): + OptzBase.__init__(self, name=name, info=info, unit=unit) + self.optz = None + self.e_str = e_str + self.code = None + self.model = model + self.owner = None + self.src = src + self.is_group = False + + def get_idx(self): + if self.is_group: + return self.owner.get_idx() + elif self.owner is None: + logger.info(f'ExpressionCalc <{self.name}> has no owner.') + return None + else: + return self.owner.idx.v + + @ensure_symbols + def parse(self): + """ + Parse the Expression. + """ + # parse the expression str + sub_map = self.om.rtn.syms.sub_map + code_expr = self.e_str + for pattern, replacement in sub_map.items(): + try: + code_expr = re.sub(pattern, replacement, code_expr) + except Exception as e: + raise Exception(f"Error in parsing expr <{self.name}>.\n{e}") + # store the parsed expression str code + self.code = code_expr + msg = f" - ExpressionCalc <{self.name}>: {self.e_str}" + logger.debug(pretty_long_message(msg, _prefix, max_length=_max_length)) + return True + + @ensure_mats_and_parsed + def evaluate(self): + """ + Evaluate the expression. + """ + msg = f" - Expression <{self.name}>: {self.code}" + logger.debug(pretty_long_message(msg, _prefix, max_length=_max_length)) + try: + local_vars = {'self': self, 'np': np, 'cp': cp} + self.optz = self._evaluate_expression(self.code, local_vars=local_vars) + except Exception as e: + raise Exception(f"Error in evaluating ExpressionCalc <{self.name}>.\n{e}") + return True + + def _evaluate_expression(self, code, local_vars=None): + """ + Helper method to evaluate the expression code. + + Parameters + ---------- + code : str + The code string representing the expression. + + Returns + ------- + cp.Expression + The evaluated cvxpy expression. + """ + return eval(code, {}, local_vars) + + @property + def v(self): + """ + Return the CVXPY expression value. + """ + if self.optz is None: + return None + else: + return self.optz.value + + @property + def e(self): + """ + Return the calculated expression value. + """ + if self.code is None: + logger.info(f"ExpressionCalc <{self.name}> is not parsed yet.") + return None + + val_map = self.om.rtn.syms.val_map + code = self.code + for pattern, replacement in val_map.items(): + try: + code = re.sub(pattern, replacement, code) + except TypeError as e: + raise TypeError(e) + + try: + logger.debug(pretty_long_message(f"Value code: {code}", + _prefix, max_length=_max_length)) + local_vars = {'self': self, 'np': np, 'cp': cp, 'val_map': val_map} + return self._evaluate_expression(code, local_vars) + except Exception as e: + logger.error(f"Error in calculating expr <{self.name}>.\n{e}") + return None diff --git a/ams/opt/expression.py b/ams/opt/expression.py new file mode 100644 index 00000000..d7d43c10 --- /dev/null +++ b/ams/opt/expression.py @@ -0,0 +1,200 @@ +""" +Module for optimization Expression. +""" +import logging + +from typing import Optional +import re + +import numpy as np + +import cvxpy as cp + +from ams.utils import pretty_long_message +from ams.shared import _prefix, _max_length + +from ams.opt import OptzBase, ensure_symbols, ensure_mats_and_parsed + +logger = logging.getLogger(__name__) + + +class Expression(OptzBase): + """ + Base class for expressions used in a routine. + + Parameters + ---------- + name : str, optional + Expression name. One should typically assigning the name directly because + it will be automatically assigned by the model. The value of ``name`` + will be the symbol name to be used in expressions. + tex_name : str, optional + LaTeX-formatted variable symbol. Defaults to the value of ``name``. + info : str, optional + Descriptive information + unit : str, optional + Unit + e_str : str, optional + Expression string + model : str, optional + Name of the owner model or group. + src : str, optional + Source expression name + vtype : type, optional + Value type + horizon : ams.routines.RParam, optional + Horizon + """ + + def __init__(self, + name: Optional[str] = None, + tex_name: Optional[str] = None, + info: Optional[str] = None, + unit: Optional[str] = None, + e_str: Optional[str] = None, + model: Optional[str] = None, + src: Optional[str] = None, + vtype: Optional[str] = float, + horizon: Optional[str] = None, + ): + OptzBase.__init__(self, name=name, info=info, unit=unit) + self.tex_name = tex_name + self.e_str = e_str + self.optz = None + self.code = None + self.model = model + self.owner = None + self.src = src + self.is_group = False + self.horizon = horizon + + def get_idx(self): + if self.is_group: + return self.owner.get_idx() + elif self.owner is None: + logger.info(f'ExpressionCalc <{self.name}> has no owner.') + return None + else: + return self.owner.idx.v + + @ensure_symbols + def parse(self): + """ + Parse the expression. + + Returns + ------- + bool + Returns True if the parsing is successful, False otherwise. + """ + sub_map = self.om.rtn.syms.sub_map + code_expr = self.e_str + for pattern, replacement in sub_map.items(): + try: + code_expr = re.sub(pattern, replacement, code_expr) + except Exception as e: + raise Exception(f"Error in parsing expr <{self.name}>.\n{e}") + self.code = code_expr + msg = f" - Expression <{self.name}>: {self.e_str}" + logger.debug(pretty_long_message(msg, _prefix, max_length=_max_length)) + return True + + @ensure_mats_and_parsed + def evaluate(self): + """ + Evaluate the expression. + + Returns + ------- + bool + Returns True if the evaluation is successful, False otherwise. + """ + msg = f" - Expression <{self.name}>: {self.code}" + logger.debug(pretty_long_message(msg, _prefix, max_length=_max_length)) + try: + local_vars = {'self': self, 'np': np, 'cp': cp, 'sub_map': self.om.rtn.syms.val_map} + self.optz = self._evaluate_expression(self.code, local_vars=local_vars) + except Exception as e: + raise Exception(f"Error in evaluating Expression <{self.name}>.\n{e}") + return True + + def _evaluate_expression(self, code, local_vars=None): + """ + Helper method to evaluate the expression code. + + Parameters + ---------- + code : str + The code string representing the expression. + + Returns + ------- + cp.Expression + The evaluated cvxpy expression. + """ + return eval(code, {}, local_vars) + + @property + def v(self): + """ + Return the CVXPY expression value. + """ + if self.optz is None: + return None + else: + return self.optz.value + + @v.setter + def v(self, value): + """ + Set the value. + """ + raise NotImplementedError('Cannot set value to an Expression.') + + @property + def shape(self): + """ + Return the shape. + """ + try: + return self.om.__dict__[self.name].shape + except KeyError: + logger.warning('Shape info is not ready before initialization.') + return None + + @property + def size(self): + """ + Return the size. + """ + if self.rtn.initialized: + return self.om.__dict__[self.name].size + else: + logger.warning(f'Routine <{self.rtn.class_name}> is not initialized yet.') + return None + + @property + def e(self): + """ + Return the calculated expression value. + """ + if self.code is None: + logger.info(f"Expression <{self.name}> is not parsed yet.") + return None + + val_map = self.om.rtn.syms.val_map + code = self.code + for pattern, replacement in val_map.items(): + try: + code = re.sub(pattern, replacement, code) + except TypeError as e: + raise TypeError(e) + + try: + logger.debug(pretty_long_message(f"Value code: {code}", + _prefix, max_length=_max_length)) + local_vars = {'self': self, 'np': np, 'cp': cp, 'val_map': val_map} + return self._evaluate_expression(code, local_vars) + except Exception as e: + logger.error(f"Error in calculating expr <{self.name}>.\n{e}") + return None diff --git a/ams/opt/objective.py b/ams/opt/objective.py new file mode 100644 index 00000000..f9695ec4 --- /dev/null +++ b/ams/opt/objective.py @@ -0,0 +1,174 @@ +""" +Module for optimization Objective. +""" +import logging + +from typing import Optional +import re + +import numpy as np + +import cvxpy as cp + +from ams.utils import pretty_long_message +from ams.shared import _prefix, _max_length + +from ams.opt import OptzBase, ensure_symbols, ensure_mats_and_parsed + +logger = logging.getLogger(__name__) + + +class Objective(OptzBase): + """ + Base class for objective functions. + + This class serves as a template for defining objective functions. Each + instance of this class represents a single objective function that can + be minimized or maximized depending on the sense ('min' or 'max'). + + Parameters + ---------- + name : str, optional + A user-defined name for the objective function. + e_str : str, optional + A mathematical expression representing the objective function. + info : str, optional + Additional informational text about the objective function. + sense : str, optional + The sense of the objective function, default to 'min'. + `min` for minimization and `max` for maximization. + + Attributes + ---------- + v : NoneType + The value of the objective function. It needs to be set through + computation. + rtn : ams.routines.Routine + The owner routine instance. + code : str + The code string for the objective function. + """ + + def __init__(self, + name: Optional[str] = None, + e_str: Optional[str] = None, + info: Optional[str] = None, + unit: Optional[str] = None, + sense: Optional[str] = 'min'): + OptzBase.__init__(self, name=name, info=info, unit=unit) + self.e_str = e_str + self.sense = sense + self.code = None + + @property + def e(self): + """ + Return the calculated objective value. + + Note that `v` should be used primarily as it is obtained + from the solver directly. + + `e` is for debugging purpose. For a successfully solved problem, + `e` should equal to `v`. However, when a problem is infeasible + or unbounded, `e` can be used to check the objective value. + """ + if self.code is None: + logger.info(f"Objective <{self.name}> is not parsed yet.") + return None + + val_map = self.om.rtn.syms.val_map + code = self.code + for pattern, replacement in val_map.items(): + try: + code = re.sub(pattern, replacement, code) + except TypeError as e: + logger.error(f"Error in parsing value for obj <{self.name}>.") + raise e + + try: + logger.debug(pretty_long_message(f"Value code: {code}", + _prefix, max_length=_max_length)) + local_vars = {'self': self, 'np': np, 'cp': cp, 'val_map': val_map} + return self._evaluate_expression(code, local_vars) + except Exception as e: + logger.error(f"Error in calculating obj <{self.name}>.\n{e}") + return None + + @property + def v(self): + """ + Return the CVXPY objective value. + """ + if self.optz is None: + return None + else: + return self.optz.value + + @v.setter + def v(self, value): + raise AttributeError("Cannot set the value of the objective function.") + + @ensure_symbols + def parse(self): + """ + Parse the objective function. + + Returns + ------- + bool + Returns True if the parsing is successful, False otherwise. + """ + # parse the expression str + sub_map = self.om.rtn.syms.sub_map + code_obj = self.e_str + for pattern, replacement, in sub_map.items(): + try: + code_obj = re.sub(pattern, replacement, code_obj) + except Exception as e: + raise Exception(f"Error in parsing obj <{self.name}>.\n{e}") + # store the parsed expression str code + self.code = code_obj + if self.sense not in ['min', 'max']: + raise ValueError(f'Objective sense {self.sense} is not supported.') + sense = 'cp.Minimize' if self.sense == 'min' else 'cp.Maximize' + self.code = f"{sense}({code_obj})" + msg = f" - Objective <{self.name}>: {self.code}" + logger.debug(pretty_long_message(msg, _prefix, max_length=_max_length)) + return True + + @ensure_mats_and_parsed + def evaluate(self): + """ + Evaluate the objective function. + + Returns + ------- + bool + Returns True if the evaluation is successful, False otherwise. + """ + logger.debug(f" - Objective <{self.name}>: {self.e_str}") + try: + local_vars = {'self': self, 'cp': cp} + self.optz = self._evaluate_expression(self.code, local_vars=local_vars) + except Exception as e: + raise Exception(f"Error in evaluating Objective <{self.name}>.\n{e}") + return True + + def _evaluate_expression(self, code, local_vars=None): + """ + Helper method to evaluate the expression code. + + Parameters + ---------- + code : str + The code string representing the expression. + + Returns + ------- + cp.Expression + The evaluated cvxpy expression. + """ + return eval(code, {}, local_vars) + + def __repr__(self): + return f"{self.class_name}: {self.name} [{self.sense.upper()}]" diff --git a/ams/opt/omodel.py b/ams/opt/omodel.py index d81f9143..e0aca8b9 100644 --- a/ams/opt/omodel.py +++ b/ams/opt/omodel.py @@ -1,930 +1,119 @@ """ -Module for optimization modeling. +Module for optimization OModel. """ import logging -from typing import Any, Optional, Union +from typing import Any from collections import OrderedDict -import re -import numpy as np - -from andes.core.common import Config from andes.utils.misc import elapsed import cvxpy as cp -from ams.utils import pretty_long_message -from ams.shared import sps - -logger = logging.getLogger(__name__) - -_prefix = r" - --------------> | " -_max_length = 80 - - -def ensure_symbols(func): - """ - Decorator to ensure that symbols are generated before parsing. - If not, it runs self.rtn.syms.generate_symbols(). - - Designed to be used on the `parse` method of the optimization elements (`OptzBase`) - and optimization model (`OModel`), i.e., `Var`, `Param`, `Constraint`, `Objective`, - and `ExpressionCalc`. - - Note: - ----- - Parsing before symbol generation can give wrong results. Ensure that symbols - are generated before calling the `parse` method. - """ +from ams.opt.optbase import ensure_symbols, ensure_mats_and_parsed - def wrapper(self, *args, **kwargs): - if not self.rtn._syms: - logger.debug(f"<{self.rtn.class_name}> symbols are not generated yet. Generating now...") - self.rtn.syms.generate_symbols() - return func(self, *args, **kwargs) - return wrapper +logger = logging.getLogger(__name__) -def ensure_mats_and_parsed(func): - """ - Decorator to ensure that system matrices are built and OModel is parsed - before evaluation. If not, it runs the necessary methods to initialize them. - - Designed to be used on the `evaluate` method of the optimization elements (`OptzBase`) - and optimization model (`OModel`), i.e., `Var`, `Param`, `Constraint`, `Objective`, - and `ExpressionCalc`. - - Note: - ----- - Evaluation before matrices building and parsing can run into errors. Ensure that - system matrices are built and OModel is parsed before calling the `evaluate` method. - """ - def wrapper(self, *args, **kwargs): - try: - if not self.rtn.system.mats.initialized: - logger.debug("System matrices are not built yet. Building now...") - self.rtn.system.mats.build() - if isinstance(self, (OptzBase, Var, Param, Constraint, Objective)): - if not self.om.parsed: - logger.debug("OModel is not parsed yet. Parsing now...") - self.om.parse() - elif isinstance(self, OModel): - if not self.parsed: - logger.debug("OModel is not parsed yet. Parsing now...") - self.parse() - except Exception as e: - logger.error(f"Error during initialization or parsing: {e}") - raise - return func(self, *args, **kwargs) - return wrapper - - -class OptzBase: +class OModelBase: """ - Base class for optimization elements, e.g., Var and Constraint. - - Parameters - ---------- - name : str, optional - Name. - info : str, optional - Descriptive information - - Attributes - ---------- - rtn : ams.routines.Routine - The owner routine instance. - - Note: - ----- - Ensure that symbols are generated before calling the `parse` method. Parsing - before symbol generation can give wrong results. + Template class for optimization models. """ - def __init__(self, - name: Optional[str] = None, - info: Optional[str] = None, - unit: Optional[str] = None, - ): - self.om = None - self.name = name - self.info = info - self.unit = unit - self.is_disabled = False - self.rtn = None - self.optz = None # corresponding optimization element - self.code = None - - @ensure_symbols - def parse(self): - """ - Parse the object. - """ - raise NotImplementedError - - @ensure_mats_and_parsed - def evaluate(self): - """ - Evaluate the object. - """ - raise NotImplementedError - - @property - def class_name(self): - """ - Return the class name - """ - return self.__class__.__name__ - - @property - def n(self): - """ - Return the number of elements. - """ - if self.owner is None: - return len(self.v) - else: - return self.owner.n - - @property - def shape(self): - """ - Return the shape. - """ - try: - return self.om.__dict__[self.name].shape - except KeyError: - logger.warning('Shape info is not ready before initialization.') - return None + def __init__(self, routine): + self.rtn = routine + self.prob = None + self.exprs = OrderedDict() + self.params = OrderedDict() + self.vars = OrderedDict() + self.constrs = OrderedDict() + self.obj = None + self.parsed = False + self.evaluated = False + self.finalized = False @property - def size(self): + def initialized(self): """ - Return the size. + Return the initialization status. """ - if self.rtn.initialized: - return self.om.__dict__[self.name].size - else: - logger.warning(f'Routine <{self.rtn.class_name}> is not initialized yet.') - return None - - -class ExpressionCalc(OptzBase): - """ - Expression for calculation. - """ - - def __init__(self, - name: Optional[str] = None, - info: Optional[str] = None, - unit: Optional[str] = None, - var: Optional[str] = None, - e_str: Optional[str] = None, - ): - OptzBase.__init__(self, name=name, info=info, unit=unit) - self.optz = None - self.var = var - self.e_str = e_str - self.code = None + return self.parsed and self.evaluated and self.finalized - @ensure_symbols - def parse(self): - """ - Parse the Expression. - """ - # parse the expression str - sub_map = self.om.rtn.syms.sub_map - code_expr = self.e_str - for pattern, replacement in sub_map.items(): - try: - code_expr = re.sub(pattern, replacement, code_expr) - except Exception as e: - raise Exception(f"Error in parsing expr <{self.name}>.\n{e}") - # store the parsed expression str code - self.code = code_expr - msg = f" - ExpressionCalc <{self.name}>: {self.e_str}" - logger.debug(pretty_long_message(msg, _prefix, max_length=_max_length)) - return True + def parse(self, force=False): + self.parsed = True + return self.parsed - @ensure_mats_and_parsed - def evaluate(self): - """ - Evaluate the expression. - """ - msg = f" - Expression <{self.name}>: {self.code}" - logger.debug(pretty_long_message(msg, _prefix, max_length=_max_length)) - try: - local_vars = {'self': self, 'np': np, 'cp': cp} - self.optz = self._evaluate_expression(self.code, local_vars=local_vars) - except Exception as e: - raise Exception(f"Error in evaluating expr <{self.name}>.\n{e}") + def _evaluate_params(self): return True - def _evaluate_expression(self, code, local_vars=None): - """ - Helper method to evaluate the expression code. - - Parameters - ---------- - code : str - The code string representing the expression. - - Returns - ------- - cp.Expression - The evaluated cvxpy expression. - """ - return eval(code, {}, local_vars) - - @property - def v(self): - """ - Return the CVXPY expression value. - """ - if self.optz is None: - return None - else: - return self.optz.value - - def __repr__(self): - return f'{self.__class__.__name__}: {self.name}' - - -class Param(OptzBase): - """ - Base class for parameters used in a routine. - - Parameters - ---------- - no_parse: bool, optional - True to skip parsing the parameter. - nonneg: bool, optional - True to set the parameter as non-negative. - nonpos: bool, optional - True to set the parameter as non-positive. - cplx: bool, optional - True to set the parameter as complex, avoiding the use of `complex`. - imag: bool, optional - True to set the parameter as imaginary. - symmetric: bool, optional - True to set the parameter as symmetric. - diag: bool, optional - True to set the parameter as diagonal. - hermitian: bool, optional - True to set the parameter as hermitian. - boolean: bool, optional - True to set the parameter as boolean. - integer: bool, optional - True to set the parameter as integer. - pos: bool, optional - True to set the parameter as positive. - neg: bool, optional - True to set the parameter as negative. - sparse: bool, optional - True to set the parameter as sparse. - """ - - def __init__(self, - name: Optional[str] = None, - info: Optional[str] = None, - unit: Optional[str] = None, - no_parse: Optional[bool] = False, - nonneg: Optional[bool] = False, - nonpos: Optional[bool] = False, - cplx: Optional[bool] = False, - imag: Optional[bool] = False, - symmetric: Optional[bool] = False, - diag: Optional[bool] = False, - hermitian: Optional[bool] = False, - boolean: Optional[bool] = False, - integer: Optional[bool] = False, - pos: Optional[bool] = False, - neg: Optional[bool] = False, - sparse: Optional[bool] = False, - ): - OptzBase.__init__(self, name=name, info=info, unit=unit) - self.no_parse = no_parse # True to skip parsing the parameter - self.sparse = sparse - - self.config = Config(name=self.class_name) # `config` that can be exported - - self.config.add(OrderedDict((('nonneg', nonneg), - ('nonpos', nonpos), - ('complex', cplx), - ('imag', imag), - ('symmetric', symmetric), - ('diag', diag), - ('hermitian', hermitian), - ('boolean', boolean), - ('integer', integer), - ('pos', pos), - ('neg', neg), - ))) - - @ensure_symbols - def parse(self): - """ - Parse the parameter. - - Returns - ------- - bool - Returns True if the parsing is successful, False otherwise. - """ - sub_map = self.om.rtn.syms.sub_map - code_param = "param(**config)" - for pattern, replacement, in sub_map.items(): - try: - code_param = re.sub(pattern, replacement, code_param) - except Exception as e: - raise Exception(f"Error in parsing param <{self.name}>.\n{e}") - self.code = code_param + def _evaluate_vars(self): return True - @ensure_mats_and_parsed - def evaluate(self): - """ - Evaluate the parameter. - """ - if self.no_parse: - return True - - config = self.config.as_dict() - try: - msg = f"Parameter <{self.name}> is set as sparse, " - msg += "but the value is not sparse." - if self.sparse: - self.v = sps.csr_matrix(self.v) - - # Create the cvxpy.Parameter object - if isinstance(self.v, np.ndarray): - self.optz = cp.Parameter(shape=self.v.shape, **config) - else: - self.optz = cp.Parameter(**config) - self.optz.value = self.v - except ValueError: - msg = f"Parameter <{self.name}> has non-numeric value, " - msg += "set `no_parse=True`." - logger.warning(msg) - self.no_parse = True - return True - except Exception as e: - raise Exception(f"Error in evaluating param <{self.name}>.\n{e}") + def _evaluate_constrs(self): return True - def update(self): - """ - Update the Parameter value. - """ - # NOTE: skip no_parse parameters - if self.optz is None: - return None - self.optz.value = self.v + def _evaluate_obj(self): return True - def __repr__(self): - return f'{self.__class__.__name__}: {self.name}' - - -class Var(OptzBase): - """ - Base class for variables used in a routine. - - When `horizon` is provided, the variable will be expanded to a matrix, - where rows are indexed by the source variable index and columns are - indexed by the horizon index. - - Parameters - ---------- - info : str, optional - Descriptive information - unit : str, optional - Unit - tex_name : str - LaTeX-formatted variable symbol. Defaults to the value of ``name``. - name : str, optional - Variable name. One should typically assigning the name directly because - it will be automatically assigned by the model. The value of ``name`` - will be the symbol name to be used in expressions. - src : str, optional - Source variable name. Defaults to the value of ``name``. - model : str, optional - Name of the owner model or group. - horizon : ams.routines.RParam, optional - Horizon idx. - nonneg : bool, optional - Non-negative variable - nonpos : bool, optional - Non-positive variable - cplx : bool, optional - Complex variable - imag : bool, optional - Imaginary variable - symmetric : bool, optional - Symmetric variable - diag : bool, optional - Diagonal variable - psd : bool, optional - Positive semi-definite variable - nsd : bool, optional - Negative semi-definite variable - hermitian : bool, optional - Hermitian variable - boolean : bool, optional - Boolean variable - integer : bool, optional - Integer variable - pos : bool, optional - Positive variable - neg : bool, optional - Negative variable - - Attributes - ---------- - a : np.ndarray - Variable address. - _v : np.ndarray - Local-storage of the variable value. - rtn : ams.routines.Routine - The owner routine instance. - """ - - def __init__(self, - name: Optional[str] = None, - tex_name: Optional[str] = None, - info: Optional[str] = None, - src: Optional[str] = None, - unit: Optional[str] = None, - model: Optional[str] = None, - shape: Optional[Union[tuple, int]] = None, - v0: Optional[str] = None, - horizon=None, - nonneg: Optional[bool] = False, - nonpos: Optional[bool] = False, - cplx: Optional[bool] = False, - imag: Optional[bool] = False, - symmetric: Optional[bool] = False, - diag: Optional[bool] = False, - psd: Optional[bool] = False, - nsd: Optional[bool] = False, - hermitian: Optional[bool] = False, - boolean: Optional[bool] = False, - integer: Optional[bool] = False, - pos: Optional[bool] = False, - neg: Optional[bool] = False, - ): - self.name = name - self.info = info - self.unit = unit - - self.tex_name = tex_name if tex_name else name - # instance of the owner Model - self.owner = None - # variable internal index inside a model (assigned in run time) - self.id = None - OptzBase.__init__(self, name=name, info=info, unit=unit) - self.src = src - self.is_group = False - self.model = model # indicate if this variable is a group variable - self.owner = None # instance of the owner model or group - self.v0 = v0 - self.horizon = horizon - self._shape = shape - self._v = None - self.a: np.ndarray = np.array([], dtype=int) - - self.config = Config(name=self.class_name) # `config` that can be exported - - self.config.add(OrderedDict((('nonneg', nonneg), - ('nonpos', nonpos), - ('complex', cplx), - ('imag', imag), - ('symmetric', symmetric), - ('diag', diag), - ('psd', psd), - ('nsd', nsd), - ('hermitian', hermitian), - ('boolean', boolean), - ('integer', integer), - ('pos', pos), - ('neg', neg), - ))) - - @property - def v(self): - """ - Return the CVXPY variable value. - """ - if self.optz is None: - return None - if self.optz.value is None: - try: - shape = self.optz.shape - return np.zeros(shape) - except AttributeError: - return None - else: - return self.optz.value - - @v.setter - def v(self, value): - if self.optz is None: - logger.info(f"Variable <{self.name}> is not initialized yet.") - else: - self.optz.value = value - - def get_idx(self): - if self.is_group: - return self.owner.get_idx() - elif self.owner is None: - logger.info(f'Variable <{self.name}> has no owner.') - return None - else: - return self.owner.idx.v - - @ensure_symbols - def parse(self): - """ - Parse the variable. - """ - sub_map = self.om.rtn.syms.sub_map - # NOTE: number of rows is the size of the source variable - if self.owner is not None: - nr = self.owner.n - nc = 0 - if self.horizon: - # NOTE: numer of columns is the horizon if exists - nc = int(self.horizon.n) - shape = (nr, nc) - else: - shape = (nr,) - elif isinstance(self._shape, int): - shape = (self._shape,) - nr = shape - nc = 0 - elif isinstance(self._shape, tuple): - shape = self._shape - nr = shape[0] - nc = shape[1] if len(shape) > 1 else 0 - else: - raise ValueError(f"Invalid shape {self._shape}.") - code_var = f"var({shape}, **config)" - logger.debug(f" - Var <{self.name}>: {self.code}") - for pattern, replacement, in sub_map.items(): - try: - code_var = re.sub(pattern, replacement, code_var) - except Exception as e: - raise Exception(f"Error in parsing var <{self.name}>.\n{e}") - self.code = code_var + def _evaluate_exprs(self): return True - @ensure_mats_and_parsed - def evaluate(self): - """ - Evaluate the variable. - - Returns - ------- - bool - Returns True if the evaluation is successful, False otherwise. - """ - # NOTE: in CVXPY, Config only allow lower case letters - config = {} # used in `self.code` - for k, v in self.config.as_dict().items(): - if k == 'psd': - config['PSD'] = v - elif k == 'nsd': - config['NSD'] = v - elif k == 'bool': - config['boolean'] = v - else: - config[k] = v - msg = f" - Var <{self.name}>: {self.code}" - logger.debug(pretty_long_message(msg, _prefix, max_length=_max_length)) - try: - local_vars = {'self': self, 'config': config, 'cp': cp} - self.optz = eval(self.code, {}, local_vars) - except Exception as e: - raise Exception(f"Error in evaluating var <{self.name}>.\n{e}") + def _evaluate_exprcs(self): return True - def __repr__(self): - return f'{self.__class__.__name__}: {self.owner.__class__.__name__}.{self.name}' - - -class Constraint(OptzBase): - """ - Base class for constraints. - - This class is used as a template for defining constraints. Each - instance of this class represents a single constraint. - - Parameters - ---------- - name : str, optional - A user-defined name for the constraint. - e_str : str, optional - A mathematical expression representing the constraint. - info : str, optional - Additional informational text about the constraint. - is_eq : str, optional - Flag indicating if the constraint is an equality constraint. False indicates - an inequality constraint in the form of `<= 0`. - - Attributes - ---------- - is_disabled : bool - Flag indicating if the constraint is disabled, False by default. - rtn : ams.routines.Routine - The owner routine instance. - is_disabled : bool, optional - Flag indicating if the constraint is disabled, False by default. - dual : float, optional - The dual value of the constraint. - code : str, optional - The code string for the constraint - """ - - def __init__(self, - name: Optional[str] = None, - e_str: Optional[str] = None, - info: Optional[str] = None, - is_eq: Optional[bool] = False, - ): - OptzBase.__init__(self, name=name, info=info) - self.e_str = e_str - self.is_eq = is_eq - self.is_disabled = False - self.dual = None - self.code = None + def evaluate(self, force=False): + self._evaluate_params() + self._evaluate_vars() + self._evaluate_exprs() + self._evaluate_constrs() + self._evaluate_obj() + self._evaluate_exprcs() + self.evaluated = True + return self.evaluated - @ensure_symbols - def parse(self): - """ - Parse the constraint. - """ - # parse the expression str - sub_map = self.om.rtn.syms.sub_map - code_constr = self.e_str - for pattern, replacement in sub_map.items(): - try: - code_constr = re.sub(pattern, replacement, code_constr) - except TypeError as e: - raise TypeError(f"Error in parsing constr <{self.name}>.\n{e}") - # parse the constraint type - code_constr += " == 0" if self.is_eq else " <= 0" - # store the parsed expression str code - self.code = code_constr - msg = f" - Constr <{self.name}>: {self.e_str}" - logger.debug(pretty_long_message(msg, _prefix, max_length=_max_length)) + def finalize(self, force=False): + self.finalized = True return True - def _evaluate_expression(self, code, local_vars=None): - """ - Helper method to evaluate the expression code. - - Parameters - ---------- - code : str - The code string representing the expression. - - Returns - ------- - cp.Expression - The evaluated cvxpy expression. - """ - return eval(code, {}, local_vars) - - @ensure_mats_and_parsed - def evaluate(self): - """ - Evaluate the constraint. - """ - msg = f" - Constr <{self.name}>: {self.code}" - logger.debug(pretty_long_message(msg, _prefix, max_length=_max_length)) - try: - local_vars = {'self': self, 'cp': cp, 'sub_map': self.om.rtn.syms.val_map} - self.optz = self._evaluate_expression(self.code, local_vars=local_vars) - except Exception as e: - raise Exception(f"Error in evaluating constr <{self.name}>.\n{e}") - - def __repr__(self): - enabled = 'OFF' if self.is_disabled else 'ON' - out = f"{self.class_name}: {self.name} [{enabled}]" - return out - - @property - def e(self): - """ - Return the calculated constraint LHS value. - Note that `v` should be used primarily as it is obtained - from the solver directly. - - `e` is for debugging purpose. For a successfully solved problem, - `e` should equal to `v`. However, when a problem is infeasible - or unbounded, `e` can be used to check the constraint LHS value. - """ - if self.code is None: - logger.info(f"Constraint <{self.name}> is not parsed yet.") - return None - - val_map = self.om.rtn.syms.val_map - code = self.code - for pattern, replacement in val_map.items(): - try: - code = re.sub(pattern, replacement, code) - except TypeError as e: - raise TypeError(e) - - try: - logger.debug(pretty_long_message(f"Value code: {code}", - _prefix, max_length=_max_length)) - local_vars = {'self': self, 'np': np, 'cp': cp, 'val_map': val_map} - return self._evaluate_expression(code, local_vars) - except Exception as e: - logger.error(f"Error in calculating constr <{self.name}>.\n{e}") - return None - - @property - def v(self): - """ - Return the CVXPY constraint LHS value. - """ - if self.optz is None: - return None - if self.optz._expr.value is None: - try: - shape = self._expr.shape - return np.zeros(shape) - except AttributeError: - return None - else: - return self.optz._expr.value - - @v.setter - def v(self, value): - raise AttributeError("Cannot set the value of the constraint.") - - -class Objective(OptzBase): - """ - Base class for objective functions. - - This class serves as a template for defining objective functions. Each - instance of this class represents a single objective function that can - be minimized or maximized depending on the sense ('min' or 'max'). - - Parameters - ---------- - name : str, optional - A user-defined name for the objective function. - e_str : str, optional - A mathematical expression representing the objective function. - info : str, optional - Additional informational text about the objective function. - sense : str, optional - The sense of the objective function, default to 'min'. - `min` for minimization and `max` for maximization. - - Attributes - ---------- - v : NoneType - The value of the objective function. It needs to be set through - computation. - rtn : ams.routines.Routine - The owner routine instance. - code : str - The code string for the objective function. - """ - - def __init__(self, - name: Optional[str] = None, - e_str: Optional[str] = None, - info: Optional[str] = None, - unit: Optional[str] = None, - sense: Optional[str] = 'min'): - OptzBase.__init__(self, name=name, info=info, unit=unit) - self.e_str = e_str - self.sense = sense - self.code = None - - @property - def e(self): - """ - Return the calculated objective value. - - Note that `v` should be used primarily as it is obtained - from the solver directly. - - `e` is for debugging purpose. For a successfully solved problem, - `e` should equal to `v`. However, when a problem is infeasible - or unbounded, `e` can be used to check the objective value. - """ - if self.code is None: - logger.info(f"Objective <{self.name}> is not parsed yet.") - return None - - val_map = self.om.rtn.syms.val_map - code = self.code - for pattern, replacement in val_map.items(): - try: - code = re.sub(pattern, replacement, code) - except TypeError as e: - logger.error(f"Error in parsing value for obj <{self.name}>.") - raise e - - try: - logger.debug(pretty_long_message(f"Value code: {code}", - _prefix, max_length=_max_length)) - local_vars = {'self': self, 'np': np, 'cp': cp, 'val_map': val_map} - return self._evaluate_expression(code, local_vars) - except Exception as e: - logger.error(f"Error in calculating obj <{self.name}>.\n{e}") - return None + def init(self, force=False): + self.parse(force) + self.evaluate(force) + self.finalize(force) + return self.initialized @property - def v(self): - """ - Return the CVXPY objective value. - """ - if self.optz is None: - return None - else: - return self.optz.value - - @v.setter - def v(self, value): - raise AttributeError("Cannot set the value of the objective function.") + def class_name(self): + return self.__class__.__name__ - @ensure_symbols - def parse(self): + def _register_attribute(self, key, value): """ - Parse the objective function. + Register a pair of attributes to OModel instance. - Returns - ------- - bool - Returns True if the parsing is successful, False otherwise. + Called within ``__setattr__``, this is where the magic happens. + Subclass attributes are automatically registered based on the variable type. """ - # parse the expression str - sub_map = self.om.rtn.syms.sub_map - code_obj = self.e_str - for pattern, replacement, in sub_map.items(): - try: - code_obj = re.sub(pattern, replacement, code_obj) - except Exception as e: - raise Exception(f"Error in parsing obj <{self.name}>.\n{e}") - # store the parsed expression str code - self.code = code_obj - if self.sense not in ['min', 'max']: - raise ValueError(f'Objective sense {self.sense} is not supported.') - sense = 'cp.Minimize' if self.sense == 'min' else 'cp.Maximize' - self.code = f"{sense}({code_obj})" - msg = f" - Objective <{self.name}>: {self.code}" - logger.debug(pretty_long_message(msg, _prefix, max_length=_max_length)) - return True + if isinstance(value, cp.Variable): + self.vars[key] = value + elif isinstance(value, cp.Constraint): + self.constrs[key] = value + elif isinstance(value, cp.Parameter): + self.params[key] = value + elif isinstance(value, cp.Expression): + self.exprs[key] = value - @ensure_mats_and_parsed - def evaluate(self): - """ - Evaluate the objective function. + def __setattr__(self, name: str, value: Any): + super().__setattr__(name, value) + self._register_attribute(name, value) - Returns - ------- - bool - Returns True if the evaluation is successful, False otherwise. - """ - logger.debug(f" - Objective <{self.name}>: {self.e_str}") - try: - local_vars = {'self': self, 'cp': cp} - self.optz = self._evaluate_expression(self.code, local_vars=local_vars) - except Exception as e: - raise Exception(f"Error in evaluating obj <{self.name}>.\n{e}") + def update(self, params): return True - def _evaluate_expression(self, code, local_vars=None): - """ - Helper method to evaluate the expression code. - - Parameters - ---------- - code : str - The code string representing the expression. - - Returns - ------- - cp.Expression - The evaluated cvxpy expression. - """ - return eval(code, {}, local_vars) - - def __repr__(self): - return f"{self.class_name}: {self.name} [{self.sense.upper()}]" + def __repr__(self) -> str: + return f'{self.rtn.class_name}.{self.__class__.__name__} at {hex(id(self))}' -class OModel: +class OModel(OModelBase): """ Base class for optimization models. @@ -937,12 +126,14 @@ class OModel: ---------- prob: cvxpy.Problem Optimization model. + exprs: OrderedDict + Expressions registry. params: OrderedDict - Parameters. + Parameters registry. vars: OrderedDict - Decision variables. + Decision variables registry. constrs: OrderedDict - Constraints. + Constraints registry. obj: Objective Objective function. initialized: bool @@ -951,25 +142,12 @@ class OModel: Flag indicating if the model is parsed. evaluated: bool Flag indicating if the model is evaluated. + finalized: bool + Flag indicating if the model is finalized. """ def __init__(self, routine): - self.rtn = routine - self.prob = None - self.params = OrderedDict() - self.vars = OrderedDict() - self.constrs = OrderedDict() - self.obj = None - self.parsed = False - self.evaluated = False - self.finalized = False - - @property - def initialized(self): - """ - Return the initialization status. - """ - return self.parsed and self.evaluated and self.finalized + OModelBase.__init__(self, routine) @ensure_symbols def parse(self, force=False): @@ -994,9 +172,14 @@ def parse(self, force=False): if self.parsed and not force: logger.debug("Model is already parsed.") return self.parsed + t, _ = elapsed() - # --- add RParams and Services as parameters --- logger.warning(f'Parsing OModel for <{self.rtn.class_name}>') + # --- add expressions --- + for key, val in self.rtn.exprs.items(): + val.parse() + + # --- add RParams and Services as parameters --- for key, val in self.rtn.params.items(): if not val.no_parse: val.parse() @@ -1009,17 +192,20 @@ def parse(self, force=False): for key, val in self.rtn.constrs.items(): val.parse() + # --- add ExpressionCalcs --- + for key, val in self.rtn.exprcs.items(): + val.parse() + # --- parse objective functions --- - if self.rtn.type != 'PF': - if self.rtn.obj is not None: - try: - self.rtn.obj.parse() - except Exception as e: - raise Exception(f"Failed to parse Objective <{self.rtn.obj.name}>.\n{e}") - else: - logger.warning(f"{self.rtn.class_name} has no objective function!") - self.parsed = False - return self.parsed + if self.rtn.obj is not None: + try: + self.rtn.obj.parse() + except Exception as e: + raise Exception(f"Failed to parse Objective <{self.rtn.obj.name}>.\n{e}") + elif self.rtn.class_name not in ['DCPF0']: + logger.warning(f"{self.rtn.class_name} has no objective function!") + self.parsed = False + return self.parsed # --- parse expressions --- for key, val in self.rtn.exprs.items(): @@ -1072,7 +258,7 @@ def _evaluate_obj(self): """ # NOTE: since we already have the attribute `obj`, # we can update it rather than setting it - if self.rtn.type != 'PF': + if self.rtn.obj is not None: self.rtn.obj.evaluate() self.obj = self.rtn.obj.optz @@ -1081,6 +267,17 @@ def _evaluate_exprs(self): Evaluate the expressions. """ for key, val in self.rtn.exprs.items(): + try: + val.evaluate() + setattr(self, key, val.optz) + except Exception as e: + raise Exception(f"Failed to evaluate Expression <{key}>.\n{e}") + + def _evaluate_exprcs(self): + """ + Evaluate the expressions. + """ + for key, val in self.rtn.exprcs.items(): try: val.evaluate() except Exception as e: @@ -1111,11 +308,13 @@ def evaluate(self, force=False): logger.warning(f"Evaluating OModel for <{self.rtn.class_name}>") t, _ = elapsed() + # NOTE: should evaluate in sequence self._evaluate_params() self._evaluate_vars() + self._evaluate_exprs() self._evaluate_constrs() self._evaluate_obj() - self._evaluate_exprs() + self._evaluate_exprcs() self.evaluated = True _, s = elapsed(t) @@ -1135,7 +334,7 @@ def finalize(self, force=False): Returns True if the finalization is successful, False otherwise. """ # NOTE: for power flow type, we skip the finalization - if self.rtn.type == 'PF': + if self.rtn.class_name in ['DCPF0']: self.finalized = True return self.finalized if self.finalized and not force: @@ -1209,6 +408,8 @@ def _register_attribute(self, key, value): self.constrs[key] = value elif isinstance(value, cp.Parameter): self.params[key] = value + elif isinstance(value, cp.Expression): + self.exprs[key] = value def __setattr__(self, name: str, value: Any): super().__setattr__(name, value) diff --git a/ams/opt/optbase.py b/ams/opt/optbase.py new file mode 100644 index 00000000..56940242 --- /dev/null +++ b/ams/opt/optbase.py @@ -0,0 +1,155 @@ +""" +Module for optimization base classes. +""" +import logging + +from typing import Optional + + +logger = logging.getLogger(__name__) + + +def ensure_symbols(func): + """ + Decorator to ensure that symbols are generated before parsing. + If not, it runs self.rtn.syms.generate_symbols(). + + Designed to be used on the `parse` method of the optimization elements (`OptzBase`) + and optimization model (`OModel`), i.e., `Var`, `Param`, `Constraint`, `Objective`, + and `ExpressionCalc`. + + Parsing before symbol generation can give wrong results. Ensure that symbols + are generated before calling the `parse` method. + """ + + def wrapper(self, *args, **kwargs): + if not self.rtn._syms: + logger.debug(f"<{self.rtn.class_name}> symbols are not generated yet. Generating now...") + self.rtn.syms.generate_symbols() + return func(self, *args, **kwargs) + return wrapper + + +def ensure_mats_and_parsed(func): + """ + Decorator to ensure that system matrices are built and the OModel is parsed + before evaluation. If not, it runs the necessary methods to initialize them. + + Designed to be used on the `evaluate` method of optimization elements (`OptzBase`) + and the optimization model (`OModel`), i.e., `Var`, `Param`, `Constraint`, `Objective`, + and `ExpressionCalc`. + + Evaluation before building matrices and parsing the OModel can lead to errors. Ensure that + system matrices are built and the OModel is parsed before calling the `evaluate` method. + """ + + def wrapper(self, *args, **kwargs): + try: + if not self.rtn.system.mats.initialized: + logger.debug("System matrices are not built yet. Building now...") + self.rtn.system.mats.build() + if isinstance(self, (OptzBase)): + if not self.om.parsed: + logger.debug("OModel is not parsed yet. Parsing now...") + self.om.parse() + else: + if not self.parsed: + logger.debug("OModel is not parsed yet. Parsing now...") + self.parse() + except Exception as e: + logger.error(f"Error during initialization or parsing: {e}") + raise e + return func(self, *args, **kwargs) + return wrapper + + +class OptzBase: + """ + Base class for optimization elements. + Ensure that symbols are generated before calling the `parse` method. Parsing + before symbol generation can lead to incorrect results. + + Parameters + ---------- + name : str, optional + Name of the optimization element. + info : str, optional + Descriptive information about the optimization element. + unit : str, optional + Unit of measurement for the optimization element. + + Attributes + ---------- + rtn : ams.routines.Routine + The owner routine instance. + """ + + def __init__(self, + name: Optional[str] = None, + info: Optional[str] = None, + unit: Optional[str] = None, + ): + self.om = None + self.name = name + self.info = info + self.unit = unit + self.is_disabled = False + self.rtn = None + self.optz = None # corresponding optimization element + self.code = None + + @ensure_symbols + def parse(self): + """ + Parse the object. + """ + raise NotImplementedError + + @ensure_mats_and_parsed + def evaluate(self): + """ + Evaluate the object. + """ + raise NotImplementedError + + @property + def class_name(self): + """ + Return the class name + """ + return self.__class__.__name__ + + @property + def n(self): + """ + Return the number of elements. + """ + if self.owner is None: + return len(self.v) + else: + return self.owner.n + + @property + def shape(self): + """ + Return the shape. + """ + try: + return self.om.__dict__[self.name].shape + except KeyError: + logger.warning('Shape info is not ready before initialization.') + return None + + @property + def size(self): + """ + Return the size. + """ + if self.rtn.initialized: + return self.om.__dict__[self.name].size + else: + logger.warning(f'Routine <{self.rtn.class_name}> is not initialized yet.') + return None + + def __repr__(self): + return f'{self.__class__.__name__}: {self.name}' diff --git a/ams/opt/param.py b/ams/opt/param.py new file mode 100644 index 00000000..f8860fe0 --- /dev/null +++ b/ams/opt/param.py @@ -0,0 +1,156 @@ +""" +Module for optimization Param. +""" +import logging + +from typing import Optional +from collections import OrderedDict +import re + +import numpy as np + +from andes.core.common import Config + +import cvxpy as cp + +from ams.shared import sps + +from ams.opt import OptzBase, ensure_symbols, ensure_mats_and_parsed + +logger = logging.getLogger(__name__) + + +class Param(OptzBase): + """ + Base class for parameters used in a routine. + + Parameters + ---------- + no_parse: bool, optional + True to skip parsing the parameter. + nonneg: bool, optional + True to set the parameter as non-negative. + nonpos: bool, optional + True to set the parameter as non-positive. + cplx: bool, optional + True to set the parameter as complex, avoiding the use of `complex`. + imag: bool, optional + True to set the parameter as imaginary. + symmetric: bool, optional + True to set the parameter as symmetric. + diag: bool, optional + True to set the parameter as diagonal. + hermitian: bool, optional + True to set the parameter as hermitian. + boolean: bool, optional + True to set the parameter as boolean. + integer: bool, optional + True to set the parameter as integer. + pos: bool, optional + True to set the parameter as positive. + neg: bool, optional + True to set the parameter as negative. + sparse: bool, optional + True to set the parameter as sparse. + """ + + def __init__(self, + name: Optional[str] = None, + info: Optional[str] = None, + unit: Optional[str] = None, + no_parse: Optional[bool] = False, + nonneg: Optional[bool] = False, + nonpos: Optional[bool] = False, + cplx: Optional[bool] = False, + imag: Optional[bool] = False, + symmetric: Optional[bool] = False, + diag: Optional[bool] = False, + hermitian: Optional[bool] = False, + boolean: Optional[bool] = False, + integer: Optional[bool] = False, + pos: Optional[bool] = False, + neg: Optional[bool] = False, + sparse: Optional[bool] = False, + ): + OptzBase.__init__(self, name=name, info=info, unit=unit) + self.no_parse = no_parse # True to skip parsing the parameter + self.sparse = sparse + + self.config = Config(name=self.class_name) # `config` that can be exported + + self.config.add(OrderedDict((('nonneg', nonneg), + ('nonpos', nonpos), + ('complex', cplx), + ('imag', imag), + ('symmetric', symmetric), + ('diag', diag), + ('hermitian', hermitian), + ('boolean', boolean), + ('integer', integer), + ('pos', pos), + ('neg', neg), + ))) + + @ensure_symbols + def parse(self): + """ + Parse the parameter. + + Returns + ------- + bool + Returns True if the parsing is successful, False otherwise. + """ + sub_map = self.om.rtn.syms.sub_map + code_param = "param(**config)" + for pattern, replacement, in sub_map.items(): + try: + code_param = re.sub(pattern, replacement, code_param) + except Exception as e: + raise Exception(f"Error in parsing param <{self.name}>.\n{e}") + self.code = code_param + return True + + @ensure_mats_and_parsed + def evaluate(self): + """ + Evaluate the parameter. + """ + if self.no_parse: + return True + + config = self.config.as_dict() + try: + msg = f"Parameter <{self.name}> is set as sparse, " + msg += "but the value is not sparse." + if self.sparse: + self.v = sps.csr_matrix(self.v) + + # Create the cvxpy.Parameter object + if isinstance(self.v, np.ndarray): + self.optz = cp.Parameter(shape=self.v.shape, **config) + else: + self.optz = cp.Parameter(**config) + self.optz.value = self.v + except ValueError: + msg = f"Parameter <{self.name}> has non-numeric value, " + msg += "set `no_parse=True`." + logger.warning(msg) + self.no_parse = True + return True + except Exception as e: + raise Exception(f"Error in evaluating Param <{self.name}>.\n{e}") + return True + + def update(self): + """ + Update the Parameter value. + """ + # NOTE: skip no_parse parameters + if self.optz is None: + return None + self.optz.value = self.v + return True + + def __repr__(self): + return f'{self.__class__.__name__}: {self.name}' diff --git a/ams/opt/var.py b/ams/opt/var.py new file mode 100644 index 00000000..ac13a8cc --- /dev/null +++ b/ams/opt/var.py @@ -0,0 +1,245 @@ +""" +Module for optimization Var. +""" +import logging + +from typing import Optional, Union +from collections import OrderedDict +import re + +import numpy as np + +from andes.core.common import Config + +import cvxpy as cp + +from ams.utils import pretty_long_message +from ams.shared import _prefix, _max_length + +from ams.opt import OptzBase, ensure_symbols, ensure_mats_and_parsed + +logger = logging.getLogger(__name__) + + +class Var(OptzBase): + """ + Base class for variables used in a routine. + + When `horizon` is provided, the variable will be expanded to a matrix, + where rows are indexed by the source variable index and columns are + indexed by the horizon index. + + Parameters + ---------- + info : str, optional + Descriptive information + unit : str, optional + Unit + tex_name : str + LaTeX-formatted variable symbol. Defaults to the value of ``name``. + name : str, optional + Variable name. One should typically assigning the name directly because + it will be automatically assigned by the model. The value of ``name`` + will be the symbol name to be used in expressions. + src : str, optional + Source variable name. Defaults to the value of ``name``. + model : str, optional + Name of the owner model or group. + horizon : ams.routines.RParam, optional + Horizon idx. + nonneg : bool, optional + Non-negative variable + nonpos : bool, optional + Non-positive variable + cplx : bool, optional + Complex variable + imag : bool, optional + Imaginary variable + symmetric : bool, optional + Symmetric variable + diag : bool, optional + Diagonal variable + psd : bool, optional + Positive semi-definite variable + nsd : bool, optional + Negative semi-definite variable + hermitian : bool, optional + Hermitian variable + boolean : bool, optional + Boolean variable + integer : bool, optional + Integer variable + pos : bool, optional + Positive variable + neg : bool, optional + Negative variable + + Attributes + ---------- + a : np.ndarray + Variable address. + _v : np.ndarray + Local-storage of the variable value. + rtn : ams.routines.Routine + The owner routine instance. + """ + + def __init__(self, + name: Optional[str] = None, + tex_name: Optional[str] = None, + info: Optional[str] = None, + src: Optional[str] = None, + unit: Optional[str] = None, + model: Optional[str] = None, + shape: Optional[Union[tuple, int]] = None, + v0: Optional[str] = None, + horizon: Optional[str] = None, + nonneg: Optional[bool] = False, + nonpos: Optional[bool] = False, + cplx: Optional[bool] = False, + imag: Optional[bool] = False, + symmetric: Optional[bool] = False, + diag: Optional[bool] = False, + psd: Optional[bool] = False, + nsd: Optional[bool] = False, + hermitian: Optional[bool] = False, + boolean: Optional[bool] = False, + integer: Optional[bool] = False, + pos: Optional[bool] = False, + neg: Optional[bool] = False, + ): + self.name = name + self.info = info + self.unit = unit + + self.tex_name = tex_name if tex_name else name + # variable internal index inside a model (assigned in run time) + self.id = None + OptzBase.__init__(self, name=name, info=info, unit=unit) + self.src = src + self.is_group = False + self.model = model # indicate if this variable is a group variable + self.owner = None # instance of the owner model or group + self.v0 = v0 + self.horizon = horizon + self._shape = shape + self._v = None + self.a: np.ndarray = np.array([], dtype=int) + + self.config = Config(name=self.class_name) # `config` that can be exported + + self.config.add(OrderedDict((('nonneg', nonneg), + ('nonpos', nonpos), + ('complex', cplx), + ('imag', imag), + ('symmetric', symmetric), + ('diag', diag), + ('psd', psd), + ('nsd', nsd), + ('hermitian', hermitian), + ('boolean', boolean), + ('integer', integer), + ('pos', pos), + ('neg', neg), + ))) + + @property + def v(self): + """ + Return the CVXPY variable value. + """ + if self.optz is None: + return None + if self.optz.value is None: + try: + shape = self.optz.shape + return np.zeros(shape) + except AttributeError: + return None + else: + return self.optz.value + + @v.setter + def v(self, value): + if self.optz is None: + logger.info(f"Variable <{self.name}> is not initialized yet.") + else: + self.optz.value = value + + def get_idx(self): + if self.is_group: + return self.owner.get_idx() + elif self.owner is None: + logger.info(f'Variable <{self.name}> has no owner.') + return None + else: + return self.owner.idx.v + + @ensure_symbols + def parse(self): + """ + Parse the variable. + """ + sub_map = self.om.rtn.syms.sub_map + # NOTE: number of rows is the size of the source variable + if self.owner is not None: + nr = self.owner.n + nc = 0 + if self.horizon: + # NOTE: numer of columns is the horizon if exists + nc = int(self.horizon.n) + shape = (nr, nc) + else: + shape = (nr,) + elif isinstance(self._shape, int): + shape = (self._shape,) + nr = shape + nc = 0 + elif isinstance(self._shape, tuple): + shape = self._shape + nr = shape[0] + nc = shape[1] if len(shape) > 1 else 0 + else: + raise ValueError(f"Invalid shape {self._shape}.") + code_var = f"var({shape}, **config)" + logger.debug(f" - Var <{self.name}>: {self.code}") + for pattern, replacement, in sub_map.items(): + try: + code_var = re.sub(pattern, replacement, code_var) + except Exception as e: + raise Exception(f"Error in parsing var <{self.name}>.\n{e}") + self.code = code_var + return True + + @ensure_mats_and_parsed + def evaluate(self): + """ + Evaluate the variable. + + Returns + ------- + bool + Returns True if the evaluation is successful, False otherwise. + """ + # NOTE: in CVXPY, Config only allow lower case letters + config = {} # used in `self.code` + for k, v in self.config.as_dict().items(): + if k == 'psd': + config['PSD'] = v + elif k == 'nsd': + config['NSD'] = v + elif k == 'bool': + config['boolean'] = v + else: + config[k] = v + msg = f" - Var <{self.name}>: {self.code}" + logger.debug(pretty_long_message(msg, _prefix, max_length=_max_length)) + try: + local_vars = {'self': self, 'config': config, 'cp': cp} + self.optz = eval(self.code, {}, local_vars) + except Exception as e: + raise Exception(f"Error in evaluating Var <{self.name}>.\n{e}") + return True + + def __repr__(self): + return f'{self.__class__.__name__}: {self.owner.__class__.__name__}.{self.name}' diff --git a/ams/pypower/routines/opf.py b/ams/pypower/routines/opf.py index a58c4e37..a855dcf6 100644 --- a/ams/pypower/routines/opf.py +++ b/ams/pypower/routines/opf.py @@ -1203,7 +1203,7 @@ def linear_constraints(self): # nnzA = nnzA + nnz(self.lin["data"].A.(self.lin.order{k})) if self.lin["N"]: - A = sp.sparse.lil_matrix((self.lin["N"], self.var["N"])) + A = sp.lil_matrix((self.lin["N"], self.var["N"])) u = inf * np.ones(self.lin["N"]) l = -u else: diff --git a/ams/report.py b/ams/report.py index ea3e62f6..abc3c0a7 100644 --- a/ams/report.py +++ b/ams/report.py @@ -180,15 +180,26 @@ def write(self): text.append(['']) row_name.append( ['Generation', 'Load']) - if rtn.type == 'ACED': + + if hasattr(rtn, 'pd'): + pd = rtn.pd.v.sum().round(6) + else: + pd = rtn.system.PQ.p0.v.sum().round(6) + if hasattr(rtn, 'qd'): + qd = rtn.qd.v.sum().round(6) + else: + qd = rtn.system.PQ.q0.v.sum().round(6) + + if rtn.type in ['ACED', 'PF']: header.append(['P (p.u.)', 'Q (p.u.)']) - Pcol = [rtn.pg.v.sum().round(6), rtn.pd.v.sum().round(6)] - Qcol = [rtn.qg.v.sum().round(6), rtn.qd.v.sum().round(6)] + Pcol = [rtn.pg.v.sum().round(6), pd] + Qcol = [rtn.qg.v.sum().round(6), qd] data.append([Pcol, Qcol]) else: header.append(['P (p.u.)']) - Pcol = [rtn.pg.v.sum().round(6), rtn.pd.v.sum().round(6)] + Pcol = [rtn.pg.v.sum().round(6), pd] data.append([Pcol]) + # --- routine data --- text.extend(text_sum) header.extend(header_sum) diff --git a/ams/routines/__init__.py b/ams/routines/__init__.py index 87615654..ad5be239 100644 --- a/ams/routines/__init__.py +++ b/ams/routines/__init__.py @@ -8,7 +8,6 @@ all_routines = OrderedDict([ ('dcpf', ['DCPF']), ('pflow', ['PFlow']), - ('pflow2', ['PFlow2']), ('cpf', ['CPF']), ('acopf', ['ACOPF']), ('dcopf', ['DCOPF']), @@ -16,6 +15,8 @@ ('rted', ['RTED', 'RTEDDG', 'RTEDES', 'RTEDVIS']), ('uc', ['UC', 'UCDG', 'UCES']), ('dopf', ['DOPF', 'DOPFVIS']), + ('pflow0', ['PFlow0']), + ('dcpf0', ['DCPF0']), ]) class_names = list_flatten(list(all_routines.values())) diff --git a/ams/routines/acopf.py b/ams/routines/acopf.py index 3a746ffb..e4403493 100644 --- a/ams/routines/acopf.py +++ b/ams/routines/acopf.py @@ -10,13 +10,13 @@ from ams.io.pypower import system2ppc from ams.core.param import RParam -from ams.routines.dcpf import DCPF -from ams.opt.omodel import Var, Constraint, Objective +from ams.routines.dcpf0 import DCPF0 +from ams.opt import Var, Constraint, Objective logger = logging.getLogger(__name__) -class ACOPF(DCPF): +class ACOPF(DCPF0): """ Standard AC optimal power flow. @@ -29,7 +29,7 @@ class ACOPF(DCPF): """ def __init__(self, system, config): - DCPF.__init__(self, system, config) + DCPF0.__init__(self, system, config) self.info = 'AC Optimal Power Flow' self.type = 'ACED' diff --git a/ams/routines/dcopf.py b/ams/routines/dcopf.py index b11a46d0..846ff0dd 100644 --- a/ams/routines/dcopf.py +++ b/ams/routines/dcopf.py @@ -5,17 +5,17 @@ import numpy as np from ams.core.param import RParam -from ams.core.service import NumOp, NumOpDual, VarSelect +from ams.core.service import NumOp, NumOpDual -from ams.routines.routine import RoutineBase +from ams.routines.dcpf import DCPFBase -from ams.opt.omodel import Var, Constraint, Objective, ExpressionCalc +from ams.opt import Constraint, Objective, ExpressionCalc, Expression logger = logging.getLogger(__name__) -class DCOPF(RoutineBase): +class DCOPF(DCPFBase): """ DC optimal power flow (DCOPF). @@ -29,7 +29,7 @@ class DCOPF(RoutineBase): """ def __init__(self, system, config): - RoutineBase.__init__(self, system, config) + DCPFBase.__init__(self, system, config) self.info = 'DC Optimal Power Flow' self.type = 'DCED' @@ -62,10 +62,6 @@ def __init__(self, system, config): indexer='gen', imodel='StaticGen', no_parse=True) # --- generator --- - self.ug = RParam(info='Gen connection status', - name='ug', tex_name=r'u_{g}', - model='StaticGen', src='u', - no_parse=True) self.ctrl = RParam(info='Gen controllability', name='ctrl', tex_name=r'c_{trl}', model='StaticGen', src='ctrl', @@ -90,24 +86,7 @@ def __init__(self, system, config): name='pmin', tex_name=r'p_{g, min}', unit='p.u.', model='StaticGen', no_parse=False,) - self.pg0 = RParam(info='Gen initial active power', - name='pg0', tex_name=r'p_{g, 0}', - unit='p.u.', model='StaticGen', - src='p0', no_parse=False,) - # --- bus --- - self.buss = RParam(info='Bus slack', - name='buss', tex_name=r'B_{us,s}', - model='Slack', src='bus', - no_parse=True,) - # --- load --- - self.upq = RParam(info='Load connection status', - name='upq', tex_name=r'u_{PQ}', - model='StaticLoad', src='u', - no_parse=True,) - self.pd = RParam(info='active demand', - name='pd', tex_name=r'p_{d}', - model='StaticLoad', src='p0', - unit='p.u.',) + # --- line --- self.ul = RParam(info='Line connection status', name='ul', tex_name=r'u_{l}', @@ -125,103 +104,39 @@ def __init__(self, system, config): name='amin', tex_name=r'\theta_{bus, min}', info='min line angle difference', no_parse=True,) - # --- shunt --- - self.gsh = RParam(info='shunt conductance', - name='gsh', tex_name=r'g_{sh}', - model='Shunt', src='g', - no_parse=True,) - # --- connection matrix --- - self.Cg = RParam(info='Gen connection matrix', - name='Cg', tex_name=r'C_{g}', - model='mats', src='Cg', - no_parse=True, sparse=True,) - self.Cl = RParam(info='Load connection matrix', - name='Cl', tex_name=r'C_{l}', - model='mats', src='Cl', - no_parse=True, sparse=True,) - self.CftT = RParam(info='Transpose of line connection matrix', - name='CftT', tex_name=r'C_{ft}^T', - model='mats', src='CftT', - no_parse=True, sparse=True,) - self.Csh = RParam(info='Shunt connection matrix', - name='Csh', tex_name=r'C_{sh}', - model='mats', src='Csh', - no_parse=True, sparse=True,) - # --- system matrix --- - self.Bbus = RParam(info='Bus admittance matrix', - name='Bbus', tex_name=r'B_{bus}', - model='mats', src='Bbus', - no_parse=True, sparse=True,) - self.Bf = RParam(info='Bf matrix', - name='Bf', tex_name=r'B_{f}', - model='mats', src='Bf', - no_parse=True, sparse=True,) - self.Pbusinj = RParam(info='Bus power injection vector', - name='Pbusinj', tex_name=r'P_{bus}^{inj}', - model='mats', src='Pbusinj', - no_parse=True,) - self.Pfinj = RParam(info='Line power injection vector', - name='Pfinj', tex_name=r'P_{f}^{inj}', - model='mats', src='Pfinj', - no_parse=True,) # --- Model Section --- - # --- generation --- - self.pg = Var(info='Gen active power', - unit='p.u.', - name='pg', tex_name=r'p_g', - model='StaticGen', src='p', - v0=self.pg0) - pglb = '-pg + mul(nctrle, pg0) + mul(ctrle, pmin)' + self.pmaxe = Expression(info='Effective pmax', + name='pmaxe', tex_name=r'p_{g, max, e}', + e_str='mul(nctrle, pg0) + mul(ctrle, pmax)', + model='StaticGen', src=None, unit='p.u.',) + self.pmine = Expression(info='Effective pmin', + name='pmine', tex_name=r'p_{g, min, e}', + e_str='mul(nctrle, pg0) + mul(ctrle, pmin)', + model='StaticGen', src=None, unit='p.u.',) self.pglb = Constraint(name='pglb', info='pg min', - e_str=pglb, is_eq=False,) - pgub = 'pg - mul(nctrle, pg0) - mul(ctrle, pmax)' + e_str='-pg + pmine', is_eq=False,) self.pgub = Constraint(name='pgub', info='pg max', - e_str=pgub, is_eq=False,) - # --- bus --- - self.aBus = Var(info='Bus voltage angle', - unit='rad', - name='aBus', tex_name=r'\theta_{bus}', - model='Bus', src='a',) - self.csb = VarSelect(info='select slack bus', - name='csb', tex_name=r'c_{sb}', - u=self.aBus, indexer='buss', - no_parse=True,) - self.sba = Constraint(info='align slack bus angle', - name='sbus', is_eq=True, - e_str='csb@aBus',) - self.pi = Var(info='nodal price', - name='pi', tex_name=r'\pi', - unit='$/p.u.', - model='Bus',) - # --- power balance --- - pb = 'Bbus@aBus + Pbusinj + Cl@(mul(upq, pd)) + Csh@gsh - Cg@pg' - self.pb = Constraint(name='pb', info='power balance', - e_str=pb, is_eq=True,) + e_str='pg - pmaxe', is_eq=False,) + # --- line flow --- - self.plf = Var(info='Line flow', - unit='p.u.', - name='plf', tex_name=r'p_{lf}', - model='Line',) self.plflb = Constraint(info='line flow lower bound', name='plflb', is_eq=False, - e_str='-Bf@aBus - Pfinj - mul(ul, rate_a)',) + e_str='-plf - mul(ul, rate_a)',) self.plfub = Constraint(info='line flow upper bound', name='plfub', is_eq=False, - e_str='Bf@aBus + Pfinj - mul(ul, rate_a)',) + e_str='plf - mul(ul, rate_a)',) self.alflb = Constraint(info='line angle difference lower bound', name='alflb', is_eq=False, e_str='-CftT@aBus + amin',) self.alfub = Constraint(info='line angle difference upper bound', name='alfub', is_eq=False, e_str='CftT@aBus - amax',) - self.plfc = ExpressionCalc(info='plf calculation', - name='plfc', var='plf', - e_str='Bf@aBus + Pfinj') # NOTE: in CVXPY, dual_variables returns a list - self.pic = ExpressionCalc(info='dual of Constraint pb', - name='pic', var='pi', - e_str='pb.dual_variables[0]') + self.pi = ExpressionCalc(info='LMP, dual of ', + name='pi', unit='$/p.u.', + model='Bus', src=None, + e_str='pb.dual_variables[0]') # --- objective --- obj = 'sum(mul(c2, pg**2))' @@ -231,13 +146,6 @@ def __init__(self, system, config): info='total cost', unit='$', sense='min', e_str=obj,) - def solve(self, **kwargs): - """ - Solve the routine optimization model. - args and kwargs go to `self.om.prob.solve()` (`cvxpy.Problem.solve()`). - """ - return self.om.prob.solve(**kwargs) - def run(self, **kwargs): """ Run the routine. @@ -284,98 +192,3 @@ def run(self, **kwargs): A custom solve method to use. """ return super().run(**kwargs) - - def _post_solve(self): - """ - Post-solve calculations. - """ - for expr in self.exprs.values(): - try: - var = getattr(self, expr.var) - var.optz.value = expr.v - logger.debug(f'Post solve: {var} = {expr.e_str}') - except AttributeError: - raise AttributeError(f'No such variable {expr.var}') - return True - - def unpack(self, **kwargs): - """ - Unpack the results from CVXPY model. - """ - # --- solver results to routine algeb --- - for _, var in self.vars.items(): - # --- copy results from routine algeb into system algeb --- - if var.model is None: # if no owner - continue - if var.src is None: # if no source - continue - else: - try: - idx = var.owner.get_idx() - except AttributeError: - idx = var.owner.idx.v - else: - pass - # NOTE: only unpack the variables that are in the model or group - try: - var.owner.set(src=var.src, idx=idx, attr='v', value=var.v) - except (KeyError, TypeError): - logger.error(f'Failed to unpack <{var}> to <{var.owner.class_name}>.') - pass - - # label the most recent solved routine - self.system.recent = self.system.routines[self.class_name] - return True - - def dc2ac(self, kloss=1.0, **kwargs): - """ - Convert the DCOPF results with ACOPF. - - Parameters - ---------- - kloss : float, optional - The loss factor for the conversion. Defaults to 1.2. - """ - exec_time = self.exec_time - if self.exec_time == 0 or self.exit_code != 0: - logger.warning(f'{self.class_name} is not executed successfully, quit conversion.') - return False - - # --- ACOPF --- - # scale up load - pq_idx = self.system.StaticLoad.get_idx() - pd0 = self.system.StaticLoad.get(src='p0', attr='v', idx=pq_idx).copy() - qd0 = self.system.StaticLoad.get(src='q0', attr='v', idx=pq_idx).copy() - self.system.StaticLoad.set(src='p0', idx=pq_idx, attr='v', value=pd0 * kloss) - self.system.StaticLoad.set(src='q0', idx=pq_idx, attr='v', value=qd0 * kloss) - - ACOPF = self.system.ACOPF - # run ACOPF - ACOPF.run() - # self.exec_time += ACOPF.exec_time - # scale load back - self.system.StaticLoad.set(src='p0', idx=pq_idx, value=pd0) - self.system.StaticLoad.set(src='q0', idx=pq_idx, value=qd0) - if not ACOPF.exit_code == 0: - logger.warning(' did not converge, conversion failed.') - # NOTE: mock results to fit interface with ANDES - self.vBus = ACOPF.vBus - self.vBus.optz.value = np.ones(self.system.Bus.n) - self.aBus.optz.value = np.zeros(self.system.Bus.n) - return False - self.pg.optz.value = ACOPF.pg.v - - # NOTE: mock results to fit interface with ANDES - self.addVars(name='vBus', - info='Bus voltage', unit='p.u.', - model='Bus', src='v',) - self.vBus.parse() - self.vBus.optz.value = ACOPF.vBus.v - self.aBus.optz.value = ACOPF.aBus.v - self.exec_time = exec_time - - # --- set status --- - self.system.recent = self - self.converted = True - logger.warning(f'<{self.class_name}> converted to AC.') - return True diff --git a/ams/routines/dcpf.py b/ams/routines/dcpf.py index d0f84735..8462bb9a 100644 --- a/ams/routines/dcpf.py +++ b/ams/routines/dcpf.py @@ -3,186 +3,205 @@ """ import logging -from andes.shared import deg2rad -from andes.utils.misc import elapsed - +from ams.opt import Var, Constraint, Expression, Objective from ams.routines.routine import RoutineBase -from ams.opt.omodel import Var -from ams.pypower import runpf -from ams.pypower.core import ppoption -from ams.io.pypower import system2ppc from ams.core.param import RParam +from ams.core.service import VarSelect logger = logging.getLogger(__name__) -class DCPF(RoutineBase): +class DCPFBase(RoutineBase): """ - DC power flow, overload the ``solve``, ``unpack``, and ``run`` methods. - - Notes - ----- - 1. DCPF is solved with PYPOWER ``runpf`` function. - 2. DCPF formulation is not complete yet, but this does not affect the - results because the data are passed to PYPOWER for solving. + Base class for DC power flow. """ def __init__(self, system, config): RoutineBase.__init__(self, system, config) - self.info = 'DC Power Flow' - self.type = 'PF' - - # --- routine data --- - self.x = RParam(info="line reactance", - name='x', tex_name='x', - unit='p.u.', - model='Line', src='x',) - self.tap = RParam(info="transformer branch tap ratio", - name='tap', tex_name=r't_{ap}', - model='Line', src='tap', - unit='float',) - self.phi = RParam(info="transformer branch phase shift in rad", - name='phi', tex_name=r'\phi', - model='Line', src='phi', - unit='radian',) + self.ug = RParam(info='Gen connection status', + name='ug', tex_name=r'u_{g}', + model='StaticGen', src='u', + no_parse=True) + self.pg0 = RParam(info='Gen initial active power', + name='pg0', tex_name=r'p_{g, 0}', + unit='p.u.', model='StaticGen', + src='p0', no_parse=False,) + # --- shunt --- + self.gsh = RParam(info='shunt conductance', + name='gsh', tex_name=r'g_{sh}', + model='Shunt', src='g', + no_parse=True,) + + self.buss = RParam(info='Bus slack', + name='buss', tex_name=r'B_{us,s}', + model='Slack', src='bus', + no_parse=True,) # --- load --- - self.pd = RParam(info='active deman', + self.pd = RParam(info='active demand', name='pd', tex_name=r'p_{d}', - unit='p.u.', - model='StaticLoad', src='p0') - # --- gen --- + model='StaticLoad', src='p0', + unit='p.u.',) + + # --- connection matrix --- + self.Cg = RParam(info='Gen connection matrix', + name='Cg', tex_name=r'C_{g}', + model='mats', src='Cg', + no_parse=True, sparse=True,) + self.Cl = RParam(info='Load connection matrix', + name='Cl', tex_name=r'C_{l}', + model='mats', src='Cl', + no_parse=True, sparse=True,) + self.CftT = RParam(info='Transpose of line connection matrix', + name='CftT', tex_name=r'C_{ft}^T', + model='mats', src='CftT', + no_parse=True, sparse=True,) + self.Csh = RParam(info='Shunt connection matrix', + name='Csh', tex_name=r'C_{sh}', + model='mats', src='Csh', + no_parse=True, sparse=True,) + + # --- system matrix --- + self.Bbus = RParam(info='Bus admittance matrix', + name='Bbus', tex_name=r'B_{bus}', + model='mats', src='Bbus', + no_parse=True, sparse=True,) + self.Bf = RParam(info='Bf matrix', + name='Bf', tex_name=r'B_{f}', + model='mats', src='Bf', + no_parse=True, sparse=True,) + self.Pbusinj = RParam(info='Bus power injection vector', + name='Pbusinj', tex_name=r'P_{bus}^{inj}', + model='mats', src='Pbusinj', + no_parse=True,) + self.Pfinj = RParam(info='Line power injection vector', + name='Pfinj', tex_name=r'P_{f}^{inj}', + model='mats', src='Pfinj', + no_parse=True,) + + # --- generation --- self.pg = Var(info='Gen active power', unit='p.u.', - name='pg', tex_name=r'p_{g}', - model='StaticGen', src='p',) + name='pg', tex_name=r'p_g', + model='StaticGen', src='p', + v0=self.pg0) # --- bus --- - self.aBus = Var(info='bus voltage angle', + self.aBus = Var(info='Bus voltage angle', unit='rad', - name='aBus', tex_name=r'a_{Bus}', + name='aBus', tex_name=r'\theta_{bus}', model='Bus', src='a',) + # --- power balance --- + pb = 'Bbus@aBus + Pbusinj + Cl@pd + Csh@gsh - Cg@pg' + self.pb = Constraint(name='pb', info='power balance', + e_str=pb, is_eq=True,) + + # --- bus --- + self.csb = VarSelect(info='select slack bus', + name='csb', tex_name=r'c_{sb}', + u=self.aBus, indexer='buss', + no_parse=True,) + self.sba = Constraint(info='align slack bus angle', + name='sbus', is_eq=True, + e_str='csb@aBus',) + # --- line flow --- - self.plf = Var(info='Line flow', - unit='p.u.', - name='plf', tex_name=r'p_{lf}', - model='Line',) + self.plf = Expression(info='Line flow', + name='plf', tex_name=r'p_{lf}', + unit='p.u.', + e_str='Bf@aBus + Pfinj', + model='Line', src=None,) - def unpack(self, res): + def solve(self, **kwargs): """ - Unpack results from PYPOWER. + Solve the routine optimization model. + args and kwargs go to `self.om.prob.solve()` (`cvxpy.Problem.solve()`). """ - system = self.system - mva = res['baseMVA'] - - # --- copy results from ppc into system algeb --- - # --- Bus --- - system.Bus.v.v = res['bus'][:, 7] # voltage magnitude - system.Bus.a.v = res['bus'][:, 8] * deg2rad # voltage angle - - # --- PV --- - system.PV.p.v = res['gen'][system.Slack.n:, 1] / mva # active power - system.PV.q.v = res['gen'][system.Slack.n:, 2] / mva # reactive power + return self.om.prob.solve(**kwargs) - # --- Slack --- - system.Slack.p.v = res['gen'][:system.Slack.n, 1] / mva # active power - system.Slack.q.v = res['gen'][:system.Slack.n, 2] / mva # reactive power - - # --- Line --- - self.plf.optz.value = res['branch'][:, 13] / mva # line flow - - # --- copy results from system algeb into routine algeb --- - for vname, var in self.vars.items(): - owner = getattr(system, var.model) # instance of owner, Model or Group - if var.src is None: # skip if no source variable is specified + def unpack(self, **kwargs): + """ + Unpack the results from CVXPY model. + """ + # --- solver Var results to routine algeb --- + for _, var in self.vars.items(): + # --- copy results from routine algeb into system algeb --- + if var.model is None: # if no owner + continue + if var.src is None: # if no source continue - elif hasattr(owner, 'group'): # if owner is a Model instance - grp = getattr(system, owner.group) - idx = grp.get_idx() - elif hasattr(owner, 'get_idx'): # if owner is a Group instance - idx = owner.get_idx() else: - msg = f"Failed to find valid source variable `{owner.class_name}.{var.src}` for " - msg += f"{self.class_name}.{vname}, skip unpacking." - logger.warning(msg) + try: + idx = var.owner.get_idx() + except AttributeError: + idx = var.owner.idx.v + + # NOTE: only unpack the variables that are in the model or group + try: + var.owner.set(src=var.src, idx=idx, attr='v', value=var.v) + except (KeyError, TypeError): + logger.error(f'Failed to unpack <{var}> to <{var.owner.class_name}>.') + + # --- solver ExpressionCalc results to routine algeb --- + for _, exprc in self.exprcs.items(): + if exprc.model is None: continue - try: - logger.debug(f"Unpacking {vname} into {owner.class_name}.{var.src}.") - var.optz.value = owner.get(src=var.src, attr='v', idx=idx) - except AttributeError: - logger.debug(f"Failed to unpack {vname} into {owner.class_name}.{var.src}.") + if exprc.src is None: continue + else: + try: + idx = exprc.owner.get_idx() + except AttributeError: + idx = exprc.owner.idx.v + else: + pass + try: + exprc.owner.set(src=exprc.src, idx=idx, attr='v', value=exprc.v) + except (KeyError, TypeError): + logger.error(f'Failed to unpack <{exprc}> to <{exprc.owner.class_name}>.') + pass + + # label the most recent solved routine self.system.recent = self.system.routines[self.class_name] return True - def solve(self, method=None): + def _post_solve(self): """ - Solve DC power flow using PYPOWER. + Post-solve calculations. """ - ppc = system2ppc(self.system) - ppopt = ppoption(PF_DC=True) + # NOTE: unpack Expressions if owner and arc are available + for expr in self.exprs.values(): + if expr.owner and expr.src: + expr.owner.set(src=expr.src, attr='v', + idx=expr.get_idx(), value=expr.v) + return True - res, sstats = runpf(casedata=ppc, ppopt=ppopt) - return res, sstats - def run(self, **kwargs): - """ - Run DC pwoer flow. - *args and **kwargs go to `self.solve()`, which are not used yet. - - Examples - -------- - >>> ss = ams.load(ams.get_case('matpower/case14.m')) - >>> ss.DCPF.run() - - Parameters - ---------- - method : str - Placeholder for future use. - - Returns - ------- - exit_code : int - Exit code of the routine. - """ - if not self.initialized: - self.init() - t0, _ = elapsed() - - res, sstats = self.solve(**kwargs) - self.converged = res['success'] - self.exit_code = 0 if res['success'] else 1 - _, s = elapsed(t0) - self.exec_time = float(s.split(' ')[0]) - n_iter = int(sstats['num_iters']) - n_iter_str = f"{n_iter} iterations " if n_iter > 1 else f"{n_iter} iteration " - if self.exit_code == 0: - msg = f"<{self.class_name}> solved in {s}, converged in " - msg += n_iter_str + f"with {sstats['solver_name']}." - logger.info(msg) - try: - self.unpack(res) - except Exception as e: - logger.error(f"Failed to unpack results from {self.class_name}.\n{e}") - return False - return True - else: - msg = f"{self.class_name} failed in " - msg += f"{int(sstats['num_iters'])} iterations with " - msg += f"{sstats['solver_name']}!" - logger.warning(msg) - return False - - def summary(self, **kwargs): - """ - # TODO: Print power flow summary. - """ - raise NotImplementedError +class DCPF(DCPFBase): + """ + DC power flow. + """ - def enable(self, name): - raise NotImplementedError + def __init__(self, system, config): + DCPFBase.__init__(self, system, config) + self.info = 'DC Power Flow' + self.type = 'PF' - def disable(self, name): - raise NotImplementedError + self.genpv = RParam(info='gen of PV', + name='genpv', tex_name=r'g_{DG}', + model='PV', src='idx', + no_parse=True,) + self.cpv = VarSelect(u=self.pg, indexer='genpv', + name='cpv', tex_name=r'C_{PV}', + info='Select PV from pg', + no_parse=True,) + + self.pvb = Constraint(name='pvb', info='PV generator', + e_str='cpv @ (pg - mul(ug, pg0))', + is_eq=True,) + + self.obj = Objective(name='obj', + info='place holder', unit='$', + sense='min', e_str='0',) diff --git a/ams/routines/dcpf0.py b/ams/routines/dcpf0.py new file mode 100644 index 00000000..d3a7403c --- /dev/null +++ b/ams/routines/dcpf0.py @@ -0,0 +1,190 @@ +""" +Power flow routines. +""" +import logging + +from andes.shared import deg2rad +from andes.utils.misc import elapsed + +from ams.routines.routine import RoutineBase +from ams.opt import Var +from ams.pypower import runpf +from ams.pypower.core import ppoption + +from ams.io.pypower import system2ppc +from ams.core.param import RParam + +logger = logging.getLogger(__name__) + + +class DCPF0(RoutineBase): + """ + DC power flow using PYPOWER. + + This class is deprecated as of version 0.9.12 and will be removed in 1.1.0. + + Notes + ----- + 1. DCPF is solved with PYPOWER ``runpf`` function. + 2. DCPF formulation is not complete yet, but this does not affect the + results because the data are passed to PYPOWER for solving. + """ + + def __init__(self, system, config): + RoutineBase.__init__(self, system, config) + self.info = 'DC Power Flow' + self.type = 'PF' + + # --- routine data --- + self.x = RParam(info="line reactance", + name='x', tex_name='x', + unit='p.u.', + model='Line', src='x',) + self.tap = RParam(info="transformer branch tap ratio", + name='tap', tex_name=r't_{ap}', + model='Line', src='tap', + unit='float',) + self.phi = RParam(info="transformer branch phase shift in rad", + name='phi', tex_name=r'\phi', + model='Line', src='phi', + unit='radian',) + + # --- load --- + self.pd = RParam(info='active deman', + name='pd', tex_name=r'p_{d}', + unit='p.u.', + model='StaticLoad', src='p0') + # --- gen --- + self.pg = Var(info='Gen active power', + unit='p.u.', + name='pg', tex_name=r'p_{g}', + model='StaticGen', src='p',) + + # --- bus --- + self.aBus = Var(info='bus voltage angle', + unit='rad', + name='aBus', tex_name=r'a_{Bus}', + model='Bus', src='a',) + + # --- line flow --- + self.plf = Var(info='Line flow', + unit='p.u.', + name='plf', tex_name=r'p_{lf}', + model='Line',) + + def unpack(self, res): + """ + Unpack results from PYPOWER. + """ + system = self.system + mva = res['baseMVA'] + + # --- copy results from ppc into system algeb --- + # --- Bus --- + system.Bus.v.v = res['bus'][:, 7] # voltage magnitude + system.Bus.a.v = res['bus'][:, 8] * deg2rad # voltage angle + + # --- PV --- + system.PV.p.v = res['gen'][system.Slack.n:, 1] / mva # active power + system.PV.q.v = res['gen'][system.Slack.n:, 2] / mva # reactive power + + # --- Slack --- + system.Slack.p.v = res['gen'][:system.Slack.n, 1] / mva # active power + system.Slack.q.v = res['gen'][:system.Slack.n, 2] / mva # reactive power + + # --- Line --- + self.plf.optz.value = res['branch'][:, 13] / mva # line flow + + # --- copy results from system algeb into routine algeb --- + for vname, var in self.vars.items(): + owner = getattr(system, var.model) # instance of owner, Model or Group + if var.src is None: # skip if no source variable is specified + continue + elif hasattr(owner, 'group'): # if owner is a Model instance + grp = getattr(system, owner.group) + idx = grp.get_idx() + elif hasattr(owner, 'get_idx'): # if owner is a Group instance + idx = owner.get_idx() + else: + msg = f"Failed to find valid source variable `{owner.class_name}.{var.src}` for " + msg += f"{self.class_name}.{vname}, skip unpacking." + logger.warning(msg) + continue + try: + logger.debug(f"Unpacking {vname} into {owner.class_name}.{var.src}.") + var.optz.value = owner.get(src=var.src, attr='v', idx=idx) + except AttributeError: + logger.debug(f"Failed to unpack {vname} into {owner.class_name}.{var.src}.") + continue + self.system.recent = self.system.routines[self.class_name] + return True + + def solve(self, method=None): + """ + Solve DC power flow using PYPOWER. + """ + ppc = system2ppc(self.system) + ppopt = ppoption(PF_DC=True) + + res, sstats = runpf(casedata=ppc, ppopt=ppopt) + return res, sstats + + def run(self, **kwargs): + """ + Run DC pwoer flow. + *args and **kwargs go to `self.solve()`, which are not used yet. + + Examples + -------- + >>> ss = ams.load(ams.get_case('matpower/case14.m')) + >>> ss.DCPF.run() + + Parameters + ---------- + method : str + Placeholder for future use. + + Returns + ------- + exit_code : int + Exit code of the routine. + """ + if not self.initialized: + self.init() + t0, _ = elapsed() + + res, sstats = self.solve(**kwargs) + self.converged = res['success'] + self.exit_code = 0 if res['success'] else 1 + _, s = elapsed(t0) + self.exec_time = float(s.split(' ')[0]) + n_iter = int(sstats['num_iters']) + n_iter_str = f"{n_iter} iterations " if n_iter > 1 else f"{n_iter} iteration " + if self.exit_code == 0: + msg = f"<{self.class_name}> solved in {s}, converged in " + msg += n_iter_str + f"with {sstats['solver_name']}." + logger.info(msg) + try: + self.unpack(res) + except Exception as e: + logger.error(f"Failed to unpack results from {self.class_name}.\n{e}") + return False + return True + else: + msg = f"{self.class_name} failed in " + msg += f"{int(sstats['num_iters'])} iterations with " + msg += f"{sstats['solver_name']}!" + logger.warning(msg) + return False + + def summary(self, **kwargs): + """ + # TODO: Print power flow summary. + """ + raise NotImplementedError + + def enable(self, name): + raise NotImplementedError + + def disable(self, name): + raise NotImplementedError diff --git a/ams/routines/dopf.py b/ams/routines/dopf.py index 2d2ef76c..1771d3cb 100644 --- a/ams/routines/dopf.py +++ b/ams/routines/dopf.py @@ -7,7 +7,7 @@ from ams.routines.dcopf import DCOPF -from ams.opt.omodel import Var, Constraint, Objective +from ams.opt import Var, Constraint, Objective class DOPF(DCOPF): diff --git a/ams/routines/ed.py b/ams/routines/ed.py index 24b4b989..944408e4 100644 --- a/ams/routines/ed.py +++ b/ams/routines/ed.py @@ -10,7 +10,7 @@ from ams.routines.rted import RTED, DGBase, ESD1Base -from ams.opt.omodel import Var, Constraint +from ams.opt import Var, Constraint logger = logging.getLogger(__name__) @@ -158,12 +158,12 @@ def __init__(self, system, config): # --- gen --- self.ctrle.u2 = self.ugt self.nctrle.u2 = self.ugt - pglb = '-pg + mul(mul(nctrle, pg0), tlv) ' - pglb += '+ mul(mul(ctrle, tlv), pmin)' - self.pglb.e_str = pglb - pgub = 'pg - mul(mul(nctrle, pg0), tlv) ' - pgub += '- mul(mul(ctrle, tlv), pmax)' - self.pgub.e_str = pgub + pmaxe = 'mul(mul(nctrle, pg0), tlv) + mul(mul(ctrle, tlv), pmax)' + self.pmaxe.e_str = pmaxe + pmine = 'mul(mul(nctrle, pg0), tlv) + mul(mul(ctrle, tlv), pmin)' + self.pmine.e_str = pmine + self.pglb.e_str = '-pg + pmine' + self.pgub.e_str = 'pg - pmaxe' self.pru.horizon = self.timeslot self.pru.info = '2D RegUp power' @@ -182,7 +182,7 @@ def __init__(self, system, config): self.alflb.e_str = '-CftT@aBus + amin@tlv' self.alfub.e_str = 'CftT@aBus - amax@tlv' - self.plfc.e_str = 'Bf@aBus + Pfinj@tlv' + self.plf.e_str = 'Bf@aBus + Pfinj@tlv' # --- power balance --- self.pb.e_str = 'Bbus@aBus + Pbusinj@tlv + Cl@pds + Csh@gsh@tlv - Cg@pg' diff --git a/ams/routines/pflow.py b/ams/routines/pflow.py index 68ec9bd2..3ae30ab1 100644 --- a/ams/routines/pflow.py +++ b/ams/routines/pflow.py @@ -1,111 +1,253 @@ """ -Power flow routines. +Power flow routines independent from PYPOWER. """ import logging +from typing import Optional, Union, Type from collections import OrderedDict -from ams.pypower import runpf +import numpy as np -from ams.io.pypower import system2ppc -from ams.pypower.core import ppoption -from ams.core.param import RParam +from andes.utils.misc import elapsed -from ams.routines.dcpf import DCPF -from ams.opt.omodel import Var +from ams.core.param import RParam +from ams.routines.routine import RoutineBase +from ams.opt import Var, Expression, Objective +from ams.interface import _to_andes_pflow, sync_adsys logger = logging.getLogger(__name__) -class PFlow(DCPF): +class PFlow(RoutineBase): """ - AC Power Flow routine. - - Notes - ----- - 1. AC pwoer flow is solved with PYPOWER ``runpf`` function. - 2. AC power flow formulation in AMS style is NOT DONE YET, - but this does not affect the results - because the data are passed to PYPOWER for solving. + Power flow analysis using ANDES PFlow routine. + + More settings can be changed via ``PFlow2._adsys.config`` and ``PFlow2._adsys.PFlow.config``. + + All generator output powers, bus voltages, and angles are included in the variable definitions. + However, not all of these are unknowns; the definitions are provided for easy access. + + References + ---------- + [1] M. L. Crow, Computational methods for electric power systems. 2015. + + [2] ANDES Documentation - Simulation and Plot. [Online]. Available: + https://docs.andes.app/en/latest/_examples/ex1.html """ def __init__(self, system, config): - DCPF.__init__(self, system, config) - self.info = "AC Power Flow" - self.type = "PF" - - self.config.add(OrderedDict((('qlim', 0), + RoutineBase.__init__(self, system, config) + self.info = 'AC Power Flow' + self.type = 'PF' + self._adsys = None + + self.config.add(OrderedDict((('tol', 1e-6), + ('max_iter', 25), + ('method', 'NR'), + ('check_conn', 1), + ('n_factorize', 4), ))) self.config.add_extra("_help", - qlim="Enforce generator q limits", + tol="convergence tolerance", + max_iter="max. number of iterations", + method="calculation method", + check_conn='check connectivity before power flow', + n_factorize="first N iterations to factorize Jacobian in dishonest method", ) self.config.add_extra("_alt", - qlim=(0, 1, 2), + tol="float", + method=("NR", "dishonest", "NK"), + check_conn=(0, 1), + max_iter=">=10", + n_factorize=">0", ) - self.qd = RParam(info="reactive power load in system base", - name="qd", tex_name=r"q_{d}", - unit="p.u.", - model="StaticLoad", src="q0",) - - # --- bus --- - self.vBus = Var(info="bus voltage magnitude", - unit="p.u.", - name="vBus", tex_name=r"v_{Bus}", - model="Bus", src="v",) - # --- gen --- - self.qg = Var(info="reactive power generation", - unit="p.u.", - name="qg", tex_name=r"q_{g}", - model="StaticGen", src="q",) - # NOTE: omit AC power flow formulation here - - def solve(self, method="newton", **kwargs): - """ - Solve the AC power flow using PYPOWER. - """ - ppc = system2ppc(self.system) - - method_map = dict(newton=1, fdxb=2, fdbx=3, gauss=4) - alg = method_map.get(method) - if alg == 4: - msg = "Gauss method is not fully tested yet, not recommended!" - logger.warning(msg) - if alg is None: - msg = f"Invalid method `{method}` for PFlow." - raise ValueError(msg) - ppopt = ppoption(PF_ALG=alg, ENFORCE_Q_LIMS=self.config.qlim, **kwargs) - - res, sstats = runpf(casedata=ppc, ppopt=ppopt) - return res, sstats + self.Bf = RParam(info='Bf matrix', + name='Bf', tex_name=r'B_{f}', + model='mats', src='Bf', + no_parse=True, sparse=True,) + self.Pfinj = RParam(info='Line power injection vector', + name='Pfinj', tex_name=r'P_{f}^{inj}', + model='mats', src='Pfinj', + no_parse=True,) + + self.pg = Var(info='Gen active power', + unit='p.u.', + name='pg', tex_name=r'p_g', + model='StaticGen', src='p') + self.qg = Var(info='Gen reactive power', + unit='p.u.', + name='qg', tex_name=r'q_g', + model='StaticGen', src='q') + self.aBus = Var(info='Bus voltage angle', + unit='rad', + name='aBus', tex_name=r'\theta_{bus}', + model='Bus', src='a',) + self.vBus = Var(info='Bus voltage magnitude', + unit='p.u.', + name='vBus', tex_name=r'V_{bus}', + model='Bus', src='v',) + self.plf = Expression(info='Line flow', + name='plf', tex_name=r'p_{lf}', + unit='p.u.', + e_str='Bf@aBus + Pfinj', + model='Line', src=None,) + + self.obj = Objective(name='obj', + info='place holder', unit='$', + sense='min', e_str='0',) + + def init(self, **kwargs): + """ + Initialize the ANDES PFlow routine. + + kwargs go to andes.system.System(). + """ + self._adsys = _to_andes_pflow(self.system, + no_output=self.system.files.no_output, + config=self.config.as_dict(), + **kwargs) + self._adsys.setup() + self.om.init() + self.initialized = True + return self.initialized + + def solve(self, **kwargs): + """ + Placeholder. + """ + return True def run(self, **kwargs): """ - Run AC power flow using PYPOWER. + Run the routine. + """ + if not self.initialized: + self.init() - Currently, four methods are supported: 'newton', 'fdxb', 'fdbx', 'gauss', - for Newton's method, fast-decoupled, XB, fast-decoupled, BX, and Gauss-Seidel, - respectively. + t0, _ = elapsed() + _ = self._adsys.PFlow.run() + self.exit_code = self._adsys.exit_code + self.converged = self.exit_code == 0 + _, s = elapsed(t0) + self.exec_time = float(s.split(" ")[0]) - Note that gauss method is not recommended because it seems to be much - more slower than the other three methods and not fully tested yet. + self.unpack() + return True - Examples - -------- - >>> ss = ams.load(ams.get_case('matpower/case14.m')) - >>> ss.PFlow.run() + def _post_solve(self): + """ + Placeholder. + """ + return True - Parameters - ---------- - force_init : bool - Force initialization. - no_code : bool - Disable showing code. - method : str - Method for solving the power flow. + def unpack(self, **kwargs): + """ + Unpack the results from ANDES PFlow routine. + """ + # TODO: maybe also include the DC devices results + sys = self.system + # --- device results --- + bus_idx = sys.Bus.idx.v + sys.Bus.set(src='v', attr='v', idx=bus_idx, + value=self._adsys.Bus.get(src='v', attr='v', idx=bus_idx)) + sys.Bus.set(src='a', attr='v', idx=bus_idx, + value=self._adsys.Bus.get(src='a', attr='v', idx=bus_idx)) + pv_idx = sys.PV.idx.v + pv_u = sys.PV.get(src='u', attr='v', idx=pv_idx) + # NOTE: for p, we should consider the online status as p0 is a param + sys.PV.set(src='p', attr='v', idx=pv_idx, + value=pv_u * sys.PV.get(src='p0', attr='v', idx=pv_idx)) + sys.PV.set(src='q', attr='v', idx=pv_idx, + value=self._adsys.PV.get(src='q', attr='v', idx=pv_idx)) + slack_idx = sys.Slack.idx.v + sys.Slack.set(src='p', attr='v', idx=slack_idx, + value=self._adsys.Slack.get(src='p', attr='v', idx=slack_idx)) + sys.Slack.set(src='q', attr='v', idx=slack_idx, + value=self._adsys.Slack.get(src='q', attr='v', idx=slack_idx)) + # --- routine results --- + self.pg.optz.value = sys.StaticGen.get(src='p', attr='v', idx=self.pg.get_idx()) + self.qg.optz.value = sys.StaticGen.get(src='q', attr='v', idx=self.qg.get_idx()) + self.aBus.optz.value = sys.Bus.get(src='a', attr='v', idx=self.aBus.get_idx()) + self.vBus.optz.value = sys.Bus.get(src='v', attr='v', idx=self.vBus.get_idx()) + return True + + def update(self, **kwargs): + """ + Placeholder. + """ + if not self.initialized: + self.init() + + sync_adsys(self.system, self._adsys) - Returns - ------- - exit_code : int - Exit code of the routine. + return True + + def enable(self, name): + raise NotImplementedError + + def disable(self, name): + raise NotImplementedError + + def addRParam(self, + name: str, + tex_name: Optional[str] = None, + info: Optional[str] = None, + src: Optional[str] = None, + unit: Optional[str] = None, + model: Optional[str] = None, + v: Optional[np.ndarray] = None, + indexer: Optional[str] = None, + imodel: Optional[str] = None,): + """ + Not supported! + """ + raise NotImplementedError + + def addService(self, + name: str, + value: np.ndarray, + tex_name: str = None, + unit: str = None, + info: str = None, + vtype: Type = None,): + """ + Not supported! + """ + raise NotImplementedError + + def addConstrs(self, + name: str, + e_str: str, + info: Optional[str] = None, + is_eq: Optional[str] = False,): + """ + Not supported! + """ + raise NotImplementedError + + def addVars(self, + name: str, + model: Optional[str] = None, + shape: Optional[Union[int, tuple]] = None, + tex_name: Optional[str] = None, + info: Optional[str] = None, + src: Optional[str] = None, + unit: Optional[str] = None, + horizon: Optional[RParam] = None, + nonneg: Optional[bool] = False, + nonpos: Optional[bool] = False, + cplx: Optional[bool] = False, + imag: Optional[bool] = False, + symmetric: Optional[bool] = False, + diag: Optional[bool] = False, + psd: Optional[bool] = False, + nsd: Optional[bool] = False, + hermitian: Optional[bool] = False, + boolean: Optional[bool] = False, + integer: Optional[bool] = False, + pos: Optional[bool] = False, + neg: Optional[bool] = False,): + """ + Not supported! """ - return super().run(**kwargs,) + raise NotImplementedError diff --git a/ams/routines/pflow0.py b/ams/routines/pflow0.py new file mode 100644 index 00000000..74f079fc --- /dev/null +++ b/ams/routines/pflow0.py @@ -0,0 +1,113 @@ +""" +Power flow routines. +""" +import logging +from collections import OrderedDict + +from ams.pypower import runpf + +from ams.io.pypower import system2ppc +from ams.pypower.core import ppoption +from ams.core.param import RParam + +from ams.routines.dcpf0 import DCPF0 +from ams.opt import Var + +logger = logging.getLogger(__name__) + + +class PFlow0(DCPF0): + """ + AC Power Flow using PYPOWER. + + This class is deprecated as of version 0.9.12 and will be removed in 1.1.0. + + Notes + ----- + 1. AC pwoer flow is solved with PYPOWER ``runpf`` function. + 2. AC power flow formulation in AMS style is NOT DONE YET, + but this does not affect the results + because the data are passed to PYPOWER for solving. + """ + + def __init__(self, system, config): + DCPF0.__init__(self, system, config) + self.info = "AC Power Flow" + self.type = "PF" + + self.config.add(OrderedDict((('qlim', 0), + ))) + self.config.add_extra("_help", + qlim="Enforce generator q limits", + ) + self.config.add_extra("_alt", + qlim=(0, 1, 2), + ) + + self.qd = RParam(info="reactive power load in system base", + name="qd", tex_name=r"q_{d}", + unit="p.u.", + model="StaticLoad", src="q0",) + + # --- bus --- + self.vBus = Var(info="bus voltage magnitude", + unit="p.u.", + name="vBus", tex_name=r"v_{Bus}", + model="Bus", src="v",) + # --- gen --- + self.qg = Var(info="reactive power generation", + unit="p.u.", + name="qg", tex_name=r"q_{g}", + model="StaticGen", src="q",) + # NOTE: omit AC power flow formulation here + + def solve(self, method="newton", **kwargs): + """ + Solve the AC power flow using PYPOWER. + """ + ppc = system2ppc(self.system) + + method_map = dict(newton=1, fdxb=2, fdbx=3, gauss=4) + alg = method_map.get(method) + if alg == 4: + msg = "Gauss method is not fully tested yet, not recommended!" + logger.warning(msg) + if alg is None: + msg = f"Invalid method `{method}` for PFlow." + raise ValueError(msg) + ppopt = ppoption(PF_ALG=alg, ENFORCE_Q_LIMS=self.config.qlim, **kwargs) + + res, sstats = runpf(casedata=ppc, ppopt=ppopt) + return res, sstats + + def run(self, **kwargs): + """ + Run AC power flow using PYPOWER. + + Currently, four methods are supported: 'newton', 'fdxb', 'fdbx', 'gauss', + for Newton's method, fast-decoupled, XB, fast-decoupled, BX, and Gauss-Seidel, + respectively. + + Note that gauss method is not recommended because it seems to be much + more slower than the other three methods and not fully tested yet. + + Examples + -------- + >>> ss = ams.load(ams.get_case('matpower/case14.m')) + >>> ss.PFlow.run() + + Parameters + ---------- + force_init : bool + Force initialization. + no_code : bool + Disable showing code. + method : str + Method for solving the power flow. + + Returns + ------- + exit_code : int + Exit code of the routine. + """ + return super().run(**kwargs,) diff --git a/ams/routines/pflow2.py b/ams/routines/pflow2.py deleted file mode 100644 index 8918aad4..00000000 --- a/ams/routines/pflow2.py +++ /dev/null @@ -1,382 +0,0 @@ -""" -Power flow routines independent from PYPOWER. -""" -import logging - -import numpy as np -from ams.core.param import RParam -from ams.core.service import NumOp, NumOpDual, VarSelect - -from ams.routines.routine import RoutineBase - -from ams.opt.omodel import Var, Constraint, Objective - - -logger = logging.getLogger(__name__) - - -class PFlow2(RoutineBase): - """ - Power flow routine. - - References - ---------- - 1. R. D. Zimmerman, C. E. Murillo-Sanchez, and R. J. Thomas, “MATPOWER: Steady-State Operations, Planning, and - Analysis Tools for Power Systems Research and Education,” IEEE Trans. Power Syst., vol. 26, no. 1, pp. 12–19, - Feb. 2011 - """ - - def __init__(self, system, config): - RoutineBase.__init__(self, system, config) - self.info = 'AC Power Flow' - self.type = 'ACED' - - # --- Mapping Section --- - # TODO: skip for now - # # --- from map --- - # self.map1.update({ - # 'ug': ('StaticGen', 'u'), - # }) - # # --- to map --- - # self.map2.update({ - # 'vBus': ('Bus', 'v0'), - # 'ug': ('StaticGen', 'u'), - # 'pg': ('StaticGen', 'p0'), - # }) - - # --- Data Section --- - # --- generator --- - self.ug = RParam(info='Gen connection status', - name='ug', tex_name=r'u_{g}', - model='StaticGen', src='u', - no_parse=True) - self.ctrl = RParam(info='Gen controllability', - name='ctrl', tex_name=r'c_{trl}', - model='StaticGen', src='ctrl', - no_parse=True) - self.ctrle = NumOpDual(info='Effective Gen controllability', - name='ctrle', tex_name=r'c_{trl, e}', - u=self.ctrl, u2=self.ug, - fun=np.multiply, no_parse=True) - self.nctrl = NumOp(u=self.ctrl, fun=np.logical_not, - name='nctrl', tex_name=r'c_{trl,n}', - info='Effective Gen uncontrollability', - no_parse=True,) - self.nctrle = NumOpDual(info='Effective Gen uncontrollability', - name='nctrle', tex_name=r'c_{trl,n,e}', - u=self.nctrl, u2=self.ug, - fun=np.multiply, no_parse=True) - self.pmax = RParam(info='Gen maximum active power', - name='pmax', tex_name=r'p_{g, max}', - unit='p.u.', model='StaticGen', - no_parse=False,) - self.pmin = RParam(info='Gen minimum active power', - name='pmin', tex_name=r'p_{g, min}', - unit='p.u.', model='StaticGen', - no_parse=False,) - self.pg0 = RParam(info='Gen initial active power', - name='pg0', tex_name=r'p_{g, 0}', - unit='p.u.', - model='StaticGen', src='p0') - self.qmax = RParam(info='Gen maximum reactive power', - name='qmax', tex_name=r'q_{g, max}', - unit='p.u.', model='StaticGen', - no_parse=False,) - self.qmin = RParam(info='Gen minimum reactive power', - name='qmin', tex_name=r'q_{g, min}', - unit='p.u.', model='StaticGen', - no_parse=False,) - self.qg0 = RParam(info='Gen initial reactive power', - name='qg0', tex_name=r'q_{g, 0}', - unit='p.u.', - model='StaticGen', src='q0') - # --- bus --- - self.buss = RParam(info='Bus slack', - name='buss', tex_name=r'B_{us,s}', - model='Slack', src='bus', - no_parse=True,) - # --- load --- - self.upq = RParam(info='Load connection status', - name='upq', tex_name=r'u_{PQ}', - model='StaticLoad', src='u', - no_parse=True,) - self.pd = RParam(info='active demand', - name='pd', tex_name=r'p_{d}', - model='StaticLoad', src='p0', - unit='p.u.',) - self.qd = RParam(info='reactive demand', - name='qd', tex_name=r'q_{d}', - model='StaticLoad', src='q0', - unit='p.u.',) - # --- line --- - self.ul = RParam(info='Line connection status', - name='ul', tex_name=r'u_{l}', - model='Line', src='u', - no_parse=True,) - self.rate_a = RParam(info='long-term flow limit', - name='rate_a', tex_name=r'R_{ATEA}', - unit='p.u.', model='Line',) - # --- line angle difference --- - self.amax = RParam(model='Line', src='amax', - name='amax', tex_name=r'\theta_{bus, max}', - info='max line angle difference', - no_parse=True,) - self.amin = RParam(model='Line', src='amin', - name='amin', tex_name=r'\theta_{bus, min}', - info='min line angle difference', - no_parse=True,) - # --- shunt --- - self.gsh = RParam(info='shunt conductance', - name='gsh', tex_name=r'g_{sh}', - model='Shunt', src='g', - no_parse=True,) - # --- connection matrix --- - self.Cg = RParam(info='Gen connection matrix', - name='Cg', tex_name=r'C_{g}', - model='mats', src='Cg', - no_parse=True, sparse=True,) - self.Cl = RParam(info='Load connection matrix', - name='Cl', tex_name=r'C_{l}', - model='mats', src='Cl', - no_parse=True, sparse=True,) - self.CftT = RParam(info='Transpose of line connection matrix', - name='CftT', tex_name=r'C_{ft}^T', - model='mats', src='CftT', - no_parse=True, sparse=True,) - self.Csh = RParam(info='Shunt connection matrix', - name='Csh', tex_name=r'C_{sh}', - model='mats', src='Csh', - no_parse=True, sparse=True,) - # --- system matrix --- - self.Bbus = RParam(info='Bus admittance matrix', - name='Bbus', tex_name=r'B_{bus}', - model='mats', src='Bbus', - no_parse=True, sparse=True,) - self.Bf = RParam(info='Bf matrix', - name='Bf', tex_name=r'B_{f}', - model='mats', src='Bf', - no_parse=True, sparse=True,) - self.Pbusinj = RParam(info='Bus power injection vector', - name='Pbusinj', tex_name=r'P_{bus}^{inj}', - model='mats', src='Pbusinj', - no_parse=True,) - self.Pfinj = RParam(info='Line power injection vector', - name='Pfinj', tex_name=r'P_{f}^{inj}', - model='mats', src='Pfinj', - no_parse=True,) - - # --- Model Section --- - # --- generation --- - self.pg = Var(info='Gen active power', - unit='p.u.', - name='pg', tex_name=r'p_g', - model='StaticGen', src='p', - v0=self.pg0) - pglb = '-pg + mul(nctrle, pg0) + mul(ctrle, pmin)' - self.pglb = Constraint(name='pglb', info='pg min', - e_str=pglb, is_eq=False,) - pgub = 'pg - mul(nctrle, pg0) - mul(ctrle, pmax)' - self.pgub = Constraint(name='pgub', info='pg max', - e_str=pgub, is_eq=False,) - self.qg = Var(info='Gen reactive power', - unit='p.u.', - name='qg', tex_name=r'q_g', - model='StaticGen', src='q',) - qglb = '-qg + mul(nctrle, qg0) + mul(ctrle, qmin)' - self.qglb = Constraint(name='qglb', info='qg min', - e_str=qglb, is_eq=False,) - qgub = 'qg - mul(nctrle, qg0) - mul(ctrle, qmax)' - self.qgub = Constraint(name='qgub', info='qg max', - e_str=qgub, is_eq=False,) - - # --- bus --- - self.aBus = Var(info='Bus voltage angle', - unit='rad', - name='aBus', tex_name=r'\theta_{bus}', - model='Bus', src='a',) - self.vBus = Var(info='Bus voltage magnitude', - unit='p.u.', - name='vBus', tex_name=r'v_{bus}', - model='Bus', src='v',) - self.csb = VarSelect(info='select slack bus', - name='csb', tex_name=r'c_{sb}', - u=self.aBus, indexer='buss', - no_parse=True,) - self.sba = Constraint(info='align slack bus angle', - name='sbus', is_eq=True, - e_str='csb@aBus',) - # --- power balance --- - pb = 'Bbus@aBus + Pbusinj + Cl@(mul(upq, pd)) + Csh@gsh - Cg@pg' - self.pb = Constraint(name='pb', info='power balance', - e_str=pb, is_eq=True,) - # --- line flow --- - self.plf = Var(info='Line flow', - unit='p.u.', - name='plf', tex_name=r'p_{lf}', - model='Line',) - self.plflb = Constraint(info='line flow lower bound', - name='plflb', is_eq=False, - e_str='-Bf@aBus - Pfinj - mul(ul, rate_a)',) - self.plfub = Constraint(info='line flow upper bound', - name='plfub', is_eq=False, - e_str='Bf@aBus + Pfinj - mul(ul, rate_a)',) - self.alflb = Constraint(info='line angle difference lower bound', - name='alflb', is_eq=False, - e_str='-CftT@aBus + amin',) - self.alfub = Constraint(info='line angle difference upper bound', - name='alfub', is_eq=False, - e_str='CftT@aBus - amax',) - - # --- objective --- - self.obj = Objective(name='obj', - info='No objective', unit='$', - sense='min', e_str='0',) - - def solve(self, **kwargs): - """ - Solve the routine optimization model. - args and kwargs go to `self.om.prob.solve()` (`cvxpy.Problem.solve()`). - """ - return self.om.prob.solve(**kwargs) - - def run(self, **kwargs): - """ - Run the routine. - - Following kwargs go to `self.init()`: `force`, `force_mats`, `force_constr`, `force_om`. - - Following kwargs go to `self.solve()`: `solver`, `verbose`, `gp`, `qcp`, `requires_grad`, - `enforce_dpp`, `ignore_dpp`, `method`, and all rest. - - Parameters - ---------- - force : bool, optional - If True, force re-initialization. Defaults to False. - force_mats : bool, optional - If True, force re-generating matrices. Defaults to False. - force_constr : bool, optional - Whether to turn on all constraints. - force_om : bool, optional - If True, force re-generating optimization model. Defaults to False. - solver: str, optional - The solver to use. For example, 'GUROBI', 'ECOS', 'SCS', or 'OSQP'. - verbose : bool, optional - Overrides the default of hiding solver output and prints logging - information describing CVXPY's compilation process. - gp : bool, optional - If True, parses the problem as a disciplined geometric program - instead of a disciplined convex program. - qcp : bool, optional - If True, parses the problem as a disciplined quasiconvex program - instead of a disciplined convex program. - requires_grad : bool, optional - Makes it possible to compute gradients of a solution with respect to Parameters - by calling problem.backward() after solving, or to compute perturbations to the variables - given perturbations to Parameters by calling problem.derivative(). - Gradients are only supported for DCP and DGP problems, not quasiconvex problems. - When computing gradients (i.e., when this argument is True), the problem must satisfy the DPP rules. - enforce_dpp : bool, optional - When True, a DPPError will be thrown when trying to solve a - non-DPP problem (instead of just a warning). - Only relevant for problems involving Parameters. Defaults to False. - ignore_dpp : bool, optional - When True, DPP problems will be treated as non-DPP, which may speed up compilation. Defaults to False. - method : function, optional - A custom solve method to use. - """ - return super().run(**kwargs) - - def _post_solve(self): - """ - Post-solve calculations. - """ - for expr in self.exprs.values(): - try: - var = getattr(self, expr.var) - var.optz.value = expr.v - logger.debug(f'Post solve: {var} = {expr.e_str}') - except AttributeError: - raise AttributeError(f'No such variable {expr.var}') - return True - - def unpack(self, **kwargs): - """ - Unpack the results from CVXPY model. - """ - # --- solver results to routine algeb --- - for _, var in self.vars.items(): - # --- copy results from routine algeb into system algeb --- - if var.model is None: # if no owner - continue - if var.src is None: # if no source - continue - else: - try: - idx = var.owner.get_idx() - except AttributeError: - idx = var.owner.idx.v - else: - pass - # NOTE: only unpack the variables that are in the model or group - try: - var.owner.set(src=var.src, idx=idx, attr='v', value=var.v) - except (KeyError, TypeError): - logger.error(f'Failed to unpack <{var}> to <{var.owner.class_name}>.') - pass - - # label the most recent solved routine - self.system.recent = self.system.routines[self.class_name] - return True - - def dc2ac(self, kloss=1.0, **kwargs): - """ - Convert the DCOPF results with ACOPF. - - Parameters - ---------- - kloss : float, optional - The loss factor for the conversion. Defaults to 1.2. - """ - exec_time = self.exec_time - if self.exec_time == 0 or self.exit_code != 0: - logger.warning(f'{self.class_name} is not executed successfully, quit conversion.') - return False - - # --- ACOPF --- - # scale up load - pq_idx = self.system.StaticLoad.get_idx() - pd0 = self.system.StaticLoad.get(src='p0', attr='v', idx=pq_idx).copy() - qd0 = self.system.StaticLoad.get(src='q0', attr='v', idx=pq_idx).copy() - self.system.StaticLoad.set(src='p0', idx=pq_idx, attr='v', value=pd0 * kloss) - self.system.StaticLoad.set(src='q0', idx=pq_idx, attr='v', value=qd0 * kloss) - - ACOPF = self.system.ACOPF - # run ACOPF - ACOPF.run() - # self.exec_time += ACOPF.exec_time - # scale load back - self.system.StaticLoad.set(src='p0', idx=pq_idx, value=pd0) - self.system.StaticLoad.set(src='q0', idx=pq_idx, value=qd0) - if not ACOPF.exit_code == 0: - logger.warning(' did not converge, conversion failed.') - # NOTE: mock results to fit interface with ANDES - self.vBus = ACOPF.vBus - self.vBus.optz.value = np.ones(self.system.Bus.n) - self.aBus.optz.value = np.zeros(self.system.Bus.n) - return False - self.pg.optz.value = ACOPF.pg.v - - # NOTE: mock results to fit interface with ANDES - self.addVars(name='vBus', - info='Bus voltage', unit='p.u.', - model='Bus', src='v',) - self.vBus.parse() - self.vBus.optz.value = ACOPF.vBus.v - self.aBus.optz.value = ACOPF.aBus.v - self.exec_time = exec_time - - # --- set status --- - self.system.recent = self - self.converted = True - logger.warning(f'<{self.class_name}> converted to AC.') - return True diff --git a/ams/routines/routine.py b/ams/routines/routine.py index 42ee82f5..1abbef33 100644 --- a/ams/routines/routine.py +++ b/ams/routines/routine.py @@ -16,7 +16,8 @@ from ams.core.symprocessor import SymProcessor from ams.core.documenter import RDocumenter from ams.core.service import RBaseService, ValueService -from ams.opt.omodel import OModel, Param, Var, Constraint, Objective, ExpressionCalc +from ams.opt import OModel +from ams.opt import Param, Var, Constraint, Objective, ExpressionCalc, Expression from ams.shared import pd @@ -51,6 +52,8 @@ class RoutineBase: Registry for Var objects. constrs : OrderedDict Registry for Constraint objects. + exprcs : OrderedDict + Registry for ExpressionCalc objects. exprs : OrderedDict Registry for Expression objects. obj : Optional[Objective] @@ -105,6 +108,7 @@ def __init__(self, system=None, config=None): self.params = OrderedDict() # Param registry self.vars = OrderedDict() # Var registry self.constrs = OrderedDict() # Constraint registry + self.exprcs = OrderedDict() # ExpressionCalc registry self.exprs = OrderedDict() # Expression registry self.obj = None # Objective self.initialized = False # initialization flag @@ -544,7 +548,7 @@ def _register_attribute(self, key, value): Called within ``__setattr__``, this is where the magic happens. Subclass attributes are automatically registered based on the variable type. """ - if isinstance(value, (Param, Var, Constraint, Objective, ExpressionCalc)): + if isinstance(value, (Param, Var, Constraint, Objective, ExpressionCalc, Expression)): value.om = self.om value.rtn = self if isinstance(value, Param): @@ -556,8 +560,11 @@ def _register_attribute(self, key, value): elif isinstance(value, Constraint): self.constrs[key] = value self.om.constrs[key] = None # cp.Constraint - elif isinstance(value, ExpressionCalc): + elif isinstance(value, Expression): self.exprs[key] = value + self.om.exprs[key] = None # cp.Expression + elif isinstance(value, ExpressionCalc): + self.exprcs[key] = value elif isinstance(value, RParam): self.rparams[key] = value elif isinstance(value, RBaseService): @@ -807,8 +814,7 @@ def addService(self, tex_name: str = None, unit: str = None, info: str = None, - vtype: Type = None, - model: str = None,): + vtype: Type = None,): """ Add `ValueService` to the routine. @@ -826,8 +832,6 @@ def addService(self, Description. vtype : Type, optional Variable type. - model : str, optional - Model name. """ item = ValueService(name=name, tex_name=tex_name, unit=unit, info=info, @@ -843,8 +847,7 @@ def addConstrs(self, name: str, e_str: str, info: Optional[str] = None, - is_eq: Optional[str] = False, - ): + is_eq: Optional[str] = False,): """ Add `Constraint` to the routine. to the routine. diff --git a/ams/routines/rted.py b/ams/routines/rted.py index c04ce86f..66aa1fe9 100644 --- a/ams/routines/rted.py +++ b/ams/routines/rted.py @@ -9,7 +9,7 @@ from ams.core.service import ZonalSum, VarSelect, NumOp, NumOpDual from ams.routines.dcopf import DCOPF -from ams.opt.omodel import Var, Constraint +from ams.opt import Var, Constraint logger = logging.getLogger(__name__) @@ -161,8 +161,8 @@ def __init__(self, system, config): self.rbu.e_str = 'gs @ mul(ug, pru) - dud' self.rbd.e_str = 'gs @ mul(ug, prd) - ddd' # RegUp/Dn reserve source - self.rru.e_str = 'mul(ug, (pg + pru)) - mul(ug, pmax)' - self.rrd.e_str = 'mul(ug, (-pg + prd)) + mul(ug, pmin)' + self.rru.e_str = 'mul(ug, (pg + pru)) - mul(ug, pmaxe)' + self.rrd.e_str = 'mul(ug, (-pg + prd)) + mul(ug, pmine)' # Gen ramping up/down self.rgu.e_str = 'mul(ug, (pg-pg0-R10))' self.rgd.e_str = 'mul(ug, (-pg+pg0-R10))' @@ -303,7 +303,7 @@ def __init__(self): name='gendg', tex_name=r'g_{DG}', model='DG', src='gen', no_parse=True,) - info = 'Ratio of DG.pge w.r.t to that of static generator', + info = 'Ratio of DG.pge w.r.t to that of static generator' self.gammapdg = RParam(name='gammapd', tex_name=r'\gamma_{p,DG}', model='DG', src='gammap', no_parse=True, info=info) diff --git a/ams/routines/uc.py b/ams/routines/uc.py index f1b29827..f3f49df3 100644 --- a/ams/routines/uc.py +++ b/ams/routines/uc.py @@ -12,7 +12,7 @@ from ams.routines.rted import RTEDBase from ams.routines.ed import SRBase, MPBase, ESD1MPBase, DGBase -from ams.opt.omodel import Var, Constraint +from ams.opt import Var, Constraint logger = logging.getLogger(__name__) @@ -153,12 +153,12 @@ def __init__(self, system, config): self.ctrle.info = 'Reshaped controllability' self.nctrle.u2 = self.tlv self.nctrle.info = 'Reshaped non-controllability' - pglb = '-pg + mul(mul(nctrl, pg0), ugd)' - pglb += '+ mul(mul(ctrl, pmin), ugd)' - self.pglb.e_str = pglb - pgub = 'pg - mul(mul(nctrl, pg0), ugd)' - pgub += '- mul(mul(ctrl, pmax), ugd)' - self.pgub.e_str = pgub + pmaxe = 'mul(mul(nctrl, pg0), ugd) + mul(mul(ctrl, pmax), ugd)' + self.pmaxe.e_str = pmaxe + pmine = 'mul(mul(ctrl, pmin), ugd) + mul(mul(nctrl, pg0), ugd)' + self.pmine.e_str = pmine + self.pglb.e_str = '-pg + pmine' + self.pgub.e_str = 'pg - pmaxe' self.ugd = Var(info='commitment decision', horizon=self.timeslot, diff --git a/ams/shared.py b/ams/shared.py index afe3c92e..92aa5978 100644 --- a/ams/shared.py +++ b/ams/shared.py @@ -23,7 +23,7 @@ pd = LazyImport('import pandas as pd') # --- an empty ANDES system --- -empty_adsys = adSystem() +empty_adsys = adSystem(autogen_stale=False) ad_models = list(empty_adsys.models.keys()) # --- NumPy constants --- @@ -31,6 +31,10 @@ inf = np.inf nan = np.nan +# --- misc constants --- +_prefix = r" - --------------> | " # NOQA +_max_length = 80 # NOQA + # NOTE: copyright year_end = datetime.now().year copyright_msg = f'Copyright (C) 2023-{year_end} Jinning Wang' diff --git a/ams/system.py b/ams/system.py index f14603b9..789655bf 100644 --- a/ams/system.py +++ b/ams/system.py @@ -229,10 +229,13 @@ def import_types(self): def _collect_group_data(self, items): """ - Set the owner for routine attributes: ``RParam``, ``Var``, and ``RBaseService``. + Set the owner for routine attributes: `RParam`, `Var`, `ExpressionCalc`, `Expression`, + and `RBaseService`. """ for item_name, item in items.items(): - if item.model in self.groups.keys(): + if item.model is None: + continue + elif item.model in self.groups.keys(): item.is_group = True item.owner = self.groups[item.model] elif item.model in self.models.keys(): @@ -271,12 +274,18 @@ def import_routines(self): type_instance = self.types[type_name] type_instance.routines[vname] = rtn # self.types[rtn.type].routines[vname] = rtn - # Collect rparams + # Collect RParams rparams = getattr(rtn, 'rparams') self._collect_group_data(rparams) - # Collect routine vars + # Collect routine Vars r_vars = getattr(rtn, 'vars') self._collect_group_data(r_vars) + # Collect ExpressionCalcs + exprc = getattr(rtn, 'exprcs') + self._collect_group_data(exprc) + # Collect Expressions + expr = getattr(rtn, 'exprs') + self._collect_group_data(expr) def import_groups(self): """ diff --git a/docs/source/examples/index.rst b/docs/source/examples/index.rst index fe929b9c..6dacab19 100644 --- a/docs/source/examples/index.rst +++ b/docs/source/examples/index.rst @@ -30,4 +30,5 @@ folder of the repository :caption: Demonstration ../_examples/demo/demo_ESD1.ipynb - ../_examples/demo/demo_AGC.ipynb \ No newline at end of file + ../_examples/demo/demo_AGC.ipynb + ../_examples/demo/demo_debug.ipynb \ No newline at end of file diff --git a/docs/source/release-notes.rst b/docs/source/release-notes.rst index b1ddb161..4f52de57 100644 --- a/docs/source/release-notes.rst +++ b/docs/source/release-notes.rst @@ -12,9 +12,17 @@ Pre-v1.0.0 v0.9.12 (202x-xx-xx) -------------------- -- Refactor OModel.initialized as a property method -- Add a demo to show using `Constraint.e` for debugging -- Fix `OModel.Param.evaluate()` when its value is a number +- Refactor ``OModel.initialized`` as a property method +- Add a demo to show using ``Constraint.e`` for debugging +- Fix ``ams.opt.omodel.Param.evaluate`` when its value is a number +- Improve ``ams.opt.omodel.ExpressionCalc`` for better performance +- Refactor module ``ams.opt`` +- Add class ``ams.opt.Expression`` +- Switch from PYPOWER to ANDES in routine ``PFlow`` +- Switch from PYPOWER to regular formulation in routine ``DCPF`` +- Refactor routines ``DCPF`` and ``DCOPF`` +- In ``RDocumenter``, set Srouce to be owner if there is no src +- Specify ``multiprocess<=0.70.16`` in requirements as 0.70.17 does not support Linux RC1 ~~~~ @@ -24,49 +32,51 @@ v0.9.11 (2024-11-14) -------------------- - Add pyproject.toml for PEP 517 and PEP 518 compliance -- Add model `Jumper` -- Fix deprecation warning related to `pandas.fillna` and `newshape` in NumPy -- Minor refactor on solvers information in the module `shared` +- Add model ``Jumper`` +- Fix deprecation warning related to ``pandas.fillna`` and ``newshape`` in NumPy +- Minor refactor on solvers information in the module ``shared`` - Change default values of minimum ON/OFF duration time of generators to be 1 and 0.5 hours -- Add parameter `uf` for enforced generator on/off status -- In servicee `LoadScale`, consider load online status -- Consider line online status in routine `ED` -- Add methods `evaluate` and `finalize` in the class `OModel` to handle optimization elements generation and assembling -- Refactor `OModel.init()` and `Routine.init()` -- Add ANDES paper as a citation file for now +- Add parameter ``uf`` for enforced generator on/off status +- In servicee ``LoadScale``, consider load online status +- Consider line online status in routine ``ED`` +- Add methods ``evaluate`` and ``finalize`` in the class ``OModel`` to handle optimization + elements generation and assembling +- Refactor ``OModel.init`` and ``Routine.init`` +- Add ANDES paper as the citation file for now - Add more routine tests for generator trip, line trip, and load trip - Add a README to overview built-in cases -- Rename methods `v2` as `e` for classes `Constraint` and `Objective` +- Rename methods ``v2`` as ``e`` for classes ``Constraint`` and ``Objective`` - Add benchmark functions -- Improve using of `eval()` in module `omodel` -- Refactor module `interop.andes` as module `interface` for simplicity +- Improve the usage of ``eval`` in module ``omodel`` +- Refactor module ``interop.andes`` as module ``interface`` for simplicity v0.9.10 (2024-09-03) -------------------- Hotfix of import issue in ``v0.9.9``. -- In module `MatProcessor`, add two parameters `permc_spec` and `use_umfpack` in function `build_ptdf` +- In module ``MatProcessor``, add two parameters ``permc_spec`` and ``use_umfpack`` in function ``build_ptdf`` - Follow RTD's deprecation of Sphinx context injection at build time - In MATPOWER conversion, set devices name as None - Skip macOS tests in azure-pipelines due to failure in fixing its configuration - Prepare to support NumPy v2.0.0, but solvers have unexpected behavior -- Improve the logic of setting `Optz` value +- Improve the logic of setting ``Optz`` value - Support NumPy v2.0.0 v0.9.9 (2024-09-02) ------------------- -**IMPORTANT NOTICE: This version has known issues and should be avoided.** +**NOTICE: This version has known issues and has been yanked on PyPI.** v0.9.8 (2024-06-18) ------------------- -- Assign `MParam.owner` when declaring -- In `MatProcessor`, improve `build_ptdf` and `build_lodf` to allow partial building and incremental building -- Add in 'cases/matpower/Benchmark.json' for benchmark with MATPOWER +- Assign ``MParam.owner`` when declaring +- In ``MatProcessor``, improve ``build_ptdf`` and ``build_lodf`` to allow partial building and + incremental building +- Add file ``cases/matpower/Benchmark.json`` for benchmark with MATPOWER - Improve known good results test -- Minor fix in `main.py` selftest part +- Minor fix in ``main.py`` selftest part - Set dependency NumPy version to be <2.0.0 to avoid CVXPY compatibility issues v0.9.7 (2024-05-24) @@ -81,30 +91,31 @@ References: Frequency Regulation," in IEEE Transactions on Smart Grid, doi: 10.1109/TSG.2024.3356948. - Fix OTDF calculation -- Add parameter `dtype='float64'` and `no_store=False` in `MatProcessor` PTDF, LODF, and OTDF calculation, to save memory -- Add placeholder parameter `Bus.type` +- Add parameter ``dtype='float64'`` and ``no_store=False`` in ``MatProcessor`` PTDF, LODF, and OTDF + calculation, to save memory +- Add placeholder parameter ``Bus.type`` v0.9.6 (2024-04-21) ------------------- -This patch release refactor and improve `MatProcessor`, where it support PTDF, LODF, +This patch release refactor and improve ``MatProcessor``, where it support PTDF, LODF, and OTDF for static analysis. The reference can be found online "PowerWorld > Web Help > Sensitivities > Line Outage Distribution Factors". - Refactor DCPF, PFlow, and ACOPF -- Add a loss factor in ``RTED.dc2ac()`` -- Add ``DCOPF.dc2ac()`` +- Add a loss factor in ``RTED.dc2ac`` +- Add ``DCOPF.dc2ac`` - Fix OModel parse status to ensure no_parsed params can be updated - Fix and rerun ex2 -- Format ``Routine.get()`` return type to be consistent with input idx type -- Remove unused ``Routine.prepare()`` -- Refactor `MatProcessor` to separate matrix building -- Add Var `plf` in `DCPF`, `PFlow`, and `ACOPF` to store the line flow -- Add `build_ptdf`, `build_lodf`, and `build_otdf` -- Fix ``Routine.get()`` to support pd.Series type idx input -- Reserve `exec_time` after ``dc2ac()`` +- Format ``Routine.get`` return type to be consistent with input idx type +- Remove unused ``Routine.prepare`` +- Refactor ``MatProcessor`` to separate matrix building +- Add Var ``plf`` in ``DCPF``, ``PFlow``, and ``ACOPF`` to store the line flow +- Add ``build_ptdf``, ``build_lodf``, and ``build_otdf`` +- Fix ``Routine.get`` to support pd.Series type idx input +- Reserve ``exec_time`` after ``dc2ac`` - Adjust kloss to fix ex2 v0.9.5 (2024-03-25) @@ -112,10 +123,10 @@ v0.9.5 (2024-03-25) - Add more plots in demo_AGC - Improve line rating adjustment -- Adjust static import sequence in `models.__init__.py` +- Adjust static import sequence in ``models.__init__.py`` - Adjust pjm5bus case line rate_a - Fix formulation of constraint line angle diff -- Align slack bus angle to zero in `DCOPF` +- Align slack bus angle to zero in ``DCOPF`` - Align StaticGen idx sequence with converted MATPOWER case - Fix several issues in MATPOWER converter @@ -130,7 +141,7 @@ v0.9.3 (2024-03-06) ------------------- - Major improvemets on demo_AGC -- Bug fix in ``RTED.dc2ac()`` +- Bug fix in ``RTED.dc2ac`` v0.9.2 (2024-03-04) ------------------- @@ -156,7 +167,7 @@ v0.9.0 (2024-02-27) - Fix ``addService``, ``addVars`` - Rename ``RoutineModel`` to ``RoutineBase`` for better naming - Fix ANDES file converter issue -- Initial release to conda-forge +- Initial release on conda-forge v0.8.5 (2024-01-31) ------------------- @@ -172,7 +183,7 @@ v0.8.4 (2024-01-30) v0.8.3 (2024-01-30) ------------------- -- Initial release to PyPI +- Initial release on PyPI v0.8.2 (2024-01-30) ------------------- diff --git a/examples/demonstration/demo_debug.ipynb b/examples/demonstration/demo_debug.ipynb index 0ca2948a..6f144393 100644 --- a/examples/demonstration/demo_debug.ipynb +++ b/examples/demonstration/demo_debug.ipynb @@ -4,9 +4,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Using LHS Value to Debug Unsolved Problems\n", + "# Using LHS Value to Debug An Infeasible Problem\n", "\n", - "This notebook aims to show how to use the left-hand side (LHS) value (`Constraint.e`) to debug unsolved problems.\n", + "This notebook aims to show how to use the left-hand side (LHS) value (`Constraint.e`) to debug an infeasible problem.\n", "\n", "The LHS value is the value of the left-hand side of the equation. It is useful to check which constraints are violated." ] @@ -140,7 +140,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Then, we can lower down the generator maximum output to create an unsolvable problem." + "Then, we can lower down the generator maximum output to create an infeasible problem." ] }, { diff --git a/examples/ex6.ipynb b/examples/ex6.ipynb index 0a68d364..b150d398 100644 --- a/examples/ex6.ipynb +++ b/examples/ex6.ipynb @@ -5,14 +5,14 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Multi-period Dispatch Simulation" + "# Multi-period Scheduling Simulation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Multi-period dispatch economic dispatch (ED) and unit commitment (UC) is also available.\n", + "Multi-period economic dispatch (ED) and unit commitment (UC) are also available.\n", "\n", "In this case, we will show a 24-hour ED simulation." ] diff --git a/requirements.txt b/requirements.txt index 5ceafe9f..6c6fe056 100644 --- a/requirements.txt +++ b/requirements.txt @@ -10,3 +10,4 @@ openpyxl andes>=1.8.7 pybind11 cvxpy +multiprocess<=0.70.16 diff --git a/tests/test_1st_system.py b/tests/test_1st_system.py index d2cc20a1..79f5fa39 100644 --- a/tests/test_1st_system.py +++ b/tests/test_1st_system.py @@ -26,7 +26,7 @@ def test_docum(self) -> None: for export in ['plain', 'rest']: docum._obj_doc(max_width=78, export=export) docum._constr_doc(max_width=78, export=export) - docum._expr_doc(max_width=78, export=export) + docum._exprc_doc(max_width=78, export=export) docum._var_doc(max_width=78, export=export) docum._service_doc(max_width=78, export=export) docum._param_doc(max_width=78, export=export) diff --git a/tests/test_cli.py b/tests/test_cli.py index 5c164d04..ad57826b 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -1,9 +1,18 @@ import unittest +import os import ams.main +import ams.cli class TestCLI(unittest.TestCase): + + def test_cli_parser(self): + ams.cli.create_parser() + + def test_cli_preamble(self): + ams.cli.preamble() + def test_main_doc(self): ams.main.doc('Bus') ams.main.doc(list_supported=True) @@ -14,3 +23,12 @@ def test_versioninfo(self): def test_misc(self): ams.main.misc(show_license=True) ams.main.misc(save_config=None, overwrite=True) + + def test_profile_run(self): + _ = ams.main.run(ams.get_case('matpower/case5.m'), + no_output=False, + profile=True,) + self.assertTrue(os.path.exists('case5_prof.prof')) + self.assertTrue(os.path.exists('case5_prof.txt')) + os.remove('case5_prof.prof') + os.remove('case5_prof.txt') diff --git a/tests/test_interop.py b/tests/test_interop.py index 11070f21..5dd8e709 100644 --- a/tests/test_interop.py +++ b/tests/test_interop.py @@ -127,9 +127,7 @@ def test_verify_pf(self): sa = to_andes(sp, setup=True, no_output=True, default_config=True, verify=False, tol=1e-3) - # NOTE: it is known that there is 1e-7~1e-6 diff in case300.m - self.assertFalse(verify_pf(amsys=sp, adsys=sa, tol=1e-6)) - self.assertTrue(verify_pf(amsys=sp, adsys=sa, tol=1e-3)) + self.assertTrue(verify_pf(amsys=sp, adsys=sa, tol=1e-7)) class TestDataExchange(unittest.TestCase): diff --git a/tests/test_known_good.py b/tests/test_known_good.py index 7f72d894..e53d1839 100644 --- a/tests/test_known_good.py +++ b/tests/test_known_good.py @@ -164,8 +164,10 @@ def test_DCPF_case118(self): Test DC power flow for case118. """ self.sp.DCPF.run() - np.testing.assert_allclose(self.sp.DCPF.aBus.v * rad2deg, - np.array(self.mpres['case118']['DCPF']['aBus']).reshape(-1), + aBus_mp = np.array(self.mpres['case118']['DCPF']['aBus']).reshape(-1) + aBus_mp -= aBus_mp[0] + np.testing.assert_allclose((self.sp.DCPF.aBus.v - self.sp.DCPF.aBus.v[0]) * rad2deg, + aBus_mp, rtol=1e-2, atol=1e-2) np.testing.assert_allclose(self.sp.DCPF.pg.v.sum() * self.sp.config.mva, diff --git a/tests/test_omodel.py b/tests/test_omodel.py index 48d24cad..97929d1c 100644 --- a/tests/test_omodel.py +++ b/tests/test_omodel.py @@ -62,3 +62,58 @@ def test_uninitialized_after_nonparametric_update(self): self.assertTrue(self.ss.DCOPF.om.initialized, "OModel should be initialized after initialization!") self.ss.DCOPF.update('ug') self.assertTrue(self.ss.DCOPF.om.initialized, "OModel should be initialized after nonparametric update!") + + +class TestOModelrepr(unittest.TestCase): + """ + Test __repr__ in module `omodel`. + """ + + def setUp(self) -> None: + self.ss = ams.load(ams.get_case("5bus/pjm5bus_demo.xlsx"), + setup=True, + default_config=True, + no_output=True, + ) + + def test_dcopf(self): + """ + Test `DCOPF.om` __repr__(). + """ + self.ss.DCOPF.om.parse() + self.assertIsInstance(self.ss.DCOPF.pg.__repr__(), str) + self.assertIsInstance(self.ss.DCOPF.pmax.__repr__(), str) + self.assertIsInstance(self.ss.DCOPF.pb.__repr__(), str) + self.assertIsInstance(self.ss.DCOPF.plflb.__repr__(), str) + self.assertIsInstance(self.ss.DCOPF.obj.__repr__(), str) + + def test_rted(self): + """ + Test `RTED.om` __repr__(). + """ + self.ss.RTED.om.parse() + self.assertIsInstance(self.ss.RTED.pru.__repr__(), str) + self.assertIsInstance(self.ss.RTED.R10.__repr__(), str) + self.assertIsInstance(self.ss.RTED.pru.__repr__(), str) + self.assertIsInstance(self.ss.RTED.rru.__repr__(), str) + self.assertIsInstance(self.ss.RTED.obj.__repr__(), str) + + def test_ed(self): + """ + Test `ED.om` __repr__(). + """ + self.ss.ED.om.parse() + self.assertIsInstance(self.ss.ED.prs.__repr__(), str) + self.assertIsInstance(self.ss.ED.R30.__repr__(), str) + self.assertIsInstance(self.ss.ED.prsb.__repr__(), str) + self.assertIsInstance(self.ss.ED.obj.__repr__(), str) + + def test_uc(self): + """ + Test `UC.om` __repr__(). + """ + self.ss.UC.om.parse() + self.assertIsInstance(self.ss.UC.prns.__repr__(), str) + self.assertIsInstance(self.ss.UC.cdp.__repr__(), str) + self.assertIsInstance(self.ss.UC.prnsb.__repr__(), str) + self.assertIsInstance(self.ss.UC.obj.__repr__(), str) diff --git a/tests/test_rtn_rted.py b/tests/test_rtn_rted.py index f112ada7..8dcd650d 100644 --- a/tests/test_rtn_rted.py +++ b/tests/test_rtn_rted.py @@ -81,7 +81,7 @@ def test_set_load(self): def test_dc2ac(self): """ - Test `RTED.init()` method. + Test `RTED.dc2ac()` method. """ self.ss.RTED.dc2ac() self.assertTrue(self.ss.RTED.converted, "AC conversion failed!")