MoochoPack_FeasibilityStepReducedStd_Strategy.cpp

Go to the documentation of this file.
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 "MoochoPack_FeasibilityStepReducedStd_Strategy.hpp"
00030 #include "MoochoPack_NLPAlgo.hpp"
00031 #include "MoochoPack_NLPAlgoState.hpp"
00032 #include "MoochoPack_Exceptions.hpp"
00033 #include "Teuchos_dyn_cast.hpp"
00034 
00035 namespace MoochoPack {
00036 
00037 FeasibilityStepReducedStd_Strategy::FeasibilityStepReducedStd_Strategy(
00038   const quasi_range_space_step_ptr_t   &quasi_range_space_step
00039   ,const qp_solver_ptr_t               &qp_solver
00040   ,const qp_tester_ptr_t               &qp_tester
00041   ,EQPObjective                        qp_objective
00042   ,EQPTesting                          qp_testing
00043   )
00044   :quasi_range_space_step_(quasi_range_space_step)
00045   ,qp_solver_(qp_solver)
00046   ,qp_tester_(qp_tester)
00047   ,qp_objective_(qp_objective)
00048   ,qp_testing_(qp_testing)
00049   ,dl_iq_(dl_name)
00050   ,du_iq_(du_name)
00051   ,current_k_(-1)
00052 {}
00053 
00054 bool FeasibilityStepReducedStd_Strategy::compute_feasibility_step(
00055   std::ostream& out, EJournalOutputLevel olevel, NLPAlgo *algo, NLPAlgoState *s
00056   ,const Vector& xo, const Vector& c_xo, VectorMutable* w
00057     )
00058 {
00059   using Teuchos::dyn_cast;
00060 
00061 /* Todo: UPdate below code!
00062 
00063   // problem dimensions
00064   const size_type
00065     n = algo->nlp().n(),
00066     m = algo->nlp().m(),
00067     r  = s->equ_decomp().size();
00068 
00069   // Compute the quasi-range space step Ywy
00070   Workspace<value_type> Ywy_ws(wss,xo.size());
00071   DVectorSlice                Ywy(&Ywy_ws[0],Ywy_ws.size());
00072   if(!quasi_range_space_step().solve_quasi_range_space_step(
00073     out,olevel,algo,s,xo,c_xo,&Ywy ))
00074     return false;
00075 
00076   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00077     out << "\n||Ywy||2     = " << DenseLinAlgPack::norm_2(Ywy);
00078     out << std::endl;
00079   }
00080   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00081     out << "\nYwy = \n" << Ywy;
00082   }
00083 
00084   //
00085   // Set the bounds on the null space QP subproblem:
00086   //
00087   // d_bounds_k.l <= (xo - x_k) + (1-eta) * Ywy + Z*wz <= d_bounds_k.u
00088   // =>
00089   // bl <= Z*wz - eta * Ywy <= bu
00090   //
00091   // where:   bl = d_bounds_k.l - (xo - x_k) - Ywy
00092   //          bu = d_bounds_k.u - (xo - x_k) - Ywy
00093   //
00094   // Above, fix the variables that are at an active bound as equalities
00095   // to maintain the same active-set.
00096   //
00097   const SparseBounds
00098     &d_bounds = d_bounds_(*s).get_k(0);
00099   const SpVectorSlice
00100     &dl = d_bounds.l,
00101     &du = d_bounds.u;
00102   const DVector
00103     &x_k = s->x().get_k(0).v();
00104   const SpVector
00105     &nu_k = s->nu().get_k(0);
00106   TEST_FOR_EXCEPT( !( nu_k.is_sorted() ) );
00107   SpVector bl(n,n), bu(n,n);
00108   sparse_bounds_itr
00109     d_bnds_itr(dl.begin(),dl.end(),dl.offset(),du.begin(),du.end(),du.offset());
00110   SpVector::const_iterator
00111     nu_itr     = nu_k.begin(),
00112     nu_end     = nu_k.end();
00113   for( ; !d_bnds_itr.at_end(); ++d_bnds_itr ) {
00114     typedef SpVectorSlice::element_type ele_t;
00115     const size_type i = d_bnds_itr.indice();
00116     while( nu_itr != nu_end && nu_itr->indice() + nu_k.offset() < i )
00117       ++nu_itr;
00118     if( nu_itr != nu_end && nu_itr->indice() + nu_k.offset() == i ) {
00119       const value_type
00120         act_bnd = nu_itr->value() > 0.0 ? d_bnds_itr.ubound() : d_bnds_itr.lbound();
00121       bl.add_element(ele_t( i, act_bnd - xo(i) + x_k(i) - Ywy(i) ));
00122       bu.add_element(ele_t( i, act_bnd - xo(i) + x_k(i) - Ywy(i) ));
00123     }
00124     else {
00125       if( d_bnds_itr.lbound() != -d_bnds_itr.big_bnd() )
00126         bl.add_element(ele_t(i,d_bnds_itr.lbound()  - xo(i) + x_k(i) - Ywy(i) ));
00127       if( d_bnds_itr.ubound() != +d_bnds_itr.big_bnd() )
00128         bu.add_element(ele_t(i, d_bnds_itr.ubound() - xo(i) + x_k(i) - Ywy(i) ));
00129     }
00130   }
00131   bl.assume_sorted(true);
00132   bu.assume_sorted(true);
00133   //
00134   // Setup the objective function for the null space QP subproblem
00135   //
00136   // 
00137   // OBJ_MIN_FULL_STEP
00138   //    min 1/2 * (Y*wy + Z*wz)'*(Y*wy + Z*wz) = 1/2*wy'*Y'*Y*wy + (Z'*Y*wy)'*wz + 1/2*wz'*(Z'*Z)*wz
00139   //    => grad = (Z'*Y*wy), Hess = Z'*Z
00140   //
00141   // OBJ_MIN_WZ
00142   //    min 1/2 * wz'*wz => grad = 0, Hess = I
00143   //
00144   // OBJ_RSQP
00145   //    min qp_grad_k'*wz + 1/2 * wz'*rHL_k*wz
00146   //    => grad = qp_grad, Hess = rHL_k
00147   //
00148   const MatrixOp
00149     &Z_k = s->Z().get_k(0);
00150   if( current_k_ != s->k() ) {
00151     if( qp_objective() != OBJ_RSQP )
00152       grad_store_.resize(n-r);
00153     if( qp_objective() == OBJ_MIN_FULL_STEP )
00154       Hess_store_.resize(n-r+1,n-r+1);
00155   }
00156   DVectorSlice grad;
00157   switch(qp_objective())
00158   {
00159       case OBJ_MIN_FULL_STEP: // grad = (Z'*Ywy), Hess = Z'*Z
00160     {
00161       grad.bind( grad_store_() );
00162       if( current_k_ != s->k() ) {
00163         // grad = (Z'*Ywy)
00164         LinAlgOpPack::V_MtV( &grad, Z_k, BLAS_Cpp::trans, Ywy );
00165         // Hess = Z'*Z
00166         DMatrixSliceSym S(Hess_store_(2,n-r+1,1,n-r),BLAS_Cpp::lower); // Must be strictly lower triangular here!
00167         Z_k.syrk( BLAS_Cpp::trans, 1.0, 0.0, &S ); // S = 1.0*Z'*Z + 0.0*S
00168         MatrixSymPosDefCholFactor
00169           *H_ptr = NULL;
00170         if( Hess_ptr_.get() == NULL || dynamic_cast<const MatrixSymPosDefCholFactor*>(Hess_ptr_.get()) == NULL )
00171           Hess_ptr_ = new MatrixSymPosDefCholFactor;
00172         H_ptr = const_cast<MatrixSymPosDefCholFactor*>(dynamic_cast<const MatrixSymPosDefCholFactor*>(Hess_ptr_.get()));
00173         TEST_FOR_EXCEPT( !( H_ptr ) ); // Should not be null!
00174         H_ptr->init_setup(
00175           &Hess_store_()  // The original matrix is stored in the lower triangular part (below diagonal)!
00176           ,NULL           // Nothing to deallocate
00177           ,n-r
00178           ,true           // maintain the original factor
00179           ,false          // don't maintain the factor
00180           ,true           // allow the factor to be computed if needed
00181           ,true           // Set the view
00182           ,1.0            // Scale the matrix by one
00183           );
00184       }
00185       break;
00186     }
00187       case OBJ_MIN_NULL_SPACE_STEP: // grad = 0, Hess = I
00188     {
00189       grad.bind( grad_store_() );
00190       MatrixSymIdent
00191         *H_ptr = NULL;
00192       if( Hess_ptr_.get() == NULL || dynamic_cast<const MatrixSymIdent*>(Hess_ptr_.get()) == NULL )
00193         Hess_ptr_ = new MatrixSymIdent;
00194       if( current_k_ != s->k() ) {
00195         H_ptr = const_cast<MatrixSymIdent*>(dynamic_cast<const MatrixSymIdent*>(Hess_ptr_.get()));
00196         TEST_FOR_EXCEPT( !( H_ptr ) ); // Should not be null!
00197         H_ptr->init_setup(n-r,1.0);
00198         grad = 0.0;
00199       }
00200       break;
00201     }
00202       case OBJ_RSQP: // grad = qp_grad, Hess = rHL_k
00203     {
00204       grad.bind( s->qp_grad().get_k(0)() );
00205       Hess_ptr_ = Hess_ptr_t( &s->rHL().get_k(0), false ); // don't delete memory!
00206       break;
00207     }
00208       defaut:
00209       TEST_FOR_EXCEPT(true); // Not a valid option
00210   }
00211 
00212   //
00213   // Solve the null space subproblem
00214   //
00215 
00216   Workspace<value_type>  wz_ws(wss,n-r),Zwz_ws(wss,n);
00217   DVectorSlice                 wz(&wz_ws[0],wz_ws.size());
00218   DVectorSlice                 Zwz(&Zwz_ws[0],Zwz_ws.size());
00219   value_type                  qp_eta      = 0;
00220 
00221   bool throw_qp_failure = false;
00222 
00223   if( bl.nz() == 0 && bu.nz() == 0 && m-r == 0 ) {
00224     //
00225     // Just solve the unconstrainted QP
00226     //
00227     // wz = - inv(Hess)*grad
00228 #ifdef _WINDOWS
00229     const MatrixFactorized &Hess = dynamic_cast<const MatrixFactorized&>(*Hess_ptr_);
00230 #else
00231     const MatrixFactorized &Hess = dyn_cast<const MatrixFactorized>(*Hess_ptr_);
00232 #endif
00233     AbstractLinAlgPack::V_InvMtV( &wz, Hess, BLAS_Cpp::no_trans, grad );
00234     DenseLinAlgPack::Vt_S(&wz,-1.0);
00235     // Zwz = Z*wz
00236     LinAlgOpPack::V_MtV( &Zwz, Z_k, BLAS_Cpp::no_trans, wz );
00237   }
00238   else {
00239 
00240     //
00241     // Set the arguments to the QP subproblem
00242     //
00243 
00244     DVectorSlice      qp_g    = grad;
00245     const MatrixOp&   qp_G    = *Hess_ptr_;
00246     const value_type    qp_etaL   = 0.0;
00247     SpVectorSlice     qp_dL(NULL,0,0,n-r);  // If nz() == 0 then no simple bounds
00248     SpVectorSlice     qp_dU(NULL,0,0,n-r);
00249     const MatrixOp    *qp_E   = NULL;
00250     BLAS_Cpp::Transp    qp_trans_E  = BLAS_Cpp::no_trans;
00251     DVectorSlice             qp_b;
00252     SpVectorSlice     qp_eL(NULL,0,0,n);
00253     SpVectorSlice     qp_eU(NULL,0,0,n);
00254     const MatrixOp    *qp_F   = NULL;
00255     BLAS_Cpp::Transp    qp_trans_F  = BLAS_Cpp::no_trans;
00256     DVectorSlice        qp_f;
00257     DVectorSlice        qp_d    = wz;
00258     SpVector        *qp_nu    = NULL;
00259     SpVector        *qp_mu    = NULL;
00260     DVectorSlice        qp_Ed;
00261     DVectorSlice        qp_lambda;
00262 
00263     SpVector _nu_wz, _nu_Dwz,   // Possible storage for multiplers for separate inequality
00264       _nu;        // constriants for wz.
00265     DVector _Dwz;       // Possible storage for D*wz computed by QP solver?
00266 
00267     //
00268     // Determine if we can use simple bounds on wz.
00269     // 
00270     // If we have a variable reduction null space matrix
00271     // (with any choice for Y) then:
00272     // 
00273     // w = Z*wz + (1-eta) * Y*wy
00274     // 
00275     // [ w(dep)   ]  = [ D ] * wz  + (1-eta) * [ Ywy(dep)   ]
00276     // [ w(indep) ]    [ I ]                   [ Ywy(indep) ]
00277     // 
00278     // For a cooridinate decomposition (Y = [ I ; 0 ]) then Ywy(indep) = 0 and
00279     // in this case the bounds on d(indep) become simple bounds on pz even
00280     // with the relaxation.
00281     // 
00282     // Otherwise, we can not use simple variable bounds and implement the
00283     // relaxation properly.
00284     // 
00285 
00286     const ZVarReductMatrix
00287       *Zvr = dynamic_cast<const ZVarReductMatrix*>( &Z_k );
00288     Range1D
00289       indep   = Zvr ? Zvr->indep() : Range1D(),
00290       dep   = Zvr ? Zvr->dep()   : Range1D();
00291 
00292     const bool
00293       use_simple_wz_bounds = ( Zvr!=NULL && DenseLinAlgPack::norm_inf(Ywy(indep))==0.0 );
00294 
00295     if( use_simple_wz_bounds ) {
00296 
00297       // Set simple bound constraints on pz
00298       qp_dL.bind( bl(indep) );
00299       qp_dU.bind( bu(indep) );
00300       qp_nu = &( _nu_wz = s->nu().get_k(0)(indep) );  // warm start?
00301     
00302       // Set general inequality constraints for D*pz
00303       qp_E = &Zvr->D();
00304       qp_b.bind( Ywy(dep) );
00305       qp_eL.bind( bl(dep) );
00306       qp_eU.bind( bu(dep) );
00307       qp_mu = &( _nu_Dwz = s->nu().get_k(0)(dep) ); // warm start?
00308       _Dwz.resize(r);
00309       qp_Ed.bind(_Dwz()); // Part of Zwz will be updated directly!
00310 
00311     }
00312     else {
00313 
00314       // Set general inequality constraints for Z*pz
00315       qp_E = &Z_k;
00316       qp_b.bind( Ywy() );
00317       qp_eL.bind( bl() );
00318       qp_eU.bind( bu() );
00319       qp_mu = &(_nu = s->nu().get_k(0));  // warm start??
00320       qp_Ed.bind(Zwz);  // Zwz
00321     }
00322 
00323     // Set the general equality constriants (if they exist)
00324     DVector q(m-r);
00325     Range1D undecomp = s->con_undecomp();
00326     if( m > r ) {
00327       TEST_FOR_EXCEPT(true); // ToDo: Implement when needed!
00328     }
00329 
00330     // Setup the rest of the arguments
00331 
00332     QPSolverRelaxed::EOutputLevel
00333       qp_olevel;
00334     switch( olevel ) {
00335         case PRINT_NOTHING:
00336         qp_olevel = QPSolverRelaxed::PRINT_NONE;
00337         break;
00338         case PRINT_BASIC_ALGORITHM_INFO:
00339         qp_olevel = QPSolverRelaxed::PRINT_BASIC_INFO;
00340         break;
00341         case PRINT_ALGORITHM_STEPS:
00342         qp_olevel = QPSolverRelaxed::PRINT_BASIC_INFO;
00343         break;
00344         case PRINT_ACTIVE_SET:
00345         qp_olevel = QPSolverRelaxed::PRINT_ITER_SUMMARY;
00346         break;
00347         case PRINT_VECTORS:
00348         qp_olevel = QPSolverRelaxed::PRINT_ITER_VECTORS;
00349         break;
00350         case PRINT_ITERATION_QUANTITIES:
00351         qp_olevel = QPSolverRelaxed::PRINT_EVERY_THING;
00352         break;
00353         default:
00354         TEST_FOR_EXCEPT(true);
00355     }
00356 
00357     //
00358     // Solve the QP
00359     // 
00360     const QPSolverStats::ESolutionType
00361       solution_type =
00362       qp_solver().solve_qp(
00363         int(olevel) == int(PRINT_NOTHING) ? NULL : &out
00364         , qp_olevel
00365         , algo->algo_cntr().check_results()
00366         ? QPSolverRelaxed::RUN_TESTS :  QPSolverRelaxed::NO_TESTS
00367         , qp_g, qp_G, qp_etaL, qp_dL, qp_dU
00368         , qp_E, qp_trans_E, qp_E ? &qp_b : NULL
00369         , qp_E ? &qp_eL : NULL, qp_E ? &qp_eU : NULL
00370         , qp_F, qp_trans_F, qp_F ? &qp_f : NULL
00371         , NULL
00372         , &qp_eta, &qp_d
00373         , qp_nu
00374         , qp_mu, qp_E ? &qp_Ed : NULL
00375         , qp_F ? &qp_lambda : NULL, NULL
00376         );
00377 
00378     //
00379     // Check the optimality conditions for the QP
00380     //
00381     std::ostringstream omsg;
00382     if(   qp_testing() == QP_TEST
00383         || ( qp_testing() == QP_TEST_DEFAULT && algo->algo_cntr().check_results() )  )
00384     {
00385       if( int(olevel) >= int(PRINT_ALGORITHM_STEPS) ) {
00386         out << "\nChecking the optimality conditions of the reduced QP subproblem ...\n";
00387       }
00388       if(!qp_tester().check_optimality_conditions(
00389         solution_type
00390         , int(olevel) == int(PRINT_NOTHING) ? NULL : &out
00391         , int(olevel) >= int(PRINT_VECTORS) ? true : false
00392         , int(olevel) >= int(PRINT_ITERATION_QUANTITIES) ? true : false
00393         , qp_g, qp_G, qp_etaL, qp_dL, qp_dU
00394         , qp_E, qp_trans_E, qp_E ? &qp_b : NULL
00395         , qp_E ? &qp_eL : NULL, qp_E ? &qp_eU : NULL
00396         , qp_F, qp_trans_F, qp_F ? &qp_f : NULL
00397         , NULL
00398         , &qp_eta, &qp_d
00399         , qp_nu
00400         , qp_mu, qp_E ? &qp_Ed : NULL
00401         , qp_F ? &qp_lambda : NULL, NULL
00402         ))
00403       {
00404         omsg << "\n*** Alert! at least one of the QP optimality conditions did not check out.\n";
00405         if(  static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00406           out << omsg.str();
00407         }
00408         throw_qp_failure = true;
00409       }
00410     }
00411 
00412     if( solution_type !=  QPSolverStats::OPTIMAL_SOLUTION ) {
00413       if(  static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00414         out << "\nCould not solve the QP!\n";
00415       }
00416       return false;
00417     }
00418 
00419     //
00420     // Set the solution
00421     //
00422     if( use_simple_wz_bounds ) {
00423       // Set Zwz
00424       Zwz(dep)   = _Dwz;
00425       Zwz(indep) = wz;
00426     }
00427     else {
00428       // Everything should already be set!
00429     }
00430 
00431     // Cut back Ywy = (1-eta) * Ywy
00432     const value_type eps = std::numeric_limits<value_type>::epsilon();
00433     if( fabs(qp_eta - 0.0) > eps ) {
00434       if(  static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00435         out
00436           << "\n*** Alert! the QP was infeasible (eta = "<<qp_eta<<").  Cutting back Ywy_k = (1.0 - eta)*Ywy  ...\n";
00437       }
00438       DenseLinAlgPack::Vt_S( &Ywy , 1.0 - qp_eta );
00439     }
00440   }
00441 
00442   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00443     out << "\n||wz||inf    = " << DenseLinAlgPack::norm_inf(wz);
00444     out << "\n||Zwz||2     = " << DenseLinAlgPack::norm_2(Zwz);
00445     if(qp_eta > 0.0) out << "\n||Ypy||2 = " << DenseLinAlgPack::norm_2(Ywy);
00446     out << std::endl;
00447   }
00448   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00449     out << "\nwz = \n" << wz;
00450     out << "\nZwz = \n" << Zwz;
00451     if(qp_eta > 0.0) out << "\nYwy = \n" << Ywy;
00452   }
00453   if( qp_eta == 1.0 ) {
00454     std::ostringstream omsg;
00455     omsg
00456       << "FeasibilityStepReducedStd_Strategy::compute_feasibility_step(...) : "
00457       << "Error, a QP relaxation parameter\n"
00458       << "of eta = " << qp_eta << " was calculated and therefore it must be assumed\n"
00459       << "that the NLP's constraints are infeasible\n"
00460       << "Throwing an InfeasibleConstraints exception!\n";
00461     if(  static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00462       out << omsg.str();
00463     }
00464     throw_qp_failure = true;
00465 //    throw InfeasibleConstraints(omsg.str());
00466   }
00467 
00468   //
00469   // Set the full step
00470   //
00471   // w = Ywy + Zwz
00472   //
00473   DenseLinAlgPack::V_VpV( w, Ywy, Zwz );
00474 
00475   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00476     out << "\n||w||inf    = " << DenseLinAlgPack::norm_inf(*w);
00477     out << std::endl;
00478   }
00479   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00480     out << "\nw = \n" << *w;
00481   }
00482 
00483   current_k_ = s->k();
00484 
00485   if( throw_qp_failure )
00486     return false;
00487 
00488 */
00489   TEST_FOR_EXCEPT(true);
00490 
00491   return true;
00492 }
00493 
00494 void FeasibilityStepReducedStd_Strategy::print_step( std::ostream& out, const std::string& L ) const
00495 {
00496   out << L << "*** Computes feasibility steps by solving a constrained QP using the range and null\n"
00497     << L << "*** space decomposition\n"
00498     << L << "begin quais-range space step: \"" << typeName(quasi_range_space_step()) << "\"\n";
00499 
00500    quasi_range_space_step().print_step( out, L + "  " );
00501 
00502   out << L << "end quasi-range space step\n"
00503     << L << "if quasi-range space step failed then\n"
00504     << L << "  this feasibility step fails also!\n"
00505     << L << "end\n"
00506     << L << "Ywy = v\n"
00507     << L << "*** Set the bounds for bl <= Z*wz <= bu\n"
00508     << L << "bl = d_bounds_k.l - (xo - x_k) - Ywy\n"
00509     << L << "bu = d_bounds_k.u - (xo - x_k) - Ywy\n"
00510     << L << "Set bl(i) = bu(i) for those nu_k(i) != 0.0\n"
00511     << L << "if (qp_objective == OBJ_MIN_FULL_STEP) and (current_k != k) then\n"
00512     << L << "  grad = Z_k'*Ywy\n"
00513     << L << "  Hess = Z_k'*Z_k\n"
00514     << L << "elseif (qp_objective == OBJ_MIN_NULL_SPACE_STEP) and (current_k != k) then\n"
00515     << L << "  grad = 0\n"
00516     << L << "  Hess = I\n"
00517     << L << "elseif (qp_objective == OBJ_RSQP) and (current_k != k) then\n"
00518     << L << "  grad = qp_grad_k\n"
00519     << L << "  Hess = rHL_k\n"
00520     << L << "end\n"
00521     << L << "if check_results == true then\n"
00522     << L << "  assert that bl and bu are valid and sorted\n"
00523     << L << "end\n"
00524     << L << "etaL = 0.0\n"
00525     << L << "*** Determine if we can use simple bounds on pz or not\n"
00526     << L << "if Z_k is a variable reduction null space matrix and norm(Ypy_k(indep),0) == 0 then\n"
00527     << L << "  use_simple_wz_bounds = true\n"
00528     << L << "else\n"
00529     << L << "  use_simple_wz_bounds = false\n"
00530     << L << "end\n"
00531     << L << "*** Setup QP arguments\n"
00532     << L << "qp_g = qp_grad_k\n"
00533     << L << "qp_G = rHL_k\n"
00534     << L << "if use_simple_wz_bounds == true then\n"
00535     << L << "  qp_dL = bl(indep),  qp_dU = bu(indep)\n"
00536     << L << "  qp_E  = Z_k.D,      qp_b  = Ywy(dep)\n"
00537     << L << "  qp_eL = bl(dep),    qp_eU = bu(dep)\n"
00538     << L << "else\n"
00539     << L << "  qp_dL = -inf,       qp_dU = +inf\n"
00540     << L << "  qp_E  = Z_k,        qp_b  = Ywy\n"
00541     << L << "  qp_eL = bl,         qp_eU = bu\n"
00542     << L << "end\n"
00543     << L << "if m > r then\n"
00544     << L << "  qp_F  = V_k,        qp_f  = c_k(undecomp) + Gc_k(undecomp)'*Ywy\n"
00545     << L << "else\n"
00546     << L << "  qp_F  = empty,      qp_f  = empty\n"
00547     << L << "end\n"
00548     << L << "Use a warm start given the active-set in nu_k\n"
00549     << L << "Solve the following QP to compute qp_d, qp_eta, qp_Ed = qp_E * qp_d\n"
00550     << L << ",qp_nu, qp_mu and qp_lambda (" << typeName(qp_solver()) << "):\n"
00551     << L << "  min    qp_g' * qp_d + 1/2 * qp_d' * qp_G * qp_d + M(eta)\n"
00552     << L << "  qp_d <: R^(n-r)\n"
00553     << L << "  s.t.\n"
00554     << L << "         etaL  <=  qp_eta\n"
00555     << L << "         qp_dL <= qp_d <= qp_dU                          [qp_nu]\n"
00556     << L << "         qp_eL <= qp_E * qp_d + (1-eta)*qp_b  <= qp_eU   [qp_mu]\n"
00557     << L << "         qp_F * d_qp + (1-eta) * qp_f = 0                [qp_lambda]\n"
00558     << L << "if (qp_teing==QP_TEST) or (fd_deriv_testing==QP_TEST_DEFAULT\n"
00559     << L << "and check_results==true) then\n"
00560     << L << "  Check the optimality conditions of the above QP\n"
00561     << L << "  if the optimality conditions do not check out then\n"
00562     << L << "    set throw_qp_failure = true\n"
00563     << L << "  end\n"
00564     << L << "end\n"
00565     << L << "*** Set the solution to the QP subproblem\n"
00566     << L << "wz  = qp_d\n"
00567     << L << "eta = qp_eta\n"
00568     << L << "if use_simple_wz_bounds == true then\n"
00569     << L << "  Zwz(dep)   = qp_Ed,  Zwz(indep) = wz\n"
00570     << L << "else\n"
00571     << L << "  Zwz = qp_Ed\n"
00572     << L << "end\n"
00573     << L << "if eta > 0 then set Ywy = (1-eta) * Ywy\n"
00574     << L << "if QP solution is suboptimal then\n"
00575     << L << "  throw_qp_failure = true\n"
00576     << L << "elseif QP solution is primal feasible (not optimal) then\n"
00577     << L << "  throw_qp_failure = primal_feasible_point_error\n"
00578     << L << "elseif QP solution is dual feasible (not optimal) then\n"
00579     << L << "  find max u s.t.\n"
00580     << L << "    d_bounds_k.l <= (xo - x) + u*(Ywy+Zwz) <= d_bounds_k.u\n"
00581     << L << "  alpha_k = u\n"
00582     << L << "  throw_qp_failure = true\n"
00583     << L << "end\n"
00584     << L << "if eta == 1.0 then\n"
00585     << L << "  The constraints are infeasible!\n"
00586     << L << "  throw_qp_failure = true\n"
00587     << L << "end\n"
00588     << L << "current_k = k\n"
00589     << L << "w = Zwz + Ywy\n"
00590     << L << "if (throw_qp_failure == true) then\n"
00591     << L << "  The feasibility step computation has failed!\n"
00592     << L << "end\n"
00593     ;
00594 }
00595 
00596 } // end namespace MoochoPack

Generated on Wed May 12 21:52:30 2010 for MOOCHO (Single Doxygen Collection) by  doxygen 1.4.7