Merge branch 'develop'

README.rst

@@ -77,7 +77,7 @@ The ESA team works hard to ensure ESA is well tested, and we strive for
 testing coverage data for ESA, using `coverage
 <https://pypi.org/project/coverage/>`__.
 
-.. table:: ESA's testing coverage as of 2022-05-20 (Git commit: bbb3f12)
+.. table:: ESA's testing coverage as of 2022-08-22 (Git commit: 4b28569)
     :widths: auto
     :align: left
 
@@ -86,7 +86,7 @@ testing coverage data for ESA, using `coverage
     +=================+===================+=================+=================+====================+
     | esa/__init__.py |                 2 |               0 |               2 |                100 |
     +-----------------+-------------------+-----------------+-----------------+--------------------+
-    | esa/saw.py      |               848 |               0 |             848 |                100 |
+    | esa/saw.py      |               920 |               0 |             920 |                100 |
     +-----------------+-------------------+-----------------+-----------------+--------------------+
 
 License

VERSION

@@ -1 +1 @@
-1.2.7
+1.2.8

docs/html/.buildinfo

@@ -1,4 +1,4 @@
 # Sphinx build info version 1
 # This file hashes the configuration used when building these files. When it is not found, a full rebuild will be done.
-config: ed1a8386bdf62a50cde5a8b4689068ae
+config: a53f9ecfe3e60fd7052d70b826a9f37e
 tags: 645f666f9bcd5a90fca523b33c5a78b7

docs/html/_static/documentation_options.js

@@ -1,6 +1,6 @@
 var DOCUMENTATION_OPTIONS = {
     URL_ROOT: document.getElementById("documentation_options").getAttribute('data-url_root'),
-    VERSION: '1.2.7',
+    VERSION: '1.2.8',
     LANGUAGE: 'None',
     COLLAPSE_INDEX: false,
     BUILDER: 'html',

docs/html/esa.html

@@ -1478,6 +1478,18 @@
 </dl>
 </dd></dl>
 
+<dl class="py method">
+<dt class="sig sig-object py" id="esa.saw.SAW.get_lodf_matrix_fast">
+<span class="sig-name descname"><span class="pre">get_lodf_matrix_fast</span></span><span class="sig-paren">(</span><span class="sig-paren">)</span><a class="headerlink" href="#esa.saw.SAW.get_lodf_matrix_fast" title="Permalink to this definition">¶</a></dt>
+<dd><p>Calculate the line outage distribution factor natively. This method should be much
+faster than the PW script command for large cases.</p>
+<dl class="field-list simple">
+<dt class="field-odd">Returns</dt>
+<dd class="field-odd"><p>A dense float matrix in the numpy array format.</p>
+</dd>
+</dl>
+</dd></dl>
+
 <dl class="py method">
 <dt class="sig sig-object py" id="esa.saw.SAW.get_power_flow_results">
 <span class="sig-name descname"><span class="pre">get_power_flow_results</span></span><span class="sig-paren">(</span><em class="sig-param"><span class="n"><span class="pre">ObjectType</span></span><span class="p"><span class="pre">:</span></span> <span class="n"><span class="pre">str</span></span></em>, <em class="sig-param"><span class="n"><span class="pre">additional_fields</span></span><span class="p"><span class="pre">:</span></span> <span class="n"><span class="pre">Union</span><span class="p"><span class="pre">[</span></span><span class="pre">None</span><span class="p"><span class="pre">,</span> </span><span class="pre">List</span><span class="p"><span class="pre">[</span></span><span class="pre">str</span><span class="p"><span class="pre">]</span></span><span class="p"><span class="pre">]</span></span></span> <span class="o"><span class="pre">=</span></span> <span class="default_value"><span class="pre">None</span></span></em><span class="sig-paren">)</span> <span class="sig-return"><span class="sig-return-icon">&#x2192;</span> <span class="sig-return-typehint"><span class="pre">Union</span><span class="p"><span class="pre">[</span></span><span class="pre">None</span><span class="p"><span class="pre">,</span> </span><span class="pre">pandas.core.frame.DataFrame</span><span class="p"><span class="pre">]</span></span></span></span><a class="headerlink" href="#esa.saw.SAW.get_power_flow_results" title="Permalink to this definition">¶</a></dt>
@@ -1502,6 +1514,45 @@
 </dl>
 </dd></dl>
 
+<dl class="py method">
+<dt class="sig sig-object py" id="esa.saw.SAW.get_ptdf_matrix_fast">
+<span class="sig-name descname"><span class="pre">get_ptdf_matrix_fast</span></span><span class="sig-paren">(</span><span class="sig-paren">)</span><a class="headerlink" href="#esa.saw.SAW.get_ptdf_matrix_fast" title="Permalink to this definition">¶</a></dt>
+<dd><p>Calculate the power transfer distribution factor natively. This method should be much
+faster than the PW script command for large cases.</p>
+<dl class="field-list simple">
+<dt class="field-odd">Returns</dt>
+<dd class="field-odd"><p>A dense float matrix in the numpy array format.</p>
+</dd>
+</dl>
+</dd></dl>
+
+<dl class="py method">
+<dt class="sig sig-object py" id="esa.saw.SAW.get_shift_factor_matrix">
+<span class="sig-name descname"><span class="pre">get_shift_factor_matrix</span></span><span class="sig-paren">(</span><em class="sig-param"><span class="n"><span class="pre">method</span></span><span class="p"><span class="pre">:</span></span> <span class="n"><span class="pre">str</span></span> <span class="o"><span class="pre">=</span></span> <span class="default_value"><span class="pre">'DC'</span></span></em><span class="sig-paren">)</span><a class="headerlink" href="#esa.saw.SAW.get_shift_factor_matrix" title="Permalink to this definition">¶</a></dt>
+<dd><p>Calculate the injection shift factor matrix using the auxiliary
+script CalculateShiftFactorsMultipleElement.</p>
+<dl class="field-list simple">
+<dt class="field-odd">Parameters</dt>
+<dd class="field-odd"><p><strong>method</strong> – The linear method to be used for the calculation. The</p>
+</dd>
+</dl>
+<p>options are AC, DC or DCPS.
+:returns: A dense float matrix in the numpy array format.</p>
+</dd></dl>
+
+<dl class="py method">
+<dt class="sig sig-object py" id="esa.saw.SAW.get_shift_factor_matrix_fast">
+<span class="sig-name descname"><span class="pre">get_shift_factor_matrix_fast</span></span><span class="sig-paren">(</span><span class="sig-paren">)</span><a class="headerlink" href="#esa.saw.SAW.get_shift_factor_matrix_fast" title="Permalink to this definition">¶</a></dt>
+<dd><p>Calculate the injection shift factor matrix directly using the incidence
+matrix and the susceptance matrix. This method should be much faster than
+the PW script command for large cases.</p>
+<dl class="field-list simple">
+<dt class="field-odd">Returns</dt>
+<dd class="field-odd"><p>A dense float matrix in the numpy array format.</p>
+</dd>
+</dl>
+</dd></dl>
+
 <dl class="py method">
 <dt class="sig sig-object py" id="esa.saw.SAW.get_shunt_admittance">
 <span class="sig-name descname"><span class="pre">get_shunt_admittance</span></span><span class="sig-paren">(</span><span class="sig-paren">)</span><a class="headerlink" href="#esa.saw.SAW.get_shunt_admittance" title="Permalink to this definition">¶</a></dt>

docs/html/genindex.html

@@ -320,15 +320,23 @@ <h2 id="G">G</h2>
       <li><a href="esa.html#esa.saw.SAW.get_key_fields_for_object_type">get_key_fields_for_object_type() (esa.saw.SAW method)</a>
 </li>
       <li><a href="esa.html#esa.saw.SAW.get_lodf_matrix">get_lodf_matrix() (esa.saw.SAW method)</a>
+</li>
+      <li><a href="esa.html#esa.saw.SAW.get_lodf_matrix_fast">get_lodf_matrix_fast() (esa.saw.SAW method)</a>
 </li>
       <li><a href="esa.html#esa.saw.SAW.get_power_flow_results">get_power_flow_results() (esa.saw.SAW method)</a>
 </li>
-      <li><a href="esa.html#esa.saw.SAW.get_shunt_admittance">get_shunt_admittance() (esa.saw.SAW method)</a>
+      <li><a href="esa.html#esa.saw.SAW.get_ptdf_matrix_fast">get_ptdf_matrix_fast() (esa.saw.SAW method)</a>
 </li>
-      <li><a href="esa.html#esa.saw.SAW.get_version_and_builddate">get_version_and_builddate() (esa.saw.SAW method)</a>
+      <li><a href="esa.html#esa.saw.SAW.get_shift_factor_matrix">get_shift_factor_matrix() (esa.saw.SAW method)</a>
+</li>
+      <li><a href="esa.html#esa.saw.SAW.get_shift_factor_matrix_fast">get_shift_factor_matrix_fast() (esa.saw.SAW method)</a>
 </li>
   </ul></td>
   <td style="width: 33%; vertical-align: top;"><ul>
+      <li><a href="esa.html#esa.saw.SAW.get_shunt_admittance">get_shunt_admittance() (esa.saw.SAW method)</a>
+</li>
+      <li><a href="esa.html#esa.saw.SAW.get_version_and_builddate">get_version_and_builddate() (esa.saw.SAW method)</a>
+</li>
       <li><a href="esa.html#esa.saw.SAW.get_ybus">get_ybus() (esa.saw.SAW method)</a>
 </li>
       <li><a href="esa.html#esa.saw.SAW.GetCaseHeader">GetCaseHeader() (esa.saw.SAW method)</a>

docs/html/overview.html

@@ -316,7 +316,7 @@ <h2>Testing Coverage<a class="headerlink" href="#testing-coverage" title="Permal
 <p>The ESA team strives to write good tests with 100% coverage. The table
 below provides the latest test coverage data for ESA.</p>
 <table class="colwidths-auto docutils align-left" id="id3">
-<caption><span class="caption-text">ESA’s testing coverage as of 2022-05-20 (Git commit: bbb3f12)</span><a class="headerlink" href="#id3" title="Permalink to this table">¶</a></caption>
+<caption><span class="caption-text">ESA’s testing coverage as of 2022-08-22 (Git commit: 4b28569)</span><a class="headerlink" href="#id3" title="Permalink to this table">¶</a></caption>
 <thead>
 <tr class="row-odd"><th class="head"><p>Name</p></th>
 <th class="head"><p>Num. Statements</p></th>
@@ -333,9 +333,9 @@ <h2>Testing Coverage<a class="headerlink" href="#testing-coverage" title="Permal
 <td><p>100</p></td>
 </tr>
 <tr class="row-odd"><td><p>esa/saw.py</p></td>
-<td><p>848</p></td>
+<td><p>920</p></td>
 <td><p>0</p></td>
-<td><p>848</p></td>
+<td><p>920</p></td>
 <td><p>100</p></td>
 </tr>
 </tbody>

docs/html/welcome.html

@@ -544,6 +544,12 @@ <h3>Credits and Thanks<a class="headerlink" href="#credits-and-thanks" title="Pe
 <p>Changes made with each ESA release are listed here. Please note that
 versions prior to 1.0.0 are not listed here, but are still available on
 <a class="reference external" href="https://pypi.org/project/esa/#history">PyPi</a>.</p>
+<section id="version-1-2-8">
+<h3>Version 1.2.8<a class="headerlink" href="#version-1-2-8" title="Permalink to this headline">¶</a></h3>
+<ul class="simple">
+<li><p>Native support on sensitivity analysis (ISF, PTDF, LODF) for better performance</p></li>
+</ul>
+</section>
 <section id="version-1-2-7">
 <h3>Version 1.2.7<a class="headerlink" href="#version-1-2-7" title="Permalink to this headline">¶</a></h3>
 <ul class="simple">

docs/rst/coverage.rst

@@ -1,4 +1,4 @@
-.. table:: ESA's testing coverage as of 2022-05-20 (Git commit: bbb3f12)
+.. table:: ESA's testing coverage as of 2022-08-22 (Git commit: 4b28569)
     :widths: auto
     :align: left
 
@@ -7,5 +7,5 @@
     +=================+===================+=================+=================+====================+
     | esa/__init__.py |                 2 |               0 |               2 |                100 |
     +-----------------+-------------------+-----------------+-----------------+--------------------+
-    | esa/saw.py      |               848 |               0 |             848 |                100 |
+    | esa/saw.py      |               920 |               0 |             920 |                100 |
     +-----------------+-------------------+-----------------+-----------------+--------------------+

docs/rst/welcome/changelog.rst

@@ -2,6 +2,11 @@ Changes made with each ESA release are listed here. Please note that
 versions prior to 1.0.0 are not listed here, but are still available on
 `PyPi <https://pypi.org/project/esa/#history>`__.
 
+Version 1.2.8
+^^^^^^^^^^^^^
+
+* Native support on sensitivity analysis (ISF, PTDF, LODF) for better performance
+
 Version 1.2.7
 ^^^^^^^^^^^^^
 

esa/__init__.py

@@ -16,4 +16,4 @@
 # Please keep the docstring above up to date with all the imports.
 from .saw import SAW, PowerWorldError, COMError, CommandNotRespectedError,\
     Error
-__version__ = "1.2.7"
+__version__ = "1.2.8"

esa/saw.py

@@ -24,6 +24,7 @@
 import numba as nb
 import pandas as pd
 from scipy.sparse import csr_matrix, coo_matrix, hstack
+from scipy import sparse as sp
 import networkx as nx
 from tqdm import trange
 import pythoncom
@@ -925,6 +926,128 @@ def get_incidence_matrix(self):
             incidence[i, row["BusNum:1"] - 1] = -1
         return incidence
 
+    def get_shift_factor_matrix(self, method: str = 'DC'):
+        """
+        Calculate the injection shift factor matrix using the auxiliary
+        script CalculateShiftFactorsMultipleElement.
+
+        :param method: The linear method to be used for the calculation. The
+        options are AC, DC or DCPS.
+        :returns: A dense float matrix in the numpy array format.
+        """
+        temp = self.pw_order
+        self.pw_order = True
+        key = self.get_key_field_list('branch')
+        fields = key + ['Selected']
+        df = self.GetParametersMultipleElement('branch', fields)
+        num_branch = df.shape[0]
+        df['Selected'] = 'YES'
+        self.change_parameters_multiple_element_df('branch', df)
+        # now run the calculation for all the seleced branches
+        self.RunScriptCommand(f"CalculateShiftFactorsMultipleElement(BRANCH,SELECTED,BUYER,"
+                              f"[SLACK],{method})")
+        isf_fields = ['MultBusTLRSens']
+        for i in range(1, num_branch):
+            isf_fields += [f"MultBusTLRSens:{i}"]
+        res = self.GetParametersMultipleElement('Bus', isf_fields).to_numpy(dtype=float)
+        self.pw_order = temp
+        return res
+
+
+    def _prepare_sensitivity(self):
+        """
+        Prepare the matrix for sensitivity analysis.
+        """
+        temp = self.pw_order
+        self.pw_order = True
+        bus = self.GetParametersMultipleElement('bus', ['BusNum', 'BusCat'])
+        br = self.GetParametersMultipleElement('branch', ['BusNum', 'BusNum:1', 'LineX'])
+        self.pw_order = temp
+        slack = bus[bus['BusCat'] == 'Slack'].index.tolist()[0]
+        noslack = bus.index.tolist()
+        noslack.remove(slack)
+        bus.reset_index(inplace=True)
+        bus.set_index('BusNum', inplace=True)
+        f = br['BusNum'].map(bus['index']).to_numpy(dtype=int)
+        t = br['BusNum:1'].map(bus['index']).to_numpy(dtype=int)
+        nl = br.shape[0]
+        nb = bus.shape[0]
+        i = np.r_[range(nl), range(nl)]
+        Cft = csr_matrix((np.r_[np.ones(nl), -np.ones(nl)], (i, np.r_[f, t])), (nl, nb))
+        x_val = br['LineX'].to_numpy(dtype=float)
+        b = 1 / x_val
+        Bf = csr_matrix((np.r_[b, -b], (i, np.r_[f, t])))
+        Bbus = Cft.T * Bf
+        Bbus[0, :] = 0
+        Bbus[:, 0] = 0
+        Bbus[0, 0] = -1
+        return Bbus, Bf, Cft, slack, noslack
+
+
+    def get_shift_factor_matrix_fast(self):
+        """
+        Calculate the injection shift factor matrix directly using the incidence
+        matrix and the susceptance matrix. This method should be much faster than
+        the PW script command for large cases.
+
+        :returns: A dense float matrix in the numpy array format.
+        """
+        Bbus, Bf, _, _, _ = self._prepare_sensitivity()
+        res = Bf * sp.linalg.inv(Bbus)
+        res[:, 0] = 0
+        return res.T.todense()
+
+
+    def get_ptdf_matrix_fast(self):
+        """
+        Calculate the power transfer distribution factor natively. This method should be much
+        faster than the PW script command for large cases.
+
+        :returns: A dense float matrix in the numpy array format.
+        """
+        Bbus, Bf, _, slack, noslack = self._prepare_sensitivity()
+        n = Bbus.shape[0]
+        noref = noslack
+        dP = np.eye(n, n)
+
+        # solve for change in voltage angles
+        dTheta = np.zeros((n, n))
+        Bref = Bbus[noslack, :][:, noref].tocsc()
+        dtheta_ref = sp.linalg.spsolve(Bref, dP[noslack, :])
+
+        dTheta[noref, :] = dtheta_ref
+
+        return Bf * dTheta
+
+
+    def get_lodf_matrix_fast(self):
+        """
+        Calculate the line outage distribution factor natively. This method should be much
+        faster than the PW script command for large cases.
+
+        :returns: A dense float matrix in the numpy array format.
+        """
+        Bbus, Bf, Cft, slack, noslack = self._prepare_sensitivity()
+        PTDF = self.get_ptdf_matrix_fast()
+        H = PTDF * Cft.T
+        numerical_zero = 1e-10
+
+        nl, nb = PTDF.shape
+        # this loop avoids the divisions by zero
+        # in those cases the LODF column should be zero
+        LODF = np.zeros((nl, nl))
+        div = 1 - H.diagonal()
+        islands = []
+        for j in range(H.shape[1]):
+            if abs(div[j]) > numerical_zero:
+                LODF[:, j] = H[:, j] / div[j]
+            else:
+                islands.append(j)
+
+        res = LODF - np.diag(np.diag(LODF)) - np.eye(nl, nl)
+        return res.T
+
+
     def change_to_temperature(self, T: Union[int, float, np.ndarray], R25=7.283, R75=8.688):
         """
         Change line resistance according to temperature.

setup.py

@@ -66,7 +66,7 @@
     # of sphinx for documentation. Coverage is necessary to keep the
     # coverage report up to date.
     extras_require={'test': ['networkx', 'coverage', 'matplotlib'],
-                    'doc': ['sphinx', 'tabulate'],
+                    'doc': ['sphinx', 'tabulate', 'sphinx_press_theme'],
                     'dev': ['pythran']},
     license='MIT',
     # TODO: Why aren't we zip safe?