source: branches/dev/Algorithm/IpQualityFunctionMuOracle.cpp @ 527

Last change on this file since 527 was 527, checked in by andreasw, 15 years ago
  • corrected CUTEr interface
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 46.9 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: IpQualityFunctionMuOracle.cpp 527 2005-09-28 03:33:46Z andreasw $
6//
7// Authors:  Carl Laird, Andreas Waechter            IBM    2004-11-12
8
9#include "IpQualityFunctionMuOracle.hpp"
10
11#ifdef HAVE_CMATH
12# include <cmath>
13#else
14# ifdef HAVE_MATH_H
15#  include <math.h>
16# else
17#  error "don't have header file for math"
18# endif
19#endif
20
21namespace Ipopt
22{
23#ifdef IP_DEBUG
24  static const Index dbg_verbosity = 0;
25#endif
26
27  QualityFunctionMuOracle::QualityFunctionMuOracle(const SmartPtr<PDSystemSolver>& pd_solver)
28      :
29      MuOracle(),
30      pd_solver_(pd_solver),
31
32      tmp_step_x_L_(NULL),
33      tmp_step_x_U_(NULL),
34      tmp_step_s_L_(NULL),
35      tmp_step_s_U_(NULL),
36      tmp_step_z_L_(NULL),
37      tmp_step_z_U_(NULL),
38      tmp_step_v_L_(NULL),
39      tmp_step_v_U_(NULL),
40
41      tmp_slack_x_L_(NULL),
42      tmp_slack_x_U_(NULL),
43      tmp_slack_s_L_(NULL),
44      tmp_slack_s_U_(NULL),
45      tmp_z_L_(NULL),
46      tmp_z_U_(NULL),
47      tmp_v_L_(NULL),
48      tmp_v_U_(NULL),
49
50      count_qf_evals_(0)
51  {
52    DBG_ASSERT(IsValid(pd_solver_));
53  }
54
55  QualityFunctionMuOracle::~QualityFunctionMuOracle()
56  {}
57
58  void QualityFunctionMuOracle::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
59  {
60    roptions->AddLowerBoundedNumberOption(
61      "sigma_max",
62      "Maximum value of the centering parameter.",
63      0.0, true, 1e2);
64    roptions->AddStringOption4(
65      "quality_function_norm_type",
66      "Norm used for components of the quality function.",
67      "2-norm-squared",
68      "1-norm", "use the 1-norm (abs sum)",
69      "2-norm-squared", "use the 2-norm squared (sum of squares)",
70      "max-norm", "use the infinity norm (max)",
71      "2-norm", "use 2-norm");
72    roptions->AddStringOption4(
73      "quality_function_centrality",
74      "The penalty term for centrality that is included in quality function.",
75      "none",
76      "none", "no penalty term is added",
77      "log", "complementarity * the log of the centrality measure",
78      "reciprocal", "complementarity * the reciprocal of the centrality measure",
79      "cubed-reciprocal", "complementarity * the reciprocal of the centrality measure cubed",
80      "This determines whether a term is added to the quality function to "
81      "penalize deviation from centrality with respect to complementarity.  The "
82      "complementarity measure here is the xi in the Loqo update rule.");
83    roptions->AddStringOption2(
84      "quality_function_balancing_term",
85      "The balancing term included in the quality function for centrality.",
86      "none",
87      "none", "no balancing term is added",
88      "cubic", "Max(0,Max(dual_inf,primal_inf)-compl)^3",
89      "This determines whether a term is added to the quality function that "
90      "penalizes situations where the complementarity is much smaller "
91      "than dual and primal infeasibilities.");
92    roptions->AddLowerBoundedIntegerOption(
93      "quality_function_max_section_steps",
94      "Maximum number of search steps during direct search procedure "
95      "determining the optimal centering parameter.",
96      0, 8,
97      "The golden section search is performed for the quality function based "
98      "mu oractle.");
99    roptions->AddBoundedNumberOption(
100      "quality_function_section_sigma_tol",
101      "Tolerance for the section search procedure determining "
102      "the optimal centering parameter (in sigma space).",
103      0.0, false, 1.0, true,
104      1e-2,
105      "The golden section search is performed for the quality function based "
106      "mu oractle.");
107    roptions->AddBoundedNumberOption(
108      "quality_function_section_qf_tol",
109      "Tolerance for the golden section search procedure determining "
110      "the optimal centering parameter (in the function value space).",
111      0.0, false, 1.0, true,
112      0e-2,
113      "The section search is performed for the quality function based mu "
114      "oractle.");
115  }
116
117
118  bool QualityFunctionMuOracle::InitializeImpl(const OptionsList& options,
119      const std::string& prefix)
120  {
121    Index enum_int;
122
123    options.GetNumericValue("sigma_max", sigma_max_, prefix);
124
125    options.GetEnumValue("quality_function_norm_type", enum_int, prefix);
126    quality_function_norm_ = NormEnum(enum_int);
127    options.GetEnumValue("quality_function_centrality", enum_int, prefix);
128    quality_function_centrality_ = CentralityEnum(enum_int);
129    options.GetEnumValue("quality_function_balancing_term", enum_int, prefix);
130    quality_function_balancing_term_ = BalancingTermEnum(enum_int);
131    options.GetIntegerValue("quality_function_max_section_steps",
132                            quality_function_max_section_steps_, prefix);
133    options.GetNumericValue("quality_function_section_sigma_tol",
134                            quality_function_section_sigma_tol_, prefix);
135    options.GetNumericValue("quality_function_section_qf_tol",
136                            quality_function_section_qf_tol_, prefix);
137
138    return true;
139  }
140
141  Number QualityFunctionMuOracle::CalculateMu(Number mu_min, Number mu_max)
142  {
143    DBG_START_METH("QualityFunctionMuOracle::CalculateMu",
144                   dbg_verbosity);
145
146    count_qf_evals_ = 0;
147
148    ///////////////////////////////////////////////////////////////////////////
149    // Reserve memory for temporary vectors used in CalculateQualityFunction //
150    ///////////////////////////////////////////////////////////////////////////
151
152    tmp_step_x_L_ = IpNLP().x_L()->MakeNew();
153    tmp_step_x_U_ = IpNLP().x_U()->MakeNew();
154    tmp_step_s_L_ = IpNLP().d_L()->MakeNew();
155    tmp_step_s_U_ = IpNLP().d_U()->MakeNew();
156    tmp_step_z_L_ = IpNLP().x_L()->MakeNew();
157    tmp_step_z_U_ = IpNLP().x_U()->MakeNew();
158    tmp_step_v_L_ = IpNLP().d_L()->MakeNew();
159    tmp_step_v_U_ = IpNLP().d_U()->MakeNew();
160
161    tmp_slack_x_L_ = IpNLP().x_L()->MakeNew();
162    tmp_slack_x_U_ = IpNLP().x_U()->MakeNew();
163    tmp_slack_s_L_ = IpNLP().d_L()->MakeNew();
164    tmp_slack_s_U_ = IpNLP().d_U()->MakeNew();
165    tmp_z_L_ = IpNLP().x_L()->MakeNew();
166    tmp_z_U_ = IpNLP().x_U()->MakeNew();
167    tmp_v_L_ = IpNLP().d_L()->MakeNew();
168    tmp_v_U_ = IpNLP().d_U()->MakeNew();
169
170    /////////////////////////////////////
171    // Compute the affine scaling step //
172    /////////////////////////////////////
173
174    Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
175                   "Solving the Primal Dual System for the affine step\n");
176    // First get the right hand side
177    SmartPtr<IteratesVector> rhs_aff = IpData().curr()->MakeNewIteratesVector(false);
178    rhs_aff->Set_x(*IpCq().curr_grad_lag_x());
179    rhs_aff->Set_s(*IpCq().curr_grad_lag_s());
180    rhs_aff->Set_y_c(*IpCq().curr_c());
181    rhs_aff->Set_y_d(*IpCq().curr_d_minus_s());
182    rhs_aff->Set_z_L(*IpCq().curr_compl_x_L());
183    rhs_aff->Set_z_U(*IpCq().curr_compl_x_U());
184    rhs_aff->Set_v_L(*IpCq().curr_compl_s_L());
185    rhs_aff->Set_v_U(*IpCq().curr_compl_s_U());
186
187    // Get space for the affine scaling step
188    SmartPtr<IteratesVector> step_aff = IpData().curr()->MakeNewIteratesVector(true);
189
190    // Now solve the primal-dual system to get the step
191    pd_solver_->Solve(-1.0, 0.0,
192                      *rhs_aff,
193                      *step_aff,
194                      false           // want accurate solution here
195                      // because we can use it to
196                      // compute the overall search
197                      // direction
198                     );
199
200    DBG_PRINT_VECTOR(2, "step_aff", *step_aff);
201
202    /////////////////////////////////////
203    // Compute the pure centering step //
204    /////////////////////////////////////
205
206    Number avrg_compl = IpCq().curr_avrg_compl();
207
208    Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
209                   "Solving the Primal Dual System for the centering step\n");
210    // First get the right hand side
211    SmartPtr<IteratesVector> rhs_cen = IpData().curr()->MakeNewIteratesVector(true);
212    rhs_cen->x_NonConst()->AddOneVector(-avrg_compl,
213                                        *IpCq().grad_kappa_times_damping_x(),
214                                        0.);
215    rhs_cen->s_NonConst()->AddOneVector(-avrg_compl,
216                                        *IpCq().grad_kappa_times_damping_s(),
217                                        0.);
218
219    rhs_cen->y_c_NonConst()->Set(0.);
220    rhs_cen->y_d_NonConst()->Set(0.);
221    rhs_cen->z_L_NonConst()->Set(avrg_compl);
222    rhs_cen->z_U_NonConst()->Set(avrg_compl);
223    rhs_cen->v_L_NonConst()->Set(avrg_compl);
224    rhs_cen->v_U_NonConst()->Set(avrg_compl);
225
226    // Get space for the centering step
227    SmartPtr<IteratesVector> step_cen = IpData().curr()->MakeNewIteratesVector(true);
228
229    // Now solve the primal-dual system to get the step
230    pd_solver_->Solve(1.0, 0.0,
231                      *rhs_cen,
232                      *step_cen,
233                      false           // want accurate solution here
234                      // because we can use it to
235                      // compute the overall search
236                      // direction
237                     );
238
239    DBG_PRINT_VECTOR(2, "step_cen", *step_cen);
240
241    // Start the timing for the quality function search here
242    IpData().TimingStats().QualityFunctionSearch.Start();
243
244    // We now compute the step for the slack variables.  This safes
245    // time, because we then don't have to do this any more for each
246    // evaluation of the quality function
247    SmartPtr<Vector> step_aff_x_L = step_aff->z_L()->MakeNew();
248    SmartPtr<Vector> step_aff_x_U = step_aff->z_U()->MakeNew();
249    SmartPtr<Vector> step_aff_s_L = step_aff->v_L()->MakeNew();
250    SmartPtr<Vector> step_aff_s_U = step_aff->v_U()->MakeNew();
251    IpNLP().Px_L()->TransMultVector(1., *step_aff->x(), 0., *step_aff_x_L);
252    IpNLP().Px_U()->TransMultVector(-1., *step_aff->x(), 0., *step_aff_x_U);
253    IpNLP().Pd_L()->TransMultVector(1., *step_aff->s(), 0., *step_aff_s_L);
254    IpNLP().Pd_U()->TransMultVector(-1., *step_aff->s(), 0., *step_aff_s_U);
255    SmartPtr<Vector> step_cen_x_L = step_cen->z_L()->MakeNew();
256    SmartPtr<Vector> step_cen_x_U = step_cen->z_U()->MakeNew();
257    SmartPtr<Vector> step_cen_s_L = step_cen->v_L()->MakeNew();
258    SmartPtr<Vector> step_cen_s_U = step_cen->v_U()->MakeNew();
259    IpNLP().Px_L()->TransMultVector(1., *step_cen->x(), 0., *step_cen_x_L);
260    IpNLP().Px_U()->TransMultVector(-1., *step_cen->x(), 0., *step_cen_x_U);
261    IpNLP().Pd_L()->TransMultVector(1., *step_cen->s(), 0., *step_cen_s_L);
262    IpNLP().Pd_U()->TransMultVector(-1., *step_cen->s(), 0., *step_cen_s_U);
263
264    Number sigma;
265
266    // First we determine whether we want to search for a value of
267    // sigma larger or smaller than 1.  For this, we estimate the
268    // slope of the quality function at sigma=1.
269    Number qf_1 = CalculateQualityFunction(1.,
270                                           *step_aff_x_L,
271                                           *step_aff_x_U,
272                                           *step_aff_s_L,
273                                           *step_aff_s_U,
274                                           *step_aff->y_c(),
275                                           *step_aff->y_d(),
276                                           *step_aff->z_L(),
277                                           *step_aff->z_U(),
278                                           *step_aff->v_L(),
279                                           *step_aff->v_U(),
280                                           *step_cen_x_L,
281                                           *step_cen_x_U,
282                                           *step_cen_s_L,
283                                           *step_cen_s_U,
284                                           *step_cen->y_c(),
285                                           *step_cen->y_d(),
286                                           *step_cen->z_L(),
287                                           *step_cen->z_U(),
288                                           *step_cen->v_L(),
289                                           *step_cen->v_U());
290
291    Number sigma_1minus = 1.-quality_function_section_sigma_tol_;
292    Number qf_1minus = CalculateQualityFunction(sigma_1minus,
293                       *step_aff_x_L,
294                       *step_aff_x_U,
295                       *step_aff_s_L,
296                       *step_aff_s_U,
297                       *step_aff->y_c(),
298                       *step_aff->y_d(),
299                       *step_aff->z_L(),
300                       *step_aff->z_U(),
301                       *step_aff->v_L(),
302                       *step_aff->v_U(),
303                       *step_cen_x_L,
304                       *step_cen_x_U,
305                       *step_cen_s_L,
306                       *step_cen_s_U,
307                       *step_cen->y_c(),
308                       *step_cen->y_d(),
309                       *step_cen->z_L(),
310                       *step_cen->z_U(),
311                       *step_cen->v_L(),
312                       *step_cen->v_U());
313
314    if (qf_1minus > qf_1) {
315      // It seems that the quality function decreases for values
316      // larger than sigma, so perform golden section search for sigma
317      // > 1.
318      Number sigma_up = Min(sigma_max_, mu_max/avrg_compl);
319      Number sigma_lo = 1.;
320      // ToDo maybe we should use different tolerances for sigma>1
321      sigma = PerformGoldenSection(sigma_up, -100., sigma_lo, qf_1,
322                                   quality_function_section_sigma_tol_,
323                                   quality_function_section_qf_tol_,
324                                   *step_aff_x_L,
325                                   *step_aff_x_U,
326                                   *step_aff_s_L,
327                                   *step_aff_s_U,
328                                   *step_aff->y_c(),
329                                   *step_aff->y_d(),
330                                   *step_aff->z_L(),
331                                   *step_aff->z_U(),
332                                   *step_aff->v_L(),
333                                   *step_aff->v_U(),
334                                   *step_cen_x_L,
335                                   *step_cen_x_U,
336                                   *step_cen_s_L,
337                                   *step_cen_s_U,
338                                   *step_cen->y_c(),
339                                   *step_cen->y_d(),
340                                   *step_cen->z_L(),
341                                   *step_cen->z_U(),
342                                   *step_cen->v_L(),
343                                   *step_cen->v_U());
344    }
345    else {
346      // Search for sigma less than 1
347
348      Number sigma_lo = mu_min/avrg_compl;
349      Number sigma_up = Max(sigma_lo, sigma_1minus);
350      // Todo: What if lower bound is less than sigma_lo? Skip search?
351      sigma = PerformGoldenSection(sigma_up, qf_1minus, sigma_lo, -100.,
352                                   quality_function_section_sigma_tol_,
353                                   quality_function_section_qf_tol_,
354                                   *step_aff_x_L,
355                                   *step_aff_x_U,
356                                   *step_aff_s_L,
357                                   *step_aff_s_U,
358                                   *step_aff->y_c(),
359                                   *step_aff->y_d(),
360                                   *step_aff->z_L(),
361                                   *step_aff->z_U(),
362                                   *step_aff->v_L(),
363                                   *step_aff->v_U(),
364                                   *step_cen_x_L,
365                                   *step_cen_x_U,
366                                   *step_cen_s_L,
367                                   *step_cen_s_U,
368                                   *step_cen->y_c(),
369                                   *step_cen->y_d(),
370                                   *step_cen->z_L(),
371                                   *step_cen->z_U(),
372                                   *step_cen->v_L(),
373                                   *step_cen->v_U());
374    }
375
376    //#define tracequalityfunction
377#ifdef tracequalityfunction
378    char fname[100];
379    sprintf(fname, "qf_values_%d.dat", IpData().iter_count());
380    FILE* fid = fopen(fname, "w");
381
382    Number sigma_1 = sigma_max_;
383    Number sigma_2 = 1e-9/avrg_compl;
384    Number sigma_trace = sigma_1;
385    while(sigma_trace > sigma_2) {
386      Number qf = CalculateQualityFunction(sigma_trace,
387                                           *step_aff_x_L,
388                                           *step_aff_x_U,
389                                           *step_aff_s_L,
390                                           *step_aff_s_U,
391                                           *step_aff->y_c(),
392                                           *step_aff->y_d(),
393                                           *step_aff->z_L(),
394                                           *step_aff->z_U(),
395                                           *step_aff->v_L(),
396                                           *step_aff->v_U(),
397                                           *step_cen_x_L,
398                                           *step_cen_x_U,
399                                           *step_cen_s_L,
400                                           *step_cen_s_U,
401                                           *step_cen->y_c(),
402                                           *step_cen->y_d(),
403                                           *step_cen->z_L(),
404                                           *step_cen->z_U(),
405                                           *step_cen->v_L(),
406                                           *step_cen->v_U());
407      fprintf(fid, "%9.2e %25.16e\n", sigma_trace, qf);
408      sigma_trace /= 1.1;
409    }
410    fclose(fid);
411#endif
412
413    // End timing of quality function search
414    IpData().TimingStats().QualityFunctionSearch.End();
415
416    Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
417                   "Sigma = %e\n", sigma);
418    Number mu = sigma*avrg_compl;
419
420    // Store the affine search direction (in case it is needed in the
421    // line search for a corrector step)
422    IpData().set_delta_aff(step_aff);
423    IpData().SetHaveAffineDeltas(true);
424
425    // Now construct the overall search direction here
426    SmartPtr<IteratesVector> step = IpData().curr()->MakeNewIteratesVector(true);
427    step->AddTwoVectors(sigma, *step_cen, 1.0, *IpData().delta_aff(), 0.0);
428
429    IpData().set_delta(step);
430    IpData().SetHaveDeltas(true);
431
432    ///////////////////////////////////////////////////////////////////////////
433    // Release memory for temporary vectors used in CalculateQualityFunction //
434    ///////////////////////////////////////////////////////////////////////////
435
436    tmp_step_x_L_ = NULL;
437    tmp_step_x_U_ = NULL;
438    tmp_step_s_L_ = NULL;
439    tmp_step_s_U_ = NULL;
440    tmp_step_z_L_ = NULL;
441    tmp_step_z_U_ = NULL;
442    tmp_step_v_L_ = NULL;
443    tmp_step_v_U_ = NULL;
444
445    tmp_slack_x_L_ = NULL;
446    tmp_slack_x_U_ = NULL;
447    tmp_slack_s_L_ = NULL;
448    tmp_slack_s_U_ = NULL;
449    tmp_z_L_ = NULL;
450    tmp_z_U_ = NULL;
451    tmp_v_L_ = NULL;
452    tmp_v_U_ = NULL;
453
454    // DELETEME
455    char ssigma[40];
456    sprintf(ssigma, " sigma=%8.2e", sigma);
457    IpData().Append_info_string(ssigma);
458    sprintf(ssigma, " qf=%d", count_qf_evals_);
459    IpData().Append_info_string(ssigma);
460    /*
461    sprintf(ssigma, " xi=%8.2e ", IpCq().curr_centrality_measure());
462    IpData().Append_info_string(ssigma);
463    if (sigma>1.) {
464      IpData().Append_info_string("LARGESIGMA");
465    }
466    */
467
468    return mu;
469  }
470
471  Number QualityFunctionMuOracle::CalculateQualityFunction
472  (Number sigma,
473   const Vector& step_aff_x_L,
474   const Vector& step_aff_x_U,
475   const Vector& step_aff_s_L,
476   const Vector& step_aff_s_U,
477   const Vector& step_aff_y_c,
478   const Vector& step_aff_y_d,
479   const Vector& step_aff_z_L,
480   const Vector& step_aff_z_U,
481   const Vector& step_aff_v_L,
482   const Vector& step_aff_v_U,
483   const Vector& step_cen_x_L,
484   const Vector& step_cen_x_U,
485   const Vector& step_cen_s_L,
486   const Vector& step_cen_s_U,
487   const Vector& step_cen_y_c,
488   const Vector& step_cen_y_d,
489   const Vector& step_cen_z_L,
490   const Vector& step_cen_z_U,
491   const Vector& step_cen_v_L,
492   const Vector& step_cen_v_U
493  )
494  {
495    DBG_START_METH("QualityFunctionMuOracle::CalculateQualityFunction",
496                   dbg_verbosity);
497    count_qf_evals_++;
498
499    Index n_dual = IpData().curr()->x()->Dim() + IpData().curr()->s()->Dim();
500    Index n_pri = IpData().curr()->y_c()->Dim() + IpData().curr()->y_d()->Dim();
501    Index n_comp = IpData().curr()->z_L()->Dim() + IpData().curr()->z_U()->Dim() +
502                   IpData().curr()->v_L()->Dim() + IpData().curr()->v_U()->Dim();
503
504    IpData().TimingStats().Task1.Start();
505    tmp_step_x_L_->AddTwoVectors(1., step_aff_x_L, sigma, step_cen_x_L, 0.);
506    tmp_step_x_U_->AddTwoVectors(1., step_aff_x_U, sigma, step_cen_x_U, 0.);
507    tmp_step_s_L_->AddTwoVectors(1., step_aff_s_L, sigma, step_cen_s_L, 0.);
508    tmp_step_s_U_->AddTwoVectors(1., step_aff_s_U, sigma, step_cen_s_U, 0.);
509    tmp_step_z_L_->AddTwoVectors(1., step_aff_z_L, sigma, step_cen_z_L, 0.);
510    tmp_step_z_U_->AddTwoVectors(1., step_aff_z_U, sigma, step_cen_z_U, 0.);
511    tmp_step_v_L_->AddTwoVectors(1., step_aff_v_L, sigma, step_cen_v_L, 0.);
512    tmp_step_v_U_->AddTwoVectors(1., step_aff_v_U, sigma, step_cen_v_U, 0.);
513    IpData().TimingStats().Task1.End();
514
515    // Compute the fraction-to-the-boundary step sizes
516    IpData().TimingStats().Task2.Start();
517    Number tau = IpData().curr_tau();
518    Number alpha_primal = IpCq().uncached_slack_frac_to_the_bound(tau,
519                          *tmp_step_x_L_,
520                          *tmp_step_x_U_,
521                          *tmp_step_s_L_,
522                          *tmp_step_s_U_);
523
524    Number alpha_dual = IpCq().uncached_dual_frac_to_the_bound(tau,
525                        *tmp_step_z_L_,
526                        *tmp_step_z_U_,
527                        *tmp_step_v_L_,
528                        *tmp_step_v_U_);
529    IpData().TimingStats().Task2.End();
530
531    Number xi = 0.; // centrality measure
532
533    IpData().TimingStats().Task1.Start();
534    tmp_slack_x_L_->AddTwoVectors(1., *IpCq().curr_slack_x_L(),
535                                  alpha_primal, *tmp_step_x_L_, 0.);
536    tmp_slack_x_U_->AddTwoVectors(1., *IpCq().curr_slack_x_U(),
537                                  alpha_primal, *tmp_step_x_U_, 0.);
538    tmp_slack_s_L_->AddTwoVectors(1., *IpCq().curr_slack_s_L(),
539                                  alpha_primal, *tmp_step_s_L_, 0.);
540    tmp_slack_s_U_->AddTwoVectors(1., *IpCq().curr_slack_s_U(),
541                                  alpha_primal, *tmp_step_s_U_, 0.);
542
543    tmp_z_L_->AddTwoVectors(1., *IpData().curr()->z_L(),
544                            alpha_dual, *tmp_step_z_L_, 0.);
545    tmp_z_U_->AddTwoVectors(1., *IpData().curr()->z_U(),
546                            alpha_dual, *tmp_step_z_U_, 0.);
547    tmp_v_L_->AddTwoVectors(1., *IpData().curr()->v_L(),
548                            alpha_dual, *tmp_step_v_L_, 0.);
549    tmp_v_U_->AddTwoVectors(1., *IpData().curr()->v_U(),
550                            alpha_dual, *tmp_step_v_U_, 0.);
551    IpData().TimingStats().Task1.End();
552
553    IpData().TimingStats().Task3.Start();
554    tmp_slack_x_L_->ElementWiseMultiply(*tmp_z_L_);
555    tmp_slack_x_U_->ElementWiseMultiply(*tmp_z_U_);
556    tmp_slack_s_L_->ElementWiseMultiply(*tmp_v_L_);
557    tmp_slack_s_U_->ElementWiseMultiply(*tmp_v_U_);
558    IpData().TimingStats().Task3.End();
559
560    DBG_PRINT_VECTOR(2, "compl_x_L", *tmp_slack_x_L_);
561    DBG_PRINT_VECTOR(2, "compl_x_U", *tmp_slack_x_U_);
562    DBG_PRINT_VECTOR(2, "compl_s_L", *tmp_slack_s_L_);
563    DBG_PRINT_VECTOR(2, "compl_s_U", *tmp_slack_s_U_);
564
565    Number dual_inf=-1.;
566    Number primal_inf=-1.;
567    Number compl_inf=-1.;
568
569    IpData().TimingStats().Task5.Start();
570    switch (quality_function_norm_) {
571      case NM_NORM_1:
572      dual_inf = (1.-alpha_dual)*(IpCq().curr_grad_lag_x()->Asum() +
573                                  IpCq().curr_grad_lag_s()->Asum());
574
575      primal_inf = (1.-alpha_primal)*(IpCq().curr_c()->Asum() +
576                                      IpCq().curr_d_minus_s()->Asum());
577
578      compl_inf = tmp_slack_x_L_->Asum() + tmp_slack_x_U_->Asum() +
579                  tmp_slack_s_L_->Asum() + tmp_slack_s_U_->Asum();
580
581      dual_inf /= n_dual;
582      if (n_pri>0) {
583        primal_inf /= n_pri;
584      }
585      DBG_ASSERT(n_comp>0);
586      compl_inf /= n_comp;
587      break;
588      case NM_NORM_2_SQUARED:
589      dual_inf =
590        pow(1.-alpha_dual, 2)*(pow(IpCq().curr_grad_lag_x()->Nrm2(), 2) +
591                               pow(IpCq().curr_grad_lag_s()->Nrm2(), 2));
592      primal_inf =
593        pow(1.-alpha_primal, 2)*(pow(IpCq().curr_c()->Nrm2(), 2) +
594                                 pow(IpCq().curr_d_minus_s()->Nrm2(), 2));
595      compl_inf =
596        pow(tmp_slack_x_L_->Nrm2(), 2) + pow(tmp_slack_x_U_->Nrm2(), 2) +
597        pow(tmp_slack_s_L_->Nrm2(), 2) + pow(tmp_slack_s_U_->Nrm2(), 2);
598
599      dual_inf /= n_dual;
600      if (n_pri>0) {
601        primal_inf /= n_pri;
602      }
603      DBG_ASSERT(n_comp>0);
604      compl_inf /= n_comp;
605      break;
606      case NM_NORM_MAX:
607      dual_inf =
608        (1.-alpha_dual)*Max(IpCq().curr_grad_lag_x()->Amax(),
609                            IpCq().curr_grad_lag_s()->Amax());
610      primal_inf =
611        (1.-alpha_primal)*Max(IpCq().curr_c()->Amax(),
612                              IpCq().curr_d_minus_s()->Amax());
613      compl_inf =
614        Max(tmp_slack_x_L_->Amax(), tmp_slack_x_U_->Amax(),
615            tmp_slack_s_L_->Amax(), tmp_slack_s_U_->Amax());
616      break;
617      case NM_NORM_2:
618      dual_inf =
619        (1.-alpha_dual)*sqrt(pow(IpCq().curr_grad_lag_x()->Nrm2(), 2) +
620                             pow(IpCq().curr_grad_lag_s()->Nrm2(), 2));
621      primal_inf =
622        (1.-alpha_primal)*sqrt(pow(IpCq().curr_c()->Nrm2(), 2) +
623                               pow(IpCq().curr_d_minus_s()->Nrm2(), 2));
624      compl_inf =
625        sqrt(pow(tmp_slack_x_L_->Nrm2(), 2) + pow(tmp_slack_x_U_->Nrm2(), 2) +
626             pow(tmp_slack_s_L_->Nrm2(), 2) + pow(tmp_slack_s_U_->Nrm2(), 2));
627
628      dual_inf /= sqrt((Number)n_dual);
629      if (n_pri>0) {
630        primal_inf /= sqrt((Number)n_pri);
631      }
632      DBG_ASSERT(n_comp>0);
633      compl_inf /= sqrt((Number)n_comp);
634      break;
635      default:
636      DBG_ASSERT(false && "Unknown value for quality_function_norm_");
637    }
638    IpData().TimingStats().Task5.End();
639
640    Number quality_function = dual_inf + primal_inf + compl_inf;
641
642    if (quality_function_centrality_!=CEN_NONE) {
643      IpData().TimingStats().Task4.Start();
644      xi = IpCq().CalcCentralityMeasure(*tmp_slack_x_L_, *tmp_slack_x_U_,
645                                        *tmp_slack_s_L_, *tmp_slack_s_U_);
646      IpData().TimingStats().Task4.End();
647    }
648    switch (quality_function_centrality_) {
649      case CEN_NONE:
650      //Nothing
651      break;
652      case CEN_LOG:
653      quality_function -= compl_inf*log(xi);
654      break;
655      case CEN_RECIPROCAL:
656      quality_function += compl_inf/xi;
657      case CEN_CUBED_RECIPROCAL:
658      quality_function += compl_inf/pow(xi,3);
659      break;
660      default:
661      DBG_ASSERT("Unknown value for quality_function_centrality_");
662    }
663
664    switch (quality_function_balancing_term_) {
665      case BT_NONE:
666      //Nothing
667      break;
668      case BT_CUBIC:
669      quality_function += pow(Max(0., Max(dual_inf,primal_inf)-compl_inf),3);
670      break;
671      default:
672      DBG_ASSERT("Unknown value for quality_function_balancing term_");
673    }
674
675    Jnlst().Printf(J_MOREDETAILED, J_BARRIER_UPDATE,
676                   "sigma = %8.2e d_inf = %18.12e p_inf = %18.12e cmpl = %18.12e q = %18.12e a_pri = %8.2e a_dual = %8.2e xi = %8.2e\n", sigma, dual_inf, primal_inf, compl_inf, quality_function, alpha_primal, alpha_dual, xi);
677
678
679    return quality_function;
680    //return compl_inf;
681  }
682
683  Number
684  QualityFunctionMuOracle::PerformGoldenSection
685  (Number sigma_up_in,
686   Number q_up,
687   Number sigma_lo_in,
688   Number q_lo,
689   Number sigma_tol,
690   Number qf_tol,
691   const Vector& step_aff_x_L,
692   const Vector& step_aff_x_U,
693   const Vector& step_aff_s_L,
694   const Vector& step_aff_s_U,
695   const Vector& step_aff_y_c,
696   const Vector& step_aff_y_d,
697   const Vector& step_aff_z_L,
698   const Vector& step_aff_z_U,
699   const Vector& step_aff_v_L,
700   const Vector& step_aff_v_U,
701   const Vector& step_cen_x_L,
702   const Vector& step_cen_x_U,
703   const Vector& step_cen_s_L,
704   const Vector& step_cen_s_U,
705   const Vector& step_cen_y_c,
706   const Vector& step_cen_y_d,
707   const Vector& step_cen_z_L,
708   const Vector& step_cen_z_U,
709   const Vector& step_cen_v_L,
710   const Vector& step_cen_v_U
711  )
712  {
713    Number sigma_up = ScaleSigma(sigma_up_in);
714    Number sigma_lo = ScaleSigma(sigma_lo_in);
715
716    Number sigma;
717    Number gfac = (3.-sqrt(5.))/2.;
718    Number sigma_mid1 = sigma_lo + gfac*(sigma_up-sigma_lo);
719    Number sigma_mid2 = sigma_lo + (1.-gfac)*(sigma_up-sigma_lo);
720
721    Number qmid1 = CalculateQualityFunction(UnscaleSigma(sigma_mid1),
722                                            step_aff_x_L,
723                                            step_aff_x_U,
724                                            step_aff_s_L,
725                                            step_aff_s_U,
726                                            step_aff_y_c,
727                                            step_aff_y_d,
728                                            step_aff_z_L,
729                                            step_aff_z_U,
730                                            step_aff_v_L,
731                                            step_aff_v_U,
732                                            step_cen_x_L,
733                                            step_cen_x_U,
734                                            step_cen_s_L,
735                                            step_cen_s_U,
736                                            step_cen_y_c,
737                                            step_cen_y_d,
738                                            step_cen_z_L,
739                                            step_cen_z_U,
740                                            step_cen_v_L,
741                                            step_cen_v_U);
742    Number qmid2 = CalculateQualityFunction(UnscaleSigma(sigma_mid2),
743                                            step_aff_x_L,
744                                            step_aff_x_U,
745                                            step_aff_s_L,
746                                            step_aff_s_U,
747                                            step_aff_y_c,
748                                            step_aff_y_d,
749                                            step_aff_z_L,
750                                            step_aff_z_U,
751                                            step_aff_v_L,
752                                            step_aff_v_U,
753                                            step_cen_x_L,
754                                            step_cen_x_U,
755                                            step_cen_s_L,
756                                            step_cen_s_U,
757                                            step_cen_y_c,
758                                            step_cen_y_d,
759                                            step_cen_z_L,
760                                            step_cen_z_U,
761                                            step_cen_v_L,
762                                            step_cen_v_U);
763
764    Index nsections = 0;
765    while ((sigma_up-sigma_lo)>=sigma_tol*sigma_up &&
766           //while ((sigma_up-sigma_lo)>=sigma_tol &&  // Note we are using the non-relative criterion here for sigma
767           (1.-Min(q_lo, q_up, qmid1, qmid2)/Max(q_lo, q_up, qmid1, qmid2))>=qf_tol &&
768           nsections<quality_function_max_section_steps_) {
769      nsections++;
770      //printf("sigma_lo=%e sigma_mid1=%e sigma_mid2=%e sigma_up=%e\n",sigma_lo, sigma_mid1, sigma_mid2, sigma_up);
771      if (qmid1 > qmid2) {
772        sigma_lo = sigma_mid1;
773        q_lo = qmid1;
774        sigma_mid1 = sigma_mid2;
775        qmid1 = qmid2;
776        sigma_mid2 = sigma_lo + (1.-gfac)*(sigma_up-sigma_lo);
777        qmid2 = CalculateQualityFunction(UnscaleSigma(sigma_mid2),
778                                         step_aff_x_L,
779                                         step_aff_x_U,
780                                         step_aff_s_L,
781                                         step_aff_s_U,
782                                         step_aff_y_c,
783                                         step_aff_y_d,
784                                         step_aff_z_L,
785                                         step_aff_z_U,
786                                         step_aff_v_L,
787                                         step_aff_v_U,
788                                         step_cen_x_L,
789                                         step_cen_x_U,
790                                         step_cen_s_L,
791                                         step_cen_s_U,
792                                         step_cen_y_c,
793                                         step_cen_y_d,
794                                         step_cen_z_L,
795                                         step_cen_z_U,
796                                         step_cen_v_L,
797                                         step_cen_v_U);
798      }
799      else {
800        sigma_up = sigma_mid2;
801        q_up = qmid2;
802        sigma_mid2 = sigma_mid1;
803        qmid2 = qmid1;
804        sigma_mid1 = sigma_lo + gfac*(sigma_up-sigma_lo);
805        qmid1 = CalculateQualityFunction(UnscaleSigma(sigma_mid1),
806                                         step_aff_x_L,
807                                         step_aff_x_U,
808                                         step_aff_s_L,
809                                         step_aff_s_U,
810                                         step_aff_y_c,
811                                         step_aff_y_d,
812                                         step_aff_z_L,
813                                         step_aff_z_U,
814                                         step_aff_v_L,
815                                         step_aff_v_U,
816                                         step_cen_x_L,
817                                         step_cen_x_U,
818                                         step_cen_s_L,
819                                         step_cen_s_U,
820                                         step_cen_y_c,
821                                         step_cen_y_d,
822                                         step_cen_z_L,
823                                         step_cen_z_U,
824                                         step_cen_v_L,
825                                         step_cen_v_U);
826      }
827    }
828
829    if ((sigma_up-sigma_lo)>=sigma_tol*sigma_up &&
830        (1.-Min(q_lo, q_up, qmid1, qmid2)/Max(q_lo, q_up, qmid1, qmid2))<qf_tol) {
831      // The qf tolerance make it stop
832      IpData().Append_info_string("qf_tol ");
833      Number qf_min = Min(q_lo, q_up, qmid1, qmid2);
834      DBG_ASSERT(qf_min>-100.);
835      if (qf_min == q_lo) {
836        sigma = sigma_lo;
837      }
838      else if (qf_min == qmid1) {
839        sigma = sigma_mid1;
840      }
841      else if (qf_min == qmid2) {
842        sigma = sigma_mid2;
843      }
844      else {
845        sigma = sigma_up;
846      }
847    }
848    else {
849      Number q;
850      if (qmid1 < qmid2) {
851        sigma = sigma_mid1;
852        q = qmid1;
853      }
854      else {
855        sigma = sigma_mid2;
856        q = qmid2;
857      }
858      if (sigma_up == ScaleSigma(sigma_up_in)) {
859        Number qtmp;
860        if (q_up<0.) {
861          qtmp = CalculateQualityFunction(UnscaleSigma(sigma_up),
862                                          step_aff_x_L,
863                                          step_aff_x_U,
864                                          step_aff_s_L,
865                                          step_aff_s_U,
866                                          step_aff_y_c,
867                                          step_aff_y_d,
868                                          step_aff_z_L,
869                                          step_aff_z_U,
870                                          step_aff_v_L,
871                                          step_aff_v_U,
872                                          step_cen_x_L,
873                                          step_cen_x_U,
874                                          step_cen_s_L,
875                                          step_cen_s_U,
876                                          step_cen_y_c,
877                                          step_cen_y_d,
878                                          step_cen_z_L,
879                                          step_cen_z_U,
880                                          step_cen_v_L,
881                                          step_cen_v_U);
882        }
883        else {
884          qtmp = q_up;
885        }
886        if (qtmp < q) {
887          sigma = sigma_up;
888          q = qtmp;
889        }
890      }
891      else if (sigma_lo == ScaleSigma(sigma_lo_in)) {
892        Number qtmp;
893        if (q_lo<0.) {
894          qtmp = CalculateQualityFunction(UnscaleSigma(sigma_lo),
895                                          step_aff_x_L,
896                                          step_aff_x_U,
897                                          step_aff_s_L,
898                                          step_aff_s_U,
899                                          step_aff_y_c,
900                                          step_aff_y_d,
901                                          step_aff_z_L,
902                                          step_aff_z_U,
903                                          step_aff_v_L,
904                                          step_aff_v_U,
905                                          step_cen_x_L,
906                                          step_cen_x_U,
907                                          step_cen_s_L,
908                                          step_cen_s_U,
909                                          step_cen_y_c,
910                                          step_cen_y_d,
911                                          step_cen_z_L,
912                                          step_cen_z_U,
913                                          step_cen_v_L,
914                                          step_cen_v_U);
915        }
916        else {
917          qtmp = q_lo;
918        }
919        if (qtmp < q) {
920          sigma = sigma_lo;
921          q = qtmp;
922        }
923      }
924    }
925
926    return UnscaleSigma(sigma);
927  }
928
929  /*
930  Number QualityFunctionMuOracle::ScaleSigma(Number sigma) {return log(sigma);}
931  Number QualityFunctionMuOracle::UnscaleSigma(Number scaled_sigma) {return exp(scaled_sigma);}
932  */
933  Number QualityFunctionMuOracle::ScaleSigma(Number sigma)
934  {
935    return sigma;
936  }
937  Number QualityFunctionMuOracle::UnscaleSigma(Number scaled_sigma)
938  {
939    return scaled_sigma;
940  }
941
942  /* AW: Tried search in the log space, but that was even worse than
943     search in unscaled space */
944  /*
945  Number
946  QualityFunctionMuOracle::PerformGoldenSectionLog
947  (Number sigma_up,
948   Number sigma_lo,
949   Number tol,
950   const Vector& step_aff_x_L,
951   const Vector& step_aff_x_U,
952   const Vector& step_aff_s_L,
953   const Vector& step_aff_s_U,
954   const Vector& step_aff_y_c,
955   const Vector& step_aff_y_d,
956   const Vector& step_aff_z_L,
957   const Vector& step_aff_z_U,
958   const Vector& step_aff_v_L,
959   const Vector& step_aff_v_U,
960   const Vector& step_cen_x_L,
961   const Vector& step_cen_x_U,
962   const Vector& step_cen_s_L,
963   const Vector& step_cen_s_U,
964   const Vector& step_cen_y_c,
965   const Vector& step_cen_y_d,
966   const Vector& step_cen_z_L,
967   const Vector& step_cen_z_U,
968   const Vector& step_cen_v_L,
969   const Vector& step_cen_v_U
970  )
971  {
972    Number log_sigma;
973    Number log_sigma_up = log(sigma_up);
974    Number log_sigma_lo = log(sigma_lo);
975
976    Number log_sigma_up_in = log_sigma_up;
977    Number log_sigma_lo_in = log_sigma_lo;
978    Number gfac = (3.-sqrt(5.))/2.;
979    Number log_sigma_mid1 = log_sigma_lo + gfac*(log_sigma_up-log_sigma_lo);
980    Number log_sigma_mid2 = log_sigma_lo + (1.-gfac)*(log_sigma_up-log_sigma_lo);
981
982    Number qmid1 = CalculateQualityFunction(exp(log_sigma_mid1),
983                                            step_aff_x_L,
984                                            step_aff_x_U,
985                                            step_aff_s_L,
986                                            step_aff_s_U,
987                                            step_aff_y_c,
988                                            step_aff_y_d,
989                                            step_aff_z_L,
990                                            step_aff_z_U,
991                                            step_aff_v_L,
992                                            step_aff_v_U,
993                                            step_cen_x_L,
994                                            step_cen_x_U,
995                                            step_cen_s_L,
996                                            step_cen_s_U,
997                                            step_cen_y_c,
998                                            step_cen_y_d,
999                                            step_cen_z_L,
1000                                            step_cen_z_U,
1001                                            step_cen_v_L,
1002                                            step_cen_v_U);
1003    Number qmid2 = CalculateQualityFunction(exp(log_sigma_mid2),
1004                                            step_aff_x_L,
1005                                            step_aff_x_U,
1006                                            step_aff_s_L,
1007                                            step_aff_s_U,
1008                                            step_aff_y_c,
1009                                            step_aff_y_d,
1010                                            step_aff_z_L,
1011                                            step_aff_z_U,
1012                                            step_aff_v_L,
1013                                            step_aff_v_U,
1014                                            step_cen_x_L,
1015                                            step_cen_x_U,
1016                                            step_cen_s_L,
1017                                            step_cen_s_U,
1018                                            step_cen_y_c,
1019                                            step_cen_y_d,
1020                                            step_cen_z_L,
1021                                            step_cen_z_U,
1022                                            step_cen_v_L,
1023                                            step_cen_v_U);
1024
1025    Index nsections = 0;
1026    while ((log_sigma_up-log_sigma_lo)>=tol*log_sigma_up && nsections<quality_function_max_section_steps_) {
1027      nsections++;
1028      if (qmid1 > qmid2) {
1029        log_sigma_lo = log_sigma_mid1;
1030        log_sigma_mid1 = log_sigma_mid2;
1031        qmid1 = qmid2;
1032        log_sigma_mid2 = log_sigma_lo + (1.-gfac)*(log_sigma_up-log_sigma_lo);
1033        qmid2 = CalculateQualityFunction(exp(log_sigma_mid2),
1034                                         step_aff_x_L,
1035                                         step_aff_x_U,
1036                                         step_aff_s_L,
1037                                         step_aff_s_U,
1038                                         step_aff_y_c,
1039                                         step_aff_y_d,
1040                                         step_aff_z_L,
1041                                         step_aff_z_U,
1042                                         step_aff_v_L,
1043                                         step_aff_v_U,
1044                                         step_cen_x_L,
1045                                         step_cen_x_U,
1046                                         step_cen_s_L,
1047                                         step_cen_s_U,
1048                                         step_cen_y_c,
1049                                         step_cen_y_d,
1050                                         step_cen_z_L,
1051                                         step_cen_z_U,
1052                                         step_cen_v_L,
1053                                         step_cen_v_U);
1054      }
1055      else {
1056        log_sigma_up = log_sigma_mid2;
1057        log_sigma_mid2 = log_sigma_mid1;
1058        qmid2 = qmid1;
1059        log_sigma_mid1 = log_sigma_lo + gfac*(log_sigma_up-log_sigma_lo);
1060        qmid1 = CalculateQualityFunction(exp(log_sigma_mid1),
1061                                         step_aff_x_L,
1062                                         step_aff_x_U,
1063                                         step_aff_s_L,
1064                                         step_aff_s_U,
1065                                         step_aff_y_c,
1066                                         step_aff_y_d,
1067                                         step_aff_z_L,
1068                                         step_aff_z_U,
1069                                         step_aff_v_L,
1070                                         step_aff_v_U,
1071                                         step_cen_x_L,
1072                                         step_cen_x_U,
1073                                         step_cen_s_L,
1074                                         step_cen_s_U,
1075                                         step_cen_y_c,
1076                                         step_cen_y_d,
1077                                         step_cen_z_L,
1078                                         step_cen_z_U,
1079                                         step_cen_v_L,
1080                                         step_cen_v_U);
1081      }
1082    }
1083
1084    Number q;
1085    if (qmid1 < qmid2) {
1086      log_sigma = log_sigma_mid1;
1087      q = qmid1;
1088    }
1089    else {
1090      log_sigma = log_sigma_mid2;
1091      q = qmid2;
1092    }
1093    if (log_sigma_up == log_sigma_up_in) {
1094      Number qtmp = CalculateQualityFunction(exp(log_sigma_up),
1095                                             step_aff_x_L,
1096                                             step_aff_x_U,
1097                                             step_aff_s_L,
1098                                             step_aff_s_U,
1099                                             step_aff_y_c,
1100                                             step_aff_y_d,
1101                                             step_aff_z_L,
1102                                             step_aff_z_U,
1103                                             step_aff_v_L,
1104                                             step_aff_v_U,
1105                                             step_cen_x_L,
1106                                             step_cen_x_U,
1107                                             step_cen_s_L,
1108                                             step_cen_s_U,
1109                                             step_cen_y_c,
1110                                             step_cen_y_d,
1111                                             step_cen_z_L,
1112                                             step_cen_z_U,
1113                                             step_cen_v_L,
1114                                             step_cen_v_U);
1115      if (qtmp < q) {
1116        log_sigma = log_sigma_up;
1117        q = qtmp;
1118      }
1119    }
1120    else if (log_sigma_lo == log_sigma_lo_in) {
1121      Number qtmp = CalculateQualityFunction(exp(log_sigma_lo),
1122                                             step_aff_x_L,
1123                                             step_aff_x_U,
1124                                             step_aff_s_L,
1125                                             step_aff_s_U,
1126                                             step_aff_y_c,
1127                                             step_aff_y_d,
1128                                             step_aff_z_L,
1129                                             step_aff_z_U,
1130                                             step_aff_v_L,
1131                                             step_aff_v_U,
1132                                             step_cen_x_L,
1133                                             step_cen_x_U,
1134                                             step_cen_s_L,
1135                                             step_cen_s_U,
1136                                             step_cen_y_c,
1137                                             step_cen_y_d,
1138                                             step_cen_z_L,
1139                                             step_cen_z_U,
1140                                             step_cen_v_L,
1141                                             step_cen_v_U);
1142      if (qtmp < q) {
1143        log_sigma = log_sigma_lo;
1144        q = qtmp;
1145      }
1146    }
1147
1148    return exp(log_sigma);
1149  }
1150  */
1151
1152
1153} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.