今天實現了《一類求解方程全部根的改進差分進化算法》（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