beeforage.m

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Simulation of social foraging of bees, an optimization method for
% finding a maximum of a function, and for resource (bee) allocation.
%
% By: K. Passino
% Version: 11/12/01
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all

% Initialize

p=2; % Set the number of dimensions of the search domain (here =2 for
     % easy visualization)
	 
S=100; % The number of potential foragers/scouts (unemployed, employed, scouts)
Ss=5; % The number of scouts (let this be a constant)

Ns=300; % Total number of foraging trips to take (sets the length
		 % of the simulation) (for convenience choose this to be divisible by 4)
		 
% Set parameters that model the collection of nectar (for one option)
%alpha=1;
%beta=10;

% Set parameter that sets the threshold for dancing for recruitment 
% to a site

delta=0.5;

% Set parameters that control the percentage of unemployed foragers that
% can be recruited in one step

gamma1=0.1*(S-Ss);
gamma2=1;

% Set the noise level that represents inaccuracies in waggle dance signaling
% and problems with a recruited forager finding the site she was recruited to.

sigma_r=2; % Set variance on noise that is added to forage site location in each 
           % dimension (another way would be to use polar coordinates and induce
		   % noise on the range and angle signaled by a bee) for a newly recruited
		   % forager

% If uncomment-out some lines below can add a feature where employed foragers will not
% necessarily perfectly make it back to the last place they were foraging.
%sigma_e=0.5; % Set variance that represents that a bee will not find exactly where they 
            % were foraging the last time that they went foraging

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Will also determine which forage sites each 
% forager is at; here, intialize
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Forage sites are (see beeforagefunc.m):
%
% Number   Site
% 1			(15,20)
% 2 		(-20,15)
% 3			(20,-10)
% 4			(-15,-15)

% Set threshold for being at a site
epsilon=10; % Note that this is small enough so that there is no overlap of regions so there
             % is no double counting of foragers at sites

foragesite=0*ones(4,Ns); % Initialize
epsilonq=0.001; % The degration factor for quality if a bee visits a forage site


%-------------------------------------------------------------------
% Nectar profile parameters option: Two equal sites, then one increases 
% (to test cross-inhibition between two sites)
Nc=0*ones(4,Ns+1); 
Nc(1,1:Ns+1)=3*ones(1,Ns+1); % One option, constant heights
Nc(2,1:Ns+1)=3*ones(1,Ns+1);
Nc(3,1:Ns+1)=0*ones(1,Ns+1);
Nc(4,1:Ns+1)=0*ones(1,Ns+1);

Nc(1,((Ns/2)+1):Ns+1)=8*ones(1,(Ns/2)+1);  % Makes site 1 much more profitable half
                                       % way through

%-------------------------------------------------------------------
% Nectar profile parameters option: To test case where forage site 1 
% becomes a poorer site, while site 2
% becomes a better site (to show reallocation)

Nc(1,1:Ns+1)=2*ones(1,Ns+1); % One option, constant heights
Nc(2,1:Ns+1)=5*ones(1,Ns+1);
Nc(1,1:(Ns/2))=5*ones(1,(Ns/2)); % Keep at 5 for half the time
Nc(2,1:(Ns/2))=2*ones(1,(Ns/2)); % Keep at 3 for half the time
Nc(3,1:Ns+1)=0*ones(1,Ns+1); % Makes site 3 nonexistent for all time
Nc(4,1:Ns+1)=0*ones(1,Ns+1); % Makes site 4 nonexistent for all time

%-------------------------------------------------------------------
% Nectar profile parameters option: Test when one site is not there the
% whole time then later "pops up"
Nc(1,1:Ns+1)=0.25*ones(1,Ns+1); % One option, constant heights
Nc(2,1:Ns+1)=0.25*ones(1,Ns+1);
Nc(3,1:Ns+1)=0.25*ones(1,Ns+1);
Nc(4,1:Ns+1)=0*ones(1,Ns+1); % Makes site 4 nonexistent for all time
Nc(4,((Ns/2)+1):Ns+1)=8*ones(1,(Ns/2)+1); % Makes site 4 the best site of all at half way through

%-------------------------------------------------------------------
% Nectar profile parameters option: Constant parameters
Nc=0*ones(4,Ns+1); 
Nc(1,1:Ns+1)=5*ones(1,Ns+1); % One option, constant heights
Nc(2,1:Ns+1)=3*ones(1,Ns+1);
Nc(3,1:Ns+1)=2*ones(1,Ns+1);
Nc(4,1:Ns+1)=1*ones(1,Ns+1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Next, plot the nectar concentration function that the bees search over.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% First, set the range of values (domain size) around the hive which is 
% situated at (0,0).

domainsize=30; % Make it a domain of size 2*domainsize in each dimension 
				% (hence, it is a square area)

x1=-domainsize:1:domainsize;   % For our function the range of values we are considering
x2=x1;

% Compute the function that we are trying to find the minimum of.

for jj=1:length(x1)
	for ii=1:length(x2)
		z(ii,jj)=beeforagefunc([x1(jj);x2(ii)],Nc(:,1));
	end
end

% First, show the actual function to be maximized

figure(1)
clf
surf(x1,x2,z);
% Use next line for generating plots to put in black and white documents.
%colormap(white);
xlabel('\theta_1');
ylabel('\theta_2');
zlabel('Jbar');
title('Nectar concentration/quality function.');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set the intial bee locations


%-------------------------------------------------------------------
% Initialization option: Initialize all foragers at hive site (0,0)
theta=0*ones(p,S,Ns); 

% Indicate number employed and unemployed

Se=0*ones(Ns,1); 
Se(1,1)=0; % Sets the number of employed foragers at the first iteration 
	       % (now set to zero, so that only scouts are sent out)

Su=0*ones(Ns,1); 
Su(1,1)=S-Ss-Se(1,1); % Sets the number of unemployed foragers 

%-------------------------------------------------------------------
% Intialization option: Set all bees foraging at either site 1 or 2
% Forage sites are (see beeforagefunc.m):
%
% Number   Site
% 1			(15,20)
% 2 		(-20,15)
% 3			(20,-10)
% 4			(-15,-15)

theta=0*ones(p,S,Ns); 

for i=1:(S-Ss-S/2)
theta(:,i,1)=[15; 20]+[10*rand-5; 10*rand-5]; % Near site 1
end

for i=(S-Ss+1-S/2):(S-Ss)
theta(:,i,1)=[-20; 15]+[10*rand-5; 10*rand-5];  % Near site 2
end

% Indicate number employed and unemployed
Su(1,1)=0;
Se(1,1)=S-Ss;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize other variables:

Jbar=0*ones(S,Ns); % Nectar concentration profile samples for each bee at each step
N=0*ones(S,Ns); % Amount of nectar retrieved by ith be on kth run
Nt=0*ones(Ns,1); % Total amount of nectar per forage trip by group
L=0*ones(S,Ns); % Dance length

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set initial scout positions
for i=(S-Ss+1):S % Set scouts at random locations on domain (no preference to location to scout)
	theta(:,i,1)=[2*domainsize*rand-domainsize; 2*domainsize*rand-domainsize];
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Next, start the main foraging loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for k=1:Ns

% For the sake of plotting, first compute the number of bees at each forage site:
% (however, there is also an option below for making the nectar concentration/quality 
% landscape dependent on the number of bees that have foraged there (so each time we 
% simulate as if some of the nectar were eaten so that the quality of the overall 
% foraging site goes down)

	for i=1:S
			
		if sqrt((theta(:,i,k)-[15; 20])'*(theta(:,i,k)-[15; 20]))<epsilon % Check forage site 1
			foragesite(1,k)=foragesite(1,k)+1; % Add one more member at site 1
		end
		if sqrt((theta(:,i,k)-[-20; 15])'*(theta(:,i,k)-[-20; 15]))<epsilon % Check forage site 2
			foragesite(2,k)=foragesite(2,k)+1; % Add one more member at site 2
		end
		if sqrt((theta(:,i,k)-[20; -10])'*(theta(:,i,k)-[20; -10]))<epsilon % Check forage site 3
			foragesite(3,k)=foragesite(3,k)+1; % Add one more member at site 3
		end
		if sqrt((theta(:,i,k)-[-15; -15])'*(theta(:,i,k)-[-15; -15]))<epsilon % Check forage site 4
			foragesite(4,k)=foragesite(4,k)+1; % Add one more member at site 4
		end
		
	end

	
%-----------------------------------------------------
% Forage (collect nectar), return, and unload