前言
之前写过在python中的多目标遗传算法了,但是很可惜的是python调用商用有限元软件比较费事,需要在ironpython的编译环境下调用,然后这个ironpython它又是个老古董,不支持pandas和numpy这些python必备的第三方库(内心是崩溃的)。没办法,被迫向matlab低头。从头学起来吧,好在有了python的基础学matlab也是比较容易的。
代码实现
这个代码水平写得比较low,毕竟我刚学matlab,多多包含。相关的理论说明在我之前的博客中给出了,在此就不赘述了,理论说明贴个链接:
直接上代码了哈。首先是两个待优化的函数。
function y=fun1(x)
%optimal function one
y=(x-1)^2;
end
function y=fun2(x)
%optimal function two
y=cos(x);
end
接着就是算法进化过程了:
%%get initial parameter
NP=70;
L=10;
Pc=0.8;
Pm=0.2;
G=30; %迭代次数
Xs=pi/2;
Xx=0;
%% generate Primary population
f=randi([0,1],NP,L);
x=zeros(1,NP*2);
tempx=zeros(1,NP*2);
f1trace=[];
f2trace=[];
%% begin to Iteration
for i=1:1:G
disp([i]);
%pause(10);
f1=zeros(1,NP*2);
f2=zeros(1,NP*2);
nf=f;
for M=1:2:NP
p=rand();
if p<Pc
q=randi([0,1],1,L);
for j=1:1:L
if q(j)==1
tempsize=nf(M+1,j);
nf(M+1,j)=nf(M,j);
nf(M,j)=tempsize;
end
end
end
end
j=1;
while j<=NP*Pm
h=randi(NP);
for k=1:1:round(NP*Pm)
g=randi(L);
nf(h,g)=~nf(h,g);
end
j=j+1;
end
newf=[f;nf]; %垂直方向合并两个矩阵
tempsize=size(newf);
for j = 1:1:tempsize(1)
U=newf(j,:);
m=0;
for k =1:1:L
m=U(k)*2^(k-2)+m;
end
x(j)=Xx+m*(Xs-Xx)/(2^(L-1));
f1(j)=func(x(j));
f2(j)=func1(x(j)); %f1,f2都是预先分配好内存了的,避免了append的使用
end
lenfs=length(f1)*2;
fs=zeros(1,lenfs);
for j=1:2:lenfs
fs(j)=f1((j+1)/2);
fs(j+1)=f2((j+1)/2);
end
fs=reshape(fs,2,lenfs/2);
ps=zeros(length(f1),length(f1));
for k =1:1:length(f1)
for j = 1:1:length(f1)
if fs(1,k)<fs(1,j) && fs(2,k)<fs(2,j)
ps(k,j)=1;
elseif fs(1,k)<fs(1,j) && fs(2,k)==fs(2,j)
ps(k,j)=1;
elseif fs(1,k)==fs(1,j) && fs(2,k)<fs(2,j)
ps(k,j)=1;
end
end
end
jishu=zeros(1,length(f1)); %问题记录,把数值大的点筛出来了的,判断语句没问题,但是结果中
for k = 1:1:length(f1) %又出现数值大的点了,这些点从何处来?
tempa=0;
for j = 1:1:length(f1)
if ps(j,k)==1
tempa=tempa+1;
end
end
jishu(k)=tempa;
end
[num,index]=sort(jishu);
uninum=unique(num); %向量去重
front=cell(1,length(uninum));
for k =1:1:length(uninum)
temmat=[];
for j =1:1:length(index)
if num(j)==uninum(k)
temmat=[temmat,index(j)]; %记录对应索引
end
end
front{k}=temmat;
end
f=[];
for j = 1:1:length(front)
if length(f)==NP
break
end
tempf=[];
for k=1:1:length(front{j})
tempf=[tempf;newf(front{j}(k),:)];
end
tempf1=[];
tempf2=[];
[widtemf,wif_]=size(tempf);
for M = 1:1:widtemf
U=tempf(M,:);
m=0;
for k =1:1:L
m=m+U(k)*(2^(k-2));
end
xtemp=Xx+m*(Xs-Xx)/(2^(L-1));
tempf1=[tempf1,func(xtemp)];
tempf2=[tempf2,func1(xtemp)];
end
[temp1sort,temp1index]=sort(tempf1);
distance=zeros(1,length(tempf1));
distance(temp1index(1))=inf;
distance(temp1index(end))=-inf;
for k=1:1:length(temp1index)-1
if abs(distance(temp1index(k)))<10
distace(temp1index(k))=(tempf1(temp1index(k+1))-tempf1(temp1index(k-1)))/(max(tempf1)-min(tempf1));
else
continue
end
end
[temp2sort,temp2index]=sort(tempf2);
distance1=zeros(1,length(tempf2));
distance1(temp2index(1))=inf;
distance1(temp2index(end))=-inf;
for k=1:1:length(temp2index)-1
if abs(distance1(temp2index(k)))<10
distace1(temp2index(k))=(tempf2(temp2index(k+1))-tempf2(temp2index(k-1)))/(max(tempf2)-min(tempf2));
else
continue
end
end
sumdis=[];
for k =1:1:length(distance)
sumdistance=distance(k)+distance1(k);
sumdis=[sumdis,sumdistance];
end
[dis_,disindex]=sort(sumdis,'descend');
for k=1:1:length(sumdis)
f=[f;tempf(disindex(k),:)];
if length(f)==NP
break
end
end
% f=reshape(f,length(f)/L,L);
f1=[];
f2=[];
[lenf,wid_]=size(f);
for j=1:1:lenf
U=f(j,:);
m=0;
for k =1:1:L
m=m+U(k)*(2^(k-2));
end
x(j)=Xx+m*(Xs-Xx)/(2^(L-1));
f1=[f1,func(x(j))];
f2=[f2,func1(x(j))];
end
end
name=['F:/单轨吊驱动电机/maxwell调用/测试/pic/','pareto',num2str(i),'.png']; %保存图片到指定文件夹并命名
plot(f1,f2,'o','MarkerSize',7);
set(gca,'FontName','Times New Roman');
grid on
title('优化前沿曲线','FontName','宋体');
xlabel('f1value')
ylabel('f2value')
saveas(gcf, name);
end
最后的话
代码的确写得比较不堪入目,但是结果是好的,功能可以保证实现。
第一代优化结果:
第30代结果:
再多说两句,之前一直用python,就不想去学matlab,限制在自己的舒服区了,之后还是要多尝试尝试,不能一直止步不前。