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