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

Last change on this file since 529 was 529, checked in by andreasw, 14 years ago
  • in Vector class, copy cached values in Copy and update them in Scal
  • perform iterative refinement only once for adaptive strategy
  • fix bugs in PDPerturbationHandler
  • avoid some overhead in CalculateQualityFunction?
  • minor changes in timing
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 24.6 KB
Line 
1// Copyright (C) 2004, 2005 International Business Machines and others.
2// All Rights Reserved.
3// This code is published under the Common Public License.
4//
5// $Id: IpAdaptiveMuUpdate.cpp 529 2005-09-29 21:12:38Z andreasw $
6//
7// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
8
9#include "IpAdaptiveMuUpdate.hpp"
10#include "IpJournalist.hpp"
11
12#ifdef HAVE_CMATH
13# include <cmath>
14#else
15# ifdef HAVE_MATH_H
16#  include <math.h>
17# else
18#  error "don't have header file for math"
19# endif
20#endif
21
22namespace Ipopt
23{
24  // ToDo Make the different globalization strategies extra classes
25
26#ifdef IP_DEBUG
27  static const Index dbg_verbosity = 0;
28#endif
29
30  AdaptiveMuUpdate::AdaptiveMuUpdate
31  (const SmartPtr<LineSearch>& line_search,
32   const SmartPtr<MuOracle>& free_mu_oracle,
33   const SmartPtr<MuOracle>& fix_mu_oracle)
34      :
35      MuUpdate(),
36      linesearch_(line_search),
37      free_mu_oracle_(free_mu_oracle),
38      fix_mu_oracle_(fix_mu_oracle),
39      filter_(2)
40  {
41    DBG_ASSERT(IsValid(linesearch_));
42    DBG_ASSERT(IsValid(free_mu_oracle_));
43    // fix_mu_oracle may be NULL
44  }
45
46  AdaptiveMuUpdate::~AdaptiveMuUpdate()
47  {}
48
49  void AdaptiveMuUpdate::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
50  {
51    roptions->AddLowerBoundedNumberOption(
52      "mu_max",
53      "Maximum value for barrier parameter.",
54      0.0, true, 1e5,
55      "This option specifies an upper bound on the barrier parameter in the "
56      "adaptive mu selection mode.");
57    roptions->AddLowerBoundedNumberOption(
58      "mu_min",
59      "Minimum value for barrier parameter.",
60      0.0, true, 1e-9,
61      "This option specifies the lower bound on the barrier parameter in the "
62      "adaptive mu selection mode. By default, it is set to "
63      "min(\"tol\",\"compl_inf_tol\")/(\"barrier_tol_factor\"+1), which "
64      "should be a very reasonable value.");
65    std::string prev_cat = roptions->RegisteringCategory();
66    roptions->SetRegisteringCategory("Undocumented");
67    roptions->AddLowerBoundedNumberOption(
68      "adaptive_mu_safeguard_factor",
69      "",
70      0.0, false, 0.0);
71    roptions->SetRegisteringCategory(prev_cat);
72
73    roptions->AddStringOption3(
74      "adaptive_mu_globalization",
75      "Globalization strategy for the adaptive mu selection mode",
76      "obj-constr-filter",
77      "kkt-error", "nonmonotone decrease of kkt-error",
78      "obj-constr-filter", "2-dim filter for objective and constraint violation",
79      "never-monotone-mode", "disables globalization",
80      "To achieve global convergence of the adaptive version, the algorithm "
81      "has to switch to the monotone mode (Fiacco-McCormick approach) when "
82      "convergence does not seem to appear.  This option sets the "
83      "criterion used to decide when to do this switch.");
84
85    roptions->AddLowerBoundedIntegerOption(
86      "adaptive_mu_kkterror_red_iters",
87      "Maximum number of iterations requiring sufficient progress.",
88      0, 4,
89      "For the \"kkt-error\" based globalization strategy, sufficient "
90      "progress must be made for \"adaptive_mu_kkterror_red_iters\" "
91      "iterations. If this number of iterations is exceeded, the "
92      "globalization strategy switches to the monotone mode.");
93
94    roptions->AddBoundedNumberOption(
95      "adaptive_mu_kkterror_red_fact",
96      "Sufficient decrease factor for \"kkt-error\" globalization strategy.",
97      0.0, true, 1.0, true,
98      0.9999,
99      "For the \"kkt-error\" based globalization strategy, the error "
100      "must decrease by this factor to be deemed sufficient decrease.");
101
102    roptions->AddBoundedNumberOption(
103      "filter_margin_fact",
104      "Factor determining width of margin for obj-constr-filter adaptive globalization strategy.",
105      0.0, true, 1.0, true,
106      1e-5,
107      "When using the adaptive globalization strategy, \"obj-constr-filter\", "
108      "sufficient progress for a filter entry is defined as "
109      "follows: (new obj) < (filter obj) - filter_margin_fact*(new "
110      "constr-voil) OR (new constr-viol) < (filter constr-viol) - "
111      "filter_margin_fact*(new constr-voil).  For the description of "
112      "the \"kkt-error-filter\" option see \"filter_max_margin\".");
113    roptions->AddLowerBoundedNumberOption(
114      "filter_max_margin",
115      "Maximum width of margin in obj-constr--filter adaptive globalization strategy.",
116      0.0, true,
117      1.0,
118      // ToDo Detailed description later
119      "");
120    roptions->AddStringOption2(
121      "adaptive_mu_restore_previous_iterate",
122      "Should the previous iterate be restored if the monotone mode is entered.",
123      "no",
124      "no", "don't restore accepted iterate",
125      "yes", "restore accepted iterate",
126      "When the globalization strategy for the adaptive barrier algorithm "
127      "switches to the monotone mode, it can either start "
128      "from the most recent iterate (no), or from the last "
129      "iterate that was accepted (yes).");
130
131    roptions->AddLowerBoundedNumberOption(
132      "adaptive_mu_monotone_init_factor",
133      "Determines the initial value of the barrier parameter when switching to the monotone mode.",
134      0.0, true, 0.8,
135      "When the globalization strategy for the adaptive barrier algorithm "
136      "switches to the monotone mode and fixed_mu_oracle is chosen as "
137      "\"average_compl\", the barrier parameter is set to the "
138      "current average complementarity times the value of "
139      "\"adaptive_mu_monotone_init_factor\".");
140
141    roptions->AddStringOption4(
142      "adaptive_mu_kkt_norm_type",
143      "Norm used for the KKT error in the adaptive mu globalization strategies.",
144      "2-norm-squared",
145      "1-norm", "use the 1-norm (abs sum)",
146      "2-norm-squared", "use the 2-norm squared (sum of squares)",
147      "max-norm", "use the infinity norm (max)",
148      "2-norm", "use 2-norm",
149      "When computing the KKT error for the globalization strategies, the "
150      "norm to be used is specified with this option. Note, this options is also used "
151      "in the QualityFunctionMuOracle.");
152
153  }
154
155  bool AdaptiveMuUpdate::InitializeImpl(const OptionsList& options,
156                                        const std::string& prefix)
157  {
158    options.GetNumericValue("mu_max", mu_max_, prefix);
159    options.GetNumericValue("tau_min", tau_min_, prefix);
160    options.GetNumericValue("adaptive_mu_safeguard_factor", adaptive_mu_safeguard_factor_, prefix);
161    options.GetNumericValue("adaptive_mu_kkterror_red_fact", refs_red_fact_, prefix);
162    options.GetIntegerValue("adaptive_mu_kkterror_red_iters", num_refs_max_, prefix);
163    Index enum_int;
164    options.GetEnumValue("adaptive_mu_globalization", enum_int, prefix);
165    adaptive_mu_globalization_ = AdaptiveMuGlobalizationEnum(enum_int);
166    options.GetNumericValue("filter_max_margin", filter_max_margin_, prefix);
167    options.GetNumericValue("filter_margin_fact", filter_margin_fact_, prefix);
168    options.GetBoolValue("adaptive_mu_restore_previous_iterate", restore_accepted_iterate_, prefix);
169
170    bool retvalue = free_mu_oracle_->Initialize(Jnlst(), IpNLP(), IpData(),
171                    IpCq(), options, prefix);
172    if (!retvalue) {
173      return retvalue;
174    }
175
176    if (IsValid(fix_mu_oracle_)) {
177      retvalue = fix_mu_oracle_->Initialize(Jnlst(), IpNLP(), IpData(),
178                                            IpCq(), options, prefix);
179      if (!retvalue) {
180        return retvalue;
181      }
182    }
183
184    options.GetNumericValue("adaptive_mu_monotone_init_factor", adaptive_mu_monotone_init_factor_, prefix);
185    options.GetNumericValue("barrier_tol_factor", barrier_tol_factor_, prefix);
186    options.GetNumericValue("mu_linear_decrease_factor", mu_linear_decrease_factor_, prefix);
187    options.GetNumericValue("mu_superlinear_decrease_power", mu_superlinear_decrease_power_, prefix);
188
189    options.GetEnumValue("quality_function_norm_type", enum_int, prefix);
190    adaptive_mu_kkt_norm_ = QualityFunctionMuOracle::NormEnum(enum_int);
191    options.GetEnumValue("quality_function_centrality", enum_int, prefix);
192    adaptive_mu_kkt_centrality_ = QualityFunctionMuOracle::CentralityEnum(enum_int);
193    options.GetEnumValue("quality_function_balancing_term", enum_int, prefix);
194    adaptive_mu_kkt_balancing_term_ = QualityFunctionMuOracle::BalancingTermEnum(enum_int);
195    options.GetNumericValue("compl_inf_tol", compl_inf_tol_, prefix);
196    if (!options.GetNumericValue("mu_min", mu_min_, prefix)) {
197      // Defer computation of the default until the scaling of the NLP
198      // is known
199      mu_min_ = -1.;
200    }
201
202    init_dual_inf_ = -1.;
203    init_primal_inf_ = -1.;
204
205    refs_vals_.clear();
206    check_if_no_bounds_ = false;
207    no_bounds_ = false;
208    filter_.Clear();
209    IpData().SetFreeMuMode(true);
210
211    accepted_point_ = NULL;
212
213    // The following lines are only here so that
214    // IpoptCalculatedQuantities::CalculateSafeSlack and the first
215    // output line have something to work with
216    IpData().Set_mu(1.);
217    IpData().Set_tau(0.);
218
219    return retvalue;
220  }
221
222  void AdaptiveMuUpdate::UpdateBarrierParameter()
223  {
224    DBG_START_METH("AdaptiveMuUpdate::UpdateBarrierParameter",
225                   dbg_verbosity);
226
227    // if min_mu_ has not been given, we now set the default (can't do
228    // that earlier, because during call of InitializeImpl, the
229    // scaling in the NLP is not yet determined
230    if (mu_min_ < 0.) {
231      mu_min_ = Min(IpData().tol(),
232                    IpNLP().NLP_scaling()->apply_obj_scaling(compl_inf_tol_))/
233                (barrier_tol_factor_+1.);
234    }
235
236    // if there are not bounds, we always return the minimum MU value
237    if (!check_if_no_bounds_) {
238      Index n_bounds = IpData().curr()->z_L()->Dim() + IpData().curr()->z_U()->Dim()
239                       + IpData().curr()->v_L()->Dim() + IpData().curr()->v_U()->Dim();
240
241      if (n_bounds==0) {
242        no_bounds_ = true;
243        IpData().Set_mu(mu_min_);
244        IpData().Set_tau(tau_min_);
245      }
246
247      check_if_no_bounds_ = true;
248    }
249
250    if (no_bounds_)
251      return;
252
253    bool tiny_step_flag = IpData().tiny_step_flag();
254    if (!IpData().FreeMuMode()) {
255      // if we are in the fixed mu mode, we need to check if the
256      // current iterate is good enough to continue with the free mode
257      bool sufficient_progress = CheckSufficientProgress();
258      if (sufficient_progress && !tiny_step_flag) {
259        Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
260                       "Switching back to free mu mode.\n");
261        IpData().SetFreeMuMode(true);
262        // Skipping Restoration phase?
263        RememberCurrentPointAsAccepted();
264      }
265      else {
266        Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
267                       "Remaining in fixed mu mode.\n");
268
269        // ToDo decide whether we want this for all options
270        Number sub_problem_error = IpCq().curr_barrier_error();
271        Number mu = IpData().curr_mu();
272        if (sub_problem_error <= barrier_tol_factor_ * mu ||
273            tiny_step_flag) {
274          // If the current barrier problem has been solved sufficiently
275          // well, decrease mu
276          // ToDo combine this code with MonotoneMuUpdate
277          Number tol = IpData().tol();
278          Number compl_inf_tol =
279            IpNLP().NLP_scaling()->apply_obj_scaling(compl_inf_tol_);
280
281          Number new_mu = Min( mu_linear_decrease_factor_*mu,
282                               pow(mu, mu_superlinear_decrease_power_) );
283          DBG_PRINT((1,"new_mu = %e, compl_inf_tol = %e tol = %e\n", new_mu, compl_inf_tol, tol));
284          new_mu = Max(new_mu,
285                       Min(compl_inf_tol, tol)/(barrier_tol_factor_+1.));
286          if (tiny_step_flag && new_mu == mu) {
287            THROW_EXCEPTION(TINY_STEP_DETECTED,
288                            "Problem solved to best possible numerical accuracy");
289          }
290          Number new_tau = Compute_tau_monotone(mu);
291          IpData().Set_mu(new_mu);
292          IpData().Set_tau(new_tau);
293          Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
294                         "Reducing mu to %24.16e in fixed mu mode. Tau becomes %24.16e\n", new_mu, new_tau);
295          linesearch_->Reset();
296        }
297      }
298    }
299    else {
300      // Here we are in the free mu mode.
301      bool sufficient_progress = CheckSufficientProgress();
302      if (linesearch_->CheckSkippedLineSearch() || tiny_step_flag ) {
303        sufficient_progress = false;
304      }
305      if (sufficient_progress) {
306        Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
307                       "Staying in free mu mode.\n");
308        RememberCurrentPointAsAccepted();
309      }
310      else {
311        IpData().SetFreeMuMode(false);
312
313        if (restore_accepted_iterate_) {
314          // Restore most recent accepted iterate to start fixed mode from
315          Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
316                         "Restoring most recent accepted point.\n");
317          SmartPtr<IteratesVector> prev_iter = accepted_point_->MakeNewContainer();
318          IpData().set_trial(prev_iter);
319          IpData().AcceptTrialPoint();
320        }
321
322        // Set the new values for mu and tau and tell the linesearch
323        // to reset its memory
324        Number mu = NewFixedMu();
325        Number tau = Compute_tau_monotone(mu);
326
327        if (tiny_step_flag && mu==IpData().curr_mu()) {
328          THROW_EXCEPTION(TINY_STEP_DETECTED,
329                          "Problem solved to best possible numerical accuracy");
330        }
331
332        IpData().Set_mu(mu);
333        IpData().Set_tau(tau);
334        Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
335                       "Switching to fixed mu mode with mu = %24.16e and tau = %24.16e.\n", mu, tau);
336        linesearch_->Reset();
337        // Skipping Restoration phase?
338      }
339    }
340
341    if (IpData().FreeMuMode()) {
342
343      // Choose the fraction-to-the-boundary parameter for the current
344      // iteration
345      // ToDo: Is curr_nlp_error really what we should use here?
346      Number tau = Max(tau_min_, 1.-IpCq().curr_nlp_error());
347      IpData().Set_tau(tau);
348
349      // Compute the new barrier parameter via the oracle
350      Number mu = free_mu_oracle_->CalculateMu(mu_min_, mu_max_);
351
352      mu = Max(mu, mu_min_);
353      Number mu_lower_safe = lower_mu_safeguard();
354      if (mu < mu_lower_safe) {
355        Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
356                       "mu = %e smaller than safeguard = %e. Increasing mu.\n",
357                       mu, mu_lower_safe);
358        mu = mu_lower_safe;
359        IpData().Append_info_string("m");
360      }
361
362      Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
363                     "Barrier parameter mu computed by oracle is %e\n",
364                     mu);
365
366      // Apply safeguards if appropriate
367      mu = Min(mu, mu_max_);
368      Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
369                     "Barrier parameter mu after safeguards is %e\n",
370                     mu);
371
372      // Set the new values
373      IpData().Set_mu(mu);
374
375      linesearch_->Reset();
376      // Uncomment the next line if the filter should not switch to
377      // the restoration phase in the free mode
378      // linesearch_->SetRigorousLineSearch(false);
379    }
380    else {
381      IpData().Append_info_string("F");
382      linesearch_->SetRigorousLineSearch(true);
383    }
384  }
385
386  bool
387  AdaptiveMuUpdate::CheckSufficientProgress()
388  {
389    bool retval = true;
390
391    switch (adaptive_mu_globalization_) {
392      case KKT_ERROR : {
393        Index num_refs = (Index)refs_vals_.size();
394        if (num_refs >= num_refs_max_) {
395          retval = false;
396          Number curr_error = quality_function_pd_system();
397          std::list<Number>::iterator iter;
398          for (iter = refs_vals_.begin(); iter != refs_vals_.end();
399               iter++) {
400            if ( curr_error <= refs_red_fact_*(*iter) ) {
401              retval = true;
402            }
403          }
404        }
405      }
406      break;
407      case FILTER_OBJ_CONSTR : {
408        /*
409               retval = filter_.Acceptable(IpCq().curr_f(),
410                                           IpCq().curr_constraint_violation());
411        */
412        // ToDo: Is curr_nlp_error really what we should use here?
413        Number curr_error = IpCq().curr_nlp_error();
414        Number margin = filter_margin_fact_*Min(filter_max_margin_, curr_error);
415        retval = filter_.Acceptable(IpCq().curr_f() + margin,
416                                    IpCq().curr_constraint_violation() + margin);
417      }
418      break;
419      case FILTER_KKT_ERROR : {
420        DBG_ASSERT("Unknown adaptive_mu_globalization value.");
421      }
422      break;
423      case NEVER_MONOTONE_MODE :
424      retval = true;
425      break;
426      default:
427      DBG_ASSERT("Unknown adaptive_mu_globalization value.");
428    }
429
430    return retval;
431  }
432
433  void
434  AdaptiveMuUpdate::RememberCurrentPointAsAccepted()
435  {
436    switch (adaptive_mu_globalization_) {
437      case KKT_ERROR : {
438        Number curr_error = quality_function_pd_system();
439        Index num_refs = (Index)refs_vals_.size();
440        if (num_refs >= num_refs_max_) {
441          refs_vals_.pop_front();
442        }
443        refs_vals_.push_back(curr_error);
444
445        if (Jnlst().ProduceOutput(J_MOREDETAILED, J_BARRIER_UPDATE)) {
446          Index num_refs = 0;
447          std::list<Number>::iterator iter;
448          for (iter = refs_vals_.begin(); iter != refs_vals_.end();
449               iter++) {
450            num_refs++;
451            Jnlst().Printf(J_MOREDETAILED, J_BARRIER_UPDATE,
452                           "pd system reference[%2d] = %.6e\n", num_refs, *iter);
453          }
454        }
455      }
456      break;
457      case FILTER_OBJ_CONSTR : {
458        /*
459               Number theta = IpCq().curr_constraint_violation();
460               filter_.AddEntry(IpCq().curr_f() - filter_margin_fact_*theta,
461                                IpCq().curr_constraint_violation() - filter_margin_fact_*theta,
462                                IpData().iter_count());
463               filter_.Print(Jnlst());
464        */
465        filter_.AddEntry(IpCq().curr_f(),
466                         IpCq().curr_constraint_violation(),
467                         IpData().iter_count());
468        filter_.Print(Jnlst());
469      }
470      break;
471      case FILTER_KKT_ERROR : {
472        DBG_ASSERT("Unknown corrector_type value.");
473      }
474      break;
475      default:
476      DBG_ASSERT("Unknown corrector_type value.");
477    }
478
479    if (restore_accepted_iterate_) {
480      // Keep pointers to this iterate so that it could be restored
481      accepted_point_ = IpData().curr();
482    }
483  }
484
485  Number
486  AdaptiveMuUpdate::Compute_tau_monotone(Number mu)
487  {
488    return Max(tau_min_, 1.-mu);
489  }
490
491  Number
492  AdaptiveMuUpdate::min_ref_val()
493  {
494    DBG_ASSERT(adaptive_mu_globalization_==KKT_ERROR);
495    Number min_ref;
496    DBG_ASSERT(refs_vals_.size()>0);
497    std::list<Number>::iterator iter = refs_vals_.begin();
498    min_ref = *iter;
499    iter++;
500    while (iter != refs_vals_.end()) {
501      min_ref = Min(min_ref, *iter);
502      iter++;
503    }
504    return min_ref;
505  }
506
507  Number
508  AdaptiveMuUpdate::max_ref_val()
509  {
510    DBG_ASSERT(adaptive_mu_globalization_==KKT_ERROR);
511    Number max_ref;
512    DBG_ASSERT(refs_vals_.size()>0);
513    std::list<Number>::iterator iter = refs_vals_.begin();
514    max_ref = *iter;
515    iter++;
516    while (iter != refs_vals_.end()) {
517      max_ref = Max(max_ref, *iter);
518      iter++;
519    }
520    return max_ref;
521  }
522
523  Number
524  AdaptiveMuUpdate::NewFixedMu()
525  {
526    Number max_ref;
527    // ToDo: Decide whether we should impose an upper bound on
528    // mu based on the smallest reference value.  For now, don't
529    // impose one.
530    max_ref = 1e20;
531    /*
532    switch (adaptive_mu_globalization_) {
533      case 1 :
534      max_ref = max_ref_val();
535      break;
536      case 2 : {
537        max_ref = 1e20;
538      }
539      break;
540      default:
541      DBG_ASSERT("Unknown corrector_type value.");
542    }
543    */
544
545    Number new_mu;
546
547    if (IsValid(fix_mu_oracle_)) {
548      new_mu = fix_mu_oracle_->CalculateMu(mu_min_, mu_max_);
549    }
550    else {
551      new_mu = adaptive_mu_monotone_init_factor_*IpCq().curr_avrg_compl();
552    }
553    new_mu = Max(new_mu, lower_mu_safeguard());
554    new_mu = Min(new_mu, 0.1 * max_ref);
555
556    new_mu = Max(new_mu, mu_min_);
557    new_mu = Min(new_mu, mu_max_);
558
559    return new_mu;
560  }
561
562  Number
563  AdaptiveMuUpdate::quality_function_pd_system()
564  {
565    Index n_dual = IpData().curr()->x()->Dim() + IpData().curr()->s()->Dim();
566    Index n_pri = IpData().curr()->y_c()->Dim() + IpData().curr()->y_d()->Dim();
567    Index n_comp = IpData().curr()->z_L()->Dim() + IpData().curr()->z_U()->Dim() +
568                   IpData().curr()->v_L()->Dim() + IpData().curr()->v_U()->Dim();
569
570    Number dual_inf=0.;
571    Number primal_inf=0.;
572    Number complty=0.;
573    switch (adaptive_mu_kkt_norm_) {
574      case QualityFunctionMuOracle::NM_NORM_1:
575      dual_inf =
576        IpCq().curr_dual_infeasibility(NORM_1);
577      primal_inf =
578        IpCq().curr_primal_infeasibility(NORM_1);
579      complty =
580        IpCq().curr_complementarity(0., NORM_1);
581      dual_inf /= (Number)n_dual;
582      DBG_ASSERT(n_pri>0 || primal_inf==0.);
583      if (n_pri>0) {
584        primal_inf /= (Number)n_pri;
585      }
586      DBG_ASSERT(n_comp>0 || complty==0.);
587      if (n_comp>0) {
588        complty /= (Number)n_comp;
589      }
590      break;
591      case QualityFunctionMuOracle::NM_NORM_2_SQUARED:
592      dual_inf =
593        IpCq().curr_dual_infeasibility(NORM_2);
594      dual_inf *= dual_inf;
595      primal_inf =
596        IpCq().curr_primal_infeasibility(NORM_2);
597      primal_inf *= primal_inf;
598      complty =
599        IpCq().curr_complementarity(0., NORM_2);
600      complty *= complty;
601      dual_inf /= (Number)n_dual;
602      DBG_ASSERT(n_pri>0 || primal_inf==0.);
603      if (n_pri>0) {
604        primal_inf /= (Number)n_pri;
605      }
606      DBG_ASSERT(n_comp>0 || complty==0.);
607      if (n_comp>0) {
608        complty /= (Number)n_comp;
609      }
610      break;
611      case QualityFunctionMuOracle::NM_NORM_MAX:
612      dual_inf =
613        IpCq().curr_dual_infeasibility(NORM_MAX);
614      primal_inf =
615        IpCq().curr_primal_infeasibility(NORM_MAX);
616      complty =
617        IpCq().curr_complementarity(0., NORM_MAX);
618      break;
619      case QualityFunctionMuOracle::NM_NORM_2:
620      dual_inf =
621        IpCq().curr_dual_infeasibility(NORM_2);
622      primal_inf =
623        IpCq().curr_primal_infeasibility(NORM_2);
624      complty =
625        IpCq().curr_complementarity(0., NORM_2);
626      dual_inf /= sqrt((Number)n_dual);
627      DBG_ASSERT(n_pri>0 || primal_inf==0.);
628      if (n_pri>0) {
629        primal_inf /= sqrt((Number)n_pri);
630      }
631      DBG_ASSERT(n_comp>0 || complty==0.);
632      if (n_comp>0) {
633        complty /= sqrt((Number)n_comp);
634      }
635      break;
636    }
637
638    Number centrality = 0.;
639    if (adaptive_mu_kkt_centrality_!=0) {
640      Number xi = IpCq().curr_centrality_measure();
641      switch (adaptive_mu_kkt_centrality_) {
642        case 1:
643        centrality = -complty*log(xi);
644        break;
645        case 2:
646        centrality = complty/xi;
647        case 3:
648        centrality = complty/pow(xi,3);
649        break;
650        default:
651        DBG_ASSERT("Unknown value for adaptive_mu_kkt_centrality_");
652      }
653    }
654
655    Number balancing_term=0.;
656    switch (adaptive_mu_kkt_balancing_term_) {
657      case 0:
658      //Nothing
659      break;
660      case 1:
661      balancing_term = pow(Max(0., Max(dual_inf,primal_inf)-complty),3);
662      break;
663      default:
664      DBG_ASSERT("Unknown value for adaptive_mu_kkt_balancing_term");
665    }
666
667    DBG_ASSERT(centrality>=0.);
668    DBG_ASSERT(balancing_term>=0);
669    Number kkt_error = primal_inf + dual_inf + complty +
670                       centrality + balancing_term;
671
672    Jnlst().Printf(J_MOREDETAILED, J_BARRIER_UPDATE,
673                   "KKT error in barrier update check:\n"
674                   "  primal infeasibility: %15.6e\n"
675                   "    dual infeasibility: %15.6e\n"
676                   "       complementarity: %15.6e\n"
677                   "            centrality: %15.6e\n"
678                   "             kkt error: %15.6e\n",
679                   primal_inf, dual_inf, complty, centrality, kkt_error);
680
681    return kkt_error;
682  }
683
684  Number
685  AdaptiveMuUpdate::lower_mu_safeguard()
686  {
687    DBG_START_METH("AdaptiveMuUpdate::lower_mu_safeguard",
688                   dbg_verbosity);
689    if (adaptive_mu_safeguard_factor_ == 0.)
690      return 0.;
691
692    Number dual_inf =
693      IpCq().curr_dual_infeasibility(NORM_1);
694    Number primal_inf =
695      IpCq().curr_primal_infeasibility(NORM_1);
696    Index n_dual = IpData().curr()->x()->Dim() + IpData().curr()->s()->Dim();
697    dual_inf /= (Number)n_dual;
698    Index n_pri = IpData().curr()->y_c()->Dim() + IpData().curr()->y_d()->Dim();
699    DBG_ASSERT(n_pri>0 || primal_inf==0.);
700    if (n_pri>0) {
701      primal_inf /= (Number)n_pri;
702    }
703
704    if (init_dual_inf_ < 0.) {
705      init_dual_inf_ = Max(1., dual_inf);
706    }
707    if (init_primal_inf_ < 0.) {
708      init_primal_inf_ = Max(1., primal_inf);
709    }
710
711    Number lower_mu_safeguard =
712      Max(adaptive_mu_safeguard_factor_ * (dual_inf/init_dual_inf_),
713          adaptive_mu_safeguard_factor_ * (primal_inf/init_primal_inf_));
714    DBG_PRINT((1,"dual_inf=%e init_dual_inf_=%e primal_inf=%e init_primal_inf_=%e\n", dual_inf, init_dual_inf_, primal_inf, init_primal_inf_));
715
716    if (adaptive_mu_globalization_==KKT_ERROR) {
717      lower_mu_safeguard = Min(lower_mu_safeguard, min_ref_val());
718    }
719
720    return lower_mu_safeguard;
721  }
722
723} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.