source: branches/dev/Algorithm/IpOrigIpoptNLP.cpp @ 414

Last change on this file since 414 was 414, checked in by andreasw, 14 years ago
  • several changes handling the NLP scaling
  • compute unscaled quantities in CalculatedQuantities? (such as errors)
  • print unscaled errors in output file and on console after solution
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 18.4 KB
Line 
1// Copyright (C) 2004, International Business Machines and others.
2// All Rights Reserved.
3// This code is published under the Common Public License.
4//
5// $Id: IpOrigIpoptNLP.cpp 414 2005-07-28 15:53:28Z andreasw $
6//
7// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
8
9#include "IpOrigIpoptNLP.hpp"
10#ifdef OLD_C_HEADERS
11# include <math.h>
12#else
13# include <cmath>
14#endif
15
16namespace Ipopt
17{
18  DBG_SET_VERBOSITY(0);
19
20  DefineIpoptType(OrigIpoptNLP);
21
22  OrigIpoptNLP::OrigIpoptNLP(const SmartPtr<const Journalist>& jnlst,
23                             const SmartPtr<NLP>& nlp,
24                             const SmartPtr<NLPScalingObject>& nlp_scaling)
25      :
26      IpoptNLP(nlp_scaling),
27      jnlst_(jnlst),
28      nlp_(nlp),
29      f_cache_(1),
30      grad_f_cache_(1),
31      c_cache_(1),
32      jac_c_cache_(1),
33      d_cache_(1),
34      jac_d_cache_(1),
35      h_cache_(1),
36      f_evals_(0),
37      grad_f_evals_(0),
38      c_evals_(0),
39      jac_c_evals_(0),
40      d_evals_(0),
41      jac_d_evals_(0),
42      h_evals_(0),
43      initialized_(false)
44  {}
45
46  OrigIpoptNLP::~OrigIpoptNLP()
47  {}
48
49  void OrigIpoptNLP::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
50  {
51    roptions->AddLowerBoundedNumberOption(
52      "bound_relax_factor",
53      "Factor for initial relaxation of the bounds.",
54      0, false,
55      1e-8,
56      "Before start of the optimization, the bounds given by the user are "
57      "relaxed.  This option determines the factor by how much.  If it "
58      "is set to zero, then this option is disabled.  (See Eqn.(35) in "
59      "implmentation paper.)");
60  }
61
62  bool OrigIpoptNLP::Initialize(const Journalist& jnlst,
63                                const OptionsList& options,
64                                const std::string& prefix)
65  {
66    options.GetNumericValue("bound_relax_factor", bound_relax_factor_, prefix);
67
68    if (!nlp_->ProcessOptions(options, prefix)) {
69      return false;
70    }
71
72    initialized_ = true;
73    return IpoptNLP::Initialize(jnlst, options, prefix);
74  }
75
76  bool OrigIpoptNLP::InitializeStructures(SmartPtr<Vector>& x,
77                                          bool init_x,
78                                          SmartPtr<Vector>& y_c,
79                                          bool init_y_c,
80                                          SmartPtr<Vector>& y_d,
81                                          bool init_y_d,
82                                          SmartPtr<Vector>& z_L,
83                                          bool init_z_L,
84                                          SmartPtr<Vector>& z_U,
85                                          bool init_z_U,
86                                          SmartPtr<Vector>& v_L,
87                                          SmartPtr<Vector>& v_U
88                                         )
89  {
90    DBG_ASSERT(initialized_);
91
92    bool retValue = nlp_->GetSpaces(x_space_, c_space_, d_space_,
93                                    x_l_space_, px_l_space_,
94                                    x_u_space_, px_u_space_,
95                                    d_l_space_, pd_l_space_,
96                                    d_u_space_, pd_u_space_,
97                                    jac_c_space_, jac_d_space_,
98                                    h_space_);
99
100    if (!retValue) {
101      return false;
102    }
103
104    NLP_scaling()->DetermineScaling(x_space_,
105                                    c_space_, d_space_,
106                                    jac_c_space_, jac_d_space_,
107                                    h_space_,
108                                    scaled_jac_c_space_, scaled_jac_d_space_,
109                                    scaled_h_space_);
110
111    ASSERT_EXCEPTION(x_space_->Dim() >= c_space_->Dim(), TOO_FEW_DOF,
112                     "Too few degrees of freedom!");
113
114    // cannot have any null pointers, want zero length vectors
115    // instead of null - this will later need to be changed for _h;
116    retValue = (IsValid(x_space_) && IsValid(c_space_) && IsValid(d_space_)
117                && IsValid(x_l_space_) && IsValid(px_l_space_)
118                && IsValid(x_u_space_) && IsValid(px_u_space_)
119                && IsValid(d_u_space_) && IsValid(pd_u_space_)
120                && IsValid(d_l_space_) && IsValid(pd_l_space_)
121                && IsValid(jac_c_space_) && IsValid(jac_d_space_)
122                && IsValid(h_space_)
123                && IsValid(scaled_jac_c_space_) && IsValid(scaled_jac_d_space_)
124                && IsValid(scaled_h_space_));
125
126    DBG_ASSERT(retValue && "Model cannot return null vector or matrix prototypes or spaces,"
127               " please return zero length vectors instead");
128
129    // Create the bounds structures
130    SmartPtr<Vector> x_L = x_l_space_->MakeNew();
131    SmartPtr<Matrix> Px_L = px_l_space_->MakeNew();
132    SmartPtr<Vector> x_U = x_u_space_->MakeNew();
133    SmartPtr<Matrix> Px_U = px_u_space_->MakeNew();
134    SmartPtr<Vector> d_L = d_l_space_->MakeNew();
135    SmartPtr<Matrix> Pd_L = pd_l_space_->MakeNew();
136    SmartPtr<Vector> d_U = d_u_space_->MakeNew();
137    SmartPtr<Matrix> Pd_U = pd_u_space_->MakeNew();
138
139    retValue = nlp_->GetBoundsInformation(*Px_L, *x_L, *Px_U, *x_U,
140                                          *Pd_L, *d_L, *Pd_U, *d_U);
141
142    if (!retValue) {
143      return false;
144    }
145
146    relax_bounds(-bound_relax_factor_, *x_L);
147    relax_bounds( bound_relax_factor_, *x_U);
148    relax_bounds(-bound_relax_factor_, *d_L);
149    relax_bounds( bound_relax_factor_, *d_U);
150
151    x_L_ = ConstPtr(x_L);
152    Px_L_ = ConstPtr(Px_L);
153    x_U_ = ConstPtr(x_U);
154    Px_U_ = ConstPtr(Px_U);
155    d_L_ = ConstPtr(d_L);
156    Pd_L_ = ConstPtr(Pd_L);
157    d_U_ = ConstPtr(d_U);
158    Pd_U_ = ConstPtr(Pd_U);
159
160    // now create and store the scaled bounds
161    x_L_ = NLP_scaling()->apply_vector_scaling_x_LU(*Px_L_, x_L_, *x_space_);
162    x_U_ = NLP_scaling()->apply_vector_scaling_x_LU(*Px_U_, x_U_, *x_space_);
163    d_L_ = NLP_scaling()->apply_vector_scaling_d_LU(*Pd_L_, d_L_, *d_space_);
164    d_U_ = NLP_scaling()->apply_vector_scaling_d_LU(*Pd_U_, d_U_, *d_space_);
165
166    // Create the iterates structures
167    x = x_space_->MakeNew();
168    y_c = c_space_->MakeNew();
169    y_d = d_space_->MakeNew();
170    z_L = x_l_space_->MakeNew();
171    z_U = x_u_space_->MakeNew();
172    v_L = d_l_space_->MakeNew();
173    v_U = d_u_space_->MakeNew();
174
175    retValue = nlp_->GetStartingPoint(GetRawPtr(x), init_x,
176                                      GetRawPtr(y_c), init_y_c,
177                                      GetRawPtr(y_d), init_y_d,
178                                      GetRawPtr(z_L), init_z_L,
179                                      GetRawPtr(z_U), init_z_U);
180
181    if (!retValue) {
182      return false;
183    }
184
185    Number obj_scal = NLP_scaling()->apply_obj_scaling(1.);
186    if (init_x) {
187      if (NLP_scaling()->have_x_scaling()) {
188        x = NLP_scaling()->apply_vector_scaling_x_NonConst(ConstPtr(x));
189      }
190    }
191    if (init_y_c) {
192      if (NLP_scaling()->have_c_scaling()) {
193        y_c = NLP_scaling()->unapply_vector_scaling_c_NonConst(ConstPtr(y_c));
194      }
195      if (obj_scal!=1.) {
196        y_c->Scal(obj_scal);
197      }
198    }
199    if (init_y_d) {
200      if (NLP_scaling()->have_d_scaling()) {
201        y_d = NLP_scaling()->unapply_vector_scaling_d_NonConst(ConstPtr(y_d));
202      }
203      if (obj_scal!=1.) {
204        y_d->Scal(obj_scal);
205      }
206    }
207    if (init_z_L) {
208      if (NLP_scaling()->have_x_scaling()) {
209        z_L = NLP_scaling()->apply_vector_scaling_x_LU_NonConst(*Px_L_, ConstPtr(z_L), *x_space_);
210      }
211      if (obj_scal!=1.) {
212        z_L->Scal(obj_scal);
213      }
214    }
215    if (init_z_U) {
216      if (NLP_scaling()->have_x_scaling()) {
217        z_U = NLP_scaling()->apply_vector_scaling_x_LU_NonConst(*Px_U_, ConstPtr(z_U), *x_space_);
218      }
219      if (obj_scal!=1.) {
220        z_U->Scal(obj_scal);
221      }
222    }
223
224    return true;
225  }
226
227  void
228  OrigIpoptNLP::relax_bounds(Number bound_relax_factor, Vector& bounds)
229  {
230    DBG_START_METH("OrigIpoptNLP::relax_bounds", dbg_verbosity);
231    if (bound_relax_factor!=0.) {
232      SmartPtr<Vector> tmp = bounds.MakeNew();
233      tmp->Copy(bounds);
234      tmp->ElementWiseAbs();
235      SmartPtr<Vector> ones = bounds.MakeNew();
236      ones->Set(1.);
237      tmp->ElementWiseMax(*ones);
238      DBG_PRINT((1, "bound_relax_factor = %e", bound_relax_factor));
239      DBG_PRINT_VECTOR(2, "tmp", *tmp);
240      DBG_PRINT_VECTOR(2, "bounds before", bounds);
241      bounds.Axpy(bound_relax_factor, *tmp);
242      DBG_PRINT_VECTOR(2, "bounds after", bounds);
243    }
244  }
245
246  Number OrigIpoptNLP::f(const Vector& x)
247  {
248    DBG_START_METH("OrigIpoptNLP::f", dbg_verbosity);
249    Number ret = 0.0;
250    DBG_PRINT((1, "x.Tag = %d\n", x.GetTag()));
251    if (!f_cache_.GetCachedResult1Dep(ret, &x)) {
252      f_evals_++;
253      SmartPtr<const Vector> unscaled_x = NLP_scaling()->unapply_vector_scaling_x(&x);
254      bool success = nlp_->Eval_f(*unscaled_x, ret);
255      ASSERT_EXCEPTION(success && FiniteNumber(ret), Eval_Error,
256                       "Error evaluating the objective function");
257      ret = NLP_scaling()->apply_obj_scaling(ret);
258      f_cache_.AddCachedResult1Dep(ret, &x);
259    }
260
261    return ret;
262  }
263
264  SmartPtr<const Vector> OrigIpoptNLP::grad_f(const Vector& x)
265  {
266    SmartPtr<Vector> unscaled_grad_f;
267    SmartPtr<const Vector> retValue;
268    if (!grad_f_cache_.GetCachedResult1Dep(retValue, &x)) {
269      grad_f_evals_++;
270      unscaled_grad_f = x_space_->MakeNew();
271
272      SmartPtr<const Vector> unscaled_x = NLP_scaling()->unapply_vector_scaling_x(&x);
273      bool success = nlp_->Eval_grad_f(*unscaled_x, *unscaled_grad_f);
274      ASSERT_EXCEPTION(success && FiniteNumber(unscaled_grad_f->Nrm2()),
275                       Eval_Error, "Error evaluating the gradient of the objective function");
276      retValue = NLP_scaling()->apply_grad_obj_scaling(ConstPtr(unscaled_grad_f));
277      grad_f_cache_.AddCachedResult1Dep(retValue, &x);
278    }
279
280    return retValue;
281  }
282
283  /** Equality constraint residual */
284  SmartPtr<const Vector> OrigIpoptNLP::c(const Vector& x)
285  {
286    SmartPtr<Vector> unscaled_c;
287    SmartPtr<const Vector> retValue;
288    if (!c_cache_.GetCachedResult1Dep(retValue, &x)) {
289      c_evals_++;
290      unscaled_c = c_space_->MakeNew();
291      SmartPtr<const Vector> unscaled_x = NLP_scaling()->unapply_vector_scaling_x(&x);
292      bool success = nlp_->Eval_c(*unscaled_x, *unscaled_c);
293      ASSERT_EXCEPTION(success && FiniteNumber(unscaled_c->Nrm2()),
294                       Eval_Error, "Error evaluating the equality constraints");
295      retValue = NLP_scaling()->apply_vector_scaling_c(ConstPtr(unscaled_c));
296      c_cache_.AddCachedResult1Dep(retValue, &x);
297    }
298
299    return retValue;
300  }
301
302  SmartPtr<const Vector> OrigIpoptNLP::d(const Vector& x)
303  {
304    DBG_START_METH("OrigIpoptNLP::d", dbg_verbosity);
305    SmartPtr<Vector> unscaled_d;
306    SmartPtr<const Vector> retValue;
307    if (!d_cache_.GetCachedResult1Dep(retValue, &x)) {
308      d_evals_++;
309      unscaled_d = d_space_->MakeNew();
310
311      DBG_PRINT_VECTOR(2, "scaled_x", x);
312      SmartPtr<const Vector> unscaled_x = NLP_scaling()->unapply_vector_scaling_x(&x);
313      bool success = nlp_->Eval_d(*unscaled_x, *unscaled_d);
314      DBG_PRINT_VECTOR(2, "unscaled_d", *unscaled_d);
315      ASSERT_EXCEPTION(success && FiniteNumber(unscaled_d->Nrm2()),
316                       Eval_Error, "Error evaluating the inequality constraints");
317      retValue = NLP_scaling()->apply_vector_scaling_d(ConstPtr(unscaled_d));
318      d_cache_.AddCachedResult1Dep(retValue, &x);
319    }
320
321    return retValue;
322  }
323
324  SmartPtr<const Matrix> OrigIpoptNLP::jac_c(const Vector& x)
325  {
326    SmartPtr<Matrix> unscaled_jac_c;
327    SmartPtr<const Matrix> retValue;
328    if (!jac_c_cache_.GetCachedResult1Dep(retValue, &x)) {
329      jac_c_evals_++;
330      unscaled_jac_c = jac_c_space_->MakeNew();
331
332      SmartPtr<const Vector> unscaled_x = NLP_scaling()->unapply_vector_scaling_x(&x);
333      bool success = nlp_->Eval_jac_c(*unscaled_x, *unscaled_jac_c);
334      ASSERT_EXCEPTION(success, Eval_Error, "Error evaluating the jacobian of the equality constraints");
335      retValue = NLP_scaling()->apply_jac_c_scaling(ConstPtr(unscaled_jac_c));
336      jac_c_cache_.AddCachedResult1Dep(retValue, &x);
337    }
338
339    return retValue;
340  }
341
342  SmartPtr<const Matrix> OrigIpoptNLP::jac_d(const Vector& x)
343  {
344    SmartPtr<Matrix> unscaled_jac_d;
345    SmartPtr<const Matrix> retValue;
346    if (!jac_d_cache_.GetCachedResult1Dep(retValue, &x)) {
347      jac_d_evals_++;
348      unscaled_jac_d = jac_d_space_->MakeNew();
349
350      SmartPtr<const Vector> unscaled_x = NLP_scaling()->unapply_vector_scaling_x(&x);
351      bool success = nlp_->Eval_jac_d(*unscaled_x, *unscaled_jac_d);
352      ASSERT_EXCEPTION(success, Eval_Error, "Error evaluating the jacobian of the inequality constraints");
353      retValue = NLP_scaling()->apply_jac_d_scaling(ConstPtr(unscaled_jac_d));
354      jac_d_cache_.AddCachedResult1Dep(retValue, &x);
355    }
356
357    return retValue;
358  }
359
360  SmartPtr<const SymMatrix> OrigIpoptNLP::h(const Vector& x,
361      Number obj_factor,
362      const Vector& yc,
363      const Vector& yd)
364  {
365    std::vector<const TaggedObject*> deps(3);
366    deps[0] = &x;
367    deps[1] = &yc;
368    deps[2] = &yd;
369    std::vector<Number> scalar_deps(1);
370    scalar_deps[0] = obj_factor;
371
372    SmartPtr<SymMatrix> unscaled_h;
373    SmartPtr<const SymMatrix> retValue;
374    if (!h_cache_.GetCachedResult(retValue, deps, scalar_deps)) {
375      h_evals_++;
376      unscaled_h = h_space_->MakeNewSymMatrix();
377
378      SmartPtr<const Vector> unscaled_x = NLP_scaling()->unapply_vector_scaling_x(&x);
379      SmartPtr<const Vector> unscaled_yc = NLP_scaling()->apply_vector_scaling_c(&yc);
380      SmartPtr<const Vector> unscaled_yd = NLP_scaling()->apply_vector_scaling_d(&yd);
381      Number scaled_obj_factor = NLP_scaling()->apply_obj_scaling(obj_factor);
382      bool success = nlp_->Eval_h(*unscaled_x, scaled_obj_factor, *unscaled_yc, *unscaled_yd, *unscaled_h);
383      ASSERT_EXCEPTION(success, Eval_Error, "Error evaluating the hessian of the lagrangian");
384      retValue = NLP_scaling()->apply_hessian_scaling(ConstPtr(unscaled_h));
385      h_cache_.AddCachedResult(retValue, deps, scalar_deps);
386    }
387
388    return retValue;
389  }
390
391  void OrigIpoptNLP::GetSpaces(SmartPtr<const VectorSpace>& x_space,
392                               SmartPtr<const VectorSpace>& c_space,
393                               SmartPtr<const VectorSpace>& d_space,
394                               SmartPtr<const VectorSpace>& x_l_space,
395                               SmartPtr<const MatrixSpace>& px_l_space,
396                               SmartPtr<const VectorSpace>& x_u_space,
397                               SmartPtr<const MatrixSpace>& px_u_space,
398                               SmartPtr<const VectorSpace>& d_l_space,
399                               SmartPtr<const MatrixSpace>& pd_l_space,
400                               SmartPtr<const VectorSpace>& d_u_space,
401                               SmartPtr<const MatrixSpace>& pd_u_space,
402                               SmartPtr<const MatrixSpace>& Jac_c_space,
403                               SmartPtr<const MatrixSpace>& Jac_d_space,
404                               SmartPtr<const SymMatrixSpace>& Hess_lagrangian_space)
405  {
406    // Make sure that we already have all the pointers
407    DBG_ASSERT(IsValid(x_space_) &&
408               IsValid(c_space_) &&
409               IsValid(d_space_) &&
410               IsValid(x_l_space_) &&
411               IsValid(px_l_space_) &&
412               IsValid(x_u_space_) &&
413               IsValid(px_u_space_) &&
414               IsValid(d_l_space_) &&
415               IsValid(pd_l_space_) &&
416               IsValid(d_u_space_) &&
417               IsValid(pd_u_space_) &&
418               IsValid(scaled_jac_c_space_) &&
419               IsValid(scaled_jac_d_space_) &&
420               IsValid(scaled_h_space_));
421
422    DBG_ASSERT(IsValid(NLP_scaling()));
423
424    x_space = x_space_;
425    c_space = c_space_;
426    d_space = d_space_;
427    x_l_space = x_l_space_;
428    px_l_space = px_l_space_;
429    x_u_space = x_u_space_;
430    px_u_space = px_u_space_;
431    d_l_space = d_l_space_;
432    pd_l_space = pd_l_space_;
433    d_u_space = d_u_space_;
434    pd_u_space = pd_u_space_;
435    Jac_c_space = scaled_jac_c_space_;
436    Jac_d_space = scaled_jac_d_space_;
437    Hess_lagrangian_space = scaled_h_space_;
438  }
439
440  void OrigIpoptNLP::FinalizeSolution(SolverReturn status,
441                                      const Vector& x, const Vector& z_L, const Vector& z_U,
442                                      const Vector& c, const Vector& d,
443                                      const Vector& y_c, const Vector& y_d,
444                                      Number obj_value)
445  {
446    // need to submit the unscaled solution back to the nlp
447    SmartPtr<const Vector> unscaled_x =
448      NLP_scaling()->unapply_vector_scaling_x(&x);
449    SmartPtr<const Vector> unscaled_c =
450      NLP_scaling()->unapply_vector_scaling_c(&c);
451    SmartPtr<const Vector> unscaled_d =
452      NLP_scaling()->unapply_vector_scaling_d(&d);
453    const Number unscaled_obj = NLP_scaling()->unapply_obj_scaling(obj_value);
454
455    SmartPtr<const Vector> unscaled_z_L;
456    SmartPtr<const Vector> unscaled_z_U;
457    SmartPtr<const Vector> unscaled_y_c;
458    SmartPtr<const Vector> unscaled_y_d;
459
460    // The objective function scaling factor also appears in the constraints
461    Number obj_unscale_factor = NLP_scaling()->unapply_obj_scaling(1.);
462    if (obj_unscale_factor!=1.) {
463      SmartPtr<Vector> tmp = NLP_scaling()->apply_vector_scaling_x_LU_NonConst(*Px_L_, &z_L, *x_space_);
464      tmp->Scal(obj_unscale_factor);
465      unscaled_z_L = ConstPtr(tmp);
466
467      tmp = NLP_scaling()->apply_vector_scaling_x_LU_NonConst(*Px_U_, &z_U, *x_space_);
468      tmp->Scal(obj_unscale_factor);
469      unscaled_z_U = ConstPtr(tmp);
470
471      tmp = NLP_scaling()->apply_vector_scaling_c_NonConst(&y_c);
472      tmp->Scal(obj_unscale_factor);
473      unscaled_y_c = ConstPtr(tmp);
474
475      tmp = NLP_scaling()->apply_vector_scaling_d_NonConst(&y_d);
476      tmp->Scal(obj_unscale_factor);
477      unscaled_y_d = ConstPtr(tmp);
478    }
479    else {
480      unscaled_z_L = NLP_scaling()->apply_vector_scaling_x_LU(*Px_L_, &z_L, *x_space_);
481      unscaled_z_U = NLP_scaling()->apply_vector_scaling_x_LU(*Px_U_, &z_U, *x_space_);
482      unscaled_y_c = NLP_scaling()->apply_vector_scaling_c(&y_c);
483      unscaled_y_d = NLP_scaling()->apply_vector_scaling_d(&y_d);
484    }
485
486    nlp_->FinalizeSolution(status, *unscaled_x,
487                           *unscaled_z_L, *unscaled_z_U,
488                           *unscaled_c, *unscaled_d,
489                           *unscaled_y_c, *unscaled_y_d,
490                           unscaled_obj);
491  }
492
493  void OrigIpoptNLP::AdjustVariableBounds(const Vector& new_x_L, const Vector& new_x_U,
494                                          const Vector& new_d_L, const Vector& new_d_U)
495  {
496    x_L_ = new_x_L.MakeNewCopy();
497    x_U_ = new_x_U.MakeNewCopy();
498    d_L_ = new_d_L.MakeNewCopy();
499    d_U_ = new_d_U.MakeNewCopy();
500  }
501
502} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.