diff --git a/.github/ISSUE_TEMPLATE/bug_report.md b/.github/ISSUE_TEMPLATE/bug_report.md new file mode 100644 index 00000000..2c7a7bd1 --- /dev/null +++ b/.github/ISSUE_TEMPLATE/bug_report.md @@ -0,0 +1,39 @@ +--- +name: Bug report +about: Create a report to help us improve +title: '' +labels: bug +assignees: cuihantao + +--- + +**Describe the bug** + +A clear and concise description of what the bug is. + +**To Reproduce** + +Update data files and code for reproducing the error. + +If you need to insert code inline, use a pair of triple backticks to format: + + ```python + Code goes here + ``` + +**Expected behavior** + +A clear and concise description of what you expected to happen. + +**Desktop (please complete the following information):** + + - OS: [e.g. Windows] + - AMS version (please paste the output from `ams misc --version`) + + +**pip packages (please paste the output from `pip list`)** + + +**Additional context** + +Add any other context about the problem here. diff --git a/.github/ISSUE_TEMPLATE/feature_request.md b/.github/ISSUE_TEMPLATE/feature_request.md new file mode 100644 index 00000000..f12c71ce --- /dev/null +++ b/.github/ISSUE_TEMPLATE/feature_request.md @@ -0,0 +1,20 @@ +--- +name: Feature request +about: Suggest an idea for this project +title: '' +labels: '' +assignees: + +--- + +**Is your feature request related to a problem? Please describe.** +A clear and concise description of what the problem is. Ex. I'm always frustrated when [...] + +**Describe the solution you'd like** +A clear and concise description of what you want to happen. + +**Describe alternatives you've considered** +A clear and concise description of any alternative solutions or features you've considered. + +**Additional context** +Add any other context or screenshots about the feature request here. diff --git a/.github/PULL_REQUEST_TEMPLATE/pull_request_template.md b/.github/PULL_REQUEST_TEMPLATE/pull_request_template.md new file mode 100644 index 00000000..7a792a6b --- /dev/null +++ b/.github/PULL_REQUEST_TEMPLATE/pull_request_template.md @@ -0,0 +1,8 @@ +Fixes # + +## Proposed Changes + + - + - + - + \ No newline at end of file diff --git a/.github/workflows/pythonapp.yml b/.github/workflows/pythonapp.yml new file mode 100644 index 00000000..a01fa9e3 --- /dev/null +++ b/.github/workflows/pythonapp.yml @@ -0,0 +1,43 @@ +name: Python application + +on: [push, pull_request] + +jobs: + build: + + runs-on: ubuntu-latest + + steps: + - uses: actions/checkout@v1 + - name: Set up Python 3.9 + uses: actions/setup-python@v1 + with: + python-version: 3.9 + - name: Install dependencies + run: | + # work around a pip issue with extras_require: https://github.com/pypa/pip/issues/8323 + python -m pip install -U "pip @ git+https://github.com/pypa/pip.git" + python -m pip install --upgrade pip + python -m pip install .[all] + python -m pip install nbmake==0.10 pytest-xdist line_profiler # add'l packages for notebook tests. + # - name: Lint with flake8 for pull requests + # if: github.event_name == 'pull_request' + # run: | + # # stop the build if there are Python syntax errors or undefined names + # flake8 . + - name: Test with pytest + run: | + pytest --log-cli-level=10 + # - name: Test notebooks. + # run: | + # pytest --log-cli-level=10 --nbmake examples --ignore=examples/verification + # - name: Build a distribution if tagged + # if: github.event_name == 'push' && startsWith(github.ref, 'refs/tags') + # run: | + # python setup.py sdist + # - name: Publish a Python distribution to PyPI if tagged + # if: github.event_name == 'push' && startsWith(github.ref, 'refs/tags') + # uses: pypa/gh-action-pypi-publish@master + # with: + # user: __token__ + # password: ${{ secrets.pypi_password }} diff --git a/.gitignore b/.gitignore index ecc52bc9..3efdd590 100644 --- a/.gitignore +++ b/.gitignore @@ -35,4 +35,7 @@ _generated icebar # example excel files cache -~$*.xlsx \ No newline at end of file +~$*.xlsx + +# conda +.conda/* \ No newline at end of file diff --git a/ams/cases/5bus/pjm5bus_uced.json b/ams/cases/5bus/pjm5bus_uced.json index 6b7baeb4..2f997f79 100644 --- a/ams/cases/5bus/pjm5bus_uced.json +++ b/ams/cases/5bus/pjm5bus_uced.json @@ -722,145 +722,169 @@ "idx": "EDT1", "u": 1.0, "name": "EDT1", - "sd": "0.793,0.0" + "sd": "0.793,0.0", + "ug": "1,1,1,1" }, { "idx": "EDT2", "u": 1.0, "name": "EDT2", - "sd": "0.756,0.0" + "sd": "0.756,0.0", + "ug": "1,1,1,1" }, { "idx": "EDT3", "u": 1.0, "name": "EDT3", - "sd": "0.723,0.0" + "sd": "0.723,0.0", + "ug": "1,1,1,1" }, { "idx": "EDT4", "u": 1.0, "name": "EDT4", - "sd": "0.708,0.0" + "sd": "0.708,0.0", + "ug": "1,1,1,1" }, { "idx": "EDT5", "u": 1.0, "name": "EDT5", - "sd": "0.7,0.0" + "sd": "0.7,0.0", + "ug": "1,1,1,1" }, { "idx": "EDT6", "u": 1.0, "name": "EDT6", - "sd": "0.706,0.0" + "sd": "0.706,0.0", + "ug": "1,1,1,1" }, { "idx": "EDT7", "u": 1.0, "name": "EDT7", - "sd": "0.75,0.0" + "sd": "0.75,0.0", + "ug": "1,1,1,1" }, { "idx": "EDT8", "u": 1.0, "name": "EDT8", - "sd": "0.802,0.0" + "sd": "0.802,0.0", + "ug": "1,1,1,1" }, { "idx": "EDT9", "u": 1.0, "name": "EDT9", - "sd": "0.828,0.0" + "sd": "0.828,0.0", + "ug": "1,1,1,1" }, { "idx": "EDT10", "u": 1.0, "name": "EDT10", - "sd": "0.851,0.0" + "sd": "0.851,0.0", + "ug": "1,1,1,1" }, { "idx": "EDT11", "u": 1.0, "name": "EDT11", - "sd": "0.874,0.0" + "sd": "0.874,0.0", + "ug": "1,1,1,1" }, { "idx": "EDT12", "u": 1.0, "name": "EDT12", - "sd": "0.898,0.0" + "sd": "0.898,0.0", + "ug": "1,1,1,1" }, { "idx": "EDT13", "u": 1.0, "name": "EDT13", - "sd": "0.919,0.0" + "sd": "0.919,0.0", + "ug": "1,1,1,1" }, { "idx": "EDT14", "u": 1.0, "name": "EDT14", - "sd": "0.947,0.0" + "sd": "0.947,0.0", + "ug": "1,1,1,1" }, { "idx": "EDT15", "u": 1.0, "name": "EDT15", - "sd": "0.97,0.0" + "sd": "0.97,0.0", + "ug": "1,1,1,1" }, { "idx": "EDT16", "u": 1.0, "name": "EDT16", - "sd": "0.987,0.0" + "sd": "0.987,0.0", + "ug": "1,1,1,1" }, { "idx": "EDT17", "u": 1.0, "name": "EDT17", - "sd": "1.0,0.0" + "sd": "1.0,0.0", + "ug": "1,1,1,1" }, { "idx": "EDT18", "u": 1.0, "name": "EDT18", - "sd": "1.0,0.0" + "sd": "1.0,0.0", + "ug": "1,1,1,1" }, { "idx": "EDT19", "u": 1.0, "name": "EDT19", - "sd": "0.991,0.0" + "sd": "0.991,0.0", + "ug": "1,1,1,1" }, { "idx": "EDT20", "u": 1.0, "name": "EDT20", - "sd": "0.956,0.0" + "sd": "0.956,0.0", + "ug": "1,1,1,1" }, { "idx": "EDT21", "u": 1.0, "name": "EDT21", - "sd": "0.93,0.0" + "sd": "0.93,0.0", + "ug": "1,1,1,1" }, { "idx": "EDT22", "u": 1.0, "name": "EDT22", - "sd": "0.905,0.0" + "sd": "0.905,0.0", + "ug": "1,1,1,1" }, { "idx": "EDT23", "u": 1.0, "name": "EDT23", - "sd": "0.849,0.0" + "sd": "0.849,0.0", + "ug": "1,1,1,1" }, { "idx": "EDT24", "u": 1.0, "name": "EDT24", - "sd": "0.784,0.0" + "sd": "0.784,0.0", + "ug": "1,1,1,1" } ], "UCTSlot": [ diff --git a/ams/cases/5bus/pjm5bus_uced.xlsx b/ams/cases/5bus/pjm5bus_uced.xlsx index ccc2bb56..2a9446ca 100644 Binary files a/ams/cases/5bus/pjm5bus_uced.xlsx and b/ams/cases/5bus/pjm5bus_uced.xlsx differ diff --git a/ams/cases/5bus/pjm5bus_uced_esd1.xlsx b/ams/cases/5bus/pjm5bus_uced_esd1.xlsx new file mode 100644 index 00000000..b0a1926d Binary files /dev/null and b/ams/cases/5bus/pjm5bus_uced_esd1.xlsx differ diff --git a/ams/cases/ieee123/ieee123_regcv1.xlsx b/ams/cases/ieee123/ieee123_regcv1.xlsx index d6237c0c..8bf6cfec 100644 Binary files a/ams/cases/ieee123/ieee123_regcv1.xlsx and b/ams/cases/ieee123/ieee123_regcv1.xlsx differ diff --git a/ams/cases/ieee14/ieee14_uced.xlsx b/ams/cases/ieee14/ieee14_uced.xlsx index 1986b6a5..dc7a723b 100644 Binary files a/ams/cases/ieee14/ieee14_uced.xlsx and b/ams/cases/ieee14/ieee14_uced.xlsx differ diff --git a/ams/cases/ieee39/ieee39_esd1.xlsx b/ams/cases/ieee39/ieee39_esd1.xlsx deleted file mode 100644 index c324dd45..00000000 Binary files a/ams/cases/ieee39/ieee39_esd1.xlsx and /dev/null differ diff --git a/ams/cases/ieee39/ieee39_uced.xlsx b/ams/cases/ieee39/ieee39_uced.xlsx index 36f33e82..1db5c73d 100644 Binary files a/ams/cases/ieee39/ieee39_uced.xlsx and b/ams/cases/ieee39/ieee39_uced.xlsx differ diff --git a/ams/cases/ieee39/ieee39_uced_esd1.xlsx b/ams/cases/ieee39/ieee39_uced_esd1.xlsx index 6555130c..eb3cb8d9 100644 Binary files a/ams/cases/ieee39/ieee39_uced_esd1.xlsx and b/ams/cases/ieee39/ieee39_uced_esd1.xlsx differ diff --git a/ams/cases/npcc/npcc_uced.xlsx b/ams/cases/npcc/npcc_uced.xlsx index c879f42f..58d63a6b 100644 Binary files a/ams/cases/npcc/npcc_uced.xlsx and b/ams/cases/npcc/npcc_uced.xlsx differ diff --git a/ams/cases/wecc/wecc_uced.xlsx b/ams/cases/wecc/wecc_uced.xlsx index 49069e3b..cf3e6a29 100644 Binary files a/ams/cases/wecc/wecc_uced.xlsx and b/ams/cases/wecc/wecc_uced.xlsx differ diff --git a/ams/core/documenter.py b/ams/core/documenter.py index 20f46716..c8f4b2e2 100644 --- a/ams/core/documenter.py +++ b/ams/core/documenter.py @@ -498,23 +498,32 @@ def _tex_pre(docm, p, tex_map): """ # NOTE: in the future, there might occur special math symbols - special_map = OrderedDict([ + map_before = OrderedDict([ + ('sum', 'SUM'), + (r'\sum', 'SUM'), + (r'\eta', 'ETA'), + (r'\gamma', 'GAMMA'), + (r'\theta', 'THETA'), + (r'\frac', 'FRAC'), + (r'\overline', 'OVERLINE'), + (r'\underline', 'UNDERLINE'), + ]) + map_post = OrderedDict([ ('SUM', r'\sum'), + ('THETA', r'\theta'), ('ETA', r'\eta'), + ('GAMMA', r'\gamma'), ('FRAC', r'\frac'), + ('OVERLINE', r'\overline'), + ('UNDERLINE', r'\underline'), ]) expr = p.e_str for pattern, replacement in tex_map.items(): - if r'\sum' in replacement: - replacement = replacement.replace(r'\sum', 'SUM') - if r'sum' in expr: - expr = expr.replace('sum', 'SUM') - if '\eta' in replacement: - replacement = replacement.replace('\eta', 'ETA') - if r'\frac' in replacement: - replacement = replacement.replace(r'\frac', 'FRAC') + for key, val in map_before.items(): + if key in replacement: + replacement = replacement.replace(key, val) if r'\p' in replacement: continue try: @@ -529,7 +538,7 @@ def _tex_pre(docm, p, tex_map): except re.error: logger.error('Remains '*' in the expression.') - for pattern, replacement in special_map.items(): + for pattern, replacement in map_post.items(): expr = expr.replace(pattern, replacement) return expr diff --git a/ams/core/matprocessor.py b/ams/core/matprocessor.py index 43d75631..a94b7fa5 100644 --- a/ams/core/matprocessor.py +++ b/ams/core/matprocessor.py @@ -2,17 +2,23 @@ Module for system matrix make. """ -import logging # NOQA -from typing import Optional # NOQA +import logging +from typing import Optional -import numpy as np # NOQA -from ams.pypower.make import makePTDF, makeBdc # NOQA -from ams.io.pypower import system2ppc # NOQA +import numpy as np + +from scipy.sparse import csr_matrix as c_sparse +from scipy.sparse import lil_matrix as l_sparse + +from ams.pypower.make import makePTDF, makeBdc +from ams.io.pypower import system2ppc + +from ams.opt.omodel import Param logger = logging.getLogger(__name__) -class MParam: +class MParam(Param): """ Class for matrix parameters built from the system. @@ -44,6 +50,7 @@ def __init__(self, unit: Optional[str] = None, v: Optional[np.ndarray] = None, ): + Param.__init__(self, name=name, info=info) self.name = name self.tex_name = tex_name if (tex_name is not None) else name self.info = info @@ -56,6 +63,14 @@ def v(self): """ Return the value of the parameter. """ + # NOTE: scipy.sparse matrix will return 2D array + # so we squeeze it here if only one row + if isinstance(self._v, (c_sparse, l_sparse)): + out = self._v.toarray() + if out.shape[0] == 1: + return np.squeeze(out) + else: + return out return self._v @property @@ -102,6 +117,9 @@ def __init__(self, system): self.Cg = MParam(name='Cg', tex_name=r'C_g', info='Generator connectivity matrix', v=None) + self.Cs = MParam(name='Cs', tex_name=r'C_s', + info='Slack connectivity matrix', + v=None) self.Cl = MParam(name='Cl', tex_name=r'Cl', info='Load connectivity matrix', v=None) @@ -122,18 +140,22 @@ def make(self): # FIXME: hard coded here gen_bus = system.StaticGen.get(src='bus', attr='v', idx=system.StaticGen.get_idx()) + slack_bus = system.Slack.get(src='bus', attr='v', + idx=system.Slack.idx.v) all_bus = system.Bus.idx.v load_bus = system.StaticLoad.get(src='bus', attr='v', idx=system.StaticLoad.get_idx()) idx_PD = system.PQ.find_idx(keys="bus", values=all_bus, allow_none=True, default=None) - self.pl._v = np.array(system.PQ.get(src='p0', attr='v', idx=idx_PD)) + self.pl._v = c_sparse(system.PQ.get(src='p0', attr='v', idx=idx_PD)) self.ql._v = np.array(system.PQ.get(src='q0', attr='v', idx=idx_PD)) + row, col = np.meshgrid(all_bus, slack_bus) + self.Cs._v = c_sparse((row == col).astype(int)) row, col = np.meshgrid(all_bus, gen_bus) - self.Cg._v = (row == col).astype(int) + self.Cg._v = c_sparse((row == col).astype(int)) row, col = np.meshgrid(all_bus, load_bus) - self.Cl._v = (row == col).astype(int) + self.Cl._v = c_sparse((row == col).astype(int)) return True diff --git a/ams/core/param.py b/ams/core/param.py index 2a88ffc3..9b900577 100644 --- a/ams/core/param.py +++ b/ams/core/param.py @@ -3,27 +3,35 @@ """ -import logging # NOQA +import logging -from typing import Optional # NOQA +from typing import Optional, Iterable -import numpy as np # NOQA -from scipy.sparse import issparse # NOQA +import numpy as np +from scipy.sparse import issparse -from andes.core.common import Config # NOQA from andes.core import BaseParam, DataParam, IdxParam, NumParam, ExtParam # NOQA from andes.models.group import GroupBase # NOQA from ams.core.var import Algeb # NOQA +from ams.opt.omodel import Param logger = logging.getLogger(__name__) -class RParam: +class RParam(Param): """ Class for parameters used in a routine. This class is developed to simplify the routine definition. + `RParm` is further used to define `Parameter` in the optimization model. + + `no_parse` is used to skip parsing the `RParam` in optimization + model. + It means that the `RParam` will not be added to the optimization model. + This is useful when the RParam contains non-numeric values, + or it is not necessary to be added to the optimization model. + Parameters ---------- name : str, optional @@ -46,6 +54,32 @@ class RParam: Indexer of the parameter. imodel : str, optional Name of the owner model or group of the indexer. + 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. + complex: bool, optional + True to set the parameter as 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. + sparsity: list, optional + Sparsity pattern of the parameter. Examples -------- @@ -78,8 +112,25 @@ def __init__(self, v: Optional[np.ndarray] = None, indexer: Optional[str] = None, imodel: Optional[str] = None, + expand_dims: Optional[int] = None, + no_parse: Optional[bool] = False, + nonneg: Optional[bool] = False, + nonpos: Optional[bool] = False, + complex: 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, + sparsity: Optional[list] = None, ): - + Param.__init__(self, nonneg=nonneg, nonpos=nonpos, + complex=complex, imag=imag, symmetric=symmetric, + diag=diag, hermitian=hermitian, boolean=boolean, + integer=integer, pos=pos, neg=neg, sparsity=sparsity) self.name = name self.tex_name = tex_name if (tex_name is not None) else name self.info = info @@ -89,6 +140,8 @@ def __init__(self, self.model = model # name of a group or model self.indexer = indexer # name of the indexer self.imodel = imodel # name of a group or model of the indexer + self.expand_dims = expand_dims + self.no_parse = no_parse self.owner = None # instance of the owner model or group self.rtn = None # instance of the owner routine self.is_ext = False # indicate if the value is set externally @@ -106,18 +159,19 @@ def v(self): - This property is a wrapper for the ``get`` method of the owner class. - The value will sort by the indexer if indexed, used for optmization modeling. """ + out = None if self.indexer is None: if self.is_ext: if issparse(self._v): - return self._v.toarray() + out = self._v.toarray() else: - return self._v + out = self._v elif self.is_group: - return self.owner.get(src=self.src, attr='v', - idx=self.owner.get_idx()) + out = self.owner.get(src=self.src, attr='v', + idx=self.owner.get_idx()) else: src_param = getattr(self.owner, self.src) - return getattr(src_param, 'v') + out = getattr(src_param, 'v') else: try: imodel = getattr(self.rtn.system, self.imodel) @@ -130,15 +184,29 @@ def v(self): except Exception as e: raise e model = getattr(self.rtn.system, self.model) - sorted_v = model.get(src=self.src, attr='v', idx=sorted_idx) - return sorted_v + out = model.get(src=self.src, attr='v', idx=sorted_idx) + if self.expand_dims is not None: + out = np.expand_dims(out, axis=self.expand_dims) + return out @property def shape(self): """ Return the shape of the parameter. """ - return self.v.shape + return np.shape(self.v) + + @property + def dtype(self): + """ + Return the data type of the parameter value. + """ + if isinstance(self.v, (str, bytes)): + return str + elif isinstance(self.v, Iterable): + return type(self.v[0]) + else: + return type(self.v) @property def n(self): @@ -158,39 +226,8 @@ def class_name(self): return self.__class__.__name__ def __repr__(self): - if self.is_ext: - span = '' - if isinstance(self.v, np.ndarray): - if 1 in self.v.shape: - if len(self.v) <= 20: - span = f', v={self.v}' - else: - if len(self.v) <= 20: - span = f', v={self.v}' - else: - span = f', v in length of {len(self.v)}' - return f'{self.__class__.__name__}: {self.name}{span}' - else: - span = '' - if 1 <= self.n <= 20: - span = f', v={self.v}' - if hasattr(self, 'vin') and (self.vin is not None): - span += f', v in length of {self.vin}' - - if isinstance(self.v, np.ndarray): - if self.v.shape[0] == 1 or self.v.ndim == 1: - if len(self.v) <= 20: - span = f', v={self.v}' - else: - span = f', v in shape({self.v.shape})' - elif isinstance(self.v, list): - if len(self.v) <= 20: - span = f', v={self.v}' - else: - span = f', v in length of {len(self.v)}' - else: - span = f', v={self.v}' - return f'{self.__class__.__name__}: {self.name} <{self.owner.__class__.__name__}>{span}' + postfix = '' if self.src is None else f'.{self.src}' + return f'{self.__class__.__name__}: {self.owner.__class__.__name__}' + postfix def get_idx(self): """ diff --git a/ams/core/service.py b/ams/core/service.py index a0547dc3..1b052710 100644 --- a/ams/core/service.py +++ b/ams/core/service.py @@ -9,11 +9,13 @@ from andes.core.service import BaseService, BackRef, RefFlatten # NOQA +from ams.opt.omodel import Param + logger = logging.getLogger(__name__) -class RBaseService(BaseService): +class RBaseService(BaseService, Param): """ Base class for services that are used in a routine. Revised from module `andes.core.service.BaseService`. @@ -32,6 +34,8 @@ class RBaseService(BaseService): Variable type. model : str, optional Model name. + no_parse: bool, optional + True to skip parsing the service. """ def __init__(self, @@ -40,9 +44,12 @@ def __init__(self, unit: str = None, info: str = None, vtype: Type = None, + no_parse: bool = False, ): - super().__init__(name=name, tex_name=tex_name, unit=unit, - info=info, vtype=vtype) + Param.__init__(self, name=name, unit=unit, info=info, + no_parse=no_parse) + BaseService.__init__(self, name=name, tex_name=tex_name, unit=unit, + info=info, vtype=vtype) self.export = False self.is_group = False self.rtn = None @@ -72,22 +79,8 @@ def class_name(self): return self.__class__.__name__ def __repr__(self): - val_str = '' - - v = self.v - - if v is None: - return f'{self.class_name}: {self.owner.class_name}.{self.name}' - elif isinstance(v, np.ndarray): - if v.shape[0] == 1: - if len(self.v) <= 20: - val_str = f', v={self.v}' - else: - val_str = f', v in shape of {self.v.shape}' - else: - val_str = f', v in shape of {self.v.shape}' - - return f'{self.class_name}: {self.rtn.class_name}.{self.name}{val_str}' + if self.name is None: + return f'{self.class_name}: {self.rtn.class_name}' else: return f'{self.class_name}: {self.rtn.class_name}.{self.name}' @@ -119,9 +112,11 @@ def __init__(self, unit: str = None, info: str = None, vtype: Type = None, + no_parse: bool = False, ): super().__init__(name=name, tex_name=tex_name, unit=unit, - info=info, vtype=vtype) + info=info, vtype=vtype, + no_parse=no_parse) self._v = value @property @@ -160,9 +155,11 @@ def __init__(self, tex_name: str = None, unit: str = None, info: str = None, - vtype: Type = None,): + vtype: Type = None, + no_parse: bool = False,): super().__init__(name=name, tex_name=tex_name, unit=unit, - info=info, vtype=vtype) + info=info, vtype=vtype, + no_parse=no_parse) self.u = u @@ -176,8 +173,6 @@ class LoadScale(ROperationService): nodal load. sd : Callable zonal load factor. - Cl: Callable - Connection matrix for Load and Bus. name : str, optional Instance name. tex_name : str, optional @@ -191,27 +186,27 @@ class LoadScale(ROperationService): def __init__(self, u: Callable, sd: Callable, - Cl: Callable, name: str = None, tex_name: str = None, unit: str = None, info: str = None, + no_parse: bool = False, ): tex_name = tex_name if tex_name is not None else u.tex_name super().__init__(name=name, tex_name=tex_name, unit=unit, - info=info, u=u,) + info=info, u=u, no_parse=no_parse) self.sd = sd - self.Cl = Cl @property def v(self): - Region = self.rtn.system.Region - yloc_zl = np.array(Region.idx2uid(self.u.v)) - PQ = self.rtn.system.PQ - p0 = PQ.get(src='p0', attr='v', idx=PQ.idx.v) - p0s = np.multiply(self.sd.v[:, yloc_zl].transpose(), - p0[:, np.newaxis]) - return np.matmul(np.linalg.pinv(self.Cl.v), p0s) + sys = self.rtn.system + u_idx = self.u.get_idx() + u_bus = self.u.owner.get(src='bus', attr='v', idx=u_idx) + u_zone = sys.Bus.get(src='zone', attr='v', idx=u_bus) + u_yloc = np.array(sys.Region.idx2uid(u_zone)) + p0s = np.multiply(self.sd.v[:, u_yloc].transpose(), + self.u.v[:, np.newaxis]) + return p0s class NumOp(ROperationService): @@ -245,6 +240,8 @@ class NumOp(ROperationService): Keyword arguments to pass to ``rfun``. expand_dims : int, optional Expand the dimensions of the output array along a specified axis. + array_out : bool, optional + Whether to force the output to be an array. """ def __init__(self, @@ -259,10 +256,12 @@ def __init__(self, rfun: Callable = None, rargs: dict = {}, expand_dims: int = None, - array_out=True): + array_out=True, + no_parse: bool = False,): tex_name = tex_name if tex_name is not None else u.tex_name super().__init__(name=name, tex_name=tex_name, unit=unit, - info=info, vtype=vtype, u=u,) + info=info, vtype=vtype, u=u, + no_parse=no_parse) self.fun = fun self.args = args self.rfun = rfun @@ -317,8 +316,8 @@ class NumExpandDim(NumOp): Description. vtype : Type, optional Variable type. - model : str, optional - Model name. + array_out : bool, optional + Whether to force the output to be an array. """ def __init__(self, @@ -330,11 +329,13 @@ def __init__(self, unit: str = None, info: str = None, vtype: Type = None, - array_out: bool = True,): + array_out: bool = True, + no_parse: bool = False,): super().__init__(name=name, tex_name=tex_name, unit=unit, info=info, vtype=vtype, u=u, fun=np.expand_dims, args=args, - array_out=array_out) + array_out=array_out, + no_parse=no_parse) self.axis = axis @property @@ -375,6 +376,8 @@ class NumOpDual(NumOp): Keyword arguments to pass to ``rfun``. expand_dims : int, optional Expand the dimensions of the output array along a specified axis. + array_out : bool, optional + Whether to force the output to be an array. """ def __init__(self, @@ -390,14 +393,16 @@ def __init__(self, rfun: Callable = None, rargs: dict = {}, expand_dims: int = None, - array_out=True): + array_out=True, + no_parse: bool = False,): tex_name = tex_name if tex_name is not None else u.tex_name super().__init__(name=name, tex_name=tex_name, unit=unit, info=info, vtype=vtype, u=u, fun=fun, args=args, rfun=rfun, rargs=rargs, expand_dims=expand_dims, - array_out=array_out) + array_out=array_out, + no_parse=no_parse) self.u2 = u2 @property @@ -437,15 +442,18 @@ def __init__(self, tex_name: str = None, unit: str = None, info: str = None, - vtype: Type = None,): + vtype: Type = None, + no_parse: bool = False,): tex_name = tex_name if tex_name is not None else u.tex_name super().__init__(name=name, tex_name=tex_name, unit=unit, info=info, vtype=vtype, u=u, u2=u2, fun=None, args=None, rfun=None, rargs=None, - expand_dims=None) + expand_dims=None, + no_parse=no_parse) if self.u.horizon is None: - msg = f'{self.class_name} {self.name}.u {self.u.name} has no horizon, likely a modeling error.' + msg = f'{self.class_name} <{self.name}>.u: <{self.u.name}> ' + msg += 'has no horizon, likely a modeling error.' logger.error(msg) @property @@ -468,7 +476,8 @@ def v(self): class NumHstack(NumOp): """ - Repeat an array along the second axis nc times + Repeat an array along the second axis nc times or the length of + reference array, using NumPy's hstack function, where nc is the column number of the reference array, ``np.hstack([u.v[:, np.newaxis] * ref.shape[1]], **kwargs)``. @@ -503,20 +512,22 @@ def __init__(self, info: str = None, vtype: Type = None, rfun: Callable = None, - rargs: dict = {}): + rargs: dict = {}, + no_parse: bool = False,): super().__init__(name=name, tex_name=tex_name, unit=unit, info=info, vtype=vtype, u=u, fun=np.hstack, args=args, - rfun=rfun, rargs=rargs) + rfun=rfun, rargs=rargs, + no_parse=no_parse) self.ref = ref @property def v0(self): nc = 1 - if hasattr(self.ref, "shape"): - nc = self.ref.shape[1] - elif isinstance(self.ref.v, (list, tuple)): + if isinstance(self.ref.v, (list, tuple)): nc = len(self.ref.v) + elif hasattr(self.ref, "shape"): + nc = self.ref.shape[1] else: raise AttributeError(f"{self.rtn.class_name}: ref {self.ref.name} has no attribute shape nor length.") return self.fun([self.u.v[:, np.newaxis]] * nc, @@ -585,11 +596,13 @@ def __init__(self, vtype: Type = None, rfun: Callable = None, rargs: dict = {}, + no_parse: bool = False, ): super().__init__(name=name, tex_name=tex_name, unit=unit, info=info, vtype=vtype, u=u, fun=None, args={}, - rfun=rfun, rargs=rargs) + rfun=rfun, rargs=rargs, + no_parse=no_parse) self.zone = zone @property @@ -619,17 +632,45 @@ 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 + power from StaticGen `pg`, following definition can be used: + ```python + class RTED: + ... + self.ce = VarSelect(u=self.pg, indexer='genE') + ... + ``` + Parameters ---------- u : Callable The input matrix variable. - idx : list - The index of the subset. + indexer: str + The name of the indexer source. + gamma : str, optional + The name of the indexer gamma. + name : str, optional + The name of the instance. + tex_name : str, optional + The TeX name for the instance. + unit : str, optional + The unit of the output. + info : str, optional + A description of the operation. + vtype : Type, optional + The variable type. + rfun : Callable, optional + Function to apply to the output of ``fun``. + rargs : dict, optional + Keyword arguments to pass to ``rfun``. + array_out : bool, optional + Whether to force the output to be an array. """ def __init__(self, u: Callable, indexer: str, + gamma: str = None, name: str = None, tex_name: str = None, unit: str = None, @@ -637,12 +678,17 @@ def __init__(self, vtype: Type = None, rfun: Callable = None, rargs: dict = {}, + array_out: bool = True, + no_parse: bool = False, **kwargs ): super().__init__(name=name, tex_name=tex_name, unit=unit, info=info, vtype=vtype, u=u, fun=None, - rfun=rfun, rargs=rargs, **kwargs) + rfun=rfun, rargs=rargs, array_out=array_out, + no_parse=no_parse, + **kwargs) self.indexer = indexer + self.gamma = gamma @property def v0(self): @@ -680,12 +726,12 @@ def v0(self): if not is_subset: raise ValueError(f'{indexer.model} contains undefined {indexer.src}, likey a data error.') - out = [1 if item in ref else 0 for item in uidx] - out = np.array(out) - if self.u.horizon is not None: - out = out[:, np.newaxis] - out = np.repeat(out, self.u.horizon.n, axis=1) - return np.array(out) + row, col = np.meshgrid(uidx, ref) + out = (row == col).astype(int) + if self.gamma: + vgamma = getattr(self.rtn, self.gamma) + out = vgamma.v[:, np.newaxis] * out + return out class VarReduction(NumOp): @@ -723,11 +769,13 @@ def __init__(self, vtype: Type = None, rfun: Callable = None, rargs: dict = {}, + no_parse: bool = False, **kwargs ): 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, **kwargs) self.fun = fun @@ -776,10 +824,12 @@ def __init__(self, vtype: Type = None, rfun: Callable = None, rargs: dict = {}, + no_parse: bool = False, ): super().__init__(name=name, tex_name=tex_name, unit=unit, info=info, vtype=vtype, - u=u, fun=None, rfun=rfun, rargs=rargs,) + u=u, fun=None, rfun=rfun, rargs=rargs, + no_parse=no_parse,) @property def v0(self): diff --git a/ams/core/symprocessor.py b/ams/core/symprocessor.py index eb5ea42e..c3e7d27b 100644 --- a/ams/core/symprocessor.py +++ b/ams/core/symprocessor.py @@ -4,10 +4,10 @@ This module is revised from ``andes.core.symprocessor``. """ -import logging # NOQA -from collections import OrderedDict # NOQA +import logging +from collections import OrderedDict -import sympy as sp # NOQA +import sympy as sp logger = logging.getLogger(__name__) @@ -65,13 +65,21 @@ def __init__(self, parent): (r' dot ', r' * '), (r'\bsum\b', f'{lang}.sum'), (r'\bvar\b', f'{lang}.Variable'), + (r'\bparam\b', f'{lang}.Parameter'), + (r'\bconst\b', f'{lang}.Constant'), (r'\bproblem\b', f'{lang}.Problem'), (r'\bmultiply\b', f'{lang}.multiply'), + (r'\bmul\b', f'{lang}.multiply'), # alias for multiply (r'\bvstack\b', f'{lang}.vstack'), (r'\bnorm\b', f'{lang}.norm'), (r'\bpos\b', f'{lang}.pos'), (r'\bpower\b', f'{lang}.power'), (r'\bsign\b', f'{lang}.sign'), + (r'\bsquare\b', f'{lang}.square'), + (r'\bquad_over_lin\b', f'{lang}.quad_over_lin'), + (r'\bdiag\b', f'{lang}.diag'), + (r'\bquad_form\b', f'{lang}.quad_form'), + (r'\bsum_squares\b', f'{lang}.sum_squares'), ]) self.tex_map = OrderedDict([ @@ -79,9 +87,13 @@ def __init__(self, parent): (r'\b(\w+)\s*\*\s*(\w+)\b', r'\1 \2'), (r'\@', r' '), (r'dot', r' '), + (r'sum_squares\((.*?)\)', r"SUM((\1))^2"), (r'multiply\(([^,]+), ([^)]+)\)', r'\1 \2'), (r'\bnp.linalg.pinv(\d+)', r'\1^{\-1}'), (r'\bpos\b', 'F^{+}'), + (r'mul\((.*?),\s*(.*?)\)', r'\1 \2'), + (r'\bmul\b\((.*?),\s*(.*?)\)', r'\1 \2'), + (r'\bsum\b', 'SUM'), ]) self.status = { @@ -123,7 +135,8 @@ def generate_symbols(self, force_generate=False): for rpname, rparam in self.parent.rparams.items(): tmp = sp.symbols(f'{rparam.name}') self.inputs_dict[rpname] = tmp - self.sub_map[rf"\b{rpname}\b"] = f'self.om.rtn.{rpname}.v' + sub_name = f'self.rtn.{rpname}.v' if rparam.no_parse else f'self.om.{rpname}' + self.sub_map[rf"\b{rpname}\b"] = sub_name self.tex_map[rf"\b{rpname}\b"] = f'{rparam.tex_name}' # Routine Services @@ -131,7 +144,8 @@ def generate_symbols(self, force_generate=False): tmp = sp.symbols(f'{service.name}') self.services_dict[sname] = tmp self.inputs_dict[sname] = tmp - self.sub_map[rf"\b{sname}\b"] = f'self.om.rtn.{sname}.v' + sub_name = f'self.rtn.{sname}.v' if service.no_parse else f'self.om.{sname}' + self.sub_map[rf"\b{sname}\b"] = sub_name self.tex_map[rf"\b{sname}\b"] = f'{service.tex_name}' # store tex names defined in `self.config` diff --git a/ams/main.py b/ams/main.py index e5a4c21b..fcffe07d 100644 --- a/ams/main.py +++ b/ams/main.py @@ -515,23 +515,17 @@ def versioninfo(): import numpy as np import cvxpy import andes + from ams.shared import INSTALLED_SOLVERS versions = {'Python': platform.python_version(), 'ams': get_versions()['version'], 'andes': andes.__version__, 'numpy': np.__version__, 'cvxpy': cvxpy.__version__, + 'solvers': ', '.join(INSTALLED_SOLVERS), } maxwidth = max([len(k) for k in versions.keys()]) - try: - import numba - except ImportError: - numba = None - - if numba is not None: - versions["numba"] = numba.__version__ - for key, val in versions.items(): print(f"{key: <{maxwidth}} {val}") diff --git a/ams/models/distributed/esd1.py b/ams/models/distributed/esd1.py index 293f3e78..baa69d9c 100644 --- a/ams/models/distributed/esd1.py +++ b/ams/models/distributed/esd1.py @@ -32,12 +32,12 @@ def __init__(self): ) self.gammap = NumParam(default=1.0, tex_name=r'\gamma_p', - info='Ratio of PVD1.pref0 w.r.t to that of static PV', + info='Ratio of ESD1.pref0 w.r.t to that of static PV', vrange='(0, 1]', ) self.gammaq = NumParam(default=1.0, tex_name=r'\gamma_q', - info='Ratio of PVD1.qref0 w.r.t to that of static PV', + info='Ratio of ESD1.qref0 w.r.t to that of static PV', vrange='(0, 1]', ) diff --git a/ams/models/timeslot.py b/ams/models/timeslot.py index c57b5290..a304e4af 100644 --- a/ams/models/timeslot.py +++ b/ams/models/timeslot.py @@ -30,22 +30,39 @@ def __init__(self, system=None, config=None): tex_name=r's_{d}', iconvert=str_list_iconv, oconvert=str_list_oconv, - vtype=float, - ) + vtype=float) class EDTSlot(TimeSlot): """ Time slot model for ED. + + `sd` is the zonal load scaling factor. + Cells in `sd` should have `nz` values seperated by comma, + where `nz` is the number of `Region` in the system. + + `ug` is the unit commitment decisions. + Cells in `ug` should have `ng` values seperated by comma, + where `ng` is the number of `StaticGen` in the system. """ def __init__(self, system=None, config=None): TimeSlot.__init__(self, system, config) + self.ug = NumParam(info='unit commitment decisions', + tex_name=r'u_{g}', + iconvert=str_list_iconv, + oconvert=str_list_oconv, + vtype=int) + class UCTSlot(TimeSlot): """ Time slot model for UC. + + `sd` is the zonal load scaling factor. + Cells in `sd` should have `nz` values seperated by comma, + where `nz` is the number of `Region` in the system. """ def __init__(self, system=None, config=None): diff --git a/ams/opt/omodel.py b/ams/opt/omodel.py index 335a3b6b..c8af43d5 100644 --- a/ams/opt/omodel.py +++ b/ams/opt/omodel.py @@ -1,23 +1,19 @@ """ -Module for optimization models. +Module for optimization modeling. """ +import logging -import logging # NOQA +from typing import Any, Optional, Union +from collections import OrderedDict +import re -from typing import Optional, Union # NOQA -from collections import OrderedDict # NOQA -import re # NOQA +import numpy as np -import numpy as np # NOQA +from andes.core.common import Config -from andes.core.common import Config # NOQA +from ams.utils import timer -from ams.core.param import RParam # NOQA -from ams.core.var import Algeb # NOQA - -from ams.utils import timer # NOQA - -import cvxpy as cp # NOQA +import cvxpy as cp logger = logging.getLogger(__name__) @@ -49,6 +45,8 @@ def __init__(self, self.info = info self.unit = unit self.is_disabled = False + self.rtn = None + self.optz = None # corresponding optimization element def parse(self): """ @@ -73,22 +71,15 @@ def n(self): else: return self.owner.n - @property - def rtn(self): - """ - Return the owner routine. - """ - return self.om.rtn - @property def shape(self): """ Return the shape. """ - if self.rtn.initialized: + try: return self.om.__dict__[self.name].shape - else: - logger.warning(f'<{self.rtn.class_name}> is not initialized yet.') + except KeyError: + logger.warning('Shape info is not ready before initialziation.') return None @property @@ -103,7 +94,114 @@ def size(self): return None -class Var(Algeb, OptzBase): +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. + complex: bool, optional + True to set the parameter as 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. + sparsity: list, optional + Sparsity pattern of the parameter. + """ + + 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, + complex: 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, + sparsity: Optional[list] = None, + ): + OptzBase.__init__(self, name=name, info=info, unit=unit) + self.no_parse = no_parse # True to skip parsing the parameter + self.sparsity = sparsity # sparsity pattern + + self.config = Config(name=self.class_name) # `config` that can be exported + + self.config.add(OrderedDict((('nonneg', nonneg), + ('nonpos', nonpos), + ('complex', complex), + ('imag', imag), + ('symmetric', symmetric), + ('diag', diag), + ('hermitian', hermitian), + ('boolean', boolean), + ('integer', integer), + ('pos', pos), + ('neg', neg), + ))) + + def parse(self): + """ + Parse the parameter. + """ + config = self.config.as_dict() # NOQA + sub_map = self.om.rtn.syms.sub_map + shape = np.shape(self.v) + code_param = f"tmp=param(shape={shape}, **config)" + for pattern, replacement, in sub_map.items(): + code_param = re.sub(pattern, replacement, code_param) + exec(code_param) + exec("self.optz=tmp") + try: + exec("self.optz.value = self.v") + except ValueError: + msg = f"Parameter <{self.name}> has non-numeric value, " + msg += "no_parse=True is applied." + logger.warning(msg) + self.no_parse = True + return False + 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}' + + +class Var(OptzBase): """ Base class for variables used in a routine. @@ -128,12 +226,6 @@ class Var(Algeb, OptzBase): Source variable name. If is None, the value of `name` will be used. model : str, optional Name of the owner model or group. - lb : str, optional - Lower bound - ub : str, optional - Upper bound - ctrl : str, optional - Controllability horizon : ams.routines.RParam, optional Horizon idx. nonneg : bool, optional @@ -181,11 +273,8 @@ def __init__(self, unit: Optional[str] = None, model: Optional[str] = None, shape: Optional[Union[tuple, int]] = None, - lb: Optional[str] = None, - ub: Optional[str] = None, - ctrl: Optional[str] = None, v0: Optional[str] = None, - horizon: Optional[RParam] = None, + horizon=None, nonneg: Optional[bool] = False, nonpos: Optional[bool] = False, complex: Optional[bool] = False, @@ -200,8 +289,6 @@ def __init__(self, pos: Optional[bool] = False, neg: Optional[bool] = False, ): - # Algeb.__init__(self, name=name, tex_name=tex_name, info=info, unit=unit) - # below info is the same as Algeb self.name = name self.info = info self.unit = unit @@ -214,14 +301,14 @@ def __init__(self, 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.lb = lb - self.ub = ub - self.ctrl = ctrl self.v0 = v0 self.horizon = horizon self._shape = shape self._v = None + self.optz_lb = None + self.optz_ub = None + self.config = Config(name=self.class_name) # `config` that can be exported self.config.add(OrderedDict((('nonneg', nonneg), @@ -248,7 +335,10 @@ def v(self): """ Return the CVXPY variable value. """ - out = self.om.vars[self.name].value if self._v is None else self._v + if self.name in self.om.vars: + out = self.om.vars[self.name].value if self._v is None else self._v + else: + out = None return out @v.setter @@ -268,7 +358,6 @@ def parse(self): """ Parse the variable. """ - om = self.om # NOQA sub_map = self.om.rtn.syms.sub_map # only used for CVXPY # NOTE: Config only allow lower case letters, do a conversion here @@ -305,45 +394,12 @@ def parse(self): code_var = f"tmp=var({shape}, **config)" for pattern, replacement, in sub_map.items(): code_var = re.sub(pattern, replacement, code_var) - exec(code_var) - exec("om.vars[self.name] = tmp") - exec(f'setattr(om, self.name, om.vars["{self.name}"])') - u_ctrl = self.ctrl.v if self.ctrl else np.ones(nr) - v0 = self.v0.v if self.v0 else np.zeros(nr) - if self.lb: - lv = self.lb.owner.get(src=self.lb.name, idx=self.get_idx(), attr='v') - u = self.lb.owner.get(src='u', idx=self.get_idx(), attr='v') - # element-wise lower bound considering online status - elv = u_ctrl * u * lv + (1 - u_ctrl) * v0 - # fit variable shape if horizon exists - elv = np.tile(elv, (nc, 1)).T if nc > 0 else elv - exec("om.constrs[self.lb.name] = tmp >= elv") - if self.ub: - uv = self.ub.owner.get(src=self.ub.name, idx=self.get_idx(), attr='v') - u = self.lb.owner.get(src='u', idx=self.get_idx(), attr='v') - # element-wise upper bound considering online status - euv = u_ctrl * u * uv + (1 - u_ctrl) * v0 - # fit variable shape if horizon exists - euv = np.tile(euv, (nc, 1)).T if nc > 0 else euv - exec("om.constrs[self.ub.name] = tmp <= euv") + exec(code_var) # build the Var object + exec("self.optz = tmp") # assign the to the optz attribute return True def __repr__(self): - if self.owner.n == 0: - span = [] - - elif 1 <= self.owner.n <= 20: - span = f'a={self.a}, v={self.v}' - - else: - span = [] - span.append(self.a[0]) - span.append(self.a[-1]) - span.append(self.a[1] - self.a[0]) - span = ':'.join([str(i) for i in span]) - span = 'a=[' + span + ']' - - return f'{self.__class__.__name__}: {self.owner.__class__.__name__}.{self.name}, {span}' + return f'{self.__class__.__name__}: {self.owner.__class__.__name__}.{self.name}' class Constraint(OptzBase): @@ -411,7 +467,6 @@ def parse(self, no_code=True): sub_map = self.om.rtn.syms.sub_map if self.is_disabled: return True - om = self.om # NOQA code_constr = self.e_str for pattern, replacement in sub_map.items(): try: @@ -420,22 +475,33 @@ def parse(self, no_code=True): logger.error(f"Error in parsing constr <{self.name}>.") raise e if self.type == 'uq': - code_constr = f'{code_constr} <= 0' + code_constr = f'tmp = {code_constr} <= 0' elif self.type == 'eq': - code_constr = f'{code_constr} == 0' + code_constr = f'tmp = {code_constr} == 0' else: raise ValueError(f'Constraint type {self.type} is not supported.') - code_constr = f'om.constrs["{self.name}"]=' + code_constr logger.debug(f"Set constrs {self.name}: {self.e_str} {'<= 0' if self.type == 'uq' else '== 0'}") if not no_code: logger.info(f"Code Constr: {code_constr}") exec(code_constr) - exec(f'setattr(om, self.name, om.constrs["{self.name}"])') + exec("self.optz = tmp") return True def __repr__(self): enabled = 'ON' if self.name in self.om.constrs else 'OFF' - return f"[{enabled}]: {self.e_str}" + out = f"[{enabled}]: {self.e_str}" + out += " =0" if self.type == 'eq' else " <=0" + return out + + @property + def v(self): + """ + Return the CVXPY constraint LHS value. + """ + if self.name in self.om.constrs: + return self.om.constrs[self.name]._expr.value + else: + return None class Objective(OptzBase): @@ -508,7 +574,6 @@ def parse(self, no_code=True): no_code : bool, optional Flag indicating if the code should be shown, True by default. """ - om = self.om # NOQA sub_map = self.om.rtn.syms.sub_map code_obj = self.e_str for pattern, replacement, in sub_map.items(): @@ -519,7 +584,7 @@ def parse(self, no_code=True): code_obj = f'cp.Maximize({code_obj})' else: raise ValueError(f'Objective sense {self.sense} is not supported.') - code_obj = 'om.obj=' + code_obj + code_obj = 'self.om.obj=' + code_obj logger.debug(f"Set obj {self.name}: {self.sense}. {self.e_str}") if not no_code: logger.info(f"Code Obj: {code_obj}") @@ -528,12 +593,16 @@ def parse(self, no_code=True): def __repr__(self): name_str = f"{self.name}=" if self.name is not None else "obj=" - value_str = f"{self.v:.4f}, " if self.v is not None else "" + try: + v = self.v + except Exception: + v = None + value_str = f"{v:.4f}, " if v is not None else "" return f"{name_str}{value_str}{self.e_str}" class OModel: - r""" + """ Base class for optimization models. Parameters @@ -545,6 +614,8 @@ class OModel: ---------- mdl: cvxpy.Problem Optimization model. + params: OrderedDict + Parameters. vars: OrderedDict Decision variables. constrs: OrderedDict @@ -555,15 +626,12 @@ class OModel: Number of decision variables. m: int Number of constraints. - - TODO: - - Add _check_attribute and _register_attribute for vars, constrs, and obj. - - Add support for user-defined vars, constrs, and obj. """ def __init__(self, routine): self.rtn = routine self.mdl = None + self.params = OrderedDict() self.vars = OrderedDict() self.constrs = OrderedDict() self.obj = None @@ -591,18 +659,46 @@ def setup(self, no_code=True, force_generate=False): """ rtn = self.rtn rtn.syms.generate_symbols(force_generate=force_generate) + # --- add RParams and Services as parameters --- + for key, val in rtn.params.items(): + if not val.no_parse: + try: + val.parse() + except Exception as e: + msg = f"Failed to parse Param <{key}>. " + msg += f"Original error: {e}" + raise Exception(msg) + setattr(self, key, val.optz) # --- add decision variables --- - for ovar in rtn.vars.values(): - ovar.parse() + for key, val in rtn.vars.items(): + try: + val.parse() + except Exception as e: + msg = f"Failed to parse Var <{key}>. " + msg += f"Original error: {e}" + raise Exception(msg) + setattr(self, key, val.optz) # --- add constraints --- - for constr in rtn.constrs.values(): - constr.parse(no_code=no_code) + for key, val in rtn.constrs.items(): + try: + val.parse(no_code=no_code) + except Exception as e: + msg = f"Failed to parse Constr <{key}>. " + msg += f"Original error: {e}" + raise Exception(msg) + setattr(self, key, val.optz) # --- parse objective functions --- if rtn.type == 'PF': # NOTE: power flow type has no objective function pass elif rtn.obj is not None: - rtn.obj.parse(no_code=no_code) + try: + rtn.obj.parse(no_code=no_code) + except Exception as e: + msg = f"Failed to parse Obj <{rtn.obj.name}>. " + msg += f"Original error: {e}" + raise Exception(msg) + # --- finalize the optimziation formulation --- code_mdl = "problem(self.obj, [constr for constr in self.constrs.values()])" for pattern, replacement in self.rtn.syms.sub_map.items(): @@ -628,3 +724,37 @@ def class_name(self): Return the class name """ return self.__class__.__name__ + + def __register_attribute(self, key, value): + """ + Register a pair of attributes to OModel instance. + + Called within ``__setattr__``, this is where the magic happens. + Subclass attributes are automatically registered based on the variable type. + """ + 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 + + def __setattr__(self, __name: str, __value: Any): + self.__register_attribute(__name, __value) + super().__setattr__(__name, __value) + + def update(self, params): + """ + Update the Parameter values. + + Parameters + ---------- + params: list + List of parameters to be updated. + """ + for param in params: + param.update() + return True + + def __repr__(self) -> str: + return f'{self.rtn.class_name}.{self.__class__.__name__} at {hex(id(self))}' diff --git a/ams/routines/__init__.py b/ams/routines/__init__.py index 29b07a21..8e2cb649 100644 --- a/ams/routines/__init__.py +++ b/ams/routines/__init__.py @@ -12,10 +12,10 @@ ('cpf', ['CPF']), ('acopf', ['ACOPF']), ('dcopf', ['DCOPF']), - ('ed', ['ED', 'ED2']), - ('rted', ['RTED', 'RTED2']), - ('uc', ['UC', 'UC2']), - ('dopf', ['LDOPF', 'LDOPF2']), + ('ed', ['ED', 'EDES']), + ('rted', ['RTED', 'RTEDES']), + ('uc', ['UC', 'UCES']), + ('dopf', ['DOPF', 'DOPFVIS']), ]) class_names = list_flatten(list(all_routines.values())) diff --git a/ams/routines/acopf.py b/ams/routines/acopf.py index 082cecd0..bfa2a97d 100644 --- a/ams/routines/acopf.py +++ b/ams/routines/acopf.py @@ -10,26 +10,12 @@ from ams.io.pypower import system2ppc # NOQA from ams.core.param import RParam # NOQA -from ams.routines.dcopf import DCOPFData # NOQA from ams.routines.dcpf import DCPFlowBase # NOQA from ams.opt.omodel import Var, Constraint, Objective # NOQA logger = logging.getLogger(__name__) -class ACOPFData(DCOPFData): - """ - ACOPF data. - """ - - def __init__(self): - DCOPFData.__init__(self) - self.ql = RParam(info='reactive power demand (system base)', - name='ql', tex_name=r'q_{l}', - model='mats', src='ql', - unit='p.u.',) - - class ACOPFBase(DCPFlowBase): """ Base class for ACOPF model. @@ -37,8 +23,6 @@ class ACOPFBase(DCPFlowBase): def __init__(self, system, config): DCPFlowBase.__init__(self, system, config) - self.info = 'AC Optimal Power Flow' - self.type = 'ACED' # NOTE: ACOPF does not receive data from dynamic self.map1 = OrderedDict() self.map2 = OrderedDict([ @@ -102,15 +86,42 @@ def run(self, force_init=False, no_code=True, **kwargs, ) -class ACOPFModel(ACOPFBase): +class ACOPF(ACOPFBase): """ - ACOPF model. + Standard AC optimal power flow. + + Notes + ----- + 1. ACOPF is solved with PYPOWER ``runopf`` function. + 2. ACOPF 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): ACOPFBase.__init__(self, system, config) self.info = 'AC Optimal Power Flow' self.type = 'ACED' + + # --- params --- + self.c2 = RParam(info='Gen cost coefficient 2', + name='c2', tex_name=r'c_{2}', + unit=r'$/(p.u.^2)', model='GCost', + indexer='gen', imodel='StaticGen', + nonneg=True) + self.c1 = RParam(info='Gen cost coefficient 1', + name='c1', tex_name=r'c_{1}', + unit=r'$/(p.u.)', model='GCost', + indexer='gen', imodel='StaticGen',) + self.c0 = RParam(info='Gen cost coefficient 0', + name='c0', tex_name=r'c_{0}', + unit=r'$', model='GCost', + indexer='gen', imodel='StaticGen', + no_parse=True) + self.ql = RParam(info='reactive power demand (system base)', + name='ql', tex_name=r'q_{l}', + model='mats', src='ql', + unit='p.u.',) # --- bus --- self.aBus = Var(info='Bus voltage angle', unit='rad', @@ -141,20 +152,3 @@ def __init__(self, system, config): info='total cost', e_str='sum(c2 * pg**2 + c1 * pg + c0)', sense='min',) - - -class ACOPF(ACOPFData, ACOPFModel): - """ - Standard AC optimal power flow. - - Notes - ----- - 1. ACOPF is solved with PYPOWER ``runopf`` function. - 2. ACOPF 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=None, config=None): - ACOPFData.__init__(self) - ACOPFModel.__init__(self, system, config) diff --git a/ams/routines/cpf.py b/ams/routines/cpf.py index 851ab069..cfff40fe 100644 --- a/ams/routines/cpf.py +++ b/ams/routines/cpf.py @@ -1,33 +1,29 @@ """ Continuous power flow routine. """ -import logging # NOQA +import logging -from ams.pypower import runcpf # NOQA +from ams.pypower import runcpf -from ams.io.pypower import system2ppc # NOQA -from ams.pypower.core import ppoption # NOQA +from ams.io.pypower import system2ppc +from ams.pypower.core import ppoption -from ams.routines.pflow import PFlowData, PFlowModel # NOQA +from ams.routines.pflow import PFlow logger = logging.getLogger(__name__) -class CPFModel(PFlowModel): +class CPF(PFlow): """ - Model for continuous power flow. + Continuous power flow. + + Still under development, not ready for use. """ def __init__(self, system, config): - PFlowModel.__init__(self, system, config) + PFlow.__init__(self, system, config) self.info = 'AC continuous power flow' self.type = 'PF' - # TODO: delete vars, constraints, and objectives - # FIXME: how? - # for v, _ in self.vars.items(): - # delattr(self, v) - # for c, _ in self.constraints.items(): - # delattr(self, c) def solve(self, method=None, **kwargs): """ @@ -67,15 +63,3 @@ def run(self, force_init=False, no_code=True, super().run(force_init=force_init, no_code=no_code, method=method, **kwargs, ) - - -class CPF(PFlowData, CPFModel): - """ - Continuous power flow. - - Still under development, not ready for use. - """ - - def __init__(self, system=None, config=None): - PFlowData.__init__(self) - CPFModel.__init__(self, system, config) diff --git a/ams/routines/dcopf.py b/ams/routines/dcopf.py index 135520af..e311630d 100644 --- a/ams/routines/dcopf.py +++ b/ams/routines/dcopf.py @@ -1,36 +1,55 @@ """ -OPF routines. +DCOPF routines. """ -import logging # NOQA +import logging -from ams.core.param import RParam # NOQA +import numpy as np +from ams.core.param import RParam +from ams.core.service import NumOp, NumOpDual -from ams.routines.routine import RoutineData, RoutineModel # NOQA +from ams.routines.routine import RoutineModel -from ams.opt.omodel import Var, Constraint, Objective # NOQA +from ams.opt.omodel import Var, Constraint, Objective logger = logging.getLogger(__name__) -class DCOPFData(RoutineData): +class DCOPFBase(RoutineModel): """ - DCOPF data. + Base class for DCOPF dispatch model. + + Overload the ``solve``, ``unpack``, and ``run`` methods. """ - def __init__(self): - RoutineData.__init__(self) + def __init__(self, system, config): + RoutineModel.__init__(self, system, config) # --- generator cost --- self.ug = RParam(info='Gen connection status', name='ug', tex_name=r'u_{g}', - model='StaticGen', src='u',) + model='StaticGen', src='u', + no_parse=True) self.ctrl = RParam(info='Gen controllability', name='ctrl', tex_name=r'c_{trl}', - model='StaticGen', src='ctrl',) + 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.c2 = RParam(info='Gen cost coefficient 2', name='c2', tex_name=r'c_{2}', unit=r'$/(p.u.^2)', model='GCost', - indexer='gen', imodel='StaticGen',) + indexer='gen', imodel='StaticGen', + nonneg=True) self.c1 = RParam(info='Gen cost coefficient 1', name='c1', tex_name=r'c_{1}', unit=r'$/(p.u.)', model='GCost', @@ -38,43 +57,67 @@ def __init__(self): self.c0 = RParam(info='Gen cost coefficient 0', name='c0', tex_name=r'c_{0}', unit=r'$', model='GCost', - indexer='gen', imodel='StaticGen',) - # --- generator limit --- - self.pmax = RParam(info='Gen maximum active power (system base)', - name='pmax', tex_name=r'p_{max}', - unit='p.u.', model='StaticGen',) - self.pmin = RParam(info='Gen minimum active power (system base)', - name='pmin', tex_name=r'p_{min}', - unit='p.u.', model='StaticGen',) - self.pg0 = RParam(info='Gen initial active power (system base)', + indexer='gen', imodel='StaticGen', + no_parse=True) + # --- generator --- + 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='p0', tex_name=r'p_{g,0}', unit='p.u.', model='StaticGen',) - self.Cg = RParam(info='connection matrix for Gen and Bus', - name='Cg', tex_name=r'C_{g}', - model='mats', src='Cg',) + # --- load --- - self.pl = RParam(info='nodal active load (system base)', - name='pl', tex_name=r'p_{l}', - model='mats', src='pl', + self.pd = RParam(info='active demand', + name='pd', tex_name=r'p_{d}', + model='StaticLoad', src='p0', unit='p.u.',) + # --- line --- + self.x = RParam(info='line reactance', + name='x', tex_name=r'x', + model='Line', src='x', + unit='p.u.', no_parse=True,) self.rate_a = RParam(info='long-term flow limit', name='rate_a', tex_name=r'R_{ATEA}', - unit='MVA', model='Line',) - self.PTDF = RParam(info='Power transfer distribution factor matrix', - name='PTDF', tex_name=r'P_{TDF}', - model='mats', src='PTDF',) + unit='MVA', model='Line') + # --- connection matrix --- + self.Cg = RParam(info='Gen connection matrix', + name='Cg', tex_name=r'C_{g}', + model='mats', src='Cg', + no_parse=True,) + self.Cs = RParam(info='Slack connection matrix', + name='Cs', tex_name=r'C_{s}', + model='mats', src='Cs', + no_parse=True,) + self.Cl = RParam(info='Load connection matrix', + name='Cl', tex_name=r'C_{l}', + model='mats', src='Cl', + no_parse=True,) + self.Cft = RParam(info='Line connection matrix', + name='Cft', tex_name=r'C_{ft}', + model='mats', src='Cft', + no_parse=True,) + self.PTDF = RParam(info='Power Transfer Distribution Factor', + name='PTDF', tex_name=r'P_{TDF}', + model='mats', src='PTDF', + no_parse=True,) -class DCOPFBase(RoutineModel): - """ - Base class for DCOPF dispatch model. + self.Cgi = NumOp(u=self.Cg, fun=np.linalg.pinv, + name='Cgi', tex_name=r'C_{g}^{-1}', + info='inverse of Cg', + no_parse=True,) + self.Cli = NumOp(u=self.Cl, fun=np.linalg.pinv, + name='Cli', tex_name=r'C_{l}^{-1}', + info='inverse of Cl', + no_parse=True,) - Overload the ``solve``, ``unpack``, and ``run`` methods. - """ - - def __init__(self, system, config): - RoutineModel.__init__(self, system, config) def solve(self, **kwargs): """ @@ -153,9 +196,12 @@ def unpack(self, **kwargs): return True -class DCOPFModel(DCOPFBase): +class DCOPF(DCOPFBase): """ - DCOPF model. + Standard DC optimal power flow (DCOPF). + + In this model, the bus injected power ``pn`` is used as internal variable + between generator output and load demand. """ def __init__(self, system, config): @@ -163,44 +209,48 @@ def __init__(self, system, config): self.info = 'DC Optimal Power Flow' self.type = 'DCED' # --- vars --- - self.pg = Var(info='Gen active power (system base)', - unit='p.u.', name='pg', src='p', - tex_name=r'p_{g}', - model='StaticGen', - lb=self.pmin, ub=self.pmax, - ctrl=self.ctrl, v0=self.pg0) - self.pn = Var(info='Bus active power injection (system base)', - unit='p.u.', name='pn', tex_name=r'p_{n}', - model='Bus',) + self.pg = Var(info='Gen active power', + unit='p.u.', + name='pg', tex_name=r'p_{g}', + model='StaticGen', src='p', + v0=self.pg0) + # NOTE: `ug*pmin` results in unexpected error + pglb = '-pg + mul(nctrle, pg0) + mul(ctrle, pmin)' + self.pglb = Constraint(name='pglb', info='pg min', + e_str=pglb, type='uq',) + pgub = 'pg - mul(nctrle, pg0) - mul(ctrle, pmax)' + self.pgub = Constraint(name='pgub', info='pg max', + e_str=pgub, type='uq',) + + self.aBus = Var(info='Bus voltage angle', + name='aBus', tex_name=r'\theta_{n}', + unit='rad', model='Bus',) + + self.plf = Var(info='Line active power', + name='plf', tex_name=r'p_{lf}', + unit='p.u.', model='Line',) + self.plflb = Constraint(name='plflb', info='Line power lower bound', + e_str='-plf - rate_a', type='uq',) + self.plfub = Constraint(name='plfub', info='Line power upper bound', + e_str='plf - rate_a', type='uq',) + # --- constraints --- self.pb = Constraint(name='pb', info='power balance', - e_str='sum(pl) - sum(pg)', + e_str='sum(pd) - sum(pg)', type='eq',) - self.pinj = Constraint(name='pinj', - info='nodal power injection', - e_str='Cg@(pn - pl) - pg', - type='eq',) - self.lub = Constraint(name='lub', info='Line limits upper bound', - e_str='PTDF @ (pn - pl) - rate_a', - type='uq',) - self.llb = Constraint(name='llb', info='Line limits lower bound', - e_str='- PTDF @ (pn - pl) - rate_a', - type='uq',) + # TODO: add eqn to get aBus + self.aref = Constraint(name='aref', type='eq', + info='reference bus angle', + e_str='Cs@aBus',) + + self.pnb = Constraint(name='pnb', type='eq', + info='nodal power injection', + e_str='PTDF@(Cgi@pg - Cli@pd) - plf',) + # --- objective --- self.obj = Objective(name='tc', info='total cost', unit='$', - e_str='sum(c2 * pg**2 + c1 * pg + ug * c0)', sense='min',) - - -class DCOPF(DCOPFData, DCOPFModel): - """ - Standard DC optimal power flow (DCOPF). - - In this model, the bus injected power ``pn`` is used as internal variable - between generator output and load demand. - """ - - def __init__(self, system, config): - DCOPFData.__init__(self) - DCOPFModel.__init__(self, system, config) + self.obj.e_str = 'sum_squares(mul(c2, pg))' \ + '+ sum(mul(c1, pg))' \ + '+ sum(mul(ug, c0))' diff --git a/ams/routines/dcpf.py b/ams/routines/dcpf.py index 2088b6e6..cd20ec27 100644 --- a/ams/routines/dcpf.py +++ b/ams/routines/dcpf.py @@ -1,30 +1,35 @@ """ Power flow routines. """ -import logging # NOQA +import logging -from andes.shared import deg2rad # NOQA -from andes.utils.misc import elapsed # NOQA +from andes.shared import deg2rad +from andes.utils.misc import elapsed -from ams.routines.routine import RoutineData, RoutineModel # NOQA -from ams.opt.omodel import Var # NOQA -from ams.pypower import runpf # NOQA -from ams.pypower.core import ppoption # NOQA +from ams.routines.routine import RoutineModel +from ams.opt.omodel import Var +from ams.pypower import runpf +from ams.pypower.core import ppoption -from ams.io.pypower import system2ppc # NOQA -from ams.core.param import RParam # NOQA +from ams.io.pypower import system2ppc +from ams.core.param import RParam logger = logging.getLogger(__name__) -class DCPFlowData(RoutineData): +class DCPFlowBase(RoutineModel): """ - Data class for power flow routines. + Base class for power flow. + + Overload the ``solve``, ``unpack``, and ``run`` methods. """ - def __init__(self): - RoutineData.__init__(self) - # --- line --- + def __init__(self, system, config): + RoutineModel.__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.', @@ -44,19 +49,6 @@ def __init__(self): unit='p.u.', model='mats', src='pl') - -class DCPFlowBase(RoutineModel): - """ - Base class for power flow. - - Overload the ``solve``, ``unpack``, and ``run`` methods. - """ - - def __init__(self, system, config): - RoutineModel.__init__(self, system, config) - self.info = 'DC Power Flow' - self.type = 'PF' - def unpack(self, res): """ Unpack results from PYPOWER. @@ -174,9 +166,15 @@ def disable(self, name): raise NotImplementedError -class DCPFlowModel(DCPFlowBase): +class DCPF(DCPFlowBase): """ - Base class for power flow model. + DC power flow. + + 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): @@ -193,19 +191,3 @@ def __init__(self, system, config): unit='p.u.', name='pg', tex_name=r'p_{g}', model='StaticGen', src='p',) - - -class DCPF(DCPFlowData, DCPFlowModel): - """ - DC power flow. - - 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=None, config=None): - DCPFlowData.__init__(self) - DCPFlowModel.__init__(self, system, config) diff --git a/ams/routines/dopf.py b/ams/routines/dopf.py index 812e9571..298e453c 100644 --- a/ams/routines/dopf.py +++ b/ams/routines/dopf.py @@ -1,67 +1,63 @@ """ Distributional optimal power flow (DOPF). """ -import numpy as np # NOQA +import numpy as np -from ams.core.param import RParam # NOQA -from ams.core.service import NumOp # NOQA +from ams.core.param import RParam -from ams.routines.dcopf import DCOPFData, DCOPFModel # NOQA +from ams.routines.dcopf import DCOPF -from ams.opt.omodel import Var, Constraint, Objective # NOQA +from ams.opt.omodel import Var, Constraint, Objective -class DOPFData(DCOPFData): +class DOPF(DCOPF): """ - Data for DOPF. + Linearzied distribution OPF, where power loss are ignored. + + Reference: + + [1] L. Bai, J. Wang, C. Wang, C. Chen, and F. Li, “Distribution Locational Marginal Pricing (DLMP) + for Congestion Management and Voltage Support,” IEEE Trans. Power Syst., vol. 33, no. 4, + pp. 4061–4073, Jul. 2018, doi: 10.1109/TPWRS.2017.2767632. """ - def __init__(self): - DCOPFData.__init__(self) - self.ql = RParam(info='reactive power demand connected to Bus (system base)', - name='ql', tex_name=r'q_{l}', unit='p.u.', - model='mats', src='ql',) + def __init__(self, system, config): + DCOPF.__init__(self, system, config) + self.info = 'Linearzied distribution OPF' + self.type = 'DED' + + # --- params --- + self.qd = RParam(info='reactive demand', + name='qd', tex_name=r'q_{d}', unit='p.u.', + model='StaticLoad', src='q0',) self.vmax = RParam(info="Bus voltage upper limit", - name='vmax', tex_name=r'v_{max}', unit='p.u.', - model='Bus', src='vmax', + name='vmax', tex_name=r'\overline{v}', + unit='p.u.', + model='Bus', src='vmax', no_parse=True, ) self.vmin = RParam(info="Bus voltage lower limit", - name='vmin', tex_name=r'v_{min}', unit='p.u.', - model='Bus', src='vmin', ) + name='vmin', tex_name=r'\underline{v}', + unit='p.u.', + model='Bus', src='vmin', no_parse=True,) self.r = RParam(info='line resistance', name='r', tex_name='r', unit='p.u.', model='Line', src='r') - self.x = RParam(info='line reactance', - name='x', tex_name='x', unit='p.u.', - model='Line', src='x', ) - self.qmax = RParam(info='generator maximum reactive power (system base)', + self.qmax = RParam(info='generator maximum reactive power', name='qmax', tex_name=r'q_{max}', unit='p.u.', model='StaticGen', src='qmax',) - self.qmin = RParam(info='generator minimum reactive power (system base)', + self.qmin = RParam(info='generator minimum reactive power', name='qmin', tex_name=r'q_{min}', unit='p.u.', model='StaticGen', src='qmin',) - self.Cft = RParam(info='connection matrix for Line and Bus', - name='Cft', tex_name=r'C_{ft}', - model='mats', src='Cft',) - - -class LDOPFModel(DCOPFModel): - """ - Linearzied distribution OPF model. - """ - - def __init__(self, system, config): - DCOPFModel.__init__(self, system, config) - self.info = 'Linearzied distribution OPF' - self.type = 'DED' # --- vars --- - self.qg = Var(info='Gen reactive power (system base)', + self.qg = Var(info='Gen reactive power', name='qg', tex_name=r'q_{g}', unit='p.u.', - model='StaticGen', src='q', - lb=self.qmin, ub=self.qmax,) - self.qn = Var(info='Bus reactive power', - name='qn', tex_name=r'q_{n}', unit='p.u.', - model='Bus',) + model='StaticGen', src='q',) + self.qglb = Constraint(name='qglb', type='uq', + info='qg min', + e_str='-qg + mul(ug, qmin)',) + self.qgub = Constraint(name='qgub', type='uq', + info='qg max', + e_str='qg - mul(ug, qmax)',) self.vsq = Var(info='square of Bus voltage', name='vsq', tex_name=r'v^{2}', unit='p.u.', @@ -75,32 +71,19 @@ def __init__(self, system, config): e_str='-vsq + vmin**2', type='uq',) - self.plf = Var(info='line active power', - name='plf', tex_name=r'p_{lf}', unit='p.u.', - model='Line',) self.qlf = Var(info='line reactive power', - name='qlf', tex_name=r'q_{lf}', unit='p.u.', - model='Line',) + name='qlf', tex_name=r'q_{lf}', + unit='p.u.', model='Line',) # --- constraints --- - self.CftT = NumOp(u=self.Cft, - fun=np.transpose, - name='CftT', - tex_name=r'C_{ft}^{T}', - info='transpose of connection matrix',) - self.pinj.e_str = 'CftT@plf - pl - pn' - self.qinj = Constraint(name='qinj', - info='node reactive power injection', - e_str='CftT@qlf - ql - qn', - type='eq',) + self.pnb.e_str = 'PTDF@(Cgi@pg - Cli@pd) - plf' self.qb = Constraint(name='qb', info='reactive power balance', - e_str='sum(ql) - sum(qg)', + e_str='sum(qd) - sum(qg)', type='eq',) - self.lvd = Constraint(name='lvd', info='line voltage drop', - e_str='Cft@vsq - (r * pl + x * qlf)', + e_str='Cft@vsq - (r * plf + x * qlf)', type='eq',) # --- objective --- @@ -112,9 +95,12 @@ def unpack(self, **kwargs): self.system.Bus.set(src='v', attr='v', value=vBus, idx=self.vsq.get_idx()) -class LDOPF(DOPFData, LDOPFModel): +class DOPFVIS(DOPF): """ - Linearzied distribution OPF, where power loss are ignored. + Linearzied distribution OPF with variables for virtual inertia and damping from from REGCV1, + where power loss are ignored. + + UNDER DEVELOPMENT! Reference: @@ -124,17 +110,9 @@ class LDOPF(DOPFData, LDOPFModel): """ def __init__(self, system, config): - DOPFData.__init__(self) - LDOPFModel.__init__(self, system, config) + DOPF.__init__(self, system, config) - -class DOPF2Data(DOPFData): - """ - Data for DOPF with PFRCost for virtual inertia and damping from REGCV1. - """ - - def __init__(self): - DOPFData.__init__(self) + # --- params --- self.cm = RParam(info='Virtual inertia cost', name='cm', src='cm', tex_name=r'c_{m}', unit=r'$/s', @@ -145,15 +123,6 @@ def __init__(self): tex_name=r'c_{d}', unit=r'$/(p.u.)', model='REGCV1Cost', indexer='reg', imodel='REGCV1',) - - -class LDOPF2Model(LDOPFModel): - """ - Linearzied distribution OPF model with VSG. - """ - - def __init__(self, system, config): - LDOPFModel.__init__(self, system, config) # --- vars --- self.M = Var(info='Emulated startup time constant (M=2H) from REGCV1', name='M', tex_name=r'M', unit='s', @@ -161,24 +130,8 @@ def __init__(self, system, config): self.D = Var(info='Emulated damping coefficient from REGCV1', name='D', tex_name=r'D', unit='p.u.', model='REGCV1',) + obj = 'sum(c2 * pg**2 + c1 * pg + ug * c0 + cm * M + cd * D)' self.obj = Objective(name='tc', info='total cost', unit='$', - e_str='sum(c2 * pg**2 + c1 * pg + ug * c0 + cm * M + cd * D)', + e_str=obj, sense='min',) - - -class LDOPF2(DOPF2Data, LDOPF2Model): - """ - Linearzied distribution OPF with variables for virtual inertia and damping from from REGCV1, - where power loss are ignored. - - Reference: - - [1] L. Bai, J. Wang, C. Wang, C. Chen, and F. Li, “Distribution Locational Marginal Pricing (DLMP) - for Congestion Management and Voltage Support,” IEEE Trans. Power Syst., vol. 33, no. 4, - pp. 4061–4073, Jul. 2018, doi: 10.1109/TPWRS.2017.2767632. - """ - - def __init__(self, system, config): - DOPF2Data.__init__(self) - LDOPF2Model.__init__(self, system, config) diff --git a/ams/routines/ed.py b/ams/routines/ed.py index 0d052bfa..88b4e8f2 100644 --- a/ams/routines/ed.py +++ b/ams/routines/ed.py @@ -1,67 +1,133 @@ """ -Real-time economic dispatch. +Economic dispatch routines. """ -import logging # NOQA -from collections import OrderedDict # NOQA -import numpy as np # NOQA +import logging +from collections import OrderedDict +import numpy as np -from ams.core.param import RParam # NOQA -from ams.core.service import (ZonalSum, NumOpDual, NumHstack, - RampSub, NumOp, LoadScale) # NOQA +from ams.core.param import RParam +from ams.core.service import (NumOpDual, NumHstack, + RampSub, NumOp, LoadScale) -from ams.routines.rted import RTEDData, ESD1Base # NOQA -from ams.routines.dcopf import DCOPFModel # NOQA +from ams.routines.rted import RTED, ESD1Base -from ams.core.service import VarSelect # NOQA -from ams.opt.omodel import Var, Constraint # NOQA +from ams.opt.omodel import Var, Constraint logger = logging.getLogger(__name__) -class EDData(RTEDData): +class SRBase: """ - Economic dispatch data. + Base class for spinning reserve. """ + + def __init__(self) -> None: + self.dsrp = RParam(info='spinning reserve requirement in percentage', + name='dsr', tex_name=r'd_{sr}', + model='SR', src='demand', + unit='%',) + self.csr = RParam(info='cost for spinning reserve', + name='csr', tex_name=r'c_{sr}', + model='SRCost', src='csr', + unit=r'$/(p.u.*h)', + indexer='gen', imodel='StaticGen',) - def __init__(self): - RTEDData.__init__(self) + self.prs = Var(name='prs', tex_name=r'p_{r,s}', + info='spinning reserve', unit='p.u.', + model='StaticGen', nonneg=True,) + + self.dsrpz = NumOpDual(u=self.pdz, u2=self.dsrp, fun=np.multiply, + name='dsrpz', tex_name=r'd_{s,r, p, z}', + info='zonal spinning reserve requirement in percentage',) + self.dsr = NumOpDual(u=self.dsrpz, u2=self.sd, fun=np.multiply, + rfun=np.transpose, + name='dsr', tex_name=r'd_{s,r,z}', + info='zonal spinning reserve requirement',) + + # NOTE: define e_str in dispatch model + self.prsb = Constraint(info='spinning reserve balance', + name='prsb', type='eq',) + self.rsr = Constraint(info='spinning reserve requirement', + name='rsr', type='uq',) + + +class MPBase: + """ + Base class for multi-period dispatch. + """ + + def __init__(self) -> None: # NOTE: Setting `ED.scale.owner` to `Horizon` will cause an error when calling `ED.scale.v`. # This is because `Horizon` is a group that only contains the model `TimeSlot`. # The `get` method of `Horizon` calls `andes.models.group.GroupBase.get` and results in an error. self.sd = RParam(info='zonal load factor for ED', name='sd', tex_name=r's_{d}', src='sd', model='EDTSlot') + self.timeslot = RParam(info='Time slot for multi-period ED', name='timeslot', tex_name=r't_{s,idx}', - src='idx', model='EDTSlot') - self.R30 = RParam(info='30-min ramp rate (system base)', + src='idx', model='EDTSlot', + no_parse=True) + + self.pdsz = NumOpDual(u=self.sd, u2=self.pdz, + fun=np.multiply, rfun=np.transpose, + name='pds', tex_name=r'p_{d,s,z}', + unit='p.u.', + info='Scaled zonal total load') + + self.tlv = NumOp(u=self.timeslot, fun=np.ones_like, + args=dict(dtype=float), + expand_dims=0, + name='tlv', tex_name=r'1_{tl}', + info='time length vector', + no_parse=True) + + self.pds = LoadScale(u=self.pd, sd=self.sd, + name='pds', tex_name=r'p_{d,s}', + info='Scaled load',) + + self.R30 = RParam(info='30-min ramp rate', name='R30', tex_name=r'R_{30}', src='R30', unit='p.u./min', model='StaticGen') - self.dsrp = RParam(info='spinning reserve requirement in percentage', - name='dsr', tex_name=r'd_{sr}', - model='SR', src='demand', - unit='%',) - self.csr = RParam(info='cost for spinning reserve', - name='csr', tex_name=r'c_{sr}', - model='SRCost', src='csr', - unit=r'$/(p.u.*h)', - indexer='gen', imodel='StaticGen',) - self.Cl = RParam(info='connection matrix for Load and Bus', - name='Cl', tex_name=r'C_{l}', - model='mats', src='Cl',) - self.zl = RParam(info='zone of load', - name='zl', tex_name=r'z_{l}', - model='StaticLoad', src='zone',) + self.Mr = RampSub(u=self.pg, name='Mr', tex_name=r'M_{r}', + info='Subtraction matrix for ramping',) + self.RR30 = NumHstack(u=self.R30, ref=self.Mr, + name='RR30', tex_name=r'R_{30,R}', + info='Repeated ramp rate',) + self.ctrl.expand_dims = 1 + self.c0.expand_dims = 1 + self.pmax.expand_dims = 1 + self.pmin.expand_dims = 1 + self.pg0.expand_dims = 1 + self.rate_a.expand_dims = 1 -class EDModel(DCOPFModel): +class ED(RTED): """ - Economic dispatch model. + DC-based multi-period economic dispatch (ED). + Dispath interval ``config.t`` (:math:`T_{cfg}`) is introduced, + 1 [Hour] by default. + + ED extends DCOPF as follows: + + 1. Vars ``pg``, ``pru``, ``prd`` are extended to 2D + + 2. 2D Vars ``rgu`` and ``rgd`` are introduced + + 3. Param ``ug`` is sourced from ``EDTSlot.ug`` as commitment decisions + + Notes + ----- + 1. Formulations has been adjusted with interval ``config.t`` + + 2. The tie-line flow is not implemented in this model. """ def __init__(self, system, config): - DCOPFModel.__init__(self, system, config) + RTED.__init__(self, system, config) + MPBase.__init__(self) + SRBase.__init__(self) self.config.add(OrderedDict((('t', 1), ))) @@ -72,142 +138,151 @@ def __init__(self, system, config): self.info = 'Economic dispatch' self.type = 'DCED' + self.ug.info = 'unit commitment decisions' + self.ug.model = 'EDTSlot' + self.ug.src = 'ug' + # self.ug.tex_name = r'u_{g}', + + self.dud.expand_dims = 1 + self.ddd.expand_dims = 1 + + # --- params --- + self.ugt = NumOp(u=self.ug, fun=np.transpose, + name='ugt', tex_name=r'u_{g}', + info='input ug transpose', + no_parse=True) + # --- vars --- # NOTE: extend pg to 2D matrix, where row is gen and col is timeslot self.pg.horizon = self.timeslot - self.pg.info = '2D power generation (system base)' + self.pg.info = '2D Gen power' + 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 + + self.plf.horizon = self.timeslot + self.plf.info = '2D Line flow' + + self.pru.horizon = self.timeslot + self.pru.info = '2D RegUp power' - self.pn.horizon = self.timeslot - self.pn.info = '2D Bus power injection (system base)' + self.prd.horizon = self.timeslot + self.prd.info = '2D RegDn power' + + self.prs.horizon = self.timeslot # --- constraints --- - # --- power balance --- - self.ds = ZonalSum(u=self.zb, zone='Region', - name='ds', tex_name=r'S_{d}', - info='Sum pl vector in shape of zone',) - self.pdz = NumOpDual(u=self.ds, u2=self.pl, - fun=np.multiply, - rfun=np.sum, rargs=dict(axis=1), - expand_dims=0, - name='pdz', tex_name=r'p_{d,z}', - unit='p.u.', info='zonal load') - self.pds = NumOpDual(u=self.sd, u2=self.pdz, - fun=np.multiply, rfun=np.transpose, - name='pds', tex_name=r'p_{d,s,t}', - unit='p.u.', info='Scaled total load as row vector') - self.gs = ZonalSum(u=self.zg, zone='Region', - name='gs', tex_name=r'S_{g}', - info='Sum Gen vars vector in shape of zone') # NOTE: Spg @ pg returns a row vector - self.pb.e_str = '- gs @ pg + pds' # power balance - - # spinning reserve - self.Rpmax = NumHstack(u=self.pmax, ref=self.timeslot, - name='Rpmax', tex_name=r'p_{max, R}', - info='Repetated pmax',) - self.Rug = NumHstack(u=self.ug, ref=self.timeslot, - name='Rug', tex_name=r'u_{g, R}', - info='Repetated ug',) - self.dsrpz = NumOpDual(u=self.pdz, u2=self.dsrp, fun=np.multiply, - name='dsrpz', tex_name=r'd_{sr, p, z}', - info='zonal spinning reserve requirement in percentage',) - self.dsr = NumOpDual(u=self.dsrpz, u2=self.sd, fun=np.multiply, - rfun=np.transpose, - name='dsr', tex_name=r'd_{sr}', - info='zonal spinning reserve requirement',) - self.sr = Constraint(name='sr', info='spinning reserve', type='uq', - e_str='-gs@multiply(Rpmax - pg, Rug) + dsr') + self.pb.e_str = '- gs @ pg + pdsz' # power balance + + self.prsb.e_str = 'mul(ugt, mul(pmax, tlv) - pg) - prs' + self.rsr.e_str = '-gs@prs + dsr' # --- bus power injection --- - self.Cli = NumOp(u=self.Cl, fun=np.linalg.pinv, - name='Cli', tex_name=r'C_{l}^{-1}', - info='inverse of Cl',) - self.Rpd = LoadScale(u=self.zl, sd=self.sd, Cl=self.Cl, - name='Rpd', tex_name=r'p_{d,R}', - info='Scaled nodal load',) - self.pinj.e_str = 'Cg @ (pn - Rpd) - pg' # power injection + self.pnb.e_str = 'PTDF@(Cgi@pg - Cli@pds) - plf' # --- line limits --- - self.RRA = NumHstack(u=self.rate_a, ref=self.timeslot, - name='RRA', tex_name=r'R_{ATEA,R}', - info='Repeated rate_a',) - self.lub.e_str = 'PTDF @ (pn - Rpd) - RRA' - self.llb.e_str = '-PTDF @ (pn - Rpd) - RRA' + self.plflb.e_str = '-plf - mul(rate_a, tlv)' + self.plfub.e_str = 'plf - mul(rate_a, tlv)' # --- ramping --- - self.Mr = RampSub(u=self.pg, name='Mr', tex_name=r'M_{r}', - info='Subtraction matrix for ramping, (ng, ng-1)',) - self.RR30 = NumHstack(u=self.R30, ref=self.Mr, - name='RR30', tex_name=r'R_{30,R}', - info='Repeated ramp rate as 2D matrix, (ng, ng-1)',) - self.rgu = Constraint(name='rgu', info='Gen ramping up', - e_str='pg @ Mr - t dot RR30', - type='uq',) - self.rgd = Constraint(name='rgd', - info='Gen ramping down', - e_str='-pg @ Mr - t dot RR30', - type='uq',) + self.rbu.e_str = 'gs@mul(ugt, pru) - mul(dud, tlv)' + self.rbd.e_str = 'gs@mul(ugt, prd) - mul(ddd, tlv)' + + self.rru.e_str = 'mul(ugt, pg + pru) - mul(pmax, tlv)' + self.rrd.e_str = 'mul(ugt, -pg + prd) - mul(pmin, tlv)' + + self.rgu.e_str = 'pg @ Mr - t dot RR30' + self.rgd.e_str = '-pg @ Mr - t dot RR30' + self.rgu0 = Constraint(name='rgu0', info='Initial gen ramping up', - e_str='pg[:, 0] - pg0 - R30', + e_str='pg[:, 0] - pg0[:, 0] - R30', type='uq',) self.rgd0 = Constraint(name='rgd0', info='Initial gen ramping down', - e_str='- pg[:, 0] + pg0 - R30', + e_str='- pg[:, 0] + pg0[:, 0] - R30', type='uq',) # --- objective --- - # NOTE: no need to fix objective function - gcost = 'sum(c2 @ (t dot pg)**2 + c1 @ (t dot pg) + ug * c0)' - rcost = ' + sum(csr * ug * (Rpmax - pg))' - self.obj.e_str = gcost + rcost + cost = 'sum(c2 @ (t dot pg)**2 + c1 @ (t dot pg))' + # constant cost + cost += '+ sum(mul(ugt, mul(c0, tlv)))' + # spinning reserve cost + cost += ' + sum(csr@prs)' + self.obj.e_str = cost + + def dc2ac(self, **kwargs): + """ + AC conversion ``dc2ac`` is not implemented yet for + multi-period dispatch. + """ + return NotImplementedError def unpack(self, **kwargs): """ - ED will not unpack results from solver into devices - because the resutls are multi-time-period. + Multi-period dispatch will not unpack results from + solver into devices. + + # TODO: unpack first period results, and allow input + # to specify which period to unpack. """ return None - -class ED(EDData, EDModel): - """ - DC-based multi-period economic dispatch (ED). - Dispath interval ``config.t`` (:math:`T_{cfg}`) is introduced, - 1 [Hour] by default. - - ED extends DCOPF as follows: - - 1. Var ``pg`` is extended to 2D - - 2. 2D Vars ``rgu`` and ``rgd`` are introduced - - Notes - ----- - 1. Formulations has been adjusted with interval ``config.t`` - - 2. The tie-line flow is not implemented in this model. - """ - - def __init__(self, system, config): - EDData.__init__(self) - EDModel.__init__(self, system, config) - # TODO: add data check # if has model ``TimeSlot``, mandatory # if has model ``Region``, optional # if ``Region``, if ``Bus`` has param ``zone``, optional, if none, auto fill -class ED2(EDData, EDModel, ESD1Base): +class ESD1MPBase(ESD1Base): + """ + Extended base class for energy storage in multi-period dispatch. + """ + + def __init__(self): + ESD1Base.__init__(self) + + self.Mre = RampSub(u=self.SOC, name='Mre', tex_name=r'M_{r,E}', + info='Subtraction matrix for SOC', + no_parse=True) + self.EnR = NumHstack(u=self.En, ref=self.Mre, + name='EnR', tex_name=r'E_{n,R}', + info='Repeated En as 2D matrix, (ng, ng-1)') + self.EtaCR = NumHstack(u=self.EtaC, ref=self.Mre, + name='EtaCR', tex_name=r'\eta_{c,R}', + info='Repeated Etac as 2D matrix, (ng, ng-1)') + self.REtaDR = NumHstack(u=self.REtaD, ref=self.Mre, + name='REtaDR', tex_name=r'R_{\eta_d,R}', + info='Repeated REtaD as 2D matrix, (ng, ng-1)') + SOCb = 'mul(EnR, SOC @ Mre) - t dot mul(EtaCR, zce[:, 1:])' + SOCb += ' + t dot mul(REtaDR, zde[:, 1:])' + self.SOCb.e_str = SOCb + + SOCb0 = 'mul(En, SOC[:, 0] - SOCinit) - t dot mul(EtaC, zce[:, 0])' + SOCb0 += ' + t dot mul(REtaD, zde[:, 0])' + self.SOCb0 = Constraint(name='SOCb', type='eq', + info='ESD1 SOC initial balance', + e_str=SOCb0,) + + self.SOCr = Constraint(name='SOCr', type='eq', + info='SOC requirement', + e_str='SOC[:, -1] - SOCinit',) + + +class EDES(ED, ESD1Base): """ ED with energy storage :ref:`ESD1`. The bilinear term in the formulation is linearized with big-M method. """ def __init__(self, system, config): - EDData.__init__(self) - EDModel.__init__(self, system, config) + ED.__init__(self, system, config) ESD1Base.__init__(self) self.config.t = 1 # dispatch interval in hour @@ -215,7 +290,11 @@ def __init__(self, system, config): self.info = 'Economic dispatch with energy storage' self.type = 'DCED' + # NOTE: extend vars to 2D self.SOC.horizon = self.timeslot - self.pge.horizon = self.timeslot - self.ued.horizon = self.timeslot - self.zue.horizon = self.timeslot + self.pce.horizon = self.timeslot + self.pde.horizon = self.timeslot + self.uce.horizon = self.timeslot + self.ude.horizon = self.timeslot + self.zce.horizon = self.timeslot + self.zde.horizon = self.timeslot diff --git a/ams/routines/pflow.py b/ams/routines/pflow.py index 1cb66170..9e689987 100644 --- a/ams/routines/pflow.py +++ b/ams/routines/pflow.py @@ -9,19 +9,29 @@ from ams.pypower.core import ppoption # NOQA from ams.core.param import RParam # NOQA -from ams.routines.dcpf import DCPFlowData, DCPFlowBase # NOQA +from ams.routines.dcpf import DCPFlowBase # NOQA from ams.opt.omodel import Var, Constraint # NOQA logger = logging.getLogger(__name__) -class PFlowData(DCPFlowData): +class PFlow(DCPFlowBase): """ 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. """ - def __init__(self): - DCPFlowData.__init__(self) + def __init__(self, system, config): + DCPFlowBase.__init__(self, system, config) + self.info = "AC Power Flow" + self.type = "PF" + self.qd = RParam( info="reactive power load in system base", name="qd", @@ -31,17 +41,6 @@ def __init__(self): model="PQ", ) - -class PFlowModel(DCPFlowBase): - """ - AC power flow model. - """ - - def __init__(self, system, config): - DCPFlowBase.__init__(self, system, config) - self.info = "AC Power Flow" - self.type = "PF" - # --- bus --- self.aBus = Var( info="bus voltage angle", @@ -133,20 +132,3 @@ def run(self, force_init=False, no_code=True, method="newton", **kwargs): method=method, **kwargs, ) - - -class PFlow(PFlowData, PFlowModel): - """ - 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. - """ - - def __init__(self, system=None, config=None): - PFlowData.__init__(self) - PFlowModel.__init__(self, system, config) diff --git a/ams/routines/routine.py b/ams/routines/routine.py index 6137358c..f891b133 100644 --- a/ams/routines/routine.py +++ b/ams/routines/routine.py @@ -2,25 +2,26 @@ Module for routine data. """ -import logging # NOQA -from typing import Optional, Union, Type, Iterable # NOQA -from collections import OrderedDict # NOQA +import logging +from typing import Optional, Union, Type, Iterable +from collections import OrderedDict -import numpy as np # NOQA +import numpy as np +import matplotlib.pyplot as plt -from andes.core import Config # NOQA +from andes.core import Config from andes.shared import deg2rad # NOQA -from andes.utils.misc import elapsed # NOQA -from ams.utils import timer # NOQA -from ams.core.param import RParam # NOQA -from ams.opt.omodel import OModel, Var, Constraint, Objective # NOQA +from andes.utils.misc import elapsed -from ams.core.symprocessor import SymProcessor # NOQA -from ams.core.documenter import RDocumenter # NOQA -from ams.core.service import RBaseService, ValueService # NOQA +from ams.core.param import RParam +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 + +from ams.shared import igraph as ig +from ams.shared import require_igraph -from ams.models.group import GroupBase # NOQA -from ams.core.model import Model # NOQA logger = logging.getLogger(__name__) @@ -31,7 +32,7 @@ class RoutineData: """ def __init__(self): - self.rparams = OrderedDict() # list out RParam in a routine + pass class RoutineModel: @@ -43,19 +44,24 @@ def __init__(self, system=None, config=None): self.system = system self.config = Config(self.class_name) self.info = None - self.tex_names = OrderedDict((('sys_f', 'f_{sys}'), - ('sys_mva', 'S_{b,sys}'), - )) + self.tex_names = OrderedDict( + ( + ("sys_f", "f_{sys}"), + ("sys_mva", "S_{b,sys}"), + ) + ) self.syms = SymProcessor(self) # symbolic processor self._syms = False # flag if symbols has been generated + self.rparams = OrderedDict() # list out RParam in a routine self.services = OrderedDict() # list out services in a routine + self.params = OrderedDict() # list out Params in a routine self.vars = OrderedDict() # list out Vars in a routine self.constrs = OrderedDict() self.obj = None self.initialized = False - self.type = 'UndefinedType' + self.type = "UndefinedType" self.docum = RDocumenter(self) # --- sync mapping --- @@ -68,17 +74,24 @@ def __init__(self, system=None, config=None): if config is not None: self.config.load(config) # TODO: these default configs might to be revised - self.config.add(OrderedDict((('sparselib', 'klu'), - ('linsolve', 0), - ))) - self.config.add_extra("_help", - sparselib="linear sparse solver name", - linsolve="solve symbolic factorization each step (enable when KLU segfaults)", - ) - self.config.add_extra("_alt", - sparselib=("klu", "umfpack", "spsolve", "cupy"), - linsolve=(0, 1), - ) + self.config.add( + OrderedDict( + ( + ("sparselib", "klu"), + ("linsolve", 0), + ) + ) + ) + self.config.add_extra( + "_help", + sparselib="linear sparse solver name", + linsolve="solve symbolic factorization each step (enable when KLU segfaults)", + ) + self.config.add_extra( + "_alt", + sparselib=("klu", "umfpack", "spsolve", "cupy"), + linsolve=(0, 1), + ) self.exec_time = 0.0 # recorded time to execute the routine in seconds # TODO: check exit_code of gurobipy or any other similiar solvers @@ -100,7 +113,8 @@ def _loc(self, src: str, idx, allow_none=False): return loc else: idx_none = [idxe for idxe in idx if idxe not in src_idx] - raise ValueError(f'Var <{self.class_name}.{src}> does not contain value with idx={idx_none[0]}') + msg = f"Var <{self.class_name}.{src}> does not contain value with idx={idx_none}" + raise ValueError(msg) def get_load(self, horizon: Union[int, str], src: str, attr: str = 'v', @@ -126,22 +140,24 @@ def get_load(self, horizon: Union[int, str], pq_zone = self.system.PQ.zone.v pq0 = self.system.PQ.get(src=src, attr=attr, idx=idx) else: - pq_zone = self.system.PQ.get(src='zone', attr='v', idx=idx) + pq_zone = self.system.PQ.get(src="zone", attr="v", idx=idx) pq0 = self.system.PQ.get(src=src, attr=attr, idx=idx) col = [all_zone.index(pq_z) for pq_z in pq_zone] mdl = self.system.__dict__[model] if mdl.n == 0: - raise ValueError(f'<{model}> does not have data, check input file.') + raise ValueError(f"<{model}> does not have data, check input file.") if factor not in mdl.__dict__.keys(): - raise ValueError(f'<{model}> does not have <{factor}>.') + raise ValueError(f"<{model}> does not have <{factor}>.") sdv = mdl.__dict__[factor].v horizon_all = mdl.idx.v try: row = horizon_all.index(horizon) - except ValueError: - raise ValueError(f'<{model}> does not have horizon with idx=<{horizon}>.') + except ValueError as e: + msg = f"<{model}> does not have horizon with idx=<{horizon}>. " + msg += f"Original error: {e}" + raise ValueError(msg) pq_factor = np.array(sdv[:, col][row, :]) pqv = np.multiply(pq0, pq_factor) return pqv @@ -163,16 +179,16 @@ def get(self, src: str, idx, attr: str = 'v', Horizon index. """ if src not in self.__dict__.keys(): - raise ValueError(f'<{src}> does not exist in <<{self.class_name}>.') + raise ValueError(f"<{src}> does not exist in <<{self.class_name}>.") item = self.__dict__[src] if not hasattr(item, attr): - raise ValueError(f'{attr} does not exist in {self.class_name}.{src}.') + raise ValueError(f"{attr} does not exist in {self.class_name}.{src}.") idx_all = item.get_idx() if idx_all is None: - raise ValueError(f'<{self.class_name}> item <{src}> has no idx.') + raise ValueError(f"<{self.class_name}> item <{src}> has no idx.") if isinstance(idx, (str, int)): idx = [idx] @@ -183,12 +199,13 @@ def get(self, src: str, idx, attr: str = 'v', loc = [idx_all.index(idxe) if idxe in idx_all else None for idxe in idx] if None in loc: idx_none = [idxe for idxe in idx if idxe not in idx_all] - raise ValueError(f'Var <{self.class_name}.{src}> does not contain value with idx={idx_none}') + msg = f"Var <{self.class_name}.{src}> does not contain value with idx={idx_none}" + raise ValueError(msg) out = getattr(item, attr)[loc] if horizon is not None: if item.horizon is None: - raise ValueError(f'horizon is not defined for {self.class_name}.{src}.') + raise ValueError(f"horizon is not defined for {self.class_name}.{src}.") horizon_all = item.horizon.get_idx() if isinstance(horizon, int): horizon = [horizon] @@ -197,16 +214,20 @@ def get(self, src: str, idx, attr: str = 'v', if isinstance(horizon, np.ndarray): horizon = horizon.tolist() if isinstance(horizon, list): - loc_h = [horizon_all.index(idxe) if idxe in horizon_all else None for idxe in horizon] + loc_h = [ + horizon_all.index(idxe) if idxe in horizon_all else None + for idxe in horizon + ] if None in loc_h: idx_none = [idxe for idxe in horizon if idxe not in horizon_all] - raise ValueError(f'Var <{self.class_name}.{src}> does not contain horizon with idx={idx_none}') + msg = f"Var <{self.class_name}.{src}> does not contain horizon with idx={idx_none}" + raise ValueError(msg) out = out[:, loc_h] if out.shape[1] == 1: out = out[:, 0] return out - def set(self, src: str, idx, attr: str = 'v', value=0.0): + def set(self, src: str, idx, attr: str = "v", value=0.0): """ Set the value of an attribute of a routine parameter. """ @@ -215,11 +236,11 @@ def set(self, src: str, idx, attr: str = 'v', value=0.0): owner = self.__dict__[src].owner return owner.set(src=src, idx=idx, attr=attr, value=value) else: - logger.info(f'Variable {self.name} has no owner.') + logger.info(f"Variable {self.name} has no owner.") # FIXME: add idx for non-grouped variables return None - def doc(self, max_width=78, export='plain'): + def doc(self, max_width=78, export="plain"): """ Retrieve routine documentation as a string. """ @@ -248,13 +269,18 @@ def _data_check(self): if rparam.owner.n == 0: no_input.append(rname) owner_list.append(rparam.owner.class_name) + # TODO: add more data config check? + if rparam.config.pos: + if not np.all(rparam.v > 0): + logger.warning(f"RParam <{rname}> should have all positive values.") if len(no_input) > 0: - logger.error(f"Following models are missing from input file: {set(owner_list)}") + msg = f"Following models are missing in input: {set(owner_list)}" + logger.warning(msg) return False # TODO: add data validation for RParam, typical range, etc. return True - def init(self, force=False, no_code=True, **kwargs): + def init(self, force=True, no_code=True, **kwargs): """ Setup optimization model. @@ -265,16 +291,16 @@ def init(self, force=False, no_code=True, **kwargs): no_code: bool Whether to show generated code. """ - # TODO: add input check, e.g., if GCost exists if not force and self.initialized: - logger.debug(f'{self.class_name} has already been initialized.') + logger.debug(f"{self.class_name} has already been initialized.") return True if self._data_check(): - logger.debug(f'{self.class_name} data check passed.') + logger.debug(f"{self.class_name} data check passed.") else: - logger.warning(f'{self.class_name} data check failed, setup may run into error!') + msg = f"{self.class_name} data check failed, setup may run into error!" + logger.warning(msg) self._constr_check() - # FIXME: build the system matrices every init might slow down the process + # FIXME: build the system matrices every init might slow down self.system.mats.make() results, elapsed_time = self.om.setup(no_code=no_code) common_msg = f"Routine <{self.class_name}> " @@ -332,9 +358,12 @@ def run(self, force_init=False, no_code=True, **kwargs): self.exit_code = self.syms.status[status] self.system.exit_code = self.exit_code _, s = elapsed(t0) - self.exec_time = float(s.split(' ')[0]) + self.exec_time = float(s.split(" ")[0]) sstats = self.om.mdl.solver_stats # solver stats - n_iter = int(sstats.num_iters) + if sstats.num_iters is None: + n_iter = -1 + else: + 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 as {status} in {s}, converged after " @@ -343,7 +372,7 @@ def run(self, force_init=False, no_code=True, **kwargs): self.unpack(**kwargs) return True else: - msg = f"{self.class_name} failed after " + msg = f"{self.class_name} failed as {status} after " msg += n_iter_str + f"using solver {sstats.solver_name}!" logger.warning(msg) return False @@ -361,7 +390,7 @@ def report(self, **kwargs): raise NotImplementedError def __repr__(self): - return f'{self.class_name} at {hex(id(self))}' + return f"{self.class_name} at {hex(id(self))}" def _ppc2ams(self): """ @@ -383,11 +412,12 @@ def _check_attribute(self, key, value): """ if key in self.__dict__: existing_keys = [] - for type in ['constrs', 'vars', 'rparams']: + for type in ["constrs", "vars", "rparams", "services"]: if type in self.__dict__: existing_keys += list(self.__dict__[type].keys()) if key in existing_keys: - logger.warning(f"{self.class_name}: redefinition of member <{key}>. Likely a modeling error.") + msg = f"Attribute <{key}> already exists in <{self.class_name}>." + logger.warning(msg) # register owner routine instance of following attributes if isinstance(value, (RBaseService)): @@ -408,8 +438,6 @@ def __setattr__(self, key, value): # NOTE: value.id is not in use yet if isinstance(value, Var): value.id = len(self.vars) - elif isinstance(value, RParam): - value.id = len(self.rparams) self._check_attribute(key, value) self._register_attribute(key, value) @@ -422,17 +450,73 @@ 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, (Var, Constraint, Objective)): + if isinstance(value, (Param, Var, Constraint, Objective)): value.om = self.om + value.rtn = self + if isinstance(value, Param): + self.params[key] = value + self.om.params[key] = None # cp.Parameter if isinstance(value, Var): self.vars[key] = value + self.om.vars[key] = None # cp.Variable elif isinstance(value, Constraint): self.constrs[key] = value + self.om.constrs[key] = None # cp.Constraint elif isinstance(value, RParam): self.rparams[key] = value elif isinstance(value, RBaseService): self.services[key] = value + def update(self, params=None, mat_make=True,): + """ + Update the values of Parameters in the optimization model. + + This method is particularly important when some `RParams` are + linked with system matrices. + In such cases, setting `mat_make=True` is necessary to rebuild + these matrices for the changes to take effect. + This is common in scenarios involving topology changes, connection statuses, + or load value modifications. + If unsure, it is advisable to use `mat_make=True` as a precautionary measure. + + Parameters + ---------- + params: Parameter, str, or list + Parameter, Parameter name, or a list of parameter names to be updated. + If None, all parameters will be updated. + mat_make: bool + True to rebuild the system matrices. Set to False to speed up the process + if no system matrices are changed. + """ + t0, _ = elapsed() + re_setup = False + # sanitize input + sparams = [] + if params is None: + sparams = [val for val in self.params.values()] + mat_make = True + elif isinstance(params, Param): + sparams = [params] + elif isinstance(params, str): + sparams = [self.params[params]] + elif isinstance(params, list): + sparams = [self.params[param] for param in params if isinstance(param, str)] + for param in params: + param.update() + for param in sparams: + if param.optz is None: # means no_parse=True + re_setup = True + break + if mat_make: + self.system.mats.make() + if re_setup: + logger.warning(f"Resetup {self.class_name} OModel due to non-parametric change.") + _, _ = self.om.setup(no_code=True) + results = self.om.update(params=sparams) + t0, s0 = elapsed(t0) + logger.debug(f"Update params in {s0}.") + return results + def __delattr__(self, name): """ Overload the delattr function to unregister attributes. @@ -443,7 +527,7 @@ def __delattr__(self, name): name of the attribute """ self._unregister_attribute(name) - if name == 'obj': + if name == "obj": self.obj = None else: super().__delattr__(name) # Call the superclass implementation @@ -675,8 +759,6 @@ def addVars(self, info: Optional[str] = None, src: Optional[str] = None, unit: Optional[str] = None, - lb: Optional[str] = None, - ub: Optional[str] = None, horizon: Optional[RParam] = None, nonneg: Optional[bool] = False, nonpos: Optional[bool] = False, @@ -749,11 +831,15 @@ def addVars(self, """ if model is None and shape is None: raise ValueError("Either model or shape must be specified.") - item = Var(name=name, tex_name=tex_name, info=info, src=src, unit=unit, - model=model, shape=shape, lb=lb, ub=ub, horizon=horizon, nonneg=nonneg, - nonpos=nonpos, complex=complex, imag=imag, symmetric=symmetric, - diag=diag, psd=psd, nsd=nsd, hermitian=hermitian, bool=bool, - integer=integer, pos=pos, neg=neg, ) + item = Var(name=name, tex_name=tex_name, + info=info, src=src, unit=unit, + model=model, shape=shape, horizon=horizon, + nonneg=nonneg, nonpos=nonpos, + complex=complex, imag=imag, + symmetric=symmetric, diag=diag, + psd=psd, nsd=nsd, hermitian=hermitian, + bool=bool, integer=integer, + pos=pos, neg=neg, ) # add the variable as an routine attribute setattr(self, name, item) @@ -766,8 +852,10 @@ def addVars(self, elif item.model in self.system.models.keys(): item.owner = self.system.models[item.model] else: - msg = f'Model indicator \'{item.model}\' of <{item.rtn.class_name}.{name}>' - msg += ' is not a model or group. Likely a modeling error.' + msg = ( + f"Model indicator '{item.model}' of <{item.rtn.class_name}.{name}>" + ) + msg += " is not a model or group. Likely a modeling error." logger.warning(msg) self._post_add_check() @@ -779,3 +867,221 @@ def _initial_guess(self): Generate initial guess for the optimization model. """ raise NotImplementedError + + @require_igraph + def igmake(self, directed=True): + """ + Build an igraph object from the system. + + Parameters + ---------- + directed: bool + Whether the graph is directed. + + Returns + ------- + igraph.Graph + An igraph object. + """ + system = self.system + edges = np.column_stack([system.Bus.idx2uid(system.Line.bus1.v), + system.Bus.idx2uid(system.Line.bus2.v)]) + g = ig.Graph(n=system.Bus.n, directed=directed, edges=edges) + return g + + @require_igraph + def igraph( + self, + input: Optional[Union[RParam, Var]] = None, + ytimes: Optional[float] = None, + decimal: Optional[int] = 6, + directed: Optional[bool] = True, + dpi: Optional[int] = 100, + figsize: Optional[tuple] = None, + adjust_bus: Optional[bool] = False, + gen_color: Optional[str] = "red", + rest_color: Optional[str] = "black", + vertex_shape: Optional[str] = "circle", + vertex_font: Optional[str] = None, + no_vertex_label: Optional[bool] = False, + vertex_label: Optional[Union[str, list]] = None, + vertex_size: Optional[float] = None, + vertex_label_size: Optional[float] = None, + vertex_label_dist: Optional[float] = 1.5, + vertex_label_angle: Optional[float] = 10.2, + edge_arrow_size: Optional[float] = None, + edge_arrow_width: Optional[float] = None, + edge_width: Optional[float] = None, + edge_align_label: Optional[bool] = True, + edge_background: Optional[str] = None, + edge_color: Optional[str] = None, + edge_curved: Optional[bool] = False, + edge_font: Optional[str] = None, + edge_label: Optional[Union[str, list]] = None, + layout: Optional[str] = "rt", + autocurve: Optional[bool] = True, + ax: Optional[plt.Axes] = None, + title: Optional[str] = None, + title_loc: Optional[str] = None, + **visual_style, + ): + """ + Plot a system uging `g.plot()` of `igraph`, with optional input. + For now, only support plotting of Bus and Line elements as input. + + Examples + -------- + >>> import ams + >>> sp = ams.load(ams.get_case('5bus/pjm5bus_uced.xlsx')) + >>> sp.DCOPF.run() + >>> sp.DCOPF.plot(input=sp.DCOPF.pn, + >>> ytimes=10, + >>> adjust_bus=True, + >>> vertex_size=10, + >>> vertex_label_size=15, + >>> vertex_label_dist=2, + >>> vertex_label_angle=90, + >>> show=False, + >>> edge_align_label=True, + >>> autocurve=True,) + + Parameters + ---------- + input: RParam or Var, optional + The variable or parameter to be plotted. + ytimes: float, optional + The scaling factor of the values. + directed: bool, optional + Whether the graph is directed. + dpi: int, optional + Dots per inch. + figsize: tuple, optional + Figure size. + adjust_bus: bool, optional + Whether to adjust the bus size. + gen_color: str, optional + Color of the generator bus. + rest_color: str, optional + Color of the rest buses. + no_vertex_label: bool, optional + Whether to show vertex labels. + vertex_shape: str, optional + Shape of the vertices. + vertex_font: str, optional + Font of the vertices. + vertex_size: float, optional + Size of the vertices. + vertex_label_size: float, optional + Size of the vertex labels. + vertex_label_dist: float, optional + Distance of the vertex labels. + vertex_label_angle: float, optional + Angle of the vertex labels. + edge_arrow_size: float, optional + Size of the edge arrows. + edge_arrow_width: float, optional + Width of the edge arrows. + edge_width: float, optional + Width of the edges. + edge_align_label: bool, optional + Whether to align the edge labels. + edge_background: str, optional + RGB colored rectangle background of the edge labels. + layout: str, optional + Layout of the graph, ['rt', 'kk', 'fr', 'drl', 'lgl', 'circle', 'grid_fr']. + autocurve: bool, optional + Whether to use autocurve. + ax: plt.Axes, optional + Matplotlib axes. + visual_style: dict, optional + Visual style, see ``igraph.plot`` for details. + + Returns + ------- + plt.Axes + Matplotlib axes. + igraph.Graph + An igraph object. + """ + + g = self.igmake(directed=directed) + + # --- visual style --- + vstyle = { + # layout style + "layout": layout, + # vertices + "vertex_shape": vertex_shape, + "vertex_font": vertex_font, + "vertex_size": vertex_size, + "vertex_label": vertex_label, + "vertex_label_size": vertex_label_size, + "vertex_label_dist": vertex_label_dist, + "vertex_label_angle": vertex_label_angle, + # edges + "edge_arrow_size": edge_arrow_size, + "edge_arrow_width": edge_arrow_width, + "edge_width": edge_width, + "edge_align_label": edge_align_label, + "edge_background": edge_background, + "edge_color": edge_color, + "edge_curved": edge_curved, + "edge_font": edge_font, + "edge_label": edge_label, + # others + **visual_style, + } + system = self.system + # bus name, will be overwritten if input is not None + vstyle["vertex_name"] = system.Bus.name.v + if vertex_label is None: + vstyle["vertex_label"] = None if no_vertex_label else system.Bus.name.v + + # bus size + gidx = system.PV.idx.v + system.Slack.idx.v + gbus = system.StaticGen.get(src="bus", attr="v", idx=gidx) + # initialize all bus size as vertex_size + bus_size = [vertex_size] * system.Bus.n + if adjust_bus and isinstance(vertex_size, (int, float)): + # adjust gen bus size using Sn + gsn = system.StaticGen.get(src="Sn", attr="v", idx=gidx) + gbsize = vertex_size * gsn / gsn.max() + gbus_dict = {bus: size for bus, size in zip(gbus, gbsize)} + for key, val in gbus_dict.items(): + bus_size[system.Bus.idx2uid(key)] = val + if isinstance(vertex_size, Iterable): + bus_size = vertex_size + vstyle["vertex_size"] = bus_size + + # bus colors + gbus_uid = system.Bus.idx2uid(gbus) + bus_uid = system.Bus.idx2uid(system.Bus.idx.v) + g.vs["label"] = system.Bus.name.v + g.vs["bus_type"] = ["gen" if bus_i in gbus_uid else "rest" for bus_i in bus_uid] + color_dict = {"gen": gen_color, "rest": rest_color} + vstyle["vertex_color"] = [color_dict[btype] for btype in g.vs["bus_type"]] + + # --- variables --- + k = ytimes if ytimes is not None else 1 + if input is not None: + if input.owner.class_name == "Bus": + logger.debug(f"Plotting <{input.name}> as vertex label.") + values = [f"${input.tex_name}$={round(k*v, decimal)}" for v in input.v] + vstyle["vertex_label"] = values + elif input.owner.class_name == "Line": + logger.debug(f"Plotting <{input.name}> as edge label.") + values = [f"${input.tex_name}$={round(k*v, decimal)}" for v in input.v] + elabel = system.Line.name.v + eout = [f"{label}" for label, ein in zip(values, elabel)] + vstyle["edge_label"] = eout + else: + logger.error(f"Unsupported input type <{input.owner.class_name}>.") + + if ax is None: + _, ax = plt.subplots(figsize=figsize, dpi=dpi) + default_name = self.class_name + if input is not None: + default_name += f"\n${input.tex_name}$" + f" [${input.unit}$]" + ax.set_title(title if title else default_name, loc=title_loc) + ig.plot(g, autocurve=autocurve, target=ax, **vstyle) + return ax, g diff --git a/ams/routines/rted.py b/ams/routines/rted.py index 9ad82601..b463830a 100644 --- a/ams/routines/rted.py +++ b/ams/routines/rted.py @@ -7,23 +7,62 @@ from ams.core.param import RParam # NOQA from ams.core.service import ZonalSum, VarSelect, NumOp, NumOpDual # NOQA -from ams.routines.dcopf import DCOPFData, DCOPFModel # NOQA +from ams.routines.dcopf import DCOPFBase, DCOPF # NOQA from ams.opt.omodel import Var, Constraint # NOQA logger = logging.getLogger(__name__) -class RTEDData(DCOPFData): +class RTEDBase(DCOPF): """ - Data for real-time economic dispatch. + Base class for real-time economic dispatch (RTED). + + Overload ``dc2ac``, ``run``. + """ + + def __init__(self, system, config): + DCOPF.__init__(self, system, config) + # --- region --- + self.zg = RParam(info='Gen zone', + name='zg', tex_name='z_{one,g}', + model='StaticGen', src='zone', + no_parse=True) + self.zd = RParam(info='Load zone', + name='zd', tex_name='z_{one,d}', + model='StaticLoad', src='zone', + no_parse=True) + self.gs = ZonalSum(u=self.zg, zone='Region', + name='gs', tex_name=r'S_{g}', + info='Sum Gen vars vector in shape of zone') + self.ds = ZonalSum(u=self.zd, zone='Region', + name='ds', tex_name=r'S_{d}', + info='Sum pd vector in shape of zone', + no_parse=True,) + self.pdz = NumOpDual(u=self.ds, u2=self.pd, + fun=np.multiply, + rfun=np.sum, rargs=dict(axis=1), + expand_dims=0, + name='pdz', tex_name=r'p_{d,z}', + unit='p.u.', info='zonal total load', + no_parse=True,) + + # --- generator --- + self.R10 = RParam(info='10-min ramp rate', + name='R10', tex_name=r'R_{10}', + model='StaticGen', src='R10', + unit='p.u./h',) + + +class SFRBase: + """ + Base class for SFR used in DCED. """ def __init__(self): - DCOPFData.__init__(self) # 1. reserve - # 1.1. reserve cost + # --- reserve cost --- self.cru = RParam(info='RegUp reserve coefficient', name='cru', tex_name=r'c_{r,u}', model='SFRCost', src='cru', @@ -32,35 +71,50 @@ def __init__(self): name='crd', tex_name=r'c_{r,d}', model='SFRCost', src='crd', unit=r'$/(p.u.)',) - # 1.2. reserve requirement + # --- reserve requirement --- self.du = RParam(info='RegUp reserve requirement in percentage', name='du', tex_name=r'd_{u}', model='SFR', src='du', - unit='%',) + unit='%', no_parse=True,) self.dd = RParam(info='RegDown reserve requirement in percentage', name='dd', tex_name=r'd_{d}', model='SFR', src='dd', - unit='%',) - self.zb = RParam(info='Bus zone', - name='zb', tex_name='z_{one,bus}', - model='Bus', src='zone', ) - self.zg = RParam(info='generator zone data', - name='zg', tex_name='z_{one,g}', - model='StaticGen', src='zone',) - # 2. generator - self.R10 = RParam(info='10-min ramp rate (system base)', - name='R10', tex_name=r'R_{10}', - model='StaticGen', src='R10', - unit='p.u./h',) + unit='%', no_parse=True,) -class RTEDModel(DCOPFModel): +class RTED(RTEDBase, SFRBase): """ - RTED model. + DC-based real-time economic dispatch (RTED). + RTED extends DCOPF with: + + 1. Param ``pg0``, which can be retrieved from dynamic simulation results. + + 2. RTED has mapping dicts to interface with ANDES. + + 3. RTED routine adds a function ``dc2ac`` to do the AC conversion using ACOPF + + 4. Vars for zonal SFR reserve: ``pru`` and ``prd``; + + 5. Param for linear cost of zonal SFR reserve ``cru`` and ``crd``; + + 6. Param for SFR requirement ``du`` and ``dd``; + + 7. Param for generator ramping: start point ``pg0`` and ramping limit ``R10``; + + The function ``dc2ac`` sets the ``vBus`` value from solved ACOPF. + Without this conversion, dynamic simulation might fail due to the gap between + DC-based dispatch results and AC-based dynamic initialization. + + Notes + ----- + 1. Formulations has been adjusted with interval ``config.t``, 5/60 [Hour] by default. + + 2. The tie-line flow has not been implemented in formulations. """ def __init__(self, system, config): - DCOPFModel.__init__(self, system, config) + RTEDBase.__init__(self, system, config) + SFRBase.__init__(self) self.config.add(OrderedDict((('t', 5/60), ))) @@ -87,28 +141,15 @@ def __init__(self, system, config): self.info = 'Real-time economic dispatch' self.type = 'DCED' - # --- service --- - self.gs = ZonalSum(u=self.zg, zone='Region', - name='gs', tex_name=r'S_{g}', - info='Sum Gen vars vector in shape of zone') - # --- vars --- - self.pru = Var(info='RegUp reserve (system base)', + self.pru = Var(info='RegUp reserve', unit='p.u.', name='pru', tex_name=r'p_{r,u}', model='StaticGen', nonneg=True,) - self.prd = Var(info='RegDn reserve (system base)', + self.prd = Var(info='RegDn reserve', unit='p.u.', name='prd', tex_name=r'p_{r,d}', model='StaticGen', nonneg=True,) + # --- constraints --- - self.ds = ZonalSum(u=self.zb, zone='Region', - name='ds', tex_name=r'S_{d}', - info='Sum pd vector in shape of zone',) - self.pdz = NumOpDual(u=self.ds, u2=self.pl, - fun=np.multiply, - rfun=np.sum, rargs=dict(axis=1), - expand_dims=0, - name='pdz', tex_name=r'p_{d,z}', - unit='p.u.', info='zonal load') self.dud = NumOpDual(u=self.pdz, u2=self.du, fun=np.multiply, rfun=np.reshape, rargs=dict(newshape=(-1,)), name='dud', tex_name=r'd_{u, d}', @@ -119,63 +160,31 @@ def __init__(self, system, config): info='zonal RegDn reserve requirement',) self.rbu = Constraint(name='rbu', type='eq', info='RegUp reserve balance', - e_str='gs @ multiply(ug, pru) - dud',) + e_str='gs @ mul(ug, pru) - dud',) self.rbd = Constraint(name='rbd', type='eq', info='RegDn reserve balance', - e_str='gs @ multiply(ug, prd) - ddd',) + e_str='gs @ mul(ug, prd) - ddd',) self.rru = Constraint(name='rru', type='uq', - info='RegUp reserve ramp', - e_str='multiply(ug, pg + pru) - pmax',) + info='RegUp reserve source', + e_str='mul(ug, pg + pru) - mul(ug, pmax)',) self.rrd = Constraint(name='rrd', type='uq', - info='RegDn reserve ramp', - e_str='multiply(ug, -pg + prd) - pmin',) + info='RegDn reserve source', + e_str='mul(ug, -pg + prd) - mul(ug, pmin)',) self.rgu = Constraint(name='rgu', type='uq', - info='ramp up limit of generator output', - e_str='multiply(ug, pg-pg0-R10)',) + info='Gen ramping up', + e_str='mul(ug, pg-pg0-R10)',) self.rgd = Constraint(name='rgd', type='uq', - info='ramp down limit of generator output', - e_str='multiply(ug, -pg+pg0-R10)',) + info='Gen ramping down', + e_str='mul(ug, -pg+pg0-R10)',) # --- objective --- self.obj.info = 'total generation and reserve cost' - # NOTE: the product of dt and pg is processed using ``dot``, because dt is a numnber - self.obj.e_str = 'sum(c2 @ (t dot pg)**2) ' + \ - '+ sum(c1 @ (t dot pg)) + ug * c0 ' + \ - '+ sum(cru * pru + crd * prd)' - - -class RTED(RTEDData, RTEDModel): - """ - DC-based real-time economic dispatch (RTED). - RTED extends DCOPF with: - - 1. Param ``pg0``, which can be retrieved from dynamic simulation results. - - 2. RTED has mapping dicts to interface with ANDES. - - 3. RTED routine adds a function ``dc2ac`` to do the AC conversion using ACOPF - - 4. Vars for zonal SFR reserve: ``pru`` and ``prd``; - - 5. Param for linear cost of zonal SFR reserve ``cru`` and ``crd``; - - 6. Param for SFR requirement ``du`` and ``dd``; - - 7. Param for generator ramping: start point ``pg0`` and ramping limit ``R10``; - - The function ``dc2ac`` sets the ``vBus`` value from solved ACOPF. - Without this conversion, dynamic simulation might fail due to the gap between - DC-based dispatch results and AC-based dynamic initialization. - - Notes - ----- - 1. Formulations has been adjusted with interval ``config.t``, 5/60 [Hour] by default. - - 2. The tie-line flow has not been implemented in formulations. - """ - - def __init__(self, system, config): - RTEDData.__init__(self) - RTEDModel.__init__(self, system, config) + # NOTE: the product of dt and pg is processed using ``dot``, + # because dt is a numnber + cost = 'sum_squares(mul(c2, pg))' + cost += '+ sum(c1 @ (t dot pg))' + cost += '+ ug * c0' # constant cost + cost += '+ sum(cru * pru + crd * prd)' # reserve cost + self.obj.e_str = cost def dc2ac(self, **kwargs): """ @@ -215,9 +224,44 @@ def dc2ac(self, **kwargs): logger.warning('RTED is converted to AC.') return True - def run(self, **kwargs): + def run(self, no_code=True, **kwargs): """ - Overload ``run()`` method. + Run the routine. + + Parameters + ---------- + no_code : bool, optional + If True, print the generated CVXPY code. Defaults to False. + + Other Parameters + ---------------- + 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. + kwargs : keywords, optional + Additional solver specific arguments. See CVXPY documentation for details. Notes ----- @@ -239,36 +283,39 @@ def __init__(self): self.En = RParam(info='Rated energy capacity', name='En', src='En', tex_name='E_n', unit='MWh', - model='ESD1',) - self.SOCmin = RParam(info='Minimum required value for SOC in limiter', - name='SOCmin', src='SOCmin', - tex_name='SOC_{min}', unit='%', - model='ESD1',) + model='ESD1', no_parse=True,) self.SOCmax = RParam(info='Maximum allowed value for SOC in limiter', name='SOCmax', src='SOCmax', - tex_name='SOC_{max}', unit='%', + tex_name=r'SOC_{max}', unit='%', model='ESD1',) - self.SOCinit = RParam(info='Initial state of charge', + self.SOCmin = RParam(info='Minimum required value for SOC in limiter', + name='SOCmin', src='SOCmin', + tex_name=r'SOC_{min}', unit='%', + model='ESD1',) + self.SOCinit = RParam(info='Initial SOC', name='SOCinit', src='SOCinit', tex_name=r'SOC_{init}', unit='%', model='ESD1',) self.EtaC = RParam(info='Efficiency during charging', name='EtaC', src='EtaC', tex_name=r'\eta_c', unit='%', - model='ESD1',) + model='ESD1', no_parse=True,) self.EtaD = RParam(info='Efficiency during discharging', name='EtaD', src='EtaD', tex_name=r'\eta_d', unit='%', - model='ESD1',) - self.genE = RParam(info='gen of ESD1', - name='genE', tex_name=r'g_{ES}', - model='ESD1', src='gen',) + model='ESD1', no_parse=True,) + self.gene = RParam(info='gen of ESD1', + name='gene', tex_name=r'g_{E}', + model='ESD1', src='gen', + no_parse=True,) + info = 'Ratio of ESD1.pge w.r.t to that of static generator', + self.gammape = RParam(name='gammape', tex_name=r'\gamma_{p,e}', + model='ESD1', src='gammap', + no_parse=True, info=info) # --- service --- self.REtaD = NumOp(name='REtaD', tex_name=r'\frac{1}{\eta_d}', u=self.EtaD, fun=np.reciprocal,) - self.REn = NumOp(name='REn', tex_name=r'\frac{1}{E_n}', - u=self.En, fun=np.reciprocal,) self.Mb = NumOp(info='10 times of max of pmax as big M', name='Mb', tex_name=r'M_{big}', u=self.pmax, fun=np.max, @@ -276,56 +323,78 @@ def __init__(self): array_out=False,) # --- vars --- - self.SOC = Var(info='ESD1 SOC', unit='%', + self.SOC = Var(info='ESD1 State of Charge', unit='%', name='SOC', tex_name=r'SOC', - model='ESD1', pos=True,) - self.ce = VarSelect(u=self.pg, indexer='genE', - name='ce', tex_name=r'C_{ES}', - info='Select pge from pg',) - self.pge = Var(info='ESD1 output power (system base)', - unit='p.u.', name='pge', tex_name=r'p_{g,ES}', - model='ESD1',) - self.ued = Var(info='ESD1 commitment decision', - name='ued', tex_name=r'u_{ES,d}', - model='ESD1', boolean=True,) - self.zue = Var(info='Aux var, :math:`z_{ue} = u_{e,d} * p_{g,ES}`', - name='zue', tex_name=r'z_{ue}', - model='ESD1', pos=True,) - - # --- constraints --- - self.cpge = Constraint(name='cpge', type='eq', - info='Select ESD1 power from StaticGen', - e_str='multiply(ce, pg) - zue',) - + model='ESD1', pos=True, + v0=self.SOCinit,) self.SOClb = Constraint(name='SOClb', type='uq', - info='ESD1 SOC lower bound', + info='SOC lower bound', e_str='-SOC + SOCmin',) self.SOCub = Constraint(name='SOCub', type='uq', - info='ESD1 SOC upper bound', + info='SOC upper bound', e_str='SOC - SOCmax',) + self.ce = VarSelect(u=self.pg, indexer='gene', + name='ce', tex_name=r'C_{E}', + info='Select zue from pg', + gamma='gammape', no_parse=True,) + self.pce = Var(info='ESD1 charging power', + unit='p.u.', name='pce', + tex_name=r'p_{c,E}', + model='ESD1', nonneg=True,) + self.pde = Var(info='ESD1 discharging power', + unit='p.u.', name='pde', + tex_name=r'p_{d,E}', + model='ESD1', nonneg=True,) + self.uce = Var(info='ESD1 charging decision', + name='uce', tex_name=r'u_{c,E}', + model='ESD1', boolean=True,) + self.ude = Var(info='ESD1 discharging decision', + name='ude', tex_name=r'u_{d,E}', + model='ESD1', boolean=True,) + self.zce = Var(name='zce', tex_name=r'z_{c,E}', + model='ESD1', nonneg=True,) + self.zce.info = 'Aux var for charging, :math:`z_{c,e}=u_{c,E}*p_{c,E}`' + self.zde = Var(name='zde', tex_name=r'z_{d,E}', + model='ESD1', nonneg=True,) + self.zde.info = 'Aux var for discharging, :math:`z_{d,e}=u_{d,E}*p_{d,E}`' - self.zclb = Constraint(name='zclb', type='uq', info='zue lower bound', - e_str='- zue + pge',) - self.zcub = Constraint(name='zcub', type='uq', info='zue upper bound', - e_str='zue - pge - Mb dot (1-ued)',) - self.zcub2 = Constraint(name='zcub2', type='uq', info='zue upper bound', - e_str='zue - Mb dot ued',) - - SOCb = 'SOC - SOCinit - t dot REn * EtaC * zue' - SOCb += '- t dot REn * REtaD * (pge - zue)' + # --- constraints --- + self.ceb = Constraint(name='ceb', type='eq', + info='Charging decision bound', + e_str='uce + ude - 1',) + self.cpe = Constraint(name='cpe', type='eq', + info='Select pce from pg', + e_str='ce @ pg - zce - zde',) + + self.zce1 = Constraint(name='zce1', type='uq', info='zce bound 1', + e_str='-zce + pce',) + self.zce2 = Constraint(name='zce2', type='uq', info='zce bound 2', + e_str='zce - pce - Mb dot (1-uce)',) + self.zce3 = Constraint(name='zce3', type='uq', info='zce bound 3', + e_str='zce - Mb dot uce',) + + self.zde1 = Constraint(name='zde1', type='uq', info='zde bound 1', + e_str='-zde + pde',) + self.zde2 = Constraint(name='zde2', type='uq', info='zde bound 2', + e_str='zde - pde - Mb dot (1-ude)',) + self.zde3 = Constraint(name='zde3', type='uq', info='zde bound 3', + e_str='zde - Mb dot ude',) + + SOCb = 'mul(En, (SOC - SOCinit)) - t dot mul(EtaC, zce)' + SOCb += '+ t dot mul(REtaD, zde)' self.SOCb = Constraint(name='SOCb', type='eq', - info='ESD1 SOC balance', e_str=SOCb,) + info='ESD1 SOC balance', + e_str=SOCb,) -class RTED2(RTEDData, RTEDModel, ESD1Base): +class RTEDES(RTED, ESD1Base): """ RTED with energy storage :ref:`ESD1`. The bilinear term in the formulation is linearized with big-M method. """ def __init__(self, system, config): - RTEDData.__init__(self) - RTEDModel.__init__(self, system, config) + RTED.__init__(self, system, config) ESD1Base.__init__(self) self.info = 'Real-time economic dispatch with energy storage' self.type = 'DCED' diff --git a/ams/routines/uc.py b/ams/routines/uc.py index 1188944a..0b2d4435 100644 --- a/ams/routines/uc.py +++ b/ams/routines/uc.py @@ -1,29 +1,105 @@ """ -Real-time economic dispatch. +Unit commitment routines. """ -import logging # NOQA -from collections import OrderedDict # NOQA -import numpy as np # NOQA -import pandas as pd # NOQA +import logging +from collections import OrderedDict +import numpy as np +import pandas as pd -from ams.core.param import RParam # NOQA -from ams.core.service import (NumOp, NumHstack, - NumOpDual, MinDur, ZonalSum) # NOQA -from ams.routines.ed import EDData, EDModel # NOQA -from ams.routines.rted import ESD1Base # NOQA +from ams.core.param import RParam +from ams.core.service import (NumOp, NumOpDual, MinDur) +from ams.routines.rted import RTEDBase +from ams.routines.ed import SRBase, MPBase, ESD1MPBase -from ams.opt.omodel import Var, Constraint # NOQA +from ams.opt.omodel import Var, Constraint logger = logging.getLogger(__name__) -class UCData(EDData): +class NSRBase: """ - UC data. + Base class for non-spinning reserve. """ - def __init__(self): - EDData.__init__(self) + def __init__(self) -> None: + self.cnsr = RParam(info='cost for non-spinning reserve', + name='cnsr', tex_name=r'c_{nsr}', + model='NSRCost', src='cnsr', + unit=r'$/(p.u.*h)', + indexer='gen', imodel='StaticGen',) + self.dnsrp = RParam(info='non-spinning reserve requirement in percentage', + name='dnsr', tex_name=r'd_{nsr}', + model='NSR', src='demand', + unit='%',) + self.prns = Var(info='non-spinning reserve', + name='prns', tex_name=r'p_{r, ns}', + model='StaticGen', nonneg=True,) + + self.dnsrpz = NumOpDual(u=self.pdz, u2=self.dnsrp, fun=np.multiply, + name='dnsrpz', tex_name=r'd_{nsr, p, z}', + info='zonal non-spinning reserve requirement in percentage',) + self.dnsr = NumOpDual(u=self.dnsrpz, u2=self.sd, fun=np.multiply, + rfun=np.transpose, + name='dnsr', tex_name=r'd_{nsr}', + info='zonal non-spinning reserve requirement', + no_parse=True,) + + # NOTE: define e_str in dispatch model + self.prnsb = Constraint(info='non-spinning reserve balance', + name='prnsb', type='eq',) + self.rnsr = Constraint(info='non-spinning reserve requirement', + name='rnsr', type='uq',) + + +class UC(RTEDBase, MPBase, SRBase, NSRBase): + """ + DC-based unit commitment (UC). + The bilinear term in the formulation is linearized with big-M method. + + Penalty for unserved load is introduced as ``config.cul`` (:math:`c_{ul, cfg}`), + 1000 [$/p.u.] by default. + + Constraints include power balance, ramping, spinning reserve, non-spinning reserve, + minimum ON/OFF duration. + The cost inludes generation cost, startup cost, shutdown cost, spinning reserve cost, + non-spinning reserve cost, and unserved energy penalty. + + Method ``_initial_guess`` is used to make initial guess for commitment decision if all + generators are online at initial. It is a simple heuristic method, which may not be optimal. + + Notes + ----- + 1. Formulations has been adjusted with interval ``config.t`` + + 3. The tie-line flow has not been implemented in formulations. + + References + ---------- + 1. Huang, Y., Pardalos, P. M., & Zheng, Q. P. (2017). Electrical power unit commitment: deterministic and + two-stage stochastic programming models and algorithms. Springer. + + 2. D. A. Tejada-Arango, S. Lumbreras, P. Sánchez-Martín and A. Ramos, "Which Unit-Commitment Formulation + is Best? A Comparison Framework," in IEEE Transactions on Power Systems, vol. 35, no. 4, pp. 2926-2936, + July 2020, doi: 10.1109/TPWRS.2019.2962024. + """ + + def __init__(self, system, config): + RTEDBase.__init__(self, system, config) + MPBase.__init__(self) + SRBase.__init__(self) + NSRBase.__init__(self) + + self.config.add(OrderedDict((('t', 1), + ('cul', 10000), + ))) + self.config.add_extra("_help", + t="time interval in hours", + cul="penalty for unserved load, $/p.u.", + ) + + self.info = 'unit commitment' + self.type = 'DCUC' + self.timeslot.model = 'UCTSlot' self.csu = RParam(info='startup cost', name='csu', tex_name=r'c_{su}', @@ -48,33 +124,32 @@ def __init__(self): self.timeslot.info = 'Time slot for multi-period UC' self.timeslot.model = 'UCTSlot' - self.cnsr = RParam(info='cost for spinning reserve', - name='cnsr', tex_name=r'c_{nsr}', - model='NSRCost', src='cnsr', - unit=r'$/(p.u.*h)', - indexer='gen', imodel='StaticGen',) - self.dnsrp = RParam(info='non-spinning reserve requirement in percentage', - name='dnsr', tex_name=r'd_{nsr}', - model='NSR', src='demand', - unit='%',) + self.ug.expand_dims = 1 + # NOTE: extend pg to 2D matrix, where row is gen and col is timeslot + self.pg.horizon = self.timeslot + self.pg.info = '2D Gen power' -class UCModel(EDModel): - """ - UC model. - """ + self.plf.horizon = self.timeslot + self.plf.info = '2D Line flow' - def __init__(self, system, config): - EDModel.__init__(self, system, config) + self.prs.horizon = self.timeslot + self.prs.info = '2D Spinning reserve' - self.config.add(OrderedDict((('cul', 1000), - ))) - self.config.add_extra("_help", - cul="penalty for unserved load, $/p.u.", - ) + self.prns.horizon = self.timeslot + self.prns.info = '2D Non-spinning reserve' - self.info = 'unit commitment' - self.type = 'DCUC' + # TODO: havn't test non-controllability? + self.ctrle.u2 = self.tlv + 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 # --- vars --- self.ugd = Var(info='commitment decision', @@ -104,17 +179,30 @@ def __init__(self, system, config): e_str='ugd @ Mr - vgd[:, 1:]',) self.actv0 = Constraint(name='actv0', type='eq', info='initial startup action', - e_str='ugd[:, 0] - ug - vgd[:, 0]',) + e_str='ugd[:, 0] - ug[:, 0] - vgd[:, 0]',) self.actw = Constraint(name='actw', type='eq', info='shutdown action', e_str='-ugd @ Mr - wgd[:, 1:]',) self.actw0 = Constraint(name='actw0', type='eq', info='initial shutdown action', - e_str='-ugd[:, 0] + ug - wgd[:, 0]',) + e_str='-ugd[:, 0] + ug[:, 0] - wgd[:, 0]',) # --- constraints --- - self.pb.e_str = '- gs @ zug + pds' # power balance - self.pb.type = 'uq' + self.pb.e_str = 'gs @ zug - pdsz' # power balance + self.pb.type = 'uq' # soft constraint + + # spinning reserve + self.prsb.e_str = 'mul(ugd, mul(pmax, tlv)) - zug - prs' + # spinning reserve requirement + self.rsr.e_str = '-gs@prs + dsr' + + # non-spinning reserve + self.prnsb.e_str = 'mul(1-ugd, mul(pmax, tlv)) - prns' + # non-spinning reserve requirement + self.rnsr.e_str = '-gs@prns + dnsr' + + # --- bus power injection --- + self.pnb.e_str = 'PTDF@(Cgi@pg - Cli@pds) - plf' # --- big M for ugd*pg --- self.Mzug = NumOp(info='10 times of max of pmax as big M for zug', @@ -129,26 +217,6 @@ def __init__(self, system, config): self.zugub2 = Constraint(name='zugub2', info='zug upper bound', type='uq', e_str='zug - Mzug dot ugd') - # --- reserve --- - # 1) non-spinning reserve - self.dnsrpz = NumOpDual(u=self.pdz, u2=self.dnsrp, fun=np.multiply, - name='dnsrpz', tex_name=r'd_{nsr, p, z}', - info='zonal non-spinning reserve requirement in percentage',) - self.dnsr = NumOpDual(u=self.dnsrpz, u2=self.sd, fun=np.multiply, - rfun=np.transpose, - name='dnsr', tex_name=r'd_{nsr}', - info='zonal non-spinning reserve requirement',) - self.nsr = Constraint(name='nsr', info='non-spinning reserve', type='uq', - e_str='-gs@(multiply((1 - ugd), Rpmax)) + dnsr') - # 2) spinning reserve - self.dsrpz = NumOpDual(u=self.pdz, u2=self.dsrp, fun=np.multiply, - name='dsrpz', tex_name=r'd_{sr, p, z}', - info='zonal spinning reserve requirement in percentage',) - self.dsr = NumOpDual(u=self.dsrpz, u2=self.sd, fun=np.multiply, - rfun=np.transpose, - name='dsr', tex_name=r'd_{sr}', - info='zonal spinning reserve requirement',) - # --- minimum ON/OFF duration --- self.Con = MinDur(u=self.pg, u2=self.td1, name='Con', tex_name=r'T_{on}', @@ -163,17 +231,13 @@ def __init__(self, system, config): name='doff', type='uq', e_str='multiply(Coff, wgd) - (1 - ugd)') - # --- penalty for unserved load --- - self.Cgi = NumOp(u=self.Cg, fun=np.linalg.pinv, - name='Cgi', tex_name=r'C_{g}^{-1}', - info='inverse of Cg',) - # --- objective --- - gcost = 'sum(c2 @ (t dot zug)**2 + c1 @ (t dot zug) + c0 * ugd)' - acost = ' + sum(csu * vgd + csd * wgd)' - srcost = ' + sum(csr @ (multiply(Rpmax, ugd) - zug))' - nsrcost = ' + sum(cnsr @ multiply((1 - ugd), Rpmax))' - dcost = ' + sum(cul dot pos(gs @ pg - pds))' + gcost = 'sum(c2 @ (t dot zug)**2 + c1 @ (t dot zug))' + gcost += '+ sum(mul(mul(ug, c0), tlv))' + acost = ' + sum(csu * vgd + csd * wgd)' # action + srcost = ' + sum(csr @ prs)' # spinning reserve + nsrcost = ' + sum(cnsr @ prns)' # non-spinning reserve + dcost = ' + sum(cul dot pos(pdsz - gs @ pg))' # unserved load self.obj.e_str = gcost + acost + srcost + nsrcost + dcost def _initial_guess(self): @@ -212,65 +276,47 @@ def _initial_guess(self): g_idx = priority[0] ug0 = 0 self.system.StaticGen.set(src='u', attr='v', idx=g_idx, value=ug0) - logger.warning(f'Turn off StaticGen {g_idx} as initial guess for commitment.') + logger.warning(f'Turn off StaticGen {g_idx} as initial commitment guess.') return True def init(self, **kwargs): self._initial_guess() - super().init(**kwargs) - - -class UC(UCData, UCModel): - """ - DC-based unit commitment (UC). - The bilinear term in the formulation is linearized with big-M method. + return super().init(**kwargs) - Penalty for unserved load is introduced as ``config.cul`` (:math:`c_{ul, cfg}`), - 1000 [$/p.u.] by default. - - Constraints include power balance, ramping, spinning reserve, non-spinning reserve, - minimum ON/OFF duration. - The cost inludes generation cost, startup cost, shutdown cost, spinning reserve cost, - non-spinning reserve cost, and unserved energy penalty. - - Method ``_initial_guess`` is used to make initial guess for commitment decision if all - generators are online at initial. It is a simple heuristic method, which may not be optimal. - - Notes - ----- - 1. Formulations has been adjusted with interval ``config.t`` - - 3. The tie-line flow has not been implemented in formulations. - - References - ---------- - 1. Huang, Y., Pardalos, P. M., & Zheng, Q. P. (2017). Electrical power unit commitment: deterministic and - two-stage stochastic programming models and algorithms. Springer. + def dc2ac(self, **kwargs): + """ + AC conversion ``dc2ac`` is not implemented yet for + multi-period dispatch. + """ + return NotImplementedError - 2. D. A. Tejada-Arango, S. Lumbreras, P. Sánchez-Martín and A. Ramos, "Which Unit-Commitment Formulation - is Best? A Comparison Framework," in IEEE Transactions on Power Systems, vol. 35, no. 4, pp. 2926-2936, - July 2020, doi: 10.1109/TPWRS.2019.2962024. - """ + def unpack(self, **kwargs): + """ + Multi-period dispatch will not unpack results from + solver into devices. - def __init__(self, system, config): - UCData.__init__(self) - UCModel.__init__(self, system, config) + # TODO: unpack first period results, and allow input + # to specify which period to unpack. + """ + return None -class UC2(UCData, UCModel, ESD1Base): +class UCES(UC, ESD1MPBase): """ UC with energy storage :ref:`ESD1`. """ def __init__(self, system, config): - UCData.__init__(self) - UCModel.__init__(self, system, config) - ESD1Base.__init__(self) + UC.__init__(self, system, config) + ESD1MPBase.__init__(self) self.info = 'unit commitment with energy storage' self.type = 'DCUC' self.SOC.horizon = self.timeslot - self.pge.horizon = self.timeslot - self.ued.horizon = self.timeslot - self.zue.horizon = self.timeslot + self.pce.horizon = self.timeslot + self.pde.horizon = self.timeslot + self.uce.horizon = self.timeslot + self.ude.horizon = self.timeslot + self.zce.horizon = self.timeslot + self.zde.horizon = self.timeslot diff --git a/ams/shared.py b/ams/shared.py new file mode 100644 index 00000000..0c5a3090 --- /dev/null +++ b/ams/shared.py @@ -0,0 +1,53 @@ +""" +Shared constants and delayed imports. + +This module is supplementary to the ``andes.shared`` module. +""" +import logging +from functools import wraps + +import cvxpy as cp + +from andes.utils.lazyimport import LazyImport + + +logger = logging.getLogger(__name__) + + +igraph = LazyImport("import igraph") + +# NOTE: copied from CVXPY documentation +MIP_SOLVERS = ['CBC', 'COPT', 'GLPK_MI', 'CPLEX', 'GUROBI', + 'MOSEK', 'SCIP', 'XPRESS', 'SCIPY'] + +INSTALLED_SOLVERS = cp.installed_solvers() + + +def require_MIP_solver(f): + """ + Decorator for functions that require MIP solver. + """ + + @wraps(f) + def wrapper(*args, **kwargs): + if not any(s in MIP_SOLVERS for s in INSTALLED_SOLVERS): + raise ModuleNotFoundError("No MIP solver is available.") + return f(*args, **kwargs) + + return wrapper + + +def require_igraph(f): + """ + Decorator for functions that require igraph. + """ + + @wraps(f) + def wrapper(*args, **kwargs): + try: + getattr(igraph, "__version__") + except AttributeError: + logger.error("Package `igraph` is not installed.") + return f(*args, **kwargs) + + return wrapper diff --git a/dev/demo/8.ESD1.ipynb b/dev/demo/8.ESD1.ipynb deleted file mode 100644 index 65be00d3..00000000 --- a/dev/demo/8.ESD1.ipynb +++ /dev/null @@ -1,248 +0,0 @@ -{ - "cells": [ - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# ESD" - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [], - "source": [ - "import numpy as np\n", - "\n", - "import andes\n", - "import ams\n", - "\n", - "import pandas as pd\n", - "\n", - "import json\n", - "\n", - "import cvxpy as cp\n", - "\n", - "import re" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "'0.6.7.post44.dev0+gc6b8426'" - ] - }, - "execution_count": 2, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "ams.__version__" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": {}, - "outputs": [], - "source": [ - "ams.config_logger(stream_level=10)" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "Input format guessed as xlsx.\n", - "Parsing input file \"/Users/jinningwang/Documents/work/ams/ams/cases/ieee39/ieee39_esd1.xlsx\"...\n", - "Input file parsed in 0.1659 seconds.\n", - "Adjust bus index to start from 0.\n", - "System set up in 0.0040 seconds.\n" - ] - } - ], - "source": [ - "sp = ams.load(ams.get_case('ieee39/ieee39_esd1.xlsx'),\n", - " setup=True,\n", - " )" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "Setup model of RTED\n", - "RTED data check passed.\n", - "- Generating symbols for RTED\n", - "Set constrs pb: sum(pd) - sum(pg) == 0\n", - "Set constrs pinj: Cg@(pn - pd) - pg == 0\n", - "Set constrs lub: PTDF @ (pn - pd) - rate_a <= 0\n", - "Set constrs llb: - PTDF @ (pn - pd) - rate_a <= 0\n", - "Set constrs rbu: gs @ pru - du == 0\n", - "Set constrs rbd: gs @ prd - dd == 0\n", - "Set constrs rru: pg + pru - pmax <= 0\n", - "Set constrs rrd: -pg + prd - pmin <= 0\n", - "Set constrs rgu: pg - pg0 - R10h <= 0\n", - "Set constrs rgd: -pg + pg0 - R10h <= 0\n", - "RTED model set up in 0.0168 seconds.\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Set parameter Username\n", - "Academic license - for non-commercial use only - expires 2024-05-21\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "RTED solved as optimal in 0.0218 seconds with exit code 0.\n" - ] - }, - { - "data": { - "text/plain": [ - "True" - ] - }, - "execution_count": 5, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "sp.RTED.run(solver=cp.GUROBI, reoptimize=True)" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "Setup model of RTED2\n", - "RTED2 data check passed.\n", - "- Generating symbols for RTED2\n", - "Set constrs pb: sum(pd) - sum(pg) == 0\n", - "Set constrs pinj: Cg@(pn - pd) - pg == 0\n", - "Set constrs lub: PTDF @ (pn - pd) - rate_a <= 0\n", - "Set constrs llb: - PTDF @ (pn - pd) - rate_a <= 0\n", - "Set constrs rbu: gs @ pru - du == 0\n", - "Set constrs rbd: gs @ prd - dd == 0\n", - "Set constrs rru: pg + pru - pmax <= 0\n", - "Set constrs rrd: -pg + prd - pmin <= 0\n", - "Set constrs rgu: pg - pg0 - R10h <= 0\n", - "Set constrs rgd: -pg + pg0 - R10h <= 0\n", - "Set constrs pges: e1s@pg - zc == 0\n", - "Set constrs SOClb: -SOC + SOCmin <= 0\n", - "Set constrs SOCub: SOC - SOCmax <= 0\n", - "Set constrs zclb: - zc + pec <= 0\n", - "Set constrs zcub: zc - pec - Mb@(1-uc) <= 0\n", - "Set constrs zcub2: zc - Mb@uc <= 0\n", - "Set constrs SOCb: SOC - SOCinit - power(dth, 1)*REn*EtaC*zc - power(dth, 1)*REn*REtaD*(pec - zc) == 0\n", - "RTED2 model set up in 0.0144 seconds.\n", - "RTED2 solved as optimal in 0.0198 seconds with exit code 0.\n" - ] - }, - { - "data": { - "text/plain": [ - "True" - ] - }, - "execution_count": 6, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "sp.RTED2.run(solver=cp.GUROBI, reoptimize=True)" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Routine: RTED\n", - "Objective: 629.8546048020892\n", - "------------------\n", - "Routine: RTED2\n", - "Objective: 629.8546048240468\n", - "ESD1 SOC level: [0.20548643]\n", - "Charging: uc=[1.]; pec=[6.58371675], zc=[6.58371675]\n" - ] - } - ], - "source": [ - "rtn = sp.RTED\n", - "\n", - "print(\"Routine: RTED\")\n", - "print(f\"Objective: {rtn.obj.v}\")\n", - "print(\"------------------\")\n", - "\n", - "rtn = sp.RTED2\n", - "\n", - "print(\"Routine: RTED2\")\n", - "print(f\"Objective: {rtn.obj.v}\")\n", - "print(f\"ESD1 SOC level: {rtn.SOC.v}\")\n", - "print(f\"Charging: uc={rtn.uc.v}; pec={rtn.pec.v}, zc={rtn.zc.v}\")" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "ams", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.9.16" - }, - "orig_nbformat": 4, - "vscode": { - "interpreter": { - "hash": "d2b3bf80176349caa68dc4a3c77bd06eaade8abc678330f7d1c813c53380e5d2" - } - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/docs/source/api.rst b/docs/source/api.rst index c3d483da..9fabf3cd 100644 --- a/docs/source/api.rst +++ b/docs/source/api.rst @@ -43,7 +43,7 @@ Routines :caption: Routines :template: autosummary/module_toctree.rst - ams.routines + ams.routines.routine Optimization diff --git a/docs/source/conf.py b/docs/source/conf.py index 770c56d5..7027c190 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -9,9 +9,8 @@ # -- General configuration --------------------------------------------------- # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration -# TODO: fix the importing error later on import ams -# import shutil +import shutil extensions = [ 'sphinx.ext.autodoc', @@ -184,11 +183,12 @@ exec(open("genmodelref.py").read()) exec(open("genroutineref.py").read()) -jupyter_execute_notebooks = "off" - # import and execute model reference generation script -# TODO: use this for routines doumentation later on -# exec(open("genroutineref.py").read()) +shutil.rmtree("_examples", ignore_errors=True) +shutil.copytree("../../examples", "_examples", ) +shutil.rmtree("_examples/demonstration") + +jupyter_execute_notebooks = "off" # sphinx-panels shouldn't add bootstrap css since the pydata-sphinx-theme # already loads it diff --git a/docs/source/examples/index.rst b/docs/source/examples/index.rst index 0ebbce0a..36303b9d 100644 --- a/docs/source/examples/index.rst +++ b/docs/source/examples/index.rst @@ -5,4 +5,16 @@ Examples .. _`development demos`: https://github.com/jinningwang/ams/tree/master/dev/demo -Refer to the development `development demos`_ for examples prior to preparing this section. \ No newline at end of file +Refer to the development `development demos`_ for examples prior to preparing this section. + +A collection of examples are presented to supplement the tutorial. The +examples below are identical to the Jupyter Notebook in the ``examples`` +folder of the repository +`here `__. + +.. toctree:: + :maxdepth: 2 + :caption: Scripting + + ../_examples/ex1.ipynb + ../_examples/ex2.ipynb diff --git a/docs/source/release-notes.rst b/docs/source/release-notes.rst index c6456808..ddf55835 100644 --- a/docs/source/release-notes.rst +++ b/docs/source/release-notes.rst @@ -9,6 +9,14 @@ The APIs before v3.0.0 are in beta and may change without prior notice. Pre-v1.0.0 ========== +v0.7.4 (2023-11-29) +------------------- + +- Refactor routins and optimization models to improve performance +- Fix routines modeling +- Add examples +- Fix built-in cases + v0.7.3 (2023-11-3) ------------------- diff --git a/examples/.placeholder b/examples/.placeholder deleted file mode 100644 index e69de29b..00000000 diff --git a/examples/demonstration/1. ESD1.ipynb b/examples/demonstration/1. ESD1.ipynb new file mode 100644 index 00000000..7e86af06 --- /dev/null +++ b/examples/demonstration/1. ESD1.ipynb @@ -0,0 +1,657 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Dispatch with Energy Storage" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this case, we will show the usage of energy storage included dispatch.\n", + "\n", + "In AMS, ``ESD1`` is an dispatch model for energy storage, which has a corresponding\n", + "dynamic model ``ESD1`` in ANDES." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import pandas as pd\n", + "\n", + "import ams\n", + "\n", + "import datetime" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Last run time: 2023-11-28 22:03:14\n", + "ams:0.7.3.post74.dev0+g47c1ae3\n" + ] + } + ], + "source": [ + "print(\"Last run time:\", datetime.datetime.now().strftime(\"%Y-%m-%d %H:%M:%S\"))\n", + "\n", + "print(f'ams:{ams.__version__}')" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "ams.config_logger(stream_level=20)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "A small-size PJM 5-bus case with ESD1 is used in this example." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Parsing input file \"/Users/jinningwang/Documents/work/ams/ams/cases/5bus/pjm5bus_uced_esd1.xlsx\"...\n", + "Input file parsed in 0.1205 seconds.\n", + "System set up in 0.0031 seconds.\n" + ] + } + ], + "source": [ + "sp = ams.load(ams.get_case('5bus/pjm5bus_uced_esd1.xlsx'),\n", + " setup=True,)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The model information can be inspected as follow." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
idxunamebusgenSngammapgammaqSOCminSOCmaxSOCinitEnEtaCEtaD
uid
0ESD1_11.0ESD1_10PV_1100.01.01.00.01.00.2100.01.01.0
\n", + "
" + ], + "text/plain": [ + " idx u name bus gen Sn gammap gammaq SOCmin SOCmax \\\n", + "uid \n", + "0 ESD1_1 1.0 ESD1_1 0 PV_1 100.0 1.0 1.0 0.0 1.0 \n", + "\n", + " SOCinit En EtaC EtaD \n", + "uid \n", + "0 0.2 100.0 1.0 1.0 " + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sp.ESD1.as_df()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "`RTEDES` is the RTED model extended to include energy storage.\n", + "\n", + "Note that mixed integer linear programming (MILP) requires\n", + "capable solvers such as Gurobi or CPLEX.\n", + "They might require license and installation.\n", + "\n", + "More details can be found at [CVXPY - Choosing a solver](https://www.cvxpy.org/tutorial/advanced/index.html#choosing-a-solver)." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Routine initialized in 0.0140 seconds.\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Set parameter Username\n", + "Academic license - for non-commercial use only - expires 2024-05-21\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "RTEDES solved as optimal in 0.0341 seconds, converged after 0 iteration using solver GUROBI.\n" + ] + }, + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sp.RTEDES.run(solver='GUROBI', reoptimize=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that, in RTED, the time interval is 5/60 [H] by default, and the\n", + "dispatch model has been adjusted accordingly." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
VarValueinfo
0uce[0.0]ESD1 charging decision
1ude[1.0]ESD1 discharging decision
2pce[0.0]ESD1 charging power
3pde[2.1]ESD1 discharging power
4SOC[0.1982]ESD1 State of Charge
5SOCinit[0.2]Initial SOC
\n", + "
" + ], + "text/plain": [ + " Var Value info\n", + "0 uce [0.0] ESD1 charging decision\n", + "1 ude [1.0] ESD1 discharging decision\n", + "2 pce [0.0] ESD1 charging power\n", + "3 pde [2.1] ESD1 discharging power\n", + "4 SOC [0.1982] ESD1 State of Charge\n", + "5 SOCinit [0.2] Initial SOC" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "RTEDESres = pd.DataFrame()\n", + "\n", + "items = [sp.RTEDES.uce, sp.RTEDES.ude,\n", + " sp.RTEDES.pce, sp.RTEDES.pde,\n", + " sp.RTEDES.SOC, sp.RTEDES.SOCinit]\n", + "\n", + "RTEDESres['Var'] = [item.name for item in items]\n", + "RTEDESres['Value'] = [item.v.round(4) for item in items]\n", + "RTEDESres['info'] = [item.info for item in items]\n", + "\n", + "RTEDESres" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Similarly, multi-period dispatch ``EDES`` and ``UCES`` are also available.\n", + "They have 1 [H] time interval by default." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Routine initialized in 0.0238 seconds.\n", + "EDES solved as optimal in 0.0349 seconds, converged after 0 iteration using solver GUROBI.\n" + ] + }, + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sp.EDES.run(solver='GUROBI', reoptimize=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
VarValueinfo
0uce[[0.0, 0.0, 0.0, 0.0]]ESD1 charging decision
1ude[[1.0, 1.0, 1.0, 1.0]]ESD1 discharging decision
2pce[[0.0, 0.0, 0.0, 0.0]]ESD1 charging power
3pde[[2.1, 2.1, 2.1, 2.1]]ESD1 discharging power
4SOC[[0.179, 0.179, 0.179, 0.179]]ESD1 State of Charge
5SOCinit[0.2]Initial SOC
\n", + "
" + ], + "text/plain": [ + " Var Value info\n", + "0 uce [[0.0, 0.0, 0.0, 0.0]] ESD1 charging decision\n", + "1 ude [[1.0, 1.0, 1.0, 1.0]] ESD1 discharging decision\n", + "2 pce [[0.0, 0.0, 0.0, 0.0]] ESD1 charging power\n", + "3 pde [[2.1, 2.1, 2.1, 2.1]] ESD1 discharging power\n", + "4 SOC [[0.179, 0.179, 0.179, 0.179]] ESD1 State of Charge\n", + "5 SOCinit [0.2] Initial SOC" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ESESres = pd.DataFrame()\n", + "\n", + "items = [sp.EDES.uce, sp.EDES.ude,\n", + " sp.EDES.pce, sp.EDES.pde,\n", + " sp.EDES.SOC, sp.EDES.SOCinit]\n", + "\n", + "ESESres['Var'] = [item.name for item in items]\n", + "ESESres['Value'] = [item.v.round(4) for item in items]\n", + "ESESres['info'] = [item.info for item in items]\n", + "\n", + "ESESres" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "All generators are online at initial, make initial guess for commitment.\n", + "Turn off StaticGen PV_1 as initial commitment guess.\n", + "Routine initialized in 0.0272 seconds.\n", + "UCES solved as optimal in 0.0369 seconds, converged after 0 iteration using solver GUROBI.\n" + ] + }, + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sp.UCES.run(solver='GUROBI', reoptimize=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
VarValueinfo
0uce[[0.0, 0.0, 0.0, 0.0]]ESD1 charging decision
1ude[[1.0, 1.0, 1.0, 1.0]]ESD1 discharging decision
2pce[[0.0, 0.0, 0.0, 0.0]]ESD1 charging power
3pde[[0.0, 0.0, 0.0, 0.0]]ESD1 discharging power
4SOC[[0.2, 0.2, 0.2, 0.2]]ESD1 State of Charge
5SOCinit[0.2]Initial SOC
\n", + "
" + ], + "text/plain": [ + " Var Value info\n", + "0 uce [[0.0, 0.0, 0.0, 0.0]] ESD1 charging decision\n", + "1 ude [[1.0, 1.0, 1.0, 1.0]] ESD1 discharging decision\n", + "2 pce [[0.0, 0.0, 0.0, 0.0]] ESD1 charging power\n", + "3 pde [[0.0, 0.0, 0.0, 0.0]] ESD1 discharging power\n", + "4 SOC [[0.2, 0.2, 0.2, 0.2]] ESD1 State of Charge\n", + "5 SOCinit [0.2] Initial SOC" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "UCESres = pd.DataFrame()\n", + "\n", + "items = [sp.UCES.uce, sp.UCES.ude,\n", + " sp.UCES.pce, sp.UCES.pde,\n", + " sp.UCES.SOC, sp.UCES.SOCinit]\n", + "\n", + "UCESres['Var'] = [item.name for item in items]\n", + "UCESres['Value'] = [item.v.round(4) for item in items]\n", + "UCESres['info'] = [item.info for item in items]\n", + "\n", + "UCESres" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "ams", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.18" + }, + "orig_nbformat": 4, + "vscode": { + "interpreter": { + "hash": "d2b3bf80176349caa68dc4a3c77bd06eaade8abc678330f7d1c813c53380e5d2" + } + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/ex1.ipynb b/examples/ex1.ipynb index d708d84f..c72d136a 100644 --- a/examples/ex1.ipynb +++ b/examples/ex1.ipynb @@ -9,12 +9,10 @@ ] }, { - "cell_type": "code", - "execution_count": 1, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "import ams" + "This example gives a \"hello world\" example to use AMS." ] }, { @@ -35,31 +33,33 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ - "import ams" + "import ams\n", + "\n", + "import datetime" ] }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 2, "metadata": {}, "outputs": [ { - "data": { - "text/plain": [ - "'0.6.5.post25.dev0+ga5fc571'" - ] - }, - "execution_count": 3, - "metadata": {}, - "output_type": "execute_result" + "name": "stdout", + "output_type": "stream", + "text": [ + "Last run time: 2023-11-28 22:06:12\n", + "ams:0.7.3.post74.dev0+g47c1ae3\n" + ] } ], "source": [ - "ams.__version__" + "print(\"Last run time:\", datetime.datetime.now().strftime(\"%Y-%m-%d %H:%M:%S\"))\n", + "\n", + "print(f'ams:{ams.__version__}')" ] }, { @@ -74,7 +74,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ @@ -114,26 +114,28 @@ "source": [ "AMS support multiple input file formats, including AMS ``.xlsx`` file, MATPOWER ``.m`` file, PYPOWER ``.py`` file, and PSS/E ``.raw`` file.\n", "\n", - "Here we use the AMS ``.xlsx`` file as an example. The source file locates at ``$HOME/ams/ams/cases/ieee14/ieee14_opf.xlsx``." + "Here we use the AMS ``.xlsx`` file as an example. The source file locates at ``$HOME/ams/ams/cases/ieee39/ieee39_uced.xlsx``." ] }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ - "Parsing input file \"/Users/jinningwang/Documents/work/ams/ams/cases/ieee14/ieee14_opf.xlsx\"...\n", - "Input file parsed in 0.4313 seconds.\n", - "System set up in 0.0021 seconds.\n" + "Parsing input file \"/Users/jinningwang/Documents/work/ams/ams/cases/5bus/pjm5bus_uced.xlsx\"...\n", + "Input file parsed in 0.0978 seconds.\n", + "System set up in 0.0032 seconds.\n" ] } ], "source": [ - "sp = ams.load(ams.get_case('ieee14/ieee14_opf.xlsx'), default_config=True)" + "sp = ams.load(ams.get_case('5bus/pjm5bus_uced.xlsx'),\n", + " setup=True,\n", + " no_output=True,)" ] }, { @@ -154,28 +156,37 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "OrderedDict([('Summary', Summary (1 device) at 0x1590a8fd0),\n", - " ('Bus', Bus (14 devices) at 0x15996c040),\n", - " ('PQ', PQ (11 devices) at 0x15996cac0),\n", - " ('PV', PV (4 devices) at 0x15996ce80),\n", - " ('Slack', Slack (1 device) at 0x1599799d0),\n", - " ('Shunt', Shunt (2 devices) at 0x159984400),\n", - " ('Line', Line (20 devices) at 0x1599848b0),\n", - " ('ESD1', ESD1 (0 devices) at 0x15998afa0),\n", - " ('Area', Area (2 devices) at 0x15999d700),\n", - " ('Region', Region (0 devices) at 0x15999dfa0),\n", - " ('SFR', SFR (0 devices) at 0x159a54790),\n", - " ('GCost', GCost (5 devices) at 0x159a54e80),\n", - " ('SFRCost', SFRCost (0 devices) at 0x159a61520)])" + "OrderedDict([('Summary', Summary (3 devices) at 0x153c7bd00),\n", + " ('Bus', Bus (5 devices) at 0x10774b160),\n", + " ('PQ', PQ (3 devices) at 0x153c8e640),\n", + " ('PV', PV (3 devices) at 0x153c8eb80),\n", + " ('Slack', Slack (1 device) at 0x153c9e610),\n", + " ('Shunt', Shunt (0 devices) at 0x153caa0d0),\n", + " ('Line', Line (7 devices) at 0x153caa580),\n", + " ('ESD1', ESD1 (0 devices) at 0x153cb8c70),\n", + " ('REGCV1', REGCV1 (0 devices) at 0x153cc63d0),\n", + " ('Area', Area (3 devices) at 0x153cc6b20),\n", + " ('Region', Region (2 devices) at 0x153cd32e0),\n", + " ('SFR', SFR (2 devices) at 0x153cd3a90),\n", + " ('SR', SR (2 devices) at 0x153ce4130),\n", + " ('NSR', NSR (2 devices) at 0x153ce4550),\n", + " ('GCost', GCost (4 devices) at 0x153ce4970),\n", + " ('SFRCost', SFRCost (4 devices) at 0x153cf1040),\n", + " ('SRCost', SRCost (4 devices) at 0x153cf15e0),\n", + " ('NSRCost', NSRCost (4 devices) at 0x153cf1a00),\n", + " ('REGCV1Cost', REGCV1Cost (0 devices) at 0x153cf1e20),\n", + " ('TimeSlot', TimeSlot (0 devices) at 0x153cfb160),\n", + " ('EDTSlot', EDTSlot (24 devices) at 0x153cfbd30),\n", + " ('UCTSlot', UCTSlot (24 devices) at 0x153d0a190)])" ] }, - "execution_count": 6, + "execution_count": 5, "metadata": {}, "output_type": "execute_result" } @@ -194,7 +205,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 6, "metadata": {}, "outputs": [ { @@ -221,15 +232,12 @@ " idx\n", " u\n", " name\n", + " bus\n", " Vn\n", + " p0\n", + " q0\n", " vmax\n", " vmin\n", - " v0\n", - " a0\n", - " xcoord\n", - " ycoord\n", - " area\n", - " zone\n", " owner\n", " \n", " \n", @@ -244,283 +252,67 @@ " \n", " \n", " \n", - " \n", - " \n", - " \n", " \n", " \n", " \n", " \n", " 0\n", - " 1\n", - " 1.0\n", - " BUS1\n", - " 69.0\n", - " 1.1\n", - " 0.9\n", - " 1.03000\n", - " 0.000000\n", - " 0\n", - " 0\n", - " 1\n", - " 1\n", - " 1\n", - " \n", - " \n", - " 1\n", - " 2\n", - " 1.0\n", - " BUS2\n", - " 69.0\n", - " 1.1\n", - " 0.9\n", - " 1.01970\n", - " -0.027981\n", - " 0\n", - " 0\n", - " 1\n", - " 1\n", - " 1\n", - " \n", - " \n", - " 2\n", - " 3\n", - " 1.0\n", - " BUS3\n", - " 69.0\n", - " 1.1\n", - " 0.9\n", - " 1.00042\n", - " -0.060097\n", - " 0\n", - " 0\n", - " 1\n", - " 1\n", - " 1\n", - " \n", - " \n", - " 3\n", - " 4\n", - " 1.0\n", - " BUS4\n", - " 69.0\n", - " 1.1\n", - " 0.9\n", - " 0.99858\n", - " -0.074721\n", - " 0\n", - " 0\n", - " 1\n", - " 1\n", - " 1\n", - " \n", - " \n", - " 4\n", - " 5\n", + " PQ_1\n", " 1.0\n", - " BUS5\n", - " 69.0\n", - " 1.1\n", - " 0.9\n", - " 1.00443\n", - " -0.064315\n", - " 0\n", - " 0\n", - " 1\n", + " PQ 1\n", " 1\n", - " 1\n", - " \n", - " \n", - " 5\n", - " 6\n", - " 1.0\n", - " BUS6\n", - " 138.0\n", - " 1.1\n", - " 0.9\n", - " 0.99871\n", - " -0.109998\n", - " 0\n", - " 0\n", - " 2\n", - " 2\n", - " 2\n", - " \n", - " \n", - " 6\n", - " 7\n", - " 1.0\n", - " BUS7\n", - " 138.0\n", - " 1.1\n", - " 0.9\n", - " 1.00682\n", - " -0.084285\n", - " 0\n", - " 0\n", - " 2\n", - " 2\n", - " 2\n", - " \n", - " \n", - " 7\n", - " 8\n", - " 1.0\n", - " BUS8\n", - " 69.0\n", + " 230.0\n", + " 3.0\n", + " 0.9861\n", " 1.1\n", " 0.9\n", - " 1.01895\n", - " -0.024339\n", - " 0\n", - " 0\n", - " 2\n", - " 2\n", - " 2\n", - " \n", - " \n", - " 8\n", - " 9\n", - " 1.0\n", - " BUS9\n", - " 138.0\n", - " 1.1\n", - " 0.9\n", - " 1.00193\n", - " -0.127502\n", - " 0\n", - " 0\n", - " 2\n", - " 2\n", - " 2\n", - " \n", - " \n", - " 9\n", - " 10\n", - " 1.0\n", - " BUS10\n", - " 138.0\n", - " 1.1\n", - " 0.9\n", - " 0.99351\n", - " -0.130202\n", - " 0\n", - " 0\n", - " 2\n", - " 2\n", - " 2\n", - " \n", - " \n", - " 10\n", - " 11\n", - " 1.0\n", - " BUS11\n", - " 138.0\n", - " 1.1\n", - " 0.9\n", - " 0.99245\n", - " -0.122948\n", - " 0\n", - " 0\n", - " 2\n", - " 2\n", - " 2\n", + " None\n", " \n", " \n", - " 11\n", - " 12\n", + " 1\n", + " PQ_2\n", " 1.0\n", - " BUS12\n", - " 138.0\n", - " 1.1\n", - " 0.9\n", - " 0.98639\n", - " -0.128934\n", - " 0\n", - " 0\n", - " 2\n", - " 2\n", + " PQ 2\n", " 2\n", - " \n", - " \n", - " 12\n", - " 13\n", - " 1.0\n", - " BUS13\n", - " 138.0\n", + " 230.0\n", + " 3.0\n", + " 0.9861\n", " 1.1\n", " 0.9\n", - " 0.98403\n", - " -0.133786\n", - " 0\n", - " 0\n", - " 2\n", - " 2\n", - " 2\n", + " None\n", " \n", " \n", - " 13\n", - " 14\n", + " 2\n", + " PQ_3\n", " 1.0\n", - " BUS14\n", - " 138.0\n", + " PQ 3\n", + " 3\n", + " 230.0\n", + " 4.0\n", + " 1.3147\n", " 1.1\n", " 0.9\n", - " 0.99063\n", - " -0.166916\n", - " 0\n", - " 0\n", - " 2\n", - " 2\n", - " 2\n", + " None\n", " \n", " \n", "\n", "" ], "text/plain": [ - " idx u name Vn vmax vmin v0 a0 xcoord ycoord \\\n", - "uid \n", - "0 1 1.0 BUS1 69.0 1.1 0.9 1.03000 0.000000 0 0 \n", - "1 2 1.0 BUS2 69.0 1.1 0.9 1.01970 -0.027981 0 0 \n", - "2 3 1.0 BUS3 69.0 1.1 0.9 1.00042 -0.060097 0 0 \n", - "3 4 1.0 BUS4 69.0 1.1 0.9 0.99858 -0.074721 0 0 \n", - "4 5 1.0 BUS5 69.0 1.1 0.9 1.00443 -0.064315 0 0 \n", - "5 6 1.0 BUS6 138.0 1.1 0.9 0.99871 -0.109998 0 0 \n", - "6 7 1.0 BUS7 138.0 1.1 0.9 1.00682 -0.084285 0 0 \n", - "7 8 1.0 BUS8 69.0 1.1 0.9 1.01895 -0.024339 0 0 \n", - "8 9 1.0 BUS9 138.0 1.1 0.9 1.00193 -0.127502 0 0 \n", - "9 10 1.0 BUS10 138.0 1.1 0.9 0.99351 -0.130202 0 0 \n", - "10 11 1.0 BUS11 138.0 1.1 0.9 0.99245 -0.122948 0 0 \n", - "11 12 1.0 BUS12 138.0 1.1 0.9 0.98639 -0.128934 0 0 \n", - "12 13 1.0 BUS13 138.0 1.1 0.9 0.98403 -0.133786 0 0 \n", - "13 14 1.0 BUS14 138.0 1.1 0.9 0.99063 -0.166916 0 0 \n", - "\n", - " area zone owner \n", - "uid \n", - "0 1 1 1 \n", - "1 1 1 1 \n", - "2 1 1 1 \n", - "3 1 1 1 \n", - "4 1 1 1 \n", - "5 2 2 2 \n", - "6 2 2 2 \n", - "7 2 2 2 \n", - "8 2 2 2 \n", - "9 2 2 2 \n", - "10 2 2 2 \n", - "11 2 2 2 \n", - "12 2 2 2 \n", - "13 2 2 2 " + " idx u name bus Vn p0 q0 vmax vmin owner\n", + "uid \n", + "0 PQ_1 1.0 PQ 1 1 230.0 3.0 0.9861 1.1 0.9 None\n", + "1 PQ_2 1.0 PQ 2 2 230.0 3.0 0.9861 1.1 0.9 None\n", + "2 PQ_3 1.0 PQ 3 3 230.0 4.0 1.3147 1.1 0.9 None" ] }, - "execution_count": 7, + "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "sp.Bus.as_df()" + "sp.PQ.as_df()" ] }, { @@ -528,27 +320,33 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Similarly, in AMS, all supported routines are registered to an OrderedDict ``routines``." + "In AMS, all supported routines are registered to an OrderedDict ``routines``." ] }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "OrderedDict([('DCPF', DCPF at 0x1590a8fa0),\n", - " ('PFlow', PFlow at 0x159a6c670),\n", - " ('ACOPF', ACOPF at 0x159b5f220),\n", - " ('DCOPF', DCOPF at 0x159b9b1c0),\n", - " ('ED', ED at 0x159b9b670),\n", - " ('RTED', RTED at 0x159b9be50),\n", - " ('UC', UC at 0x159bb1b50)])" + "OrderedDict([('DCPF', DCPF at 0x153c7bcd0),\n", + " ('PFlow', PFlow at 0x153d0adc0),\n", + " ('CPF', CPF at 0x153d193d0),\n", + " ('ACOPF', ACOPF at 0x153d19a30),\n", + " ('DCOPF', DCOPF at 0x153d31370),\n", + " ('ED', ED at 0x153d45040),\n", + " ('EDES', EDES at 0x153d6d100),\n", + " ('RTED', RTED at 0x153d93220),\n", + " ('RTEDES', RTEDES at 0x153da2430),\n", + " ('UC', UC at 0x153db6f10),\n", + " ('UCES', UCES at 0x16816ce20),\n", + " ('DOPF', DOPF at 0x16819e490),\n", + " ('DOPFVIS', DOPFVIS at 0x1681b35b0)])" ] }, - "execution_count": 8, + "execution_count": 7, "metadata": {}, "output_type": "execute_result" } @@ -562,522 +360,317 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### Solve Power Flow" + "### Solve an Routine" ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "In AMS, the AC power flow and DC power flow are solved by PYPOWER ``runpf()`` and ``rundcpf()`` functions, respectively." + "Before solving an routine, we need to initialize it first.\n", + "Here Real-time Economic Dispatch (RTED) is used as an example." ] }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ - "Setup model for DCPF\n", - "DCPF has no objective function.\n", - "DCPF model set up in 0.0011 seconds.\n", - "PYPOWER Version 5.1.4, 27-June-2018\n", - " -- DC Power Flow\n", - "\n", - "DCPF completed in 0.0033 seconds with exit code 0.\n" + "Routine initialized in 0.0130 seconds.\n" ] }, { "data": { "text/plain": [ - "0" + "True" ] }, - "execution_count": 9, + "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "sp.DCPF.run()" + "sp.RTED.init()" ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "Inspect the generator power outputs and bus angles." - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "array([0.4 , 0.4 , 0.3 , 0.35 , 0.787])" - ] - }, - "execution_count": 10, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "sp.DCPF.pg.v" + "Then, one can solve it by calling ``run()``.\n", + "Here, argument `solver` can be passed to specify the solver to use, such as `solver='ECOS'`.\n", + "\n", + "Installed solvers can be listed by ``ams.shared.INSTALLED_SOLVERS``,\n", + "and more detailes of solver can be found at [CVXPY-Choosing a solver](https://www.cvxpy.org/tutorial/advanced/index.html#choosing-a-solver)." ] }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "array([ 0. , -0.02887178, -0.0641022 , -0.07743542, -0.06670094,\n", - " -0.11334928, -0.08549152, -0.02403816, -0.12824676, -0.13308983,\n", - " -0.12681192, -0.1339144 , -0.13779051, -0.16285214])" + "['CLARABEL',\n", + " 'CVXOPT',\n", + " 'ECOS',\n", + " 'ECOS_BB',\n", + " 'GLPK',\n", + " 'GLPK_MI',\n", + " 'GUROBI',\n", + " 'OSQP',\n", + " 'SCIPY',\n", + " 'SCS']" ] }, - "execution_count": 11, + "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "sp.DCPF.aBus.v" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Solve the AC power flow." + "ams.shared.INSTALLED_SOLVERS" ] }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ - "Setup model for PFlow\n", - "PFlow has no objective function.\n", - "PFlow model set up in 0.0012 seconds.\n", - "PYPOWER Version 5.1.4, 27-June-2018\n", - "\n", - "Newton's method power flow converged in 3 iterations.\n", - "\n", - "PFlow completed in 0.0091 seconds with exit code 0.\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " -- AC Power Flow (Newton)\n", - "\n" + "RTED solved as optimal in 0.0193 seconds, converged after 9 iterations using solver ECOS.\n" ] }, { "data": { "text/plain": [ - "0" + "True" ] }, - "execution_count": 12, + "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "sp.PFlow.run()" + "sp.RTED.run(solver='ECOS')" ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "Inspect the generator power outputs, bus angles, and bus voltages." + "The solved results are stored in each variable itself.\n", + "For example, the solved power generation of ten generators\n", + "are stored in ``pg.v``." ] }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "array([0.4 , 0.4 , 0.3 , 0.35 , 0.81427213])" + "array([2.1, 5.2, 0.7, 2. ])" ] }, - "execution_count": 13, + "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "sp.PFlow.pg.v" + "sp.RTED.pg.v" ] }, { - "cell_type": "code", - "execution_count": 14, + "cell_type": "markdown", "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "array([ 0. , -0.03078883, -0.06173455, -0.07696511, -0.06707345,\n", - " -0.11262148, -0.08526268, -0.02687732, -0.12646404, -0.12942483,\n", - " -0.12356406, -0.13042896, -0.13475262, -0.16547667])" - ] - }, - "execution_count": 14, - "metadata": {}, - "output_type": "execute_result" - } - ], "source": [ - "sp.PFlow.aBus.v" + "Here, ``get_idx()`` can be used to get the index of a variable." ] }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "array([1.03 , 1.03 , 1.01 , 1.01140345, 1.01725551,\n", - " 1.03 , 1.0224715 , 1.03 , 1.02176879, 1.01554207,\n", - " 1.01911515, 1.01740699, 1.01445023, 1.0163402 ])" + "['PV_1', 'PV_3', 'PV_5', 'Slack_4']" ] }, - "execution_count": 15, + "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "sp.PFlow.vBus.v" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Solve ACOPF" + "sp.RTED.pg.get_idx()" ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "In AMS, the ACOPF is solved by PYPOWER ``runopf()`` function." + "Part of the solved results can be accessed with given indices." ] }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 13, "metadata": {}, "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "Setup model for ACOPF\n", - "ACOPF model set up in 0.0015 seconds.\n", - "PYPOWER Version 5.1.4, 27-June-2018\n", - " -- AC Optimal Power Flow\n", - "\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Python Interior Point Solver - PIPS, Version 1.0, 07-Feb-2011\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "ACOPF completed in 0.2583 seconds with exit code 0.\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Converged!\n" - ] - }, { "data": { "text/plain": [ - "0" + "array([2.1, 5.2])" ] }, - "execution_count": 16, + "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "sp.ACOPF.run()" + "sp.RTED.get(src='pg', attr='v', idx=['PV_1', 'PV_3'])" ] }, { - "cell_type": "code", - "execution_count": 17, + "cell_type": "markdown", "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "array([0.49996655, 0.10000065, 0.10000064, 0.10000065, 1.49809555])" - ] - }, - "execution_count": 17, - "metadata": {}, - "output_type": "execute_result" - } - ], "source": [ - "sp.ACOPF.pg.v" + "All Vars are listed in an OrderedDict ``vars``." ] }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "array([0.14990061, 0.14999993, 0.09997174, 0.09993129, 0.08059299])" + "OrderedDict([('pg', Var: StaticGen.pg),\n", + " ('aBus', Var: Bus.aBus),\n", + " ('plf', Var: Line.plf),\n", + " ('pru', Var: StaticGen.pru),\n", + " ('prd', Var: StaticGen.prd)])" ] }, - "execution_count": 18, + "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "sp.ACOPF.qg.v" + "sp.RTED.vars" ] }, { - "cell_type": "code", - "execution_count": 19, + "cell_type": "markdown", "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "array([ 0. , -0.04662694, -0.11452421, -0.11585264, -0.09988879,\n", - " -0.18205907, -0.16148797, -0.14627229, -0.19500474, -0.19788944,\n", - " -0.19253337, -0.19878038, -0.20322846, -0.23228802])" - ] - }, - "execution_count": 19, - "metadata": {}, - "output_type": "execute_result" - } - ], "source": [ - "sp.ACOPF.aBus.v" + "The Objective value can be accessed with ``obj.v``." ] }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "array([1.09999682, 1.08133094, 1.05126009, 1.05401909, 1.06111871,\n", - " 1.05714675, 1.06795632, 1.08056893, 1.06468318, 1.05600708,\n", - " 1.05321342, 1.04626446, 1.04476777, 1.055634 ])" + "2287.0833333698793" ] }, - "execution_count": 20, + "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "sp.ACOPF.vBus.v" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Solve DCOPF" + "sp.RTED.obj.v" ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "In AMS, DCOPF and other routines are modeled and solved with CVXPY." - ] - }, - { - "cell_type": "code", - "execution_count": 21, - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "Setup model of DCOPF\n", - "DCOPF model set up in 0.0017 seconds.\n", - "DCOPF solved as optimal in 0.0055 seconds with exit code 0.\n" - ] - } - ], - "source": [ - "sp.DCOPF.run()" + "Similarly, all Constrs are listed in an OrderedDict ``constrs``,\n", + "and the expression values can also be accessed." ] }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "array([0.5 , 0.5 , 0.3685, 0.3685, 0.5 ])" + "OrderedDict([('pglb', [ON]: -pg + mul(nctrle, pg0) + mul(ctrle, pmin) <=0),\n", + " ('pgub', [ON]: pg - mul(nctrle, pg0) - mul(ctrle, pmax) <=0),\n", + " ('plflb', [ON]: -plf - rate_a <=0),\n", + " ('plfub', [ON]: plf - rate_a <=0),\n", + " ('pb', [ON]: sum(pd) - sum(pg) =0),\n", + " ('aref', [ON]: Cs@aBus =0),\n", + " ('pnb', [ON]: PTDF@(Cgi@pg - Cli@pd) - plf =0),\n", + " ('rbu', [ON]: gs @ mul(ug, pru) - dud =0),\n", + " ('rbd', [ON]: gs @ mul(ug, prd) - ddd =0),\n", + " ('rru', [ON]: mul(ug, pg + pru) - mul(ug, pmax) <=0),\n", + " ('rrd', [ON]: mul(ug, -pg + prd) - mul(ug, pmin) <=0),\n", + " ('rgu', [ON]: mul(ug, pg-pg0-R10) <=0),\n", + " ('rgd', [ON]: mul(ug, -pg+pg0-R10) <=0)])" ] }, - "execution_count": 22, + "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "sp.DCOPF.pg.v" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "One can specify the solver by passing the solver name to the ``solver`` argument of ``run()``." - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The default CVXPV solver ``ECOS`` might be slow or incapable of solving large-scale problems." - ] - }, - { - "cell_type": "code", - "execution_count": 23, - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "Parsing input file \"/Users/jinningwang/Documents/work/ams/ams/cases/matpower/case_ACTIVSg2000.m\"...\n", - "Input file parsed in 0.4876 seconds.\n", - "The bus index is not continuous, adjusted automatically.\n", - "System set up in 0.8002 seconds.\n" - ] - } - ], - "source": [ - "sp = ams.load(ams.get_case('matpower/case_ACTIVSg2000.m'), default_config=True)" + "sp.RTED.constrs" ] }, { - "cell_type": "code", - "execution_count": 24, - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "Setup model of DCOPF\n", - "DCOPF model set up in 0.0094 seconds.\n", - "DCOPF solved as optimal in 29.8392 seconds with exit code 0.\n" - ] - } - ], - "source": [ - "sp.DCOPF.run(solver='ECOS')" - ] - }, - { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "One can specify the solver by passing the solver name to the ``solver`` argument of ``run()``." + "One can also inspect the Constr values." ] }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 17, "metadata": {}, "outputs": [ { - "name": "stderr", - "output_type": "stream", - "text": [ - "DCOPF model set up in 0.0129 seconds.\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Set parameter Username\n", - "Academic license - for non-commercial use only - expires 2024-05-21\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "DCOPF solved as optimal in 0.7355 seconds with exit code 0.\n" - ] + "data": { + "text/plain": [ + "array([ -997.9 , -997.0349, -1002.9651, -997. ])" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" } ], "source": [ - "sp.DCOPF.run(solver='GUROBI')" + "sp.RTED.rgu.v" ] } ], @@ -1097,7 +690,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.16" + "version": "3.9.18" }, "orig_nbformat": 4, "vscode": { diff --git a/examples/ex2.ipynb b/examples/ex2.ipynb new file mode 100644 index 00000000..839a1956 --- /dev/null +++ b/examples/ex2.ipynb @@ -0,0 +1,1222 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Manipulate the Dispatch Model" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This example shows how to manipulate the model,\n", + "such as trip a generator, change the load, etc." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import ams\n", + "\n", + "import datetime" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Last run time: 2023-11-28 22:06:37\n", + "ams:0.7.3.post74.dev0+g47c1ae3\n" + ] + } + ], + "source": [ + "print(\"Last run time:\", datetime.datetime.now().strftime(\"%Y-%m-%d %H:%M:%S\"))\n", + "\n", + "print(f'ams:{ams.__version__}')" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "ams.config_logger(stream_level=20)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Run Simulations" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Load Case" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Parsing input file \"/Users/jinningwang/Documents/work/ams/ams/cases/5bus/pjm5bus_uced.xlsx\"...\n", + "Input file parsed in 0.1106 seconds.\n", + "System set up in 0.0034 seconds.\n" + ] + } + ], + "source": [ + "sp = ams.load(ams.get_case('5bus/pjm5bus_uced.xlsx'),\n", + " setup=True,\n", + " no_output=True,)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The system load are defined in model `PQ`." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
idxunamebusVnp0q0vmaxvminowner
uid
0PQ_11.0PQ 11230.03.00.98611.10.9None
1PQ_21.0PQ 22230.03.00.98611.10.9None
2PQ_31.0PQ 33230.04.01.31471.10.9None
\n", + "
" + ], + "text/plain": [ + " idx u name bus Vn p0 q0 vmax vmin owner\n", + "uid \n", + "0 PQ_1 1.0 PQ 1 1 230.0 3.0 0.9861 1.1 0.9 None\n", + "1 PQ_2 1.0 PQ 2 2 230.0 3.0 0.9861 1.1 0.9 None\n", + "2 PQ_3 1.0 PQ 3 3 230.0 4.0 1.3147 1.1 0.9 None" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sp.PQ.as_df()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For simplicity, PQ load is typically marked as `pd` in RTED." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([3., 3., 4.])" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sp.RTED.pd.v" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Run Simulation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "RTED can be solved and one can inspect the results as discussed in\n", + "previous example." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Routine initialized in 0.0088 seconds.\n", + "RTED solved as optimal in 0.0170 seconds, converged after 9 iterations using solver ECOS.\n" + ] + }, + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sp.RTED.run(solver='ECOS')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Power generation (`pg`) and line flow (`plf`) can be accessed as follows." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([2.1, 5.2, 0.7, 2. ])" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sp.RTED.pg.v" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 0.70595331, 0.68616798, 0.00192539, -1.58809337, 0.61190663,\n", + " -0.70192539, 0.70595331])" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sp.RTED.plf.v" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Change Load" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The load values can be manipulated in the model `PQ`." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sp.PQ.set(src='p0', attr='v', idx=['PQ_1', 'PQ_2'], value=[3.2, 3.2])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "According parameters need to be updated to make the changes effective.\n", + "If not sure which parameters need to be updated, one can use\n", + "``update()`` to update all parameters." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sp.RTED.update('pd')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "After manipulation, the routined can be solved again." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "RTED solved as optimal in 0.0016 seconds, converged after 8 iterations using solver ECOS.\n" + ] + }, + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sp.RTED.run(solver='ECOS')" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([2.1 , 5.19999999, 1.10000002, 1.99999999])" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sp.RTED.pg.v" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Trip a Generator" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can see that there are three PV generators in the system." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
idxunameSnVnbusbusrp0q0pmax...Qc2minQc2maxRagcR10R30Rqapfpg0td1td2
uid
0PV_11.0Alta100.0230.00None1.00000.02.1...0.00.0999.0999.0999.0999.00.00.00.50.0
1PV_31.0Solitude100.0230.02None3.23490.05.2...0.00.0999.0999.0999.0999.00.00.00.50.0
2PV_51.0Brighton100.0230.04None4.66510.06.0...0.00.0999.0999.0999.0999.00.00.00.50.0
\n", + "

3 rows × 33 columns

\n", + "
" + ], + "text/plain": [ + " idx u name Sn Vn bus busr p0 q0 pmax ... \\\n", + "uid ... \n", + "0 PV_1 1.0 Alta 100.0 230.0 0 None 1.0000 0.0 2.1 ... \n", + "1 PV_3 1.0 Solitude 100.0 230.0 2 None 3.2349 0.0 5.2 ... \n", + "2 PV_5 1.0 Brighton 100.0 230.0 4 None 4.6651 0.0 6.0 ... \n", + "\n", + " Qc2min Qc2max Ragc R10 R30 Rq apf pg0 td1 td2 \n", + "uid \n", + "0 0.0 0.0 999.0 999.0 999.0 999.0 0.0 0.0 0.5 0.0 \n", + "1 0.0 0.0 999.0 999.0 999.0 999.0 0.0 0.0 0.5 0.0 \n", + "2 0.0 0.0 999.0 999.0 999.0 999.0 0.0 0.0 0.5 0.0 \n", + "\n", + "[3 rows x 33 columns]" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sp.PV.as_df()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "`PV_1` is tripped by setting its connection status `u` to 0." + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sp.StaticGen.set(src='u', attr='v', idx='PV_1', value=0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In AMS, some parameters are defiend as constants in the numerical optimization model\n", + "to follow the CVXPY DCP and DPP rules.\n", + "Once non-parametric parameters are changed, the model needs to be resetup\n", + "to make the changes effective.\n", + "\n", + "More details can be found at [CVXPY - Disciplined Convex Programming](https://www.cvxpy.org/tutorial/dcp/index.html#disciplined-convex-programming)." + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Resetup RTED OModel due to non-parametric change.\n" + ] + }, + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sp.RTED.update()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Then we can resolve the model." + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "RTED solved as optimal in 0.0171 seconds, converged after 8 iterations using solver ECOS.\n" + ] + }, + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sp.RTED.run(solver='ECOS')" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([0. , 5.2 , 3.20000001, 2. ])" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sp.RTED.pg.v" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 0.79973461, 0.56088965, -2.16035886, -1.60053079, 0.39946921,\n", + " -1.03964115, 0.79973461])" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sp.RTED.plf.v" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Trip a Line" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can inspect the `Line` model to check the system topology." + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
idxunamebus1bus2SnfnVn1Vn2r...tapphirate_arate_brate_cownerxcoordycoordaminamax
uid
001.0Line 1-201100.060.0230.0230.00.00281...1.00.0999.0999.0999.0NoneNoneNone-6.2831856.283185
111.0Line 1-403100.060.0230.0230.00.00304...1.00.0999.0999.0999.0NoneNoneNone-6.2831856.283185
221.0Line 1-504100.060.0230.0230.00.00064...1.00.0999.0999.0999.0NoneNoneNone-6.2831856.283185
331.0Line 2-312100.060.0230.0230.00.00108...1.00.0999.0999.0999.0NoneNoneNone-6.2831856.283185
441.0Line 3-423100.060.0230.0230.00.00297...1.00.0999.0999.0999.0NoneNoneNone-6.2831856.283185
551.0Line 4-534100.060.0230.0230.00.00297...1.00.0999.0999.0999.0NoneNoneNone-6.2831856.283185
661.0Line 1-2 (2)01100.060.0230.0230.00.00281...1.00.0999.0999.0999.0NoneNoneNone-6.2831856.283185
\n", + "

7 rows × 28 columns

\n", + "
" + ], + "text/plain": [ + " idx u name bus1 bus2 Sn fn Vn1 Vn2 r \\\n", + "uid \n", + "0 0 1.0 Line 1-2 0 1 100.0 60.0 230.0 230.0 0.00281 \n", + "1 1 1.0 Line 1-4 0 3 100.0 60.0 230.0 230.0 0.00304 \n", + "2 2 1.0 Line 1-5 0 4 100.0 60.0 230.0 230.0 0.00064 \n", + "3 3 1.0 Line 2-3 1 2 100.0 60.0 230.0 230.0 0.00108 \n", + "4 4 1.0 Line 3-4 2 3 100.0 60.0 230.0 230.0 0.00297 \n", + "5 5 1.0 Line 4-5 3 4 100.0 60.0 230.0 230.0 0.00297 \n", + "6 6 1.0 Line 1-2 (2) 0 1 100.0 60.0 230.0 230.0 0.00281 \n", + "\n", + " ... tap phi rate_a rate_b rate_c owner xcoord ycoord amin \\\n", + "uid ... \n", + "0 ... 1.0 0.0 999.0 999.0 999.0 None None None -6.283185 \n", + "1 ... 1.0 0.0 999.0 999.0 999.0 None None None -6.283185 \n", + "2 ... 1.0 0.0 999.0 999.0 999.0 None None None -6.283185 \n", + "3 ... 1.0 0.0 999.0 999.0 999.0 None None None -6.283185 \n", + "4 ... 1.0 0.0 999.0 999.0 999.0 None None None -6.283185 \n", + "5 ... 1.0 0.0 999.0 999.0 999.0 None None None -6.283185 \n", + "6 ... 1.0 0.0 999.0 999.0 999.0 None None None -6.283185 \n", + "\n", + " amax \n", + "uid \n", + "0 6.283185 \n", + "1 6.283185 \n", + "2 6.283185 \n", + "3 6.283185 \n", + "4 6.283185 \n", + "5 6.283185 \n", + "6 6.283185 \n", + "\n", + "[7 rows x 28 columns]" + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sp.Line.as_df()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here line `1` is tripped by setting its connection status `u` to 0." + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 21, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sp.Line.set(src='u', attr='v', idx=0, value=0)" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Resetup RTED OModel due to non-parametric change.\n" + ] + }, + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sp.RTED.update()" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "RTED solved as optimal in 0.0190 seconds, converged after 8 iterations using solver ECOS.\n" + ] + }, + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 23, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sp.RTED.run(solver='ECOS')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here we can see the tripped line has no flow." + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([-0. , 0.70423831, -2.03964421, -1.8645941 , 0.13540589,\n", + " -1.16035581, 1.3354059 ])" + ] + }, + "execution_count": 24, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sp.RTED.plf.v" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "ams", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.18" + }, + "orig_nbformat": 4, + "vscode": { + "interpreter": { + "hash": "d2b3bf80176349caa68dc4a3c77bd06eaade8abc678330f7d1c813c53380e5d2" + } + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/requirements-extra.txt b/requirements-extra.txt index 9a60f48f..48a96c77 100644 --- a/requirements-extra.txt +++ b/requirements-extra.txt @@ -8,17 +8,13 @@ # Only one `#` is allowed per line. Lines starting with `#` are ignored. pytest==7.0.1 # dev -flake8 # dev -sphinx # dev, doc -pydata-sphinx-theme # dev, doc -numpydoc # dev, doc -sphinx-copybutton # dev, doc -sphinx-panels # dev, doc -pydata-sphinx-theme # dev, doc -myst-parser==0.15.2 # dev, doc -myst-nb # dev, doc -graphviz # doc -pydeps # doc -scip # opt -pyscipopt # opt -gurobipy # opt \ No newline at end of file +flake8 # dev +igraph # dev +sphinx # dev, doc +pydata-sphinx-theme # dev, doc +numpydoc # dev, doc +sphinx-copybutton # dev, doc +sphinx-panels # dev, doc +pydata-sphinx-theme # dev, doc +myst-parser==0.15.2 # dev, doc +myst-nb # dev, doc \ No newline at end of file diff --git a/requirements.txt b/requirements.txt index 78f263c5..fad7a13c 100644 --- a/requirements.txt +++ b/requirements.txt @@ -5,4 +5,5 @@ matplotlib psutil openpyxl andes -cvxpy \ No newline at end of file +cvxpy +cvxopt \ No newline at end of file diff --git a/tests/test_case.py b/tests/test_case.py index 625fdfd4..a9d8be54 100644 --- a/tests/test_case.py +++ b/tests/test_case.py @@ -237,8 +237,8 @@ def test_ieee39_esd1_init(self): default_config=True, no_output=True, ) - ss.ED2.init() - ss.UC2.init() + ss.EDES.init() + ss.UCES.init() - self.assertEqual(ss.ED2.exit_code, 0, "Exit code is not 0.") - self.assertEqual(ss.UC2.exit_code, 0, "Exit code is not 0.") + self.assertEqual(ss.EDES.exit_code, 0, "Exit code is not 0.") + self.assertEqual(ss.UCES.exit_code, 0, "Exit code is not 0.") diff --git a/tests/test_omodel.py b/tests/test_omodel.py new file mode 100644 index 00000000..5ead7cb5 --- /dev/null +++ b/tests/test_omodel.py @@ -0,0 +1,64 @@ +import unittest +import numpy as np + +import ams +import cvxpy as cp + + +def require_MIP_solver(f): + """ + Decorator for skipping tests that require MIP solver. + """ + def wrapper(*args, **kwargs): + all_solvers = cp.installed_solvers() + mip_solvers = ['CBC', 'COPT', 'GLPK_MI', 'CPLEX', 'GUROBI', + 'MOSEK', 'SCIP', 'XPRESS', 'SCIPY'] + if any(s in mip_solvers for s in all_solvers): + pass + else: + raise unittest.SkipTest("MIP solver is not available.") + return f(*args, **kwargs) + return wrapper + + +class TestOMdel(unittest.TestCase): + """ + Test methods of `OModel`. + """ + def setUp(self) -> None: + self.ss = ams.load(ams.get_case("ieee39/ieee39_uced_esd1.xlsx"), + default_config=True, + no_output=True, + ) + + def test_var_access_brfore_solve(self): + """ + Test `Var` access before solve. + """ + self.ss.DCOPF.init() + self.assertIsNone(self.ss.DCOPF.pg.v) + + def test_var_access_after_solve(self): + """ + Test `Var` access after solve. + """ + self.ss.DCOPF.run() + np.testing.assert_equal(self.ss.DCOPF.pg.v, + self.ss.StaticGen.get(src='p', attr='v', + idx=self.ss.DCOPF.pg.get_idx())) + + def test_constr_access_brfore_solve(self): + """ + Test `Constr` access before solve. + """ + self.ss.DCOPF.init(force=True) + np.testing.assert_equal(self.ss.DCOPF.plflb.v, None) + + def test_constr_access_after_solve(self): + """ + Test `Constr` access after solve. + """ + self.ss.DCOPF.run() + self.assertIsInstance(self.ss.DCOPF.plflb.v, np.ndarray) + + # NOTE: add Var, Constr add functions diff --git a/tests/test_routine.py b/tests/test_routine.py index 6da09bfe..7d317bc1 100644 --- a/tests/test_routine.py +++ b/tests/test_routine.py @@ -1,29 +1,37 @@ import unittest +from functools import wraps import numpy as np import ams -import cvxpy as cp +try: + from ams.shared import igraph + getattr(igraph, '__version__') + HAVE_IGRAPH = True +except (ImportError, AttributeError): + HAVE_IGRAPH = False -def require_MIP_solver(f): + +def require_igraph(f): """ - Decorator for skipping tests that require MIP solver. + Decorator for functions that require igraph. """ - def wrapper(*args, **kwargs): - all_solvers = cp.installed_solvers() - mip_solvers = ['CBC', 'COPT', 'GLPK_MI', 'CPLEX', 'GUROBI', - 'MOSEK', 'SCIP', 'XPRESS', 'SCIPY'] - if any(s in mip_solvers for s in all_solvers): - pass - else: - raise unittest.SkipTest("MIP solver is not available.") - return f(*args, **kwargs) + + @wraps(f) + def wrapper(*args, **kwds): + try: + getattr(igraph, '__version__') + except AttributeError: + raise ModuleNotFoundError("igraph needs to be manually installed.") + + return f(*args, **kwds) + return wrapper class TestRoutineMethods(unittest.TestCase): """ - Test methods of Routine. + Test methods of `Routine`. """ def setUp(self) -> None: self.ss = ams.load(ams.get_case("ieee39/ieee39_uced_esd1.xlsx"), @@ -44,47 +52,101 @@ def test_routine_get(self): Test `Routine.get()` method. """ - # get a rparam value + # get an rparam value np.testing.assert_equal(self.ss.DCOPF.get('ug', 'PV_30'), 1) - # before solving, vars values are not available - self.assertRaises(KeyError, self.ss.DCOPF.get, 'pg', 'PV_30') - - self.ss.DCOPF.run(solver='OSQP') + # get an unpacked var value + self.ss.DCOPF.run() self.assertEqual(self.ss.DCOPF.exit_code, 0, "Exit code is not 0.") np.testing.assert_equal(self.ss.DCOPF.get('pg', 'PV_30', 'v'), self.ss.StaticGen.get('p', 'PV_30', 'v')) - def test_RTED(self): - """ - Test `RTED.run()`. - """ - self.ss.RTED.run(solver='OSQP') - self.assertEqual(self.ss.RTED.exit_code, 0, "Exit code is not 0.") +class TestRoutineInit(unittest.TestCase): + """ + Test solving routines. + """ + def setUp(self) -> None: + self.s1 = ams.load(ams.get_case("ieee39/ieee39_uced_esd1.xlsx"), + default_config=True, + no_output=True, + ) + self.s2 = ams.load(ams.get_case("ieee123/ieee123_regcv1.xlsx"), + default_config=True, + no_output=True, + ) - @require_MIP_solver - def test_RTED2(self): + def test_Init(self): """ - Test `RTED2.run()`. + Test `routine.init()`. """ + # NOTE: for DED, using ieee123 as a test case + for rtn in self.s1.routines.values(): + if not rtn._data_check(): + rtn = getattr(self.s2, rtn.class_name) + if not rtn._data_check(): + continue # TODO: here should be a warning? + self.assertTrue(rtn.init(force=True), + f"{rtn.class_name} initialization failed!") + + +@unittest.skipUnless(HAVE_IGRAPH, "igraph not available") +class TestRoutineGraph(unittest.TestCase): + """ + Test routine graph. + """ - self.ss.RTED2.run() - self.assertEqual(self.ss.RTED2.exit_code, 0, "Exit code is not 0.") - - def test_ED(self): + def test_5bus_graph(self): """ - Test `ED.run()`. + Test routine graph of PJM 5-bus system. """ - - self.ss.ED.run(solver='OSQP') - self.assertEqual(self.ss.ED.exit_code, 0, "Exit code is not 0.") - - @require_MIP_solver - def test_ED2(self): + ss = ams.load(ams.get_case("5bus/pjm5bus_uced.xlsx"), + default_config=True, + no_output=True, + ) + _, g = ss.DCOPF.igraph() + self.assertGreaterEqual(np.min(g.degree()), 1) + + def test_ieee14_graph(self): """ - Test `ED2.run()`. + Test routine graph of IEEE 14-bus system. """ - - self.ss.ED2.run() - self.assertEqual(self.ss.ED2.exit_code, 0, "Exit code is not 0.") + ss = ams.load(ams.get_case("ieee14/ieee14_uced.xlsx"), + default_config=True, + no_output=True, + ) + _, g = ss.DCOPF.igraph() + self.assertGreaterEqual(np.min(g.degree()), 1) + + def test_ieee39_graph(self): + """ + Test routine graph of IEEE 39-bus system. + """ + ss = ams.load(ams.get_case("ieee39/ieee39_uced_esd1.xlsx"), + default_config=True, + no_output=True, + ) + _, g = ss.DCOPF.igraph() + self.assertGreaterEqual(np.min(g.degree()), 1) + + def test_npcc_graph(self): + """ + Test routine graph of NPCC 140-bus system. + """ + ss = ams.load(ams.get_case("npcc/npcc_uced.xlsx"), + default_config=True, + no_output=True, + ) + _, g = ss.DCOPF.igraph() + self.assertGreaterEqual(np.min(g.degree()), 1) + + def test_wecc_graph(self): + """ + Test routine graph of WECC 179-bus system. + """ + ss = ams.load(ams.get_case("wecc/wecc_uced.xlsx"), + default_config=True, + no_output=True, + ) + _, g = ss.DCOPF.igraph() + self.assertGreaterEqual(np.min(g.degree()), 1) diff --git a/tests/test_service.py b/tests/test_service.py new file mode 100644 index 00000000..44f37e04 --- /dev/null +++ b/tests/test_service.py @@ -0,0 +1,126 @@ +import unittest +import numpy as np + +from scipy.sparse import csr_matrix as c_sparse # NOQA +from scipy.sparse import lil_matrix as l_sparse # NOQA + +import ams +from ams.core.matprocessor import MatProcessor, MParam +from ams.core.service import NumOp, NumOpDual, ZonalSum + + +class TestMatProcessor(unittest.TestCase): + """ + Test functionality of MatProcessor. + """ + + def setUp(self) -> None: + self.ss = ams.load( + ams.get_case("ieee39/ieee39_uced_esd1.xlsx"), + default_config=True, + no_output=True, + ) + self.nR = self.ss.Region.n + self.nB = self.ss.Bus.n + self.nl = self.ss.Line.n + + self.mats = MatProcessor(self.ss) + self.mats.make() + + def test_MParam(self): + """ + Test `MParam`. + """ + one_vec = MParam(v=c_sparse(np.ones(self.ss.Bus.n))) + # check if `_v` is `c_sparse` instance + self.assertIsInstance(one_vec._v, c_sparse) + # check if `v` is 1D-array + self.assertEqual(one_vec.v.shape, (self.ss.Bus.n,)) + + def test_PTDF(self): + """ + Test `PTDF`. + """ + self.assertIsInstance(self.mats.PTDF.v, np.ndarray) + self.assertEqual(self.mats.PTDF.v.shape, (self.nl, self.nB)) + + def test_Cft(self): + """ + Test `Cft`. + """ + self.assertIsInstance(self.mats.Cft._v, (c_sparse, l_sparse)) + self.assertIsInstance(self.mats.Cft.v, np.ndarray) + self.assertEqual(self.mats.Cft.v.max(), 1) + + def test_pl(self): + """ + Test `pl`. + """ + self.assertEqual(self.mats.pl.v.ndim, 1) + + +class TestService(unittest.TestCase): + """ + Test functionality of Services. + """ + + def setUp(self) -> None: + self.ss = ams.load(ams.get_case("ieee39/ieee39_uced_esd1.xlsx"), + default_config=True, + no_output=True,) + self.nR = self.ss.Region.n + self.nB = self.ss.Bus.n + self.nL = self.ss.Line.n + self.nD = self.ss.StaticLoad.n # number of static loads + + def test_NumOp_norfun(self): + """ + Test `NumOp` without return function. + """ + CftT = NumOp(u=self.ss.mats.Cft, fun=np.transpose) + np.testing.assert_array_equal(CftT.v.transpose(), self.ss.mats.Cft.v) + + def test_NumOp_rfun(self): + """ + Test `NumOp` with return function. + """ + CftTT = NumOp(u=self.ss.mats.Cft, fun=np.transpose, rfun=np.transpose) + np.testing.assert_array_equal(CftTT.v, self.ss.mats.Cft.v) + + def test_NumOp_ArrayOut(self): + """ + Test `NumOp` non-array output. + """ + M = NumOp(u=self.ss.PV.pmax, + fun=np.max, + rfun=np.dot, rargs=dict(b=10), + array_out=True,) + M2 = NumOp(u=self.ss.PV.pmax, + fun=np.max, + rfun=np.dot, rargs=dict(b=10), + array_out=False,) + self.assertIsInstance(M.v, np.ndarray) + self.assertIsInstance(M2.v, (int, float)) + + def test_NumOpDual(self): + """ + Test `NumOpDual`. + """ + p_vec = MParam(v=self.ss.mats.pl.v) + one_vec = MParam(v=np.ones(self.ss.Bus.n)) + p_sum = NumOpDual(u=p_vec, u2=one_vec, + fun=np.multiply, rfun=np.sum) + self.assertEqual(p_sum.v, self.ss.PQ.p0.v.sum()) + + def test_ZonalSum(self): + """ + Test `ZonalSum`. + """ + ds = ZonalSum(u=self.ss.RTED.zd, zone="Region", + name="ds", tex_name=r"S_{d}", + info="Sum pl vector in shape of zone",) + ds.rtn = self.ss.RTED + # check if the shape is correct + np.testing.assert_array_equal(ds.v.shape, (self.nR, self.nD)) + # check if the values are correct + self.assertTrue(np.all(ds.v.sum(axis=1) <= np.array([self.nB, self.nB])))