mxanthus_swarm_opt.m

一个用MATLAB编写的优化控制工具箱

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Bacterial swarm foraging optimization - M. xanthus model
%   
%   Kevin Passino
%   Version: 8/1/00
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear  all           % Initialize memory

flag=0; % If flag=0 indicates that will have nutrients, cell-cell attraction, and slime trails
        % If flag=1 indicates that will have no (zero) nutrients but have cell-cell attraction and slime trails
		% If flag=2 indicates that will have nutrients, no cell-cell attraction and but have slime trails

p=2;                 % Dimension of the search space (for bacteria p=2 or 3)

S=50;	% The number of bacteria in the population

Nc=100; % Number of chemotactic steps per bacteria lifetime, assumed
        % for convenience to be the same for every bacteria in the population
		% (keep even due to how plotting is done below, or divisible by 10 also)

% Define the domain of the optimization domain
thetamin=[0; 0];  % Set edges of region want to search in
thetamax=[30;30];

% Define a landscape (grid) for the bacteria to move on (assume that they stay only at grid points)
% This grid is defined for the domain of the function to be o

gridstep=0.1;

% Next find the grid point values and the number of grid points on each dimension

for i=1:p
	xt(:,i)=(thetamin(i,1):gridstep:thetamax(i,1))';
	Ngrid(i)=length(xt(:,i));
end

% Define a landscape on which the bacteria will move, and the nutrient/noxious substance map.  
% This is a discrete grid

nutrients=0*ones(Ngrid);

% 

temp=1;
if flag==1
	temp=0;
end

for ii=1:Ngrid(1)
	for jj=1:Ngrid(2)
		nutrients(jj,ii)=0.01*rand+temp*mxanthus_nutrients([xt(ii,1); xt(jj,2)]);  % Note odd index order on nutrients
											% - needed for the way Matlab plots in 3d
											% Add the rand since food is not uniform
	end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the nutrient map:

figure(1)
clf
surf(xt(:,1),xt(:,2),nutrients);
colormap(jet)
% Use next line for generating plots to put in black and white documents.
%colormap(white);
xlabel('x=\theta_1');
ylabel('y=\theta_2');
zlabel('z=J');
title('Nutrient concentration (valleys=food, peaks=noxious)');

%%%%%%%%%%%
% Define a grid on which trails are stored

slimedepth=0.1;  % Each cell secrets slime, this specifies the depth (amount)
slimeevaporationrate=0.98;  % Set rate at which slime will evaporate

slimetrails=0*ones(Ngrid); % This sets up a matrix that will be used to store the trails. 
slimetrailspart=slimetrails; % For use in plotting
% The values at each point on the grid will be either 0 (no trail) or -slimedepth which
% represents that a bacteria has been at this location before; however, will also evaporate


% Initial population

Pfront(:,:,:)=0*ones(p,S,Nc);  % First, allocate needed memory
Pback(:,:,:)=0*ones(p,S,Nc);  % First, allocate needed memory
Pindexfront(:,:,:)=0*ones(p,S,Nc);  % First, allocate needed memory
Pindexback(:,:,:)=0*ones(p,S,Nc);  % First, allocate needed memory


% Randomly place on domain for starting locations (setting cell orientations also): 

for m=1:S
	frontxindices(m)=round(Ngrid(1)*rand); % Generate random indices on each axis, for cell front
	frontxindices(m)=min(frontxindices(m),Ngrid(1));  % If put index beyond limit, then put it at the limit
	frontxindices(m)=max(frontxindices(m),1);  % If put index too low, put at lower limit
	frontyindices(m)=round(Ngrid(2)*rand);
	frontyindices(m)=min(frontyindices(m),Ngrid(1));  % If put index beyond limit, then put it at the limit
	frontyindices(m)=max(frontyindices(m),1);  % If put index too low, put at lower limit
end

Delta1=(2*round(rand(S,1))-1)'; % Generate random back locations (on diagonals), could be out of bounds but will be fixed below
backxindices=frontxindices+Delta1; 

for m=1:S
	backxindices(m)=min(backxindices(m),Ngrid(1));  % If put index beyond limit, then put it at the limit
	backxindices(m)=max(backxindices(m),1);  % If put index too low, put at lower limit
end

Delta2=(2*round(rand(S,1))-1)';
backyindices=frontyindices+Delta2;

for m=1:S
	backyindices(m)=min(backyindices(m),Ngrid(2));  % If put index beyond limit, then put it at the limit
	backyindices(m)=max(backyindices(m),1);  % If put index too low, put at lower limit
end

% Next, set bacteria locations and slime trails:

for m=1:S
	Pfront(:,m,1)=[xt(frontxindices(m),1); xt(frontyindices(m),2)]; % Place randomly on grid, fronts
	Pback(:,m,1)=[xt(backxindices(m),1); xt(backyindices(m),2)]; % Back locations
	Pindexfront(:,m,1)=[frontxindices(m);frontyindices(m)]; % Also save the index of the bacterium (front, where it senses)
	Pindexback(:,m,1)=[backxindices(m);backyindices(m)];  % Save indices of back of bacteria
	slimetrails(frontyindices(m),frontxindices(m))=-slimedepth; % Sets initial slime 
	slimetrails(backyindices(m),backxindices(m))=-slimedepth; 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot initial bacteria locations 

figure(2) 
clf
if flag~=1
	contour(xt(:,1),xt(:,2),nutrients,25)
	colormap(jet)
end
hold on
plot(squeeze(Pfront(1,:,1)),squeeze(Pfront(2,:,1)),'bx')
plot(squeeze(Pback(1,:,1)),squeeze(Pback(2,:,1)),'bo') 
xlabel('\theta_1');
ylabel('\theta_2');
title('Bacteria initial locations (x=front, o=back)')
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%---------------------------------
% Swarm optimization loop: 
%---------------------------------

for j=1:Nc
	
j  % Print swarm step to show progress

	for i=1:S  % For each bacterium
		
		% Compute the sensed values for the cell, for all its neighbors
		
		% First, set indices of neighboring locations to the head of the bacterium:
		
		n=[Pindexfront(1,i,j); Pindexfront(2,i,j)]+[0;1]; % north location to bacterium
		s=[Pindexfront(1,i,j); Pindexfront(2,i,j)]+[0;-1]; % south location to bacterium
		e=[Pindexfront(1,i,j); Pindexfront(2,i,j)]+[1;0]; % east location to bacterium
		w=[Pindexfront(1,i,j); Pindexfront(2,i,j)]+[-1;0]; % west location to bacterium
		
		ne=[Pindexfront(1,i,j); Pindexfront(2,i,j)]+[1;1]; % northeast location to bacterium
		nw=[Pindexfront(1,i,j); Pindexfront(2,i,j)]+[-1;1]; % northwest location to bacterium
		se=[Pindexfront(1,i,j); Pindexfront(2,i,j)]+[1;-1]; % southeast location to bacterium
		sw=[Pindexfront(1,i,j); Pindexfront(2,i,j)]+[-1;-1]; % southwest location to bacterium
		
		% Next, need to keep in domain by keeping the indices in range (keep fronts in range, and
		% fronts follow backs so they will stay in range):
		
		n=min(n,[Ngrid(1);Ngrid(2)]);  % If incremented index beyond limit, then put it at the limit
		n=max(n,[1;1]);  % If decremented index too low, put at lower limit
		s=min(s,[Ngrid(1);Ngrid(2)]);  % If incremented index beyond limit, then put it at the limit
		s=max(s,[1;1]);  % If decremented index too low, put at lower limit
		e=min(e,[Ngrid(1);Ngrid(2)]);  % If incremented index beyond limit, then put it at the limit
		e=max(e,[1;1]);  % If decremented index too low, put at lower limit
		w=min(w,[Ngrid(1);Ngrid(2)]);  % If incremented index beyond limit, then put it at the limit
		w=max(w,[1;1]);  % If decremented index too low, put at lower limit
		ne=min(ne,[Ngrid(1);Ngrid(2)]);  % If incremented index beyond limit, then put it at the limit
		ne=max(ne,[1;1]);  % If decremented index too low, put at lower limit
		nw=min(nw,[Ngrid(1);Ngrid(2)]);  % If incremented index beyond limit, then put it at the limit
		nw=max(nw,[1;1]);  % If decremented index too low, put at lower limit
		se=min(se,[Ngrid(1);Ngrid(2)]);  % If incremented index beyond limit, then put it at the limit
		se=max(se,[1;1]);  % If decremented index too low, put at lower limit
		sw=min(sw,[Ngrid(1);Ngrid(2)]);  % If incremented index beyond limit, then put it at the limit
		sw=max(sw,[1;1]);  % If decremented index too low, put at lower limit

		% Compute the nutrient concentration:
		
		Jn=nutrients(n(2,1),n(1,1)); 
		Js=nutrients(s(2,1),s(1,1)); 
		Je=nutrients(e(2,1),e(1,1)); 
		Jw=nutrients(w(2,1),w(1,1)); 
		Jne=nutrients(ne(2,1),ne(1,1)); 
		Jnw=nutrients(nw(2,1),nw(1,1)); 
		Jse=nutrients(se(2,1),se(1,1)); 
		Jsw=nutrients(sw(2,1),sw(1,1)); 

		% Add the cell-to-cell attractant:
		
		Jn=Jn+mxanthus_attract_func([xt(n(1,1),1);xt(n(2,1),2)],Pfront(:,:,j),S,flag);
		Js=Js+mxanthus_attract_func([xt(s(1,1),1);xt(s(2,1),2)],Pfront(:,:,j),S,flag);
		Je=Je+mxanthus_attract_func([xt(e(1,1),1);xt(e(2,1),2)],Pfront(:,:,j),S,flag);
		Jw=Jw+mxanthus_attract_func([xt(w(1,1),1);xt(w(2,1),2)],Pfront(:,:,j),S,flag);
		Jne=Jne+mxanthus_attract_func([xt(ne(1,1),1);xt(ne(2,1),2)],Pfront(:,:,j),S,flag);
		Jnw=Jnw+mxanthus_attract_func([xt(nw(1,1),1);xt(nw(2,1),2)],Pfront(:,:,j),S,flag);
		Jse=Jse+mxanthus_attract_func([xt(se(1,1),1);xt(se(2,1),2)],Pfront(:,:,j),S,flag);
		Jsw=Jsw+mxanthus_attract_func([xt(sw(1,1),1);xt(sw(2,1),2)],Pfront(:,:,j),S,flag);
		
		
		% Add the slime-trail effects:
		
		Jn=Jn+slimetrails(n(2,1),n(1,1)); % n=north
		Js=Js+slimetrails(s(2,1),s(1,1)); % s=south
		Je=Je+slimetrails(e(2,1),e(1,1)); % e=east
		Jw=Jw+slimetrails(w(2,1),w(1,1)); % w=west
		Jne=Jne+slimetrails(ne(2,1),ne(1,1)); % ne=northeast, etc.
		Jnw=Jnw+slimetrails(nw(2,1),nw(1,1)); 
		Jse=Jse+slimetrails(se(2,1),se(1,1)); 
		Jsw=Jsw+slimetrails(sw(2,1),sw(1,1)); 
		

		% There are 8 possible locations of the back of the cell, if we consider the front
		% to be at the center of a grid.  Consider each in sequence, and make moves (note
		% that we cannot move out of bounds):
		
		% Back of cell at N location:
		
		if Pindexback(:,i,j)==n
			Jmove=min([Jw,Jsw,Js,Jse,Je]); % Find the minimum cost of the three possible forward locations+left/right
			if Jmove==Jw % best to move in W direction
				Pindexfront(:,i,j+1)=w; 
			elseif Jmove==Jsw % best to move in SW direction
				Pindexfront(:,i,j+1)=sw;