第十二章：现代优化算法_x=sj(:,1:2:8);x=x(:)-CSDN博客

本文链接：https://blog.csdn.net/qq_41000485/article/details/99697398

12.1 模拟退火算法

12.1.1 算法简介

理论上，降温过程要足够缓慢，要使得在每一温度下达到热平衡
要确定在每一温度下状态转换的结束准则
选择初始温度和确定某个可行解的邻域的方法要适当

12.1.2 应用举例

解空间
目标函数
新解的产生
代价函数差
接受准则
降温
结束状态

clc, clear
sj0=load('sj.txt');    
x=sj0(:,[1:2:8]);x=x(:);
y=sj0(:,[2:2:8]);y=y(:);
sj=[x y]; d1=[70,40]; 
sj=[d1;sj;d1]; sj=sj*pi/180; %角度化成弧度
d=zeros(102); %举例矩阵d初始化
for i=1:101
   for j=i+1:102
d(i,j)=6370*acos(cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*cos(sj(j,2))+sin(sj(i,2))*sin(sj(j,2)));
   end
end
d=d+d';
path=[];long=inf; %巡航路径及长度初始化
rand('state',sum(clock));  %初始化随机发生器
for j=1:1000  %求较好的初始解
    path0=[1 1+randperm(100),102]; temp=0;
    for i=1:101
        temp=temp+d(path0(i),path0(i+1));
    end
    if temp<long
        path=path0; long=temp;
    end
end
e=0.1^30;L=20000;at=0.999;T=1;
for k=1:L  %退火过程
c=2+floor(100*rand(1,2));  %产生新解
c=sort(c); c1=c(1);c2=c(2);
  %计算代价函数值的增量
df=d(path(c1-1),path(c2))+d(path(c1),path(c2+1))-d(path(c1-1),path(c1))-d(path(c2),path(c2+1));
  if df<0 %接受准则
  path=[path(1:c1-1),path(c2:-1:c1),path(c2+1:102)]; long=long+df;
  elseif exp(-df/T)>rand
  path=[path(1:c1-1),path(c2:-1:c1),path(c2+1:102)]; long=long+df;
  end
  T=T*at;
   if T<e
       break;
   end
end
path, long %输出巡航路径及路径长度
xx=sj(path,1);yy=sj(path,2);
plot(xx,yy,'-*') %画出巡航路径

12.2 遗传算法

12.2.1 简介

根据具体问题确定可行解域，确定一种编码方法，能用数值串或字符串表示可行解域的每一解
对每一解应有一个度量好坏的依据，它用一函数表示，叫做适应函数，一般由目标函数构成
确定进化参数群体规模M，交叉概率p，变异概率 $p_m$ ，进化终止条件

12.2.2 模型及算法

编码策略
初始种群：改良圈算法
目标函数
交叉操作：单点交叉
变异操作
选择

clc,clear
sj0=load('sj.txt');       
x=sj0(:,1:2:8); x=x(:);
y=sj0(:,2:2:8); y=y(:);
sj=[x y]; d1=[70,40]; 
sj=[d1;sj;d1]; sj=sj*pi/180; 
d=zeros(102); %距离矩阵d的初始值
for i=1:101
  for j=i+1:102
  d(i,j)=6370*acos(cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*cos(sj(j,2))+sin(sj(i,2))*sin(sj(j,2)));
  end
end
d=d+d'; w=50; g=100; %w为种群的个数，g为进化的代数
rand('state',sum(clock)); 
for k=1:w  %通过改良圈算法选取初始种群
    c=randperm(100); %产生1，...，100的一个全排列  
    c1=[1,c+1,102]; %生成初始解
    for t=1:102 %该层循环是修改圈
        flag=0; %修改圈退出标志
    for m=1:100
      for n=m+2:101
        if d(c1(m),c1(n))+d(c1(m+1),c1(n+1))<d(c1(m),c1(m+1))+d(c1(n),c1(n+1))
           c1(m+1:n)=c1(n:-1:m+1);  flag=1; %修改圈
        end
      end
    end
   if flag==0
      J(k,c1)=1:102; break %记录下较好的解并退出当前层循环
   end
   end
end
J(:,1)=0; J=J/102; %把整数序列转换成[0,1]区间上的实数，即转换成染色体编码
for k=1:g  %¸该层循环进行遗传算法的操作
    A=J; %交配产生子代A的初始染色体
    c=randperm(w); %产生下面交叉操作的染色体对
    for i=1:2:w  
        F=2+floor(100*rand(1)); %产生交叉操作的地址
        temp=A(c(i),[F:102]); %中间变量的保存值
        A(c(i),[F:102])=A(c(i+1),[F:102]); %交叉操作
        A(c(i+1),F:102)=temp;  
    end
    by=[];  %为了防止下面产生空地址，这里先初始化
while ~length(by)
    by=find(rand(1,w)<0.1); %产生变异操作的地址
end
B=A(by,:); %计算变异操作的初始染色体
for j=1:length(by)
   bw=sort(2+floor(100*rand(1,3)));  %产生变异操作的3个地址
   B(j,:)=B(j,[1:bw(1)-1,bw(2)+1:bw(3),bw(1):bw(2),bw(3)+1:102]); %交换位置
end
   G=[J;A;B]; %父代和子代种群合在一起
   [SG,ind1]=sort(G,2); %把染色体翻译成1，...，102的序列ind1
   num=size(G,1); long=zeros(1,num); %路径长度的初始值
   for j=1:num
       for i=1:101
           long(j)=long(j)+d(ind1(j,i),ind1(j,i+1)); %计算每条路径长度
       end
   end
     [slong,ind2]=sort(long); %对路径长度从小到大排序
     J=G(ind2(1:w),:); %精选前w个较短的路径对应的染色体
end
path=ind1(ind2(1),:), flong=slong(1)  %解的路径及路径长度
xx=sj(path,1);yy=sj(path,2);
plot(xx,yy,'-o') %画出路径

12.3 改进的遗传算法

12.3.1 模型及算法

两点改进
7. 交叉操作：改进型交叉
8. 变异操作

tic  %计时开始
clc,clear
sj0=load('sj.txt');       
x=sj0(:,1:2:8); x=x(:);
y=sj0(:,2:2:8); y=y(:);
sj=[x y]; d1=[70,40]; 
sj=[d1;sj;d1]; sj=sj*pi/180;  
d=zeros(102); 
for i=1:101
  for j=i+1:102
  d(i,j)=6370*acos(cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*cos(sj(j,2))+sin(sj(i,2))*sin(sj(j,2)));
  end
end
d=d+d'; w=50; g=100; 
rand('state',sum(clock)); 
for k=1:w  
    c=randperm(100); 
    c1=[1,c+1,102]; 
    for t=1:102 
        flag=0; 
    for m=1:100
      for n=m+2:101
        if d(c1(m),c1(n))+d(c1(m+1),c1(n+1))<d(c1(m),c1(m+1))+d(c1(n),c1(n+1))
           c1(m+1:n)=c1(n:-1:m+1);  flag=1; 
        end
      end
    end
   if flag==0
      J(k,c1)=1:102; break %记录下较好的解并退出当前循环
   end
   end
end
J(:,1)=0; J=J/102; 
for k=1:g  
    A=J; 
    for i=1:2:w
        ch1(1)=rand; %混沌序列的初始值
        for j=2:50
            ch1(j)=4*ch1(j-1)*(1-ch1(j-1)); %产生混沌序列
        end
        ch1=2+floor(100*ch1); 
        temp=A(i,ch1); 
        A(i,ch1)=A(i+1,ch1); 
        A(i+1,ch1)=temp;
    end
    by=[];  
while ~length(by)
    by=find(rand(1,w)<0.1); 
end
num1=length(by); B=J(by,:); 
ch2=rand;  %产生混沌序列的初始值
for t=2:2*num1 
       ch2(t)=4*ch2(t-1)*(1-ch2(t-1)); %产生混沌序列
end
for j=1:num1
   bw=sort(2+floor(100*rand(1,2)));  
   B(j,bw)=ch2([j,j+1]);
end
   G=[J;A;B]; 
   [SG,ind1]=sort(G,2); 
   num2=size(G,1); long=zeros(1,num2); 
   for j=1:num2
       for i=1:101
           long(j)=long(j)+d(ind1(j,i),ind1(j,i+1)); 
       end
   end
     [slong,ind2]=sort(long); 
     J=G(ind2(1:w),:); 
end
path=ind1(ind2(1),:), flong=slong(1)  
toc  %计时结束
xx=sj(path,1);yy=sj(path,2);
plot(xx,yy,'-o')

12.4 Matlab遗传算法工具

12.4.1 命令行形式调用遗传算法ga

适应度函数

function y=ycfun1(x);   %x为行向量
c1=[2 3 1]; c2=[3 1 0];
y= c1* x' + c2* x'.^2; y=-y;

非线性约束函数

function [f,g]=ycfun2(x);
f=[x(1)+2*x(1)^2+x(2)+2*x(2)^2+x(3)-10
   x(1)+x(1)^2+x(2)+x(2)^2-x(3)-50
   2*x(1)+x(1)^2+2*x(2)+x(3)-40];
g=x(1)^2+x(3)-2;

主函数

clc, clear
a=[-1 -2 0;-1 0 0];b=[-1;0];
[x,y]=ga(@ycfun1,3,a,b,[],[],[],[],@ycfun2);
x, y=-y

12.4.2 通过GUI使用遗传算法

使用optimtool命令
输入适应度函数——欲求最小值的目标函数
变量个数——适应度函数的输入向量长度

12.4.3 直接搜索工具

function ex12_3
a=[-1 -2 0;-1 0 0];b=[-1;0];
[x,y]=patternsearch(@fun1,rand(1,3),a,b,[],[],[],[],@fun2);  %初始值必须为行向量
x,y=-y
%定义目标函数
function y=fun1(x);   %x为行向量
c1=[2 3 1]; c2=[3 1 0];
y= c1* x' + c2* x'.^2; y=-y;
%定义非线性约束函数
function [f,g]=fun2(x);
f=[x(1)+2*x(1)^2+x(2)+2*x(2)^2+x(3)-10
   x(1)+x(1)^2+x(2)+x(2)^2-x(3)-50
   2*x(1)+x(1)^2+2*x(2)+x(3)-40];
g=x(1)^2+x(3)-2;

（直接搜索算法和每次的计算结果也是不一样的）