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 17 18from scipy.sparse import coo_matrix 19from pyomo.common.deprecation import deprecated 20from pyomo.common.tempfiles import TempfileManager 21from pyomo.opt import WriterFactory 22import pyomo.core.base as pyo 23from pyomo.common.collections import ComponentMap 24from pyomo.common.env import CtypesEnviron 25from ..sparse.block_matrix import BlockMatrix 26from pyomo.contrib.pynumero.interfaces.ampl_nlp import AslNLP 27from pyomo.contrib.pynumero.interfaces.nlp import NLP 28from .external_grey_box import ExternalGreyBoxBlock 29 30 31__all__ = ['PyomoNLP'] 32 33# TODO: There are todos in the code below 34class PyomoNLP(AslNLP): 35 def __init__(self, pyomo_model, nl_file_options=None): 36 """ 37 Pyomo nonlinear program interface 38 39 Parameters 40 ---------- 41 pyomo_model: pyomo.environ.ConcreteModel 42 Pyomo concrete model 43 """ 44 TempfileManager.push() 45 try: 46 # get the temp file names for the nl file 47 nl_file = TempfileManager.create_tempfile( 48 suffix='pynumero.nl') 49 50 # The current AmplInterface code only supports a single 51 # objective function Therefore, we throw an error if there 52 # is not one (and only one) active objective function. This 53 # is better than adding a dummy objective that the user does 54 # not know about (since we do not have a good place to 55 # remove this objective later) 56 # 57 # TODO: extend the AmplInterface and the AslNLP to correctly 58 # handle this 59 # 60 # This currently addresses issue #1217 61 objectives = list(pyomo_model.component_data_objects( 62 ctype=pyo.Objective, active=True, descend_into=True)) 63 if len(objectives) != 1: 64 raise NotImplementedError( 65 'The ASL interface and PyomoNLP in PyNumero currently ' 66 'only support single objective problems. Deactivate ' 67 'any extra objectives you may have, or add a dummy ' 68 'objective (f(x)=0) if you have a square problem ' 69 '(found %s objectives).' % (len(objectives),)) 70 self._objective = objectives[0] 71 72 # write the nl file for the Pyomo model and get the symbolMap 73 if nl_file_options is None: 74 nl_file_options = dict() 75 fname, symbolMap = WriterFactory('nl')( 76 pyomo_model, nl_file, lambda x:True, nl_file_options) 77 78 # create component maps from vardata to idx and condata to idx 79 self._vardata_to_idx = vdidx = ComponentMap() 80 self._condata_to_idx = cdidx = ComponentMap() 81 82 # TODO: Are these names totally consistent? 83 for name, obj in symbolMap.bySymbol.items(): 84 if name[0] == 'v': 85 vdidx[obj()] = int(name[1:]) 86 elif name[0] == 'c': 87 cdidx[obj()] = int(name[1:]) 88 89 # The NL writer advertises the external function libraries 90 # through the PYOMO_AMPLFUNC environment variable; merge it 91 # with any preexisting AMPLFUNC definitions 92 amplfunc = "\n".join( 93 val for val in ( 94 os.environ.get('AMPLFUNC', ''), 95 os.environ.get('PYOMO_AMPLFUNC', ''), 96 ) if val) 97 with CtypesEnviron(AMPLFUNC=amplfunc): 98 super(PyomoNLP, self).__init__(nl_file) 99 100 # keep pyomo model in cache 101 self._pyomo_model = pyomo_model 102 103 # Create ComponentMap corresponding to equality constraint indices 104 # This must be done after the call to super-init. 105 full_to_equality = self._con_full_eq_map 106 equality_mask = self._con_full_eq_mask 107 self._condata_to_eq_idx = ComponentMap( 108 (con, full_to_equality[i]) 109 for con, i in self._condata_to_idx.items() 110 if equality_mask[i] 111 ) 112 full_to_inequality = self._con_full_ineq_map 113 inequality_mask = self._con_full_ineq_mask 114 self._condata_to_ineq_idx = ComponentMap( 115 (con, full_to_inequality[i]) 116 for con, i in self._condata_to_idx.items() 117 if inequality_mask[i] 118 ) 119 120 finally: 121 # delete the nl file 122 TempfileManager.pop() 123 124 125 def pyomo_model(self): 126 """ 127 Return optimization model 128 """ 129 return self._pyomo_model 130 131 def get_pyomo_objective(self): 132 """ 133 Return an instance of the active objective function on the Pyomo model. 134 (there can be only one) 135 """ 136 return self._objective 137 138 def get_pyomo_variables(self): 139 """ 140 Return an ordered list of the Pyomo VarData objects in 141 the order corresponding to the primals 142 """ 143 # ToDo: is there a more efficient way to do this 144 idx_to_vardata = {i:v for v,i in self._vardata_to_idx.items()} 145 return [idx_to_vardata[i] for i in range(len(idx_to_vardata))] 146 147 def get_pyomo_constraints(self): 148 """ 149 Return an ordered list of the Pyomo ConData objects in 150 the order corresponding to the primals 151 """ 152 # ToDo: is there a more efficient way to do this 153 idx_to_condata = {i:v for v,i in self._condata_to_idx.items()} 154 return [idx_to_condata[i] for i in range(len(idx_to_condata))] 155 156 def get_pyomo_equality_constraints(self): 157 """ 158 Return an ordered list of the Pyomo ConData objects in 159 the order corresponding to the equality constraints. 160 """ 161 idx_to_condata = {i: c for c, i in 162 self._condata_to_eq_idx.items()} 163 return [idx_to_condata[i] for i in range(len(idx_to_condata))] 164 165 def get_pyomo_inequality_constraints(self): 166 """ 167 Return an ordered list of the Pyomo ConData objects in 168 the order corresponding to the inequality constraints. 169 """ 170 idx_to_condata = {i: c for c, i in 171 self._condata_to_ineq_idx.items()} 172 return [idx_to_condata[i] for i in range(len(idx_to_condata))] 173 174 @deprecated(msg='This method has been replaced with primals_names', version='6.0.0.dev0', remove_in='6.0') 175 def variable_names(self): 176 return self.primals_names() 177 178 def primals_names(self): 179 """ 180 Return an ordered list of the Pyomo variable 181 names in the order corresponding to the primals 182 """ 183 pyomo_variables = self.get_pyomo_variables() 184 return [v.getname(fully_qualified=True) for v in pyomo_variables] 185 186 def constraint_names(self): 187 """ 188 Return an ordered list of the Pyomo constraint 189 names in the order corresponding to internal constraint order 190 """ 191 pyomo_constraints = self.get_pyomo_constraints() 192 return [v.getname(fully_qualified=True) for v in pyomo_constraints] 193 194 def equality_constraint_names(self): 195 """ 196 Return an ordered list of the Pyomo ConData names in 197 the order corresponding to the equality constraints. 198 """ 199 equality_constraints = self.get_pyomo_equality_constraints() 200 return [v.getname(fully_qualified=True) for v in equality_constraints] 201 202 def inequality_constraint_names(self): 203 """ 204 Return an ordered list of the Pyomo ConData names in 205 the order corresponding to the inequality constraints. 206 """ 207 inequality_constraints = self.get_pyomo_inequality_constraints() 208 return [v.getname(fully_qualified=True) for v in inequality_constraints] 209 210 def get_primal_indices(self, pyomo_variables): 211 """ 212 Return the list of indices for the primals 213 corresponding to the list of Pyomo variables provided 214 215 Parameters 216 ---------- 217 pyomo_variables : list of Pyomo Var or VarData objects 218 """ 219 assert isinstance(pyomo_variables, list) 220 var_indices = [] 221 for v in pyomo_variables: 222 if v.is_indexed(): 223 for vd in v.values(): 224 var_id = self._vardata_to_idx[vd] 225 var_indices.append(var_id) 226 else: 227 var_id = self._vardata_to_idx[v] 228 var_indices.append(var_id) 229 return var_indices 230 231 def get_constraint_indices(self, pyomo_constraints): 232 """ 233 Return the list of indices for the constraints 234 corresponding to the list of Pyomo constraints provided 235 236 Parameters 237 ---------- 238 pyomo_constraints : list of Pyomo Constraint or ConstraintData objects 239 """ 240 assert isinstance(pyomo_constraints, list) 241 con_indices = [] 242 for c in pyomo_constraints: 243 if c.is_indexed(): 244 for cd in c.values(): 245 con_id = self._condata_to_idx[cd] 246 con_indices.append(con_id) 247 else: 248 con_id = self._condata_to_idx[c] 249 con_indices.append(con_id) 250 return con_indices 251 252 def get_equality_constraint_indices(self, constraints): 253 """ 254 Return the list of equality indices for the constraints 255 corresponding to the list of Pyomo constraints provided. 256 257 Parameters 258 ---------- 259 constraints : list of Pyomo Constraints or ConstraintData objects 260 """ 261 indices = [] 262 for c in constraints: 263 if c.is_indexed(): 264 for cd in c.values(): 265 con_eq_idx = self._condata_to_eq_idx[cd] 266 indices.append(con_eq_idx) 267 else: 268 con_eq_idx = self._condata_to_eq_idx[c] 269 indices.append(con_eq_idx) 270 return indices 271 272 def get_inequality_constraint_indices(self, constraints): 273 """ 274 Return the list of inequality indices for the constraints 275 corresponding to the list of Pyomo constraints provided. 276 277 Parameters 278 ---------- 279 constraints : list of Pyomo Constraints or ConstraintData objects 280 """ 281 indices = [] 282 for c in constraints: 283 if c.is_indexed(): 284 for cd in c.values(): 285 con_ineq_idx = self._condata_to_ineq_idx[cd] 286 indices.append(con_ineq_idx) 287 else: 288 con_ineq_idx = self._condata_to_ineq_idx[c] 289 indices.append(con_ineq_idx) 290 return indices 291 292 # overloaded from NLP 293 def get_obj_scaling(self): 294 obj = self.get_pyomo_objective() 295 scaling_suffix = self._pyomo_model.component('scaling_factor') 296 if scaling_suffix and scaling_suffix.ctype is pyo.Suffix: 297 if obj in scaling_suffix: 298 return scaling_suffix[obj] 299 return 1.0 300 return None 301 302 # overloaded from NLP 303 def get_primals_scaling(self): 304 scaling_suffix = self._pyomo_model.component('scaling_factor') 305 if scaling_suffix and scaling_suffix.ctype is pyo.Suffix: 306 primals_scaling = np.ones(self.n_primals()) 307 for i,v in enumerate(self.get_pyomo_variables()): 308 if v in scaling_suffix: 309 primals_scaling[i] = scaling_suffix[v] 310 return primals_scaling 311 return None 312 313 # overloaded from NLP 314 def get_constraints_scaling(self): 315 scaling_suffix = self._pyomo_model.component('scaling_factor') 316 if scaling_suffix and scaling_suffix.ctype is pyo.Suffix: 317 constraints_scaling = np.ones(self.n_constraints()) 318 for i,c in enumerate(self.get_pyomo_constraints()): 319 if c in scaling_suffix: 320 constraints_scaling[i] = scaling_suffix[c] 321 return constraints_scaling 322 return None 323 324 def extract_subvector_grad_objective(self, pyomo_variables): 325 """Compute the gradient of the objective and return the entries 326 corresponding to the given Pyomo variables 327 328 Parameters 329 ---------- 330 pyomo_variables : list of Pyomo Var or VarData objects 331 """ 332 grad_obj = self.evaluate_grad_objective() 333 return grad_obj[self.get_primal_indices(pyomo_variables)] 334 335 def extract_subvector_constraints(self, pyomo_constraints): 336 """ 337 Return the values of the constraints 338 corresponding to the list of Pyomo constraints provided 339 340 Parameters 341 ---------- 342 pyomo_constraints : list of Pyomo Constraint or ConstraintData objects 343 """ 344 residuals = self.evaluate_constraints() 345 return residuals[self.get_constraint_indices(pyomo_constraints)] 346 347 def extract_submatrix_jacobian(self, pyomo_variables, pyomo_constraints): 348 """ 349 Return the submatrix of the jacobian that corresponds to the list 350 of Pyomo variables and list of Pyomo constraints provided 351 352 Parameters 353 ---------- 354 pyomo_variables : list of Pyomo Var or VarData objects 355 pyomo_constraints : list of Pyomo Constraint or ConstraintData objects 356 """ 357 jac = self.evaluate_jacobian() 358 primal_indices = self.get_primal_indices(pyomo_variables) 359 constraint_indices = self.get_constraint_indices(pyomo_constraints) 360 row_mask = np.isin(jac.row, constraint_indices) 361 col_mask = np.isin(jac.col, primal_indices) 362 submatrix_mask = row_mask & col_mask 363 submatrix_irows = np.compress(submatrix_mask, jac.row) 364 submatrix_jcols = np.compress(submatrix_mask, jac.col) 365 submatrix_data = np.compress(submatrix_mask, jac.data) 366 367 # ToDo: this is expensive - have to think about how to do this with numpy 368 row_submatrix_map = {j:i for i,j in enumerate(constraint_indices)} 369 for i, v in enumerate(submatrix_irows): 370 submatrix_irows[i] = row_submatrix_map[v] 371 372 col_submatrix_map = {j:i for i,j in enumerate(primal_indices)} 373 for i, v in enumerate(submatrix_jcols): 374 submatrix_jcols[i] = col_submatrix_map[v] 375 376 return coo_matrix((submatrix_data, (submatrix_irows, submatrix_jcols)), shape=(len(constraint_indices), len(primal_indices))) 377 378 def extract_submatrix_hessian_lag(self, pyomo_variables_rows, pyomo_variables_cols): 379 """ 380 Return the submatrix of the hessian of the lagrangian that 381 corresponds to the list of Pyomo variables provided 382 383 Parameters 384 ---------- 385 pyomo_variables_rows : list of Pyomo Var or VarData objects 386 List of Pyomo Var or VarData objects corresponding to the desired rows 387 pyomo_variables_cols : list of Pyomo Var or VarData objects 388 List of Pyomo Var or VarData objects corresponding to the desired columns 389 """ 390 hess_lag = self.evaluate_hessian_lag() 391 primal_indices_rows = self.get_primal_indices(pyomo_variables_rows) 392 primal_indices_cols = self.get_primal_indices(pyomo_variables_cols) 393 row_mask = np.isin(hess_lag.row, primal_indices_rows) 394 col_mask = np.isin(hess_lag.col, primal_indices_cols) 395 submatrix_mask = row_mask & col_mask 396 submatrix_irows = np.compress(submatrix_mask, hess_lag.row) 397 submatrix_jcols = np.compress(submatrix_mask, hess_lag.col) 398 submatrix_data = np.compress(submatrix_mask, hess_lag.data) 399 400 # ToDo: this is expensive - have to think about how to do this with numpy 401 submatrix_map = {j:i for i,j in enumerate(primal_indices_rows)} 402 for i, v in enumerate(submatrix_irows): 403 submatrix_irows[i] = submatrix_map[v] 404 405 submatrix_map = {j:i for i,j in enumerate(primal_indices_cols)} 406 for i, v in enumerate(submatrix_jcols): 407 submatrix_jcols[i] = submatrix_map[v] 408 409 return coo_matrix((submatrix_data, (submatrix_irows, submatrix_jcols)), shape=(len(primal_indices_rows), len(primal_indices_cols))) 410 411 def load_state_into_pyomo(self, bound_multipliers=None): 412 primals = self.get_primals() 413 variables = self.get_pyomo_variables() 414 for var, val in zip(variables, primals): 415 var.set_value(val) 416 m = self.pyomo_model() 417 model_suffixes = dict( 418 pyo.suffix.active_import_suffix_generator(m)) 419 if 'dual' in model_suffixes: 420 duals = self.get_duals() 421 constraints = self.get_pyomo_constraints() 422 model_suffixes['dual'].clear() 423 model_suffixes['dual'].update( 424 zip(constraints, duals)) 425 if 'ipopt_zL_out' in model_suffixes: 426 model_suffixes['ipopt_zL_out'].clear() 427 if bound_multipliers is not None: 428 model_suffixes['ipopt_zL_out'].update( 429 zip(variables, bound_multipliers[0])) 430 if 'ipopt_zU_out' in model_suffixes: 431 model_suffixes['ipopt_zU_out'].clear() 432 if bound_multipliers is not None: 433 model_suffixes['ipopt_zU_out'].update( 434 zip(variables, bound_multipliers[1])) 435 436# TODO: look for the [:-i] when i might be zero 437class PyomoGreyBoxNLP(NLP): 438 def __init__(self, pyomo_model): 439 # store all the greybox custom block data objects 440 greybox_components = [] 441 try: 442 # We support Pynumero's ExternalGreyBoxBlock modeling 443 # objects. We need to find them and convert them to Blocks 444 # before calling the NL writer so that the attached Vars get 445 # picked up by the writer. 446 for greybox in pyomo_model.component_objects( 447 ExternalGreyBoxBlock, descend_into=True): 448 greybox.parent_block().reclassify_component_type( 449 greybox, pyo.Block) 450 greybox_components.append(greybox) 451 452 self._pyomo_model = pyomo_model 453 self._pyomo_nlp = PyomoNLP(pyomo_model) 454 455 finally: 456 # Restore the ctypes of the ExternalGreyBoxBlock components 457 for greybox in greybox_components: 458 greybox.parent_block().reclassify_component_type( 459 greybox, ExternalGreyBoxBlock) 460 461 # get the greybox block data objects 462 greybox_data = [] 463 for greybox in greybox_components: 464 greybox_data.extend(data for data in greybox.values() 465 if data.active) 466 467 if len(greybox_data) > 1: 468 raise NotImplementedError("The PyomoGreyBoxModel interface has not" 469 " been tested with Pyomo models that contain" 470 " more than one ExternalGreyBoxBlock. Currently," 471 " only a single block is supported.") 472 473 if self._pyomo_nlp.n_primals() == 0: 474 raise ValueError( 475 "No variables were found in the Pyomo part of the model." 476 " PyomoGreyBoxModel requires at least one variable" 477 " to be active in a Pyomo objective or constraint") 478 479 # check that the greybox model supports what we would expect 480 # TODO: add support for models that do not provide jacobians 481 """ 482 for data in greybox_data: 483 c = data._ex_model.model_capabilities() 484 if (c.n_equality_constraints() > 0 \ 485 and not c.supports_jacobian_equality_constraints) \ 486 or (c.n_equality_constraints() > 0 \ 487 and not c.supports_jacobian_equality_constraints) 488 raise NotImplementedError('PyomoGreyBoxNLP does not support models' 489 ' without explicit Jacobian support') 490 """ 491 492 # number of additional variables required - they are in the 493 # greybox models but not included in the NL file 494 self._n_greybox_primals = 0 495 496 # number of residuals (equality constraints + output constraints 497 # coming from the grey box models 498 self._greybox_primals_names = [] 499 self._greybox_constraints_names = [] 500 501 # Update the primal index map with any variables in the 502 # greybox models that do not otherwise appear in the NL 503 # and capture some other book keeping items 504 n_primals = self._pyomo_nlp.n_primals() 505 greybox_primals = [] 506 self._vardata_to_idx = ComponentMap(self._pyomo_nlp._vardata_to_idx) 507 for data in greybox_data: 508 # check that none of the inputs / outputs are fixed 509 for v in data.inputs.values(): 510 if v.fixed: 511 raise NotImplementedError('Found a grey box model input that is fixed: {}.' 512 ' This interface does not currently support fixed' 513 ' variables. Please add a constraint instead.' 514 ''.format(v.getname(fully_qualified=True))) 515 for v in data.outputs.values(): 516 if v.fixed: 517 raise NotImplementedError('Found a grey box model output that is fixed: {}.' 518 ' This interface does not currently support fixed' 519 ' variables. Please add a constraint instead.' 520 ''.format(v.getname(fully_qualified=True))) 521 522 block_name = data.getname() 523 for nm in data._ex_model.equality_constraint_names(): 524 self._greybox_constraints_names.append('{}.{}'.format(block_name, nm)) 525 for nm in data._ex_model.output_names(): 526 self._greybox_constraints_names.append('{}.{}_con'.format(block_name, nm)) 527 528 for var in data.component_data_objects(pyo.Var): 529 if var not in self._vardata_to_idx: 530 # there is a variable in the greybox block that 531 # is not in the NL - append this to the end 532 self._vardata_to_idx[var] = n_primals 533 n_primals += 1 534 greybox_primals.append(var) 535 self._greybox_primals_names.append(var.getname(fully_qualified=True)) 536 self._n_greybox_primals = len(greybox_primals) 537 self._greybox_primal_variables = greybox_primals 538 539 # Configure the primal and dual data caches 540 self._greybox_primals_lb = np.zeros(self._n_greybox_primals) 541 self._greybox_primals_ub = np.zeros(self._n_greybox_primals) 542 self._init_greybox_primals = np.zeros(self._n_greybox_primals) 543 for i, var in enumerate(greybox_primals): 544 if var.value is not None: 545 self._init_greybox_primals[i] = var.value 546 self._greybox_primals_lb[i] = -np.inf if var.lb is None else var.lb 547 self._greybox_primals_ub[i] = np.inf if var.ub is None else var.ub 548 self._greybox_primals_lb.flags.writeable = False 549 self._greybox_primals_ub.flags.writeable = False 550 551 self._greybox_primals = self._init_greybox_primals.copy() 552 553 # data member to store the cached greybox constraints and jacobian 554 self._cached_greybox_constraints = None 555 self._cached_greybox_jac = None 556 557 # create the helper objects 558 con_offset = self._pyomo_nlp.n_constraints() 559 self._external_greybox_helpers = [] 560 for data in greybox_data: 561 h = _ExternalGreyBoxModelHelper(data, self._vardata_to_idx, con_offset) 562 self._external_greybox_helpers.append(h) 563 con_offset += h.n_residuals() 564 565 self._n_greybox_constraints = con_offset - self._pyomo_nlp.n_constraints() 566 assert len(self._greybox_constraints_names) == self._n_greybox_constraints 567 568 # make sure the primal values get to the greybox models 569 self.set_primals(self.get_primals()) 570 571 # If any part of the problem is scaled (i.e., obj, primals, 572 # or any of the constraints for any of the external grey boxes), 573 # then we want scaling factors for everything (defaulting to 574 # ones for any of the missing factors). 575 # This code builds all the scaling factors, with the defaults, 576 # but then sets them all to None if *no* scaling factors are provided 577 # for any of the pieces. (inefficient, but only done once) 578 need_scaling = False 579 self._obj_scaling = self._pyomo_nlp.get_obj_scaling() 580 if self._obj_scaling is None: 581 self._obj_scaling = 1.0 582 else: 583 need_scaling = True 584 585 self._primals_scaling = np.ones(self.n_primals()) 586 scaling_suffix = self._pyomo_nlp._pyomo_model.component('scaling_factor') 587 if scaling_suffix and scaling_suffix.ctype is pyo.Suffix: 588 need_scaling = True 589 for i,v in enumerate(self.get_pyomo_variables()): 590 if v in scaling_suffix: 591 self._primals_scaling[i] = scaling_suffix[v] 592 593 self._constraints_scaling = [] 594 pyomo_nlp_scaling = self._pyomo_nlp.get_constraints_scaling() 595 if pyomo_nlp_scaling is None: 596 pyomo_nlp_scaling = np.ones(self._pyomo_nlp.n_constraints()) 597 else: 598 need_scaling = True 599 self._constraints_scaling.append(pyomo_nlp_scaling) 600 601 for h in self._external_greybox_helpers: 602 tmp_scaling = h.get_residual_scaling() 603 if tmp_scaling is None: 604 tmp_scaling = np.ones(h.n_residuals()) 605 else: 606 need_scaling = True 607 self._constraints_scaling.append(tmp_scaling) 608 609 if need_scaling: 610 self._constraints_scaling = np.concatenate(self._constraints_scaling) 611 else: 612 self._obj_scaling = None 613 self._primals_scaling = None 614 self._constraints_scaling = None 615 616 # might want the user to be able to specify these at some point 617 self._init_greybox_duals = np.zeros(self._n_greybox_constraints) 618 self._init_greybox_primals.flags.writeable = False 619 self._init_greybox_duals.flags.writeable = False 620 self._greybox_duals = self._init_greybox_duals.copy() 621 622 # compute the jacobian for the external greybox models 623 # to get some of the statistics 624 self._evaluate_greybox_jacobians_and_cache_if_necessary() 625 self._nnz_greybox_jac = len(self._cached_greybox_jac.data) 626 627 # make sure the duals get to the external greybox models 628 self.set_duals(self.get_duals()) 629 630 # compute the hessian for the external greybox models 631 # to get some of the statistics 632 try: 633 self._evaluate_greybox_hessians_and_cache_if_necessary() 634 self._nnz_greybox_hess = len(self._cached_greybox_hess.data) 635 except (AttributeError, NotImplementedError): 636 self._nnz_greybox_hess = None 637 638 def _invalidate_greybox_primals_cache(self): 639 self._greybox_constraints_cached = False 640 self._greybox_jac_cached = False 641 self._greybox_hess_cached = False 642 643 def _invalidate_greybox_duals_cache(self): 644 self._greybox_hess_cached = False 645 646 # overloaded from NLP 647 def n_primals(self): 648 return self._pyomo_nlp.n_primals() + self._n_greybox_primals 649 650 # overloaded from NLP 651 def n_constraints(self): 652 return self._pyomo_nlp.n_constraints() + self._n_greybox_constraints 653 654 # overloaded from ExtendedNLP 655 def n_eq_constraints(self): 656 return self._pyomo_nlp.n_eq_constraints() + self._n_greybox_constraints 657 658 # overloaded from ExtendedNLP 659 def n_ineq_constraints(self): 660 return self._pyomo_nlp.n_ineq_constraints() 661 662 # overloaded from NLP 663 def nnz_jacobian(self): 664 return self._pyomo_nlp.nnz_jacobian() + self._nnz_greybox_jac 665 666 # overloaded from AslNLP 667 def nnz_jacobian_eq(self): 668 return self._pyomo_nlp.nnz_jacobian_eq() + self._nnz_greybox_jac 669 670 # overloaded from NLP 671 def nnz_hessian_lag(self): 672 return self._pyomo_nlp.nnz_hessian_lag() + self._nnz_greybox_hess 673 674 # overloaded from NLP 675 def primals_lb(self): 676 return np.concatenate((self._pyomo_nlp.primals_lb(), 677 self._greybox_primals_lb, 678 )) 679 680 # overloaded from NLP 681 def primals_ub(self): 682 return np.concatenate(( 683 self._pyomo_nlp.primals_ub(), 684 self._greybox_primals_ub, 685 )) 686 687 # overloaded from NLP 688 def constraints_lb(self): 689 return np.concatenate(( 690 self._pyomo_nlp.constraints_lb(), 691 np.zeros(self._n_greybox_constraints, dtype=np.float64), 692 )) 693 694 # overloaded from NLP 695 def constraints_ub(self): 696 return np.concatenate(( 697 self._pyomo_nlp.constraints_ub(), 698 np.zeros(self._n_greybox_constraints, dtype=np.float64), 699 )) 700 701 # overloaded from NLP 702 def init_primals(self): 703 return np.concatenate(( 704 self._pyomo_nlp.init_primals(), 705 self._init_greybox_primals, 706 )) 707 708 # overloaded from NLP 709 def init_duals(self): 710 return np.concatenate(( 711 self._pyomo_nlp.init_duals(), 712 self._init_greybox_duals, 713 )) 714 715 # overloaded from ExtendedNLP 716 def init_duals_eq(self): 717 return np.concatenate(( 718 self._pyomo_nlp.init_duals_eq(), 719 self._init_greybox_duals, 720 )) 721 722 # overloaded from NLP / Extended NLP 723 def create_new_vector(self, vector_type): 724 if vector_type == 'primals': 725 return np.zeros(self.n_primals(), dtype=np.float64) 726 elif vector_type == 'constraints' or vector_type == 'duals': 727 return np.zeros(self.n_constraints(), dtype=np.float64) 728 elif vector_type == 'eq_constraints' or vector_type == 'duals_eq': 729 return np.zeros(self.n_eq_constraints(), dtype=np.float64) 730 elif vector_type == 'ineq_constraints' or vector_type == 'duals_ineq': 731 return np.zeros(self.n_ineq_constraints(), dtype=np.float64) 732 else: 733 raise RuntimeError('Called create_new_vector with an unknown vector_type') 734 735 # overloaded from NLP 736 def set_primals(self, primals): 737 self._invalidate_greybox_primals_cache() 738 739 # set the primals on the "pyomo" part of the nlp 740 self._pyomo_nlp.set_primals( 741 primals[:self._pyomo_nlp.n_primals()]) 742 743 # copy the values for the greybox primals 744 np.copyto(self._greybox_primals, primals[self._pyomo_nlp.n_primals():]) 745 746 for external in self._external_greybox_helpers: 747 external.set_primals(primals) 748 749 # overloaded from AslNLP 750 def get_primals(self): 751 # return the value of the primals that the pyomo 752 # part knows about as well as any extra values that 753 # are only in the greybox part 754 return np.concatenate(( 755 self._pyomo_nlp.get_primals(), 756 self._greybox_primals, 757 )) 758 759 # overloaded from NLP 760 def set_duals(self, duals): 761 self._invalidate_greybox_duals_cache() 762 763 # set the duals for the pyomo part of the nlp 764 self._pyomo_nlp.set_duals( 765 duals[:self._pyomo_nlp.n_constraints()]) 766 767 # set the duals for the greybox part of the nlp 768 np.copyto(self._greybox_duals, duals[self._pyomo_nlp.n_constraints():]) 769 770 # set the duals in the helpers for the hessian computation 771 for h in self._external_greybox_helpers: 772 h.set_duals(duals) 773 774 # overloaded from NLP 775 def get_duals(self): 776 # return the duals for the pyomo part of the nlp 777 # concatenated with the greybox part 778 return np.concatenate(( 779 self._pyomo_nlp.get_duals(), 780 self._greybox_duals, 781 )) 782 783 # overloaded from ExtendedNLP 784 def set_duals_eq(self, duals): 785 raise NotImplementedError('set_duals_eq not implemented for PyomoGreyBoxNLP') 786 # we think the code below is correct, but it has not yet been tested 787 """ 788 #self._invalidate_greybox_duals_cache() 789 790 # set the duals for the pyomo part of the nlp 791 self._pyomo_nlp.set_duals_eq( 792 duals[:self._pyomo_nlp.n_equality_constraints()]) 793 794 # set the duals for the greybox part of the nlp 795 np.copyto(self._greybox_duals, duals[self._pyomo_nlp.n_equality_constraints():]) 796 # set the duals in the helpers for the hessian computation 797 for h in self._external_greybox_helpers: 798 h.set_duals_eq(duals) 799 """ 800 801 # TODO: Implement set_duals_ineq 802 803 # overloaded from NLP 804 def get_duals_eq(self): 805 raise NotImplementedError('get_duals_eq not implemented for PyomoGreyBoxNLP') 806 """ 807 # return the duals for the pyomo part of the nlp 808 # concatenated with the greybox part 809 return np.concatenate(( 810 self._pyomo_nlp.get_duals_eq(), 811 self._greybox_duals, 812 )) 813 """ 814 815 # overloaded from NLP 816 def set_obj_factor(self, obj_factor): 817 # objective is owned by the pyomo model 818 self._pyomo_nlp.set_obj_factor(obj_factor) 819 820 # overloaded from NLP 821 def get_obj_factor(self): 822 # objective is owned by the pyomo model 823 return self._pyomo_nlp.get_obj_factor() 824 825 # overloaded from NLP 826 def get_obj_scaling(self): 827 return self._obj_scaling 828 829 # overloaded from NLP 830 def get_primals_scaling(self): 831 return self._primals_scaling 832 833 # overloaded from NLP 834 def get_constraints_scaling(self): 835 return self._constraints_scaling 836 837 # overloaded from NLP 838 def evaluate_objective(self): 839 # objective is owned by the pyomo model 840 return self._pyomo_nlp.evaluate_objective() 841 842 # overloaded from NLP 843 def evaluate_grad_objective(self, out=None): 844 # objective is owned by the pyomo model 845 return np.concatenate(( 846 self._pyomo_nlp.evaluate_grad_objective(out), 847 np.zeros(self._n_greybox_primals))) 848 849 def _evaluate_greybox_constraints_and_cache_if_necessary(self): 850 if self._greybox_constraints_cached: 851 return 852 853 self._cached_greybox_constraints = np.concatenate(tuple( 854 external.evaluate_residuals() 855 for external in self._external_greybox_helpers)) 856 self._greybox_constraints_cached = True 857 858 # overloaded from NLP 859 def evaluate_constraints(self, out=None): 860 self._evaluate_greybox_constraints_and_cache_if_necessary() 861 862 if out is not None: 863 if not isinstance(out, np.ndarray) \ 864 or out.size != self.n_constraints(): 865 raise RuntimeError( 866 'Called evaluate_constraints with an invalid' 867 ' "out" argument - should take an ndarray of ' 868 'size {}'.format(self.n_constraints())) 869 870 # call on the pyomo part of the nlp 871 self._pyomo_nlp.evaluate_constraints( 872 out[:-self._n_greybox_constraints]) 873 874 # call on the greybox part of the nlp 875 np.copyto(out[-self._n_greybox_constraints:], 876 self._cached_greybox_constraints) 877 return out 878 879 else: 880 # concatenate the pyomo and external constraint residuals 881 return np.concatenate(( 882 self._pyomo_nlp.evaluate_constraints(), 883 self._cached_greybox_constraints, 884 )) 885 886 # overloaded from ExtendedNLP 887 def evaluate_eq_constraints(self, out=None): 888 raise NotImplementedError('Not yet implemented for PyomoGreyBoxNLP') 889 """ 890 self._evaluate_greybox_constraints_and_cache_if_necessary() 891 892 if out is not None: 893 if not isinstance(out, np.ndarray) \ 894 or out.size != self.n_eq_constraints(): 895 raise RuntimeError( 896 'Called evaluate_eq_constraints with an invalid' 897 ' "out" argument - should take an ndarray of ' 898 'size {}'.format(self.n_eq_constraints())) 899 self._pyomo_nlp.evaluate_eq_constraints( 900 out[:-self._n_greybox_constraints]) 901 np.copyto(out[-self._n_greybox_constraints:], self._cached_greybox_constraints) 902 return out 903 else: 904 return np.concatenate(( 905 self._pyomo_nlp.evaluate_eq_constraints(), 906 self._cached_greybox_constraints, 907 )) 908 """ 909 910 def _evaluate_greybox_jacobians_and_cache_if_necessary(self): 911 if self._greybox_jac_cached: 912 return 913 914 jac = BlockMatrix(len(self._external_greybox_helpers), 1) 915 for i, external in enumerate(self._external_greybox_helpers): 916 jac.set_block(i, 0, external.evaluate_jacobian()) 917 self._cached_greybox_jac = jac.tocoo() 918 self._greybox_jac_cached = True 919 920 # overloaded from NLP 921 def evaluate_jacobian(self, out=None): 922 self._evaluate_greybox_jacobians_and_cache_if_necessary() 923 924 if out is not None: 925 if ( not isinstance(out, coo_matrix) 926 or out.shape[0] != self.n_constraints() 927 or out.shape[1] != self.n_primals() 928 or out.nnz != self.nnz_jacobian() ): 929 raise RuntimeError( 930 'evaluate_jacobian called with an "out" argument' 931 ' that is invalid. This should be a coo_matrix with' 932 ' shape=({},{}) and nnz={}' 933 .format(self.n_constraints(), self.n_primals(), 934 self.nnz_jacobian())) 935 n_pyomo_constraints = self.n_constraints() - self._n_greybox_constraints 936 937 # to avoid an additional copy, we pass in a slice (numpy view) of the underlying 938 # data, row, and col that we were handed to be populated in evaluate_jacobian 939 self._pyomo_nlp.evaluate_jacobian( 940 out=coo_matrix((out.data[:-self._nnz_greybox_jac], 941 (out.row[:-self._nnz_greybox_jac], 942 out.col[:-self._nnz_greybox_jac])), 943 shape=(n_pyomo_constraints, self._pyomo_nlp.n_primals()))) 944 np.copyto(out.data[-self._nnz_greybox_jac:], 945 self._cached_greybox_jac.data) 946 return out 947 else: 948 base = self._pyomo_nlp.evaluate_jacobian() 949 base = coo_matrix((base.data, (base.row, base.col)), 950 shape=(base.shape[0], self.n_primals())) 951 952 jac = BlockMatrix(2,1) 953 jac.set_block(0, 0, base) 954 jac.set_block(1, 0, self._cached_greybox_jac) 955 return jac.tocoo() 956 957 # TODO: Doesn't this need a "shape" specification? 958 #return coo_matrix(( 959 # np.concatenate((base.data, self._cached_greybox_jac.data)), 960 # ( np.concatenate((base.row, self._cached_greybox_jac.row)), 961 # np.concatenate((base.col, self._cached_greybox_jac.col)) ) 962 #)) 963 964 # overloaded from ExtendedNLP 965 """ 966 def evaluate_jacobian_eq(self, out=None): 967 raise NotImplementedError() 968 self._evaluate_greybox_jacobians_and_cache_if_necessary() 969 970 if out is not None: 971 if ( not isinstance(out, coo_matrix) 972 or out.shape[0] != self.n_eq_constraints() 973 or out.shape[1] != self.n_primals() 974 or out.nnz != self.nnz_jacobian_eq() ): 975 raise RuntimeError( 976 'evaluate_jacobian called with an "out" argument' 977 ' that is invalid. This should be a coo_matrix with' 978 ' shape=({},{}) and nnz={}' 979 .format(self.n_eq_constraints(), self.n_primals(), 980 self.nnz_jacobian_eq())) 981 self._pyomo_nlp.evaluate_jacobian_eq( 982 coo_matrix((out.data[:-self._nnz_greybox_jac], 983 (out.row[:-self._nnz_greybox_jac], 984 out.col[:-self._nnz_greybox_jac]))) 985 ) 986 np.copyto(out.data[-self._nnz_greybox_jac], 987 self._cached_greybox_jac.data) 988 return out 989 else: 990 base = self._pyomo_nlp.evaluate_jacobian_eq() 991 # TODO: Doesn't this need a "shape" specification? 992 return coo_matrix(( 993 np.concatenate((base.data, self._cached_greybox_jac.data)), 994 ( np.concatenate((base.row, self._cached_greybox_jac.row)), 995 np.concatenate((base.col, self._cached_greybox_jac.col)) ) 996 )) 997 """ 998 999 # overloaded from NLP 1000 def _evaluate_greybox_hessians_and_cache_if_necessary(self): 1001 if self._greybox_hess_cached: 1002 return 1003 1004 data = list() 1005 irow = list() 1006 jcol = list() 1007 for external in self._external_greybox_helpers: 1008 hess = external.evaluate_hessian() 1009 data.append(hess.data) 1010 irow.append(hess.row) 1011 jcol.append(hess.col) 1012 1013 data = np.concatenate(data) 1014 irow = np.concatenate(irow) 1015 jcol = np.concatenate(jcol) 1016 1017 self._cached_greybox_hess = coo_matrix( (data, (irow,jcol)), shape=(self.n_primals(), self.n_primals())) 1018 self._greybox_hess_cached = True 1019 1020 def evaluate_hessian_lag(self, out=None): 1021 self._evaluate_greybox_hessians_and_cache_if_necessary() 1022 if out is not None: 1023 if ( not isinstance(out, coo_matrix) 1024 or out.shape[0] != self.n_primals() 1025 or out.shape[1] != self.n_primals() 1026 or out.nnz != self.nnz_hessian_lag() ): 1027 raise RuntimeError( 1028 'evaluate_hessian_lag called with an "out" argument' 1029 ' that is invalid. This should be a coo_matrix with' 1030 ' shape=({},{}) and nnz={}' 1031 .format(self.n_primals(), self.n_primals(), 1032 self.nnz_hessian())) 1033 # to avoid an additional copy, we pass in a slice (numpy view) of the underlying 1034 # data, row, and col that we were handed to be populated in evaluate_hessian_lag 1035 # the coo_matrix is simply a holder of the data, row, and col structures 1036 self._pyomo_nlp.evaluate_hessian_lag( 1037 out=coo_matrix((out.data[:-self._nnz_greybox_hess], 1038 (out.row[:-self._nnz_greybox_hess], 1039 out.col[:-self._nnz_greybox_hess])), 1040 shape=(self._pyomo_nlp.n_primals(), self._pyomo_nlp.n_primals()))) 1041 np.copyto(out.data[-self._nnz_greybox_hess:], 1042 self._cached_greybox_hess.data) 1043 return out 1044 else: 1045 hess = self._pyomo_nlp.evaluate_hessian_lag() 1046 data = np.concatenate((hess.data, self._cached_greybox_hess.data)) 1047 row = np.concatenate((hess.row, self._cached_greybox_hess.row)) 1048 col = np.concatenate((hess.col, self._cached_greybox_hess.col)) 1049 hess = coo_matrix((data, (row, col)), shape=(self.n_primals(), self.n_primals())) 1050 return hess 1051 1052 # overloaded from NLP 1053 def report_solver_status(self, status_code, status_message): 1054 raise NotImplementedError('Todo: implement this') 1055 1056 @deprecated(msg='This method has been replaced with primals_names', version='6.0.0.dev0', remove_in='6.0') 1057 def variable_names(self): 1058 return self.primals_names() 1059 1060 def primals_names(self): 1061 names = list(self._pyomo_nlp.variable_names()) 1062 names.extend(self._greybox_primals_names) 1063 return names 1064 1065 def constraint_names(self): 1066 names = list(self._pyomo_nlp.constraint_names()) 1067 names.extend(self._greybox_constraints_names) 1068 return names 1069 1070 def pyomo_model(self): 1071 """ 1072 Return optimization model 1073 """ 1074 return self._pyomo_model 1075 1076 def get_pyomo_objective(self): 1077 """ 1078 Return an instance of the active objective function on the Pyomo model. 1079 (there can be only one) 1080 """ 1081 return self._pyomo_nlp.get_pyomo_objective() 1082 1083 def get_pyomo_variables(self): 1084 """ 1085 Return an ordered list of the Pyomo VarData objects in 1086 the order corresponding to the primals 1087 """ 1088 return self._pyomo_nlp.get_pyomo_variables() + \ 1089 self._greybox_primal_variables 1090 1091 def get_pyomo_constraints(self): 1092 """ 1093 Return an ordered list of the Pyomo ConData objects in 1094 the order corresponding to the primals 1095 """ 1096 # FIXME: what do we return for the external block constraints? 1097 # return self._pyomo_nlp.get_pyomo_constraints() 1098 raise NotImplementedError( 1099 "returning list of all constraints when using an external " 1100 "model is TBD") 1101 1102 def load_state_into_pyomo(self, bound_multipliers=None): 1103 primals = self.get_primals() 1104 variables = self.get_pyomo_variables() 1105 for var, val in zip(variables, primals): 1106 var.set_value(val) 1107 m = self.pyomo_model() 1108 model_suffixes = dict( 1109 pyo.suffix.active_import_suffix_generator(m)) 1110 if 'dual' in model_suffixes: 1111 model_suffixes['dual'].clear() 1112 # Until we sort out how to return the duals for the external 1113 # block (implied) constraints, I am disabling *all* duals 1114 # 1115 # duals = self.get_duals() 1116 # constraints = self.get_pyomo_constraints() 1117 # model_suffixes['dual'].update( 1118 # zip(constraints, duals)) 1119 if 'ipopt_zL_out' in model_suffixes: 1120 model_suffixes['ipopt_zL_out'].clear() 1121 if bound_multipliers is not None: 1122 model_suffixes['ipopt_zL_out'].update( 1123 zip(variables, bound_multipliers[0])) 1124 if 'ipopt_zU_out' in model_suffixes: 1125 model_suffixes['ipopt_zU_out'].clear() 1126 if bound_multipliers is not None: 1127 model_suffixes['ipopt_zU_out'].update( 1128 zip(variables, bound_multipliers[1])) 1129 1130 1131class _ExternalGreyBoxModelHelper(object): 1132 def __init__(self, ex_grey_box_block, vardata_to_idx, con_offset): 1133 """This helper takes an ExternalGreyBoxModel and provides the residual, 1134 Jacobian, and Hessian computation mapped to the correct variable space. 1135 1136 The ExternalGreyBoxModel provides an interface that supports 1137 equality constraints (pure residuals) and output equations. Let 1138 u be the inputs, o be the outputs, and x be the full set of 1139 primal variables from the entire pyomo_nlp. 1140 1141 With this, the ExternalGreyBoxModel provides the residual 1142 computations w_eq(u), and w_o(u), as well as the Jacobians, 1143 Jw_eq(u), and Jw_o(u). This helper provides h(x)=0, where h(x) = 1144 [h_eq(x); h_o(x)-o] and h_eq(x)=w_eq(Pu*x), and 1145 h_o(x)=w_o(Pu*x), and Pu is a mapping from the full primal 1146 variables "x" to the inputs "u". 1147 1148 It also provides the Jacobian of h w.r.t. x. 1149 J_h(x) = [Jw_eq(Pu*x); Jw_o(Pu*x)-Po*x] 1150 where Po is a mapping from the full primal variables "x" to the 1151 outputs "o". 1152 1153 """ 1154 self._block = ex_grey_box_block 1155 self._ex_model = ex_grey_box_block.get_external_model() 1156 self._n_primals = len(vardata_to_idx) 1157 n_inputs = len(self._block.inputs) 1158 assert n_inputs == self._ex_model.n_inputs() 1159 n_eq_constraints = self._ex_model.n_equality_constraints() 1160 n_outputs = len(self._block.outputs) 1161 assert n_outputs == self._ex_model.n_outputs() 1162 1163 # store the map of input indices (0 .. n_inputs) to 1164 # the indices in the full primals vector 1165 self._inputs_to_primals_map = np.fromiter( 1166 (vardata_to_idx[v] for v in self._block.inputs.values()), 1167 dtype=np.int64, count=n_inputs) 1168 1169 # store the map of output indices (0 .. n_outputs) to 1170 # the indices in the full primals vector 1171 self._outputs_to_primals_map = np.fromiter( 1172 (vardata_to_idx[v] for v in self._block.outputs.values()), 1173 dtype=np.int64, count=n_outputs) 1174 1175 if self._ex_model.n_outputs() == 0 and \ 1176 self._ex_model.n_equality_constraints() == 0: 1177 raise ValueError( 1178 'ExternalGreyBoxModel has no equality constraints ' 1179 'or outputs. It must have at least one or both.') 1180 1181 self._ex_eq_duals_to_full_map = None 1182 if n_eq_constraints > 0: 1183 self._ex_eq_duals_to_full_map = \ 1184 list(range(con_offset, con_offset + n_eq_constraints)) 1185 1186 self._ex_output_duals_to_full_map = None 1187 if n_outputs > 0: 1188 self._ex_output_duals_to_full_map = \ 1189 list(range(con_offset + n_eq_constraints, con_offset + n_eq_constraints + n_outputs)) 1190 1191 # we need to change the column indices in the jacobian 1192 # from the 0..n_inputs provided by the external model 1193 # to the indices corresponding to the full Pyomo model 1194 # so we create that here 1195 self._eq_jac_primal_jcol = None 1196 self._outputs_jac_primal_jcol = None 1197 self._additional_output_entries_irow = None 1198 self._additional_output_entries_jcol = None 1199 self._additional_output_entries_data = None 1200 self._con_offset = con_offset 1201 self._eq_hess_jcol = None 1202 self._eq_hess_irow = None 1203 self._output_hess_jcol = None 1204 self._output_hess_irow = None 1205 1206 def set_primals(self, primals): 1207 # map the full primals "x" to the inputs "u" and set 1208 # the values on the external model 1209 input_values = primals[self._inputs_to_primals_map] 1210 self._ex_model.set_input_values(input_values) 1211 1212 # map the full primals "x" to the outputs "o" and 1213 # store a vector of the current output values to 1214 # use when evaluating residuals 1215 self._output_values = primals[self._outputs_to_primals_map] 1216 1217 def set_duals(self, duals): 1218 # map the full duals to the duals for the equality 1219 # and the output constraints 1220 if self._ex_eq_duals_to_full_map is not None: 1221 duals_eq = duals[self._ex_eq_duals_to_full_map] 1222 self._ex_model.set_equality_constraint_multipliers(duals_eq) 1223 1224 if self._ex_output_duals_to_full_map is not None: 1225 duals_outputs = duals[self._ex_output_duals_to_full_map] 1226 self._ex_model.set_output_constraint_multipliers(duals_outputs) 1227 1228 def n_residuals(self): 1229 return self._ex_model.n_equality_constraints() \ 1230 + self._ex_model.n_outputs() 1231 1232 def get_residual_scaling(self): 1233 eq_scaling = self._ex_model.get_equality_constraint_scaling_factors() 1234 output_con_scaling = self._ex_model.get_output_constraint_scaling_factors() 1235 if eq_scaling is None and output_con_scaling is None: 1236 return None 1237 if eq_scaling is None: 1238 eq_scaling = np.ones(self._ex_model.n_equality_constraints()) 1239 if output_con_scaling is None: 1240 output_con_scaling = np.ones(self._ex_model.n_outputs()) 1241 1242 return np.concatenate(( 1243 eq_scaling, 1244 output_con_scaling)) 1245 1246 def evaluate_residuals(self): 1247 # evalute the equality constraints and the output equations 1248 # and return a single vector of residuals 1249 # returns residual for h(x)=0, where h(x) = [h_eq(x); h_o(x)-o] 1250 resid_list = [] 1251 if self._ex_model.n_equality_constraints() > 0: 1252 resid_list.append(self._ex_model.evaluate_equality_constraints()) 1253 1254 if self._ex_model.n_outputs() > 0: 1255 computed_output_values = self._ex_model.evaluate_outputs() 1256 output_resid = computed_output_values - self._output_values 1257 resid_list.append(output_resid) 1258 1259 return np.concatenate(resid_list) 1260 1261 def evaluate_jacobian(self): 1262 # compute the jacobian of h(x) w.r.t. x 1263 # J_h(x) = [Jw_eq(Pu*x); Jw_o(Pu*x)-Po*x] 1264 1265 # Jw_eq(x) 1266 eq_jac = None 1267 if self._ex_model.n_equality_constraints() > 0: 1268 eq_jac = self._ex_model.evaluate_jacobian_equality_constraints() 1269 if self._eq_jac_primal_jcol is None: 1270 # The first time through, we won't have created the 1271 # mapping of external primals ('u') to the full space 1272 # primals ('x') 1273 self._eq_jac_primal_jcol = self._inputs_to_primals_map[ 1274 eq_jac.col] 1275 # map the columns from the inputs "u" back to the full primals "x" 1276 eq_jac = coo_matrix( 1277 (eq_jac.data, (eq_jac.row, self._eq_jac_primal_jcol)), 1278 (eq_jac.shape[0], self._n_primals)) 1279 1280 outputs_jac = None 1281 if self._ex_model.n_outputs() > 0: 1282 outputs_jac = self._ex_model.evaluate_jacobian_outputs() 1283 1284 row = outputs_jac.row 1285 # map the columns from the inputs "u" back to the full primals "x" 1286 if self._outputs_jac_primal_jcol is None: 1287 # The first time through, we won't have created the 1288 # mapping of external outputs ('o') to the full space 1289 # primals ('x') 1290 self._outputs_jac_primal_jcol = self._inputs_to_primals_map[ 1291 outputs_jac.col] 1292 1293 # We also need tocreate the irow, jcol, nnz structure for the 1294 # output variable portion of h(u)-o=0 1295 self._additional_output_entries_irow = np.asarray(range(self._ex_model.n_outputs())) 1296 self._additional_output_entries_jcol = self._outputs_to_primals_map 1297 self._additional_output_entries_data = -1.0*np.ones(self._ex_model.n_outputs()) 1298 1299 col = self._outputs_jac_primal_jcol 1300 data = outputs_jac.data 1301 1302 # add the additional entries for the -Po*x portion of the jacobian 1303 row = np.concatenate((row, self._additional_output_entries_irow)) 1304 col = np.concatenate((col, self._additional_output_entries_jcol)) 1305 data = np.concatenate((data, self._additional_output_entries_data)) 1306 outputs_jac = coo_matrix( 1307 (data, (row, col)), 1308 shape=(outputs_jac.shape[0], self._n_primals)) 1309 1310 jac = None 1311 if eq_jac is not None: 1312 if outputs_jac is not None: 1313 # create a jacobian with both Jw_eq and Jw_o 1314 jac = BlockMatrix(2,1) 1315 jac.name = 'external model jacobian' 1316 jac.set_block(0,0,eq_jac) 1317 jac.set_block(1,0,outputs_jac) 1318 else: 1319 assert self._ex_model.n_outputs() == 0 1320 assert self._ex_model.n_equality_constraints() > 0 1321 # only need the Jacobian with Jw_eq (there are not 1322 # output equations) 1323 jac = eq_jac 1324 else: 1325 assert outputs_jac is not None 1326 assert self._ex_model.n_outputs() > 0 1327 assert self._ex_model.n_equality_constraints() == 0 1328 # only need the Jacobian with Jw_o (there are no equalities) 1329 jac = outputs_jac 1330 1331 return jac 1332 1333 def evaluate_hessian(self): 1334 # compute the portion of the Hessian of the Lagrangian 1335 # of h(x) w.r.t. x 1336 # H_h(x) = sum_i y_eq_i * Hw_eq(Pu*x) +_ sum_k y_o_j * Hw_o(Pu*x)] 1337 1338 data_list = [] 1339 irow_list = [] 1340 jcol_list = [] 1341 1342 # Hw_eq(x) 1343 eq_hess = None 1344 if self._ex_model.n_equality_constraints() > 0: 1345 eq_hess = self._ex_model.evaluate_hessian_equality_constraints() 1346 if self._eq_hess_jcol is None: 1347 # first time through, let's also check that it is lower triangular 1348 if np.any(eq_hess.row < eq_hess.col): 1349 raise ValueError('ExternalGreyBoxModel must return lower ' 1350 'triangular portion of the Hessian only') 1351 1352 # The first time through, we won't have created the 1353 # mapping of external primals ('u') to the full space 1354 # primals ('x') 1355 self._eq_hess_irow = row = self._inputs_to_primals_map[eq_hess.row] 1356 self._eq_hess_jcol = col = self._inputs_to_primals_map[eq_hess.col] 1357 1358 # mapping may have made this not lower triangular 1359 mask = col > row 1360 row[mask], col[mask] = col[mask], row[mask] 1361 1362 data_list.append(eq_hess.data) 1363 irow_list.append(self._eq_hess_irow) 1364 jcol_list.append(self._eq_hess_jcol) 1365 1366 if self._ex_model.n_outputs() > 0: 1367 output_hess = self._ex_model.evaluate_hessian_outputs() 1368 if self._output_hess_irow is None: 1369 # first time through, let's also check that it is lower triangular 1370 if np.any(output_hess.row < output_hess.col): 1371 raise ValueError('ExternalGreyBoxModel must return lower ' 1372 'triangular portion of the Hessian only') 1373 1374 # The first time through, we won't have created the 1375 # mapping of external outputs ('o') to the full space 1376 # primals ('x') 1377 self._output_hess_irow = row = self._inputs_to_primals_map[output_hess.row] 1378 self._output_hess_jcol = col = self._inputs_to_primals_map[output_hess.col] 1379 1380 # mapping may have made this not lower triangular 1381 mask = col > row 1382 row[mask], col[mask] = col[mask], row[mask] 1383 1384 data_list.append(output_hess.data) 1385 irow_list.append(self._output_hess_irow) 1386 jcol_list.append(self._output_hess_jcol) 1387 1388 data = np.concatenate(data_list) 1389 irow = np.concatenate(irow_list) 1390 jcol = np.concatenate(jcol_list) 1391 hess = coo_matrix( (data, (irow, jcol)), (self._n_primals, self._n_primals) ) 1392 1393 return hess 1394 1395