diff --git a/mpisppy/extensions/sensi_rho.py b/mpisppy/extensions/sensi_rho.py new file mode 100644 index 00000000..880bc5ab --- /dev/null +++ b/mpisppy/extensions/sensi_rho.py @@ -0,0 +1,243 @@ +############################################################################### +# mpi-sppy: MPI-based Stochastic Programming in PYthon +# +# Copyright (c) 2024, Lawrence Livermore National Security, LLC, Alliance for +# Sustainable Energy, LLC, The Regents of the University of California, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md for +# full copyright and license information. +############################################################################### + +import numpy as np + +import pyomo.environ as pyo +from pyomo.contrib.pynumero.linalg.scipy_interface import ScipyLU + +import mpisppy.extensions.extension +import mpisppy.MPI as MPI +from mpisppy.utils.kkt.interface import InteriorPointInterface + + +class SensiRho(mpisppy.extensions.extension.Extension): + """ + Rho determination algorithm using nonant sensitivities + """ + + def __init__(self, ph): + self.ph = ph + + self.multiplier = 1.0 + + if ( + "sensi_rho_options" in ph.options + and "multiplier" in ph.options["sensi_rho_options"] + ): + self.multiplier = ph.options["sensi_rho_options"]["multiplier"] + + @staticmethod + def _compute_rho_min_max(ph, npop, mpiop, start): + local_nodenames = [] + local_xmaxmin = {} + global_xmaxmin = {} + + for k, s in ph.local_scenarios.items(): + nlens = s._mpisppy_data.nlens + for node in s._mpisppy_node_list: + if node.name not in local_nodenames: + ndn = node.name + local_nodenames.append(ndn) + nlen = nlens[ndn] + + local_xmaxmin[ndn] = start * np.ones(nlen, dtype="d") + global_xmaxmin[ndn] = np.zeros(nlen, dtype="d") + + for k, s in ph.local_scenarios.items(): + nlens = s._mpisppy_data.nlens + rho = s._mpisppy_model.rho + for node in s._mpisppy_node_list: + ndn = node.name + xmaxmin = local_xmaxmin[ndn] + + xmaxmin_partial = np.fromiter( + (rho[ndn,i]._value for i, _ in enumerate(node.nonant_vardata_list)), + dtype="d", + count=nlens[ndn], + ) + xmaxmin = npop(xmaxmin, xmaxmin_partial) + local_xmaxmin[ndn] = xmaxmin + + for nodename in local_nodenames: + ph.comms[nodename].Allreduce( + [local_xmaxmin[nodename], MPI.DOUBLE], + [global_xmaxmin[nodename], MPI.DOUBLE], + op=mpiop, + ) + + xmaxmin_dict = {} + for ndn, global_xmaxmin_dict in global_xmaxmin.items(): + for i, v in enumerate(global_xmaxmin_dict): + xmaxmin_dict[ndn, i] = v + + return xmaxmin_dict + + @staticmethod + def _compute_rho_avg(ph): + local_nodenames = [] + local_avg = {} + global_avg = {} + + for k, s in ph.local_scenarios.items(): + nlens = s._mpisppy_data.nlens + rho = s._mpisppy_model.rho + for node in s._mpisppy_node_list: + if node.name not in local_nodenames: + ndn = node.name + local_nodenames.append(ndn) + nlen = nlens[ndn] + + local_avg[ndn] = np.zeros(nlen, dtype="d") + global_avg[ndn] = np.zeros(nlen, dtype="d") + + for k, s in ph.local_scenarios.items(): + nlens = s._mpisppy_data.nlens + rho = s._mpisppy_model.rho + for node in s._mpisppy_node_list: + ndn = node.name + + local_rhos = np.fromiter( + (rho[ndn,i]._value for i, _ in enumerate(node.nonant_vardata_list)), + dtype="d", + count=nlens[ndn], + ) + # print(f"{k=}, {local_rhos=}, {s._mpisppy_probability=}, {s._mpisppy_data.prob_coeff[ndn]=}") + # TODO: is this the right thing, or should it be s._mpisppy_probability? + local_rhos *= s._mpisppy_data.prob_coeff[ndn] + + local_avg[ndn] += local_rhos + + for nodename in local_nodenames: + ph.comms[nodename].Allreduce( + [local_avg[nodename], MPI.DOUBLE], + [global_avg[nodename], MPI.DOUBLE], + op=MPI.SUM, + ) + + rhoavg_dict = {} + for ndn, global_rhoavg_dict in global_avg.items(): + for i, v in enumerate(global_rhoavg_dict): + rhoavg_dict[ndn, i] = v + + return rhoavg_dict + + @staticmethod + def _compute_rho_max(ph): + return SensiRho._compute_rho_min_max(ph, np.maximum, MPI.MAX, -np.inf) + + @staticmethod + def _compute_rho_min(ph): + return SensiRho._compute_rho_min_max(ph, np.minimum, MPI.MIN, np.inf) + + def pre_iter0(self): + pass + + def post_iter0(self): + ph = self.ph + + # first, solve the subproblems with Ipopt, + # and gather sensitivity information + ipopt = pyo.SolverFactory("ipopt") + nonant_sensis = {} + for k, s in ph.local_subproblems.items(): + solution_cache = pyo.ComponentMap() + for var in s.component_data_objects(pyo.Var): + solution_cache[var] = var._value + relax_int = pyo.TransformationFactory('core.relax_integer_vars') + relax_int.apply_to(s) + + assert hasattr(s, "_relaxed_integer_vars") + + # add the needed suffixes / remove later + s.ipopt_zL_out = pyo.Suffix(direction=pyo.Suffix.IMPORT) + s.ipopt_zU_out = pyo.Suffix(direction=pyo.Suffix.IMPORT) + s.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT_EXPORT) + + results = ipopt.solve(s) + pyo.assert_optimal_termination(results) + + kkt_builder = InteriorPointInterface(s) + kkt_builder.set_barrier_parameter(1e-9) + kkt_builder.set_bounds_relaxation_factor(1e-8) + #rhs = kkt_builder.evaluate_primal_dual_kkt_rhs() + #print(f"{rhs}") + #print(f"{rhs.flatten()}") + kkt = kkt_builder.evaluate_primal_dual_kkt_matrix() + + # print(f"{kkt=}") + # could do better than SuperLU + kkt_lu = ScipyLU() + # always regularize equality constraints + kkt_builder.regularize_equality_gradient(kkt=kkt, coef=-1e-8, copy_kkt=False) + kkt_lu.do_numeric_factorization(kkt, raise_on_error=True) + + grad_vec = np.zeros(kkt.shape[1]) + grad_vec[0:kkt_builder._nlp.n_primals()] = kkt_builder._nlp.evaluate_grad_objective() + + grad_vec_kkt_inv = kkt_lu._lu.solve(grad_vec, "T") + + for scenario_name in s.scen_list: + nonant_sensis[scenario_name] = {} + rho = ph.local_scenarios[scenario_name]._mpisppy_model.rho + for ndn_i, v in ph.local_scenarios[scenario_name]._mpisppy_data.nonant_indices.items(): + var_idx = kkt_builder._nlp._vardata_to_idx[v] + + y_vec = np.zeros(kkt.shape[0]) + y_vec[var_idx] = 1.0 + + x_denom = y_vec.T @ kkt_lu._lu.solve(y_vec) + x = (-1 / x_denom) + e_x = x * y_vec + + sensitivity = grad_vec_kkt_inv @ -e_x + # print(f"df/d{v.name}: {sensitivity:.2e}, ∂f/∂{v.name}: {grad_vec[var_idx]:.2e}, " + # f"rho {v.name}: {ph.local_scenarios[scenario_name]._mpisppy_model.rho[ndn_i]._value:.2e}, ", + # f"value: {v._value:.2e}" + # ) + + rho[ndn_i]._value = abs(sensitivity) + + relax_int.apply_to(s, options={"undo":True}) + assert not hasattr(s, "_relaxed_integer_vars") + del s.ipopt_zL_out + del s.ipopt_zU_out + del s.dual + for var, val in solution_cache.items(): + var._value = val + + for s in ph.local_scenarios.values(): + xbars = s._mpisppy_model.xbars + for ndn_i, rho in s._mpisppy_model.rho.items(): + nv = s._mpisppy_data.nonant_indices[ndn_i] # var_data object + rho._value = rho._value / max(1, abs(nv._value - xbars[ndn_i]._value)) + rho._value *= self.multiplier + # if ph.cylinder_rank == 0: + # print(f"{s.name=}, {nv.name=}, {rho.value=}") + + rhoavg = self._compute_rho_avg(ph) + for s in ph.local_scenarios.values(): + xbars = s._mpisppy_model.xbars + for ndn_i, rho in s._mpisppy_model.rho.items(): + rho._value = rhoavg[ndn_i] + # if ph.cylinder_rank == 0: + # nv = s._mpisppy_data.nonant_indices[ndn_i] # var_data object + # print(f"{s.name=}, {nv.name=}, {rho.value=}") + + if ph.cylinder_rank == 0: + print("Rho values updated by SensiRho Extension") + + def miditer(self): + pass + + def enditer(self): + pass + + def post_everything(self): + pass diff --git a/mpisppy/generic_cylinders.py b/mpisppy/generic_cylinders.py index ebd009f2..70ab324c 100644 --- a/mpisppy/generic_cylinders.py +++ b/mpisppy/generic_cylinders.py @@ -72,6 +72,7 @@ def _parse_args(m): cfg.reduced_costs_args() cfg.sep_rho_args() cfg.coeff_rho_args() + cfg.sensi_rho_args() cfg.parse_command_line(f"mpi-sppy for {cfg.module_name}") return cfg @@ -99,7 +100,7 @@ def _name_lists(module, cfg): def _do_decomp(module, cfg, scenario_creator, scenario_creator_kwargs, scenario_denouement): rho_setter = module._rho_setter if hasattr(module, '_rho_setter') else None if cfg.default_rho is None and rho_setter is None: - if cfg.sep_rho or cfg.coeff_rho: + if cfg.sep_rho or cfg.coeff_rho or cfg.sensi_rho: cfg.default_rho = 1 else: raise RuntimeError("No rho_setter so a default must be specified via --default-rho") @@ -167,6 +168,9 @@ def _do_decomp(module, cfg, scenario_creator, scenario_creator_kwargs, scenario_ if cfg.coeff_rho: vanilla.add_coeff_rho(hub_dict, cfg) + + if cfg.sensi_rho: + vanilla.add_sensi_rho(hub_dict, cfg) if len(ext_classes) != 0: hub_dict['opt_kwargs']['extensions'] = MultiExtension @@ -215,6 +219,9 @@ def _do_decomp(module, cfg, scenario_creator, scenario_creator_kwargs, scenario_ vanilla.add_sep_rho(ph_ob_spoke, cfg) if cfg.coeff_rho: vanilla.add_coeff_rho(ph_ob_spoke, cfg) + if cfg.sensi_rho: + vanilla.add_sensi_rho(ph_ob_spoke, cfg) + # subgradient outer bound spoke if cfg.subgradient: @@ -227,6 +234,8 @@ def _do_decomp(module, cfg, scenario_creator, scenario_creator_kwargs, scenario_ vanilla.add_sep_rho(subgradient_spoke, cfg) if cfg.coeff_rho: vanilla.add_coeff_rho(subgradient_spoke, cfg) + if cfg.sensi_rho: + vanilla.add_sensi_rho(subgradient_spoke, cfg) # xhat shuffle bound spoke if cfg.xhatshuffle: diff --git a/mpisppy/utils/cfg_vanilla.py b/mpisppy/utils/cfg_vanilla.py index e09e8e67..59ba3d86 100644 --- a/mpisppy/utils/cfg_vanilla.py +++ b/mpisppy/utils/cfg_vanilla.py @@ -42,6 +42,7 @@ from mpisppy.extensions.reduced_costs_fixer import ReducedCostsFixer from mpisppy.extensions.sep_rho import SepRho from mpisppy.extensions.coeff_rho import CoeffRho +from mpisppy.extensions.sensi_rho import SensiRho from mpisppy.utils.wxbarreader import WXBarReader from mpisppy.utils.wxbarwriter import WXBarWriter @@ -209,6 +210,10 @@ def add_coeff_rho(hub_dict, cfg): hub_dict = extension_adder(hub_dict,CoeffRho) hub_dict["opt_kwargs"]["options"]["coeff_rho_options"] = {"multiplier" : cfg.coeff_rho_multiplier} +def add_sensi_rho(hub_dict, cfg): + hub_dict = extension_adder(hub_dict,SensiRho) + hub_dict["opt_kwargs"]["options"]["sensi_rho_options"] = {"multiplier" : cfg.sensi_rho_multiplier} + def add_cross_scenario_cuts(hub_dict, cfg, ): diff --git a/mpisppy/utils/config.py b/mpisppy/utils/config.py index 6955be40..0e72b1f1 100644 --- a/mpisppy/utils/config.py +++ b/mpisppy/utils/config.py @@ -421,6 +421,17 @@ def sep_rho_args(self): default=1.0) + def sensi_rho_args(self): + self.add_to_config("sensi_rho", + description="have a SensiRho extension", + domain=bool, + default=False) + self.add_to_config("sensi_rho_multiplier", + description="multiplier for SensiRho (default 1.0)", + domain=float, + default=1.0) + + def coeff_rho_args(self): self.add_to_config("coeff_rho", description="have a CoeffRho extension", diff --git a/mpisppy/utils/kkt/LICENSE.md b/mpisppy/utils/kkt/LICENSE.md new file mode 100644 index 00000000..5f6adcaf --- /dev/null +++ b/mpisppy/utils/kkt/LICENSE.md @@ -0,0 +1,185 @@ +BSD 3-Clause License + +Copyright Notice +================= + +Copyright 2020 National Technology & Engineering Solutions of Sandia, LLC +(NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. +Government retains certain rights in this software. + +License Notice +================= + +This software is distributed under the Revised BSD License (see +below). Parapint also leverages a variety of third-party software +packages, which have separate licensing policies. + +Revised BSD License +------------------- + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions +are met: + +* Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. +* Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. +* Neither the name of Sandia National Laboratories, nor the names of + its contributors may be used to endorse or promote products derived from + this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED +TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + +Third-Party Libraries +===================== + +Parapint links to or includes source code from the following +third-party libraries: + + +Pyomo +----- + +Copyright 2017 National Technology and Engineering Solutions of +Sandia, LLC . Under the terms of Contract DE-NA0003525 with National +Technology and Engineering Solutions of Sandia, LLC , the U.S. +Government retains certain rights in this software. + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions +are met: + +* Redistributions of source code must retain the above copyright notice, +this list of conditions and the following disclaimer. + +* Redistributions in binary form must reproduce the above copyright +notice, this list of conditions and the following disclaimer in the +documentation and/or other materials provided with the distribution. + +* Neither the name of the Sandia National Laboratories nor the names of +its contributors may be used to endorse or promote products derived from +this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED +TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + +NumPy +----- +Copyright (c) 2005-2020, NumPy Developers. +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of the NumPy Developers nor the names of any + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + +SciPy +----- +Copyright (c) 2001-2002 Enthought, Inc. 2003-2019, SciPy Developers. +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions +are met: + +1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + +2. Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + +3. Neither the name of the copyright holder nor the names of its + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + +MPI for Python +-------------- +Copyright (c) 2019, Lisandro Dalcin. +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions +are met: + +* Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + +* Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. \ No newline at end of file diff --git a/mpisppy/utils/kkt/interface.py b/mpisppy/utils/kkt/interface.py new file mode 100644 index 00000000..63913160 --- /dev/null +++ b/mpisppy/utils/kkt/interface.py @@ -0,0 +1,446 @@ +############################################################################### +# mpi-sppy: MPI-based Stochastic Programming in PYthon +# +# Copyright (c) 2024, Lawrence Livermore National Security, LLC, Alliance for +# Sustainable Energy, LLC, The Regents of the University of California, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md for +# full copyright and license information. +############################################################################### +############################################################################### +# This file was originally part of parapint, available: +# https://github.com/sandialabs/parapint/ +# See LICENSE.md in this directory for full copyright and license information. +############################################################################### + +from pyomo.contrib.pynumero.interfaces import pyomo_nlp +from pyomo.contrib.pynumero.sparse import BlockMatrix, BlockVector +import numpy as np +import scipy.sparse +from pyomo.common.timing import HierarchicalTimer + +class InteriorPointInterface: + def __init__(self, pyomo_model): + self._nlp = pyomo_nlp.PyomoNLP(pyomo_model, nl_file_options={'skip_trivial_constraints': True}) + + self.bounds_relaxation_factor = 0 + + self._slacks = self.init_slacks() + + # set the init_duals_primals_lb/ub from ipopt_zL_out, ipopt_zU_out if available + # need to compress them as well and initialize the duals_primals_lb/ub + self._init_duals_primals_lb, self._init_duals_primals_ub =\ + self._get_full_duals_primals_bounds() + self._init_duals_primals_lb[np.isneginf(self.primals_lb())] = 0 + self._init_duals_primals_ub[np.isinf(self.primals_ub())] = 0 + self._duals_primals_lb = self._init_duals_primals_lb.copy() + self._duals_primals_ub = self._init_duals_primals_ub.copy() + + # set the init_duals_slacks_lb/ub from the init_duals_ineq + # need to be compressed and set according to their sign + # (-) value indicates it the upper is active, while (+) indicates + # that lower is active + self._init_duals_slacks_lb = self._nlp.init_duals_ineq().copy() + self._init_duals_slacks_lb[self._init_duals_slacks_lb < 0] = 0 + self._init_duals_slacks_ub = self._nlp.init_duals_ineq().copy() + self._init_duals_slacks_ub[self._init_duals_slacks_ub > 0] = 0 + self._init_duals_slacks_ub *= -1.0 + + self._duals_slacks_lb = self._init_duals_slacks_lb.copy() + self._duals_slacks_ub = self._init_duals_slacks_ub.copy() + + self._delta_primals = None + self._delta_slacks = None + self._delta_duals_eq = None + self._delta_duals_ineq = None + self._barrier = None + + def get_bounds_relaxation_factor(self) -> float: + return self.bounds_relaxation_factor + + def set_bounds_relaxation_factor(self, val: float): + self.bounds_relaxation_factor = val + + def n_primals(self): + return self._nlp.n_primals() + + def nnz_hessian_lag(self): + return self._nlp.nnz_hessian_lag() + + def set_obj_factor(self, obj_factor): + self._nlp.set_obj_factor(obj_factor) + + def get_obj_factor(self): + return self._nlp.get_obj_factor() + + def n_eq_constraints(self): + return self._nlp.n_eq_constraints() + + def n_ineq_constraints(self): + return self._nlp.n_ineq_constraints() + + def nnz_jacobian_eq(self): + return self._nlp.nnz_jacobian_eq() + + def nnz_jacobian_ineq(self): + return self._nlp.nnz_jacobian_ineq() + + def init_primals(self): + primals = self._nlp.init_primals() + return primals + + def init_slacks(self): + slacks = self._nlp.evaluate_ineq_constraints() + return slacks + + def init_duals_eq(self): + return self._nlp.init_duals_eq() + + def init_duals_ineq(self): + return self._nlp.init_duals_ineq() + + def init_duals_primals_lb(self): + return self._init_duals_primals_lb + + def init_duals_primals_ub(self): + return self._init_duals_primals_ub + + def init_duals_slacks_lb(self): + return self._init_duals_slacks_lb + + def init_duals_slacks_ub(self): + return self._init_duals_slacks_ub + + def set_primals(self, primals): + self._nlp.set_primals(primals) + + def set_slacks(self, slacks): + self._slacks = slacks + + def set_duals_eq(self, duals): + self._nlp.set_duals_eq(duals) + + def set_duals_ineq(self, duals): + self._nlp.set_duals_ineq(duals) + + def set_duals_primals_lb(self, duals): + self._duals_primals_lb = duals + + def set_duals_primals_ub(self, duals): + self._duals_primals_ub = duals + + def set_duals_slacks_lb(self, duals): + self._duals_slacks_lb = duals + + def set_duals_slacks_ub(self, duals): + self._duals_slacks_ub = duals + + def get_primals(self): + return self._nlp.get_primals() + + def get_slacks(self): + return self._slacks + + def get_duals_eq(self): + return self._nlp.get_duals_eq() + + def get_duals_ineq(self): + return self._nlp.get_duals_ineq() + + def get_duals_primals_lb(self): + return self._duals_primals_lb + + def get_duals_primals_ub(self): + return self._duals_primals_ub + + def get_duals_slacks_lb(self): + return self._duals_slacks_lb + + def get_duals_slacks_ub(self): + return self._duals_slacks_ub + + def primals_lb(self): + lbs = self._nlp.primals_lb() + if self.bounds_relaxation_factor == 0: + return lbs + eye = np.ones(lbs.size) + lbs_mod = lbs - self.bounds_relaxation_factor * np.max(np.array([eye, np.abs(lbs)]), axis=0) + return lbs_mod + + def primals_ub(self): + ubs = self._nlp.primals_ub() + if self.bounds_relaxation_factor == 0: + return ubs + eye = np.ones(ubs.size) + ubs_mod = ubs + self.bounds_relaxation_factor * np.max(np.array([eye, np.abs(ubs)]), axis=0) + return ubs_mod + + def ineq_lb(self): + lbs = self._nlp.ineq_lb() + if self.bounds_relaxation_factor == 0: + return lbs + eye = np.ones(lbs.size) + lbs_mod = lbs - self.bounds_relaxation_factor * np.max(np.array([eye, np.abs(lbs)]), axis=0) + return lbs_mod + + def ineq_ub(self): + ubs = self._nlp.ineq_ub() + if self.bounds_relaxation_factor == 0: + return ubs + eye = np.ones(ubs.size) + ubs_mod = ubs + self.bounds_relaxation_factor * np.max(np.array([eye, np.abs(ubs)]), axis=0) + return ubs_mod + + def set_barrier_parameter(self, barrier): + self._barrier = barrier + + def pyomo_nlp(self): + return self._nlp + + def evaluate_primal_dual_kkt_matrix(self, timer=None): + if timer is None: + timer = HierarchicalTimer() + timer.start('eval hess') + hess_block = self._nlp.evaluate_hessian_lag() + timer.stop('eval hess') + timer.start('eval jac') + jac_eq = self._nlp.evaluate_jacobian_eq() + jac_ineq = self._nlp.evaluate_jacobian_ineq() + timer.stop('eval jac') + + duals_primals_lb = self._duals_primals_lb + duals_primals_ub = self._duals_primals_ub + duals_slacks_lb = self._duals_slacks_lb + duals_slacks_ub = self._duals_slacks_ub + primals = self._nlp.get_primals() + + timer.start('hess block') + data = (duals_primals_lb/(primals - self.primals_lb()) + + duals_primals_ub/(self.primals_ub() - primals)) + n = self._nlp.n_primals() + indices = np.arange(n) + hess_block.row = np.concatenate([hess_block.row, indices]) + hess_block.col = np.concatenate([hess_block.col, indices]) + hess_block.data = np.concatenate([hess_block.data, data]) + timer.stop('hess block') + + timer.start('slack block') + data = (duals_slacks_lb/(self._slacks - self.ineq_lb()) + + duals_slacks_ub/(self.ineq_ub() - self._slacks)) + n = self._nlp.n_ineq_constraints() + indices = np.arange(n) + slack_block = scipy.sparse.coo_matrix((data, (indices, indices)), shape=(n, n)) + timer.stop('slack block') + + timer.start('regularization block') + eq_reg_blk = scipy.sparse.identity(self._nlp.n_eq_constraints(), format='coo') + eq_reg_blk.data.fill(0) + ineq_reg_blk = scipy.sparse.identity(self._nlp.n_ineq_constraints(), format='coo') + ineq_reg_blk.data.fill(0) + timer.stop('regularization block') + + timer.start('set block') + kkt = BlockMatrix(4, 4) + kkt.set_block(0, 0, hess_block) + kkt.set_block(1, 1, slack_block) + kkt.set_block(2, 0, jac_eq) + kkt.set_block(0, 2, jac_eq.transpose()) + kkt.set_block(3, 0, jac_ineq) + kkt.set_block(0, 3, jac_ineq.transpose()) + kkt.set_block(3, 1, -scipy.sparse.identity( + self._nlp.n_ineq_constraints(), + format='coo')) + kkt.set_block(1, 3, -scipy.sparse.identity( + self._nlp.n_ineq_constraints(), + format='coo')) + kkt.set_block(2, 2, eq_reg_blk) + kkt.set_block(3, 3, ineq_reg_blk) + timer.stop('set block') + return kkt + + def evaluate_primal_dual_kkt_rhs(self, timer=None): + if timer is None: + timer = HierarchicalTimer() + timer.start('eval grad obj') + grad_obj = self.get_obj_factor() * self.evaluate_grad_objective() + timer.stop('eval grad obj') + timer.start('eval jac') + jac_eq = self._nlp.evaluate_jacobian_eq() + jac_ineq = self._nlp.evaluate_jacobian_ineq() + timer.stop('eval jac') + timer.start('eval cons') + eq_resid = self._nlp.evaluate_eq_constraints() + ineq_resid = self._nlp.evaluate_ineq_constraints() - self._slacks + timer.stop('eval cons') + + timer.start('grad_lag_primals') + grad_lag_primals = (grad_obj + + jac_eq.transpose() * self._nlp.get_duals_eq() + + jac_ineq.transpose() * self._nlp.get_duals_ineq() - + self._barrier / (self._nlp.get_primals() - self.primals_lb()) + + self._barrier / (self.primals_ub() - self._nlp.get_primals())) + timer.stop('grad_lag_primals') + + timer.start('grad_lag_slacks') + grad_lag_slacks = (-self._nlp.get_duals_ineq() - + self._barrier / (self._slacks - self.ineq_lb()) + + self._barrier / (self.ineq_ub() - self._slacks)) + timer.stop('grad_lag_slacks') + + rhs = BlockVector(4) + rhs.set_block(0, grad_lag_primals) + rhs.set_block(1, grad_lag_slacks) + rhs.set_block(2, eq_resid) + rhs.set_block(3, ineq_resid) + rhs = -rhs + return rhs + + def set_primal_dual_kkt_solution(self, sol): + self._delta_primals = sol.get_block(0) + self._delta_slacks = sol.get_block(1) + self._delta_duals_eq = sol.get_block(2) + self._delta_duals_ineq = sol.get_block(3) + + def get_delta_primals(self): + return self._delta_primals + + def get_delta_slacks(self): + return self._delta_slacks + + def get_delta_duals_eq(self): + return self._delta_duals_eq + + def get_delta_duals_ineq(self): + return self._delta_duals_ineq + + def get_delta_duals_primals_lb(self): + res = (((self._barrier - self._duals_primals_lb * self._delta_primals) / + (self._nlp.get_primals() - self.primals_lb())) - + self._duals_primals_lb) + return res + + def get_delta_duals_primals_ub(self): + res = (((self._barrier + self._duals_primals_ub * self._delta_primals) / + (self.primals_ub() - self._nlp.get_primals())) - + self._duals_primals_ub) + return res + + def get_delta_duals_slacks_lb(self): + res = (((self._barrier - self._duals_slacks_lb * self._delta_slacks) / + (self._slacks - self.ineq_lb())) - + self._duals_slacks_lb) + return res + + def get_delta_duals_slacks_ub(self): + res = (((self._barrier + self._duals_slacks_ub * self._delta_slacks) / + (self.ineq_ub() - self._slacks)) - + self._duals_slacks_ub) + return res + + def evaluate_objective(self): + return self._nlp.evaluate_objective() + + def evaluate_eq_constraints(self): + return self._nlp.evaluate_eq_constraints() + + def evaluate_ineq_constraints(self): + return self._nlp.evaluate_ineq_constraints() + + def evaluate_grad_objective(self): + return self._nlp.evaluate_grad_objective() + + def evaluate_jacobian_eq(self): + return self._nlp.evaluate_jacobian_eq() + + def evaluate_jacobian_ineq(self): + return self._nlp.evaluate_jacobian_ineq() + + def regularize_equality_gradient(self, kkt, coef, copy_kkt=True): + # Not technically regularizing the equality gradient ... + # Replace this with a regularize_diagonal_block function? + # Then call with kkt matrix and the value of the perturbation? + + # Use a constant perturbation to regularize the equality constraint + # gradient + if copy_kkt: + kkt = kkt.copy() + reg_coef = coef + eq_ptb = (reg_coef * + scipy.sparse.identity(self._nlp.n_eq_constraints(), + format='coo')) + ineq_ptb = (reg_coef * + scipy.sparse.identity(self._nlp.n_ineq_constraints(), + format='coo')) + + kkt.set_block(2, 2, eq_ptb) + kkt.set_block(3, 3, ineq_ptb) + return kkt + + def regularize_hessian(self, kkt, coef, copy_kkt=True): + if copy_kkt: + kkt = kkt.copy() + + hess = kkt.get_block(0, 0) + ptb = coef * scipy.sparse.identity(self._nlp.n_primals(), format='coo') + hess += ptb + kkt.set_block(0, 0, hess) + return kkt + + def _get_full_duals_primals_bounds(self): + full_duals_primals_lb = None + full_duals_primals_ub = None + # Check in case _nlp was constructed as an AmplNLP (from an nl file) + if (hasattr(self._nlp, 'pyomo_model') and + hasattr(self._nlp, 'get_pyomo_variables')): + pyomo_model = self._nlp.pyomo_model() + pyomo_variables = self._nlp.get_pyomo_variables() + if hasattr(pyomo_model,'ipopt_zL_out'): + zL_suffix = pyomo_model.ipopt_zL_out + full_duals_primals_lb = np.empty(self._nlp.n_primals()) + for i,v in enumerate(pyomo_variables): + if v in zL_suffix: + full_duals_primals_lb[i] = zL_suffix[v] + + if hasattr(pyomo_model,'ipopt_zU_out'): + zU_suffix = pyomo_model.ipopt_zU_out + full_duals_primals_ub = np.empty(self._nlp.n_primals()) + for i,v in enumerate(pyomo_variables): + if v in zU_suffix: + full_duals_primals_ub[i] = zU_suffix[v] + + if full_duals_primals_lb is None: + full_duals_primals_lb = np.ones(self._nlp.n_primals()) + + if full_duals_primals_ub is None: + full_duals_primals_ub = np.ones(self._nlp.n_primals()) + + return full_duals_primals_lb, full_duals_primals_ub + + def load_primals_into_pyomo_model(self): + if not isinstance(self._nlp, pyomo_nlp.PyomoNLP): + raise RuntimeError('Can only load primals into a pyomo model if a pyomo model was used in the constructor.') + + pyomo_variables = self._nlp.get_pyomo_variables() + primals = self._nlp.get_primals() + for i, v in enumerate(pyomo_variables): + v.value = primals[i] + + def pyomo_model(self): + return self._nlp.pyomo_model() + + def get_pyomo_variables(self): + return self._nlp.get_pyomo_variables() + + def get_pyomo_constraints(self): + return self._nlp.get_pyomo_constraints() + + def variable_names(self): + return self._nlp.variable_names() + + def constraint_names(self): + return self._nlp.constraint_names() + + def get_primal_indices(self, pyomo_variables): + return self._nlp.get_primal_indices(pyomo_variables) + + def get_constraint_indices(self, pyomo_constraints): + return self._nlp.get_constraint_indices(pyomo_constraints)