00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "BonOuterDescription.hpp"
00012 #include "BonOsiTMINLPInterface.hpp"
00013
00014 namespace Bonmin{
00015
00016
00017
00018
00019
00020 static inline
00021 bool cleanNnz(double &value, double colLower, double colUpper,
00022 double rowLower, double rowUpper, double colsol,
00023 double & lb, double &ub, double tiny, double veryTiny)
00024 {
00025 if(fabs(value)>= tiny) return 1;
00026
00027 if(fabs(value)<veryTiny) return 0;
00028
00029
00030 double infty = 1e20;
00031 bool colUpBounded = colUpper < 10000;
00032 bool colLoBounded = colLower > -10000;
00033 bool rowNotLoBounded = rowLower <= - infty;
00034 bool rowNotUpBounded = rowUpper >= infty;
00035 bool pos = value > 0;
00036
00037 if(colLoBounded && pos && rowNotUpBounded) {
00038 lb += value * (colsol - colLower);
00039 return 0;
00040 }
00041 else
00042 if(colLoBounded && !pos && rowNotLoBounded) {
00043 ub += value * (colsol - colLower);
00044 return 0;
00045 }
00046 else
00047 if(colUpBounded && !pos && rowNotUpBounded) {
00048 lb += value * (colsol - colUpper);
00049 return 0;
00050 }
00051 else
00052 if(colUpBounded && pos && rowNotLoBounded) {
00053 ub += value * (colsol - colUpper);
00054 return 0;
00055 }
00056
00057 if(pos) value = tiny;
00058 else
00059 value = - tiny;
00060 return 1;
00061 }
00062
00063
00067 void getMyOuterApproximation(
00068 OsiTMINLPInterface &si, OsiCuts &cs, int ind,
00069 const double * x, int getObj, const double * x2, double theta,
00070 bool global) {
00071 int n, m, nnz_jac_g, nnz_h_lag;
00072 Ipopt::TNLP::IndexStyleEnum index_style;
00073 TMINLP2TNLP* problem = si.problem();
00074 problem->get_nlp_info(n, m, nnz_jac_g, nnz_h_lag, index_style);
00075
00076 double g_i = 0;
00077 problem->eval_gi(n, x, 1, ind, g_i);
00078 vector<int> jCol(n);
00079 int nnz;
00080 problem->eval_grad_gi(n, x, 0, ind, nnz, jCol(), NULL);
00081 vector<double> jValues(nnz);
00082 problem->eval_grad_gi(n, x, 0, ind, nnz, NULL, jValues());
00083
00084 CoinPackedVector cut;
00085 double lb;
00086 double ub;
00087
00088 const double * rowLower = si.getRowLower();
00089 const double * rowUpper = si.getRowUpper();
00090 const double * colLower = si.getColLower();
00091 const double * colUpper = si.getColUpper();
00092 const double * duals = si.getRowPrice() + 2 * n;
00093 double infty = si.getInfinity();
00094 double nlp_infty = infty;
00095
00096 int rowIdx = ind;
00097
00098 if (rowLower[rowIdx] > -nlp_infty)
00099 lb = rowLower[rowIdx] - g_i;
00100 else
00101 lb = -infty;
00102 if (rowUpper[rowIdx] < nlp_infty)
00103 ub = rowUpper[rowIdx] - g_i;
00104 else
00105 ub = infty;
00106 if (rowLower[rowIdx] > -infty && rowUpper[rowIdx] < infty) {
00107 if (duals[rowIdx] >= 0)
00108 lb = -infty;
00109 if (duals[rowIdx] <= 0)
00110 ub = infty;
00111 }
00112
00113 double tiny = 1e-08;
00114 double veryTiny = 1e-20;
00115
00116 for (int i = 0; i < nnz; i++) {
00117 if(index_style == Ipopt::TNLP::FORTRAN_STYLE) jCol[i]--;
00118 const int &colIdx = jCol[i];
00119
00120 if (cleanNnz(jValues[i], colLower[colIdx], colUpper[colIdx],
00121 rowLower[rowIdx], rowUpper[rowIdx], x[colIdx], lb, ub,
00122 tiny, veryTiny)) {
00123 cut.insert(colIdx, jValues[i]);
00124 if (lb > -infty)
00125 lb += jValues[i] * x[colIdx];
00126 if (ub < infty)
00127 ub += jValues[i] * x[colIdx];
00128 }
00129 }
00130
00131
00132 bool add = true;
00133
00134 if (x2 != NULL) {
00135 double rhs = cut.dotProduct(x2);
00136 double violation = 0.;
00137 if (ub < infty)
00138 violation = std::max(violation, fabs(rhs - ub));
00139 if (lb > -infty)
00140 violation = std::max(violation, fabs(lb - rhs));
00141 if (violation < theta) {
00142 add = false;
00143 }
00144 }
00145 OsiRowCut newCut;
00146
00147
00148 if (add) {
00149 if (global) {
00150 newCut.setGloballyValidAsInteger(1);
00151 }
00152
00153 newCut.setLb(lb);
00154 newCut.setUb(ub);
00155 newCut.setRow(cut);
00156
00157 cs.insert(newCut);
00158 }
00159 }
00160
00161 void addOuterDescription(OsiTMINLPInterface &nlp, OsiSolverInterface &si,
00162 const double * x, int nbAp, bool getObj) {
00163 int n;
00164 int m;
00165 int nnz_jac_g;
00166 int nnz_h_lag;
00167 Ipopt::TNLP::IndexStyleEnum index_style;
00168
00169 TMINLP2TNLP* problem = nlp.problem();
00170 problem->get_nlp_info(n, m, nnz_jac_g, nnz_h_lag, index_style);
00171
00172 const double * colLower = nlp.getColLower();
00173 const double * colUpper = nlp.getColUpper();
00174 const Bonmin::TMINLP::VariableType* variableType = problem->var_types();
00175 vector<Ipopt::TNLP::LinearityType> constTypes(m);
00176 problem->get_constraints_linearity(m, constTypes());
00177
00178 int OuterDesc = 0;
00179 if (OuterDesc == 0) {
00180 OsiCuts cs;
00181
00182 double * p = CoinCopyOfArray(nlp.getColLower(), n);
00183 double * pp = CoinCopyOfArray(nlp.getColLower(), n);
00184 double * up = CoinCopyOfArray(nlp.getColUpper(), n);
00185
00186 std::vector<int> nbG(m, 0);
00187
00188 std::vector<double> step(n);
00189
00190 for (int i = 0; i < n; i++) {
00191
00192 if (colUpper[i] > 1e08) {
00193 up[i] = 0;
00194 }
00195
00196 if (colUpper[i] > 1e08 || colLower[i] < -1e08 || (variableType[i]
00197 == TMINLP::BINARY) || (variableType[i] == TMINLP::INTEGER)) {
00198 step[i] = 0;
00199 } else
00200 step[i] = (up[i] - colLower[i]) / (nbAp);
00201
00202 if (colLower[i] < -1e08) {
00203 p[i] = 0;
00204 pp[i] = 0;
00205 }
00206 }
00207 vector<double> g_p(m);
00208 vector<double> g_pp(m);
00209 for (int i = 0; (i < m); i++) {
00210 if(constTypes[i] != Ipopt::TNLP::NON_LINEAR) continue;
00211 getMyOuterApproximation(nlp, cs, i, p, 0, NULL, 10000, true);
00212 }
00213 for (int j = 1; j <= nbAp; j++) {
00214
00215 for (int i = 0; i < n; i++) {
00216 pp[i] += step[i];
00217 }
00218
00219
00220 problem->eval_g(n, p, 1, m, g_p());
00221 problem->eval_g(n, pp, 1, m, g_pp());
00222 double diff = 0;
00223 int varInd = 0;
00224 for (int i = 0; (i < m); i++) {
00225 if(constTypes[i] != Ipopt::TNLP::NON_LINEAR) continue;
00226 if (varInd == n - 1)
00227 varInd = 0;
00228 diff = std::abs(g_p[i] - g_pp[i]);
00229
00230 if (nbG[i] < nbAp && diff ) {
00231 getMyOuterApproximation(nlp, cs, i, pp, 0, NULL, 10000, true);
00232 p[varInd] = pp[varInd];
00233 nbG[i]++;
00234 }
00235 varInd++;
00236 }
00237 }
00238 for (int i = 0; i < m ; i++) {
00239 if(constTypes[i] != Ipopt::TNLP::NON_LINEAR) continue;
00240 getMyOuterApproximation(nlp, cs, i, up, 0, NULL, 10000, true);
00241 }
00242
00243 si.applyCuts(cs);
00244 delete [] p;
00245 delete [] pp;
00246 delete [] up;
00247 }
00248
00249 }
00250
00251 }
00252