MOOCHO (Single Doxygen Collection) Version of the Day
MoochoPack_MeritFunc_PenaltyParamUpdateGuts_AddedStep.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 // 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 <ostream>
00043 #include <typeinfo>
00044 
00045 #include "MoochoPack_MeritFunc_PenaltyParamUpdateGuts_AddedStep.hpp"
00046 #include "MoochoPack_moocho_algo_conversion.hpp"
00047 #include "IterationPack_print_algorithm_step.hpp"
00048 #include "ConstrainedOptPack_MeritFuncNLP.hpp"
00049 #include "ConstrainedOptPack_MeritFuncPenaltyParam.hpp"
00050 #include "ConstrainedOptPack_MeritFuncNLPDirecDeriv.hpp"
00051 #include "AbstractLinAlgPack_Vector.hpp"
00052 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00053 
00054 namespace {
00055 template< class T >
00056 inline
00057 T my_max( const T& v1, const T& v2 ) { return v1 > v2 ? v1 : v2; }
00058 } // end namespace
00059 
00060 namespace MoochoPack {
00061 
00062 MeritFunc_PenaltyParamUpdateGuts_AddedStep::MeritFunc_PenaltyParamUpdateGuts_AddedStep(
00063   value_type                     small_mu
00064   ,value_type                    mult_factor
00065   ,value_type                    kkt_near_sol
00066   )
00067   :near_solution_(false)
00068   ,small_mu_(small_mu)
00069   ,mult_factor_(mult_factor)
00070   ,kkt_near_sol_(kkt_near_sol)
00071 {}
00072 
00073 bool MeritFunc_PenaltyParamUpdateGuts_AddedStep::do_step(
00074   Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00075   ,poss_type assoc_step_poss
00076   )
00077 {
00078   NLPAlgo &algo = rsqp_algo(_algo);
00079   NLPAlgoState  &s    = algo.rsqp_state();
00080   NLP     &nlp  = algo.nlp();
00081   
00082   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00083   std::ostream& out = algo.track().journal_out();
00084 
00085   // print step header.
00086   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00087     using IterationPack::print_algorithm_step;
00088     print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
00089   }
00090 
00091   const size_type
00092     //n  = nlp.n(),
00093     m  = nlp.m();
00094   IterQuantityAccess<MeritFuncNLP>
00095     &merit_func_nlp_iq = s.merit_func_nlp();
00096 
00097   if( !merit_func_nlp_iq.updated_k(0) ) {
00098     const int merit_func_k_last_updated = merit_func_nlp_iq.last_updated();
00099     if( merit_func_k_last_updated != IterQuantity::NONE_UPDATED ) {
00100       MeritFuncNLP
00101         &merit_func_nlp_k_last = merit_func_nlp_iq.get_k(merit_func_k_last_updated);
00102       merit_func_nlp_iq.set_k(0) = merit_func_nlp_k_last;
00103     }
00104     else {
00105       merit_func_nlp_iq.set_k(0); // Just use default constructor
00106     }
00107     MeritFuncNLP
00108       &merit_func_nlp_k = merit_func_nlp_iq.get_k(0);
00109     MeritFuncPenaltyParam
00110       *param = dynamic_cast<MeritFuncPenaltyParam*>(&merit_func_nlp_k);
00111     TEUCHOS_TEST_FOR_EXCEPTION(
00112       !param, std::logic_error
00113       ,"MeritFunc_PenaltyParamUpdateGuts_AddedStep::do_step(...), Error "
00114       << "The class " << typeName(merit_func_nlp_k) << " does not support the "
00115       << "MeritFuncPenaltyParam iterface" );
00116     MeritFuncNLPDirecDeriv
00117       *direc_deriv = dynamic_cast<MeritFuncNLPDirecDeriv*>(&merit_func_nlp_k);
00118     TEUCHOS_TEST_FOR_EXCEPTION(
00119       !direc_deriv, std::logic_error
00120       ,"MeritFunc_PenaltyParamUpdateGuts_AddedStep::do_step(...), Error "
00121       << "The class " << typeName(merit_func_nlp_k) << " does not support the "
00122       << "MeritFuncNLPDirecDeriv iterface" );
00123     value_type  new_mu = 0.0;
00124     value_type  min_mu = 0.0;
00125     if ( this->min_mu(s,&min_mu) ) {
00126       // Update the penalty parameter as defined in the fortran rSQP code (EXACT2())
00127       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00128         out << "\nUpdate the penalty parameter...\n";
00129       }
00130       value_type
00131         mu_km1 = param->mu(),
00132         mult_fact = (1.0 + mult_factor_);
00133       if(near_solution_) {
00134         if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00135           out << "\nNear solution, forcing mu_k >= mu_km1...\n";
00136         }
00137         new_mu = my_max( my_max( mu_km1, mult_fact * min_mu ), small_mu_ );
00138       }
00139       else {
00140         if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00141           out << "\nNot near solution, allowing reduction in mu ...\n";
00142         }
00143         new_mu =  my_max(
00144           (3.0 * mu_km1 + min_mu) / 4.0 
00145           , my_max( mult_fact * min_mu, small_mu_ )
00146           ); 
00147         value_type
00148           kkt_error = s.opt_kkt_err().get_k(0) + s.feas_kkt_err().get_k(0);
00149         if(kkt_error <= kkt_near_sol_) {
00150           if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00151             out << "\nkkt_error = " << kkt_error << " <= kkt_near_sol = "
00152               << kkt_near_sol_ << std::endl
00153               << "Switching to forcing mu_k >= mu_km1 in the future\n";
00154           }
00155           near_solution_ = true;
00156         }
00157       }
00158     }
00159     else {
00160       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00161         out << "\nDon't have the info to update penalty parameter so just use the last updated...\n";
00162       }
00163       new_mu = param->mu();
00164     }
00165     // Set the penalty parameter
00166     param->mu( new_mu );
00167     // In addition also compute the directional derivative
00168     direc_deriv->calc_deriv(
00169       s.Gf().get_k(0)
00170       ,m  ? &s.c().get_k(0) : NULL
00171       ,NULL   // h
00172       ,NULL   // hl
00173       ,NULL   // hu
00174       ,s.d().get_k(0)
00175       );
00176     
00177     if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00178       out << "\nmu = " << new_mu << "\n";
00179     }
00180   }
00181   return true;
00182 }
00183 
00184 void MeritFunc_PenaltyParamUpdateGuts_AddedStep::print_step( const Algorithm& algo
00185   , poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
00186   , std::ostream& out, const std::string& L ) const
00187 {
00188   out
00189     << L << "*** Update the penalty parameter for the merit function to ensure\n"
00190     << L << "*** a descent direction a directional derivatieve.\n"
00191     << L << "*** phi is a merit function object that uses the penalty parameter mu.\n"
00192     << L << "default: near_solution = false\n"
00193     << L << "         small_mu = " << small_mu_ << std::endl
00194     << L << "         mult_factor = " << mult_factor_ << std::endl
00195     << L << "         kkt_near_sol = " << kkt_near_sol_ << std::endl
00196     << L << "if merit_func_nlp_k is not already updated then\n"
00197     << L << "    if some merit_func_nlp_k(?) has been udpated then\n"
00198     << L << "        merit_func_nlp_k = merit_func_nlp_k(last_udpated)\n"
00199     << L << "    else\n"
00200     << L << "        merit_func_nlp_k = default construction\n"
00201     << L << "    end\n"
00202     << L << "    if merit_func_nlp_k does not support MeritFuncPenaltyParam throw excpetion\n"
00203     << L << "    if merit_func_nlp_k does not support MeritFuncNLPDirecDeriv throw excpetion\n"
00204     ;
00205               print_min_mu_step( out, L + "    " ); 
00206   out
00207     << L << "   mu_new = merit_func_nlp_k.mu()\n"
00208     << L << "   if update_mu == true then\n"
00209     << L << "       mu_last = merit_func_nlp_k.mu()\n"
00210     << L << "       mult_fact = 1.0 + mult_factor\n"
00211     << L << "       if near_solution == true\n"
00212     << L << "           mu_new = max( max( mu_last, mult_fact*min_mu ), small_mu )\n"
00213     << L << "       else\n"
00214     << L << "           mu_new = max(   ( 3.0 * mu_last + min_mu ) / 4.0\n"
00215     << L << "                           , max( mult_fact * min_mu , small_mu ) )\n"
00216     << L << "           kkt_error = opt_kkt_err_k + feas_kkt_err_k\n"
00217     << L << "           if kkt_error <= kkt_near_sol then\n"
00218     << L << "               near_solution = true\n"
00219     << L << "           end\n"
00220     << L << "       end\n"
00221     << L << "   else\n"
00222     << L << "       mu_new = merit_func_nlp_k.mu()\n"
00223     << L << "   end\n"
00224     << L << "   merit_func_nlp_k..mu(mu_new)\n"
00225     << L << "   merit_func_nlp_k.calc_deriv(Gf_k,c_k,h_k,hl,hu,d_k)\n"
00226     << L << "end\n"
00227     ;
00228 }
00229 
00230 // Overridden from MeritFunc_PenaltyParamUpdate_AddedStep
00231 
00232 void MeritFunc_PenaltyParamUpdateGuts_AddedStep::small_mu( value_type small_mu )
00233 {
00234   small_mu_ = small_mu;
00235 }
00236 
00237 value_type MeritFunc_PenaltyParamUpdateGuts_AddedStep::small_mu() const
00238 {
00239   return small_mu_;
00240 }
00241 
00242 void MeritFunc_PenaltyParamUpdateGuts_AddedStep::mult_factor( value_type mult_factor )
00243 {
00244   mult_factor_ = mult_factor;
00245 }
00246 
00247 value_type MeritFunc_PenaltyParamUpdateGuts_AddedStep::mult_factor() const
00248 {
00249   return mult_factor_;
00250 }
00251 
00252 void MeritFunc_PenaltyParamUpdateGuts_AddedStep::kkt_near_sol( value_type kkt_near_sol )
00253 {
00254   kkt_near_sol_ = kkt_near_sol;
00255 }
00256 
00257 value_type MeritFunc_PenaltyParamUpdateGuts_AddedStep::kkt_near_sol() const
00258 {
00259   return kkt_near_sol_;
00260 }
00261 
00262 } // end namespace MoochoPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines