00001
00020 #include "IpoptSolver.h"
00021 #include "IpIpoptApplication.hpp"
00022 #include "CommonUtil.h"
00023
00024 using std::cout;
00025 using std::endl;
00026 using std::ostringstream;
00027 using namespace Ipopt;
00028
00029
00030 IpoptSolver::IpoptSolver() {
00031 osrlwriter = new OSrLWriter();
00032 ipoptErrorMsg = "";
00033 }
00034
00035 IpoptSolver::~IpoptSolver() {
00036 #ifdef DEBUG
00037 cout << "inside IpoptSolver destructor" << endl;
00038 #endif
00039 delete osrlwriter;
00040 osrlwriter = NULL;
00041 #ifdef DEBUG
00042 cout << "leaving IpoptSolver destructor" << endl;
00043 #endif
00044 }
00045
00046
00047 bool IpoptSolver::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
00048 Index& nnz_h_lag, IndexStyleEnum& index_style)
00049 {
00050 if(osinstance->getObjectiveNumber() <= 0) throw ErrorClass("Ipopt NEEDS AN OBJECTIVE FUNCTION");
00051
00052 n = osinstance->getVariableNumber();
00053
00054 m = osinstance->getConstraintNumber();
00055 cout << "number variables !!!!!!!!!!!!!!!!!!!!!!!!!!!" << n << endl;
00056 cout << "number constraints !!!!!!!!!!!!!!!!!!!!!!!!!!!" << m << endl;
00057 try{
00058 osinstance->initForAlgDiff( );
00059 }
00060 catch(const ErrorClass& eclass){
00061 ipoptErrorMsg = eclass.errormsg;
00062 throw;
00063 }
00064
00065 osinstance->bUseExpTreeForFunEval = true;
00066 std::cout << "Call sparse jacobian" << std::endl;
00067 SparseJacobianMatrix *sparseJacobian = NULL;
00068 try{
00069 sparseJacobian = osinstance->getJacobianSparsityPattern();
00070 }
00071 catch(const ErrorClass& eclass){
00072 ipoptErrorMsg = eclass.errormsg;
00073 throw;
00074 }
00075 std::cout << "Done calling sparse jacobian" << std::endl;
00076 nnz_jac_g = sparseJacobian->valueSize;
00077
00078
00079
00080 if( (osinstance->getNumberOfNonlinearExpressions() == 0) && (osinstance->getNumberOfQuadraticTerms() == 0) ) {
00081 cout << "This is a linear program" << endl;
00082 nnz_h_lag = 0;
00083 }
00084 else{
00085 std::cout << "Get Lagrangain Hessian Sparsity Pattern " << std::endl;
00086 SparseHessianMatrix *sparseHessian = osinstance->getLagrangianHessianSparsityPattern();
00087 std::cout << "Done Getting Lagrangain Hessian Sparsity Pattern " << std::endl;
00088 nnz_h_lag = sparseHessian->hessDimension;
00089 }
00090
00091
00092 index_style = TNLP::C_STYLE;
00093
00095
00096 return true;
00097 }
00098
00099
00100 bool IpoptSolver::get_bounds_info(Index n, Number* x_l, Number* x_u,
00101 Index m, Number* g_l, Number* g_u){
00102 int i;
00103 double * mdVarLB = osinstance->getVariableLowerBounds();
00104 std::cout << "GET BOUNDS INFORMATION FOR IPOPT !!!!!!!!!!!!!!!!!" << std::endl;
00105
00106 double * mdVarUB = osinstance->getVariableUpperBounds();
00107
00108 for(i = 0; i < n; i++){
00109 x_l[ i] = mdVarLB[ i];
00110 x_u[ i] = mdVarUB[ i];
00111
00112
00113 }
00114
00115
00116
00117
00118
00119
00120 double * mdConLB = osinstance->getConstraintLowerBounds();
00121
00122 double * mdConUB = osinstance->getConstraintUpperBounds();
00123
00124 for(int i = 0; i < m; i++){
00125 g_l[ i] = mdConLB[ i];
00126 g_u[ i] = mdConUB[ i];
00127
00128
00129 }
00130 return true;
00131 }
00132
00133
00134
00135 bool IpoptSolver::get_starting_point(Index n, bool init_x, Number* x,
00136 bool init_z, Number* z_L, Number* z_U, Index m, bool init_lambda,
00137 Number* lambda) {
00138
00139
00140
00141 assert(init_x == true);
00142 assert(init_z == false);
00143 assert(init_lambda == false);
00144 int i;
00145
00146 double *mdXInit = osinstance->getVariableInitialValues();
00147 if( mdXInit != NULL) {
00148 for(i = 0; i < n; i++){
00149 x[ i] = mdXInit[ i];
00150
00151 }
00152 }
00153 else{
00154 for(i = 0; i < n; i++){
00155 x[ i] = 1.7171;
00156 }
00157 }
00158
00159 return true;
00160 }
00161
00162
00163 bool IpoptSolver::eval_f(Index n, const Number* x, bool new_x, Number& obj_value){
00164 try{
00165 obj_value = osinstance->calculateAllObjectiveFunctionValues( const_cast<double*>(x), NULL, NULL, new_x, 0 )[ 0];
00166 }
00167 catch(const ErrorClass& eclass){
00168 ipoptErrorMsg = eclass.errormsg;
00169 throw;
00170 }
00171 if( CommonUtil::ISOSNAN( (double)obj_value) ) return false;
00172 return true;
00173 }
00174
00175 bool IpoptSolver::eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f){
00176 int i;
00177 double *objGrad;
00178 try{
00179 objGrad = osinstance->calculateAllObjectiveFunctionGradients( const_cast<double*>(x), NULL, NULL, new_x, 1)[ 0];
00180 }
00181 catch(const ErrorClass& eclass){
00182 ipoptErrorMsg = eclass.errormsg;
00183 throw;
00184 }
00185 for(i = 0; i < n; i++){
00186 grad_f[ i] = objGrad[ i];
00187 }
00188 return true;
00189 }
00190
00191
00192 bool IpoptSolver::eval_g(Index n, const Number* x, bool new_x, Index m, Number* g) {
00193 try{
00194 double *conVals = osinstance->calculateAllConstraintFunctionValues( const_cast<double*>(x), NULL, NULL, new_x, 0 );
00195 int i;
00196 for(i = 0; i < m; i++){
00197 if( CommonUtil::ISOSNAN( (double)conVals[ i] ) ) return false;
00198 g[i] = conVals[ i] ;
00199 }
00200 return true;
00201 }
00202 catch(const ErrorClass& eclass){
00203 ipoptErrorMsg = eclass.errormsg;
00204 throw;
00205 }
00206 }
00207
00208
00209
00210 bool IpoptSolver::eval_jac_g(Index n, const Number* x, bool new_x,
00211 Index m, Index nele_jac, Index* iRow, Index *jCol,
00212 Number* values){
00213 SparseJacobianMatrix *sparseJacobian;
00214 if (values == NULL) {
00215
00216 cout << "n: " << n << endl;
00217 cout << "m: " << m << endl;
00218 cout << "nele_jac: " << nele_jac << endl;
00219
00220 try{
00221 sparseJacobian = osinstance->getJacobianSparsityPattern();
00222 }
00223 catch(const ErrorClass& eclass){
00224 ipoptErrorMsg = eclass.errormsg;
00225 throw;
00226 }
00227 int i = 0;
00228 int k, idx;
00229 for(idx = 0; idx < m; idx++){
00230 for(k = *(sparseJacobian->starts + idx); k < *(sparseJacobian->starts + idx + 1); k++){
00231 iRow[i] = idx;
00232 jCol[i] = *(sparseJacobian->indexes + k);
00233
00234
00235 i++;
00236 }
00237 }
00238 }
00239 else {
00240
00241 try{
00242 sparseJacobian = osinstance->calculateAllConstraintFunctionGradients( const_cast<double*>(x), NULL, NULL, new_x, 1);
00243 }
00244 catch(const ErrorClass& eclass){
00245 ipoptErrorMsg = eclass.errormsg;
00246 throw;
00247 }
00248
00249 for(int i = 0; i < nele_jac; i++){
00250 values[ i] = sparseJacobian->values[i];
00251
00252
00253
00254 }
00255 }
00256 return true;
00257 }
00258
00259
00260 bool IpoptSolver::eval_h(Index n, const Number* x, bool new_x,
00261 Number obj_factor, Index m, const Number* lambda,
00262 bool new_lambda, Index nele_hess, Index* iRow,
00263 Index* jCol, Number* values){
00264
00266 SparseHessianMatrix *sparseHessian;
00267
00268 int i;
00269 if (values == NULL) {
00270
00271
00272 try{
00273 sparseHessian = osinstance->getLagrangianHessianSparsityPattern( );
00274 }
00275 catch(const ErrorClass& eclass){
00276 ipoptErrorMsg = eclass.errormsg;
00277 throw;
00278 }
00279 cout << "got structure of HESSIAN !!!!!!!!!!!!!!!!!!!!!!!!!! " << endl;
00280 for(i = 0; i < nele_hess; i++){
00281 iRow[i] = *(sparseHessian->hessColIdx + i);
00282 jCol[i] = *(sparseHessian->hessRowIdx + i);
00283
00284
00285 }
00286 }
00287 else {
00288
00289
00290 double* objMultipliers = new double[1];
00291 objMultipliers[0] = obj_factor;
00292 try{
00293 sparseHessian = osinstance->calculateLagrangianHessian( const_cast<double*>(x), objMultipliers, (double*)lambda , new_x, 2);
00294 }
00295 catch(const ErrorClass& eclass){
00296 ipoptErrorMsg = eclass.errormsg;
00297 throw;
00298 }
00299 for(i = 0; i < nele_hess; i++){
00300 values[ i] = *(sparseHessian->hessValues + i);
00301 }
00302 }
00304 return true;
00305 }
00306
00307 bool IpoptSolver::get_scaling_parameters(Number& obj_scaling,
00308 bool& use_x_scaling, Index n,
00309 Number* x_scaling,
00310 bool& use_g_scaling, Index m,
00311 Number* g_scaling){
00312 if( osinstance->instanceData->objectives->obj[ 0]->maxOrMin.compare("min") != 0){
00313 obj_scaling = -1;
00314 }
00315 else obj_scaling = 1;
00316 use_x_scaling = false;
00317 use_g_scaling = false;
00318 return true;
00319 }
00320
00321 void IpoptSolver::finalize_solution(SolverReturn status,
00322 Index n, const Number* x, const Number* z_L, const Number* z_U,
00323 Index m, const Number* g, const Number* lambda,
00324 Number obj_value,
00325 const IpoptData* ip_data,
00326 IpoptCalculatedQuantities* ip_cq)
00327 {
00328
00329
00330
00331 printf("\n\nSolution of the primal variables, x\n");
00332 for (Index i=0; i<n; i++) {
00333 printf("x[%d] = %e\n", i, x[i]);
00334 }
00335
00336 printf("\n\nSolution of the bound multipliers, z_L and z_U\n");
00337 for (Index i=0; i<n; i++) {
00338 printf("z_L[%d] = %e\n", i, z_L[i]);
00339 }
00340 for (Index i=0; i<n; i++) {
00341 printf("z_U[%d] = %e\n", i, z_U[i]);
00342 }
00343
00344 printf("\n\nObjective value\n");
00345 printf("f(x*) = %e\n", obj_value);
00346 int solIdx = 0;
00347 ostringstream outStr;
00348 double* mdObjValues = new double[1];
00349 std::string message = "Ipopt solver finishes to the end.";
00350 std::string solutionDescription = "";
00351
00352 try{
00353
00354 if(osresult->setServiceName( "Ipopt solver service") != true)
00355 throw ErrorClass("OSResult error: setServiceName");
00356 if(osresult->setInstanceName( osinstance->getInstanceName()) != true)
00357 throw ErrorClass("OSResult error: setInstanceName");
00358
00359
00360
00361
00362
00363 if(osresult->setVariableNumber( osinstance->getVariableNumber()) != true)
00364 throw ErrorClass("OSResult error: setVariableNumer");
00365 if(osresult->setObjectiveNumber( 1) != true)
00366 throw ErrorClass("OSResult error: setObjectiveNumber");
00367 if(osresult->setConstraintNumber( osinstance->getConstraintNumber()) != true)
00368 throw ErrorClass("OSResult error: setConstraintNumber");
00369 if(osresult->setSolutionNumber( 1) != true)
00370 throw ErrorClass("OSResult error: setSolutionNumer");
00371
00372
00373 if(osresult->setGeneralMessage( message) != true)
00374 throw ErrorClass("OSResult error: setGeneralMessage");
00375
00376 switch( status){
00377 case SUCCESS:
00378 solutionDescription = "SUCCESS[IPOPT]: Algorithm terminated successfully at a locally optimal point, satisfying the convergence tolerances.";
00379 osresult->setSolutionStatus(solIdx, "locallyOptimal", solutionDescription);
00380 osresult->setPrimalVariableValues(solIdx, const_cast<double*>(x));
00381 mdObjValues[0] = obj_value;
00382 osresult->setObjectiveValues(solIdx, mdObjValues);
00383 break;
00384 case MAXITER_EXCEEDED:
00385 solutionDescription = "MAXITER_EXCEEDED[IPOPT]: Maximum number of iterations exceeded.";
00386 osresult->setSolutionStatus(solIdx, "stoppedByLimit", solutionDescription);
00387 osresult->setPrimalVariableValues(solIdx, const_cast<double*>(x));
00388 mdObjValues[0] = obj_value;
00389 osresult->setObjectiveValues(solIdx, mdObjValues);
00390 break;
00391 case STOP_AT_TINY_STEP:
00392 solutionDescription = "STOP_AT_TINY_STEP[IPOPT]: Algorithm proceeds with very little progress.";
00393 osresult->setSolutionStatus(solIdx, "stoppedByLimit", solutionDescription);
00394 osresult->setPrimalVariableValues(solIdx, const_cast<double*>(x));
00395 mdObjValues[0] = obj_value;
00396 osresult->setObjectiveValues(solIdx, mdObjValues);
00397 break;
00398 case STOP_AT_ACCEPTABLE_POINT:
00399 solutionDescription = "STOP_AT_ACCEPTABLE_POINT[IPOPT]: Algorithm stopped at a point that was converged, not to _desired_ tolerances, but to _acceptable_ tolerances";
00400 osresult->setSolutionStatus(solIdx, "locallyOptimal", solutionDescription);
00401 osresult->setPrimalVariableValues(solIdx, const_cast<double*>(x));
00402 mdObjValues[0] = obj_value;
00403 osresult->setObjectiveValues(solIdx, mdObjValues);
00404 break;
00405 case LOCAL_INFEASIBILITY:
00406 solutionDescription = "LOCAL_INFEASIBILITY[IPOPT]: Algorithm converged to a point of local infeasibility. Problem may be infeasible.";
00407 osresult->setSolutionStatus(solIdx, "infeasible", solutionDescription);
00408 break;
00409 case USER_REQUESTED_STOP:
00410 solutionDescription = "USER_REQUESTED_STOP[IPOPT]: The user call-back function intermediate_callback returned false, i.e., the user code requested a premature termination of the optimization.";
00411 osresult->setSolutionStatus(solIdx, "error", solutionDescription);
00412 break;
00413 case DIVERGING_ITERATES:
00414 solutionDescription = "DIVERGING_ITERATES[IPOPT]: It seems that the iterates diverge.";
00415 osresult->setSolutionStatus(solIdx, "unbounded", solutionDescription);
00416 break;
00417 case RESTORATION_FAILURE:
00418 solutionDescription = "RESTORATION_FAILURE[IPOPT]: Restoration phase failed, algorithm doesn't know how to proceed.";
00419 osresult->setSolutionStatus(solIdx, "error", solutionDescription);
00420 break;
00421 case ERROR_IN_STEP_COMPUTATION:
00422 solutionDescription = "ERROR_IN_STEP_COMPUTATION[IPOPT]: An unrecoverable error occurred while IPOPT tried to compute the search direction.";
00423 osresult->setSolutionStatus(solIdx, "error", solutionDescription);
00424 break;
00425 case INVALID_NUMBER_DETECTED:
00426 solutionDescription = "INVALID_NUMcatBER_DETECTED[IPOPT]: Algorithm received an invalid number (such as NaN or Inf) from the NLP; see also option check_derivatives_for_naninf.";
00427 osresult->setSolutionStatus(solIdx, "error", solutionDescription);
00428 break;
00429 case INTERNAL_ERROR:
00430 solutionDescription = "INTERNAL_ERROR[IPOPT]: An unknown internal error occurred. Please contact the IPOPT authors through the mailing list.";
00431 osresult->setSolutionStatus(solIdx, "error", solutionDescription);
00432 break;
00433 default:
00434 solutionDescription = "OTHER[IPOPT]: other unknown solution status from Ipopt solver";
00435 osresult->setSolutionStatus(solIdx, "other", solutionDescription);
00436 }
00437
00438 osresult->setGeneralStatusType("success");
00439 osrl = osrlwriter->writeOSrL( osresult);
00440
00441 }
00442 catch(const ErrorClass& eclass){
00443 osresult->setGeneralMessage( eclass.errormsg);
00444 osresult->setGeneralStatusType( "error");
00445 osrl = osrlwriter->writeOSrL( osresult);
00446 throw ;
00447 }
00449 }
00450
00451
00452
00453 void IpoptSolver::solve() throw (ErrorClass) {
00454 try{
00455 OSiLReader* osilreader = NULL;
00456 osresult = new OSResult();
00457 if(osil.length() == 0 && osinstance == NULL) throw ErrorClass("there is no instance");
00458 clock_t start, finish;
00459 double duration;
00460 start = clock();
00461 if(osinstance == NULL){
00462 osilreader = new OSiLReader();
00463 osinstance = osilreader->readOSiL( &osil);
00464 }
00465 OSiLWriter osilwriter;
00466
00467 if(osinstance->getVariableNumber() <= 0)throw ErrorClass("Ipopt requires decision variables");
00468 finish = clock();
00469 duration = (double) (finish - start) / CLOCKS_PER_SEC;
00470 cout << "Parsing took (seconds): "<< duration << endl;
00471
00472
00473
00474
00475
00476
00477 SmartPtr<TNLP> nlp = this;
00478
00479
00480
00481 SmartPtr<IpoptApplication> app = new IpoptApplication();
00482
00483
00484
00485 app->Options()->SetNumericValue("tol", 1e-9);
00486 app->Options()->SetStringValue("mu_strategy", "adaptive");
00487 app->Options()->SetStringValue("output_file", "ipopt.out");
00488 app->Options()->SetStringValue("check_derivatives_for_naninf", "yes");
00489
00490 if(osinstance->getObjectiveNumber() <= 0) throw ErrorClass("Ipopt NEEDS AN OBJECTIVE FUNCTION");
00491 if( (osinstance->getNumberOfNonlinearExpressions() == 0) && (osinstance->getNumberOfQuadraticTerms() == 0) )
00492 app->Options()->SetStringValue("hessian_approximation", "limited-memory");
00493
00494 if( osinstance->instanceData->objectives->obj[ 0]->maxOrMin.compare("min") != 0){
00495 app->Options()->SetStringValue("nlp_scaling_method", "user-scaling");
00496 }
00497
00498 std::cout << "Call Ipopt Initialize" << std::endl;
00499 app->Initialize();
00500 std::cout << "Finished Ipopt Initialize" << std::endl;
00501
00502 std::cout << "Call Ipopt Optimize" << std::endl;
00503 ApplicationReturnStatus status = app->OptimizeTNLP(nlp);
00504 std::cout << "Finish Ipopt Optimize" << std::endl;
00505 if (status != Solve_Succeeded) {
00506 throw ErrorClass("Ipopt FAILED TO SOLVE THE PROBLEM: " + ipoptErrorMsg);
00507 }
00508 delete osilreader;
00509 osilreader = NULL;
00510 }
00511 catch(const ErrorClass& eclass){
00512 osresult->setGeneralMessage( eclass.errormsg);
00513 osresult->setGeneralStatusType( "error");
00514 osrl = osrlwriter->writeOSrL( osresult);
00515 throw ErrorClass( osrl) ;
00516 }
00517 }
00518
00519
00520 void IpoptSolver::dataEchoCheck(){
00521
00522 int i;
00523
00524
00525 cout << "This is problem: " << osinstance->getInstanceName() << endl;
00526 cout << "The problem source is: " << osinstance->getInstanceSource() << endl;
00527 cout << "The problem description is: " << osinstance->getInstanceDescription() << endl;
00528 cout << "number of variables = " << osinstance->getVariableNumber() << endl;
00529 cout << "number of Rows = " << osinstance->getConstraintNumber() << endl;
00530
00531
00532 if(osinstance->getVariableNumber() > 0){
00533 for(i = 0; i < osinstance->getVariableNumber(); i++){
00534 if(osinstance->getVariableNames() != NULL) cout << "variable Names " << osinstance->getVariableNames()[ i] << endl;
00535 if(osinstance->getVariableTypes() != NULL) cout << "variable Types " << osinstance->getVariableTypes()[ i] << endl;
00536 if(osinstance->getVariableLowerBounds() != NULL) cout << "variable Lower Bounds " << osinstance->getVariableLowerBounds()[ i] << endl;
00537 if(osinstance->getVariableUpperBounds() != NULL) cout << "variable Upper Bounds " << osinstance->getVariableUpperBounds()[i] << endl;
00538 }
00539 }
00540
00541
00542 if(osinstance->getVariableNumber() > 0 || osinstance->instanceData->objectives->obj != NULL || osinstance->instanceData->objectives->numberOfObjectives > 0){
00543 if( osinstance->getObjectiveMaxOrMins()[0] == "min") cout << "problem is a minimization" << endl;
00544 else cout << "problem is a maximization" << endl;
00545 for(i = 0; i < osinstance->getVariableNumber(); i++){
00546 cout << "OBJ COEFFICIENT = " << osinstance->getDenseObjectiveCoefficients()[0][i] << endl;
00547 }
00548 }
00549
00550 if(osinstance->getConstraintNumber() > 0){
00551 for(i = 0; i < osinstance->getConstraintNumber(); i++){
00552 if(osinstance->getConstraintNames() != NULL) cout << "row name = " << osinstance->getConstraintNames()[i] << endl;
00553 if(osinstance->getConstraintLowerBounds() != NULL) cout << "row lower bound = " << osinstance->getConstraintLowerBounds()[i] << endl;
00554 if(osinstance->getConstraintUpperBounds() != NULL) cout << "row upper bound = " << osinstance->getConstraintUpperBounds()[i] << endl;
00555 }
00556 }
00557
00558
00559 cout << endl;
00560 cout << "number of nonzeros = " << osinstance->getLinearConstraintCoefficientNumber() << endl;
00561 for(i = 0; i <= osinstance->getVariableNumber(); i++){
00562 cout << "Start Value = " << osinstance->getLinearConstraintCoefficientsInColumnMajor()->starts[ i] << endl;
00563 }
00564 cout << endl;
00565 for(i = 0; i < osinstance->getLinearConstraintCoefficientNumber(); i++){
00566 cout << "Index Value = " << osinstance->getLinearConstraintCoefficientsInColumnMajor()->indexes[i] << endl;
00567 cout << "Nonzero Value = " << osinstance->getLinearConstraintCoefficientsInColumnMajor()->values[i] << endl;
00568 }
00569
00570
00571 cout << "number of qterms = " << osinstance->getNumberOfQuadraticTerms() << endl;
00572 for(int i = 0; i < osinstance->getNumberOfQuadraticTerms(); i++){
00573 cout << "Row Index = " << osinstance->getQuadraticTerms()->rowIndexes[i] << endl;
00574 cout << "Var Index 1 = " << osinstance->getQuadraticTerms()->varOneIndexes[ i] << endl;
00575 cout << "Var Index 2 = " << osinstance->getQuadraticTerms()->varTwoIndexes[ i] << endl;
00576 cout << "Coefficient = " << osinstance->getQuadraticTerms()->coefficients[ i] << endl;
00577 }
00578 }
00579
00580
00581
00582
00583