source: branches/dev/Algorithm/IpPDPerturbationHandler.cpp @ 517

Last change on this file since 517 was 517, checked in by andreasw, 15 years ago
  • enabled warm start option for using same application and also for solving problems with identical structure
  • minor cosmetic changes
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.2 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 517 2005-09-13 22:51:57Z 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,
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      bool retval = get_deltas_for_wrong_inertia(delta_x, delta_s,
179                    delta_c, delta_d);
180      Jnlst().Printf(J_ERROR, J_MAIN, "Cannot determine appropriate modification of Hessian matrix.\nThis can happen if there are NaN or Inf numbers in the user-provided Hessian.\n");
181      ASSERT_EXCEPTION(retval, INTERNAL_ABORT,
182                       "get_deltas_for_wrong_inertia returns false.");
183    }
184    else {
185      delta_x = 0.;
186      delta_s = delta_x;
187    }
188
189    delta_x_curr_ = delta_x;
190    delta_s_curr_ = delta_s;
191    delta_c_curr_ = delta_c;
192    delta_d_curr_ = delta_d;
193
194    IpData().Set_info_regu_x(delta_x);
195
196    return true;
197  }
198
199  bool
200  PDPerturbationHandler::PerturbForSingularity(
201    Number& delta_x, Number& delta_s,
202    Number& delta_c, Number& delta_d)
203  {
204    DBG_START_METH("PDPerturbationHandler::PerturbForSingularity",
205                   dbg_verbosity);
206
207    bool retval;
208
209    // Check for structural degeneracy
210    if (hess_degenerate_ == NOT_YET_DETERMINED ||
211        jac_degenerate_ == NOT_YET_DETERMINED) {
212      switch (test_status_) {
213        case TEST_DELTA_C_EQ_0_DELTA_X_EQ_0:
214        DBG_ASSERT(delta_x_curr_ == 0. && delta_c_curr_ == 0.);
215        // in this case we haven't tried anything for this matrix yet
216        if (jac_degenerate_ == NOT_YET_DETERMINED) {
217          delta_d_curr_ = delta_c_curr_ = delta_cd_val_;
218          test_status_ = TEST_DELTA_C_GT_0_DELTA_X_EQ_0;
219        }
220        else {
221          DBG_ASSERT(hess_degenerate_ == NOT_YET_DETERMINED);
222          retval = get_deltas_for_wrong_inertia(delta_x, delta_s,
223                                                delta_c, delta_d);
224          ASSERT_EXCEPTION(retval, INTERNAL_ABORT,
225                           "get_deltas_for_wrong_inertia returns false.");
226          DBG_ASSERT(delta_c == 0. && delta_d == 0.);
227          test_status_ = TEST_DELTA_C_EQ_0_DELTA_X_GT_0;
228        }
229        break;
230        case TEST_DELTA_C_GT_0_DELTA_X_EQ_0:
231        DBG_ASSERT(delta_x_curr_ == 0. && delta_c_curr_ > 0.);
232        DBG_ASSERT(jac_degenerate_ == NOT_YET_DETERMINED);
233        delta_d_curr_ = delta_c_curr_ = 0.;
234        retval = get_deltas_for_wrong_inertia(delta_x, delta_s,
235                                              delta_c, delta_d);
236        ASSERT_EXCEPTION(retval, INTERNAL_ABORT,
237                         "get_deltas_for_wrong_inertia returns false.");
238        DBG_ASSERT(delta_c == 0. && delta_d == 0.);
239        test_status_ = TEST_DELTA_C_EQ_0_DELTA_X_GT_0;
240        break;
241        case TEST_DELTA_C_EQ_0_DELTA_X_GT_0:
242        DBG_ASSERT(delta_x_curr_ > 0. && delta_c_curr_ == 0.);
243        delta_d_curr_ = delta_c_curr_ = delta_cd_val_;
244        retval = get_deltas_for_wrong_inertia(delta_x, delta_s,
245                                              delta_c, delta_d);
246        ASSERT_EXCEPTION(retval, INTERNAL_ABORT,
247                         "get_deltas_for_wrong_inertia returns false.");
248        test_status_ = TEST_DELTA_C_GT_0_DELTA_X_GT_0;
249        break;
250        case TEST_DELTA_C_GT_0_DELTA_X_GT_0:
251        retval = get_deltas_for_wrong_inertia(delta_x, delta_s,
252                                              delta_c, delta_d);
253        ASSERT_EXCEPTION(retval, INTERNAL_ABORT,
254                         "get_deltas_for_wrong_inertia returns false.");
255        break;
256        case NO_TEST:
257        DBG_ASSERT(false && "we should not get here.");
258      }
259    }
260    else {
261      if (delta_c_curr_ > 0.) {
262        // If we already used a perturbation for the constraints, we do
263        // the same thing as if we were encountering negative curvature
264        retval = get_deltas_for_wrong_inertia(delta_x, delta_s,
265                                              delta_c, delta_d);
266        ASSERT_EXCEPTION(retval, INTERNAL_ABORT,
267                         "get_deltas_for_wrong_inertia returns false.");
268      }
269      else {
270        // Otherwise we now perturb the lower right corner
271        delta_d_curr_ = delta_c_curr_ = delta_cd_val_;
272
273        // ToDo - also perturb Hessian?
274        IpData().Append_info_string("L");
275      }
276    }
277
278    delta_x = delta_x_curr_;
279    delta_s = delta_s_curr_;
280    delta_c = delta_c_curr_;
281    delta_d = delta_d_curr_;
282
283    IpData().Set_info_regu_x(delta_x);
284
285    return true;
286  }
287
288  bool
289  PDPerturbationHandler::get_deltas_for_wrong_inertia(
290    Number& delta_x, Number& delta_s,
291    Number& delta_c, Number& delta_d)
292  {
293    if (delta_x_curr_ == 0.) {
294      if (delta_x_last_ == 0.) {
295        delta_x_curr_ = delta_xs_init_;
296      }
297      else {
298        delta_x_curr_ = Max(delta_xs_min_,
299                            delta_x_last_*delta_xs_dec_fact_);
300      }
301    }
302    else {
303      if (delta_x_last_ == 0. || 1e5*delta_x_last_<delta_x_curr_) {
304        delta_x_curr_ = delta_xs_first_inc_fact_*delta_x_curr_;
305      }
306      else {
307        delta_x_curr_ = delta_xs_inc_fact_*delta_x_curr_;
308      }
309    }
310    if (delta_x_curr_ > delta_xs_max_) {
311      // Give up trying to solve the linear system
312      return false;
313    }
314
315    delta_s_curr_ = delta_x_curr_;
316
317    delta_x = delta_x_curr_;
318    delta_s = delta_s_curr_;
319    delta_c = delta_c_curr_;
320    delta_d = delta_d_curr_;
321
322    IpData().Set_info_regu_x(delta_x);
323
324    return true;
325  }
326
327  bool
328  PDPerturbationHandler::PerturbForWrongInertia(
329    Number& delta_x, Number& delta_s,
330    Number& delta_c, Number& delta_d)
331  {
332    DBG_START_METH("PDPerturbationHandler::PerturbForWrongInertia",
333                   dbg_verbosity);
334
335    // Check if we can conclude that components of the system are
336    // structurally degenerate (we only get here if the most recent
337    // perturbation for a test did not result in a singular system)
338    finalize_test();
339
340    return get_deltas_for_wrong_inertia(delta_x, delta_s, delta_c, delta_d);
341  }
342
343  void
344  PDPerturbationHandler::CurrentPerturbation(
345    Number& delta_x, Number& delta_s,
346    Number& delta_c, Number& delta_d)
347  {
348    delta_x = delta_x_curr_;
349    delta_s = delta_s_curr_;
350    delta_c = delta_c_curr_;
351    delta_d = delta_d_curr_;
352  }
353
354  void
355  PDPerturbationHandler::finalize_test()
356  {
357    switch (test_status_) {
358      case NO_TEST:
359      return;
360      break;
361      case TEST_DELTA_C_EQ_0_DELTA_X_EQ_0:
362      if (hess_degenerate_ == NOT_YET_DETERMINED &&
363          jac_degenerate_ == NOT_YET_DETERMINED) {
364        hess_degenerate_ = NOT_DEGENERATE;
365        jac_degenerate_ = NOT_DEGENERATE;
366        IpData().Append_info_string("Nhj ");
367      }
368      else if (hess_degenerate_ == NOT_YET_DETERMINED) {
369        hess_degenerate_ = NOT_DEGENERATE;
370        IpData().Append_info_string("Nh ");
371      }
372      else if (jac_degenerate_ == NOT_YET_DETERMINED) {
373        jac_degenerate_ = NOT_DEGENERATE;
374        IpData().Append_info_string("Nj ");
375      }
376      break;
377      case TEST_DELTA_C_GT_0_DELTA_X_EQ_0:
378      if (hess_degenerate_ == NOT_YET_DETERMINED) {
379        hess_degenerate_ = NOT_DEGENERATE;
380        IpData().Append_info_string("Nh ");
381      }
382      if (jac_degenerate_ == NOT_YET_DETERMINED) {
383        degen_iters_++;
384        if (degen_iters_ >= degen_iters_max_) {
385          jac_degenerate_ = DEGENERATE;
386          IpData().Append_info_string("Dj ");
387        }
388        IpData().Append_info_string("L");
389      }
390      break;
391      case TEST_DELTA_C_EQ_0_DELTA_X_GT_0:
392      if (jac_degenerate_ == NOT_YET_DETERMINED) {
393        jac_degenerate_ = NOT_DEGENERATE;
394        IpData().Append_info_string("Nj ");
395      }
396      if (hess_degenerate_ == NOT_YET_DETERMINED) {
397        degen_iters_++;
398        if (degen_iters_ >= degen_iters_max_) {
399          hess_degenerate_ = DEGENERATE;
400          IpData().Append_info_string("Dh ");
401        }
402      }
403      break;
404      case TEST_DELTA_C_GT_0_DELTA_X_GT_0:
405      degen_iters_++;
406      if (degen_iters_ >= degen_iters_max_) {
407        hess_degenerate_ = DEGENERATE;
408        jac_degenerate_ = DEGENERATE;
409        IpData().Append_info_string("Dhj ");
410      }
411      IpData().Append_info_string("L");
412      break;
413    }
414  }
415
416} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.