Skip to content

Commit a1410b8

Browse files
quantresearch1quantresearch1
authored
Merge branch 'main' into feature/3381_quadratic_obj_highs
2 parents a5317fa + 926c91a commit a1410b8

File tree

10 files changed

+625
-58
lines changed

10 files changed

+625
-58
lines changed

‎pyomo/contrib/parmest/examples/semibatch/semibatch.py‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -284,6 +284,20 @@ def label_model(self):
284284

285285
m = self.model
286286

287+
m.experiment_outputs = Suffix(direction=Suffix.LOCAL)
288+
m.experiment_outputs.update(
289+
(m.Ca[t], self.data["Ca_meas"][f"{t}"]) for t in m.measT
290+
)
291+
m.experiment_outputs.update(
292+
(m.Cb[t], self.data["Cb_meas"][f"{t}"]) for t in m.measT
293+
)
294+
m.experiment_outputs.update(
295+
(m.Cc[t], self.data["Cc_meas"][f"{t}"]) for t in m.measT
296+
)
297+
m.experiment_outputs.update(
298+
(m.Tr[t], self.data["Tr_meas"][f"{t}"]) for t in m.measT
299+
)
300+
287301
m.unknown_parameters = Suffix(direction=Suffix.LOCAL)
288302
m.unknown_parameters.update(
289303
(k, ComponentUID(k)) for k in [m.k1, m.k2, m.E1, m.E2]

‎pyomo/contrib/parmest/parmest.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -284,13 +284,13 @@ def __init__(
284284
try:
285285
outputs = [k.name for k, v in model.experiment_outputs.items()]
286286
except:
287-
RuntimeError(
287+
raise RuntimeError(
288288
'Experiment list model does not have suffix ' + '"experiment_outputs".'
289289
)
290290
try:
291291
params = [k.name for k, v in model.unknown_parameters.items()]
292292
except:
293-
RuntimeError(
293+
raise RuntimeError(
294294
'Experiment list model does not have suffix ' + '"unknown_parameters".'
295295
)
296296

‎pyomo/contrib/parmest/tests/test_parmest.py‎

Lines changed: 94 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -80,6 +80,39 @@ def SSE(model):
8080
exp_list, obj_function=SSE, solver_options=solver_options, tee=True
8181
)
8282

83+
def test_parmest_exception(self):
84+
"""
85+
Test the exception raised by parmest when the "experiment_outputs"
86+
attribute is not defined in the model
87+
"""
88+
from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import (
89+
RooneyBieglerExperiment,
90+
)
91+
92+
# create an instance of the RooneyBieglerExperiment class
93+
# without the "experiment_outputs" attribute
94+
class RooneyBieglerExperimentException(RooneyBieglerExperiment):
95+
def label_model(self):
96+
m = self.model
97+
98+
# add the unknown parameters
99+
m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL)
100+
m.unknown_parameters.update(
101+
(k, pyo.ComponentUID(k)) for k in [m.asymptote, m.rate_constant]
102+
)
103+
104+
# create an experiment list
105+
exp_list = []
106+
for i in range(self.data.shape[0]):
107+
exp_list.append(RooneyBieglerExperimentException(self.data.loc[i, :]))
108+
109+
# check the exception raised by parmest due to not defining
110+
# the "experiment_outputs"
111+
with self.assertRaises(RuntimeError) as context:
112+
parmest.Estimator(exp_list, obj_function="SSE", tee=True)
113+
114+
self.assertIn("experiment_outputs", str(context.exception))
115+
83116
def test_theta_est(self):
84117
objval, thetavals = self.pest.theta_est()
85118

@@ -816,6 +849,22 @@ def label_model(self):
816849

817850
m = self.model
818851

852+
if isinstance(self.data, pd.DataFrame):
853+
meas_time_points = self.data.index
854+
else:
855+
meas_time_points = list(self.data["ca"].keys())
856+
857+
m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
858+
m.experiment_outputs.update(
859+
(m.ca[t], self.data["ca"][t]) for t in meas_time_points
860+
)
861+
m.experiment_outputs.update(
862+
(m.cb[t], self.data["cb"][t]) for t in meas_time_points
863+
)
864+
m.experiment_outputs.update(
865+
(m.cc[t], self.data["cc"][t]) for t in meas_time_points
866+
)
867+
819868
m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL)
820869
m.unknown_parameters.update(
821870
(k, pyo.ComponentUID(k)) for k in [m.k1, m.k2]
@@ -885,6 +934,51 @@ def get_labeled_model(self):
885934
self.m_df = ABC_model(data_df)
886935
self.m_dict = ABC_model(data_dict)
887936

937+
# create an instance of the ReactorDesignExperimentDAE class
938+
# without the "unknown_parameters" attribute
939+
class ReactorDesignExperimentException(ReactorDesignExperimentDAE):
940+
def label_model(self):
941+
942+
m = self.model
943+
944+
if isinstance(self.data, pd.DataFrame):
945+
meas_time_points = self.data.index
946+
else:
947+
meas_time_points = list(self.data["ca"].keys())
948+
949+
m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
950+
m.experiment_outputs.update(
951+
(m.ca[t], self.data["ca"][t]) for t in meas_time_points
952+
)
953+
m.experiment_outputs.update(
954+
(m.cb[t], self.data["cb"][t]) for t in meas_time_points
955+
)
956+
m.experiment_outputs.update(
957+
(m.cc[t], self.data["cc"][t]) for t in meas_time_points
958+
)
959+
960+
# create an experiment list without the "unknown_parameters" attribute
961+
exp_list_df_no_params = [ReactorDesignExperimentException(data_df)]
962+
exp_list_dict_no_params = [ReactorDesignExperimentException(data_dict)]
963+
964+
self.exp_list_df_no_params = exp_list_df_no_params
965+
self.exp_list_dict_no_params = exp_list_dict_no_params
966+
967+
def test_parmest_exception(self):
968+
"""
969+
Test the exception raised by parmest when the "unknown_parameters"
970+
attribute is not defined in the model
971+
"""
972+
with self.assertRaises(RuntimeError) as context:
973+
parmest.Estimator(self.exp_list_df_no_params, obj_function="SSE")
974+
975+
self.assertIn("unknown_parameters", str(context.exception))
976+
977+
with self.assertRaises(RuntimeError) as context:
978+
parmest.Estimator(self.exp_list_dict_no_params, obj_function="SSE")
979+
980+
self.assertIn("unknown_parameters", str(context.exception))
981+
888982
def test_dataformats(self):
889983
obj1, theta1 = self.pest_df.theta_est()
890984
obj2, theta2 = self.pest_dict.theta_est()

‎pyomo/contrib/sensitivity_toolbox/pynumero.py‎

Lines changed: 129 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,129 @@
1+
# ___________________________________________________________________________
2+
#
3+
# Pyomo: Python Optimization Modeling Objects
4+
# Copyright (c) 2008-2025
5+
# National Technology and Engineering Solutions of Sandia, LLC
6+
# Under the terms of Contract DE-NA0003525 with National Technology and
7+
# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
8+
# rights in this software.
9+
# This software is distributed under the 3-clause BSD License.
10+
# ___________________________________________________________________________
11+
12+
from pyomo.common.dependencies import numpy as np
13+
from pyomo.common.dependencies import scipy, attempt_import
14+
15+
import pyomo.environ as pyo
16+
17+
# Use attempt_import here due to unguarded NumPy import in this file
18+
nlp = attempt_import('pyomo.contrib.pynumero.interfaces.pyomo_nlp')[0]
19+
from pyomo.common.collections import ComponentSet, ComponentMap
20+
21+
22+
def _coo_reorder_cols(mat, remap):
23+
"""Change the order of columns in a COO matrix. The main use of this is
24+
to reorder variables in the Jacobian matrix. This changes the matrix in
25+
place.
26+
27+
Parameters
28+
----------
29+
mat: scipy.sparse.coo_matrix
30+
Reorder the columns of this matrix
31+
remap: dict
32+
dictionary where keys are old column and value is new column, if a column
33+
doesn't move, it doesn't need to be included.
34+
35+
Returns
36+
-------
37+
NoneType
38+
None
39+
"""
40+
for i in range(len(mat.data)):
41+
try:
42+
mat.col[i] = remap[mat.col[i]]
43+
except KeyError:
44+
pass # it's fine if we don't move a col in remap
45+
46+
47+
def get_dsdp_dfdp(model, theta):
48+
"""Calculate the derivatives of the state variables (s) with respect to
49+
parameters (p) (ds/dp), and the derivative of the objective function (f)
50+
with respect to p (df/dp). The number of parameters in theta should be the
51+
same as the number of degrees of freedom.
52+
53+
Parameters
54+
----------
55+
model: pyomo.environ.Block | pyomo.contrib.pynumero.interfaces.PyomoNLP
56+
Model to calculate sensitivity on. To retain the cached objects in
57+
the pynumero interface, create a PyomoNLP first and pass it to this function.
58+
theta: list
59+
A list of parameters as pyomo.environ.VarData, the number of parameters
60+
should be equal to the degrees of freedom.
61+
62+
Returns
63+
-------
64+
scipy.sparse.csc_matrix, csc_matrix, ComponentMap, ComponentMap
65+
ds/dp (ns by np), df/dp (1 by np), row map, column map.
66+
The column map maps Pyomo variables p to columns and the
67+
row map maps Pyomo variables s to rows.
68+
"""
69+
# Create a Pynumero NLP and get Jacobian
70+
if isinstance(model, nlp.PyomoNLP):
71+
m2 = model
72+
else:
73+
m2 = nlp.PyomoNLP(model)
74+
J = m2.evaluate_jacobian_eq()
75+
v_list = m2.get_pyomo_variables()
76+
# Map variables to columns in J
77+
mv_map = {id(v): i for i, v in enumerate(v_list)}
78+
s_list = list(ComponentSet(v_list) - ComponentSet(theta))
79+
ns = len(s_list)
80+
np = len(theta)
81+
col_remap = {mv_map[id(v)]: i for i, v in enumerate(s_list + theta)}
82+
_coo_reorder_cols(J, remap=col_remap)
83+
J = J.tocsc()
84+
dB = -(
85+
J
86+
@ scipy.sparse.vstack(
87+
(scipy.sparse.coo_matrix((ns, np)), scipy.sparse.identity(np))
88+
).tocsc()
89+
)
90+
# Calculate sensitivity matrix
91+
dsdp = scipy.sparse.linalg.spsolve(J[:, range(ns)], dB)
92+
# Get a map of state vars to columns
93+
s_map = {id(v): i for i, v in enumerate(s_list)}
94+
# Get the outputs we are interested in from the list of output vars
95+
column_map = ComponentMap([(v, i) for i, v in enumerate(theta)])
96+
row_map = ComponentMap([(v, i) for i, v in enumerate(s_list)])
97+
dfdx = scipy.sparse.coo_matrix(m2.evaluate_grad_objective())
98+
_coo_reorder_cols(dfdx, remap=col_remap)
99+
dfdx = dfdx.tocsc()
100+
dfdp = dfdx[0, :ns] @ dsdp + dfdx[0, ns:]
101+
# return sensitivity of the outputs to p and component maps
102+
return dsdp, dfdp, row_map, column_map
103+
104+
105+
def get_dydp(y_list, dsdp, row_map):
106+
"""Reduce the sensitivity matrix from get_dsdp_dfdp to only
107+
a specified set of state variables of interest.
108+
109+
Parameters
110+
----------
111+
y_list: list
112+
A list of state variables of interest (a subset of s)
113+
dsdp: csc_matrix
114+
A sensitivity matrix calculated by get_dsdp_dfdp
115+
row_map: ComponentMap
116+
A row map from get_dsdp_dfdp
117+
118+
Returns
119+
-------
120+
csc_matrix, ComponentMap
121+
dy/dp and a new row map with only y variables
122+
123+
"""
124+
new_row_map = ComponentMap()
125+
for i, v in enumerate(y_list):
126+
new_row_map[v] = i
127+
rows = [row_map[v] for v in y_list]
128+
dydp = dsdp[rows, :]
129+
return dydp, new_row_map

‎pyomo/contrib/sensitivity_toolbox/tests/test_pynumero.py‎

Lines changed: 107 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,107 @@
1+
# ___________________________________________________________________________
2+
#
3+
# Pyomo: Python Optimization Modeling Objects
4+
# Copyright (c) 2008-2025
5+
# National Technology and Engineering Solutions of Sandia, LLC
6+
# Under the terms of Contract DE-NA0003525 with National Technology and
7+
# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
8+
# rights in this software.
9+
# This software is distributed under the 3-clause BSD License.
10+
# ___________________________________________________________________________
11+
12+
13+
import pyomo.common.unittest as unittest
14+
from pyomo.common.dependencies import numpy as np, numpy_available
15+
from pyomo.common.dependencies import scipy_available, attempt_import
16+
17+
import pyomo.environ as pyo
18+
19+
# Use attempt_import here due to unguarded NumPy import in this file
20+
nlp = attempt_import('pyomo.contrib.pynumero.interfaces.pyomo_nlp')[0]
21+
import pyomo.contrib.sensitivity_toolbox.pynumero as pnsens
22+
from pyomo.contrib.pynumero.asl import AmplInterface
23+
24+
if not scipy_available or not numpy_available:
25+
raise unittest.SkipTest("scipy or numpy is not available")
26+
27+
if not AmplInterface.available():
28+
raise unittest.SkipTest("Pynumero needs the ASL extension to run NLP tests")
29+
30+
31+
class TestSeriesData(unittest.TestCase):
32+
def test_dsdp_dfdp_pyomo(self):
33+
m = pyo.ConcreteModel()
34+
m.x1 = pyo.Var(initialize=200)
35+
m.x2 = pyo.Var(initialize=5)
36+
m.p1 = pyo.Var(initialize=10)
37+
m.p2 = pyo.Var(initialize=5)
38+
m.obj = pyo.Objective(
39+
expr=m.x1 * m.p1 + m.x2 * m.x2 * m.p2 + m.p1 * m.p2, sense=pyo.minimize
40+
)
41+
m.c1 = pyo.Constraint(expr=m.x1 == 2 * m.p1**2)
42+
m.c2 = pyo.Constraint(expr=m.x2 == m.p2)
43+
theta = [m.p1, m.p2]
44+
45+
dsdp, dfdp, rmap, cmap = pnsens.get_dsdp_dfdp(m, theta)
46+
47+
# Since x1 = p1
48+
np.testing.assert_almost_equal(dsdp[rmap[m.x1], cmap[m.p1]], 40.0)
49+
np.testing.assert_almost_equal(dsdp[rmap[m.x1], cmap[m.p2]], 0.0)
50+
np.testing.assert_almost_equal(dsdp[rmap[m.x2], cmap[m.p1]], 0.0)
51+
np.testing.assert_almost_equal(dsdp[rmap[m.x2], cmap[m.p2]], 1.0)
52+
53+
# if x1 = 2 * p1 and x2 = p2 then
54+
# df/dp1 = 6 p1**2 + p2 = 45.0
55+
# df/dp2 = 3 p2 + p1 = 85.0
56+
np.testing.assert_almost_equal(dfdp[0, cmap[m.p1]], 605.0)
57+
np.testing.assert_almost_equal(dfdp[0, cmap[m.p2]], 85.0)
58+
59+
def test_dsdp_dfdp_pyomo_nlp(self):
60+
m = pyo.ConcreteModel()
61+
m.x1 = pyo.Var(initialize=200)
62+
m.x2 = pyo.Var(initialize=5)
63+
m.p1 = pyo.Var(initialize=10)
64+
m.p2 = pyo.Var(initialize=5)
65+
m.obj = pyo.Objective(
66+
expr=m.x1 * m.p1 + m.x2 * m.x2 * m.p2 + m.p1 * m.p2, sense=pyo.minimize
67+
)
68+
m.c1 = pyo.Constraint(expr=m.x1 == 2 * m.p1**2)
69+
m.c2 = pyo.Constraint(expr=m.x2 == m.p2)
70+
theta = [m.p1, m.p2]
71+
72+
m2 = nlp.PyomoNLP(m)
73+
dsdp, dfdp, rmap, cmap = pnsens.get_dsdp_dfdp(m2, theta)
74+
75+
# Since x1 = p1
76+
assert dsdp.shape == (2, 2)
77+
np.testing.assert_almost_equal(dsdp[rmap[m.x1], cmap[m.p1]], 40.0)
78+
np.testing.assert_almost_equal(dsdp[rmap[m.x1], cmap[m.p2]], 0.0)
79+
np.testing.assert_almost_equal(dsdp[rmap[m.x2], cmap[m.p1]], 0.0)
80+
np.testing.assert_almost_equal(dsdp[rmap[m.x2], cmap[m.p2]], 1.0)
81+
82+
# if x1 = 2 * p1 and x2 = p2 then
83+
# df/dp1 = 6 p1**2 + p2 = 45.0
84+
# df/dp2 = 3 p2 + p1 = 85.0
85+
assert dfdp.shape == (1, 2)
86+
np.testing.assert_almost_equal(dfdp[0, cmap[m.p1]], 605.0)
87+
np.testing.assert_almost_equal(dfdp[0, cmap[m.p2]], 85.0)
88+
89+
def test_dydp_pyomo(self):
90+
m = pyo.ConcreteModel()
91+
m.x1 = pyo.Var(initialize=200)
92+
m.x2 = pyo.Var(initialize=5)
93+
m.p1 = pyo.Var(initialize=10)
94+
m.p2 = pyo.Var(initialize=5)
95+
m.obj = pyo.Objective(
96+
expr=m.x1 * m.p1 + m.x2 * m.x2 * m.p2 + m.p1 * m.p2, sense=pyo.minimize
97+
)
98+
m.c1 = pyo.Constraint(expr=m.x1 == 2 * m.p1**2)
99+
m.c2 = pyo.Constraint(expr=m.x2 == m.p2)
100+
theta = [m.p1, m.p2]
101+
dsdp, dfdp, rmap, cmap = pnsens.get_dsdp_dfdp(m, theta)
102+
103+
dydp, rmap = pnsens.get_dydp([m.x2], dsdp, rmap)
104+
105+
np.testing.assert_almost_equal(dydp[rmap[m.x2], cmap[m.p1]], 0.0)
106+
np.testing.assert_almost_equal(dydp[rmap[m.x2], cmap[m.p2]], 1.0)
107+
assert dydp.shape == (1, 2)

0 commit comments

Comments
 (0)