set i starting / 1*3 / j ending / 1*4 / k set to derive iterations / 1*64 /; table alternative(i,j) distance between i and j 1 2 3 4 1 10 12 0 0 2 0 0 8 10 3 20 0 0 5; alias (i,ii,iii),(j,jj,jjj),(k,kk); variable objective some objective variable; binary variable y(i,j) binary variable; parameter report1(k,i,j) solution report when i = 1 report2(k,ii,jj) solution report when i = 2 report3(k,iii,jjj) solution report when i = 3 sol(k) feasible solution record ; report1(k,i,j) = 0; report2(k,ii,jj) = 0; report3(k,iii,jjj) = 0; equation c1 objective function c2 logical constraint only one alternative can be selected for each i, cut(kk,i,j,ii,jj,iii,jjj) cutting constraints; c1.. objective =e= sum((i,j),alternative(i,j) *y(i,j)); c2(i).. sum(j, y(i,j))=e=1; cut(kk,i,j,ii,jj,iii,jjj)$(report1(kk,i,j) and report2(kk,ii,jj)and report3(kk,iii,jjj)).. y(i,j) + y(ii,jj) + y(iii,jjj) =l= 2; model mymodel /all/; option limrow = 100; option optcr = 0; solve mymodel maximize objective using mip; loop (k, report1(k,i,j) $(ord(i) eq 1)=y.l(i,j); report2(k,ii,jj)$(ord(ii) eq 2)=y.l(ii,jj); report3(k,iii,jjj)$(ord(iii) eq 3)=y.l(iii,jjj); sol(k) = objective.l; solve mymodel maximize objective using mip; display report1, report2,report3, sol; abort$(mymodel.modelstat <> 1) " model not normally completed", mymodel.modelstat; ); display report1, report2,report3, sol;