ConstrainedOptPack_QPSolverRelaxedQPOPTSOL.cpp

00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include <assert.h>
00030 
00031 #include <vector>
00032 
00033 #include "ConstrainedOptPack_QPSolverRelaxedQPOPTSOL.hpp"
00034 #include "AbstractLinAlgPack_SpVectorOp.hpp"
00035 #include "AbstractLinAlgPack_MatrixSymOp.hpp"
00036 #include "AbstractLinAlgPack_EtaVector.hpp"
00037 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
00038 #include "AbstractLinAlgPack_SortByDescendingAbsValue.hpp"
00039 #include "AbstractLinAlgPack_sparse_bounds.hpp"
00040 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
00041 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00042 #include "ProfileHackPack_profile_hack.hpp"
00043 
00044 // /////////////////////////////////////////////////////////////////
00045 //
00046 // This subclass uses a relaxation of the equality and inequality
00047 // constraints.  The mapping to the arguments of QPOPT or QPSOL
00048 // is done as follows.
00049 //
00050 //  QP formulation:
00051 //  ---------------
00052 //
00053 //  min          g'*d + 1/2*d'*G*d + (eta + 1/2*eta^2)*M
00054 //  d <: R^n
00055 //         
00056 //  s.t.
00057 //               etaL <=  eta
00058 //               dL   <=  d                       <= dU
00059 //               eL   <=  op(E)*d - b*eta         <= eU
00060 //                        op(F)*d + (1 - eta) * f  = 0
00061 //
00062 //  Rearranged to :
00063 //  ---------------
00064 //
00065 //  min          [ g', M ] * [  d  ] + 1/2 * [ d', eta ] * [ G  0 ] * [  d  ]
00066 //                           [ eta ]                       [ 0  M ]   [ eta ]
00067 //
00068 //  s.t.         [  bL  ]    [   I  ,  0 ]              [ dU  ]
00069 //               [ etaL ] <= [   0  ,  1 ] * [  d  ] <= [ inf ]
00070 //               [  eL  ]    [ op(E), -b ]   [ eta ]    [ eU  ]
00071 //               [  -f  ]    [ op(F), -f ]              [ -f  ]
00072 //
00073 //  Which maps to the QPSOL interface which is:
00074 //  -------------------------------------------
00075 //
00076 //  min           CVEC' * X + 1/2 * X'* H * X
00077 //
00078 //  s.t.          BL <= [    X   ] <= BU
00079 //                      [  A * X ]
00080 //
00081 //  Which gives us:
00082 //
00083 //  X    = [ d; eta ] 
00084 //  CVEC = [ g; M ]
00085 //  H    = [ G, 0; 0, M ]
00086 //  BL   = [ bL, etaL, eL, -f ]
00087 //  BU   = [ bU, inf,  eU, -f ]
00088 //  A    = [ op(E), -b; op(F), -f ]
00089 //
00090 
00091 namespace LinAlgOpPack {
00092   using AbstractLinAlgPack::Vp_StV;
00093   using AbstractLinAlgPack::Mp_StM;
00094   using AbstractLinAlgPack::Vp_StMtV;
00095 }
00096 
00097 // ///////////////////////////////////////
00098 // Members for QPSolverRelaxedQPOPTSOL
00099 
00100 namespace ConstrainedOptPack {
00101 
00102 QPSolverRelaxedQPOPTSOL::QPSolverRelaxedQPOPTSOL()
00103   :N_(0)
00104   ,bigM_(1e+10)
00105   ,use_as_bigM_(1e+10)
00106   ,G_(NULL)
00107 {}
00108 
00109 QPSolverRelaxedQPOPTSOL::~QPSolverRelaxedQPOPTSOL()
00110 {
00111   this->release_memory();
00112 }
00113 
00114 const MatrixOp* QPSolverRelaxedQPOPTSOL::G() const
00115 {
00116   return G_;
00117 }
00118 
00119 value_type QPSolverRelaxedQPOPTSOL::use_as_bigM() const
00120 {
00121   return use_as_bigM_;
00122 }
00123 
00124 // Overridden from QPSolverRelaxed
00125 
00126 QPSolverStats
00127 QPSolverRelaxedQPOPTSOL::get_qp_stats() const
00128 {
00129   return qp_stats_;
00130 }
00131 
00132 void QPSolverRelaxedQPOPTSOL::release_memory()
00133 {
00134   // Todo: resize to zero all the matrices and vectors
00135 }
00136 
00137 QPSolverStats::ESolutionType
00138 QPSolverRelaxedQPOPTSOL::imp_solve_qp(
00139     std::ostream* out, EOutputLevel olevel, ERunTests test_what
00140     ,const Vector& g, const MatrixSymOp& G
00141     ,value_type etaL
00142     ,const Vector* dL, const Vector* dU
00143     ,const MatrixOp* E, BLAS_Cpp::Transp trans_E, const Vector* b
00144     ,const Vector* eL, const Vector* eU
00145     ,const MatrixOp* F, BLAS_Cpp::Transp trans_F, const Vector* f
00146     ,value_type* obj_d
00147     ,value_type* eta, VectorMutable* d
00148     ,VectorMutable* nu
00149     ,VectorMutable* mu, VectorMutable* Ed
00150     ,VectorMutable* lambda, VectorMutable* Fd
00151   )
00152 {
00153 
00154   using AbstractLinAlgPack::VectorDenseEncap;
00155   using AbstractLinAlgPack::VectorDenseMutableEncap;
00156 
00157 #ifdef PROFILE_HACK_ENABLED
00158   ProfileHackPack::ProfileTiming profile_timing( "QPSolverRelaxedQPOPTSOL::imp_solve_qp(...)" );
00159 #endif
00160 
00161   const size_type n = d->dim();
00162   const value_type inf = this->infinite_bound();
00163   
00164   //
00165   // Map to the input arguments for QPOPT or QPSOL
00166   //
00167 
00168   // N
00169   N_ = n + 1; // With relaxation
00170 
00171   // NCLIN
00172   n_inequ_bnds_ = ( E ? AbstractLinAlgPack::num_bounded(*eL,*eU,inf) : 0 );
00173   NCLIN_ = n_inequ_bnds_ + (F ? f->dim() : 0);
00174 
00175   // A, BL, BU
00176   A_.resize(NCLIN_,N_);
00177   BL_.resize(N_+NCLIN_);
00178   BU_.resize(N_+NCLIN_);
00179   if(dL) {
00180     VectorDenseEncap dL_de(*dL);
00181     BL_(1,n) = dL_de();
00182   }
00183   else {
00184     BL_(1,n) = -inf;
00185   }
00186   if(dU) {
00187     VectorDenseEncap dU_de(*dU);
00188     BU_(1,n) = dU_de();
00189   }
00190   else {
00191     BU_(1,n) = -inf;
00192   }
00193   BL_(N_) = etaL;
00194   BU_(N_) = +inf;
00195   TEST_FOR_EXCEPTION(
00196     E!=NULL, std::logic_error
00197     ,"Error, the QPOPT/QPSOL wrapper has not been updated for general inequalities yet!"
00198     );
00199 /* ToDo: Update this when needed!
00200   if( E ) {
00201     i_inequ_bnds_.resize(n_inequ_bnds_);
00202     if( n_inequ_bnds_ < b->dim() ) {
00203       // Initialize BL, BU, and A for sparse bounds on general inequalities
00204       //
00205       // read iterators
00206       AbstractLinAlgPack::sparse_bounds_itr
00207         eLU_itr( eL->begin(), eL->end(), eL->offset()
00208              , eU->begin(), eU->end(), eU->offset(), inf );
00209       // written iterators
00210       DVector::iterator
00211         BL_itr    = BL_.begin() + N_,
00212         BU_itr    = BU_.begin() + N_;
00213       ibnds_t::iterator
00214         ibnds_itr = i_inequ_bnds_.begin();
00215       // loop
00216       for(size_type i = 1; i <= n_inequ_bnds_; ++i, ++eLU_itr, ++ibnds_itr ) {
00217         TEST_FOR_EXCEPT( !( !eLU_itr.at_end() ) );
00218         const size_type k      = eLU_itr.indice();
00219         *BL_itr++              = eLU_itr.lbound();
00220         *BU_itr++              = eLU_itr.ubound();
00221         *ibnds_itr             = k;  // Only for my record
00222         // Add the corresponding row of [ op(E), -b ] to A
00223         // y == A.row(i)
00224         // y(1,n) = op(E')*e_k
00225         DVectorSlice y = A_.row(i);
00226         AbstractLinAlgPack::EtaVector e_k(k,eL->dim());
00227         LinAlgOpPack::V_MtV( &y(1,n), *E, BLAS_Cpp::trans_not(trans_E), e_k() ); // op(E')*e_k
00228         // y(n+1) = -b(k)
00229         y(n+1) = -(*b)(k);
00230       }
00231     }
00232     else {
00233       // Initialize BL, BU and A for dense bounds on general inequalities
00234       //
00235       // Initialize BL(N+1:N+n_inequ_bnds), BU(N+1:N+n_inequ_bnds)
00236       // and i_inequ_bnds_ = identity (only for my record, not used by QPKWIK)
00237       AbstractLinAlgPack::sparse_bounds_itr
00238         eLU_itr( eL->begin(), eL->end(), eL->offset()
00239              , eU->begin(), eU->end(), eU->offset(), inf );
00240       DVector::iterator
00241         BL_itr    = BL_.begin() + N_,
00242         BU_itr    = BU_.begin() + N_;
00243       ibnds_t::iterator
00244         ibnds_itr = i_inequ_bnds_.begin();
00245       for(size_type i = 1; i <= n_inequ_bnds_; ++i, ++eLU_itr, ++ibnds_itr ) {
00246         TEST_FOR_EXCEPT( !( !eLU_itr.at_end() ) );
00247         const size_type k      = eLU_itr.indice();
00248         *BL_itr++              = eLU_itr.lbound();
00249         *BU_itr++              = eLU_itr.ubound();
00250         *ibnds_itr             = k;  // Only for my record
00251       }
00252       // A(1:n_inequ_bnds,1:n) = op(E)
00253       LinAlgOpPack::assign( &A_(1,n_inequ_bnds_,1,n), *E, trans_E );
00254       // A(1:n_inequ_bnds,n+1) = -b
00255       LinAlgOpPack::V_StV( &A_.col(n+1)(1,n_inequ_bnds_), -1.0, *b );
00256     }
00257   }
00258 */
00259   TEST_FOR_EXCEPTION(
00260     F!=NULL, std::logic_error
00261     ,"Error, the QPOPT/QPSOL wrapper has not been updated for general equalities yet!"
00262     );
00263 /* ToDo: Update this when needed!
00264   if( F ) {
00265     // BL(N+n_inequ_bnds+1:N+NCLIN) = -f
00266     LinAlgOpPack::V_StV( &BL_(N_+n_inequ_bnds_+1,N_+NCLIN_), -1.0, *f );
00267     // BU(N+n_inequ_bnds+1:N+NCLIN) = -f
00268     LinAlgOpPack::V_StV( &BU_(N_+n_inequ_bnds_+1,N_+NCLIN_), -1.0, *f );
00269     // A(n_inequ_bnds+1:NCLIN,1:n) = op(F)
00270     LinAlgOpPack::assign( &A_(n_inequ_bnds_+1,NCLIN_,1,n), *F, trans_F );
00271     // A(n_inequ_bnds+1:NCLIN,n+1) = -f
00272     LinAlgOpPack::V_StV( &A_.col(n+1)(n_inequ_bnds_+1,NCLIN_), -1.0, *f );
00273   }
00274 */
00275   
00276   // CVEC
00277   CVEC_.resize(N_);
00278   CVEC_(1,n) = VectorDenseEncap(g)();
00279   CVEC_(n+1) = bigM_;
00280 
00281   // HESS
00282   G_ = &G; // That's all there is to it!
00283 
00284   // ISTATE
00285   ISTATE_.resize(N_+NCLIN_);
00286   std::fill( ISTATE_.begin(), ISTATE_.end(), 0 ); // cold start!
00287   ISTATE_[n] = 1; // Make eta >= etaL active
00288 
00289   // X
00290   X_.resize(N_);
00291   X_(1,n) = VectorDenseEncap(*d)();
00292   X_(n+1) = *eta;
00293 
00294   // AX
00295   // will be resized by QPOPT but not QPSOL
00296 
00297   // CLAMBDA
00298   CLAMDA_.resize(N_+NCLIN_);
00299 
00300   // LIWORK, IWORK
00301   LIWORK_ = liwork(N_,NCLIN_);
00302   if(static_cast<f_int>(IWORK_.size()) < LIWORK_) IWORK_.resize(LIWORK_);
00303 
00304   // LWORK, WORK
00305   LWORK_ = lrwork(N_,NCLIN_);
00306   if(static_cast<f_int>(WORK_.size()) < LWORK_) WORK_.resize(LWORK_);
00307 
00308   // We need to initialize some warm start information if
00309   // it was given by the user!
00310   bool warm_start = false;
00311   if( (nu && nu->nz()) || (mu && mu->nz() ) ) {
00312     // Let's try a warm start
00313     if(nu) {
00314       VectorDenseEncap nu_de(*nu);
00315       for(int i = 1; i <= n; ++i ) {
00316         if( nu_de()(i) < 0.0 )
00317           ISTATE_[i-1] = 1; // Lower bound is active
00318         else if( nu_de()(i) > 0.0 )
00319           ISTATE_[i-1] = 2; // Upper bound is active
00320       }
00321     }
00322     TEST_FOR_EXCEPTION(
00323       mu!=NULL, std::logic_error
00324       ,"Error, the QPOPT/QPSOL wrapper has not been updated for general inequalities yet!"
00325       );
00326 /* ToDo: Update below when needed!
00327     if(mu) {
00328       const SpVectorSlice::difference_type o = mu->offset();
00329       for( SpVectorSlice::const_iterator itr = mu->begin(); itr != mu->end(); ++itr ) {
00330         if( itr->value() < 0.0 )
00331           ISTATE_[ itr->indice() + o + n ] = 1; // Lower bound is active
00332         else if( itr->value() > 0.0 )
00333           ISTATE_[ itr->indice() + o + n ] = 2; // Upper bound is active
00334       }
00335     }
00336 */
00337     warm_start = true;
00338   }
00339 
00340   //
00341   // Solve the QP using QPOPT or QPSOL
00342   //
00343 
00344   const EInform inform_return = call_qp_solver(warm_start);
00345 
00346   //
00347   // Map from the output from QPOPT or QPSOL
00348   //
00349 
00350   // d
00351   {
00352     VectorDenseMutableEncap d_de(*d);
00353     d_de() = X_(1,n);
00354   }
00355   
00356   // eta
00357   *eta = X_(n+1);
00358 
00359   // obj_d
00360   if(obj_d)
00361     *obj_d = OBJ_ - (*eta) * bigM_ - 0.5 * (*eta)*(*eta) * use_as_bigM_; 
00362 
00363   // nu
00364   if(nu) {
00365     VectorDenseMutableEncap nu_de(*nu);
00366     nu_de() = 0.0;
00367     ISTATE_t::const_iterator
00368       istate_itr = ISTATE_.begin();
00369     DVector::const_iterator
00370       clamda_itr = CLAMDA_.begin();
00371     for( size_type i = 1; i <= n; ++i, ++istate_itr, ++clamda_itr ) {
00372       const f_int state = *istate_itr;
00373       switch(state) {
00374         case -2: // The lower bound is violated by more than feas_tol
00375         case -1: // The upper bound is violated by more than feas_tol
00376           // What do we do?
00377           break;
00378         case 0: // Within bounds by more than feas_tol
00379           break;
00380         case 1: // lower bound is active
00381         case 2: // upper bound is active
00382         case 3: // the bounds are equal and are satisfied
00383           nu_de()(i) = -(*clamda_itr); // Different sign convention
00384           break;
00385         case 4: // Temporaraly fixed at current value
00386           // What do we do?
00387           break;
00388         default:
00389           TEST_FOR_EXCEPT(true); // Should not get here!
00390       }
00391     }
00392   }
00393   
00394   // mu
00395   TEST_FOR_EXCEPTION(
00396     n_inequ_bnds_!=0, std::logic_error
00397     ,"Error, the QPOPT/QPSOL wrapper has not been updated for general inequalities yet!"
00398     );
00399 /* ToDo: Update below when needed!
00400   if( n_inequ_bnds_ ) {
00401     mu->resize(b->dim(),n_inequ_bnds_);
00402     typedef SpVector::element_type ele_t;
00403     ISTATE_t::const_iterator
00404       istate_itr = ISTATE_.begin() + N_;
00405     DVector::const_iterator
00406       clamda_itr = CLAMDA_.begin() + N_;
00407     ibnds_t::const_iterator
00408       bnd_itr = i_inequ_bnds_.begin();
00409     for( size_type k = 1; k <= n_inequ_bnds_; ++k, ++istate_itr, ++clamda_itr, ++bnd_itr )
00410     {
00411       const f_int state = *istate_itr;
00412       const size_type j = *bnd_itr;
00413       switch(state) {
00414           case -2: // The lower bound is violated by more than feas_tol
00415           case -1: // The upper bound is violated by more than feas_tol
00416           // What do we do?
00417           break;
00418           case 0: // Within bounds by more than feas_tol
00419           break;
00420           case 1: // lower bound is active
00421           case 2: // upper bound is active
00422           case 3: // the bounds are equal and are satisfied
00423           mu->add_element(ele_t(j,-(*clamda_itr))); // Different sign!
00424           break;
00425           case 4: // Temporaraly fixed at current value
00426           // What do we do?
00427           break;
00428           default:
00429           TEST_FOR_EXCEPT(true); // Should not get here!
00430       }
00431     }
00432     mu->assume_sorted(true);
00433   }
00434   else if(E) {
00435     mu->resize( eL->dim(), 0 );
00436   }
00437 */
00438 
00439   TEST_FOR_EXCEPTION(
00440     F!=NULL, std::logic_error
00441     ,"Error, the QPOPT/QPSOL wrapper has not been updated for general equalities yet!"
00442     );
00443 /* ToDo: Update this when needed!
00444 
00445   // lambda
00446   if( F ) {
00447     LinAlgOpPack::V_StV( lambda, -1.0, CLAMDA_(N_+n_inequ_bnds_+1,N_+NCLIN_) );
00448     // Validate istate
00449     ISTATE_t::const_iterator
00450       istate_itr = ISTATE_.begin() + N_ + n_inequ_bnds_;
00451     for( size_type k = 1; k <= f->dim(); ++k, ++istate_itr ) {
00452       TEST_FOR_EXCEPT( !(  *istate_itr == 3  ) );
00453     }
00454   }
00455 
00456   // Ed, Fd
00457   if( E && AX_.size() && eL->dim() == n_inequ_bnds_ ) {
00458     if( Ed ) { // Ed = AX + b*eta
00459       *Ed = AX_(1,n_inequ_bnds_);
00460       if( *eta > 0.0 )
00461         LinAlgOpPack::Vp_StV( Ed, *eta, *b );
00462     }
00463     if( Fd ) { // Fd = AX + f*eta
00464       *Fd = AX_(n_inequ_bnds_+1,NCLIN_);
00465       if( *eta > 0.0 )
00466         LinAlgOpPack::Vp_StV( Fd, *eta, *f );
00467     }
00468   }
00469   else {
00470     if(Ed)
00471       LinAlgOpPack::V_MtV( Ed, *E, trans_E, *d );
00472     if(Fd)
00473       LinAlgOpPack::V_MtV( Fd, *F, trans_F, *d );
00474   }
00475 
00476 */
00477 
00478   //
00479   // Setup the QP statistics
00480   //
00481 
00482   QPSolverStats::ESolutionType solution_type  = QPSolverStats::SOLUTION_TYPE_NOT_KNOWN;
00483   QPSolverStats::EConvexity    convexity_type = QPSolverStats::CONVEXITY_NOT_KNOWN;
00484   switch(inform_return) {
00485       case STRONG_LOCAL_MIN:
00486       solution_type  = QPSolverStats::OPTIMAL_SOLUTION;
00487       convexity_type =  QPSolverStats::CONVEX;
00488       break;
00489       case WEAK_LOCAL_MIN:
00490       solution_type  = QPSolverStats::OPTIMAL_SOLUTION;
00491       convexity_type =  QPSolverStats::NONCONVEX;
00492       break;
00493       case MAX_ITER_EXCEEDED:
00494       solution_type  = QPSolverStats::PRIMAL_FEASIBLE_POINT;
00495       convexity_type =  QPSolverStats::CONVEXITY_NOT_KNOWN;
00496       break;
00497       case OTHER_ERROR:
00498       solution_type  = QPSolverStats::SUBOPTIMAL_POINT;
00499       convexity_type =  QPSolverStats::CONVEXITY_NOT_KNOWN;
00500       break;
00501   }
00502   qp_stats_.set_stats(
00503     solution_type, convexity_type
00504     ,ITER_, QPSolverStats::NOT_KNOWN, QPSolverStats::NOT_KNOWN
00505     ,warm_start, *eta > 0.0 );
00506 
00507   return qp_stats_.solution_type();
00508 }
00509 
00510 } // end namespace ConstrainedOptPack

Generated on Wed May 12 21:51:09 2010 for ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization by  doxygen 1.4.7