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