source: branches/dev/Algorithm/IpPDPerturbationHandler.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: 14.3 KB
Line 
1// Copyright (C) 2005 International Business Machines and others.
2// All Rights Reserved.
3// This code is published under the Common Public License.
4//
5// $Id: IpPDPerturbationHandler.cpp 529 2005-09-29 21:12:38Z andreasw $
6//
7// Authors:  Carl Laird, Andreas Waechter              IBM    2005-08-04
8
9#include "IpPDPerturbationHandler.hpp"
10
11namespace Ipopt
12{
13#ifdef IP_DEBUG
14  static const Index dbg_verbosity = 0;
15#endif
16
17  PDPerturbationHandler::PDPerturbationHandler()
18      :
19      delta_xs_first_inc_fact_(100.),
20      delta_xs_inc_fact_(8.),
21      delta_xs_dec_fact_(1./3.),
22      delta_xs_init_(1e-4),
23      delta_cd_val_(1e-9),
24      always_perturb_cd_(false),
25      reset_last_(false),
26      degen_iters_max_(3)
27  {}
28
29  void
30  PDPerturbationHandler::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
31  {
32    roptions->AddLowerBoundedNumberOption(
33      "max_hessian_perturbation",
34      "Maximum value of regularization parameter for handling negative curvature.",
35      0, true,
36      1e40, // ToDo make this 1e20 ?
37      "In order to guarantee that the search directions are indeed proper "
38      "descent directions, Ipopt requires that the inertia of the "
39      "(augmented) linear system for the step computation has the "
40      "correct number of negative and positive eigenvalues. The idea "
41      "is that this guides the algorithm away from maximizers and makes "
42      "Ipopt more likely converge to first order optimal points that "
43      "are minimizers. If the inertia is not correct, a multiple of the "
44      "identity matrix is added to the Hessian of the Lagrangian in the "
45      "augmented system. This parameter gives the maximum value of the "
46      "regularization parameter. If a regularization of that size is "
47      "not enough, the algorithm skips this iteration and goes to the "
48      "restoration phase. (This is delta_w^max in the implementation paper.)");
49    roptions->AddLowerBoundedNumberOption(
50      "min_hessian_perturbation",
51      "Smallest perturbation of the Hessian block.",
52      0., false, 1e-20,
53      "The size of the perturbation of the Hessian block is never selected "
54      "smaller than this value, unless no perturbation is necessary. (This "
55      "is delta_w^min in implementation paper.)");
56    roptions->AddLowerBoundedNumberOption(
57      "perturb_inc_fact_first",
58      "Increase factor for x-s perturbation for very first perturbation.",
59      1., true, 100.,
60      "The factor by which the perturbation is increased when a trial value "
61      "was not sufficient - this value is used for the computation of the "
62      "very first perturbation and allows a different value for for the first "
63      "perturbation than that used for the remaining perturbations. "
64      "(This is bar_kappa_w^+ in the implementation paper.)");
65    roptions->AddLowerBoundedNumberOption(
66      "perturb_inc_fact",
67      "Increase factor for x-s perturbation.",
68      1., true, 8.,
69      "The factor by which the perturbation is increased when a trial value "
70      "was not sufficient - this value is used for the computation of "
71      "all perturbations except for the first. "
72      "(This is kappa_w^+ in the implementation paper.)");
73    roptions->AddBoundedNumberOption(
74      "perturb_dec_fact",
75      "Decrease factor for x-s perturbation.",
76      0., true, 1., true, 1./3.,
77      "The factor by which the perturbation is decreased when a trial value "
78      "is deduced from the size of the most recent successful perturbation. "
79      "(This is kappa_w^- in the implementation paper.)");
80    roptions->AddLowerBoundedNumberOption(
81      "first_hessian_perturbation",
82      "Size of first x-s perturbation tried.",
83      0., true, 1e-4,
84      "The first value tried for the x-s perturbation in the inertia "
85      "correction scheme."
86      "(This is delta_0 in the implementation paper.)");
87    roptions->AddLowerBoundedNumberOption(
88      "jacobian_regularization",
89      "Size of the regularization for rank-deficient constraint Jacobians.",
90      0., false, 1e-9,
91      "(This is delta_c in the implementation paper - assuming that kappa_c=0.)");
92  }
93
94  bool PDPerturbationHandler::InitializeImpl(const OptionsList& options,
95      const std::string& prefix)
96  {
97    options.GetNumericValue("max_hessian_perturbation", delta_xs_max_, prefix);
98    options.GetNumericValue("min_hessian_perturbation", delta_xs_min_, prefix);
99    options.GetNumericValue("perturb_inc_fact_first", delta_xs_first_inc_fact_, prefix);
100    options.GetNumericValue("perturb_inc_fact", delta_xs_inc_fact_, prefix);
101    options.GetNumericValue("perturb_dec_fact", delta_xs_dec_fact_, prefix);
102    options.GetNumericValue("first_hessian_perturbation", delta_xs_init_, prefix);
103    options.GetNumericValue("jacobian_regularization", delta_cd_val_, prefix);
104
105    hess_degenerate_ = NOT_YET_DETERMINED;
106    jac_degenerate_ = NOT_YET_DETERMINED;
107    degen_iters_ = 0;
108
109    delta_x_curr_ = 0.;
110    delta_s_curr_ = 0.;
111    delta_c_curr_ = 0.;
112    delta_d_curr_ = 0.;
113    delta_x_last_ = 0.;
114    delta_s_last_ = 0.;
115    delta_c_last_ = 0.;
116    delta_d_last_ = 0.;
117
118    test_status_ = NO_TEST;
119
120    return true;
121  }
122
123  bool
124  PDPerturbationHandler::ConsiderNewSystem(Number& delta_x, Number& delta_s,
125      Number& delta_c, Number& delta_d)
126  {
127    DBG_START_METH("PDPerturbationHandler::ConsiderNewSystem",dbg_verbosity);
128
129    // Check if we can conclude that some components of the system are
130    // structurally degenerate
131    finalize_test();
132
133    // Store the perturbation from the previous matrix
134    if (reset_last_) {
135      delta_x_last_ = delta_x_curr_;
136      delta_s_last_ = delta_s_curr_;
137      delta_c_last_ = delta_c_curr_;
138      delta_d_last_ = delta_d_curr_;
139    }
140    else {
141      if (delta_x_curr_ > 0.) {
142        delta_x_last_ = delta_x_curr_;
143      }
144      if (delta_s_curr_ > 0.) {
145        delta_s_last_ = delta_s_curr_;
146      }
147      if (delta_c_curr_ > 0.) {
148        delta_c_last_ = delta_c_curr_;
149      }
150      if (delta_d_curr_ > 0.) {
151        delta_d_last_ = delta_d_curr_;
152      }
153    }
154
155    DBG_ASSERT((hess_degenerate_ != NOT_YET_DETERMINED ||
156                jac_degenerate_ != DEGENERATE) &&
157               (jac_degenerate_ != NOT_YET_DETERMINED ||
158                hess_degenerate_ != DEGENERATE));
159
160    if (hess_degenerate_ == NOT_YET_DETERMINED ||
161        jac_degenerate_ == NOT_YET_DETERMINED) {
162      test_status_ = TEST_DELTA_C_EQ_0_DELTA_X_EQ_0;
163    }
164    else {
165      test_status_ = NO_TEST;
166    }
167
168    if (jac_degenerate_ == DEGENERATE) {
169      delta_c = delta_c_curr_ = delta_cd_val_;
170      IpData().Append_info_string("l");
171    }
172    else {
173      delta_c = delta_c_curr_ = 0.;
174    }
175    delta_d = delta_d_curr_ = delta_c;
176
177    if (hess_degenerate_ == DEGENERATE) {
178      delta_x_curr_ = 0.;
179      delta_s_curr_ = 0.;
180      bool retval = get_deltas_for_wrong_inertia(delta_x, delta_s,
181                    delta_c, delta_d);
182      ASSERT_EXCEPTION(retval, INTERNAL_ABORT,
183                       "get_deltas_for_wrong_inertia returns false for degenerate Hessian before any further modification.");
184    }
185    else {
186      delta_x = 0.;
187      delta_s = delta_x;
188    }
189
190    delta_x_curr_ = delta_x;
191    delta_s_curr_ = delta_s;
192    delta_c_curr_ = delta_c;
193    delta_d_curr_ = delta_d;
194
195    IpData().Set_info_regu_x(delta_x);
196
197    return true;
198  }
199
200  bool
201  PDPerturbationHandler::PerturbForSingularity(
202    Number& delta_x, Number& delta_s,
203    Number& delta_c, Number& delta_d)
204  {
205    DBG_START_METH("PDPerturbationHandler::PerturbForSingularity",
206                   dbg_verbosity);
207
208    bool retval;
209
210    // Check for structural degeneracy
211    if (hess_degenerate_ == NOT_YET_DETERMINED ||
212        jac_degenerate_ == NOT_YET_DETERMINED) {
213      switch (test_status_) {
214        case TEST_DELTA_C_EQ_0_DELTA_X_EQ_0:
215        DBG_ASSERT(delta_x_curr_ == 0. && delta_c_curr_ == 0.);
216        // in this case we haven't tried anything for this matrix yet
217        if (jac_degenerate_ == NOT_YET_DETERMINED) {
218          delta_d_curr_ = delta_c_curr_ = delta_cd_val_;
219          test_status_ = TEST_DELTA_C_GT_0_DELTA_X_EQ_0;
220        }
221        else {
222          DBG_ASSERT(hess_degenerate_ == NOT_YET_DETERMINED);
223          retval = get_deltas_for_wrong_inertia(delta_x, delta_s,
224                                                delta_c, delta_d);
225          ASSERT_EXCEPTION(retval, INTERNAL_ABORT,
226                           "get_deltas_for_wrong_inertia returns false.");
227          DBG_ASSERT(delta_c == 0. && delta_d == 0.);
228          test_status_ = TEST_DELTA_C_EQ_0_DELTA_X_GT_0;
229        }
230        break;
231        case TEST_DELTA_C_GT_0_DELTA_X_EQ_0:
232        DBG_ASSERT(delta_x_curr_ == 0. && delta_c_curr_ > 0.);
233        DBG_ASSERT(jac_degenerate_ == NOT_YET_DETERMINED);
234        delta_d_curr_ = delta_c_curr_ = 0.;
235        retval = get_deltas_for_wrong_inertia(delta_x, delta_s,
236                                              delta_c, delta_d);
237        ASSERT_EXCEPTION(retval, INTERNAL_ABORT,
238                         "get_deltas_for_wrong_inertia returns false.");
239        DBG_ASSERT(delta_c == 0. && delta_d == 0.);
240        test_status_ = TEST_DELTA_C_EQ_0_DELTA_X_GT_0;
241        break;
242        case TEST_DELTA_C_EQ_0_DELTA_X_GT_0:
243        DBG_ASSERT(delta_x_curr_ > 0. && delta_c_curr_ == 0.);
244        delta_d_curr_ = delta_c_curr_ = delta_cd_val_;
245        retval = get_deltas_for_wrong_inertia(delta_x, delta_s,
246                                              delta_c, delta_d);
247        ASSERT_EXCEPTION(retval, INTERNAL_ABORT,
248                         "get_deltas_for_wrong_inertia returns false.");
249        test_status_ = TEST_DELTA_C_GT_0_DELTA_X_GT_0;
250        break;
251        case TEST_DELTA_C_GT_0_DELTA_X_GT_0:
252        retval = get_deltas_for_wrong_inertia(delta_x, delta_s,
253                                              delta_c, delta_d);
254        ASSERT_EXCEPTION(retval, INTERNAL_ABORT,
255                         "get_deltas_for_wrong_inertia returns false.");
256        break;
257        case NO_TEST:
258        DBG_ASSERT(false && "we should not get here.");
259      }
260    }
261    else {
262      if (delta_c_curr_ > 0.) {
263        // If we already used a perturbation for the constraints, we do
264        // the same thing as if we were encountering negative curvature
265        retval = get_deltas_for_wrong_inertia(delta_x, delta_s,
266                                              delta_c, delta_d);
267        if (!retval) {
268          Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
269                         "Can't get_deltas_for_wrong_inertia for delta_x_curr_ = %e and delta_c_curr_ = %e\n",
270                         delta_x_curr_, delta_c_curr_);
271          return false;
272        }
273      }
274      else {
275        // Otherwise we now perturb the lower right corner
276        delta_d_curr_ = delta_c_curr_ = delta_cd_val_;
277
278        // ToDo - also perturb Hessian?
279        IpData().Append_info_string("L");
280      }
281    }
282
283    delta_x = delta_x_curr_;
284    delta_s = delta_s_curr_;
285    delta_c = delta_c_curr_;
286    delta_d = delta_d_curr_;
287
288    IpData().Set_info_regu_x(delta_x);
289
290    return true;
291  }
292
293  bool
294  PDPerturbationHandler::get_deltas_for_wrong_inertia(
295    Number& delta_x, Number& delta_s,
296    Number& delta_c, Number& delta_d)
297  {
298    if (delta_x_curr_ == 0.) {
299      if (delta_x_last_ == 0.) {
300        delta_x_curr_ = delta_xs_init_;
301      }
302      else {
303        delta_x_curr_ = Max(delta_xs_min_,
304                            delta_x_last_*delta_xs_dec_fact_);
305      }
306    }
307    else {
308      if (delta_x_last_ == 0. || 1e5*delta_x_last_<delta_x_curr_) {
309        delta_x_curr_ = delta_xs_first_inc_fact_*delta_x_curr_;
310      }
311      else {
312        delta_x_curr_ = delta_xs_inc_fact_*delta_x_curr_;
313      }
314    }
315    if (delta_x_curr_ > delta_xs_max_) {
316      // Give up trying to solve the linear system
317      return false;
318    }
319
320    delta_s_curr_ = delta_x_curr_;
321
322    delta_x = delta_x_curr_;
323    delta_s = delta_s_curr_;
324    delta_c = delta_c_curr_;
325    delta_d = delta_d_curr_;
326
327    IpData().Set_info_regu_x(delta_x);
328
329    return true;
330  }
331
332  bool
333  PDPerturbationHandler::PerturbForWrongInertia(
334    Number& delta_x, Number& delta_s,
335    Number& delta_c, Number& delta_d)
336  {
337    DBG_START_METH("PDPerturbationHandler::PerturbForWrongInertia",
338                   dbg_verbosity);
339
340    // Check if we can conclude that components of the system are
341    // structurally degenerate (we only get here if the most recent
342    // perturbation for a test did not result in a singular system)
343    finalize_test();
344
345    return get_deltas_for_wrong_inertia(delta_x, delta_s, delta_c, delta_d);
346  }
347
348  void
349  PDPerturbationHandler::CurrentPerturbation(
350    Number& delta_x, Number& delta_s,
351    Number& delta_c, Number& delta_d)
352  {
353    delta_x = delta_x_curr_;
354    delta_s = delta_s_curr_;
355    delta_c = delta_c_curr_;
356    delta_d = delta_d_curr_;
357  }
358
359  void
360  PDPerturbationHandler::finalize_test()
361  {
362    switch (test_status_) {
363      case NO_TEST:
364      return;
365      break;
366      case TEST_DELTA_C_EQ_0_DELTA_X_EQ_0:
367      if (hess_degenerate_ == NOT_YET_DETERMINED &&
368          jac_degenerate_ == NOT_YET_DETERMINED) {
369        hess_degenerate_ = NOT_DEGENERATE;
370        jac_degenerate_ = NOT_DEGENERATE;
371        IpData().Append_info_string("Nhj ");
372      }
373      else if (hess_degenerate_ == NOT_YET_DETERMINED) {
374        hess_degenerate_ = NOT_DEGENERATE;
375        IpData().Append_info_string("Nh ");
376      }
377      else if (jac_degenerate_ == NOT_YET_DETERMINED) {
378        jac_degenerate_ = NOT_DEGENERATE;
379        IpData().Append_info_string("Nj ");
380      }
381      break;
382      case TEST_DELTA_C_GT_0_DELTA_X_EQ_0:
383      if (hess_degenerate_ == NOT_YET_DETERMINED) {
384        hess_degenerate_ = NOT_DEGENERATE;
385        IpData().Append_info_string("Nh ");
386      }
387      if (jac_degenerate_ == NOT_YET_DETERMINED) {
388        degen_iters_++;
389        if (degen_iters_ >= degen_iters_max_) {
390          jac_degenerate_ = DEGENERATE;
391          IpData().Append_info_string("Dj ");
392        }
393        IpData().Append_info_string("L");
394      }
395      break;
396      case TEST_DELTA_C_EQ_0_DELTA_X_GT_0:
397      if (jac_degenerate_ == NOT_YET_DETERMINED) {
398        jac_degenerate_ = NOT_DEGENERATE;
399        IpData().Append_info_string("Nj ");
400      }
401      if (hess_degenerate_ == NOT_YET_DETERMINED) {
402        degen_iters_++;
403        if (degen_iters_ >= degen_iters_max_) {
404          hess_degenerate_ = DEGENERATE;
405          IpData().Append_info_string("Dh ");
406        }
407      }
408      break;
409      case TEST_DELTA_C_GT_0_DELTA_X_GT_0:
410      degen_iters_++;
411      if (degen_iters_ >= degen_iters_max_) {
412        hess_degenerate_ = DEGENERATE;
413        jac_degenerate_ = DEGENERATE;
414        IpData().Append_info_string("Dhj ");
415      }
416      IpData().Append_info_string("L");
417      break;
418    }
419  }
420
421} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.