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# ______________________________________________________________________________ 10from pyomo.environ import ( 11 Param, Var, Block, ComponentMap, Objective, Constraint, 12 ConstraintList, Suffix, value, ComponentUID, 13) 14 15from pyomo.common.sorting import sorted_robust 16from pyomo.core.expr.current import ExpressionReplacementVisitor 17 18from pyomo.common.modeling import unique_component_name 19from pyomo.common.deprecation import deprecated 20from pyomo.common.tempfiles import TempfileManager 21from pyomo.opt import SolverFactory, SolverStatus 22from pyomo.contrib.sensitivity_toolbox.k_aug import K_augInterface, InTempDir 23import logging 24import os 25import io 26import shutil 27from pyomo.common.dependencies import ( 28 numpy as np, numpy_available 29 ) 30from pyomo.common.dependencies import scipy, scipy_available 31 32logger = logging.getLogger('pyomo.contrib.sensitivity_toolbox') 33 34@deprecated("The sipopt function has been deprecated. Use the sensitivity_calculation() " 35 "function with method='sipopt' to access this functionality.", 36 logger='pyomo.contrib.sensitivity_toolbox', 37 version='6.1') 38def sipopt(instance, paramSubList, perturbList, 39 cloneModel=True, tee=False, keepfiles=False, 40 streamSoln=False): 41 m = sensitivity_calculation('sipopt', instance, paramSubList, perturbList, 42 cloneModel, tee, keepfiles, solver_options=None) 43 44 return m 45 46@deprecated("The kaug function has been deprecated. Use the sensitivity_calculation() " 47 "function with method='k_aug' to access this functionality.", 48 logger='pyomo.contrib.sensitivity_toolbox', 49 version='6.1') 50def kaug(instance, paramSubList, perturbList, 51 cloneModel=True, tee=False, keepfiles=False, solver_options=None, 52 streamSoln=False): 53 m = sensitivity_calculation('k_aug', instance, paramSubList, perturbList, 54 cloneModel, tee, keepfiles, solver_options) 55 56 return m 57 58_SIPOPT_SUFFIXES = { 59 'sens_state_0': Suffix.EXPORT, 60 # ^ Not sure what this suffix does -RBP 61 'sens_state_1': Suffix.EXPORT, 62 'sens_state_value_1': Suffix.EXPORT, 63 'sens_init_constr': Suffix.EXPORT, 64 65 'sens_sol_state_1': Suffix.IMPORT, 66 'sens_sol_state_1_z_L': Suffix.IMPORT, 67 'sens_sol_state_1_z_U': Suffix.IMPORT, 68 } 69 70_K_AUG_SUFFIXES = { 71 'ipopt_zL_out': Suffix.IMPORT, 72 'ipopt_zU_out': Suffix.IMPORT, 73 'ipopt_zL_in': Suffix.EXPORT, 74 'ipopt_zU_in': Suffix.EXPORT, 75 'dual': Suffix.IMPORT_EXPORT, 76 'dcdp': Suffix.EXPORT, 77 'DeltaP': Suffix.EXPORT, 78 } 79 80def _add_sensitivity_suffixes(block): 81 suffix_dict = {} 82 suffix_dict.update(_SIPOPT_SUFFIXES) 83 suffix_dict.update(_K_AUG_SUFFIXES) 84 for name, direction in suffix_dict.items(): 85 if block.component(name) is None: 86 # Only add suffix if it doesn't already exist. 87 # If something of this name does already exist, just 88 # assume it is the proper suffix and move on. 89 block.add_component(name, Suffix(direction=direction)) 90 91class _NotAnIndex(object): 92 pass 93 94def _generate_component_items(components): 95 if type(components) not in {list, tuple}: 96 components = (components,) 97 for comp in components: 98 if comp.is_indexed(): 99 for idx in sorted_robust(comp): 100 yield idx, comp[idx] 101 else: 102 yield _NotAnIndex, comp 103 104def sensitivity_calculation(method, instance, paramList, perturbList, 105 cloneModel=True, tee=False, keepfiles=False, solver_options=None): 106 """This function accepts a Pyomo ConcreteModel, a list of parameters, and 107 their corresponding perturbation list. The model is then augmented with 108 dummy constraints required to call sipopt or k_aug to get an approximation 109 of the perturbed solution. 110 111 Parameters 112 ---------- 113 method: string 114 'sipopt' or 'k_aug' 115 instance: Block 116 pyomo block or model object 117 paramSubList: list 118 list of mutable parameters or fixed variables 119 perturbList: list 120 list of perturbed parameter values 121 cloneModel: bool, optional 122 indicator to clone the model. If set to False, the original 123 model will be altered 124 tee: bool, optional 125 indicator to stream solver log 126 keepfiles: bool, optional 127 preserve solver interface files 128 solver_options: dict, optional 129 Provides options to the solver (also the name of an attribute) 130 131 Returns 132 ------- 133 The model that was manipulated by the sensitivity interface 134 135 """ 136 137 sens = SensitivityInterface(instance, clone_model=cloneModel) 138 sens.setup_sensitivity(paramList) 139 140 m = sens.model_instance 141 142 if method not in {"k_aug", "sipopt"}: 143 raise ValueError("Only methods 'k_aug' and 'sipopt' are supported'") 144 145 if method == 'k_aug': 146 k_aug = SolverFactory('k_aug', solver_io='nl') 147 dot_sens = SolverFactory('dot_sens', solver_io='nl') 148 ipopt = SolverFactory('ipopt', solver_io='nl') 149 150 k_aug_interface = K_augInterface(k_aug=k_aug, dot_sens=dot_sens) 151 152 ipopt.solve(m, tee=tee) 153 m.ipopt_zL_in.update(m.ipopt_zL_out) #: important! 154 m.ipopt_zU_in.update(m.ipopt_zU_out) #: important! 155 156 k_aug.options['dsdp_mode'] = "" #: sensitivity mode! 157 k_aug_interface.k_aug(m, tee=tee) 158 159 sens.perturb_parameters(perturbList) 160 161 if method == 'sipopt': 162 ipopt_sens = SolverFactory('ipopt_sens', solver_io='nl') 163 ipopt_sens.options['run_sens'] = 'yes' 164 165 # Send the model to ipopt_sens and collect the solution 166 results = ipopt_sens.solve(m, keepfiles=keepfiles, tee=tee) 167 168 elif method == 'k_aug': 169 dot_sens.options["dsdp_mode"] = "" 170 k_aug_interface.dot_sens(m, tee=tee) 171 172 return m 173 174def get_dsdp(model, theta_names, theta, tee=False): 175 """This function calculates gradient vector of the variables 176 with respect to the parameters (theta_names). 177 178 e.g) min f: p1*x1+ p2*(x2^2) + p1*p2 179 s.t c1: x1 + x2 = p1 180 c2: x2 + x3 = p2 181 0 <= x1, x2, x3 <= 10 182 p1 = 10 183 p2 = 5 184 the function retuns dx/dp and dp/dp, and column orders. 185 186 The following terms are used to define the output dimensions: 187 Ncon = number of constraints 188 Nvar = number of variables (Nx + Ntheta) 189 Nx = number of decision (primal) variables 190 Ntheta = number of uncertain parameters. 191 192 Parameters 193 ---------- 194 model: Pyomo ConcreteModel 195 model should include an objective function 196 theta_names: list of strings 197 List of Var names 198 theta: dict 199 Estimated parameters e.g) from parmest 200 tee: bool, optional 201 Indicates that ef solver output should be teed 202 203 Returns 204 ------- 205 dsdp: scipy.sparse.csr.csr_matrix 206 Ntheta by Nvar size sparse matrix. A Jacobian matrix of the 207 (decision variables, parameters) with respect to parameters 208 (theta_names). number of rows = len(theta_name), number of 209 columns = len(col) 210 col: list 211 List of variable names 212 """ 213 # Get parameters from names. In SensitivityInterface, we expect 214 # these to be parameters on the original model. 215 param_list = [] 216 for name in theta_names: 217 comp = model.find_component(name) 218 if comp is None: 219 raise RuntimeError("Cannot find component %s on model" % name) 220 if comp.ctype is Var: 221 # If theta_names correspond to Vars in the model, these vars 222 # need to be fixed. 223 comp.fix() 224 param_list.append(comp) 225 226 sens = SensitivityInterface(model, clone_model=True) 227 m = sens.model_instance 228 229 # Setup model and calculate sensitivity matrix with k_aug 230 sens.setup_sensitivity(param_list) 231 k_aug = K_augInterface() 232 k_aug.k_aug(m, tee=tee) 233 234 # Write row and col files in a temp dir, then immediately 235 # read into a Python data structure. 236 nl_data = {} 237 with InTempDir(): 238 base_fname = "col_row" 239 nl_file = ".".join((base_fname, "nl")) 240 row_file = ".".join((base_fname, "row")) 241 col_file = ".".join((base_fname, "col")) 242 m.write(nl_file, io_options={"symbolic_solver_labels": True}) 243 for fname in [nl_file, row_file, col_file]: 244 with open(fname, "r") as fp: 245 nl_data[fname] = fp.read() 246 247 # Create more useful data structures from strings 248 dsdp = np.fromstring(k_aug.data["dsdp_in_.in"], sep="\n\t") 249 col = nl_data[col_file].strip("\n").split("\n") 250 row = nl_data[row_file].strip("\n").split("\n") 251 252 dsdp = dsdp.reshape((len(theta_names), int(len(dsdp)/len(theta_names)))) 253 dsdp = dsdp[:len(theta_names), :len(col)] 254 255 col = [i for i in col if sens.get_default_block_name() not in i] 256 dsdp_out = np.zeros((len(theta_names),len(col))) 257 for i in range(len(theta_names)): 258 for j in range(len(col)): 259 if sens.get_default_block_name() not in col[j]: 260 dsdp_out[i,j] = -dsdp[i, j] # e.g) k_aug dsdp returns -dx1/dx1 = -1.0 261 262 return scipy.sparse.csr_matrix(dsdp_out), col 263 264 265def get_dfds_dcds(model, theta_names, tee=False, solver_options=None): 266 """This function calculates gradient vector of the objective function 267 and constraints with respect to the variables and parameters. 268 269 e.g) min f: p1*x1+ p2*(x2^2) + p1*p2 270 s.t c1: x1 + x2 = p1 271 c2: x2 + x3 = p2 272 0 <= x1, x2, x3 <= 10 273 p1 = 10 274 p2 = 5 275 - Variables = (x1, x2, x3, p1, p2) 276 - Fix p1 and p2 with estimated values 277 278 The following terms are used to define the output dimensions: 279 Ncon = number of constraints 280 Nvar = number of variables (Nx + Ntheta) 281 Nx = number of decision (primal) variables 282 Ntheta = number of uncertain parameters. 283 284 Parameters 285 ---------- 286 model: Pyomo ConcreteModel 287 model should include an objective function 288 theta_names: list of strings 289 List of Var names 290 tee: bool, optional 291 Indicates that ef solver output should be teed 292 solver_options: dict, optional 293 Provides options to the solver (also the name of an attribute) 294 295 Returns 296 ------- 297 gradient_f: numpy.ndarray 298 Length Nvar array. A gradient vector of the objective function 299 with respect to the (decision variables, parameters) at the optimal 300 solution 301 gradient_c: scipy.sparse.csr.csr_matrix 302 Ncon by Nvar size sparse matrix. A Jacobian matrix of the 303 constraints with respect to the (decision variables, parameters) 304 at the optimal solution. Each row contains [column number, 305 row number, and value], column order follows variable order in col 306 and index starts from 1. Note that it follows k_aug. 307 If no constraint exists, return [] 308 col: list 309 Size Nvar list of variable names 310 row: list 311 Size Ncon+1 list of constraints and objective function names. 312 The final element is the objective function name. 313 line_dic: dict 314 column numbers of the theta_names in the model. Index starts from 1 315 316 Raises 317 ------ 318 RuntimeError 319 When ipopt or k_aug or dotsens is not available 320 Exception 321 When ipopt fails 322 """ 323 # Create the solver plugin using the ASL interface 324 ipopt = SolverFactory('ipopt',solver_io='nl') 325 if solver_options is not None: 326 ipopt.options = solver_options 327 k_aug = SolverFactory('k_aug',solver_io='nl') 328 if not ipopt.available(False): 329 raise RuntimeError('ipopt is not available') 330 if not k_aug.available(False): 331 raise RuntimeError('k_aug is not available') 332 333 # Declare Suffixes 334 _add_sensitivity_suffixes(model) 335 336 # K_AUG SUFFIXES 337 model.dof_v = Suffix(direction=Suffix.EXPORT) #: SUFFIX FOR K_AUG 338 model.rh_name = Suffix(direction=Suffix.IMPORT) #: SUFFIX FOR K_AUG AS WELL 339 k_aug.options["print_kkt"] = "" 340 341 results = ipopt.solve(model,tee=tee) 342 343 # Raise exception if ipopt fails 344 if (results.solver.status == SolverStatus.warning): 345 raise Exception(results.solver.Message) 346 347 for o in model.component_objects(Objective, active=True): 348 f_mean = value(o) 349 model.ipopt_zL_in.update(model.ipopt_zL_out) 350 model.ipopt_zU_in.update(model.ipopt_zU_out) 351 352 # run k_aug 353 k_aug_interface = K_augInterface(k_aug=k_aug) 354 k_aug_interface.k_aug(model, tee=tee) #: always call k_aug AFTER ipopt. 355 356 nl_data = {} 357 with InTempDir(): 358 base_fname = "col_row" 359 nl_file = ".".join((base_fname, "nl")) 360 row_file = ".".join((base_fname, "row")) 361 col_file = ".".join((base_fname, "col")) 362 model.write(nl_file, io_options={"symbolic_solver_labels": True}) 363 for fname in [nl_file, row_file, col_file]: 364 with open(fname, "r") as fp: 365 nl_data[fname] = fp.read() 366 367 col = nl_data[col_file].strip("\n").split("\n") 368 row = nl_data[row_file].strip("\n").split("\n") 369 370 # get the column numbers of "parameters" 371 line_dic = {name: col.index(name) for name in theta_names} 372 373 grad_f_file = os.path.join("GJH", "gradient_f_print.txt") 374 grad_f_string = k_aug_interface.data[grad_f_file] 375 gradient_f = np.fromstring(grad_f_string, sep="\n\t") 376 col = [i for i in col if SensitivityInterface.get_default_block_name() not in i] 377 378 grad_c_file = os.path.join("GJH", "A_print.txt") 379 grad_c_string = k_aug_interface.data[grad_c_file] 380 gradient_c = np.fromstring(grad_c_string, sep="\n\t") 381 382 # Jacobian file is in "COO format," i.e. an nnz-by-3 array. 383 # Reshape to a numpy array that matches this format. 384 gradient_c = gradient_c.reshape((-1, 3)) 385 386 num_constraints = len(row)-1 # Objective is included as a row 387 if num_constraints > 0 : 388 row_idx = gradient_c[:,1]-1 389 col_idx = gradient_c[:,0]-1 390 data = gradient_c[:,2] 391 gradient_c = scipy.sparse.csr_matrix((data, (row_idx, col_idx)), 392 shape=(num_constraints, len(col))) 393 else: 394 gradient_c = np.array([]) 395 396 return gradient_f, gradient_c, col,row, line_dic 397 398 399def line_num(file_name, target): 400 """This function returns the line number that contains 'target' in the 401 file_name. This function identifies constraints that have variables 402 in theta_names. 403 404 Parameters 405 ---------- 406 file_name: string 407 file includes the variable order (i.e. col file) 408 target: string 409 variable name to check 410 411 Returns 412 ------- 413 count: int 414 line number of target in the file 415 416 Raises 417 ------ 418 Exception 419 When file does not include target 420 """ 421 with open(file_name) as f: 422 count = int(1) 423 for line in f: 424 if line.strip() == target: 425 return int(count) 426 count += 1 427 raise Exception(file_name + " does not include "+target) 428 429 430class SensitivityInterface(object): 431 432 def __init__(self, instance, clone_model=True): 433 """ Constructor clones model if necessary and attaches 434 to this object. 435 """ 436 self._original_model = instance 437 438 if clone_model: 439 # Note that we are not "cloning" the user's parameters 440 # or perturbations. 441 self.model_instance = instance.clone() 442 else: 443 self.model_instance = instance 444 445 @classmethod 446 def get_default_block_name(self): 447 return '_SENSITIVITY_TOOLBOX_DATA' 448 449 @staticmethod 450 def get_default_var_name(name): 451 #return '_'.join(('sens_var', name)) 452 return name 453 454 @staticmethod 455 def get_default_param_name(name): 456 #return '_'.join(('sens_param', name)) 457 return name 458 459 def _process_param_list(self, paramList): 460 # We need to translate the components in paramList into 461 # components in our possibly cloned model. 462 orig = self._original_model 463 instance = self.model_instance 464 if orig is not instance: 465 paramList = list( 466 ComponentUID(param, context=orig).find_component_on(instance) 467 for param in paramList 468 ) 469 return paramList 470 471 def _add_data_block(self, existing_block=None): 472 # If a sensitivity block already exists, and we have not done 473 # any expression replacement, we delete the old block, re-fix the 474 # sensitivity variables, and start again. 475 # 476 # Don't do this in the constructor as we could want to call 477 # the constructor once, then perform multiple sensitivity 478 # calculations with the same model instance. 479 if existing_block is not None: 480 if (hasattr(existing_block, '_has_replaced_expressions') and 481 not existing_block._has_replaced_expressions): 482 for var, _, _, _ in existing_block._sens_data_list: 483 # Re-fix variables that the previous block was 484 # treating as parameters. 485 var.fix() 486 self.model_instance.del_component(existing_block) 487 else: 488 msg = ("Re-using sensitivity interface is not supported " 489 "when calculating sensitivity for mutable parameters. " 490 "Used fixed vars instead if you want to do this." 491 ) 492 raise RuntimeError(msg) 493 494 # Add a block to keep track of model components necessary for this 495 # sensitivity calculation. 496 block = Block() 497 self.model_instance.add_component(self.get_default_block_name(), block) 498 self.block = block 499 500 # If the user tells us they will perturb a set of parameters, we will 501 # need to replace these parameters in the user's model's constraints. 502 # This affects what we can do with the model later, so we add a flag. 503 block._has_replaced_expressions = False 504 505 # This is the main data structure for keeping track of "sensitivity 506 # vars" and their associated params. It will be a list of tuples of 507 # (vardata, paramdata, list_index, comp_index) 508 # where: 509 # vardata is the "sensitivity variable" data object, 510 # paramdata is the associated parameter, 511 # list_index is its index in the user-provided list, and 512 # comp_index is its index in the component provided by the user. 513 block._sens_data_list = [] 514 515 # This will hold the user-provided list of 516 # variables and/or parameters to perturb 517 block._paramList = None 518 519 # This will hold any constraints where we have replaced 520 # parameters with variables. 521 block.constList = ConstraintList() 522 523 return block 524 525 def _add_sensitivity_data(self, param_list): 526 block = self.block 527 sens_data_list = block._sens_data_list 528 for i, comp in enumerate(param_list): 529 if comp.ctype is Param: 530 parent = comp.parent_component() 531 if not parent.mutable: 532 raise ValueError( 533 "Parameters within paramList must be mutable. " 534 "Got %s, which is not mutable." % comp.name 535 ) 536 # Add a Var: 537 if comp.is_indexed(): 538 d = {k: value(comp[k]) for k in comp.index_set()} 539 var = Var(comp.index_set(), initialize=d) 540 else: 541 d = value(comp) 542 var = Var(initialize=d) 543 name = self.get_default_var_name(parent.local_name) 544 name = unique_component_name(block, name) 545 block.add_component(name, var) 546 547 if comp.is_indexed(): 548 sens_data_list.extend( 549 (var[idx], param, i, idx) 550 for idx, param in _generate_component_items(comp) 551 ) 552 else: 553 sens_data_list.append((var, comp, i, _NotAnIndex)) 554 555 elif comp.ctype is Var: 556 parent = comp.parent_component() 557 for _, data in _generate_component_items(comp): 558 if not data.fixed: 559 raise ValueError( 560 "Specified \"parameter\" variables must be " 561 "fixed. Got %s, which is not fixed." 562 % comp.name 563 ) 564 # Add a Param: 565 if comp.is_indexed(): 566 d = {k: value(comp[k]) for k in comp.index_set()} 567 param = Param(comp.index_set(), mutable=True, initialize=d) 568 else: 569 d = value(comp) 570 param = Param(mutable=True, initialize=d) 571 name = self.get_default_param_name(parent.local_name) 572 name = unique_component_name(block, name) 573 block.add_component(name, param) 574 575 if comp.is_indexed(): 576 sens_data_list.extend( 577 (var, param[idx], i, idx) 578 for idx, var in _generate_component_items(comp) 579 ) 580 else: 581 sens_data_list.append((comp, param, i, _NotAnIndex)) 582 583 def _replace_parameters_in_constraints(self, variableSubMap): 584 instance = self.model_instance 585 block = self.block 586 # Visitor that we will use to replace user-provided parameters 587 # in the objective and the constraints. 588 param_replacer = ExpressionReplacementVisitor( 589 substitute=variableSubMap, 590 remove_named_expressions=True, 591 ) 592 # TODO: Flag to ExpressionReplacementVisitor to only replace 593 # named expressions if a node has been replaced within that 594 # expression. 595 596 new_old_comp_map = ComponentMap() 597 598 # clone Objective, add to Block, and update any Expressions 599 for obj in list(instance.component_data_objects(Objective, 600 active=True, 601 descend_into=True)): 602 tempName = unique_component_name(block, obj.local_name) 603 new_expr = param_replacer.dfs_postorder_stack(obj.expr) 604 block.add_component(tempName, Objective(expr=new_expr)) 605 new_old_comp_map[block.component(tempName)] = obj 606 obj.deactivate() 607 608 # clone Constraints, add to Block, and update any Expressions 609 # 610 # Unfortunate that this deactivates and replaces constraints 611 # even if they don't contain the parameters. 612 # 613 old_con_list = list(instance.component_data_objects(Constraint, 614 active=True, descend_into=True)) 615 last_idx = 0 616 for con in old_con_list: 617 if (con.equality or con.lower is None or con.upper is None): 618 new_expr = param_replacer.dfs_postorder_stack(con.expr) 619 block.constList.add(expr=new_expr) 620 last_idx += 1 621 new_old_comp_map[block.constList[last_idx]] = con 622 else: 623 # Constraint must be a ranged inequality, break into 624 # separate constraints 625 new_body = param_replacer.dfs_postorder_stack(con.body) 626 new_lower = param_replacer.dfs_postorder_stack(con.lower) 627 new_upper = param_replacer.dfs_postorder_stack(con.upper) 628 629 # Add constraint for lower bound 630 block.constList.add(expr=(new_lower <= new_body)) 631 last_idx += 1 632 new_old_comp_map[block.constList[last_idx]] = con 633 634 # Add constraint for upper bound 635 block.constList.add(expr=(new_body <= new_upper)) 636 last_idx += 1 637 new_old_comp_map[block.constList[last_idx]] = con 638 con.deactivate() 639 640 return new_old_comp_map 641 642 def setup_sensitivity(self, paramList): 643 """ 644 """ 645 instance = self.model_instance 646 paramList = self._process_param_list(paramList) 647 648 existing_block = instance.component(self.get_default_block_name()) 649 block = self._add_data_block(existing_block=existing_block) 650 block._sens_data_list = [] 651 block._paramList = paramList 652 653 self._add_sensitivity_data(paramList) 654 sens_data_list = block._sens_data_list 655 656 for var, _, _, _ in sens_data_list: 657 # This unfixes all variables, not just those the user added. 658 var.unfix() 659 660 # Map used to replace user-provided parameters. 661 variableSubMap = dict((id(param), var) 662 for var, param, list_idx, _ in sens_data_list 663 if paramList[list_idx].ctype is Param) 664 665 if variableSubMap: 666 # We now replace the provided parameters in the user's 667 # expressions. Only do this if we have to, i.e. the 668 # user provided some parameters rather than all vars. 669 block._replaced_map = \ 670 self._replace_parameters_in_constraints(variableSubMap) 671 672 # Assume that we just replaced some params 673 block._has_replaced_expressions = True 674 675 block.paramConst = ConstraintList() 676 for var, param, _, _ in sens_data_list: 677 #block.paramConst.add(param - var == 0) 678 block.paramConst.add(var - param == 0) 679 680 # Declare Suffixes 681 _add_sensitivity_suffixes(instance) 682 683 for i, (var, _, _, _) in enumerate(sens_data_list): 684 idx = i + 1 685 con = block.paramConst[idx] 686 687 # sipopt 688 instance.sens_state_0[var] = idx 689 instance.sens_state_1[var] = idx 690 instance.sens_init_constr[con] = idx 691 692 # k_aug 693 instance.dcdp[con] = idx 694 695 696 def perturb_parameters(self, perturbList): 697 """ 698 """ 699 # Note that entries of perturbList need not be components 700 # of the cloned model. All we need are the values. 701 instance = self.model_instance 702 sens_data_list = self.block._sens_data_list 703 paramConst = self.block.paramConst 704 705 if len(self.block._paramList) != len(perturbList): 706 raise ValueError( 707 "Length of paramList argument does not equal " 708 "length of perturbList") 709 710 for i, (var, param, list_idx, comp_idx) in enumerate(sens_data_list): 711 con = paramConst[i+1] 712 if comp_idx is _NotAnIndex: 713 ptb = value(perturbList[list_idx]) 714 else: 715 try: 716 ptb = value(perturbList[list_idx][comp_idx]) 717 except TypeError: 718 # If the user provided a scalar value to perturb 719 # an indexed component. 720 ptb = value(perturbList[list_idx]) 721 722 # sipopt 723 instance.sens_state_value_1[var] = ptb 724 725 # k_aug 726 #instance.DeltaP[con] = value(ptb - var) 727 instance.DeltaP[con] = value(var - ptb) 728 # FIXME: ^ This is incorrect. DeltaP should be (ptb - current). 729 # But at least one test doesn't pass unless I use (current - ptb). 730