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

Last change on this file since 567 was 567, checked in by andreasw, 15 years ago

adapted for Portland compilers

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