prompt: 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()
  
  ```
  
golden_diff: 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])])
  
verification_info: {"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"}]}