1 // Copyright (C) 2005, 2006 International Business Machines and others. 2 // All Rights Reserved. 3 // This code is published under the Common Public License. 4 // 5 // $Id: IpWarmStartIterateInitializer.cpp 759 2006-07-07 03:07:08Z andreasw $ 6 // 7 // Authors: Carl Laird, Andreas Waechter IBM 2005-04-01 8 9 #include "IpWarmStartIterateInitializer.hpp" 10 #include "IpDefaultIterateInitializer.hpp" 11 12 // ToDo make independent of DenseVector 13 #include "IpDenseVector.hpp" 14 15 #ifdef HAVE_CMATH 16 # include <cmath> 17 #else 18 # ifdef HAVE_MATH_H 19 # include <math.h> 20 # else 21 # error "don't have header file for math" 22 # endif 23 #endif 24 25 namespace SimTKIpopt 26 { 27 #ifdef IP_DEBUG 28 static const Index dbg_verbosity = 0; 29 #endif 30 WarmStartIterateInitializer()31 WarmStartIterateInitializer::WarmStartIterateInitializer() 32 : 33 IterateInitializer() 34 {} 35 RegisterOptions(SmartPtr<RegisteredOptions> roptions)36 void WarmStartIterateInitializer::RegisterOptions(SmartPtr<RegisteredOptions> roptions) 37 { 38 roptions->AddLowerBoundedNumberOption( 39 "warm_start_bound_push", 40 "same as bound_push for the regular initializer.", 41 0.0, true, 1e-3); 42 roptions->AddBoundedNumberOption( 43 "warm_start_bound_frac", 44 "same as bound_frac for the regular initializer.", 45 0.0, true, 0.5, false, 1e-3); 46 roptions->AddLowerBoundedNumberOption( 47 "warm_start_mult_bound_push", 48 "same as mult_bound_push for the regular initializer.", 49 0.0, true, 1e-3); 50 roptions->AddNumberOption( 51 "warm_start_mult_init_max", 52 "Maximum initial value for the equality multipliers.", 53 1e6); 54 roptions->AddNumberOption( 55 "warm_start_target_mu", 56 "Unsupported!", 57 0e-3); 58 roptions->AddStringOption2( 59 "warm_start_entire_iterate", 60 "Tells algorithm whether to use the GetWarmStartIterate method in the NLP.", 61 "no", 62 "no", "call GetStartingPoint in the NLP", 63 "yes", "call GetWarmStartIterate in the NLP", 64 ""); 65 } 66 InitializeImpl(const OptionsList & options,const std::string & prefix)67 bool WarmStartIterateInitializer::InitializeImpl(const OptionsList& options, 68 const std::string& prefix) 69 { 70 options.GetNumericValue("warm_start_bound_push", 71 warm_start_bound_push_, prefix); 72 options.GetNumericValue("warm_start_bound_frac", 73 warm_start_bound_frac_, prefix); 74 options.GetNumericValue("warm_start_mult_bound_push", 75 warm_start_mult_bound_push_, prefix); 76 options.GetNumericValue("warm_start_mult_init_max", 77 warm_start_mult_init_max_, prefix); 78 options.GetNumericValue("warm_start_target_mu", 79 warm_start_target_mu_, prefix); 80 options.GetBoolValue("warm_start_entire_iterate", 81 warm_start_entire_iterate_, prefix); 82 83 return true; 84 } 85 SetInitialIterates()86 bool WarmStartIterateInitializer::SetInitialIterates() 87 { 88 DBG_START_METH("WarmStartIterateInitializer::SetInitialIterates", 89 dbg_verbosity); 90 91 // Get the starting values provided by the NLP and store them 92 // in the ip_data current fields. 93 94 SmartPtr<IteratesVector> init_vec; 95 bool have_iterate = false; 96 97 if (warm_start_entire_iterate_) { 98 IpData().InitializeDataStructures(IpNLP(), false, false, false, 99 false, false); 100 101 init_vec = IpData().curr()->MakeNewIteratesVector(true); 102 103 have_iterate = IpNLP().GetWarmStartIterate(*init_vec); 104 105 if (!have_iterate) { 106 Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, 107 "Tried to obtain entire warm start iterate from NLP, but it returned false.\n"); 108 IpData().Append_info_string("NW"); 109 } 110 111 // Make sure given bounds are respected 112 if (have_iterate && warm_start_mult_init_max_>0.) { 113 SmartPtr<Vector> y_c = init_vec->create_new_y_c_copy(); 114 SmartPtr<Vector> tmp = y_c->MakeNew(); 115 tmp->Set(warm_start_mult_init_max_); 116 y_c->ElementWiseMin(*tmp); 117 tmp->Set(-warm_start_mult_init_max_); 118 y_c->ElementWiseMax(*tmp); 119 120 SmartPtr<Vector> y_d = init_vec->create_new_y_d_copy(); 121 tmp = y_d->MakeNew(); 122 tmp->Set(warm_start_mult_init_max_); 123 y_d->ElementWiseMin(*tmp); 124 tmp->Set(-warm_start_mult_init_max_); 125 y_d->ElementWiseMax(*tmp); 126 127 SmartPtr<Vector> z_L = init_vec->create_new_z_L_copy(); 128 tmp = z_L->MakeNew(); 129 tmp->Set(warm_start_mult_init_max_); 130 z_L->ElementWiseMin(*tmp); 131 132 SmartPtr<Vector> z_U = init_vec->create_new_z_U_copy(); 133 tmp = z_U->MakeNew(); 134 tmp->Set(warm_start_mult_init_max_); 135 z_U->ElementWiseMin(*tmp); 136 137 SmartPtr<Vector> v_L = init_vec->create_new_v_L_copy(); 138 tmp = v_L->MakeNew(); 139 tmp->Set(warm_start_mult_init_max_); 140 v_L->ElementWiseMin(*tmp); 141 142 SmartPtr<Vector> v_U = init_vec->create_new_v_U_copy(); 143 tmp = v_U->MakeNew(); 144 tmp->Set(warm_start_mult_init_max_); 145 v_U->ElementWiseMin(*tmp); 146 } 147 } 148 149 if (!have_iterate) { 150 151 ///////////////////////////////////////////////////////////////////// 152 // Initialize primal variables // 153 ///////////////////////////////////////////////////////////////////// 154 155 // Get the initial values for x, y_c, y_d, z_L, z_U, 156 IpData().InitializeDataStructures(IpNLP(), true, true, true, true, true); 157 158 IpData().curr()->x()->Print(Jnlst(), J_VECTOR, J_INITIALIZATION, 159 "user-provided x"); 160 IpData().curr()->y_c()->Print(Jnlst(), J_VECTOR, J_INITIALIZATION, 161 "user-provided y_c"); 162 IpData().curr()->y_d()->Print(Jnlst(), J_VECTOR, J_INITIALIZATION, 163 "user-provided y_d"); 164 IpData().curr()->z_L()->Print(Jnlst(), J_VECTOR, J_INITIALIZATION, 165 "user-provided z_L"); 166 IpData().curr()->z_U()->Print(Jnlst(), J_VECTOR, J_INITIALIZATION, 167 "user-provided z_U"); 168 if (Jnlst().ProduceOutput(J_MOREVECTOR, J_INITIALIZATION)) { 169 IpCq().curr_d()->Print(Jnlst(), J_MOREVECTOR, J_INITIALIZATION, 170 "d at user-provided x"); 171 } 172 173 SmartPtr<Vector> tmp; 174 175 init_vec = IpData().curr()->MakeNewContainer(); 176 177 // If requested, make sure that the multipliers are not too large 178 if (warm_start_mult_init_max_>0.) { 179 SmartPtr<Vector> y_c = init_vec->create_new_y_c_copy(); 180 tmp = y_c->MakeNew(); 181 tmp->Set(warm_start_mult_init_max_); 182 y_c->ElementWiseMin(*tmp); 183 tmp->Set(-warm_start_mult_init_max_); 184 y_c->ElementWiseMax(*tmp); 185 186 SmartPtr<Vector> y_d = init_vec->create_new_y_d_copy(); 187 tmp = y_d->MakeNew(); 188 tmp->Set(warm_start_mult_init_max_); 189 y_d->ElementWiseMin(*tmp); 190 tmp->Set(-warm_start_mult_init_max_); 191 y_d->ElementWiseMax(*tmp); 192 193 SmartPtr<Vector> z_L = init_vec->create_new_z_L_copy(); 194 tmp = z_L->MakeNew(); 195 tmp->Set(warm_start_mult_init_max_); 196 z_L->ElementWiseMin(*tmp); 197 198 SmartPtr<Vector> z_U = init_vec->create_new_z_U_copy(); 199 tmp = z_U->MakeNew(); 200 tmp->Set(warm_start_mult_init_max_); 201 z_U->ElementWiseMin(*tmp); 202 } 203 204 // Get the initial values for v_L and v_U out of y_d 205 SmartPtr<Vector> v_L = init_vec->create_new_v_L(); 206 IpNLP().Pd_L()->TransMultVector(-1., *init_vec->y_d(), 0., *v_L); 207 tmp = v_L->MakeNew(); 208 tmp->Set(warm_start_mult_bound_push_); 209 v_L->ElementWiseMax(*tmp); 210 211 SmartPtr<Vector> v_U = init_vec->create_new_v_U(); 212 IpNLP().Pd_U()->TransMultVector(1., *init_vec->y_d(), 0., *v_U); 213 tmp = v_U->MakeNew(); 214 tmp->Set(warm_start_mult_bound_push_); 215 v_U->ElementWiseMax(*tmp); 216 217 // Initialize slack variables 218 init_vec->Set_s(*IpCq().curr_d()); 219 } 220 221 // Make the corrected values current (and initialize s) 222 IpData().set_trial(init_vec); 223 IpData().AcceptTrialPoint(); 224 225 // Now apply the target mu heuristic if required 226 if (warm_start_target_mu_>0.) { 227 SmartPtr<const Vector> new_x; 228 SmartPtr<const Vector> new_z_L; 229 230 SmartPtr<const IteratesVector> curr = IpData().curr(); 231 process_target_mu(1., *curr->x(), *IpCq().curr_slack_x_L(), 232 *curr->z_L(), *IpNLP().Px_L(), 233 new_x, new_z_L); 234 SmartPtr<const Vector> new_s; 235 SmartPtr<const Vector> new_v_L; 236 process_target_mu(1., *curr->s(), *IpCq().curr_slack_s_L(), 237 *curr->v_L(), *IpNLP().Pd_L(), 238 new_s, new_v_L); 239 240 // Set the trial pointers to new_x and new_s. The process_target_mu 241 // methods below create new vectors in new_x and new_s and do not alter 242 // the existing ones. 243 init_vec->Set_x(*new_x); 244 init_vec->Set_s(*new_s); 245 IpData().set_trial(init_vec); 246 247 SmartPtr<const Vector> new_z_U; 248 process_target_mu(-1., *IpData().trial()->x(), *IpCq().trial_slack_x_U(), 249 *IpData().curr()->z_U(), *IpNLP().Px_U(), 250 new_x, new_z_U); 251 SmartPtr<const Vector> new_v_U; 252 process_target_mu(-1., *IpData().trial()->s(), *IpCq().trial_slack_s_U(), 253 *IpData().curr()->v_U(), *IpNLP().Pd_U(), 254 new_s, new_v_U); 255 256 // Now submit the full modified point 257 init_vec->Set_x(*new_x); 258 init_vec->Set_s(*new_s); 259 // y_c and y_d currently contain a copy of curr()->y_c... 260 // we set them back to the actual pointer to reuse the tags 261 init_vec->Set_y_c(*IpData().curr()->y_c()); 262 init_vec->Set_y_d(*IpData().curr()->y_d()); 263 init_vec->Set_z_L(*new_z_L); 264 init_vec->Set_z_U(*new_z_U); 265 init_vec->Set_v_L(*new_v_L); 266 init_vec->Set_v_U(*new_v_U); 267 IpData().set_trial(init_vec); 268 IpData().AcceptTrialPoint(); 269 270 // We need to call this to make sure that we don't get an error 271 // message because at some point a slack became too small 272 IpCq().ResetAdjustedTrialSlacks(); 273 } 274 275 SmartPtr<const Vector> new_x; 276 SmartPtr<const Vector> new_s; 277 // Push the primal x variables 278 DefaultIterateInitializer::push_variables(Jnlst(), 279 warm_start_bound_push_, 280 warm_start_bound_frac_, 281 "x", 282 *IpData().curr()->x(), 283 new_x, 284 *IpNLP().x_L(), 285 *IpNLP().x_U(), 286 *IpNLP().Px_L(), 287 *IpNLP().Px_U()); 288 289 // Push the primal s variables 290 DefaultIterateInitializer::push_variables(Jnlst(), 291 warm_start_bound_push_, 292 warm_start_bound_frac_, 293 "s", 294 *IpData().curr()->s(), 295 new_s, 296 *IpNLP().d_L(), 297 *IpNLP().d_U(), 298 *IpNLP().Pd_L(), 299 *IpNLP().Pd_U()); 300 301 // Push the multipliers 302 SmartPtr<Vector> new_z_L = IpData().curr()->z_L()->MakeNewCopy(); 303 SmartPtr<Vector> tmp = IpData().curr()->z_L()->MakeNew(); 304 tmp->Set(warm_start_mult_bound_push_); 305 new_z_L->ElementWiseMax(*tmp); 306 307 SmartPtr<Vector> new_z_U = IpData().curr()->z_U()->MakeNewCopy(); 308 tmp = IpData().curr()->z_U()->MakeNew(); 309 tmp->Set(warm_start_mult_bound_push_); 310 new_z_U->ElementWiseMax(*tmp); 311 312 SmartPtr<Vector> new_v_L = IpData().curr()->v_L()->MakeNewCopy(); 313 tmp = IpData().curr()->v_L()->MakeNew(); 314 tmp->Set(warm_start_mult_bound_push_); 315 new_v_L->ElementWiseMax(*tmp); 316 317 SmartPtr<Vector> new_v_U = IpData().curr()->v_U()->MakeNewCopy(); 318 tmp = IpData().curr()->v_U()->MakeNew(); 319 tmp->Set(warm_start_mult_bound_push_); 320 new_v_U->ElementWiseMax(*tmp); 321 322 // Make sure the new variables are current 323 init_vec = IpData().curr()->MakeNewContainer(); 324 init_vec->Set_x(*new_x); 325 init_vec->Set_s(*new_s); 326 init_vec->Set_z_L(*new_z_L); 327 init_vec->Set_z_U(*new_z_U); 328 init_vec->Set_v_L(*new_v_L); 329 init_vec->Set_v_U(*new_v_U); 330 IpData().set_trial(init_vec); 331 IpData().AcceptTrialPoint(); 332 333 IpData().curr()->x()->Print(Jnlst(), J_VECTOR, J_INITIALIZATION, 334 "initial x"); 335 IpData().curr()->s()->Print(Jnlst(), J_VECTOR, J_INITIALIZATION, 336 "initial s"); 337 IpData().curr()->y_c()->Print(Jnlst(), J_VECTOR, J_INITIALIZATION, 338 "initial y_c"); 339 IpData().curr()->y_d()->Print(Jnlst(), J_VECTOR, J_INITIALIZATION, 340 "initial y_d"); 341 IpData().curr()->z_L()->Print(Jnlst(), J_VECTOR, J_INITIALIZATION, 342 "initial z_L"); 343 IpData().curr()->z_U()->Print(Jnlst(), J_VECTOR, J_INITIALIZATION, 344 "initial z_U"); 345 IpData().curr()->v_L()->Print(Jnlst(), J_VECTOR, J_INITIALIZATION, 346 "initial v_L"); 347 IpData().curr()->v_U()->Print(Jnlst(), J_VECTOR, J_INITIALIZATION, 348 "initial v_U"); 349 if (Jnlst().ProduceOutput(J_MOREVECTOR, J_INITIALIZATION)) { 350 IpCq().curr_slack_x_L()->Print(Jnlst(), J_MOREVECTOR, J_INITIALIZATION, 351 "initial slack_x_L"); 352 IpCq().curr_slack_x_U()->Print(Jnlst(), J_MOREVECTOR, J_INITIALIZATION, 353 "initial slack_x_U"); 354 IpCq().curr_slack_s_L()->Print(Jnlst(), J_MOREVECTOR, J_INITIALIZATION, 355 "initial slack_s_L"); 356 IpCq().curr_slack_s_U()->Print(Jnlst(), J_MOREVECTOR, J_INITIALIZATION, 357 "initial slack_s_U"); 358 } 359 360 return true; 361 } 362 process_target_mu(Number factor,const Vector & curr_vars,const Vector & curr_slacks,const Vector & curr_mults,const Matrix & P,SmartPtr<const Vector> & ret_vars,SmartPtr<const Vector> & ret_mults)363 void WarmStartIterateInitializer::process_target_mu(Number factor, 364 const Vector& curr_vars, 365 const Vector& curr_slacks, 366 const Vector& curr_mults, 367 const Matrix& P, 368 SmartPtr<const Vector>& ret_vars, 369 SmartPtr<const Vector>& ret_mults) 370 { 371 SmartPtr<Vector> new_slacks = curr_slacks.MakeNewCopy(); 372 SmartPtr<Vector> new_mults = curr_mults.MakeNewCopy(); 373 adapt_to_target_mu(*new_slacks, *new_mults, warm_start_target_mu_); 374 new_slacks->Axpy(-1, curr_slacks); // this is now correction step 375 SmartPtr<Vector> new_vars = curr_vars.MakeNew(); 376 new_vars->Copy(curr_vars); 377 P.MultVector(factor, *new_slacks, 1., *new_vars); 378 379 ret_vars = ConstPtr(new_vars); 380 ret_mults = ConstPtr(new_mults); 381 } 382 adapt_to_target_mu(Vector & new_s,Vector & new_z,Number target_mu)383 void WarmStartIterateInitializer::adapt_to_target_mu(Vector& new_s, 384 Vector& new_z, 385 Number target_mu) 386 { 387 DBG_ASSERT(new_s.Dim() == new_z.Dim()); 388 389 DenseVector* dnew_s = dynamic_cast<DenseVector*>(&new_s); 390 assert(dnew_s); 391 DenseVector* dnew_z = dynamic_cast<DenseVector*>(&new_z); 392 assert(dnew_z); 393 Number* values_s = dnew_s->Values(); 394 Number* values_z = dnew_z->Values(); 395 396 for (Index i=0; i<new_s.Dim(); i++) { 397 if (values_s[i] > 1e4*values_z[i]) { 398 values_z[i] = target_mu/values_s[i]; 399 if (values_z[i]>values_s[i]) { 400 values_s[i] = values_z[i] = sqrt(target_mu); 401 } 402 } 403 else if (values_z[i] > 1e4*values_s[i]) { 404 values_s[i] = target_mu/values_z[i]; 405 if (values_s[i]>values_z[i]) { 406 values_s[i] = values_z[i] = sqrt(target_mu); 407 } 408 } 409 else { 410 values_s[i] = values_z[i] = sqrt(target_mu); 411 } 412 } 413 } 414 415 } // namespace Ipopt 416