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

Last change on this file since 529 was 529, checked in by andreasw, 15 years ago
  • in Vector class, copy cached values in Copy and update them in Scal
  • perform iterative refinement only once for adaptive strategy
  • fix bugs in PDPerturbationHandler
  • avoid some overhead in CalculateQualityFunction?
  • minor changes in timing
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.1 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: IpProbingMuOracle.cpp 529 2005-09-29 21:12:38Z andreasw $
6//
7// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
8
9#include "IpProbingMuOracle.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  ProbingMuOracle::ProbingMuOracle(const SmartPtr<PDSystemSolver>& pd_solver)
28      :
29      MuOracle(),
30      pd_solver_(pd_solver)
31  {
32    DBG_ASSERT(IsValid(pd_solver_));
33  }
34
35  ProbingMuOracle::~ProbingMuOracle()
36  {}
37
38  void ProbingMuOracle::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
39  {
40    // None to register...
41  }
42
43  bool ProbingMuOracle::InitializeImpl(const OptionsList& options,
44                                       const std::string& prefix)
45  {
46    options.GetNumericValue("sigma_max", sigma_max_, prefix);
47
48    return true;
49  }
50
51  Number ProbingMuOracle::CalculateMu(Number mu_min, Number mu_max)
52  {
53    DBG_START_METH("ProbingMuOracle::CalculateMu",
54                   dbg_verbosity);
55
56    /////////////////////////////////////
57    // Compute the affine scaling step //
58    /////////////////////////////////////
59
60    Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
61                   "Solving the Primal Dual System for the affine step\n");
62    // First get the right hand side
63    SmartPtr<IteratesVector> rhs = IpData().curr()->MakeNewContainer();
64
65    rhs->Set_x(*IpCq().curr_grad_lag_x());
66    rhs->Set_s(*IpCq().curr_grad_lag_s());
67    rhs->Set_y_c(*IpCq().curr_c());
68    rhs->Set_y_d(*IpCq().curr_d_minus_s());
69    rhs->Set_z_L(*IpCq().curr_compl_x_L());
70    rhs->Set_z_U(*IpCq().curr_compl_x_U());
71    rhs->Set_v_L(*IpCq().curr_compl_s_L());
72    rhs->Set_v_U(*IpCq().curr_compl_s_U());
73
74    // Get space for the affine scaling step
75    SmartPtr<IteratesVector> step = rhs->MakeNewIteratesVector(true);
76
77    // Now solve the primal-dual system to get the affine step.  We
78    // allow a somewhat inexact solution here
79    bool allow_inexact = true;
80    pd_solver_->Solve(-1.0, 0.0,
81                      *rhs,
82                      *step,
83                      allow_inexact
84                     );
85
86    DBG_PRINT_VECTOR(2, "step", *step);
87
88    /////////////////////////////////////////////////////////////
89    // Use Mehrotra's formula to compute the barrier parameter //
90    /////////////////////////////////////////////////////////////
91
92    // First compute the fraction-to-the-boundary step sizes
93    Number alpha_primal_aff = IpCq().primal_frac_to_the_bound(1.0,
94                              *step->x(),
95                              *step->s());
96
97    Number alpha_dual_aff = IpCq().dual_frac_to_the_bound(1.0,
98                            *step->z_L(),
99                            *step->z_U(),
100                            *step->v_L(),
101                            *step->v_U());
102
103    Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
104                   "  The affine maximal step sizes are\n"
105                   "   alpha_primal_aff = %23.16e\n"
106                   "   alpha_dual_aff = %23.16e\n",
107                   alpha_primal_aff,
108                   alpha_dual_aff);
109
110    // now compute the average complementarity at the affine step
111    // ToDo shoot for mu_min instead of 0?
112    Number mu_aff = CalculateAffineMu(alpha_primal_aff, alpha_dual_aff, *step);
113    Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
114                   "  The average complementariy at the affine step is %23.16e\n",
115                   mu_aff);
116
117    // get the current average complementarity
118    Number mu_curr = IpCq().curr_avrg_compl();
119    Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
120                   "  The average complementariy at the current point is %23.16e\n",
121                   mu_curr);
122    DBG_ASSERT(mu_curr>0.);
123
124    // Apply Mehrotra's rule
125    Number sigma = pow((mu_aff/mu_curr),3);
126    // Make sure, sigma is not too large
127    sigma = Min(sigma, sigma_max_);
128
129    Number mu = sigma*mu_curr;
130
131    // Store the affine search direction (in case it is needed in the
132    // line search for a corrector step)
133    IpData().set_delta_aff(step);
134    IpData().SetHaveAffineDeltas(true);
135
136    // DELETEME
137    char ssigma[40];
138    sprintf(ssigma, " sigma=%8.2e", sigma);
139    IpData().Append_info_string(ssigma);
140    sprintf(ssigma, " xi=%8.2e ", IpCq().curr_centrality_measure());
141    IpData().Append_info_string(ssigma);
142
143    return Max(Min(mu, mu_max), mu_min);
144  }
145
146  Number ProbingMuOracle::CalculateAffineMu
147  (
148    Number alpha_primal,
149    Number alpha_dual,
150    const IteratesVector& step)
151  {
152    // Get the current values of the slack variables and bound multipliers
153    SmartPtr<const Vector> slack_x_L = IpCq().curr_slack_x_L();
154    SmartPtr<const Vector> slack_x_U = IpCq().curr_slack_x_U();
155    SmartPtr<const Vector> slack_s_L = IpCq().curr_slack_s_L();
156    SmartPtr<const Vector> slack_s_U = IpCq().curr_slack_s_U();
157
158    SmartPtr<const Vector> z_L = IpData().curr()->z_L();
159    SmartPtr<const Vector> z_U = IpData().curr()->z_U();
160    SmartPtr<const Vector> v_L = IpData().curr()->v_L();
161    SmartPtr<const Vector> v_U = IpData().curr()->v_U();
162
163    SmartPtr<Vector> tmp_slack;
164    SmartPtr<Vector> tmp_mult;
165    SmartPtr<const Matrix> P;
166    Index ncomp = 0;
167    Number sum =0.;
168
169    // For each combination of slack and multiplier, compute the new
170    // values and their dot products.
171
172    // slack_x_L
173    if (slack_x_L->Dim()>0) {
174      ncomp += slack_x_L->Dim();
175
176      P = IpNLP().Px_L();
177      tmp_slack = slack_x_L->MakeNew();
178      tmp_slack->Copy(*slack_x_L);
179      P->TransMultVector(alpha_primal, *step.x(), 1.0, *tmp_slack);
180
181      tmp_mult = z_L->MakeNew();
182      tmp_mult->Copy(*z_L);
183      tmp_mult->Axpy(alpha_dual, *step.z_L());
184
185      sum += tmp_slack->Dot(*tmp_mult);
186    }
187
188    // slack_x_U
189    if (slack_x_U->Dim()>0) {
190      ncomp += slack_x_U->Dim();
191
192      P = IpNLP().Px_U();
193      tmp_slack = slack_x_U->MakeNew();
194      tmp_slack->Copy(*slack_x_U);
195      P->TransMultVector(-alpha_primal, *step.x(), 1.0, *tmp_slack);
196
197      tmp_mult = z_U->MakeNew();
198      tmp_mult->Copy(*z_U);
199      tmp_mult->Axpy(alpha_dual, *step.z_U());
200
201      sum += tmp_slack->Dot(*tmp_mult);
202    }
203
204    // slack_s_L
205    if (slack_s_L->Dim()>0) {
206      ncomp += slack_s_L->Dim();
207
208      P = IpNLP().Pd_L();
209      tmp_slack = slack_s_L->MakeNew();
210      tmp_slack->Copy(*slack_s_L);
211      P->TransMultVector(alpha_primal, *step.s(), 1.0, *tmp_slack);
212
213      tmp_mult = v_L->MakeNew();
214      tmp_mult->Copy(*v_L);
215      tmp_mult->Axpy(alpha_dual, *step.v_L());
216
217      sum += tmp_slack->Dot(*tmp_mult);
218    }
219
220    // slack_s_U
221    if (slack_s_U->Dim()>0) {
222      ncomp += slack_s_U->Dim();
223
224      P = IpNLP().Pd_U();
225      tmp_slack = slack_s_U->MakeNew();
226      tmp_slack->Copy(*slack_s_U);
227      P->TransMultVector(-alpha_primal, *step.s(), 1.0, *tmp_slack);
228
229      tmp_mult = v_U->MakeNew();
230      tmp_mult->Copy(*v_U);
231      tmp_mult->Axpy(alpha_dual, *step.v_U());
232
233      sum += tmp_slack->Dot(*tmp_mult);
234    }
235
236    DBG_ASSERT(ncomp>0);
237
238    return sum/((Number)ncomp);
239  }
240
241} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.