reliabilityweib.m

基于蒙特卡罗法的复杂系统可靠性仿真威布尔分布

%基于蒙特卡罗法的复杂系统可靠性仿真威布尔分布
clear
clc

%可靠性仿真
ybsl=1000;
for i=1:ybsl
%电源系统
%     T1=[exprnd(1200),exprnd(3000),exprnd(1550),normrnd(3000,220),exprnd(1400)];
%     T2=[exprnd(1200),exprnd(3000),exprnd(1550),normrnd(3000,220),exprnd(1400)];
%     TW1=exprnd(2000);
%     if TW1>=min(T1)
%        TZ=min(T1)+min(T2);
%     else
%         TZ=min(T1);
%     end
%     TK1=[exprnd(1600),exprnd(1600)];
%     TK2=[exprnd(1600),exprnd(1600)];
%     T3=[exprnd(1300),max(TK1),exprnd(1500),exprnd(1500),normrnd(2500,180)];
%     T4=[TZ,exprnd(1000),max(TK2),normrnd(2800,200)];
%     TW2=exprnd(1800);
%     if TW2>=min(T3)
%         T(i)=min(T3)+min(T4);
%     else
%         T(i)=min(T3);
%     end
    T1=[exprnd(1200),exprnd(3000),exprnd(1550),normrnd(3000,220),exprnd(1400)];
    T2=[exprnd(1200),exprnd(3000),exprnd(1550),normrnd(3000,220),exprnd(1400)];

    TZ=max([min(T1),min(T2)]);
    TK1=[exprnd(1600),exprnd(1600)];
    TK2=[exprnd(1600),exprnd(1600)];
    T3=[exprnd(1300),max(TK1),exprnd(1500),exprnd(1500),normrnd(2500,180)];
    T4=[TZ,exprnd(1000),max(TK2),normrnd(2800,200)];

    T(i)=max([min(T3),min(T4)]);
%二单元旁联系统
%     t1=exprnd(120);
%     t2=exprnd(100);
%     tw=exprnd(200);
%     if(tw>=t1)
%         T(i)=t1+t2;
%     else
%         T(i)=t1;
%     end
end

Y=1:ybsl;          %样本序号
XX=sort(T);
YY=Y/ybsl;
% ybjz=mean(XX);
% ybfc=std(XX);
%参数估计*******************************

[wphat,wpci] =weibfit(XX,0.05)

%初始值设定，n编码种群数 ,daishu迭代次数，pc交叉概率，pm变异概率，
n=100;
pc=0.6;
pm=0.001;


%指数分布*******************************************************************
m=20;                     %m编码长度 
x1=round(rand(n,m));
x2=round(rand(n,m));

%大循环开始*******************************************
eewc=1;
daishu=2;
pjz(1)=0;

while eewc>1.0000e-05
  for i=1:n
  y1(i)=0;
  y2(i)=0;
    for j=1:m   
    y1(i)=y1(i)+2^(m-j)*x1(i,j);
    y2(i)=y2(i)+2^(m-j)*x2(i,j);
    end
  end

yweib1=(wpci(2,1)-wpci(1,1))*y1/(2^m-1)+wpci(1,1);
yweib2=(wpci(2,2)-wpci(1,2))*y1/(2^m-1)+wpci(1,2);

%选择**************************
for i=1:n 
    f(i)=1/sum((weibcdf(XX,yweib1(i),yweib2(i))-YY).^2);
end

sumf=sum(f);

for i=1:n
    p=rand;
    j=1;
    fs=f(1);
     wheel=sumf*p;
    while ((fs<=wheel)&(j<=n))
        j=j+1;
        fs=fs+f(j);
    end
    sel1(i,:)=x1(j,:);
    sel2(i,:)=x2(j,:);
end 
%赌盘规则
%选择结束************************

%单点交叉************************
xx1=randperm(n);
xx2=randperm(n);
for i=1:2:n-1
    jc1=rand;
    if(jc1<pc)
       crossdot1=round((m-1)*rand)+1;
       for j=crossdot1:m
           X1=sel1(xx1(i),j);
           sel1(xx1(i),j)=sel1(xx1(i+1),j);   
           sel1(xx1(i+1),j)=X1;
       end
    end
    jc2=rand;
    if(jc2<pc)
       crossdot2=round((m-1)*rand)+1;
       for j=crossdot2:m
           X2=sel2(xx2(i),j);
           sel2(xx2(i),j)=sel2(xx2(i+1),j);   
           sel2(xx2(i+1),j)=X2;
       end
    end
end
%交叉结束************************

% 变异*************************
for i=1:n
    for j=1:m
        my1=rand;
        if my1<pm       
          if sel1(i,j)==1
             sel1(i,j)=0;
          else 
             sel1(i,j)=1;
          end
        end
        my2=rand;
        if my2<pm
          if sel2(i,j)==1
             sel2(i,j)=0;
          else 
             sel2(i,j)=1;
          end
        end
    end
end
%变异结束************************

%保存最优个体********************
fm=0;  
   for i=1:n
       if fm<f(i)
           fm=f(i);
           sel1(1,:)=x1(i,:); 
           sel2(1,:)=x2(i,:);
       end
   end
%保存最优个体********************

%绘图******************************
maxf(daishu)=0;  
   for i=1:n
       if maxf(daishu)<f(i)
           maxf(daishu)=f(i);
       end
   end
pjz(daishu)=sumf/n;

% hold on;
% plot(daishu,pjz(daishu),'b-');
% plot(daishu,maxf(daishu),'r.');  

eewc=abs(pjz(daishu)-pjz(daishu-1));
daishu=daishu+1;
%绘图******************************
x1=sel1;   
x2=sel2;
end
%大循环结束*******************************************

%最优点(yezy,zezy),函数值为fezy***************
weibulzy=0;
for i=1:n
    if weibulzy<f(i)
        weibulzy=f(i);
        yweib1zy=yweib1(i);
        yweib2zy=yweib2(i);
    end
end

yweib1zy
yweib2zy

%求得最优解*********求得最优解***************  

% %解析表达式
% t=0:20:1000;
% RR=exp(-t/120)+2.5*(exp(-t/100)-exp(-t/75));
% R1=1-weibcdf(t,yweib1zy,yweib2zy);
% R2=1-weibcdf(t,wphat(1),wphat(2));
% figure(1)
% hold on
% plot(t,RR,'r-');
% plot(t,R1,'b.');
% figure(2)
% hold on
% j=0;
% for i=1:40
%     j=j+25;
%     htXX(i)=XX(j);
% end
% htYY=25:25:ybsl;
% plot(t,RR,'r-');
% plot(htXX,1-htYY/ybsl,'bo')
% 
% figure(3)
% hold on
% plot(t,abs(RR-R1),'r-');
% plot(t,abs(RR-R2),'b-.');