MoochoPack_TangentialStepWithInequStd_Step.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 <math.h>
00030 
00031 #include <ostream>
00032 #include <sstream>
00033 
00034 #include "MoochoPack_TangentialStepWithInequStd_Step.hpp"
00035 #include "MoochoPack_moocho_algo_conversion.hpp"
00036 #include "MoochoPack_Exceptions.hpp"
00037 #include "IterationPack_print_algorithm_step.hpp"
00038 #include "ConstrainedOptPack_MatrixIdentConcat.hpp"
00039 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00040 #include "AbstractLinAlgPack_VectorOut.hpp"
00041 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
00042 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00043 #include "Teuchos_dyn_cast.hpp"
00044 
00045 namespace {
00046 template< class T >
00047 inline
00048 T my_max( const T& v1, const T& v2 ) { return v1 > v2 ? v1 : v2; }
00049 template< class T >
00050 inline
00051 T my_min( const T& v1, const T& v2 ) { return v1 < v2 ? v1 : v2; }
00052 } // end namespace
00053 
00054 namespace MoochoPack {
00055 
00056 TangentialStepWithInequStd_Step::TangentialStepWithInequStd_Step(
00057   const qp_solver_ptr_t       &qp_solver
00058   ,const qp_tester_ptr_t      &qp_tester
00059   ,value_type                 warm_start_frac
00060   ,EQPTesting                 qp_testing
00061   ,bool                       primal_feasible_point_error
00062   ,bool                       dual_feasible_point_error
00063   )
00064   :qp_solver_(qp_solver)
00065   ,qp_tester_(qp_tester)
00066   ,warm_start_frac_(warm_start_frac)
00067   ,qp_testing_(qp_testing)
00068   ,primal_feasible_point_error_(primal_feasible_point_error)
00069   ,dual_feasible_point_error_(dual_feasible_point_error)
00070   ,dl_iq_(dl_name)
00071   ,du_iq_(du_name)
00072 {}
00073 
00074 bool TangentialStepWithInequStd_Step::do_step(
00075   Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00076   ,poss_type assoc_step_poss
00077   )
00078 {
00079   namespace mmp = MemMngPack;
00080   using Teuchos::dyn_cast;
00081   using ::fabs;
00082   using LinAlgOpPack::Vt_S;
00083   using LinAlgOpPack::V_VpV;
00084   using LinAlgOpPack::V_VmV;
00085   using LinAlgOpPack::Vp_StV;
00086   using LinAlgOpPack::Vp_V;
00087   using LinAlgOpPack::V_StV;
00088   using LinAlgOpPack::V_MtV;
00089 //  using ConstrainedOptPack::min_abs;
00090   using AbstractLinAlgPack::max_near_feas_step;
00091   typedef VectorMutable::vec_mut_ptr_t   vec_mut_ptr_t;
00092 
00093   NLPAlgo &algo = rsqp_algo(_algo);
00094   NLPAlgoState &s = algo.rsqp_state();
00095   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00096   EJournalOutputLevel ns_olevel = algo.algo_cntr().null_space_journal_output_level();
00097   std::ostream &out = algo.track().journal_out();
00098   //const bool check_results = algo.algo_cntr().check_results();
00099 
00100   // print step header.
00101   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00102     using IterationPack::print_algorithm_step;
00103     print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
00104   }
00105 
00106   // problem dimensions
00107   const size_type
00108     //n  = algo.nlp().n(),
00109     m  = algo.nlp().m(),
00110     r  = s.equ_decomp().size();
00111 
00112   // Get the iteration quantity container objects
00113   IterQuantityAccess<value_type>
00114     &alpha_iq = s.alpha(),
00115     &zeta_iq  = s.zeta(),
00116     &eta_iq   = s.eta();
00117   IterQuantityAccess<VectorMutable>
00118     &dl_iq      = dl_iq_(s),
00119     &du_iq      = du_iq_(s),
00120     &nu_iq      = s.nu(),
00121     *c_iq       = m > 0  ? &s.c()       : NULL,
00122     *lambda_iq  = m > 0  ? &s.lambda()  : NULL,
00123     &rGf_iq     = s.rGf(),
00124     &w_iq       = s.w(),
00125     &qp_grad_iq = s.qp_grad(),
00126     &py_iq      = s.py(),
00127     &pz_iq      = s.pz(),
00128     &Ypy_iq     = s.Ypy(),
00129     &Zpz_iq     = s.Zpz();
00130   IterQuantityAccess<MatrixOp>
00131     &Z_iq   = s.Z(),
00132     //*Uz_iq  = (m > r)  ? &s.Uz() : NULL,
00133     *Uy_iq  = (m > r)  ? &s.Uy() : NULL;
00134   IterQuantityAccess<MatrixSymOp>
00135     &rHL_iq = s.rHL();
00136   IterQuantityAccess<ActSetStats>
00137     &act_set_stats_iq = act_set_stats_(s);
00138   
00139   // Accessed/modified/updated (just some)
00140   VectorMutable  *Ypy_k = (m ? &Ypy_iq.get_k(0) : NULL);
00141   const MatrixOp  &Z_k   = Z_iq.get_k(0);
00142   VectorMutable  &pz_k  = pz_iq.set_k(0);
00143   VectorMutable  &Zpz_k = Zpz_iq.set_k(0);
00144 
00145   // Comupte qp_grad which is an approximation to rGf + Z'*HL*Y*py
00146 
00147   // qp_grad = rGf
00148   VectorMutable
00149     &qp_grad_k = ( qp_grad_iq.set_k(0) = rGf_iq.get_k(0) );
00150 
00151   // qp_grad += zeta * w
00152   if( w_iq.updated_k(0) ) {
00153     if(zeta_iq.updated_k(0))
00154       Vp_StV( &qp_grad_k, zeta_iq.get_k(0), w_iq.get_k(0) );
00155     else
00156       Vp_V( &qp_grad_k, w_iq.get_k(0) );
00157   }
00158 
00159   // Set the bounds for:
00160   //   dl <= Z*pz + Y*py <= du  ->  dl - Ypy <= Z*pz <= du - Ypz
00161 
00162   vec_mut_ptr_t
00163     bl = s.space_x().create_member(),
00164     bu = s.space_x().create_member();
00165 
00166   if(m) {
00167     // bl = dl_k - Ypy_k
00168     V_VmV( bl.get(), dl_iq.get_k(0), *Ypy_k );
00169     // bu = du_k - Ypy_k
00170     V_VmV( bu.get(), du_iq.get_k(0), *Ypy_k );
00171   }
00172   else {
00173     *bl = dl_iq.get_k(0);
00174     *bu = du_iq.get_k(0);
00175   }
00176 
00177   // Print out the QP bounds for the constraints
00178   if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00179     out << "\nqp_grad_k = \n" << qp_grad_k;
00180   }
00181   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00182     out << "\nbl = \n" << *bl;
00183     out << "\nbu = \n" << *bu;
00184   }
00185 
00186   //
00187   // Determine if we should perform a warm start or not.
00188   //
00189   bool do_warm_start = false;
00190   if( act_set_stats_iq.updated_k(-1) ) {
00191     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00192       out << "\nDetermining if the QP should use a warm start ...\n";
00193     }
00194     // We need to see if we should preform a warm start for the next iteration
00195     ActSetStats &stats = act_set_stats_iq.get_k(-1);
00196     const size_type
00197       num_active = stats.num_active(),
00198       num_adds   = stats.num_adds(),
00199       num_drops  = stats.num_drops();
00200     const value_type
00201       frac_same
00202       = ( num_adds == ActSetStats::NOT_KNOWN || num_active == 0
00203         ? 0.0
00204         : my_max(((double)(num_active)-num_adds-num_drops) / num_active, 0.0 ) );
00205     do_warm_start = ( num_active > 0 && frac_same >= warm_start_frac() );
00206     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00207       out << "\nnum_active = " << num_active;
00208       if( num_active ) {
00209         out << "\nmax(num_active-num_adds-num_drops,0)/(num_active) = "
00210           << "max("<<num_active<<"-"<<num_adds<<"-"<<num_drops<<",0)/("<<num_active<<") = "
00211           << frac_same;
00212         if( do_warm_start )
00213           out << " >= ";
00214         else
00215           out << " < ";
00216         out << "warm_start_frac = " << warm_start_frac();
00217       }
00218       if( do_warm_start )
00219         out << "\nUse a warm start this time!\n";
00220       else
00221         out << "\nDon't use a warm start this time!\n";
00222     }
00223   }
00224 
00225   // Use active set from last iteration as an estimate for current active set
00226   // if we are to use a warm start.
00227   // 
00228   // ToDo: If the selection of dependent and independent variables changes
00229   // then you will have to adjust this or not perform a warm start at all!
00230   if( do_warm_start ) {
00231     nu_iq.set_k(0,-1);
00232   }
00233   else {
00234     nu_iq.set_k(0) = 0.0; // No guess of the active set
00235   }
00236   VectorMutable  &nu_k  = nu_iq.get_k(0);
00237 
00238   //
00239   // Setup the reduced QP subproblem
00240   //
00241   // The call to the QP is setup for the more flexible call to the QPSolverRelaxed
00242   // interface to deal with the three independent variabilities: using simple
00243   // bounds for pz or not, general inequalities included or not, and extra equality
00244   // constraints included or not.
00245   // If this method of calling the QP solver were not used then 4 separate
00246   // calls to solve_qp(...) would have to be included to handle the four possible
00247   // QP formulations.
00248   //
00249 
00250   // The numeric arguments for the QP solver (in the nomenclatrue of QPSolverRelaxed)
00251 
00252   const value_type  qp_bnd_inf = NLP::infinite_bound();
00253 
00254   const Vector            &qp_g       = qp_grad_k;
00255   const MatrixSymOp       &qp_G       = rHL_iq.get_k(0);
00256   const value_type        qp_etaL     = 0.0;
00257   vec_mut_ptr_t           qp_dL       = Teuchos::null;
00258   vec_mut_ptr_t           qp_dU       = Teuchos::null;
00259   Teuchos::RefCountPtr<const MatrixOp>
00260                           qp_E        = Teuchos::null;
00261   BLAS_Cpp::Transp        qp_trans_E  = BLAS_Cpp::no_trans;
00262   vec_mut_ptr_t           qp_b        = Teuchos::null;
00263   vec_mut_ptr_t           qp_eL       = Teuchos::null;
00264   vec_mut_ptr_t           qp_eU       = Teuchos::null;
00265   Teuchos::RefCountPtr<const MatrixOp>
00266                           qp_F        = Teuchos::null;
00267   BLAS_Cpp::Transp        qp_trans_F  = BLAS_Cpp::no_trans;
00268   vec_mut_ptr_t           qp_f        = Teuchos::null;
00269   value_type              qp_eta      = 0.0;
00270   VectorMutable           &qp_d       = pz_k;  // pz_k will be updated directly!
00271   vec_mut_ptr_t           qp_nu       = Teuchos::null;
00272   vec_mut_ptr_t           qp_mu       = Teuchos::null;
00273   vec_mut_ptr_t           qp_Ed       = Teuchos::null;
00274   vec_mut_ptr_t           qp_lambda   = Teuchos::null;
00275 
00276   //
00277   // Determine if we can use simple bounds on pz.
00278   // 
00279   // If we have a variable-reduction null-space matrix
00280   // (with any choice for Y) then:
00281   // 
00282   // d = Z*pz + (1-eta) * Y*py
00283   // 
00284   // [ d(var_dep)   ]  = [ D ] * pz  + (1-eta) * [ Ypy(var_dep)   ]
00285   // [ d(var_indep) ]    [ I ]                   [ Ypy(var_indep) ]
00286   // 
00287   // For a cooridinate decomposition (Y = [ I ; 0 ]) then Ypy(var_indep) == 0.0 and
00288   // in this case the bounds on d(var_indep) become simple bounds on pz even
00289   // with the relaxation.
00290   // 
00291   // Otherwise, we can not use simple variable bounds and implement the
00292   // relaxation properly.
00293   // 
00294 
00295   const MatrixIdentConcat
00296     *Zvr = dynamic_cast<const MatrixIdentConcat*>( &Z_k );
00297   const Range1D
00298     var_dep   = Zvr ? Zvr->D_rng() : Range1D::Invalid,
00299     var_indep = Zvr ? Zvr->I_rng() : Range1D();
00300 
00301   const value_type Ypy_indep_norm_inf = ( m ? Ypy_k->sub_view(var_indep)->norm_inf() : 0.0);
00302 
00303   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
00304     out
00305       << "\nDetermine if we can use simple bounds on pz ...\n"
00306       << "    m = " << m << std::endl
00307       << "    dynamic_cast<const MatrixIdentConcat*>(&Z_k) = " << Zvr << std::endl
00308       << "    ||Ypy_k(var_indep)||inf = " << Ypy_indep_norm_inf << std::endl;
00309 
00310   const bool
00311     use_simple_pz_bounds = ( m == 0 || ( Zvr != NULL && Ypy_indep_norm_inf == 0.0 ) ),
00312     bounded_var_dep      = ( m > 0 && num_bounded( *bl->sub_view(var_dep), *bu->sub_view(var_dep), qp_bnd_inf ) );
00313 
00314   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
00315     out
00316       << (use_simple_pz_bounds
00317           ? "\nUsing simple bounds on pz ...\n"
00318           : "\nUsing bounds on full Z*pz ...\n")
00319       << (bounded_var_dep
00320           ? "\nThere are finite bounds on dependent variables.  Adding extra inequality constrints for D*pz ...\n" 
00321           : "\nThere are no finite bounds on dependent variables.  There will be no extra inequality constraints added on D*pz ...\n" ) ;
00322 
00323   if( use_simple_pz_bounds ) {
00324     // Set simple bound constraints on pz
00325     qp_dL = bl->sub_view(var_indep);
00326     qp_dU = bu->sub_view(var_indep);
00327     qp_nu = nu_k.sub_view(var_indep); // nu_k(var_indep) will be updated directly!
00328     if( m && bounded_var_dep ) {
00329       // Set general inequality constraints for D*pz
00330       qp_E   = Teuchos::rcp(&Zvr->D(),false);
00331       qp_b   = Ypy_k->sub_view(var_dep);
00332       qp_eL  = bl->sub_view(var_dep);
00333       qp_eU  = bu->sub_view(var_dep);
00334       qp_mu  = nu_k.sub_view(var_dep);  // nu_k(var_dep) will be updated directly!
00335       qp_Ed  = Zpz_k.sub_view(var_dep); // Zpz_k(var_dep) will be updated directly!
00336     }
00337     else {
00338       // Leave these as NULL since there is no extra general inequality constraints
00339     }
00340   }
00341   else if( !use_simple_pz_bounds ) { // ToDo: Leave out parts for unbounded dependent variables!
00342     // There are no simple bounds! (leave qp_dL, qp_dU and qp_nu as null)
00343     // Set general inequality constraints for Z*pz
00344     qp_E   = Teuchos::rcp(&Z_k,false);
00345     qp_b   = Teuchos::rcp(Ypy_k,false);
00346     qp_eL  = bl;
00347     qp_eU  = bu;
00348     qp_mu  = Teuchos::rcp(&nu_k,false);
00349     qp_Ed  = Teuchos::rcp(&Zpz_k,false); // Zpz_k will be updated directly!
00350   }
00351   else {
00352     TEST_FOR_EXCEPT(true);
00353   }
00354 
00355   // Set the general equality constriants (if they exist)
00356   Range1D equ_undecomp = s.equ_undecomp();
00357   if( m > r && m > 0 ) {
00358     // qp_f = Uy_k * py_k + c_k(equ_undecomp)
00359     qp_f = s.space_c().sub_space(equ_undecomp)->create_member();
00360     V_MtV( qp_f.get(), Uy_iq->get_k(0), BLAS_Cpp::no_trans, py_iq.get_k(0) );
00361     Vp_V( qp_f.get(), *c_iq->get_k(0).sub_view(equ_undecomp) );
00362     // Must resize for the undecomposed constriants if it has not already been
00363     qp_F       = Teuchos::rcp(&Uy_iq->get_k(0),false);
00364     qp_lambda  = lambda_iq->set_k(0).sub_view(equ_undecomp); // lambda_k(equ_undecomp), will be updated directly!
00365   }
00366 
00367   // Setup the rest of the arguments
00368 
00369   QPSolverRelaxed::EOutputLevel
00370     qp_olevel;
00371   switch( olevel ) {
00372     case PRINT_NOTHING:
00373       qp_olevel = QPSolverRelaxed::PRINT_NONE;
00374       break;
00375     case PRINT_BASIC_ALGORITHM_INFO:
00376       qp_olevel = QPSolverRelaxed::PRINT_NONE;
00377       break;
00378     case PRINT_ALGORITHM_STEPS:
00379       qp_olevel = QPSolverRelaxed::PRINT_BASIC_INFO;
00380       break;
00381     case PRINT_ACTIVE_SET:
00382       qp_olevel = QPSolverRelaxed::PRINT_ITER_SUMMARY;
00383       break;
00384     case PRINT_VECTORS:
00385       qp_olevel = QPSolverRelaxed::PRINT_ITER_VECTORS;
00386       break;
00387     case PRINT_ITERATION_QUANTITIES:
00388       qp_olevel = QPSolverRelaxed::PRINT_EVERY_THING;
00389       break;
00390     default:
00391       TEST_FOR_EXCEPT(true);
00392   }
00393   // ToDo: Set print options so that only vectors matrices etc
00394   // are only printed in the null space
00395 
00396   //
00397   // Solve the QP
00398   // 
00399   qp_solver().infinite_bound(qp_bnd_inf);
00400   const QPSolverStats::ESolutionType
00401     solution_type =
00402     qp_solver().solve_qp(
00403       int(olevel) == int(PRINT_NOTHING) ? NULL : &out
00404       ,qp_olevel
00405       ,( algo.algo_cntr().check_results()
00406          ? QPSolverRelaxed::RUN_TESTS :  QPSolverRelaxed::NO_TESTS )
00407       ,qp_g, qp_G, qp_etaL, qp_dL.get(), qp_dU.get()
00408       ,qp_E.get(), qp_trans_E, qp_E.get() ? qp_b.get() : NULL
00409       ,qp_E.get() ? qp_eL.get() : NULL, qp_E.get() ? qp_eU.get() : NULL
00410       ,qp_F.get(), qp_trans_F, qp_F.get() ? qp_f.get() : NULL
00411       ,NULL // obj_d
00412       ,&qp_eta, &qp_d
00413       ,qp_nu.get()
00414       ,qp_mu.get(), qp_E.get() ? qp_Ed.get() : NULL
00415       ,qp_F.get() ? qp_lambda.get() : NULL
00416       ,NULL // qp_Fd
00417       );
00418 
00419   //
00420   // Check the optimality conditions for the QP
00421   //
00422   std::ostringstream omsg;
00423   bool throw_qp_failure = false;
00424   if(   qp_testing() == QP_TEST
00425     || ( qp_testing() == QP_TEST_DEFAULT && algo.algo_cntr().check_results() )  )
00426   {
00427     if( int(olevel) >= int(PRINT_ALGORITHM_STEPS) ) {
00428       out << "\nChecking the optimality conditions of the reduced QP subproblem ...\n";
00429     }
00430     if(!qp_tester().check_optimality_conditions(
00431       solution_type,qp_solver().infinite_bound()
00432       ,int(olevel) == int(PRINT_NOTHING) ? NULL : &out
00433       ,int(olevel) >= int(PRINT_VECTORS) ? true : false
00434       ,int(olevel) >= int(PRINT_ITERATION_QUANTITIES) ? true : false
00435       ,qp_g, qp_G, qp_etaL, qp_dL.get(), qp_dU.get()
00436       ,qp_E.get(), qp_trans_E, qp_E.get() ? qp_b.get() : NULL
00437       ,qp_E.get() ? qp_eL.get() : NULL, qp_E.get() ? qp_eU.get() : NULL
00438       ,qp_F.get(), qp_trans_F, qp_F.get() ? qp_f.get() : NULL
00439       ,NULL // obj_d
00440       ,&qp_eta, &qp_d
00441       ,qp_nu.get()
00442       ,qp_mu.get(), qp_E.get() ? qp_Ed.get() : NULL
00443       ,qp_F.get() ? qp_lambda.get() : NULL
00444       ,NULL // qp_Fd
00445       ))
00446     {
00447       omsg << "\n*** Alert! at least one of the QP optimality conditions did not check out.\n";
00448       if(  static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00449         out << omsg.str();
00450       }
00451       throw_qp_failure = true;
00452     }
00453   }
00454 
00455   //
00456   // Set the solution
00457   //
00458   if( !use_simple_pz_bounds ) {
00459     // Everything is already updated!
00460   }
00461   else if( use_simple_pz_bounds ) {
00462     // Just have to set Zpz_k(var_indep) = pz_k
00463     *Zpz_k.sub_view(var_indep) = pz_k;
00464     if( m && !bounded_var_dep ) {
00465       // Must compute Zpz_k(var_dep) = D*pz
00466       LinAlgOpPack::V_MtV( &*Zpz_k.sub_view(var_dep), Zvr->D(), BLAS_Cpp::no_trans, pz_k );
00467       // ToDo: Remove the compuation of Zpz here unless you must
00468     }
00469   }
00470   else {
00471     TEST_FOR_EXCEPT(true);
00472   }
00473 
00474   // Set the solution statistics
00475   qp_solver_stats_(s).set_k(0) = qp_solver().get_qp_stats();
00476 
00477   // Cut back Ypy_k = (1-eta) * Ypy_k
00478   const value_type eps = std::numeric_limits<value_type>::epsilon();
00479   if( fabs(qp_eta - 0.0) > eps ) {
00480     if(  static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00481       out
00482         << "\n*** Alert! the QP was infeasible (eta = "<<qp_eta<<").  Cutting back Ypy_k = (1.0 - eta)*Ypy_k  ...\n";
00483     }
00484     Vt_S( Ypy_k, 1.0 - qp_eta );
00485   }
00486 
00487   // eta_k
00488   eta_iq.set_k(0) = qp_eta;
00489 
00490   //
00491   // Modify the solution if we have to!
00492   // 
00493   switch(solution_type) {
00494     case QPSolverStats::OPTIMAL_SOLUTION:
00495       break;  // we are good!
00496     case QPSolverStats::PRIMAL_FEASIBLE_POINT:
00497     {
00498       omsg
00499         << "\n*** Alert! the returned QP solution is PRIMAL_FEASIBLE_POINT but not optimal!\n";
00500       if( primal_feasible_point_error() )
00501         omsg
00502           << "\n*** primal_feasible_point_error == true, this is an error!\n";
00503       if(  static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00504         out << omsg.str();
00505       }
00506       throw_qp_failure = primal_feasible_point_error();
00507       break;
00508     } 
00509     case QPSolverStats::DUAL_FEASIBLE_POINT:
00510     {
00511       omsg
00512         << "\n*** Alert! the returned QP solution is DUAL_FEASIBLE_POINT"
00513         << "\n*** but not optimal so we cut back the step ...\n";
00514       if( dual_feasible_point_error() )
00515         omsg
00516           << "\n*** dual_feasible_point_error == true, this is an error!\n";
00517       if(  static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00518         out << omsg.str();
00519       }
00520       // Cut back the step to fit in the bounds
00521       // 
00522       // dl <= u*(Ypy_k+Zpz_k) <= du
00523       //
00524       vec_mut_ptr_t
00525         zero  = s.space_x().create_member(0.0),
00526         d_tmp = s.space_x().create_member();
00527       V_VpV( d_tmp.get(), *Ypy_k, Zpz_k );
00528       const std::pair<value_type,value_type>
00529         u_steps = max_near_feas_step( *zero, *d_tmp, dl_iq.get_k(0), du_iq.get_k(0), 0.0 );
00530       const value_type
00531         u = my_min( u_steps.first, 1.0 ); // largest positive step size
00532       alpha_iq.set_k(0) = u;
00533       if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00534         out << "\nFinding u s.t. dl <= u*(Ypy_k+Zpz_k) <= du\n"
00535           << "max step length u = " << u << std::endl
00536           << "alpha_k = u = " << alpha_iq.get_k(0) << std::endl;
00537       }
00538       throw_qp_failure = dual_feasible_point_error();
00539       break;
00540     } 
00541     case QPSolverStats::SUBOPTIMAL_POINT:
00542     {
00543       omsg
00544         << "\n*** Alert!, the returned QP solution is SUBOPTIMAL_POINT!\n";
00545       if(  static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00546         out << omsg.str();
00547       }
00548       throw_qp_failure = true;
00549       break;
00550     }
00551     default:
00552       TEST_FOR_EXCEPT(true);  // should not happen!
00553   }
00554 
00555   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00556     out << "\n||pz_k||inf    = " << s.pz().get_k(0).norm_inf()
00557       << "\nnu_k.nz()      = " << s.nu().get_k(0).nz()
00558       << "\nmax(|nu_k(i)|) = " << s.nu().get_k(0).norm_inf()
00559 //      << "\nmin(|nu_k(i)|) = " << min_abs( s.nu().get_k(0)() )
00560       ;
00561     if( m > r ) out << "\n||lambda_k(undecomp)||inf = " << s.lambda().get_k(0).norm_inf();
00562     out << "\n||Zpz_k||2     = " << s.Zpz().get_k(0).norm_2()
00563       ;
00564     if(qp_eta > 0.0) out << "\n||Ypy||2 = " << s.Ypy().get_k(0).norm_2();
00565     out << std::endl;
00566   }
00567 
00568   if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00569     out << "\npz_k = \n" << s.pz().get_k(0);
00570   }
00571 
00572   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00573     out << "\nnu_k = \n" << s.nu().get_k(0);
00574     if( m > r )
00575       out << "\nlambda_k(equ_undecomp) = \n" << *s.lambda().get_k(0).sub_view(equ_undecomp);
00576     out << "\nZpz_k = \n" << s.Zpz().get_k(0);
00577     if(qp_eta > 0.0) out << "\nYpy = \n" << s.Ypy().get_k(0);
00578   }
00579 
00580   if( qp_eta == 1.0 ) {
00581     omsg
00582       << "TangentialStepWithInequStd_Step::do_step(...) : Error, a QP relaxation parameter\n"
00583       << "of eta = " << qp_eta << " was calculated and therefore it must be assumed\n"
00584       << "that the NLP's constraints are infeasible\n"
00585       << "Throwing an InfeasibleConstraints exception!\n";
00586     if(  static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00587       out << omsg.str();
00588     }
00589     throw InfeasibleConstraints(omsg.str());
00590   }
00591 
00592   if( throw_qp_failure )
00593     throw QPFailure( omsg.str(), qp_solver().get_qp_stats() );
00594 
00595   return true;
00596 }
00597 
00598 void TangentialStepWithInequStd_Step::print_step(
00599   const Algorithm& algo, poss_type step_poss, IterationPack::EDoStepType type
00600   ,poss_type assoc_step_poss, std::ostream& out, const std::string& L
00601   ) const
00602 {
00603   out
00604     << L << "*** Calculate the null-space step by solving a constrained QP\n"
00605     << L << "qp_grad_k = rGf_k\n"
00606     << L << "if w_k is updated then set qp_grad_k = qp_grad_k + zeta_k * w_k\n"
00607     << L << "bl = dl_k - Ypy_k\n"
00608     << L << "bu = du_k - Ypy_k\n"
00609     << L << "etaL = 0.0\n"
00610     << L << "*** Determine if we can use simple bounds on pz or not\n"
00611     << L << "if m==0 or (Z_k is a variable reduction null space matrix and ||Ypy_k(var_indep)||inf == 0) then\n"
00612     << L << "  use_simple_pz_bounds = true\n"
00613     << L << "else\n"
00614     << L << "  use_simple_pz_bounds = false\n"
00615     << L << "end\n"
00616     << L << "*** Setup QP arguments\n"
00617     << L << "qp_g = qp_grad_k\n"
00618     << L << "qp_G = rHL_k\n"
00619     << L << "if (use_simple_pz_bounds == true) then\n"
00620     << L << "  qp_dL = bl(var_indep),  qp_dU = bu(var_indep)\n"
00621     << L << "  if (m > 0) then\n"
00622     << L << "    qp_E  = Z_k.D,          qp_b  = Ypy_k(var_dep)\n"
00623     << L << "    qp_eL = bl(var_dep),    qp_eU = bu(var_dep)\n"
00624     << L << "  end\n"
00625     << L << "elseif (use_simple_pz_bounds == false) then\n"
00626     << L << "  qp_dL = -inf,           qp_dU = +inf\n"
00627     << L << "  qp_E  = Z_k,            qp_b  = Ypy_k\n"
00628     << L << "  qp_eL = bl,             qp_eU = bu\n"
00629     << L << "end\n"
00630     << L << "if (m > r) then\n"
00631     << L << "  qp_F  = V_k,            qp_f  = Uy_k*py_k + c_k(equ_undecomp)\n"
00632     << L << "else\n"
00633     << L << "  qp_F  = empty,          qp_f  = empty\n"
00634     << L << "end\n"
00635     << L << "Given active-set statistics (act_set_stats_km1)\n"
00636     << L << "  frac_same = max(num_active-num_adds-num_drops,0)/(num_active)\n"
00637     << L << "Use a warm start when frac_same >= warm_start_frac\n"
00638     << L << "Solve the following QP to compute qp_d, qp_eta, qp_Ed = qp_E * qp_d\n"
00639     << L << ",qp_nu, qp_mu and qp_lambda (" << typeid(qp_solver()).name() << "):\n"
00640     << L << "  min  qp_g' * qp_d + 1/2 * qp_d' * qp_G * qp_d + M(eta)\n"
00641     << L << "  qp_d <: R^(n-r)\n"
00642     << L << "  s.t.\n"
00643     << L << "       etaL  <=  eta\n"
00644     << L << "       qp_dL <= qp_d                         <= qp_dU   [qp_nu]\n"
00645     << L << "       qp_eL <= qp_E * qp_d + (1-eta)*qp_b   <= qp_eU   [qp_mu]\n"
00646     << L << "                qp_F * d_qp + (1-eta) * qp_f  = 0       [qp_lambda]\n"
00647     << L << "if (qp_testing==QP_TEST) or (fd_deriv_testing==QP_TEST_DEFAULT\n"
00648     << L << "and check_results==true) then\n"
00649     << L << "  Check the optimality conditions of the above QP\n"
00650     << L << "  if the optimality conditions do not check out then\n"
00651     << L << "    set throw_qp_failure = true\n"
00652     << L << "  end\n"
00653     << L << "end\n"
00654     << L << "*** Set the solution to the QP subproblem\n"
00655     << L << "pz_k  = qp_d\n"
00656     << L << "eta_k = qp_eta\n"
00657     << L << "if (use_simple_pz_bounds == true) then\n"
00658     << L << "  nu_k(var_dep)  = qp_mu,  nu_k(var_indep)  = qp_nu\n"
00659     << L << "  Zpz_k(var_dep) = qp_Ed,  Zpz_k(var_indep) = pz_k\n"
00660     << L << "elseif (use_simple_pz_bounds == false) then\n"
00661     << L << "  nu_k  = qp_mu\n"
00662     << L << "  Zpz_k = qp_Ed\n"
00663     << L << "end\n"
00664     << L << "if m > r then\n"
00665     << L << "  lambda_k(equ_undecomp) = qp_lambda\n"
00666     << L << "end\n"
00667     << L << "if (eta_k > 0) then set Ypy_k = (1-eta_k) * Ypy_k\n"
00668     << L << "if QP solution is suboptimal then\n"
00669     << L << "  throw_qp_failure = true\n"
00670     << L << "elseif QP solution is primal feasible (not optimal) then\n"
00671     << L << "  throw_qp_failure = primal_feasible_point_error\n"
00672     << L << "elseif QP solution is dual feasible (not optimal) then\n"
00673     << L << "  find max u s.t.\n"
00674     << L << "    dl_k <= u*(Ypy_k+Zpz_k) <= du_k\n"
00675     << L << "  alpha_k = u\n"
00676     << L << "  throw_qp_failure = dual_feasible_point_error\n"
00677     << L << "end\n"
00678     << L << "if (eta_k == 1.0) then\n"
00679     << L << "  The constraints are infeasible!\n"
00680     << L << "  throw InfeasibleConstraints(...)\n"
00681     << L << "end\n"
00682     << L << "if (throw_qp_failure == true) then\n"
00683     << L << "  throw QPFailure(...)\n"
00684     << L << "end\n"
00685     ;
00686 }
00687 
00688 } // end namespace MoochoPack

Generated on Thu Sep 18 12:35:18 2008 for MOOCHO (Single Doxygen Collection) by doxygen 1.3.9.1