Merge pull request #3655 from eslickj/pynum_sens

Add independent variable rows to sensitivity toolbox Pynumero-based sensitivity calcs

Author: blnicho

pyomo/contrib/sensitivity_toolbox/pynumero.py

@@ -62,7 +62,7 @@ def get_dsdp_dfdp(model, theta):
     Returns
     -------
     scipy.sparse.csc_matrix, csc_matrix, ComponentMap, ComponentMap
-        ds/dp (ns by np), df/dp (1 by np), row map, column map.
+        ds/dp (ns + np by np), df/dp (1 by np), row map, column map.
         The column map maps Pyomo variables p to columns and the
         row map maps Pyomo variables s to rows.
     """
@@ -98,6 +98,15 @@ def get_dsdp_dfdp(model, theta):
     _coo_reorder_cols(dfdx, remap=col_remap)
     dfdx = dfdx.tocsc()
     dfdp = dfdx[0, :ns] @ dsdp + dfdx[0, ns:]
+
+    # Add independent variables to the rows as the end
+    dsdp = scipy.sparse.vstack(
+        [dsdp, scipy.sparse.identity(len(column_map), format="csc")]
+    )
+    n = len(row_map)
+    for p, i in column_map.items():
+        row_map[p] = i + n
+
     # return sensitivity of the outputs to p and component maps
     return dsdp, dfdp, row_map, column_map
 

pyomo/contrib/sensitivity_toolbox/tests/test_pynumero.py

@@ -73,11 +73,15 @@ def test_dsdp_dfdp_pyomo_nlp(self):
         dsdp, dfdp, rmap, cmap = pnsens.get_dsdp_dfdp(m2, theta)
 
         # Since x1 = p1
-        assert dsdp.shape == (2, 2)
+        assert dsdp.shape == (4, 2)
         np.testing.assert_almost_equal(dsdp[rmap[m.x1], cmap[m.p1]], 40.0)
         np.testing.assert_almost_equal(dsdp[rmap[m.x1], cmap[m.p2]], 0.0)
         np.testing.assert_almost_equal(dsdp[rmap[m.x2], cmap[m.p1]], 0.0)
         np.testing.assert_almost_equal(dsdp[rmap[m.x2], cmap[m.p2]], 1.0)
+        np.testing.assert_almost_equal(dsdp[rmap[m.p1], cmap[m.p2]], 0.0)
+        np.testing.assert_almost_equal(dsdp[rmap[m.p2], cmap[m.p2]], 1.0)
+        np.testing.assert_almost_equal(dsdp[rmap[m.p1], cmap[m.p1]], 1.0)
+        np.testing.assert_almost_equal(dsdp[rmap[m.p2], cmap[m.p1]], 0.0)
 
         # if x1 = 2 * p1 and x2 = p2 then
         #   df/dp1 = 6 p1**2 + p2 = 45.0