55 #include "ConstrainedOptPack/src/VectorWithNorms.h"
56 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixOp.hpp"
57 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_SpVectorClass.hpp"
58 #include "AbstractLinAlgPack/src/max_near_feas_step.hpp"
64 namespace LinAlgOpPack {
68 namespace MoochoPack {
71 const direct_ls_sqp_ptr_t& direct_ls_sqp
72 ,
const merit_func_ptr_t& merit_func
73 ,
const feasibility_step_ptr_t& feasibility_step
74 ,
const direct_ls_newton_ptr_t& direct_ls_newton
76 ,ENewtonOutputLevel newton_olevel
80 ,EForcedConstrReduction forced_constr_reduction
85 :direct_ls_sqp_(direct_ls_sqp)
86 ,merit_func_(merit_func)
87 ,feasibility_step_(feasibility_step)
88 ,direct_ls_newton_(direct_ls_newton)
90 ,newton_olevel_(newton_olevel)
91 ,constr_norm_threshold_(constr_norm_threshold)
92 ,constr_incr_ratio_(constr_incr_ratio)
93 ,after_k_iter_(after_k_iter)
94 ,forced_constr_reduction_(forced_constr_reduction)
95 ,forced_reduct_ratio_(forced_reduct_ratio)
96 ,max_step_ratio_(max_step_ratio)
97 ,max_newton_iter_(max_newton_iter)
102 ,poss_type assoc_step_poss
696 const Algorithm& algo
698 , std::ostream&
out,
const std::string& L )
const
700 out << L <<
"*** Calculate a second order correction when near solution.\n"
701 << L <<
"*** If we can compute a correction w then perform a curved\n"
702 << L <<
"*** line search along x_kp1 = x_k + alpha_k * d_k + alpha_k^2 * w.\n"
703 << L <<
"default: eta = " << eta() << std::endl
704 << L <<
" constr_norm_threshold = " << constr_norm_threshold() << std::endl
705 << L <<
" constr_incr_ratio = " << constr_incr_ratio() << std::endl
706 << L <<
" after_k_iter = " << after_k_iter() << std::endl
707 << L <<
" forced_constr_reduction = " << (forced_constr_reduction()==
CONSTR_LESS_X_D
708 ?
"CONSTR_LESS_X_D\n"
709 :
"CONSTR_LESS_X\n" )
710 << L <<
" forced_reduct_ratio = " << forced_reduct_ratio() << std::endl
711 << L <<
" max_step_ratio = " << max_step_ratio() << std::endl
712 << L <<
" max_newton_iter = " << max_newton_iter() << std::endl
713 << L <<
"begin definition of NLP merit function phi.value(f(x),c(x)):\n";
715 merit_func().print_merit_func( out, L +
" " );
717 out << L <<
"end definition\n"
718 << L <<
"Dphi_k = phi.deriv()\n"
719 << L <<
"if Dphi_k >= 0 then\n"
720 << L <<
" throw line_search_failure\n"
722 << L <<
"phi_kp1 = phi_k.value(f_kp1,c_kp1)\n"
723 << L <<
"phi_k = phi.value(f_k,c_k)\n"
724 << L <<
"if (norm(c_k,inf) <= constr_norm_threshold)\n"
725 << L <<
"and (norm(c_kp1,inf)/(norm(c_k,inf)+small_num) <= constr_incr_ratio)\n"
726 << L <<
"and (k >= after_k_iter) then\n"
727 << L <<
"considering_correction = ( (norm(c_k,inf) <= constr_norm_threshold)\n"
728 << L <<
" and (norm(c_kp1,inf)/(1.0 + norm(c_k,inf)) <= constr_incr_ratio)\n"
729 << L <<
" and (k >= after_k_iter) )\n"
730 << L <<
"chose_point = false\n"
731 << L <<
"if phi_kp1 < phi_k + eta * Dphi_k then\n"
732 << L <<
" chose_point = true\n"
734 << L <<
"if (considering_correction == true) and (chose_point == false) then\n"
735 << L <<
" considered_correction = true\n"
736 << L <<
" begin definition of c(x) merit function phi_c.value(c(x)):\n";
741 out << L <<
" end definition\n"
742 << L <<
" xdw = x_kp1;\n"
743 << L <<
" phi_c_x = phi_c.value(c_k);\n"
744 << L <<
" phi_c_xd = phi_c.value(c_kp1);\n"
745 << L <<
" phi_c_xdw = phi_c_xd;\n"
746 << L <<
" phi_c_xdww = phi_c_xdw;\n"
747 << L <<
" if phi_c_xd < forced_reduct_ratio * phi_c_x then\n"
748 << L <<
" *** There is no need to perform a correction.\n"
749 << L <<
" use_correction = false;\n"
751 << L <<
" *** Lets try to compute a correction by performing\n"
752 << L <<
" *** a series of newton steps to compute local steps w\n"
753 << L <<
" for newton_i = 1...max_newton_itr\n"
754 << L <<
" begin feasibility step calculation: \"" <<
typeName(feasibility_step()) <<
"\"\n";
756 feasibility_step().print_step( out, L +
" " );
758 out << L <<
" end feasibility step calculation\n"
759 << L <<
" Find the largest positive step u where the slightly\n"
760 << L <<
" relaxed variable bounds:\n"
761 << L <<
" xl - delta <= xdw + u * w <= xu + delta\n"
762 << L <<
" are strictly satisfied (where delta = max_var_bounds_viol).\n"
763 << L <<
" a = min(1,u);\n"
764 << L <<
" step_ratio = norm(w,inf)/norm(d,inf);\n"
765 << L <<
" a = min(a,max_step_ratio/step_ratio);\n"
766 << L <<
" Perform line search for phi_c.value(c(xdww = xdw+a*w))\n"
767 << L <<
" starting from a and compute:\n"
769 << L <<
" xdww = xdw + a * w,\n"
770 << L <<
" phi_c_xdww = phi_c.value(c(xdww))\n"
771 << L <<
" print summary information;\n"
772 << L <<
" if line search failed then\n"
773 << L <<
" use_correction = false;\n"
774 << L <<
" exit for loop;\n"
776 << L <<
" *** Determine if this is sufficent reduction in c(x) error\n"
777 << L <<
" if forced_constr_reduction == CONSTR_LESS_X_D then\n"
778 << L <<
" good_correction = (phi_c.value(c(xdww)) < forced_reduct_ratio*phi_c_xd);\n"
779 << L <<
" else if forced_constr_reduction == CONSTR_LESS_X then\n"
780 << L <<
" good_correction = (phi_c.value(c(xdww)) < forced_reduct_ratio*phi_c_x);\n"
782 << L <<
" if good_correction == true then\n"
783 << L <<
" w = xdww - (x_k+d_k);\n"
784 << L <<
" use_correction = true;\n"
785 << L <<
" exit for loop;\n"
787 << L <<
" *** This is not sufficient reduction is c(x) error yet.\n"
788 << L <<
" xdw = xdww;\n"
789 << L <<
" phi_c_xdw = phi_c_xdww;\n"
791 << L <<
" if use_correction == false then\n"
792 << L <<
" if forced_constr_reduction == CONSTR_LESS_X_D then\n"
793 << L <<
" *** Dam! We could not find a point phi_c_xdww < phi_c_xd.\n"
794 << L <<
" *** Perhaps Gc_k does not give a descent direction for phi_c!\n"
795 << L <<
" else if forced_constr_reduction == CONSTR_LESS_X then\n"
796 << L <<
" if phi_c_dww < phi_c_xd then\n"
797 << L <<
" *** Accept this correction anyway.\n"
798 << L <<
" use_correction = true\n"
800 << L <<
" *** Dam! we could not find any reduction in phi_c so\n"
801 << L <<
" *** Perhaps Gc_k does not give a descent direction for phi_c!\n"
806 << L <<
"if chose_point == false then\n"
807 << L <<
" if use_correction == true then\n"
808 << L <<
" Perform line search along x_kp1 = x_k + alpha_k * d_k + alpha_k^2 * w\n"
810 << L <<
" Perform line search along x_kp1 = x_k + alpha_k * d_k\n"
812 << L <<
" begin direct line search : \"" <<
typeName(direct_ls_sqp()) <<
"\"\n";
814 direct_ls_sqp().print_algorithm( out, L +
" " );
817 << L <<
" end direct line search\n"
818 << L <<
" if maximum number of linesearch iterations are exceeded then\n"
819 << L <<
" throw line_search_failure\n"
void print_merit_func(std::ostream &out, const std::string &leading_str) const
std::string typeName(const T &t)
bool do_step(Algorithm &algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss)
void Vp_StMtV(VectorMutable *v_lhs, value_type alpha, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta=1.0)
v_lhs = alpha * op(M_rhs1) * v_rhs2 + beta * v_lhs (BLAS xGEMV)
A merit function for the square of the constriant values.
AbstractLinAlgPack::value_type value_type
void print_step(const Algorithm &algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss, std::ostream &out, const std::string &leading_str) const
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
LineSearch2ndOrderCorrect_Step(const direct_ls_sqp_ptr_t &direct_ls_sqp=NULL, const merit_func_ptr_t &merit_func=NULL, const feasibility_step_ptr_t &feasibility_step=NULL, const direct_ls_newton_ptr_t &direct_ls_newton=0, value_type eta=1.0e-4, ENewtonOutputLevel newton_olevel=PRINT_USE_DEFAULT, value_type constr_norm_threshold=1.0, value_type constr_incr_ratio=10.0, int after_k_iter=0, EForcedConstrReduction forced_constr_reduction=CONSTR_LESS_X, value_type forced_reduct_ratio=1.0, value_type max_step_ratio=1.0, int max_newton_iter=3)