genetic_adaptive_ind.m

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

      
% Next, make the parent chromosome

    parent_chrom(:,pop_member)=pop(:,member_count);
	end
   end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reproduce section (i.e., make off-spring - "children")
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% In the approach below each individual gets to mate and
% they randomly pick someone else (not themselves) in the 
% mating pool to mate with.  Resulting children are 
% composed of a combination of genetic material of their
% parents.  If elitism is on, when the elite member gets
% a chance to mate they do not; they are simply copied
% over to the next generation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for parent_number1 = 1:POP_SIZE,    % Crossover (parent_number1 is the 
									% individual who gets to mate
	if ELITISM ==1 & parent_number1==bestmember(k) % If elitism on, and 
											 % have the elite member
		child(:,parent_number1)=parent_chrom(:,parent_number1);
	else
	   parent_number2=parent_number1;  % Initialize who the mate is
	   while parent_number2 == parent_number1   % Iterate until find 
		                                        % a mate other than
												% yourself
         parent_number2 = rand*POP_SIZE; % Choose parent number 2
		 								 % randomly (a random mate)
         parent_number2 = parent_number2-rem(parent_number2,1)+1;     
	   end
    if CROSS_PROB > rand           % If true then crossover occurs

         site = rand*CHROM_LENGTH;   % Choose site for crossover
         site = site-rem(site,1)+1;  % and make it a valid integer 
         							 % number for a site

% The next two lines form the child by the swapping of genetic
% material between the parents

child(1:site,parent_number1)=parent_chrom(1:site,parent_number1);

child(site+1:CHROM_LENGTH,parent_number1)=...
     parent_chrom(site+1:CHROM_LENGTH,parent_number2);

      else                           % No crossover occurs

% Copy non-crossovered chromosomes into next generation
% In this case we simply take one parent and make them
% the child.

child(:,parent_number1)=parent_chrom(:,parent_number1);

	 end
   end  % End the "if ELITISM..." statement
end  % End "for parent_number1=..." loop

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mutate children.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Here, we mutate to a different allele
% with a probability MUTAT_PROB
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   for pop_member= 1:POP_SIZE,

	if ELITISM ==1 & pop_member==bestmember(k)  % If elitism on, and 
											 % have the elite member
		child(:,pop_member)=child(:,pop_member); % Do not mutate
											% the elite member
	else
	  for site = 1:CHROM_LENGTH,

         if MUTAT_PROB > rand        % If true then mutate  
            rand_gene=rand*10; 		 % Creat a random gene
            
            % If it is the same as the one already there then 
            % generate another random allele in the alphabet
            
            while child(site,pop_member) == rand_gene-rem(rand_gene,1),
               rand_gene=rand*10;
            end;
            
            % If it is not the same one, then mutate
            
            child(site,pop_member)=rand_gene-rem(rand_gene,1);
            
            % If takes a value of 10 (which it cannot
            % mutate to) then try again (this is a very low probability
			% event (most random number generators generate numbers
			% on the *closed* interval [0,1] and this is why this line
			% is included).
            
            if rand_gene == 10
                 site+site-1;
            end
         end  % End "if MUTAT_PROB > rand ... 
	   end  % End for site... loop
	 end  % End "if ELITISM..."
   end  % End for pop_member loop

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create the next generation 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   pop=child;               % Create next generation (children 
   							% become parents)
							
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Next, we have to convert the population (pop) to the base-10 
% representation (called trait) so that we can check if the traits
% all still lie in the proper ranges specified by HIGHTRAIT and LOWTRAIT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for pop_member = 1:POP_SIZE

	for current_trait = 1:NUM_TRAITS,
			 
trait(current_trait,pop_member,k)=0; % Initialize variables
place_pointer=1;

% Change each of the coded traits on the chromosomes into base-10 traits: 
% For each gene on the current_trait past the sign digit but before the
% next trait find its real number amount and hence after finishing
% the next loop trait(current_trait,pop_member,k) will be the base-10
% number representing the trait

for gene=TRAIT_START(current_trait)+1:TRAIT_START(current_trait+1)-1,
 place=DECIMAL(current_trait)-place_pointer;
 trait(current_trait,pop_member,k)=...
 trait(current_trait,pop_member,k)+...
    (pop(gene,pop_member))*10^place;
 place_pointer=place_pointer+1;
end

% Determine sign of the traits and fix 
% trait(current_trait,pop_member,k) so that it has the right sign:

if pop(TRAIT_START(current_trait),pop_member) < 5
   trait(current_trait,pop_member,k)=...
      -trait(current_trait,pop_member,k);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fix bad traits (i.e., ones that are out of the range 
% specified by HIGHTRAIT and LOWTRAIT) by saturation at the extremes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if trait(current_trait,pop_member,k)>HIGHTRAIT(current_trait)

                     % The trait has went higher than the upper
                     % bound so let the trait equal to the
                     % HIGHTRAIT bound.

trait(current_trait,pop_member,k)=HIGHTRAIT(current_trait);

% Now consider the other case:

elseif trait(current_trait,pop_member,k)<LOWTRAIT(current_trait)

                     % The trait has went lower than the lower 
                     % bound so let the trait equal to the
                     % LOWTRAIT bound

trait(current_trait,pop_member,k)=LOWTRAIT(current_trait);
                  
		end


	end % Ends "for current_trait=..." loop
         
end    % Ends "for pop_member=..." loop
 
	
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% End main GA operations
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
else  % Hold estimates constant if have not gotten enough data
	  % to be able to compute the fitness function (hence, adaptation
	  % does not start until after N steps)

	trait(:,:,k)=trait(:,:,k-1);
	
	% Pick the parameters to be the same (here, just pick it to be the first
	% population member - since they are all the same initially it does not
	% matter which one you pick)
	
	thetaalpha(k)=trait(1,1,k-1);  
	thetabeta(k)=trait(2,1,k-1);  
	
end

	% Next, find the estimates of the plant dynamics (elite member)
	
	betahat(k)=thetabeta(k);
	alphahat(k)=thetaalpha(k)*h(k);

	% Store the estimate of the plant dynamics
	
	hhat(k)=alphahat(k-1)+betahat(k-1)*u(k-1);
	
	% Next, use the certainty equivalence controller
	
	u(k)=(1/(betahat(k)))*(-alphahat(k)+r(k+1));
		
	% Define some parameters to be used in the plant
	
	A(k)=abs(abar*h(k)+bbar);

	alpha(k)=h(k)-T*dbar*sqrt(2*g*h(k))/A(k);
	beta(k)=T*cbar/A(k);
	
end


%%%%%%%%
% Plot the results

figure(1)
clf
subplot(211)
plot(time,h,'k-',time,ref,'k--')
grid
ylabel('Liquid height, h')
title('Liquid level h and reference input r')

subplot(212)
plot(time,u,'k-')
grid
title('Tank input, u')
xlabel('Time, k')
axis([min(time) max(time) -50 50])


%%%%%%%%
figure(2)
clf
subplot(311)
plot(time,h,'k-',time,hhat,'k--')
grid
title('Liquid level h and estimate of h')

subplot(312)
plot(time,alpha,'k-',time,alphahat,'k--')
grid
title('Plant nonlinearity \alpha and its estimate')

subplot(313)
plot(time,beta,'k-',time,betahat,'k--')
grid
xlabel('Time, k')
title('Plant nonlinearity \beta and its estimate')


%%%%%%%%
figure(3)
clf
subplot(211)
plot(time,Jbaravg,'k-')
grid
title('Average fitness')

subplot(212)
plot(time,bestmember,'k-')
grid
xlabel('Time, k')
T=num2str(POP_SIZE);
T=strcat('Index of best member in population of size=',T);
title(T)



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