求助，为什么我的梯度算法，迭代出来的结果不对
这里是要计算一个网络的优化问题其中xij(t),hij(t),pi(t)是变量，xij(t)是原始变量；hij(t),pi(t)是拉格朗日松弛因子，将有约束的问题转换成无约束
因为我刚刚接触matlab以前都没有用过，有点赶鸭子上架，第一次编写，可能写的很垃圾，高手们帮我看看
我也没有分成多个m文件，就是一个，希望快点实现功能而已，结果死活不对的。。。
我也没有设定迭代的最小值，我只是简单的希望能在迭代200次以后，看看结果，不过为什么我的结果是严格单调的呢，奇怪的。。
帮我看看吧，谢谢了
我写的很烂的源文件，帮我提提意见吧。。
function L=my_grad_mincost(x0,p0,h0)
% input-> x0=[; ; ;] p0=[; ;] h0=[; ; ;]
%循环参数和网络参数设置
format long
maxiter = 200; % 定义循环次数
m=1;
R(1,1)=2;R(2,1)=0;R(3,1)=0;R(4,1)=0;R(5,1)=0;R(6,1)=-2;R(7,1)=0;R(1,2)=2;R(2,2)=0;R(3,2)=0;R(4,2)=0;R(5,2)=0;R(6,2)=0;R(7,2)=-2;%节点到目的节点的流量
link=;%链路存在为1，不存在为0
link_des(:,:,1)=link;
link_des(3,7,1)=0;
link_des(5,7,1)=0;
link_des(:,:,2)=link;
link_des(2,6,2)=0;
link_des(5,6,2)=0;
%链路的代价权值
a =;
%变量定义
syms x121 x131 x241 x261 x341 x451 x561;
syms x122 x132 x242 x342 x372 x452 x572;
syms h121 h131 h241 h261 h341 h451 h561;
syms h122 h132 h242 h342 h372 h452 h572;
syms p11 p12 p21 p22 p31 p32 p41 p42 p51 p52 p61 p72;
x=sym(ones(7,7,2));
for ii=1:7
for jj=1:7
for kk=1:2
x(ii,jj,kk)=['x' num2str(ii) num2str(jj) num2str(kk)];
end
end
end
x=x.*link_des;
p=sym(ones(7,2));
for ii=1:7
for jj=1:2
p(ii,jj)=['p' num2str(ii) num2str(jj)];
end
end
p(6,2)=0;
p(7,1)=0;
h=sym(ones(7,7,2));
for ii=1:7
for jj=1:7
for kk=1:2
h(ii,jj,kk)=['h' num2str(ii) num2str(jj) num2str(kk)];
end
end
end
h=h.*link_des;
%目标函数
X = zeros(7,7);
X=x(:,:,1)+x(:,:,2);
X=X.*X;
for ii=1:7
for kk=1:2
Lcap(ii,kk)=p(ii,kk)*(sum(x(ii,:,kk))-sum(x(:,ii,kk))-R(ii,kk))
end
end
Lh=h.*x;
L2=Lh(:,:,1)+Lh(:,:,2);
L=(sum(sum(a.*X)))-sum(sum(Lcap))-(sum(sum(L2)));
%计算 h p x各个变量的梯度
grad_h=sym(zeros(7,7,2));
for ii=1:7
for jj=1:7
for kk=1:2
if h(ii,jj,kk)~=0
grad_h(ii,jj,kk)=diff(L,h(ii,jj,kk));
end
end
end
end
grad_p=sym(zeros(7,2));
for ii=1:7
for kk=1:2
if p(ii,kk)~=0
grad_p(ii,kk)=diff(L,p(ii,kk));
end
end
end
grad_x=sym(zeros(7,7,2));
for ii=1:7
for jj=1:7
for kk=1:2
if x(ii,jj,kk)~=0
grad_x(ii,jj,kk)=diff(L,x(ii,jj,kk));
end
end
end
end
%计算 h p x各个变量变化的步长
step_h=sym(zeros(7,7,2));
for ii=1:7
for jj=1:7
for kk=1:2
if h(ii,jj,kk)~=0
step_h(ii,jj,kk)=0.01;
end
end
end
end
step_p=sym(zeros(7,2));
for ii=1:7
for kk=1:2
if p(ii,kk)~=0
step_p(ii,kk)=0.01;
end
end
end
step_x=sym(zeros(7,7,2));
for ii=1:7
for jj=1:7
for kk=1:2
if x(ii,jj,kk)~=0
step_x(ii,jj,kk)=0.01;
end
end
end
end
%迭代循环过程
while(m<=maxiter)
gradp=subs(grad_p,{x121,x131,x241,x261,x341,x451,x561,x122,x132,x242,x342,x372,x452,x572},{x0(1,2,1),x0(1,3,1),x0(2,4,1),x0(2,6,1),x0(3,4,1),x0(4,5,1),x0(5,6,1),x0(1,2,2),x0(1,3,2),x0(2,4,2),x0(3,4,2),x0(3,7,2),x0(4,5,2),x0(5,7,2)});
gradh1=subs(grad_h(:,:,1),{x121,x131,x241,x261,x341,x451,x561,x122,x132,x242,x342,x372,x452,x572},{x0(1,2,1),x0(1,3,1),x0(2,4,1),x0(2,6,1),x0(3,4,1),x0(4,5,1),x0(5,6,1),x0(1,2,2),x0(1,3,2),x0(2,4,2),x0(3,4,2),x0(3,7,2),x0(4,5,2),x0(5,7,2)});
gradh2=subs(grad_h(:,:,2),{x121,x131,x241,x261,x341,x451,x561,x122,x132,x242,x342,x372,x452,x572},{x0(1,2,1),x0(1,3,1),x0(2,4,1),x0(2,6,1),x0(3,4,1),x0(4,5,1),x0(5,6,1),x0(1,2,2),x0(1,3,2),x0(2,4,2),x0(3,4,2),x0(3,7,2),x0(4,5,2),x0(5,7,2)});
gradh(:,:,1)=gradh1;
gradh(:,:,2)=gradh2;
gradx1=subs(grad_x(:,:,1),{x121,x131,x241,x261,x341,x451,x561,x122,x132,x242,x342,x372,x452,x572,h121,h131,h241,h261,h341,h451,h561,h122,h132,h242,h342,h372,h452,h572,p11,p12,p21,p22,p31,p32,p41,p42,p51,p52,p61,p72},{x0(1,2,1),x0(1,3,1),x0(2,4,1),x0(2,6,1),x0(3,4,1),x0(4,5,1),x0(5,6,1),x0(1,2,2),x0(1,3,2),x0(2,4,2),x0(3,4,2),x0(3,7,2),x0(4,5,2),x0(5,7,2),h0(1,2,1),h0(1,3,1),h0(2,4,1),h0(2,6,1),h0(3,4,1),h0(4,5,1),h0(5,6,1),h0(1,2,2),h0(1,3,2),h0(2,4,2),h0(3,4,2),h0(3,7,2),h0(4,5,2),h0(5,7,2),p0(1,1),p0(1,2),p0(2,1),p0(2,2),p0(3,1),p0(3,2),p0(4,1),p0(4,2),p0(5,1),p0(5,2),p0(6,1),p0(7,2)});
gradx2=subs(grad_x(:,:,2),{x121,x131,x241,x261,x341,x451,x561,x122,x132,x242,x342,x372,x452,x572,h121,h131,h241,h261,h341,h451,h561,h122,h132,h242,h342,h372,h452,h572,p11,p12,p21,p22,p31,p32,p41,p42,p51,p52,p61,p72},{x0(1,2,1),x0(1,3,1),x0(2,4,1),x0(2,6,1),x0(3,4,1),x0(4,5,1),x0(5,6,1),x0(1,2,2),x0(1,3,2),x0(2,4,2),x0(3,4,2),x0(3,7,2),x0(4,5,2),x0(5,7,2),h0(1,2,1),h0(1,3,1),h0(2,4,1),h0(2,6,1),h0(3,4,1),h0(4,5,1),h0(5,6,1),h0(1,2,2),h0(1,3,2),h0(2,4,2),h0(3,4,2),h0(3,7,2),h0(4,5,2),h0(5,7,2),p0(1,1),p0(1,2),p0(2,1),p0(2,2),p0(3,1),p0(3,2),p0(4,1),p0(4,2),p0(5,1),p0(5,2),p0(6,1),p0(7,2)});
gradx(:,:,1)=gradx1;
gradx(:,:,2)=gradx2;
stepp=subs(step_p,{p11,p12,p21,p22,p31,p32,p41,p42,p51,p52,p61,p72},{p0(1,1),p0(1,2),p0(2,1),p0(2,2),p0(3,1),p0(3,2),p0(4,1),p0(4,2),p0(5,1),p0(5,2),p0(6,1),p0(7,2)});
steph1=subs(step_h(:,:,1),{h121,h131,h241,h261,h341,h451,h561,h122,h132,h242,h342,h372,h452,h572},{h0(1,2,1),h0(1,3,1),h0(2,4,1),h0(2,6,1),h0(3,4,1),h0(4,5,1),h0(5,6,1),h0(1,2,2),h0(1,3,2),h0(2,4,2),h0(3,4,2),h0(3,7,2),h0(4,5,2),h0(5,7,2)});
steph2=subs(step_h(:,:,2),{h121,h131,h241,h261,h341,h451,h561,h122,h132,h242,h342,h372,h452,h572},{h0(1,2,1),h0(1,3,1),h0(2,4,1),h0(2,6,1),h0(3,4,1),h0(4,5,1),h0(5,6,1),h0(1,2,2),h0(1,3,2),h0(2,4,2),h0(3,4,2),h0(3,7,2),h0(4,5,2),h0(5,7,2)});
steph(:,:,1)=steph1;
steph(:,:,2)=steph2;
stepx1=subs(step_x(:,:,1),{x121,x131,x241,x261,x341,x451,x561,x122,x132,x242,x342,x372,x452,x572,h121,h131,h241,h261,h341,h451,h561,h122,h132,h242,h342,h372,h452,h572,p11,p12,p21,p22,p31,p32,p41,p42,p51,p52,p61,p72},{x0(1,2,1),x0(1,3,1),x0(2,4,1),x0(2,6,1),x0(3,4,1),x0(4,5,1),x0(5,6,1),x0(1,2,2),x0(1,3,2),x0(2,4,2),x0(3,4,2),x0(3,7,2),x0(4,5,2),x0(5,7,2),h0(1,2,1),h0(1,3,1),h0(2,4,1),h0(2,6,1),h0(3,4,1),h0(4,5,1),h0(5,6,1),h0(1,2,2),h0(1,3,2),h0(2,4,2),h0(3,4,2),h0(3,7,2),h0(4,5,2),h0(5,7,2),p0(1,1),p0(1,2),p0(2,1),p0(2,2),p0(3,1),p0(3,2),p0(4,1),p0(4,2),p0(5,1),p0(5,2),p0(6,1),p0(7,2)});
stepx2=subs(step_x(:,:,2),{x121,x131,x241,x261,x341,x451,x561,x122,x132,x242,x342,x372,x452,x572,h121,h131,h241,h261,h341,h451,h561,h122,h132,h242,h342,h372,h452,h572,p11,p12,p21,p22,p31,p32,p41,p42,p51,p52,p61,p72},{x0(1,2,1),x0(1,3,1),x0(2,4,1),x0(2,6,1),x0(3,4,1),x0(4,5,1),x0(5,6,1),x0(1,2,2),x0(1,3,2),x0(2,4,2),x0(3,4,2),x0(3,7,2),x0(4,5,2),x0(5,7,2),h0(1,2,1),h0(1,3,1),h0(2,4,1),h0(2,6,1),h0(3,4,1),h0(4,5,1),h0(5,6,1),h0(1,2,2),h0(1,3,2),h0(2,4,2),h0(3,4,2),h0(3,7,2),h0(4,5,2),h0(5,7,2),p0(1,1),p0(1,2),p0(2,1),p0(2,2),p0(3,1),p0(3,2),p0(4,1),p0(4,2),p0(5,1),p0(5,2),p0(6,1),p0(7,2)});
stepx(:,:,1)=stepx1;
stepx(:,:,2)=stepx2;
for ii=1:7
for kk=1:2
p(ii,kk) = p0(ii,kk)+stepp(ii,kk)*gradp(ii,kk);
end
end
for ii=1:7
for jj=1:7
for kk=1:2
if double(h0(ii,jj,kk)) > 0
h(ii,jj,kk) = h0(ii,jj,kk)+steph(ii,jj,kk)*gradh(ii,jj,kk);
else
if double(gradh(ii,jj,kk))>0
h(ii,jj,kk) = h0(ii,jj,kk)+steph(ii,jj,kk)*gradh(ii,jj,kk);
else
h(ii,jj,kk) = h0(ii,jj,kk);
end
end
end
end
end
for ii=1:7
for jj=1:7
for kk=1:2
x(ii,jj,kk) = x0(ii,jj,kk)+stepx(ii,jj,kk)*gradx(ii,jj,kk);
end
end
end
%此次更新的x,h,p值作为下次更新的初值x0,h0,p0
x0=x;
p0=p;
h0=h;
La=subs(L,{x121,x131,x241,x261,x341,x451,x561,x122,x132,x242,x342,x372,x452,x572,h121,h131,h241,h261,h341,h451,h561,h122,h132,h242,h342,h372,h452,h572,p11,p12,p21,p22,p31,p32,p41,p42,p51,p52,p61,p72},{x0(1,2,1),x0(1,3,1),x0(2,4,1),x0(2,6,1),x0(3,4,1),x0(4,5,1),x0(5,6,1),x0(1,2,2),x0(1,3,2),x0(2,4,2),x0(3,4,2),x0(3,7,2),x0(4,5,2),x0(5,7,2),h0(1,2,1),h0(1,3,1),h0(2,4,1),h0(2,6,1),h0(3,4,1),h0(4,5,1),h0(5,6,1),h0(1,2,2),h0(1,3,2),h0(2,4,2),h0(3,4,2),h0(3,7,2),h0(4,5,2),h0(5,7,2),p0(1,1),p0(1,2),p0(2,1),p0(2,2),p0(3,1),p0(3,2),p0(4,1),p0(4,2),p0(5,1),p0(5,2),p0(6,1),p0(7,2)});
iterstep(m,:)=La;
iterstep_x(:,:,:,m)=x0;
iterstep_gradp(:,:,m)=gradp;
iterstep_gradh(:,:,:,m)=gradh;
iterstep_gradx(:,:,:,m)=gradx;
m=m+1;
x=sym(ones(7,7,2));
for ii=1:7
for jj=1:7
for kk=1:2
x(ii,jj,kk)=['x' num2str(ii) num2str(jj) num2str(kk)];
end
end
end
x=x.*link_des;
p=sym(ones(7,2));
for ii=1:7
for jj=1:2
p(ii,jj)=['p' num2str(ii) num2str(jj)];
end
end
p(6,2)=0;
p(7,1)=0;
h=sym(ones(7,7,2));
for ii=1:7
for jj=1:7
for kk=1:2
h(ii,jj,kk)=['h' num2str(ii) num2str(jj) num2str(kk)];
end
end
end
h=h.*link_des;
end
save La
m;
iterstep;
回复 楼主 的帖子
这个是子函数？调用的时候没有出现语法错误吗？这种错的话，只有你自己知道了
页:
[1]