pso2.m

这是从matlab上down的源文件。上面显示调试通过

function [c,a] = tsp020(a)
iworth = find(a(:,3)>0);
if isempty(iworth), c=[1,1]; return; end
mfreight = mean( a(iworth,3) ) / 15;
igood=find( a(:,4)>a(1,4)/15 | a(:,3)>=mfreight );
a = a(igood,:); n=size(a,1);
ichk = 2:n;
d = abs( a(ichk,1)-a(1) + 1i*(a(ichk,2)-a(1,2)) );
iworth = find( a(ichk,3)>0  ); ifuel = find( a(ichk,4)>0  );
gas = a(ichk(ifuel),4); afuel = sum( gas );
inf1a = any( a(1,4) >= d(ifuel) );
inf1b = any( a(1,4) >= d(iworth) );
inf1  = (inf1a | inf1b);
inf2 = any( a(1,4)+afuel >= 2*d(iworth) );
if ~( inf1 & inf2 ) , c=[1 1]; return; end
if n < 6
   c = tsp20(a);
else
   c = robsolver(a);
end
c=igood(c);
function [p,a] = tsp20(a)
npts = size(a,1);
x = a(:,1); y = a(:,2); freight = a(:,3); fuel = a(:,4);
availableFreight = sum(freight);
M = Inf;
   D = M+zeros(npts);
   for count1=1:npts,
      for count2=1:count1-1,
         x1 = x(count1); y1 = y(count1); x2 = x(count2); y2 = y(count2);
         d=sqrt((x1-x2)^2+(y1-y2)^2);
         D(count1,count2)=d - fuel(count1);
         D(count2,count1)=d - fuel(count2);
      end;
   end;
nPrb = 1;
S(1).D = D;
optc = Inf; p = [1];
while nPrb > 0
   D = S(1).D;
   [iz, unfeas] = bghungar(-D);
   S(1) = []; nPrb = nPrb - 1;
   if ~unfeas
      [pnew, nArc, df] = cycle( iz, D );
      g = cumsum( df );
      c = g(nArc)^2;
      if c < optc
         if any( g>0 ) | nArc < npts
            pLast = pnew(nArc);
            iPts = 1:npts;
            iArc = 1; iNext = iz(1);
            m1 = 0;
            while m1 < nArc
               m1 = m1 + 1;
               C = D;
               jPrev = pLast; jArc = 1; jNext = iz(1);
               m2 = 0;
               while m2 < m1-1
                  m2 = m2 + 1;
                  jNon = iPts( iPts ~= jNext & iPts ~= jPrev );
                  C(jArc,jNon) = M;
                  jPrev = jArc; jArc = jNext; jNext = iz(jArc);
               end 
               C(iArc,iNext) = M;
               nPrb = nPrb + 1; S(nPrb).D = C;
               iArc = iNext; iNext = iz(iArc);
            end 
if 1
            while nArc > 1
               nArc = nArc - 1;
               pnew = pnew(1:nArc); c = g(nArc)^2;
               if ~any( g(1:nArc)<0 ) ...
               & 0.95 < sum( a(pnew,3) )/availableFreight
                  p=pnew; optc=c; break;
               end
            end
end 
         
         else
         
            p=pnew; optc=c;
         end
      end 
   end 
end 
p = [p,1];
function [cyc, nArc, df] = cycle( iz, D );
   nArc = 1;
   iArc = 1; iNext = iz(1);
   cyc(nArc) = 1; df(nArc) = D(iArc,iNext);
   while iNext ~= 1
      nArc = nArc + 1;
      cyc(nArc) = iNext;
      iArc = iNext; iNext = iz(iArc);
      df(nArc) = D(iArc,iNext);
   end 
function [iz, unfeas, D] = bghungar(C);
unfeas = 0;
[m,n] = size(C);
[r,ii] = max(C,[],1); C1 = ones(m,1)*r - C;
[c,jj] = min(C1,[],2); D = C1 - c*ones(1,n);
iz = zeros(m,1);
if any(isnan(D)), unfeas = 1; return; end
for j = 1:m
   kk = find( ~D(:,j) );
   while ~isempty(kk),
      i=kk(1); if ~iz(i), iz(i) = j; break; else, kk(1) = []; end
   end
end 
izPrev = zeros(1,m);
keepSets = 0;
ii = find( iz );
while length( ii ) < m
if ~keepSets
   keepSets = 0;
   
   rr = 1:m; cc = 1:n;
   
   
   
   zz = iz(ii);
   
   
   cc( zz ) = [];
end 
while 1
   if any(isnan(D)), unfeas = 2; break; end
jz = 0;
for i = rr
   kk = find( ~D(i,cc) );
   if ~isempty(kk), jz = i; kz = cc(kk(1)); break; end
end 
if jz
   if iz(jz)
      kz = iz(jz); cc = [cc,kz]; zz(zz==kz) = []; rr(rr==jz) = [];
   else
      jz0 = jz; kz0 = kz; iz0 = iz;
      while 1
         iz(jz) = kz;
         rr1 = [1:jz-1, jz+1:m]; next = find( iz(rr1) == kz );
         if isempty(next), break; end
         jz = rr1(next(1));
         iz(jz) = 0;
         cc1 = [1:kz-1, kz+1:n]; next = find( ~D(jz,cc1) );
         if isempty(next), break; end
         kz = cc1(next(1));
      end
      if any( ismember( iz', izPrev,  'rows') )
         unfeas = 3; 
         iz = iz0;
         zz = [zz;kz0]; cc(cc==kz0) = []; 
      else
        izPrev = [ izPrev; iz' ];
      end
      break;
   end
else
   p = min( min( D(rr,cc) ) );
   if p >= Inf, unfeas = 4; break; end
   D(rr,:) = D(rr,:) - p;
   D(:,zz) = D(:,zz) + p;
end
end 
if unfeas==3,
   unfeas = 0; keepSets = 1;
else
   keepSets = 0;
end
if unfeas, break; end
ii = find( iz );
end 
function c =  robsolver(a)
iworth=find(a(:,3));
mfreight=sum(a(iworth,3))/length(iworth);
ia=find( a(:,4)>a(1,4)/14 | a(:,3)>=mfreight*.0713);
a=a(ia,:);
xy = (a(:,1)-a(1)) + 1i*(a(:,2)-a(1,2));
if ~any(abs(xy(2:end)) <= a(1,4))
  
  c=[1 1];
  return
end
n=size(a,1);
da=a(:,4)-a(:,3);
X = xy(:, ones(n,1));
dist = abs(X-X.');
dist0=dist;
I=1;
c=ones(1,n+1);
g=0;
for k=2:n
	g=g+a(I,4);
    dist(:,I)=inf;
    [Y,I]=min(dist(I,:));
	g=g-Y;
	if g<0
		k=k-1;
		break;
	end
    c(k) = I;
end
c=c(1:k+1);
g=ones(k,1);
g(1)=a(1,4)-dist0(1,c(2));
for ig=2:k	
	g(ig)=g(ig-1)+a(c(ig),4)-dist0(c(ig),c(ig+1));
end
if g(end)<0
	if k<3
		c=[1 1];
	elseif g(end)<dist0(c(k))
		
		g=g(1:k-2)+a(c(2:k-1),4)-dist0(c(2:k-1),1);
		i=find(g>=0);
		if isempty(i)
			c=[1 1];
		else
			i=i(end)+2;
			c=c(1:i);
			c(i)=1;
		end
	else
		c = c(1:k+1);
		c(end)=1;
	end
end
score1 = sum(da(c(2:end-1)));
if sum(a(c,3))/sum(a(:,3))<0.95
	if n<20
		c2=turbodiesel(a,dist0);
	else
		c2=diesel2(a,dist0);
	end
	score2 = sum(da(c2(2:end-1)));
	if(score2<score1)
		c=c2;
		score1=score2;
	end
	c3 = petrol3(a,dist0);
	score3 = sum(da(c3(2:end-1)));
	if(score3<score1)
		c=c3;
		score1=score3;
	end
end
da1=da;
da1(c)=0;
if any(da1<0)
	fail=0;
	while ~fail 
		[c,fail] = addGas(a,c,dist0);
		if fail 
			break 
		end;
	end
	fail=0;
	while ~fail 
		[c,fail] = addFreight(a,c,dist0);
		if fail 
			break 
		end;
	end
end
[c,score1]=optim(a,c,score1,xy,da);
if da(1)<score1
	
	c=[1 1];
end
c=ia(c);
function [c,fail] = addGas(a,c,dist)
fail=1;
gasv=a(:,4);
gasv(c)=0;
c=c(:);
Nc = length(c);
fIdx = find(gasv>0);
if isempty(fIdx), return; end 
N = length(fIdx);
d=ones(Nc-1,1);
for i=1:Nc-1	
	d(i)=dist(c(i),c(i+1));
end
eg = cumsum(a(c(1:Nc-1),4)-d);
dispgas=zeros(Nc-1,1);
dispgas(Nc-1)=eg(Nc-1);
for counter=Nc-2:-1:1
	dispgas(counter)=min(dispgas(counter+1),eg(counter));
end
ad = dist(c,fIdx);
h = dispgas(:,ones(1,N)) - (ad(1:Nc-1,:) + ad(2:Nc,:) - d(:,ones(1,N)));
proximity = max(h);
[ignore,idx] = sort(-a(fIdx,4).*proximity(:));
fIdx = fIdx(idx);
for n=1:N,
	
	i=fIdx(n);
	
	ad = dist(c,i);
	ed = ad(1:Nc-1) + ad(2:Nc) - d;  
if(any(dispgas>=ed))   
	h=dispgas-ed;
	[ignore,closest] = max(h);     
	
	if ed(closest)<a(fIdx(n),4)
		if h(closest)>=0
			fail=0;
			d1=dist(c(closest),i);
			d2=dist(i,c(closest+1));
			d0=dist(c(closest),c(closest+1));
			c  = [c(1:closest); i; c(closest+1:Nc)];
			d=[d(1:closest-1);d1;d2;d(closest+1:Nc-1)];
			eg = cumsum(a(c(1:Nc),4)-d);
			dispgas=zeros(Nc,1);
			dispgas(Nc)=eg(Nc);
			for counter=Nc-1:-1:1
				dispgas(counter)=min(dispgas(counter+1),eg(counter));
			end
			Nc=Nc+1;
		end
	end
end
end
function [c,fail] = addFreight(a,c,dist)
fail=1;
c=c(:);
Nc = length(c);
a(c,3)=0;
fIdx = find(a(:,3)>0);
if isempty(fIdx), return; end 
N = length(fIdx);
d=ones(Nc-1,1);
for i=1:Nc-1	
	d(i)=dist(c(i),c(i+1));
end
eg = cumsum(a(c(1:Nc-1),4)-d);
dispgas=zeros(Nc-1,1);
dispgas(Nc-1)=eg(Nc-1);
for counter=Nc-2:-1:1
    dispgas(counter)=min(dispgas(counter+1),eg(counter));
end
ad = dist(c,fIdx);
h = dispgas(:,ones(1,N)) - (ad(1:Nc-1,:) + ad(2:Nc,:) - d(:,ones(1,N)));
proximity = max(h);
[ignore,idx] = sort(-a(fIdx,3).*proximity(:));
fIdx = fIdx(idx);
for n=1:N
	
	i=fIdx(n);
	
	ad = dist(c,i);
	ed = ad(1:Nc-1) + ad(2:Nc) - d;     
if(any(dispgas>=ed))
	h=dispgas-ed;
	[ignore,closest] = max(h);     
	
	if h(closest)>=0
		fail=0;
		d1=dist(c(closest),i);