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

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

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