source: branches/dev/Algorithm/IpFilterLSAcceptor.cpp @ 542

Last change on this file since 542 was 542, checked in by andreasw, 14 years ago
  • cleaned up line search to allow for alternative globalization scheme
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 27.2 KB
RevLine 
[438]1// Copyright (C) 2004, 2005 International Business Machines and others.
[2]2// All Rights Reserved.
3// This code is published under the Common Public License.
4//
5// $Id: IpFilterLSAcceptor.cpp 542 2005-10-13 22:43:08Z andreasw $
6//
7// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
[542]8//           Andreas Waechter                 IBM    2005-10-13
9//               derived file from IpFilterLineSearch.cpp
[2]10
[542]11#include "IpFilterLSAcceptor.hpp"
[2]12#include "IpJournalist.hpp"
13#include "IpRestoPhase.hpp"
[387]14#include "IpAlgTypes.hpp"
[2]15
[438]16#ifdef HAVE_CMATH
17# include <cmath>
[2]18#else
[438]19# ifdef HAVE_MATH_H
20#  include <math.h>
21# else
22#  error "don't have header file for math"
23# endif
24#endif
25
[2]26namespace Ipopt
27{
28
[465]29#ifdef IP_DEBUG
30  static const Index dbg_verbosity = 0;
31#endif
[2]32
[542]33  FilterLSAcceptor::FilterLSAcceptor(const SmartPtr<PDSystemSolver>& pd_solver)
[2]34      :
[422]35      filter_(2),
[542]36      pd_solver_(pd_solver)
[177]37  {
[542]38    DBG_START_FUN("FilterLSAcceptor::FilterLSAcceptor",
[224]39                  dbg_verbosity);
[177]40  }
[2]41
[542]42  FilterLSAcceptor::~FilterLSAcceptor()
[177]43  {
[542]44    DBG_START_FUN("FilterLSAcceptor::~FilterLSAcceptor()",
[224]45                  dbg_verbosity);
[177]46  }
[2]47
[542]48  void FilterLSAcceptor::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
[2]49  {
[336]50    roptions->AddLowerBoundedNumberOption(
51      "theta_max_fact",
52      "Determines upper bound for constraint violation in the filter.",
53      0.0, true, 1e4,
54      "The algorithmic parameter theta_max is determined as theta_max_fact "
55      "times the maximum of 1 and the constraint violation at initial point.  "
56      "Any point with a constraint violation larger than theta_max is "
57      "unacceptable to the filter (see Eqn. (21) in implementation paper).");
58    roptions->AddLowerBoundedNumberOption(
59      "theta_min_fact",
[490]60      "Determines constraint violation threshold in the switching rule.",
[336]61      0.0, true, 1e-4,
[490]62      "The algorithmic parameter theta_min is determined as theta_min_fact "
[336]63      "times the maximum of 1 and the constraint violation at initial point.  "
[490]64      "The switching rules treats an iteration as an h-type iteration whenever "
[336]65      "the current constraint violation is larger than theta_min (see "
66      "paragraph before Eqn. (19) in implementation paper).");
67    roptions->AddBoundedNumberOption(
68      "eta_phi",
69      "Relaxation factor in the Armijo condition.",
[402]70      0.0, true, 0.5, true, 1e-8,
[336]71      "(See Eqn. (20) in implementation paper)");
72    roptions->AddLowerBoundedNumberOption(
[490]73      "delta", "Multiplier for constraint violation in the switching rule.",
[336]74      0.0, true, 1.0,
[493]75      "(See Eqn. (19) in implementation paper.)");
[336]76    roptions->AddLowerBoundedNumberOption(
77      "s_phi",
[490]78      "Exponent for linear barrier function model in the switching rule.",
[336]79      1.0, true, 2.3,
[493]80      "(See Eqn. (19) in implementation paper.)");
[336]81    roptions->AddLowerBoundedNumberOption(
82      "s_theta",
[490]83      "Exponent for current constraint violation in the switching rule.",
[336]84      1.0, true, 1.1,
[493]85      "(See Eqn. (19) in implementation paper.)");
[336]86    roptions->AddBoundedNumberOption(
87      "gamma_phi",
[490]88      "Relaxation factor in the filter margin for the barrier function.",
[336]89      0.0, true, 1.0, true, 1e-8,
[493]90      "(See Eqn. (18a) in implementation paper.)");
[336]91    roptions->AddBoundedNumberOption(
92      "gamma_theta",
[490]93      "Relaxation factor in the filter margin for the constraint violation.",
[336]94      0.0, true, 1.0, true, 1e-5,
[493]95      "(See Eqn. (18b) in implementation paper.)");
[336]96    roptions->AddBoundedNumberOption(
97      "alpha_min_frac",
[490]98      "Safety factor for the minimal step size (before switching to restoration phase).",
[336]99      0.0, true, 1.0, true, 0.05,
[493]100      "(This is gamma_alpha in Eqn. (20) in implementation paper.)");
[336]101    roptions->AddLowerBoundedIntegerOption(
102      "max_soc",
[490]103      "Maximum number of second order correction trial steps at each iteration.",
[336]104      0, 4,
[490]105      "Choosing 0 disables the second order "
[336]106      "corrections. (This is p^{max} of Step A-5.9 of "
107      "Algorithm A in implementation paper.)");
108    roptions->AddLowerBoundedNumberOption(
109      "kappa_soc",
[490]110      "Factor in the sufficient reduction rule for second order correction.",
[336]111      0.0, true, 0.99,
[490]112      "This option determines how much a second order correction step must reduce the "
[336]113      "constraint violation so that further correction steps are attempted.  "
114      "(See Step A-5.9 of Algorithm A in implementation paper.)");
115    roptions->AddLowerBoundedNumberOption(
116      "obj_max_inc",
[490]117      "Determines the upper bound on the acceptable increase of barrier objective function.",
[336]118      1.0, true, 5.0,
[493]119      "Trial points are rejected if they lead to an increase in the "
[490]120      "barrier objective function by more than obj_max_inc orders "
121      "of magnitude.");
[493]122
[336]123    roptions->AddStringOption3(
124      "corrector_type",
[490]125      "The type of corrector steps that should be taken.",
[336]126      "none",
127      "none", "no corrector",
128      "affine", "corrector step towards mu=0",
129      "primal-dual", "corrector step towards current mu",
[490]130      "If \"mu_strategy\" is \"adaptive\", this option determines "
131      "what kind of corrector steps should be tried.");
[2]132
[336]133    roptions->AddStringOption2(
134      "skip_corr_if_neg_curv",
135      "Skip the corrector step in negative curvature iteration.",
136      "yes",
137      "no", "don't skip",
138      "yes", "skip",
[490]139      "The corrector step is not tried if negative curvature has been "
140      "encountered during the computation of the search direction in "
141      "the current iteration. This option is only used if \"mu_strategy\" is "
142      "\"adaptive\".");
[2]143
[336]144    roptions->AddStringOption2(
145      "skip_corr_in_monotone_mode",
146      "Skip the corrector step during monotone barrier parameter mode.",
147      "yes",
148      "no", "don't skip",
149      "yes", "skip",
150      "The corrector step is not tried if the algorithm is currently in the "
[490]151      "monotone mode (see also option \"barrier_strategy\")."
152      "This option is only used if \"mu_strategy\" is \"adaptive\".");
[2]153
[336]154    roptions->AddLowerBoundedNumberOption(
155      "corrector_compl_avrg_red_fact",
[490]156      "Complementarity tolerance factor for accepting corrector step.",
[336]157      0.0, true, 1.0,
[490]158      "This option determines the factor by which complementarity is allowed to increase "
[336]159      "for a corrector step to be accepted.");
[310]160  }
[2]161
[542]162  bool FilterLSAcceptor::InitializeImpl(const OptionsList& options,
[310]163                                        const std::string& prefix)
164  {
165    options.GetNumericValue("theta_max_fact", theta_max_fact_, prefix);
166    options.GetNumericValue("theta_min_fact", theta_min_fact_, prefix);
[433]167    ASSERT_EXCEPTION(theta_min_fact_ < theta_max_fact_, OPTION_INVALID,
[321]168                     "Option \"theta_min_fact\": This value must be larger than 0 and less than theta_max_fact.");
[310]169    options.GetNumericValue("eta_phi", eta_phi_, prefix);
170    options.GetNumericValue("delta", delta_, prefix);
171    options.GetNumericValue("s_phi", s_phi_, prefix);
172    options.GetNumericValue("s_theta", s_theta_, prefix);
173    options.GetNumericValue("gamma_phi", gamma_phi_, prefix);
174    options.GetNumericValue("gamma_theta", gamma_theta_, prefix);
175    options.GetNumericValue("alpha_min_frac", alpha_min_frac_, prefix);
176    options.GetIntegerValue("max_soc", max_soc_, prefix);
[2]177    if (max_soc_>0) {
[433]178      ASSERT_EXCEPTION(IsValid(pd_solver_), OPTION_INVALID,
[542]179                       "Option \"max_soc\": This option is non-negative, but no linear solver for computing the SOC given to FilterLSAcceptor object.");
[2]180    }
[310]181    options.GetNumericValue("kappa_soc", kappa_soc_, prefix);
182    options.GetNumericValue("obj_max_inc", obj_max_inc_, prefix);
[314]183    Index enum_int;
184    options.GetEnumValue("corrector_type", enum_int, prefix);
185    corrector_type_ = CorrectorTypeEnum(enum_int);
[310]186    options.GetBoolValue("skip_corr_if_neg_curv", skip_corr_if_neg_curv_, prefix);
[336]187    options.GetBoolValue("skip_corr_in_monotone_mode", skip_corr_in_monotone_mode_, prefix);
[310]188    options.GetNumericValue("corrector_compl_avrg_red_fact", corrector_compl_avrg_red_fact_, prefix);
[401]189
[520]190    theta_min_ = -1.;
191    theta_max_ = -1.;
[152]192
[273]193    Reset();
194
[542]195    return true;
[2]196  }
197
[542]198  void FilterLSAcceptor::InitThisLineSearch(bool in_watchdog)
[2]199  {
[542]200    DBG_START_METH("FilterLSAcceptor::InitThisLineSearch",
[2]201                   dbg_verbosity);
202
[273]203    // Set the values for the reference point
[542]204    if (!in_watchdog) {
[273]205      reference_theta_ = IpCq().curr_constraint_violation();
206      reference_barr_ = IpCq().curr_barrier_obj();
207      reference_gradBarrTDelta_ = IpCq().curr_gradBarrTDelta();
208    }
209    else {
[336]210      reference_theta_ = watchdog_theta_;
211      reference_barr_ = watchdog_barr_;
212      reference_gradBarrTDelta_ = watchdog_gradBarrTDelta_;
[273]213    }
214    filter_.Print(Jnlst());
215  }
216
[542]217  bool FilterLSAcceptor::IsFtype(Number alpha_primal_test)
[2]218  {
[542]219    DBG_START_METH("FilterLSAcceptor::IsFtype",
[2]220                   dbg_verbosity);
[412]221    DBG_ASSERT(reference_theta_>0. || reference_gradBarrTDelta_ < 0.0);
[273]222    return (reference_gradBarrTDelta_ < 0.0 &&
223            alpha_primal_test*pow(-reference_gradBarrTDelta_,s_phi_) >
224            delta_*pow(reference_theta_,s_theta_));
[2]225  }
226
[542]227  void FilterLSAcceptor::AugmentFilter()
[2]228  {
[542]229    DBG_START_METH("FilterLSAcceptor::AugmentFilter",
[2]230                   dbg_verbosity);
231
[273]232    Number phi_add = reference_barr_ - gamma_phi_*reference_theta_;
233    Number theta_add = (1.-gamma_theta_)*reference_theta_;
[2]234
235    filter_.AddEntry(phi_add, theta_add, IpData().iter_count());
236  }
237
238  bool
[542]239  FilterLSAcceptor::CheckAcceptabilityOfTrialPoint(Number alpha_primal_test)
[2]240  {
[542]241    DBG_START_METH("FilterLSAcceptor::CheckAcceptabilityOfTrialPoint",
[2]242                   dbg_verbosity);
243
244    bool accept;
245
246    // First compute the barrier function and constraint violation at the
247    // current iterate and the trial point
248
249    Number trial_theta = IpCq().trial_constraint_violation();
250    // Check if constraint violation is becoming too large
251    if (theta_max_ < 0.0) {
[273]252      theta_max_ = theta_max_fact_*Max(1.0, reference_theta_);
[2]253    }
254    if (theta_min_ < 0.0) {
[273]255      theta_min_ = theta_min_fact_*Max(1.0, reference_theta_);
[2]256    }
257
258    if (theta_max_>0 && trial_theta>theta_max_) {
[488]259      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
260                     "trial_theta = %e is larger than theta_max = %e\n",
261                     trial_theta, theta_max_);
262      IpData().Append_info_string("Tmax");
[2]263      return false;
264    }
265
266    Number trial_barr = IpCq().trial_barrier_obj();
[438]267    DBG_ASSERT(IsFiniteNumber(trial_barr));
[2]268
269    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
270                   "Checking acceptability for trial step size alpha_primal_test=%13.6e:\n", alpha_primal_test);
271    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
[273]272                   "  New values of barrier function     = %23.16e  (reference %23.16e):\n", trial_barr, reference_barr_);
[2]273    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
[273]274                   "  New values of constraint violation = %23.16e  (reference %23.16e):\n", trial_theta, reference_theta_);
[2]275
276    // Check if point is acceptable w.r.t current iterate
[542]277    if (alpha_primal_test>0. && IsFtype(alpha_primal_test) &&
278        reference_theta_ <= theta_min_) {
[2]279      // Armijo condition for the barrier function has to be satisfied
280      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Checking Armijo Condition...\n");
281      accept = ArmijoHolds(alpha_primal_test);
282    }
283    else {
284      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Checking sufficient reduction...\n");
285      accept = IsAcceptableToCurrentIterate(trial_barr, trial_theta);
286    }
287
288    if (!accept) {
289      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Failed...\n");
290      return accept;
291    }
292    else {
293      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Succeeded...\n");
294    }
295
296    // Now check if that pair is acceptable to the filter
297    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Checking filter acceptability...\n");
298    accept = IsAcceptableToCurrentFilter(trial_barr, trial_theta);
299    if (!accept) {
300      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Failed...\n");
301      return accept;
302    }
303    else {
304      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Succeeded...\n");
305    }
306
307    return accept;
308  }
309
[542]310  bool FilterLSAcceptor::ArmijoHolds(Number alpha_primal_test)
[2]311  {
[401]312    /*
313    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
314                   "ArmijoHolds test with trial_barr = %25.16e reference_barr = %25.16e\n        alpha_primal_test = %25.16e reference_gradBarrTDelta = %25.16e\n", IpCq().trial_barrier_obj(), reference_barr_,alpha_primal_test,reference_gradBarrTDelta_);
315    */
[273]316    return Compare_le(IpCq().trial_barrier_obj()-reference_barr_,
317                      eta_phi_*alpha_primal_test*reference_gradBarrTDelta_,
318                      reference_barr_);
[2]319  }
320
[542]321  Number FilterLSAcceptor::CalculateAlphaMin()
[2]322  {
323    Number gBD = IpCq().curr_gradBarrTDelta();
324    Number curr_theta = IpCq().curr_constraint_violation();
325    Number alpha_min = gamma_theta_;
326
327    if (gBD < 0) {
328      alpha_min = Min( gamma_theta_,
329                       gamma_phi_*curr_theta/(-gBD));
330      if (curr_theta <= theta_min_) {
331        alpha_min = Min( alpha_min,
332                         delta_*pow(curr_theta,s_theta_)/pow(-gBD,s_phi_)
333                       );
334      }
335    }
336
337    return alpha_min_frac_*alpha_min;
338  }
339
[542]340  bool FilterLSAcceptor::IsAcceptableToCurrentIterate(Number trial_barr,
[224]341      Number trial_theta,
342      bool called_from_restoration /*=false*/) const
[2]343  {
[542]344    DBG_START_METH("FilterLSAcceptor::IsAcceptableToCurrentIterate",
[2]345                   dbg_verbosity);
[11]346
347    // Check if the barrier objective function is increasing to
348    // rapidly (according to option obj_max_inc)
[273]349    if (!called_from_restoration && trial_barr > reference_barr_) {
[11]350      Number basval = 1.;
[273]351      if (fabs(reference_barr_)>10.) {
352        basval = log10(fabs(reference_barr_));
[11]353      }
[273]354      if (log10(trial_barr-reference_barr_)>obj_max_inc_+basval) {
355        Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Rejecting trial point because barrier objective function increasing too rapidly (from %27.15e to %27.15e)\n",reference_barr_,trial_barr);
[11]356        return false;
357      }
358    }
359
[273]360    DBG_PRINT((1,"trial_barr  = %e reference_barr  = %e\n", trial_barr, reference_barr_));
361    DBG_PRINT((1,"trial_theta = %e reference_theta = %e\n", trial_theta, reference_theta_));
362    return (Compare_le(trial_theta, (1.-gamma_theta_)*reference_theta_, reference_theta_)
363            || Compare_le(trial_barr-reference_barr_, -gamma_phi_*reference_theta_, reference_barr_));
[2]364  }
365
[542]366  bool FilterLSAcceptor::IsAcceptableToCurrentFilter(Number trial_barr, Number trial_theta) const
[2]367  {
368    return filter_.Acceptable(trial_barr, trial_theta);
369  }
370
[542]371  bool FilterLSAcceptor::Compare_le(Number lhs, Number rhs, Number BasVal)
[2]372  {
[542]373    DBG_START_FUN("FilterLSAcceptor::Compare_le",
[2]374                  dbg_verbosity);
[65]375    DBG_PRINT((1,"lhs = %27.16e rhs = %27.16e  BasVal = %27.16e\n",lhs,rhs,BasVal));
[2]376
[542]377    Number mach_eps = std::numeric_limits<Number>::epsilon();
378    return (lhs - rhs <= 1e10*mach_eps*fabs(BasVal));
[264]379  }
380
[542]381  void FilterLSAcceptor::StartWatchDog()
[273]382  {
[542]383    DBG_START_FUN("FilterLSAcceptor::StartWatchDog", dbg_verbosity);
[273]384
[336]385    watchdog_theta_ = IpCq().curr_constraint_violation();
386    watchdog_barr_ = IpCq().curr_barrier_obj();
387    watchdog_gradBarrTDelta_ = IpCq().curr_gradBarrTDelta();
[273]388  }
389
[542]390  void FilterLSAcceptor::StopWatchDog()
[273]391  {
[542]392    DBG_START_FUN("FilterLSAcceptor::StopWatchDog", dbg_verbosity);
[273]393
[336]394    reference_theta_ = watchdog_theta_;
395    reference_barr_ = watchdog_barr_;
396    reference_gradBarrTDelta_ = watchdog_gradBarrTDelta_;
[273]397  }
398
[542]399  void FilterLSAcceptor::Reset()
[2]400  {
[542]401    DBG_START_FUN("FilterLSAcceptor::Reset", dbg_verbosity);
[273]402
[2]403    filter_.Clear();
404  }
405
[49]406  bool
[542]407  FilterLSAcceptor::TrySecondOrderCorrection(
[56]408    Number alpha_primal_test,
409    Number& alpha_primal,
[305]410    SmartPtr<IteratesVector>& actual_delta)
[2]411  {
[542]412    DBG_START_METH("FilterLSAcceptor::TrySecondOrderCorrection",
[49]413                   dbg_verbosity);
[2]414
[542]415    if (max_soc_==0) {
416      return false;
417    }
418
[49]419    bool accept = false;
420    Index count_soc = 0;
421
422    Number theta_soc_old = 0.;
[264]423    Number theta_trial = IpCq().trial_constraint_violation();
[49]424    Number alpha_primal_soc = alpha_primal;
425
426    SmartPtr<Vector> c_soc = IpCq().curr_c()->MakeNew();
427    SmartPtr<Vector> dms_soc = IpCq().curr_d_minus_s()->MakeNew();
428    c_soc->Copy(*IpCq().curr_c());
429    dms_soc->Copy(*IpCq().curr_d_minus_s());
[264]430    while (count_soc<max_soc_ && !accept &&
431           (count_soc==0 || theta_trial<=kappa_soc_*theta_soc_old) ) {
432      theta_soc_old = theta_trial;
[49]433
434      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
[56]435                     "Trying second order correction number %d\n",
436                     count_soc+1);
[49]437
438      // Compute SOC constraint violation
[198]439      c_soc->AddOneVector(1.0, *IpCq().trial_c(), alpha_primal_soc);
440      dms_soc->AddOneVector(1.0, *IpCq().trial_d_minus_s(), alpha_primal_soc);
[49]441
442      // Compute the SOC search direction
[305]443      SmartPtr<IteratesVector> delta_soc = actual_delta->MakeNewIteratesVector(true);
444      SmartPtr<IteratesVector> rhs = actual_delta->MakeNewContainer();
445      rhs->Set_x(*IpCq().curr_grad_lag_with_damping_x());
446      rhs->Set_s(*IpCq().curr_grad_lag_with_damping_s());
447      rhs->Set_y_c(*c_soc);
448      rhs->Set_y_d(*dms_soc);
449      rhs->Set_z_L(*IpCq().curr_relaxed_compl_x_L());
450      rhs->Set_z_U(*IpCq().curr_relaxed_compl_x_U());
451      rhs->Set_v_L(*IpCq().curr_relaxed_compl_s_L());
452      rhs->Set_v_U(*IpCq().curr_relaxed_compl_s_U());
[307]453      pd_solver_->Solve(-1.0, 0.0, *rhs, *delta_soc, true);
[49]454
455      // Compute step size
456      alpha_primal_soc =
[56]457        IpCq().primal_frac_to_the_bound(IpData().curr_tau(),
[305]458                                        *delta_soc->x(),
459                                        *delta_soc->s());
[49]460
461      // Check if trial point is acceptable
462      try {
[56]463        // Compute the primal trial point
[321]464        IpData().SetTrialPrimalVariablesFromStep(alpha_primal_soc, *delta_soc->x(), *delta_soc->s());
[49]465
[56]466        // in acceptance tests, use original step size!
467        accept = CheckAcceptabilityOfTrialPoint(alpha_primal_test);
[49]468      }
469      catch(IpoptNLP::Eval_Error& e) {
[498]470        e.ReportException(Jnlst(), J_DETAILED);
[56]471        Jnlst().Printf(J_WARNING, J_MAIN, "Warning: SOC step rejected due to evaluation error\n");
[498]472        IpData().Append_info_string("e");
[56]473        accept = false;
[416]474        // There is no point in continuing SOC procedure
475        break;
[49]476      }
477
478      if (accept) {
[56]479        Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Second order correction step accepted with %d corrections.\n", count_soc+1);
480        // Accept all SOC quantities
481        alpha_primal = alpha_primal_soc;
[321]482        actual_delta = delta_soc;
[49]483      }
484      else {
[56]485        count_soc++;
486        theta_trial = IpCq().trial_constraint_violation();
[49]487      }
488    }
489    return accept;
[2]490  }
491
[49]492  bool
[542]493  FilterLSAcceptor::TryCorrector(
[56]494    Number alpha_primal_test,
495    Number& alpha_primal,
[305]496    SmartPtr<IteratesVector>& actual_delta)
[49]497  {
[542]498    if (corrector_type_==NO_CORRECTOR ||
499        (skip_corr_if_neg_curv_ && IpData().info_regu_x()!=0.) ||
500        (skip_corr_in_monotone_mode_ && !IpData().FreeMuMode())) {
501      return false;
502    }
503
504    DBG_START_METH("FilterLSAcceptor::TryCorrector",
[49]505                   dbg_verbosity);
506
[305]507    Index n_bounds = IpData().curr()->z_L()->Dim() + IpData().curr()->z_U()->Dim()
508                     + IpData().curr()->v_L()->Dim() + IpData().curr()->v_U()->Dim();
[264]509    if (n_bounds==0) {
510      // Nothing to be done
511      return false;
512    }
513
[529]514    IpData().TimingStats().TryCorrector().Start();
[523]515
[49]516    bool accept = false;
517
518    // Compute the corrector step based on corrector_type parameter
[305]519    // create a new iterates vector and allocate space for all the entries
520    SmartPtr<IteratesVector> delta_corr = actual_delta->MakeNewIteratesVector(true);
[49]521
522    switch (corrector_type_) {
[336]523      case AFFINE_CORRECTOR : {
[65]524        // 1: Standard MPC corrector
[49]525
[295]526        if (!IpData().HaveAffineDeltas()) {
527          Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
528                         "Solving the Primal Dual System for the affine step\n");
529          // First get the right hand side
[321]530          SmartPtr<IteratesVector> rhs_aff = delta_corr->MakeNewContainer();
531
[305]532          rhs_aff->Set_x(*IpCq().curr_grad_lag_x());
533          rhs_aff->Set_s(*IpCq().curr_grad_lag_s());
534          rhs_aff->Set_y_c(*IpCq().curr_c());
535          rhs_aff->Set_y_d(*IpCq().curr_d_minus_s());
536          rhs_aff->Set_z_L(*IpCq().curr_compl_x_L());
[321]537          rhs_aff->Set_z_U(*IpCq().curr_compl_x_U());
538          rhs_aff->Set_v_L(*IpCq().curr_compl_s_L());
539          rhs_aff->Set_v_U(*IpCq().curr_compl_s_U());
[295]540
[321]541          // create a new iterates vector (with allocated space)
542          // for the affine scaling step
543          SmartPtr<IteratesVector> step_aff = delta_corr->MakeNewIteratesVector(true);
[295]544
545          // Now solve the primal-dual system to get the step
[307]546          pd_solver_->Solve(-1.0, 0.0, *rhs_aff, *step_aff, false);
[295]547
[321]548          DBG_PRINT_VECTOR(2, "step_aff", *step_aff);
[295]549
[321]550          IpData().set_delta_aff(step_aff);
[295]551          IpData().SetHaveAffineDeltas(true);
552        }
553
[224]554        DBG_ASSERT(IpData().HaveAffineDeltas());
[65]555
[321]556        const SmartPtr<const IteratesVector> delta_aff = IpData().delta_aff();
[49]557
[321]558        delta_corr->Copy(*actual_delta);
[49]559
[321]560        // create a rhs vector and allocate entries
561        SmartPtr<IteratesVector> rhs = actual_delta->MakeNewIteratesVector(true);
562
[305]563        rhs->x_NonConst()->Set(0.);
564        rhs->s_NonConst()->Set(0.);
565        rhs->y_c_NonConst()->Set(0.);
566        rhs->y_d_NonConst()->Set(0.);
567        IpNLP().Px_L()->TransMultVector(-1., *delta_aff->x(), 0., *rhs->z_L_NonConst());
568        rhs->z_L_NonConst()->ElementWiseMultiply(*delta_aff->z_L());
569        IpNLP().Px_U()->TransMultVector(1., *delta_aff->x(), 0., *rhs->z_U_NonConst());
570        rhs->z_U_NonConst()->ElementWiseMultiply(*delta_aff->z_U());
571        IpNLP().Pd_L()->TransMultVector(-1., *delta_aff->s(), 0., *rhs->v_L_NonConst());
572        rhs->v_L_NonConst()->ElementWiseMultiply(*delta_aff->v_L());
573        IpNLP().Pd_U()->TransMultVector(1., *delta_aff->s(), 0., *rhs->v_U_NonConst());
574        rhs->v_U_NonConst()->ElementWiseMultiply(*delta_aff->v_U());
[49]575
[307]576        pd_solver_->Solve(1.0, 1.0, *rhs, *delta_corr, true);
[49]577
[321]578        DBG_PRINT_VECTOR(2, "delta_corr", *delta_corr);
[49]579      }
580      break;
[336]581      case PRIMAL_DUAL_CORRECTOR : {
[65]582        // 2: Second order correction for primal-dual step to
583        // primal-dual mu
584
[321]585        delta_corr->Copy(*actual_delta);
[65]586
[321]587        // allocate space for the rhs
588        SmartPtr<IteratesVector> rhs = actual_delta->MakeNewIteratesVector(true);
589
[305]590        rhs->x_NonConst()->Set(0.);
591        rhs->s_NonConst()->Set(0.);
592        rhs->y_c_NonConst()->Set(0.);
593        rhs->y_d_NonConst()->Set(0.);
[65]594
595        Number mu = IpData().curr_mu();
596        SmartPtr<Vector> tmp;
597
[305]598        rhs->z_L_NonConst()->Copy(*IpCq().curr_slack_x_L());
599        IpNLP().Px_L()->TransMultVector(-1., *actual_delta->x(),
600                                        -1., *rhs->z_L_NonConst());
601        tmp = actual_delta->z_L()->MakeNew();
602        tmp->AddTwoVectors(1., *IpData().curr()->z_L(), 1., *actual_delta->z_L(), 0.);
603        rhs->z_L_NonConst()->ElementWiseMultiply(*tmp);
604        rhs->z_L_NonConst()->AddScalar(mu);
[65]605
[305]606        rhs->z_U_NonConst()->Copy(*IpCq().curr_slack_x_U());
607        IpNLP().Px_U()->TransMultVector(1., *actual_delta->x(),
608                                        -1., *rhs->z_U_NonConst());
609        tmp = actual_delta->z_U()->MakeNew();
610        tmp->AddTwoVectors(1., *IpData().curr()->z_U(), 1., *actual_delta->z_U(), 0.);
611        rhs->z_U_NonConst()->ElementWiseMultiply(*tmp);
612        rhs->z_U_NonConst()->AddScalar(mu);
[65]613
[305]614        rhs->v_L_NonConst()->Copy(*IpCq().curr_slack_s_L());
615        IpNLP().Pd_L()->TransMultVector(-1., *actual_delta->s(),
616                                        -1., *rhs->v_L_NonConst());
617        tmp = actual_delta->v_L()->MakeNew();
618        tmp->AddTwoVectors(1., *IpData().curr()->v_L(), 1., *actual_delta->v_L(), 0.);
619        rhs->v_L_NonConst()->ElementWiseMultiply(*tmp);
620        rhs->v_L_NonConst()->AddScalar(mu);
[65]621
[305]622        rhs->v_U_NonConst()->Copy(*IpCq().curr_slack_s_U());
623        IpNLP().Pd_U()->TransMultVector(1., *actual_delta->s(),
624                                        -1., *rhs->v_U_NonConst());
625        tmp = actual_delta->v_U()->MakeNew();
626        tmp->AddTwoVectors(1., *IpData().curr()->v_U(), 1., *actual_delta->v_U(), 0.);
627        rhs->v_U_NonConst()->ElementWiseMultiply(*tmp);
628        rhs->v_U_NonConst()->AddScalar(mu);
[65]629
[321]630        DBG_PRINT_VECTOR(2, "rhs", *rhs);
[65]631
[307]632        pd_solver_->Solve(1.0, 1.0, *rhs, *delta_corr, true);
[65]633
[321]634        DBG_PRINT_VECTOR(2, "delta_corr", *delta_corr);
[65]635      }
636      break;
[56]637      default:
[49]638      DBG_ASSERT("Unknown corrector_type value.");
639    }
640
641    // Compute step size
642    Number alpha_primal_corr =
643      IpCq().primal_frac_to_the_bound(IpData().curr_tau(),
[305]644                                      *delta_corr->x(),
645                                      *delta_corr->s());
646    // Set the primal trial point
647    IpData().SetTrialPrimalVariablesFromStep(alpha_primal_corr, *delta_corr->x(), *delta_corr->s());
[49]648
649    // Check if we want to not even try the filter criterion
650    Number alpha_dual_max =
651      IpCq().dual_frac_to_the_bound(IpData().curr_tau(),
[305]652                                    *delta_corr->z_L(), *delta_corr->z_U(),
653                                    *delta_corr->v_L(), *delta_corr->v_U());
[49]654
[305]655    IpData().SetTrialBoundMultipliersFromStep(alpha_dual_max, *delta_corr->z_L(), *delta_corr->z_U(), *delta_corr->v_L(), *delta_corr->v_U());
[49]656
657    Number trial_avrg_compl = IpCq().trial_avrg_compl();
658    Number curr_avrg_compl = IpCq().curr_avrg_compl();
659    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
[56]660                   "avrg_compl(curr) = %e, avrg_compl(trial) = %e\n",
661                   curr_avrg_compl, trial_avrg_compl);
[336]662    if (corrector_type_==AFFINE_CORRECTOR &&
[65]663        trial_avrg_compl>=corrector_compl_avrg_red_fact_*curr_avrg_compl) {
[49]664      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
[65]665                     "Rejecting corrector step, because trial complementarity is too large.\n" );
[529]666      IpData().TimingStats().TryCorrector().End();
[49]667      return false;
668    }
669
670    // Check if trial point is acceptable
671    try {
672      // in acceptance tests, use original step size!
673      accept = CheckAcceptabilityOfTrialPoint(alpha_primal_test);
674    }
675    catch(IpoptNLP::Eval_Error& e) {
[498]676      e.ReportException(Jnlst(), J_DETAILED);
[336]677      Jnlst().Printf(J_WARNING, J_MAIN,
678                     "Warning: Corrector step rejected due to evaluation error\n");
[498]679      IpData().Append_info_string("e");
[49]680      accept = false;
681    }
682
683    if (accept) {
684      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
[56]685                     "Corrector step accepted with alpha_primal = %e\n",
686                     alpha_primal_corr);
[49]687      // Accept all SOC quantities
688      alpha_primal = alpha_primal_corr;
[305]689      actual_delta = delta_corr;
[49]690
691      if (Jnlst().ProduceOutput(J_MOREVECTOR, J_MAIN)) {
[56]692        Jnlst().Printf(J_MOREVECTOR, J_MAIN,
693                       "*** Accepted corrector for Iteration: %d\n",
694                       IpData().iter_count());
[430]695        delta_corr->Print(Jnlst(), J_MOREVECTOR, J_MAIN, "delta_corr");
[49]696      }
697    }
698
[529]699    IpData().TimingStats().TryCorrector().End();
[49]700    return accept;
701  }
702
[542]703  char FilterLSAcceptor::UpdateForNextIteration(Number alpha_primal_test)
[11]704  {
[542]705    char info_alpha_primal_char;
706    // Augment the filter if required
707    if (!IsFtype(alpha_primal_test) ||
708        !ArmijoHolds(alpha_primal_test)) {
709      AugmentFilter();
710      info_alpha_primal_char = 'h';
[11]711    }
[542]712    else {
713      info_alpha_primal_char = 'f';
714    }
715    return info_alpha_primal_char;
[11]716  }
717
[542]718  void FilterLSAcceptor::PrepareRestoPhaseStart()
[265]719  {
[542]720    AugmentFilter();
[265]721  }
722
[416]723
[2]724} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.