求助，为什么我的梯度算法，迭代出来的结果不对

dudumaotc

这里是要计算一个网络的优化问题
其中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=[0 1 1 0 0 0 0;0 0 0 1 0 1 0;0 0 0 1 0 0 1;0 0 0 0 1 0 0;0 0 0 0 0 1 1;0 0 0 0 0 0 0;0 0 0 0 0 0 0];%链路存在为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 =[0,0.4611,0.3050,0,0,0,0;0,0,0,0.4983,0,0.6288,0;0,0,0,0.2140,0,0,0.0439;0,0,0,0,0.4399,0,0;0,0,0,0,0,0.6072,0.3127;0,0,0,0,0,0,0;0,0,0,0,0,0,0];

%变量定义
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;

sigma665

这个是子函数？

调用的时候没有出现语法错误吗？这种错的话，只有你自己知道了