matlab差分算法
今天实现了《一类求解方程全部根的改进差分进化算法》（by 宁桂英，周永权），虽然最后的实现结果并没有文中分析的那么好，但是本文依然是给了一个求解多项式全部实根的基本思路。思路是对的，利用了代数原理。
求解全部根的理论还是很有必要说一下的。就是利用了多项式综合除法，在matlab中可以采用deconv(A,B)直接实现。同时为了确定多项式方程根的范围，还采用了代数方程根的分布理论，个人觉得这两点是值得借鉴的一种方法。
% 首先定义常量，包括最大迭代次数、搜索范围、个体维度、缩放因子等。程序如下
function DE
close all
clc
maxIteration=1000;%最大迭代次数
Generation=0;%进化代数，或者当前迭代代数
Xmax=30;%搜索上界，可以根据需要改为向量形式
Xmin=-30;%搜索下界
Dim=30;%个体维数
NP=50;%population size,种群规模
F=0.5;%scaling factor 缩放因子
CR=0.3;%crossover rate 交叉概率
FunIndex=4;%测试方程索引
mutationStrategy=1;%变异策略
crossStrategy=1;%交叉策略
% 步骤1：对应算法中Step 1，即初始化
X=(Xmax-Xmin)*rand(NP,Dim)+Xmin;%X：行代表个体i，列代表i的维度j
% 步骤2：对应算法中Step 2:
%step2 mutation,crossover,selection
while Generation<maxIteration
%求bestX
for i=1:NP
fitnessX(i)=testFun(X(i,:),FunIndex);%fitnessX表示X的适应值
end
[fitnessbestX,indexbestX]=min(fitnessX);%fitnessbestX最优适应值
bestX=X(indexbestX,:);%bestX表示最优值对应的位置
%%
%step2.1 mutation
%mutationStrategy=1：DE/rand/1,
%mutationStrategy=2：DE/best/1,
%mutationStrategy=3：DE/rand-to-best/1,
%mutationStrategy=4：DE/best/2,
%mutationStrategy=5：DE/rand/2,
%产生为每一个个体Xi,G 产生一个变异向量Vi,G。 G代表进化代数
V=mutation(X,bestX,F,mutationStrategy);
%%
%step2.2 crossover
%crossStrategy=1:binomial crossover
%crossStrategy=2:Exponential crossover
%产生为每一个个体Xi,G 产生一个交叉向量Ui,G。 G代表进化代数
U=crossover(X,V,CR,crossStrategy);
%%
%step2.3 selection
for i=1:NP
fitnessU(i)=testFun(U(i,:),FunIndex);
if fitnessU(i)<=fitnessX(i)
X(i,:)=U(i,:);
fitnessX(i)=fitnessU(i);
if fitnessU(i)<fitnessbestX
bestX=U(i,:);
fitnessbestX=fitnessU(i);
end
end
end
%%
Generation=Generation+1;
bestfitnessG(Generation)=fitnessbestX;
end
% 步骤3：显示结果
%画图
%plot(bestfitnessG);
optValue=num2str(fitnessbestX);
Location=num2str(bestX);
disp(strcat('the optimal value','=',optValue));
disp(strcat('the best location','=',Location));
end
%变异向量用函数mutation(X,bestX,F,mutationStrategy)
function V=mutation(X,bestX,F,mutationStrategy)
NP=length(X);
for i=1:NP
%在[1 NP]中产生nrandI个互不相等的随机数，且与i皆不相等
nrandI=5;
r=randi([1,NP],1,nrandI);%产生一个1*nrandI的伪随机标准分布，范围是1：NP
for j=1:nrandI %保证随机数互不相等。如果出现随机数相等的情况，则sum(equlr)>nrandi;
equalr(j)=sum(r==r(j));
end
equali=sum(r==i); %保证产生的随机数与i不相等,如果相等的话，则equali>0;则两者合并在一起可得，当且仅当sum(equalr)+equali=nradi时，产生的随机数是符合条件；任一条件不满足，则sum(equalr)+equali》nrandi
equalval=sum(equalr)+equali;
while(equalval>nrandI)
r=randi([1,NP],1,nrandI);
for j=1:nrandI
equalr(j)=sum(r==r(j));
end
equali=sum(r==i);
equalval=sum(equalr)+equali;
end
switch mutationStrategy
case 1
%mutationStrategy=1:DE/rand/1;
V(i,:)=X(r(1),:)+F*(X(r(2),:)-X(r(3),:));
case 2
%mutationStrategy=2:DE/best/1;
V(i,:)=bestX+F*(X(r(1),:)-X(r(2),:));
case 3
%mutationStrategy=3:DE/rand-to-best/1;
V(i,:)=X(i,:)+F*(bestX-X(i,:))+F*(X(r(1),:)-X(r(2),:));
case 4
%mutationStrategy=4:DE/best/2;
V(i,:)=bestX+F*(X(r(1),:)-X(r(2),:))+F*(X(r(3),:)-X(r(4),:));
case 5
%mutationStrategy=5:DE/rand/2;
V(i,:)=X(r(1),:)+F*(X(r(2),:)-X(r(3),:))+F*(X(r(4),:)-X(r(5),:));
otherwise
error('没有所指定的变异策略，请重新设定mutationStrategy的值');
end
end
end
%交叉函数，根据算法中提供的两种交叉方法编写，即binomial crossover和 %二项式交叉和指数交叉
%Exponential crossover
function U=crossover(X,V,CR,crossStrategy)
%V为产生的变异向量
[NP,Dim]=size(X);
switch crossStrategy
%crossStrategy=1:binomial crossover
case 1
for i=1:NP
jRand=floor(rand*Dim);%由于jRand要在[0,1)*Dim中取值，故而用floor
for j=1:Dim
if rand<CR||j==jRand
U(i,j)=V(i,j);
else
U(i,j)=X(i,j);
end
end
end
%crossStrategy=2:Exponential crossover
case 2
for i=1:NP
j=floor(rand*Dim);%由于j在[0,1)*Dim中取值，故而用floor
L=0;
U=X;
U(i,j)=V(i,j);
j=mod(j+1,D);
L=L+1;
while(rand<CR&&L<Dim)
U(i,j)=V(i,j);
j=mod(j+1,D);
L=L+1;
end
end
otherwise
error('没有所指定的交叉策略，请重新设定crossStrategy的值');
end
end
%测试函数，可以根据需要添加相应的程序
function y=testFun(x,index)
%x代表参数，index代表测试的函数的选择
%该测试函数为通用测试函数，可以移植
%目录
% 函数名 位置 最优值
%1.Sphere 0 0
%2.Camel 多个
%3.Rosenbrock
switch index
case 1 %Sphere函数
y=sum(x.^2);
case 2 %Camel函数,Dim只能取2
if length(x)>2
error('x的维度超出了2');
end
xx=x(1);yy=x(2);y=(4-2.1*xx^2+xx^4/3)*xx^2+xx*yy+(-4+4*yy^2)*yy^2;
case 3 %Rosenbrock函数
y=0;
for i=2:length(x)
y=y+100*(x(i)-x(i-1)^2)^2+(x(i-1)-1)^2;
end
case 4
y=sum((x-1).^2);
otherwise
disp('no such function, please choose another');
end
end
% %%该函数用修正的差分进化算法求解多项式的全部根
function DE_to_polynomial
clc
close all
popsize=20;%种群规模
% F=0.5;%缩放因子
CR=0.5;%交叉概率 控制参数是三个
Gnow=1; %种群代数
Gmax=500;
Dim=1;
% Xmax=10;
% Xmin=-10;
R=11;
A=[1 0 1 -10 -1 0 -1 10];
len=length(A);
for i=len:-1:1
r(i)=R.^(len-i);
end
A=A.*r;
nA=length(A);
JGoal=zeros(nA,1);
XGoal=zeros(nA,1);
beta=0.3;
for n=1:nA
% step2: initialize
Gnow=1;
a=2*rand(popsize,Dim)-1;
b=2*rand(popsize,Dim)-1;
X=complex(a,b);
if n>1
B=[1,XGoal(n,1)];
A=deconv(A,B);
end
while Gnow<Gmax
lamda=Gmax/(Gnow+Gmax);
sita=beta*(exp(lamda)-1); %这里采用的自适应变异算子
% step3: choose the well-fitness seeds 即1/2原则
if Gnow==1
Jt1=fitness_sort(X,A);
end
% step 4:mutation operation
L=mutation(X(1:popsize/2,:),sita);
%step 5:cross
V=crossover(X,L,CR);
%step6:choose
X=selection(X,V,A);
% 步骤7：终止检验
Jt2=fitness_sort(X,A);
eps=1e-5;
Jbest(Gnow)=Jt2(1);
fitbestX(Gnow)=X(1);
if Jbest(Gnow)<eps
JGoal(n,1)=Jbest(Gnow);
XGoal(n,1)=fitbestX(Gnow);
break;
end
JGoal(n,1)=Jbest(Gnow);
XGoal(n,1)=fitbestX(Gnow);
Gnow=Gnow+1;
end
end
JGoal
XGoal*11
end
%% 变异操作
function L=mutation(X,sita)
Xbest=X(1,:);%由于进行排序之后，较小适应度的X排在前面
[NP,Dim]=size(X);
for i=1:NP
%接下来只需要产生si个不相同的r ，并且r与i不相等
numr=4;
r=randi([1,NP],1,numr);%产生一个1*nrandI的伪随机标准分布，范围是1：NP
for j=1:numr
equalr(j)= sum(r==r(j));
end
equali=sum(r==i);
equalvalue=sum(equalr(j))+equali;
while equalvalue>numr
r=randi([1,NP],1,numr);%产生一个1*nrandI的伪随机标准分布，范围是1：NP
for j=1:numr
equalr(j)= sum(r==r(j));
end
equali=sum(r==i);
equalvalue=sum(equalr(j))+equali;
end
% 变异策略,将popsize/2个个体生成为popsize个个体
u(i,:)=Xbest+sita*(X(r(1),:)-X(r(2),:));
w(i,:)=(X(r(3),:)-X(r(4),:))/2;
end
k=1:NP;
L(k,:)=u(k,:);
k=NP+1:2*NP;
L(k,:)=w(k-NP,:);
end
%% 交叉操作
function V=crossover(X,L,CR)
[popsize,Dim]=size(X);
for i=1:popsize
jrand=floor(rand*Dim);
for j=1:Dim
if rand<CR ||j==jrand
V(i,j)=L(i,j);
else
V(i,j)=X(i,j);
end
end
end
end
%% 选择操作
function T=selection(X,V,A)
[popsize,Dim]=size(X);
for i=1:popsize
Jv=fitness(V(i,:),A);
Jx=fitness(X(i,:),A);
if Jv<Jx
T(i,:)=V(i,:);
else
T(i,:)=X(i,:);
end
end
end
% %% 步骤7：终止检验
% function termination(X)
% eps=1e-6;%终止此次求根的精度要求
% if
%
%
%
% end
%% 用于测试的多项式方程
function y=testfunc(x,A)
% 测试向量
% A=[1 0 1 -10 -1 0 -1 10];
len=length(A);
for i=len:-1:1
vect(i)=x.^i-1;
end
y=A*vect';
end
function y=func(x,A)
R=11;%多项式根的变化范围在R的圆内
% A=[1 4 1 -10 -1 0 -1 10]/R;
% syms x;
len=length(A);
for i=len:-1:1
vect(i)=x.^(len-i);
end
y=A*vect';
end
%% 计算各个函数的适应度
function J=fitness(x,A)
y=func(x,A);
J=abs(y);
end
%% 计算种群中每个个体的适应度并排序，利用二分之一规则选取个体，其中J 值越小，说明适应度越好
function J=fitness_sort(X,A)
[popsize,Dim]=size(X);
for i=1:popsize
J(i)=fitness(X(i),A);
end
[J,index]=sort(J);
X=X(index);%按照升序排列的X
end
