1 // Copyright (C) 2005, 2007 International Business Machines and others. 2 // All Rights Reserved. 3 // This code is published under the Eclipse Public License. 4 // 5 // $Id$ 6 // 7 // Authors: Carl Laird, Andreas Waechter IBM 2005-08-04 8 9 #include "IpCGPerturbationHandler.hpp" 10 #include "IpCGPenaltyData.hpp" 11 #include "IpCGPenaltyCq.hpp" 12 13 #ifdef HAVE_CMATH 14 # include <cmath> 15 #else 16 # ifdef HAVE_MATH_H 17 # include <math.h> 18 # else 19 # error "don't have header file for math" 20 # endif 21 #endif 22 23 namespace Ipopt 24 { 25 #if COIN_IPOPT_VERBOSITY > 0 26 static const Index dbg_verbosity = 0; 27 #endif 28 CGPerturbationHandler()29 CGPerturbationHandler::CGPerturbationHandler() 30 : 31 PDPerturbationHandler(), 32 reset_last_(false), 33 degen_iters_max_(3) 34 {} 35 36 void RegisterOptions(SmartPtr<RegisteredOptions> roptions)37 CGPerturbationHandler::RegisterOptions(SmartPtr<RegisteredOptions> roptions) 38 {} 39 InitializeImpl(const OptionsList & options,const std::string & prefix)40 bool CGPerturbationHandler::InitializeImpl(const OptionsList& options, 41 const std::string& prefix) 42 { 43 options.GetNumericValue("max_hessian_perturbation", delta_xs_max_, prefix); 44 options.GetNumericValue("min_hessian_perturbation", delta_xs_min_, prefix); 45 options.GetNumericValue("perturb_inc_fact_first", delta_xs_first_inc_fact_, prefix); 46 options.GetNumericValue("perturb_inc_fact", delta_xs_inc_fact_, prefix); 47 options.GetNumericValue("perturb_dec_fact", delta_xs_dec_fact_, prefix); 48 options.GetNumericValue("first_hessian_perturbation", delta_xs_init_, prefix); 49 options.GetNumericValue("jacobian_regularization_value", delta_cd_val_, prefix); 50 options.GetNumericValue("jacobian_regularization_exponent", delta_cd_exp_, prefix); 51 options.GetBoolValue("perturb_always_cd", perturb_always_cd_, prefix); 52 // The following option has been registered from CGSearchDirCalc 53 options.GetNumericValue("penalty_max", penalty_max_, prefix); 54 // The following option has been registered from CGPenaltyLSAccepter 55 options.GetNumericValue("mult_diverg_feasibility_tol", mult_diverg_feasibility_tol_, prefix); 56 hess_degenerate_ = NOT_YET_DETERMINED; 57 if (!perturb_always_cd_) { 58 jac_degenerate_ = NOT_YET_DETERMINED; 59 } 60 else { 61 jac_degenerate_ = NOT_DEGENERATE; 62 } 63 degen_iters_ = 0; 64 65 delta_x_curr_ = 0.; 66 delta_s_curr_ = 0.; 67 delta_c_curr_ = 0.; 68 delta_d_curr_ = 0.; 69 delta_x_last_ = 0.; 70 delta_s_last_ = 0.; 71 delta_c_last_ = 0.; 72 delta_d_last_ = 0.; 73 74 test_status_ = NO_TEST; 75 76 return PDPerturbationHandler::InitializeImpl(options, prefix); 77 } 78 79 bool ConsiderNewSystem(Number & delta_x,Number & delta_s,Number & delta_c,Number & delta_d)80 CGPerturbationHandler::ConsiderNewSystem(Number& delta_x, Number& delta_s, 81 Number& delta_c, Number& delta_d) 82 { 83 DBG_START_METH("CGPerturbationHandler::ConsiderNewSystem",dbg_verbosity); 84 85 // Check if we can conclude that some components of the system are 86 // structurally degenerate 87 finalize_test(); 88 // If the current iterate is restored from a previous iteration, 89 // initialize perturbationhandler data 90 if (CGPenData().restor_iter() == IpData().iter_count()) { 91 hess_degenerate_ = NOT_YET_DETERMINED; 92 jac_degenerate_ = NOT_YET_DETERMINED; 93 degen_iters_ = 0; 94 hess_degenerate_ = NOT_DEGENERATE; 95 jac_degenerate_ = NOT_DEGENERATE; 96 delta_x_curr_ = 0.; 97 delta_s_curr_ = 0.; 98 delta_c_curr_ = 0.; 99 delta_d_curr_ = 0.; 100 delta_x_last_ = 0.; 101 delta_s_last_ = 0.; 102 delta_c_last_ = 0.; 103 delta_d_last_ = 0.; 104 test_status_ = NO_TEST; 105 } 106 107 // Store the perturbation from the previous matrix 108 if (reset_last_) { 109 delta_x_last_ = delta_x_curr_; 110 delta_s_last_ = delta_s_curr_; 111 delta_c_last_ = delta_c_curr_; 112 delta_d_last_ = delta_d_curr_; 113 } 114 else { 115 if (delta_x_curr_ > 0.) { 116 delta_x_last_ = delta_x_curr_; 117 } 118 if (delta_s_curr_ > 0.) { 119 delta_s_last_ = delta_s_curr_; 120 } 121 if (delta_c_curr_ > 0.) { 122 delta_c_last_ = delta_c_curr_; 123 } 124 if (delta_d_curr_ > 0.) { 125 delta_d_last_ = delta_d_curr_; 126 } 127 } 128 129 DBG_ASSERT((hess_degenerate_ != NOT_YET_DETERMINED || 130 jac_degenerate_ != DEGENERATE) && 131 (jac_degenerate_ != NOT_YET_DETERMINED || 132 hess_degenerate_ != DEGENERATE)); 133 134 if (hess_degenerate_ == NOT_YET_DETERMINED || 135 jac_degenerate_ == NOT_YET_DETERMINED) { 136 if (!perturb_always_cd_ 137 || CGPenCq().curr_cg_pert_fact() < delta_cd() 138 || !CGPenData().NeverTryPureNewton()) { 139 test_status_ = TEST_DELTA_C_EQ_0_DELTA_X_EQ_0; 140 } 141 else { 142 test_status_ = TEST_DELTA_C_GT_0_DELTA_X_EQ_0; 143 } 144 } 145 else { 146 test_status_ = NO_TEST; 147 } 148 149 Number pert_fact = CGPenCq().curr_cg_pert_fact(); 150 if (jac_degenerate_ == DEGENERATE 151 || CGPenData().NeverTryPureNewton() 152 || perturb_always_cd_) { 153 Number mach_eps = std::numeric_limits<Number>::epsilon(); 154 if (pert_fact < 100.*mach_eps 155 && jac_degenerate_ == DEGENERATE) { 156 delta_c = delta_c_curr_ = 100.*mach_eps; 157 } 158 else { 159 delta_c = delta_c_curr_ = pert_fact; 160 } 161 } 162 else { 163 delta_c = delta_c_curr_ = 0.; 164 } 165 CGPenData().SetCurrPenaltyPert(delta_c); 166 167 delta_d = delta_d_curr_ = delta_c; 168 169 if (hess_degenerate_ == DEGENERATE) { 170 delta_x_curr_ = 0.; 171 delta_s_curr_ = 0.; 172 bool retval = get_deltas_for_wrong_inertia(delta_x, delta_s, 173 delta_c, delta_d); 174 if (!retval) { 175 return false; 176 } 177 } 178 else { 179 delta_x = 0.; 180 delta_s = delta_x; 181 } 182 183 delta_x_curr_ = delta_x; 184 delta_s_curr_ = delta_s; 185 delta_c_curr_ = delta_c; 186 delta_d_curr_ = delta_d; 187 188 IpData().Set_info_regu_x(delta_x); 189 190 get_deltas_for_wrong_inertia_called_ = false; 191 192 return true; 193 } 194 195 bool PerturbForSingularity(Number & delta_x,Number & delta_s,Number & delta_c,Number & delta_d)196 CGPerturbationHandler::PerturbForSingularity( 197 Number& delta_x, Number& delta_s, 198 Number& delta_c, Number& delta_d) 199 { 200 DBG_START_METH("CGPerturbationHandler::PerturbForSingularity", 201 dbg_verbosity); 202 203 bool retval; 204 205 // Check for structural degeneracy 206 if (hess_degenerate_ == NOT_YET_DETERMINED || 207 jac_degenerate_ == NOT_YET_DETERMINED) { 208 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 209 "Degeneracy test for hess_degenerate_ = %d and jac_degenerate_ = %d\n test_status_ = %d\n", 210 hess_degenerate_, jac_degenerate_, test_status_); 211 switch (test_status_) { 212 case TEST_DELTA_C_EQ_0_DELTA_X_EQ_0: 213 DBG_ASSERT(delta_x_curr_ == 0. && delta_c_curr_ == 0.); 214 // in this case we haven't tried anything for this matrix yet 215 if (jac_degenerate_ == NOT_YET_DETERMINED) { 216 delta_d_curr_ = delta_c_curr_ = delta_cd(); 217 test_status_ = TEST_DELTA_C_GT_0_DELTA_X_EQ_0; 218 } 219 else { 220 DBG_ASSERT(hess_degenerate_ == NOT_YET_DETERMINED); 221 retval = get_deltas_for_wrong_inertia(delta_x, delta_s, 222 delta_c, delta_d); 223 if (!retval) { 224 return false; 225 } 226 DBG_ASSERT(delta_c == 0. && delta_d == 0.); 227 test_status_ = TEST_DELTA_C_EQ_0_DELTA_X_GT_0; 228 } 229 break; 230 case TEST_DELTA_C_GT_0_DELTA_X_EQ_0: 231 DBG_ASSERT(delta_x_curr_ == 0. && delta_c_curr_ > 0.); 232 DBG_ASSERT(jac_degenerate_ == NOT_YET_DETERMINED); 233 //if (!perturb_always_cd_) { 234 delta_d_curr_ = delta_c_curr_ = 235 Max(delta_cd(), CGPenCq().curr_cg_pert_fact()); 236 //delta_d_curr_ = delta_c_curr_ = 237 // Max(delta_cd(), CGPenCq().curr_cg_pert_fact()); 238 if (delta_d_curr_ < delta_cd()) { 239 test_status_ = TEST_DELTA_C_EQ_0_DELTA_X_GT_0; 240 } 241 else { 242 test_status_ = TEST_DELTA_C_GT_0_DELTA_X_GT_0; 243 } 244 retval = get_deltas_for_wrong_inertia(delta_x, delta_s, 245 delta_c, delta_d); 246 if (!retval) { 247 return false; 248 } 249 /* DBG_ASSERT(delta_c == 0. && delta_d == 0.); */ 250 test_status_ = TEST_DELTA_C_EQ_0_DELTA_X_GT_0; 251 //} 252 /* 253 else { 254 retval = get_deltas_for_wrong_inertia(delta_x, delta_s, 255 delta_c, delta_d); 256 if (!retval) { 257 return false; 258 } 259 DBG_ASSERT(delta_c > 0. && delta_d > 0.); 260 test_status_ = TEST_DELTA_C_GT_0_DELTA_X_GT_0; 261 }*/ 262 break; 263 case TEST_DELTA_C_EQ_0_DELTA_X_GT_0: 264 DBG_ASSERT(delta_x_curr_ > 0. && delta_c_curr_ == 0.); 265 delta_d_curr_ = delta_c_curr_ = Max(delta_cd(), CGPenCq().curr_cg_pert_fact()); 266 //delta_d_curr_ = delta_c_curr_ = CGPenCq().curr_cg_pert_fact(); 267 retval = get_deltas_for_wrong_inertia(delta_x, delta_s, 268 delta_c, delta_d); 269 if (!retval) { 270 return false; 271 } 272 test_status_ = TEST_DELTA_C_GT_0_DELTA_X_GT_0; 273 break; 274 case TEST_DELTA_C_GT_0_DELTA_X_GT_0: 275 retval = get_deltas_for_wrong_inertia(delta_x, delta_s, 276 delta_c, delta_d); 277 if (!retval) { 278 return false; 279 } 280 break; 281 case NO_TEST: 282 DBG_ASSERT(false && "we should not get here."); 283 } 284 } 285 else { 286 if (delta_c_curr_ > 0. || get_deltas_for_wrong_inertia_called_) { 287 // If we already used a perturbation for the constraints, we do 288 // the same thing as if we were encountering negative curvature 289 retval = get_deltas_for_wrong_inertia(delta_x, delta_s, 290 delta_c, delta_d); 291 if (!retval) { 292 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 293 "Can't get_deltas_for_wrong_inertia for delta_x_curr_ = %e and delta_c_curr_ = %e\n", 294 delta_x_curr_, delta_c_curr_); 295 return false; 296 } 297 } 298 else { 299 // Otherwise we now perturb the lower right corner 300 delta_d_curr_ = delta_c_curr_ = delta_cd(); 301 302 // ToDo - also perturb Hessian? 303 IpData().Append_info_string("L"); 304 Number curr_inf = IpCq().curr_primal_infeasibility(NORM_2); 305 if (!CGPenData().NeverTryPureNewton() 306 && curr_inf > mult_diverg_feasibility_tol_) { 307 Number penalty = CGPenCq().compute_curr_cg_penalty_scale(); 308 penalty = Min(penalty_max_, Max(penalty, 309 CGPenData().curr_kkt_penalty())); 310 CGPenData().Set_kkt_penalty(penalty); 311 Number mach_pro = std::numeric_limits<Number>::epsilon(); 312 delta_d_curr_ = delta_c_curr_ = 313 Max(1e3*mach_pro,Max(CGPenCq().curr_cg_pert_fact(),delta_cd())); 314 IpData().Append_info_string("u"); 315 } 316 } 317 } 318 319 delta_x = delta_x_curr_; 320 delta_s = delta_s_curr_; 321 delta_c = delta_c_curr_; 322 delta_d = delta_d_curr_; 323 324 IpData().Set_info_regu_x(delta_x); 325 326 return true; 327 } 328 329 bool get_deltas_for_wrong_inertia(Number & delta_x,Number & delta_s,Number & delta_c,Number & delta_d)330 CGPerturbationHandler::get_deltas_for_wrong_inertia( 331 Number& delta_x, Number& delta_s, 332 Number& delta_c, Number& delta_d) 333 { 334 if (delta_x_curr_ == 0.) { 335 if (delta_x_last_ == 0.) { 336 delta_x_curr_ = delta_xs_init_; 337 } 338 else { 339 delta_x_curr_ = Max(delta_xs_min_, 340 delta_x_last_*delta_xs_dec_fact_); 341 } 342 } 343 else { 344 if (delta_x_last_ == 0. || 1e5*delta_x_last_<delta_x_curr_) { 345 delta_x_curr_ = delta_xs_first_inc_fact_*delta_x_curr_; 346 } 347 else { 348 delta_x_curr_ = delta_xs_inc_fact_*delta_x_curr_; 349 } 350 } 351 if (delta_x_curr_ > delta_xs_max_) { 352 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 353 "delta_x perturbation is becoming too large: %e\n", 354 delta_x_curr_); 355 delta_x_last_ = 0.; 356 delta_s_last_ = 0.; 357 IpData().Append_info_string("dx"); 358 return false; 359 } 360 361 delta_s_curr_ = delta_x_curr_; 362 363 delta_x = delta_x_curr_; 364 delta_s = delta_s_curr_; 365 delta_c = delta_c_curr_; 366 delta_d = delta_d_curr_; 367 368 IpData().Set_info_regu_x(delta_x); 369 370 get_deltas_for_wrong_inertia_called_ = true; 371 372 return true; 373 } 374 375 bool PerturbForWrongInertia(Number & delta_x,Number & delta_s,Number & delta_c,Number & delta_d)376 CGPerturbationHandler::PerturbForWrongInertia( 377 Number& delta_x, Number& delta_s, 378 Number& delta_c, Number& delta_d) 379 { 380 DBG_START_METH("CGPerturbationHandler::PerturbForWrongInertia", 381 dbg_verbosity); 382 383 // Check if we can conclude that components of the system are 384 // structurally degenerate (we only get here if the most recent 385 // perturbation for a test did not result in a singular system) 386 finalize_test(); 387 388 bool retval = get_deltas_for_wrong_inertia(delta_x, delta_s, delta_c, delta_d); 389 if (!retval && delta_c==0.) { 390 DBG_ASSERT(delta_d == 0.); 391 delta_c_curr_ = delta_cd(); 392 delta_d_curr_ = delta_c_curr_; 393 delta_x_curr_ = 0.; 394 delta_s_curr_ = 0.; 395 test_status_ = NO_TEST; 396 if (hess_degenerate_ == DEGENERATE) { 397 hess_degenerate_ = NOT_YET_DETERMINED; 398 } 399 retval = get_deltas_for_wrong_inertia(delta_x, delta_s, delta_c, delta_d); 400 } 401 return retval; 402 } 403 404 void CurrentPerturbation(Number & delta_x,Number & delta_s,Number & delta_c,Number & delta_d)405 CGPerturbationHandler::CurrentPerturbation( 406 Number& delta_x, Number& delta_s, 407 Number& delta_c, Number& delta_d) 408 { 409 delta_x = delta_x_curr_; 410 delta_s = delta_s_curr_; 411 delta_c = delta_c_curr_; 412 delta_d = delta_d_curr_; 413 } 414 415 Number delta_cd()416 CGPerturbationHandler::delta_cd() 417 { 418 return delta_cd_val_ * pow(IpData().curr_mu(), delta_cd_exp_); 419 } 420 421 void finalize_test()422 CGPerturbationHandler::finalize_test() 423 { 424 switch (test_status_) { 425 case NO_TEST: 426 return; 427 case TEST_DELTA_C_EQ_0_DELTA_X_EQ_0: 428 if (hess_degenerate_ == NOT_YET_DETERMINED && 429 jac_degenerate_ == NOT_YET_DETERMINED) { 430 hess_degenerate_ = NOT_DEGENERATE; 431 jac_degenerate_ = NOT_DEGENERATE; 432 IpData().Append_info_string("Nhj "); 433 } 434 else if (hess_degenerate_ == NOT_YET_DETERMINED) { 435 hess_degenerate_ = NOT_DEGENERATE; 436 IpData().Append_info_string("Nh "); 437 } 438 else if (jac_degenerate_ == NOT_YET_DETERMINED) { 439 jac_degenerate_ = NOT_DEGENERATE; 440 IpData().Append_info_string("Nj "); 441 } 442 break; 443 case TEST_DELTA_C_GT_0_DELTA_X_EQ_0: 444 if (hess_degenerate_ == NOT_YET_DETERMINED) { 445 hess_degenerate_ = NOT_DEGENERATE; 446 IpData().Append_info_string("Nh "); 447 } 448 if (jac_degenerate_ == NOT_YET_DETERMINED) { 449 degen_iters_++; 450 if (degen_iters_ >= degen_iters_max_) { 451 jac_degenerate_ = DEGENERATE; 452 IpData().Append_info_string("Dj "); 453 } 454 IpData().Append_info_string("L"); 455 } 456 break; 457 case TEST_DELTA_C_EQ_0_DELTA_X_GT_0: 458 if (jac_degenerate_ == NOT_YET_DETERMINED) { 459 jac_degenerate_ = NOT_DEGENERATE; 460 IpData().Append_info_string("Nj "); 461 } 462 if (hess_degenerate_ == NOT_YET_DETERMINED) { 463 degen_iters_++; 464 if (degen_iters_ >= degen_iters_max_) { 465 hess_degenerate_ = DEGENERATE; 466 IpData().Append_info_string("Dh "); 467 } 468 } 469 break; 470 case TEST_DELTA_C_GT_0_DELTA_X_GT_0: 471 degen_iters_++; 472 if (degen_iters_ >= degen_iters_max_) { 473 hess_degenerate_ = DEGENERATE; 474 jac_degenerate_ = DEGENERATE; 475 IpData().Append_info_string("Dhj "); 476 } 477 IpData().Append_info_string("L"); 478 break; 479 } 480 } 481 482 } // namespace Ipopt 483