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

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