|
|
全国赛，本人用matlab求解过程中发现matlab所求的结果中没有满足约束条件，最近本人测试了大部分matlab优化函数，得到惊人的发现：
matlab在求解时约束条件得不到满足，当方程是病态时matlab很可能出错。
最近我查阅有关资料，本人试图利用构造罚函数来解决该问题，具体算法正在设想之中，同时本人编制了分支定界函数，不知哪位有兴趣一起研究。
以下是本人编的模拟退火法解决垃圾运输问题的程序：
clear
%+++++++++++++++++++++++++++++++++模拟退火算法+++++++++++++++++++++动态规划法++++++++++++++++++++++++++
tf=0.7;%终温
t0=10;%初始温度
t1=t0;
m=[1.5 1.5 0.55 1.2 1.30 0.85 1.2 2.3 1.4 1.5 1.1 2.7 1.8 1.8 1.4 1.5 0.8 1.5 0.8 0.6 1.3 1.8 1.4 1.6 1.6 1.0 2.0 1.0 2.1 1.2 1.9 1.2 1.6 1.2 1.5 1.3 0];
x=[ 3 1 5 4 3 0 7 9 10 14 17 14 12 10 19 2 6 11 15 7 17 21 27 15 15 20 21 24 25 28 5 22 25 9 9 30 0];
y=[2 5 4 7 11 8 9 6 2 0 3 6 9 12 9 16 18 17 12 14 16 0 9 19 14 17 13 20 16 18 12 5 7 20 15 12 0];
p1=0.4;
p2=1.8;
%最大运量
max_no=6;
%max_no=8;
fv=40;
init=100000;
n=length(x);
fstart=zeros(1,n);
path=zeros(1,n);
cost=0;
time=0;
xfstart=zeros(1,n);
xpath=zeros(1,n);
xcost=0;
xtime=0;
fmin=init;
fpath=zeros(1,n);
fst=zeros(1,n);
ftime=0
p_next=zeros(n,n)
g=zeros(1,n);
d=ones(n,n)*init;
l=10*n;
%+++++++++++++++++++++++++++建立连接图与路径++++++++++++++++++++++++++++++++++++
for i=1:n
for j=1:n
if i~=j&x(i)>=x(j)&y(i)>=y(j)
d(i,j)=abs(x(i)-x(j))+abs(y(i)-y(j));
end
end
end
%++++++++++++++++++++++++++++++以下用计算状态量+++++++++++++++++++++++++++++
for j=1:n
g(j)=(x(j)+y(j))*p2*m(j); %取第j点增加的费用
end
for i=1:n
[v,path(i)]=min(d(i, );
end
clear v;
for i=1:n
temp=0;
for j=1:n
if x(i)>=x(j)&y(i)>=y(j)&i~=j
temp=temp+1;
p_next(i,temp)=j;
end
end
end
%+++++++++++++++++++++++++++++++++++++++++打印图形++++++++++++++++++++++++++++++++++++++++++++++
%path=[37 37 37 37 37 37 4 3 1 37 9 8 7 31 13 6 16 35 37 37 25 10 33 18 19 21 37 26 27 29 5 37 32 17 20 23 37 ];
% for i=1:n
% if path(i)==37
% path(i)=i;
% end
% end
% for i=1:n-1
% a= [ x(i) , x(path(i))];
% b= [ y(i) , y(path(i))];
% hold on
% plot(a,b)
% grid on
% end
% for i=1:n
% hold on
% plot(x(i),y(i),'r*')
% grid on
% end
% break
%---------------------------------------------------------------------------------------------------------------------
for i=1:n
j=p_next(i,1);
k=1;
while j~=0
a= [ x(i) , x(p_next(i,k))];
b= [ y(i) , y(p_next(i,k))];
hold on
plot(a,b)
grid on
k=k+1;
if k==37 break;end
j=p_next(i,k);
end
end
for i=1:n
hold on
plot(x(i),y(i),'r*')
grid on
end
beark
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++计算初始值+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[path,fstart,time,cost]=myfun(d,m,g,path,p1,p2,max_no,x,y,fv)%调用动态规划函数；
%*************************************************************************************************************************************************************
% %打印出使图图形
% for i=1:n
% a= [ x(i) , x(path(i))];
% b= [ y(i) , y(path(i))];
% hold on
% plot(a,b)
% grid on
% end
% for i=1:n
% hold on
% plot(x(i),y(i),'r*')
% grid on
% end
% %******************************************************************************************************
while t1>tf
for k=1:l
% *************************************************随机交换*****************************************
xpath=change1(path,p_next,n)
%++++++++++++++++++++++++++++++++++++++++++++++++++++计算交换值+++++++++++++++++++++++++++++++++++++++++++++++++++++++
[xpath,xfstart,xtime,xcost]=myfun(d,m,g,xpath,p1,p2,max_no,x,y,fv)
if xcost==0
break
xcost=fmin;
end
%***********************************************************************************************************************
df=cost-xcost;
if df>0
path=xpath;
fstart=xfstart;
time=xtime;
cost=xcost;
if cost==0
break;
end
if cost<fmin
fmin=cost;
fpath=path;
ftime=time;
fst=fstart;
end
else
p=rand;%随机发生器
if p<exp(df/t1)
path=xpath;
fstart=xfstart;
time=xtime;
cost=xcost;
end;
end;break
end
t1=0.87*t1;
break
end
function [xpath,fstart,time,cost]=myfun(d,m,g,xpath,p1,p2,max_no,x,y,fv)
n=length(xpath);
fstart=zeros(1,n);
start=ones(1,n);
jj=0;temp=0;i=0;
f=zeros(n,n);
fw=zeros(n,n);
cost=0;dist=0;time=0.000;c=0;
ck=ones(1,n);count=0
for i=1:n-1
j =xpath(i);
start(j)=0;
end
%************************************************************************
while count<37 %全部分配完毕
temp=0;
t=n;
for k=1:n-1
if start(k)==1&ck(k)==1 %如是起始结点开始计算
if temp<d(k,n)
temp=d(k,n);
t=k;
end
end
end
i=t ;
if i==n
break;
end
jj=jj+1;
fstart(jj)=i;
f(i,i)=p1*d(i,n)+p2*d(i,n)*m(i);
count=count+1; % the start point and the end point cost 每一点到原点的距离空载时的费用
ck(i)=0;
fw(i,i)=m(i) ; %开始装载量
%下面是递规
k=xpath(i);
t=i;
while k~=n
if fw(i,t)+m(k) < max_no & ck(k)>0 %未装满,且当前还有垃圾
f(i,k)=f(i,t)+g(k); %动态规法求以i点为起始点以k点为最后铲点的运费
fw(i,k)=fw(i,t)+m(k) ; %在结点k时的载重
xpath(t)=k;
t=k;
ck(k)=0;
count=count+1;
start(k)=0;
end
k=xpath(k); %递规调用
end
f(i,n)=f(i,t);
xpath(t)=n;
%++++++++++++++++++++++++++
start=ones(1,n);
for i=1:n
j =xpath(i);
start(j)=0;%第j个不是起始点
end
end
%-------++++++++++++++++++++++---------------------------分配完毕------------------------------------------
for i=1:n
if fstart(i)~=0
tt=fstart(i);
cost=f(tt,n)+cost;
end
end
i=1;
while fstart(i)~=0
dist=x(fstart(i))+y(fstart(i))+dist;
i=i+1;
end
time=2*dist/40+(n-1)/6;
% %***********************************************************************
% i=1;value=0
% while fstart(i)~=0
% i
% c=fstart(i)
% value=value+(x(c)+y(c))*(p1+p2*m(i));
% while c~=n
% c= xpath(c)
% value=value+(x(c)+y(c))*p2*m(i);
% end
% i=i+1;
% end
function [xpath]=change1(path,p_next,n)
path=[37 37 37 37 37 37 4 3 1 37 9 8 7 31 13 6 16 35 37 37 25 10 33 18 19 21 37 26 27 29 5 37 32 17 20 23 37 ];
%path(24)=18;
for i=1:n
p=rand;
p=p*3
p=round(p)
if p==0
p=1;
end
[u,v]=sort(p_next(i,);
if u(p)~=0
path(i)=u(p);
end
end
path(37)=37;
xpath=path;
%******************************************************************
作者：逆风飞
QQ:176466480
CUMT 机自01-2班
[此贴子已经被作者于2003-12-6 22:50:57编辑过]
|
|
发表于 2003-12-6 20:02:35
楼主