Skip to content

Commit 7022bfd

Browse files
jeanpaulwatsonjeanpaulwatson
authored
Merge pull request #433 from DLWoodruff/nonant_sensitivities
move sensitivity code to utils from the extension that uses it
2 parents c4560bf + cdc7b8c commit 7022bfd

File tree

3 files changed

+105
-73
lines changed

3 files changed

+105
-73
lines changed

Diff for: ‎examples/generic_cylinders.bash

+10
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,16 @@
33

44
SOLVER="cplex"
55

6+
echo "^^^ netdes sensi-rho ^^^"
7+
cd netdes
8+
mpiexec -np 3 python -m mpi4py ../../mpisppy/generic_cylinders.py --module-name netdes --netdes-data-path ./data --instance-name network-10-20-L-01 --solver-name ${SOLVER} --max-iterations 10 --max-solver-threads 4 --default-rho 1 --lagrangian --xhatshuffle --rel-gap 0.01 --sensi-rho
9+
cd ..
10+
11+
echo "^^^ farmer sensi-rho ^^^"
12+
mpiexec -np 3 python -m mpi4py ../mpisppy/generic_cylinders.py --module-name farmer/farmer --num-scens 3 --solver-name ${SOLVER} --max-iterations 10 --max-solver-threads 4 --default-rho 1 --lagrangian --xhatshuffle --rel-gap 0.01 --sensi-rho
13+
14+
exit
15+
616
# sslp EF
717
echo "^^^ sslp ef ^^^"
818
cd sslp

Diff for: ‎mpisppy/extensions/sensi_rho.py

+5-73
Original file line numberDiff line numberDiff line change
@@ -9,12 +9,10 @@
99

1010
import numpy as np
1111

12-
import pyomo.environ as pyo
13-
from pyomo.contrib.pynumero.linalg.scipy_interface import ScipyLU
1412

1513
import mpisppy.extensions.extension
1614
import mpisppy.MPI as MPI
17-
from mpisppy.utils.kkt.interface import InteriorPointInterface
15+
from mpisppy.utils.nonant_sensitivities import nonant_sensitivies
1816

1917

2018
class SensiRho(mpisppy.extensions.extension.Extension):
@@ -142,81 +140,15 @@ def pre_iter0(self):
142140
def post_iter0(self):
143141
ph = self.ph
144142

145-
# first, solve the subproblems with Ipopt,
146-
# and gather sensitivity information
147-
ipopt = pyo.SolverFactory("ipopt")
148-
nonant_sensis = {}
143+
nonant_sensis = dict() # dict of dicts [s][ndn_i]
149144
for k, s in ph.local_subproblems.items():
150-
solution_cache = pyo.ComponentMap()
151-
for var in s.component_data_objects(pyo.Var):
152-
solution_cache[var] = var._value
153-
relax_int = pyo.TransformationFactory('core.relax_integer_vars')
154-
relax_int.apply_to(s)
155-
156-
assert hasattr(s, "_relaxed_integer_vars")
157-
158-
# add the needed suffixes / remove later
159-
s.ipopt_zL_out = pyo.Suffix(direction=pyo.Suffix.IMPORT)
160-
s.ipopt_zU_out = pyo.Suffix(direction=pyo.Suffix.IMPORT)
161-
s.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT_EXPORT)
162-
163-
results = ipopt.solve(s)
164-
pyo.assert_optimal_termination(results)
165-
166-
kkt_builder = InteriorPointInterface(s)
167-
kkt_builder.set_barrier_parameter(1e-9)
168-
kkt_builder.set_bounds_relaxation_factor(1e-8)
169-
#rhs = kkt_builder.evaluate_primal_dual_kkt_rhs()
170-
#print(f"{rhs}")
171-
#print(f"{rhs.flatten()}")
172-
kkt = kkt_builder.evaluate_primal_dual_kkt_matrix()
173-
174-
# print(f"{kkt=}")
175-
# could do better than SuperLU
176-
kkt_lu = ScipyLU()
177-
# always regularize equality constraints
178-
kkt_builder.regularize_equality_gradient(kkt=kkt, coef=-1e-8, copy_kkt=False)
179-
kkt_lu.do_numeric_factorization(kkt, raise_on_error=True)
180-
181-
grad_vec = np.zeros(kkt.shape[1])
182-
grad_vec[0:kkt_builder._nlp.n_primals()] = kkt_builder._nlp.evaluate_grad_objective()
183-
184-
grad_vec_kkt_inv = kkt_lu._lu.solve(grad_vec, "T")
185-
186-
for scenario_name in s.scen_list:
187-
nonant_sensis[scenario_name] = {}
188-
rho = ph.local_scenarios[scenario_name]._mpisppy_model.rho
189-
for ndn_i, v in ph.local_scenarios[scenario_name]._mpisppy_data.nonant_indices.items():
190-
var_idx = kkt_builder._nlp._vardata_to_idx[v]
191-
192-
y_vec = np.zeros(kkt.shape[0])
193-
y_vec[var_idx] = 1.0
194-
195-
x_denom = y_vec.T @ kkt_lu._lu.solve(y_vec)
196-
x = (-1 / x_denom)
197-
e_x = x * y_vec
198-
199-
sensitivity = grad_vec_kkt_inv @ -e_x
200-
# print(f"df/d{v.name}: {sensitivity:.2e}, ∂f/∂{v.name}: {grad_vec[var_idx]:.2e}, "
201-
# f"rho {v.name}: {ph.local_scenarios[scenario_name]._mpisppy_model.rho[ndn_i]._value:.2e}, ",
202-
# f"value: {v._value:.2e}"
203-
# )
204-
205-
rho[ndn_i]._value = abs(sensitivity)
206-
207-
relax_int.apply_to(s, options={"undo":True})
208-
assert not hasattr(s, "_relaxed_integer_vars")
209-
del s.ipopt_zL_out
210-
del s.ipopt_zU_out
211-
del s.dual
212-
for var, val in solution_cache.items():
213-
var._value = val
214-
145+
nonant_sensis[s] = nonant_sensitivies(s, ph)
146+
215147
for s in ph.local_scenarios.values():
216148
xbars = s._mpisppy_model.xbars
217149
for ndn_i, rho in s._mpisppy_model.rho.items():
218150
nv = s._mpisppy_data.nonant_indices[ndn_i] # var_data object
219-
rho._value = rho._value / max(1, abs(nv._value - xbars[ndn_i]._value))
151+
rho._value = abs(nonant_sensis[s][ndn_i]) / max(1, abs(nv._value - xbars[ndn_i]._value))
220152
rho._value *= self.multiplier
221153
# if ph.cylinder_rank == 0:
222154
# print(f"{s.name=}, {nv.name=}, {rho.value=}")

Diff for: ‎mpisppy/utils/nonant_sensitivities.py

+90
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,90 @@
1+
###############################################################################
2+
# mpi-sppy: MPI-based Stochastic Programming in PYthon
3+
#
4+
# Copyright (c) 2024, Lawrence Livermore National Security, LLC, Alliance for
5+
# Sustainable Energy, LLC, The Regents of the University of California, et al.
6+
# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md for
7+
# full copyright and license information.
8+
###############################################################################
9+
10+
import numpy as np
11+
12+
import pyomo.environ as pyo
13+
from pyomo.contrib.pynumero.linalg.scipy_interface import ScipyLU
14+
15+
from mpisppy.utils.kkt.interface import InteriorPointInterface
16+
17+
def nonant_sensitivies(s, ph):
18+
""" Compute the sensitivities of noants (w.r.t. the Lagrangian for s)
19+
Args:
20+
s: (Pyomo ConcreteModel): the scenario
21+
ph: (PHBase Object): to deal with bundles (that are not proper)
22+
Returns:
23+
nonant_sensis (dict): [ndn_i]: sensitivity for the Var
24+
"""
25+
26+
# first, solve the subproblems with Ipopt,
27+
# and gather sensitivity information
28+
ipopt = pyo.SolverFactory("ipopt")
29+
solution_cache = pyo.ComponentMap()
30+
for var in s.component_data_objects(pyo.Var):
31+
solution_cache[var] = var._value
32+
relax_int = pyo.TransformationFactory('core.relax_integer_vars')
33+
relax_int.apply_to(s)
34+
35+
assert hasattr(s, "_relaxed_integer_vars")
36+
37+
# add the needed suffixes / remove later
38+
s.ipopt_zL_out = pyo.Suffix(direction=pyo.Suffix.IMPORT)
39+
s.ipopt_zU_out = pyo.Suffix(direction=pyo.Suffix.IMPORT)
40+
s.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT_EXPORT)
41+
42+
results = ipopt.solve(s)
43+
pyo.assert_optimal_termination(results)
44+
45+
kkt_builder = InteriorPointInterface(s)
46+
kkt_builder.set_barrier_parameter(1e-9)
47+
kkt_builder.set_bounds_relaxation_factor(1e-8)
48+
#rhs = kkt_builder.evaluate_primal_dual_kkt_rhs()
49+
#print(f"{rhs}")
50+
#print(f"{rhs.flatten()}")
51+
kkt = kkt_builder.evaluate_primal_dual_kkt_matrix()
52+
53+
# print(f"{kkt=}")
54+
# could do better than SuperLU
55+
kkt_lu = ScipyLU()
56+
# always regularize equality constraints
57+
kkt_builder.regularize_equality_gradient(kkt=kkt, coef=-1e-8, copy_kkt=False)
58+
kkt_lu.do_numeric_factorization(kkt, raise_on_error=True)
59+
60+
grad_vec = np.zeros(kkt.shape[1])
61+
grad_vec[0:kkt_builder._nlp.n_primals()] = kkt_builder._nlp.evaluate_grad_objective()
62+
63+
grad_vec_kkt_inv = kkt_lu._lu.solve(grad_vec, "T")
64+
65+
nonant_sensis = dict()
66+
# bundles?
67+
for scenario_name in s.scen_list:
68+
for ndn_i, v in ph.local_scenarios[scenario_name]._mpisppy_data.nonant_indices.items():
69+
var_idx = kkt_builder._nlp._vardata_to_idx[v]
70+
71+
y_vec = np.zeros(kkt.shape[0])
72+
y_vec[var_idx] = 1.0
73+
74+
x_denom = y_vec.T @ kkt_lu._lu.solve(y_vec)
75+
x = (-1 / x_denom)
76+
e_x = x * y_vec
77+
78+
sensitivity = grad_vec_kkt_inv @ -e_x
79+
#rho[ndn_i]._value = abs(sensitivity)
80+
nonant_sensis[ndn_i] = sensitivity
81+
82+
relax_int.apply_to(s, options={"undo":True})
83+
assert not hasattr(s, "_relaxed_integer_vars")
84+
del s.ipopt_zL_out
85+
del s.ipopt_zU_out
86+
del s.dual
87+
for var, val in solution_cache.items():
88+
var._value = val
89+
90+
return nonant_sensis

0 commit comments

Comments
 (0)