clear all;clc
%输入磨损数据 DeltaV
load('abrasion_loss.mat')
guancexulie=[]
for i=1:120
if(0<=DeltaV(1,i)&&DeltaV(1,i)<=0.01275)
guancexulie(i)=1;
elseif (0.01275<DeltaV(1,i)&&DeltaV(1,i)<=0.02545)
guancexulie(i)=2;
elseif (0.02545<DeltaV(1,i)&&DeltaV(1,i)<=0.048915)
guancexulie(i)=3;
elseif (0.048915<DeltaV(1,i)&&DeltaV(1,i)<=0.06815)
guancexulie(i)=4;
else (0.06815<DeltaV(1,i)&&DeltaV(1,i)<=1)
guancexulie(i)=5;
end
end
% 状态转移矩阵
A=[0.6 0.2 0.1 0.1;
0 0.6 0.2 0.2;
0 0 0.8 0.2;
0 0 0 1];
% 输出观测矩阵
B=[0.8 0.15 0.03 0.01 0.01;
0 0.8 0.15 0.025 0.025;
0 0 0.8 0.15 0.05;
0 0 0 0.8 0.2];
% 输出概率
Pi=[1 0 0 0];
Alpha=zeros(4,120);
Beta=zeros(4,120);
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%发射概率
P=[];
for i=1:4
for t=1:120
c=guancexulie(1,t);
P(i,t)=B(i,c);
end
end
%前向概率
for i=1:4 %初始化Alpha的第一列
Alpha(i,1)=Pi(i)*P(i,1);
end
for t=1:119
for j = 1:4
sum1=0;
for i =1:4
sum1=sum1+Alpha(i,t)*A(i,j);
end
Alpha(j,t+1)=sum1*P(j,t+1);
end
end
%后向概率
for i=1:4
Beta(i,120)=1;%初始化
end
%计算Beta
for t=119:-1:1
for i=1:4
sum2=0;
for j=1:4
sum2=sum2+P(j,t+1)*Beta(j,t+1)*A(i,j);
end
Beta(i,t)=sum2;
end
end
houxianggailv=0;
for i=1:4
houxianggailv=houxianggailv+Pi(i)*P(i,1)*Beta(i,1);
end
Alpha
Beta
L=Alpha(1,120)+Alpha(2,120)+Alpha(3,120)+Alpha(4,120)
T=120;
m=4;
%其他参数
Gamma=zeros(m,T);
g=zeros(1,m);
Epsilon=zeros(m,m,T-1);
H=zeros(m,m);
%计算Gamma/g
for i=1:m
for t=1:T
Gamma(i,t)=Alpha(i,t)*Beta(i,t)/L;
end
end
%for i=1:m
% sum1=0;
% for t=1:T-1
% sum1=sum1+Gamma(i,t);
% end
% g(i)=sum1;
%end
for i=1:m
for t=1:T-1
g(i)=g(i)+Gamma(i,t);
end
end
%计算Epsilon/H
for t=1:T-1
for i=1:m
for j=1:m
Epsilon(i,j,t)=Alpha(i,t)*A(i,j)*Beta(j,t+1)*P(j,t+1)/L;
end
end
end
for i=1:m
for j=1:m
for t=1:T-1;
H(i,j)=H(i,j)+Epsilon(i,j,t);
end
end
end
Gamma
g
Epsilon
H
Pi_new=Gamma(:,1)'
O = guancexulie;
A_size = size(A) ;
B_size = size(B) ;
O_size = size(O) ;
N = A_size(1,1) ;%状态集个数
M = A_size(1,2) ;
Y = B_size(1,2); %B的列数，即观测集个数
K = O_size(1,2);
A_1 = zeros();
f = sum(Gamma,2);%计算Gamma的每一行之和，表示在观测时间序列的所有时刻都处于某一个状态的概率和
for i = 1:M
for j = 1:N
A_1(i,j) = sum(sum(Epsilon(i,j,:))) / (f(i,1) - Gamma(i,120));
end
end
for i=1:m
for j=1:m
A_new(i,j)=H(i,j)/g(i);
end
end
Pi_1 = zeros();
for i = 1:M
Pi_1(i,1) = Gamma(i,1);
end
B_1 = zeros();
G = zeros();
c = 0;
for y = 1:Y % Y输出集个数，即红、白
for j= 1:N % N状态个数
for t = 1:K % K输出个数
if O(1,t) == y % 用于比较输出的观测值是否由对应的状态值输出，即条件ot = vk;
c = c + 1;
G(c,1) = Gamma(j,t);
end
end
B_1(j,y) = sum(G) / f(j,1); % f为Gamma每一行的和 即我定义的g
c = 0;
G = zeros();
end
end
A_1
A=A_1
Pi=Pi_new
%A_new =
%
% 0.9512 0.0488 0 0
% 0 0.9592 0.0408 0
% 0 0 0.9633 0.0458
% 0 0 0 1.0000
%%A
n = 100;
Q =A-eye(4);
[l,l] = size(Q);
A = eye(l) + 10 * Q;
[l,m] = size(A);%+++++
a = [1 zeros(1,m-1)]; %++++++
t = 0:0.5:n;
A
R_lisan = lisan_kkd_jisuan( a,A_1,n );
R_lianxu = lianxu_kkd_jisuan( Q,t );
load('mckkd.mat');
load('lskkd');
% plot(0:0.1:10,R_lisan,'k:',t/10,R_lianxu,'b-',0:0.1:9.8,KKD,'v','MarkerEdgeColor','r','MarkerFaceColor','r','MarkerFaceColor','r');
%
% legend('Discontinue method [X]','Continue method [X]','Monte Carlo method');
% legend boxoff;
% axis([0 n/10 0 1.0]);
% xlabel('t','FontName','Times New Roman','FontSize',10);
% ylabel('Reliability','FontName','Times New Roman','FontSize',10);
% set(gca,'ytick',[0:0.2:1.0]);
% set(gca,'FontName','Times New Roman','FontSize',10,'LineWidth',1);
% set(gcf,'Units','centimeters','Position',[10 7 12 6.75]);
% set(gca,'box','off');
%%
%复现可靠度计算
fuxiankkd=[];
Al=[];
for h=1:100
Pi=Pi*A_1^(h-1)
for i=1:4 %初始化Alpha的第一列%
Al(i,1)=Pi(i)*P(i,h);
end
for t=1:20
for j = 1:4
sum1=0;
for i =1:4
sum1=sum1+Al(i,t)*A_1(i,j);
end
Al(j,t+1)=sum1*P(j,t+1);
end
end
fuxiankkd(h)=sum(Al(:,20))
end
kkkkkkkkkkkkkdddddd=log(fuxiankkd)
A
DELTA=Stationary_Distribution(A)
%%
plot(0:1:99,kkkkkkkkkkkkkdddddd,'k:','MarkerEdgeColor','r','MarkerFaceColor','r','MarkerFaceColor','r');
评论0