请问下列优化问题为什么只迭代一次，我找了很久都找不到错误，希望了解fmincon的朋友能够批评指正，谢谢
function sansan()
clear;

options = optimset('LargeScale','off','Display','iter','TOLFUN',1e-30,'TOLX',1e-30);
x0=[50;1.35;1.37];   % Starting guess
lb=[1;0.95;0.97];
ub=[5;1.15;1.17];
[x,fval,exitflag,output]=fmincon(@myfun,x0,[],[],[],[],lb,ub,[],options) %后面两个向量是变量的LB和UB

function f=myfun(x)
u=1;
for b=1:1:50;
s=((b-1)*x(1)*1000):(b*x(1)*1000);
for m=0:100;
F1=((exp(-2.5*10^-5*s).* (1.896*10^-5*s).^m)/ gamma(m+1)).*normcdf((1250-(8.4823*10^-3*s+100*m))./sqrt(3.6*10^-7*s.^2+400*m)).^2.*normcdf((1270-(8.4823*10^-3*s+100*m))./sqrt(3.6*10^-7*s.^2+400*m)).^2;
Z1= trapz(s,F1);%积分
Zrecord1(m+1)=Z1;%用于记录每个m值所对应的积的值,之所以Zrecord后加1，因为matlab不允许序号为零的数组元素
end
Ztotal(u)=sum(Zrecord1);%分母在这
u=u+1;
end
for k=2:51  %时间间隔前的倍数
     for j=1:50   %打击次数
     a=-2.5*10^-5;
     b=1.896*10^-5;
     c=8.4823*10^-3;
     d=3.6*10^-7;
     K1a(j,k)=exp(a*(k-1).*x(1)*1000)* (b*(k-1).*x(1)*1000).^j./ gamma(j+1).* normcdf((1000*x(2)-(c*(k-1).*x(1)*1000+100*j))./sqrt(d*((k-1).*x(1)*1000)^2+400*j)).^2.*normcdf((1000*x(3)-(c*(k-1).*x(1)*1000+100*j))./sqrt(d*((k-1).*x(1)*1000)^2+400*j)).^2;
     K1b(j,k)=exp(a*(k-1).*x(1)*1000)* (b*(k-1).*x(1)*1000).^j./ gamma(j+1).* normcdf((1250-(c*(k-1).*x(1)*1000+100*j))./sqrt(d*((k-1).*x(1)*1000)^2+400*j)).^2.*normcdf((1270-(c*(k-1).*x(1)*1000+100*j))./sqrt(d*((k-1).*x(1)*1000)^2+400*j)).^2;
%给定t,在固定打击j下，t前倍数k的值
  end
  K2a(k)=1-sum(K1a(:,k));          %给定t,t前倍数k的值，j求和
  K2a(1)=0;
  K2b(k)=1-sum(K1b(:,k));          %给定t,t前倍数k的值，j求和
  K2b(1)=0;
end

for k=2:50;
K3a(k)=-x(1)*1000*K2a(k);
end
K3=sum(K3a)+50*x(1)*1000*K2a(51);  %t前倍数累加， EK

for k=2:50;
K4a(k)=-K2a(k);
end
K4=sum(K4a)+50*K2a(51);        %t前倍数累加， EN

for k=2:50;
K5(k)=-x(1)*1000*K2b(k).*K2b(k+1);
end                             %t前倍数累加， EP1
K5=sum(K5(k));
for k=2:50;
K6(k)=x(1)*1000*K2b(k).^2;
end
K6=sum(K6(k));
                            %t前倍数累加， EP2
for k=2:50;
K7(k)=K2b(k).*(x(1)*1000-Ztotal(k-1));
end
K7=sum(K7(k));                            %t前倍数累加， EP3
for k=2:50;
K8(k)=-K2b(k).*(x(1)*1000-Ztotal(k));    %t前倍数累加， EP4
end         %t前倍数累加， EP4
K8=sum(K8(k));
K9=K5+K6+K7+K8;       %downtime expected value
f=(2000*K4+5000*K9+27000)./K3;