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

Last change on this file since 542 was 542, checked in by andreasw, 14 years ago
  • cleaned up line search to allow for alternative globalization scheme
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 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 542 2005-10-13 22:43:08Z 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 line search should not switch to
377      // the restoration phase in the free mode
378
379      // linesearch_->SetRigorousLineSearch(false);
380    }
381    else {
382      IpData().Append_info_string("F");
383      linesearch_->SetRigorousLineSearch(true);
384    }
385  }
386
387  bool
388  AdaptiveMuUpdate::CheckSufficientProgress()
389  {
390    bool retval = true;
391
392    switch (adaptive_mu_globalization_) {
393      case KKT_ERROR : {
394        Index num_refs = (Index)refs_vals_.size();
395        if (num_refs >= num_refs_max_) {
396          retval = false;
397          Number curr_error = quality_function_pd_system();
398          std::list<Number>::iterator iter;
399          for (iter = refs_vals_.begin(); iter != refs_vals_.end();
400               iter++) {
401            if ( curr_error <= refs_red_fact_*(*iter) ) {
402              retval = true;
403            }
404          }
405        }
406      }
407      break;
408      case FILTER_OBJ_CONSTR : {
409        /*
410               retval = filter_.Acceptable(IpCq().curr_f(),
411                                           IpCq().curr_constraint_violation());
412        */
413        // ToDo: Is curr_nlp_error really what we should use here?
414        Number curr_error = IpCq().curr_nlp_error();
415        Number margin = filter_margin_fact_*Min(filter_max_margin_, curr_error);
416        retval = filter_.Acceptable(IpCq().curr_f() + margin,
417                                    IpCq().curr_constraint_violation() + margin);
418      }
419      break;
420      case FILTER_KKT_ERROR : {
421        DBG_ASSERT("Unknown adaptive_mu_globalization value.");
422      }
423      break;
424      case NEVER_MONOTONE_MODE :
425      retval = true;
426      break;
427      default:
428      DBG_ASSERT("Unknown adaptive_mu_globalization value.");
429    }
430
431    return retval;
432  }
433
434  void
435  AdaptiveMuUpdate::RememberCurrentPointAsAccepted()
436  {
437    switch (adaptive_mu_globalization_) {
438      case KKT_ERROR : {
439        Number curr_error = quality_function_pd_system();
440        Index num_refs = (Index)refs_vals_.size();
441        if (num_refs >= num_refs_max_) {
442          refs_vals_.pop_front();
443        }
444        refs_vals_.push_back(curr_error);
445
446        if (Jnlst().ProduceOutput(J_MOREDETAILED, J_BARRIER_UPDATE)) {
447          Index num_refs = 0;
448          std::list<Number>::iterator iter;
449          for (iter = refs_vals_.begin(); iter != refs_vals_.end();
450               iter++) {
451            num_refs++;
452            Jnlst().Printf(J_MOREDETAILED, J_BARRIER_UPDATE,
453                           "pd system reference[%2d] = %.6e\n", num_refs, *iter);
454          }
455        }
456      }
457      break;
458      case FILTER_OBJ_CONSTR : {
459        /*
460               Number theta = IpCq().curr_constraint_violation();
461               filter_.AddEntry(IpCq().curr_f() - filter_margin_fact_*theta,
462                                IpCq().curr_constraint_violation() - filter_margin_fact_*theta,
463                                IpData().iter_count());
464               filter_.Print(Jnlst());
465        */
466        filter_.AddEntry(IpCq().curr_f(),
467                         IpCq().curr_constraint_violation(),
468                         IpData().iter_count());
469        filter_.Print(Jnlst());
470      }
471      break;
472      case FILTER_KKT_ERROR : {
473        DBG_ASSERT("Unknown corrector_type value.");
474      }
475      break;
476      default:
477      DBG_ASSERT("Unknown corrector_type value.");
478    }
479
480    if (restore_accepted_iterate_) {
481      // Keep pointers to this iterate so that it could be restored
482      accepted_point_ = IpData().curr();
483    }
484  }
485
486  Number
487  AdaptiveMuUpdate::Compute_tau_monotone(Number mu)
488  {
489    return Max(tau_min_, 1.-mu);
490  }
491
492  Number
493  AdaptiveMuUpdate::min_ref_val()
494  {
495    DBG_ASSERT(adaptive_mu_globalization_==KKT_ERROR);
496    Number min_ref;
497    DBG_ASSERT(refs_vals_.size()>0);
498    std::list<Number>::iterator iter = refs_vals_.begin();
499    min_ref = *iter;
500    iter++;
501    while (iter != refs_vals_.end()) {
502      min_ref = Min(min_ref, *iter);
503      iter++;
504    }
505    return min_ref;
506  }
507
508  Number
509  AdaptiveMuUpdate::max_ref_val()
510  {
511    DBG_ASSERT(adaptive_mu_globalization_==KKT_ERROR);
512    Number max_ref;
513    DBG_ASSERT(refs_vals_.size()>0);
514    std::list<Number>::iterator iter = refs_vals_.begin();
515    max_ref = *iter;
516    iter++;
517    while (iter != refs_vals_.end()) {
518      max_ref = Max(max_ref, *iter);
519      iter++;
520    }
521    return max_ref;
522  }
523
524  Number
525  AdaptiveMuUpdate::NewFixedMu()
526  {
527    Number max_ref;
528    // ToDo: Decide whether we should impose an upper bound on
529    // mu based on the smallest reference value.  For now, don't
530    // impose one.
531    max_ref = 1e20;
532    /*
533    switch (adaptive_mu_globalization_) {
534      case 1 :
535      max_ref = max_ref_val();
536      break;
537      case 2 : {
538        max_ref = 1e20;
539      }
540      break;
541      default:
542      DBG_ASSERT("Unknown corrector_type value.");
543    }
544    */
545
546    Number new_mu;
547
548    if (IsValid(fix_mu_oracle_)) {
549      new_mu = fix_mu_oracle_->CalculateMu(mu_min_, mu_max_);
550    }
551    else {
552      new_mu = adaptive_mu_monotone_init_factor_*IpCq().curr_avrg_compl();
553    }
554    new_mu = Max(new_mu, lower_mu_safeguard());
555    new_mu = Min(new_mu, 0.1 * max_ref);
556
557    new_mu = Max(new_mu, mu_min_);
558    new_mu = Min(new_mu, mu_max_);
559
560    return new_mu;
561  }
562
563  Number
564  AdaptiveMuUpdate::quality_function_pd_system()
565  {
566    Index n_dual = IpData().curr()->x()->Dim() + IpData().curr()->s()->Dim();
567    Index n_pri = IpData().curr()->y_c()->Dim() + IpData().curr()->y_d()->Dim();
568    Index n_comp = IpData().curr()->z_L()->Dim() + IpData().curr()->z_U()->Dim() +
569                   IpData().curr()->v_L()->Dim() + IpData().curr()->v_U()->Dim();
570
571    Number dual_inf=0.;
572    Number primal_inf=0.;
573    Number complty=0.;
574    switch (adaptive_mu_kkt_norm_) {
575      case QualityFunctionMuOracle::NM_NORM_1:
576      dual_inf =
577        IpCq().curr_dual_infeasibility(NORM_1);
578      primal_inf =
579        IpCq().curr_primal_infeasibility(NORM_1);
580      complty =
581        IpCq().curr_complementarity(0., NORM_1);
582      dual_inf /= (Number)n_dual;
583      DBG_ASSERT(n_pri>0 || primal_inf==0.);
584      if (n_pri>0) {
585        primal_inf /= (Number)n_pri;
586      }
587      DBG_ASSERT(n_comp>0 || complty==0.);
588      if (n_comp>0) {
589        complty /= (Number)n_comp;
590      }
591      break;
592      case QualityFunctionMuOracle::NM_NORM_2_SQUARED:
593      dual_inf =
594        IpCq().curr_dual_infeasibility(NORM_2);
595      dual_inf *= dual_inf;
596      primal_inf =
597        IpCq().curr_primal_infeasibility(NORM_2);
598      primal_inf *= primal_inf;
599      complty =
600        IpCq().curr_complementarity(0., NORM_2);
601      complty *= complty;
602      dual_inf /= (Number)n_dual;
603      DBG_ASSERT(n_pri>0 || primal_inf==0.);
604      if (n_pri>0) {
605        primal_inf /= (Number)n_pri;
606      }
607      DBG_ASSERT(n_comp>0 || complty==0.);
608      if (n_comp>0) {
609        complty /= (Number)n_comp;
610      }
611      break;
612      case QualityFunctionMuOracle::NM_NORM_MAX:
613      dual_inf =
614        IpCq().curr_dual_infeasibility(NORM_MAX);
615      primal_inf =
616        IpCq().curr_primal_infeasibility(NORM_MAX);
617      complty =
618        IpCq().curr_complementarity(0., NORM_MAX);
619      break;
620      case QualityFunctionMuOracle::NM_NORM_2:
621      dual_inf =
622        IpCq().curr_dual_infeasibility(NORM_2);
623      primal_inf =
624        IpCq().curr_primal_infeasibility(NORM_2);
625      complty =
626        IpCq().curr_complementarity(0., NORM_2);
627      dual_inf /= sqrt((Number)n_dual);
628      DBG_ASSERT(n_pri>0 || primal_inf==0.);
629      if (n_pri>0) {
630        primal_inf /= sqrt((Number)n_pri);
631      }
632      DBG_ASSERT(n_comp>0 || complty==0.);
633      if (n_comp>0) {
634        complty /= sqrt((Number)n_comp);
635      }
636      break;
637    }
638
639    Number centrality = 0.;
640    if (adaptive_mu_kkt_centrality_!=0) {
641      Number xi = IpCq().curr_centrality_measure();
642      switch (adaptive_mu_kkt_centrality_) {
643        case 1:
644        centrality = -complty*log(xi);
645        break;
646        case 2:
647        centrality = complty/xi;
648        case 3:
649        centrality = complty/pow(xi,3);
650        break;
651        default:
652        DBG_ASSERT("Unknown value for adaptive_mu_kkt_centrality_");
653      }
654    }
655
656    Number balancing_term=0.;
657    switch (adaptive_mu_kkt_balancing_term_) {
658      case 0:
659      //Nothing
660      break;
661      case 1:
662      balancing_term = pow(Max(0., Max(dual_inf,primal_inf)-complty),3);
663      break;
664      default:
665      DBG_ASSERT("Unknown value for adaptive_mu_kkt_balancing_term");
666    }
667
668    DBG_ASSERT(centrality>=0.);
669    DBG_ASSERT(balancing_term>=0);
670    Number kkt_error = primal_inf + dual_inf + complty +
671                       centrality + balancing_term;
672
673    Jnlst().Printf(J_MOREDETAILED, J_BARRIER_UPDATE,
674                   "KKT error in barrier update check:\n"
675                   "  primal infeasibility: %15.6e\n"
676                   "    dual infeasibility: %15.6e\n"
677                   "       complementarity: %15.6e\n"
678                   "            centrality: %15.6e\n"
679                   "             kkt error: %15.6e\n",
680                   primal_inf, dual_inf, complty, centrality, kkt_error);
681
682    return kkt_error;
683  }
684
685  Number
686  AdaptiveMuUpdate::lower_mu_safeguard()
687  {
688    DBG_START_METH("AdaptiveMuUpdate::lower_mu_safeguard",
689                   dbg_verbosity);
690    if (adaptive_mu_safeguard_factor_ == 0.)
691      return 0.;
692
693    Number dual_inf =
694      IpCq().curr_dual_infeasibility(NORM_1);
695    Number primal_inf =
696      IpCq().curr_primal_infeasibility(NORM_1);
697    Index n_dual = IpData().curr()->x()->Dim() + IpData().curr()->s()->Dim();
698    dual_inf /= (Number)n_dual;
699    Index n_pri = IpData().curr()->y_c()->Dim() + IpData().curr()->y_d()->Dim();
700    DBG_ASSERT(n_pri>0 || primal_inf==0.);
701    if (n_pri>0) {
702      primal_inf /= (Number)n_pri;
703    }
704
705    if (init_dual_inf_ < 0.) {
706      init_dual_inf_ = Max(1., dual_inf);
707    }
708    if (init_primal_inf_ < 0.) {
709      init_primal_inf_ = Max(1., primal_inf);
710    }
711
712    Number lower_mu_safeguard =
713      Max(adaptive_mu_safeguard_factor_ * (dual_inf/init_dual_inf_),
714          adaptive_mu_safeguard_factor_ * (primal_inf/init_primal_inf_));
715    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_));
716
717    if (adaptive_mu_globalization_==KKT_ERROR) {
718      lower_mu_safeguard = Min(lower_mu_safeguard, min_ref_val());
719    }
720
721    return lower_mu_safeguard;
722  }
723
724} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.