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

/* 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 1

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 2

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 3

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 4