量子遗传算法求解背包问题程序

%qgan=100;%群体规模g=100;%进化代数m=50;w=[80,82,85,70,72,70,66,50,55,25,50,55,40,48,50,32,22,60,30,32,40,38,35,32,25,28,30,22,50,30, 45,30,60,50,20,65,20,25,30,10,20,25,15,10,10,10,4,4,2,1]; %物品重量p=[220,208,198,192,108,108,165,162,160,158,155,130,125,122,120,118,115,110,105,101,100,100 ,98,96,95,90,88,82,80,77,75,73,72,70,69,66,65,63,60,58,56,50,30,20,15,10,8,5,3,1]; %物品价值%clfclearglobal m n;%全局变量，m为染色体串长，即背包问题中的物品数量，n为群体规模m=input('please input chromsome length m=:');%输入串长w=[80,82,85,70,72,70,66,50,55,25,50,55,40,48,50,32,22,60,30,32,40,38,35,32,25,28,30,22,50,30, 45,30,60,50,20,65,20,25,30,10,20,25,15,10,10,10,4,4,2,1]p=[220,208,198,192,180,180,165,162,160,158,155,130,125,122,120,118,115,110,105,101,100,100 ,98,96,95,90,88,82,80,77,75,73,72,70,69,66,65,63,60,58,56,50,30,20,15,10,8,5,3,1]C=1000;save mwpC m w p C%knapsackclfclearglobal m n;%全局变量，m为染色体串长，即背包问题中的物品数量，n为群体规模m=input('please input chromsome length m=:');%输入串长for i=1:mw(i)=1+rand()*9;%物品重量p(i)=w(i)+5; %物品价值endC=sum(w)/2;%限制重量w,psave mwpc m w p C %保留数据，重复试验使用相同的数据%初始化群体，规模1，染色体位数10，%n=input('please input population size n=:');%群体规模%g=input('please input max-generation g=:');%进化代数for i=1:ma(i)=1/sqrt(2);%‘0’态系数b(i)=1/sqrt(2); %‘1’态系数end%MAX=zeros(number,g) %保持的最高适应度值%BEST=zeros(number,m)%保持的问题最优解q=zeros(2,m,n);%定义群体染色体for j=1:nfor i=1:mq(:,i,j)=[a(i),b(i)]';%单个染色体。

即q(1,i,j)为第j个染色体的第i位的‘0’态系数，%q(1,i,j)为第j个染色体的第i位的‘1’态系数endendfunction[q]=initialize(n,m)t=0;while t<g %进化循环t=t+1;observe;repair;evaluate;update;%量子门更新，产生下一代的量子态染色体store;endendfor t=1:gmaxqga10(t)=mean(MAX(:,t));aveqga10(t)=mean(ave(:,t));endsave data/data maxqga10 aveqga10for number=1:30plot(MAX(number,:));endplot(maxqga10,'r-');hold on;plot(aveqga10);%observe%x=zeros(n,m)function[x]=observe(n,m,q)n=100;m=50;q=zeros(2,m,n);%定义群体染色体for j=1:n %共有n个个体for i=1:m %每个个体为m位r=rand();if r>q(1,i,j)^2 %如果r>a(i,j)^2,则该位二进制串置为1，即取该物品x(j,i)=1;else x(j,i)=0;% 该为二进制串置为0，即不取该物品endendendplot(x);%repair 修改超重的问题解，即选择的物品重量不能超过限重Cfunction[overfiled]=repair(n,m,k,x,C)overfiled=0;%不超重n=100;m=50;x=zeros(n,m);w=[80,82,85,70,72,70,66,50,55,25,50,55,40,48,50,32,22,60,30,32,40,38,35,32,25,28,30,22,50,30, 45,30,60,50,20,65,20,25,30,10,20,25,15,10,10,10,4,4,2,1]; %物品重量c=1000;C=sum(w)/2;%限制重量for j=1:nif sum(x(j,:)*w')>C %超重overfiled=1;%超重符号endwhile overfiledk=fix(1+rand()*(m-1));%选择其中一个物品放弃,fix是求最接近0的整数x(j,k)=0;if sum(x(j,:)*w')<=C%不超重了overfiled=0;endendwhile ~overfiled %不超重k=fix(1+rand()*(m-1));%尽可能的再多选一个物品x(j,k)=1;if sum(x(j,:)*w')>C %超重了overfiled=1;endendx(j,k)=0;%将刚才选择后导致超重的那个物品丢弃x(j,:);endplot(k,fix(k));hold on;plot(k);plot(x);%evaluate 评估%fit=zeros(1,n);function[f,v]=observe(n,m,p)n=100;m=50;x=zeros(n,m);p=[220,208,198,192,108,108,165,162,160,158,155,130,125,122,120,118,115,110,105,101,100,100 ,98,96,95,90,88,82,80,77,75,73,72,70,69,66,65,63,60,58,56,50,30,20,15,10,8,5,3,1]; %物品价值for j=1:n %n个个体for i=1:m %m位。

即m个物品fit(j)=sum(x(j,:)*p');%问题解的适应度，即选择物品的总价值endend%for j=n:-1:2%if fit(j)>fit(j-1)% t=fit(j);%fit(j)=fit(j-1);% fit(j-1)=t;% end%end[f,v]=max(fit);%f位fit中的最大值，v为最大fit的位置，即本代最优解的对应序号g=100;t=1:g;if t==0ave0(n)=mean(fit);else ave(n,t)=mean(fit);endplot(f,t);hold on;plot(f);% update 量子门更新策略function[sign,BEST,angle,q]=update(n,m,t,z)n=100;m=50;g=100;t=1:g;for j=1:n %n个个体，依次更新if t==1if fit(j)>=MAX0(number);sign=1;else sign=0;endelseif fit(j)>=MAX(number,t-1)%本代的最优解比上代保持的最优解sign=1;else sign=0;endi=0;while i<m %n位i=i+1;if x(j,i)==0 %新解此位状态为‘0’，即第j个个体的第i位为0if BEST(number,i)==0 %上代保持的最优解此位为0angle=0;%不旋转elseif sign==0 %上代保持的最优解此位为‘1’，且新解不如上代保持的最优解适应度值高angle=0;%不旋转elseif q(1,i,j)*q(2,i,j)>0 %上代保持的最优解此位为‘1’，且新解比上代保持的最优解适应度%值高，且新解的‘0’态、‘1’态系数同号，即在一、三象限angle=-0.05*pi;%顺时针旋转，使下一代系数更靠近‘0’态elseif q(2,i,j)==0 %上代保持的最优解此位为‘1’，且新解比上代保持的最优解适应度值高，%且新解在实轴上，即b(i)即‘0’态，不旋转angle=0;elseangle=0.05*pi;%新解比上代保持的最优解适应度高，还包括两种情况，%一是新解在虚轴上，即a(i)处于‘1’态，此时无论顺时针还是逆时针旋转%都可以，使下个状态更靠近‘0’态，本程序使用顺指针逆时针旋转，%第二种情况是在二四象限，此时必须逆时针旋转endelseif BEST(n,i)==0 %新解状态为1，且保持的最优解状态为0if sign==0 % 新解适应度值小于保持的最优解适应度值if q(1,i,j)*q(2,i,j)>0%一三象限，所以顺时针旋转angle=-0.01*pi;elseif q(2,i,j)==0 %新解适应度值小于保持的最优值，且新解在实轴上，不旋转？？angle=0;else %两种情况。

一是新解在虚轴上，即在‘1’态，顺时针逆时针%旋转都可，本程序使用逆时针旋转，二是在三四象限，逆时针旋转angle=0.01*pi;endelseif q(1,i,j)*q(2,i,j)>0 %新解为1，保持的最优解为0，且新解适应度值比保持的最优解高，%且在一三象限，所以逆时针旋转angle=0.025*pi;elseif q(1,i,j)==0 %新解为1，保持的最优解为0，且新解适应度值比保持的最优解高，新解在虚轴上，不旋转angle=0;else %新解为1，保持的最优解为0，且新解适应度值比保持的最优解高，两种情况%一是新解在虚轴上，顺时针逆时针旋转都可，本程序使用顺时针旋转，二是在三四象限，顺时针旋转angle=-0.025*pi;endelseif sign==0%新解和保持的最优解都是1，且新解适应度值小于保持最优解if q(1,i,j)*q(2,i,j)>0 %一三象限，逆时针旋转？angle=0.005*pi;elseif q(1,i,j)==0 %新解在虚轴上angle=0;else %两种情况。