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