/* project.c version trim 2: alternative use of - - PDF document

Ralf Östermark 1

/* project.c version trim 2: alternative use of MILP_MACHINE/SIMPLEX/LP_SOLVE to solve the Trim Loss problem of Pörn, Harjunkoski & Westerlund This file contains the following problem specific functions: (1) evaluator() (2) accelerator() (3) gradient() (4) hessian() (5) pre_processor() (6) post_processor() Testing the powerful MILP-machine of GHA */ /* Global Definitions */ #define MILP_MACHINE 1 #define SIMPLEX 2 #define LP_SOLVE 3 #include "project.h" #include "rpa_proj.h" extern MINLP_ptr *MINLP_PROB; extern GHA_ptr *GHA_PROB; #include "rpa_proj.c" #define FREE_PROJECT free_dmatrix(A,0,0);free_dvector(b,0);free_dvector(c,0);\ free_dvector(x,0);free_dvector(RC,0);free_ivector(IRTYPE,0);\ free_ivector(INTEGERS,0);free_dvector(LB,0);free_dvector(UB,0);free_ive ctor(base,0); void resize_SHAREX(int m,int n,DVECTOR b,DVECTOR c,DVECTOR x,IVECTOR IRTYPE, DVECTOR LB,DVECTOR UB,DMATRIX A,MINLP_ptr *MINLP_PROB) { int i; MINLP_PROB->H = resize_dmatrix(MINLP_PROB->H,0,(long)(n- 1),0,(long)(n-1)); MINLP_PROB->A = resize_dmatrix(MINLP_PROB->A,0,(long)(m- 1),0,(long)(n-1)); MINLP_PROB->b = resize_dvector(MINLP_PROB->b,0,(long)(m-1)); MINLP_PROB->IRTYPE = resize_ivector(MINLP_PROB->IRTYPE,0,(long)(m-1)); MINLP_PROB->c = resize_dvector(MINLP_PROB->c,0,(long)(n-1)); MINLP_PROB->x = resize_dvector(MINLP_PROB->x,0,(long)(n-1)); MINLP_PROB->x_best = resize_dvector(MINLP_PROB->x_best,0,(long)(n-1)); MINLP_PROB->LB = resize_dvector(MINLP_PROB->LB,0,(long)(n-1)); MINLP_PROB->UB = resize_dvector(MINLP_PROB->UB,0,(long)(n-1)); MINLP_PROB->LB_best = resize_dvector(MINLP_PROB->LB_best,0,(long)(n- 1)); MINLP_PROB->UB_best = resize_dvector(MINLP_PROB->UB_best,0,(long)(n- 1));

Ralf Östermark 2

memcpy(MINLP_PROB->b,b,m*sizeof(double)); memcpy(MINLP_PROB->c,c,n*sizeof(double)); memcpy(MINLP_PROB->x,x,n*sizeof(double)); memcpy(MINLP_PROB->IRTYPE,IRTYPE,m*sizeof(int)); memcpy(MINLP_PROB->LB,LB,n*sizeof(double)); memcpy(MINLP_PROB->UB,UB,n*sizeof(double)); for(i=0;i<m;i++) { memcpy(MINLP_PROB->A[i],A[i],n*sizeof(double)); } /* end of i-loop */ } /* end of resize_SHAREX) */ void trim_debug(int phase,int m,int n,DVECTOR b,DVECTOR c, DVECTOR LB,DVECTOR UB,IVECTOR IRTYPE,DMATRIX A) { FILE *fid; int i,j; if(phase EQ 1) fid=fopen("trim_debug.out","w"); if(phase EQ 2) fid=fopen("trim_debug.out","a"); fprintf(fid,"phase %d:\r\n",phase); fprintf(fid,"b:\r\n");fflush(fid); for(i=0;i<m;i++) fprintf(fid,"%7.2f ",b[i]);fflush(fid);fprintf(fid,"\r\n"); fprintf(fid,"c:\r\n"); for(i=0;i<n;i++) fprintf(fid,"%7.2f ",c[i]);fflush(fid);fprintf(fid,"\r\n"); fprintf(fid,"LB:\r\n"); for(i=0;i<n;i++) fprintf(fid,"%7.2f ",LB[i]);fflush(fid);fprintf(fid,"\r\n"); fprintf(fid,"UB:\r\n"); for(i=0;i<n;i++) fprintf(fid,"%7.2e ",UB[i]);fflush(fid);fprintf(fid,"\r\n"); fprintf(fid,"IRTYPE:\r\n");fflush(fid); for(i=0;i<m;i++) fprintf(fid,"%d ",IRTYPE[i]);fprintf(fid,"\r\n");fflush(fid); fprintf(fid,"A:\r\n");fflush(fid); for(i=0;i<m;i++) { for(j=0;j<n;j++) fprintf(fid,"%7.2f ",A[i][j]);fprintf(fid,"\r\n");fflush(fid); } fclose(fid); } /* end of trim_debug() */ void print_model(int m,int n,double fx,DVECTOR w,DVECTOR x,DVECTOR RC,DMATRIX A, DVECTOR b,DVECTOR c,IVECTOR IRTYPE) { int i,j; SETTINGS_ptr *TASK = GHA_PROB->TASK; printf("\r\nvector POP: "); for(i=0;i < *(TASK->NVAR);i++) printf("%3.1f ",w[i]); printf("\r\nobjective = %f at vector x: ",fx); for(i=0;i<n+m;i++) printf("%3.1f ",x[i+1]); printf("\r\nvector RC: "); for(i=0;i<n+m;i++) printf("%3.1f ",RC[i+1]); printf("\r\nvector c: "); for(i=0;i<n;i++) printf("%3.1f ",c[i]); printf("\r\nvector b: ");

Ralf Östermark 3

for(i=0;i<m;i++) printf("%3.1f ",b[i]); printf("\r\nmatrix A,vector b and IRTYPE:\r\n"); for(i=0;i<m;i++) { for(j=0;j<n;j++) { printf("%3.0f ",A[i][j]); } printf("%4.3f %1d \r\n",b[i],IRTYPE[i]); } } /* end of print_model */ /* end of Global Definitions */ void evaluator(GHA_ptr *GHA_PROB,DVECTOR w,double Penalty,double *gf,double *F,double *Dev) { /* The structure of w: {x[i],i,j=0,...,4;yj,zj} */ int ACCELERATE,i,GA_t=GHA_PROB->INT_STATUS[3]; DVECTOR x1;DVECTOR x2;DVECTOR x3;DVECTOR x4; DVECTOR z;DVECTOR y; double n_o[4] = {15,28,21,30}; double width[4] = {290,315,350,455}; double BEST_gf = GHA_PROB->REAL_STATUS[1]; double BEST_Dev = GHA_PROB->REAL_STATUS[3]; ACCELERATE = GHA_PROB->INT_STATUS[0]; GHA_PROB->REAL_STATUS[4] += 1.0; /* Add one for calling evaluator() */ z=&(w[0]);x1=&(w[4]);x2=&(w[8]);x3=&(w[12]);x4=&(w[16]);y=&(w[20]); *F = 0.0;*gf = 0.0;*Dev = 0.0; if(ACCELERATE EQ TRUE) {} for(i=0;i<4;i++) {*F += z[i]+0.1*(i+1)*y[i];} for(i=0;i<4;i++) { *Dev += max(0.0,width[0]*x1[i]+width[1]*x2[i]+width[2]*x3[i]+width[3]*x4[i]- 1850*y[i]); *Dev += max(0.0,1750*y[i]- (width[0]*x1[i]+width[1]*x2[i]+width[2]*x3[i]+width[3]*x4[i])); *Dev += max(0.0,x1[i]+x2[i]+x3[i]+x4[i]-5*y[i]); *Dev += max(0.0,y[i]-z[i]); *Dev += max(0.0,z[i]-30*y[i]); } *Dev += max(0.0,n_o[0]-z[0]*x1[0]-z[1]*x1[1]-z[2]*x1[2]-z[3]*x1[3]); *Dev += max(0.0,n_o[1]-z[0]*x2[0]-z[1]*x2[1]-z[2]*x2[2]-z[3]*x2[3]); *Dev += max(0.0,n_o[2]-z[0]*x3[0]-z[1]*x3[1]-z[2]*x3[2]-z[3]*x3[3]); *Dev += max(0.0,n_o[3]-z[0]*x4[0]-z[1]*x4[1]-z[2]*x4[2]-z[3]*x4[3]); *gf = *F+Penalty*(*Dev); if(((BEST_Dev LE ZEROLIM) AND (BEST_gf LE 19.61)) OR ((GA_t EQ 0) AND (*gf LE 19.61))) STOP_FLIP=TRUE; } void accelerator(GHA_ptr *GHA_PROB,int FIX_Flip,int ROW,double Penalty,DVECTOR LOWER, DVECTOR UPPER,DVECTOR x_IN,DVECTOR x_OUT) { /* IRTYPEs: = -> EQ_rh (0); < -> LE_rh (1); > -> GE_rh (2)

Ralf Östermark 4

NOTE: DEBUG_MODE=0 -> no debug; 1 -> prGHA_PROB->INT_solution; 2 -> full analysis structure of POP: {zj,x[i],i,j=1,...,4,yj} structure of x (simplex formulation 1): {zj,yj,dj given x[i],i,j=1,...,4} structure of x (simplex formulation 2): {x[i],i,j=1,...,4,dj given zj,yj} */ int SOLVER,i,j,m,n,IVARS,DEBUG_MODE=FALSE,LOOP=1,MAXLOOP=1,SIGN=1,MILP_STAT US,s_OFFSET; double simplex_1=1.0,simplex_2=2.0,DF=0.0; DMATRIX A;DVECTOR b;DVECTOR c;DVECTOR x;DVECTOR RC;IVECTOR IRTYPE;IVECTOR base; DVECTOR z;DVECTOR y;DVECTOR x1;DVECTOR x2;DVECTOR x3;DVECTOR x4; IVECTOR INTEGERS;DVECTOR LB;DVECTOR UB; double n_o[4] = {15,28,21,30}; double width[4] = {290,315,350,455}; SETTINGS_ptr *TASK = GHA_PROB->TASK; GHA_PROB->REAL_STATUS[5] += 1.0; /* Add one for calling accelerator() */ SOLVER=MILP_MACHINE;if(SOLVER EQ SIMPLEX) SIGN=-1; z=&(x_IN[0]); for(j=0;j<4;j++) { x_OUT[j] = min(z[j],30); x_OUT[j+20] = (z[j]>0 ? 1.0:0.0); } for(j=4;j < *(TASK->NVAR)-4;j++) x_OUT[j] = x_IN[j]; if(*(TASK->GA_RUNS) EQ 0) DEBUG_MODE=TRUE; if((FIX_Flip > -2) OR (DEBUG_MODE EQ TRUE)) { while((dabs(DF-(simplex_1+simplex_2)) > EPS) AND (LOOP LE MAXLOOP)) { LOOP++; DF=simplex_1+simplex_2; m=24;n=8+8+16;IVARS=8; /* we introduce 8+16 slacks for the GE_rh+LE_rh constraints */ s_OFFSET = 8; /* the slack variable section begins here */ A=dmatrixCALLOC(0,m-1,0,n-1);b=dvectorCALLOC(0,m- 1);c=dvectorCALLOC(0,n-1); x=dvectorCALLOC(0,m+n);RC=dvectorCALLOC(0,m+n);LB=dvectorCALLOC(0,n- 1);UB=dvectorCALLOC(0,n-1); IRTYPE=ivectorCALLOC(0,m-1);INTEGERS=ivectorCALLOC(0,IVARS- 1);base=ivectorCALLOC(0,m-1); for(i=0;i<4;i++) { A[i][4+i]=1850;A[4+i][4+i]=1750; A[i][8+i]=A[8+i][4+i]=A[20+i][4+i]=1.0; A[8+i][i] = -1;A[12+i][4+i] = -30; } for(j=0;j<IVARS;j++) {INTEGERS[j]=j+1;} for(j=8;j<16;j++) {b[j] = 0.0;} for(j=16;j<20;j++) {b[j] = n_o[j-16];} for(j=20;j<24;j++) {b[j] = 1.0;} for(j=0;j<4;j++) {IRTYPE[j]=IRTYPE[j+16]=GE_rh;} for(j=4;j<16;j++) {IRTYPE[j]=LE_rh;} for(j=20;j<24;j++) {IRTYPE[j]=LE_rh;}

Ralf Östermark 5

for(j=0;j<4;j++) {c[j] = 1;c[j+4] = 0.1*(j+1);} for(j=s_OFFSET;j<n;j++) {c[j] = 5;} for(j=0;j<n;j++) {c[j] = SIGN*c[j];} for(j=0;j<4;j++) { UB[j]=UPPER[j];LB[j]=LOWER[j];UB[j+4]=UPPER[j+20];LB[j+4]=LOWER[j+20]; } for(j=8;j<n;j++) {UB[j]=100000.0;LB[j]=0.0;} x1=&(x_IN[4]);x2=&(x_IN[8]);x3=&(x_IN[12]);x4=&(x_IN[16]); for(i=0;i<4;i++) { b[i] = width[0]*x1[i]+width[1]*x2[i]+width[2]*x3[i]+width[3]*x4[i]; b[i+4] = width[0]*x1[i]+width[1]*x2[i]+width[2]*x3[i]+width[3]*x4[i]; } for(j=0;j<4;j++) {A[16][j]=x1[j];A[17][j]=x2[j];A[18][j]=x3[j];A[19][j]=x4[j];} /* add slacks */ for(i=0;i<m;i++) { if(IRTYPE[i] EQ LE_rh) A[i][s_OFFSET+i] = -1.0; if(IRTYPE[i] EQ GE_rh) A[i][s_OFFSET+i] = 1.0; } /* end of i-loop */ switch(SOLVER) { case MILP_MACHINE: MINLP_PROB->INITIALIZED = FALSE; MINLP_PROB->TREE_COUNTER = 0.0; MINLP_PROB->ACTIVE_TREE = 0.0; MINLP_PROB->m = m; MINLP_PROB->n = n; MINLP_PROB->n_i = IVARS; MINLP_PROB->n_x = n-IVARS; MINLP_PROB->n_d = 0; resize_SHAREX(m,n,b,c,x,IRTYPE,LB,UB,A,MINLP_PROB); if(DEBUG_MODE EQ TRUE) trim_debug(1,m,n,MINLP_PROB->b,MINLP_PROB- >c,MINLP_PROB->LB, MINLP_PROB->UB,MINLP_PROB- >IRTYPE,MINLP_PROB->A); MILP_STATUS = milp_caller(MINLP_PROB); if(MINLP_PROB->INITIALIZED EQ TRUE) { for(j=0;j<4;j++) { x_OUT[j] = MINLP_PROB->x[j]; /* extract z-vector */ x_OUT[j+20] = MINLP_PROB->x[j+4]; /* extract y-vector */ } } /* end of if() */ break; case LP_SOLVE: for(j=0;j<n;j++) { if(LB[j]>UB[j]) { printf("accelerator(): Lower-Upper conflict in variable %d. %20.19f > %20.19f\r\n",j,LB[j],UB[j]); exit(-1); } }

Ralf Östermark 6

simplex_1=lp_solve(m,n,A,b,c,LB,UB,IRTYPE,IVARS,INTEGERS,x,RC,DEBUG_MOD E,&MILP_STATUS,base); if(MILP_STATUS EQ 0) { for(j=0;j<4;j++) { x_OUT[j] = x[j]; /* extract z-vector */ x_OUT[j+20] = x[j+4]; /* extract y-vector */ } } break; case SIMPLEX: simplex_1=simplex(m,n,A,b,c,IRTYPE,x,RC); for(j=0;j<4;j++) { x_OUT[j] = x[j+1]; /* extract z-vector */ x_OUT[j+20] = x[j+5]; /* extract y-vector */ } if(DEBUG_MODE EQ TRUE) {print_model(m,n,simplex_1,x_OUT,x,RC,A,b,c,IRTYPE);} break; } /* end of switch() */ FREE_PROJECT m=16;n=16+8+8;IVARS=16; /* we introduce 8+8 slacks for the GE_rh+LE_rh constraints */ s_OFFSET = 16; /* the slack variable section begins here */ A=dmatrixCALLOC(0,m-1,0,n-1);b=dvectorCALLOC(0,m- 1);c=dvectorCALLOC(0,n-1); x=dvectorCALLOC(0,m+n);RC=dvectorCALLOC(0,m+n);IRTYPE=ivectorCALLOC(0,m

• 1);base=ivectorCALLOC(0,m-1);

INTEGERS=ivectorCALLOC(0,IVARS-1);LB=dvectorCALLOC(0,n- 1);UB=dvectorCALLOC(0,n-1); for(j=0;j<s_OFFSET;j++) c[j] = 0.0; for(j=s_OFFSET;j<n;j++) c[j] = SIGN; for(j=0;j<4;j++) { IRTYPE[j]=LE_rh;IRTYPE[j+4]=GE_rh;IRTYPE[j+8]=LE_rh;IRTYPE[j+12]=GE_rh; } for(i=0;i<4;i++) { for(j=0;j<4;j++) { A[i][4*j+i] = A[4+i][4*j+i] = width[j]; A[8+i][4*j+i]=1; } } /* end of i-loop */ for(j=0;j<IVARS;j++) INTEGERS[j]=j+1; for(j=0;j<16;j++) {UB[j]=UPPER[j+4];LB[j]=LOWER[j+4];} for(j=16;j<n;j++) {UB[j]=100000.0;LB[j]=0.0;} z =&(x_OUT[0]);y=&(x_OUT[20]); for(i=0;i<4;i++) { for(j=0;j<4;j++) {A[12+i][4*i+j]=x_OUT[j];} } for(j=0;j<4;j++) { b[j]=1850*y[j];b[j+4]=1750*y[j];b[j+8]=5*y[j]; b[j+12]=n_o[j]; } /* add slacks */

Ralf Östermark 7

for(i=0;i<m;i++) { if(IRTYPE[i] EQ LE_rh) A[i][s_OFFSET+i] = -1.0; if(IRTYPE[i] EQ GE_rh) A[i][s_OFFSET+i] = 1.0; } /* end of i-loop */ switch(SOLVER) { case MILP_MACHINE: MINLP_PROB->INITIALIZED = FALSE; MINLP_PROB->TREE_COUNTER = 0.0; MINLP_PROB->ACTIVE_TREE = 0.0; MINLP_PROB->m = m; MINLP_PROB->n = n; MINLP_PROB->n_i = IVARS; MINLP_PROB->n_x = n-IVARS; MINLP_PROB->n_d = 0; resize_SHAREX(m,n,b,c,x,IRTYPE,LB,UB,A,MINLP_PROB); if(DEBUG_MODE EQ TRUE) trim_debug(2,m,n,MINLP_PROB->b,MINLP_PROB- >c,MINLP_PROB->LB, MINLP_PROB->UB,MINLP_PROB- >IRTYPE,MINLP_PROB->A); MILP_STATUS = milp_caller(MINLP_PROB); if(MINLP_PROB->INITIALIZED EQ TRUE) { for(j=0;j<s_OFFSET;j++) {x_OUT[j+4] = MINLP_PROB->x[j];} /* extract x-vector */ } break; case LP_SOLVE: simplex_2=lp_solve(m,n,A,b,c,LB,UB,IRTYPE,IVARS,INTEGERS,x,RC,DEBUG_MOD E,&MILP_STATUS,base); if(MILP_STATUS EQ 0) { for(j=0;j<s_OFFSET;j++) {x_OUT[j+4] = x[j];} /* extract x-vector */ } break; case SIMPLEX: simplex_2=simplex(m,n,A,b,c,IRTYPE,x,RC); for(j=0;j<s_OFFSET;j++) {x_OUT[j+4] = x[j+1];} /* extract x-vector */ if(DEBUG_MODE EQ TRUE) {print_model(m,n,simplex_2,x_OUT,x,RC,A,b,c,IRTYPE);} break; } /* end of switch() */ memcpy(x_IN,x_OUT,*(TASK->NVAR)*sizeof(double)); FREE_PROJECT } /* end of while{} */ if(DEBUG_MODE EQ FALSE) { bound_POP(GHA_PROB,x_OUT,0,*(TASK->NVAR),LOWER,UPPER); for(j=0;j<*(TASK->INTEGERS)+(*(TASK->FIXED_INTEGERS));j++) { x_OUT[j]=round_integer(x_OUT[j]); } } /* end of if(DEBUG_MODE) */ if(ROW GE 0) GHA_PROB->POP->FLAGS[ROW] = FALSE; } /* end of if(FIX_Flip) */ } /* end of accelerator() */

Ralf Östermark 8

void gradient(GHA_ptr *GHA_PROB,DVECTOR w,int n_w,double Penalty,DVECTOR d,int n_d,IVECTOR pos) { /* NOTE: FIXED_h=1,VARIABLE_h=2 */ int h_TYPE=2; num_gradient(GHA_PROB,h_TYPE,w,n_w,Penalty,d,n_d,pos); } /* end of gradient */ void hessian(GHA_ptr *GHA_PROB,DVECTOR w,int n_w,double Penalty,DMATRIX G,int n_G,IVECTOR pos) { /* NOTE: FIXED_h=1,VARIABLE_h=2 */ int h_TYPE=2; num_hessian(GHA_PROB,h_TYPE,w,n_w,Penalty,G,n_G,pos); } /* end of hessian() */ void pre_processor(GHA_ptr *GHA_PROB,struct SOLUTION *START) { int PHASE = GHA_PROB->INT_STATUS[17]; IVECTOR ALLOCATION_SWITCH = &(GHA_PROB->INT_PROTOCOL[30]); SETTINGS_ptr *TASK = GHA_PROB->TASK; ARCHITECTURE_ptr *MESHINFO = GHA_PROB->MESHINFO; if((PHASE EQ 1) AND (*(TASK->ACTIVE_SYSTEM) EQ 0) AND (*ALLOCATION_SWITCH EQ FALSE)) { *ALLOCATION_SWITCH = TRUE; /* initlp_();*/ } /* end of if() */ if((MESHINFO->NODE_INDEX GE MESHINFO->MESH_SIZE/3) AND (MESHINFO->NODE_INDEX < 2*MESHINFO->MESH_SIZE/3)) { /* Adjust TASK for Node */ *(TASK->GET_HESSIAN)=TRUE; *(TASK->GA_SIBS)=32; *(TASK->LINEMETHOD)=2; *(TASK->T_INCR)=5; *(TASK->F_INCR)=2; *(TASK->T_SLOPE)=1; *(TASK->RND_STARTPOINTS)=1000; *(TASK->RND_ACCELERATOR)=1; *(TASK->CROSSOVER)=3; } if(MESHINFO->NODE_INDEX GE 2*MESHINFO->MESH_SIZE/3) { *(TASK->GET_HESSIAN)=FALSE; *(TASK->GA_SIBS)=32; *(TASK->LINEMETHOD)=1; *(TASK->F_INCR)=2; *(TASK->T_INCR)=3; *(TASK->T_SLOPE)=1; *(TASK->RND_STARTPOINTS)=2000; *(TASK->RND_ACCELERATOR)=1; *(TASK->CROSSOVER)=4; *(TASK->MUTATION)=2; *(TASK->GET_RCOND)=FALSE; } } /* end of pre_processor */ void post_processor(GHA_ptr *GHA_PROB,struct SOLUTION *BEST) { double Penalty;

Ralf Östermark 9

int PHASE = GHA_PROB->INT_STATUS[17]; IVECTOR SUPER_FLAG = &(GHA_PROB->INT_STATUS[16]); CONTROL_ptr *GHA_SYS = GHA_PROB->GHA_SYS; Penalty = GHA_PROB->REAL_STATUS[0]; if((*SUPER_FLAG EQ TRUE) AND (PHASE EQ GHA_SYS->SYSTEM_CALLS)){ deallocate_MINLP(MINLP_PROB); } /* end of if() */ } /* end of post_processor() */