00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "CouenneCutGenerator.hpp"
00012
00013 #include "CouenneTypes.hpp"
00014 #include "CouenneExprMul.hpp"
00015 #include "CouenneExprTrilinear.hpp"
00016 #include "CouenneProblem.hpp"
00017 #include "CouenneExprAux.hpp"
00018
00019 #include <vector>
00020
00021
00022 using namespace Couenne;
00023
00024 #define EPSILONT 1.e-6
00025
00026
00027
00029 void permutation3(int **ind,int *ibnd)
00030 {
00031 ind[0][0] = ibnd[0]; ind[0][1] = ibnd[1]; ind[0][2] = ibnd[2];
00032 ind[1][0] = ibnd[0]; ind[1][1] = ibnd[2]; ind[1][2] = ibnd[1];
00033 ind[2][0] = ibnd[1]; ind[2][1] = ibnd[0]; ind[2][2] = ibnd[2];
00034 ind[3][0] = ibnd[1]; ind[3][1] = ibnd[2]; ind[3][2] = ibnd[0];
00035 ind[4][0] = ibnd[2]; ind[4][1] = ibnd[0]; ind[4][2] = ibnd[1];
00036 ind[5][0] = ibnd[2]; ind[5][1] = ibnd[1]; ind[5][2] = ibnd[0];
00037 }
00038
00039
00041 void TriLinCuts (double *vlb, double *vub, int *varIndices,
00042 std::vector <std::vector <int> > &cutIndices,
00043 std::vector <std::vector <double> > &cutCoeff,
00044 std::vector <double> &cutLb,
00045 std::vector <double> &cutUb) {
00046
00047
00048 int v1, v2, v3, v4;
00049
00050 int defcons_size = 20;
00051
00052 double *bnd = new double [12];
00053
00054 CouNumber cf = 1.;
00055
00056 int **ind;
00057 ind = new int*[6];
00058 for(int i=0; i<6; i++) {
00059 ind[i] = new int[6];
00060 }
00061
00062 int *ibnd;
00063 ibnd = new int[3];
00064 ibnd[0] = varIndices[0]; ibnd[1] = varIndices[1]; ibnd[2] = varIndices[2];
00065 #ifdef DEBUG
00066 std::cout << "ibnd[0] =" << ibnd[0] << " ibnd[1] =" << ibnd[1] << " ibnd[2] =" << ibnd[2] << std::endl;
00067 std::cout << "vlb[ibnd[0]] =" << vlb[ibnd[0]] << " vub[ibnd[0]] =" << vub[ibnd[0]] << std::endl;
00068 std::cout << "vlb[ibnd[1]] =" << vlb[ibnd[1]] << " vub[ibnd[1]] =" << vub[ibnd[1]] << std::endl;
00069 std::cout << "vlb[ibnd[2]] =" << vlb[ibnd[2]] << " vub[ibnd[2]] =" << vub[ibnd[2]] << std::endl;
00070 #endif
00071
00072
00073 permutation3(ind,ibnd);
00074
00075 int i, flag=0, idx=0;
00076 i = 0;
00077 while(i < 6 && flag == 0) {
00078 if(vlb[ind[i][0]] >=0 && vlb[ind[i][1]] >=0 && vlb[ind[i][2]] <=0 && vub[ind[i][2]] >=0) {
00079 idx = i;
00080 flag = 1;
00081 }
00082 i++;
00083 }
00084 i = 0;
00085 while(i < 6 && flag == 0) {
00086 if(vlb[ind[i][0]] >=0 && vlb[ind[i][1]] <=0 && vlb[ind[i][2]] <=0
00087 && vub[ind[i][1]] >=0 && vub[ind[i][2]] >=0) {
00088 idx = i;
00089 flag = 2;
00090 }
00091 i++;
00092 }
00093 i = 0;
00094 while(i < 6 && flag == 0) {
00095 if(vlb[ind[i][0]] <=0 && vlb[ind[i][1]] <=0 && vlb[ind[i][2]] <=0
00096 && vub[ind[i][0]] >=0 && vub[ind[i][1]] >=0 && vub[ind[i][2]] >=0) {
00097 idx = i;
00098 flag = 3;
00099 }
00100 i++;
00101 }
00102 i = 0;
00103 while(i < 6 && flag == 0) {
00104 if(vlb[ind[i][0]] >=0 && vlb[ind[i][1]] <=0 && vlb[ind[i][2]] <=0
00105 && vub[ind[i][1]] >=0 && vub[ind[i][2]] <=0) {
00106 idx = i;
00107 flag = 4;
00108 }
00109 i++;
00110 }
00111 i = 0;
00112 while(i < 6 && flag == 0) {
00113 if(vlb[ind[i][0]] <=0 && vlb[ind[i][1]] <=0 && vlb[ind[i][2]] <=0
00114 && vub[ind[i][0]] >=0 && vub[ind[i][1]] >=0 && vub[ind[i][2]] <=0) {
00115 idx = i;
00116 flag = 5;
00117 }
00118 i++;
00119 }
00120 i = 0;
00121 while(i < 6 && flag == 0) {
00122 if(vlb[ind[i][0]] <=0 && vub[ind[i][0]] >=0 && vub[ind[i][1]] <=0 && vub[ind[i][2]] <=0) {
00123 idx = i;
00124 flag = 6;
00125 }
00126 i++;
00127 }
00128 i = 0;
00129 while(i < 6 && flag == 0) {
00130 if(vlb[ind[i][0]] >=0 && vlb[ind[i][1]] >=0 && vlb[ind[i][2]] >=0) {
00131 idx = i;
00132 flag = 7;
00133 }
00134 i++;
00135 }
00136 i = 0;
00137 while(i < 6 && flag == 0) {
00138 if(vlb[ind[i][0]] >=0 && vlb[ind[i][1]] >=0 && vlb[ind[i][2]] <=0 && vub[ind[i][2]] <=0) {
00139 idx = i;
00140 flag = 8;
00141 }
00142 i++;
00143 }
00144 i = 0;
00145 while(i < 6 && flag == 0) {
00146 if(vlb[ind[i][0]] >=0 && vlb[ind[i][1]] <=0 && vub[ind[i][1]] <=0
00147 && vlb[ind[i][2]] <=0 && vub[ind[i][2]] <=0) {
00148 idx = i;
00149 flag = 9;
00150 }
00151 i++;
00152 }
00153 i = 0;
00154 while(i < 6 && flag == 0) {
00155 if(vlb[ind[i][0]] <=0 && vub[ind[i][0]] <=0 && vlb[ind[i][1]] <=0 && vub[ind[i][1]] <=0
00156 && vlb[ind[i][2]] <=0 && vub[ind[i][2]] <=0) {
00157 idx = i;
00158 flag = 10;
00159 }
00160 i++;
00161 }
00162
00163 if (flag==0) {
00164 std::cout << "ERROR: case not implemented" << std::endl;
00165 exit(0);
00166 }
00167
00168
00169 v1 = ind[idx][0];
00170 v2 = ind[idx][1];
00171 v3 = ind[idx][2];
00172 v4 = varIndices [3];
00173
00174
00175 double xL1 = cf*vlb[v1]; double xU1 = cf*vub[v1];
00176 double xL2 = vlb[v2]; double xU2 = vub[v2];
00177 double xL3 = vlb[v3]; double xU3 = vub[v3];
00178
00179 #define prepareVectors(a) { \
00180 \
00181 int size = (int) (cutIndices.size ()); \
00182 \
00183 for (int i=0; i<a; i++) { \
00184 \
00185 cutIndices. push_back (std::vector <int> ()); \
00186 cutCoeff. push_back (std::vector <double> ()); \
00187 \
00188 cutLb. push_back (-COUENNE_INFINITY); \
00189 cutUb. push_back ( COUENNE_INFINITY); \
00190 \
00191 for (int j=0; j<4; j++) { \
00192 \
00193 cutIndices [size+i].push_back (-1); \
00194 cutCoeff [size+i].push_back (0.); \
00195 } \
00196 } \
00197 }
00198
00199
00200
00201
00202 if(flag == 1) {
00203 #ifdef DEBUG
00204 std::cout << " -- case 1 --" << std::endl;
00205 #endif
00206
00207 double theta = xL1*xU2*xU3 - xU1*xU2*xL3 - xL1*xL2*xU3 + xU1*xL2*xU3;
00208 double theta1 = xU1*xL2*xL3 - xU1*xU2*xU3 - xL1*xL2*xL3 + xL1*xU2*xL3;
00209
00210 defcons_size = 12;
00211
00212 prepareVectors(defcons_size);
00213
00214 for(int ii = 0; ii < defcons_size; ii++) {
00215
00216 cutIndices [ii][0] = v1;
00217 cutIndices [ii][1] = v2;
00218 cutIndices [ii][2] = v3;
00219 cutIndices [ii][3] = v4;
00220
00221 cutCoeff [ii][3] = 1.;
00222 }
00223
00224 cutCoeff [0][0] = -xU2*xU3; cutCoeff [0][1] = -xU1*xU3; cutCoeff [0][2] = -xU1*xU2; bnd[0] = - 2.*xU1*xU2*xU3;
00225 cutCoeff [1][0] = -xU2*xL3; cutCoeff [1][1] = -xL1*xU3; cutCoeff [1][2] = -xL1*xU2; bnd[1] = - xL1*xU2*xL3 - xL1*xU2*xU3;
00226 cutCoeff [2][0] = -xU2*xL3; cutCoeff [2][1] = -xL1*xL3; cutCoeff [2][2] = -xL1*xL2; bnd[2] = - xL1*xU2*xL3 - xL1*xL2*xL3;
00227 cutCoeff [3][0] = -xL2*xU3; cutCoeff [3][1] = -xU1*xL3; cutCoeff [3][2] = -xU1*xL2; bnd[3] = - xU1*xL2*xU3 - xU1*xL2*xL3;
00228 cutCoeff [4][0] = -xL2*xL3; cutCoeff [4][1] = -xU1*xL3; cutCoeff [4][2] = -xL1*xL2; bnd[4] = - xU1*xL2*xL3 - xL1*xL2*xL3;
00229 cutCoeff [5][0] = -xL2*xU3; cutCoeff [5][1] = -xL1*xU3; cutCoeff [5][2] = -(theta/(xU3-xL3));
00230 bnd[5] = (-(theta*xL3)/(xU3-xL3)) - xL1*xU2*xU3 - xU1*xL2*xU3 + xU1*xU2*xL3;
00231
00232 cutCoeff [6][0] = -xU2*xL3; cutCoeff [6][1] = -xU1*xL3; cutCoeff [6][2] = -xU1*xU2; bnd[6] = - 2.*xU1*xU2*xL3;
00233 cutCoeff [7][0] = -xL2*xL3; cutCoeff [7][1] = -xU1*xU3; cutCoeff [7][2] = -xU1*xL2; bnd[7] = - xU1*xL2*xU3 - xU1*xL2*xL3;
00234 cutCoeff [8][0] = -xU2*xU3; cutCoeff [8][1] = -xL1*xU3; cutCoeff [8][2] = -xL1*xL2; bnd[8] = - xL1*xU2*xU3 - xL1*xL2*xU3;
00235 cutCoeff [9][0] = -xU2*xU3; cutCoeff [9][1] = -xL1*xL3; cutCoeff [9][2] = -xL1*xU2; bnd[9] = - xL1*xU2*xU3 - xL1*xU2*xL3;
00236 cutCoeff [10][0] = -xL2*xU3; cutCoeff [10][1] = -xU1*xU3; cutCoeff [10][2] = -xL1*xL2; bnd[10] = - xU1*xL2*xU3 - xL1*xL2*xU3;
00237 cutCoeff [11][0] = -xL2*xL3; cutCoeff [11][1] = -xL1*xL3; cutCoeff [11][2] = -(theta1/(xL3-xU3)); cutCoeff [11][3] = 1.;
00238 bnd[11] = (-(theta1*xU3)/(xL3-xU3)) - xU1*xL2*xL3 - xL1*xU2*xL3 + xU1*xU2*xU3;
00239 }
00240
00241
00242
00243
00244 if(flag == 2) {
00245 #ifdef DEBUG
00246 std::cout << " -- case 2 --" << std::endl;
00247 #endif
00248
00249 defcons_size = 12;
00250
00251 prepareVectors (defcons_size);
00252
00253
00254 ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3;
00255 permutation3(ind,ibnd);
00256 int i, flagg=0, idx=0;
00257 i = 0;
00258 while(i < 6 && flagg == 0) {
00259 if(vub[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]]
00260 <= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]])
00261 {
00262 idx = i;
00263 flagg = 1;
00264 }
00265 i++;
00266 }
00267 if (flagg==0) {
00268 std::cout << "ERROR!!!" << std::endl; exit(0);
00269 }
00270 v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
00271
00272 double xL1(cf*vlb[v1]); double xU1(cf*vub[v1]);
00273 double xL2(vlb[v2]); double xU2(vub[v2]);
00274 double xL3(vlb[v3]); double xU3(vub[v3]);
00275
00276 for(int ii = 0; ii < defcons_size; ii++) {
00277
00278 cutIndices [ii][0] = v1;
00279 cutIndices [ii][1] = v2;
00280 cutIndices [ii][2] = v3;
00281 cutIndices [ii][3] = v4;
00282 cutCoeff [ii][3] = 1.;
00283 }
00284
00285 double theta1 = xL1*xL2*xL3 - xU1*xU2*xL3 - xL1*xL2*xU3 + xU1*xL2*xU3;
00286 double theta2 = xU1*xL2*xU3 - xU1*xU2*xL3 - xL1*xL2*xU3 + xL1*xU2*xU3;
00287
00288 cutCoeff [0][0] = -xU2*xU3; cutCoeff [0][1] = -xU1*xU3; cutCoeff [0][2] = -xU1*xU2; bnd[0] = - 2.*xU1*xU2*xU3;
00289 cutCoeff [1][0] = -xL2*xL3; cutCoeff [1][1] = -xU1*xL3; cutCoeff [1][2] = -xU1*xL2; bnd[1] = - 2.*xU1*xL2*xL3;
00290 cutCoeff [2][0] = -xU2*xL3; cutCoeff [2][1] = -xL1*xU3; cutCoeff [2][2] = -xL1*xU2; bnd[2] = - xL1*xU2*xL3 - xL1*xU2*xU3;
00291 cutCoeff [3][0] = -xU2*xL3; cutCoeff [3][1] = -xL1*xL3; cutCoeff [3][2] = -xL1*xL2; bnd[3] = - xL1*xU2*xL3 - xL1*xL2*xL3;
00292 cutCoeff [4][0] = -xL2*xU3; cutCoeff [4][1] = -(theta1/(xL2-xU2)); cutCoeff [4][2] = -xL1*xL2;
00293 bnd[4] = (-(theta1*xU2)/(xL2-xU2)) - xL1*xL2*xL3 - xU1*xL2*xU3 + xU1*xU2*xL3;
00294 cutCoeff [5][0] = -xL2*xU3; cutCoeff [5][1] = -xL1*xU3; cutCoeff [5][2] = -(theta2/(xU3-xL3));
00295 bnd[5] = (-(theta2*xL3)/(xU3-xL3)) - xU1*xL2*xU3 - xL1*xU2*xU3 + xU1*xU2*xL3;
00296
00297 if ( vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vub[v3]
00298 >= vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vub[v2]*vub[v3]) {
00299
00300 double theta1c = xU1*xL2*xL3 - xU1*xU2*xU3 - xL1*xL2*xL3 + xL1*xU2*xL3;
00301 double theta2c = xU1*xL2*xL3 - xU1*xU2*xU3 - xL1*xL2*xL3 + xL1*xL2*xU3;
00302
00303 cutCoeff [6][0] = -xL2*xU3; cutCoeff [6][1] = -xU1*xU3; cutCoeff [6][2] = -xU1*xL2; bnd[6] = - 2.*xU1*xL2*xU3;
00304 cutCoeff [7][0] = -xU2*xL3; cutCoeff [7][1] = -xU1*xL3; cutCoeff [7][2] = -xU1*xU2; bnd[7] = - 2.*xU1*xU2*xL3;
00305 cutCoeff [8][0] = -xU2*xU3; cutCoeff [8][1] = -xL1*xU3; cutCoeff [8][2] = -xL1*xL2; bnd[8] = - xL1*xU2*xU3 - xL1*xL2*xU3;
00306 cutCoeff [9][0] = -xU2*xU3; cutCoeff [9][1] = -xL1*xL3; cutCoeff [9][2] = -xL1*xU2; bnd[9] = - xL1*xU2*xU3 - xL1*xU2*xL3;
00307 cutCoeff [10][0] = -xL2*xL3; cutCoeff [10][1] = -xL1*xL3; cutCoeff [10][2] = -(theta1c/(xL3-xU3));
00308 bnd[10] = (-(theta1c*xU3)/(xL3-xU3)) - xU1*xL2*xL3 - xL1*xU2*xL3 + xU1*xU2*xU3;
00309 cutCoeff [11][0] = -xL2*xL3; cutCoeff [11][1] = -(theta2c/(xL2-xU2)); cutCoeff [11][2] = -xL1*xL2;
00310 bnd[11] = (-(theta2c*xU2)/(xL2-xU2)) - xU1*xL2*xL3 - xL1*xL2*xU3 + xU1*xU2*xU3;
00311
00312 } else {
00313 double theta1c = xU1*xU2*xU3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xL2*xU3;
00314 double theta2c = xU1*xU2*xU3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xU2*xL3;
00315
00316 cutCoeff [6][0] = -xL2*xU3; cutCoeff [6][1] = -xU1*xU3; cutCoeff [6][2] = -xU1*xL2; bnd[6] = - 2.*xU1*xL2*xU3;
00317 cutCoeff [7][0] = -xU2*xL3; cutCoeff [7][1] = -xU1*xL3; cutCoeff [7][2] = -xU1*xU2; bnd[7] = - 2.*xU1*xU2*xL3;
00318 cutCoeff [8][0] = -xL2*xL3; cutCoeff [8][1] = -xL1*xL3; cutCoeff [8][2] = -xL1*xU2; bnd[8] = - xL1*xL2*xL3 - xL1*xU2*xL3;
00319 cutCoeff [9][0] = -xL2*xL3; cutCoeff [9][1] = -xL1*xU3; cutCoeff [9][2] = -xL1*xL2; bnd[9] = - xL1*xL2*xL3 - xL1*xL2*xU3;
00320 cutCoeff [10][0] = -xU2*xU3; cutCoeff [10][1] = -xL1*xU3; cutCoeff [10][2] = -(theta1c/(xU3-xL3));
00321 bnd[10] = (-(theta1c*xL3)/(xU3-xL3)) - xU1*xU2*xU3 - xL1*xL2*xU3 + xU1*xL2*xL3;
00322 cutCoeff [11][0] = -xU2*xU3; cutCoeff [11][1] = -(theta2c/(xU2-xL2)); cutCoeff [11][2] = -xL1*xU2;
00323 bnd[11] = (-(theta2c*xL2)/(xU2-xL2)) - xU1*xU2*xU3 - xL1*xU2*xL3 + xU1*xL2*xL3;
00324 }
00325
00326 }
00327
00328
00329
00330
00331 if(flag == 3) {
00332 #ifdef DEBUG
00333 std::cout << " -- case 3 --" << std::endl;
00334 #endif
00335
00336 int last;
00337
00338 if(vub[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3] + vlb[v1]*vub[v2]*vub[v3]
00339 <= vlb[v1]*vlb[v2]*vlb[v3] + 2.*vub[v1]*vub[v2]*vub[v3] &&
00340 vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3]
00341 <= vlb[v1]*vub[v2]*vub[v3] + 2.*vub[v1]*vlb[v2]*vlb[v3] &&
00342 vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vlb[v3] + vlb[v1]*vub[v2]*vub[v3]
00343 <= vub[v1]*vlb[v2]*vub[v3] + 2.*vlb[v1]*vub[v2]*vlb[v3] &&
00344 vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3] + vlb[v1]*vub[v2]*vub[v3]
00345 <= vub[v1]*vub[v2]*vlb[v3] + 2.*vlb[v1]*vlb[v2]*vub[v3] ) {
00346
00347 double theta3x = 0.5*(xL1*xU2*xU3 + xL1*xL2*xL3 - xU1*xU2*xL3 - xU1*xL2*xU3)/(xL1-xU1);
00348 double theta3y = 0.5*(xU1*xL2*xU3 + xL1*xL2*xL3 - xU1*xU2*xL3 - xL1*xU2*xU3)/(xL2-xU2);
00349 double theta3z = 0.5*(xU1*xU2*xL3 + xL1*xL2*xL3 - xU1*xL2*xU3 - xL1*xU2*xU3)/(xL3-xU3);
00350 double theta3c = xL1*xL2*xL3 - theta3x*xL1 - theta3y*xL2 - theta3z*xL3;
00351
00352 defcons_size = 5;
00353 prepareVectors (defcons_size);
00354
00355 for(int ii = 0; ii < 5; ii++) {
00356
00357 cutIndices [ii][0] = v1;
00358 cutIndices [ii][1] = v2;
00359 cutIndices [ii][2] = v3;
00360 cutIndices [ii][3] = v4;
00361
00362 cutCoeff [ii][3] = 1.;
00363 }
00364
00365 cutCoeff [0][0] = -xU2*xL3; cutCoeff [0][1] = -xL1*xL3; cutCoeff [0][2] = -xL1*xU2; bnd[0] = - 2.*xL1*xU2*xL3;
00366 cutCoeff [1][0] = -xL2*xU3; cutCoeff [1][1] = -xL1*xU3; cutCoeff [1][2] = -xL1*xL2; bnd[1] = - 2.*xL1*xL2*xU3;
00367 cutCoeff [2][0] = -xU2*xU3; cutCoeff [2][1] = -xU1*xU3; cutCoeff [2][2] = -xU1*xU2; bnd[2] = - 2.*xU1*xU2*xU3;
00368 cutCoeff [3][0] = -xL2*xL3; cutCoeff [3][1] = -xU1*xL3; cutCoeff [3][2] = -xU1*xL2; bnd[3] = - 2.*xU1*xL2*xL3;
00369 cutCoeff [4][0] = -theta3x; cutCoeff [4][1] = -theta3y; cutCoeff [4][2] = -theta3z; bnd[4] = theta3c;
00370
00371 last=4;
00372
00373 } else if (vub[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3] + vlb[v1]*vub[v2]*vub[v3]
00374 >= vlb[v1]*vlb[v2]*vlb[v3] + 2.*vub[v1]*vub[v2]*vub[v3]) {
00375 #ifdef DEBUG
00376 std::cout << "else if " << std::endl;
00377 #endif
00378
00379 double theta1 = xU1*xU2*xL3 - xL1*xL2*xL3 - xU1*xU2*xU3 + xU1*xL2*xU3;
00380 double theta2 = xU1*xL2*xU3 - xL1*xL2*xL3 - xU1*xU2*xU3 + xL1*xU2*xU3;
00381 double theta3 = xU1*xU2*xL3 - xL1*xL2*xL3 - xU1*xU2*xU3 + xL1*xU2*xU3;
00382
00383 defcons_size = 6;
00384 prepareVectors (defcons_size);
00385
00386 for(int ii = 0; ii < 6; ii++) {
00387
00388 cutIndices [ii][0] = v1;
00389 cutIndices [ii][1] = v2;
00390 cutIndices [ii][2] = v3;
00391 cutIndices [ii][3] = v4;
00392
00393 cutCoeff [ii][3] = 1.;
00394 }
00395
00396 cutCoeff [0][0] = -xU2*xL3; cutCoeff [0][1] = -xL1*xL3; cutCoeff [0][2] = -xL1*xU2; bnd[0] = - 2.*xL1*xU2*xL3;
00397 cutCoeff [1][0] = -xL2*xU3; cutCoeff [1][1] = -xL1*xU3; cutCoeff [1][2] = -xL1*xL2; bnd[1] = - 2.*xL1*xL2*xU3;
00398 cutCoeff [2][0] = -xL2*xL3; cutCoeff [2][1] = -xU1*xL3; cutCoeff [2][2] = -xU1*xL2; bnd[2] = - 2.*xU1*xL2*xL3;
00399 cutCoeff [3][0] = -(theta1/(xU1-xL1)); cutCoeff [3][1] = -xU1*xU3; cutCoeff [3][2] = -xU1*xU2;
00400 bnd[3] = (-(theta1*xL1)/(xU1-xL1)) - xU1*xU2*xL3 - xU1*xL2*xU3 + xL1*xL2*xL3;
00401 cutCoeff [4][0] = -xU2*xU3; cutCoeff [4][1] = -xU1*xU3; cutCoeff [4][2] = -(theta2/(xU3-xL3));
00402 bnd[4] = (-(theta2*xL3)/(xU3-xL3)) - xU1*xL2*xU3 - xL1*xU2*xU3 + xL1*xL2*xL3;
00403 cutCoeff [5][0] = -xU2*xU3; cutCoeff [5][1] = -(theta3/(xU2-xL2)); cutCoeff [5][2] = -xU1*xU2;
00404 bnd[5] = (-(theta3*xL2)/(xU2-xL2)) - xU1*xU2*xL3 - xL1*xU2*xU3 + xL1*xL2*xL3;
00405
00406 last=5;
00407
00408 } else {
00409 #ifdef DEBUG
00410 std::cout << "else " << std::endl;
00411 #endif
00412
00413 ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3;
00414 permutation3(ind,ibnd);
00415 int i, flagg=0, idx=0;
00416 i = 0;
00417 while(i < 6 && flagg == 0) {
00418 if (((vlb[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]]
00419 >= vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]] + 2.*vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]]) ))
00420 {
00421 idx = i;
00422 flagg = 1;
00423 }
00424 i++;
00425 }
00426 if (flagg==0) {
00427 std::cout << "ERROR!!!" << std::endl; exit(0);
00428 }
00429 v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
00430
00431 double xL1(cf*vlb[v1]); double xU1(cf*vub[v1]);
00432 double xL2(vlb[v2]); double xU2(vub[v2]);
00433 double xL3(vlb[v3]); double xU3(vub[v3]);
00434
00435
00436
00437
00438 defcons_size = 6;
00439 prepareVectors (defcons_size);
00440
00441 for(int ii = 0; ii < 6; ii++) {
00442
00443 cutIndices [ii][0] = v1;
00444 cutIndices [ii][1] = v2;
00445 cutIndices [ii][2] = v3;
00446 cutIndices [ii][3] = v4;
00447
00448 cutCoeff [ii][3] = 1.;
00449 }
00450
00451 double theta1 = xL1*xL2*xL3 - xL1*xU2*xU3 - xU1*xL2*xL3 + xU1*xL2*xU3;
00452 double theta2 = xU1*xU2*xL3 - xL1*xU2*xU3 - xU1*xL2*xL3 + xU1*xL2*xU3;
00453 double theta3 = xL1*xL2*xL3 - xL1*xU2*xU3 - xU1*xL2*xL3 + xU1*xU2*xL3;
00454
00455 cutCoeff [0][0] = -xU2*xL3; cutCoeff [0][1] = -xL1*xL3; cutCoeff [0][2] = -xL1*xU2; bnd[0] = - 2.*xL1*xU2*xL3;
00456 cutCoeff [1][0] = -xL2*xU3; cutCoeff [1][1] = -xL1*xU3; cutCoeff [1][2] = -xL1*xL2; bnd[1] = - 2.*xL1*xL2*xU3;
00457 cutCoeff [2][0] = -xU2*xU3; cutCoeff [2][1] = -xU1*xU3; cutCoeff [2][2] = -xU1*xU2; bnd[2] = - 2.*xU1*xU2*xU3;
00458 cutCoeff [3][0] = -xL2*xL3; cutCoeff [3][1] = -(theta1/(xL2-xU2)); cutCoeff [3][2] = -xU1*xL2;
00459 bnd[3] = (-(theta1*xU2)/(xL2-xU2)) - xL1*xL2*xL3 - xU1*xL2*xU3 + xL1*xU2*xU3;
00460 cutCoeff [4][0] = -(theta2/(xU1-xL1)); cutCoeff [4][1] = -xU1*xL3; cutCoeff [4][2] = -xU1*xL2;
00461 bnd[4] = (-(theta2*xL1)/(xU1-xL1)) - xU1*xU2*xL3 - xU1*xL2*xU3 + xL1*xU2*xU3;
00462 cutCoeff [5][0] = -xL2*xL3; cutCoeff [5][1] = -xU1*xL3; cutCoeff [5][2] = -(theta3/(xL3-xU3));
00463 bnd[5] = (-(theta3*xU3)/(xL3-xU3)) - xL1*xL2*xL3 - xU1*xU2*xL3 + xL1*xU2*xU3;
00464
00465 last=5;
00466 }
00467
00468 if(vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3]
00469 >= vub[v1]*vub[v2]*vub[v3] + 2.*vlb[v1]*vlb[v2]*vlb[v3] &&
00470 vlb[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3] + vub[v1]*vub[v2]*vub[v3]
00471 >= vub[v1]*vlb[v2]*vlb[v3] + 2.*vlb[v1]*vub[v2]*vub[v3] &&
00472 vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3] + vub[v1]*vub[v2]*vub[v3]
00473 >= vlb[v1]*vub[v2]*vlb[v3] + 2.*vub[v1]*vlb[v2]*vub[v3] &&
00474 vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vub[v2]*vlb[v3] + vub[v1]*vub[v2]*vub[v3]
00475 >= vlb[v1]*vlb[v2]*vub[v3] + 2.*vub[v1]*vub[v2]*vlb[v3] ) {
00476
00477 #ifdef DEBUG
00478 std::cout << "2 - if " << std::endl;
00479 #endif
00480 double theta3x = 0.5*(xU1*xL2*xL3 + xU1*xU2*xU3 - xL1*xL2*xU3 - xL1*xU2*xL3)/(xU1-xL1);
00481 double theta3y = 0.5*(xL1*xU2*xL3 + xU1*xU2*xU3 - xL1*xL2*xU3 - xU1*xL2*xL3)/(xU2-xL2);
00482 double theta3z = 0.5*(xL1*xL2*xU3 + xU1*xU2*xU3 - xL1*xU2*xL3 - xU1*xL2*xL3)/(xU3-xL3);
00483 double theta3c = xU1*xU2*xU3 - theta3x*xU1 - theta3y*xU2 - theta3z*xU3;
00484
00485 prepareVectors (5);
00486
00487 for(int ii = last+1; ii <= last+5; ii++) {
00488
00489 cutIndices [ii][0] = v1;
00490 cutIndices [ii][1] = v2;
00491 cutIndices [ii][2] = v3;
00492 cutIndices [ii][3] = v4;
00493 cutCoeff [ii][3] = 1.;
00494 }
00495
00496 cutCoeff [last+1][0] = -xL2*xL3; cutCoeff [last+1][1] = -xL1*xL3; cutCoeff [last+1][2] = -xL1*xL2; bnd[last+1] = - 2.*xL1*xL2*xL3;
00497 cutCoeff [last+2][0] = -xU2*xL3; cutCoeff [last+2][1] = -xU1*xL3; cutCoeff [last+2][2] = -xU1*xU2; bnd[last+2] = - 2.*xU1*xU2*xL3;
00498 cutCoeff [last+3][0] = -xL2*xU3; cutCoeff [last+3][1] = -xU1*xU3; cutCoeff [last+3][2] = -xU1*xL2; bnd[last+3] = - 2.*xU1*xL2*xU3;
00499 cutCoeff [last+4][0] = -xU2*xU3; cutCoeff [last+4][1] = -xL1*xU3; cutCoeff [last+4][2] = -xL1*xU2; bnd[last+4] = - 2.*xL1*xU2*xU3;
00500 cutCoeff [last+5][0] = -theta3x; cutCoeff [last+5][1] = -theta3y; cutCoeff [last+5][2] = -theta3z; bnd[last+5] = theta3c;
00501
00502 defcons_size = last+6;
00503
00504 } else if (vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3]
00505 <= vub[v1]*vub[v2]*vub[v3] + 2.*vlb[v1]*vlb[v2]*vlb[v3]) {
00506 #ifdef DEBUG
00507 std::cout << "2 - else if" << std::endl;
00508 #endif
00509
00510 double theta1 = xU1*xL2*xL3 - xU1*xU2*xU3 - xL1*xL2*xL3 + xL1*xU2*xL3;
00511 double theta2 = xL1*xL2*xU3 - xU1*xU2*xU3 - xL1*xL2*xL3 + xL1*xU2*xL3;
00512 double theta3 = xL1*xL2*xU3 - xU1*xU2*xU3 - xL1*xL2*xL3 + xU1*xL2*xL3;
00513
00514 prepareVectors (6);
00515
00516 for(int ii = last+1; ii <= last+6; ii++) {
00517
00518 cutIndices [ii][0] = v1;
00519 cutIndices [ii][1] = v2;
00520 cutIndices [ii][2] = v3;
00521 cutIndices [ii][3] = v4;
00522
00523 cutCoeff [ii][3] = 1.;
00524 }
00525
00526 cutCoeff [last+1][0] = -xU2*xL3; cutCoeff [last+1][1] = -xU1*xL3; cutCoeff [last+1][2] = -xU1*xU2; bnd[last+1] = - 2.*xU1*xU2*xL3;
00527 cutCoeff [last+2][0] = -xL2*xU3; cutCoeff [last+2][1] = -xU1*xU3; cutCoeff [last+2][2] = -xU1*xL2; bnd[last+2] = - 2.*xU1*xL2*xU3;
00528 cutCoeff [last+3][0] = -xU2*xU3; cutCoeff [last+3][1] = -xL1*xU3; cutCoeff [last+3][2] = -xL1*xU2; bnd[last+3] = - 2.*xL1*xU2*xU3;
00529 cutCoeff [last+4][0] = -xL2*xL3; cutCoeff [last+4][1] = -xL1*xL3; cutCoeff [last+4][2] = -(theta1/(xL3-xU3));
00530 bnd[last+4] = (-(theta1*xU3)/(xL3-xU3)) - xU1*xL2*xL3 - xL1*xU2*xL3 + xU1*xU2*xU3;
00531 cutCoeff [last+5][0] = -(theta2/(xL1-xU1)); cutCoeff [last+5][1] = -xL1*xL3; cutCoeff [last+5][2] = -xL1*xL2;
00532 bnd[last+5] = (-(theta2*xU1)/(xL1-xU1)) - xL1*xL2*xU3 - xL1*xU2*xL3 + xU1*xU2*xU3;
00533 cutCoeff [last+6][0] = - xL2*xL3; cutCoeff [last+6][1] = -(theta3/(xL2-xU2)); cutCoeff [last+6][2] = -xL1*xL2;
00534 bnd[last+6] = (-(theta3*xU2)/(xL2-xU2)) - xL1*xL2*xU3 - xU1*xL2*xL3 + xU1*xU2*xU3;
00535
00536 defcons_size = last+7;
00537
00538 } else
00539
00540 {
00541 #ifdef DEBUG
00542 std::cout << "2 - another else if" << std::endl;
00543 std::cout << "v1 = " << v1 << " v2 =" << v2 << " v3 =" << v3 << std::endl;
00544 #endif
00545
00546
00547 permutation3(ind,ibnd);
00548 int i, flagg=0, idx=0;
00549 i = 0;
00550 while(i < 6 && flagg == 0) {
00551 if (vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]] + vub[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
00552 <= vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + 2.*vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]])
00553 {
00554 idx = i;
00555 flagg = 1;
00556 }
00557 i++;
00558 }
00559 if (flagg==0) {
00560 std::cout << "ERROR!!!" << std::endl; exit(0);
00561 }
00562 v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
00563
00564 double xL1(cf*vlb[v1]); double xU1(cf*vub[v1]);
00565 double xL2(vlb[v2]); double xU2(vub[v2]);
00566 double xL3(vlb[v3]); double xU3(vub[v3]);
00567
00568 double theta1 = xL1*xL2*xU3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xU2*xL3;
00569 double theta2 = xU1*xU2*xU3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xU2*xL3;
00570 double theta3 = xL1*xL2*xU3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xU1*xU2*xU3;
00571
00572 prepareVectors (6);
00573
00574 for(int ii = last+1; ii <= last+6; ii++) {
00575
00576 cutIndices [ii][0] = v1;
00577 cutIndices [ii][1] = v2;
00578 cutIndices [ii][2] = v3;
00579 cutIndices [ii][3] = v4;
00580
00581 cutCoeff [ii][3] = 1.;
00582 }
00583
00584 cutCoeff [last+1][0] = -xL2*xL3; cutCoeff [last+1][1] = -xL1*xL3; cutCoeff [last+1][2] = -xL1*xL2; bnd[last+1] = - 2.*xL1*xL2*xL3;
00585 cutCoeff [last+2][0] = -xU2*xL3; cutCoeff [last+2][1] = -xU1*xL3; cutCoeff [last+2][2] = -xU1*xU2; bnd[last+2] = - 2.*xU1*xU2*xL3;
00586 cutCoeff [last+3][0] = -xL2*xU3; cutCoeff [last+3][1] = -xU1*xU3; cutCoeff [last+3][2] = -xU1*xL2; bnd[last+3] = - 2.*xU1*xL2*xU3;
00587 cutCoeff [last+4][0] = -(theta1/(xL1-xU1)); cutCoeff [last+4][1] = -xL1*xU3; cutCoeff [last+4][2] = -xL1*xU2;
00588 bnd[last+4] = (-(theta1*xU1)/(xL1-xU1)) - xL1*xL2*xU3 - xL1*xU2*xL3 + xU1*xL2*xL3;
00589 cutCoeff [last+5][0] = -xU2*xU3; cutCoeff [last+5][1] = -(theta2/(xU2-xL2)); cutCoeff [last+5][2] = -xL1*xU2;
00590 bnd[last+5] = (-(theta2*xL2)/(xU2-xL2)) - xU1*xU2*xU3 - xL1*xU2*xL3 + xU1*xL2*xL3;
00591 cutCoeff [last+6][0] = -xU2*xU3; cutCoeff [last+6][1] = -xL1*xU3; cutCoeff [last+6][2] = -(theta3/(xU3-xL3));
00592 bnd[last+6] = (-(theta3*xL3)/(xU3-xL3)) - xL1*xL2*xU3 - xU1*xU2*xU3 + xU1*xL2*xL3;
00593
00594 defcons_size = last+7;
00595 }
00596
00597 }
00598
00599
00600
00601
00602 if(flag == 4) {
00603 #ifdef DEBUG
00604 std::cout << " -- case 4 --" << std::endl;
00605 #endif
00606
00607 double theta = xU1*xL2*xU3 - xU1*xU2*xL3 - xL1*xL2*xU3 + xL1*xL2*xL3;
00608 double theta1 = xL1*xU2*xL3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xU1*xU2*xU3;
00609
00610 defcons_size = 12;
00611
00612 prepareVectors (defcons_size);
00613
00614 for(int ii = 0; ii < defcons_size; ii++) {
00615 cutIndices [ii][0] = v1;
00616 cutIndices [ii][1] = v2;
00617 cutIndices [ii][2] = v3;
00618 cutIndices [ii][3] = v4;
00619 cutCoeff [ii][3] = 1.;
00620 }
00621
00622 cutCoeff [0][0] = -xL2*xL3; cutCoeff [0][1] = -xU1*xL3; cutCoeff [0][2] = -xU1*xL2; bnd[0] = - 2.*xU1*xL2*xL3;
00623 cutCoeff [1][0] = -xU2*xL3; cutCoeff [1][1] = -xL1*xU3; cutCoeff [1][2] = -xL1*xU2; bnd[1] = - xL1*xU2*xL3 - xL1*xU2*xU3;
00624 cutCoeff [2][0] = -xU2*xL3; cutCoeff [2][1] = -xL1*xL3; cutCoeff [2][2] = -xL1*xL2; bnd[2] = - xL1*xU2*xL3 - xL1*xL2*xL3;
00625 cutCoeff [3][0] = -xL2*xU3; cutCoeff [3][1] = -xU1*xU3; cutCoeff [3][2] = -xU1*xU2; bnd[3] = - xU1*xL2*xU3 - xU1*xU2*xU3;
00626 cutCoeff [4][0] = -xU2*xU3; cutCoeff [4][1] = -xL1*xU3; cutCoeff [4][2] = -xU1*xU2; bnd[4] = - xL1*xU2*xU3 - xU1*xU2*xU3;
00627 cutCoeff [5][0] = -xL2*xU3; cutCoeff [5][1] = -(theta/(xL2-xU2)); cutCoeff [5][2] = -xL1*xL2;
00628 bnd[5] = (-(theta*xU2)/(xL2-xU2)) - xU1*xL2*xU3 - xL1*xL2*xL3 + xU1*xU2*xL3;
00629
00630 cutCoeff [6][0] = -xU2*xL3; cutCoeff [6][1] = -xU1*xL3; cutCoeff [6][2] = -xU1*xU2; bnd[6] = - 2.*xU1*xU2*xL3;
00631 cutCoeff [7][0] = -xU2*xU3; cutCoeff [7][1] = -xU1*xU3; cutCoeff [7][2] = -xU1*xL2; bnd[7] = - xU1*xU2*xU3 - xU1*xL2*xU3;
00632 cutCoeff [8][0] = -xL2*xL3; cutCoeff [8][1] = -xL1*xL3; cutCoeff [8][2] = -xL1*xU2; bnd[8] = - xL1*xU2*xL3 - xL1*xL2*xL3;
00633 cutCoeff [9][0] = -xL2*xL3; cutCoeff [9][1] = -xL1*xU3; cutCoeff [9][2] = -xL1*xL2; bnd[9] = - xL1*xL2*xU3 - xL1*xL2*xL3;
00634 cutCoeff [10][0] = -xL2*xU3; cutCoeff [10][1] = -xL1*xU3; cutCoeff [10][2] = -xU1*xL2; bnd[10] = - xL1*xL2*xU3 - xU1*xL2*xU3;
00635 cutCoeff [11][0] = -xU2*xU3; cutCoeff [11][1] = -(theta1/(xU2-xL2)); cutCoeff [11][2] = -xL1*xU2;
00636 bnd[11] = (-(theta1*xL2)/(xU2-xL2)) - xL1*xU2*xL3 - xU1*xU2*xU3 + xU1*xL2*xL3;
00637
00638 }
00639
00640
00641
00642
00643 if(flag == 5) {
00644 #ifdef DEBUG
00645 std::cout << " -- case 5 --" << std::endl;
00646 std::cout << "v1 = " << v1 << " v2 =" << v2 << " v3 =" << v3 << std::endl;
00647 #endif
00648
00649 defcons_size = 12;
00650 prepareVectors (defcons_size);
00651
00652
00653 ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3;
00654 ind[0][0] = ibnd[0]; ind[0][1] = ibnd[1]; ind[0][2] = ibnd[2];
00655 ind[1][0] = ibnd[1]; ind[1][1] = ibnd[0]; ind[1][2] = ibnd[2];
00656 int i, flagg=0, idx=0;
00657 i = 0;
00658 while(i < 2 && flagg == 0) {
00659 if(vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
00660 >= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]])
00661 {
00662 idx = i;
00663 flagg = 1;
00664 }
00665 i++;
00666 }
00667 if (flagg==0) {
00668 std::cout << "ERROR!!!" << std::endl; exit(0);
00669 }
00670 v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
00671 #ifdef DEBUG
00672 std::cout << "v1 = " << v1 << " v2 =" << v2 << " v3 =" << v3 << std::endl;
00673 #endif
00674
00675 double xL1 = cf*vlb[v1]; double xU1 = cf*vub[v1];
00676 double xL2 = vlb[v2]; double xU2 = vub[v2];
00677 double xL3 = vlb[v3]; double xU3 = vub[v3];
00678
00679 for(int ii = 0; ii < defcons_size; ii++) {
00680
00681 cutIndices [ii][0] = v1;
00682 cutIndices [ii][1] = v2;
00683 cutIndices [ii][2] = v3;
00684 cutIndices [ii][3] = v4;
00685 cutCoeff [ii][3] = 1.;
00686 }
00687
00688 if(vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vub[v3]
00689 <= vub[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3]) {
00690 #ifdef DEBUG
00691 std::cout << " -- 5 if --" << std::endl;
00692 #endif
00693
00694 double theta1 = xU1*xU2*xL3 - xL1*xL2*xL3 - xU1*xU2*xU3 + xU1*xL2*xU3;
00695 double theta2 = xU1*xU2*xL3 - xL1*xL2*xL3 - xU1*xU2*xU3 + xL1*xU2*xU3;
00696
00697 cutCoeff [0][0] = -xU2*xL3; cutCoeff [0][1] = -xL1*xL3; cutCoeff [0][2] = -xL1*xU2; bnd[0] = - 2.*xL1*xU2*xL3;
00698 cutCoeff [1][0] = -xL2*xL3; cutCoeff [1][1] = -xU1*xL3; cutCoeff [1][2] = -xU1*xL2; bnd[1] = - 2.*xU1*xL2*xL3;
00699 cutCoeff [2][0] = -xU2*xU3; cutCoeff [2][1] = -xL1*xU3; cutCoeff [2][2] = -xL1*xL2; bnd[2] = - xL1*xU2*xU3 - xL1*xL2*xU3;
00700 cutCoeff [3][0] = -xL2*xU3; cutCoeff [3][1] = -xU1*xU3; cutCoeff [3][2] = -xL1*xL2; bnd[3] = - xU1*xL2*xU3 - xL1*xL2*xU3;
00701 cutCoeff [4][0] = -(theta1/(xU1-xL1)); cutCoeff [4][1] = -xU1*xU3; cutCoeff [4][2] = -xU1*xU2;
00702 bnd[4] = (-(theta1*xL1)/(xU1-xL1)) - xU1*xU2*xL3 - xU1*xL2*xU3 + xL1*xL2*xL3;
00703 cutCoeff [5][0] = -xU2*xU3; cutCoeff [5][1] = -(theta2/(xU2-xL2)); cutCoeff [5][2] = -xU1*xU2;
00704 bnd[5] = (-(theta2*xL2)/(xU2-xL2)) - xU1*xU2*xL3 - xL1*xU2*xU3 + xL1*xL2*xL3;
00705
00706 } else {
00707 #ifdef DEBUG
00708 std::cout << " -- 5 else --" << std::endl;
00709 #endif
00710 double theta1 = xL1*xL2*xL3 - xU1*xU2*xL3 - xL1*xL2*xU3 + xL1*xU2*xU3;
00711 double theta2 = xL1*xL2*xL3 - xU1*xU2*xL3 - xL1*xL2*xU3 + xU1*xL2*xU3;
00712
00713 cutCoeff [0][0] = -xU2*xL3; cutCoeff [0][1] = -xL1*xL3; cutCoeff [0][2] = -xL1*xU2; bnd[0] = - 2.*xL1*xU2*xL3;
00714 cutCoeff [1][0] = -xL2*xL3; cutCoeff [1][1] = -xU1*xL3; cutCoeff [1][2] = -xU1*xL2; bnd[1] = - 2.*xU1*xL2*xL3;
00715 cutCoeff [2][0] = -xL2*xU3; cutCoeff [2][1] = -xU1*xU3; cutCoeff [2][2] = -xU1*xU2; bnd[2] = - xU1*xL2*xU3 - xU1*xU2*xU3;
00716 cutCoeff [3][0] = -xU2*xU3; cutCoeff [3][1] = -xL1*xU3; cutCoeff [3][2] = -xU1*xU2; bnd[3] = - xL1*xU2*xU3 - xU1*xU2*xU3;
00717 cutCoeff [4][0] = -(theta1/(xL1-xU1)); cutCoeff [4][1] = -xL1*xU3; cutCoeff [4][2] = -xL1*xL2;
00718 bnd[4] = (-(theta1*xU1)/(xL1-xU1)) - xL1*xL2*xL3 - xL1*xU2*xU3 + xU1*xU2*xL3;
00719 cutCoeff [5][0] = -xL2*xU3; cutCoeff [5][1] = -(theta2/(xL2-xU2)); cutCoeff [5][2] = -xL1*xL2;
00720 bnd[5] = (-(theta2*xU2)/(xL2-xU2)) - xL1*xL2*xL3 - xU1*xL2*xU3 + xU1*xU2*xL3;
00721 }
00722
00723 double theta1c = xL1*xU2*xL3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xL2*xU3;
00724 double theta2c = xU1*xU2*xU3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xU2*xL3;
00725
00726 cutCoeff [6][0] = -xL2*xL3; cutCoeff [6][1] = -xL1*xL3; cutCoeff [6][2] = -xL1*xL2; bnd[6] = - 2.*xL1*xL2*xL3;
00727 cutCoeff [7][0] = -xU2*xL3; cutCoeff [7][1] = -xU1*xL3; cutCoeff [7][2] = -xU1*xU2; bnd[7] = - 2.*xU1*xU2*xL3;
00728 cutCoeff [8][0] = -xL2*xU3; cutCoeff [8][1] = -xL1*xU3; cutCoeff [8][2] = -xU1*xL2; bnd[8] = - xL1*xL2*xU3 - xU1*xL2*xU3;
00729 cutCoeff [9][0] = -xU2*xU3; cutCoeff [9][1] = -xU1*xU3; cutCoeff [9][2] = -xU1*xL2; bnd[9] = - xU1*xU2*xU3 - xU1*xL2*xU3;
00730 cutCoeff [10][0] = -(theta1c/(xL1-xU1)); cutCoeff [10][1] = -xL1*xU3; cutCoeff [10][2] = -xL1*xU2;
00731 bnd[10] = (-(theta1c*xU1)/(xL1-xU1)) - xL1*xU2*xL3 - xL1*xL2*xU3 + xU1*xL2*xL3;
00732 cutCoeff [11][0] = -xU2*xU3; cutCoeff [11][1] = -(theta2c/(xU2-xL2)); cutCoeff [11][2] = -xL1*xU2;
00733 bnd[11] = (-(theta2c*xL2)/(xU2-xL2)) - xU1*xU2*xU3 - xL1*xU2*xL3 + xU1*xL2*xL3;
00734
00735 }
00736
00737
00738
00739
00740 if(flag == 6) {
00741 #ifdef DEBUG
00742 std::cout << " -- case 6 --" << std::endl;
00743 #endif
00744
00745 double theta = xU1*xU2*xL3 - xL1*xL2*xL3 - xU1*xU2*xU3 + xU1*xL2*xU3;
00746 double theta1 = xL1*xU2*xL3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xL2*xU3;
00747
00748 defcons_size = 12;
00749 prepareVectors (defcons_size);
00750
00751 for (int ii = 0; ii < defcons_size; ii++) {
00752
00753 cutIndices [ii][0] = v1;
00754 cutIndices [ii][1] = v2;
00755 cutIndices [ii][2] = v3;
00756 cutIndices [ii][3] = v4;
00757
00758 cutCoeff [ii][3] = 1.;
00759 }
00760
00761 cutCoeff [0][0] = -xL2*xL3; cutCoeff [0][1] = -xU1*xL3; cutCoeff [0][2] = -xU1*xL2; bnd[0] = - 2.*xU1*xL2*xL3;
00762 cutCoeff [1][0] = -xU2*xU3; cutCoeff [1][1] = -xL1*xL3; cutCoeff [1][2] = -xL1*xU2; bnd[1] = - xL1*xU2*xU3 - xL1*xU2*xL3;
00763 cutCoeff [2][0] = -xU2*xU3; cutCoeff [2][1] = -xL1*xU3; cutCoeff [2][2] = -xL1*xL2; bnd[2] = - xL1*xU2*xU3 - xL1*xL2*xU3;
00764 cutCoeff [3][0] = -xL2*xU3; cutCoeff [3][1] = -xU1*xU3; cutCoeff [3][2] = -xL1*xL2; bnd[3] = - xU1*xL2*xU3 - xL1*xL2*xU3;
00765 cutCoeff [4][0] = -xU2*xL3; cutCoeff [4][1] = -xL1*xL3; cutCoeff [4][2] = -xU1*xU2; bnd[4] = - xL1*xU2*xL3 - xU1*xU2*xL3;
00766 cutCoeff [5][0] = -(theta/(xU1-xL1)); cutCoeff [5][1] = -xU1*xU3; cutCoeff [5][2] = -xU1*xU2;
00767 bnd[5] = (-(theta*xL1)/(xU1-xL1)) - xU1*xU2*xL3 - xU1*xL2*xU3 + xL1*xL2*xL3;
00768
00769 cutCoeff [6][0] = -xL2*xL3; cutCoeff [6][1] = -xL1*xL3; cutCoeff [6][2] = -xL1*xL2; bnd[6] = - 2.*xL1*xL2*xL3;
00770 cutCoeff [7][0] = -xL2*xU3; cutCoeff [7][1] = -xL1*xU3; cutCoeff [7][2] = -xU1*xL2; bnd[7] = - xL1*xL2*xU3 - xU1*xL2*xU3;
00771 cutCoeff [8][0] = -xU2*xU3; cutCoeff [8][1] = -xU1*xU3; cutCoeff [8][2] = -xU1*xL2; bnd[8] = - xU1*xU2*xU3 - xU1*xL2*xU3;
00772 cutCoeff [9][0] = -xU2*xU3; cutCoeff [9][1] = -xU1*xL3; cutCoeff [9][2] = -xU1*xU2; bnd[9] = - xU1*xU2*xU3 - xU1*xU2*xL3;
00773 cutCoeff [10][0] = -xU2*xL3; cutCoeff [10][1] = -xU1*xL3; cutCoeff [10][2] = -xL1*xU2; bnd[10] = - xL1*xU2*xL3 - xU1*xU2*xL3;
00774 cutCoeff [11][0] = -(theta1/(xL1-xU1)); cutCoeff [11][1] = -xL1*xU3; cutCoeff [11][2] = -xL1*xU2;
00775 bnd[11] = (-(theta1*xU1)/(xL1-xU1)) - xL1*xU2*xL3 - xL1*xL2*xU3 + xU1*xL2*xL3;
00776 }
00777
00778
00779
00780
00781 if(flag == 7) {
00782 #ifdef DEBUG
00783 std::cout << " -- case 7 --" << std::endl;
00784 #endif
00785
00786 defcons_size = 12;
00787 prepareVectors (defcons_size);
00788
00789 if ((vlb[v1]<=EPSILONT && vlb[v2]<=EPSILONT && vlb[v3]<=EPSILONT) ||
00790 (vlb[v1]==vlb[v2] && vlb[v1]==vlb[v3] && vub[v1]==vub[v2] && vub[v1]==vub[v3])) {
00791 #ifdef DEBUG
00792 std::cout << " -- epsilonT --" << std::endl;
00793 #endif
00794 } else {
00795
00796 ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3;
00797 permutation3(ind,ibnd);
00798 int i, flagg=0, idx=0;
00799 i = 0;
00800 while(i < 6 && flagg == 0) {
00801 if(vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
00802 <= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]] &&
00803 vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
00804 <= vub[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]])
00805 {
00806 idx = i;
00807 flagg = 1;
00808 }
00809 i++;
00810 }
00811 if (flagg==0) {
00812 std::cout << "ERROR!!!" << std::endl; exit(0);
00813 }
00814 }
00815 v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
00816
00817 double xL1 = cf*vlb[v1]; double xU1 = cf*vub[v1];
00818 double xL2 = vlb[v2]; double xU2 = vub[v2];
00819 double xL3 = vlb[v3]; double xU3 = vub[v3];
00820
00821
00822
00823
00824
00825
00826 double theta1 = xU1*xU2*xL3 - xL1*xU2*xU3 - xU1*xL2*xL3 + xU1*xL2*xU3;
00827 double theta2 = xL1*xL2*xU3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xU2*xL3;
00828
00829 for(int ii = 0; ii < defcons_size; ii++) {
00830 cutIndices [ii][0] = v1;
00831 cutIndices [ii][1] = v2;
00832 cutIndices [ii][2] = v3;
00833 cutIndices [ii][3] = v4;
00834
00835 cutCoeff [ii][3] = 1.;
00836 }
00837
00838 cutCoeff [0][0] = -xL2*xL3; cutCoeff [0][1] = -xL1*xL3; cutCoeff [0][2] = -xL1*xL2; bnd[0] = - 2.*xL1*xL2*xL3;
00839 cutCoeff [1][0] = -xU2*xU3; cutCoeff [1][1] = -xU1*xU3; cutCoeff [1][2] = -xU1*xU2; bnd[1] = - 2.*xU1*xU2*xU3;
00840 cutCoeff [2][0] = -xL2*xU3; cutCoeff [2][1] = -xL1*xU3; cutCoeff [2][2] = -xU1*xL2; bnd[2] = - xL1*xL2*xU3 - xU1*xL2*xU3;
00841 cutCoeff [3][0] = -xU2*xL3; cutCoeff [3][1] = -xU1*xL3; cutCoeff [3][2] = -xL1*xU2; bnd[3] = - xU1*xU2*xL3 - xL1*xU2*xL3;
00842 cutCoeff [4][0] = -(theta1/(xU1-xL1)); cutCoeff [4][1] = -xU1*xL3; cutCoeff [4][2] = -xU1*xL2;
00843 bnd[4] = (-(theta1*xL1)/(xU1-xL1)) - xU1*xU2*xL3 - xU1*xL2*xU3 + xL1*xU2*xU3;
00844 cutCoeff [5][0] = -(theta2/(xL1-xU1)); cutCoeff [5][1] = -xL1*xU3; cutCoeff [5][2] = -xL1*xU2;
00845 bnd[5] = (-(theta2*xU1)/(xL1-xU1)) - xL1*xL2*xU3 - xL1*xU2*xL3 + xU1*xL2*xL3;
00846
00847 cutCoeff [6][0] = -xL2*xL3; cutCoeff [6][1] = -xU1*xL3; cutCoeff [6][2] = -xU1*xU2; bnd[6] = - xU1*xU2*xL3 - xU1*xL2*xL3;
00848 cutCoeff [7][0] = -xU2*xL3; cutCoeff [7][1] = -xL1*xL3; cutCoeff [7][2] = -xU1*xU2; bnd[7] = - xU1*xU2*xL3 - xL1*xU2*xL3;
00849 cutCoeff [8][0] = -xL2*xL3; cutCoeff [8][1] = -xU1*xU3; cutCoeff [8][2] = -xU1*xL2; bnd[8] = - xU1*xL2*xU3 - xU1*xL2*xL3;
00850 cutCoeff [9][0] = -xU2*xU3; cutCoeff [9][1] = -xL1*xL3; cutCoeff [9][2] = -xL1*xU2; bnd[9] = - xL1*xU2*xU3 - xL1*xU2*xL3;
00851 cutCoeff [10][0] = -xL2*xU3; cutCoeff [10][1] = -xU1*xU3; cutCoeff [10][2] = -xL1*xL2; bnd[10] = - xU1*xL2*xU3 - xL1*xL2*xU3;
00852 cutCoeff [11][0] = -xU2*xU3; cutCoeff [11][1] = -xL1*xU3; cutCoeff [11][2] = -xL1*xL2; bnd[11] = - xL1*xU2*xU3 - xL1*xL2*xU3;
00853 }
00854
00855
00856
00857
00858 if(flag == 8) {
00859 #ifdef DEBUG
00860 std::cout << " -- case 8 --" << std::endl;
00861 #endif
00862
00863 defcons_size = 12;
00864 prepareVectors (defcons_size);
00865
00866
00867 ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3;
00868 permutation3(ind,ibnd);
00869 int i, flagg=0, idx=0;
00870 i = 0;
00871 while(i < 6 && flagg == 0) {
00872 if(vub[ind[i][2]] <=0 &&
00873 ((vlb[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
00874 >= vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]] &&
00875 vlb[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
00876 >= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]]) ||
00877 (vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
00878 >= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]])) )
00879 {
00880 idx = i;
00881 flagg = 1;
00882 }
00883 i++;
00884 }
00885 if (flagg==0) {
00886 std::cout << "ERROR!!!" << std::endl; exit(0);
00887 }
00888 v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
00889
00890 double xL1(cf*vlb[v1]); double xU1(cf*vub[v1]);
00891 double xL2(vlb[v2]); double xU2(vub[v2]);
00892 double xL3(vlb[v3]); double xU3(vub[v3]);
00893
00894 for(int ii = 0; ii < defcons_size; ii++) {
00895
00896 cutIndices [ii][0] = v1;
00897 cutIndices [ii][1] = v2;
00898 cutIndices [ii][2] = v3;
00899 cutIndices [ii][3] = v4;
00900 cutCoeff [ii][3] = 1.;
00901 }
00902
00903
00904
00905 cutCoeff [0][0] = -xU2*xL3; cutCoeff [0][1] = -xL1*xL3; cutCoeff [0][2] = -xL1*xL2; bnd[0] = - xL1*xU2*xL3 - xL1*xL2*xL3;
00906 cutCoeff [1][0] = -xU2*xL3; cutCoeff [1][1] = -xL1*xU3; cutCoeff [1][2] = -xL1*xU2; bnd[1] = - xL1*xU2*xL3 - xL1*xU2*xU3;
00907 cutCoeff [2][0] = -xL2*xU3; cutCoeff [2][1] = -xU1*xL3; cutCoeff [2][2] = -xU1*xL2; bnd[2] = - xU1*xL2*xU3 - xU1*xL2*xL3;
00908 cutCoeff [3][0] = -xL2*xU3; cutCoeff [3][1] = -xU1*xU3; cutCoeff [3][2] = -xU1*xU2; bnd[3] = - xU1*xL2*xU3 - xU1*xU2*xU3;
00909 cutCoeff [4][0] = -xL2*xL3; cutCoeff [4][1] = -xU1*xL3; cutCoeff [4][2] = -xL1*xL2; bnd[4] = - xU1*xL2*xL3 - xL1*xL2*xL3;
00910 cutCoeff [5][0] = -xU2*xU3; cutCoeff [5][1] = -xL1*xU3; cutCoeff [5][2] = -xU1*xU2; bnd[5] = - xU1*xU2*xU3 - xL1*xU2*xU3;
00911
00912 if(vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vub[v3]
00913 >= vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vub[v2]*vub[v3] &&
00914 vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vub[v3]
00915 >= vlb[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3]) {
00916
00917 double theta1c = xU1*xL2*xL3 - xU1*xU2*xU3 - xL1*xL2*xL3 + xL1*xU2*xL3;
00918 double theta2c = xU1*xL2*xU3 - xL1*xL2*xL3 - xU1*xU2*xU3 + xL1*xU2*xU3;
00919
00920 cutCoeff [6][0] = -xL2*xU3; cutCoeff [6][1] = -xL1*xU3; cutCoeff [6][2] = -xL1*xL2; bnd[6] = - 2.*xL1*xL2*xU3;
00921 cutCoeff [7][0] = -xU2*xL3; cutCoeff [7][1] = -xU1*xL3; cutCoeff [7][2] = -xU1*xU2; bnd[7] = - 2.*xU1*xU2*xL3;
00922 cutCoeff [8][0] = -xL2*xL3; cutCoeff [8][1] = -xU1*xU3; cutCoeff [8][2] = -xU1*xL2; bnd[8] = - xU1*xL2*xU3 - xU1*xL2*xL3;
00923 cutCoeff [9][0] = -xU2*xU3; cutCoeff [9][1] = -xL1*xL3; cutCoeff [9][2] = -xL1*xU2; bnd[9] = - xL1*xU2*xU3 - xL1*xU2*xL3;
00924 cutCoeff [10][0] = -xL2*xL3; cutCoeff [10][1] = -xL1*xL3; cutCoeff [10][2] = -(theta1c/(xL3-xU3));
00925 bnd[10] = (-(theta1c*xU3)/(xL3-xU3)) - xU1*xL2*xL3 - xL1*xU2*xL3 + xU1*xU2*xU3;
00926 cutCoeff [11][0] = -xU2*xU3; cutCoeff [11][1] = -xU1*xU3; cutCoeff [11][2] = -(theta2c/(xU3-xL3));
00927 bnd[11] = (-(theta2c*xL3)/(xU3-xL3)) - xU1*xL2*xU3 - xL1*xU2*xU3 + xL1*xL2*xL3;
00928
00929 } else if(vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vub[v2]*vub[v3]
00930 >= vlb[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3]) {
00931
00932 double theta1c = xL1*xL2*xL3 - xL1*xU2*xU3 - xU1*xL2*xL3 + xU1*xL2*xU3;
00933 double theta2c = xL1*xU2*xL3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xU1*xU2*xU3;
00934
00935 cutCoeff [6][0] = -xL2*xU3; cutCoeff [6][1] = -xL1*xU3; cutCoeff [6][2] = -xL1*xL2; bnd[6] = - 2.*xL1*xL2*xU3;
00936 cutCoeff [7][0] = -xU2*xL3; cutCoeff [7][1] = -xU1*xL3; cutCoeff [7][2] = -xU1*xU2; bnd[7] = - 2.*xU1*xU2*xL3;
00937 cutCoeff [8][0] = -xL2*xL3; cutCoeff [8][1] = -xL1*xL3; cutCoeff [8][2] = -xL1*xU2; bnd[8] = - xL1*xL2*xL3 - xL1*xU2*xL3;
00938 cutCoeff [9][0] = -xU2*xU3; cutCoeff [9][1] = -xU1*xU3; cutCoeff [9][2] = -xU1*xL2; bnd[9] = - xU1*xL2*xU3 - xU1*xU2*xU3;
00939 cutCoeff [10][0] = -xL2*xL3; cutCoeff [10][1] = -(theta1c/(xL2-xU2)); cutCoeff [10][2] = -xU1*xL2;
00940 bnd[10] = (-(theta1c*xU2)/(xL2-xU2)) - xL1*xL2*xL3 - xU1*xL2*xU3 + xL1*xU2*xU3;
00941 cutCoeff [11][0] = -xU2*xU3; cutCoeff [11][1] = -(theta2c/(xU2-xL2)); cutCoeff [11][2] = -xL1*xU2;
00942 bnd[11] = (-(theta2c*xL2)/(xU2-xL2)) - xL1*xU2*xL3 - xU1*xU2*xU3 + xU1*xL2*xL3;
00943 }
00944
00945 }
00946
00947
00948
00949
00950 if(flag == 9) {
00951 #ifdef DEBUG
00952 std::cout << " -- case 9 --" << std::endl;
00953 #endif
00954
00955 defcons_size = 12;
00956 prepareVectors (defcons_size);
00957
00958
00959 ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3;
00960 permutation3(ind,ibnd);
00961 int i, flagg=0, idx=0;
00962 i = 0;
00963 while(i < 6 && flagg == 0) {
00964 if(vlb[ind[i][0]] >=0 &&
00965 ((vlb[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
00966 <= vub[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]] &&
00967 vlb[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
00968 <= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]]) ||
00969 (vub[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]]
00970 <= vlb[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]] &&
00971 vub[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]]
00972 <= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]]) ))
00973 {
00974 idx = i;
00975 flagg = 1;
00976 }
00977 i++;
00978 }
00979 if (flagg==0) {
00980 std::cout << "ERROR!!!" << std::endl; exit(0);
00981 }
00982 v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
00983
00984 double xL1(cf*vlb[v1]); double xU1(cf*vub[v1]);
00985 double xL2(vlb[v2]); double xU2(vub[v2]);
00986 double xL3(vlb[v3]); double xU3(vub[v3]);
00987
00988 for (int ii = 0; ii < defcons_size; ii++) {
00989
00990 cutIndices [ii][0] = v1;
00991 cutIndices [ii][1] = v2;
00992 cutIndices [ii][2] = v3;
00993 cutIndices [ii][3] = v4;
00994 cutCoeff [ii][3] = 1.;
00995 }
00996
00997
00998 if(vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vub[v3]
00999 <= vub[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3] &&
01000 vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vub[v3]
01001 <= vlb[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3]) {
01002
01003 double theta1 = xL1*xL2*xU3 - xU1*xU2*xU3 - xL1*xL2*xL3 + xL1*xU2*xL3;
01004 double theta2 = xU1*xL2*xU3 - xL1*xL2*xL3 - xU1*xU2*xU3 + xU1*xU2*xL3;
01005
01006 cutCoeff [0][0] = -xU2*xU3; cutCoeff [0][1] = -xL1*xU3; cutCoeff [0][2] = -xL1*xU2; bnd[0] = - 2.*xL1*xU2*xU3;
01007 cutCoeff [1][0] = -xL2*xL3; cutCoeff [1][1] = -xU1*xL3; cutCoeff [1][2] = -xU1*xL2; bnd[1] = - 2.*xU1*xL2*xL3;
01008 cutCoeff [2][0] = -xL2*xU3; cutCoeff [2][1] = -xU1*xU3; cutCoeff [2][2] = -xL1*xL2; bnd[2] = - xU1*xL2*xU3 - xL1*xL2*xU3;
01009 cutCoeff [3][0] = -xU2*xL3; cutCoeff [3][1] = -xL1*xL3; cutCoeff [3][2] = -xU1*xU2; bnd[3] = - xU1*xU2*xL3 - xL1*xU2*xL3;
01010 cutCoeff [4][0] = -(theta1/(xL1-xU1)); cutCoeff [4][1] = -xL1*xL3; cutCoeff [4][2] = -xL1*xL2;
01011 bnd[4] = (-(theta1*xU1)/(xL1-xU1)) - xL1*xL2*xU3 - xL1*xU2*xL3 + xU1*xU2*xU3;
01012 cutCoeff [5][0] = -(theta2/(xU1-xL1)); cutCoeff [5][1] = -xU1*xU3; cutCoeff [5][2] = -xU1*xU2;
01013 bnd[5] = (-(theta2*xL1)/(xU1-xL1)) - xU1*xL2*xU3 - xU1*xU2*xL3 + xL1*xL2*xL3;
01014
01015 } else if(vub[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3]
01016 <= vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vub[v3] &&
01017 vub[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3]
01018 <= vlb[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3]) {
01019
01020 double theta1 = xU1*xU2*xU3 - xL1*xL2*xU3 - xU1*xU2*xL3 + xL1*xU2*xL3;
01021 double theta2 = xL1*xL2*xL3 - xU1*xU2*xL3 - xL1*xL2*xU3 + xU1*xL2*xU3;
01022
01023 cutCoeff [0][0] = -xU2*xU3; cutCoeff [0][1] = -xL1*xU3; cutCoeff [0][2] = -xL1*xU2; bnd[0] = - 2.*xL1*xU2*xU3;
01024 cutCoeff [1][0] = -xL2*xL3; cutCoeff [1][1] = -xU1*xL3; cutCoeff [1][2] = -xU1*xL2; bnd[1] = - 2.*xU1*xL2*xL3;
01025 cutCoeff [2][0] = -xL2*xU3; cutCoeff [2][1] = -xU1*xU3; cutCoeff [2][2] = -xU1*xU2; bnd[2] = - xU1*xL2*xU3 - xU1*xU2*xU3;
01026 cutCoeff [3][0] = -xU2*xL3; cutCoeff [3][1] = -xL1*xL3; cutCoeff [3][2] = -xL1*xL2; bnd[3] = - xL1*xL2*xL3 - xL1*xU2*xL3;
01027 cutCoeff [4][0] = -xU2*xL3; cutCoeff [4][1] = -(theta1/(xU2-xL2)); cutCoeff [4][2] = -xU1*xU2;
01028 bnd[4] = (-(theta1*xL2)/(xU2-xL2)) - xU1*xU2*xU3 - xL1*xU2*xL3 + xL1*xL2*xU3;
01029 cutCoeff [5][0] = -xL2*xU3; cutCoeff [5][1] = -(theta2/(xL2-xU2)); cutCoeff [5][2] = -xL1*xL2;
01030 bnd[5] = (-(theta2*xU2)/(xL2-xU2)) - xL1*xL2*xL3 - xU1*xL2*xU3 + xU1*xU2*xL3;
01031 }
01032
01033 cutCoeff [6][0] = -xL2*xL3; cutCoeff [6][1] = -xL1*xU3; cutCoeff [6][2] = -xL1*xL2; bnd[6] = - xL1*xL2*xU3 - xL1*xL2*xL3;
01034 cutCoeff [7][0] = -xL2*xL3; cutCoeff [7][1] = -xL1*xL3; cutCoeff [7][2] = -xL1*xU2; bnd[7] = - xL1*xU2*xL3 - xL1*xL2*xL3;
01035 cutCoeff [8][0] = -xL2*xU3; cutCoeff [8][1] = -xL1*xU3; cutCoeff [8][2] = -xU1*xL2; bnd[8] = - xL1*xL2*xU3 - xU1*xL2*xU3;
01036 cutCoeff [9][0] = -xU2*xU3; cutCoeff [9][1] = -xU1*xU3; cutCoeff [9][2] = -xU1*xL2; bnd[9] = - xU1*xU2*xU3 - xU1*xL2*xU3;
01037 cutCoeff [10][0] = -xU2*xU3; cutCoeff [10][1] = -xU1*xL3; cutCoeff [10][2] = -xU1*xU2; bnd[10] = - xU1*xU2*xU3 - xU1*xU2*xL3;
01038 cutCoeff [11][0] = -xU2*xL3; cutCoeff [11][1] = -xU1*xL3; cutCoeff [11][2] = -xL1*xU2; bnd[11] = - xL1*xU2*xL3 - xU1*xU2*xL3;
01039
01040 }
01041
01042
01043
01044
01045 if(flag == 10) {
01046 #ifdef DEBUG
01047 std::cout << " -- case 10 --" << std::endl;
01048 #endif
01049
01050 defcons_size = 12;
01051 prepareVectors (defcons_size);
01052
01053 ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3;
01054 permutation3(ind,ibnd);
01055 int i, flagg=0, idx=0;
01056 i = 0;
01057 while(i < 6 && flagg == 0) {
01058 if(vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
01059 >= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]] &&
01060 vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
01061 >= vub[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]])
01062 {
01063 idx = i;
01064 flagg = 1;
01065 }
01066 i++;
01067 }
01068 if (flagg==0) {
01069 std::cout << "ERROR!!!" << std::endl; exit(0);
01070 }
01071 v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
01072
01073 double xL1(cf*vlb[v1]); double xU1(cf*vub[v1]);
01074 double xL2(vlb[v2]); double xU2(vub[v2]);
01075 double xL3(vlb[v3]); double xU3(vub[v3]);
01076
01077 for(int ii = 0; ii < defcons_size; ii++) {
01078
01079 cutIndices [ii][0] = v1;
01080 cutIndices [ii][1] = v2;
01081 cutIndices [ii][2] = v3;
01082 cutIndices [ii][3] = v4;
01083 cutCoeff [ii][3] = 1.;
01084 }
01085
01086
01087 cutCoeff [0][0] = -xL2*xL3; cutCoeff [0][1] = -xU1*xL3; cutCoeff [0][2] = -xU1*xU2; bnd[0] = - xU1*xU2*xL3 - xU1*xL2*xL3;
01088 cutCoeff [1][0] = -xU2*xU3; cutCoeff [1][1] = -xL1*xL3; cutCoeff [1][2] = -xL1*xU2; bnd[1] = - xL1*xU2*xU3 - xL1*xU2*xL3;
01089 cutCoeff [2][0] = -xU2*xL3; cutCoeff [2][1] = -xL1*xL3; cutCoeff [2][2] = -xU1*xU2; bnd[2] = - xU1*xU2*xL3 - xL1*xU2*xL3;
01090 cutCoeff [3][0] = -xU2*xU3; cutCoeff [3][1] = -xL1*xU3; cutCoeff [3][2] = -xL1*xL2; bnd[3] = - xL1*xU2*xU3 - xL1*xL2*xU3;
01091 cutCoeff [4][0] = -xL2*xU3; cutCoeff [4][1] = -xU1*xU3; cutCoeff [4][2] = -xL1*xL2; bnd[4] = - xU1*xL2*xU3 - xL1*xL2*xU3;
01092 cutCoeff [5][0] = -xL2*xL3; cutCoeff [5][1] = -xU1*xU3; cutCoeff [5][2] = -xU1*xL2; bnd[5] = - xU1*xL2*xU3 - xU1*xL2*xL3;
01093
01094
01095
01096
01097
01098
01099 double theta1c = xL1*xU2*xL3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xL2*xU3;
01100 double theta2c = xU1*xU2*xL3 - xL1*xU2*xU3 - xU1*xL2*xL3 + xU1*xL2*xU3;
01101
01102 cutCoeff [6][0] = -xL2*xL3; cutCoeff [6][1] = -xL1*xL3; cutCoeff [6][2] = -xL1*xL2; bnd[6] = - 2.*xL1*xL2*xL3;
01103 cutCoeff [7][0] = -xU2*xU3; cutCoeff [7][1] = -xU1*xU3; cutCoeff [7][2] = -xU1*xU2; bnd[7] = - 2.*xU1*xU2*xU3;
01104 cutCoeff [8][0] = -xL2*xU3; cutCoeff [8][1] = -xL1*xU3; cutCoeff [8][2] = -xU1*xL2; bnd[8] = - xL1*xL2*xU3 - xU1*xL2*xU3;
01105 cutCoeff [9][0] = -xU2*xL3; cutCoeff [9][1] = -xU1*xL3; cutCoeff [9][2] = -xL1*xU2; bnd[9] = - xL1*xU2*xL3 - xU1*xU2*xL3;
01106 cutCoeff [10][0] = -(theta1c/(xL1-xU1)); cutCoeff [10][1] = -xL1*xU3; cutCoeff [10][2] = -xL1*xU2;
01107 bnd[10] = (-(theta1c*xU1)/(xL1-xU1)) - xL1*xU2*xL3 - xL1*xL2*xU3 + xU1*xL2*xL3;
01108 cutCoeff [11][0] = -(theta2c/(xU1-xL1)); cutCoeff [11][1] = -xU1*xL3; cutCoeff [11][2] = -xU1*xL2;
01109 bnd[11] = (-(theta2c*xL1)/(xU1-xL1)) - xU1*xU2*xL3 - xU1*xL2*xU3 + xL1*xU2*xU3;
01110
01111 }
01112
01113
01114
01115
01116 for (int ii = 0; ii < defcons_size; ii++) {
01117
01118
01119
01120 if (ii < defcons_size/2) {
01121
01122 if (cf > 0) {cutLb[ii] = bnd[ii]; cutUb[ii] = COUENNE_INFINITY;}
01123 if (cf < 0) {cutLb[ii] = -COUENNE_INFINITY; cutUb[ii] = bnd[ii];}
01124
01125 } else {
01126
01127 if (cf > 0) {cutLb[ii] = -COUENNE_INFINITY; cutUb[ii] = bnd[ii];}
01128 if (cf < 0) {cutLb[ii] = bnd[ii]; cutUb[ii] = COUENNE_INFINITY;}
01129 }
01130 #ifdef DEBUG
01131 std::cout << ii << ") cutLb =" << cutLb[ii] << " " << "cutUb = " << cutUb[ii] << std::endl;
01132 #endif
01133 }
01134
01135 for (int i=0; i<6; i++)
01136 delete [] ind [i];
01137
01138 delete [] ibnd;
01139 delete [] bnd;
01140 delete [] ind;
01141 }
01142
01143
01144
01145 void exprTrilinear::generateCuts (expression *w,
01146 OsiCuts &cs, const CouenneCutGenerator *cg,
01147 t_chg_bounds *chg, int wind,
01148 CouNumber lbw, CouNumber ubw) {
01149
01150 expression **args = w -> Image () -> ArgList ();
01151
01152 int *varInd = new int [4];
01153
01154 for (int i=0; i<3; i++)
01155 varInd [i] = args [i] -> Index ();
01156
01157 varInd [3] = w -> Index ();
01158
01159 std::vector <std::vector <int> > cutIndices;
01160 std::vector <std::vector <double> > cutCoeff;
01161 std::vector <double> cutLb, cutUb;
01162
01163 TriLinCuts (cg -> Problem () -> Lb (),
01164 cg -> Problem () -> Ub (),
01165 varInd,
01166 cutIndices, cutCoeff,
01167 cutLb, cutUb);
01168
01169
01170 assert (cutIndices.size () == cutCoeff.size () &&
01171 cutIndices.size () == cutLb.size () &&
01172 cutIndices.size () == cutUb.size ());
01173
01174
01175
01176 for (int i = (int) cutIndices.size (); i--;) {
01177
01178 int
01179 size = (int) cutIndices [i].size (),
01180 *ind = new int [size];
01181
01182 double *coe = new double [size];
01183
01184 std::copy (cutIndices [i].begin (), cutIndices [i].end (), ind);
01185 std::copy (cutCoeff [i].begin (), cutCoeff [i].end (), coe);
01186
01187 OsiRowCut cut (cutLb [i], cutUb [i], 4, 4, ind, coe);
01188
01189
01190 delete [] ind;
01191 delete [] coe;
01192
01193 if (cg -> Problem () -> bestSol ()) {
01194
01195
01196
01197
01198 double *sol = cg -> Problem () -> bestSol ();
01199 const double
01200 *lb = cg -> Problem () -> Lb (),
01201 *ub = cg -> Problem () -> Ub ();
01202
01203 int nVars = cg -> Problem () -> nVars ();
01204
01205 bool optIn = true;
01206
01207 for (int i=0; i< nVars; i++)
01208 if ((sol [i] < lb [i] - COUENNE_EPS) ||
01209 (sol [i] > ub [i] + COUENNE_EPS)) {
01210 optIn = false;
01211 break;
01212 }
01213
01214 if (optIn) {
01215
01216 for (unsigned int i=0; i<cutIndices.size (); i++) {
01217
01218 double chs = 0.;
01219
01220 for (unsigned int j=0; j<cutIndices[i].size(); j++)
01221 chs += cutCoeff [i] [j] * sol [cutIndices [i] [j]];
01222
01223 if ((chs < cutLb [i] - COUENNE_EPS) ||
01224 (chs > cutUb [i] + COUENNE_EPS)) {
01225
01226 printf ("cut %d violates optimum:\n", i);
01227
01228 if (cutLb [i] > -COUENNE_INFINITY) printf ("%g <= ", cutLb [i]);
01229 for (unsigned int j=0; j<cutIndices[i].size(); j++) printf ("%+g x%d ", cutCoeff [i] [j], cutIndices [i] [j]); printf ("\n = ");
01230 for (unsigned int j=0; j<cutIndices[i].size(); j++) printf ("%+g *%g ", cutCoeff [i] [j], sol [cutIndices [i] [j]]); printf ("\n = ");
01231 for (unsigned int j=0; j<cutIndices[i].size(); j++) printf ("%+g ", cutCoeff [i] [j] * sol [cutIndices [i] [j]]); printf (" = %g", chs);
01232 if (cutUb [i] < COUENNE_INFINITY) printf (" <= %g", cutUb [i]);
01233 printf ("\n");
01234
01235 }
01236 }
01237 }
01238 }
01239
01240 cs.insert (cut);
01241 }
01242
01243 delete [] varInd;
01244 }