We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
bug in pynumero example: sensitivity.py
It seems like the [sensitivity.py](https://github.com/Pyomo/pyomo/blob/master/pyomo/contrib/pynumero/examples/sensitivity.py) code is not calculating the sensitivity correctly. The sensitivity based x_new is a bit different from the result obtained by solving the perturbed NLP in the code. The cause of this seems to be the creation of `Np`

```
Np = BlockMatrix(2, 1)
Np[0, 0] = nlp.extract_submatrix_hessian_lag(pyomo_variables_rows=nlp.get_pyomo_variables(), pyomo_variables_cols=[m.eta1, m.eta2])
Np[1, 0] = nlp.extract_submatrix_jacobian(pyomo_variables=[m.eta1, m.eta2], pyomo_constraints=nlp.get_pyomo_constraints())
```
It seems that in this example Np should be 9x2 matrix with all zeros except elements (8, 1) and (9, 2) which should be one (based on eq 13 in this [reference paper](http://www.optimization-online.org/DB_FILE/2011/04/3008.pdf)). 

I have a [jupyter notebook](https://github.com/notesofdabbler/learn_pyomo/blob/master/wk_pynumero/wk_sensitivity.ipynb) that recreates this example and the proposed modification which seems to fix the issue.

--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `pyomo/contrib/pynumero/examples/sensitivity.py`
Content:
```
1 #  ___________________________________________________________________________
2 #
3 #  Pyomo: Python Optimization Modeling Objects
4 #  Copyright 2017 National Technology and Engineering Solutions of Sandia, LLC
5 #  Under the terms of Contract DE-NA0003525 with National Technology and
6 #  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
7 #  rights in this software.
8 #  This software is distributed under the 3-clause BSD License.
9 #  ___________________________________________________________________________
10 import pyomo.environ as pyo
11 from pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP
12 from pyomo.contrib.pynumero.sparse import BlockSymMatrix, BlockMatrix, BlockVector
13 from scipy.sparse import identity
14 from scipy.sparse.linalg import spsolve
15 import numpy as np
16 
17 
18 def create_model(eta1, eta2):
19     model = pyo.ConcreteModel()
20     # variables
21     model.x1 = pyo.Var(initialize=0.15)
22     model.x2 = pyo.Var(initialize=0.15)
23     model.x3 = pyo.Var(initialize=0.0)
24     # parameters
25     model.eta1 = pyo.Var()
26     model.eta2 = pyo.Var()
27 
28     model.nominal_eta1 = pyo.Param(initialize=eta1, mutable=True)
29     model.nominal_eta2 = pyo.Param(initialize=eta2, mutable=True)
30 
31     # constraints + objective
32     model.const1 = pyo.Constraint(expr=6*model.x1+3*model.x2+2*model.x3 - model.eta1 == 0)
33     model.const2 = pyo.Constraint(expr=model.eta2*model.x1+model.x2-model.x3-1 == 0)
34     model.cost = pyo.Objective(expr=model.x1**2 + model.x2**2 + model.x3**2)
35     model.consteta1 = pyo.Constraint(expr=model.eta1 == model.nominal_eta1)
36     model.consteta2 = pyo.Constraint(expr=model.eta2 == model.nominal_eta2)
37 
38     return model
39 
40 def compute_init_lam(nlp, x=None, lam_max=1e3):
41     if x is None:
42         x = nlp.init_primals()
43     else:
44         assert x.size == nlp.n_primals()
45     nlp.set_primals(x)
46 
47     assert nlp.n_ineq_constraints() == 0, "only supported for equality constrained nlps for now"
48 
49     nx = nlp.n_primals()
50     nc = nlp.n_constraints()
51 
52     # create Jacobian
53     jac = nlp.evaluate_jacobian()
54 
55     # create gradient of objective
56     df = nlp.evaluate_grad_objective()
57 
58     # create KKT system
59     kkt = BlockSymMatrix(2)
60     kkt[0, 0] = identity(nx)
61     kkt[1, 0] = jac
62 
63     zeros = np.zeros(nc)
64     rhs = BlockVector([-df, zeros])
65 
66     flat_kkt = kkt.tocoo().tocsc()
67     flat_rhs = rhs.flatten()
68 
69     sol = spsolve(flat_kkt, flat_rhs)
70     return sol[nlp.n_primals() : nlp.n_primals() + nlp.n_constraints()]
71 
72 #################################################################
73 m = create_model(4.5, 1.0)
74 opt = pyo.SolverFactory('ipopt')
75 results = opt.solve(m, tee=True)
76 
77 #################################################################
78 nlp = PyomoNLP(m)
79 x = nlp.init_primals()
80 y = compute_init_lam(nlp, x=x)
81 nlp.set_primals(x)
82 nlp.set_duals(y)
83 
84 J = nlp.evaluate_jacobian()
85 H = nlp.evaluate_hessian_lag()
86 
87 M = BlockSymMatrix(2)
88 M[0, 0] = H
89 M[1, 0] = J
90 
91 Np = BlockMatrix(2, 1)
92 Np[0, 0] = nlp.extract_submatrix_hessian_lag(pyomo_variables_rows=nlp.get_pyomo_variables(), pyomo_variables_cols=[m.eta1, m.eta2])
93 Np[1, 0] = nlp.extract_submatrix_jacobian(pyomo_variables=[m.eta1, m.eta2], pyomo_constraints=nlp.get_pyomo_constraints())
94 
95 ds = spsolve(M.tocsc(), Np.tocsc())
96 print(nlp.variable_names())
97 
98 #################################################################
99 
100 p0 = np.array([pyo.value(m.nominal_eta1), pyo.value(m.nominal_eta2)])
101 p = np.array([4.45, 1.05])
102 dp = p - p0
103 dx = ds.dot(dp)[0:nlp.n_primals()]
104 new_x = x + dx
105 print(new_x)
106 
107 #################################################################
108 m = create_model(4.45, 1.05)
109 opt = pyo.SolverFactory('ipopt')
110 results = opt.solve(m, tee=True)
111 nlp = PyomoNLP(m)
112 print(nlp.init_primals())
113 
114 
115
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/pyomo/contrib/pynumero/examples/sensitivity.py b/pyomo/contrib/pynumero/examples/sensitivity.py
--- a/pyomo/contrib/pynumero/examples/sensitivity.py
+++ b/pyomo/contrib/pynumero/examples/sensitivity.py
@@ -81,34 +81,39 @@
 nlp.set_primals(x)
 nlp.set_duals(y)
 
-J = nlp.evaluate_jacobian()
-H = nlp.evaluate_hessian_lag()
+J = nlp.extract_submatrix_jacobian(pyomo_variables=[m.x1, m.x2, m.x3], pyomo_constraints=[m.const1, m.const2])
+H = nlp.extract_submatrix_hessian_lag(pyomo_variables_rows=[m.x1, m.x2, m.x3], pyomo_variables_cols=[m.x1, m.x2, m.x3])
 
 M = BlockSymMatrix(2)
 M[0, 0] = H
 M[1, 0] = J
 
 Np = BlockMatrix(2, 1)
-Np[0, 0] = nlp.extract_submatrix_hessian_lag(pyomo_variables_rows=nlp.get_pyomo_variables(), pyomo_variables_cols=[m.eta1, m.eta2])
-Np[1, 0] = nlp.extract_submatrix_jacobian(pyomo_variables=[m.eta1, m.eta2], pyomo_constraints=nlp.get_pyomo_constraints())
+Np[0, 0] = nlp.extract_submatrix_hessian_lag(pyomo_variables_rows=[m.x1, m.x2, m.x3], pyomo_variables_cols=[m.eta1, m.eta2])
+Np[1, 0] = nlp.extract_submatrix_jacobian(pyomo_variables=[m.eta1, m.eta2], pyomo_constraints=[m.const1, m.const2])
 
-ds = spsolve(M.tocsc(), Np.tocsc())
-print(nlp.variable_names())
+ds = spsolve(M.tocsc(), -Np.tocsc())
 
+print("ds:\n", ds.todense())
 #################################################################
 
 p0 = np.array([pyo.value(m.nominal_eta1), pyo.value(m.nominal_eta2)])
 p = np.array([4.45, 1.05])
 dp = p - p0
-dx = ds.dot(dp)[0:nlp.n_primals()]
-new_x = x + dx
-print(new_x)
+dx = ds.dot(dp)[0:3]
+x_indices = nlp.get_primal_indices([m.x1, m.x2, m.x3])
+x_names = np.array(nlp.variable_names())
+new_x = x[x_indices] + dx
+print("dp:", dp)
+print("dx:", dx)
+print("Variable names: \n",x_names[x_indices])
+print("Sensitivity based x:\n", new_x)
 
 #################################################################
 m = create_model(4.45, 1.05)
 opt = pyo.SolverFactory('ipopt')
-results = opt.solve(m, tee=True)
+results = opt.solve(m, tee=False)
 nlp = PyomoNLP(m)
-print(nlp.init_primals())
-
+new_x = nlp.init_primals()
+print("NLP based x:\n", new_x[nlp.get_primal_indices([m.x1, m.x2, m.x3])])

{"golden_diff": "diff --git a/pyomo/contrib/pynumero/examples/sensitivity.py b/pyomo/contrib/pynumero/examples/sensitivity.py\n--- a/pyomo/contrib/pynumero/examples/sensitivity.py\n+++ b/pyomo/contrib/pynumero/examples/sensitivity.py\n@@ -81,34 +81,39 @@\n nlp.set_primals(x)\n nlp.set_duals(y)\n \n-J = nlp.evaluate_jacobian()\n-H = nlp.evaluate_hessian_lag()\n+J = nlp.extract_submatrix_jacobian(pyomo_variables=[m.x1, m.x2, m.x3], pyomo_constraints=[m.const1, m.const2])\n+H = nlp.extract_submatrix_hessian_lag(pyomo_variables_rows=[m.x1, m.x2, m.x3], pyomo_variables_cols=[m.x1, m.x2, m.x3])\n \n M = BlockSymMatrix(2)\n M[0, 0] = H\n M[1, 0] = J\n \n Np = BlockMatrix(2, 1)\n-Np[0, 0] = nlp.extract_submatrix_hessian_lag(pyomo_variables_rows=nlp.get_pyomo_variables(), pyomo_variables_cols=[m.eta1, m.eta2])\n-Np[1, 0] = nlp.extract_submatrix_jacobian(pyomo_variables=[m.eta1, m.eta2], pyomo_constraints=nlp.get_pyomo_constraints())\n+Np[0, 0] = nlp.extract_submatrix_hessian_lag(pyomo_variables_rows=[m.x1, m.x2, m.x3], pyomo_variables_cols=[m.eta1, m.eta2])\n+Np[1, 0] = nlp.extract_submatrix_jacobian(pyomo_variables=[m.eta1, m.eta2], pyomo_constraints=[m.const1, m.const2])\n \n-ds = spsolve(M.tocsc(), Np.tocsc())\n-print(nlp.variable_names())\n+ds = spsolve(M.tocsc(), -Np.tocsc())\n \n+print(\"ds:\\n\", ds.todense())\n #################################################################\n \n p0 = np.array([pyo.value(m.nominal_eta1), pyo.value(m.nominal_eta2)])\n p = np.array([4.45, 1.05])\n dp = p - p0\n-dx = ds.dot(dp)[0:nlp.n_primals()]\n-new_x = x + dx\n-print(new_x)\n+dx = ds.dot(dp)[0:3]\n+x_indices = nlp.get_primal_indices([m.x1, m.x2, m.x3])\n+x_names = np.array(nlp.variable_names())\n+new_x = x[x_indices] + dx\n+print(\"dp:\", dp)\n+print(\"dx:\", dx)\n+print(\"Variable names: \\n\",x_names[x_indices])\n+print(\"Sensitivity based x:\\n\", new_x)\n \n #################################################################\n m = create_model(4.45, 1.05)\n opt = pyo.SolverFactory('ipopt')\n-results = opt.solve(m, tee=True)\n+results = opt.solve(m, tee=False)\n nlp = PyomoNLP(m)\n-print(nlp.init_primals())\n-\n+new_x = nlp.init_primals()\n+print(\"NLP based x:\\n\", new_x[nlp.get_primal_indices([m.x1, m.x2, m.x3])])\n", "issue": "bug in pynumero example: sensitivity.py\nIt seems like the [sensitivity.py](https://github.com/Pyomo/pyomo/blob/master/pyomo/contrib/pynumero/examples/sensitivity.py) code is not calculating the sensitivity correctly. The sensitivity based x_new is a bit different from the result obtained by solving the perturbed NLP in the code. The cause of this seems to be the creation of `Np`\r\n\r\n```\r\nNp = BlockMatrix(2, 1)\r\nNp[0, 0] = nlp.extract_submatrix_hessian_lag(pyomo_variables_rows=nlp.get_pyomo_variables(), pyomo_variables_cols=[m.eta1, m.eta2])\r\nNp[1, 0] = nlp.extract_submatrix_jacobian(pyomo_variables=[m.eta1, m.eta2], pyomo_constraints=nlp.get_pyomo_constraints())\r\n```\r\nIt seems that in this example Np should be 9x2 matrix with all zeros except elements (8, 1) and (9, 2) which should be one (based on eq 13 in this [reference paper](http://www.optimization-online.org/DB_FILE/2011/04/3008.pdf)). \r\n\r\nI have a [jupyter notebook](https://github.com/notesofdabbler/learn_pyomo/blob/master/wk_pynumero/wk_sensitivity.ipynb) that recreates this example and the proposed modification which seems to fix the issue.\n", "before_files": [{"content": "#  ___________________________________________________________________________\n#\n#  Pyomo: Python Optimization Modeling Objects\n#  Copyright 2017 National Technology and Engineering Solutions of Sandia, LLC\n#  Under the terms of Contract DE-NA0003525 with National Technology and\n#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain\n#  rights in this software.\n#  This software is distributed under the 3-clause BSD License.\n#  ___________________________________________________________________________\nimport pyomo.environ as pyo\nfrom pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP\nfrom pyomo.contrib.pynumero.sparse import BlockSymMatrix, BlockMatrix, BlockVector\nfrom scipy.sparse import identity\nfrom scipy.sparse.linalg import spsolve\nimport numpy as np\n\n\ndef create_model(eta1, eta2):\n    model = pyo.ConcreteModel()\n    # variables\n    model.x1 = pyo.Var(initialize=0.15)\n    model.x2 = pyo.Var(initialize=0.15)\n    model.x3 = pyo.Var(initialize=0.0)\n    # parameters\n    model.eta1 = pyo.Var()\n    model.eta2 = pyo.Var()\n\n    model.nominal_eta1 = pyo.Param(initialize=eta1, mutable=True)\n    model.nominal_eta2 = pyo.Param(initialize=eta2, mutable=True)\n\n    # constraints + objective\n    model.const1 = pyo.Constraint(expr=6*model.x1+3*model.x2+2*model.x3 - model.eta1 == 0)\n    model.const2 = pyo.Constraint(expr=model.eta2*model.x1+model.x2-model.x3-1 == 0)\n    model.cost = pyo.Objective(expr=model.x1**2 + model.x2**2 + model.x3**2)\n    model.consteta1 = pyo.Constraint(expr=model.eta1 == model.nominal_eta1)\n    model.consteta2 = pyo.Constraint(expr=model.eta2 == model.nominal_eta2)\n\n    return model\n\ndef compute_init_lam(nlp, x=None, lam_max=1e3):\n    if x is None:\n        x = nlp.init_primals()\n    else:\n        assert x.size == nlp.n_primals()\n    nlp.set_primals(x)\n\n    assert nlp.n_ineq_constraints() == 0, \"only supported for equality constrained nlps for now\"\n\n    nx = nlp.n_primals()\n    nc = nlp.n_constraints()\n\n    # create Jacobian\n    jac = nlp.evaluate_jacobian()\n\n    # create gradient of objective\n    df = nlp.evaluate_grad_objective()\n\n    # create KKT system\n    kkt = BlockSymMatrix(2)\n    kkt[0, 0] = identity(nx)\n    kkt[1, 0] = jac\n\n    zeros = np.zeros(nc)\n    rhs = BlockVector([-df, zeros])\n\n    flat_kkt = kkt.tocoo().tocsc()\n    flat_rhs = rhs.flatten()\n\n    sol = spsolve(flat_kkt, flat_rhs)\n    return sol[nlp.n_primals() : nlp.n_primals() + nlp.n_constraints()]\n\n#################################################################\nm = create_model(4.5, 1.0)\nopt = pyo.SolverFactory('ipopt')\nresults = opt.solve(m, tee=True)\n\n#################################################################\nnlp = PyomoNLP(m)\nx = nlp.init_primals()\ny = compute_init_lam(nlp, x=x)\nnlp.set_primals(x)\nnlp.set_duals(y)\n\nJ = nlp.evaluate_jacobian()\nH = nlp.evaluate_hessian_lag()\n\nM = BlockSymMatrix(2)\nM[0, 0] = H\nM[1, 0] = J\n\nNp = BlockMatrix(2, 1)\nNp[0, 0] = nlp.extract_submatrix_hessian_lag(pyomo_variables_rows=nlp.get_pyomo_variables(), pyomo_variables_cols=[m.eta1, m.eta2])\nNp[1, 0] = nlp.extract_submatrix_jacobian(pyomo_variables=[m.eta1, m.eta2], pyomo_constraints=nlp.get_pyomo_constraints())\n\nds = spsolve(M.tocsc(), Np.tocsc())\nprint(nlp.variable_names())\n\n#################################################################\n\np0 = np.array([pyo.value(m.nominal_eta1), pyo.value(m.nominal_eta2)])\np = np.array([4.45, 1.05])\ndp = p - p0\ndx = ds.dot(dp)[0:nlp.n_primals()]\nnew_x = x + dx\nprint(new_x)\n\n#################################################################\nm = create_model(4.45, 1.05)\nopt = pyo.SolverFactory('ipopt')\nresults = opt.solve(m, tee=True)\nnlp = PyomoNLP(m)\nprint(nlp.init_primals())\n\n\n", "path": "pyomo/contrib/pynumero/examples/sensitivity.py"}], "after_files": [{"content": "#  ___________________________________________________________________________\n#\n#  Pyomo: Python Optimization Modeling Objects\n#  Copyright 2017 National Technology and Engineering Solutions of Sandia, LLC\n#  Under the terms of Contract DE-NA0003525 with National Technology and\n#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain\n#  rights in this software.\n#  This software is distributed under the 3-clause BSD License.\n#  ___________________________________________________________________________\nimport pyomo.environ as pyo\nfrom pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP\nfrom pyomo.contrib.pynumero.sparse import BlockSymMatrix, BlockMatrix, BlockVector\nfrom scipy.sparse import identity\nfrom scipy.sparse.linalg import spsolve\nimport numpy as np\n\n\ndef create_model(eta1, eta2):\n    model = pyo.ConcreteModel()\n    # variables\n    model.x1 = pyo.Var(initialize=0.15)\n    model.x2 = pyo.Var(initialize=0.15)\n    model.x3 = pyo.Var(initialize=0.0)\n    # parameters\n    model.eta1 = pyo.Var()\n    model.eta2 = pyo.Var()\n\n    model.nominal_eta1 = pyo.Param(initialize=eta1, mutable=True)\n    model.nominal_eta2 = pyo.Param(initialize=eta2, mutable=True)\n\n    # constraints + objective\n    model.const1 = pyo.Constraint(expr=6*model.x1+3*model.x2+2*model.x3 - model.eta1 == 0)\n    model.const2 = pyo.Constraint(expr=model.eta2*model.x1+model.x2-model.x3-1 == 0)\n    model.cost = pyo.Objective(expr=model.x1**2 + model.x2**2 + model.x3**2)\n    model.consteta1 = pyo.Constraint(expr=model.eta1 == model.nominal_eta1)\n    model.consteta2 = pyo.Constraint(expr=model.eta2 == model.nominal_eta2)\n\n    return model\n\ndef compute_init_lam(nlp, x=None, lam_max=1e3):\n    if x is None:\n        x = nlp.init_primals()\n    else:\n        assert x.size == nlp.n_primals()\n    nlp.set_primals(x)\n\n    assert nlp.n_ineq_constraints() == 0, \"only supported for equality constrained nlps for now\"\n\n    nx = nlp.n_primals()\n    nc = nlp.n_constraints()\n\n    # create Jacobian\n    jac = nlp.evaluate_jacobian()\n\n    # create gradient of objective\n    df = nlp.evaluate_grad_objective()\n\n    # create KKT system\n    kkt = BlockSymMatrix(2)\n    kkt[0, 0] = identity(nx)\n    kkt[1, 0] = jac\n\n    zeros = np.zeros(nc)\n    rhs = BlockVector([-df, zeros])\n\n    flat_kkt = kkt.tocoo().tocsc()\n    flat_rhs = rhs.flatten()\n\n    sol = spsolve(flat_kkt, flat_rhs)\n    return sol[nlp.n_primals() : nlp.n_primals() + nlp.n_constraints()]\n\n#################################################################\nm = create_model(4.5, 1.0)\nopt = pyo.SolverFactory('ipopt')\nresults = opt.solve(m, tee=True)\n\n#################################################################\nnlp = PyomoNLP(m)\nx = nlp.init_primals()\ny = compute_init_lam(nlp, x=x)\nnlp.set_primals(x)\nnlp.set_duals(y)\n\nJ = nlp.extract_submatrix_jacobian(pyomo_variables=[m.x1, m.x2, m.x3], pyomo_constraints=[m.const1, m.const2])\nH = nlp.extract_submatrix_hessian_lag(pyomo_variables_rows=[m.x1, m.x2, m.x3], pyomo_variables_cols=[m.x1, m.x2, m.x3])\n\nM = BlockSymMatrix(2)\nM[0, 0] = H\nM[1, 0] = J\n\nNp = BlockMatrix(2, 1)\nNp[0, 0] = nlp.extract_submatrix_hessian_lag(pyomo_variables_rows=[m.x1, m.x2, m.x3], pyomo_variables_cols=[m.eta1, m.eta2])\nNp[1, 0] = nlp.extract_submatrix_jacobian(pyomo_variables=[m.eta1, m.eta2], pyomo_constraints=[m.const1, m.const2])\n\nds = spsolve(M.tocsc(), -Np.tocsc())\n\nprint(\"ds:\\n\", ds.todense())\n#################################################################\n\np0 = np.array([pyo.value(m.nominal_eta1), pyo.value(m.nominal_eta2)])\np = np.array([4.45, 1.05])\ndp = p - p0\ndx = ds.dot(dp)[0:3]\nx_indices = nlp.get_primal_indices([m.x1, m.x2, m.x3])\nx_names = np.array(nlp.variable_names())\nnew_x = x[x_indices] + dx\nprint(\"dp:\", dp)\nprint(\"dx:\", dx)\nprint(\"Variable names: \\n\",x_names[x_indices])\nprint(\"Sensitivity based x:\\n\", new_x)\n\n#################################################################\nm = create_model(4.45, 1.05)\nopt = pyo.SolverFactory('ipopt')\nresults = opt.solve(m, tee=False)\nnlp = PyomoNLP(m)\nnew_x = nlp.init_primals()\nprint(\"NLP based x:\\n\", new_x[nlp.get_primal_indices([m.x1, m.x2, m.x3])])\n\n", "path": "pyomo/contrib/pynumero/examples/sensitivity.py"}]}