小编刚开始看书中的概念感觉讲得云里雾里,于是小编直接看书中的代码,边看代码边理解书中所讲的内容。
量子遗传算法与标准的遗传算法的区别有两点:
编码方式不同,把遗传算法的染色体编码操作符换为量子比特概率幅表示
种群进化方式不同,标准的量子遗传算法采用量子旋转门更新种群,而不是采用标准遗传算法的交叉和变异操作。
代码求解的问题是一个复杂二元函数求最值问题
使用量子遗传算法求解这样一个二元函数问题的第一步是如何进行编码?
标准遗传算法求解这个问题的编码一般会采用二进制编码,而量子遗传算法会采用量子比特编码的方式进行编码,不过依然会将量子比特编码转换成二进制编码,最终转换成实数(十进制)。
01 | 染色体编码
小编在这里简单介绍一下二进制编码,以代码中20位二进制编码为例,比如说我们想把[1,3]这个区间任意一个实数用20位二进制编码的形式表现出来,怎么办呢?
以一个5位二进制串为例01101,这个二进制串对应的实数为
如果二进制串为11111,则对应的实数为
所以用20位二进制串编码,能得到最大的实数为。
这时候我们对二进制串进行等比例缩放,缩放到[1,3]这个区间内,假设这20位二进制串每一位上的数字用b来表达(b要么等于1,要么等于0),则这20位二进制串表达为,则所得到的实数为
所以任意一个二进制串映射在[1,3]这个区间内所能得到的实数为
二进制转换成十进制的代码如下
function X=bin2decFun(x,lenchrom,bound)
%% 二进制转化成十进制
% 输入 x:二进制编码
% lenchrom:各变量的二进制位数
% bound:各变量的范围
% 输出 X:十进制数
M=length(lenchrom);
n=1;
X=zeros(1,M);
for i=1:M
for j=lenchrom(i)-1:-1:0
X(i)=X(i)+x(n).*2.^j;
n=n+1;
end
end
X=bound(:,1)'+X./(2.^lenchrom-1).*(bound(:,2)-bound(:,1))';
我们现在只是知道如何将二进制转换成十进制,但是不知道很关键的一步,量子遗传算法如何用量子比特编码的方式进行编码,有了这一步才会将量子比特编码转换成二进制编码,然后才有上述将二进制编码转换成实数(十进制)的步骤。
常规遗传算法求解上述二元函数的编码形式为,因为有两个变量,所以需要两端长度为20的二进制串。
量子比特编码采用一种叠加的形式进行编码,就书中代码而言,一个染色体的量子比特编码形式为
其中表示第i个量子比特的概率幅(这是术语,可忽略不看),并且满足
接下来详细介绍一下量子比特编码:
量子比特(qutbit)是一个充当信息存储单元的物理介质的双态量子系统,是定义在二维复向量空间中的一个单位向量,该空间由一对特定的标准正交基{| 0>,| 1>}组成。因此,它可以同时处于两个量子态的叠加态中。如,其中0和1分别表示自旋向下态和自旋向上态,“|〉”为一种量子态的表示,α,β是两个复数称几率幅对。可看成量子处于自旋向下态的概率,而可看成量子处于自旋向上态的概率。
在量子遗传算法中,每个自变量的分量可用一个量子基因来表示,如可用
表示第i个量子基因,而每个量子基因又是由若干位(k位)量子比特组成,每个量子比特具有两个状态,即分别以一定的概率和取到0态和1态.因此,每个多维自变量(n维)就可以用n个量子基因组成.显然,若问题维数为n,则每个量子个体由n个量子基因组成,每个量子个体可用如下的量子编码形式编码
其中为第t代第j个个体,k为每个自变量的分量所用的量子比特数.
量子比特编码就是采用上述编码形式,一般来说将初始设置为
下述代码是种群初始化的代码
function chrom=InitPop(M,N)
%% 初始化种群-量子比特编码
% M:为种群大小×2,(α和β)
% N:为量子比特编码长度
for i=1:M
for j=1:N
chrom(i,j)=1/sqrt(2);
end
end
虽然我们知道量子遗传算法如何用量子比特编码的方式进行编码,但是我们不知道如何将量子比特编码转换成二进制编码。
转换过程中,随机产生一个[0,1]区间的数,若它大于概率幅(或)的平方,则将对应的二进制编码赋值为1,否则为0
子比特编码转换成二进制编码的代码如下
function binary=collapse(chrom)
%% 对种群实施一次测量 得到二进制编码
% 输入chrom :为量子比特编码
% 输出binary:二进制编码
[M,N]=size(chrom); %得到种群大小 和编码长度
M=M/2; % 种群大小
binary=zeros(M,N); %二进制编码大小初始化
for i=1:M
for j=1:N
pick=rand; %产生【0,1】随机数
if pick>(chrom(2.*i-1,j)^2) % 随机数大于α的平方
binary(i,j)=1;
else
binary(i,j)=0;
end
end
end
经过上述一系列操作,可以将量子比特编码转换成实数,所以求单个染色体目标函数值的代码如下:
function [Y,X]=Objfunction(x,lenchrom)
%% 目标函数
% 输入 x:二进制编码
% lenchrom:各变量的二进制位数
% 输出 Y:目标值
% X:十进制数
bound=[-3.0 12.1;4.1 5.8]; % 函数自变量的范围
%% 将binary数组转化成十进制数组
X=bin2decFun(x,lenchrom,bound);
%% 计算适应度-函数值
Y=sin(4*pi*X(1))*X(1)+sin(20*pi*X(2))*X(2);
求种群适应度值的代码如下
function [fitness,X]=FitnessFunction(binary,lenchrom)
%% 适应度函数
% 输入 binary:二进制编码
% lenchrom:各变量的二进制位数
% 输出 fitness:适应度
% X:十进制数(待优化参数)
sizepop=size(binary,1);
fitness=zeros(1,sizepop);
num=size(lenchrom,2);
X=zeros(sizepop,num);
for i=1:sizepop
[fitness(i),X(i,:)]=Objfunction(binary(i,:),lenchrom); % 使用目标函数计算适应度
end
02 | 种群进化
在了解量子遗传算法如何用量子比特编码的方式进行编码以后,还需知道如何将种群进化,量子遗传算法采用量子门更新的方法实现种群的进化。
我们通过上述的学习知道,量子遗传算法中一个染色体的编码形式为
所以在进化的过程中,上述需要按照某种方式变化,进而朝着理想的方向进化。
书中提到使用量子旋转门(术语,可忽略)
实现的更新,更新方式如下:
但是
可以看出变换之后的值仍为1。
但是公式中的是怎么确定的?具体做法如下:
为当前染色体的第i位(这里的染色体是指转换为二进制编码后的染色体);为当前的最优染色体的第i位;f(x)为适应度函数;为旋转方向,为旋转角大小,和的值根据下表确定,则中的
下面解释一下上表怎么使用?
将当前个体的适应度值与该种群中最优个体适应度值进行比较,如果,则使向着有利于出现的方向演化;反之如果,则使向着有利于出现的方向演化。
种群进化的代码如下
function chrom=Qgate(chrom,fitness,best,binary)
%% 量子旋转门调整策略
% 输入 chrom:更新前的量子比特编码
% fitness:适应度值
% best:当前种群中最优个体
% binary:二进制编码
% 输出 chrom:更新后的量子比特编码
sizepop=size(chrom,1)/2;
lenchrom=size(binary,2);
for i=1:sizepop
for j=1:lenchrom
A=chrom(2*i-1,j); % α
B=chrom(2*i,j); % β
x=binary(i,j);
b=best.binary(j);
if ((x==0)&(b==0))||((x==1)&(b==1))
delta=0; % delta为旋转角的大小
s=0; % s为旋转角的符号,即旋转方向
elseif (x==0)&(b==1)&(fitness(i)<best.fitness)
delta=0.01*pi;
if A*B>0
s=1;
elseif A*B<0
s=-1;
elseif A==0
s=0;
elseif B==0
s=sign(randn);
end
elseif (x==0)&(b==1)&(fitness(i)>=best.fitness)
delta=0.01*pi;
if A*B>0
s=-1;
elseif A*B<0
s=1;
elseif A==0
s=sign(randn);
elseif B==0
s=0;
end
elseif (x==1)&(b==0)&(fitness(i)<best.fitness)
delta=0.01*pi;
if A*B>0
s=-1;
elseif A*B<0
s=1;
elseif A==0
s=sign(randn);
elseif B==0
s=0;
end
elseif (x==1)&(b==0)&(fitness(i)>=best.fitness)
delta=0.01*pi;
if A*B>0
s=1;
elseif A*B<0
s=-1;
elseif A==0
s=0;
elseif B==0
s=sign(randn);
end
end
e=s*delta; % e为旋转角
U=[cos(e) -sin(e);sin(e) cos(e)]; % 量子旋转门
y=U*[A B]'; % y为更新后的量子位
chrom(2*i-1,j)=y(1);
chrom(2*i,j)=y(2);
end
end
03 | 量子遗传算法流程图
在了解上述关键步骤后,量子遗传算法流程图如下所示
主函数QuantumMain.m代码如下
clc;
clear all;
close all;
%----------------参数设置-----------------------
MAXGEN=200; % 最大遗传代数
sizepop=40; % 种群大小
lenchrom=[20 20]; % 每个变量的二进制长度
trace=zeros(1,MAXGEN);
%--------------------------------------------------------------------------
% 最佳个体 记录其适应度值、十进制值、二进制编码、量子比特编码
best=struct('fitness',0,'X',[],'binary',[],'chrom',[]);
%% 初始化种群
chrom=InitPop(sizepop*2,sum(lenchrom));
%% 对种群实施一次测量 得到二进制编码
binary=collapse(chrom);
%% 求种群个体的适应度值,和对应的十进制值
[fitness,X]=FitnessFunction(binary,lenchrom); % 使用目标函数计算适应度
%% 记录最佳个体到best
[best.fitness bestindex]=max(fitness); % 找出最大值
best.binary=binary(bestindex,:);
best.chrom=chrom([2*bestindex-1:2*bestindex],:);
best.X=X(bestindex,:);
trace(1)=best.fitness;
fprintf('%d\n',1)
%% 进化
for gen=2:MAXGEN
fprintf('%d\n',gen) %提示进化代数
%% 对种群实施一次测量
binary=collapse(chrom);
%% 计算适应度
[fitness,X]=FitnessFunction(binary,lenchrom);
%% 量子旋转门
chrom=Qgate(chrom,fitness,best,binary);
[newbestfitness,newbestindex]=max(fitness); % 找到最佳值
% 记录最佳个体到best
if newbestfitness>best.fitness
best.fitness=newbestfitness;
best.binary=binary(newbestindex,:);
best.chrom=chrom([2*newbestindex-1:2*newbestindex],:);
best.X=X(newbestindex,:);
end
trace(gen)=best.fitness;
end
%% 画进化曲线
plot(1:MAXGEN,trace);
title('进化过程');
xlabel('进化代数');
ylabel('每代的最佳适应度');
%% 显示优化结果
disp(['最优解X:',num2str(best.X)])
disp(['最大值Y:',num2str(best.fitness)]);
求解结果为x1=11.6281,x2=5.72503,最大值为17.3441