source: branches/dev/Algorithm/IpFilterLineSearch.cpp @ 529

Last change on this file since 529 was 529, checked in by andreasw, 14 years ago
  • in Vector class, copy cached values in Copy and update them in Scal
  • perform iterative refinement only once for adaptive strategy
  • fix bugs in PDPerturbationHandler
  • avoid some overhead in CalculateQualityFunction?
  • minor changes in timing
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 62.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: IpFilterLineSearch.cpp 529 2005-09-29 21:12:38Z andreasw $
6//
7// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
8
9#include "IpFilterLineSearch.hpp"
10#include "IpJournalist.hpp"
11#include "IpRestoPhase.hpp"
12#include "IpAlgTypes.hpp"
13
14#ifdef HAVE_CMATH
15# include <cmath>
16#else
17# ifdef HAVE_MATH_H
18#  include <math.h>
19# else
20#  error "don't have header file for math"
21# endif
22#endif
23
24#ifdef HAVE_CCTYPE
25# include <cctype>
26#else
27# ifdef HAVE_CTYPE_H
28#  include <ctype.h>
29# else
30#  error "don't have header file for ctype"
31# endif
32#endif
33
34namespace Ipopt
35{
36
37#ifdef IP_DEBUG
38  static const Index dbg_verbosity = 0;
39#endif
40
41  FilterLineSearch::FilterLineSearch(const SmartPtr<RestorationPhase>& resto_phase,
42                                     const SmartPtr<PDSystemSolver>& pd_solver,
43                                     const SmartPtr<ConvergenceCheck>& conv_check)
44      :
45      LineSearch(),
46      filter_(2),
47      resto_phase_(resto_phase),
48      pd_solver_(pd_solver),
49      conv_check_(conv_check)
50  {
51    DBG_START_FUN("FilterLineSearch::FilterLineSearch",
52                  dbg_verbosity);
53  }
54
55  FilterLineSearch::~FilterLineSearch()
56  {
57    DBG_START_FUN("FilterLineSearch::~FilterLineSearch()",
58                  dbg_verbosity);
59  }
60
61  void FilterLineSearch::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
62  {
63    roptions->AddLowerBoundedNumberOption(
64      "theta_max_fact",
65      "Determines upper bound for constraint violation in the filter.",
66      0.0, true, 1e4,
67      "The algorithmic parameter theta_max is determined as theta_max_fact "
68      "times the maximum of 1 and the constraint violation at initial point.  "
69      "Any point with a constraint violation larger than theta_max is "
70      "unacceptable to the filter (see Eqn. (21) in implementation paper).");
71    roptions->AddLowerBoundedNumberOption(
72      "theta_min_fact",
73      "Determines constraint violation threshold in the switching rule.",
74      0.0, true, 1e-4,
75      "The algorithmic parameter theta_min is determined as theta_min_fact "
76      "times the maximum of 1 and the constraint violation at initial point.  "
77      "The switching rules treats an iteration as an h-type iteration whenever "
78      "the current constraint violation is larger than theta_min (see "
79      "paragraph before Eqn. (19) in implementation paper).");
80    roptions->AddBoundedNumberOption(
81      "eta_phi",
82      "Relaxation factor in the Armijo condition.",
83      0.0, true, 0.5, true, 1e-8,
84      "(See Eqn. (20) in implementation paper)");
85    roptions->AddLowerBoundedNumberOption(
86      "delta", "Multiplier for constraint violation in the switching rule.",
87      0.0, true, 1.0,
88      "(See Eqn. (19) in implementation paper.)");
89    roptions->AddLowerBoundedNumberOption(
90      "s_phi",
91      "Exponent for linear barrier function model in the switching rule.",
92      1.0, true, 2.3,
93      "(See Eqn. (19) in implementation paper.)");
94    roptions->AddLowerBoundedNumberOption(
95      "s_theta",
96      "Exponent for current constraint violation in the switching rule.",
97      1.0, true, 1.1,
98      "(See Eqn. (19) in implementation paper.)");
99    roptions->AddBoundedNumberOption(
100      "gamma_phi",
101      "Relaxation factor in the filter margin for the barrier function.",
102      0.0, true, 1.0, true, 1e-8,
103      "(See Eqn. (18a) in implementation paper.)");
104    roptions->AddBoundedNumberOption(
105      "gamma_theta",
106      "Relaxation factor in the filter margin for the constraint violation.",
107      0.0, true, 1.0, true, 1e-5,
108      "(See Eqn. (18b) in implementation paper.)");
109    roptions->AddBoundedNumberOption(
110      "alpha_min_frac",
111      "Safety factor for the minimal step size (before switching to restoration phase).",
112      0.0, true, 1.0, true, 0.05,
113      "(This is gamma_alpha in Eqn. (20) in implementation paper.)");
114    roptions->AddBoundedNumberOption(
115      "alpha_red_factor",
116      "Fractional reduction of the trial step size in the backtracking line search.",
117      0.0, true, 1.0, true, 0.5,
118      "At every step of the backtracking line search, the trial step size is "
119      "reduced by this factor.");
120    roptions->AddLowerBoundedIntegerOption(
121      "max_soc",
122      "Maximum number of second order correction trial steps at each iteration.",
123      0, 4,
124      "Choosing 0 disables the second order "
125      "corrections. (This is p^{max} of Step A-5.9 of "
126      "Algorithm A in implementation paper.)");
127    roptions->AddLowerBoundedNumberOption(
128      "kappa_soc",
129      "Factor in the sufficient reduction rule for second order correction.",
130      0.0, true, 0.99,
131      "This option determines how much a second order correction step must reduce the "
132      "constraint violation so that further correction steps are attempted.  "
133      "(See Step A-5.9 of Algorithm A in implementation paper.)");
134    roptions->AddLowerBoundedNumberOption(
135      "obj_max_inc",
136      "Determines the upper bound on the acceptable increase of barrier objective function.",
137      1.0, true, 5.0,
138      "Trial points are rejected if they lead to an increase in the "
139      "barrier objective function by more than obj_max_inc orders "
140      "of magnitude.");
141
142    std::string prev_category = roptions->RegisteringCategory();
143    roptions->SetRegisteringCategory("Undocumented");
144    roptions->AddStringOption2(
145      "magic_steps",
146      "Enables magic steps.",
147      "no",
148      "no", "don't take magic steps",
149      "yes", "take magic steps",
150      "DOESN'T REALLY WORK YET!");
151    roptions->SetRegisteringCategory(prev_category);
152    roptions->AddStringOption3(
153      "corrector_type",
154      "The type of corrector steps that should be taken.",
155      "none",
156      "none", "no corrector",
157      "affine", "corrector step towards mu=0",
158      "primal-dual", "corrector step towards current mu",
159      "If \"mu_strategy\" is \"adaptive\", this option determines "
160      "what kind of corrector steps should be tried.");
161
162    roptions->AddStringOption2(
163      "skip_corr_if_neg_curv",
164      "Skip the corrector step in negative curvature iteration.",
165      "yes",
166      "no", "don't skip",
167      "yes", "skip",
168      "The corrector step is not tried if negative curvature has been "
169      "encountered during the computation of the search direction in "
170      "the current iteration. This option is only used if \"mu_strategy\" is "
171      "\"adaptive\".");
172
173    roptions->AddStringOption2(
174      "skip_corr_in_monotone_mode",
175      "Skip the corrector step during monotone barrier parameter mode.",
176      "yes",
177      "no", "don't skip",
178      "yes", "skip",
179      "The corrector step is not tried if the algorithm is currently in the "
180      "monotone mode (see also option \"barrier_strategy\")."
181      "This option is only used if \"mu_strategy\" is \"adaptive\".");
182
183    roptions->AddStringOption2(
184      "accept_every_trial_step",
185      "Always accept the frist trial step.",
186      "no",
187      "no", "don't arbitrarily accept the full step",
188      "yes", "always accept the full step",
189      "Setting this option to \"yes\" essentially disables the line search "
190      "and makes the algorithm take aggressive steps.");
191
192    roptions->AddStringOption7(
193      "alpha_for_y",
194      "Method to determine the step size for constraint multipliers.",
195      "primal",
196      "primal", "use primal step size",
197      "bound_mult", "use step size for the bound multipliers",
198      "min", "use the min of primal and bound multipliers",
199      "max", "use the max of primal and bound multipliers",
200      "full", "take a full step of size one",
201      "min_dual_infeas", "choose step size minimizing new dual infeasibility",
202      "safe_min_dual_infeas", "like \"min_dual_infeas\", but safeguarded by \"min\" and \"max\"",
203      "This option determines how the step size (alpha_y) will be calculated when updating the "
204      "constraint multipliers.");
205
206    roptions->AddLowerBoundedNumberOption(
207      "corrector_compl_avrg_red_fact",
208      "Complementarity tolerance factor for accepting corrector step.",
209      0.0, true, 1.0,
210      "This option determines the factor by which complementarity is allowed to increase "
211      "for a corrector step to be accepted.");
212
213    roptions->AddStringOption2(
214      "expect_infeasible_problem",
215      "Enable heuristics to quickly detect an infeasible problem.",
216      "no",
217      "no", "the problem probably be feasible",
218      "yes", "the problem has a good chance to be infeasible",
219      "This options is meant to activate heuristics that may speed up the "
220      "infeasibility determination if you expect that there is a good chance for the problem to be "
221      "infeasible.  In the filter line search procedure, the restoration "
222      "phase is called more qucikly than usually, and more reduction in "
223      "the constraint violation is enforced. If the problem is square, this "
224      "option is enabled automatically.");
225    roptions->AddLowerBoundedNumberOption(
226      "expect_infeasible_problem_ctol",
227      "Threshold for disabling \"expect_infeasible_problem\" option.",
228      0.0, false, 1e-3,
229      "If the constraint violation becomes smaller than this threshold, "
230      "the \"expect_infeasible_problem\" heuristics in the filter line "
231      "search are disabled. If the problem is square, this options is set to "
232      "0.");
233    roptions->AddLowerBoundedNumberOption(
234      "soft_resto_pderror_reduction_factor",
235      "Required reduction in primal-dual error in the soft restoration phase.",
236      0.0, false, (1.0 - 1e-4),
237      "The soft restoration phase attempts to reduce the "
238      "primal-dual error with regular steps. If the damped "
239      "primal-dual step (damped only to satisfy the "
240      "fraction-to-the-boundary rule) is not decreasing the primal-dual error "
241      "by at least this factor, then the regular restoration phase is called. "
242      "Choosing \"0\" here disables the soft "
243      "restoration phase.");
244    roptions->AddStringOption2(
245      "start_with_resto",
246      "Tells algorithm to switch to restoration phase in first iteration.",
247      "no",
248      "no", "don't force start in restoration phase",
249      "yes", "force start in restoration phase",
250      "Setting this option to \"yes\" forces the algorithm to switch to the "
251      "feasibility restoration phase in the first iteration. If the initial "
252      "point is feasible, the algorithm will abort with a failure.");
253    roptions->AddLowerBoundedNumberOption(
254      "tiny_step_tol",
255      "Tolerance for detecting numerically insignificant steps.",
256      0.0, false, 10.0*std::numeric_limits<double>::epsilon(),
257      "If the search direction in the primal variables (x and s) is, in "
258      "relative terms for each component, less than this value, the "
259      "algorithm accepts the full step without line search.  If this happens "
260      "repeatedly, the algorithm will terminate with a corresponding exit "
261      "message. The default value is 10 times machine precision.");
262    roptions->AddLowerBoundedIntegerOption(
263      "watchdog_shortened_iter_trigger",
264      "Number of shortened iterations that trigger the watchdog.",
265      0, 10,
266      "If the number of iterations in which the backtracking line search "
267      "did not accept the first trial point exceedes this number, the "
268      "watchdog procedure is activated.  Choosing \"0\" here disables the "
269      "watchdog procedure.");
270    roptions->AddLowerBoundedIntegerOption(
271      "watchdog_trial_iter_max",
272      "Maximum number of watchdog iterations.",
273      1, 3,
274      "This option determines the number of trial iterations "
275      "allowed before the watchdog "
276      "procedure is aborted and the algorithm returns to the stored point.");
277  }
278
279  bool FilterLineSearch::InitializeImpl(const OptionsList& options,
280                                        const std::string& prefix)
281  {
282    options.GetNumericValue("theta_max_fact", theta_max_fact_, prefix);
283    options.GetNumericValue("theta_min_fact", theta_min_fact_, prefix);
284    ASSERT_EXCEPTION(theta_min_fact_ < theta_max_fact_, OPTION_INVALID,
285                     "Option \"theta_min_fact\": This value must be larger than 0 and less than theta_max_fact.");
286    options.GetNumericValue("eta_phi", eta_phi_, prefix);
287    options.GetNumericValue("delta", delta_, prefix);
288    options.GetNumericValue("s_phi", s_phi_, prefix);
289    options.GetNumericValue("s_theta", s_theta_, prefix);
290    options.GetNumericValue("gamma_phi", gamma_phi_, prefix);
291    options.GetNumericValue("gamma_theta", gamma_theta_, prefix);
292    options.GetNumericValue("alpha_min_frac", alpha_min_frac_, prefix);
293    options.GetNumericValue("alpha_red_factor", alpha_red_factor_, prefix);
294    options.GetIntegerValue("max_soc", max_soc_, prefix);
295    if (max_soc_>0) {
296      ASSERT_EXCEPTION(IsValid(pd_solver_), OPTION_INVALID,
297                       "Option \"max_soc\": This option is non-negative, but no linear solver for computing the SOC given to FilterLineSearch object.");
298    }
299    options.GetNumericValue("kappa_soc", kappa_soc_, prefix);
300    options.GetNumericValue("obj_max_inc", obj_max_inc_, prefix);
301    options.GetBoolValue("magic_steps", magic_steps_, prefix);
302    Index enum_int;
303    options.GetEnumValue("corrector_type", enum_int, prefix);
304    corrector_type_ = CorrectorTypeEnum(enum_int);
305    options.GetBoolValue("skip_corr_if_neg_curv", skip_corr_if_neg_curv_, prefix);
306    options.GetBoolValue("skip_corr_in_monotone_mode", skip_corr_in_monotone_mode_, prefix);
307    options.GetBoolValue("accept_every_trial_step", accept_every_trial_step_, prefix);
308    options.GetEnumValue("alpha_for_y", enum_int, prefix);
309    alpha_for_y_ = AlphaForYEnum(enum_int);
310    options.GetNumericValue("corrector_compl_avrg_red_fact", corrector_compl_avrg_red_fact_, prefix);
311    options.GetNumericValue("expect_infeasible_problem_ctol", expect_infeasible_problem_ctol_, prefix);
312    options.GetBoolValue("expect_infeasible_problem", expect_infeasible_problem_, prefix);
313
314    options.GetBoolValue("start_with_resto", start_with_resto_, prefix);
315
316    bool retvalue = true;
317    if (IsValid(resto_phase_)) {
318      retvalue = resto_phase_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
319                                          options, prefix);
320    }
321
322    options.GetNumericValue("soft_resto_pderror_reduction_factor",
323                            soft_resto_pderror_reduction_factor_, prefix);
324    options.GetNumericValue("tiny_step_tol", tiny_step_tol_, prefix);
325    options.GetIntegerValue("watchdog_trial_iter_max", watchdog_trial_iter_max_, prefix);
326    options.GetIntegerValue("watchdog_shortened_iter_trigger", watchdog_shortened_iter_trigger_, prefix);
327
328    rigorous_ = true;
329    skipped_line_search_ = false;
330    tiny_step_last_iteration_ = false;
331    theta_min_ = -1.;
332    theta_max_ = -1.;
333
334    Reset();
335
336    count_successive_shortened_steps_ = 0;
337
338    acceptable_iterate_ = NULL;
339
340    return retvalue;
341  }
342
343  void FilterLineSearch::FindAcceptableTrialPoint()
344  {
345    DBG_START_METH("FilterLineSearch::FindAcceptableTrialPoint",
346                   dbg_verbosity);
347    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
348                   "--> Starting filter line search in iteration %d <--\n",
349                   IpData().iter_count());
350
351    // If the problem is square, we want to enable the
352    // expect_infeasible_problem option automatically so that the
353    // restoration phase is entered soon
354    if (IpCq().IsSquareProblem()) {
355      expect_infeasible_problem_ = true;
356      expect_infeasible_problem_ctol_ = 0.;
357    }
358
359    // Store current iterate if the optimality error is on acceptable
360    // level to restored if things fail later
361    if (CurrentIsAcceptable()) {
362      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
363                     "Storing current iterate as backup acceptable point.\n");
364      StoreAcceptablePoint();
365    }
366
367    // First assume that line search will find an acceptable trial point
368    skipped_line_search_ = false;
369
370    // Set the values for the reference point
371    if (!in_watchdog_) {
372      reference_theta_ = IpCq().curr_constraint_violation();
373      reference_barr_ = IpCq().curr_barrier_obj();
374      reference_gradBarrTDelta_ = IpCq().curr_gradBarrTDelta();
375    }
376    else {
377      reference_theta_ = watchdog_theta_;
378      reference_barr_ = watchdog_barr_;
379      reference_gradBarrTDelta_ = watchdog_gradBarrTDelta_;
380    }
381
382    // Get the search directions (this will store the actual search
383    // direction, possibly including higher order corrections)
384    SmartPtr<IteratesVector> actual_delta = IpData().delta()->MakeNewContainer();
385
386    bool goto_resto = false;
387    if (actual_delta->Asum() == 0. ) {
388      // In this case, a search direction could not be computed, and
389      // we should immediately go to the restoration phase.  ToDo: Cue
390      // off of a something else than the norm of the search direction
391      goto_resto = true;
392    }
393
394    if (start_with_resto_) {
395      // If the use requested to start with the restoration phase,
396      // skip the lin e search and do exactly that.  Reset the flag so
397      // that this happens only once.
398      goto_resto = true;
399      start_with_resto_= false;
400    }
401
402    bool accept = false;
403    bool corr_taken = false;
404    bool soc_taken = false;
405    Index n_steps = 0;
406    Number alpha_primal = 0.;
407
408    // Check if search direction becomes too small
409    // ToDo: move this into place independent of this particular line search?
410    bool tiny_step = (!goto_resto && DetectTinyStep());
411
412    if (in_watchdog_ && (goto_resto || tiny_step)) {
413      // If the step could not be computed or is too small and the
414      // watchdog is active, stop the watch dog and resume everything
415      // from reference point
416      StopWatchDog(actual_delta);
417      goto_resto = false;
418      tiny_step = false;
419    }
420
421    // Check if we want to wake up the watchdog
422    if (watchdog_shortened_iter_trigger_ > 0 &&
423        !in_watchdog_ && !goto_resto && !tiny_step &&
424        !in_soft_resto_phase_ && !expect_infeasible_problem_ &&
425        watchdog_shortened_iter_ >= watchdog_shortened_iter_trigger_) {
426      StartWatchDog();
427    }
428
429    if (tiny_step) {
430      alpha_primal =
431        IpCq().curr_primal_frac_to_the_bound(IpData().curr_tau());
432      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
433                     "Tiny step detected. Use step size alpha = %e unchecked\n",
434                     alpha_primal);
435      IpData().SetTrialPrimalVariablesFromStep(alpha_primal, *IpData().delta()->x(), *IpData().delta()->s());
436      IpData().Set_info_ls_count(0);
437
438      if (tiny_step_last_iteration_) {
439        IpData().Set_info_alpha_primal_char('T');
440        IpData().Set_tiny_step_flag(true);
441      }
442      else {
443        IpData().Set_info_alpha_primal_char('t');
444      }
445
446      tiny_step_last_iteration_ = true;
447      accept = true;
448    }
449    else {
450      tiny_step_last_iteration_ = false;
451    }
452
453    if (!goto_resto && !tiny_step) {
454
455      if (in_soft_resto_phase_) {
456        bool satisfies_original_filter = false;
457        // ToDo use tiny_step in TrySoftRestoStep?
458        accept = TrySoftRestoStep(actual_delta,
459                                  satisfies_original_filter);
460        if (accept) {
461          IpData().Set_info_alpha_primal_char('s');
462          if (satisfies_original_filter) {
463            in_soft_resto_phase_ = false;
464            IpData().Set_info_alpha_primal_char('S');
465          }
466        }
467      }
468      else {
469        bool done = false;
470        bool skip_first_trial_point = false;
471        bool evaluation_error;
472        while (!done) {
473          accept = DoBacktrackingLineSearch(skip_first_trial_point,
474                                            alpha_primal,
475                                            corr_taken,
476                                            soc_taken,
477                                            n_steps,
478                                            evaluation_error,
479                                            actual_delta);
480          DBG_PRINT((1, "evaluation_error = %d\n", evaluation_error));
481          if (in_watchdog_) {
482            if (accept) {
483              in_watchdog_ = false;
484              IpData().Append_info_string("W");
485              Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
486                             "Watch dog procedure successful!\n");
487              done = true;
488            }
489            else {
490              watchdog_trial_iter_++;
491              if (evaluation_error ||
492                  watchdog_trial_iter_ > watchdog_trial_iter_max_) {
493                StopWatchDog(actual_delta);
494                skip_first_trial_point = true;
495              }
496              else {
497                done = true;
498                accept = true;
499              }
500            }
501          }
502          else {
503            done = true;
504          }
505        }
506      } /* else: if (in_soft_resto_phase_) { */
507    } /* if (!goto_resto && !tiny_step) { */
508
509    // If line search has been aborted because the step size becomes too small,
510    // go to the restoration phase
511    if (!accept) {
512      // If we are not asked to do a rigorous line search, do no call
513      // the restoration phase.
514      if (!rigorous_) {
515        Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Skipping call of restoration phase...\n");
516        skipped_line_search_=true;
517      }
518      else {
519        // Check if we should start the soft restoration phase
520        if (!in_soft_resto_phase_ && soft_resto_pderror_reduction_factor_>0.
521            && !goto_resto && !expect_infeasible_problem_) {
522          Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
523                         "--> Starting soft restoration phase <--\n");
524          // Augment the filter with the current point
525          AugmentFilter();
526
527          // Try the current search direction for the soft restoration phase
528          bool satisfies_original_filter;
529          accept = TrySoftRestoStep(actual_delta,
530                                    satisfies_original_filter);
531          // If it has been accepted: If the original filter is also
532          // satisfied, we can just take that step and continue with
533          // the regular algorithm, otherwise we stay in the soft
534          // restoration phase
535          if (accept) {
536            if (satisfies_original_filter) {
537              IpData().Set_info_alpha_primal_char('S');
538            }
539            else {
540              in_soft_resto_phase_ = true;
541              IpData().Set_info_alpha_primal_char('s');
542            }
543          }
544        }
545
546        if (!accept) {
547          if (!in_soft_resto_phase_) {
548            // Augment the filter with the current point if we are
549            // already in the soft restoration phase, this has been
550            // done earlier
551            AugmentFilter();
552          }
553          if (CurrentIsAcceptable()) {
554            THROW_EXCEPTION(ACCEPTABLE_POINT_REACHED,
555                            "Restoration phase called at acceptable point.");
556          }
557
558          if (!IsValid(resto_phase_)) {
559            THROW_EXCEPTION(IpoptException, "No Restoration Phase given to this Filter Line Search Object!");
560          }
561          // ToDo make the 1e-2 below a parameter?
562          if (IpCq().curr_constraint_violation()<=
563              1e-2*IpData().tol()) {
564            bool found_acceptable = RestoreAcceptablePoint();
565            if (found_acceptable) {
566              THROW_EXCEPTION(ACCEPTABLE_POINT_REACHED,
567                              "Restoration phase called at almost feasible point, but acceptable point could be restored.\n");
568            }
569            else {
570              // ToDo does that happen too often?
571              THROW_EXCEPTION(RESTORATION_FAILED,
572                              "Restoration phase called, but point is almost feasible.");
573            }
574          }
575
576          // Set the info fields for the first output line in the
577          // restoration phase which reflects why the restoration phase
578          // was called
579          IpData().Set_info_alpha_primal(alpha_primal);
580          IpData().Set_info_alpha_dual(0.);
581          IpData().Set_info_alpha_primal_char('R');
582          IpData().Set_info_ls_count(n_steps+1);
583
584          accept = resto_phase_->PerformRestoration();
585          if (!accept) {
586            bool found_acceptable = RestoreAcceptablePoint();
587            if (found_acceptable) {
588              THROW_EXCEPTION(ACCEPTABLE_POINT_REACHED,
589                              "Restoration phase failed, but acceptable point could be restore.\n");
590            }
591            else {
592              THROW_EXCEPTION(RESTORATION_FAILED,
593                              "Failed restoration phase!!!");
594            }
595          }
596          count_successive_shortened_steps_ = 0;
597          if (expect_infeasible_problem_) {
598            expect_infeasible_problem_ = false;
599          }
600          in_soft_resto_phase_ = false;
601          watchdog_shortened_iter_ = 0;
602        }
603      }
604    }
605    else if (!in_soft_resto_phase_ || tiny_step) {
606      // we didn't do the restoration phase and are now updating the
607      // dual variables of the trial point
608      Number alpha_dual_max =
609        IpCq().dual_frac_to_the_bound(IpData().curr_tau(),
610                                      *actual_delta->z_L(), *actual_delta->z_U(),
611                                      *actual_delta->v_L(), *actual_delta->v_U());
612
613      PerformDualStep(alpha_primal, alpha_dual_max, actual_delta);
614
615      if (n_steps==0) {
616        // accepted this if a full step was
617        // taken
618        count_successive_shortened_steps_ = 0;
619        watchdog_shortened_iter_ = 0;
620      }
621      else {
622        count_successive_shortened_steps_++;
623        watchdog_shortened_iter_++;
624      }
625
626      if (expect_infeasible_problem_ &&
627          IpCq().curr_constraint_violation() <= expect_infeasible_problem_ctol_) {
628        Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
629                       "Constraint violation is with %e less than expect_infeasible_problem_ctol.\nDisable expect_infeasible_problem_heuristic.\n", IpCq().curr_constraint_violation());
630        expect_infeasible_problem_ = false;
631      }
632    }
633  }
634
635  bool FilterLineSearch::DoBacktrackingLineSearch(bool skip_first_trial_point,
636      Number& alpha_primal,
637      bool& corr_taken,
638      bool& soc_taken,
639      Index& n_steps,
640      bool& evaluation_error,
641      SmartPtr<IteratesVector>& actual_delta)
642  {
643    evaluation_error = false;
644    bool accept = false;
645
646    DBG_START_METH("FilterLineSearch::DoBacktrackingLineSearch",
647                   dbg_verbosity);
648
649    // Compute primal fraction-to-the-boundary value
650    Number alpha_primal_max =
651      IpCq().primal_frac_to_the_bound(IpData().curr_tau(),
652                                      *actual_delta->x(),
653                                      *actual_delta->s());
654
655    // Compute smallest step size allowed
656    Number alpha_min = alpha_primal_max;
657    if (!in_watchdog_) {
658      alpha_min = CalculateAlphaMin();
659    }
660    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
661                   "minimal step size ALPHA_MIN = %E\n", alpha_min);
662
663    // Start line search from maximal step size
664    alpha_primal = alpha_primal_max;
665
666    // Step size used in ftype and armijo tests
667    Number alpha_primal_test = alpha_primal;
668    if (in_watchdog_) {
669      alpha_primal_test = watchdog_alpha_primal_test_;
670    }
671
672    if (skip_first_trial_point) {
673      alpha_primal *= alpha_red_factor_;
674    }
675
676    filter_.Print(Jnlst());
677
678    if (corrector_type_!=NO_CORRECTOR && !skip_first_trial_point &&
679        (!skip_corr_if_neg_curv_ || IpData().info_regu_x()==0.) &&
680        (!skip_corr_in_monotone_mode_ || IpData().FreeMuMode()) ) {
681      // Before we do the actual backtracking line search for the
682      // regular primal-dual search direction, let's see if a step
683      // including a higher-order correctior is already acceptable
684      accept = TryCorrector(alpha_primal_test,
685                            alpha_primal,
686                            actual_delta);
687    }
688    if (accept) {
689      corr_taken = true;
690    }
691
692    if (!accept) {
693      // Loop over decreaseing step sizes until acceptable point is
694      // found or until step size becomes too small
695
696      while (alpha_primal>alpha_min ||
697             n_steps == 0) { // always allow the "full" step if it is
698        // acceptable (even if alpha_primal<=alpha_min)
699        Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
700                       "Starting checks for alpha (primal) = %8.2e\n",
701                       alpha_primal);
702
703        try {
704          // Compute the primal trial point
705          IpData().SetTrialPrimalVariablesFromStep(alpha_primal, *actual_delta->x(), *actual_delta->s());
706
707          if (magic_steps_) {
708            PerformMagicStep();
709          }
710
711          // If it is acceptable, stop the search
712          alpha_primal_test = alpha_primal;
713          accept = CheckAcceptabilityOfTrialPoint(alpha_primal_test);
714        }
715        catch(IpoptNLP::Eval_Error& e) {
716          e.ReportException(Jnlst(), J_DETAILED);
717          Jnlst().Printf(J_WARNING, J_LINE_SEARCH,
718                         "Warning: Cutting back alpha due to evaluation error\n");
719          IpData().Append_info_string("e");
720          accept = false;
721          evaluation_error = true;
722        }
723
724        if (accept) {
725          break;
726        }
727
728        if (in_watchdog_) {
729          break;
730        }
731
732        // Decide if we want to go to the restoration phase in a
733        // short cut to check if the problem is infeasible
734        if (expect_infeasible_problem_) {
735          if (count_successive_shortened_steps_>=5) {
736            break;
737          }
738        }
739
740        // try second order correction step if the function could
741        // be evaluated
742        // DoTo: check if we want to do SOC when watchdog is active
743        if (!evaluation_error) {
744          Number theta_curr = IpCq().curr_constraint_violation();
745          Number theta_trial = IpCq().trial_constraint_violation();
746          if (alpha_primal==alpha_primal_max &&       // i.e. first trial point
747              theta_curr<=theta_trial && max_soc_>0) {
748            // Try second order correction
749            accept = TrySecondOrderCorrection(alpha_primal_test,
750                                              alpha_primal,
751                                              actual_delta);
752          }
753          if (accept) {
754            soc_taken = true;
755            break;
756          }
757        }
758
759        // Point is not yet acceptable, try a shorter one
760        alpha_primal *= alpha_red_factor_;
761        n_steps++;
762      }
763    } /* if (!accept) */
764
765    char info_alpha_primal_char;
766    if (!accept && in_watchdog_) {
767      info_alpha_primal_char = 'w';
768    }
769    else {
770      // Augment the filter if required
771      if (!IsFtype(alpha_primal_test) ||
772          !ArmijoHolds(alpha_primal_test)) {
773        AugmentFilter();
774        info_alpha_primal_char = 'h';
775      }
776      else {
777        info_alpha_primal_char = 'f';
778      }
779    }
780    if (soc_taken) {
781      info_alpha_primal_char = toupper(info_alpha_primal_char);
782    }
783    IpData().Set_info_alpha_primal_char(info_alpha_primal_char);
784    IpData().Set_info_ls_count(n_steps+1);
785    if (corr_taken) {
786      IpData().Append_info_string("C");
787    }
788
789    return accept;
790  }
791
792  bool FilterLineSearch::IsFtype(Number alpha_primal_test)
793  {
794    DBG_START_METH("FilterLineSearch::IsFtype",
795                   dbg_verbosity);
796    DBG_ASSERT(reference_theta_>0. || reference_gradBarrTDelta_ < 0.0);
797    return (reference_gradBarrTDelta_ < 0.0 &&
798            alpha_primal_test*pow(-reference_gradBarrTDelta_,s_phi_) >
799            delta_*pow(reference_theta_,s_theta_));
800  }
801
802  void FilterLineSearch::AugmentFilter()
803  {
804    DBG_START_METH("FilterLineSearch::AugmentFilter",
805                   dbg_verbosity);
806
807    Number phi_add = reference_barr_ - gamma_phi_*reference_theta_;
808    Number theta_add = (1.-gamma_theta_)*reference_theta_;
809
810    filter_.AddEntry(phi_add, theta_add, IpData().iter_count());
811  }
812
813  bool
814  FilterLineSearch::CheckAcceptabilityOfTrialPoint(Number alpha_primal_test)
815  {
816    DBG_START_METH("FilterLineSearch::CheckAcceptabilityOfTrialPoint",
817                   dbg_verbosity);
818
819    if (accept_every_trial_step_) {
820
821      // We call the evaluation at the trial point here, so that an
822      // exception will the thrown if there are problem during the
823      // evaluation of the functions (in that case, we want to further
824      // reduce the step size
825      /* Number trial_barr = */ IpCq().trial_barrier_obj();
826      /* Number trial_theta = */
827      IpCq().trial_constraint_violation();
828      return true;
829    }
830
831    bool accept;
832
833    // First compute the barrier function and constraint violation at the
834    // current iterate and the trial point
835
836    Number trial_theta = IpCq().trial_constraint_violation();
837    // Check if constraint violation is becoming too large
838    if (theta_max_ < 0.0) {
839      theta_max_ = theta_max_fact_*Max(1.0, reference_theta_);
840    }
841    if (theta_min_ < 0.0) {
842      theta_min_ = theta_min_fact_*Max(1.0, reference_theta_);
843    }
844
845    if (theta_max_>0 && trial_theta>theta_max_) {
846      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
847                     "trial_theta = %e is larger than theta_max = %e\n",
848                     trial_theta, theta_max_);
849      IpData().Append_info_string("Tmax");
850      return false;
851    }
852
853    Number trial_barr = IpCq().trial_barrier_obj();
854    DBG_ASSERT(IsFiniteNumber(trial_barr));
855
856    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
857                   "Checking acceptability for trial step size alpha_primal_test=%13.6e:\n", alpha_primal_test);
858    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
859                   "  New values of barrier function     = %23.16e  (reference %23.16e):\n", trial_barr, reference_barr_);
860    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
861                   "  New values of constraint violation = %23.16e  (reference %23.16e):\n", trial_theta, reference_theta_);
862
863    // Check if point is acceptable w.r.t current iterate
864    if (IsFtype(alpha_primal_test) && reference_theta_ <= theta_min_) {
865      // Armijo condition for the barrier function has to be satisfied
866      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Checking Armijo Condition...\n");
867      accept = ArmijoHolds(alpha_primal_test);
868    }
869    else {
870      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Checking sufficient reduction...\n");
871      accept = IsAcceptableToCurrentIterate(trial_barr, trial_theta);
872    }
873
874    if (!accept) {
875      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Failed...\n");
876      return accept;
877    }
878    else {
879      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Succeeded...\n");
880    }
881
882    // Now check if that pair is acceptable to the filter
883    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Checking filter acceptability...\n");
884    accept = IsAcceptableToCurrentFilter(trial_barr, trial_theta);
885    if (!accept) {
886      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Failed...\n");
887      return accept;
888    }
889    else {
890      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Succeeded...\n");
891    }
892
893    return accept;
894  }
895
896  bool FilterLineSearch::ArmijoHolds(Number alpha_primal_test)
897  {
898    /*
899    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
900                   "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_);
901    */
902    return Compare_le(IpCq().trial_barrier_obj()-reference_barr_,
903                      eta_phi_*alpha_primal_test*reference_gradBarrTDelta_,
904                      reference_barr_);
905  }
906
907  Number FilterLineSearch::CalculateAlphaMin()
908  {
909    Number gBD = IpCq().curr_gradBarrTDelta();
910    Number curr_theta = IpCq().curr_constraint_violation();
911    Number alpha_min = gamma_theta_;
912
913    if (gBD < 0) {
914      alpha_min = Min( gamma_theta_,
915                       gamma_phi_*curr_theta/(-gBD));
916      if (curr_theta <= theta_min_) {
917        alpha_min = Min( alpha_min,
918                         delta_*pow(curr_theta,s_theta_)/pow(-gBD,s_phi_)
919                       );
920      }
921    }
922
923    return alpha_min_frac_*alpha_min;
924  }
925
926  bool FilterLineSearch::IsAcceptableToCurrentIterate(Number trial_barr,
927      Number trial_theta,
928      bool called_from_restoration /*=false*/) const
929  {
930    DBG_START_METH("FilterLineSearch::IsAcceptableToCurrentIterate",
931                   dbg_verbosity);
932
933    // Check if the barrier objective function is increasing to
934    // rapidly (according to option obj_max_inc)
935    if (!called_from_restoration && trial_barr > reference_barr_) {
936      Number basval = 1.;
937      if (fabs(reference_barr_)>10.) {
938        basval = log10(fabs(reference_barr_));
939      }
940      if (log10(trial_barr-reference_barr_)>obj_max_inc_+basval) {
941        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);
942        return false;
943      }
944    }
945
946    DBG_PRINT((1,"trial_barr  = %e reference_barr  = %e\n", trial_barr, reference_barr_));
947    DBG_PRINT((1,"trial_theta = %e reference_theta = %e\n", trial_theta, reference_theta_));
948    return (Compare_le(trial_theta, (1.-gamma_theta_)*reference_theta_, reference_theta_)
949            || Compare_le(trial_barr-reference_barr_, -gamma_phi_*reference_theta_, reference_barr_));
950  }
951
952  bool FilterLineSearch::IsAcceptableToCurrentFilter(Number trial_barr, Number trial_theta) const
953  {
954    return filter_.Acceptable(trial_barr, trial_theta);
955  }
956
957  bool FilterLineSearch::Compare_le(Number lhs, Number rhs, Number BasVal)
958  {
959    DBG_START_FUN("FilterLineSearch::Compare_le",
960                  dbg_verbosity);
961    DBG_PRINT((1,"lhs = %27.16e rhs = %27.16e  BasVal = %27.16e\n",lhs,rhs,BasVal));
962    // ToDo: Comparison based on machine precision
963    return (lhs - rhs <= 1e-15*fabs(BasVal));
964  }
965
966  bool FilterLineSearch::TrySoftRestoStep(SmartPtr<IteratesVector>& actual_delta,
967                                          bool &satisfies_original_filter)
968  {
969    DBG_START_FUN("FilterLineSearch::TrySoftRestoStep", dbg_verbosity);
970
971    satisfies_original_filter = false;
972
973    // ToDo: Need to decide if we want to try a corrector step first
974
975    // Compute the maximal step sizes (we use identical step sizes for
976    // primal and dual variables
977    Number alpha_primal_max =
978      IpCq().primal_frac_to_the_bound(IpData().curr_tau(),
979                                      *actual_delta->x(),
980                                      *actual_delta->s());
981    Number alpha_dual_max =
982      IpCq().dual_frac_to_the_bound(IpData().curr_tau(),
983                                    *actual_delta->z_L(),
984                                    *actual_delta->z_U(),
985                                    *actual_delta->v_L(),
986                                    *actual_delta->v_U());
987    Number alpha_max =  Min(alpha_primal_max, alpha_dual_max);
988
989    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
990                   "Trying soft restoration phase step with step length %13.6e\n",
991                   alpha_max);
992
993    // Set the trial point
994    IpData().SetTrialPrimalVariablesFromStep(alpha_max, *actual_delta->x(), *actual_delta->s());
995    PerformDualStep(alpha_max, alpha_max,
996                    actual_delta);
997
998    // Check if that point is acceptable with respect to the current
999    // original filter
1000
1001    Number trial_barr;
1002    Number trial_theta;
1003    try {
1004      trial_barr = IpCq().trial_barrier_obj();
1005      trial_theta = IpCq().trial_constraint_violation();
1006    }
1007    catch(IpoptNLP::Eval_Error& e) {
1008      e.ReportException(Jnlst(), J_DETAILED);
1009      Jnlst().Printf(J_WARNING, J_LINE_SEARCH,
1010                     "Warning: Evaluation error during soft restoration phase step.\n");
1011      IpData().Append_info_string("e");
1012      return false;
1013    }
1014    if (theta_max_<=0 || trial_theta<=theta_max_) {
1015      if (IsAcceptableToCurrentIterate(trial_barr, trial_theta)) {
1016        Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
1017                       "  Trial step acceptable with respect to original filter.\n");
1018        satisfies_original_filter = true;
1019        return true;
1020      }
1021    }
1022
1023    // Evaluate the optimality error at the new point
1024    Number mu = .0;
1025    if (!IpData().FreeMuMode()) {
1026      mu = IpData().curr_mu();
1027    }
1028    Number trial_pderror;
1029    Number curr_pderror;
1030    try {
1031      trial_pderror = IpCq().trial_primal_dual_system_error(mu);
1032      curr_pderror = IpCq().curr_primal_dual_system_error(mu);
1033    }
1034    catch(IpoptNLP::Eval_Error& e) {
1035      e.ReportException(Jnlst(), J_DETAILED);
1036      Jnlst().Printf(J_WARNING, J_LINE_SEARCH,
1037                     "Warning: Evaluation error during soft restoration phase step.\n");
1038      IpData().Append_info_string("e");
1039      return false;
1040    }
1041
1042    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
1043                   "  Primal-dual error at current point:  %23.16e\n", curr_pderror);
1044    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
1045                   "  Primal-dual error at trial point  :  %23.16e\n", trial_pderror);
1046    // Check if there is sufficient reduction in the optimality error
1047    if (trial_pderror <= soft_resto_pderror_reduction_factor_*curr_pderror) {
1048      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
1049                     "  Trial step accepted.\n");
1050      return true;
1051    }
1052
1053    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
1054                   "  Trial step rejected.\n");
1055    return false;
1056  }
1057
1058  void FilterLineSearch::StartWatchDog()
1059  {
1060    DBG_START_FUN("FilterLineSearch::StartWatchDog", dbg_verbosity);
1061
1062    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
1063                   "Starting Watch Dog\n");
1064
1065    in_watchdog_ = true;
1066    watchdog_iterate_ = IpData().curr();
1067    watchdog_delta_ = IpData().delta();
1068    watchdog_trial_iter_ = 0;
1069    watchdog_alpha_primal_test_ =
1070      IpCq().curr_primal_frac_to_the_bound(IpData().curr_tau());
1071    watchdog_theta_ = IpCq().curr_constraint_violation();
1072    watchdog_barr_ = IpCq().curr_barrier_obj();
1073    watchdog_gradBarrTDelta_ = IpCq().curr_gradBarrTDelta();
1074  }
1075
1076  void FilterLineSearch::StopWatchDog(SmartPtr<IteratesVector>& actual_delta)
1077  {
1078    DBG_START_FUN("FilterLineSearch::StopWatchDog", dbg_verbosity);
1079
1080    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
1081                   "Stopping Watch Dog\n");
1082
1083    IpData().Append_info_string("w");
1084
1085    in_watchdog_ = false;
1086
1087    // Reset all fields in IpData to reference point
1088    SmartPtr<IteratesVector> old_trial = watchdog_iterate_->MakeNewContainer();
1089    IpData().set_trial(old_trial);
1090    IpData().AcceptTrialPoint();
1091    actual_delta = watchdog_delta_->MakeNewContainer();
1092    IpData().SetHaveAffineDeltas(false);
1093
1094    // reset the stored watchdog iterates
1095    watchdog_iterate_ = NULL;
1096    watchdog_delta_ = NULL;
1097
1098    watchdog_shortened_iter_ = 0;
1099
1100    reference_theta_ = watchdog_theta_;
1101    reference_barr_ = watchdog_barr_;
1102    reference_gradBarrTDelta_ = watchdog_gradBarrTDelta_;
1103  }
1104
1105  void FilterLineSearch::Reset()
1106  {
1107    DBG_START_FUN("FilterLineSearch::Reset", dbg_verbosity);
1108    in_soft_resto_phase_ = false;
1109
1110    // Inactivate the watchdog and release all stored data
1111    in_watchdog_ = false;
1112    watchdog_iterate_ = NULL;
1113    watchdog_delta_ = NULL;
1114    watchdog_shortened_iter_ = 0;
1115
1116    filter_.Clear();
1117  }
1118
1119  void FilterLineSearch::PerformDualStep(Number alpha_primal,
1120                                         Number alpha_dual,
1121                                         SmartPtr<IteratesVector>& delta)
1122  {
1123    DBG_START_FUN("FilterLineSearch::PerformDualStep", dbg_verbosity);
1124
1125    // set the bound multipliers from the step
1126    IpData().SetTrialBoundMultipliersFromStep(alpha_dual, *delta->z_L(), *delta->z_U(), *delta->v_L(), *delta->v_U());
1127
1128    Number alpha_y=-1.;
1129    switch (alpha_for_y_) {
1130      case PRIMAL_ALPHA_FOR_Y:
1131      alpha_y = alpha_primal;
1132      break;
1133      case DUAL_ALPHA_FOR_Y:
1134      alpha_y = alpha_dual;
1135      break;
1136      case MIN_ALPHA_FOR_Y:
1137      alpha_y = Min(alpha_dual, alpha_primal);
1138      break;
1139      case MAX_ALPHA_FOR_Y:
1140      alpha_y = Max(alpha_dual, alpha_primal);
1141      break;
1142      case FULL_STEP_FOR_Y:
1143      alpha_y = 1;
1144      break;
1145      case MIN_DUAL_INFEAS_ALPHA_FOR_Y:
1146      case SAFE_MIN_DUAL_INFEAS_ALPHA_FOR_Y:
1147      // Here we compute the step size for y so that the dual
1148      // infeasibility is minimized along delta_y
1149
1150      // compute the dual infeasibility at new point with old y
1151      SmartPtr<IteratesVector> temp_trial
1152      = IpData().trial()->MakeNewContainer();
1153      temp_trial->Set_y_c(*IpData().curr()->y_c());
1154      temp_trial->Set_y_d(*IpData().curr()->y_d());
1155      IpData().set_trial(temp_trial);
1156      SmartPtr<const Vector> dual_inf_x = IpCq().trial_grad_lag_x();
1157      SmartPtr<const Vector> dual_inf_s = IpCq().trial_grad_lag_s();
1158
1159      SmartPtr<Vector> new_jac_times_delta_y =
1160        IpData().curr()->x()->MakeNew();
1161      new_jac_times_delta_y->AddTwoVectors(1., *IpCq().trial_jac_cT_times_vec(*delta->y_c()),
1162                                           1., *IpCq().trial_jac_dT_times_vec(*delta->y_d()),
1163                                           0.);
1164
1165      Number a = pow(new_jac_times_delta_y->Nrm2(), 2.) +
1166                 pow(delta->y_d()->Nrm2(), 2.);
1167      Number b = dual_inf_x->Dot(*new_jac_times_delta_y) -
1168                 dual_inf_s->Dot(*delta->y_d());
1169
1170      Number alpha = - b/a;
1171
1172      if (alpha_for_y_==SAFE_MIN_DUAL_INFEAS_ALPHA_FOR_Y) {
1173        alpha_y = Min(Max(alpha_primal, alpha_dual), Max(alpha, Min(alpha_primal, alpha_dual)));
1174      }
1175      else {
1176        alpha_y = Min(1., Max(0., alpha));
1177      }
1178      break;
1179    }
1180
1181    // Set the eq multipliers from the step now that alpha_y
1182    // has been calculated.
1183    DBG_PRINT((1, "alpha_y = %e\n", alpha_y));
1184    DBG_PRINT_VECTOR(2, "delta_y_c", *delta->y_c());
1185    DBG_PRINT_VECTOR(2, "delta_y_d", *delta->y_d());
1186    IpData().SetTrialEqMultipliersFromStep(alpha_y, *delta->y_c(), *delta->y_d());
1187
1188    // Set some information for iteration summary output
1189    IpData().Set_info_alpha_primal(alpha_primal);
1190    IpData().Set_info_alpha_dual(alpha_dual);
1191  }
1192
1193  bool
1194  FilterLineSearch::TrySecondOrderCorrection(
1195    Number alpha_primal_test,
1196    Number& alpha_primal,
1197    SmartPtr<IteratesVector>& actual_delta)
1198  {
1199    DBG_START_METH("FilterLineSearch::TrySecondOrderCorrection",
1200                   dbg_verbosity);
1201
1202    bool accept = false;
1203    Index count_soc = 0;
1204
1205    Number theta_soc_old = 0.;
1206    Number theta_trial = IpCq().trial_constraint_violation();
1207    Number alpha_primal_soc = alpha_primal;
1208
1209    SmartPtr<Vector> c_soc = IpCq().curr_c()->MakeNew();
1210    SmartPtr<Vector> dms_soc = IpCq().curr_d_minus_s()->MakeNew();
1211    c_soc->Copy(*IpCq().curr_c());
1212    dms_soc->Copy(*IpCq().curr_d_minus_s());
1213    while (count_soc<max_soc_ && !accept &&
1214           (count_soc==0 || theta_trial<=kappa_soc_*theta_soc_old) ) {
1215      theta_soc_old = theta_trial;
1216
1217      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
1218                     "Trying second order correction number %d\n",
1219                     count_soc+1);
1220
1221      // Compute SOC constraint violation
1222      c_soc->AddOneVector(1.0, *IpCq().trial_c(), alpha_primal_soc);
1223      dms_soc->AddOneVector(1.0, *IpCq().trial_d_minus_s(), alpha_primal_soc);
1224
1225      // Compute the SOC search direction
1226      SmartPtr<IteratesVector> delta_soc = actual_delta->MakeNewIteratesVector(true);
1227      SmartPtr<IteratesVector> rhs = actual_delta->MakeNewContainer();
1228      rhs->Set_x(*IpCq().curr_grad_lag_with_damping_x());
1229      rhs->Set_s(*IpCq().curr_grad_lag_with_damping_s());
1230      rhs->Set_y_c(*c_soc);
1231      rhs->Set_y_d(*dms_soc);
1232      rhs->Set_z_L(*IpCq().curr_relaxed_compl_x_L());
1233      rhs->Set_z_U(*IpCq().curr_relaxed_compl_x_U());
1234      rhs->Set_v_L(*IpCq().curr_relaxed_compl_s_L());
1235      rhs->Set_v_U(*IpCq().curr_relaxed_compl_s_U());
1236      pd_solver_->Solve(-1.0, 0.0, *rhs, *delta_soc, true);
1237
1238      // Compute step size
1239      alpha_primal_soc =
1240        IpCq().primal_frac_to_the_bound(IpData().curr_tau(),
1241                                        *delta_soc->x(),
1242                                        *delta_soc->s());
1243
1244      // Check if trial point is acceptable
1245      try {
1246        // Compute the primal trial point
1247        IpData().SetTrialPrimalVariablesFromStep(alpha_primal_soc, *delta_soc->x(), *delta_soc->s());
1248
1249        // in acceptance tests, use original step size!
1250        accept = CheckAcceptabilityOfTrialPoint(alpha_primal_test);
1251      }
1252      catch(IpoptNLP::Eval_Error& e) {
1253        e.ReportException(Jnlst(), J_DETAILED);
1254        Jnlst().Printf(J_WARNING, J_MAIN, "Warning: SOC step rejected due to evaluation error\n");
1255        IpData().Append_info_string("e");
1256        accept = false;
1257        // There is no point in continuing SOC procedure
1258        break;
1259      }
1260
1261      if (accept) {
1262        Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Second order correction step accepted with %d corrections.\n", count_soc+1);
1263        // Accept all SOC quantities
1264        alpha_primal = alpha_primal_soc;
1265        actual_delta = delta_soc;
1266      }
1267      else {
1268        count_soc++;
1269        theta_trial = IpCq().trial_constraint_violation();
1270      }
1271    }
1272    return accept;
1273  }
1274
1275  bool
1276  FilterLineSearch::TryCorrector(
1277    Number alpha_primal_test,
1278    Number& alpha_primal,
1279    SmartPtr<IteratesVector>& actual_delta)
1280  {
1281    DBG_START_METH("FilterLineSearch::TryCorrector",
1282                   dbg_verbosity);
1283
1284    Index n_bounds = IpData().curr()->z_L()->Dim() + IpData().curr()->z_U()->Dim()
1285                     + IpData().curr()->v_L()->Dim() + IpData().curr()->v_U()->Dim();
1286    if (n_bounds==0) {
1287      // Nothing to be done
1288      return false;
1289    }
1290
1291    IpData().TimingStats().TryCorrector().Start();
1292
1293    bool accept = false;
1294
1295    // Compute the corrector step based on corrector_type parameter
1296    // create a new iterates vector and allocate space for all the entries
1297    SmartPtr<IteratesVector> delta_corr = actual_delta->MakeNewIteratesVector(true);
1298
1299    switch (corrector_type_) {
1300      case AFFINE_CORRECTOR : {
1301        // 1: Standard MPC corrector
1302
1303        if (!IpData().HaveAffineDeltas()) {
1304          Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
1305                         "Solving the Primal Dual System for the affine step\n");
1306          // First get the right hand side
1307          SmartPtr<IteratesVector> rhs_aff = delta_corr->MakeNewContainer();
1308
1309          rhs_aff->Set_x(*IpCq().curr_grad_lag_x());
1310          rhs_aff->Set_s(*IpCq().curr_grad_lag_s());
1311          rhs_aff->Set_y_c(*IpCq().curr_c());
1312          rhs_aff->Set_y_d(*IpCq().curr_d_minus_s());
1313          rhs_aff->Set_z_L(*IpCq().curr_compl_x_L());
1314          rhs_aff->Set_z_U(*IpCq().curr_compl_x_U());
1315          rhs_aff->Set_v_L(*IpCq().curr_compl_s_L());
1316          rhs_aff->Set_v_U(*IpCq().curr_compl_s_U());
1317
1318          // create a new iterates vector (with allocated space)
1319          // for the affine scaling step
1320          SmartPtr<IteratesVector> step_aff = delta_corr->MakeNewIteratesVector(true);
1321
1322          // Now solve the primal-dual system to get the step
1323          pd_solver_->Solve(-1.0, 0.0, *rhs_aff, *step_aff, false);
1324
1325          DBG_PRINT_VECTOR(2, "step_aff", *step_aff);
1326
1327          IpData().set_delta_aff(step_aff);
1328          IpData().SetHaveAffineDeltas(true);
1329        }
1330
1331        DBG_ASSERT(IpData().HaveAffineDeltas());
1332
1333        const SmartPtr<const IteratesVector> delta_aff = IpData().delta_aff();
1334
1335        delta_corr->Copy(*actual_delta);
1336
1337        // create a rhs vector and allocate entries
1338        SmartPtr<IteratesVector> rhs = actual_delta->MakeNewIteratesVector(true);
1339
1340        rhs->x_NonConst()->Set(0.);
1341        rhs->s_NonConst()->Set(0.);
1342        rhs->y_c_NonConst()->Set(0.);
1343        rhs->y_d_NonConst()->Set(0.);
1344        IpNLP().Px_L()->TransMultVector(-1., *delta_aff->x(), 0., *rhs->z_L_NonConst());
1345        rhs->z_L_NonConst()->ElementWiseMultiply(*delta_aff->z_L());
1346        IpNLP().Px_U()->TransMultVector(1., *delta_aff->x(), 0., *rhs->z_U_NonConst());
1347        rhs->z_U_NonConst()->ElementWiseMultiply(*delta_aff->z_U());
1348        IpNLP().Pd_L()->TransMultVector(-1., *delta_aff->s(), 0., *rhs->v_L_NonConst());
1349        rhs->v_L_NonConst()->ElementWiseMultiply(*delta_aff->v_L());
1350        IpNLP().Pd_U()->TransMultVector(1., *delta_aff->s(), 0., *rhs->v_U_NonConst());
1351        rhs->v_U_NonConst()->ElementWiseMultiply(*delta_aff->v_U());
1352
1353        pd_solver_->Solve(1.0, 1.0, *rhs, *delta_corr, true);
1354
1355        DBG_PRINT_VECTOR(2, "delta_corr", *delta_corr);
1356      }
1357      break;
1358      case PRIMAL_DUAL_CORRECTOR : {
1359        // 2: Second order correction for primal-dual step to
1360        // primal-dual mu
1361
1362        delta_corr->Copy(*actual_delta);
1363
1364        // allocate space for the rhs
1365        SmartPtr<IteratesVector> rhs = actual_delta->MakeNewIteratesVector(true);
1366
1367        rhs->x_NonConst()->Set(0.);
1368        rhs->s_NonConst()->Set(0.);
1369        rhs->y_c_NonConst()->Set(0.);
1370        rhs->y_d_NonConst()->Set(0.);
1371
1372        Number mu = IpData().curr_mu();
1373        SmartPtr<Vector> tmp;
1374
1375        rhs->z_L_NonConst()->Copy(*IpCq().curr_slack_x_L());
1376        IpNLP().Px_L()->TransMultVector(-1., *actual_delta->x(),
1377                                        -1., *rhs->z_L_NonConst());
1378        tmp = actual_delta->z_L()->MakeNew();
1379        tmp->AddTwoVectors(1., *IpData().curr()->z_L(), 1., *actual_delta->z_L(), 0.);
1380        rhs->z_L_NonConst()->ElementWiseMultiply(*tmp);
1381        rhs->z_L_NonConst()->AddScalar(mu);
1382
1383        rhs->z_U_NonConst()->Copy(*IpCq().curr_slack_x_U());
1384        IpNLP().Px_U()->TransMultVector(1., *actual_delta->x(),
1385                                        -1., *rhs->z_U_NonConst());
1386        tmp = actual_delta->z_U()->MakeNew();
1387        tmp->AddTwoVectors(1., *IpData().curr()->z_U(), 1., *actual_delta->z_U(), 0.);
1388        rhs->z_U_NonConst()->ElementWiseMultiply(*tmp);
1389        rhs->z_U_NonConst()->AddScalar(mu);
1390
1391        rhs->v_L_NonConst()->Copy(*IpCq().curr_slack_s_L());
1392        IpNLP().Pd_L()->TransMultVector(-1., *actual_delta->s(),
1393                                        -1., *rhs->v_L_NonConst());
1394        tmp = actual_delta->v_L()->MakeNew();
1395        tmp->AddTwoVectors(1., *IpData().curr()->v_L(), 1., *actual_delta->v_L(), 0.);
1396        rhs->v_L_NonConst()->ElementWiseMultiply(*tmp);
1397        rhs->v_L_NonConst()->AddScalar(mu);
1398
1399        rhs->v_U_NonConst()->Copy(*IpCq().curr_slack_s_U());
1400        IpNLP().Pd_U()->TransMultVector(1., *actual_delta->s(),
1401                                        -1., *rhs->v_U_NonConst());
1402        tmp = actual_delta->v_U()->MakeNew();
1403        tmp->AddTwoVectors(1., *IpData().curr()->v_U(), 1., *actual_delta->v_U(), 0.);
1404        rhs->v_U_NonConst()->ElementWiseMultiply(*tmp);
1405        rhs->v_U_NonConst()->AddScalar(mu);
1406
1407        DBG_PRINT_VECTOR(2, "rhs", *rhs);
1408
1409        pd_solver_->Solve(1.0, 1.0, *rhs, *delta_corr, true);
1410
1411        DBG_PRINT_VECTOR(2, "delta_corr", *delta_corr);
1412      }
1413      break;
1414      default:
1415      DBG_ASSERT("Unknown corrector_type value.");
1416    }
1417
1418    // Compute step size
1419    Number alpha_primal_corr =
1420      IpCq().primal_frac_to_the_bound(IpData().curr_tau(),
1421                                      *delta_corr->x(),
1422                                      *delta_corr->s());
1423    // Set the primal trial point
1424    IpData().SetTrialPrimalVariablesFromStep(alpha_primal_corr, *delta_corr->x(), *delta_corr->s());
1425
1426    // Check if we want to not even try the filter criterion
1427    Number alpha_dual_max =
1428      IpCq().dual_frac_to_the_bound(IpData().curr_tau(),
1429                                    *delta_corr->z_L(), *delta_corr->z_U(),
1430                                    *delta_corr->v_L(), *delta_corr->v_U());
1431
1432    IpData().SetTrialBoundMultipliersFromStep(alpha_dual_max, *delta_corr->z_L(), *delta_corr->z_U(), *delta_corr->v_L(), *delta_corr->v_U());
1433
1434    Number trial_avrg_compl = IpCq().trial_avrg_compl();
1435    Number curr_avrg_compl = IpCq().curr_avrg_compl();
1436    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
1437                   "avrg_compl(curr) = %e, avrg_compl(trial) = %e\n",
1438                   curr_avrg_compl, trial_avrg_compl);
1439    if (corrector_type_==AFFINE_CORRECTOR &&
1440        trial_avrg_compl>=corrector_compl_avrg_red_fact_*curr_avrg_compl) {
1441      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
1442                     "Rejecting corrector step, because trial complementarity is too large.\n" );
1443      IpData().TimingStats().TryCorrector().End();
1444      return false;
1445    }
1446
1447    // Check if trial point is acceptable
1448    try {
1449      // in acceptance tests, use original step size!
1450      accept = CheckAcceptabilityOfTrialPoint(alpha_primal_test);
1451    }
1452    catch(IpoptNLP::Eval_Error& e) {
1453      e.ReportException(Jnlst(), J_DETAILED);
1454      Jnlst().Printf(J_WARNING, J_MAIN,
1455                     "Warning: Corrector step rejected due to evaluation error\n");
1456      IpData().Append_info_string("e");
1457      accept = false;
1458    }
1459
1460    if (accept) {
1461      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
1462                     "Corrector step accepted with alpha_primal = %e\n",
1463                     alpha_primal_corr);
1464      // Accept all SOC quantities
1465      alpha_primal = alpha_primal_corr;
1466      actual_delta = delta_corr;
1467
1468      if (Jnlst().ProduceOutput(J_MOREVECTOR, J_MAIN)) {
1469        Jnlst().Printf(J_MOREVECTOR, J_MAIN,
1470                       "*** Accepted corrector for Iteration: %d\n",
1471                       IpData().iter_count());
1472        delta_corr->Print(Jnlst(), J_MOREVECTOR, J_MAIN, "delta_corr");
1473      }
1474    }
1475
1476    IpData().TimingStats().TryCorrector().End();
1477    return accept;
1478  }
1479
1480  void
1481  FilterLineSearch::PerformMagicStep()
1482  {
1483    DBG_START_METH("FilterLineSearch::PerformMagicStep",
1484                   2);//dbg_verbosity);
1485
1486    DBG_PRINT((1,"Incoming barr = %e and constrviol %e\n",
1487               IpCq().trial_barrier_obj(),
1488               IpCq().trial_constraint_violation()));
1489    DBG_PRINT_VECTOR(2, "s in", *IpData().trial()->s());
1490    DBG_PRINT_VECTOR(2, "d minus s in", *IpCq().trial_d_minus_s());
1491    DBG_PRINT_VECTOR(2, "slack_s_L in", *IpCq().trial_slack_s_L());
1492    DBG_PRINT_VECTOR(2, "slack_s_U in", *IpCq().trial_slack_s_U());
1493
1494    SmartPtr<const Vector> d_L = IpNLP().d_L();
1495    SmartPtr<const Matrix> Pd_L = IpNLP().Pd_L();
1496    SmartPtr<Vector> delta_s_magic_L = d_L->MakeNew();
1497    delta_s_magic_L->Set(0.);
1498    SmartPtr<Vector> tmp = d_L->MakeNew();
1499    Pd_L->TransMultVector(1., *IpCq().trial_d_minus_s(), 0., *tmp);
1500    delta_s_magic_L->ElementWiseMax(*tmp);
1501
1502    SmartPtr<const Vector> d_U = IpNLP().d_U();
1503    SmartPtr<const Matrix> Pd_U = IpNLP().Pd_U();
1504    SmartPtr<Vector> delta_s_magic_U = d_U->MakeNew();
1505    delta_s_magic_U->Set(0.);
1506    tmp = d_U->MakeNew();
1507    Pd_U->TransMultVector(1., *IpCq().trial_d_minus_s(), 0., *tmp);
1508    delta_s_magic_U->ElementWiseMin(*tmp);
1509
1510    SmartPtr<Vector> delta_s_magic = IpData().trial()->s()->MakeNew();
1511    Pd_L->MultVector(1., *delta_s_magic_L, 0., *delta_s_magic);
1512    Pd_U->MultVector(1., *delta_s_magic_U, 1., *delta_s_magic);
1513    delta_s_magic_L = NULL; // free memory
1514    delta_s_magic_U = NULL; // free memory
1515
1516    // Now find those entries with both lower and upper bounds, there
1517    // the step is too large
1518    // ToDo this should only be done if there are inequality
1519    // constraints with two bounds
1520    // also this can be done in a smaller space (d_L or d_U whichever
1521    // is smaller)
1522    tmp = delta_s_magic->MakeNew();
1523    tmp->Copy(*IpData().trial()->s());
1524    Pd_L->MultVector(1., *d_L, -2., *tmp);
1525    Pd_U->MultVector(1., *d_U, 1., *tmp);
1526    SmartPtr<Vector> tmp2 = tmp->MakeNew();
1527    tmp2->Copy(*tmp);
1528    tmp2->ElementWiseAbs();
1529    tmp->Axpy(-2., *delta_s_magic);
1530    tmp->ElementWiseAbs();
1531    // now, tmp2 = |d_L + d_u - 2*s| and tmp = |d_L + d_u - 2*(s+Delta s)|
1532    // we want to throw out those for which tmp2 > tmp
1533    tmp->Axpy(-1., *tmp2);
1534    tmp->ElementWiseSgn();
1535    tmp2->Set(0.);
1536    tmp2->ElementWiseMax(*tmp);
1537    tmp = d_L->MakeNew();
1538    Pd_L->TransMultVector(1., *tmp2, 0., *tmp);
1539    Pd_L->MultVector(1., *tmp, 0., *tmp2);
1540    tmp = d_U->MakeNew();
1541    Pd_U->TransMultVector(1., *tmp2, 0., *tmp);
1542    Pd_U->MultVector(1., *tmp, 0., *tmp2);
1543    DBG_PRINT_VECTOR(2, "tmp indicator", *tmp2)
1544    // tmp2 now is one for those entries with both bounds, for which
1545    // no step should be taken
1546
1547    tmp = delta_s_magic->MakeNew();
1548    tmp->Copy(*delta_s_magic);
1549    tmp->ElementWiseMultiply(*tmp2);
1550    delta_s_magic->Axpy(-1., *tmp);
1551
1552    Number delta_s_magic_max = delta_s_magic->Amax();
1553    Number mach_eps = std::numeric_limits<Number>::epsilon();
1554    if (delta_s_magic_max>0.) {
1555      if (delta_s_magic_max > 10*mach_eps*IpData().trial()->s()->Amax()) {
1556        IpData().Append_info_string("M");
1557        Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Magic step with max-norm %.6e taken.\n", delta_s_magic->Amax());
1558        delta_s_magic->Print(Jnlst(), J_MOREVECTOR, J_LINE_SEARCH,
1559                             "delta_s_magic");
1560      }
1561
1562      // now finally compute the new overall slacks
1563      delta_s_magic->Axpy(1., *IpData().trial()->s());
1564      SmartPtr<IteratesVector> trial = IpData().trial()->MakeNewContainer();
1565      trial->Set_s(*delta_s_magic);
1566
1567      // also update the set in the dual variables
1568
1569
1570      IpData().set_trial(trial);
1571    }
1572
1573    DBG_PRINT((1,"Outgoing barr = %e and constrviol %e\n", IpCq().trial_barrier_obj(), IpCq().trial_constraint_violation()));
1574    DBG_PRINT_VECTOR(2, "s out", *IpData().trial()->s());
1575    DBG_PRINT_VECTOR(2, "d minus s out", *IpCq().trial_d_minus_s());
1576    DBG_PRINT_VECTOR(2, "slack_s_L out", *IpCq().trial_slack_s_L());
1577    DBG_PRINT_VECTOR(2, "slack_s_U out", *IpCq().trial_slack_s_U());
1578  }
1579
1580  bool
1581  FilterLineSearch::DetectTinyStep()
1582  {
1583    DBG_START_METH("FilterLineSearch::DetectTinyStep",
1584                   dbg_verbosity);
1585
1586    Number max_step_x;
1587    Number max_step_s;
1588
1589    if (tiny_step_tol_==0.)
1590      return false;
1591
1592    // ToDo try to find more efficient implementation
1593
1594    SmartPtr<Vector> tmp = IpData().curr()->x()->MakeNew();
1595    tmp->Copy(*IpData().curr()->x());
1596    tmp->ElementWiseAbs();
1597    tmp->AddScalar(1.);
1598
1599    SmartPtr<Vector> tmp2 = IpData().curr()->x()->MakeNew();
1600    tmp2->Copy(*IpData().delta()->x());
1601    tmp2->ElementWiseDivide(*tmp);
1602    max_step_x = tmp2->Amax();
1603    Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
1604                   "Relative step size for delta_x = %e\n",
1605                   max_step_x);
1606    if (max_step_x > tiny_step_tol_)
1607      return false;
1608
1609    tmp = IpData().curr()->s()->MakeNew();
1610    tmp->Copy(*IpData().curr()->s());
1611    tmp->ElementWiseAbs();
1612    tmp->AddScalar(1.);
1613
1614    tmp2 = IpData().curr()->s()->MakeNew();
1615    tmp2->Copy(*IpData().delta()->s());
1616    tmp2->ElementWiseDivide(*tmp);
1617    max_step_s = tmp2->Amax();
1618    Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
1619                   "Relative step size for delta_s = %e\n",
1620                   max_step_s);
1621    if (max_step_s > tiny_step_tol_)
1622      return false;
1623
1624    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
1625                   "Tiny step of relative size %e detected.\n",
1626                   Max(max_step_x, max_step_s));
1627
1628    return true;
1629  }
1630
1631  bool FilterLineSearch::CurrentIsAcceptable()
1632  {
1633    return (IsValid(conv_check_) &&
1634            conv_check_->CurrentIsAcceptable());
1635  }
1636
1637  void FilterLineSearch::StoreAcceptablePoint()
1638  {
1639    DBG_START_METH("FilterLineSearch::StoreAcceptablePoint",
1640                   dbg_verbosity);
1641
1642    acceptable_iterate_ = IpData().curr();
1643  }
1644
1645  bool FilterLineSearch::RestoreAcceptablePoint()
1646  {
1647    DBG_START_METH("FilterLineSearch::RestoreAcceptablePoint",
1648                   dbg_verbosity);
1649
1650    if (!IsValid(acceptable_iterate_)) {
1651      return false;
1652    }
1653
1654    SmartPtr<IteratesVector> prev_iterate = acceptable_iterate_->MakeNewContainer();
1655    IpData().set_trial(prev_iterate);
1656    IpData().AcceptTrialPoint();
1657
1658    return true;
1659  }
1660
1661
1662} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.