CppAD: A C++ Algorithmic Differentiation Package
20130918
|
00001 /* $Id$ */ 00002 /* -------------------------------------------------------------------------- 00003 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 Bradley M. Bell 00004 00005 CppAD is distributed under multiple licenses. This distribution is under 00006 the terms of the 00007 Eclipse Public License Version 1.0. 00008 00009 A copy of this license is included in the COPYING file of this distribution. 00010 Please visit http://www.coin-or.org/CppAD/ for information on other licenses. 00011 -------------------------------------------------------------------------- */ 00012 # include <limits> 00013 00014 # include "cppad_ipopt_nlp.hpp" 00015 # include "sparse_map2vec.hpp" 00016 # include "jac_g_map.hpp" 00017 # include "hes_fg_map.hpp" 00018 # include "vec_fun_pattern.hpp" 00019 # include "fun_record.hpp" 00020 00021 00022 /// If 0 tracing is off, otherwise tracing is on. 00023 # define CPPAD_IPOPT_NLP_TRACE 0 00024 00025 # if CPPAD_IPOPT_NLP_TRACE 00026 # include <cstdio> 00027 # endif 00028 00029 // --------------------------------------------------------------------------- 00030 namespace cppad_ipopt { 00031 // --------------------------------------------------------------------------- 00032 /*! 00033 \{ 00034 \file cppad_ipopt_nlp.cpp 00035 \brief Member functions for the cppad_ipopt_nlp class. 00036 */ 00037 00038 00039 /*! 00040 Constructor for the \ref Nonlinear_Programming_Problem. 00041 00042 \param n 00043 dimension of the domain space for f(x) and g(x). 00044 00045 \param m 00046 dimension of the range space for g(x) 00047 00048 \param x_i 00049 initial value of x during the optimization procedure (size n). 00050 00051 \param x_l 00052 lower limit for x (size n). 00053 00054 \param x_u 00055 upper limit for x (size n). 00056 00057 \param g_l 00058 lower limit for g(x) (size m). 00059 00060 \param g_u 00061 upper limit for g(x) (size m). 00062 00063 \param fg_info 00064 pointer to base class version of derived class object used to get 00065 information about the user's representation for f(x) and g(x). 00066 (The object pointed to must not be deleted before this cppad_ipopt_nlp object). 00067 00068 \param solution 00069 pointer to object where final results are stored. 00070 (The object pointed to must not be deleted before this cppad_ipopt_nlp object). 00071 00072 \par Constants 00073 The following values are set by the constructor and are \c const 00074 or effectively \c const; i.e., they are set by the constructor and should 00075 not be changed: 00076 \verbatim 00077 n_, m_, x_i_, x_l_, x_u_, g_l_, g_u_, K_, L_, p_, q_, retape_, 00078 pattern_jac_r_, pattern_hes_r_, index_jac_g_, index_hes_fg_, 00079 nnz_jac_g_, iRow_jac_g_, jCol_jac_g_, 00080 nnz_h_lag_, iRow_h_lag_, jCol_h_lag_, 00081 \endverbatim 00082 In addition, the function calls <tt>fg_info->set_n(n)</tt> 00083 and <tt>fg_info->set_m(m)</tt> are used to set the values of \c n 00084 and \c m in \c fg_info. 00085 00086 \par Variables 00087 The following arrays have fixed size which is set during this constructor: 00088 00089 \li \c tape_ok_ has size \c K_. It is initialized as true for indices 00090 \c k such that <tt>retape[k]</tt> is false. 00091 00092 \li \c r_fun_ has size \c K_. It is initilaize with the default 00093 \c ADFun constructor. Then, for indices \c k such that 00094 <tt>retape[k]</tt> is false, the operation sequence corresponding 00095 to \f$ r_k (u) \f$ is stored in <tt>r_fun_[k]</tt>. 00096 00097 \li \c I_ has size equal to the maximum of <tt>p[k]</tt> w.r.t \c k. 00098 00099 \li \c J_ has size equal to the maximum of <tt>q[k]</tt> w.r.t \c k. 00100 00101 \par NDEBUG 00102 If the preprocessor symbol \c NEBUG is not defined, 00103 certain of the assumptions about the function calls of the form 00104 \verbatim 00105 fg_info->index(k, ell, I, J) 00106 \endverbatim 00107 are checked to make sure they hold. 00108 */ 00109 cppad_ipopt_nlp::cppad_ipopt_nlp( 00110 size_t n , 00111 size_t m , 00112 const NumberVector &x_i , 00113 const NumberVector &x_l , 00114 const NumberVector &x_u , 00115 const NumberVector &g_l , 00116 const NumberVector &g_u , 00117 cppad_ipopt_fg_info* fg_info , 00118 cppad_ipopt_solution* solution ) 00119 : n_ ( n ), 00120 m_ ( m ), 00121 x_i_ ( x_i ), 00122 x_l_ ( x_l ), 00123 x_u_ ( x_u ), 00124 g_l_ ( g_l ), 00125 g_u_ ( g_u ), 00126 fg_info_ ( fg_info ) , 00127 solution_ (solution) , 00128 infinity_ ( std::numeric_limits<Number>::infinity() ) 00129 { size_t k; 00130 00131 // set information needed in cppad_ipopt_fg_info 00132 fg_info_->set_n(n); 00133 fg_info_->set_m(m); 00134 00135 // get information from derived class version of fg_info 00136 K_ = fg_info_->number_functions(); 00137 L_.resize(K_); 00138 p_.resize(K_); 00139 q_.resize(K_); 00140 r_fun_.resize(K_); 00141 retape_.resize(K_); 00142 tape_ok_.resize(K_); 00143 pattern_jac_r_.resize(K_); 00144 pattern_hes_r_.resize(K_); 00145 size_t max_p = 0; 00146 size_t max_q = 0; 00147 for(k = 0; k < K_; k++) 00148 { L_[k] = fg_info_->number_terms(k); 00149 p_[k] = fg_info_->range_size(k); 00150 q_[k] = fg_info_->domain_size(k); 00151 retape_[k] = fg_info_->retape(k); 00152 max_p = std::max(max_p, p_[k]); 00153 max_q = std::max(max_q, q_[k]); 00154 pattern_jac_r_[k].resize( p_[k] * q_[k] ); 00155 pattern_hes_r_[k].resize( q_[k] * q_[k] ); 00156 } 00157 I_.resize(max_p); 00158 J_.resize(max_q); 00159 # ifndef NDEBUG 00160 size_t i, j, ell; 00161 // check for valid range and domain indices 00162 for(k = 0; k < K_; k++) for(ell = 0; ell < L_[k]; ell++) 00163 { 00164 for( i = 0; i < p_[k]; i++) 00165 I_[i] = m+1; // an invalid range index 00166 for( j = 0; j < q_[k]; j++) 00167 J_[j] = n; // an invalid domain index 00168 fg_info_->index(k, ell, I_, J_); 00169 for( i = 0; i < p_[k]; i++) if( I_[i] > m ) 00170 { std::cerr << "k=" << k << ", ell=" << ell 00171 << ", I[" << i << "]=" << I_[i] << std::endl; 00172 CPPAD_ASSERT_KNOWN( I_[i] <= m, 00173 "cppad_ipopt_nlp: invalid value in index vector I" 00174 ); 00175 } 00176 for( j = 0; j < q_[k]; j++) if( J_[j] >= n ) 00177 { std::cerr << "k=" << k << ", ell=" << ell 00178 << ", J[" << j << "]=" << J_[j] << std::endl; 00179 CPPAD_ASSERT_KNOWN( J_[j] < n, 00180 "cppad_ipopt_nlp: invalid value in index vector J" 00181 ); 00182 } 00183 } 00184 # endif 00185 // record r[k] for functions that do not need retaping 00186 for(k = 0; k < K_; k++) 00187 { tape_ok_[k] = false; 00188 if( ! retape_[k] ) 00189 { // Operation sequence does not depend on value 00190 // of u so record it once here in the constructor. 00191 fg_info_->index(k, 0, I_, J_); 00192 fun_record( 00193 fg_info_ , // inputs 00194 k , 00195 p_ , 00196 q_ , 00197 n_ , 00198 x_i_ , 00199 J_ , 00200 r_fun_ // output 00201 ); 00202 // take time to optimize because only recording once 00203 r_fun_[k].optimize(); 00204 // ok and will stay that way 00205 tape_ok_[k] = true; 00206 } 00207 } 00208 00209 // compute a sparsity patterns for each r_k (u) 00210 vec_fun_pattern( 00211 K_, p_, q_, retape_, r_fun_, // inputs 00212 pattern_jac_r_, pattern_hes_r_ // outputs 00213 ); 00214 00215 // mapping from (i,j) to Ipopt sparsity index for Jacobian of g 00216 jac_g_map( 00217 fg_info_, m_, n_, K_, L_, p_, q_, pattern_jac_r_, // inputs 00218 I_, J_, // work 00219 index_jac_g_ // outputs 00220 ); 00221 00222 // mapping from (i,j) to Ipopt sparsity index for Hessian of Lagragian 00223 hes_fg_map( 00224 fg_info_, m_, n_, K_, L_, p_, q_, pattern_hes_r_, // inputs 00225 I_, J_, // work 00226 index_hes_fg_ // outputs 00227 ); 00228 00229 // Compute Ipopt sparsity structure for Jacobian of g 00230 sparse_map2vec( 00231 index_jac_g_, // inputs 00232 nnz_jac_g_, iRow_jac_g_, jCol_jac_g_ // outputs 00233 ); 00234 00235 // Compute Ipopt sparsity structure for Hessian of Lagragian 00236 sparse_map2vec( 00237 index_hes_fg_, // inputs 00238 nnz_h_lag_, iRow_h_lag_, jCol_h_lag_ // outputs 00239 ); 00240 } 00241 00242 /// The destructor takes no special action. 00243 cppad_ipopt_nlp::~cppad_ipopt_nlp() 00244 {} 00245 00246 /*! 00247 Return dimension information about optimization problem. 00248 00249 \param[out] n 00250 is set to the value \c n_. 00251 00252 \param[out] m 00253 is set to the value \c m_. 00254 00255 \param[out] nnz_jac_g 00256 is set to the value of \c nnz_jac_g_. 00257 00258 \param[out] nnz_h_lag 00259 is set to the vlaue of \c nnz_h_lag_. 00260 00261 \param[out] index_style 00262 is set to C_STYLE; i.e., zeoro based indexing is used in the 00263 information passed to Ipopt. 00264 */ 00265 bool cppad_ipopt_nlp::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g, 00266 Index& nnz_h_lag, IndexStyleEnum& index_style) 00267 { 00268 n = n_; 00269 m = m_; 00270 nnz_jac_g = nnz_jac_g_; 00271 nnz_h_lag = nnz_h_lag_; 00272 00273 // use the fortran index style for row/col entries 00274 index_style = C_STYLE; 00275 00276 return true; 00277 } 00278 00279 /*! 00280 Return bound information about optimization problem. 00281 00282 \param[in] n 00283 is the dimension of the domain space for f(x) and g(x); i.e., 00284 it must be equal to \c n_. 00285 00286 \param[out] x_l 00287 is a vector of size \c n. 00288 The input value of its elements does not matter. 00289 On output, it is a copy of the lower bound for \f$ x \f$; i.e., 00290 \c x_l_. 00291 00292 \param[out] x_u 00293 is a vector of size \c n. 00294 The input value of its elements does not matter. 00295 On output, it is a copy of the upper bound for \f$ x \f$; i.e., 00296 \c x_u_. 00297 00298 \param[in] m 00299 is the dimension of the range space for g(x). i.e., 00300 it must be equal to \c m_. 00301 00302 \param[out] g_l 00303 is a vector of size \c m. 00304 The input value of its elements does not matter. 00305 On output, it is a copy of the lower bound for \f$ g(x) \f$; i.e., \c g_l_. 00306 00307 \param[out] g_u 00308 is a vector of size \c m. 00309 The input value of its elements does not matter. 00310 On output, it is a copy of the upper bound for \f$ g(x) \f$; i.e, \c g_u_. 00311 */ 00312 bool cppad_ipopt_nlp::get_bounds_info(Index n, Number* x_l, Number* x_u, 00313 Index m, Number* g_l, Number* g_u) 00314 { size_t i, j; 00315 // here, the n and m we gave IPOPT in get_nlp_info are passed back 00316 CPPAD_ASSERT_UNKNOWN(size_t(n) == n_); 00317 CPPAD_ASSERT_UNKNOWN(size_t(m) == m_); 00318 00319 // pass back bounds 00320 for(j = 0; j < n_; j++) 00321 { x_l[j] = x_l_[j]; 00322 x_u[j] = x_u_[j]; 00323 } 00324 for(i = 0; i < m_; i++) 00325 { g_l[i] = g_l_[i]; 00326 g_u[i] = g_u_[i]; 00327 } 00328 00329 return true; 00330 } 00331 00332 /*! 00333 Return initial x value where optimiation is started. 00334 00335 \param[in] n 00336 must be equal to the domain dimension for f(x) and g(x); i.e., 00337 it must be equal to \c n_. 00338 00339 \param[in] init_x 00340 must be equal to true. 00341 00342 \param[out] x 00343 is a vector of size \c n. 00344 The input value of its elements does not matter. 00345 On output, it is a copy of the initial value for \f$ x \f$; i.e. \c x_i_. 00346 00347 \param[in] init_z 00348 must be equal to false. 00349 00350 \param z_L 00351 is not used. 00352 00353 \param z_U 00354 is not used. 00355 00356 \param[in] m 00357 must be equal to the range dimension for g(x); i.e., 00358 it must be equal to \c m_. 00359 00360 \param init_lambda 00361 must be equal to false. 00362 00363 \param lambda 00364 is not used. 00365 */ 00366 bool cppad_ipopt_nlp::get_starting_point(Index n, bool init_x, Number* x, 00367 bool init_z, Number* z_L, Number* z_U, 00368 Index m, bool init_lambda, 00369 Number* lambda) 00370 { size_t j; 00371 00372 CPPAD_ASSERT_UNKNOWN(size_t(n) == n_ ); 00373 CPPAD_ASSERT_UNKNOWN(size_t(m) == m_ ); 00374 CPPAD_ASSERT_UNKNOWN(init_x == true); 00375 CPPAD_ASSERT_UNKNOWN(init_z == false); 00376 CPPAD_ASSERT_UNKNOWN(init_lambda == false); 00377 00378 for(j = 0; j < n_; j++) 00379 x[j] = x_i_[j]; 00380 00381 return true; 00382 } 00383 00384 /*! 00385 Evaluate the objective fucntion f(x). 00386 00387 \param[in] n 00388 is the dimension of the argument space for f(x); i.e., must be equal \c n_. 00389 00390 \param[in] x 00391 is a vector of size \c n containing the point at which to evaluate 00392 the function f(x). 00393 00394 \param[in] new_x 00395 is false if the previous call to any one of the 00396 \ref Evaluation_Methods used the same value for \c x. 00397 00398 \param[out] obj_value 00399 is the value of the objective f(x) at this value of \c x. 00400 00401 \return 00402 The return value is always true; see \ref Evaluation_Methods. 00403 00404 \par Efficiency 00405 This routine could be more efficient 00406 (for certain when when L[k] > 1 and retape[k] is true) 00407 if the users also provided a version 00408 of the function <tt>fg_info->eval_r(k, u)</tt> where \c u was of type 00409 \c NumberVector. 00410 */ 00411 bool cppad_ipopt_nlp::eval_f( 00412 Index n, const Number* x, bool new_x, Number& obj_value 00413 ) 00414 { 00415 CPPAD_ASSERT_UNKNOWN(size_t(n) == n_ ); 00416 00417 size_t iobj, j, k, ell; 00418 00419 // initialize summation 00420 obj_value = 0.; 00421 00422 // update tape_ok_ flag 00423 for(k = 0; k < K_; k++) 00424 { if( retape_[k] && (new_x || L_[k] > 1) ) 00425 tape_ok_[k] = false; 00426 } 00427 00428 for(k = 0; k < K_; k++) for(ell = 0; ell < L_[k]; ell++) 00429 { fg_info_->index(k, ell, I_, J_); 00430 for(iobj = 0; iobj < p_[k]; iobj++) if( I_[iobj] == 0 ) 00431 { if( ! tape_ok_[k] ) 00432 { // Record r_k for value of u corresponding to x 00433 fun_record( 00434 fg_info_ , // inputs 00435 k , 00436 p_ , 00437 q_ , 00438 n_ , 00439 x , 00440 J_ , 00441 r_fun_ // output 00442 ); 00443 tape_ok_[k] = ! (retape_[k] || L_[k] > 1); 00444 } 00445 NumberVector u(q_[k]); 00446 NumberVector r(p_[k]); 00447 for(j = 0; j < q_[k]; j++) 00448 { CPPAD_ASSERT_UNKNOWN( J_[j] < n_ ); 00449 u[j] = x[ J_[j] ]; 00450 } 00451 r = r_fun_[k].Forward(0, u); 00452 obj_value += r[iobj]; 00453 } 00454 } 00455 # if CPPAD_IPOPT_NLP_TRACE 00456 using std::printf; 00457 for(j = 0; j < n_; j++) 00458 printf("cppad_ipopt_nlp::eval_f::x[%d] = %20.14g\n", j, x[j]); 00459 printf("cppad_ipopt_nlp::eval_f::obj_value = %20.14g\n", obj_value); 00460 # endif 00461 # ifndef NDEBUG 00462 CPPAD_ASSERT_KNOWN( 00463 (-infinity_ < obj_value) && (obj_value < infinity_), 00464 "cppad_ipopt_nlp::eval_f:: objective value is not finite" 00465 ); 00466 # endif 00467 return true; 00468 } 00469 00470 /*! 00471 Evaluate the gradient of f(x). 00472 00473 \param[in] n 00474 is the dimension of the argument space for f(x); i.e., must be equal \c n_. 00475 00476 \param[in] x 00477 has a vector of size \c n containing the point at which to evaluate 00478 the gradient of f(x). 00479 00480 \param[in] new_x 00481 is false if the previous call to any one of the 00482 \ref Evaluation_Methods used the same value for \c x. 00483 00484 \param[out] grad_f 00485 is a vector of size \c n. 00486 The input value of its elements does not matter. 00487 The output value of its elements is the gradient of f(x) 00488 at this value of. 00489 00490 \return 00491 The return value is always true; see \ref Evaluation_Methods. 00492 */ 00493 bool cppad_ipopt_nlp::eval_grad_f( 00494 Index n, const Number* x, bool new_x, Number* grad_f 00495 ) 00496 { CPPAD_ASSERT_UNKNOWN(size_t(n) == n_ ); 00497 00498 size_t iobj, i, j, k, ell; 00499 00500 // initialize summation 00501 for(j = 0; j < n_; j++) 00502 grad_f[j] = 0.; 00503 00504 // update tape_ok_ flag 00505 for(k = 0; k < K_; k++) 00506 { if( retape_[k] && (new_x || L_[k] > 1) ) 00507 tape_ok_[k] = false; 00508 } 00509 00510 for(k = 0; k < K_; k++) for(ell = 0; ell < L_[k]; ell++) 00511 { fg_info_->index(k, ell, I_, J_); 00512 for(iobj = 0; iobj < p_[k]; iobj++) if( I_[iobj] == 0 ) 00513 { if( ! tape_ok_[k] ) 00514 { // Record r_k for value of u corresponding to x 00515 fun_record( 00516 fg_info_ , // inputs 00517 k , 00518 p_ , 00519 q_ , 00520 n_ , 00521 x , 00522 J_ , 00523 r_fun_ // output 00524 ); 00525 tape_ok_[k] = ! (retape_[k] || L_[k] > 1); 00526 } 00527 NumberVector u(q_[k]); 00528 NumberVector w(p_[k]); 00529 NumberVector r_grad(q_[k]); 00530 for(j = 0; j < q_[k]; j++) 00531 { CPPAD_ASSERT_UNKNOWN( J_[j] < n_ ); 00532 u[j] = x[ J_[j] ]; 00533 } 00534 r_fun_[k].Forward(0, u); 00535 for(i = 0; i < p_[k]; i++) 00536 w[i] = 0.; 00537 w[iobj] = 1.; 00538 r_grad = r_fun_[k].Reverse(1, w); 00539 for(j = 0; j < q_[k]; j++) 00540 { CPPAD_ASSERT_UNKNOWN( J_[j] < n_ ); 00541 grad_f[ J_[j] ] += r_grad[j]; 00542 } 00543 } 00544 } 00545 # if CPPAD_IPOPT_NLP_TRACE 00546 using std::printf; 00547 for(j = 0; j < n_; j++) printf( 00548 "cppad_ipopt_nlp::eval_grad_f::x[%d] = %20.14g\n", j, x[j] 00549 ); 00550 for(j = 0; j < n_; j++) printf( 00551 "cppad_ipopt_nlp::eval_grad_f::grad_f[%d] = %20.14g\n", j, grad_f[j] 00552 ); 00553 # endif 00554 # ifndef NDEBUG 00555 for(j = 0; j < n_; j++) CPPAD_ASSERT_KNOWN( 00556 (-infinity_ < grad_f[j]) && (grad_f[j] < infinity_), 00557 "cppad_ipopt_nlp::grad_f:: gradient of objective is not finite" 00558 ); 00559 # endif 00560 return true; 00561 } 00562 00563 /*! 00564 Evaluate the function g(x). 00565 00566 \param[in] n 00567 is the dimension of the argument space for g(x); i.e., must be equal \c n_. 00568 00569 \param[in] x 00570 has a vector of size \c n containing the point at which to evaluate 00571 the constraint function g(x). 00572 00573 \param[in] new_x 00574 is false if the previous call to any one of the 00575 \ref Evaluation_Methods used the same value for \c x. 00576 00577 \param[in] m 00578 is the dimension of the range space for g(x); i.e., must be equal to \c m_. 00579 00580 \param[out] g 00581 is a vector of size \c m. 00582 The input value of its elements does not matter. 00583 The output value of its elements is 00584 the value of the function g(x) at this value of \c x. 00585 00586 \return 00587 The return value is always true; see \ref Evaluation_Methods. 00588 */ 00589 bool cppad_ipopt_nlp::eval_g( 00590 Index n, const Number* x, bool new_x, Index m, Number* g 00591 ) 00592 { CPPAD_ASSERT_UNKNOWN(size_t(n) == n_ ); 00593 00594 size_t i, j, k, ell; 00595 00596 // initialize summation 00597 for(i = 0; i < m_; i++) 00598 g[i] = 0.; 00599 00600 // update tape_ok_ flag 00601 for(k = 0; k < K_; k++) 00602 { if( retape_[k] && (new_x || L_[k] > 1) ) 00603 tape_ok_[k] = false; 00604 } 00605 00606 for(k = 0; k < K_; k++) for(ell = 0; ell < L_[k]; ell++) 00607 { fg_info_->index(k, ell, I_, J_); 00608 if( ! tape_ok_[k] ) 00609 { // Record r_k for value of u corresponding to x 00610 fun_record( 00611 fg_info_ , // inputs 00612 k , 00613 p_ , 00614 q_ , 00615 n_ , 00616 x , 00617 J_ , 00618 r_fun_ // output 00619 ); 00620 } 00621 tape_ok_[k] = ! (retape_[k] || L_[k] > 1); 00622 NumberVector u(q_[k]); 00623 NumberVector r(p_[k]); 00624 for(j = 0; j < q_[k]; j++) 00625 { CPPAD_ASSERT_UNKNOWN( J_[j] < n_ ); 00626 u[j] = x[ J_[j] ]; 00627 } 00628 r = r_fun_[k].Forward(0, u); 00629 for(i = 0; i < p_[k]; i++) 00630 { CPPAD_ASSERT_UNKNOWN( I_[i] <= m_ ); 00631 if( I_[i] >= 1 ) 00632 g[ I_[i] - 1 ] += r[i]; 00633 } 00634 } 00635 # if CPPAD_IPOPT_NLP_TRACE 00636 using std::printf; 00637 for(j = 0; j < n_; j++) 00638 printf("cppad_ipopt_nlp::eval_g::x[%d] = %20.14g\n", j, x[j]); 00639 for(i = 0; i < m_; i++) 00640 printf("cppad_ipopt_nlp::eval_g::g[%d] = %20.14g\n", i, g[i]); 00641 # endif 00642 # ifndef NDEBUG 00643 for(i = 0; i < m_; i++) CPPAD_ASSERT_KNOWN( 00644 (-infinity_ < g[i]) && (g[i] < infinity_), 00645 "cppad_ipopt_nlp::eval_g:: not all constraints are not finite" 00646 ); 00647 # endif 00648 return true; 00649 } 00650 00651 /*! 00652 Evaluate the Jacobian of g(x). 00653 00654 \param[in] n 00655 is the dimension of the argument space for g(x); i.e., must be equal \c n_. 00656 00657 \param x 00658 if \c values is not \c NULL, 00659 \c x is a vector of size \c n containing the point at which to evaluate 00660 the gradient of g(x). 00661 00662 \param[in] new_x 00663 is false if the previous call to any one of the 00664 \ref Evaluation_Methods used the same value for \c x. 00665 00666 \param[in] m 00667 is the dimension of the range space for g(x); i.e., must be equal to \c m_. 00668 00669 \param[in] nele_jac 00670 is the number of possibly non-zero elements in the Jacobian of g(x); 00671 i.e., must be equal to \c nnz_jac_g_. 00672 00673 \param iRow 00674 if \c values is not \c NULL, \c iRow is not defined. 00675 if \c values is \c NULL, \c iRow 00676 is a vector with size \c nele_jac. 00677 The input value of its elements does not matter. 00678 On output, 00679 For <tt>k = 0 , ... , nele_jac-1, iRow[k]</tt> is the 00680 base zero row index for the 00681 k-th possibly non-zero entry in the Jacobian of g(x). 00682 00683 \param jCol 00684 if \c values is not \c NULL, \c jCol is not defined. 00685 if \c values is \c NULL, \c jCol 00686 is a vector with size \c nele_jac. 00687 The input value of its elements does not matter. 00688 On output, 00689 For <tt>k = 0 , ... , nele_jac-1, jCol[k]</tt> is the 00690 base zero column index for the 00691 k-th possibly non-zero entry in the Jacobian of g(x). 00692 00693 \param values 00694 if \c values is not \c NULL, \c values 00695 is a vector with size \c nele_jac. 00696 The input value of its elements does not matter. 00697 On output, 00698 For <tt>k = 0 , ... , nele_jac-1, values[k]</tt> is the 00699 value for the 00700 k-th possibly non-zero entry in the Jacobian of g(x). 00701 00702 \return 00703 The return value is always true; see \ref Evaluation_Methods. 00704 */ 00705 bool cppad_ipopt_nlp::eval_jac_g(Index n, const Number* x, bool new_x, 00706 Index m, Index nele_jac, Index* iRow, Index *jCol, 00707 Number* values) 00708 { CPPAD_ASSERT_UNKNOWN(size_t(m) == m_ ); 00709 CPPAD_ASSERT_UNKNOWN(size_t(n) == n_ ); 00710 CPPAD_ASSERT_UNKNOWN( size_t(nele_jac) == nnz_jac_g_ ); 00711 00712 size_t i, j, k, ell, l; 00713 std::map<size_t,size_t>::iterator index_ij; 00714 00715 00716 if (values == NULL) 00717 { for(k = 0; k < nnz_jac_g_; k++) 00718 { iRow[k] = iRow_jac_g_[k]; 00719 jCol[k] = jCol_jac_g_[k]; 00720 } 00721 return true; 00722 } 00723 00724 // initialize summation 00725 l = nnz_jac_g_; 00726 while(l--) 00727 values[l] = 0.; 00728 00729 // update tape_ok_ flag 00730 for(k = 0; k < K_; k++) 00731 { if( retape_[k] && (new_x || L_[k] > 1) ) 00732 tape_ok_[k] = false; 00733 } 00734 00735 for(k = 0; k < K_; k++) for(ell = 0; ell < L_[k]; ell++) 00736 { fg_info_->index(k, ell, I_, J_); 00737 if( ! tape_ok_[k] ) 00738 { // Record r_k for value of u corresponding to x 00739 fun_record( 00740 fg_info_ , // inputs 00741 k , 00742 p_ , 00743 q_ , 00744 n_ , 00745 x , 00746 J_ , 00747 r_fun_ // output 00748 ); 00749 } 00750 tape_ok_[k] = ! (retape_[k] || L_[k] > 1); 00751 NumberVector u(q_[k]); 00752 NumberVector jac_r(p_[k] * q_[k]); 00753 for(j = 0; j < q_[k]; j++) 00754 { CPPAD_ASSERT_UNKNOWN( J_[j] < n_ ); 00755 u[j] = x[ J_[j] ]; 00756 } 00757 if( retape_[k] ) 00758 jac_r = r_fun_[k].Jacobian(u); 00759 else jac_r = r_fun_[k].SparseJacobian(u, pattern_jac_r_[k]); 00760 for(i = 0; i < p_[k]; i++) if( I_[i] != 0 ) 00761 { CPPAD_ASSERT_UNKNOWN( I_[i] <= m_ ); 00762 for(j = 0; j < q_[k]; j++) 00763 { index_ij = index_jac_g_[I_[i]-1].find(J_[j]); 00764 if( index_ij != index_jac_g_[I_[i]-1].end() ) 00765 { l = index_ij->second; 00766 values[l] += jac_r[i * q_[k] + j]; 00767 } 00768 else CPPAD_ASSERT_UNKNOWN( 00769 jac_r[i * q_[k] + j] == 0. 00770 ); 00771 } 00772 } 00773 } 00774 # ifndef NDEBUG 00775 for(l = 0; l < nnz_jac_g_; l++) CPPAD_ASSERT_KNOWN( 00776 (-infinity_ < values[l]) && (values[l] < infinity_), 00777 "cppad_ipopt_nlp::eval_jac_g:: a component of " 00778 "gradient of g is not finite" 00779 ); 00780 # endif 00781 00782 return true; 00783 } 00784 00785 /*! 00786 Evaluate the Hessian of the Lagragian 00787 00788 \section The_Hessian_of_the_Lagragian The Hessian of the Lagragian 00789 The Hessian of the Lagragian is defined as 00790 \f[ 00791 H(x, \sigma, \lambda ) 00792 = 00793 \sigma \nabla^2 f(x) + \sum_{i=0}^{m-1} \lambda_i \nabla^2 g(x)_i 00794 \f] 00795 00796 \param[in] n 00797 is the dimension of the argument space for g(x); i.e., must be equal \c n_. 00798 00799 \param x 00800 if \c values is not \c NULL, \c x 00801 is a vector of size \c n containing the point at which to evaluate 00802 the Hessian of the Lagrangian. 00803 00804 \param[in] new_x 00805 is false if the previous call to any one of the 00806 \ref Evaluation_Methods used the same value for \c x. 00807 00808 \param[in] obj_factor 00809 the value \f$ \sigma \f$ multiplying the Hessian of 00810 f(x) in the expression for \ref The_Hessian_of_the_Lagragian. 00811 00812 \param[in] m 00813 is the dimension of the range space for g(x); i.e., must be equal to \c m_. 00814 00815 \param[in] lambda 00816 if \c values is not \c NULL, \c lambda 00817 is a vector of size \c m specifing the value of \f$ \lambda \f$ 00818 in the expression for \ref The_Hessian_of_the_Lagragian. 00819 00820 \param[in] new_lambda 00821 is true if the previous call to \c eval_h had the same value for 00822 \c lambda and false otherwise. 00823 (Not currently used.) 00824 00825 \param[in] nele_hess 00826 is the number of possibly non-zero elements in the Hessian of the Lagragian; 00827 i.e., must be equal to \c nnz_h_lag_. 00828 00829 \param iRow 00830 if \c values is not \c NULL, \c iRow is not defined. 00831 if \c values is \c NULL, \c iRow 00832 is a vector with size \c nele_jac. 00833 The input value of its elements does not matter. 00834 On output, 00835 For <tt>k = 0 , ... , nele_jac-1, iRow[k]</tt> is the 00836 base zero row index for the 00837 k-th possibly non-zero entry in the Hessian of the Lagragian. 00838 00839 \param jCol 00840 if \c values is not \c NULL, \c jCol is not defined. 00841 if \c values is \c NULL, \c jCol 00842 is a vector with size \c nele_jac. 00843 The input value of its elements does not matter. 00844 On output, 00845 For <tt>k = 0 , ... , nele_jac-1, jCol[k]</tt> is the 00846 base zero column index for the 00847 k-th possibly non-zero entry in the Hessian of the Lagragian. 00848 00849 \param values 00850 if \c values is not \c NULL, it 00851 is a vector with size \c nele_jac. 00852 The input value of its elements does not matter. 00853 On output, 00854 For <tt>k = 0 , ... , nele_jac-1, values[k]</tt> is the 00855 value for the 00856 k-th possibly non-zero entry in the Hessian of the Lagragian. 00857 00858 \return 00859 The return value is always true; see \ref Evaluation_Methods. 00860 */ 00861 bool cppad_ipopt_nlp::eval_h(Index n, const Number* x, bool new_x, 00862 Number obj_factor, Index m, const Number* lambda, 00863 bool new_lambda, Index nele_hess, Index* iRow, 00864 Index* jCol, Number* values) 00865 { CPPAD_ASSERT_UNKNOWN(size_t(m) == m_ ); 00866 CPPAD_ASSERT_UNKNOWN(size_t(n) == n_ ); 00867 00868 size_t i, j, k, ell, l; 00869 std::map<size_t,size_t>::iterator index_ij; 00870 00871 if (values == NULL) 00872 { for(k = 0; k < nnz_h_lag_; k++) 00873 { iRow[k] = iRow_h_lag_[k]; 00874 jCol[k] = jCol_h_lag_[k]; 00875 } 00876 return true; 00877 } 00878 00879 // initialize summation 00880 l = nnz_h_lag_; 00881 while(l--) 00882 values[l] = 0.; 00883 00884 // update tape_ok_ flag 00885 for(k = 0; k < K_; k++) 00886 { if( retape_[k] && (new_x || L_[k] > 1) ) 00887 tape_ok_[k] = false; 00888 } 00889 00890 for(k = 0; k < K_; k++) for(ell = 0; ell < L_[k]; ell++) 00891 { fg_info_->index(k, ell, I_, J_); 00892 bool in_use = false; 00893 for(i = 0; i < p_[k]; i++) 00894 { if( I_[i] == 0 ) 00895 in_use |= obj_factor > 0.; 00896 else in_use |= lambda[ I_[i] - 1 ] > 0; 00897 } 00898 if( in_use ) 00899 { 00900 if( ! tape_ok_[k] ) 00901 { // Record r_k for value of u corresponding to x 00902 fun_record( 00903 fg_info_ , // inputs 00904 k , 00905 p_ , 00906 q_ , 00907 n_ , 00908 x , 00909 J_ , 00910 r_fun_ // output 00911 ); 00912 tape_ok_[k] = ! (retape_[k] || L_[k] > 1); 00913 } 00914 NumberVector w(p_[k]); 00915 NumberVector r_hes(q_[k] * q_[k]); 00916 NumberVector u(q_[k]); 00917 for(j = 0; j < q_[k]; j++) 00918 { CPPAD_ASSERT_UNKNOWN( J_[j] < n_ ); 00919 u[j] = x[ J_[j] ]; 00920 } 00921 for(i = 0; i < p_[k]; i++) 00922 { CPPAD_ASSERT_UNKNOWN( I_[i] <= m_ ); 00923 if( I_[i] == 0 ) 00924 w[i] = obj_factor; 00925 else w[i] = lambda[ I_[i] - 1 ]; 00926 } 00927 if( retape_[k] ) 00928 r_hes = r_fun_[k].Hessian(u, w); 00929 else r_hes = r_fun_[k].SparseHessian( 00930 u, w, pattern_hes_r_[k] 00931 ); 00932 for(i = 0; i < q_[k]; i++) for(j = 0; j < q_[k]; j++) 00933 if( J_[j] <= J_[i] ) 00934 { index_ij = index_hes_fg_[J_[i]].find(J_[j]); 00935 if( index_ij != index_hes_fg_[J_[i]].end() ) 00936 { l = index_ij->second; 00937 values[l] += r_hes[i * q_[k] + j]; 00938 } 00939 else CPPAD_ASSERT_UNKNOWN( 00940 r_hes[i * q_[k] + j] == 0. 00941 ); 00942 } 00943 } 00944 } 00945 # ifndef NDEBUG 00946 for(l = 0; l < nnz_h_lag_; l++) CPPAD_ASSERT_KNOWN( 00947 (-infinity_ < values[l]) && (values[l] < infinity_), 00948 "cppad_ipopt_nlp::eval_h:: a component of " 00949 "Hessian of Lagragian is not finite" 00950 ); 00951 # endif 00952 return true; 00953 } 00954 00955 /*! 00956 Pass solution information from Ipopt to users solution structure. 00957 00958 \param[in] status 00959 is value that the Ipopt solution status 00960 which gets mapped to a correponding value for 00961 \n 00962 <tt>solution_->status</tt>. 00963 00964 \param[in] n 00965 is the dimension of the domain space for f(x) and g(x); i.e., 00966 it must be equal to \c n_. 00967 00968 \param[in] x 00969 is a vector with size \c n specifing the final solution. 00970 \n 00971 <tt>solution_->x</tt> is set to be a vector with size \c n 00972 and to have the same element values. 00973 00974 \param[in] z_L 00975 is a vector with size \c n specifing the Lagragian multipliers for the 00976 constraint \f$ x^l \leq x \f$. 00977 \n 00978 <tt>solution_->z_l</tt> is set to be a vector with size \c n 00979 and to have the same element values. 00980 00981 \param[in] z_U 00982 is a vector with size \c n specifing the Lagragian multipliers for the 00983 constraint \f$ x \leq x^u \f$. 00984 \n 00985 <tt>solution_->z_u</tt> is set to be a vector with size \c n 00986 and to have the same element values. 00987 00988 \param[in] m 00989 is the dimension of the range space for g(x). i.e., 00990 it must be equal to \c m_. 00991 00992 \param[in] g 00993 is a vector with size \c m containing the value of the constraint function 00994 g(x) at the final solution for \c x. 00995 \n 00996 <tt>solution_->g</tt> is set to be a vector with size \c m 00997 and to have the same element values. 00998 00999 \param[in] lambda 01000 is a vector with size \c m specifing the Lagragian multipliers for the 01001 constraints \f$ g^l \leq g(x) \leq g^u \f$. 01002 \n 01003 <tt>solution_->lambda</tt> is set to be a vector with size \c m 01004 and to have the same element values. 01005 01006 \param[in] obj_value 01007 is the value of the objective function f(x) at the final solution for \c x. 01008 \n 01009 <tt>solution_->obj_value</tt> is set to have the same value. 01010 01011 \param[in] ip_data 01012 is unspecified (by Ipopt) and hence not used. 01013 01014 \param[in] ip_cq 01015 is unspecified (by Ipopt) and hence not used. 01016 01017 \par solution_[out] 01018 the pointer \c solution_ , which is equal to the pointer \c solution 01019 in the constructor for \c cppad_ipopt_nlp, 01020 is used to set output values (see documentation above). 01021 */ 01022 void cppad_ipopt_nlp::finalize_solution( 01023 Ipopt::SolverReturn status , 01024 Index n , 01025 const Number* x , 01026 const Number* z_L , 01027 const Number* z_U , 01028 Index m , 01029 const Number* g , 01030 const Number* lambda , 01031 Number obj_value , 01032 const Ipopt::IpoptData* ip_data , 01033 Ipopt::IpoptCalculatedQuantities* ip_cq 01034 ) 01035 { size_t i, j; 01036 01037 CPPAD_ASSERT_UNKNOWN(size_t(n) == n_ ); 01038 CPPAD_ASSERT_UNKNOWN(size_t(m) == m_ ); 01039 01040 switch(status) 01041 { // convert status from Ipopt enum to cppad_ipopt_solution enum 01042 case Ipopt::SUCCESS: 01043 solution_->status = 01044 cppad_ipopt_solution::success; 01045 break; 01046 01047 case Ipopt::MAXITER_EXCEEDED: 01048 solution_->status = 01049 cppad_ipopt_solution::maxiter_exceeded; 01050 break; 01051 01052 case Ipopt::STOP_AT_TINY_STEP: 01053 solution_->status = 01054 cppad_ipopt_solution::stop_at_tiny_step; 01055 break; 01056 01057 case Ipopt::STOP_AT_ACCEPTABLE_POINT: 01058 solution_->status = 01059 cppad_ipopt_solution::stop_at_acceptable_point; 01060 break; 01061 01062 case Ipopt::LOCAL_INFEASIBILITY: 01063 solution_->status = 01064 cppad_ipopt_solution::local_infeasibility; 01065 break; 01066 01067 case Ipopt::USER_REQUESTED_STOP: 01068 solution_->status = 01069 cppad_ipopt_solution::user_requested_stop; 01070 break; 01071 01072 case Ipopt::DIVERGING_ITERATES: 01073 solution_->status = 01074 cppad_ipopt_solution::diverging_iterates; 01075 break; 01076 01077 case Ipopt::RESTORATION_FAILURE: 01078 solution_->status = 01079 cppad_ipopt_solution::restoration_failure; 01080 break; 01081 01082 case Ipopt::ERROR_IN_STEP_COMPUTATION: 01083 solution_->status = 01084 cppad_ipopt_solution::error_in_step_computation; 01085 break; 01086 01087 case Ipopt::INVALID_NUMBER_DETECTED: 01088 solution_->status = 01089 cppad_ipopt_solution::invalid_number_detected; 01090 break; 01091 01092 case Ipopt::INTERNAL_ERROR: 01093 solution_->status = 01094 cppad_ipopt_solution::internal_error; 01095 break; 01096 01097 default: 01098 solution_->status = 01099 cppad_ipopt_solution::unknown; 01100 } 01101 01102 solution_->x.resize(n_); 01103 solution_->z_l.resize(n_); 01104 solution_->z_u.resize(n_); 01105 for(j = 0; j < n_; j++) 01106 { solution_->x[j] = x[j]; 01107 solution_->z_l[j] = z_L[j]; 01108 solution_->z_u[j] = z_U[j]; 01109 } 01110 solution_->g.resize(m_); 01111 solution_->lambda.resize(m_); 01112 for(i = 0; i < m_; i++) 01113 { solution_->g[i] = g[i]; 01114 solution_->lambda[i] = lambda[i]; 01115 } 01116 solution_->obj_value = obj_value; 01117 return; 01118 } 01119 01120 // This routine is defined, but not yet used 01121 // (trying to figure out a problem with Ipopt-3.9.1 and dismod4). 01122 bool cppad_ipopt_nlp::intermediate_callback( 01123 Ipopt::AlgorithmMode mode, 01124 Index iter, 01125 Number obj_value, 01126 Number inf_pr, 01127 Number inf_du, 01128 Number mu, 01129 Number d_norm, 01130 Number regularization_size, 01131 Number alpha_du, 01132 Number alpha_pr, 01133 Index ls_trials, 01134 const Ipopt::IpoptData* ip_data, 01135 Ipopt::IpoptCalculatedQuantities* ip_cq 01136 ) 01137 { 01138 // std::cout << "intermediate_callback called" << std::endl; 01139 return true; 01140 } 01141 01142 // --------------------------------------------------------------------------- 01143 } // end namespace cppad_ipopt 01144 // ---------------------------------------------------------------------------