/home/coin/SVN-release/OS-2.4.2/OS/src/OSSolverInterfaces/OSIpoptSolver.cpp

Go to the documentation of this file.
00001 /* $Id: OSIpoptSolver.cpp 4292 2011-09-21 05:47:18Z kmartin $ */
00018 //#define DEBUG
00019 
00020 #include "OSIpoptSolver.h"
00021 #include "OSGeneral.h"
00022 #include "OSParameters.h"
00023 #include "OSMathUtil.h"
00024 #include "CoinFinite.hpp"
00025 
00026 
00027 using std::cout;
00028 using std::endl;
00029 using std::ostringstream;
00030 using namespace Ipopt;
00031 
00032 
00033 IpoptSolver::IpoptSolver()
00034 {
00035     osrlwriter = new OSrLWriter();
00036     osresult = new OSResult();
00037     m_osilreader = NULL;
00038     m_osolreader = NULL;
00039     ipoptErrorMsg = new std::string("");
00040 }
00041 
00042 IpoptSolver::~IpoptSolver()
00043 {
00044 #ifdef DEBUG
00045     cout << "inside IpoptSolver destructor" << endl;
00046 #endif
00047     if(m_osilreader != NULL) delete m_osilreader;
00048     m_osilreader = NULL;
00049     if(m_osolreader != NULL) delete m_osolreader;
00050     m_osolreader = NULL;
00051     delete osresult;
00052     osresult = NULL;
00053     delete osrlwriter;
00054     osrlwriter = NULL;
00055     //delete osinstance;
00056     //osinstance = NULL;
00057     delete ipoptErrorMsg;
00058 #ifdef DEBUG
00059     cout << "leaving IpoptSolver destructor" << endl;
00060 #endif
00061 }
00062 
00063 // returns the size of the problem
00064 bool IpoptProblem::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
00065                                 Index& nnz_h_lag, IndexStyleEnum& index_style)
00066 {
00067     try
00068     {
00069         //if(osinstance->getObjectiveNumber() <= 0) throw ErrorClass("Ipopt NEEDS AN OBJECTIVE FUNCTION");
00070         if( (osinstance->getNumberOfIntegerVariables() + osinstance->getNumberOfBinaryVariables()) > 0 )
00071             throw ErrorClass("Ipopt does not solve integer programs -- please try Bonmin or Couenne");
00072         // number of variables
00073         n = osinstance->getVariableNumber();
00074         // number of constraints
00075         m = osinstance->getConstraintNumber();
00076 #ifdef DEBUG
00077         cout << "number variables  !!!!!!!!!!!!!!!!!!!!!!!!!!!" << n << endl;
00078         cout << "number constraints  !!!!!!!!!!!!!!!!!!!!!!!!!!!" << m << endl;
00079 #endif
00080         try
00081         {
00082             osinstance->initForAlgDiff( );
00083         }
00084         catch(const ErrorClass& eclass)
00085         {
00086 #ifdef DEBUG
00087             cout << "error in OSIpoptSolver, line 78:\n" << eclass.errormsg << endl;
00088 #endif
00089             *ipoptErrorMsg = eclass.errormsg;
00090             throw;
00091         }
00092         // use the OS Expression tree for function evaluations instead of CppAD
00093         osinstance->bUseExpTreeForFunEval = true;
00094         //std::cout << "Call sparse jacobian" << std::endl;
00095         SparseJacobianMatrix *sparseJacobian = NULL;
00096         try
00097         {
00098             sparseJacobian = osinstance->getJacobianSparsityPattern();
00099         }
00100         catch(const ErrorClass& eclass)
00101         {
00102 #ifdef DEBUG
00103             cout << "error in OSIpoptSolver, line 91:\n" << eclass.errormsg << endl;
00104 #endif
00105             *ipoptErrorMsg = eclass.errormsg;
00106             throw;
00107         }
00108         //std::cout << "Done calling sparse jacobian" << std::endl;
00109         if (sparseJacobian != NULL)
00110         {
00111             nnz_jac_g = sparseJacobian->valueSize;
00112         }
00113         else
00114         {
00115             nnz_jac_g = 0;
00116         }
00117 
00118 #ifdef DEBUG
00119         cout << "nnz_jac_g  !!!!!!!!!!!!!!!!!!!!!!!!!!!" << nnz_jac_g << endl;
00120 #endif
00121         // nonzeros in upper hessian
00122 
00123         if( (osinstance->getNumberOfNonlinearExpressions() == 0) && (osinstance->getNumberOfQuadraticTerms() == 0) )
00124         {
00125 #ifdef DEBUG
00126             cout << "This is a linear program"  << endl;
00127 #endif
00128             nnz_h_lag = 0;
00129         }
00130         else
00131         {
00132             //std::cout << "Get Lagrangain Hessian Sparsity Pattern " << std::endl;
00133             SparseHessianMatrix *sparseHessian = osinstance->getLagrangianHessianSparsityPattern();
00134             //std::cout << "Done Getting Lagrangain Hessian Sparsity Pattern " << std::endl;
00135             if(sparseHessian != NULL)
00136             {
00137                 nnz_h_lag = sparseHessian->hessDimension;
00138             }
00139             else
00140             {
00141                 nnz_h_lag = 0;
00142             }
00143         }
00144 #ifdef DEBUG
00145         cout << "print nnz_h_lag (OSIpoptSolver.cpp)" << endl;
00146         cout << "nnz_h_lag  !!!!!!!!!!!!!!!!!!!!!!!!!!!" << nnz_h_lag << endl;
00147         cout << "set index_style (OSIpoptSolver.cpp)" << endl;
00148 #endif
00149         // use the C style indexing (0-based)
00150         index_style = TNLP::C_STYLE;
00151 #ifdef DEBUG
00152         cout << "return from get_nlp_info (OSIpoptSolver.cpp)" << nnz_h_lag << endl;
00153 #endif
00154 
00156 
00157         return true;
00158     }
00159     catch(const ErrorClass& eclass)
00160     {
00161 
00162         *ipoptErrorMsg = eclass.errormsg;
00163         throw;
00164     }
00165 
00166 }//get_nlp_info
00167 
00168 
00169 bool  IpoptProblem::get_bounds_info(Index n, Number* x_l, Number* x_u,
00170                                     Index m, Number* g_l, Number* g_u)
00171 {
00172     int i;
00173     double * mdVarLB = osinstance->getVariableLowerBounds();
00174     //std::cout << "GET BOUNDS INFORMATION FOR IPOPT !!!!!!!!!!!!!!!!!" << std::endl;
00175     // variables upper bounds
00176     double * mdVarUB = osinstance->getVariableUpperBounds();
00177 
00178     for(i = 0; i < n; i++)
00179     {
00180         x_l[ i] = mdVarLB[ i];
00181         x_u[ i] = mdVarUB[ i];
00182         //cout << "x_l !!!!!!!!!!!!!!!!!!!!!!!!!!!" << x_l[i] << endl;
00183         //cout << "x_u !!!!!!!!!!!!!!!!!!!!!!!!!!!" << x_u[i] << endl;
00184     }
00185     // Ipopt interprets any number greater than nlp_upper_bound_inf as
00186     // infinity. The default value of nlp_upper_bound_inf and nlp_lower_bound_inf
00187     // is 1e19 and can be changed through ipopt options.
00188     // e.g. g_u[0] = 2e19;
00189 
00190     //constraint lower bounds
00191     double * mdConLB = osinstance->getConstraintLowerBounds();
00192     //constraint upper bounds
00193     double * mdConUB = osinstance->getConstraintUpperBounds();
00194 
00195     for(int i = 0; i < m; i++)
00196     {
00197         g_l[ i] = mdConLB[ i];
00198         g_u[ i] = mdConUB[ i];
00199         //cout << "lower !!!!!!!!!!!!!!!!!!!!!!!!!!!" << g_l[i] << endl;
00200         //cout << "upper !!!!!!!!!!!!!!!!!!!!!!!!!!!" << g_u[i] << endl;
00201     }
00202     return true;
00203 }//get_bounds_info
00204 
00205 
00206 // returns the initial point for the problem
00207 bool IpoptProblem::get_starting_point(Index n, bool init_x, Number* x,
00208                                       bool init_z, Number* z_L, Number* z_U, Index m, bool init_lambda,
00209                                       Number* lambda)
00210 {
00211     // Here, we assume we only have starting values for x, if you code
00212     // your own NLP, you can provide starting values for the dual variables
00213     // if you wish
00214     assert(init_x == true);
00215     assert(init_z == false);
00216     assert(init_lambda == false);
00217     int i, m1, n1;
00218 
00219 #ifdef DEBUG
00220     cout << "get initial values !!!!!!!!!!!!!!!!!!!!!!!!!! " << endl;
00221 #endif
00222 
00223     //now set initial values
00224 #ifdef DEBUG
00225     cout << "get number of initial values !!!!!!!!!!!!!!!!!!!!!!!!!! " << endl;
00226 #endif
00227     int k;
00228     if (osoption != NULL)
00229         m1 = osoption->getNumberOfInitVarValues();
00230     else
00231         m1 = 0;
00232 #ifdef DEBUG
00233     cout << "number of variables initialed: " << m1 << endl;
00234 #endif
00235 
00236     n1 = osinstance->getVariableNumber();
00237     bool* initialed;
00238     initialed = new bool[n1];
00239 #ifdef DEBUG
00240     cout << "number of variables in total: " << n1 << endl;
00241 #endif
00242 
00243     for(k = 0; k < n1; k++)
00244         initialed[k] = false;
00245 
00246     if (m1 > 0)
00247     {
00248 #ifdef DEBUG
00249         cout << "get initial values " << endl;
00250 #endif
00251 
00252         InitVarValue**  initVarVector = osoption->getInitVarValuesSparse();
00253 #ifdef DEBUG
00254         cout << "done " << endl;
00255 #endif
00256         try
00257         {
00258             double initval;
00259             for(k = 0; k < m1; k++)
00260             {
00261                 i = initVarVector[k]->idx;
00262                 if (initVarVector[k]->idx > n1)
00263                     throw ErrorClass ("Illegal index value in variable initialization");
00264 
00265                 initval = initVarVector[k]->value;
00266                 if (osinstance->instanceData->variables->var[k]->ub == OSDBL_MAX)
00267                 {
00268                     if (osinstance->instanceData->variables->var[k]->lb > initval)
00269                         throw ErrorClass ("Initial value outside of bounds");
00270                 }
00271                 else if (osinstance->instanceData->variables->var[k]->lb == -OSDBL_MAX)
00272                 {
00273                     if (osinstance->instanceData->variables->var[k]->ub < initval)
00274                         throw ErrorClass ("Initial value outside of bounds");
00275                 }
00276                 else
00277                 {
00278                     if ((osinstance->instanceData->variables->var[k]->lb > initval) ||
00279                             (osinstance->instanceData->variables->var[k]->ub < initval))
00280                         throw ErrorClass ("Initial value outside of bounds");
00281                 }
00282 
00283                 x[initVarVector[k]->idx] = initval;
00284                 initialed[initVarVector[k]->idx] = true;
00285             }
00286         }
00287         catch(const ErrorClass& eclass)
00288         {
00289             cout << "Error in IpoptProblem::get_starting_point (OSIpoptSolver.cpp, line 247)";
00290             cout << endl << endl << endl;
00291         }
00292     }  //  end if (m1 > 0)
00293 
00294     double default_initval;
00295     default_initval = 1.7171;
00296     //default_initval = 0;
00297 
00298     for(k = 0; k < n1; k++)
00299     {
00300         if (!initialed[k])
00301             if (osinstance->instanceData->variables->var[k]->ub == OSDBL_MAX)
00302                 if (osinstance->instanceData->variables->var[k]->lb <= default_initval)
00303                     x[k] = default_initval;
00304                 else
00305                     x[k] = osinstance->instanceData->variables->var[k]->lb;
00306             else if (osinstance->instanceData->variables->var[k]->lb == -OSDBL_MAX)
00307                 if (osinstance->instanceData->variables->var[k]->ub >= default_initval)
00308                     x[k] = default_initval;
00309                 else
00310                     x[k] = osinstance->instanceData->variables->var[k]->ub;
00311             else if ((osinstance->instanceData->variables->var[k]->lb <= default_initval) &&
00312                      (osinstance->instanceData->variables->var[k]->ub >= default_initval))
00313                 x[k] = default_initval;
00314             else if (osinstance->instanceData->variables->var[k]->lb > default_initval)
00315                 x[k] = osinstance->instanceData->variables->var[k]->lb;
00316             else
00317                 x[k] = osinstance->instanceData->variables->var[k]->ub;
00318     }
00319 
00320 #ifdef DEBUG
00321     for(i = 0; i < n1; i++)
00322     {
00323         std::cout << "INITIAL VALUE !!!!!!!!!!!!!!!!!!!!  " << x[ i] << std::endl;
00324     }
00325 #endif
00326     //make sure objvalue is initialized
00327     osinstance->calculateAllObjectiveFunctionValues( x, true);
00328     delete[] initialed;
00329 
00330     return true;
00331 }//get_starting_point
00332 
00333 // returns the value of the objective function
00334 bool IpoptProblem::eval_f(Index n, const Number* x, bool new_x, Number& obj_value)
00335 {
00336     try
00337     {
00338         if(osinstance->getObjectiveNumber() > 0)
00339         {
00340             //the following is a kludge for ipopt, new_x does not seem to get initilized if there are no constraints.
00341             //if(osinstance->getConstraintNumber() <= 0) new_x = true;
00342             if(new_x == false) obj_value  = osinstance->calculateAllObjectiveFunctionValues( const_cast<double*>(x), false)[ 0];
00343             else obj_value = osinstance->calculateAllObjectiveFunctionValues( const_cast<double*>(x), NULL, NULL, true, 0 )[ 0];
00344             //if( CoinIsnan( (double)obj_value) ) return false;
00345             //if( CoinIsnan( obj_value ) )return false;
00346         }
00347 
00348     }
00349     catch(const ErrorClass& eclass)
00350     {
00351 
00352         *ipoptErrorMsg = eclass.errormsg;
00353         throw;
00354     }
00355 
00356     return true;
00357 }
00358 
00359 bool IpoptProblem::eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f)
00360 {
00361     int i;
00362     double *objGrad = NULL;
00363     if(osinstance->getObjectiveNumber() > 0)
00364     {
00365         try
00366         {
00367             //objGrad = osinstance->calculateAllObjectiveFunctionGradients( const_cast<double*>(x), NULL, NULL,  new_x, 1)[ 0];
00368             objGrad = osinstance->calculateObjectiveFunctionGradient( const_cast<double*>(x), NULL, NULL, -1,  new_x, 1);
00369         }
00370         catch(const ErrorClass& eclass)
00371         {
00372 #ifdef DEBUG
00373             cout << "error in OSIpoptSolver, line 314:\n" << eclass.errormsg << endl;
00374 #endif
00375             *ipoptErrorMsg = eclass.errormsg;
00376             throw;
00377         }
00378         for(i = 0; i < n; i++)
00379         {
00380             grad_f[ i]  = objGrad[ i];
00381         }
00382     }
00383     return true;
00384 }//eval_grad_f
00385 
00386 // return the value of the constraints: g(x)
00387 bool IpoptProblem::eval_g(Index n, const Number* x, bool new_x, Index m, Number* g)
00388 {
00389     try
00390     {
00391         double *conVals = osinstance->calculateAllConstraintFunctionValues( const_cast<double*>(x), NULL, NULL, new_x, 0 );
00392         int i;
00393         for(i = 0; i < m; i++)
00394         {
00395             if( CoinIsnan( (double)conVals[ i] ) ) return false;
00396             g[i] = conVals[ i]  ;
00397         }
00398         return true;
00399     }
00400     catch(const ErrorClass& eclass)
00401     {
00402 #ifdef DEBUG
00403         cout << "error in OSIpoptSolver, line 338:\n" << eclass.errormsg << endl;
00404 #endif
00405         *ipoptErrorMsg = eclass.errormsg;
00406         throw;
00407     }
00408 }//eval_g
00409 
00410 
00411 // return the structure or values of the jacobian
00412 bool IpoptProblem::eval_jac_g(Index n, const Number* x, bool new_x,
00413                               Index m, Index nele_jac, Index* iRow, Index *jCol,
00414                               Number* values)
00415 {
00416     SparseJacobianMatrix *sparseJacobian;
00417     if (values == NULL)
00418     {
00419         // return the values of the jacobian of the constraints
00420         //cout << "n: " << n << endl;
00421         //cout << "m: " << m << endl;
00422         //cout << "nele_jac: " <<  nele_jac << endl;
00423         // return the structure of the jacobian
00424         try
00425         {
00426             sparseJacobian = osinstance->getJacobianSparsityPattern();
00427         }
00428         catch(const ErrorClass& eclass)
00429         {
00430 #ifdef DEBUG
00431             cout << "error in OSIpoptSolver, line 362:\n" << eclass.errormsg << endl;
00432 #endif
00433             *ipoptErrorMsg =  eclass.errormsg;
00434             throw;
00435         }
00436         int i = 0;
00437         int k, idx;
00438         for(idx = 0; idx < m; idx++)
00439         {
00440             for(k = *(sparseJacobian->starts + idx); k < *(sparseJacobian->starts + idx + 1); k++)
00441             {
00442                 iRow[i] = idx;
00443                 jCol[i] = *(sparseJacobian->indexes + k);
00444                 //cout << "ROW IDX  !!!!!!!!!!!!!!!!!!!!!!!!!!!"  << iRow[i] << endl;
00445                 //cout << "COL IDX  !!!!!!!!!!!!!!!!!!!!!!!!!!!"  << jCol[i] << endl;
00446                 i++;
00447             }
00448         }
00449     }
00450     else
00451     {
00452         //std::cout << "EVALUATING JACOBIAN" << std::endl;
00453         try
00454         {
00455             sparseJacobian = osinstance->calculateAllConstraintFunctionGradients( const_cast<double*>(x), NULL, NULL,  new_x, 1);
00456         }
00457         catch(const ErrorClass& eclass)
00458         {
00459 #ifdef DEBUG
00460             cout << "error in OSIpoptSolver, line 386:\n" << eclass.errormsg << endl;
00461 #endif
00462             *ipoptErrorMsg = eclass.errormsg;
00463             throw;
00464         }
00465         //osinstance->getIterateResults( (double*)x, 0.0, NULL, -1, new_x,  1);
00466         for(int i = 0; i < nele_jac; i++)
00467         {
00468             values[ i] = sparseJacobian->values[i];
00469             //values[ i] = osinstance->m_mdJacValue[ i];
00470             //cout << "values[i]:!!!!!!!!!!!!  " <<  values[ i] << endl;
00471             //cout << "m_mdJacValue[ i]:!!!!!!!!!!!!  " <<  osinstance->m_mdJacValue[ i] << endl;
00472         }
00473     }
00474     return true;
00475 }//eval_jac_g
00476 
00477 //return the structure or values of the hessian
00478 bool IpoptProblem::eval_h(Index n, const Number* x, bool new_x,
00479                           Number obj_factor, Index m, const Number* lambda,
00480                           bool new_lambda, Index nele_hess, Index* iRow,
00481                           Index* jCol, Number* values)
00482 {
00483 
00485     SparseHessianMatrix *sparseHessian;
00486 
00487     int i;
00488     if (values == NULL)
00489     {
00490         // return the structure. This is a symmetric matrix, fill the lower left triangle only.
00491         //cout << "get structure of HESSIAN !!!!!!!!!!!!!!!!!!!!!!!!!! "  << endl;
00492         try
00493         {
00494             sparseHessian = osinstance->getLagrangianHessianSparsityPattern( );
00495         }
00496         catch(const ErrorClass& eclass)
00497         {
00498 
00499             *ipoptErrorMsg = eclass.errormsg;
00500             throw;
00501         }
00502 //              cout << "got structure of HESSIAN !!!!!!!!!!!!!!!!!!!!!!!!!! "  << endl;
00503         for(i = 0; i < nele_hess; i++)
00504         {
00505             iRow[i] = *(sparseHessian->hessColIdx + i);
00506             jCol[i] = *(sparseHessian->hessRowIdx + i);
00507             //cout << "ROW HESS IDX  !!!!!!!!!!!!!!!!!!!!!!!!!!!"  << iRow[i] << endl;
00508             //cout << "COL HESS IDX  !!!!!!!!!!!!!!!!!!!!!!!!!!!"  << jCol[i] << endl;
00509         }
00510     }
00511     else
00512     {
00513         //std::cout << "EVALUATING HESSIAN" << std::endl;
00514         // return the values. This is a symmetric matrix, fill the lower left triangle only
00515         double* objMultipliers = new double[1];
00516         objMultipliers[0] = obj_factor;
00517         try
00518         {
00519             sparseHessian = osinstance->calculateLagrangianHessian( const_cast<double*>(x), objMultipliers, const_cast<double*>(lambda) ,  new_x, 2);
00520             delete[]  objMultipliers;
00521         }
00522         catch(const ErrorClass& eclass)
00523         {
00524 #ifdef DEBUG
00525             cout << "error in OSIpoptSolver, line 444:\n" << eclass.errormsg << endl;
00526 #endif
00527             *ipoptErrorMsg = eclass.errormsg;
00528             delete[]  objMultipliers;
00529             throw;
00530         }
00531         for(i = 0; i < nele_hess; i++)
00532         {
00533             values[ i]  = *(sparseHessian->hessValues + i);
00534         }
00535     }
00537     return true;
00538 }//eval_h
00539 
00540 bool IpoptProblem::get_scaling_parameters(Number& obj_scaling,
00541         bool& use_x_scaling, Index n,
00542         Number* x_scaling,
00543         bool& use_g_scaling, Index m,
00544         Number* g_scaling)
00545 {
00546     if( osinstance->instanceData->objectives->obj[ 0]->maxOrMin.compare("min") != 0)
00547     {
00548         obj_scaling = -1;
00549     }
00550     else obj_scaling = 1;
00551     use_x_scaling = false;
00552     use_g_scaling = false;
00553     return true;
00554 }//get_scaling_parameters
00555 
00556 void IpoptProblem::finalize_solution(SolverReturn status,
00557                                      Index n, const Number* x, const Number* z_L, const Number* z_U,
00558                                      Index m, const Number* g, const Number* lambda, Number obj_value,
00559                                      const IpoptData* ip_data, IpoptCalculatedQuantities* ip_cq)
00560 {
00561     // here is where we would store the solution to variables, or write to a file, etc
00562     // so we could use the solution.
00563     // For this example, we write the solution to the console
00564     if(osinstance->getObjectiveNumber() > 0)
00565         obj_value = osinstance->calculateAllObjectiveFunctionValues( const_cast<double*>(x), true)[ 0];
00566     OSrLWriter *osrlwriter ;
00567     osrlwriter = new OSrLWriter();
00568 #ifdef DEBUG
00569     printf("\n\nSolution of the primal variables, x\n");
00570     for (Index i=0; i<n; i++)
00571     {
00572         printf("x[%d] = %e\n", i, x[i]);
00573     }
00574 
00575     printf("\n\nSolution of the bound multipliers, z_L and z_U\n");
00576     for (Index i=0; i<n; i++)
00577     {
00578         printf("z_L[%d] = %e\n", i, z_L[i]);
00579     }
00580     for (Index i=0; i<n; i++)
00581     {
00582         printf("z_U[%d] = %e\n", i, z_U[i]);
00583     }
00584 #endif
00585     if(osinstance->getObjectiveNumber() > 0) printf("\nObjective value f(x*) = %e\n", obj_value  );
00586 
00587     int solIdx = 0;
00588     int numberOfOtherVariableResults;
00589     int otherIdx;
00590     int numCon = osinstance->getConstraintNumber();
00591     ostringstream outStr;
00592 
00593 
00594     //make sure the sign on the dual values is correct
00595     double *dualValue = NULL;
00596     std::string *rcost = NULL;
00597     int* idx = NULL;
00598     double*  mdObjValues = NULL;
00599     if(osinstance->getObjectiveNumber() > 0)
00600     {
00601         mdObjValues = new double[1];
00602     }
00603 
00604     std::string message = "Ipopt solver finishes to the end.";
00605     std::string solutionDescription = "";
00606 
00607     try
00608     {
00609         // resultHeader infomration
00610         if(osresult->setSolverInvoked( "COIN-OR Ipopt") != true)
00611             throw ErrorClass("OSResult error: setSolverInvoked");
00612         if(osresult->setServiceName( OSgetVersionInfo()) != true)
00613             throw ErrorClass("OSResult error: setServiceName");
00614         if(osresult->setInstanceName(  osinstance->getInstanceName()) != true)
00615             throw ErrorClass("OSResult error: setInstanceName");
00616 
00617         //if(osresult->setJobID( osoption->jobID) != true)
00618         //      throw ErrorClass("OSResult error: setJobID");
00619 
00620         // set basic problem parameters
00621         if(osresult->setVariableNumber( osinstance->getVariableNumber()) != true)
00622             throw ErrorClass("OSResult error: setVariableNumer");
00623         if(osresult->setObjectiveNumber( 1) != true)
00624             throw ErrorClass("OSResult error: setObjectiveNumber");
00625         if(osresult->setConstraintNumber( osinstance->getConstraintNumber()) != true)
00626             throw ErrorClass("OSResult error: setConstraintNumber");
00627         if(osresult->setSolutionNumber(  1) != true)
00628             throw ErrorClass("OSResult error: setSolutionNumer");
00629 
00630 
00631         if(osresult->setGeneralMessage( message) != true)
00632             throw ErrorClass("OSResult error: setGeneralMessage");
00633 
00634         switch( status)
00635         {
00636         case SUCCESS:
00637             solutionDescription = "SUCCESS[IPOPT]: Algorithm terminated normally at a locally optimal point, satisfying the convergence tolerances.";
00638             osresult->setSolutionStatus(solIdx,  "locallyOptimal", solutionDescription);
00639 
00640 
00641             if(osinstance->getVariableNumber() > 0) osresult->setPrimalVariableValuesDense(solIdx, const_cast<double*>(x));
00642 
00643             if(numCon > 0)
00644             {
00645                 dualValue = new double[ numCon];
00646                 for (Index i=0; i < numCon; i++)
00647                 {
00648                     dualValue[ i] = -lambda[ i];
00649                 }
00650                 //osresult->setDualVariableValuesDense(solIdx, const_cast<double*>(lambda));
00651                 osresult->setDualVariableValuesDense(solIdx, dualValue);
00652             }
00653 
00654             if(osinstance->getObjectiveNumber() > 0)
00655             {
00656                 mdObjValues[0] = obj_value ;
00657                 osresult->setObjectiveValuesDense(solIdx, mdObjValues);
00658             }
00659             // set other
00660 
00661             if(osinstance->getVariableNumber() > 0)
00662             {
00663                 numberOfOtherVariableResults = 1;
00664                 osresult->setNumberOfOtherVariableResults(solIdx, numberOfOtherVariableResults);
00665                 rcost = new std::string[ osinstance->getVariableNumber()];
00666                 idx = new int[ osinstance->getVariableNumber()];
00667                 for (Index i = 0; i < n; i++)
00668                 {
00669                     rcost[ i] =  os_dtoa_format( z_L[i] - z_U[i]);
00670                     idx[ i] = i;
00671                 }
00672                 otherIdx = 0;
00673                 osresult->setAnOtherVariableResultSparse(solIdx, otherIdx, "reduced costs", "", "the variable reduced costs",
00674                         idx, rcost, osinstance->getVariableNumber());
00675             }
00676             //set dual values on variable upper and lower bounds
00677             /*
00678             for (Index i = 0; i < n; i++) {
00679                 rcost[ i] =  os_dtoa_format( z_L[i]);
00680                 idx[ i] = i;
00681              }
00682             otherIdx = 0;
00683             osresult->setAnOtherVariableResultSparse(solIdx, otherIdx, "varL", "", "Lagrange Multiplier on the Variable Lower Bound",  idx,  rcost,  osinstance->getVariableNumber() );
00684 
00685 
00686             for (Index i = 0; i < n; i++) {
00687                 rcost[ i] =  os_dtoa_format( z_U[i]);
00688              }
00689             otherIdx = 1;
00690             osresult->setAnOtherVariableResultSparse(solIdx, otherIdx, "varU", "", "Lagrange Multiplier on the Variable Upper Bound", idx,  rcost, osinstance->getVariableNumber() );
00691             */
00692             if(osinstance->getVariableNumber() > 0) delete[] rcost;
00693             if(osinstance->getVariableNumber() > 0) delete[] idx;
00694             if(osinstance->getConstraintNumber() > 0) delete[] dualValue;
00695 
00696             // done with other
00697 
00698 
00699             break;
00700         case MAXITER_EXCEEDED:
00701             solutionDescription = "MAXITER_EXCEEDED[IPOPT]: Maximum number of iterations exceeded.";
00702             osresult->setSolutionStatus(solIdx,  "stoppedByLimit", solutionDescription);
00703             if(x != NULL) osresult->setPrimalVariableValuesDense(solIdx, const_cast<double*>(x));
00704             if(lambda != NULL) osresult->setDualVariableValuesDense(solIdx, const_cast<double*>( lambda));
00705             if(osinstance->getObjectiveNumber() > 0)
00706             {
00707                 mdObjValues[0] = obj_value ;
00708                 osresult->setObjectiveValuesDense(solIdx, mdObjValues);
00709             }
00710             break;
00711         case STOP_AT_TINY_STEP:
00712             solutionDescription = "STOP_AT_TINY_STEP[IPOPT]: Algorithm proceeds with very little progress.";
00713             osresult->setSolutionStatus(solIdx,  "stoppedByLimit", solutionDescription);
00714             if(x != NULL) osresult->setPrimalVariableValuesDense(solIdx, const_cast<double*>( x));
00715             if(lambda != NULL) osresult->setDualVariableValuesDense(solIdx, const_cast<double*>( lambda));
00716             if(osinstance->getObjectiveNumber() > 0)
00717             {
00718                 mdObjValues[0] = obj_value ;
00719                 osresult->setObjectiveValuesDense(solIdx, mdObjValues);
00720             }
00721             break;
00722         case STOP_AT_ACCEPTABLE_POINT:
00723             solutionDescription = "STOP_AT_ACCEPTABLE_POINT[IPOPT]: Algorithm stopped at a point that was converged, not to _desired_ tolerances, but to _acceptable_ tolerances";
00724             osresult->setSolutionStatus(solIdx,  "IpoptAccetable", solutionDescription);
00725             if(lambda != NULL) osresult->setDualVariableValuesDense(solIdx, const_cast<double*>( lambda));
00726             if(x != NULL)osresult->setPrimalVariableValuesDense(solIdx, const_cast<double*>(x));
00727             if(osinstance->getObjectiveNumber() > 0)
00728             {
00729                 mdObjValues[0] = obj_value ;
00730                 osresult->setObjectiveValuesDense(solIdx, mdObjValues);
00731             }
00732             break;
00733         case LOCAL_INFEASIBILITY:
00734             solutionDescription = "LOCAL_INFEASIBILITY[IPOPT]: Algorithm converged to a point of local infeasibility. Problem may be infeasible.";
00735             osresult->setSolutionStatus(solIdx,  "infeasible", solutionDescription);
00736             if( osinstance->getVariableNumber() == 0) osresult->setSolutionMessage(solIdx, "Warning: this problem has zero decision variables!");
00737             break;
00738         case USER_REQUESTED_STOP:
00739             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.";
00740             osresult->setSolutionStatus(solIdx,  "error", solutionDescription);
00741             if( osinstance->getVariableNumber() == 0) osresult->setSolutionMessage(solIdx, "Warning: this problem has zero decision variables!");
00742             break;
00743         case DIVERGING_ITERATES:
00744             solutionDescription = "DIVERGING_ITERATES[IPOPT]: It seems that the iterates diverge.";
00745             osresult->setSolutionStatus(solIdx,  "unbounded", solutionDescription);
00746             if( osinstance->getVariableNumber() == 0) osresult->setSolutionMessage(solIdx, "Warning: this problem has zero decision variables!");
00747             break;
00748         case RESTORATION_FAILURE:
00749             solutionDescription = "RESTORATION_FAILURE[IPOPT]: Restoration phase failed, algorithm doesn't know how to proceed.";
00750             osresult->setSolutionStatus(solIdx,  "error", solutionDescription);
00751             if( osinstance->getVariableNumber() == 0) osresult->setSolutionMessage(solIdx, "Warning: this problem has zero decision variables!");
00752             break;
00753         case ERROR_IN_STEP_COMPUTATION:
00754             solutionDescription = "ERROR_IN_STEP_COMPUTATION[IPOPT]: An unrecoverable error occurred while IPOPT tried to compute the search direction.";
00755             osresult->setSolutionStatus(solIdx,  "error", solutionDescription);
00756             if( osinstance->getVariableNumber() == 0) osresult->setSolutionMessage(solIdx, "Warning: this problem has zero decision variables!");
00757             break;
00758         case INVALID_NUMBER_DETECTED:
00759             solutionDescription = "INVALID_NUMBER_DETECTED[IPOPT]: Algorithm received an invalid number (such as NaN or Inf) from the NLP; see also option check_derivatives_for_naninf.";
00760             osresult->setSolutionStatus(solIdx,  "error", solutionDescription);
00761             if( osinstance->getVariableNumber() == 0) osresult->setSolutionMessage(solIdx, "Warning: this problem has zero decision variables!");
00762             break;
00763         case INTERNAL_ERROR:
00764             solutionDescription = "INTERNAL_ERROR[IPOPT]: An unknown internal error occurred. Please contact the IPOPT authors through the mailing list.";
00765             osresult->setSolutionStatus(solIdx,  "error", solutionDescription);
00766             if( osinstance->getVariableNumber() == 0) osresult->setSolutionMessage(solIdx, "Warning: this problem has zero decision variables!");
00767             break;
00768         default:
00769             solutionDescription = "OTHER[IPOPT]: other unknown solution status from Ipopt solver";
00770             osresult->setSolutionStatus(solIdx,  "other", solutionDescription);
00771             if( osinstance->getVariableNumber() == 0) osresult->setSolutionMessage(solIdx, "Warning: this problem has zero decision variables!");
00772         }
00773         osresult->setGeneralStatusType("normal");
00774         delete osrlwriter;
00775         if(osinstance->getObjectiveNumber() > 0)
00776         {
00777             delete[] mdObjValues;
00778         }
00779         osrlwriter = NULL;
00780 
00781     }
00782     catch(const ErrorClass& eclass)
00783     {
00784 #ifdef DEBUG
00785         cout << "error in OSIpoptSolver, line 636:\n" << eclass.errormsg << endl;
00786 #endif
00787         osresult->setGeneralMessage( eclass.errormsg);
00788         osresult->setGeneralStatusType( "error");
00789         std::string osrl = osrlwriter->writeOSrL( osresult);
00790         delete osrlwriter;
00791         osrlwriter = NULL;
00792         throw ErrorClass(  osrl) ;
00793         if(osinstance->getObjectiveNumber() > 0)
00794         {
00795             delete[] mdObjValues;
00796         }
00797         mdObjValues = NULL;
00798     }
00800 }
00801 
00802 
00803 void IpoptSolver::setSolverOptions() throw (ErrorClass)
00804 {
00805     try
00806     {
00807         //if(osinstance->getObjectiveNumber() <= 0) throw ErrorClass("Ipopt NEEDS AN OBJECTIVE FUNCTION");
00808         this->bSetSolverOptions = true;
00809         /* set the default options */
00810         //app->Options()->SetNumericValue("tol", 1e-9);
00811         app->Options()->SetIntegerValue("print_level", 0);
00812         app->Options()->SetIntegerValue("max_iter", 20000);
00813         app->Options()->SetNumericValue("bound_relax_factor", 0, true, true);
00814         app->Options()->SetStringValue("mu_strategy", "adaptive", true, true);
00815         //app->Options()->SetStringValue("output_file", "ipopt.out");
00816         app->Options()->SetStringValue("check_derivatives_for_naninf", "yes");
00817         // hessian constant for an LP
00818         if( (osinstance-> getNumberOfNonlinearExpressions() <= 0) && (osinstance->getNumberOfQuadraticTerms() <= 0) )
00819         {
00820             app->Options()->SetStringValue("hessian_constant", "yes", true, true);
00821         }
00822         if(osinstance->getObjectiveNumber() > 0)
00823         {
00824             if( osinstance->instanceData->objectives->obj[ 0]->maxOrMin.compare("min") != 0)
00825             {
00826                 app->Options()->SetStringValue("nlp_scaling_method", "user-scaling");
00827             }
00828         }
00829         /* end of the default options, now get options from OSoL */
00830 
00831 
00832         if(osoption == NULL && osol.length() > 0)
00833         {
00834             m_osolreader = new OSoLReader();
00835             osoption = m_osolreader->readOSoL( osol);
00836         }
00837 
00838         if( osoption != NULL  &&  osoption->getNumberOfSolverOptions() > 0 )
00839         {
00840 //                      std::cout << "number of solver options "  <<  osoption->getNumberOfSolverOptions() << std::endl;
00841             std::vector<SolverOption*> optionsVector;
00842             optionsVector = osoption->getSolverOptions( "ipopt",true);
00843             char *pEnd;
00844             int i;
00845             int num_ipopt_options = optionsVector.size();
00846             for(i = 0; i < num_ipopt_options; i++)
00847             {
00848 //                              std::cout << "ipopt solver option  "  << optionsVector[ i]->name << std::endl;
00849                 if(optionsVector[ i]->type == "numeric" )
00850                 {
00851 //                                      std::cout << "FOUND A NUMERIC OPTION  "  <<  os_strtod( optionsVector[ i]->value.c_str(), &pEnd ) << std::endl;
00852                     app->Options()->SetNumericValue(optionsVector[ i]->name, os_strtod( optionsVector[ i]->value.c_str(), &pEnd ) );
00853                 }
00854                 else if(optionsVector[ i]->type == "integer" )
00855                 {
00856 //                                      std::cout << "FOUND AN INTEGER OPTION  "  << atoi( optionsVector[ i]->value.c_str() ) << std::endl;
00857                     app->Options()->SetIntegerValue(optionsVector[ i]->name, atoi( optionsVector[ i]->value.c_str() ) );
00858                 }
00859                 else if(optionsVector[ i]->type == "string" )
00860                 {
00861 //                                      std::cout << "FOUND A STRING OPTION  "  <<  optionsVector[ i]->value.c_str() << std::endl;
00862                     app->Options()->SetStringValue(optionsVector[ i]->name, optionsVector[ i]->value);
00863                 }
00864             }
00865         }
00866     }
00867     catch(const ErrorClass& eclass)
00868     {
00869         std::cout << "THERE IS AN ERROR" << std::endl;
00870         osresult->setGeneralMessage( eclass.errormsg);
00871         osresult->setGeneralStatusType( "error");
00872         osrl = osrlwriter->writeOSrL( osresult);
00873         throw ErrorClass( osrl) ;
00874     }
00875 }//end setSolverOptions()
00876 
00877 
00878 
00879 void IpoptSolver::buildSolverInstance() throw (ErrorClass)
00880 {
00881     try
00882     {
00883 
00884         if(osil.length() == 0 && osinstance == NULL) throw ErrorClass("there is no instance");
00885         if(osinstance == NULL)
00886         {
00887             m_osilreader = new OSiLReader();
00888             osinstance = m_osilreader->readOSiL( osil);
00889         }
00890         // Create a new instance of your nlp
00891         nlp = new IpoptProblem( osinstance, osoption, osresult, ipoptErrorMsg);
00892         app = new IpoptApplication();
00893         this->bCallbuildSolverInstance = true;
00894     }
00895     catch(const ErrorClass& eclass)
00896     {
00897 #ifdef DEBUG
00898         cout << "error in OSIpoptSolver, line 722:\n" << eclass.errormsg << endl;
00899 #endif
00900         std::cout << "THERE IS AN ERROR" << std::endl;
00901         osresult->setGeneralMessage( eclass.errormsg);
00902         osresult->setGeneralStatusType( "error");
00903         osrl = osrlwriter->writeOSrL( osresult);
00904         throw ErrorClass( osrl) ;
00905     }
00906 }//end buildSolverInstance()
00907 
00908 
00909 
00910 void IpoptSolver::solve() throw (ErrorClass)
00911 {
00912     if( this->bCallbuildSolverInstance == false) buildSolverInstance();
00913     if( this->bSetSolverOptions == false) setSolverOptions();
00914     try
00915     {
00916         //clock_t start, finish;
00917         //double duration;
00918         //start = clock();
00919         //OSiLWriter osilwriter;
00920         //cout << osilwriter.writeOSiL( osinstance) << endl;
00921         //if(osinstance->getVariableNumber() <= 0)throw ErrorClass("Ipopt requires decision variables");
00922         //finish = clock();
00923         //duration = (double) (finish - start) / CLOCKS_PER_SEC;
00924         //cout << "Parsing took (seconds): " << duration << endl;
00925         //dataEchoCheck();
00926         /***************now the ipopt invokation*********************/
00927         // see if we have a linear program
00928         //if(osinstance->getObjectiveNumber() <= 0) throw ErrorClass("Ipopt NEEDS AN OBJECTIVE FUNCTION");
00929         // Intialize the IpoptApplication and process the options
00930 //              std::cout << "Call Ipopt Initialize" << std::endl;
00931         app->Initialize();
00932 //              std::cout << "Finished Ipopt Initialize" << std::endl;
00933         //nlp->osinstance = this->osinstance;
00934         // Ask Ipopt to solve the problem
00935 //              std::cout << "Call Ipopt Optimize" << std::endl;
00936         ApplicationReturnStatus status = app->OptimizeTNLP( nlp);
00937 //              std::cout << "Finish Ipopt Optimize" << std::endl;
00938         osrl = osrlwriter->writeOSrL( osresult);
00939 //              std::cout << "Finish writing the osrl" << std::endl;
00940         //if (status != Solve_Succeeded) {
00941         if (status < -2)
00942         {
00943             throw ErrorClass("Ipopt FAILED TO SOLVE THE PROBLEM: " + *ipoptErrorMsg);
00944         }
00945     }
00946     catch(const ErrorClass& eclass)
00947     {
00948 #ifdef DEBUG
00949         cout << "error in OSIpoptSolver, line 775:\n" << eclass.errormsg << endl;
00950 #endif
00951         osresult->setGeneralMessage( eclass.errormsg);
00952         osresult->setGeneralStatusType( "error");
00953         osrl = osrlwriter->writeOSrL( osresult);
00954         throw ErrorClass( osrl) ;
00955     }
00956 }//solve
00957 
00958 
00959 void IpoptSolver::dataEchoCheck()
00960 {
00961 
00962     int i;
00963 
00964     // print out problem parameters
00965     cout << "This is problem:  " << osinstance->getInstanceName() << endl;
00966     cout << "The problem source is:  " << osinstance->getInstanceSource() << endl;
00967     cout << "The problem description is:  " << osinstance->getInstanceDescription() << endl;
00968     cout << "number of variables = " << osinstance->getVariableNumber() << endl;
00969     cout << "number of Rows = " << osinstance->getConstraintNumber() << endl;
00970 
00971     // print out the variable information
00972     if(osinstance->getVariableNumber() > 0)
00973     {
00974         for(i = 0; i < osinstance->getVariableNumber(); i++)
00975         {
00976             if(osinstance->getVariableNames() != NULL) cout << "variable Names  " << osinstance->getVariableNames()[ i]  << endl;
00977             if(osinstance->getVariableTypes() != NULL) cout << "variable Types  " << osinstance->getVariableTypes()[ i]  << endl;
00978             if(osinstance->getVariableLowerBounds() != NULL) cout << "variable Lower Bounds  " << osinstance->getVariableLowerBounds()[ i]  << endl;
00979             if(osinstance->getVariableUpperBounds() != NULL) cout << "variable Upper Bounds  " <<  osinstance->getVariableUpperBounds()[i] << endl;
00980         }
00981     }
00982 
00983     // print out objective function information
00984     if(osinstance->getVariableNumber() > 0 || osinstance->instanceData->objectives->obj != NULL || osinstance->instanceData->objectives->numberOfObjectives > 0)
00985     {
00986         if( osinstance->getObjectiveMaxOrMins()[0] == "min")  cout <<  "problem is a minimization" << endl;
00987         else cout <<  "problem is a maximization" << endl;
00988         for(i = 0; i < osinstance->getVariableNumber(); i++)
00989         {
00990             cout << "OBJ COEFFICIENT =  " <<  osinstance->getDenseObjectiveCoefficients()[0][i] << endl;
00991         }
00992     }
00993     // print out constraint information
00994     if(osinstance->getConstraintNumber() > 0)
00995     {
00996         for(i = 0; i < osinstance->getConstraintNumber(); i++)
00997         {
00998             if(osinstance->getConstraintNames() != NULL) cout << "row name = " << osinstance->getConstraintNames()[i] <<  endl;
00999             if(osinstance->getConstraintLowerBounds() != NULL) cout << "row lower bound = " << osinstance->getConstraintLowerBounds()[i] <<  endl;
01000             if(osinstance->getConstraintUpperBounds() != NULL) cout << "row upper bound = " << osinstance->getConstraintUpperBounds()[i] <<  endl;
01001         }
01002     }
01003 
01004     // print out linear constraint data
01005     cout << endl;
01006     cout << "number of nonzeros =  " << osinstance->getLinearConstraintCoefficientNumber() << endl;
01007     for(i = 0; i <= osinstance->getVariableNumber(); i++)
01008     {
01009         cout << "Start Value =  " << osinstance->getLinearConstraintCoefficientsInColumnMajor()->starts[ i] << endl;
01010     }
01011     cout << endl;
01012     for(i = 0; i < osinstance->getLinearConstraintCoefficientNumber(); i++)
01013     {
01014         cout << "Index Value =  " << osinstance->getLinearConstraintCoefficientsInColumnMajor()->indexes[i] << endl;
01015         cout << "Nonzero Value =  " << osinstance->getLinearConstraintCoefficientsInColumnMajor()->values[i] << endl;
01016     }
01017 
01018     // print out quadratic data
01019     cout << "number of qterms =  " <<  osinstance->getNumberOfQuadraticTerms() << endl;
01020     for(int i = 0; i <  osinstance->getNumberOfQuadraticTerms(); i++)
01021     {
01022         cout << "Row Index = " <<  osinstance->getQuadraticTerms()->rowIndexes[i] << endl;
01023         cout << "Var Index 1 = " << osinstance->getQuadraticTerms()->varOneIndexes[ i] << endl;
01024         cout << "Var Index 2 = " << osinstance->getQuadraticTerms()->varTwoIndexes[ i] << endl;
01025         cout << "Coefficient = " << osinstance->getQuadraticTerms()->coefficients[ i] << endl;
01026     }
01027 } // end dataEchoCheck
01028 
01029 
01030 IpoptProblem::IpoptProblem(OSInstance *osinstance_,  OSOption *osoption_, OSResult *osresult_, std::string* ipoptErrorMsg_)
01031 {
01032     osinstance = osinstance_;
01033     osoption = osoption_;
01034     osresult = osresult_;
01035     ipoptErrorMsg = ipoptErrorMsg_;
01036 }
01037 
01038 IpoptProblem::~IpoptProblem()
01039 {
01040 
01041 }
01042 
01043 
01044 

Generated on Wed Nov 30 03:04:24 2011 by  doxygen 1.4.7