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