source: branches/dev/Algorithm/IpBacktrackingLineSearch.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: 38.4 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 542 2005-10-13 22:43:08Z 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 qucikly 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  }
173
174  bool BacktrackingLineSearch::InitializeImpl(const OptionsList& options,
175      const std::string& prefix)
176  {
177    options.GetNumericValue("alpha_red_factor", alpha_red_factor_, prefix);
178    options.GetBoolValue("magic_steps", magic_steps_, prefix);
179    options.GetBoolValue("accept_every_trial_step", accept_every_trial_step_, prefix);
180    Index enum_int;
181    options.GetEnumValue("alpha_for_y", enum_int, prefix);
182    alpha_for_y_ = AlphaForYEnum(enum_int);
183    options.GetNumericValue("expect_infeasible_problem_ctol", expect_infeasible_problem_ctol_, prefix);
184    options.GetBoolValue("expect_infeasible_problem", expect_infeasible_problem_, prefix);
185
186    options.GetBoolValue("start_with_resto", start_with_resto_, prefix);
187
188    options.GetNumericValue("tiny_step_tol", tiny_step_tol_, prefix);
189    options.GetIntegerValue("watchdog_trial_iter_max", watchdog_trial_iter_max_, prefix);
190    options.GetIntegerValue("watchdog_shortened_iter_trigger", watchdog_shortened_iter_trigger_, prefix);
191    options.GetNumericValue("soft_resto_pderror_reduction_factor",
192                            soft_resto_pderror_reduction_factor_, prefix);
193
194    bool retvalue = true;
195    if (IsValid(resto_phase_)) {
196      if (!resto_phase_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
197                                    options, prefix)) {
198        return false;
199      }
200      if (!acceptor_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
201                                 options, prefix)) {
202        return false;
203      }
204    }
205
206    rigorous_ = true;
207    skipped_line_search_ = false;
208    tiny_step_last_iteration_ = false;
209
210    Reset();
211
212    count_successive_shortened_steps_ = 0;
213
214    acceptable_iterate_ = NULL;
215
216    return retvalue;
217  }
218
219  void BacktrackingLineSearch::FindAcceptableTrialPoint()
220  {
221    DBG_START_METH("BacktrackingLineSearch::FindAcceptableTrialPoint",
222                   dbg_verbosity);
223    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
224                   "--> Starting filter line search in iteration %d <--\n",
225                   IpData().iter_count());
226
227    // If the problem is square, we want to enable the
228    // expect_infeasible_problem option automatically so that the
229    // restoration phase is entered soon
230    if (IpCq().IsSquareProblem()) {
231      expect_infeasible_problem_ = true;
232      expect_infeasible_problem_ctol_ = 0.;
233    }
234
235    // Store current iterate if the optimality error is on acceptable
236    // level to restored if things fail later
237    if (CurrentIsAcceptable()) {
238      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
239                     "Storing current iterate as backup acceptable point.\n");
240      StoreAcceptablePoint();
241    }
242
243    // First assume that line search will find an acceptable trial point
244    skipped_line_search_ = false;
245
246    // Initialize the acceptor for this backtracking line search
247    acceptor_->InitThisLineSearch(in_watchdog_);
248
249    // Get the search directions (this will store the actual search
250    // direction, possibly including higher order corrections)
251    SmartPtr<IteratesVector> actual_delta =
252      IpData().delta()->MakeNewContainer();
253
254    bool goto_resto = false;
255    if (actual_delta->Asum() == 0. ) {
256      // In this case, a search direction could not be computed, and
257      // we should immediately go to the restoration phase.  ToDo: Cue
258      // off of a something else than the norm of the search direction
259      goto_resto = true;
260    }
261
262    if (start_with_resto_) {
263      // If the user requested to start with the restoration phase,
264      // skip the line search and do exactly that.  Reset the flag so
265      // that this happens only once.
266      goto_resto = true;
267      start_with_resto_= false;
268    }
269
270    bool accept = false;
271    bool corr_taken = false;
272    bool soc_taken = false;
273    Index n_steps = 0;
274    Number alpha_primal = 0.;
275
276    // Check if search direction becomes too small
277    bool tiny_step = (!goto_resto && DetectTinyStep());
278
279    if (in_watchdog_ && (goto_resto || tiny_step)) {
280      // If the step could not be computed or is too small and the
281      // watchdog is active, stop the watch dog and resume everything
282      // from reference point
283      StopWatchDog(actual_delta);
284      goto_resto = false;
285      tiny_step = false;
286    }
287
288    // Check if we want to wake up the watchdog
289    if (watchdog_shortened_iter_trigger_ > 0 &&
290        !in_watchdog_ && !goto_resto && !tiny_step &&
291        !in_soft_resto_phase_ && !expect_infeasible_problem_ &&
292        watchdog_shortened_iter_ >= watchdog_shortened_iter_trigger_) {
293      StartWatchDog();
294    }
295
296    // Handle the situation of a tiny step
297    if (tiny_step) {
298      alpha_primal =
299        IpCq().curr_primal_frac_to_the_bound(IpData().curr_tau());
300      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
301                     "Tiny step detected. Use step size alpha = %e unchecked\n",
302                     alpha_primal);
303      IpData().SetTrialPrimalVariablesFromStep(alpha_primal, *IpData().delta()->x(), *IpData().delta()->s());
304      IpData().Set_info_ls_count(0);
305
306      if (tiny_step_last_iteration_) {
307        IpData().Set_info_alpha_primal_char('T');
308        IpData().Set_tiny_step_flag(true);
309      }
310      else {
311        IpData().Set_info_alpha_primal_char('t');
312      }
313
314      tiny_step_last_iteration_ = true;
315      accept = true;
316    }
317    else {
318      tiny_step_last_iteration_ = false;
319    }
320
321    if (!goto_resto && !tiny_step) {
322
323      if (in_soft_resto_phase_) {
324        // If we are currently in the soft restoration phase, continue
325        // that way, and switch back if enough progress is made to the
326        // original criterion (e.g., the filter)
327        bool satisfies_original_criterion = false;
328        // ToDo use tiny_step in TrySoftRestoStep?
329        accept = TrySoftRestoStep(actual_delta,
330                                  satisfies_original_criterion);
331        if (accept) {
332          IpData().Set_info_alpha_primal_char('s');
333          if (satisfies_original_criterion) {
334            in_soft_resto_phase_ = false;
335            IpData().Set_info_alpha_primal_char('S');
336          }
337        }
338      }
339      else {
340        // Start the backtracking line search
341        bool done = false;
342        bool skip_first_trial_point = false;
343        bool evaluation_error;
344        while (!done) {
345          accept = DoBacktrackingLineSearch(skip_first_trial_point,
346                                            alpha_primal,
347                                            corr_taken,
348                                            soc_taken,
349                                            n_steps,
350                                            evaluation_error,
351                                            actual_delta);
352          DBG_PRINT((1, "evaluation_error = %d\n", evaluation_error));
353          if (in_watchdog_) {
354            if (accept) {
355              in_watchdog_ = false;
356              IpData().Append_info_string("W");
357              Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
358                             "Watch dog procedure successful!\n");
359              done = true;
360            }
361            else {
362              watchdog_trial_iter_++;
363              if (evaluation_error ||
364                  watchdog_trial_iter_ > watchdog_trial_iter_max_) {
365                StopWatchDog(actual_delta);
366                skip_first_trial_point = true;
367              }
368              else {
369                done = true;
370                accept = true;
371              }
372            }
373          }
374          else {
375            done = true;
376          }
377        }
378      } /* else: if (in_soft_resto_phase_) { */
379    } /* if (!goto_resto && !tiny_step) { */
380
381    // If line search has been aborted because the step size becomes
382    // too small, go to the restoration phase or continue with soft
383    // restoration phase
384    if (!accept) {
385      // If we are not asked to do a rigorous line search, do no call
386      // the restoration phase.
387      if (!rigorous_) {
388        Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Skipping call of restoration phase...\n");
389        skipped_line_search_=true;
390      }
391      else {
392        // Check if we should start the soft restoration phase
393        if (!in_soft_resto_phase_
394            && !goto_resto && !expect_infeasible_problem_) {
395          Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
396                         "--> Starting soft restoration phase <--\n");
397          // Prepare the restoration phase, e.g., augment the filter
398          // with the current point.
399          acceptor_->PrepareRestoPhaseStart();
400
401          // Try the current search direction for the soft restoration phase
402          bool satisfies_original_criterion;
403          accept = TrySoftRestoStep(actual_delta,
404                                    satisfies_original_criterion);
405          // If it has been accepted: If the original criterion is also
406          // satisfied, we can just take that step and continue with
407          // the regular algorithm, otherwise we stay in the soft
408          // restoration phase
409          if (accept) {
410            if (satisfies_original_criterion) {
411              IpData().Set_info_alpha_primal_char('S');
412            }
413            else {
414              in_soft_resto_phase_ = true;
415              IpData().Set_info_alpha_primal_char('s');
416            }
417          }
418        }
419
420        if (!accept) {
421          // Go to the restoration phase
422          if (!in_soft_resto_phase_) {
423            // Prepare the restoration phase, e.g., augment the filter
424            // with the current point. If we are already in the soft
425            // restoration phase, this has been done earlier
426            acceptor_->PrepareRestoPhaseStart();
427          }
428          if (CurrentIsAcceptable()) {
429            THROW_EXCEPTION(ACCEPTABLE_POINT_REACHED,
430                            "Restoration phase called at acceptable point.");
431          }
432
433          if (!IsValid(resto_phase_)) {
434            THROW_EXCEPTION(IpoptException, "No Restoration Phase given to this Filter Line Search Object!");
435          }
436          // ToDo make the 1e-2 below a parameter?
437          if (IpCq().curr_constraint_violation()<=
438              1e-2*IpData().tol()) {
439            bool found_acceptable = RestoreAcceptablePoint();
440            if (found_acceptable) {
441              THROW_EXCEPTION(ACCEPTABLE_POINT_REACHED,
442                              "Restoration phase called at almost feasible point, but acceptable point could be restored.\n");
443            }
444            else {
445              // ToDo does that happen too often?
446              THROW_EXCEPTION(RESTORATION_FAILED,
447                              "Restoration phase called, but point is almost feasible.");
448            }
449          }
450
451          // Set the info fields for the first output line in the
452          // restoration phase which reflects why the restoration phase
453          // was called
454          IpData().Set_info_alpha_primal(alpha_primal);
455          IpData().Set_info_alpha_dual(0.);
456          IpData().Set_info_alpha_primal_char('R');
457          IpData().Set_info_ls_count(n_steps+1);
458
459          accept = resto_phase_->PerformRestoration();
460          if (!accept) {
461            bool found_acceptable = RestoreAcceptablePoint();
462            if (found_acceptable) {
463              THROW_EXCEPTION(ACCEPTABLE_POINT_REACHED,
464                              "Restoration phase failed, but acceptable point could be restore.\n");
465            }
466            else {
467              THROW_EXCEPTION(RESTORATION_FAILED,
468                              "Failed restoration phase!!!");
469            }
470          }
471          count_successive_shortened_steps_ = 0;
472          if (expect_infeasible_problem_) {
473            expect_infeasible_problem_ = false;
474          }
475          in_soft_resto_phase_ = false;
476          watchdog_shortened_iter_ = 0;
477        }
478      }
479    }
480    else if (!in_soft_resto_phase_ || tiny_step) {
481      // we didn't do the restoration phase and are now updating the
482      // dual variables of the trial point
483      Number alpha_dual_max =
484        IpCq().dual_frac_to_the_bound(IpData().curr_tau(),
485                                      *actual_delta->z_L(), *actual_delta->z_U(),
486                                      *actual_delta->v_L(), *actual_delta->v_U());
487
488      PerformDualStep(alpha_primal, alpha_dual_max, actual_delta);
489
490      if (n_steps==0) {
491        // accepted this if a full step was
492        // taken
493        count_successive_shortened_steps_ = 0;
494        watchdog_shortened_iter_ = 0;
495      }
496      else {
497        count_successive_shortened_steps_++;
498        watchdog_shortened_iter_++;
499      }
500
501      if (expect_infeasible_problem_ &&
502          IpCq().curr_constraint_violation() <= expect_infeasible_problem_ctol_) {
503        Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
504                       "Constraint violation is with %e less than expect_infeasible_problem_ctol.\nDisable expect_infeasible_problem_heuristic.\n", IpCq().curr_constraint_violation());
505        expect_infeasible_problem_ = false;
506      }
507    }
508  }
509
510  bool BacktrackingLineSearch::DoBacktrackingLineSearch(bool skip_first_trial_point,
511      Number& alpha_primal,
512      bool& corr_taken,
513      bool& soc_taken,
514      Index& n_steps,
515      bool& evaluation_error,
516      SmartPtr<IteratesVector>& actual_delta)
517  {
518    evaluation_error = false;
519    bool accept = false;
520
521    DBG_START_METH("BacktrackingLineSearch::DoBacktrackingLineSearch",
522                   dbg_verbosity);
523
524    // Compute primal fraction-to-the-boundary value
525    Number alpha_primal_max =
526      IpCq().primal_frac_to_the_bound(IpData().curr_tau(),
527                                      *actual_delta->x(),
528                                      *actual_delta->s());
529
530    // Compute smallest step size allowed
531    Number alpha_min = alpha_primal_max;
532    if (!in_watchdog_) {
533      alpha_min = acceptor_->CalculateAlphaMin();
534    }
535    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
536                   "minimal step size ALPHA_MIN = %E\n", alpha_min);
537
538    // Start line search from maximal step size
539    alpha_primal = alpha_primal_max;
540
541    // Step size used in ftype and armijo tests
542    Number alpha_primal_test = alpha_primal;
543    if (in_watchdog_) {
544      alpha_primal_test = watchdog_alpha_primal_test_;
545    }
546
547    if (skip_first_trial_point) {
548      alpha_primal *= alpha_red_factor_;
549    }
550
551    if (!skip_first_trial_point) {
552      // Before we do the actual backtracking line search for the
553      // regular primal-dual search direction, let's see if a step
554      // including a higher-order correctior is already acceptable
555      accept = acceptor_->TryCorrector(alpha_primal_test,
556                                       alpha_primal,
557                                       actual_delta);
558    }
559    if (accept) {
560      corr_taken = true;
561    }
562
563    if (!accept) {
564      // Loop over decreaseing step sizes until acceptable point is
565      // found or until step size becomes too small
566
567      while (alpha_primal>alpha_min ||
568             n_steps == 0) { // always allow the "full" step if it is
569        // acceptable (even if alpha_primal<=alpha_min)
570        Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
571                       "Starting checks for alpha (primal) = %8.2e\n",
572                       alpha_primal);
573
574        try {
575          // Compute the primal trial point
576          IpData().SetTrialPrimalVariablesFromStep(alpha_primal, *actual_delta->x(), *actual_delta->s());
577
578          if (magic_steps_) {
579            PerformMagicStep();
580          }
581
582          // If it is acceptable, stop the search
583          alpha_primal_test = alpha_primal;
584          if (accept_every_trial_step_) {
585            // We call the evaluation at the trial point here, so that an
586            // exception will the thrown if there are problem during the
587            // evaluation of the functions (in that case, we want to further
588            // reduce the step size
589            IpCq().trial_barrier_obj();
590            IpCq().trial_constraint_violation();
591            accept = true;
592          }
593          else {
594            accept = acceptor_->CheckAcceptabilityOfTrialPoint(alpha_primal_test);
595          }
596        }
597        catch(IpoptNLP::Eval_Error& e) {
598          e.ReportException(Jnlst(), J_DETAILED);
599          Jnlst().Printf(J_WARNING, J_LINE_SEARCH,
600                         "Warning: Cutting back alpha due to evaluation error\n");
601          IpData().Append_info_string("e");
602          accept = false;
603          evaluation_error = true;
604        }
605
606        if (accept) {
607          break;
608        }
609
610        if (in_watchdog_) {
611          break;
612        }
613
614        // Decide if we want to go to the restoration phase in a
615        // short cut to check if the problem is infeasible
616        if (expect_infeasible_problem_) {
617          if (count_successive_shortened_steps_>=5) {
618            break;
619          }
620        }
621
622        // try second order correction step if the function could
623        // be evaluated
624        // DoTo: check if we want to do SOC when watchdog is active
625        if (!evaluation_error) {
626          Number theta_curr = IpCq().curr_constraint_violation();
627          Number theta_trial = IpCq().trial_constraint_violation();
628          if (alpha_primal==alpha_primal_max &&       // i.e. first trial point
629              theta_curr<=theta_trial) {
630            // Try second order correction
631            accept = acceptor_->TrySecondOrderCorrection(alpha_primal_test,
632                     alpha_primal,
633                     actual_delta);
634          }
635          if (accept) {
636            soc_taken = true;
637            break;
638          }
639        }
640
641        // Point is not yet acceptable, try a shorter one
642        alpha_primal *= alpha_red_factor_;
643        n_steps++;
644      }
645    } /* if (!accept) */
646
647    char info_alpha_primal_char;
648    if (!accept && in_watchdog_) {
649      info_alpha_primal_char = 'w';
650    }
651    else if (accept) {
652      info_alpha_primal_char =
653        acceptor_->UpdateForNextIteration(alpha_primal_test);
654    }
655    if (soc_taken) {
656      info_alpha_primal_char = toupper(info_alpha_primal_char);
657    }
658    IpData().Set_info_alpha_primal_char(info_alpha_primal_char);
659    IpData().Set_info_ls_count(n_steps+1);
660    if (corr_taken) {
661      IpData().Append_info_string("C");
662    }
663
664    return accept;
665  }
666
667  void BacktrackingLineSearch::StartWatchDog()
668  {
669    DBG_START_FUN("BacktrackingLineSearch::StartWatchDog", dbg_verbosity);
670
671    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
672                   "Starting Watch Dog\n");
673
674    in_watchdog_ = true;
675    watchdog_iterate_ = IpData().curr();
676    watchdog_delta_ = IpData().delta();
677    watchdog_trial_iter_ = 0;
678    watchdog_alpha_primal_test_ =
679      IpCq().curr_primal_frac_to_the_bound(IpData().curr_tau());
680
681    acceptor_->StartWatchDog();
682  }
683
684  void BacktrackingLineSearch::StopWatchDog(SmartPtr<IteratesVector>& actual_delta)
685  {
686    DBG_START_FUN("BacktrackingLineSearch::StopWatchDog", dbg_verbosity);
687
688    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
689                   "Stopping Watch Dog\n");
690
691    IpData().Append_info_string("w");
692
693    in_watchdog_ = false;
694
695    // Reset all fields in IpData to reference point
696    SmartPtr<IteratesVector> old_trial = watchdog_iterate_->MakeNewContainer();
697    IpData().set_trial(old_trial);
698    IpData().AcceptTrialPoint();
699    actual_delta = watchdog_delta_->MakeNewContainer();
700    IpData().SetHaveAffineDeltas(false);
701
702    // reset the stored watchdog iterates
703    watchdog_iterate_ = NULL;
704    watchdog_delta_ = NULL;
705
706    watchdog_shortened_iter_ = 0;
707
708    acceptor_->StopWatchDog();
709  }
710
711  void BacktrackingLineSearch::Reset()
712  {
713    DBG_START_FUN("BacktrackingLineSearch::Reset", dbg_verbosity);
714    in_soft_resto_phase_ = false;
715
716    // Inactivate the watchdog and release all stored data
717    in_watchdog_ = false;
718    watchdog_iterate_ = NULL;
719    watchdog_delta_ = NULL;
720    watchdog_shortened_iter_ = 0;
721
722    acceptor_->Reset();
723  }
724
725  void BacktrackingLineSearch::PerformDualStep(Number alpha_primal,
726      Number alpha_dual,
727      SmartPtr<IteratesVector>& delta)
728  {
729    DBG_START_FUN("BacktrackingLineSearch::PerformDualStep", dbg_verbosity);
730
731    // set the bound multipliers from the step
732    IpData().SetTrialBoundMultipliersFromStep(alpha_dual, *delta->z_L(), *delta->z_U(), *delta->v_L(), *delta->v_U());
733
734    Number alpha_y=-1.;
735    switch (alpha_for_y_) {
736      case PRIMAL_ALPHA_FOR_Y:
737      alpha_y = alpha_primal;
738      break;
739      case DUAL_ALPHA_FOR_Y:
740      alpha_y = alpha_dual;
741      break;
742      case MIN_ALPHA_FOR_Y:
743      alpha_y = Min(alpha_dual, alpha_primal);
744      break;
745      case MAX_ALPHA_FOR_Y:
746      alpha_y = Max(alpha_dual, alpha_primal);
747      break;
748      case FULL_STEP_FOR_Y:
749      alpha_y = 1;
750      break;
751      case MIN_DUAL_INFEAS_ALPHA_FOR_Y:
752      case SAFE_MIN_DUAL_INFEAS_ALPHA_FOR_Y:
753      // Here we compute the step size for y so that the dual
754      // infeasibility is minimized along delta_y
755
756      // compute the dual infeasibility at new point with old y
757      SmartPtr<IteratesVector> temp_trial
758      = IpData().trial()->MakeNewContainer();
759      temp_trial->Set_y_c(*IpData().curr()->y_c());
760      temp_trial->Set_y_d(*IpData().curr()->y_d());
761      IpData().set_trial(temp_trial);
762      SmartPtr<const Vector> dual_inf_x = IpCq().trial_grad_lag_x();
763      SmartPtr<const Vector> dual_inf_s = IpCq().trial_grad_lag_s();
764
765      SmartPtr<Vector> new_jac_times_delta_y =
766        IpData().curr()->x()->MakeNew();
767      new_jac_times_delta_y->AddTwoVectors(1., *IpCq().trial_jac_cT_times_vec(*delta->y_c()),
768                                           1., *IpCq().trial_jac_dT_times_vec(*delta->y_d()),
769                                           0.);
770
771      Number a = pow(new_jac_times_delta_y->Nrm2(), 2.) +
772                 pow(delta->y_d()->Nrm2(), 2.);
773      Number b = dual_inf_x->Dot(*new_jac_times_delta_y) -
774                 dual_inf_s->Dot(*delta->y_d());
775
776      Number alpha = - b/a;
777
778      if (alpha_for_y_==SAFE_MIN_DUAL_INFEAS_ALPHA_FOR_Y) {
779        alpha_y = Min(Max(alpha_primal, alpha_dual), Max(alpha, Min(alpha_primal, alpha_dual)));
780      }
781      else {
782        alpha_y = Min(1., Max(0., alpha));
783      }
784      break;
785    }
786
787    // Set the eq multipliers from the step now that alpha_y
788    // has been calculated.
789    DBG_PRINT((1, "alpha_y = %e\n", alpha_y));
790    DBG_PRINT_VECTOR(2, "delta_y_c", *delta->y_c());
791    DBG_PRINT_VECTOR(2, "delta_y_d", *delta->y_d());
792    IpData().SetTrialEqMultipliersFromStep(alpha_y, *delta->y_c(), *delta->y_d());
793
794    // Set some information for iteration summary output
795    IpData().Set_info_alpha_primal(alpha_primal);
796    IpData().Set_info_alpha_dual(alpha_dual);
797  }
798
799  void
800  BacktrackingLineSearch::PerformMagicStep()
801  {
802    DBG_START_METH("BacktrackingLineSearch::PerformMagicStep",
803                   2);//dbg_verbosity);
804
805    DBG_PRINT((1,"Incoming barr = %e and constrviol %e\n",
806               IpCq().trial_barrier_obj(),
807               IpCq().trial_constraint_violation()));
808    DBG_PRINT_VECTOR(2, "s in", *IpData().trial()->s());
809    DBG_PRINT_VECTOR(2, "d minus s in", *IpCq().trial_d_minus_s());
810    DBG_PRINT_VECTOR(2, "slack_s_L in", *IpCq().trial_slack_s_L());
811    DBG_PRINT_VECTOR(2, "slack_s_U in", *IpCq().trial_slack_s_U());
812
813    SmartPtr<const Vector> d_L = IpNLP().d_L();
814    SmartPtr<const Matrix> Pd_L = IpNLP().Pd_L();
815    SmartPtr<Vector> delta_s_magic_L = d_L->MakeNew();
816    delta_s_magic_L->Set(0.);
817    SmartPtr<Vector> tmp = d_L->MakeNew();
818    Pd_L->TransMultVector(1., *IpCq().trial_d_minus_s(), 0., *tmp);
819    delta_s_magic_L->ElementWiseMax(*tmp);
820
821    SmartPtr<const Vector> d_U = IpNLP().d_U();
822    SmartPtr<const Matrix> Pd_U = IpNLP().Pd_U();
823    SmartPtr<Vector> delta_s_magic_U = d_U->MakeNew();
824    delta_s_magic_U->Set(0.);
825    tmp = d_U->MakeNew();
826    Pd_U->TransMultVector(1., *IpCq().trial_d_minus_s(), 0., *tmp);
827    delta_s_magic_U->ElementWiseMin(*tmp);
828
829    SmartPtr<Vector> delta_s_magic = IpData().trial()->s()->MakeNew();
830    Pd_L->MultVector(1., *delta_s_magic_L, 0., *delta_s_magic);
831    Pd_U->MultVector(1., *delta_s_magic_U, 1., *delta_s_magic);
832    delta_s_magic_L = NULL; // free memory
833    delta_s_magic_U = NULL; // free memory
834
835    // Now find those entries with both lower and upper bounds, there
836    // the step is too large
837    // ToDo this should only be done if there are inequality
838    // constraints with two bounds
839    // also this can be done in a smaller space (d_L or d_U whichever
840    // is smaller)
841    tmp = delta_s_magic->MakeNew();
842    tmp->Copy(*IpData().trial()->s());
843    Pd_L->MultVector(1., *d_L, -2., *tmp);
844    Pd_U->MultVector(1., *d_U, 1., *tmp);
845    SmartPtr<Vector> tmp2 = tmp->MakeNew();
846    tmp2->Copy(*tmp);
847    tmp2->ElementWiseAbs();
848    tmp->Axpy(-2., *delta_s_magic);
849    tmp->ElementWiseAbs();
850    // now, tmp2 = |d_L + d_u - 2*s| and tmp = |d_L + d_u - 2*(s+Delta s)|
851    // we want to throw out those for which tmp2 > tmp
852    tmp->Axpy(-1., *tmp2);
853    tmp->ElementWiseSgn();
854    tmp2->Set(0.);
855    tmp2->ElementWiseMax(*tmp);
856    tmp = d_L->MakeNew();
857    Pd_L->TransMultVector(1., *tmp2, 0., *tmp);
858    Pd_L->MultVector(1., *tmp, 0., *tmp2);
859    tmp = d_U->MakeNew();
860    Pd_U->TransMultVector(1., *tmp2, 0., *tmp);
861    Pd_U->MultVector(1., *tmp, 0., *tmp2);
862    DBG_PRINT_VECTOR(2, "tmp indicator", *tmp2)
863    // tmp2 now is one for those entries with both bounds, for which
864    // no step should be taken
865
866    tmp = delta_s_magic->MakeNew();
867    tmp->Copy(*delta_s_magic);
868    tmp->ElementWiseMultiply(*tmp2);
869    delta_s_magic->Axpy(-1., *tmp);
870
871    Number delta_s_magic_max = delta_s_magic->Amax();
872    Number mach_eps = std::numeric_limits<Number>::epsilon();
873    if (delta_s_magic_max>0.) {
874      if (delta_s_magic_max > 10*mach_eps*IpData().trial()->s()->Amax()) {
875        IpData().Append_info_string("M");
876        Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Magic step with max-norm %.6e taken.\n", delta_s_magic->Amax());
877        delta_s_magic->Print(Jnlst(), J_MOREVECTOR, J_LINE_SEARCH,
878                             "delta_s_magic");
879      }
880
881      // now finally compute the new overall slacks
882      delta_s_magic->Axpy(1., *IpData().trial()->s());
883      SmartPtr<IteratesVector> trial = IpData().trial()->MakeNewContainer();
884      trial->Set_s(*delta_s_magic);
885
886      // also update the set in the dual variables
887
888
889      IpData().set_trial(trial);
890    }
891
892    DBG_PRINT((1,"Outgoing barr = %e and constrviol %e\n", IpCq().trial_barrier_obj(), IpCq().trial_constraint_violation()));
893    DBG_PRINT_VECTOR(2, "s out", *IpData().trial()->s());
894    DBG_PRINT_VECTOR(2, "d minus s out", *IpCq().trial_d_minus_s());
895    DBG_PRINT_VECTOR(2, "slack_s_L out", *IpCq().trial_slack_s_L());
896    DBG_PRINT_VECTOR(2, "slack_s_U out", *IpCq().trial_slack_s_U());
897  }
898
899  bool BacktrackingLineSearch::TrySoftRestoStep(SmartPtr<IteratesVector>& actual_delta,
900      bool &satisfies_original_criterion)
901  {
902    DBG_START_FUN("FilterLSAcceptor::TrySoftRestoStep", dbg_verbosity);
903
904    if (soft_resto_pderror_reduction_factor_==0.) {
905      return false;
906    }
907
908    satisfies_original_criterion = false;
909
910    // ToDo: Need to decide if we want to try a corrector step first
911
912    // Compute the maximal step sizes (we use identical step sizes for
913    // primal and dual variables
914    Number alpha_primal_max =
915      IpCq().primal_frac_to_the_bound(IpData().curr_tau(),
916                                      *actual_delta->x(),
917                                      *actual_delta->s());
918    Number alpha_dual_max =
919      IpCq().dual_frac_to_the_bound(IpData().curr_tau(),
920                                    *actual_delta->z_L(),
921                                    *actual_delta->z_U(),
922                                    *actual_delta->v_L(),
923                                    *actual_delta->v_U());
924    Number alpha_max =  Min(alpha_primal_max, alpha_dual_max);
925
926    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
927                   "Trying soft restoration phase step with step length %13.6e\n",
928                   alpha_max);
929
930    // Set the trial point
931    IpData().SetTrialPrimalVariablesFromStep(alpha_max, *actual_delta->x(), *actual_delta->s());
932    PerformDualStep(alpha_max, alpha_max,
933                    actual_delta);
934
935    // Check if that point is acceptable with respect to the current
936    // original filter
937
938    Number trial_barr;
939    Number trial_theta;
940    try {
941      trial_barr = IpCq().trial_barrier_obj();
942      trial_theta = IpCq().trial_constraint_violation();
943    }
944    catch(IpoptNLP::Eval_Error& e) {
945      e.ReportException(Jnlst(), J_DETAILED);
946      Jnlst().Printf(J_WARNING, J_LINE_SEARCH,
947                     "Warning: Evaluation error during soft restoration phase step.\n");
948      IpData().Append_info_string("e");
949      return false;
950    }
951    if (acceptor_->CheckAcceptabilityOfTrialPoint(0.)) {
952      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
953                     "  Trial step acceptable with respect to original backtracking globalization.\n");
954      satisfies_original_criterion = true;
955      return true;
956    }
957
958    // Evaluate the optimality error at the new point
959    Number mu = .0;
960    if (!IpData().FreeMuMode()) {
961      mu = IpData().curr_mu();
962    }
963    Number trial_pderror;
964    Number curr_pderror;
965    try {
966      trial_pderror = IpCq().trial_primal_dual_system_error(mu);
967      curr_pderror = IpCq().curr_primal_dual_system_error(mu);
968    }
969    catch(IpoptNLP::Eval_Error& e) {
970      e.ReportException(Jnlst(), J_DETAILED);
971      Jnlst().Printf(J_WARNING, J_LINE_SEARCH,
972                     "Warning: Evaluation error during soft restoration phase step.\n");
973      IpData().Append_info_string("e");
974      return false;
975    }
976
977    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
978                   "  Primal-dual error at current point:  %23.16e\n", curr_pderror);
979    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
980                   "  Primal-dual error at trial point  :  %23.16e\n", trial_pderror);
981    // Check if there is sufficient reduction in the optimality error
982    if (trial_pderror <= soft_resto_pderror_reduction_factor_*curr_pderror) {
983      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
984                     "  Trial step accepted.\n");
985      return true;
986    }
987
988    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
989                   "  Trial step rejected.\n");
990    return false;
991  }
992
993  bool
994  BacktrackingLineSearch::DetectTinyStep()
995  {
996    DBG_START_METH("BacktrackingLineSearch::DetectTinyStep",
997                   dbg_verbosity);
998
999    Number max_step_x;
1000    Number max_step_s;
1001
1002    if (tiny_step_tol_==0.)
1003      return false;
1004
1005    // ToDo try to find more efficient implementation
1006
1007    SmartPtr<Vector> tmp = IpData().curr()->x()->MakeNew();
1008    tmp->Copy(*IpData().curr()->x());
1009    tmp->ElementWiseAbs();
1010    tmp->AddScalar(1.);
1011
1012    SmartPtr<Vector> tmp2 = IpData().curr()->x()->MakeNew();
1013    tmp2->Copy(*IpData().delta()->x());
1014    tmp2->ElementWiseDivide(*tmp);
1015    max_step_x = tmp2->Amax();
1016    Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
1017                   "Relative step size for delta_x = %e\n",
1018                   max_step_x);
1019    if (max_step_x > tiny_step_tol_)
1020      return false;
1021
1022    tmp = IpData().curr()->s()->MakeNew();
1023    tmp->Copy(*IpData().curr()->s());
1024    tmp->ElementWiseAbs();
1025    tmp->AddScalar(1.);
1026
1027    tmp2 = IpData().curr()->s()->MakeNew();
1028    tmp2->Copy(*IpData().delta()->s());
1029    tmp2->ElementWiseDivide(*tmp);
1030    max_step_s = tmp2->Amax();
1031    Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
1032                   "Relative step size for delta_s = %e\n",
1033                   max_step_s);
1034    if (max_step_s > tiny_step_tol_)
1035      return false;
1036
1037    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
1038                   "Tiny step of relative size %e detected.\n",
1039                   Max(max_step_x, max_step_s));
1040
1041    return true;
1042  }
1043
1044  bool BacktrackingLineSearch::CurrentIsAcceptable()
1045  {
1046    return (IsValid(conv_check_) &&
1047            conv_check_->CurrentIsAcceptable());
1048  }
1049
1050  void BacktrackingLineSearch::StoreAcceptablePoint()
1051  {
1052    DBG_START_METH("BacktrackingLineSearch::StoreAcceptablePoint",
1053                   dbg_verbosity);
1054
1055    acceptable_iterate_ = IpData().curr();
1056  }
1057
1058  bool BacktrackingLineSearch::RestoreAcceptablePoint()
1059  {
1060    DBG_START_METH("BacktrackingLineSearch::RestoreAcceptablePoint",
1061                   dbg_verbosity);
1062
1063    if (!IsValid(acceptable_iterate_)) {
1064      return false;
1065    }
1066
1067    SmartPtr<IteratesVector> prev_iterate = acceptable_iterate_->MakeNewContainer();
1068    IpData().set_trial(prev_iterate);
1069    IpData().AcceptTrialPoint();
1070
1071    return true;
1072  }
1073
1074
1075} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.