source: branches/dev/Algorithm/IpIpoptAlg.cpp @ 526

Last change on this file since 526 was 526, checked in by andreasw, 15 years ago
  • (minor) changes in the quality function search
  • added more timing statistics
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 24.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: IpIpoptAlg.cpp 526 2005-09-27 23:03:20Z andreasw $
6//
7// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
8
9#include "IpIpoptAlg.hpp"
10#include "IpJournalist.hpp"
11#include "IpRestoPhase.hpp"
12
13namespace Ipopt
14{
15#ifdef IP_DEBUG
16  static const Index dbg_verbosity = 0;
17#endif
18
19  static bool message_printed = false;
20
21  static void print_message(const Journalist& jnlst)
22  {
23    jnlst.Printf(J_INSUPPRESSIBLE, J_MAIN,
24                 "\n******************************************************************************\n"
25                 "This program contains Ipopt, a library for large-scale nonlinear optimization.\n"
26                 " Ipopt is released as open source code under the Common Public License (CPL).\n"
27                 "         For more information visit http://projects.coin-or.org/Ipopt\n"
28                 "******************************************************************************\n\n");
29    message_printed = true;
30  }
31
32  IpoptAlgorithm::IpoptAlgorithm(const SmartPtr<PDSystemSolver>& pd_solver,
33                                 const SmartPtr<LineSearch>& line_search,
34                                 const SmartPtr<MuUpdate>& mu_update,
35                                 const SmartPtr<ConvergenceCheck>& conv_check,
36                                 const SmartPtr<IterateInitializer>& iterate_initializer,
37                                 const SmartPtr<IterationOutput>& iter_output)
38      :
39      pd_solver_(pd_solver),
40      line_search_(line_search),
41      mu_update_(mu_update),
42      conv_check_(conv_check),
43      iterate_initializer_(iterate_initializer),
44      iter_output_(iter_output)
45  {
46    DBG_START_METH("IpoptAlgorithm::IpoptAlgorithm",
47                   dbg_verbosity);
48    DBG_ASSERT(IsValid(pd_solver_));
49    DBG_ASSERT(IsValid(line_search_));
50    DBG_ASSERT(IsValid(mu_update_));
51    DBG_ASSERT(IsValid(conv_check_));
52    DBG_ASSERT(IsValid(iterate_initializer_));
53    DBG_ASSERT(IsValid(iter_output_));
54  }
55
56  IpoptAlgorithm::~IpoptAlgorithm()
57  {
58    DBG_START_METH("IpoptAlgorithm::~IpoptAlgorithm()",
59                   dbg_verbosity);
60  }
61
62  void IpoptAlgorithm::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
63  {
64    roptions->SetRegisteringCategory("Line Search");
65    roptions->AddLowerBoundedNumberOption(
66      "kappa_sigma",
67      "Factor limiting the deviation of dual variables from primal estimates.",
68      0, true, 1e10,
69      "If the dual variables deviate from their primal estimates, a correction "
70      "is performed. (See Eqn. (16) in the implementation paper.) "
71      "Setting the to value less than one disables the correction.");
72  }
73
74  bool IpoptAlgorithm::InitializeImpl(const OptionsList& options,
75                                      const std::string& prefix)
76  {
77    DBG_START_METH("IpoptAlgorithm::InitializeImpl",
78                   dbg_verbosity);
79
80    // Read the IpoptAlgorithm options
81    // Initialize the Data object
82    bool retvalue = IpData().Initialize(Jnlst(),
83                                        options, prefix);
84    ASSERT_EXCEPTION(retvalue, FAILED_INITIALIZATION,
85                     "the IpIpoptData object failed to initialize.");
86
87    // Initialize the CQ object
88    retvalue = IpCq().Initialize(Jnlst(),
89                                 options, prefix);
90    ASSERT_EXCEPTION(retvalue, FAILED_INITIALIZATION,
91                     "the IpIpoptCalculatedQuantities object failed to initialize.");
92
93    // Initialize the CQ object
94    retvalue = IpNLP().Initialize(Jnlst(),
95                                  options, prefix);
96    ASSERT_EXCEPTION(retvalue, FAILED_INITIALIZATION,
97                     "the IpIpoptNLP object failed to initialize.");
98
99    // Initialize all the strategies
100    retvalue = iterate_initializer_->Initialize(Jnlst(), IpNLP(), IpData(),
101               IpCq(), options, prefix);
102    ASSERT_EXCEPTION(retvalue, FAILED_INITIALIZATION,
103                     "the iterate_initializer strategy failed to initialize.");
104
105    retvalue = mu_update_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
106                                      options, prefix);
107    ASSERT_EXCEPTION(retvalue, FAILED_INITIALIZATION,
108                     "the mu_update strategy failed to initialize.");
109
110    retvalue = pd_solver_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
111                                      options, prefix);
112    ASSERT_EXCEPTION(retvalue, FAILED_INITIALIZATION,
113                     "the pd_solver strategy failed to initialize.");
114
115    retvalue = line_search_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
116                                        options,prefix);
117    ASSERT_EXCEPTION(retvalue, FAILED_INITIALIZATION,
118                     "the line_search strategy failed to initialize.");
119
120    retvalue = conv_check_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
121                                       options, prefix);
122    ASSERT_EXCEPTION(retvalue, FAILED_INITIALIZATION,
123                     "the conv_check strategy failed to initialize.");
124
125    retvalue = iter_output_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
126                                        options, prefix);
127    ASSERT_EXCEPTION(retvalue, FAILED_INITIALIZATION,
128                     "the iter_output strategy failed to initialize.");
129
130    options.GetNumericValue("kappa_sigma", kappa_sigma_, prefix);
131
132    if (prefix=="resto.") {
133      skip_print_problem_stats_ = true;
134    }
135    else {
136      skip_print_problem_stats_ = false;
137    }
138
139    return true;
140  }
141
142  SolverReturn IpoptAlgorithm::Optimize()
143  {
144    DBG_START_METH("IpoptAlgorithm::Optimize", dbg_verbosity);
145
146    // Start measuring CPU time
147    IpData().TimingStats().OverallAlgorithm.Start();
148
149    if (!message_printed) {
150      print_message(Jnlst());
151    }
152
153    IpData().TimingStats().InitializeIterates.Start();
154    // Initialize the iterates
155    InitializeIterates();
156    IpData().TimingStats().InitializeIterates.End();
157
158    if (!skip_print_problem_stats_) {
159      IpData().TimingStats().PrintProblemStatistics.Start();
160      PrintProblemStatistics();
161      IpData().TimingStats().PrintProblemStatistics.End();
162    }
163
164    IpData().TimingStats().CheckConvergence.Start();
165    ConvergenceCheck::ConvergenceStatus conv_status
166    = conv_check_->CheckConvergence();
167    IpData().TimingStats().CheckConvergence.End();
168
169    try {
170      // main loop
171      while (conv_status == ConvergenceCheck::CONTINUE) {
172        // Set the Hessian Matrix
173        IpData().TimingStats().ActualizeHessian.Start();
174        ActualizeHessian();
175        IpData().TimingStats().ActualizeHessian.End();
176
177        // do all the output for this iteration
178        IpData().TimingStats().OutputIteration.Start();
179        OutputIteration();
180        IpData().ResetInfo();
181        IpData().TimingStats().OutputIteration.End();
182
183        // update the barrier parameter if necessary
184        IpData().TimingStats().UpdateBarrierParameter.Start();
185        UpdateBarrierParameter();
186        IpData().TimingStats().UpdateBarrierParameter.End();
187
188        // solve the primal-dual system to get the full step
189        IpData().TimingStats().ComputeSearchDirection.Start();
190        ComputeSearchDirection();
191        IpData().TimingStats().ComputeSearchDirection.End();
192
193        // Compute the new iterate
194        IpData().TimingStats().ComputeAcceptableTrialPoint.Start();
195        ComputeAcceptableTrialPoint();
196        IpData().TimingStats().ComputeAcceptableTrialPoint.End();
197
198        // Accept the new iterate
199        IpData().TimingStats().AcceptTrialPoint.Start();
200        AcceptTrialPoint();
201        IpData().TimingStats().AcceptTrialPoint.End();
202
203        IpData().Set_iter_count(IpData().iter_count()+1);
204
205        IpData().TimingStats().CheckConvergence.Start();
206        conv_status  = conv_check_->CheckConvergence();
207        IpData().TimingStats().CheckConvergence.End();
208      }
209
210      IpData().TimingStats().OutputIteration.Start();
211      OutputIteration();
212      IpData().TimingStats().OutputIteration.End();
213
214      IpData().TimingStats().OverallAlgorithm.End();
215
216      if (conv_status == ConvergenceCheck::CONVERGED) {
217        return SUCCESS;
218      }
219      if (conv_status == ConvergenceCheck::CONVERGED_TO_ACCEPTABLE_POINT) {
220        return STOP_AT_ACCEPTABLE_POINT;
221      }
222      else if (conv_status == ConvergenceCheck::MAXITER_EXCEEDED) {
223        return MAXITER_EXCEEDED;
224      }
225    }
226    catch(TINY_STEP_DETECTED& exc) {
227      exc.ReportException(Jnlst(), J_DETAILED);
228      IpData().TimingStats().UpdateBarrierParameter.End();
229      IpData().TimingStats().OverallAlgorithm.End();
230      return STOP_AT_TINY_STEP;
231    }
232    catch(ACCEPTABLE_POINT_REACHED& exc) {
233      exc.ReportException(Jnlst(), J_DETAILED);
234      IpData().TimingStats().ComputeAcceptableTrialPoint.End();
235      IpData().TimingStats().OverallAlgorithm.End();
236      return STOP_AT_ACCEPTABLE_POINT;
237    }
238    catch(LOCALLY_INFEASIBLE& exc) {
239      exc.ReportException(Jnlst(), J_DETAILED);
240      IpData().TimingStats().ComputeAcceptableTrialPoint.EndIfStarted();
241      IpData().TimingStats().CheckConvergence.EndIfStarted();
242      IpData().TimingStats().OverallAlgorithm.End();
243      return LOCAL_INFEASIBILITY;
244    }
245    catch(RESTORATION_FAILED& exc) {
246      exc.ReportException(Jnlst(), J_DETAILED);
247      IpData().TimingStats().ComputeAcceptableTrialPoint.End();
248      IpData().TimingStats().OverallAlgorithm.End();
249      return RESTORATION_FAILURE;
250    }
251    catch(FEASIBILITY_PROBLEM_SOLVED& exc) {
252      exc.ReportException(Jnlst(), J_DETAILED);
253      IpData().TimingStats().ComputeAcceptableTrialPoint.End();
254      IpData().TimingStats().OverallAlgorithm.End();
255      return SUCCESS;
256    }
257    catch(INTERNAL_ABORT& exc) {
258      exc.ReportException(Jnlst());
259      IpData().TimingStats().OverallAlgorithm.End();
260      return INTERNAL_ERROR;
261    }
262
263    DBG_ASSERT(false && "Unknown return code in the algorithm");
264
265    IpData().TimingStats().OverallAlgorithm.End();
266    return INTERNAL_ERROR;
267  }
268
269  void IpoptAlgorithm::ActualizeHessian()
270  {
271    // At this point, just compute the exact Hessian
272    IpData().Set_W(IpCq().curr_exact_hessian());
273  }
274
275
276  void IpoptAlgorithm::UpdateBarrierParameter()
277  {
278    Jnlst().Printf(J_DETAILED, J_MAIN, "\n**************************************************\n");
279    Jnlst().Printf(J_DETAILED, J_MAIN, "*** Update Barrier Parameter for Iteration %d:", IpData().iter_count());
280    Jnlst().Printf(J_DETAILED, J_MAIN, "\n**************************************************\n\n");
281    mu_update_->UpdateBarrierParameter();
282
283    Jnlst().Printf(J_DETAILED, J_MAIN, "Barrier Parameter: %e\n", IpData().curr_mu());
284
285  }
286
287  void IpoptAlgorithm::ComputeSearchDirection()
288  {
289    DBG_START_METH("IpoptAlgorithm::ComputeSearchDirection", dbg_verbosity);
290
291    Jnlst().Printf(J_DETAILED, J_MAIN,
292                   "\n**************************************************\n");
293    Jnlst().Printf(J_DETAILED, J_MAIN,
294                   "*** Solving the Primal Dual System for Iteration %d:",
295                   IpData().iter_count());
296    Jnlst().Printf(J_DETAILED, J_MAIN,
297                   "\n**************************************************\n\n");
298
299    if (IpData().HaveDeltas()) {
300      Jnlst().Printf(J_DETAILED, J_MAIN,
301                     "No need to compute search direction - it has already been computed.\n");
302      return;
303    }
304
305    SmartPtr<IteratesVector> rhs = IpData().curr()->MakeNewContainer();
306    rhs->Set_x(*IpCq().curr_grad_lag_with_damping_x());
307    rhs->Set_s(*IpCq().curr_grad_lag_with_damping_s());
308    rhs->Set_y_c(*IpCq().curr_c());
309    rhs->Set_y_d(*IpCq().curr_d_minus_s());
310    rhs->Set_z_L(*IpCq().curr_relaxed_compl_x_L());
311    rhs->Set_z_U(*IpCq().curr_relaxed_compl_x_U());
312    rhs->Set_v_L(*IpCq().curr_relaxed_compl_s_L());
313    rhs->Set_v_U(*IpCq().curr_relaxed_compl_s_U());
314
315    DBG_PRINT_VECTOR(2, "rhs", *rhs);
316
317    // Get space for the search direction
318    SmartPtr<IteratesVector> delta = IpData().curr()->MakeNewIteratesVector(true);
319
320    pd_solver_->Solve(-1.0, 0.0, *rhs, *delta);
321
322    // Store the search directions in the IpData object
323    IpData().set_delta(delta);
324
325    Jnlst().Printf(J_MOREVECTOR, J_MAIN,
326                   "*** Step Calculated for Iteration: %d\n",
327                   IpData().iter_count());
328    IpData().delta()->Print(Jnlst(), J_MOREVECTOR, J_MAIN, "delta");
329  }
330
331  void IpoptAlgorithm::ComputeAcceptableTrialPoint()
332  {
333    Jnlst().Printf(J_DETAILED, J_MAIN,
334                   "\n**************************************************\n");
335    Jnlst().Printf(J_DETAILED, J_MAIN,
336                   "*** Finding Acceptable Trial Point for Iteration %d:",
337                   IpData().iter_count());
338    Jnlst().Printf(J_DETAILED, J_MAIN,
339                   "\n**************************************************\n\n");
340    line_search_->FindAcceptableTrialPoint();
341  }
342
343  void IpoptAlgorithm::OutputIteration()
344  {
345    iter_output_->WriteOutput();
346  }
347
348  void IpoptAlgorithm::InitializeIterates()
349  {
350    DBG_START_METH("IpoptAlgorithm::InitializeIterates", dbg_verbosity);
351
352    iterate_initializer_->SetInitialIterates();
353  }
354
355  void IpoptAlgorithm::AcceptTrialPoint()
356  {
357    DBG_START_METH("IpoptAlgorithm::AcceptTrialPoint", dbg_verbosity);
358    // If the line search didn't determine a new acceptable trial
359    // point, do not accept a new iterate
360    if (line_search_->CheckSkippedLineSearch()) {
361      Jnlst().Printf(J_SUMMARY, J_MAIN,
362                     "Line search didn't find acceptable trial point.\n");
363      return;
364    }
365
366    // Adjust the bounds if necessary
367    Index adjusted_slacks = IpCq().AdjustedTrialSlacks();
368    DBG_PRINT((1, "adjusted_slacks = %d\n", adjusted_slacks));
369    if (adjusted_slacks>0) {
370      IpCq().ResetAdjustedTrialSlacks();
371      if (adjusted_slacks==1) {
372        Jnlst().Printf(J_SUMMARY, J_MAIN,
373                       "%d Slack too small, adjusting variable bound\n",
374                       adjusted_slacks);
375      }
376      else {
377        Jnlst().Printf(J_SUMMARY, J_MAIN,
378                       "%d Slacks too small, adjusting variable bounds\n",
379                       adjusted_slacks);
380      }
381      if (Jnlst().ProduceOutput(J_VECTOR, J_MAIN)) {
382        IpNLP().x_L()->Print(Jnlst(), J_VECTOR, J_MAIN, "old_x_L");
383        IpNLP().x_U()->Print(Jnlst(), J_VECTOR, J_MAIN, "old_x_U");
384        IpNLP().d_L()->Print(Jnlst(), J_VECTOR, J_MAIN, "old_d_L");
385        IpNLP().d_U()->Print(Jnlst(), J_VECTOR, J_MAIN, "old_d_U");
386      }
387
388      SmartPtr<Vector> new_x_l = IpNLP().x_L()->MakeNew();
389      IpNLP().Px_L()->TransMultVector(1.0, *IpData().trial()->x(),
390                                      0.0, *new_x_l);
391      new_x_l->Axpy(-1.0, *IpCq().trial_slack_x_L());
392
393      SmartPtr<Vector> new_x_u = IpNLP().x_U()->MakeNew();
394      IpNLP().Px_U()->TransMultVector(1.0, *IpData().trial()->x(),
395                                      0.0, *new_x_u);
396      new_x_u->Axpy(1.0, *IpCq().trial_slack_x_U());
397
398      SmartPtr<Vector> new_d_l = IpNLP().d_L()->MakeNew();
399      IpNLP().Pd_L()->TransMultVector(1.0, *IpData().trial()->s(),
400                                      0.0, *new_d_l);
401      new_d_l->Axpy(-1.0, *IpCq().trial_slack_s_L());
402
403      SmartPtr<Vector> new_d_u = IpNLP().d_U()->MakeNew();
404      IpNLP().Pd_U()->TransMultVector(1.0, *IpData().trial()->s(),
405                                      0.0, *new_d_u);
406      new_d_u->Axpy(1.0, *IpCq().trial_slack_s_U());
407
408      IpNLP().AdjustVariableBounds(*new_x_l, *new_x_u, *new_d_l, *new_d_u);
409
410      if (Jnlst().ProduceOutput(J_VECTOR, J_MAIN)) {
411        IpNLP().x_L()->Print(Jnlst(), J_VECTOR, J_MAIN, "new_x_L");
412        IpNLP().x_U()->Print(Jnlst(), J_VECTOR, J_MAIN, "new_x_U");
413        IpNLP().d_L()->Print(Jnlst(), J_VECTOR, J_MAIN, "new_d_L");
414        IpNLP().d_U()->Print(Jnlst(), J_VECTOR, J_MAIN, "new_d_U");
415      }
416
417    }
418
419    // Make sure that bound multipliers are not too far from \mu * S^{-1}
420    // (see kappa_sigma in paper)
421    bool corrected = false;
422    Number max_correction;
423    SmartPtr<const Vector> new_z_L;
424    max_correction = correct_bound_multiplier(
425                       *IpData().trial()->z_L(),
426                       *IpCq().trial_slack_x_L(),
427                       *IpCq().trial_compl_x_L(),
428                       new_z_L);
429    if (max_correction>0.) {
430      Jnlst().Printf(J_DETAILED, J_MAIN,
431                     "Some value in z_L becomes too large - maximal correction = %8.2e\n",
432                     max_correction);
433      corrected = true;
434    }
435    SmartPtr<const Vector> new_z_U;
436    max_correction = correct_bound_multiplier(
437                       *IpData().trial()->z_U(),
438                       *IpCq().trial_slack_x_U(),
439                       *IpCq().trial_compl_x_U(),
440                       new_z_U);
441    if (max_correction>0.) {
442      Jnlst().Printf(J_DETAILED, J_MAIN,
443                     "Some value in z_U becomes too large - maximal correction = %8.2e\n",
444                     max_correction);
445      corrected = true;
446    }
447    SmartPtr<const Vector> new_v_L;
448    max_correction = correct_bound_multiplier(
449                       *IpData().trial()->v_L(),
450                       *IpCq().trial_slack_s_L(),
451                       *IpCq().trial_compl_s_L(),
452                       new_v_L);
453    if (max_correction>0.) {
454      Jnlst().Printf(J_DETAILED, J_MAIN,
455                     "Some value in v_L becomes too large - maximal correction = %8.2e\n",
456                     max_correction);
457      corrected = true;
458    }
459    SmartPtr<const Vector> new_v_U;
460    max_correction = correct_bound_multiplier(
461                       *IpData().trial()->v_U(),
462                       *IpCq().trial_slack_s_U(),
463                       *IpCq().trial_compl_s_U(),
464                       new_v_U);
465    if (max_correction>0.) {
466      Jnlst().Printf(J_DETAILED, J_MAIN,
467                     "Some value in v_U becomes too large - maximal correction = %8.2e\n",
468                     max_correction);
469      corrected = true;
470    }
471    SmartPtr<IteratesVector> trial = IpData().trial()->MakeNewContainer();
472    trial->Set_bound_mult(*new_z_L, *new_z_U, *new_v_L, *new_v_U);
473    IpData().set_trial(trial);
474
475    if (corrected) {
476      IpData().Append_info_string("z");
477    }
478
479    // Accept the step
480    IpData().AcceptTrialPoint();
481  }
482
483  void IpoptAlgorithm::PrintProblemStatistics()
484  {
485    if (!Jnlst().ProduceOutput(J_SUMMARY, J_STATISTICS)) {
486      // nothing to print
487      return;
488    }
489
490    SmartPtr<const Vector> x = IpData().curr()->x();
491    SmartPtr<const Vector> x_L = IpNLP().x_L();
492    SmartPtr<const Vector> x_U = IpNLP().x_U();
493    SmartPtr<const Matrix> Px_L = IpNLP().Px_L();
494    SmartPtr<const Matrix> Px_U = IpNLP().Px_U();
495
496    Index nx_tot, nx_only_lower, nx_both, nx_only_upper;
497    calc_number_of_bounds(*IpData().curr()->x(), *IpNLP().x_L(), *IpNLP().x_U(),
498                          *IpNLP().Px_L(), *IpNLP().Px_U(),
499                          nx_tot, nx_only_lower, nx_both, nx_only_upper);
500
501    Index ns_tot, ns_only_lower, ns_both, ns_only_upper;
502    calc_number_of_bounds(*IpData().curr()->s(), *IpNLP().d_L(), *IpNLP().d_U(),
503                          *IpNLP().Pd_L(), *IpNLP().Pd_U(),
504                          ns_tot, ns_only_lower, ns_both, ns_only_upper);
505
506    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
507                   "Total number of variables............................: %8d\n",nx_tot);
508    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
509                   "                     variables with only lower bounds: %8d\n",
510                   nx_only_lower);
511    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
512                   "                variables with lower and upper bounds: %8d\n",nx_both);
513    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
514                   "                     variables with only upper bounds: %8d\n",
515                   nx_only_upper);
516    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
517                   "Total number of equality constraints.................: %8d\n",
518                   IpData().curr()->y_c()->Dim());
519    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
520                   "Total number of inequality constraints...............: %8d\n",ns_tot);
521    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
522                   "        inequality constraints with only lower bounds: %8d\n",
523                   ns_only_lower);
524    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
525                   "   inequality constraints with lower and upper bounds: %8d\n",ns_both);
526    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
527                   "        inequality constraints with only upper bounds: %8d\n\n",
528                   ns_only_upper);
529  }
530
531  void IpoptAlgorithm::calc_number_of_bounds(
532    const Vector& x,
533    const Vector& x_L,
534    const Vector& x_U,
535    const Matrix& Px_L,
536    const Matrix& Px_U,
537    Index& n_tot,
538    Index& n_only_lower,
539    Index& n_both,
540    Index& n_only_upper)
541  {
542    DBG_START_METH("IpoptAlgorithm::calc_number_of_bounds",
543                   dbg_verbosity);
544
545    n_tot = x.Dim();
546
547    SmartPtr<Vector> tmpx = x.MakeNew();
548    SmartPtr<Vector> tmpxL = x_L.MakeNew();
549    SmartPtr<Vector> tmpxU = x_U.MakeNew();
550
551    tmpxL->Set(-1.);
552    tmpxU->Set(2.);
553    Px_L.MultVector(1.0, *tmpxL, 0.0, *tmpx);
554    Px_U.MultVector(1.0, *tmpxU, 1.0, *tmpx);
555    // Now, x has elements
556    //  -1 : if component has only lower bound
557    //   0 : if component has no bound
558    //   1 : if component has both lower and upper bound
559    //   2 : if component has only upper bound
560    DBG_PRINT_VECTOR(2, "x-indicator", *tmpx);
561
562    SmartPtr<Vector> tmpx0 = x.MakeNew();
563    tmpx0->Set(0.);
564
565    SmartPtr<Vector> tmpx2 = x.MakeNew();
566    tmpx2->Set(-1.0);
567    tmpx2->Axpy(1.0, *tmpx);
568    tmpx2->ElementWiseMax(*tmpx0); // tmpx2 is now 1 in those
569    // components with only upper bounds
570    n_only_upper = (Index)tmpx2->Asum();
571
572    tmpx->Axpy(-2., *tmpx2);       // now make all those entries for
573    // only upper bounds zero in tmpx
574
575    tmpx2->Copy(*tmpx);
576    tmpx2->ElementWiseMax(*tmpx0); // tmpx2 is now 1 in those
577    // components with both bounds
578    n_both = (Index)tmpx2->Asum();
579
580    tmpx->Axpy(-1., *tmpx2);
581    tmpx->ElementWiseMin(*tmpx);   // tmpx is now -1 in those with only
582    // lower bounds
583    n_only_lower = (Index)tmpx->Asum();
584
585  }
586
587  Number IpoptAlgorithm::correct_bound_multiplier(
588    const Vector& trial_z,
589    const Vector& trial_slack,
590    const Vector& trial_compl,
591    SmartPtr<const Vector>& new_trial_z)
592  {
593    DBG_START_METH("IpoptAlgorithm::CorrectBoundMultiplier",
594                   dbg_verbosity);
595
596    if (kappa_sigma_<1. || trial_z.Dim()==0) {
597      new_trial_z = &trial_z;
598      return 0.;
599    }
600
601    // We choose as barrier parameter to be used either the current
602    // algorithmic barrier parameter (if we are not in the free mode),
603    // or the average complementarity (at the trial point)
604    Number mu;
605    if (IpData().FreeMuMode()) {
606      mu = IpCq().trial_avrg_compl();
607      mu = Min(mu, 1e3);
608    }
609    else {
610      mu = IpData().curr_mu();
611    }
612    DBG_PRINT((1,"mu = %8.2e\n", mu));
613    DBG_PRINT_VECTOR(2, "trial_z", trial_z);
614
615    // First check quickly if anything need to be corrected, using the
616    // trial complementarity directly.  Here, Amax is the same as Max
617    // (and we use Amax because that can be used later)
618    if (trial_compl.Amax() <= kappa_sigma_*mu &&
619        trial_compl.Min() >= 1./kappa_sigma_*mu) {
620      new_trial_z = &trial_z;
621      return 0.;
622    }
623
624    SmartPtr<Vector> one_over_s = trial_z.MakeNew();
625    one_over_s->Copy(trial_slack);
626    one_over_s->ElementWiseReciprocal();
627
628    SmartPtr<Vector> step_z = trial_z.MakeNew();
629    step_z->AddTwoVectors(kappa_sigma_*mu, *one_over_s, -1., trial_z, 0.);
630
631    DBG_PRINT_VECTOR(2, "step_z", *step_z);
632
633    Number max_correction_up = Max(0., -step_z->Min());
634    if (max_correction_up>0.) {
635      SmartPtr<Vector> tmp = trial_z.MakeNew();
636      tmp->Set(0.);
637      step_z->ElementWiseMin(*tmp);
638      tmp->AddTwoVectors(1., trial_z, 1., *step_z, 0.);
639      new_trial_z = GetRawPtr(tmp);
640    }
641    else {
642      new_trial_z = &trial_z;
643    }
644
645    step_z->AddTwoVectors(1./kappa_sigma_*mu, *one_over_s, -1., *new_trial_z, 0.);
646
647    Number max_correction_low = Max(0., step_z->Max());
648    if (max_correction_low>0.) {
649      SmartPtr<Vector> tmp = trial_z.MakeNew();
650      tmp->Set(0.);
651      step_z->ElementWiseMax(*tmp);
652      tmp->AddTwoVectors(1., *new_trial_z, 1., *step_z, 0.);
653      new_trial_z = GetRawPtr(tmp);
654    }
655
656    DBG_PRINT_VECTOR(2, "new_trial_z", *new_trial_z);
657
658    return Max(max_correction_up, max_correction_low);
659  }
660
661} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.