end

end

%第四步

D(n)=0;

for i=1:r

for j=1:s

D(n)=D(n)+Pa(i)*Pba(i,j,n)*d(i,j); % 隶D(n)

end

end

R(n)=0;

for i=1:r

for j=1:s

if(Pba(i,j,n)~=0)

R(n)=R(n)+Pa(i)*Pba(i,j,n)*log2(Pba(i,j,n)/Pb(j,n)); %求R(n)

end

end

end

%判断差别是否在允许的精度范围之内

if(abs(R(n)-R(n-1))<=10^(-7)) %R(n)精度判断

if(abs(D(n)-D(n-1))<=10^(-7)) %D(n)精度判断

break;

end

end

n=n+1; %内层循环次数加1

end %内层循环结束

%第五步

S(m+1)=S(m)+0.5;

%第六步

if(abs(R(n)<10^(-7)))

if(m<=10)

disp('此时S的值为：'),disp(S(m));

error('初始拉氏乘子S取得大了，请取小些！');

%判断拉氏乘子S的初始值是否合适

else

[k,l,q]=size(pba);

Pba=pba(:,:,q);

Rmin=min(RS)

Dmax=max(DS);

Smax=S(m-1);

break;

end

end