source: branches/dev-penalty/Algorithm/IpCGPenaltyLSAcceptor.cpp @ 551

Last change on this file since 551 was 551, checked in by andreasw, 15 years ago
  • improvement to the CG-penatly implememtation
  • Property svn:keywords set to Author Date Id Revision
File size: 17.0 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: IpCGPenaltyLSAcceptor.cpp 551 2005-10-27 00:31:28Z andreasw $
6//
7// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
8//           Andreas Waechter                 IBM    2005-10-13
9//               derived file from IpFilterLineSearch.cpp
10
11#include "IpCGPenaltyLSAcceptor.hpp"
12#include "IpJournalist.hpp"
13#include "IpRestoPhase.hpp"
14#include "IpAlgTypes.hpp"
15
16#ifdef HAVE_CMATH
17# include <cmath>
18#else
19# ifdef HAVE_MATH_H
20#  include <math.h>
21# else
22#  error "don't have header file for math"
23# endif
24#endif
25
26namespace Ipopt
27{
28
29#ifdef IP_DEBUG
30  static const Index dbg_verbosity = 0;
31#endif
32
33  CGPenaltyLSAcceptor::CGPenaltyLSAcceptor(const SmartPtr<PDSystemSolver>& pd_solver)
34      :
35      pd_solver_(pd_solver)
36  {
37    DBG_START_FUN("CGPenaltyLSAcceptor::CGPenaltyLSAcceptor",
38                  dbg_verbosity);
39  }
40
41  CGPenaltyLSAcceptor::~CGPenaltyLSAcceptor()
42  {
43    DBG_START_FUN("CGPenaltyLSAcceptor::~CGPenaltyLSAcceptor()",
44                  dbg_verbosity);
45  }
46
47  void CGPenaltyLSAcceptor::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
48  {
49    roptions->AddBoundedNumberOption(
50      "eta_penalty",
51      "Relaxation factor in the Armijo condition for the penalty function.",
52      0.0, true, 0.5, true, 1e-8);
53    roptions->AddLowerBoundedNumberOption(
54      "penalty_update_infeasibility_tol",
55      "Threshold for infeasibility in penalty parameter update test.",
56      0.0, true, 1e-9,
57      "If the new constraint violation is smaller than this tolerance, the "
58      "penalty parameter is not increased.");
59    roptions->AddLowerBoundedNumberOption(
60      "eta_min",
61      "LIFENG WRITES THIS.",
62      0.0, true, 1e-2,
63      "");
64    roptions->AddLowerBoundedNumberOption(
65      "penalty_update_compl_tol",
66      "LIFENG WRITES THIS.",
67      0.0, true, 0.2,
68      "");
69    roptions->AddLowerBoundedNumberOption(
70      "chi_hat",
71      "LIFENG WRITES THIS.",
72      0.0, true, 2.,
73      "");
74    roptions->AddLowerBoundedNumberOption(
75      "chi_tilde",
76      "LIFENG WRITES THIS.",
77      0.0, true, 5.,
78      "");
79    roptions->AddLowerBoundedNumberOption(
80      "chi_cup",
81      "LIFENG WRITES THIS.",
82      0.0, true, 1.5,
83      "");
84    roptions->AddLowerBoundedNumberOption(
85      "gamma_hat",
86      "LIFENG WRITES THIS.",
87      0.0, true, 0.04,
88      "");
89    roptions->AddLowerBoundedNumberOption(
90      "gamma_tilde",
91      "LIFENG WRITES THIS.",
92      0.0, true, 4.,
93      "");
94    roptions->AddLowerBoundedNumberOption(
95      "penalty_max",
96      "LIFENG WRITES THIS.",
97      0.0, true, 1e20,
98      "");
99    roptions->AddLowerBoundedNumberOption(
100      "epsilon_c",
101      "LIFENG WRITES THIS.",
102      0.0, true, 1e-2,
103      "");
104  }
105
106  bool CGPenaltyLSAcceptor::InitializeImpl(const OptionsList& options,
107      const std::string& prefix)
108  {
109    options.GetNumericValue("eta_penalty", eta_penalty_, prefix);
110    options.GetNumericValue("penalty_update_infeasibility_tol",
111                            penalty_update_infeasibility_tol_, prefix);
112    options.GetNumericValue("eta_min", eta_min_, prefix);
113    options.GetNumericValue("penalty_update_compl_tol",
114                            penalty_update_compl_tol_, prefix);
115    options.GetNumericValue("chi_hat", chi_hat_, prefix);
116    options.GetNumericValue("chi_tilde", chi_tilde_, prefix);
117    options.GetNumericValue("chi_cup", chi_cup_, prefix);
118    options.GetNumericValue("gamma_hat", gamma_hat_, prefix);
119    options.GetNumericValue("gamma_tilde", gamma_tilde_, prefix);
120    options.GetNumericValue("penalty_max", penalty_max_, prefix);
121    options.GetNumericValue("epsilon_c", epsilon_c_, prefix);
122    // The following two option is registered by FilterLSAcceptor
123    options.GetIntegerValue("max_soc", max_soc_, prefix);
124    if (max_soc_>0) {
125      ASSERT_EXCEPTION(IsValid(pd_solver_), OPTION_INVALID,
126                       "Option \"max_soc\": This option is non-negative, but no linear solver for computing the SOC given to FilterLSAcceptor object.");
127    }
128    options.GetNumericValue("kappa_soc", kappa_soc_, prefix);
129
130    Reset();
131
132    counter_penalty_updates_ = 0;
133    curr_eta_ = -1.;
134    IpData().SetPenaltyUninitialized();
135
136    return true;
137  }
138
139  void CGPenaltyLSAcceptor::InitThisLineSearch(bool in_watchdog)
140  {
141    DBG_START_METH("CGPenaltyLSAcceptor::InitThisLineSearch",
142                   dbg_verbosity);
143
144    // Set the values for the reference point
145    if (!in_watchdog) {
146      reference_penalty_function_ = IpCq().curr_penalty_function();
147      if (IpData().HaveCgFastDeltas()) {
148        // use the fast step
149        reference_direct_deriv_penalty_function_ =
150          IpCq().curr_fast_direct_deriv_penalty_function();
151      }
152      else {
153        reference_direct_deriv_penalty_function_ =
154          IpCq().curr_direct_deriv_penalty_function();
155      }
156    }
157    else {
158      reference_penalty_function_ = watchdog_penalty_function_;
159      reference_direct_deriv_penalty_function_ =
160        watchdog_direct_deriv_penalty_function_;
161    }
162  }
163
164  bool
165  CGPenaltyLSAcceptor::CheckAcceptabilityOfTrialPoint(Number alpha_primal_test)
166  {
167    DBG_START_METH("CGPenaltyLSAcceptor::CheckAcceptabilityOfTrialPoint",
168                   dbg_verbosity);
169
170    bool accept;
171
172    Number trial_penalty_function = IpCq().trial_penalty_function();
173    DBG_ASSERT(IsFiniteNumber(trial_penalty_function));
174
175    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
176                   "Checking acceptability for trial step size alpha_primal_test=%13.6e:\n", alpha_primal_test);
177    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
178                   "  New values of penalty function     = %23.16e  (reference %23.16e):\n", trial_penalty_function, reference_penalty_function_);
179    if (Jnlst().ProduceOutput(J_DETAILED, J_LINE_SEARCH)) {
180      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
181                     "curr_barr  = %23.16e curr_inf  = %23.16e\n",
182                     IpCq().curr_barrier_obj(),
183                     IpCq().curr_primal_infeasibility(NORM_2));
184      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
185                     "trial_barr = %23.16e trial_inf = %23.16e\n",
186                     IpCq().trial_barrier_obj(),
187                     IpCq().trial_primal_infeasibility(NORM_2));
188    }
189
190    // Now check the Armijo condition
191    accept = Compare_le(trial_penalty_function-reference_penalty_function_,
192                        eta_penalty_*alpha_primal_test*reference_direct_deriv_penalty_function_,
193                        reference_penalty_function_);
194
195    return accept;
196  }
197
198  Number CGPenaltyLSAcceptor::CalculateAlphaMin()
199  {
200    // ToDo For now we just return zero
201    return 0.;
202  }
203
204  bool CGPenaltyLSAcceptor::Compare_le(Number lhs, Number rhs, Number BasVal)
205  {
206    DBG_START_FUN("CGPenaltyLSAcceptor::Compare_le",
207                  dbg_verbosity);
208    DBG_PRINT((1,"lhs = %27.16e rhs = %27.16e  BasVal = %27.16e\n",lhs,rhs,BasVal));
209
210    Number mach_eps = std::numeric_limits<Number>::epsilon();
211    return (lhs - rhs <= 10.*mach_eps*fabs(BasVal));
212  }
213
214  void CGPenaltyLSAcceptor::StartWatchDog()
215  {
216    DBG_START_FUN("CGPenaltyLSAcceptor::StartWatchDog", dbg_verbosity);
217
218    watchdog_penalty_function_ = IpCq().curr_penalty_function();
219    watchdog_direct_deriv_penalty_function_ =
220      IpCq().curr_direct_deriv_penalty_function();
221    watchdog_delta_cgpen_ = IpData().delta_cgpen();
222  }
223
224  void CGPenaltyLSAcceptor::StopWatchDog()
225  {
226    DBG_START_FUN("CGPenaltyLSAcceptor::StopWatchDog", dbg_verbosity);
227
228    reference_penalty_function_ = watchdog_penalty_function_;
229    reference_direct_deriv_penalty_function_ =
230      watchdog_direct_deriv_penalty_function_;
231    IpData().set_delta_cgpen(watchdog_delta_cgpen_);
232    watchdog_delta_cgpen_ = NULL;
233  }
234
235  void CGPenaltyLSAcceptor::Reset()
236  {
237    DBG_START_FUN("CGPenaltyLSAcceptor::Reset", dbg_verbosity);
238    /*
239    counter_penalty_updates_ = 0;
240    curr_eta_ = -1.;
241    IpData().SetPenaltyUninitialized();
242    */
243  }
244
245  bool
246  CGPenaltyLSAcceptor::TrySecondOrderCorrection(
247    Number alpha_primal_test,
248    Number& alpha_primal,
249    SmartPtr<IteratesVector>& actual_delta)
250  {
251    DBG_START_METH("CGPenaltyLSAcceptor::TrySecondOrderCorrection",
252                   dbg_verbosity);
253
254    if (max_soc_==0) {
255      return false;
256    }
257
258    bool accept = false;
259    Index count_soc = 0;
260
261    Number theta_soc_old = 0.;
262    Number theta_trial = IpCq().trial_constraint_violation();
263    Number alpha_primal_soc = alpha_primal;
264
265    // delta_y_c and delta_y_d are the steps used in the right hand
266    // side for the SOC step
267    SmartPtr<const Vector> delta_y_c = IpData().delta()->y_c();
268    SmartPtr<const Vector> delta_y_d = IpData().delta()->y_d();
269
270    SmartPtr<Vector> c_soc = IpCq().curr_c()->MakeNewCopy();
271    SmartPtr<Vector> dms_soc = IpCq().curr_d_minus_s()->MakeNewCopy();
272
273    while (count_soc<max_soc_ && !accept &&
274           (count_soc==0 || theta_trial<=kappa_soc_*theta_soc_old) ) {
275      theta_soc_old = theta_trial;
276
277      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
278                     "Trying second order correction number %d\n",
279                     count_soc+1);
280
281      // Compute SOC constraint violation
282      Number c_over_r = IpCq().curr_cg_pert_fact();
283      c_soc->AddTwoVectors(1.0, *IpCq().trial_c(),
284                           -c_over_r, *delta_y_c,
285                           alpha_primal_soc);
286      dms_soc->AddTwoVectors(1.0, *IpCq().trial_d_minus_s(),
287                             -c_over_r, *delta_y_d,
288                             alpha_primal_soc);
289
290      // Compute the SOC search direction
291      SmartPtr<IteratesVector> delta_soc =
292        actual_delta->MakeNewIteratesVector(true);
293      SmartPtr<IteratesVector> rhs = actual_delta->MakeNewContainer();
294      rhs->Set_x(*IpCq().curr_grad_lag_with_damping_x());
295      rhs->Set_s(*IpCq().curr_grad_lag_with_damping_s());
296      rhs->Set_y_c(*c_soc);
297      rhs->Set_y_d(*dms_soc);
298      rhs->Set_z_L(*IpCq().curr_relaxed_compl_x_L());
299      rhs->Set_z_U(*IpCq().curr_relaxed_compl_x_U());
300      rhs->Set_v_L(*IpCq().curr_relaxed_compl_s_L());
301      rhs->Set_v_U(*IpCq().curr_relaxed_compl_s_U());
302      pd_solver_->Solve(-1.0, 0.0, *rhs, *delta_soc, true);
303
304      // Update the delta_y_c and delta_y_d vectors in case we do
305      // additional SOC steps
306      delta_y_c = ConstPtr(delta_soc->y_c());
307      delta_y_d = ConstPtr(delta_soc->y_d());
308
309      // Compute step size
310      alpha_primal_soc =
311        IpCq().primal_frac_to_the_bound(IpData().curr_tau(),
312                                        *delta_soc->x(),
313                                        *delta_soc->s());
314
315      // Check if trial point is acceptable
316      try {
317        // Compute the primal trial point
318        IpData().SetTrialPrimalVariablesFromStep(alpha_primal_soc, *delta_soc->x(), *delta_soc->s());
319
320        // in acceptance tests, use original step size!
321        accept = CheckAcceptabilityOfTrialPoint(alpha_primal_test);
322      }
323      catch(IpoptNLP::Eval_Error& e) {
324        e.ReportException(Jnlst(), J_DETAILED);
325        Jnlst().Printf(J_WARNING, J_MAIN, "Warning: SOC step rejected due to evaluation error\n");
326        IpData().Append_info_string("e");
327        accept = false;
328        // There is no point in continuing SOC procedure
329        break;
330      }
331
332      if (accept) {
333        Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Second order correction step accepted with %d corrections.\n", count_soc+1);
334        // Accept all SOC quantities
335        alpha_primal = alpha_primal_soc;
336        actual_delta = delta_soc;
337      }
338      else {
339        count_soc++;
340        theta_trial = IpCq().trial_constraint_violation();
341      }
342    }
343    return accept;
344  }
345
346  bool
347  CGPenaltyLSAcceptor::TryCorrector(
348    Number alpha_primal_test,
349    Number& alpha_primal,
350    SmartPtr<IteratesVector>& actual_delta)
351  {
352    DBG_START_METH("CGPenaltyLSAcceptor::TryCorrector",
353                   dbg_verbosity);
354
355    return false;
356  }
357
358  char CGPenaltyLSAcceptor::UpdateForNextIteration(Number alpha_primal_test)
359  {
360    char info_alpha_primal_char;
361
362    if (curr_eta_<0.) {
363      // We need to initialize the eta tolerance
364      curr_eta_ = Max(eta_min_, Min(gamma_tilde_,
365                                    gamma_hat_*IpCq().curr_nlp_error()));
366    }
367
368    // Check if the penalty parameter is to be increased
369    Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
370                   "Starting tests for penalty parameter update:\n");
371
372    // We use the new infeasibility here...
373    Number trial_inf = IpCq().trial_primal_infeasibility(NORM_2);
374    Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
375                   "trial infeasibility = %8.2\n", trial_inf);
376    bool increase = (trial_inf >= penalty_update_infeasibility_tol_);
377    if (!increase) {
378      info_alpha_primal_char='i';
379    }
380
381    if (increase) {
382      Number max_step = Max(IpData().delta_cgpen()->x()->Amax(),
383                            IpData().delta_cgpen()->s()->Amax());
384      Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
385                     "Max norm of step = %8.2\n", max_step);
386      increase = (max_step <= curr_eta_);
387      if (!increase) {
388        info_alpha_primal_char='d';
389      }
390    }
391
392    // Lifeng: Should we use the new complementarity here?  If so, I
393    // have to restructure BacktrackingLineSearch
394    Number mu = IpData().curr_mu();
395    if (increase) {
396      Number min_compl = mu;
397      Number max_compl = mu;
398      if (IpNLP().x_L()->Dim()>0) {
399        SmartPtr<const Vector> compl_x_L = IpCq().curr_compl_x_L();
400        min_compl = Min(min_compl, compl_x_L->Min());
401        max_compl = Max(max_compl, compl_x_L->Max());
402      }
403      if (IpNLP().x_U()->Dim()>0) {
404        SmartPtr<const Vector> compl_x_U = IpCq().curr_compl_x_U();
405        min_compl = Min(min_compl, compl_x_U->Min());
406        max_compl = Max(max_compl, compl_x_U->Max());
407      }
408      if (IpNLP().d_L()->Dim()>0) {
409        SmartPtr<const Vector> compl_s_L = IpCq().curr_compl_s_L();
410        min_compl = Min(min_compl, compl_s_L->Min());
411        max_compl = Max(max_compl, compl_s_L->Max());
412      }
413      if (IpNLP().d_U()->Dim()>0) {
414        SmartPtr<const Vector> compl_s_U = IpCq().curr_compl_s_U();
415        min_compl = Min(min_compl, compl_s_U->Min());
416        max_compl = Max(max_compl, compl_s_U->Max());
417      }
418      Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
419                     "Minimal compl = %8.2\n", min_compl);
420      Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
421                     "Maximal compl = %8.2\n", max_compl);
422      increase = (min_compl >= mu*penalty_update_compl_tol_ &&
423                  max_compl <= mu/penalty_update_compl_tol_);
424      if (!increase) {
425        info_alpha_primal_char='c';
426      }
427    }
428
429    // Lifeng: Here I'm using the information from the current step
430    // and the current infeasibility
431    if (increase) {
432      SmartPtr<Vector> vec = IpData().curr()->y_c()->MakeNewCopy();
433      vec->AddTwoVectors(1., *IpData().delta_cgpen()->y_c(),
434                         -1./IpCq().curr_cg_pert_fact(), *IpCq().curr_c(),
435                         1.);
436      Number omega_test = vec->Amax();
437      Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
438                     "omega_test for c = %8.2\n", omega_test);
439      increase = (omega_test < curr_eta_);
440      if (increase) {
441        SmartPtr<Vector> vec = IpData().curr()->y_d()->MakeNewCopy();
442        vec->AddTwoVectors(1., *IpData().delta()->y_d(),
443                           -1./IpCq().curr_cg_pert_fact(), *IpCq().curr_d_minus_s(),
444                           1.);
445        omega_test = vec->Amax();
446        Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
447                       "omega_test for d = %8.2\n", omega_test);
448        increase = (omega_test < curr_eta_);
449      }
450      if (!increase) {
451        info_alpha_primal_char='m';
452      }
453    }
454
455    if (increase) {
456      // Ok, now we should increase the penalty parameter
457      counter_penalty_updates_++;
458
459      // Update the eta tolerance
460      curr_eta_ = Max(eta_min_, curr_eta_/2.);
461      Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
462                     "Updating eta to = %8.2\n", curr_eta_);
463      Number penalty = IpData().curr_penalty();
464      Number y_full_step_max;
465      SmartPtr<Vector> vec = IpData().curr()->y_c()->MakeNew();
466      vec->AddTwoVectors(1., *IpData().curr()->y_c(),
467                         1., *IpData().delta_cgpen()->y_c(), 0.);
468      y_full_step_max = vec->Amax();
469      vec = IpData().curr()->y_d()->MakeNew();
470      vec->AddTwoVectors(1., *IpData().curr()->y_d(),
471                         1., *IpData().delta_cgpen()->y_d(), 0.);
472      y_full_step_max = Max(y_full_step_max, vec->Amax());
473      if (IpCq().curr_primal_infeasibility(NORM_2) >= epsilon_c_) {
474        penalty = Max(chi_hat_*penalty, y_full_step_max + 1.);
475        info_alpha_primal_char = 'l';
476      }
477      else {
478        penalty = Max(chi_tilde_*penalty, chi_cup_*y_full_step_max);
479        info_alpha_primal_char = 's';
480      }
481      if (penalty > penalty_max_) {
482        THROW_EXCEPTION(IpoptException, "Penalty parameter becomes too large.");
483      }
484      IpData().Set_penalty(penalty);
485
486      char spen[40];
487      sprintf(spen, " penalty=%8.2e", penalty);
488      IpData().Append_info_string(spen);
489    }
490
491    return info_alpha_primal_char;
492  }
493
494  void CGPenaltyLSAcceptor::PrepareRestoPhaseStart()
495  {
496    DBG_ASSERT(false && "PrepareRestoPhaseStart not yet implemented");
497  }
498
499
500} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.