1 // Copyright (C) 2008, 2010 International Business Machines and others. 2 // All Rights Reserved. 3 // This code is published under the Eclipse Public License. 4 // 5 // $Id: IpTNLP.hpp 1235 2008-05-22 14:38:40Z andreasw $ 6 // 7 // Authors: Andreas Waechter IBM 2008-08-25 8 9 #include "IpNLPBoundsRemover.hpp" 10 #include "IpCompoundVector.hpp" 11 #include "IpCompoundMatrix.hpp" 12 #include "IpCompoundSymMatrix.hpp" 13 #include "IpIdentityMatrix.hpp" 14 #include "IpTransposeMatrix.hpp" 15 #include "IpDenseVector.hpp" 16 #include "IpZeroMatrix.hpp" 17 18 namespace Ipopt 19 { 20 #if COIN_IPOPT_VERBOSITY > 0 21 static const Index dbg_verbosity = 0; 22 #endif 23 NLPBoundsRemover(NLP & nlp,bool allow_twosided_inequalities)24 NLPBoundsRemover::NLPBoundsRemover(NLP& nlp, 25 bool allow_twosided_inequalities /* = false */) 26 : 27 nlp_(&nlp), 28 allow_twosided_inequalities_(allow_twosided_inequalities) 29 {} 30 31 bool GetSpaces(SmartPtr<const VectorSpace> & x_space,SmartPtr<const VectorSpace> & c_space,SmartPtr<const VectorSpace> & d_space,SmartPtr<const VectorSpace> & x_l_space,SmartPtr<const MatrixSpace> & px_l_space,SmartPtr<const VectorSpace> & x_u_space,SmartPtr<const MatrixSpace> & px_u_space,SmartPtr<const VectorSpace> & d_l_space,SmartPtr<const MatrixSpace> & pd_l_space,SmartPtr<const VectorSpace> & d_u_space,SmartPtr<const MatrixSpace> & pd_u_space,SmartPtr<const MatrixSpace> & Jac_c_space,SmartPtr<const MatrixSpace> & Jac_d_space,SmartPtr<const SymMatrixSpace> & Hess_lagrangian_space)32 NLPBoundsRemover::GetSpaces(SmartPtr<const VectorSpace>& x_space, 33 SmartPtr<const VectorSpace>& c_space, 34 SmartPtr<const VectorSpace>& d_space, 35 SmartPtr<const VectorSpace>& x_l_space, 36 SmartPtr<const MatrixSpace>& px_l_space, 37 SmartPtr<const VectorSpace>& x_u_space, 38 SmartPtr<const MatrixSpace>& px_u_space, 39 SmartPtr<const VectorSpace>& d_l_space, 40 SmartPtr<const MatrixSpace>& pd_l_space, 41 SmartPtr<const VectorSpace>& d_u_space, 42 SmartPtr<const MatrixSpace>& pd_u_space, 43 SmartPtr<const MatrixSpace>& Jac_c_space, 44 SmartPtr<const MatrixSpace>& Jac_d_space, 45 SmartPtr<const SymMatrixSpace>& Hess_lagrangian_space) 46 { 47 DBG_START_METH("NLPBoundsRemover::GetSpaces", dbg_verbosity); 48 SmartPtr<const VectorSpace> d_space_orig; 49 SmartPtr<const VectorSpace> x_l_space_orig; 50 SmartPtr<const MatrixSpace> px_l_space_orig; 51 SmartPtr<const VectorSpace> x_u_space_orig; 52 SmartPtr<const MatrixSpace> px_u_space_orig; 53 SmartPtr<const VectorSpace> d_l_space_orig; 54 SmartPtr<const MatrixSpace> pd_l_space_orig; 55 SmartPtr<const VectorSpace> d_u_space_orig; 56 SmartPtr<const MatrixSpace> pd_u_space_orig; 57 SmartPtr<const MatrixSpace> Jac_d_space_orig; 58 59 bool retval = nlp_->GetSpaces(x_space, c_space, d_space_orig, 60 x_l_space_orig, px_l_space_orig, 61 x_u_space_orig, px_u_space_orig, 62 d_l_space_orig, pd_l_space_orig, 63 d_u_space_orig, pd_u_space_orig, 64 Jac_c_space, Jac_d_space_orig, 65 Hess_lagrangian_space); 66 if (!retval) { 67 return retval; 68 } 69 // Keep a copy of the expansion matrices for the x bounds 70 Px_l_orig_ = px_l_space_orig->MakeNew(); 71 Px_u_orig_ = px_u_space_orig->MakeNew(); 72 73 // create the new d_space 74 Index total_dim = d_space_orig->Dim() + x_l_space_orig->Dim() + 75 x_u_space_orig->Dim(); 76 SmartPtr<CompoundVectorSpace> d_space_new = 77 new CompoundVectorSpace(3, total_dim); 78 d_space_new->SetCompSpace(0, *d_space_orig); 79 d_space_new->SetCompSpace(1, *x_l_space_orig); 80 d_space_new->SetCompSpace(2, *x_u_space_orig); 81 d_space = GetRawPtr(d_space_new); 82 83 // create the new (emply) x_l and x_u spaces, and also the 84 // corresponding projection matrix spaces 85 x_l_space = new DenseVectorSpace(0); 86 x_u_space = new DenseVectorSpace(0); 87 px_l_space = new ZeroMatrixSpace(x_space->Dim(),0); 88 px_u_space = new ZeroMatrixSpace(x_space->Dim(),0); 89 90 // create the new d_l and d_u vector spaces 91 total_dim = d_l_space_orig->Dim() + x_l_space_orig->Dim(); 92 SmartPtr<CompoundVectorSpace> d_l_space_new = 93 new CompoundVectorSpace(2, total_dim); 94 d_l_space_new->SetCompSpace(0, *d_l_space_orig); 95 d_l_space_new->SetCompSpace(1, *x_l_space_orig); 96 d_l_space = GetRawPtr(d_l_space_new); 97 total_dim = d_u_space_orig->Dim() + x_u_space_orig->Dim(); 98 SmartPtr<CompoundVectorSpace> d_u_space_new = 99 new CompoundVectorSpace(2, total_dim); 100 d_u_space_new->SetCompSpace(0, *d_u_space_orig); 101 d_u_space_new->SetCompSpace(1, *x_u_space_orig); 102 d_u_space = GetRawPtr(d_u_space_new); 103 104 // create the new d_l and d_u projection matrix spaces 105 Index total_rows = d_space_orig->Dim() + x_l_space_orig->Dim() + 106 x_u_space_orig->Dim(); 107 Index total_cols = d_l_space_orig->Dim() + x_l_space_orig->Dim(); 108 SmartPtr<CompoundMatrixSpace> pd_l_space_new = 109 new CompoundMatrixSpace(3, 2, total_rows, total_cols); 110 pd_l_space_new->SetBlockRows(0, d_space_orig->Dim()); 111 pd_l_space_new->SetBlockRows(1, x_l_space_orig->Dim()); 112 pd_l_space_new->SetBlockRows(2, x_u_space_orig->Dim()); 113 pd_l_space_new->SetBlockCols(0, d_l_space_orig->Dim()); 114 pd_l_space_new->SetBlockCols(1, x_l_space_orig->Dim()); 115 pd_l_space_new->SetCompSpace(0, 0, *pd_l_space_orig, true); 116 SmartPtr<const MatrixSpace> identity_space = 117 new IdentityMatrixSpace(x_l_space_orig->Dim()); 118 pd_l_space_new->SetCompSpace(1, 1, *identity_space, true); 119 pd_l_space = GetRawPtr(pd_l_space_new); 120 121 total_cols = d_u_space_orig->Dim() + x_u_space_orig->Dim(); 122 SmartPtr<CompoundMatrixSpace> pd_u_space_new = 123 new CompoundMatrixSpace(3, 2, total_rows, total_cols); 124 pd_u_space_new->SetBlockRows(0, d_space_orig->Dim()); 125 pd_u_space_new->SetBlockRows(1, x_l_space_orig->Dim()); 126 pd_u_space_new->SetBlockRows(2, x_u_space_orig->Dim()); 127 pd_u_space_new->SetBlockCols(0, d_u_space_orig->Dim()); 128 pd_u_space_new->SetBlockCols(1, x_u_space_orig->Dim()); 129 pd_u_space_new->SetCompSpace(0, 0, *pd_u_space_orig, true); 130 identity_space = new IdentityMatrixSpace(x_u_space_orig->Dim()); 131 pd_u_space_new->SetCompSpace(2, 1, *identity_space, true); 132 pd_u_space = GetRawPtr(pd_u_space_new); 133 134 // Jacobian for inequalities matrix space 135 total_rows = d_space_orig->Dim() + x_l_space_orig->Dim() + 136 x_u_space_orig->Dim(); 137 total_cols = x_space->Dim(); 138 SmartPtr<CompoundMatrixSpace> Jac_d_space_new = 139 new CompoundMatrixSpace(3, 1, total_rows, total_cols); 140 Jac_d_space_new->SetBlockRows(0, d_space_orig->Dim()); 141 Jac_d_space_new->SetBlockRows(1, x_l_space_orig->Dim()); 142 Jac_d_space_new->SetBlockRows(2, x_u_space_orig->Dim()); 143 Jac_d_space_new->SetBlockCols(0, x_space->Dim()); 144 Jac_d_space_new->SetCompSpace(0, 0, *Jac_d_space_orig); 145 SmartPtr<MatrixSpace> trans_px_l_space_orig = 146 new TransposeMatrixSpace(GetRawPtr(px_l_space_orig)); 147 Jac_d_space_new->SetCompSpace(1, 0, *trans_px_l_space_orig, true); 148 SmartPtr<MatrixSpace> trans_px_u_space_orig = 149 new TransposeMatrixSpace(GetRawPtr(px_u_space_orig)); 150 Jac_d_space_new->SetCompSpace(2, 0, *trans_px_u_space_orig, true); 151 Jac_d_space = GetRawPtr(Jac_d_space_new); 152 153 // We keep the original d_space around in order to be able to do 154 // the sanity check later 155 d_space_orig_ = d_space_orig; 156 157 return true; 158 } 159 160 bool GetBoundsInformation(const Matrix & Px_L,Vector & x_L,const Matrix & Px_U,Vector & x_U,const Matrix & Pd_L,Vector & d_L,const Matrix & Pd_U,Vector & d_U)161 NLPBoundsRemover::GetBoundsInformation(const Matrix& Px_L, 162 Vector& x_L, 163 const Matrix& Px_U, 164 Vector& x_U, 165 const Matrix& Pd_L, 166 Vector& d_L, 167 const Matrix& Pd_U, 168 Vector& d_U) 169 { 170 const CompoundMatrix* comp_pd_l = 171 static_cast<const CompoundMatrix*>(&Pd_L); 172 DBG_ASSERT(dynamic_cast<const CompoundMatrix*>(&Pd_L)); 173 SmartPtr<const Matrix> pd_l_orig = comp_pd_l->GetComp(0,0); 174 175 const CompoundMatrix* comp_pd_u = 176 static_cast<const CompoundMatrix*>(&Pd_U); 177 DBG_ASSERT(dynamic_cast<const CompoundMatrix*>(&Pd_U)); 178 SmartPtr<const Matrix> pd_u_orig = comp_pd_u->GetComp(0,0); 179 180 CompoundVector* comp_d_l = static_cast<CompoundVector*>(&d_L); 181 DBG_ASSERT(dynamic_cast<CompoundVector*>(&d_L)); 182 SmartPtr<Vector> d_l_orig = comp_d_l->GetCompNonConst(0); 183 SmartPtr<Vector> x_l_orig = comp_d_l->GetCompNonConst(1); 184 185 CompoundVector* comp_d_u = static_cast<CompoundVector*>(&d_U); 186 DBG_ASSERT(dynamic_cast<CompoundVector*>(&d_U)); 187 SmartPtr<Vector> d_u_orig = comp_d_u->GetCompNonConst(0); 188 SmartPtr<Vector> x_u_orig = comp_d_u->GetCompNonConst(1); 189 190 // Here we do a santiy check to make sure that no inequality 191 // constraint has two non-infite bounds. 192 if (d_space_orig_->Dim()>0 && !allow_twosided_inequalities_) { 193 SmartPtr<Vector> d = d_space_orig_->MakeNew(); 194 SmartPtr<Vector> tmp = d_l_orig->MakeNew(); 195 tmp->Set(1.); 196 pd_l_orig->MultVector(1., *tmp, 0., *d); 197 tmp = d_u_orig->MakeNew(); 198 tmp->Set(1.); 199 pd_u_orig->MultVector(1., *tmp, 1., *d); 200 Number dmax = d->Amax(); 201 ASSERT_EXCEPTION(dmax==1., INVALID_NLP, "In NLPBoundRemover, an inequality with both lower and upper bounds was detected"); 202 Number dmin = d->Min(); 203 ASSERT_EXCEPTION(dmin==1., INVALID_NLP, "In NLPBoundRemover, an inequality with without bounds was detected."); 204 } 205 206 bool retval = 207 nlp_->GetBoundsInformation(*Px_l_orig_, *x_l_orig, *Px_u_orig_, 208 *x_u_orig, *pd_l_orig, *d_l_orig, 209 *pd_u_orig, *d_u_orig); 210 return retval; 211 } 212 213 bool GetStartingPoint(SmartPtr<Vector> x,bool need_x,SmartPtr<Vector> y_c,bool need_y_c,SmartPtr<Vector> y_d,bool need_y_d,SmartPtr<Vector> z_L,bool need_z_L,SmartPtr<Vector> z_U,bool need_z_U)214 NLPBoundsRemover::GetStartingPoint(SmartPtr<Vector> x, 215 bool need_x, 216 SmartPtr<Vector> y_c, 217 bool need_y_c, 218 SmartPtr<Vector> y_d, 219 bool need_y_d, 220 SmartPtr<Vector> z_L, 221 bool need_z_L, 222 SmartPtr<Vector> z_U, 223 bool need_z_U) 224 { 225 SmartPtr<Vector> y_d_orig; 226 SmartPtr<Vector> z_L_orig; 227 SmartPtr<Vector> z_U_orig; 228 if (need_y_d) { 229 CompoundVector* comp_y_d = static_cast<CompoundVector*>(GetRawPtr(y_d)); 230 DBG_ASSERT(dynamic_cast<CompoundVector*>(GetRawPtr(y_d))); 231 y_d_orig = comp_y_d->GetCompNonConst(0); 232 z_L_orig = comp_y_d->GetCompNonConst(1); 233 z_U_orig = comp_y_d->GetCompNonConst(2); 234 } 235 bool retval = 236 nlp_->GetStartingPoint(x, need_x, y_c, need_y_c, y_d_orig, need_y_d, 237 z_L_orig, need_y_d, z_U_orig, need_y_d); 238 return retval; 239 } 240 241 bool Eval_d(const Vector & x,Vector & d)242 NLPBoundsRemover::Eval_d(const Vector& x, Vector& d) 243 { 244 CompoundVector* comp_d = static_cast<CompoundVector*>(&d); 245 DBG_ASSERT(dynamic_cast<CompoundVector*>(&d)); 246 SmartPtr<Vector> d_orig = comp_d->GetCompNonConst(0); 247 248 bool retval = nlp_->Eval_d(x, *d_orig); 249 if (retval) { 250 SmartPtr<Vector> x_L = comp_d->GetCompNonConst(1); 251 SmartPtr<Vector> x_U = comp_d->GetCompNonConst(2); 252 Px_l_orig_->TransMultVector(1., x, 0., *x_L); 253 Px_u_orig_->TransMultVector(1., x, 0., *x_U); 254 } 255 return retval; 256 } 257 258 bool Eval_jac_d(const Vector & x,Matrix & jac_d)259 NLPBoundsRemover::Eval_jac_d(const Vector& x, Matrix& jac_d) 260 { 261 CompoundMatrix* comp_jac_d = static_cast<CompoundMatrix*>(&jac_d); 262 DBG_ASSERT(dynamic_cast<CompoundMatrix*>(&jac_d)); 263 SmartPtr<const MatrixSpace> jac_d_space = comp_jac_d->OwnerSpace(); 264 const CompoundMatrixSpace* comp_jac_d_space = 265 static_cast<const CompoundMatrixSpace*>(GetRawPtr(jac_d_space)); 266 DBG_ASSERT(dynamic_cast<const CompoundMatrixSpace*>(GetRawPtr(jac_d_space))); 267 SmartPtr<Matrix> jac_d_orig = comp_jac_d_space->GetCompSpace(0,0)->MakeNew(); 268 bool retval = nlp_->Eval_jac_d(x, *jac_d_orig); 269 if (retval) { 270 comp_jac_d->SetComp(0, 0, *jac_d_orig); 271 } 272 return retval; 273 } 274 275 bool Eval_h(const Vector & x,Number obj_factor,const Vector & yc,const Vector & yd,SymMatrix & h)276 NLPBoundsRemover::Eval_h(const Vector& x, Number obj_factor, 277 const Vector& yc, const Vector& yd, SymMatrix& h) 278 { 279 const CompoundVector* comp_yd = static_cast<const CompoundVector*>(&yd); 280 DBG_ASSERT(dynamic_cast<const CompoundVector*>(&yd)); 281 SmartPtr<const Vector> yd_orig = comp_yd->GetComp(0); 282 283 bool retval = nlp_->Eval_h(x, obj_factor, yc, *yd_orig, h); 284 return retval; 285 } 286 287 void FinalizeSolution(SolverReturn status,const Vector & x,const Vector & z_L,const Vector & z_U,const Vector & c,const Vector & d,const Vector & y_c,const Vector & y_d,Number obj_value,const IpoptData * ip_data,IpoptCalculatedQuantities * ip_cq)288 NLPBoundsRemover::FinalizeSolution(SolverReturn status, 289 const Vector& x, const Vector& z_L, 290 const Vector& z_U, 291 const Vector& c, const Vector& d, 292 const Vector& y_c, const Vector& y_d, 293 Number obj_value, 294 const IpoptData* ip_data, 295 IpoptCalculatedQuantities* ip_cq) 296 { 297 const CompoundVector* comp_d = static_cast<const CompoundVector*>(&d); 298 DBG_ASSERT(dynamic_cast<const CompoundVector*>(&d)); 299 SmartPtr<const Vector> d_orig = comp_d->GetComp(0); 300 301 const CompoundVector* comp_y_d = static_cast<const CompoundVector*>(&y_d); 302 DBG_ASSERT(dynamic_cast<const CompoundVector*>(&y_d)); 303 SmartPtr<const Vector> y_d_orig = comp_y_d->GetComp(0); 304 SmartPtr<const Vector> z_L_orig = comp_y_d->GetComp(1); 305 SmartPtr<const Vector> z_U_orig = comp_y_d->GetComp(2); 306 307 SmartPtr<Vector> z_L_new = z_L_orig->MakeNewCopy(); 308 z_L_new->Scal(-1.); 309 310 nlp_->FinalizeSolution(status, x, *z_L_new, *z_U_orig, c, *d_orig, 311 y_c, *y_d_orig, obj_value, ip_data, ip_cq); 312 } 313 314 void GetScalingParameters(const SmartPtr<const VectorSpace> x_space,const SmartPtr<const VectorSpace> c_space,const SmartPtr<const VectorSpace> d_space,Number & obj_scaling,SmartPtr<Vector> & x_scaling,SmartPtr<Vector> & c_scaling,SmartPtr<Vector> & d_scaling) const315 NLPBoundsRemover::GetScalingParameters( 316 const SmartPtr<const VectorSpace> x_space, 317 const SmartPtr<const VectorSpace> c_space, 318 const SmartPtr<const VectorSpace> d_space, 319 Number& obj_scaling, 320 SmartPtr<Vector>& x_scaling, 321 SmartPtr<Vector>& c_scaling, 322 SmartPtr<Vector>& d_scaling) const 323 { 324 const CompoundVectorSpace* comp_d_space = 325 static_cast<const CompoundVectorSpace*>(GetRawPtr(d_space)); 326 DBG_ASSERT(dynamic_cast<const CompoundVectorSpace*>(GetRawPtr(d_space))); 327 SmartPtr<const VectorSpace> d_space_orig = comp_d_space->GetCompSpace(0); 328 329 SmartPtr<Vector> d_scaling_orig; 330 nlp_->GetScalingParameters(x_space, c_space, d_space_orig, obj_scaling, 331 x_scaling, c_scaling, d_scaling_orig); 332 333 if (IsValid(x_scaling) || IsValid(d_scaling_orig)) { 334 335 SmartPtr<CompoundVector> comp_d_scaling = 336 comp_d_space->MakeNewCompoundVector(); 337 338 SmartPtr<Vector> xL_scaling = comp_d_scaling->GetCompNonConst(1); 339 SmartPtr<Vector> xU_scaling = comp_d_scaling->GetCompNonConst(2); 340 if (IsValid(x_scaling)) { 341 Px_l_orig_->TransMultVector(1., *x_scaling, 0., *xL_scaling); 342 Px_u_orig_->TransMultVector(1., *x_scaling, 0., *xU_scaling); 343 } 344 else { 345 xL_scaling->Set(1.); 346 xU_scaling->Set(1.); 347 } 348 349 if (IsValid(d_scaling_orig)) { 350 comp_d_scaling->SetComp(0, *d_scaling_orig); 351 } 352 else { 353 comp_d_scaling->GetCompNonConst(0)->Set(1.); 354 } 355 356 d_scaling = GetRawPtr(comp_d_scaling); 357 } 358 else { 359 d_scaling = NULL; 360 } 361 } 362 363 } 364 365