source: branches/dev/Algorithm/IpBacktrackingLineSearch.cpp @ 546

Last change on this file since 546 was 546, checked in by andreasw, 15 years ago
  • added missing svn:keyword for a few files
  • increased factor for theta reduction in restoration convergence check from 1e1 to 1e2
  • trigger restoration phase in expect_infeasible_problem mode if multipliers become larger than 1e8
  • revert to regular restoration phase if more than max_soft_resto_iters (new option) iterations are taken in soft restoration phase
  • made quality_function default for adaptive barrier strategy
  • added new print level to avoid one line of output per iteration
  • fixed TNLPAdapter, so that only some parts of the problem (objective, constraints, variables) are scaled if desired
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 39.3 KB
Line 
1// Copyright (C) 2004, 2005 International Business Machines and others.
2// All Rights Reserved.
3// This code is published under the Common Public License.
4//
5// $Id: IpBacktrackingLineSearch.cpp 546 2005-10-20 22:40:15Z andreasw $
6//
7// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
8//           Andreas Waechter                 IBM    2005-10-13
9//               derived file from IpFilterLineSearch.cpp
10
11#include "IpBacktrackingLineSearch.hpp"
12#include "IpJournalist.hpp"
13#include "IpRestoPhase.hpp"
14#include "IpAlgTypes.hpp"
15
16#ifdef HAVE_CMATH
17# include <cmath>
18#else
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
26#ifdef HAVE_CCTYPE
27# include <cctype>
28#else
29# ifdef HAVE_CTYPE_H
30#  include <ctype.h>
31# else
32#  error "don't have header file for ctype"
33# endif
34#endif
35
36namespace Ipopt
37{
38
39#ifdef IP_DEBUG
40  static const Index dbg_verbosity = 0;
41#endif
42
43  BacktrackingLineSearch::BacktrackingLineSearch(
44    const SmartPtr<BacktrackingLSAcceptor>& acceptor,
45    const SmartPtr<RestorationPhase>& resto_phase,
46    const SmartPtr<ConvergenceCheck>& conv_check)
47      :
48      LineSearch(),
49      acceptor_(acceptor),
50      resto_phase_(resto_phase),
51      conv_check_(conv_check)
52  {
53    DBG_START_FUN("BacktrackingLineSearch::BacktrackingLineSearch",
54                  dbg_verbosity);
55    DBG_ASSERT(IsValid(acceptor_));
56  }
57
58  BacktrackingLineSearch::~BacktrackingLineSearch()
59  {
60    DBG_START_FUN("BacktrackingLineSearch::~BacktrackingLineSearch()",
61                  dbg_verbosity);
62  }
63
64  void BacktrackingLineSearch::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
65  {
66    roptions->AddBoundedNumberOption(
67      "alpha_red_factor",
68      "Fractional reduction of the trial step size in the backtracking line search.",
69      0.0, true, 1.0, true, 0.5,
70      "At every step of the backtracking line search, the trial step size is "
71      "reduced by this factor.");
72
73    std::string prev_category = roptions->RegisteringCategory();
74    roptions->SetRegisteringCategory("Undocumented");
75    roptions->AddStringOption2(
76      "magic_steps",
77      "Enables magic steps.",
78      "no",
79      "no", "don't take magic steps",
80      "yes", "take magic steps",
81      "DOESN'T REALLY WORK YET!");
82    roptions->SetRegisteringCategory(prev_category);
83
84    roptions->AddStringOption2(
85      "accept_every_trial_step",
86      "Always accept the frist trial step.",
87      "no",
88      "no", "don't arbitrarily accept the full step",
89      "yes", "always accept the full step",
90      "Setting this option to \"yes\" essentially disables the line search "
91      "and makes the algorithm take aggressive steps.");
92
93    roptions->AddStringOption7(
94      "alpha_for_y",
95      "Method to determine the step size for constraint multipliers.",
96      "primal",
97      "primal", "use primal step size",
98      "bound_mult", "use step size for the bound multipliers",
99      "min", "use the min of primal and bound multipliers",
100      "max", "use the max of primal and bound multipliers",
101      "full", "take a full step of size one",
102      "min_dual_infeas", "choose step size minimizing new dual infeasibility",
103      "safe_min_dual_infeas", "like \"min_dual_infeas\", but safeguarded by \"min\" and \"max\"",
104      "This option determines how the step size (alpha_y) will be calculated when updating the "
105      "constraint multipliers.");
106
107    roptions->AddStringOption2(
108      "expect_infeasible_problem",
109      "Enable heuristics to quickly detect an infeasible problem.",
110      "no",
111      "no", "the problem probably be feasible",
112      "yes", "the problem has a good chance to be infeasible",
113      "This options is meant to activate heuristics that may speed up the "
114      "infeasibility determination if you expect that there is a good chance for the problem to be "
115      "infeasible.  In the filter line search procedure, the restoration "
116      "phase is called more quickly than usually, and more reduction in "
117      "the constraint violation is enforced. If the problem is square, this "
118      "option is enabled automatically.");
119    roptions->AddLowerBoundedNumberOption(
120      "expect_infeasible_problem_ctol",
121      "Threshold for disabling \"expect_infeasible_problem\" option.",
122      0.0, false, 1e-3,
123      "If the constraint violation becomes smaller than this threshold, "
124      "the \"expect_infeasible_problem\" heuristics in the filter line "
125      "search are disabled. If the problem is square, this options is set to "
126      "0.");
127    roptions->AddStringOption2(
128      "start_with_resto",
129      "Tells algorithm to switch to restoration phase in first iteration.",
130      "no",
131      "no", "don't force start in restoration phase",
132      "yes", "force start in restoration phase",
133      "Setting this option to \"yes\" forces the algorithm to switch to the "
134      "feasibility restoration phase in the first iteration. If the initial "
135      "point is feasible, the algorithm will abort with a failure.");
136    roptions->AddLowerBoundedNumberOption(
137      "tiny_step_tol",
138      "Tolerance for detecting numerically insignificant steps.",
139      0.0, false, 10.0*std::numeric_limits<double>::epsilon(),
140      "If the search direction in the primal variables (x and s) is, in "
141      "relative terms for each component, less than this value, the "
142      "algorithm accepts the full step without line search.  If this happens "
143      "repeatedly, the algorithm will terminate with a corresponding exit "
144      "message. The default value is 10 times machine precision.");
145    roptions->AddLowerBoundedIntegerOption(
146      "watchdog_shortened_iter_trigger",
147      "Number of shortened iterations that trigger the watchdog.",
148      0, 10,
149      "If the number of iterations in which the backtracking line search "
150      "did not accept the first trial point exceedes this number, the "
151      "watchdog procedure is activated.  Choosing \"0\" here disables the "
152      "watchdog procedure.");
153    roptions->AddLowerBoundedIntegerOption(
154      "watchdog_trial_iter_max",
155      "Maximum number of watchdog iterations.",
156      1, 3,
157      "This option determines the number of trial iterations "
158      "allowed before the watchdog "
159      "procedure is aborted and the algorithm returns to the stored point.");
160
161    roptions->AddLowerBoundedNumberOption(
162      "soft_resto_pderror_reduction_factor",
163      "Required reduction in primal-dual error in the soft restoration phase.",
164      0.0, false, (1.0 - 1e-4),
165      "The soft restoration phase attempts to reduce the "
166      "primal-dual error with regular steps. If the damped "
167      "primal-dual step (damped only to satisfy the "
168      "fraction-to-the-boundary rule) is not decreasing the primal-dual error "
169      "by at least this factor, then the regular restoration phase is called. "
170      "Choosing \"0\" here disables the soft "
171      "restoration phase.");
172    roptions->AddLowerBoundedIntegerOption(
173      "max_soft_resto_iters",
174      "Maximum number of iterations performed successively in soft restoration phase.",
175      0, 10,
176      "If the soft restoration phase is performed for more than so many "
177      "iteratins in a row, the regular restoration phase is called.");
178  }
179
180  bool BacktrackingLineSearch::InitializeImpl(const OptionsList& options,
181      const std::string& prefix)
182  {
183    options.GetNumericValue("alpha_red_factor", alpha_red_factor_, prefix);
184    options.GetBoolValue("magic_steps", magic_steps_, prefix);
185    options.GetBoolValue("accept_every_trial_step", accept_every_trial_step_, prefix);
186    Index enum_int;
187    options.GetEnumValue("alpha_for_y", enum_int, prefix);
188    alpha_for_y_ = AlphaForYEnum(enum_int);
189    options.GetNumericValue("expect_infeasible_problem_ctol", expect_infeasible_problem_ctol_, prefix);
190    options.GetBoolValue("expect_infeasible_problem", expect_infeasible_problem_, prefix);
191
192    options.GetBoolValue("start_with_resto", start_with_resto_, prefix);
193
194    options.GetNumericValue("tiny_step_tol", tiny_step_tol_, prefix);
195    options.GetIntegerValue("watchdog_trial_iter_max", watchdog_trial_iter_max_, prefix);
196    options.GetIntegerValue("watchdog_shortened_iter_trigger", watchdog_shortened_iter_trigger_, prefix);
197    options.GetNumericValue("soft_resto_pderror_reduction_factor",
198                            soft_resto_pderror_reduction_factor_, prefix);
199    options.GetIntegerValue("max_soft_resto_iters", max_soft_resto_iters_,
200                            prefix);
201
202    bool retvalue = true;
203    if (IsValid(resto_phase_)) {
204      if (!resto_phase_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
205                                    options, prefix)) {
206        return false;
207      }
208      if (!acceptor_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
209                                 options, prefix)) {
210        return false;
211      }
212    }
213
214    rigorous_ = true;
215    skipped_line_search_ = false;
216    tiny_step_last_iteration_ = false;
217
218    Reset();
219
220    count_successive_shortened_steps_ = 0;
221
222    acceptable_iterate_ = NULL;
223
224    return retvalue;
225  }
226
227  void BacktrackingLineSearch::FindAcceptableTrialPoint()
228  {
229    DBG_START_METH("BacktrackingLineSearch::FindAcceptableTrialPoint",
230                   dbg_verbosity);
231    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
232                   "--> Starting filter line search in iteration %d <--\n",
233                   IpData().iter_count());
234
235    // If the problem is square, we want to enable the
236    // expect_infeasible_problem option automatically so that the
237    // restoration phase is entered soon
238    if (IpCq().IsSquareProblem()) {
239      expect_infeasible_problem_ = true;
240      expect_infeasible_problem_ctol_ = 0.;
241    }
242
243    // Store current iterate if the optimality error is on acceptable
244    // level to restored if things fail later
245    if (CurrentIsAcceptable()) {
246      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
247                     "Storing current iterate as backup acceptable point.\n");
248      StoreAcceptablePoint();
249    }
250
251    // First assume that line search will find an acceptable trial point
252    skipped_line_search_ = false;
253
254    // Initialize the acceptor for this backtracking line search
255    acceptor_->InitThisLineSearch(in_watchdog_);
256
257    // Get the search directions (this will store the actual search
258    // direction, possibly including higher order corrections)
259    SmartPtr<IteratesVector> actual_delta =
260      IpData().delta()->MakeNewContainer();
261
262    bool goto_resto = false;
263    if (actual_delta->Asum() == 0. ) {
264      // In this case, a search direction could not be computed, and
265      // we should immediately go to the restoration phase.  ToDo: Cue
266      // off of a something else than the norm of the search direction
267      goto_resto = true;
268    }
269
270    if (start_with_resto_) {
271      // If the user requested to start with the restoration phase,
272      // skip the line search and do exactly that.  Reset the flag so
273      // that this happens only once.
274      goto_resto = true;
275      start_with_resto_= false;
276    }
277
278    if (expect_infeasible_problem_ &&
279        Max(IpData().curr()->y_c()->Amax(),
280            IpData().curr()->y_d()->Amax()) > 1e8) {
281      goto_resto = true;
282    }
283
284    bool accept = false;
285    bool corr_taken = false;
286    bool soc_taken = false;
287    Index n_steps = 0;
288    Number alpha_primal = 0.;
289
290    // Check if search direction becomes too small
291    bool tiny_step = (!goto_resto && DetectTinyStep());
292
293    if (in_watchdog_ && (goto_resto || tiny_step)) {
294      // If the step could not be computed or is too small and the
295      // watchdog is active, stop the watch dog and resume everything
296      // from reference point
297      StopWatchDog(actual_delta);
298      goto_resto = false;
299      tiny_step = false;
300    }
301
302    // Check if we want to wake up the watchdog
303    if (watchdog_shortened_iter_trigger_ > 0 &&
304        !in_watchdog_ && !goto_resto && !tiny_step &&
305        !in_soft_resto_phase_ && !expect_infeasible_problem_ &&
306        watchdog_shortened_iter_ >= watchdog_shortened_iter_trigger_) {
307      StartWatchDog();
308    }
309
310    // Handle the situation of a tiny step
311    if (tiny_step) {
312      alpha_primal =
313        IpCq().curr_primal_frac_to_the_bound(IpData().curr_tau());
314      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
315                     "Tiny step detected. Use step size alpha = %e unchecked\n",
316                     alpha_primal);
317      IpData().SetTrialPrimalVariablesFromStep(alpha_primal, *IpData().delta()->x(), *IpData().delta()->s());
318      IpData().Set_info_ls_count(0);
319
320      if (tiny_step_last_iteration_) {
321        IpData().Set_info_alpha_primal_char('T');
322        IpData().Set_tiny_step_flag(true);
323      }
324      else {
325        IpData().Set_info_alpha_primal_char('t');
326      }
327
328      tiny_step_last_iteration_ = true;
329      accept = true;
330    }
331    else {
332      tiny_step_last_iteration_ = false;
333    }
334
335    if (!goto_resto && !tiny_step) {
336
337      if (in_soft_resto_phase_) {
338        soft_resto_counter_++;
339        if (soft_resto_counter_ > max_soft_resto_iters_) {
340          accept = false;
341        }
342        else {
343          // If we are currently in the soft restoration phase, continue
344          // that way, and switch back if enough progress is made to the
345          // original criterion (e.g., the filter)
346          bool satisfies_original_criterion = false;
347          // ToDo use tiny_step in TrySoftRestoStep?
348          accept = TrySoftRestoStep(actual_delta,
349                                    satisfies_original_criterion);
350          if (accept) {
351            IpData().Set_info_alpha_primal_char('s');
352            if (satisfies_original_criterion) {
353              in_soft_resto_phase_ = false;
354              soft_resto_counter_ = 0;
355              IpData().Set_info_alpha_primal_char('S');
356            }
357          }
358        }
359      }
360      else {
361        // Start the backtracking line search
362        bool done = false;
363        bool skip_first_trial_point = false;
364        bool evaluation_error;
365        while (!done) {
366          accept = DoBacktrackingLineSearch(skip_first_trial_point,
367                                            alpha_primal,
368                                            corr_taken,
369                                            soc_taken,
370                                            n_steps,
371                                            evaluation_error,
372                                            actual_delta);
373          DBG_PRINT((1, "evaluation_error = %d\n", evaluation_error));
374          if (in_watchdog_) {
375            if (accept) {
376              in_watchdog_ = false;
377              IpData().Append_info_string("W");
378              Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
379                             "Watch dog procedure successful!\n");
380              done = true;
381            }
382            else {
383              watchdog_trial_iter_++;
384              if (evaluation_error ||
385                  watchdog_trial_iter_ > watchdog_trial_iter_max_) {
386                StopWatchDog(actual_delta);
387                skip_first_trial_point = true;
388              }
389              else {
390                done = true;
391                accept = true;
392              }
393            }
394          }
395          else {
396            done = true;
397          }
398        }
399      } /* else: if (in_soft_resto_phase_) { */
400    } /* if (!goto_resto && !tiny_step) { */
401
402    // If line search has been aborted because the step size becomes
403    // too small, go to the restoration phase or continue with soft
404    // restoration phase
405    if (!accept) {
406      // If we are not asked to do a rigorous line search, do no call
407      // the restoration phase.
408      if (!rigorous_) {
409        Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Skipping call of restoration phase...\n");
410        skipped_line_search_=true;
411      }
412      else {
413        // Check if we should start the soft restoration phase
414        if (!in_soft_resto_phase_
415            && !goto_resto && !expect_infeasible_problem_) {
416          Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
417                         "--> Starting soft restoration phase <--\n");
418          // Prepare the restoration phase, e.g., augment the filter
419          // with the current point.
420          acceptor_->PrepareRestoPhaseStart();
421
422          // Try the current search direction for the soft restoration phase
423          bool satisfies_original_criterion;
424          accept = TrySoftRestoStep(actual_delta,
425                                    satisfies_original_criterion);
426          // If it has been accepted: If the original criterion is also
427          // satisfied, we can just take that step and continue with
428          // the regular algorithm, otherwise we stay in the soft
429          // restoration phase
430          if (accept) {
431            if (satisfies_original_criterion) {
432              IpData().Set_info_alpha_primal_char('S');
433            }
434            else {
435              in_soft_resto_phase_ = true;
436              IpData().Set_info_alpha_primal_char('s');
437            }
438          }
439        }
440
441        if (!accept) {
442          // Go to the restoration phase
443          if (!in_soft_resto_phase_) {
444            // Prepare the restoration phase, e.g., augment the filter
445            // with the current point. If we are already in the soft
446            // restoration phase, this has been done earlier
447            acceptor_->PrepareRestoPhaseStart();
448          }
449          if (CurrentIsAcceptable()) {
450            THROW_EXCEPTION(ACCEPTABLE_POINT_REACHED,
451                            "Restoration phase called at acceptable point.");
452          }
453
454          if (!IsValid(resto_phase_)) {
455            THROW_EXCEPTION(IpoptException, "No Restoration Phase given to this Filter Line Search Object!");
456          }
457          // ToDo make the 1e-2 below a parameter?
458          if (IpCq().curr_constraint_violation()<=
459              1e-2*IpData().tol()) {
460            bool found_acceptable = RestoreAcceptablePoint();
461            if (found_acceptable) {
462              THROW_EXCEPTION(ACCEPTABLE_POINT_REACHED,
463                              "Restoration phase called at almost feasible point, but acceptable point could be restored.\n");
464            }
465            else {
466              // ToDo does that happen too often?
467              THROW_EXCEPTION(RESTORATION_FAILED,
468                              "Restoration phase called, but point is almost feasible.");
469            }
470          }
471
472          // Set the info fields for the first output line in the
473          // restoration phase which reflects why the restoration phase
474          // was called
475          IpData().Set_info_alpha_primal(alpha_primal);
476          IpData().Set_info_alpha_dual(0.);
477          IpData().Set_info_alpha_primal_char('R');
478          IpData().Set_info_ls_count(n_steps+1);
479
480          accept = resto_phase_->PerformRestoration();
481          if (!accept) {
482            bool found_acceptable = RestoreAcceptablePoint();
483            if (found_acceptable) {
484              THROW_EXCEPTION(ACCEPTABLE_POINT_REACHED,
485                              "Restoration phase failed, but acceptable point could be restore.\n");
486            }
487            else {
488              THROW_EXCEPTION(RESTORATION_FAILED,
489                              "Failed restoration phase!!!");
490            }
491          }
492          count_successive_shortened_steps_ = 0;
493          if (expect_infeasible_problem_) {
494            expect_infeasible_problem_ = false;
495          }
496          in_soft_resto_phase_ = false;
497          soft_resto_counter_ = 0;
498          watchdog_shortened_iter_ = 0;
499        }
500      }
501    }
502    else if (!in_soft_resto_phase_ || tiny_step) {
503      // we didn't do the restoration phase and are now updating the
504      // dual variables of the trial point
505      Number alpha_dual_max =
506        IpCq().dual_frac_to_the_bound(IpData().curr_tau(),
507                                      *actual_delta->z_L(), *actual_delta->z_U(),
508                                      *actual_delta->v_L(), *actual_delta->v_U());
509
510      PerformDualStep(alpha_primal, alpha_dual_max, actual_delta);
511
512      if (n_steps==0) {
513        // accepted this if a full step was
514        // taken
515        count_successive_shortened_steps_ = 0;
516        watchdog_shortened_iter_ = 0;
517      }
518      else {
519        count_successive_shortened_steps_++;
520        watchdog_shortened_iter_++;
521      }
522
523      if (expect_infeasible_problem_ &&
524          IpCq().curr_constraint_violation() <= expect_infeasible_problem_ctol_) {
525        Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
526                       "Constraint violation is with %e less than expect_infeasible_problem_ctol.\nDisable expect_infeasible_problem_heuristic.\n", IpCq().curr_constraint_violation());
527        expect_infeasible_problem_ = false;
528      }
529    }
530  }
531
532  bool BacktrackingLineSearch::DoBacktrackingLineSearch(bool skip_first_trial_point,
533      Number& alpha_primal,
534      bool& corr_taken,
535      bool& soc_taken,
536      Index& n_steps,
537      bool& evaluation_error,
538      SmartPtr<IteratesVector>& actual_delta)
539  {
540    evaluation_error = false;
541    bool accept = false;
542
543    DBG_START_METH("BacktrackingLineSearch::DoBacktrackingLineSearch",
544                   dbg_verbosity);
545
546    // Compute primal fraction-to-the-boundary value
547    Number alpha_primal_max =
548      IpCq().primal_frac_to_the_bound(IpData().curr_tau(),
549                                      *actual_delta->x(),
550                                      *actual_delta->s());
551
552    // Compute smallest step size allowed
553    Number alpha_min = alpha_primal_max;
554    if (!in_watchdog_) {
555      alpha_min = acceptor_->CalculateAlphaMin();
556    }
557    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
558                   "minimal step size ALPHA_MIN = %E\n", alpha_min);
559
560    // Start line search from maximal step size
561    alpha_primal = alpha_primal_max;
562
563    // Step size used in ftype and armijo tests
564    Number alpha_primal_test = alpha_primal;
565    if (in_watchdog_) {
566      alpha_primal_test = watchdog_alpha_primal_test_;
567    }
568
569    if (skip_first_trial_point) {
570      alpha_primal *= alpha_red_factor_;
571    }
572
573    if (!skip_first_trial_point) {
574      // Before we do the actual backtracking line search for the
575      // regular primal-dual search direction, let's see if a step
576      // including a higher-order correctior is already acceptable
577      accept = acceptor_->TryCorrector(alpha_primal_test,
578                                       alpha_primal,
579                                       actual_delta);
580    }
581    if (accept) {
582      corr_taken = true;
583    }
584
585    if (!accept) {
586      // Loop over decreaseing step sizes until acceptable point is
587      // found or until step size becomes too small
588
589      while (alpha_primal>alpha_min ||
590             n_steps == 0) { // always allow the "full" step if it is
591        // acceptable (even if alpha_primal<=alpha_min)
592        Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
593                       "Starting checks for alpha (primal) = %8.2e\n",
594                       alpha_primal);
595
596        try {
597          // Compute the primal trial point
598          IpData().SetTrialPrimalVariablesFromStep(alpha_primal, *actual_delta->x(), *actual_delta->s());
599
600          if (magic_steps_) {
601            PerformMagicStep();
602          }
603
604          // If it is acceptable, stop the search
605          alpha_primal_test = alpha_primal;
606          if (accept_every_trial_step_) {
607            // We call the evaluation at the trial point here, so that an
608            // exception will the thrown if there are problem during the
609            // evaluation of the functions (in that case, we want to further
610            // reduce the step size
611            IpCq().trial_barrier_obj();
612            IpCq().trial_constraint_violation();
613            accept = true;
614          }
615          else {
616            accept = acceptor_->CheckAcceptabilityOfTrialPoint(alpha_primal_test);
617          }
618        }
619        catch(IpoptNLP::Eval_Error& e) {
620          e.ReportException(Jnlst(), J_DETAILED);
621          Jnlst().Printf(J_WARNING, J_LINE_SEARCH,
622                         "Warning: Cutting back alpha due to evaluation error\n");
623          IpData().Append_info_string("e");
624          accept = false;
625          evaluation_error = true;
626        }
627
628        if (accept) {
629          break;
630        }
631
632        if (in_watchdog_) {
633          break;
634        }
635
636        // Decide if we want to go to the restoration phase in a
637        // short cut to check if the problem is infeasible
638        if (expect_infeasible_problem_) {
639          if (count_successive_shortened_steps_>=5) {
640            break;
641          }
642        }
643
644        // try second order correction step if the function could
645        // be evaluated
646        // DoTo: check if we want to do SOC when watchdog is active
647        if (!evaluation_error) {
648          Number theta_curr = IpCq().curr_constraint_violation();
649          Number theta_trial = IpCq().trial_constraint_violation();
650          if (alpha_primal==alpha_primal_max &&       // i.e. first trial point
651              theta_curr<=theta_trial) {
652            // Try second order correction
653            accept = acceptor_->TrySecondOrderCorrection(alpha_primal_test,
654                     alpha_primal,
655                     actual_delta);
656          }
657          if (accept) {
658            soc_taken = true;
659            break;
660          }
661        }
662
663        // Point is not yet acceptable, try a shorter one
664        alpha_primal *= alpha_red_factor_;
665        n_steps++;
666      }
667    } /* if (!accept) */
668
669    char info_alpha_primal_char='?';
670    if (!accept && in_watchdog_) {
671      info_alpha_primal_char = 'w';
672    }
673    else if (accept) {
674      info_alpha_primal_char =
675        acceptor_->UpdateForNextIteration(alpha_primal_test);
676    }
677    if (soc_taken) {
678      info_alpha_primal_char = toupper(info_alpha_primal_char);
679    }
680    IpData().Set_info_alpha_primal_char(info_alpha_primal_char);
681    IpData().Set_info_ls_count(n_steps+1);
682    if (corr_taken) {
683      IpData().Append_info_string("C");
684    }
685
686    return accept;
687  }
688
689  void BacktrackingLineSearch::StartWatchDog()
690  {
691    DBG_START_FUN("BacktrackingLineSearch::StartWatchDog", dbg_verbosity);
692
693    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
694                   "Starting Watch Dog\n");
695
696    in_watchdog_ = true;
697    watchdog_iterate_ = IpData().curr();
698    watchdog_delta_ = IpData().delta();
699    watchdog_trial_iter_ = 0;
700    watchdog_alpha_primal_test_ =
701      IpCq().curr_primal_frac_to_the_bound(IpData().curr_tau());
702
703    acceptor_->StartWatchDog();
704  }
705
706  void BacktrackingLineSearch::StopWatchDog(SmartPtr<IteratesVector>& actual_delta)
707  {
708    DBG_START_FUN("BacktrackingLineSearch::StopWatchDog", dbg_verbosity);
709
710    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
711                   "Stopping Watch Dog\n");
712
713    IpData().Append_info_string("w");
714
715    in_watchdog_ = false;
716
717    // Reset all fields in IpData to reference point
718    SmartPtr<IteratesVector> old_trial = watchdog_iterate_->MakeNewContainer();
719    IpData().set_trial(old_trial);
720    IpData().AcceptTrialPoint();
721    actual_delta = watchdog_delta_->MakeNewContainer();
722    IpData().SetHaveAffineDeltas(false);
723
724    // reset the stored watchdog iterates
725    watchdog_iterate_ = NULL;
726    watchdog_delta_ = NULL;
727
728    watchdog_shortened_iter_ = 0;
729
730    acceptor_->StopWatchDog();
731  }
732
733  void BacktrackingLineSearch::Reset()
734  {
735    DBG_START_FUN("BacktrackingLineSearch::Reset", dbg_verbosity);
736    in_soft_resto_phase_ = false;
737    soft_resto_counter_ = 0;
738
739    // Inactivate the watchdog and release all stored data
740    in_watchdog_ = false;
741    watchdog_iterate_ = NULL;
742    watchdog_delta_ = NULL;
743    watchdog_shortened_iter_ = 0;
744
745    acceptor_->Reset();
746  }
747
748  void BacktrackingLineSearch::PerformDualStep(Number alpha_primal,
749      Number alpha_dual,
750      SmartPtr<IteratesVector>& delta)
751  {
752    DBG_START_FUN("BacktrackingLineSearch::PerformDualStep", dbg_verbosity);
753
754    // set the bound multipliers from the step
755    IpData().SetTrialBoundMultipliersFromStep(alpha_dual, *delta->z_L(), *delta->z_U(), *delta->v_L(), *delta->v_U());
756
757    Number alpha_y=-1.;
758    switch (alpha_for_y_) {
759      case PRIMAL_ALPHA_FOR_Y:
760      alpha_y = alpha_primal;
761      break;
762      case DUAL_ALPHA_FOR_Y:
763      alpha_y = alpha_dual;
764      break;
765      case MIN_ALPHA_FOR_Y:
766      alpha_y = Min(alpha_dual, alpha_primal);
767      break;
768      case MAX_ALPHA_FOR_Y:
769      alpha_y = Max(alpha_dual, alpha_primal);
770      break;
771      case FULL_STEP_FOR_Y:
772      alpha_y = 1;
773      break;
774      case MIN_DUAL_INFEAS_ALPHA_FOR_Y:
775      case SAFE_MIN_DUAL_INFEAS_ALPHA_FOR_Y:
776      // Here we compute the step size for y so that the dual
777      // infeasibility is minimized along delta_y
778
779      // compute the dual infeasibility at new point with old y
780      SmartPtr<IteratesVector> temp_trial
781      = IpData().trial()->MakeNewContainer();
782      temp_trial->Set_y_c(*IpData().curr()->y_c());
783      temp_trial->Set_y_d(*IpData().curr()->y_d());
784      IpData().set_trial(temp_trial);
785      SmartPtr<const Vector> dual_inf_x = IpCq().trial_grad_lag_x();
786      SmartPtr<const Vector> dual_inf_s = IpCq().trial_grad_lag_s();
787
788      SmartPtr<Vector> new_jac_times_delta_y =
789        IpData().curr()->x()->MakeNew();
790      new_jac_times_delta_y->AddTwoVectors(1., *IpCq().trial_jac_cT_times_vec(*delta->y_c()),
791                                           1., *IpCq().trial_jac_dT_times_vec(*delta->y_d()),
792                                           0.);
793
794      Number a = pow(new_jac_times_delta_y->Nrm2(), 2.) +
795                 pow(delta->y_d()->Nrm2(), 2.);
796      Number b = dual_inf_x->Dot(*new_jac_times_delta_y) -
797                 dual_inf_s->Dot(*delta->y_d());
798
799      Number alpha = - b/a;
800
801      if (alpha_for_y_==SAFE_MIN_DUAL_INFEAS_ALPHA_FOR_Y) {
802        alpha_y = Min(Max(alpha_primal, alpha_dual), Max(alpha, Min(alpha_primal, alpha_dual)));
803      }
804      else {
805        alpha_y = Min(1., Max(0., alpha));
806      }
807      break;
808    }
809
810    // Set the eq multipliers from the step now that alpha_y
811    // has been calculated.
812    DBG_PRINT((1, "alpha_y = %e\n", alpha_y));
813    DBG_PRINT_VECTOR(2, "delta_y_c", *delta->y_c());
814    DBG_PRINT_VECTOR(2, "delta_y_d", *delta->y_d());
815    IpData().SetTrialEqMultipliersFromStep(alpha_y, *delta->y_c(), *delta->y_d());
816
817    // Set some information for iteration summary output
818    IpData().Set_info_alpha_primal(alpha_primal);
819    IpData().Set_info_alpha_dual(alpha_dual);
820  }
821
822  void
823  BacktrackingLineSearch::PerformMagicStep()
824  {
825    DBG_START_METH("BacktrackingLineSearch::PerformMagicStep",
826                   2);//dbg_verbosity);
827
828    DBG_PRINT((1,"Incoming barr = %e and constrviol %e\n",
829               IpCq().trial_barrier_obj(),
830               IpCq().trial_constraint_violation()));
831    DBG_PRINT_VECTOR(2, "s in", *IpData().trial()->s());
832    DBG_PRINT_VECTOR(2, "d minus s in", *IpCq().trial_d_minus_s());
833    DBG_PRINT_VECTOR(2, "slack_s_L in", *IpCq().trial_slack_s_L());
834    DBG_PRINT_VECTOR(2, "slack_s_U in", *IpCq().trial_slack_s_U());
835
836    SmartPtr<const Vector> d_L = IpNLP().d_L();
837    SmartPtr<const Matrix> Pd_L = IpNLP().Pd_L();
838    SmartPtr<Vector> delta_s_magic_L = d_L->MakeNew();
839    delta_s_magic_L->Set(0.);
840    SmartPtr<Vector> tmp = d_L->MakeNew();
841    Pd_L->TransMultVector(1., *IpCq().trial_d_minus_s(), 0., *tmp);
842    delta_s_magic_L->ElementWiseMax(*tmp);
843
844    SmartPtr<const Vector> d_U = IpNLP().d_U();
845    SmartPtr<const Matrix> Pd_U = IpNLP().Pd_U();
846    SmartPtr<Vector> delta_s_magic_U = d_U->MakeNew();
847    delta_s_magic_U->Set(0.);
848    tmp = d_U->MakeNew();
849    Pd_U->TransMultVector(1., *IpCq().trial_d_minus_s(), 0., *tmp);
850    delta_s_magic_U->ElementWiseMin(*tmp);
851
852    SmartPtr<Vector> delta_s_magic = IpData().trial()->s()->MakeNew();
853    Pd_L->MultVector(1., *delta_s_magic_L, 0., *delta_s_magic);
854    Pd_U->MultVector(1., *delta_s_magic_U, 1., *delta_s_magic);
855    delta_s_magic_L = NULL; // free memory
856    delta_s_magic_U = NULL; // free memory
857
858    // Now find those entries with both lower and upper bounds, there
859    // the step is too large
860    // ToDo this should only be done if there are inequality
861    // constraints with two bounds
862    // also this can be done in a smaller space (d_L or d_U whichever
863    // is smaller)
864    tmp = delta_s_magic->MakeNew();
865    tmp->Copy(*IpData().trial()->s());
866    Pd_L->MultVector(1., *d_L, -2., *tmp);
867    Pd_U->MultVector(1., *d_U, 1., *tmp);
868    SmartPtr<Vector> tmp2 = tmp->MakeNew();
869    tmp2->Copy(*tmp);
870    tmp2->ElementWiseAbs();
871    tmp->Axpy(-2., *delta_s_magic);
872    tmp->ElementWiseAbs();
873    // now, tmp2 = |d_L + d_u - 2*s| and tmp = |d_L + d_u - 2*(s+Delta s)|
874    // we want to throw out those for which tmp2 > tmp
875    tmp->Axpy(-1., *tmp2);
876    tmp->ElementWiseSgn();
877    tmp2->Set(0.);
878    tmp2->ElementWiseMax(*tmp);
879    tmp = d_L->MakeNew();
880    Pd_L->TransMultVector(1., *tmp2, 0., *tmp);
881    Pd_L->MultVector(1., *tmp, 0., *tmp2);
882    tmp = d_U->MakeNew();
883    Pd_U->TransMultVector(1., *tmp2, 0., *tmp);
884    Pd_U->MultVector(1., *tmp, 0., *tmp2);
885    DBG_PRINT_VECTOR(2, "tmp indicator", *tmp2)
886    // tmp2 now is one for those entries with both bounds, for which
887    // no step should be taken
888
889    tmp = delta_s_magic->MakeNew();
890    tmp->Copy(*delta_s_magic);
891    tmp->ElementWiseMultiply(*tmp2);
892    delta_s_magic->Axpy(-1., *tmp);
893
894    Number delta_s_magic_max = delta_s_magic->Amax();
895    Number mach_eps = std::numeric_limits<Number>::epsilon();
896    if (delta_s_magic_max>0.) {
897      if (delta_s_magic_max > 10*mach_eps*IpData().trial()->s()->Amax()) {
898        IpData().Append_info_string("M");
899        Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Magic step with max-norm %.6e taken.\n", delta_s_magic->Amax());
900        delta_s_magic->Print(Jnlst(), J_MOREVECTOR, J_LINE_SEARCH,
901                             "delta_s_magic");
902      }
903
904      // now finally compute the new overall slacks
905      delta_s_magic->Axpy(1., *IpData().trial()->s());
906      SmartPtr<IteratesVector> trial = IpData().trial()->MakeNewContainer();
907      trial->Set_s(*delta_s_magic);
908
909      // also update the set in the dual variables
910
911
912      IpData().set_trial(trial);
913    }
914
915    DBG_PRINT((1,"Outgoing barr = %e and constrviol %e\n", IpCq().trial_barrier_obj(), IpCq().trial_constraint_violation()));
916    DBG_PRINT_VECTOR(2, "s out", *IpData().trial()->s());
917    DBG_PRINT_VECTOR(2, "d minus s out", *IpCq().trial_d_minus_s());
918    DBG_PRINT_VECTOR(2, "slack_s_L out", *IpCq().trial_slack_s_L());
919    DBG_PRINT_VECTOR(2, "slack_s_U out", *IpCq().trial_slack_s_U());
920  }
921
922  bool BacktrackingLineSearch::TrySoftRestoStep(SmartPtr<IteratesVector>& actual_delta,
923      bool &satisfies_original_criterion)
924  {
925    DBG_START_FUN("FilterLSAcceptor::TrySoftRestoStep", dbg_verbosity);
926
927    if (soft_resto_pderror_reduction_factor_==0.) {
928      return false;
929    }
930
931    satisfies_original_criterion = false;
932
933    // ToDo: Need to decide if we want to try a corrector step first
934
935    // Compute the maximal step sizes (we use identical step sizes for
936    // primal and dual variables
937    Number alpha_primal_max =
938      IpCq().primal_frac_to_the_bound(IpData().curr_tau(),
939                                      *actual_delta->x(),
940                                      *actual_delta->s());
941    Number alpha_dual_max =
942      IpCq().dual_frac_to_the_bound(IpData().curr_tau(),
943                                    *actual_delta->z_L(),
944                                    *actual_delta->z_U(),
945                                    *actual_delta->v_L(),
946                                    *actual_delta->v_U());
947    Number alpha_max =  Min(alpha_primal_max, alpha_dual_max);
948
949    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
950                   "Trying soft restoration phase step with step length %13.6e\n",
951                   alpha_max);
952
953    // Set the trial point
954    IpData().SetTrialPrimalVariablesFromStep(alpha_max, *actual_delta->x(), *actual_delta->s());
955    PerformDualStep(alpha_max, alpha_max,
956                    actual_delta);
957
958    // Check if that point is acceptable with respect to the current
959    // original filter
960
961    Number trial_barr;
962    Number trial_theta;
963    try {
964      trial_barr = IpCq().trial_barrier_obj();
965      trial_theta = IpCq().trial_constraint_violation();
966    }
967    catch(IpoptNLP::Eval_Error& e) {
968      e.ReportException(Jnlst(), J_DETAILED);
969      Jnlst().Printf(J_WARNING, J_LINE_SEARCH,
970                     "Warning: Evaluation error during soft restoration phase step.\n");
971      IpData().Append_info_string("e");
972      return false;
973    }
974    if (acceptor_->CheckAcceptabilityOfTrialPoint(0.)) {
975      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
976                     "  Trial step acceptable with respect to original backtracking globalization.\n");
977      satisfies_original_criterion = true;
978      return true;
979    }
980
981    // Evaluate the optimality error at the new point
982    Number mu = .0;
983    if (!IpData().FreeMuMode()) {
984      mu = IpData().curr_mu();
985    }
986    Number trial_pderror;
987    Number curr_pderror;
988    try {
989      trial_pderror = IpCq().trial_primal_dual_system_error(mu);
990      curr_pderror = IpCq().curr_primal_dual_system_error(mu);
991    }
992    catch(IpoptNLP::Eval_Error& e) {
993      e.ReportException(Jnlst(), J_DETAILED);
994      Jnlst().Printf(J_WARNING, J_LINE_SEARCH,
995                     "Warning: Evaluation error during soft restoration phase step.\n");
996      IpData().Append_info_string("e");
997      return false;
998    }
999
1000    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
1001                   "  Primal-dual error at current point:  %23.16e\n", curr_pderror);
1002    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
1003                   "  Primal-dual error at trial point  :  %23.16e\n", trial_pderror);
1004    // Check if there is sufficient reduction in the optimality error
1005    if (trial_pderror <= soft_resto_pderror_reduction_factor_*curr_pderror) {
1006      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
1007                     "  Trial step accepted.\n");
1008      return true;
1009    }
1010
1011    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
1012                   "  Trial step rejected.\n");
1013    return false;
1014  }
1015
1016  bool
1017  BacktrackingLineSearch::DetectTinyStep()
1018  {
1019    DBG_START_METH("BacktrackingLineSearch::DetectTinyStep",
1020                   dbg_verbosity);
1021
1022    Number max_step_x;
1023    Number max_step_s;
1024
1025    if (tiny_step_tol_==0.)
1026      return false;
1027
1028    // ToDo try to find more efficient implementation
1029
1030    SmartPtr<Vector> tmp = IpData().curr()->x()->MakeNew();
1031    tmp->Copy(*IpData().curr()->x());
1032    tmp->ElementWiseAbs();
1033    tmp->AddScalar(1.);
1034
1035    SmartPtr<Vector> tmp2 = IpData().curr()->x()->MakeNew();
1036    tmp2->Copy(*IpData().delta()->x());
1037    tmp2->ElementWiseDivide(*tmp);
1038    max_step_x = tmp2->Amax();
1039    Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
1040                   "Relative step size for delta_x = %e\n",
1041                   max_step_x);
1042    if (max_step_x > tiny_step_tol_)
1043      return false;
1044
1045    tmp = IpData().curr()->s()->MakeNew();
1046    tmp->Copy(*IpData().curr()->s());
1047    tmp->ElementWiseAbs();
1048    tmp->AddScalar(1.);
1049
1050    tmp2 = IpData().curr()->s()->MakeNew();
1051    tmp2->Copy(*IpData().delta()->s());
1052    tmp2->ElementWiseDivide(*tmp);
1053    max_step_s = tmp2->Amax();
1054    Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
1055                   "Relative step size for delta_s = %e\n",
1056                   max_step_s);
1057    if (max_step_s > tiny_step_tol_)
1058      return false;
1059
1060    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
1061                   "Tiny step of relative size %e detected.\n",
1062                   Max(max_step_x, max_step_s));
1063
1064    return true;
1065  }
1066
1067  bool BacktrackingLineSearch::CurrentIsAcceptable()
1068  {
1069    return (IsValid(conv_check_) &&
1070            conv_check_->CurrentIsAcceptable());
1071  }
1072
1073  void BacktrackingLineSearch::StoreAcceptablePoint()
1074  {
1075    DBG_START_METH("BacktrackingLineSearch::StoreAcceptablePoint",
1076                   dbg_verbosity);
1077
1078    acceptable_iterate_ = IpData().curr();
1079  }
1080
1081  bool BacktrackingLineSearch::RestoreAcceptablePoint()
1082  {
1083    DBG_START_METH("BacktrackingLineSearch::RestoreAcceptablePoint",
1084                   dbg_verbosity);
1085
1086    if (!IsValid(acceptable_iterate_)) {
1087      return false;
1088    }
1089
1090    SmartPtr<IteratesVector> prev_iterate = acceptable_iterate_->MakeNewContainer();
1091    IpData().set_trial(prev_iterate);
1092    IpData().AcceptTrialPoint();
1093
1094    return true;
1095  }
1096
1097
1098} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.