beeforage.m

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

%-----------------------------------------------------

Nt(k,1)=0; % Intialize total amount of nectar accumulator variable

for i=1:S
	if theta(:,i,k)==[0; 0]
		Jbar(i,k)=0; % An unemployed forager does not sample the nectar profile
		N(i,k)=0; 
	else
	Jbar(i,k)=beeforagefunc(theta(:,i,k),Nc(:,k)); % Sample nectar nutrient/quality landscape
	N(i,k)=Jbar(i,k); % Nectar brought back by bee i
%	N(i,k)=min([alpha*Jbar(i,k), beta]); % Nectar brought back by bee i (where can saturate assessment capability)
	Nt(k,1)=Nt(k,1)+N(i,k); % Accumulate nectar
    end
end

%-------------------------------------------------------
% An optional feature: change foraging landscape based on 
% the number of foragers that visit it
%-------------------------------------------------------

for i=1:4
	Nc(i,k+1)=max([Nc(i,k)-epsilonq*foragesite(i,k), 0]); % Remove some epsilon amount of quality/concentration
	                                                    % for each be at foragesite i at time k
end

%-----------------------------------------------------
% Determine dance length, and create unemployed foragers
%-----------------------------------------------------

for i=1:S
	if theta(:,i,k)==[0; 0]
		L(i,k)=0; % An unemployed forager does not dance
	else
	L(i,k)=max([N(i,k)-(delta/(Se(k,1)+Ss))*Nt(k,1), 0]); % Bee i determines dance length
	end

end

%-----------------------------------------------------
% Dance, recruit, and set positions of foragers/scouts
%-----------------------------------------------------

theta(:,:,k+1)=theta(:,:,k); % Set new forage sites to be the same as the previous ones (implies that 
                             % employed foragers return precisely to same site as last run); however,
							 % will modify this below also...

Su(k+1,1)=0;

% Create unemployed foragers, based on how good they foraged at the last step (this sets the number 
% of unemployed foragers before  they watch dances)
for i=1:(S-Ss)
	if L(i,k)==0
		theta(:,i,k+1)=[0;0]; % Here, make it unemployed for the next iteration
		Su(k+1,1)=Su(k+1,1)+1; % Change counts of numbers of unemployed and employed foragers
% Next, shows an option where you could add another feature to the program:
%	else % It is an employed forager and so set where it will forage at in the next step 
		 % (note that since sigma_e is chosen to be small the bee will find its way back close
		 % to the previous profitable site that it had found)
%		theta(:,i,k+1)=theta(:,i,k+1)++[sqrt(sigma_e)*randn; sqrt(sigma_e)*randn];
	end
end

Se(k+1,1)=S-Ss-Su(k+1,1);

% Place the scouts (scouts will be held in last Ss positions in theta):

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

Lscaled(k)=gamma1*sum(L(:,k))/(1+gamma2*sum(L(:,k))); % Used below

if Su(k+1,1)>0  % If there are some unemployed foragers
	
Nr=min([floor(Lscaled(k)), Su(k+1,1)]);

for j=1:Nr % Recruit some percentage of unemployed foragers
	
	% For some unemployed forager
	for i=1:(S-Ss)
		if theta(:,i,k+1)==[0; 0]
			jstar=i; % Finds an unemployed forager 
			break % Found an unemployed forager so break loop
		end
	end

	% Recruit the unemployed forager by some employed one, in a proportional manner to how
	% long the dancing of that employed forager is (like fitness-proportionate selection for
	% genetic algorithms, here call it dance-length proportionate recruitment)
    pointer=rand*sum(L(:,k)); % This makes the pointer for the roulette wheel.  
    bee_index=1;              % Initialization
    total=L(1,k);

    while total < pointer,                % This spins the wheel to the 
       									  % pointer and finds the 
       									  % bee there - which is
       									  % identified by bee_count - and 
										  % this is the recruiter bee
    bee_index=bee_index+1;
    total=total+L(bee_index,k);  % Notice that if L(bee_index,k)=0 (unemployed forager) then it will
								 % add nothing to total so it will skip that index (so this makes
								 % it so that unemployed foragers will only follow employed foragers
								 % or scouts.
    end  % At the end bee_index is the index of the recruiter bee for the unemployed forager jstar

	% Recruit to forage site
	theta(:,jstar,k+1)=theta(:,bee_index,k+1)+[sqrt(sigma_r)*randn; sqrt(sigma_r)*randn]; 
	    % Unemployed forager follows dance of employed forager, but not perfectly
	
	Su(k+1,1)=Su(k+1,1)-1; % Change counts of numbers of unemployed and employed foragers
	Se(k+1,1)=Se(k+1,1)+1;
end

end % End if there are some unemployed foragers loop

end % End of main loop



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


figure(2) 
clf
contour(x1,x2,z,25);
grid
% Use next line for generating plots to put in black and white documents.
%colormap(white);
hold on
for i=1:(S-Ss) % Plot the forager positions
	xx=squeeze(theta(1,i,:));
	yy=squeeze(theta(2,i,:));
	siteplot=plot(xx,yy,'k.');
	set(siteplot,'MarkerSize',8);
end
for i=(S-Ss+1):S % Plot the scout positions
	xx=squeeze(theta(1,i,:));
	yy=squeeze(theta(2,i,:));
	siteplot=plot(xx,yy,'r.');
	set(siteplot,'MarkerSize',4);
end
xlabel('x_1=\theta_1');
ylabel('x_2=\theta_2');
title('Nectar concentration (contour plot) and forage sites');
hold off

time=0:Ns;

figure(3)
clf
plot(time,Se,'k-',time,Su,'b--')
xlabel('Iteration, k')
ylabel('S_e, S_u')
title('Number of employed (-) and unemployed (--) foragers')

time2=0:Ns-1;

figure(4)
clf
plot(time2,Nt,'b-')
xlabel('Iteration, k')
ylabel('N_t')
title('Total amount of nectar collected')


figure(5)
clf
plot(time2,Lscaled,'b-')
xlabel('Iteration, k')
ylabel('N_t')
title('Scaled total dance length')


figure(6)
clf
subplot(221)
plot(time2,foragesite(1,:))
axis([min(time2) max(time2) 0 S])
xlabel('Iternation, k')
ylabel('F_1')
title('Number of bees at forage site 1')
subplot(222)
plot(time2,foragesite(2,:))
axis([min(time2) max(time2) 0 S])
xlabel('Iternation, k')
ylabel('F_2')
title('Number of bees at forage site 2')
subplot(223)
plot(time2,foragesite(3,:))
axis([min(time2) max(time2) 0 S])
xlabel('Iternation, k')
ylabel('F_3')
title('Number of bees at forage site 3')
subplot(224)
plot(time2,foragesite(4,:))
axis([min(time2) max(time2) 0 S])
xlabel('Iternation, k')
ylabel('F_4')
title('Number of bees at forage site 4')

figure(7)
clf
subplot(221)
plot(time,Nc(1,:))
axis([min(time) max(time) 0 10])
xlabel('Iternation, k')
ylabel('N_c_1')
title('Nectar concentration/quality at forage site 1')
subplot(222)
plot(time,Nc(2,:))
axis([min(time) max(time) 0 10])
xlabel('Iternation, k')
ylabel('N_c_2')
title('Nectar concentration/quality at forage site 2')
subplot(223)
plot(time,Nc(3,:))
axis([min(time) max(time) 0 10])
xlabel('Iternation, k')
ylabel('N_c_3')
title('Nectar concentration/quality at forage site 3')
subplot(224)
plot(time,Nc(4,:))
axis([min(time) max(time) 0 10])
xlabel('Iternation, k')
ylabel('N_c_4')
title('Nectar concentration/quality at forage site 4')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End of program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%