source: branches/dev/Algorithm/IpRestoIpoptNLP.cpp @ 416

Last change on this file since 416 was 416, checked in by andreasw, 14 years ago
  • revised handling of "acceptable level of accuracy" (now in ConvCheck?)
  • fixed uncaught evaluation error exceptions
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 26.8 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: IpRestoIpoptNLP.cpp 416 2005-07-29 19:11:41Z andreasw $
6//
7// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
8
9#include "IpRestoIpoptNLP.hpp"
10#include "IpIdentityMatrix.hpp"
11#include "IpSumSymMatrix.hpp"
12#include "IpSumMatrix.hpp"
13#include "IpNLPScaling.hpp"
14
15#ifdef OLD_C_HEADERS
16#include <math.h>
17#else
18#include <cmath>
19#endif
20
21namespace Ipopt
22{
23  DBG_SET_VERBOSITY(0);
24
25  DefineIpoptType(RestoIpoptNLP);
26
27  RestoIpoptNLP::RestoIpoptNLP(IpoptNLP& orig_ip_nlp,
28                               IpoptData& orig_ip_data,
29                               IpoptCalculatedQuantities& orig_ip_cq)
30      :
31      IpoptNLP(new NoNLPScalingObject()),
32      orig_ip_nlp_(&orig_ip_nlp),
33      orig_ip_data_(&orig_ip_data),
34      orig_ip_cq_(&orig_ip_cq),
35      eta_factor_(1.0),
36      eta_mu_exponent_(0.5),
37      rho_(1000.)
38  {}
39
40  RestoIpoptNLP::~RestoIpoptNLP()
41  {}
42
43  void RestoIpoptNLP::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
44  {
45    roptions->AddStringOption2(
46      "evaluate_orig_obj_at_resto_trial",
47      "Determines if the original objective function should be evalutated at restoration phase trial points.",
48      "yes",
49      "no", "skip evaluation",
50      "yes", "evaluate at every trial point",
51      "Setting this option to true makes the restoration phase algorithm "
52      "evaluate the objective function of the original problem at every trial "
53      "point encountered during the restoration phase, even if this value is "
54      "not required.  In this way, it is guaranteed that the original "
55      "objective function can be evaluated without error at all accepted "
56      "iterates; otherwise the algorithm might fail at a point where the "
57      "restoration phase accepts an iterate that is good for the restoration "
58      "phase problem, but not the original problem.  On the other hand, if "
59      "the evaluation of the original objective is expensive, this might be "
60      "costly.");
61  }
62
63  bool RestoIpoptNLP::Initialize(const Journalist& jnlst,
64                                 const OptionsList& options,
65                                 const std::string& prefix)
66  {
67    options.GetBoolValue("evaluate_orig_obj_at_resto_trial",
68                         evaluate_orig_obj_at_resto_trial_, prefix);
69
70    initialized_ = true;
71    return IpoptNLP::Initialize(jnlst, options, prefix);
72  }
73
74  bool RestoIpoptNLP::InitializeStructures(SmartPtr<Vector>& x,
75      bool init_x,
76      SmartPtr<Vector>& y_c,
77      bool init_y_c,
78      SmartPtr<Vector>& y_d,
79      bool init_y_d,
80      SmartPtr<Vector>& z_L,
81      bool init_z_L,
82      SmartPtr<Vector>& z_U,
83      bool init_z_U,
84      SmartPtr<Vector>& v_L,
85      SmartPtr<Vector>& v_U
86                                          )
87  {
88    DBG_START_METH("RestoIpoptNLP::InitializeStructures", 0);
89    DBG_ASSERT(initialized_);
90    ///////////////////////////////////////////////////////////
91    // Get the vector/matrix spaces for the original problem //
92    ///////////////////////////////////////////////////////////
93
94    SmartPtr<const VectorSpace> orig_x_space;
95    SmartPtr<const VectorSpace> orig_c_space;
96    SmartPtr<const VectorSpace> orig_d_space;
97    SmartPtr<const VectorSpace> orig_x_l_space;
98    SmartPtr<const MatrixSpace> orig_px_l_space;
99    SmartPtr<const VectorSpace> orig_x_u_space;
100    SmartPtr<const MatrixSpace> orig_px_u_space;
101    SmartPtr<const VectorSpace> orig_d_l_space;
102    SmartPtr<const MatrixSpace> orig_pd_l_space;
103    SmartPtr<const VectorSpace> orig_d_u_space;
104    SmartPtr<const MatrixSpace> orig_pd_u_space;
105    SmartPtr<const MatrixSpace> orig_jac_c_space;
106    SmartPtr<const MatrixSpace> orig_jac_d_space;
107    SmartPtr<const SymMatrixSpace> orig_h_space;
108
109    orig_ip_nlp_->GetSpaces(orig_x_space, orig_c_space, orig_d_space,
110                            orig_x_l_space, orig_px_l_space,
111                            orig_x_u_space, orig_px_u_space,
112                            orig_d_l_space, orig_pd_l_space,
113                            orig_d_u_space, orig_pd_u_space,
114                            orig_jac_c_space, orig_jac_d_space,
115                            orig_h_space);
116
117    // Create the restoration phase problem vector/matrix spaces, based
118    // on the original spaces (pretty inconvenient with all the
119    // matrix spaces, isn't it?!?)
120
121    DBG_PRINT((1, "Creating the x_space_\n"));
122    // vector x
123    Index total_dim = orig_x_space->Dim() + 2*orig_c_space->Dim()
124                      //orig  + orig_d_l_space->Dim() + orig_d_u_space->Dim();
125                      + 2*orig_d_space->Dim();
126    x_space_ = new CompoundVectorSpace(5, total_dim);
127    x_space_->SetCompSpace(0, *orig_x_space);
128    x_space_->SetCompSpace(1, *orig_c_space); // n_c
129    x_space_->SetCompSpace(2, *orig_c_space); // p_c
130    x_space_->SetCompSpace(3, *orig_d_space); // n_d
131    x_space_->SetCompSpace(4, *orig_d_space); // p_d
132    //orig    x_space_->SetCompSpace(3, *orig_d_l_space); // n_d
133    //orig    x_space_->SetCompSpace(4, *orig_d_u_space); // p_d
134
135    DBG_PRINT((1, "Setting the c_space_\n"));
136    // vector c
137    c_space_ = orig_c_space;
138
139    DBG_PRINT((1, "Setting the d_space_\n"));
140    // vector d
141    d_space_ = orig_d_space;
142
143    DBG_PRINT((1, "Creating the x_l_space_\n"));
144    // vector x_L
145    total_dim = orig_x_l_space->Dim() + 2*orig_c_space->Dim()
146                //orig      + orig_d_l_space->Dim() + orig_d_u_space->Dim();
147                + 2*orig_d_space->Dim();
148    x_l_space_ = new CompoundVectorSpace(5, total_dim);
149    x_l_space_->SetCompSpace(0, *orig_x_l_space);
150    x_l_space_->SetCompSpace(1, *orig_c_space); // n_c >=0
151    x_l_space_->SetCompSpace(2, *orig_c_space); // p_c >=0
152    x_l_space_->SetCompSpace(3, *orig_d_space); // n_d >=0
153    x_l_space_->SetCompSpace(4, *orig_d_space); // p_d >=0
154    //orig    x_l_space_->SetCompSpace(3, *orig_d_l_space); // n_d >=0
155    //orig    x_l_space_->SetCompSpace(4, *orig_d_u_space); // p_d >=0
156
157    DBG_PRINT((1, "Setting the x_u_space_\n"));
158    // vector x_U
159    x_u_space_ = orig_x_u_space;
160
161    DBG_PRINT((1, "Creating the px_l_space_\n"));
162    // matrix px_l
163    Index total_rows = orig_x_space->Dim() + 2*orig_c_space->Dim()
164                       //orig      + orig_d_l_space->Dim() + orig_d_u_space->Dim();
165                       + 2*orig_d_space->Dim();
166    Index total_cols = orig_x_l_space->Dim() + 2*orig_c_space->Dim()
167                       //orig      + orig_d_l_space->Dim() + orig_d_u_space->Dim();
168                       + 2*orig_d_space->Dim();
169    px_l_space_ = new CompoundMatrixSpace(5, 5, total_rows, total_cols);
170    px_l_space_->SetBlockRows(0, orig_x_space->Dim());
171    px_l_space_->SetBlockRows(1, orig_c_space->Dim());
172    px_l_space_->SetBlockRows(2, orig_c_space->Dim());
173    px_l_space_->SetBlockRows(3, orig_d_space->Dim());
174    px_l_space_->SetBlockRows(4, orig_d_space->Dim());
175    //orig    px_l_space_->SetBlockRows(3, orig_d_l_space->Dim());
176    //orig    px_l_space_->SetBlockRows(4, orig_d_u_space->Dim());
177    px_l_space_->SetBlockCols(0, orig_x_l_space->Dim());
178    px_l_space_->SetBlockCols(1, orig_c_space->Dim());
179    px_l_space_->SetBlockCols(2, orig_c_space->Dim());
180    px_l_space_->SetBlockCols(3, orig_d_space->Dim());
181    px_l_space_->SetBlockCols(4, orig_d_space->Dim());
182    //orig    px_l_space_->SetBlockCols(3, orig_d_l_space->Dim());
183    //orig    px_l_space_->SetBlockCols(4, orig_d_u_space->Dim());
184
185    px_l_space_->SetCompSpace(0, 0, *orig_px_l_space);
186    // now setup the identity matrix
187    // This could be changed to be something like...
188    // px_l_space_->SetBlockToIdentity(1,1,1.0);
189    // px_l_space_->SetBlockToIdentity(2,2,other_factor);
190    // ... etc with some simple changes to the CompoundMatrixSpace
191    // to allow this (space should auto create the matrices)
192    //
193    // for now, we use the new feature and set the true flag for this block
194    // to say that the matrices should be auto_allocated
195    SmartPtr<const MatrixSpace> identity_mat_space_nc
196    = new IdentityMatrixSpace(orig_c_space->Dim());
197    px_l_space_->SetCompSpace(1, 1, *identity_mat_space_nc, true);
198    px_l_space_->SetCompSpace(2, 2, *identity_mat_space_nc, true);
199    //orig    SmartPtr<const MatrixSpace> identity_mat_space_nd_l
200    //orig    = new IdentityMatrixSpace(orig_d_l_space->Dim());
201    SmartPtr<const MatrixSpace> identity_mat_space_nd
202    = new IdentityMatrixSpace(orig_d_space->Dim());
203    //orig    px_l_space_->SetCompSpace(3, 3, *identity_mat_space_nd_l, true);
204    px_l_space_->SetCompSpace(3, 3, *identity_mat_space_nd, true);
205    //orig    SmartPtr<const MatrixSpace> identity_mat_space_nd_u
206    //orig      = new IdentityMatrixSpace(orig_d_u_space->Dim());
207    //orig    px_l_space_->SetCompSpace(4, 4, *identity_mat_space_nd_u, true);
208    px_l_space_->SetCompSpace(4, 4, *identity_mat_space_nd, true);
209
210    DBG_PRINT((1, "Creating the px_u_space_\n"));
211    // matrix px_u    px_u_space_->SetBlockRows(0, orig_x_space->Dim());
212
213    total_rows = orig_x_space->Dim() + 2*orig_c_space->Dim()
214                 //orig      + orig_d_l_space->Dim() + orig_d_u_space->Dim();
215                 + 2*orig_d_space->Dim();
216    total_cols = orig_x_u_space->Dim();
217    DBG_PRINT((1, "total_rows = %d, total_cols = %d\n",total_rows, total_cols));
218    px_u_space_ = new CompoundMatrixSpace(5, 1, total_rows, total_cols);
219    px_u_space_->SetBlockRows(0, orig_x_space->Dim());
220    px_u_space_->SetBlockRows(1, orig_c_space->Dim());
221    px_u_space_->SetBlockRows(2, orig_c_space->Dim());
222    px_u_space_->SetBlockRows(3, orig_d_space->Dim());
223    px_u_space_->SetBlockRows(4, orig_d_space->Dim());
224    //orig    px_u_space_->SetBlockRows(3, orig_d_l_space->Dim());
225    //orig    px_u_space_->SetBlockRows(4, orig_d_u_space->Dim());
226    px_u_space_->SetBlockCols(0, orig_x_u_space->Dim());
227
228    px_u_space_->SetCompSpace(0, 0, *orig_px_u_space);
229    // other matrices are zero'ed out
230
231    // vector d_L
232    d_l_space_ = orig_d_l_space;
233
234    // vector d_U
235    d_u_space_ = orig_d_u_space;
236
237    // matrix pd_L
238    pd_l_space_ = orig_pd_l_space;
239
240    // matrix pd_U
241    pd_u_space_ = orig_pd_u_space;
242
243    DBG_PRINT((1, "Creating the jac_c_space_\n"));
244    // matrix jac_c
245    total_rows = orig_c_space->Dim();
246    total_cols = orig_x_space->Dim() + 2*orig_c_space->Dim()
247                 + 2*orig_d_space->Dim();
248    //orig      + orig_d_l_space->Dim() + orig_d_u_space->Dim();
249    jac_c_space_ = new CompoundMatrixSpace(1, 5, total_rows, total_cols);
250    jac_c_space_->SetBlockRows(0, orig_c_space->Dim());
251    jac_c_space_->SetBlockCols(0, orig_x_space->Dim());
252    jac_c_space_->SetBlockCols(1, orig_c_space->Dim());
253    jac_c_space_->SetBlockCols(2, orig_c_space->Dim());
254    jac_c_space_->SetBlockCols(3, orig_d_space->Dim());
255    jac_c_space_->SetBlockCols(4, orig_d_space->Dim());
256    //orig    jac_c_space_->SetBlockCols(3, orig_d_l_space->Dim());
257    //orig    jac_c_space_->SetBlockCols(4, orig_d_u_space->Dim());
258
259    jac_c_space_->SetCompSpace(0, 0, *orig_jac_c_space);
260    jac_c_space_->SetCompSpace(0, 1, *identity_mat_space_nc, true);
261    jac_c_space_->SetCompSpace(0, 2, *identity_mat_space_nc, true);
262    // remaining blocks are zero'ed
263
264    DBG_PRINT((1, "Creating the jac_d_space_\n"));
265    // matrix jac_d
266    total_rows = orig_d_space->Dim();
267    total_cols = orig_x_space->Dim() + 2*orig_c_space->Dim()
268                 + 2*orig_d_space->Dim();
269    //orig      + orig_d_l_space->Dim() + orig_d_u_space->Dim();
270    jac_d_space_ = new CompoundMatrixSpace(1, 5, total_rows, total_cols);
271    jac_d_space_->SetBlockRows(0, orig_d_space->Dim());
272    jac_d_space_->SetBlockCols(0, orig_x_space->Dim());
273    jac_d_space_->SetBlockCols(1, orig_c_space->Dim());
274    jac_d_space_->SetBlockCols(2, orig_c_space->Dim());
275    jac_d_space_->SetBlockCols(3, orig_d_space->Dim());
276    jac_d_space_->SetBlockCols(4, orig_d_space->Dim());
277    //orig    jac_d_space_->SetBlockCols(3, orig_d_l_space->Dim());
278    //orig    jac_d_space_->SetBlockCols(4, orig_d_u_space->Dim());
279
280    jac_d_space_->SetCompSpace(0, 0, *orig_jac_d_space);
281    // Blocks (0,1) and (0,2) are zero'ed out
282    jac_d_space_->SetCompSpace(0, 3, *identity_mat_space_nd, true);
283    jac_d_space_->SetCompSpace(0, 4, *identity_mat_space_nd, true);
284    //orig    jac_d_space_->SetCompSpace(0, 3, *orig_pd_l_space, true);
285    //orig    SmartPtr<SumMatrixSpace> sum_pd_u
286    //orig    = new SumMatrixSpace(orig_d_space->Dim(), orig_d_u_space->Dim(), 1);
287    //orig    jac_d_space_->SetCompSpace(0, 4, *sum_pd_u, true);
288
289    DBG_PRINT((1, "Creating the h_space_\n"));
290    // matrix h
291    total_dim = orig_x_space->Dim() + 2*orig_c_space->Dim()
292                + 2*orig_d_space->Dim();
293    //orig      + orig_d_l_space->Dim() + orig_d_u_space->Dim();
294    h_space_ = new CompoundSymMatrixSpace(5, total_dim);
295    h_space_->SetBlockDim(0, orig_x_space->Dim());
296    h_space_->SetBlockDim(1, orig_c_space->Dim());
297    h_space_->SetBlockDim(2, orig_c_space->Dim());
298    h_space_->SetBlockDim(3, orig_d_space->Dim());
299    h_space_->SetBlockDim(4, orig_d_space->Dim());
300    //orig    h_space_->SetBlockDim(3, orig_d_l_space->Dim());
301    //orig    h_space_->SetBlockDim(4, orig_d_u_space->Dim());
302
303    SmartPtr<const MatrixSpace> sumsym_mat_space =
304      new SumSymMatrixSpace(orig_x_space->Dim(), 2);
305    h_space_->SetCompSpace(0, 0, *sumsym_mat_space, true);
306    // All remaining blocks are zero'ed out
307
308    SmartPtr<const MatrixSpace> scaled_jac_c_space;
309    SmartPtr<const MatrixSpace> scaled_jac_d_space;
310    SmartPtr<const SymMatrixSpace> scaled_h_space;
311    NLP_scaling()->DetermineScaling(GetRawPtr(x_space_),
312                                    c_space_, d_space_,
313                                    GetRawPtr(jac_c_space_),
314                                    GetRawPtr(jac_d_space_),
315                                    GetRawPtr(h_space_),
316                                    scaled_jac_c_space, scaled_jac_d_space,
317                                    scaled_h_space);
318    // For now we assume that no scaling is done inside the NLP_Scaling
319    DBG_ASSERT(scaled_jac_c_space == jac_c_space_);
320    DBG_ASSERT(scaled_jac_d_space == jac_d_space_);
321    DBG_ASSERT(scaled_h_space == h_space_);
322
323    ///////////////////////////
324    // Create the bound data //
325    ///////////////////////////
326
327    // x_L
328    x_L_ = x_l_space_->MakeNewCompoundVector();
329    x_L_->SetComp(0, *orig_ip_nlp_->x_L()); // x >= x_L
330    x_L_->GetCompNonConst(1)->Set(0.0); // n_c >= 0
331    x_L_->GetCompNonConst(2)->Set(0.0); // p_c >= 0
332    x_L_->GetCompNonConst(3)->Set(0.0); // n_d >= 0
333    x_L_->GetCompNonConst(4)->Set(0.0); // p_d >= 0
334    DBG_PRINT_VECTOR(2,"resto_x_L", *x_L_);
335
336    // x_U
337    x_U_ = orig_ip_nlp_->x_U();
338
339    // d_L
340    d_L_ = orig_ip_nlp_->d_L();
341
342    // d_U
343    d_U_ = orig_ip_nlp_->d_U();
344
345    // Px_L
346    Px_L_ = px_l_space_->MakeNewCompoundMatrix();
347    Px_L_->SetComp(0, 0, *orig_ip_nlp_->Px_L());
348    // Identities are auto-created (true flag passed into SetCompSpace)
349
350    // Px_U
351    Px_U_ = px_u_space_->MakeNewCompoundMatrix();
352    Px_U_->SetComp(0, 0, *orig_ip_nlp_->Px_U());
353    // Remaining matrices will be zero'ed out
354
355    // Pd_L
356    Pd_L_ = orig_ip_nlp_->Pd_L();
357
358    // Pd_U
359    Pd_U_ = orig_ip_nlp_->Pd_U();
360
361    /////////////////////////////////////////////////////////////////////////
362    // Create and initialize the vectors for the restoration phase problem //
363    /////////////////////////////////////////////////////////////////////////
364
365    // Vector x
366    SmartPtr<CompoundVector> comp_x = x_space_->MakeNewCompoundVector();
367    if (init_x) {
368      comp_x->GetCompNonConst(0)->Copy(*orig_ip_data_->curr()->x());
369      comp_x->GetCompNonConst(1)->Set(1.0);
370      comp_x->GetCompNonConst(2)->Set(1.0);
371      comp_x->GetCompNonConst(3)->Set(1.0);
372      comp_x->GetCompNonConst(4)->Set(1.0);
373    }
374    x = GetRawPtr(comp_x);
375
376    // Vector y_c
377    y_c = c_space_->MakeNew();
378    if (init_y_c) {
379      y_c->Set(0.0);  // ToDo
380    }
381
382    // Vector y_d
383    y_d = d_space_->MakeNew();
384    if (init_y_d) {
385      y_d->Set(0.0);
386    }
387
388    // Vector z_L
389    z_L = x_l_space_->MakeNew();
390    if (init_z_L) {
391      z_L->Set(1.0);
392    }
393
394    // Vector z_U
395    z_U = x_u_space_->MakeNew();
396    if (init_z_U) {
397      z_U->Set(1.0);
398    }
399
400    // Vector v_L
401    v_L = d_l_space_->MakeNew();
402
403    // Vector v_U
404    v_U = d_u_space_->MakeNew();
405
406    // Initialize other data needed by the restoration nlp.  x_ref is
407    // the point to reference to which we based the regularization
408    // term
409    x_ref_ = orig_x_space->MakeNew();
410    x_ref_->Copy(*orig_ip_data_->curr()->x());
411
412    SmartPtr<DiagMatrixSpace> DR_x_space
413    = new DiagMatrixSpace(orig_x_space->Dim());
414    dr_x_ = orig_x_space->MakeNew();
415    dr_x_->Set(1.0);
416    SmartPtr<Vector> tmp = dr_x_->MakeNew();
417    tmp->Copy(*x_ref_);
418    dr_x_->ElementWiseMax(*tmp);
419    tmp->Scal(-1.);
420    dr_x_->ElementWiseMax(*tmp);
421    dr_x_->ElementWiseReciprocal();
422    DBG_PRINT_VECTOR(2, "dr_x_", *dr_x_);
423    DR_x_ = DR_x_space->MakeNewDiagMatrix();
424    DR_x_->SetDiag(*dr_x_);
425
426    return true;
427  }
428
429  Number RestoIpoptNLP::f(const Vector& x, Number mu)
430  {
431    DBG_START_METH("RestoIpoptNLP::f",
432                   dbg_verbosity);
433    Number ret = 0.0;
434    // rho*(pcTe + ncTe + pdT*e + ndT*e) + eta/2*||Dr*(x-xr)||_2^2
435    const CompoundVector* c_vec = dynamic_cast<const CompoundVector*>(&x);
436    DBG_ASSERT(c_vec);
437    SmartPtr<const Vector> x_only = c_vec->GetComp(0);
438    ret = x.Sum() - x_only->Sum();
439    DBG_PRINT((1,"xdiff sum = %e\n",ret));
440    ret = rho_ * ret;
441    DBG_PRINT((1,"rho_ = %e\n",rho_));
442
443    SmartPtr<Vector> x_diff = x_only->MakeNew();
444    x_diff->Copy(*x_only);
445    x_diff->Axpy(-1.0, *x_ref_);
446    DBG_PRINT_VECTOR(2,"x_ref",*x_ref_);
447    x_diff->ElementWiseMultiply(*dr_x_);
448    Number ret2 = x_diff->Nrm2();
449    DBG_PRINT((1,"Eta = %e\n",Eta(mu)));
450    ret2 = Eta(mu)/2.0*ret2*ret2;
451
452    ret += ret2;
453
454    // We evaluate also the objective function for the original
455    // problem here.  This might be wasteful, but it will detect if
456    // the original objective function cannot be evaluated at the
457    // trial point in the restoration phase
458    if (evaluate_orig_obj_at_resto_trial_) {
459      Number orig_f = orig_ip_nlp_->f(*x_only);
460    }
461
462    return ret;
463  }
464
465  SmartPtr<const Vector> RestoIpoptNLP::grad_f(const Vector& x, Number mu)
466  {
467    SmartPtr<Vector> retPtr = x.MakeNew();
468    // Scale the p's and n's by rho (Scale all, take out the x part later)
469    retPtr->Set(rho_);
470
471    const CompoundVector* c_vec_in = dynamic_cast<const CompoundVector*>(&x);
472    SmartPtr<const Vector> x_only_in = c_vec_in->GetComp(0);
473
474    CompoundVector* c_vec = dynamic_cast<CompoundVector*>(GetRawPtr(retPtr));
475    DBG_ASSERT(c_vec);
476    SmartPtr<Vector> x_only = c_vec->GetCompNonConst(0);
477    x_only->Copy(*x_only_in);
478    x_only->Axpy(-1.0, *x_ref_);
479    x_only->ElementWiseMultiply(*dr_x_);
480    x_only->Scal(Eta(mu));
481
482    return ConstPtr(retPtr);
483  }
484
485  SmartPtr<const Vector> RestoIpoptNLP::c(const Vector& x)
486  {
487    const CompoundVector* c_vec = dynamic_cast<const CompoundVector*>(&x);
488    SmartPtr<const Vector> x_only = c_vec->GetComp(0);
489    SmartPtr<const Vector> nc_only = c_vec->GetComp(1);
490    SmartPtr<const Vector> pc_only = c_vec->GetComp(2);
491
492    SmartPtr<const Vector> orig_c = orig_ip_nlp_->c(*x_only);
493    SmartPtr<Vector> retPtr = c_space_->MakeNew();
494    retPtr->Copy(*orig_c);
495    retPtr->Axpy(1.0, *nc_only);
496    retPtr->Axpy(-1.0, *pc_only);
497
498    return GetRawPtr(retPtr);
499  }
500
501
502  SmartPtr<const Vector> RestoIpoptNLP::d(const Vector& x)
503  {
504    const CompoundVector* c_vec = dynamic_cast<const CompoundVector*>(&x);
505    SmartPtr<const Vector> x_only = c_vec->GetComp(0);
506    SmartPtr<const Vector> nd_only = c_vec->GetComp(3);
507    SmartPtr<const Vector> pd_only = c_vec->GetComp(4);
508
509    SmartPtr<const Vector> orig_d = orig_ip_nlp_->d(*x_only);
510    SmartPtr<Vector> retPtr = d_space_->MakeNew();
511    retPtr->Copy(*orig_d);
512    retPtr->Axpy(1., *nd_only);
513    retPtr->Axpy(-1., *pd_only);
514#ifdef orig
515
516    SmartPtr<Vector> tmp = orig_d->MakeNew();
517    orig_ip_nlp_->Pd_L()->MultVector(1.0, *nd_only, 0.0, *tmp);
518    retPtr->Axpy(1.0, *tmp);
519    orig_ip_nlp_->Pd_U()->MultVector(1.0, *pd_only, 0.0, *tmp);
520    retPtr->Axpy(-1.0, *tmp);
521#endif
522
523    return GetRawPtr(retPtr);
524  }
525
526  SmartPtr<const Matrix> RestoIpoptNLP::jac_c(const Vector& x)
527  {
528
529    // Here, we set the (0,0) block with the values from the
530    // original jac_c and set the factor for the -I (jac w.r.t. p_c)
531
532    // get out the x_only part
533    const CompoundVector* c_vec = dynamic_cast<const CompoundVector*>(&x);
534    DBG_ASSERT(c_vec);
535    SmartPtr<const Vector> x_only = c_vec->GetComp(0);
536
537    // calculate the jacobian for the original problem
538    SmartPtr<const Matrix> jac_c_only = orig_ip_nlp_->jac_c(*x_only);
539
540    // Create the new compound matrix
541    // The zero parts remain NULL, the identities are created from the matrix
542    // space (since auto_allocate was set to true in SetCompSpace)
543    SmartPtr<CompoundMatrix> retPtr = jac_c_space_->MakeNewCompoundMatrix();
544
545    // set the (0,0) block to the original jacobian
546    retPtr->SetComp(0,0,*jac_c_only);
547
548    // we currently do not have a default factor in the matrix spaces
549    // so we need to set the factor on the identity (jacobian of the
550    // restoration c w.r.t. p_c is -I)
551    // This could easily be changed to include special processing
552    // for identities in the CompoundMatrixSpace (and a factor)
553    SmartPtr<Matrix> jac_c_pc_mat = retPtr->GetCompNonConst(0,2);
554    IdentityMatrix* jac_c_pc = dynamic_cast<IdentityMatrix*>(GetRawPtr(jac_c_pc_mat));
555    DBG_ASSERT(jac_c_pc);
556    jac_c_pc->SetFactor(-1.0);
557
558    return GetRawPtr(retPtr);
559  }
560
561  SmartPtr<const Matrix> RestoIpoptNLP::jac_d(const Vector& x)
562  {
563    // Here, we set the (0,0) block with the values from the
564    // original jac_d and set the factor for the -I (jac w.r.t. p_d)
565
566    // get out the x_only part
567    const CompoundVector* c_vec = dynamic_cast<const CompoundVector*>(&x);
568    DBG_ASSERT(c_vec);
569    SmartPtr<const Vector> x_only = c_vec->GetComp(0);
570
571    // calculate the jacobian for the original problem
572    SmartPtr<const Matrix> jac_d_only = orig_ip_nlp_->jac_d(*x_only);
573
574    // Create the new compound matrix
575    // The zero parts remain NULL, the identities are created from the matrix
576    // space (since auto_allocate was set to true in SetCompSpace)
577    SmartPtr<CompoundMatrix> retPtr = jac_d_space_->MakeNewCompoundMatrix();
578
579    // Set the block for the original jacobian
580    retPtr->SetComp(0,0,*jac_d_only);
581
582    // (0,1) and (0,2) blocks are zero (NULL)
583
584    // set the factor for the identity matrix for the pd variables
585    // (likr in jac_c)
586    SmartPtr<Matrix> jac_d_pd_mat = retPtr->GetCompNonConst(0,4);
587    IdentityMatrix* jac_d_pd = dynamic_cast<IdentityMatrix*>(GetRawPtr(jac_d_pd_mat));
588    DBG_ASSERT(jac_d_pd);
589    jac_d_pd->SetFactor(-1.0);
590
591#ifdef orig
592    // Jacobian of resto d w.r.t. n_d is Pd_L
593    retPtr->SetComp(0, 3, *orig_ip_nlp_->Pd_L());
594
595    // Change the matrix factor to -1 for Pd_U
596    // Jacobian of the resto d w.r.t. p_d is -Pd_U
597    SmartPtr<Matrix> jac_d_pd_mat = retPtr->GetCompNonConst(0,4);
598    SumMatrix* jac_d_pd_sum = dynamic_cast<SumMatrix*>(GetRawPtr(jac_d_pd_mat));
599    DBG_ASSERT(jac_d_pd_sum);
600    jac_d_pd_sum->SetTerm(0, -1.0, *orig_ip_nlp_->Pd_U());
601#endif
602
603    return GetRawPtr(retPtr);
604  }
605
606  SmartPtr<const SymMatrix> RestoIpoptNLP::h(const Vector& x,
607      Number obj_factor,
608      const Vector& yc,
609      const Vector& yd,
610      Number mu)
611  {
612    // Here, we use a SumSymMatrix for the (0,0) block of the
613    // Hessian. We need to set this to the hessian of the restoration
614    // problem, which is the hessian of the objective from the restoration
615    // problem + the constraint only part of the hessian from the original
616    // problem
617    // All other blocks are zero'ed (NULL)
618
619    // get the x_only part
620    const CompoundVector* c_vec = dynamic_cast<const CompoundVector*>(&x);
621    DBG_ASSERT(c_vec);
622    SmartPtr<const Vector> x_only = c_vec->GetComp(0);
623
624    // yc and yd should not be compound vectors
625
626    // calculate the original hessian
627    SmartPtr<const SymMatrix> h_con_orig = orig_ip_nlp_->h(*x_only, 0.0, yc, yd);
628
629    // Create the new compound matrix
630    // The SumSymMatrix is auto_allocated
631    SmartPtr<CompoundSymMatrix> retPtr = h_space_->MakeNewCompoundSymMatrix();
632
633    // Set the entries in the SumSymMatrix
634    SmartPtr<Matrix> h_sum_mat = retPtr->GetCompNonConst(0,0);
635    SmartPtr<SumSymMatrix> h_sum = dynamic_cast<SumSymMatrix*>(GetRawPtr(h_sum_mat));
636    h_sum->SetTerm(0, 1.0, *h_con_orig);
637    h_sum->SetTerm(1, obj_factor*Eta(mu), *DR_x_);
638
639    return GetRawPtr(retPtr);
640  }
641
642  void RestoIpoptNLP::GetSpaces(SmartPtr<const VectorSpace>& x_space,
643                                SmartPtr<const VectorSpace>& c_space,
644                                SmartPtr<const VectorSpace>& d_space,
645                                SmartPtr<const VectorSpace>& x_l_space,
646                                SmartPtr<const MatrixSpace>& px_l_space,
647                                SmartPtr<const VectorSpace>& x_u_space,
648                                SmartPtr<const MatrixSpace>& px_u_space,
649                                SmartPtr<const VectorSpace>& d_l_space,
650                                SmartPtr<const MatrixSpace>& pd_l_space,
651                                SmartPtr<const VectorSpace>& d_u_space,
652                                SmartPtr<const MatrixSpace>& pd_u_space,
653                                SmartPtr<const MatrixSpace>& Jac_c_space,
654                                SmartPtr<const MatrixSpace>& Jac_d_space,
655                                SmartPtr<const SymMatrixSpace>& Hess_lagrangian_space)
656  {
657    x_space = GetRawPtr(x_space_);
658    c_space = GetRawPtr(c_space_);
659    d_space = GetRawPtr(d_space_);
660    x_l_space = GetRawPtr(x_l_space_);
661    px_l_space = GetRawPtr(px_l_space_);
662    x_u_space = GetRawPtr(x_u_space_);
663    px_u_space = GetRawPtr(px_u_space_);
664    d_l_space = GetRawPtr(d_l_space_);
665    pd_l_space = GetRawPtr(pd_l_space_);
666    d_u_space = GetRawPtr(d_u_space_);
667    pd_u_space = GetRawPtr(pd_u_space_);
668    Jac_c_space = GetRawPtr(jac_c_space_);
669    Jac_d_space = GetRawPtr(jac_d_space_);
670    Hess_lagrangian_space = GetRawPtr(h_space_);
671  }
672
673  Number RestoIpoptNLP::Eta(Number mu) const
674  {
675    return eta_factor_ * pow(mu, eta_mu_exponent_);
676  }
677
678  void RestoIpoptNLP::AdjustVariableBounds(const Vector& new_x_L, const Vector& new_x_U,
679      const Vector& new_d_L, const Vector& new_d_U)
680  {
681
682    const CompoundVector* comp_new_x_L =
683      dynamic_cast<const CompoundVector*>(&new_x_L);
684    DBG_ASSERT(comp_new_x_L);
685
686    SmartPtr<const Vector> new_orig_x_L = comp_new_x_L->GetComp(0);
687
688    // adapt bounds for the original NLP
689    orig_ip_nlp_->AdjustVariableBounds(*new_orig_x_L, new_x_U, new_d_L, new_d_U);
690
691    // adapt bounds for the p and n variables
692    SmartPtr<const Vector> new_nc_L = comp_new_x_L->GetComp(1);
693    SmartPtr<const Vector> new_pc_L = comp_new_x_L->GetComp(2);
694    SmartPtr<const Vector> new_nd_L = comp_new_x_L->GetComp(3);
695    SmartPtr<const Vector> new_pd_L = comp_new_x_L->GetComp(4);
696
697    x_L_->GetCompNonConst(1)->Copy(*new_nc_L);
698    x_L_->GetCompNonConst(2)->Copy(*new_pc_L);
699    x_L_->GetCompNonConst(3)->Copy(*new_nd_L);
700    x_L_->GetCompNonConst(4)->Copy(*new_pd_L);
701
702  }
703
704} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.