最优化原理及应用

“最优化原理及应用”2008200388姚远1．用C语言，或者Matlab,或者Fortran等编写一个完整的SimulatedAnnealing算法和Genetic算法的优化程序。解:此题采用Matlab语言编写一个完整的SA算法优化程序。在该程序中选用的代价函数为：，初始的C0=1000，每一个阶段的Lk选为20，接受概率设为0.6，迭代的终止条件为e<0.00001(e=)。初始值的选取范围为，每次迭代的扰动=6。C=C0/k。的趋势如下列图所示：经过SA算法得出的结果为：x1=5.6682f(x)=-3.8854。程序如下：%退火算法clearallclcC0=1000;x0=20*rand(1,1)-10; %取初始值k=0;Lk=20;F=sin((x0-1.5)^2)+(x0-6)^2-3; %代价函数delta_x=6; %扰动e=1;epsilon=0.6; %接受概率i=1;while(e>0.00001)k=k+1;C=C0/k;for(i=1:Lk)w=2*rand(1,1)-1;x1=x0+w*delta_x; %产生一个x1F1=x1^4-x1^3-15*x1^2+1;delta_f=F1-F;e=abs(delta_f);if(F1<F)x0=x1;F=F1;X(i)=x1;i=i+1;elseprob=exp(-delta_f/C);if(prob>epsilon)x0=x1;F=F1;X(i)=i+1;endendendendx0F下面为采用遗传算法的优化程序。该程序中的代价函数与上面的SA算法所用的一致。设定变量的二进制码链长度为10，基因库中的二进制码链个数为10。自变量的取值区间为，设定遗传算法的迭代次数为500。在每次迭代保存基因的选择中分别采用了均值法和轮盘赌的方法，在保存、交换和异化的过程中每次都将最好的基因保存在基因库的最后一行，即精英操作。经过GA算法得出的结果为：x1=5.6696f(x)=-3.8851。从图中我们可以看出SA和GA算法均找到了该代价函数的最小值点。在运行的过程中GA的速度要明显快于SA。程序如下：%遗传算法%F=sin((x-1.5)^2)+(x-6)^2-3;最优化问题的代价函数%Q=1000-(sin((x-1.5)^2)+(x-6)^2-3);遗传算法中定义的fitnessclcclearallnum=10;length=10;total=2^(length)-1;min=-5;%取值区间的长度max=5;G=100;genome=round(rand(num,length));%定义自变量的随机二进制编码，长度为lengthfor(k=1:G)X=zeros(1,num);%二进制转化为十进制for(i=1:num)for(j=length:-1:1)X(i)=X(i)+genome(i,j)*2^(length-j);endX(i)=X(i)/total*(max-min)+min;endfor(i=1:num)%计算各变量的fitnessQ(i)=1000-(sin((X(i)-1.5)^2)+(X(i)-6)^2-3);end[order,location]=sort(Q);bestgene=location(num);%最好的基因位置BEST(k,:)=genome(bestgene,:);%最好的基因%%保存%q=sum(Q);%vector=(Q/q)*num;%vector=floor(vector);%将品质因数大于均值的gene选出来,保存到下次%j=1;%for(i=1:num)%if(vector(i)==1)%geneupdate(j,:)=genome(i,:);%j=j+1;%end%end%geneupdate(num,:)=BEST(k,:);%保存%使用轮盘赌的筛选方法q=sum(Q);selectvariable=ceil(q)*rand(1,1);j=1;while(1)a=0;i=1;a=a+Q(i);while(1)if(a<selectvariable)i=i+1;a=a+Q(i);elsegeneupdate(j,:)=genome(i,:);j=j+1;break;endendif(j==num)break;endendgeneupdate(num,:)=BEST(k,:);%交换kk=1;probc=0.5;for(i=1:2:num-1)variable1=rand(1,1);variable2=floor((length-1)*rand(1,1)+1);%找出截断点if(variable1<probc)for(j=variable2:length)a(kk)=geneupdate(i,j);b(kk)=geneupdate(i+1,j);kk=kk+1;endkk=1;for(j=variable2:length)geneupdate(i,j)=b(kk);geneupdate(i+1,j)=a(kk);kk=kk+1;endendendgeneupdate(num,:)=BEST(k,:);genome=geneupdate;%异化probt=0.1;for(i=1:num)for(n=1:length)variable3=rand(1,1);if(variable3<probt)if(genome(i,n)==1)genome(i,n)=0;elsegenome(i,n)=1;endendendendgeneupdate(num,:)=BEST(k,:);genome=geneupdate;endx=0;for(j=length:-1:1)x=x+genome(num,j)*2^(length-j);endx=x/total*(max-min)+min;xF=sin((x-1.5)^2)+(x-6)^2-3;F利用模拟退火和遗传算法求函数的最小值点，，x的精度控制在0.001范围内。解：在此题中所用的代价函数为：，初始的C0=1000，每一阶段的Lk取为20，接受概率设为0.6，迭代的终止条件为e<0.000001(e=)。初始值x0的选取范围为，每次迭代的扰动=2，C=C0/k。的趋势如下列图所示：经过SA算法得出的结果为：x1=0.2505f(x)=-1.1250。由原函数我们可以容易地得出：f(x)的最小值为0.25。经过SA得到的最优解与真实值的差值为0.0005，满足题目所给的精度要求。程序如下：%退火算法clearallclcC0=1000;x0=(2*rand(1,1)-1)*0.5;k=0;Lk=20;F=2*x0^2-x0-1;delta_x=2;e=1;epsilon=0.6;while(e>0.000001)k=k+1;C=C0/k;for(i=1:Lk)w=2*rand(1,1)-1;x1=x0+w*delta_x;F1=2*x1^2-x1-1;delta_f=F1-F;e=abs(delta_f);if(F1<F)x0=x1;F=F1;elseprob=exp(-delta_f/C);if(prob>epsilon)x0=x1;F=F1;endendendendx0FX=-1:0.1:1;y=2*X.^2-X-1;plot(X,y);gridon;set(gcf,'color','w');下面利用GA算法来计算f(x)的最小值，为了到达题目中给定的精度要求，这里设定二进制码链的长度为11，基因库中二进制码链的个数为10，自变量的取值区间定为：。遗传算法的迭代次数为100。在每次迭代保存基因的选择中分别采用了均值法和轮盘赌的方法，在保存、交换和异化的过程中每次都将最好的基因保存在基因库的最后一行，即精英操作。经过GA算法得出的结果为：x1=0.2496f(x)=-1.1250。满足题目中所给的精度要求。程序如下：%遗传算法%F=2*x^2-x-1;最优化问题的代价函数%Q=1000-(2*x^2-x-1);遗传算法中定义的fitnessclcclearallnum=10;length=11;total=2^(length)-1;min=-1;%取值区间的长度max=1;G=100;genome=round(rand(num,length));%定义自变量的随机二进制编码，长度为lengthfor(k=1:G)X=zeros(1,num);%二进制转化为十进制for(i=1:num)for(j=length:-1:1)X(i)=X(i)+genome(i,j)*2^(length-j);endX(i)=X(i)/total*(max-min)+min;endfor(i=1:num)%计算各变量的fitnessQ(i)=1000-(2*X(i)^2-X(i)-1);end[order,location]=sort(Q);bestgene=location(num);%最好的基因位置BEST(k,:)=genome(bestgene,:);%最好的基因%%保存%q=sum(Q);%vector=(Q/q)*num;%vector=floor(vector);%将品质因数大于均值的gene选出来,保存到下次%j=1;%for(i=1:num)%if(vector(i)==1)%geneupdate(j,:)=genome(i,:);%j=j+1;%end%end%geneupdate(num,:)=BEST(k,:);%保存%使用轮盘赌的筛选方法q=sum(Q);selectvariable=ceil(q)*rand(1,1);j=1;while(1)a=0;i=1;a=a+Q(i);while(1)if(a<selectvariable)i=i+1;a=a+Q(i);elsegeneupdate(j,:)=genome(i,:);j=j+1;break;endendif(j==num)break;endendgeneupdate(num,:)=BEST(k,:);%交换kk=1;probc=0.5;for(i=1:2:num-1)variable1=rand(1,1);variable2=floor((length-1)*rand(1,1)+1);if(variable1<probc)for(j=variable2:length)a(kk)=geneupdate(i,j);b(kk)=geneupdate(i+1,j);kk=kk+1;endkk=1;for(j=variable2:length)geneupdate(i,j)=b(kk);geneupdate(i+1,j)=a(kk);kk=kk+1;endendendgeneupdate(num,:)=BEST(k,:);genome=geneupdate;%异化probt=0.1;for(i=1:num)for(n=1:length)variable3=rand(1,1);%n=floor(7*rand(1,1)+1);if(variable3<probt)if(genome(i,n)==1)genome(i,n)=0;elsegenome(i,n)=1;endendendendgeneupdate(num,:)=BEST(k,:);genome=geneupdate;endx=0;for(j=length:-1:1)x=x+genome(num,j)*2^(length-j);endx=x/total*(max-min)+min;xF=2*x^2-x-1;F利用模拟退火和遗传算法求函数的最小值点，，x的精度控制在0.001范围内。解：首先采用SA算法，取C0=1000，每阶段的Lk=100，接受概率取为0.5，设定迭代次数K=100，为了加快收敛，令C0=C0/log(k)。x，y的扰动都定为0.1。从给定的代价函数我们容易看出其最小值点是〔0，0〕。经过SA运算可以得到优化后的结果为：x1=2.2114e-004，y1=3.1916e-004。x，y的收敛趋势如下列图所示：结果到达了题目中的精度要求。程序如下：%SA退火算法计算F=2*x^2+y^2的最小值clearallclcC0=1000;x0=0.5*(2*rand(1,1)-1);y0=0.5*(2*rand(1,1)-1);k=1;Lk=1000;F=2*x0^2+y0^2;Fdelta_x=0.1;delta_y=0.1;epsilon=0.8;K=100;I=1;for(j=1:K)k=k+1;C0=C0/log(k);for(i=1:Lk)x1=x0+(2*rand(1,1)-1)*delta_x;y1=y0+(2*rand(1,1)-1)*delta_y;F1=2*x1^2+y1^2;delta_f=F1-F;if(F1<F)x0=x1;y0=y1;F=F1;temp1(I)=x1;temp2(I)=y1;I=I+1;elseprob=exp(-delta_f/C0);if(prob>epsilon)x0=x1;y0=y1;F=F1;temp1(I)=x1;temp2(I)=y1;I=I+1;endendendendplot(temp1,'r');gridon;set(gcf,'color','w');holdonplot(temp2,'g');gridon;set(gcf,'color','w');x0y0F采用GA算法，为了到达题中所给的精度要求，所以取二进制码链的长度为22。其中前11位表示x，后11位表示y。基因库里二进制码链的个数为25个。遗传迭代200次。采用轮盘赌和均值筛选两种方式，进行精英操作。得出的结果为：x1=4.8852e-004，y1=4.8852e-004。将每次迭代的精英代码对应的x，y值画出来可以得到下列图：得到的结果到达了题目中所要求的精度要求。程序如下：%应用遗传算法计算F=2*x^2+y^2的最小值clcclearallnum=25;length_x=11;length_y=11;length=length_x+length_y;total_x=2^(length_x)-1;total_y=2^(length_y)-1;min=-1;%取值区间的长度max=1;G=200;k=1;genome=round(rand(num,length));%自变量定义的随机二进制编码，长度为lengthfor(i=1:G)%二进制转化为十进制X=zeros(1,num);Y=zeros(1,num);for(i=1:num)for(j=length_x:-1:1)X(i)=X(i)+genome(i,j)*2^(length_x-j);endX(i)=X(i)/total_x*(max-min)+min;endfor(i=1:num)for(j=length:-1:length-length_x+1)Y(i)=Y(i)+genome(i,j)*2^(length-j);endY(i)=Y(i)/total_y*(max-min)+min;endfor(i=1:num)%计算各变量的fitnessQ(i)=1000-(2*X(i)^2+Y(i)^2);end[order,location]=sort(Q);bestgene=location(num);%最好的基因位置BEST(k,:)=genome(bestgene,:);%最好的基因%保存q=sum(Q);vector=Q/q*num;vector=floor(vector);%将品质因数大于均值的gene选出来,保存到下次j=1;for(i=1:num)if(vector(i)==1)geneupdate(j,:)=genome(i,:);j=j+1;endendgeneupdate(num,:)=BEST(k,:);%交换kk=1;probc=0.5;for(i=1:2:num-1)variable1=rand(1,1);variable2=floor((length-1)*rand(1,1)+1);if(variable1<probc)for(j=variable2:length)a(kk)=geneupdate(i,j);b(kk)=geneupdate(i+1,j);kk=kk+1;endkk=1;for(j=variable2:length)geneupdate(i,j)=b(kk);geneupdate(i+1,j)=a(kk);kk=kk+1;endendendgeneupdate(num,:)=BEST(k,:);genome=geneupdate;%异化probt=0.2;for(i=1:num)for(n=1:length)variable3=rand(1,1);%n=floor(7*rand(1,1)+1);if(variable3<probt)if(genome(i,n)==1)genome(i,n)=0;elsegenome(i,n)=1;endendendendgeneupdate(num,:)=BEST(k,:);genome=geneupdate;k=k+1;endx=0;for(j=length_x:-1:1)x=x+genome(num,j)*2^(length_x-j);endx=x/total_x*(max-min)+min;xy=0;for(j=length:-1:length-length_x+1)y=y+genome(num,j)*2^(length-j);endy=y/total_x*(max-min)+min;y%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%X=zeros(1,G);Y=zeros(1,G);for(i=1:G)for(j=length_x:-1:1)X(i)=X(i)+BEST(i,j)*2^(length_x-j);endX(i)=X(i)/total_x*(max-min)+min;endfor(i=1:G)for(j=length:-1:length-length_x+1)Y(i)=Y(i)+BEST(i,j)*2^(length-j);endY(i)=Y(i)/total_y*(max-min)+min;endplot(X,'g');gridon;set(gcf,'color','w');holdonplot(Y,'r');gridon;set(gcf,'color','w');4．推销员问题是很有名的最正确化问题。假定一个推销员必须访问某区域里的N个村庄，他从第一号村庄出发走访其他N-1个村庄。如果村庄i和村庄j之间的距离为dij,那么他如何选择路线才能使他走的距离最近？给出其最正确化解决方法。解：在该题中假设有10个村庄，代价函数定为走遍所有村庄的距离和，固定一个起始点，运用SA算法对其进行优化。得到以下的结果：没有优化之前的路径：一号村庄的坐标为〔46，26〕。总长度为：554.6568优化之后的路径：总长度为：237.0437程序如下：%利用SA算法实现推销员问题clearallclcC0=500000000;Lk=2000;N=10;%position=round(100*rand(2,N)); %随机给出10个村庄的位置坐标position=[46,4,61,92,11,32,13,96,42,71;26,27,77,44,99,93,20,32,92,79;]; %给出10个村庄的固定位置F0=0;F1=0;k=1;epsilon=0.5;posorder0=1:N;posorder1=1:N;for(i=1:N-1)F0=F0+sqrt((position(1,posorder0(i+1))-position(1,posorder0(i)))^2+(position(2,posorder0(i+1))-position(2,posorder0(i)))^2);endF0 %代价函数%循环开始的地方while(1)k=k+1;C0=C0/log(k);if(C0<eps)break;endfor(i=1:Lk)posorder1=posorder0;%随机交换除第一个村庄之外的其余村庄的位置p=ceil((N-1)*rand(1,1)+1);while(p==1)%防止p=1，即防止交换第一个村庄p=ceil((N-1)*rand(1,1)+1);endq=ceil((N-1)*rand(1,1)+1);while(q==1)q=ceil((N-1)*rand(1,1)+1);endtemp=posorder1(p);posorder1(p)=posorder1(q);posorder1(q)=temp;F1=0;for(i=1:N-1)F1=F1+sqrt((position(1,posorder1(i+1))-position(1,posorder1(i)))^2+(position(2,posorder1(i+1))-position(2,posorder1(i)))^2);endif(F1<F0)F0=F1;posorder0=posorder1;elseifexp((F0-F1)/C0)>epsilonF0=F1;posorder0=posorder1;endendendendF0for(i=1:N)position1(:,i)=position(:,posorder0(i));endfigure(1)set(gcf,'color','w');plot(position1(1,:),position1(2,:));holdon;plot(position1(1,:),position1(2,:),'*');title=('优化后的行走路线');figure(2)set(gcf,'color','w');plot(position(1,:),position(2,:));holdonplot(position(1,:),position(2,:),'*');title=('没有经过优化的路线');5．如下递推公式：是有名的WidrowLMS算法递推公式。其中Wn是权向量，Xn是输入向量，μ是常数，ξn是误差函数。该公式的推导是一个典型的二次方程式加随机编程处理最正确化问题，请给出该公式的完整推导。〔参阅B.Widrow的’Adaptivesignalprocessing’教材和龚耀寰的“自适应原理及应用”或其他有关自适应信号处理书籍〕解：由各变量的定义，，最小均方误差MSE=E[]，，根据递推公式：，所以，将代入上式，并对两边求数学期望，可以得到：，又因为，，所以可以求得。6．基阵方向图设计是阵设计的根本问题，对于等间隔线列阵，已经证明切比雪夫加权是最正确权系数，其最正确化意义在于对于相同次瓣级，主瓣宽度最窄；而对于固定主瓣宽度次瓣级最低〔参见Urick“工程水声原理”〕。假设有一6元等间隔线列阵，试利用模拟退火或遗传算法求其最正确加权系数。并画出有关波束图*由切比雪夫加权法求得的最正确加权系数为[0.3,0.69,1,1,0.69,0.3],可利用该系数求出方向图函数G0(θ),θ为空间角，范围[-90,90],构成代价函数f=∑|[G(θ)－G0(θ)]|和品质函数F=C-f。解：应用SA退火算法求一个6元等间隔线列阵的最正确加权系数，假设6元的等间隔线列阵相