source: branches/dev/Algorithm/IpIpoptAlg.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: 24.4 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 529 2005-09-29 21:12:38Z 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    try {
154      IpData().TimingStats().InitializeIterates().Start();
155      // Initialize the iterates
156      InitializeIterates();
157      IpData().TimingStats().InitializeIterates().End();
158
159      if (!skip_print_problem_stats_) {
160        IpData().TimingStats().PrintProblemStatistics().Start();
161        PrintProblemStatistics();
162        IpData().TimingStats().PrintProblemStatistics().End();
163      }
164
165      IpData().TimingStats().CheckConvergence().Start();
166      ConvergenceCheck::ConvergenceStatus conv_status
167      = conv_check_->CheckConvergence();
168      IpData().TimingStats().CheckConvergence().End();
169
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    bool improve_solution = false;
300    if (IpData().HaveDeltas()) {
301      /*
302      Jnlst().Printf(J_DETAILED, J_MAIN,
303                     "No need to compute search direction - it has already been computed.\n");
304      return;
305      */
306      improve_solution = true;
307    }
308
309    SmartPtr<IteratesVector> rhs = IpData().curr()->MakeNewContainer();
310    rhs->Set_x(*IpCq().curr_grad_lag_with_damping_x());
311    rhs->Set_s(*IpCq().curr_grad_lag_with_damping_s());
312    rhs->Set_y_c(*IpCq().curr_c());
313    rhs->Set_y_d(*IpCq().curr_d_minus_s());
314    rhs->Set_z_L(*IpCq().curr_relaxed_compl_x_L());
315    rhs->Set_z_U(*IpCq().curr_relaxed_compl_x_U());
316    rhs->Set_v_L(*IpCq().curr_relaxed_compl_s_L());
317    rhs->Set_v_U(*IpCq().curr_relaxed_compl_s_U());
318
319    DBG_PRINT_VECTOR(2, "rhs", *rhs);
320
321    // Get space for the search direction
322    SmartPtr<IteratesVector> delta =
323      IpData().curr()->MakeNewIteratesVector(true);
324
325    if (improve_solution) {
326      // We can probably avoid copying and scaling...
327      delta->AddOneVector(-1., *IpData().delta(), 0.);
328    }
329
330    bool allow_inexact = false;
331    pd_solver_->Solve(-1.0, 0.0, *rhs, *delta, allow_inexact,
332                      improve_solution);
333
334    // Store the search directions in the IpData object
335    IpData().set_delta(delta);
336
337    Jnlst().Printf(J_MOREVECTOR, J_MAIN,
338                   "*** Step Calculated for Iteration: %d\n",
339                   IpData().iter_count());
340    IpData().delta()->Print(Jnlst(), J_MOREVECTOR, J_MAIN, "delta");
341  }
342
343  void IpoptAlgorithm::ComputeAcceptableTrialPoint()
344  {
345    Jnlst().Printf(J_DETAILED, J_MAIN,
346                   "\n**************************************************\n");
347    Jnlst().Printf(J_DETAILED, J_MAIN,
348                   "*** Finding Acceptable Trial Point for Iteration %d:",
349                   IpData().iter_count());
350    Jnlst().Printf(J_DETAILED, J_MAIN,
351                   "\n**************************************************\n\n");
352    line_search_->FindAcceptableTrialPoint();
353  }
354
355  void IpoptAlgorithm::OutputIteration()
356  {
357    iter_output_->WriteOutput();
358  }
359
360  void IpoptAlgorithm::InitializeIterates()
361  {
362    DBG_START_METH("IpoptAlgorithm::InitializeIterates", dbg_verbosity);
363
364    iterate_initializer_->SetInitialIterates();
365  }
366
367  void IpoptAlgorithm::AcceptTrialPoint()
368  {
369    DBG_START_METH("IpoptAlgorithm::AcceptTrialPoint", dbg_verbosity);
370    // If the line search didn't determine a new acceptable trial
371    // point, do not accept a new iterate
372    if (line_search_->CheckSkippedLineSearch()) {
373      Jnlst().Printf(J_SUMMARY, J_MAIN,
374                     "Line search didn't find acceptable trial point.\n");
375      return;
376    }
377
378    // Adjust the bounds if necessary
379    Index adjusted_slacks = IpCq().AdjustedTrialSlacks();
380    DBG_PRINT((1, "adjusted_slacks = %d\n", adjusted_slacks));
381    if (adjusted_slacks>0) {
382      IpCq().ResetAdjustedTrialSlacks();
383      if (adjusted_slacks==1) {
384        Jnlst().Printf(J_SUMMARY, J_MAIN,
385                       "%d Slack too small, adjusting variable bound\n",
386                       adjusted_slacks);
387      }
388      else {
389        Jnlst().Printf(J_SUMMARY, J_MAIN,
390                       "%d Slacks too small, adjusting variable bounds\n",
391                       adjusted_slacks);
392      }
393      if (Jnlst().ProduceOutput(J_VECTOR, J_MAIN)) {
394        IpNLP().x_L()->Print(Jnlst(), J_VECTOR, J_MAIN, "old_x_L");
395        IpNLP().x_U()->Print(Jnlst(), J_VECTOR, J_MAIN, "old_x_U");
396        IpNLP().d_L()->Print(Jnlst(), J_VECTOR, J_MAIN, "old_d_L");
397        IpNLP().d_U()->Print(Jnlst(), J_VECTOR, J_MAIN, "old_d_U");
398      }
399
400      SmartPtr<Vector> new_x_l = IpNLP().x_L()->MakeNew();
401      IpNLP().Px_L()->TransMultVector(1.0, *IpData().trial()->x(),
402                                      0.0, *new_x_l);
403      new_x_l->Axpy(-1.0, *IpCq().trial_slack_x_L());
404
405      SmartPtr<Vector> new_x_u = IpNLP().x_U()->MakeNew();
406      IpNLP().Px_U()->TransMultVector(1.0, *IpData().trial()->x(),
407                                      0.0, *new_x_u);
408      new_x_u->Axpy(1.0, *IpCq().trial_slack_x_U());
409
410      SmartPtr<Vector> new_d_l = IpNLP().d_L()->MakeNew();
411      IpNLP().Pd_L()->TransMultVector(1.0, *IpData().trial()->s(),
412                                      0.0, *new_d_l);
413      new_d_l->Axpy(-1.0, *IpCq().trial_slack_s_L());
414
415      SmartPtr<Vector> new_d_u = IpNLP().d_U()->MakeNew();
416      IpNLP().Pd_U()->TransMultVector(1.0, *IpData().trial()->s(),
417                                      0.0, *new_d_u);
418      new_d_u->Axpy(1.0, *IpCq().trial_slack_s_U());
419
420      IpNLP().AdjustVariableBounds(*new_x_l, *new_x_u, *new_d_l, *new_d_u);
421
422      if (Jnlst().ProduceOutput(J_VECTOR, J_MAIN)) {
423        IpNLP().x_L()->Print(Jnlst(), J_VECTOR, J_MAIN, "new_x_L");
424        IpNLP().x_U()->Print(Jnlst(), J_VECTOR, J_MAIN, "new_x_U");
425        IpNLP().d_L()->Print(Jnlst(), J_VECTOR, J_MAIN, "new_d_L");
426        IpNLP().d_U()->Print(Jnlst(), J_VECTOR, J_MAIN, "new_d_U");
427      }
428
429    }
430
431    // Make sure that bound multipliers are not too far from \mu * S^{-1}
432    // (see kappa_sigma in paper)
433    bool corrected = false;
434    Number max_correction;
435    SmartPtr<const Vector> new_z_L;
436    max_correction = correct_bound_multiplier(
437                       *IpData().trial()->z_L(),
438                       *IpCq().trial_slack_x_L(),
439                       *IpCq().trial_compl_x_L(),
440                       new_z_L);
441    if (max_correction>0.) {
442      Jnlst().Printf(J_DETAILED, J_MAIN,
443                     "Some value in z_L becomes too large - maximal correction = %8.2e\n",
444                     max_correction);
445      corrected = true;
446    }
447    SmartPtr<const Vector> new_z_U;
448    max_correction = correct_bound_multiplier(
449                       *IpData().trial()->z_U(),
450                       *IpCq().trial_slack_x_U(),
451                       *IpCq().trial_compl_x_U(),
452                       new_z_U);
453    if (max_correction>0.) {
454      Jnlst().Printf(J_DETAILED, J_MAIN,
455                     "Some value in z_U becomes too large - maximal correction = %8.2e\n",
456                     max_correction);
457      corrected = true;
458    }
459    SmartPtr<const Vector> new_v_L;
460    max_correction = correct_bound_multiplier(
461                       *IpData().trial()->v_L(),
462                       *IpCq().trial_slack_s_L(),
463                       *IpCq().trial_compl_s_L(),
464                       new_v_L);
465    if (max_correction>0.) {
466      Jnlst().Printf(J_DETAILED, J_MAIN,
467                     "Some value in v_L becomes too large - maximal correction = %8.2e\n",
468                     max_correction);
469      corrected = true;
470    }
471    SmartPtr<const Vector> new_v_U;
472    max_correction = correct_bound_multiplier(
473                       *IpData().trial()->v_U(),
474                       *IpCq().trial_slack_s_U(),
475                       *IpCq().trial_compl_s_U(),
476                       new_v_U);
477    if (max_correction>0.) {
478      Jnlst().Printf(J_DETAILED, J_MAIN,
479                     "Some value in v_U becomes too large - maximal correction = %8.2e\n",
480                     max_correction);
481      corrected = true;
482    }
483    SmartPtr<IteratesVector> trial = IpData().trial()->MakeNewContainer();
484    trial->Set_bound_mult(*new_z_L, *new_z_U, *new_v_L, *new_v_U);
485    IpData().set_trial(trial);
486
487    if (corrected) {
488      IpData().Append_info_string("z");
489    }
490
491    // Accept the step
492    IpData().AcceptTrialPoint();
493  }
494
495  void IpoptAlgorithm::PrintProblemStatistics()
496  {
497    if (!Jnlst().ProduceOutput(J_SUMMARY, J_STATISTICS)) {
498      // nothing to print
499      return;
500    }
501
502    SmartPtr<const Vector> x = IpData().curr()->x();
503    SmartPtr<const Vector> x_L = IpNLP().x_L();
504    SmartPtr<const Vector> x_U = IpNLP().x_U();
505    SmartPtr<const Matrix> Px_L = IpNLP().Px_L();
506    SmartPtr<const Matrix> Px_U = IpNLP().Px_U();
507
508    Index nx_tot, nx_only_lower, nx_both, nx_only_upper;
509    calc_number_of_bounds(*IpData().curr()->x(), *IpNLP().x_L(), *IpNLP().x_U(),
510                          *IpNLP().Px_L(), *IpNLP().Px_U(),
511                          nx_tot, nx_only_lower, nx_both, nx_only_upper);
512
513    Index ns_tot, ns_only_lower, ns_both, ns_only_upper;
514    calc_number_of_bounds(*IpData().curr()->s(), *IpNLP().d_L(), *IpNLP().d_U(),
515                          *IpNLP().Pd_L(), *IpNLP().Pd_U(),
516                          ns_tot, ns_only_lower, ns_both, ns_only_upper);
517
518    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
519                   "Total number of variables............................: %8d\n",nx_tot);
520    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
521                   "                     variables with only lower bounds: %8d\n",
522                   nx_only_lower);
523    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
524                   "                variables with lower and upper bounds: %8d\n",nx_both);
525    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
526                   "                     variables with only upper bounds: %8d\n",
527                   nx_only_upper);
528    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
529                   "Total number of equality constraints.................: %8d\n",
530                   IpData().curr()->y_c()->Dim());
531    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
532                   "Total number of inequality constraints...............: %8d\n",ns_tot);
533    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
534                   "        inequality constraints with only lower bounds: %8d\n",
535                   ns_only_lower);
536    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
537                   "   inequality constraints with lower and upper bounds: %8d\n",ns_both);
538    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
539                   "        inequality constraints with only upper bounds: %8d\n\n",
540                   ns_only_upper);
541  }
542
543  void IpoptAlgorithm::calc_number_of_bounds(
544    const Vector& x,
545    const Vector& x_L,
546    const Vector& x_U,
547    const Matrix& Px_L,
548    const Matrix& Px_U,
549    Index& n_tot,
550    Index& n_only_lower,
551    Index& n_both,
552    Index& n_only_upper)
553  {
554    DBG_START_METH("IpoptAlgorithm::calc_number_of_bounds",
555                   dbg_verbosity);
556
557    n_tot = x.Dim();
558
559    SmartPtr<Vector> tmpx = x.MakeNew();
560    SmartPtr<Vector> tmpxL = x_L.MakeNew();
561    SmartPtr<Vector> tmpxU = x_U.MakeNew();
562
563    tmpxL->Set(-1.);
564    tmpxU->Set(2.);
565    Px_L.MultVector(1.0, *tmpxL, 0.0, *tmpx);
566    Px_U.MultVector(1.0, *tmpxU, 1.0, *tmpx);
567    // Now, x has elements
568    //  -1 : if component has only lower bound
569    //   0 : if component has no bound
570    //   1 : if component has both lower and upper bound
571    //   2 : if component has only upper bound
572    DBG_PRINT_VECTOR(2, "x-indicator", *tmpx);
573
574    SmartPtr<Vector> tmpx0 = x.MakeNew();
575    tmpx0->Set(0.);
576
577    SmartPtr<Vector> tmpx2 = x.MakeNew();
578    tmpx2->Set(-1.0);
579    tmpx2->Axpy(1.0, *tmpx);
580    tmpx2->ElementWiseMax(*tmpx0); // tmpx2 is now 1 in those
581    // components with only upper bounds
582    n_only_upper = (Index)tmpx2->Asum();
583
584    tmpx->Axpy(-2., *tmpx2);       // now make all those entries for
585    // only upper bounds zero in tmpx
586
587    tmpx2->Copy(*tmpx);
588    tmpx2->ElementWiseMax(*tmpx0); // tmpx2 is now 1 in those
589    // components with both bounds
590    n_both = (Index)tmpx2->Asum();
591
592    tmpx->Axpy(-1., *tmpx2);
593    tmpx->ElementWiseMin(*tmpx);   // tmpx is now -1 in those with only
594    // lower bounds
595    n_only_lower = (Index)tmpx->Asum();
596
597  }
598
599  Number IpoptAlgorithm::correct_bound_multiplier(
600    const Vector& trial_z,
601    const Vector& trial_slack,
602    const Vector& trial_compl,
603    SmartPtr<const Vector>& new_trial_z)
604  {
605    DBG_START_METH("IpoptAlgorithm::CorrectBoundMultiplier",
606                   dbg_verbosity);
607
608    if (kappa_sigma_<1. || trial_z.Dim()==0) {
609      new_trial_z = &trial_z;
610      return 0.;
611    }
612
613    // We choose as barrier parameter to be used either the current
614    // algorithmic barrier parameter (if we are not in the free mode),
615    // or the average complementarity (at the trial point)
616    Number mu;
617    if (IpData().FreeMuMode()) {
618      mu = IpCq().trial_avrg_compl();
619      mu = Min(mu, 1e3);
620    }
621    else {
622      mu = IpData().curr_mu();
623    }
624    DBG_PRINT((1,"mu = %8.2e\n", mu));
625    DBG_PRINT_VECTOR(2, "trial_z", trial_z);
626
627    // First check quickly if anything need to be corrected, using the
628    // trial complementarity directly.  Here, Amax is the same as Max
629    // (and we use Amax because that can be used later)
630    if (trial_compl.Amax() <= kappa_sigma_*mu &&
631        trial_compl.Min() >= 1./kappa_sigma_*mu) {
632      new_trial_z = &trial_z;
633      return 0.;
634    }
635
636    SmartPtr<Vector> one_over_s = trial_z.MakeNew();
637    one_over_s->Copy(trial_slack);
638    one_over_s->ElementWiseReciprocal();
639
640    SmartPtr<Vector> step_z = trial_z.MakeNew();
641    step_z->AddTwoVectors(kappa_sigma_*mu, *one_over_s, -1., trial_z, 0.);
642
643    DBG_PRINT_VECTOR(2, "step_z", *step_z);
644
645    Number max_correction_up = Max(0., -step_z->Min());
646    if (max_correction_up>0.) {
647      SmartPtr<Vector> tmp = trial_z.MakeNew();
648      tmp->Set(0.);
649      step_z->ElementWiseMin(*tmp);
650      tmp->AddTwoVectors(1., trial_z, 1., *step_z, 0.);
651      new_trial_z = GetRawPtr(tmp);
652    }
653    else {
654      new_trial_z = &trial_z;
655    }
656
657    step_z->AddTwoVectors(1./kappa_sigma_*mu, *one_over_s, -1., *new_trial_z, 0.);
658
659    Number max_correction_low = Max(0., step_z->Max());
660    if (max_correction_low>0.) {
661      SmartPtr<Vector> tmp = trial_z.MakeNew();
662      tmp->Set(0.);
663      step_z->ElementWiseMax(*tmp);
664      tmp->AddTwoVectors(1., *new_trial_z, 1., *step_z, 0.);
665      new_trial_z = GetRawPtr(tmp);
666    }
667
668    DBG_PRINT_VECTOR(2, "new_trial_z", *new_trial_z);
669
670    return Max(max_correction_up, max_correction_low);
671  }
672
673} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.