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

Last change on this file since 412 was 412, checked in by andreasw, 14 years ago

added debug check to make sure direction is always indeed descent direction

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