1 // Copyright (C) 2004, 2011 International Business Machines and others. 2 // All Rights Reserved. 3 // This code is published under the Eclipse Public License. 4 // 5 // Authors: Carl Laird, Andreas Waechter IBM 2004-08-13 6 7 #ifndef __IPIPOPTCALCULATEDQUANTITIES_HPP__ 8 #define __IPIPOPTCALCULATEDQUANTITIES_HPP__ 9 10 #include "IpSmartPtr.hpp" 11 #include "IpCachedResults.hpp" 12 13 #include <string> 14 15 namespace Ipopt 16 { 17 class IpoptNLP; 18 class IpoptData; 19 class Vector; 20 class Matrix; 21 class SymMatrix; 22 class Journalist; 23 class OptionsList; 24 class RegisteredOptions; 25 26 /** Norm types */ 27 enum ENormType 28 { 29 NORM_1 = 0, 30 NORM_2, 31 NORM_MAX 32 }; 33 34 /** Base class for additional calculated quantities that is special 35 * to a particular type of algorithm, such as the CG penalty 36 * function, or using iterative linear solvers. 37 * 38 * The regular IpoptCalculatedQuantities object should be given 39 * a derivation of this base class when it is created. 40 */ 41 class IPOPTLIB_EXPORT IpoptAdditionalCq: public ReferencedObject 42 { 43 public: 44 /**@name Constructors/Destructors */ 45 //@{ 46 /** Default Constructor */ IpoptAdditionalCq()47 IpoptAdditionalCq() 48 { } 49 50 /** Destructor */ ~IpoptAdditionalCq()51 virtual ~IpoptAdditionalCq() 52 { } 53 //@} 54 55 /** This method is called to initialize the global algorithmic 56 * parameters. 57 * 58 * The parameters are taken from the OptionsList object. 59 */ 60 virtual bool Initialize( 61 const Journalist& jnlst, 62 const OptionsList& options, 63 const std::string& prefix 64 ) = 0; 65 66 private: 67 /**@name Default Compiler Generated Methods 68 * (Hidden to avoid implicit creation/calling). 69 * 70 * These methods are not implemented and 71 * we do not want the compiler to implement 72 * them for us, so we declare them private 73 * and do not define them. This ensures that 74 * they will not be implicitly created/called. 75 */ 76 //@{ 77 /** Copy Constructor */ 78 IpoptAdditionalCq( 79 const IpoptAdditionalCq&); 80 81 /** Default Assignment Operator */ 82 void operator=( 83 const IpoptAdditionalCq&); 84 //@} 85 }; 86 87 /** Class for all IPOPT specific calculated quantities. */ 88 class IPOPTLIB_EXPORT IpoptCalculatedQuantities: public ReferencedObject 89 { 90 public: 91 92 /**@name Constructors/Destructors */ 93 //@{ 94 /** Constructor */ 95 IpoptCalculatedQuantities( 96 const SmartPtr<IpoptNLP>& ip_nlp, 97 const SmartPtr<IpoptData>& ip_data 98 ); 99 /** Destructor */ 100 virtual ~IpoptCalculatedQuantities(); 101 //@} 102 103 /** Method for setting pointer for additional calculated 104 * quantities. 105 * 106 * This needs to be called before Initialized. 107 */ SetAddCq(SmartPtr<IpoptAdditionalCq> add_cq)108 void SetAddCq( 109 SmartPtr<IpoptAdditionalCq> add_cq 110 ) 111 { 112 DBG_ASSERT(!HaveAddCq()); 113 add_cq_ = add_cq; 114 } 115 116 /** Method detecting if additional object for calculated 117 * quantities has already been set 118 */ HaveAddCq()119 bool HaveAddCq() 120 { 121 return IsValid(add_cq_); 122 } 123 124 /** This method must be called to initialize the global 125 * algorithmic parameters. 126 * 127 * The parameters are taken from the OptionsList object. 128 */ 129 bool Initialize( 130 const Journalist& jnlst, 131 const OptionsList& options, 132 const std::string& prefix 133 ); 134 135 /** @name Slacks */ 136 //@{ 137 /** Slacks for x_L (at current iterate) */ 138 SmartPtr<const Vector> curr_slack_x_L(); 139 /** Slacks for x_U (at current iterate) */ 140 SmartPtr<const Vector> curr_slack_x_U(); 141 /** Slacks for s_L (at current iterate) */ 142 SmartPtr<const Vector> curr_slack_s_L(); 143 /** Slacks for s_U (at current iterate) */ 144 SmartPtr<const Vector> curr_slack_s_U(); 145 /** Slacks for x_L (at trial point) */ 146 SmartPtr<const Vector> trial_slack_x_L(); 147 /** Slacks for x_U (at trial point) */ 148 SmartPtr<const Vector> trial_slack_x_U(); 149 /** Slacks for s_L (at trial point) */ 150 SmartPtr<const Vector> trial_slack_s_L(); 151 /** Slacks for s_U (at trial point) */ 152 SmartPtr<const Vector> trial_slack_s_U(); 153 /** Indicating whether or not we "fudged" the slacks */ 154 Index AdjustedTrialSlacks(); 155 /** Reset the flags for "fudged" slacks */ 156 void ResetAdjustedTrialSlacks(); 157 //@} 158 159 /** @name Objective function */ 160 //@{ 161 /** Value of objective function (at current point) */ 162 virtual Number curr_f(); 163 /** Unscaled value of the objective function (at the current point) */ 164 virtual Number unscaled_curr_f(); 165 /** Value of objective function (at trial point) */ 166 virtual Number trial_f(); 167 /** Unscaled value of the objective function (at the trial point) */ 168 virtual Number unscaled_trial_f(); 169 /** Gradient of objective function (at current point) */ 170 SmartPtr<const Vector> curr_grad_f(); 171 /** Gradient of objective function (at trial point) */ 172 SmartPtr<const Vector> trial_grad_f(); 173 //@} 174 175 /** @name Barrier Objective Function */ 176 //@{ 177 /** Barrier Objective Function Value 178 * (at current iterate with current mu) 179 */ 180 virtual Number curr_barrier_obj(); 181 /** Barrier Objective Function Value 182 * (at trial point with current mu) 183 */ 184 virtual Number trial_barrier_obj(); 185 186 /** Gradient of barrier objective function with respect to x 187 * (at current point with current mu) 188 */ 189 SmartPtr<const Vector> curr_grad_barrier_obj_x(); 190 /** Gradient of barrier objective function with respect to s 191 * (at current point with current mu) 192 */ 193 SmartPtr<const Vector> curr_grad_barrier_obj_s(); 194 195 /** Gradient of the damping term with respect to x (times kappa_d) */ 196 SmartPtr<const Vector> grad_kappa_times_damping_x(); 197 /** Gradient of the damping term with respect to s (times kappa_d) */ 198 SmartPtr<const Vector> grad_kappa_times_damping_s(); 199 //@} 200 201 /** @name Constraints */ 202 //@{ 203 /** c(x) (at current point) */ 204 SmartPtr<const Vector> curr_c(); 205 /** unscaled c(x) (at current point) */ 206 SmartPtr<const Vector> unscaled_curr_c(); 207 /** c(x) (at trial point) */ 208 SmartPtr<const Vector> trial_c(); 209 /** unscaled c(x) (at trial point) */ 210 SmartPtr<const Vector> unscaled_trial_c(); 211 /** d(x) (at current point) */ 212 SmartPtr<const Vector> curr_d(); 213 /** unscaled d(x) (at current point) */ 214 SmartPtr<const Vector> unscaled_curr_d(); 215 /** d(x) (at trial point) */ 216 SmartPtr<const Vector> trial_d(); 217 /** d(x) - s (at current point) */ 218 SmartPtr<const Vector> curr_d_minus_s(); 219 /** d(x) - s (at trial point) */ 220 SmartPtr<const Vector> trial_d_minus_s(); 221 /** Jacobian of c (at current point) */ 222 SmartPtr<const Matrix> curr_jac_c(); 223 /** Jacobian of c (at trial point) */ 224 SmartPtr<const Matrix> trial_jac_c(); 225 /** Jacobian of d (at current point) */ 226 SmartPtr<const Matrix> curr_jac_d(); 227 /** Jacobian of d (at trial point) */ 228 SmartPtr<const Matrix> trial_jac_d(); 229 /** Product of Jacobian (evaluated at current point) of C 230 * transpose with general vector 231 */ 232 SmartPtr<const Vector> curr_jac_cT_times_vec( 233 const Vector& vec 234 ); 235 /** Product of Jacobian (evaluated at trial point) of C 236 * transpose with general vector 237 */ 238 SmartPtr<const Vector> trial_jac_cT_times_vec( 239 const Vector& vec 240 ); 241 /** Product of Jacobian (evaluated at current point) of D 242 * transpose with general vector 243 */ 244 SmartPtr<const Vector> curr_jac_dT_times_vec( 245 const Vector& vec 246 ); 247 /** Product of Jacobian (evaluated at trial point) of D 248 * transpose with general vector 249 */ 250 SmartPtr<const Vector> trial_jac_dT_times_vec( 251 const Vector& vec 252 ); 253 /** Product of Jacobian (evaluated at current point) of C 254 * transpose with current y_c 255 */ 256 SmartPtr<const Vector> curr_jac_cT_times_curr_y_c(); 257 /** Product of Jacobian (evaluated at trial point) of C 258 * transpose with trial y_c 259 */ 260 SmartPtr<const Vector> trial_jac_cT_times_trial_y_c(); 261 /** Product of Jacobian (evaluated at current point) of D 262 * transpose with current y_d 263 */ 264 SmartPtr<const Vector> curr_jac_dT_times_curr_y_d(); 265 /** Product of Jacobian (evaluated at trial point) of D 266 * transpose with trial y_d 267 */ 268 SmartPtr<const Vector> trial_jac_dT_times_trial_y_d(); 269 /** Product of Jacobian (evaluated at current point) of C 270 * with general vector 271 */ 272 SmartPtr<const Vector> curr_jac_c_times_vec( 273 const Vector& vec 274 ); 275 /** Product of Jacobian (evaluated at current point) of D 276 * with general vector 277 */ 278 SmartPtr<const Vector> curr_jac_d_times_vec( 279 const Vector& vec 280 ); 281 /** Constraint Violation (at current iterate). 282 * 283 * This value should 284 * be used in the line search, and not curr_primal_infeasibility(). 285 * What type of norm is used depends on constr_viol_normtype 286 */ 287 virtual Number curr_constraint_violation(); 288 /** Constraint Violation (at trial point). 289 * 290 * This value should 291 * be used in the line search, and not curr_primal_infeasibility(). 292 * What type of norm is used depends on constr_viol_normtype 293 */ 294 virtual Number trial_constraint_violation(); 295 /** Real constraint violation in a given norm (at current iterate). 296 * 297 * This considers the inequality constraints without slacks. 298 */ 299 virtual Number curr_nlp_constraint_violation( 300 ENormType NormType 301 ); 302 /** Unscaled real constraint violation in a given norm (at current iterate). 303 * 304 * This considers the inequality constraints without slacks. 305 */ 306 virtual Number unscaled_curr_nlp_constraint_violation( 307 ENormType NormType 308 ); 309 /** Unscaled real constraint violation in a given norm (at trial iterate). 310 * 311 * This considers the inequality constraints without slacks. 312 */ 313 virtual Number unscaled_trial_nlp_constraint_violation( 314 ENormType NormType 315 ); 316 //@} 317 318 /** @name Hessian matrices */ 319 //@{ 320 /** exact Hessian at current iterate (uncached) */ 321 SmartPtr<const SymMatrix> curr_exact_hessian(); 322 //@} 323 324 /** @name primal-dual error and its components */ 325 //@{ 326 /** x-part of gradient of Lagrangian function (at current point) */ 327 SmartPtr<const Vector> curr_grad_lag_x(); 328 /** x-part of gradient of Lagrangian function (at trial point) */ 329 SmartPtr<const Vector> trial_grad_lag_x(); 330 /** s-part of gradient of Lagrangian function (at current point) */ 331 SmartPtr<const Vector> curr_grad_lag_s(); 332 /** s-part of gradient of Lagrangian function (at trial point) */ 333 SmartPtr<const Vector> trial_grad_lag_s(); 334 /** x-part of gradient of Lagrangian function (at current point) 335 * including linear damping term 336 */ 337 SmartPtr<const Vector> curr_grad_lag_with_damping_x(); 338 /** s-part of gradient of Lagrangian function (at current point) 339 * including linear damping term 340 */ 341 SmartPtr<const Vector> curr_grad_lag_with_damping_s(); 342 /** Complementarity for x_L (for current iterate) */ 343 SmartPtr<const Vector> curr_compl_x_L(); 344 /** Complementarity for x_U (for current iterate) */ 345 SmartPtr<const Vector> curr_compl_x_U(); 346 /** Complementarity for s_L (for current iterate) */ 347 SmartPtr<const Vector> curr_compl_s_L(); 348 /** Complementarity for s_U (for current iterate) */ 349 SmartPtr<const Vector> curr_compl_s_U(); 350 /** Complementarity for x_L (for trial iterate) */ 351 SmartPtr<const Vector> trial_compl_x_L(); 352 /** Complementarity for x_U (for trial iterate) */ 353 SmartPtr<const Vector> trial_compl_x_U(); 354 /** Complementarity for s_L (for trial iterate) */ 355 SmartPtr<const Vector> trial_compl_s_L(); 356 /** Complementarity for s_U (for trial iterate) */ 357 SmartPtr<const Vector> trial_compl_s_U(); 358 /** Relaxed complementarity for x_L (for current iterate and current mu) */ 359 SmartPtr<const Vector> curr_relaxed_compl_x_L(); 360 /** Relaxed complementarity for x_U (for current iterate and current mu) */ 361 SmartPtr<const Vector> curr_relaxed_compl_x_U(); 362 /** Relaxed complementarity for s_L (for current iterate and current mu) */ 363 SmartPtr<const Vector> curr_relaxed_compl_s_L(); 364 /** Relaxed complementarity for s_U (for current iterate and current mu) */ 365 SmartPtr<const Vector> curr_relaxed_compl_s_U(); 366 367 /** Primal infeasibility in a given norm (at current iterate). */ 368 virtual Number curr_primal_infeasibility( 369 ENormType NormType 370 ); 371 /** Primal infeasibility in a given norm (at trial point) */ 372 virtual Number trial_primal_infeasibility( 373 ENormType NormType 374 ); 375 376 /** Dual infeasibility in a given norm (at current iterate) */ 377 virtual Number curr_dual_infeasibility( 378 ENormType NormType 379 ); 380 /** Dual infeasibility in a given norm (at trial iterate) */ 381 virtual Number trial_dual_infeasibility( 382 ENormType NormType 383 ); 384 /** Unscaled dual infeasibility in a given norm (at current iterate) */ 385 virtual Number unscaled_curr_dual_infeasibility( 386 ENormType NormType 387 ); 388 389 /** Complementarity (for all complementarity conditions together) 390 * in a given norm (at current iterate) 391 */ 392 virtual Number curr_complementarity( 393 Number mu, 394 ENormType NormType 395 ); 396 /** Complementarity (for all complementarity conditions together) 397 * in a given norm (at trial iterate) 398 */ 399 virtual Number trial_complementarity( 400 Number mu, 401 ENormType NormType 402 ); 403 /** Complementarity (for all complementarity conditions together) 404 * in a given norm (at current iterate) without NLP scaling. 405 */ 406 virtual Number unscaled_curr_complementarity( 407 Number mu, 408 ENormType NormType 409 ); 410 411 /** Centrality measure (in spirit of the -infinity-neighborhood. */ 412 Number CalcCentralityMeasure( 413 const Vector& compl_x_L, 414 const Vector& compl_x_U, 415 const Vector& compl_s_L, 416 const Vector& compl_s_U 417 ); 418 /** Centrality measure at current point */ 419 virtual Number curr_centrality_measure(); 420 421 /** Total optimality error for the original NLP at the current 422 * iterate, using scaling factors based on multipliers. 423 * 424 * Note that here the constraint violation is measured without slacks 425 * (nlp_constraint_violation) 426 */ 427 virtual Number curr_nlp_error(); 428 /** Total optimality error for the original NLP at the current 429 * iterate, but using no scaling based on multipliers, and no 430 * scaling for the NLP. 431 * 432 * Note that here the constraint violation 433 * is measured without slacks (nlp_constraint_violation) 434 */ 435 virtual Number unscaled_curr_nlp_error(); 436 437 /** Total optimality error for the barrier problem at the 438 * current iterate, using scaling factors based on multipliers. 439 */ 440 virtual Number curr_barrier_error(); 441 442 /** Norm of the primal-dual system for a given mu (at current iterate). 443 * 444 * The norm is defined as the sum of the 1-norms of 445 * dual infeasibiliy, primal infeasibility, and complementarity, 446 * all divided by the number of elements of the vectors of which 447 * the norm is taken. 448 */ 449 virtual Number curr_primal_dual_system_error( 450 Number mu 451 ); 452 /** Norm of the primal-dual system for a given mu (at trial iterate). 453 * 454 * The norm is defined as the sum of the 1-norms of 455 * dual infeasibiliy, primal infeasibility, and complementarity, 456 * all divided by the number of elements of the vectors of which 457 * the norm is taken. 458 */ 459 virtual Number trial_primal_dual_system_error( 460 Number mu 461 ); 462 //@} 463 464 /** @name Computing fraction-to-the-boundary step sizes */ 465 //@{ 466 /** Fraction to the boundary from (current) primal variables x and s 467 * for a given step 468 */ 469 Number primal_frac_to_the_bound( 470 Number tau, 471 const Vector& delta_x, 472 const Vector& delta_s 473 ); 474 /** Fraction to the boundary from (current) primal variables x and s 475 * for internal (current) step 476 */ 477 Number curr_primal_frac_to_the_bound( 478 Number tau 479 ); 480 /** Fraction to the boundary from (current) dual variables z and v 481 * for a given step 482 */ 483 Number dual_frac_to_the_bound( 484 Number tau, 485 const Vector& delta_z_L, 486 const Vector& delta_z_U, 487 const Vector& delta_v_L, 488 const Vector& delta_v_U 489 ); 490 /** Fraction to the boundary from (current) dual variables z and v 491 * for a given step, without caching 492 */ 493 Number uncached_dual_frac_to_the_bound( 494 Number tau, 495 const Vector& delta_z_L, 496 const Vector& delta_z_U, 497 const Vector& delta_v_L, 498 const Vector& delta_v_U 499 ); 500 /** Fraction to the boundary from (current) dual variables z and v 501 * for internal (current) step 502 */ 503 Number curr_dual_frac_to_the_bound( 504 Number tau 505 ); 506 /** Fraction to the boundary from (current) slacks for a given 507 * step in the slacks. 508 * 509 * Usually, one will use the 510 * primal_frac_to_the_bound method to compute the primal fraction 511 * to the boundary step size, but if it is cheaper to provide the 512 * steps in the slacks directly (e.g. when the primal step sizes 513 * are only temporary), the this method is more efficient. This 514 * method does not cache computations. 515 */ 516 Number uncached_slack_frac_to_the_bound( 517 Number tau, 518 const Vector& delta_x_L, 519 const Vector& delta_x_U, 520 const Vector& delta_s_L, 521 const Vector& delta_s_U 522 ); 523 //@} 524 525 /** @name Sigma matrices */ 526 //@{ 527 SmartPtr<const Vector> curr_sigma_x(); 528 SmartPtr<const Vector> curr_sigma_s(); 529 //@} 530 531 /** average of current values of the complementarities */ 532 Number curr_avrg_compl(); 533 /** average of trial values of the complementarities */ 534 Number trial_avrg_compl(); 535 536 /** inner_product of current barrier obj. fn. gradient with 537 * current search direction 538 */ 539 Number curr_gradBarrTDelta(); 540 541 /** Compute the norm of a specific type of a set of vectors (uncached) */ 542 Number 543 CalcNormOfType( 544 ENormType NormType, 545 std::vector<SmartPtr<const Vector> > vecs 546 ); 547 548 /** Compute the norm of a specific type of two vectors (uncached) */ 549 Number 550 CalcNormOfType( 551 ENormType NormType, 552 const Vector& vec1, 553 const Vector& vec2 554 ); 555 556 /** Norm type used for calculating constraint violation */ constr_viol_normtype() const557 ENormType constr_viol_normtype() const 558 { 559 return constr_viol_normtype_; 560 } 561 562 /** Method returning true if this is a square problem */ 563 bool IsSquareProblem() const; 564 565 /** Method returning the IpoptNLP object. 566 * 567 * This should only be used with care! 568 */ GetIpoptNLP()569 SmartPtr<IpoptNLP>& GetIpoptNLP() 570 { 571 return ip_nlp_; 572 } 573 AdditionalCq()574 IpoptAdditionalCq& AdditionalCq() 575 { 576 DBG_ASSERT(IsValid(add_cq_)); 577 return *add_cq_; 578 } 579 580 /** Called by IpoptType to register the options */ 581 static void RegisterOptions( 582 SmartPtr<RegisteredOptions> roptions 583 ); 584 585 private: 586 /**@name Default Compiler Generated Methods 587 * (Hidden to avoid implicit creation/calling). 588 * 589 * These methods are not implemented and 590 * we do not want the compiler to implement 591 * them for us, so we declare them private 592 * and do not define them. This ensures that 593 * they will not be implicitly created/called. 594 */ 595 //@{ 596 /** Default Constructor */ 597 IpoptCalculatedQuantities(); 598 599 /** Copy Constructor */ 600 IpoptCalculatedQuantities( 601 const IpoptCalculatedQuantities& 602 ); 603 604 /** Default Assignment Operator */ 605 void operator=( 606 const IpoptCalculatedQuantities& 607 ); 608 //@} 609 610 /** @name Pointers for easy access to data and NLP information */ 611 //@{ 612 /** Ipopt NLP object */ 613 SmartPtr<IpoptNLP> ip_nlp_; 614 /** Ipopt Data object */ 615 SmartPtr<IpoptData> ip_data_; 616 /** Chen-Goldfarb specific calculated quantities */ 617 SmartPtr<IpoptAdditionalCq> add_cq_; 618 //@} 619 620 /** @name Algorithmic Parameters that can be set through the 621 * options list. 622 * 623 * Those parameters are initialize by calling the Initialize method. 624 */ 625 //@{ 626 /** Parameter in formula for computing overall primal-dual 627 * optimality error 628 */ 629 Number s_max_; 630 /** Weighting factor for the linear damping term added to the 631 * barrier objective function. 632 */ 633 Number kappa_d_; 634 /** fractional movement allowed in bounds */ 635 Number slack_move_; 636 /** Norm type to be used when calculating the constraint violation */ 637 ENormType constr_viol_normtype_; 638 /** Flag indicating whether the TNLP with identical structure has 639 * already been solved before. 640 */ 641 bool warm_start_same_structure_; 642 /** Desired value of the barrier parameter */ 643 Number mu_target_; 644 //@} 645 646 /** @name Caches for slacks */ 647 //@{ 648 CachedResults<SmartPtr<Vector> > curr_slack_x_L_cache_; 649 CachedResults<SmartPtr<Vector> > curr_slack_x_U_cache_; 650 CachedResults<SmartPtr<Vector> > curr_slack_s_L_cache_; 651 CachedResults<SmartPtr<Vector> > curr_slack_s_U_cache_; 652 CachedResults<SmartPtr<Vector> > trial_slack_x_L_cache_; 653 CachedResults<SmartPtr<Vector> > trial_slack_x_U_cache_; 654 CachedResults<SmartPtr<Vector> > trial_slack_s_L_cache_; 655 CachedResults<SmartPtr<Vector> > trial_slack_s_U_cache_; 656 Index num_adjusted_slack_x_L_; 657 Index num_adjusted_slack_x_U_; 658 Index num_adjusted_slack_s_L_; 659 Index num_adjusted_slack_s_U_; 660 //@} 661 662 /** @name Cached for objective function stuff */ 663 //@{ 664 CachedResults<Number> curr_f_cache_; 665 CachedResults<Number> trial_f_cache_; 666 CachedResults<SmartPtr<const Vector> > curr_grad_f_cache_; 667 CachedResults<SmartPtr<const Vector> > trial_grad_f_cache_; 668 //@} 669 670 /** @name Caches for barrier function stuff */ 671 //@{ 672 CachedResults<Number> curr_barrier_obj_cache_; 673 CachedResults<Number> trial_barrier_obj_cache_; 674 CachedResults<SmartPtr<const Vector> > curr_grad_barrier_obj_x_cache_; 675 CachedResults<SmartPtr<const Vector> > curr_grad_barrier_obj_s_cache_; 676 CachedResults<SmartPtr<const Vector> > grad_kappa_times_damping_x_cache_; 677 CachedResults<SmartPtr<const Vector> > grad_kappa_times_damping_s_cache_; 678 //@} 679 680 /** @name Caches for constraint stuff */ 681 //@{ 682 CachedResults<SmartPtr<const Vector> > curr_c_cache_; 683 CachedResults<SmartPtr<const Vector> > trial_c_cache_; 684 CachedResults<SmartPtr<const Vector> > curr_d_cache_; 685 CachedResults<SmartPtr<const Vector> > trial_d_cache_; 686 CachedResults<SmartPtr<const Vector> > curr_d_minus_s_cache_; 687 CachedResults<SmartPtr<const Vector> > trial_d_minus_s_cache_; 688 CachedResults<SmartPtr<const Matrix> > curr_jac_c_cache_; 689 CachedResults<SmartPtr<const Matrix> > trial_jac_c_cache_; 690 CachedResults<SmartPtr<const Matrix> > curr_jac_d_cache_; 691 CachedResults<SmartPtr<const Matrix> > trial_jac_d_cache_; 692 CachedResults<SmartPtr<const Vector> > curr_jac_cT_times_vec_cache_; 693 CachedResults<SmartPtr<const Vector> > trial_jac_cT_times_vec_cache_; 694 CachedResults<SmartPtr<const Vector> > curr_jac_dT_times_vec_cache_; 695 CachedResults<SmartPtr<const Vector> > trial_jac_dT_times_vec_cache_; 696 CachedResults<SmartPtr<const Vector> > curr_jac_c_times_vec_cache_; 697 CachedResults<SmartPtr<const Vector> > curr_jac_d_times_vec_cache_; 698 CachedResults<Number> curr_constraint_violation_cache_; 699 CachedResults<Number> trial_constraint_violation_cache_; 700 CachedResults<Number> curr_nlp_constraint_violation_cache_; 701 CachedResults<Number> unscaled_curr_nlp_constraint_violation_cache_; 702 CachedResults<Number> unscaled_trial_nlp_constraint_violation_cache_; 703 //@} 704 705 /** Cache for the exact Hessian */ 706 CachedResults<SmartPtr<const SymMatrix> > curr_exact_hessian_cache_; 707 708 /** @name Components of primal-dual error */ 709 //@{ 710 CachedResults<SmartPtr<const Vector> > curr_grad_lag_x_cache_; 711 CachedResults<SmartPtr<const Vector> > trial_grad_lag_x_cache_; 712 CachedResults<SmartPtr<const Vector> > curr_grad_lag_s_cache_; 713 CachedResults<SmartPtr<const Vector> > trial_grad_lag_s_cache_; 714 CachedResults<SmartPtr<const Vector> > curr_grad_lag_with_damping_x_cache_; 715 CachedResults<SmartPtr<const Vector> > curr_grad_lag_with_damping_s_cache_; 716 CachedResults<SmartPtr<const Vector> > curr_compl_x_L_cache_; 717 CachedResults<SmartPtr<const Vector> > curr_compl_x_U_cache_; 718 CachedResults<SmartPtr<const Vector> > curr_compl_s_L_cache_; 719 CachedResults<SmartPtr<const Vector> > curr_compl_s_U_cache_; 720 CachedResults<SmartPtr<const Vector> > trial_compl_x_L_cache_; 721 CachedResults<SmartPtr<const Vector> > trial_compl_x_U_cache_; 722 CachedResults<SmartPtr<const Vector> > trial_compl_s_L_cache_; 723 CachedResults<SmartPtr<const Vector> > trial_compl_s_U_cache_; 724 CachedResults<SmartPtr<const Vector> > curr_relaxed_compl_x_L_cache_; 725 CachedResults<SmartPtr<const Vector> > curr_relaxed_compl_x_U_cache_; 726 CachedResults<SmartPtr<const Vector> > curr_relaxed_compl_s_L_cache_; 727 CachedResults<SmartPtr<const Vector> > curr_relaxed_compl_s_U_cache_; 728 CachedResults<Number> curr_primal_infeasibility_cache_; 729 CachedResults<Number> trial_primal_infeasibility_cache_; 730 CachedResults<Number> curr_dual_infeasibility_cache_; 731 CachedResults<Number> trial_dual_infeasibility_cache_; 732 CachedResults<Number> unscaled_curr_dual_infeasibility_cache_; 733 CachedResults<Number> curr_complementarity_cache_; 734 CachedResults<Number> trial_complementarity_cache_; 735 CachedResults<Number> curr_centrality_measure_cache_; 736 CachedResults<Number> curr_nlp_error_cache_; 737 CachedResults<Number> unscaled_curr_nlp_error_cache_; 738 CachedResults<Number> curr_barrier_error_cache_; 739 CachedResults<Number> curr_primal_dual_system_error_cache_; 740 CachedResults<Number> trial_primal_dual_system_error_cache_; 741 //@} 742 743 /** @name Caches for fraction to the boundary step sizes */ 744 //@{ 745 CachedResults<Number> primal_frac_to_the_bound_cache_; 746 CachedResults<Number> dual_frac_to_the_bound_cache_; 747 //@} 748 749 /** @name Caches for sigma matrices */ 750 //@{ 751 CachedResults<SmartPtr<const Vector> > curr_sigma_x_cache_; 752 CachedResults<SmartPtr<const Vector> > curr_sigma_s_cache_; 753 //@} 754 755 /** Cache for average of current complementarity */ 756 CachedResults<Number> curr_avrg_compl_cache_; 757 /** Cache for average of trial complementarity */ 758 CachedResults<Number> trial_avrg_compl_cache_; 759 760 /** Cache for grad barrier obj. fn inner product with step */ 761 CachedResults<Number> curr_gradBarrTDelta_cache_; 762 763 /** @name Indicator vectors required for the linear damping terms 764 * to handle unbounded solution sets. */ 765 //@{ 766 /** Indicator vector for selecting the elements in x that have 767 * only lower bounds. 768 */ 769 SmartPtr<Vector> dampind_x_L_; 770 /** Indicator vector for selecting the elements in x that have 771 * only upper bounds. 772 */ 773 SmartPtr<Vector> dampind_x_U_; 774 /** Indicator vector for selecting the elements in s that have 775 * only lower bounds. 776 */ 777 SmartPtr<Vector> dampind_s_L_; 778 /** Indicator vector for selecting the elements in s that have 779 * only upper bounds. 780 */ 781 SmartPtr<Vector> dampind_s_U_; 782 //@} 783 784 /** @name Temporary vectors for intermediate calculations. 785 * 786 * We keep these around to avoid unnecessarily many new allocations 787 * of Vectors. 788 */ 789 //@{ 790 SmartPtr<Vector> tmp_x_; 791 SmartPtr<Vector> tmp_s_; 792 SmartPtr<Vector> tmp_c_; 793 SmartPtr<Vector> tmp_d_; 794 SmartPtr<Vector> tmp_x_L_; 795 SmartPtr<Vector> tmp_x_U_; 796 SmartPtr<Vector> tmp_s_L_; 797 SmartPtr<Vector> tmp_s_U_; 798 799 /** Accessor methods for the temporary vectors */ 800 Vector& Tmp_x(); 801 Vector& Tmp_s(); 802 Vector& Tmp_c(); 803 Vector& Tmp_d(); 804 Vector& Tmp_x_L(); 805 Vector& Tmp_x_U(); 806 Vector& Tmp_s_L(); 807 Vector& Tmp_s_U(); 808 //@} 809 810 /** flag indicating if Initialize method has been called (for 811 * debugging) 812 */ 813 bool initialize_called_; 814 815 /** @name Auxiliary functions */ 816 //@{ 817 /** Compute new vector containing the slack to a lower bound 818 * (uncached) 819 */ 820 SmartPtr<Vector> CalcSlack_L( 821 const Matrix& P, 822 const Vector& x, 823 const Vector& x_bound 824 ); 825 /** Compute new vector containing the slack to a upper bound 826 * (uncached) 827 */ 828 SmartPtr<Vector> CalcSlack_U( 829 const Matrix& P, 830 const Vector& x, 831 const Vector& x_bound 832 ); 833 /** Compute barrier term at given point 834 * (uncached) 835 */ 836 Number CalcBarrierTerm( 837 Number mu, 838 const Vector& slack_x_L, 839 const Vector& slack_x_U, 840 const Vector& slack_s_L, 841 const Vector& slack_s_U 842 ); 843 844 /** Compute complementarity for slack / multiplier pair */ 845 SmartPtr<const Vector> CalcCompl( 846 const Vector& slack, 847 const Vector& mult 848 ); 849 850 /** Compute fraction to the boundary parameter for lower and upper bounds */ 851 Number CalcFracToBound( 852 const Vector& slack_L, 853 Vector& tmp_L, 854 const Matrix& P_L, 855 const Vector& slack_U, 856 Vector& tmp_U, 857 const Matrix& P_U, 858 const Vector& delta, 859 Number tau 860 ); 861 862 /** Compute the scaling factors for the optimality error. */ 863 void ComputeOptimalityErrorScaling( 864 const Vector& y_c, 865 const Vector& y_d, 866 const Vector& z_L, 867 const Vector& z_U, 868 const Vector& v_L, 869 const Vector& v_U, 870 Number s_max, 871 Number& s_d, 872 Number& s_c 873 ); 874 875 /** Check if slacks are becoming too small. 876 * 877 * If slacks are becoming too small, they are change. 878 * 879 * @return number of corrected slacks 880 */ 881 Index CalculateSafeSlack( 882 SmartPtr<Vector>& slack, 883 const SmartPtr<const Vector>& bound, 884 const SmartPtr<const Vector>& curr_point, 885 const SmartPtr<const Vector>& multiplier 886 ); 887 888 /** Computes the indicator vectors that can be used to filter out 889 * those entries in the slack_... variables, that correspond to 890 * variables with only lower and upper bounds. 891 * 892 * This is required for the linear damping term in the barrier 893 * objective function to handle unbounded solution sets. 894 */ 895 void ComputeDampingIndicators( 896 SmartPtr<const Vector>& dampind_x_L, 897 SmartPtr<const Vector>& dampind_x_U, 898 SmartPtr<const Vector>& dampind_s_L, 899 SmartPtr<const Vector>& dampind_s_U 900 ); 901 902 /** Check if we are in the restoration phase. 903 * 904 * @return true, if the ip_nlp is of the type RestoIpoptNLP 905 * 906 * @todo We probably want to 907 * handle this more elegant and don't have an explicit dependency 908 * here. Now I added this because otherwise the caching doesn't 909 * work properly since the restoration phase objective function 910 * depends on the current barrier parameter. 911 */ 912 bool in_restoration_phase(); 913 914 //@} 915 }; 916 917 } // namespace Ipopt 918 919 #endif 920