%改进模拟退火求解圆packing
function f=TM
%清屏
clc;
%清内存变量
clear;
circle=[0,0,10.0,100.0;
0,0,11.0,121.0;
0,0,12.0,144.0;
0,0,11.5,132.25;
0,0,9.5,90.25;
0,0,8.5,72.25;
0,0,10.5,110.25];
dim=size(circle,1);
maxradius=max(circle(:,3))*3;
r=-maxradius+2.*maxradius.*rand(dim,2);
r=[r,zeros(dim,2)];
circle=circle+r;
figure
rectangle('Position',[-maxradius,-maxradius,2*maxradius,2*maxradius],'Curvature',[1,1],'EdgeColor','b');
for k=1:7
rectangle('Position',[circle(k,1)-circle(k,3),circle(k,2)-circle(k,3),2*circle(k,3),2*circle(k,3)],'Curvature',[1,1],'EdgeColor','b');
text(circle(k,1),circle(k,2),num2str(k));
end
xlim([-maxradius,maxradius])
ylim([-maxradius,maxradius])
axis equal
fai=1e-5*3.4;
T_max=100;
T_min=0.01;
iter_max=10;
T=T_max;
a=1e-10;
R=36;
while T>=T_min
iter=1;
while iter
X0=circle;%记录生成初始格局
EX0=E(circle,dim,R,fai);%计算当前目标值E(X0)
b=(-maxradius+maxradius*rand(1,2));
circle(Ci,1:2)=b;
Xnew=circle(Ci,1:2);%计算每一个小圆的相对挤压痛苦度REDPi(i=1,2,3…n),其中最大痛苦度的圆记为Ci,
%按照置圆策略,将挑选出来的Ci的圆心放置在大圆布局的一个空点上,产生一个新的格局Xnew,
EXnew=E(circle,dim,R,fai);%计算新格局的能量EXnew
if EXnew
E(circle,dim,maxradius,fai)=EXnew;
circle=Xnew;
break;
elseif (exp((EXnew-EX0)/(T))>random(0,1));
EX0=EXnew;
X0=Xnew;
circle=X0;
else
EX0=E(circle,dim,R,fai);%计算目标值E(X0)
circle=X0;
end
iter=iter+1;
%每次停0.1秒
pause(0.1)
%清除上次画的图形
cla
%重绘图形
rectangle('Position',[-maxradius,-maxradius,2*maxradius,2*maxradius],'Curvature',[1,1],'EdgeColor','b');
for k=1:7
rectangle('Position',[circle(k,1)-circle(k,3),circle(k,2)-circle(k,3),2*circle(k,3),2*circle(k,3)],'Curvature',[1,1],'EdgeColor','b');
text(circle(k,1),circle(k,2),num2str(k));
disp(['迭代次数' num2str(iter) ',整体势能:' num2str(E(circle,dim,maxradius,fai)) ',静不平衡量:' num2str(J(circle,dim))]);
end
end
T=T*0.9;
disp(['迭代次数' num2str(iter) ',整体势能:' num2str(E(circle,dim,maxradius,fai)) ',静不平衡量:' num2str(J(circle,dim)) ',最大外包罗圆半径:' num2str(find_max_radius(circle,dim))]);
end
end
%求i,j两个圆之间的干涉量
function f=dij(circle,i,j)
f=circle(i,3) + circle(j,3) - sqrt( (circle(i,1)-circle(j,1)).^2 + (circle(i,2)-circle(j,2)).^2 );
if (f<0)
f=0;
end
end
%求第i个圆与外包罗圆之间的干涉量
function f=d0i(circle,i,R)
f=sqrt( circle(i,1).^2 + circle(i,2).^2 ) + circle(i,3) - R;
if (f<0)
f=0;
end
end
%整体势能
function f=E(circle,dim,R,fai)
f=0;
for i=1:dim
for j=1:dim
if(i~=i)
f=f+dij(circle,i,j).^2;
end
end
end
f=f./2;
for i=1:dim
f=f+d0i(circle,i,R).^2;
end
f=f+fai*(J(circle,dim).^2);
end
%求静不平衡量
function f=J(circle,dim)
tx=0;
ty=0;
for i=1:dim
tx=tx+circle(i,4).*circle(i,1);
ty=ty+circle(i,4).*circle(i,2);
end
tx=tx.^2;
ty=ty.^2;
f=sqrt(tx+ty);
end
%求最大包罗圆半径
function f=find_max_radius(circle,dim)
maxradius=dist(circle(1,1),circle(1,2),0,0)+circle(1,3);
for i=2:dim
t=dist(circle(i,1),circle(i,2),0,0)+circle(i,3);
if(t>maxradius)
maxradius=t;
end
end
f=maxradius;
end
%求两点间的距离
function f=dist(x1,y1,x2,y2)
f=sqrt((x1-x2).^2+(y1-y2).^2);
end
%标记最大痛苦度圆Ci
function f=Ci(circle,dim)
maxredp=0;
Ci=0;
for i=1:dim
r=Ui/(circle(i,3).^2);
if r>maxredp;
maxredp=r;
Ci=i;
end
f=Ci;
end
end
%第i个圆饼具有的挤压弹性势能
function f=Ui(circle,dim)
f=0;
for i=1:dim
for j=0:dim
if(i~=i)
f=f+dij(circle,i,j).^2;
end
end
end
end
主函数无法调用最后两个function,单步调试时发现在主函数中定义的dim,circle在子函数中无法使用。
错误提示为:??? Input argument "dim" is undefined.
Error in ==> TM>Ci at 136
for i=1:dim
Error in ==> TM at 41
circle(Ci,1:2)=b;
自己是matlab小白,靠百度代码意思学习,希望能得到帮助,谢谢