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 11from abc import ABCMeta, abstractmethod 12from pyomo.contrib.pynumero.interfaces import pyomo_nlp, ampl_nlp 13from pyomo.contrib.pynumero.sparse import BlockMatrix, BlockVector 14import numpy as np 15import scipy.sparse 16from pyomo.common.timing import HierarchicalTimer 17 18 19class BaseInteriorPointInterface(object, metaclass=ABCMeta): 20 @abstractmethod 21 def n_primals(self): 22 pass 23 24 @abstractmethod 25 def nnz_hessian_lag(self): 26 pass 27 28 @abstractmethod 29 def primals_lb(self): 30 pass 31 32 @abstractmethod 33 def primals_ub(self): 34 pass 35 36 @abstractmethod 37 def init_primals(self): 38 pass 39 40 @abstractmethod 41 def set_primals(self, primals): 42 pass 43 44 @abstractmethod 45 def get_primals(self): 46 pass 47 48 @abstractmethod 49 def get_obj_factor(self): 50 pass 51 52 @abstractmethod 53 def set_obj_factor(self, obj_factor): 54 pass 55 56 @abstractmethod 57 def evaluate_objective(self): 58 pass 59 60 @abstractmethod 61 def evaluate_grad_objective(self): 62 pass 63 64 @abstractmethod 65 def n_eq_constraints(self): 66 pass 67 68 @abstractmethod 69 def n_ineq_constraints(self): 70 pass 71 72 @abstractmethod 73 def nnz_jacobian_eq(self): 74 pass 75 76 @abstractmethod 77 def nnz_jacobian_ineq(self): 78 pass 79 80 @abstractmethod 81 def ineq_lb(self): 82 pass 83 84 @abstractmethod 85 def ineq_ub(self): 86 pass 87 88 @abstractmethod 89 def init_duals_eq(self): 90 pass 91 92 @abstractmethod 93 def init_duals_ineq(self): 94 pass 95 96 @abstractmethod 97 def set_duals_eq(self, duals_eq): 98 pass 99 100 @abstractmethod 101 def set_duals_ineq(self, duals_ineq): 102 pass 103 104 @abstractmethod 105 def get_duals_eq(self): 106 pass 107 108 @abstractmethod 109 def get_duals_ineq(self): 110 pass 111 112 @abstractmethod 113 def evaluate_eq_constraints(self): 114 pass 115 116 @abstractmethod 117 def evaluate_ineq_constraints(self): 118 pass 119 120 @abstractmethod 121 def evaluate_jacobian_eq(self): 122 pass 123 124 @abstractmethod 125 def evaluate_jacobian_ineq(self): 126 pass 127 128 @abstractmethod 129 def init_slacks(self): 130 pass 131 132 @abstractmethod 133 def init_duals_primals_lb(self): 134 pass 135 136 @abstractmethod 137 def init_duals_primals_ub(self): 138 pass 139 140 @abstractmethod 141 def init_duals_slacks_lb(self): 142 pass 143 144 @abstractmethod 145 def init_duals_slacks_ub(self): 146 pass 147 148 @abstractmethod 149 def set_slacks(self, slacks): 150 pass 151 152 @abstractmethod 153 def set_duals_primals_lb(self, duals): 154 pass 155 156 @abstractmethod 157 def set_duals_primals_ub(self, duals): 158 pass 159 160 @abstractmethod 161 def set_duals_slacks_lb(self, duals): 162 pass 163 164 @abstractmethod 165 def set_duals_slacks_ub(self, duals): 166 pass 167 168 @abstractmethod 169 def get_slacks(self): 170 pass 171 172 @abstractmethod 173 def get_duals_primals_lb(self): 174 pass 175 176 @abstractmethod 177 def get_duals_primals_ub(self): 178 pass 179 180 @abstractmethod 181 def get_duals_slacks_lb(self): 182 pass 183 184 @abstractmethod 185 def get_duals_slacks_ub(self): 186 pass 187 188 @abstractmethod 189 def set_barrier_parameter(self, barrier): 190 pass 191 192 @abstractmethod 193 def evaluate_primal_dual_kkt_matrix(self, timer=None): 194 pass 195 196 @abstractmethod 197 def evaluate_primal_dual_kkt_rhs(self, timer=None): 198 pass 199 200 @abstractmethod 201 def set_primal_dual_kkt_solution(self, sol): 202 pass 203 204 @abstractmethod 205 def get_delta_primals(self): 206 pass 207 208 @abstractmethod 209 def get_delta_slacks(self): 210 pass 211 212 @abstractmethod 213 def get_delta_duals_eq(self): 214 pass 215 216 @abstractmethod 217 def get_delta_duals_ineq(self): 218 pass 219 220 @abstractmethod 221 def get_delta_duals_primals_lb(self): 222 pass 223 224 @abstractmethod 225 def get_delta_duals_primals_ub(self): 226 pass 227 228 @abstractmethod 229 def get_delta_duals_slacks_lb(self): 230 pass 231 232 @abstractmethod 233 def get_delta_duals_slacks_ub(self): 234 pass 235 236 def regularize_equality_gradient(self, kkt, coef, copy_kkt=True): 237 raise RuntimeError( 238 'Equality gradient regularization is necessary but no ' 239 'function has been implemented for doing so.') 240 241 def regularize_hessian(self, kkt, coef, copy_kkt=True): 242 raise RuntimeError( 243 'Hessian of Lagrangian regularization is necessary but no ' 244 'function has been implemented for doing so.') 245 246 247class InteriorPointInterface(BaseInteriorPointInterface): 248 def __init__(self, pyomo_model): 249 if type(pyomo_model) is str: 250 # Assume argument is the name of an nl file 251 self._nlp = ampl_nlp.AmplNLP(pyomo_model) 252 else: 253 self._nlp = pyomo_nlp.PyomoNLP(pyomo_model) 254 self._slacks = self.init_slacks() 255 256 # set the init_duals_primals_lb/ub from ipopt_zL_out, ipopt_zU_out if available 257 # need to compress them as well and initialize the duals_primals_lb/ub 258 self._init_duals_primals_lb, self._init_duals_primals_ub =\ 259 self._get_full_duals_primals_bounds() 260 self._init_duals_primals_lb[np.isneginf(self._nlp.primals_lb())] = 0 261 self._init_duals_primals_ub[np.isinf(self._nlp.primals_ub())] = 0 262 self._duals_primals_lb = self._init_duals_primals_lb.copy() 263 self._duals_primals_ub = self._init_duals_primals_ub.copy() 264 265 # set the init_duals_slacks_lb/ub from the init_duals_ineq 266 # need to be compressed and set according to their sign 267 # (-) value indicates it the upper is active, while (+) indicates 268 # that lower is active 269 self._init_duals_slacks_lb = self._nlp.init_duals_ineq().copy() 270 self._init_duals_slacks_lb[self._init_duals_slacks_lb < 0] = 0 271 self._init_duals_slacks_ub = self._nlp.init_duals_ineq().copy() 272 self._init_duals_slacks_ub[self._init_duals_slacks_ub > 0] = 0 273 self._init_duals_slacks_ub *= -1.0 274 275 self._duals_slacks_lb = self._init_duals_slacks_lb.copy() 276 self._duals_slacks_ub = self._init_duals_slacks_ub.copy() 277 278 self._delta_primals = None 279 self._delta_slacks = None 280 self._delta_duals_eq = None 281 self._delta_duals_ineq = None 282 self._barrier = None 283 284 def n_primals(self): 285 return self._nlp.n_primals() 286 287 def nnz_hessian_lag(self): 288 return self._nlp.nnz_hessian_lag() 289 290 def set_obj_factor(self, obj_factor): 291 self._nlp.set_obj_factor(obj_factor) 292 293 def get_obj_factor(self): 294 return self._nlp.get_obj_factor() 295 296 def n_eq_constraints(self): 297 return self._nlp.n_eq_constraints() 298 299 def n_ineq_constraints(self): 300 return self._nlp.n_ineq_constraints() 301 302 def nnz_jacobian_eq(self): 303 return self._nlp.nnz_jacobian_eq() 304 305 def nnz_jacobian_ineq(self): 306 return self._nlp.nnz_jacobian_ineq() 307 308 def init_primals(self): 309 primals = self._nlp.init_primals() 310 return primals 311 312 def init_slacks(self): 313 slacks = self._nlp.evaluate_ineq_constraints() 314 return slacks 315 316 def init_duals_eq(self): 317 return self._nlp.init_duals_eq() 318 319 def init_duals_ineq(self): 320 return self._nlp.init_duals_ineq() 321 322 def init_duals_primals_lb(self): 323 return self._init_duals_primals_lb 324 325 def init_duals_primals_ub(self): 326 return self._init_duals_primals_ub 327 328 def init_duals_slacks_lb(self): 329 return self._init_duals_slacks_lb 330 331 def init_duals_slacks_ub(self): 332 return self._init_duals_slacks_ub 333 334 def set_primals(self, primals): 335 self._nlp.set_primals(primals) 336 337 def set_slacks(self, slacks): 338 self._slacks = slacks 339 340 def set_duals_eq(self, duals): 341 self._nlp.set_duals_eq(duals) 342 343 def set_duals_ineq(self, duals): 344 self._nlp.set_duals_ineq(duals) 345 346 def set_duals_primals_lb(self, duals): 347 self._duals_primals_lb = duals 348 349 def set_duals_primals_ub(self, duals): 350 self._duals_primals_ub = duals 351 352 def set_duals_slacks_lb(self, duals): 353 self._duals_slacks_lb = duals 354 355 def set_duals_slacks_ub(self, duals): 356 self._duals_slacks_ub = duals 357 358 def get_primals(self): 359 return self._nlp.get_primals() 360 361 def get_slacks(self): 362 return self._slacks 363 364 def get_duals_eq(self): 365 return self._nlp.get_duals_eq() 366 367 def get_duals_ineq(self): 368 return self._nlp.get_duals_ineq() 369 370 def get_duals_primals_lb(self): 371 return self._duals_primals_lb 372 373 def get_duals_primals_ub(self): 374 return self._duals_primals_ub 375 376 def get_duals_slacks_lb(self): 377 return self._duals_slacks_lb 378 379 def get_duals_slacks_ub(self): 380 return self._duals_slacks_ub 381 382 def primals_lb(self): 383 return self._nlp.primals_lb() 384 385 def primals_ub(self): 386 return self._nlp.primals_ub() 387 388 def ineq_lb(self): 389 return self._nlp.ineq_lb() 390 391 def ineq_ub(self): 392 return self._nlp.ineq_ub() 393 394 def set_barrier_parameter(self, barrier): 395 self._barrier = barrier 396 397 def pyomo_nlp(self): 398 return self._nlp 399 400 def evaluate_primal_dual_kkt_matrix(self, timer=None): 401 if timer is None: 402 timer = HierarchicalTimer() 403 timer.start('eval hess') 404 hessian = self._nlp.evaluate_hessian_lag() 405 timer.stop('eval hess') 406 timer.start('eval jac') 407 jac_eq = self._nlp.evaluate_jacobian_eq() 408 jac_ineq = self._nlp.evaluate_jacobian_ineq() 409 timer.stop('eval jac') 410 411 duals_primals_lb = self._duals_primals_lb 412 duals_primals_ub = self._duals_primals_ub 413 duals_slacks_lb = self._duals_slacks_lb 414 duals_slacks_ub = self._duals_slacks_ub 415 primals = self._nlp.get_primals() 416 417 timer.start('hess block') 418 data = (duals_primals_lb/(primals - self._nlp.primals_lb()) + 419 duals_primals_ub/(self._nlp.primals_ub() - primals)) 420 n = self._nlp.n_primals() 421 indices = np.arange(n) 422 hess_block = scipy.sparse.coo_matrix((data, (indices, indices)), shape=(n, n)) 423 hess_block += hessian 424 timer.stop('hess block') 425 426 timer.start('slack block') 427 data = (duals_slacks_lb/(self._slacks - self._nlp.ineq_lb()) + 428 duals_slacks_ub/(self._nlp.ineq_ub() - self._slacks)) 429 n = self._nlp.n_ineq_constraints() 430 indices = np.arange(n) 431 slack_block = scipy.sparse.coo_matrix((data, (indices, indices)), shape=(n, n)) 432 timer.stop('slack block') 433 434 timer.start('set block') 435 kkt = BlockMatrix(4, 4) 436 kkt.set_block(0, 0, hess_block) 437 kkt.set_block(1, 1, slack_block) 438 kkt.set_block(2, 0, jac_eq) 439 kkt.set_block(0, 2, jac_eq.transpose()) 440 kkt.set_block(3, 0, jac_ineq) 441 kkt.set_block(0, 3, jac_ineq.transpose()) 442 kkt.set_block(3, 1, -scipy.sparse.identity( 443 self._nlp.n_ineq_constraints(), 444 format='coo')) 445 kkt.set_block(1, 3, -scipy.sparse.identity( 446 self._nlp.n_ineq_constraints(), 447 format='coo')) 448 timer.stop('set block') 449 return kkt 450 451 def evaluate_primal_dual_kkt_rhs(self, timer=None): 452 if timer is None: 453 timer = HierarchicalTimer() 454 timer.start('eval grad obj') 455 grad_obj = self.get_obj_factor() * self.evaluate_grad_objective() 456 timer.stop('eval grad obj') 457 timer.start('eval jac') 458 jac_eq = self._nlp.evaluate_jacobian_eq() 459 jac_ineq = self._nlp.evaluate_jacobian_ineq() 460 timer.stop('eval jac') 461 timer.start('eval cons') 462 eq_resid = self._nlp.evaluate_eq_constraints() 463 ineq_resid = self._nlp.evaluate_ineq_constraints() - self._slacks 464 timer.stop('eval cons') 465 466 timer.start('grad_lag_primals') 467 grad_lag_primals = (grad_obj + 468 jac_eq.transpose() * self._nlp.get_duals_eq() + 469 jac_ineq.transpose() * self._nlp.get_duals_ineq() - 470 self._barrier / (self._nlp.get_primals() - self._nlp.primals_lb()) + 471 self._barrier / (self._nlp.primals_ub() - self._nlp.get_primals())) 472 timer.stop('grad_lag_primals') 473 474 timer.start('grad_lag_slacks') 475 grad_lag_slacks = (-self._nlp.get_duals_ineq() - 476 self._barrier / (self._slacks - self._nlp.ineq_lb()) + 477 self._barrier / (self._nlp.ineq_ub() - self._slacks)) 478 timer.stop('grad_lag_slacks') 479 480 rhs = BlockVector(4) 481 rhs.set_block(0, grad_lag_primals) 482 rhs.set_block(1, grad_lag_slacks) 483 rhs.set_block(2, eq_resid) 484 rhs.set_block(3, ineq_resid) 485 rhs = -rhs 486 return rhs 487 488 def set_primal_dual_kkt_solution(self, sol): 489 self._delta_primals = sol.get_block(0) 490 self._delta_slacks = sol.get_block(1) 491 self._delta_duals_eq = sol.get_block(2) 492 self._delta_duals_ineq = sol.get_block(3) 493 494 def get_delta_primals(self): 495 return self._delta_primals 496 497 def get_delta_slacks(self): 498 return self._delta_slacks 499 500 def get_delta_duals_eq(self): 501 return self._delta_duals_eq 502 503 def get_delta_duals_ineq(self): 504 return self._delta_duals_ineq 505 506 def get_delta_duals_primals_lb(self): 507 res = (((self._barrier - self._duals_primals_lb * self._delta_primals) / 508 (self._nlp.get_primals() - self._nlp.primals_lb())) - 509 self._duals_primals_lb) 510 return res 511 512 def get_delta_duals_primals_ub(self): 513 res = (((self._barrier + self._duals_primals_ub * self._delta_primals) / 514 (self._nlp.primals_ub() - self._nlp.get_primals())) - 515 self._duals_primals_ub) 516 return res 517 518 def get_delta_duals_slacks_lb(self): 519 res = (((self._barrier - self._duals_slacks_lb * self._delta_slacks) / 520 (self._slacks - self._nlp.ineq_lb())) - 521 self._duals_slacks_lb) 522 return res 523 524 def get_delta_duals_slacks_ub(self): 525 res = (((self._barrier + self._duals_slacks_ub * self._delta_slacks) / 526 (self._nlp.ineq_ub() - self._slacks)) - 527 self._duals_slacks_ub) 528 return res 529 530 def evaluate_objective(self): 531 return self._nlp.evaluate_objective() 532 533 def evaluate_eq_constraints(self): 534 return self._nlp.evaluate_eq_constraints() 535 536 def evaluate_ineq_constraints(self): 537 return self._nlp.evaluate_ineq_constraints() 538 539 def evaluate_grad_objective(self): 540 return self._nlp.evaluate_grad_objective() 541 542 def evaluate_jacobian_eq(self): 543 return self._nlp.evaluate_jacobian_eq() 544 545 def evaluate_jacobian_ineq(self): 546 return self._nlp.evaluate_jacobian_ineq() 547 548 def regularize_equality_gradient(self, kkt, coef, copy_kkt=True): 549 # Not technically regularizing the equality gradient ... 550 # Replace this with a regularize_diagonal_block function? 551 # Then call with kkt matrix and the value of the perturbation? 552 553 # Use a constant perturbation to regularize the equality constraint 554 # gradient 555 if copy_kkt: 556 kkt = kkt.copy() 557 reg_coef = coef 558 ptb = (reg_coef * 559 scipy.sparse.identity(self._nlp.n_eq_constraints(), 560 format='coo')) 561 562 kkt.set_block(2, 2, ptb) 563 return kkt 564 565 def regularize_hessian(self, kkt, coef, copy_kkt=True): 566 if copy_kkt: 567 kkt = kkt.copy() 568 569 hess = kkt.get_block(0, 0) 570 ptb = coef * scipy.sparse.identity(self._nlp.n_primals(), format='coo') 571 hess += ptb 572 kkt.set_block(0, 0, hess) 573 return kkt 574 575 def _get_full_duals_primals_bounds(self): 576 full_duals_primals_lb = None 577 full_duals_primals_ub = None 578 # Check in case _nlp was constructed as an AmplNLP (from an nl file) 579 if (hasattr(self._nlp, 'pyomo_model') and 580 hasattr(self._nlp, 'get_pyomo_variables')): 581 pyomo_model = self._nlp.pyomo_model() 582 pyomo_variables = self._nlp.get_pyomo_variables() 583 if hasattr(pyomo_model,'ipopt_zL_out'): 584 zL_suffix = pyomo_model.ipopt_zL_out 585 full_duals_primals_lb = np.empty(self._nlp.n_primals()) 586 for i,v in enumerate(pyomo_variables): 587 if v in zL_suffix: 588 full_duals_primals_lb[i] = zL_suffix[v] 589 590 if hasattr(pyomo_model,'ipopt_zU_out'): 591 zU_suffix = pyomo_model.ipopt_zU_out 592 full_duals_primals_ub = np.empty(self._nlp.n_primals()) 593 for i,v in enumerate(pyomo_variables): 594 if v in zU_suffix: 595 full_duals_primals_ub[i] = zU_suffix[v] 596 597 if full_duals_primals_lb is None: 598 full_duals_primals_lb = np.ones(self._nlp.n_primals()) 599 600 if full_duals_primals_ub is None: 601 full_duals_primals_ub = np.ones(self._nlp.n_primals()) 602 603 return full_duals_primals_lb, full_duals_primals_ub 604 605 def load_primals_into_pyomo_model(self): 606 if not isinstance(self._nlp, pyomo_nlp.PyomoNLP): 607 raise RuntimeError('Can only load primals into a pyomo model if a pyomo model was used in the constructor.') 608 609 pyomo_variables = self._nlp.get_pyomo_variables() 610 primals = self._nlp.get_primals() 611 for i, v in enumerate(pyomo_variables): 612 v.value = primals[i] 613 614 def pyomo_model(self): 615 return self._nlp.pyomo_model() 616 617 def get_pyomo_variables(self): 618 return self._nlp.get_pyomo_variables() 619 620 def get_pyomo_constraints(self): 621 return self._nlp.get_pyomo_constraints() 622 623 def variable_names(self): 624 return self._nlp.variable_names() 625 626 def constraint_names(self): 627 return self._nlp.constraint_names() 628 629 def get_primal_indices(self, pyomo_variables): 630 return self._nlp.get_primal_indices(pyomo_variables) 631 632 def get_constraint_indices(self, pyomo_constraints): 633 return self._nlp.get_constraint_indices(pyomo_constraints) 634