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""" 11This module defines the classes that provide an NLP interface based on 12the Ampl Solver Library (ASL) implementation 13""" 14 15import os 16import numpy as np 17import logging 18 19from scipy.sparse import coo_matrix, identity 20from pyomo.common.deprecation import deprecated 21import pyomo.core.base as pyo 22from pyomo.common.collections import ComponentMap 23from pyomo.contrib.pynumero.sparse.block_matrix import BlockMatrix 24from pyomo.contrib.pynumero.sparse.block_vector import BlockVector 25from pyomo.contrib.pynumero.interfaces.nlp import NLP 26from pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP 27from pyomo.contrib.pynumero.interfaces.utils import make_lower_triangular_full, CondensedSparseSummation 28from pyomo.contrib.pynumero.interfaces.external_grey_box import ExternalGreyBoxBlock 29from pyomo.contrib.pynumero.interfaces.nlp_projections import ProjectedNLP 30 31# Todo: make some of the numpy arise not writable from __init__ 32class PyomoNLPWithGreyBoxBlocks(NLP): 33 def __init__(self, pyomo_model): 34 super(PyomoNLPWithGreyBoxBlocks,self).__init__() 35 36 # get the list of all grey box blocks and build _ExternalGreyBoxAsNLP objects 37 greybox_components = [] 38 # build a map from the names to the variable data objects 39 # this is done over *all* variables in active blocks, even 40 # if they are not included in this model 41 self._pyomo_model_var_names_to_datas = None 42 try: 43 # We support Pynumero's ExternalGreyBoxBlock modeling 44 # objects that are provided through ExternalGreyBoxBlock objects 45 # We reclassify these as Pyomo Block objects before building the 46 # PyomoNLP object to expose any variables on the block to 47 # the underlying Pyomo machinery 48 for greybox in pyomo_model.component_objects( 49 ExternalGreyBoxBlock, descend_into=True): 50 greybox.parent_block().reclassify_component_type( 51 greybox, pyo.Block) 52 greybox_components.append(greybox) 53 54 # store the pyomo model 55 self._pyomo_model = pyomo_model 56 # build a PyomoNLP object (will include the "pyomo" 57 # part of the model only) 58 self._pyomo_nlp = PyomoNLP(pyomo_model) 59 self._pyomo_model_var_names_to_datas = \ 60 {v.getname(fully_qualified=True):v for v in pyomo_model.component_data_objects(ctype=pyo.Var, descend_into=True)} 61 self._pyomo_model_constraint_names_to_datas = \ 62 {c.getname(fully_qualified=True):c for c in pyomo_model.component_data_objects(ctype=pyo.Constraint, descend_into=True)} 63 64 finally: 65 # Restore the ctypes of the ExternalGreyBoxBlock components 66 for greybox in greybox_components: 67 greybox.parent_block().reclassify_component_type( 68 greybox, ExternalGreyBoxBlock) 69 70 if self._pyomo_nlp.n_primals() == 0: 71 raise ValueError( 72 "No variables were found in the Pyomo part of the model." 73 " PyomoGreyBoxModel requires at least one variable" 74 " to be active in a Pyomo objective or constraint") 75 76 # build the list of NLP wrappers for the greybox objects 77 greybox_nlps = [] 78 fixed_vars = [] 79 for greybox in greybox_components: 80 # iterate through the data objects if component is indexed 81 for data in greybox.values(): 82 if data.active: 83 # check that no variables are fixed 84 fixed_vars.extend(v for v in data.inputs.values() if v.fixed) 85 fixed_vars.extend(v for v in data.outputs.values() if v.fixed) 86 greybox_nlp = _ExternalGreyBoxAsNLP(data) 87 greybox_nlps.append(greybox_nlp) 88 89 if fixed_vars: 90 logging.getLogger(__name__).error('PyomoNLPWithGreyBoxBlocks found fixed variables for the' 91 ' inputs and/or outputs of an ExternalGreyBoxBlock. This' 92 ' is not currently supported. The fixed variables were:\n\t' 93 + '\n\t'.join(f.getname(fully_qualified=True) for f in fixed_vars) 94 ) 95 raise NotImplementedError('PyomoNLPWithGreyBoxBlocks does not support fixed inputs or outputs') 96 97 # let's build up the union of all the primal variables names 98 # RBP: Why use names here? Why not just ComponentSet of all 99 # data objects? 100 primals_names = set(self._pyomo_nlp.primals_names()) 101 for gbnlp in greybox_nlps: 102 primals_names.update(gbnlp.primals_names()) 103 104 # sort the names for consistency run to run 105 self._n_primals = len(primals_names) 106 self._primals_names = primals_names = sorted(primals_names) 107 self._pyomo_model_var_datas = [self._pyomo_model_var_names_to_datas[nm] for nm in self._primals_names] 108 109 # get the names of all the constraints 110 self._constraint_names = list(self._pyomo_nlp.constraint_names()) 111 self._constraint_datas = [self._pyomo_model_constraint_names_to_datas.get(nm) for nm in self._constraint_names] 112 for gbnlp in greybox_nlps: 113 self._constraint_names.extend(gbnlp.constraint_names()) 114 self._constraint_datas.extend([(gbnlp._block, nm) for nm in gbnlp.constraint_names()]) 115 self._n_constraints = len(self._constraint_names) 116 117 self._has_hessian_support = True 118 for nlp in greybox_nlps: 119 if not nlp.has_hessian_support(): 120 self._has_hessian_support = False 121 122 # wrap all the nlp objects with projected nlp objects 123 self._pyomo_nlp = ProjectedNLP(self._pyomo_nlp, primals_names) 124 for i,gbnlp in enumerate(greybox_nlps): 125 greybox_nlps[i] = ProjectedNLP(greybox_nlps[i], primals_names) 126 127 # build a list of all the nlps in order 128 self._nlps = nlps = [self._pyomo_nlp] 129 nlps.extend(greybox_nlps) 130 131 # build the primal and dual inits and lb, ub vectors 132 self._init_primals = self._pyomo_nlp.init_primals() 133 self._primals_lb = self._pyomo_nlp.primals_lb() 134 self._primals_ub = self._pyomo_nlp.primals_ub() 135 for gbnlp in greybox_nlps: 136 local = gbnlp.init_primals() 137 mask = ~np.isnan(local) 138 self._init_primals[mask] = local[mask] 139 140 local = gbnlp.primals_lb() 141 mask = ~np.isnan(local) 142 self._primals_lb[mask] = np.maximum(self._primals_lb[mask], local[mask]) 143 144 local = gbnlp.primals_ub() 145 mask = ~np.isnan(local) 146 self._primals_ub[mask] = np.minimum(self._primals_ub[mask], local[mask]) 147 148 # all the nan's should be gone (every primal should be initialized) 149 if np.any(np.isnan(self._init_primals)) \ 150 or np.any(np.isnan(self._primals_lb)) \ 151 or np.any(np.isnan(self._primals_ub)): 152 raise ValueError('NaN values found in initialization of primals or' 153 ' primals_lb or primals_ub in _PyomoNLPWithGreyBoxBlocks.') 154 155 self._init_duals = BlockVector(len(nlps)) 156 self._dual_values_blockvector = BlockVector(len(nlps)) 157 self._constraints_lb = BlockVector(len(nlps)) 158 self._constraints_ub = BlockVector(len(nlps)) 159 for i,nlp in enumerate(nlps): 160 self._init_duals.set_block(i, nlp.init_duals()) 161 self._constraints_lb.set_block(i, nlp.constraints_lb()) 162 self._constraints_ub.set_block(i, nlp.constraints_ub()) 163 self._dual_values_blockvector.set_block(i, np.nan*np.zeros(nlp.n_constraints())) 164 self._init_duals = self._init_duals.flatten() 165 self._constraints_lb = self._constraints_lb.flatten() 166 self._constraints_ub = self._constraints_ub.flatten() 167 # verify that there are no nans in the init_duals 168 if np.any(np.isnan(self._init_duals)) \ 169 or np.any(np.isnan(self._constraints_lb)) \ 170 or np.any(np.isnan(self._constraints_ub)): 171 raise ValueError('NaN values found in initialization of duals or' 172 ' constraints_lb or constraints_ub in' 173 ' _PyomoNLPWithGreyBoxBlocks.') 174 175 self._primal_values = np.nan*np.ones(self._n_primals) 176 # set the values of the primals and duals to make sure initial 177 # values get all the way through to the underlying models 178 self.set_primals(self._init_primals) 179 self.set_duals(self._init_duals) 180 assert not np.any(np.isnan(self._primal_values)) 181 assert not np.any(np.isnan(self._dual_values_blockvector)) 182 183 # if any of the problem is scaled (i.e., one or more of primals, 184 # constraints, or objective), then we want scaling factors for 185 # all of them (defaulted to 1) 186 need_scaling = False 187 # objective is owned by self._pyomo_nlp, not in any of the greybox models 188 self._obj_scaling = self._pyomo_nlp.get_obj_scaling() 189 if self._obj_scaling is None: 190 self._obj_scaling = 1.0 191 else: 192 need_scaling = True 193 194 self._primals_scaling = np.ones(self.n_primals()) 195 scaling_suffix = pyomo_model.component('scaling_factor') 196 if scaling_suffix and scaling_suffix.ctype is pyo.Suffix: 197 need_scaling = True 198 for i,v in enumerate(self._pyomo_model_var_datas): 199 if v in scaling_suffix: 200 self._primals_scaling[i] = scaling_suffix[v] 201 202 self._constraints_scaling = BlockVector(len(nlps)) 203 for i,nlp in enumerate(nlps): 204 local_constraints_scaling = nlp.get_constraints_scaling() 205 if local_constraints_scaling is None: 206 self._constraints_scaling.set_block(i, np.ones(nlp.n_constraints())) 207 else: 208 self._constraints_scaling.set_block(i, local_constraints_scaling) 209 need_scaling = True 210 if need_scaling: 211 self._constraints_scaling = self._constraints_scaling.flatten() 212 else: 213 self._obj_scaling = None 214 self._primals_scaling = None 215 self._constraints_scaling = None 216 217 # compute the jacobian and the hessian to get nnz 218 jac = self.evaluate_jacobian() 219 self._nnz_jacobian = len(jac.data) 220 221 self._sparse_hessian_summation = None 222 self._nnz_hessian_lag = None 223 if self._has_hessian_support: 224 hess = self.evaluate_hessian_lag() 225 self._nnz_hessian_lag = len(hess.data) 226 227 # overloaded from NLP 228 def n_primals(self): 229 return self._n_primals 230 231 # overloaded from NLP 232 def primals_names(self): 233 return self._primals_names 234 235 # overloaded from NLP 236 def n_constraints(self): 237 return self._n_constraints 238 239 # overloaded from NLP 240 def constraint_names(self): 241 return self._constraint_names 242 243 # overloaded from NLP 244 def nnz_jacobian(self): 245 return self._nnz_jacobian 246 247 # overloaded from NLP 248 def nnz_hessian_lag(self): 249 return self._nnz_hessian_lag 250 251 # overloaded from NLP 252 def primals_lb(self): 253 return self._primals_lb 254 255 # overloaded from NLP 256 def primals_ub(self): 257 return self._primals_ub 258 259 # overloaded from NLP 260 def constraints_lb(self): 261 return self._constraints_lb 262 263 # overloaded from NLP 264 def constraints_ub(self): 265 return self._constraints_ub 266 267 # overloaded from NLP 268 def init_primals(self): 269 return self._init_primals 270 271 # overloaded from NLP 272 def init_duals(self): 273 return self._init_duals 274 275 # overloaded from NLP / Extended NLP 276 def create_new_vector(self, vector_type): 277 if vector_type == 'primals': 278 return np.zeros(self.n_primals(), dtype=np.float64) 279 elif vector_type == 'constraints' or vector_type == 'duals': 280 return np.zeros(self.n_constraints(), dtype=np.float64) 281 else: 282 raise RuntimeError('Called create_new_vector with an unknown vector_type') 283 284 # overloaded from NLP 285 def set_primals(self, primals): 286 np.copyto(self._primal_values, primals) 287 for nlp in self._nlps: 288 nlp.set_primals(primals) 289 290 # overloaded from AslNLP 291 def get_primals(self): 292 return np.copy(self._primal_values) 293 294 # overloaded from NLP 295 def set_duals(self, duals): 296 self._dual_values_blockvector.copyfrom(duals) 297 for i,nlp in enumerate(self._nlps): 298 nlp.set_duals(self._dual_values_blockvector.get_block(i)) 299 300 # overloaded from NLP 301 def get_duals(self): 302 return self._dual_values_blockvector.flatten() 303 304 # overloaded from NLP 305 def set_obj_factor(self, obj_factor): 306 # objective is owned by the pyomo model 307 self._pyomo_nlp.set_obj_factor(obj_factor) 308 309 # overloaded from NLP 310 def get_obj_factor(self): 311 # objective is owned by the pyomo model 312 return self._pyomo_nlp.get_obj_factor() 313 314 # overloaded from NLP 315 def get_obj_scaling(self): 316 return self._obj_scaling 317 318 # overloaded from NLP 319 def get_primals_scaling(self): 320 return self._primals_scaling 321 322 # overloaded from NLP 323 def get_constraints_scaling(self): 324 return self._constraints_scaling 325 326 # overloaded from NLP 327 def evaluate_objective(self): 328 # objective is owned by the pyomo model 329 return self._pyomo_nlp.evaluate_objective() 330 331 # overloaded from NLP 332 def evaluate_grad_objective(self, out=None): 333 return self._pyomo_nlp.evaluate_grad_objective(out=out) 334 335 # overloaded from NLP 336 def evaluate_constraints(self, out=None): 337 # todo: implement the "out" version more efficiently 338 ret = BlockVector(len(self._nlps)) 339 for i,nlp in enumerate(self._nlps): 340 ret.set_block(i, nlp.evaluate_constraints()) 341 342 if out is not None: 343 ret.copyto(out) 344 return out 345 346 return ret.flatten() 347 348 # overloaded from NLP 349 def evaluate_jacobian(self, out=None): 350 ret = BlockMatrix(len(self._nlps),1) 351 for i,nlp in enumerate(self._nlps): 352 ret.set_block(i, 0, nlp.evaluate_jacobian()) 353 ret = ret.tocoo() 354 355 if out is not None: 356 assert np.array_equal(ret.row, out.row) 357 assert np.array_equal(ret.col, out.col) 358 np.copyto(out.data, ret.data) 359 return out 360 return ret 361 362 def evaluate_hessian_lag(self, out=None): 363 list_of_hessians = [nlp.evaluate_hessian_lag() for nlp in self._nlps] 364 if self._sparse_hessian_summation is None: 365 # This is assuming that the nonzero structures of Hessians 366 # do not change 367 self._sparse_hessian_summation = CondensedSparseSummation(list_of_hessians) 368 ret = self._sparse_hessian_summation.sum(list_of_hessians) 369 370 if out is not None: 371 assert np.array_equal(ret.row, out.row) 372 assert np.array_equal(ret.col, out.col) 373 np.copyto(out.data, ret.data) 374 return out 375 return ret 376 377 def report_solver_status(self, status_code, status_message): 378 raise NotImplementedError('This is not yet implemented.') 379 380 def load_state_into_pyomo(self, bound_multipliers=None): 381 # load the values of the primals into the pyomo 382 primals = self.get_primals() 383 for value,vardata in zip(primals, self._pyomo_model_var_datas): 384 vardata.set_value(value) 385 386 # get the active suffixes 387 m = self._pyomo_model 388 model_suffixes = dict( 389 pyo.suffix.active_import_suffix_generator(m)) 390 391 # we need to correct the sign of the multipliers based on whether or 392 # not we are minimizing or maximizing - this is done in the ASL interface 393 # for ipopt, but does not appear to be done in cyipopt. 394 obj_sign = 1.0 395 objs = list(m.component_data_objects(ctype=pyo.Objective, descend_into=True)) 396 assert len(objs) == 1 397 if objs[0].sense == pyo.maximize: 398 obj_sign = -1.0 399 400 if 'dual' in model_suffixes: 401 model_suffixes['dual'].clear() 402 dual_values = self._dual_values_blockvector.flatten() 403 for value,t in zip(dual_values, self._constraint_datas): 404 if type(t) is tuple: 405 model_suffixes['dual'].setdefault(t[0], {})[t[1]] = -obj_sign*value 406 else: 407 # t is a constraint data 408 model_suffixes['dual'][t] = -obj_sign*value 409 410 if 'ipopt_zL_out' in model_suffixes: 411 model_suffixes['ipopt_zL_out'].clear() 412 if bound_multipliers is not None: 413 model_suffixes['ipopt_zL_out'].update( 414 zip(self._pyomo_model_var_datas, obj_sign*bound_multipliers[0])) 415 if 'ipopt_zU_out' in model_suffixes: 416 model_suffixes['ipopt_zU_out'].clear() 417 if bound_multipliers is not None: 418 model_suffixes['ipopt_zU_out'].update( 419 zip(self._pyomo_model_var_datas, -obj_sign*bound_multipliers[1])) 420 421def _default_if_none(value, default): 422 if value is None: 423 return default 424 return value 425 426class _ExternalGreyBoxAsNLP(NLP): 427 """ 428 This class takes an ExternalGreyBoxModel and makes it look 429 like an NLP so it can be used with other interfaces. Currently, 430 the ExternalGreyBoxModel supports constraints only (no objective), 431 so some of the methods are not appropriate and raise exceptions 432 """ 433 def __init__(self, external_grey_box_block): 434 self._block = external_grey_box_block 435 self._ex_model = external_grey_box_block.get_external_model() 436 n_inputs = len(self._block.inputs) 437 assert n_inputs == self._ex_model.n_inputs() 438 n_eq_constraints = self._ex_model.n_equality_constraints() 439 n_outputs = len(self._block.outputs) 440 assert n_outputs == self._ex_model.n_outputs() 441 442 if self._ex_model.n_outputs() == 0 and \ 443 self._ex_model.n_equality_constraints() == 0: 444 raise ValueError( 445 'ExternalGreyBoxModel has no equality constraints ' 446 'or outputs. To use _ExternalGreyBoxAsNLP, it must' 447 ' have at least one or both.') 448 449 # create the list of primals and constraint names 450 # primals will be ordered inputs, followed by outputs 451 self._primals_names = \ 452 [self._block.inputs[k].getname(fully_qualified=True) \ 453 for k in self._block.inputs] 454 self._primals_names.extend( 455 [self._block.outputs[k].getname(fully_qualified=True) \ 456 for k in self._block.outputs] 457 ) 458 n_primals = len(self._primals_names) 459 460 prefix = self._block.getname(fully_qualified=True) 461 self._constraint_names = \ 462 ['{}.{}'.format(prefix, nm) \ 463 for nm in self._ex_model.equality_constraint_names()] 464 output_var_names = \ 465 [self._block.outputs[k].getname(fully_qualified=False) \ 466 for k in self._block.outputs] 467 self._constraint_names.extend( 468 ['{}.output_constraints[{}]'.format(prefix, nm) \ 469 for nm in self._ex_model.output_names()]) 470 471 # create the numpy arrays of bounds on the primals 472 self._primals_lb = BlockVector(2) 473 self._primals_ub = BlockVector(2) 474 self._init_primals = BlockVector(2) 475 lb = np.nan*np.zeros(n_inputs) 476 ub = np.nan*np.zeros(n_inputs) 477 init_primals = np.nan*np.zeros(n_inputs) 478 for i,k in enumerate(self._block.inputs): 479 lb[i] = _default_if_none(self._block.inputs[k].lb, -np.inf) 480 ub[i] = _default_if_none(self._block.inputs[k].ub, np.inf) 481 init_primals[i] = _default_if_none(self._block.inputs[k].value, 0.0) 482 self._primals_lb.set_block(0,lb) 483 self._primals_ub.set_block(0,ub) 484 self._init_primals.set_block(0, init_primals) 485 486 lb = np.nan*np.zeros(n_outputs) 487 ub = np.nan*np.zeros(n_outputs) 488 init_primals = np.nan*np.zeros(n_outputs) 489 for i,k in enumerate(self._block.outputs): 490 lb[i] = _default_if_none(self._block.outputs[k].lb, -np.inf) 491 ub[i] = _default_if_none(self._block.outputs[k].ub, np.inf) 492 init_primals[i] = _default_if_none(self._block.outputs[k].value, 0.0) 493 self._primals_lb.set_block(1,lb) 494 self._primals_ub.set_block(1,ub) 495 self._init_primals.set_block(1,init_primals) 496 self._primals_lb = self._primals_lb.flatten() 497 self._primals_ub = self._primals_ub.flatten() 498 self._init_primals = self._init_primals.flatten() 499 500 # create a numpy array to store the values of the primals 501 self._primal_values = np.copy(self._init_primals) 502 # make sure the values are passed through to other objects 503 self.set_primals(self._init_primals) 504 505 # create the numpy arrays for the duals and initial values 506 # for now, initialize the duals to zero 507 self._init_duals = np.zeros(self.n_constraints(), dtype=np.float64) 508 # create the numpy arrays to store the dual variables 509 self._dual_values = np.copy(self._init_duals) 510 # make sure the values are passed through to other objects 511 self.set_duals(self._init_duals) 512 513 # create the numpy arrays for bounds on the constraints 514 # for now all of these are equalities 515 self._constraints_lb = np.zeros(self.n_constraints(), dtype=np.float64) 516 self._constraints_ub = np.zeros(self.n_constraints(), dtype=np.float64) 517 518 # do we have hessian support 519 self._has_hessian_support = True 520 if self._ex_model.n_equality_constraints() > 0 \ 521 and not hasattr(self._ex_model, 'evaluate_hessian_equality_constraints'): 522 self._has_hessian_support = False 523 if self._ex_model.n_outputs() > 0 \ 524 and not hasattr(self._ex_model, 'evaluate_hessian_outputs'): 525 self._has_hessian_support = False 526 527 self._nnz_jacobian = None 528 self._nnz_hessian_lag = None 529 self._cached_constraint_residuals = None 530 self._cached_jacobian = None 531 self._cached_hessian = None 532 533 def n_primals(self): 534 return len(self._primals_names) 535 536 def primals_names(self): 537 return list(self._primals_names) 538 539 def n_constraints(self): 540 return len(self._constraint_names) 541 542 def constraint_names(self): 543 return list(self._constraint_names) 544 545 def nnz_jacobian(self): 546 if self._nnz_jacobian is None: 547 J = self.evaluate_jacobian() 548 self._nnz_jacobian = len(J.data) 549 return self._nnz_jacobian 550 551 def nnz_hessian_lag(self): 552 if self._nnz_hessian_lag is None: 553 H = self.evaluate_hessian_lag() 554 self._nnz_hessian_lag = len(H.data) 555 return self._nnz_hessian_lag 556 557 def primals_lb(self): 558 return np.copy(self._primals_lb) 559 560 def primals_ub(self): 561 return np.copy(self._primals_ub) 562 563 def constraints_lb(self): 564 return np.copy(self._constraints_lb) 565 566 def constraints_ub(self): 567 return np.copy(self._constraints_ub) 568 569 def init_primals(self): 570 return np.copy(self._init_primals) 571 572 def init_duals(self): 573 return np.copy(self._init_duals) 574 575 def create_new_vector(self, vector_type): 576 if vector_type == 'primals': 577 return np.zeros(self.n_primals(), dtype=np.float64) 578 elif vector_type == 'constraints' or vector_type == 'duals': 579 return np.zeros(self.n_constraints(), dtype=np.float64) 580 581 def _cache_invalidate_primals(self): 582 self._cached_constraint_residuals = None 583 self._cached_jacobian = None 584 self._cached_hessian = None 585 586 def set_primals(self, primals): 587 self._cache_invalidate_primals() 588 assert len(primals) == self.n_primals() 589 np.copyto(self._primal_values, primals) 590 self._ex_model.set_input_values(primals[:self._ex_model.n_inputs()]) 591 592 def get_primals(self): 593 return np.copy(self._primal_values) 594 595 def _cache_invalidate_duals(self): 596 self._cached_hessian = None 597 598 def set_duals(self, duals): 599 self._cache_invalidate_duals() 600 assert len(duals) == self.n_constraints() 601 np.copyto(self._dual_values, duals) 602 if self._ex_model.n_equality_constraints() > 0: 603 self._ex_model.set_equality_constraint_multipliers( 604 self._dual_values[:self._ex_model.n_equality_constraints()] 605 ) 606 if self._ex_model.n_outputs() > 0: 607 self._ex_model.set_output_constraint_multipliers( 608 self._dual_values[self._ex_model.n_equality_constraints():] 609 ) 610 611 def get_duals(self): 612 return np.copy(self._dual_values) 613 614 def set_obj_factor(self, obj_factor): 615 raise NotImplementedError('_ExternalGreyBoxAsNLP does not support objectives') 616 617 def get_obj_factor(self): 618 raise NotImplementedError('_ExternalGreyBoxAsNLP does not support objectives') 619 620 def get_obj_scaling(self): 621 raise NotImplementedError('_ExternalGreyBoxAsNLP does not support objectives') 622 623 def get_primals_scaling(self): 624 raise NotImplementedError( 625 '_ExternalGreyBoxAsNLP does not support scaling of primals ' 626 'directly. This should be handled at a higher level using ' 627 'suffixes on the Pyomo variables.' 628 ) 629 630 def get_constraints_scaling(self): 631 # todo: would this be better with block vectors 632 scaling = np.ones(self.n_constraints(), dtype=np.float64) 633 scaled = False 634 if self._ex_model.n_equality_constraints() > 0: 635 eq_scaling = self._ex_model.get_equality_constraint_scaling_factors() 636 if eq_scaling is not None: 637 scaling[:self._ex_model.n_equality_constraints()] = eq_scaling 638 scaled = True 639 if self._ex_model.n_outputs() > 0: 640 output_scaling = self._ex_model.get_output_constraint_scaling_factors() 641 if output_scaling is not None: 642 scaling[self._ex_model.n_equality_constraints():] = output_scaling 643 scaled = True 644 if scaled: 645 return scaling 646 return None 647 648 def evaluate_objective(self): 649 # todo: Should we return 0 here? 650 raise NotImplementedError('_ExternalGreyBoxNLP does not support objectives') 651 652 def evaluate_grad_objective(self, out=None): 653 # todo: Should we return 0 here? 654 raise NotImplementedError('_ExternalGreyBoxNLP does not support objectives') 655 656 def _evaluate_constraints_if_necessary_and_cache(self): 657 if self._cached_constraint_residuals is None: 658 c = BlockVector(2) 659 if self._ex_model.n_equality_constraints() > 0: 660 c.set_block(0, self._ex_model.evaluate_equality_constraints()) 661 else: 662 c.set_block(0, np.zeros(0, dtype=np.float64)) 663 if self._ex_model.n_outputs() > 0: 664 output_values = self._primal_values[self._ex_model.n_inputs():] 665 c.set_block(1, self._ex_model.evaluate_outputs() - output_values) 666 else: 667 c.set_block(1,np.zeros(0, dtype=np.float64)) 668 self._cached_constraint_residuals = c.flatten() 669 670 def evaluate_constraints(self, out=None): 671 self._evaluate_constraints_if_necessary_and_cache() 672 if out is not None: 673 assert len(out) == self.n_constraints() 674 np.copyto(out, self._cached_constraint_residuals) 675 return out 676 677 return np.copy(self._cached_constraint_residuals) 678 679 def _evaluate_jacobian_if_necessary_and_cache(self): 680 if self._cached_jacobian is None: 681 jac = BlockMatrix(2,2) 682 jac.set_row_size(0,self._ex_model.n_equality_constraints()) 683 jac.set_row_size(1,self._ex_model.n_outputs()) 684 jac.set_col_size(0,self._ex_model.n_inputs()) 685 jac.set_col_size(1,self._ex_model.n_outputs()) 686 687 if self._ex_model.n_equality_constraints() > 0: 688 jac.set_block(0,0,self._ex_model.evaluate_jacobian_equality_constraints()) 689 if self._ex_model.n_outputs() > 0: 690 jac.set_block(1,0,self._ex_model.evaluate_jacobian_outputs()) 691 jac.set_block(1,1,-1.0*identity(self._ex_model.n_outputs())) 692 693 self._cached_jacobian = jac.tocoo() 694 695 def evaluate_jacobian(self, out=None): 696 self._evaluate_jacobian_if_necessary_and_cache() 697 if out is not None: 698 jac = self._cached_jacobian 699 assert np.array_equal(jac.row, out.row) 700 assert np.array_equal(jac.col, out.col) 701 np.copyto(out.data, jac.data) 702 return out 703 704 return self._cached_jacobian.copy() 705 706 def _evaluate_hessian_if_necessary_and_cache(self): 707 if self._cached_hessian is None: 708 hess = BlockMatrix(2,2) 709 hess.set_row_size(0,self._ex_model.n_inputs()) 710 hess.set_row_size(1,self._ex_model.n_outputs()) 711 hess.set_col_size(0,self._ex_model.n_inputs()) 712 hess.set_col_size(1,self._ex_model.n_outputs()) 713 714 # get the hessian w.r.t. the equality constraints 715 eq_hess = None 716 if self._ex_model.n_equality_constraints() > 0: 717 eq_hess = self._ex_model.evaluate_hessian_equality_constraints() 718 # let's check that it is lower triangular 719 if np.any(eq_hess.row < eq_hess.col): 720 raise ValueError('ExternalGreyBoxModel must return lower ' 721 'triangular portion of the Hessian only') 722 723 eq_hess = make_lower_triangular_full(eq_hess) 724 725 output_hess = None 726 if self._ex_model.n_outputs() > 0: 727 output_hess = self._ex_model.evaluate_hessian_outputs() 728 # let's check that it is lower triangular 729 if np.any(output_hess.row < output_hess.col): 730 raise ValueError('ExternalGreyBoxModel must return lower ' 731 'triangular portion of the Hessian only') 732 733 output_hess = make_lower_triangular_full(output_hess) 734 735 input_hess = None 736 if eq_hess is not None and output_hess is not None: 737 # we may want to make this more efficient 738 row = np.concatenate((eq_hess.row, output_hess.row)) 739 col = np.concatenate((eq_hess.col, output_hess.col)) 740 data = np.concatenate((eq_hess.data, output_hess.data)) 741 742 assert eq_hess.shape == output_hess.shape 743 input_hess = coo_matrix( (data, (row,col)), shape=eq_hess.shape) 744 elif eq_hess is not None: 745 input_hess = eq_hess 746 elif output_hess is not None: 747 input_hess = output_hess 748 assert input_hess is not None # need equality or outputs or both 749 750 hess.set_block(0,0,input_hess) 751 self._cached_hessian = hess.tocoo() 752 753 def has_hessian_support(self): 754 return self._has_hessian_support 755 756 def evaluate_hessian_lag(self, out=None): 757 if not self._has_hessian_support: 758 raise NotImplementedError( 759 'Hessians not supported for all of the external grey box' 760 ' models. Therefore, Hessians are not supported overall.' 761 ) 762 763 self._evaluate_hessian_if_necessary_and_cache() 764 if out is not None: 765 hess = self._cached_hessian 766 assert np.array_equal(hess.row, out.row) 767 assert np.array_equal(hess.col, out.col) 768 np.copyto(out.data, hess.data) 769 return out 770 771 return self._cached_hessian.copy() 772 773 def report_solver_status(self, status_code, status_message): 774 raise NotImplementedError('report_solver_status not implemented') 775