基于模拟落体碰撞的优化算法——本科毕业论文

之前本科做了个优化算法，文章也没发出去，这里就把算法贴出来。

1 算法核心思想

碰撞过程的机械能损耗导致落体运动走向结束。

2 方程

2.1 落体过程

简单，不再赘述。强调一点，利用时间定步长的算法，算法时间复杂度可能很高（未考虑归一化）。

2.2 碰撞过程

设定了个速度恢复系数，包括切向的和法向的，是对现实的简单抽象处理。本算法的核心部分——耗能。

2.3 碰壁

类似于碰撞，只不过不考虑速度损失，这是对所有限制条件的描述方式。

3 验证

结果：利用9个常用的验证函数进行了验证，结果正确，精度较低（相对于遗传算法），效率底下。

优势（相对于遗传算法）：对于波动很大而且很复杂的函数，同样得到局部最优解，本算法的结果明显优于简单遗传算法。

4 改进

I 对于碰撞，不计算地面的法向量，只通过各个自由度的一次比较，确定地面变化趋势。然后碰撞之后的方向变为所谓的“45度仰望天空”，即碰撞后方向固定为几个中选一个。

II 对于落体，采用变步长，时间步长的确定为当前“向上”的速度除以一个数（N，比如2），那么大概2N次计算能够判定是不是落地（每次都判断是否落地或者碰壁）。

5 改进结果

选取100个用于碰撞的个体进行两个自变量函数的测试，结果精度10^-5左右（高度上限1，下限0，另外两个自由度的长度级别基本和高度的相同），运行耗时和简单遗传算法处于相同状态。算法复杂度O（1）。

可能会在评论里给论文的网盘地址！

代码（MATLAB，针对的优化函数sin(x)*sin(y)/(x*y)   x,y [-10,10]）：

1 主函数（年代久远，参数记得个大概）

nind=input('The number of ping-pong:');%个体数目，建议100个
cts=input('The collision times:');%碰撞次数，建议100
g=input('The acceleration of gravity:');%重力加速度，建议1~10
ratio=input('The collision coefficient of kinetic energy:');%碰撞动能恢复系数，建议0.9以上
af=input('The amplification factor:');%目标函数于中间过程中的放大系数，此处为1比较好
value0=input('The initial height:');%置空高度，这里可以设为1，也可以更高
lb=[-10 -10]';
ub=[10 10]';
tic;
value0=af*value0;
position0=20*rand(2,nind)-10;
smallvalueindex=find(af*objfun(position0)>=value0);
position0(:,smallvalueindex)=[];
ppn=size(position0,2);%ping-pong个数
collidingtimes=0;%碰撞计数
v0=zeros(3,1);%初速度
positionall=zeros(2,ppn,cts);
valueall=zeros(ppn,cts);
for i=1:ppn
    position1=position0(:,i);
    value1=value0;
    v1=v0;
while collidingtimes<cts
    [position2,v2]=p2pdelh(position1,v1,value1,lb,ub,g,af);
    value1=objfun(position2);
    positionall(:,i,collidingtimes+1)=position2;
    valueall(i,collidingtimes+1)=value1;
    value1=af*value1;
    v3=v2vfun2(position2,v2,ratio);
    position1=position2;
    v1=v3;
    collidingtimes=collidingtimes+1;
end
collidingtimes=0;
end
positionall2=positionall(:,:,cts);
valuemin=min(valueall);
plot(1:cts,-valuemin,'k-');
xlabel('collision times');
ylabel('max value');
grid on;
title('Collision Procedure');
valuelast=valueall(:,cts);
[minmin index]=min(valuelast);
index=index(1);
minmin=-minmin;
toc;
fprintf('The number of ping-pong after filter according to the initial height is:%f.\n',ppn);
fprintf('The best position is:\n [%f %f].\n',positionall2(1,index),positionall2(2,index));
fprintf('The maxmum value is:%f.\n',minmin);
fprintf('This procedure spends %f seconds.\n',toc);

2 落体

function [pc,vc]=p2pdelh(position1,v1,value1,lb,ub,g,af)
%依据竖向速度变步长
ts=max([v1(3)/(2*g) 0.01]);
position1=position1+ts*v1(1:2);
value1=value1+v1(3)*ts-1/2*g*ts^2;
v1(3)=v1(3)-g*ts;
t=ts;
while value1>af*objfun(position1)
    indextoosmall=find(position1<lb);
    indextoobig=find(position1>ub);
    if ~isempty(indextoosmall)||~isempty(indextoobig)
        dn=zeros(3,1);
        dn(indextoosmall)=1;
        dn(indextoobig)=-1;
        dn=dn/sqrt(sum(dn.^2));
        v1n=sum(v1.*dn)*dn;
        v1=v1-2*v1n;
    end
    t=t+ts;
    position1=position1+ts*v1(1:2);
    value1=value1+v1(3)*ts-1/2*g*ts^2;
    v1(3)=v1(3)-g*ts;
end
pc=position1;
vc=v1;
end

3 碰撞

function v2 = v2vfun2(position,v1,ratio)
%Velocity which is a column vector changes when colliding.
%position is the colliding point.
%n is amplification factor.
fundiff=zeros(2,1);
for i=1:2
    position1=position;
    position1(i)=position(i)+0.01;
    fundiff(i)=objfun(position1)-objfun(position);
end
a= fundiff>=0;
dv=ones(3,1);
dv(a)=-1;
dv=1/sqrt(3)*dv;
nv=sqrt(ratio*sum(v1.^2)/3);
v2=dv*nv;
end