DISTRIBUCION DE WEIBULL POR EL METODO DE MATLAB
Es el nombre abreviado de “MATriz LABoratory”. Es un programa para realizar cálculos
numéricos con vectores y matrices, y por tanto se puede trabajar también con números
escalares (tanto reales como complejos), con cadenas de caracteres y con otras estructuras de
información más complejas.
Matlab es un lenguaje de alto rendimiento para cálculos técnicos, es al mismo tiempo un
entorno y un lenguaje de programación. Uno de sus puntos fuertes es que permite construir
nuestras propias herramientas reutilizables. Podemos crear fácilmente nuestras propias
funciones y programas especiales (conocidos como M-archivos) en código Matlab, los
podemos agrupar en Toolbox (también llamadas librerías): colección especializada de
Marchivos para trabajar en clases particulares de problemas.
Matlab, a parte del cálculo matricial y álgebra lineal, también puede manejar polinomios,
funciones, ecuaciones diferenciales ordinarias, gráficos
1. MARCO TEORICO
DISTRIBUCION DE WEIBULL
La función de distribución de Weibull es un modelo estadístico que representa la probabilidad
de fallo después de un tiempo t (R (t)) en función del tiempo transcurrido o de una variable
análoga. O dicho de otra manera, R (t) es la probabilidad de que los componentes de un
conjunto sobrevivan hasta el momento t.
Esta función de probabilidad de fallo o función de fiabilidad R (t), viene dada por:
Donde los parámetros que definen la función son:
t : Tiempo entre fallas.
β : Es el parámetro de forma (Adimensional).
η : Es el parámetro de escala o tiempo característico (En unidades de tiempo).
γ : Es el parámetro de localización (En unidades de tiempo).
La función distribución acumulativa F (t) es el complemento de la función confiabilidad, la
función distribución acumulativa se puede interpretar como la Probabilidad de Falla y se
define de la siguiente manera:
2. DEDUCCIÓN DE LA ECUACIÓN LINEAL DE REGRESIÓN
Debido a que se desconoce los valores β, η e γ, se debe linealizar las curvas, es decir usar el
método de regresión lineal, este método permitirá obtener un polinomio que linealizará la
distribución de Weibull y permitirá estimar los parámetros β, η e γ, y tenemos los siguientes
pasos:
PASO I. Tenemos la Función acumulativa de Weibull:
PAS O II. Haciendo unos arreglos algebraicos tenemos:
PAS O III. Aplicando logaritmos naturales a ambos miembros:
PAS O IV. Nuevamente Aplicando logaritmos naturales a ambos miembros:
PAS O V. Luego queda al final la ecuación linealizada:
Por lo tanto la expresión representa una ecuación lineal de la forma:
Donde comparando tenemos:
3. CALCULO DE LOS PARAMETROS DE WEIBULL
Se tiene la siguiente tabla de recopilación de Tiempos entre Fallas obtenida por personal del
Camión, equipo HT044, modelo KOM 930E-4SE.
4. METODO DEL RANGO MEDIANO DE BERNARD PARA EL CALCULO DE F (t i)
A continuación, se explica un método (no el único, aunque sí uno de los más utilizados) que
permite obtener estimaciones puntuales de la función de distribución empírica, F (t i), a
partir de las observaciones, ti. Considerando que t1< t2<…tn son los “n” tiempos de fallo
observados. Para i = 1,2,…,n se puede obtener la segunda componente del punto (t,F(t i))
tomando para tiempos de fallo mayor a 20:
Para el caso del Camión, equipo HT044, modelo KOM 930E-4SE se tiene los siguientes valores:
n = 54
i = {1, 2,3,…, 38, 39, 40, 41, 42, 43,…..52, 53, 54 }
De lo expuesto se tiene definido las coordenadas con un valor de (γ = 0):
Abscisas(X) = Ln(t)
Ordenadas (Y) = Ln(Ln(1/(1-F (t))))
Una vez calculado los valores de las constantes “a” y “b” podemos calcular los siguientes
valores:
β=a
η = exp(-b/ β)
“Existen diferentes métodos para calcular los parámetros de WEIBULL, en este Manual se
detallará paso a paso el desarrollo del script en el SOFTWARE MATLAB además de la
comparación con los métodos: GRAFICO.
METODO DE MATLAB
Paso numero 1
Colocamos en primer lugar las funciones clc y clear para que cada vez que se ejecute el edit
en primer lugar borre los resultados que podrían haber en el Command Windows (clc ) y
también que limpie el Workspace ( clear )
clc, clear;
Paso numero 2
Para colocar un comentario en el edit de MATLAB es necesario antecederle el símbolo
“%” por ejemplo:
% Cálculo de los Parámetros de Weibull y la confiabilidad de una SISTEMA
DE REFRIGERACION DE UN MOTOR DE CAMION
% t = tiempo entre fallas (TBF)
Paso numero 3
Creamos una matriz fila “t” para introducir nuestros datos:
t=[0.86,1.14,1.35,1.45,2.96,4.72,6.06,6.79,7.43,7.63,7.74,8.53,8.85,11.75
,12.22,13.35,13.90,15.44,18.57,26.95,35.88,37.63,37.68,42.28,43.01,48.33,
49.21,54.49,58.88,60.96,72.40,79.57,87.07,92.04,94.71,95.08,97.28,103.51,
107.71,117.85,121.37,133.16,135.47,137.26,138.56,154.45,183.23,188.20,218
.10,221.35,277.37,284.56,353.55,417.48]*60;
t=sort(t);
Paso numero 4
Usamos la función s ort para ordenar nuestros datos de menor a mayor, es decir creamos
otra matriz con los datos iniciales pero ahora ordenados de forma ascendente.
t = sort (t);
Paso numero 5
Creamos una matriz para los valores de la fórmula de Bernard, para ello usamos la función
length (t). Length, hace el conteo del número de elementos de una matriz.
Entonces podemos calcularla siguiente matriz [1: leng th(t)] ; esta expresión nos indica que
se tomaran datos desde 1 hasta el valor resultante del length(t) ;
posteriormente calculamos el ajuste mediante el método de Bernard.
%FORMULA DE AJUSTE MEDIANTE EL MÉTODO DE RANGOS
MEDIANOS DE BENARD
F = ([1: length (t)]-0.3)/(length(t)+0.4);
Paso numero 6
Definimos las coordenadas mediante la función logaritmo natural ( log ):
% Definimos las coordenadas de "x" e "y"
x=log (t);
y=log (log (1. / (1-F)));
Paso numero 7
En el paso siguiente usaremos la función polyfi t para generar una recta.
Polyfit (x, y, n); esta función usa los valores de “x” e “y” para generar un polinomio de grado
“n”. Pero como nosotros buscamos una recta haremos que “n=1”, es así que se obtendrán dos
valores que vendrían a ser los coeficientes de la recta Y=aX+b es decir obtenemos una matriz
fila como sigue [a b] (linealizando la curva) % Linealización de la curva
Pol=polyfit (x, y, 1); %Se obtiene los valores pol= [pol(1) pol(2)]
coeficientes de la recta(y=pol(1)x + pol(2))
Para determinar los valores de beta y eta se quieren de los valores que nos da la función
polyfit(x, y, 1) anteriormente hallados como se mencionó [a b] para ello “a” será pol (1) y
“b” será pol (2).
%Se hallan los valores de "beta" y "eta"
beta = Pol (1);
eta = exp ((Pol(2)/(-Pol(1))));
Paso numero 8
Luego de generar la recta con ayuda de la función polyfit (x , y, 1) evaluaremos los valores
de x en esta recta, para ello nos ayudamos de la función polyv al (P ol, x) y obtendremos los
valores de la ordenada para esta recta.
%Se evalúa los valores de X en "pol" para obtener la ordenada de la recta
yy = polyval(Pol,x);
Paso numero 9
Para poder colocar algún texto o datos en el Command Windows usaremos la función disp.
disp(' '); coloca un espacio en blanco
d i s p ( 'TIEMPO ENTRE FA LLA S D EL C A RG A DOR FRONTA L C A T MOD
ELO
994F (5029)); el texto escrito entre apostrofes aparecerá en el Command Windows.
disp([(1:length(t))',t']); si se colocan matrices también aparecerán en el Command
Windows.
disp (' ');
disp('TIEMPO ENTRE FALLAS DEL CARGADOR FRONTAL CAT
MODELO 994F (5029)');
disp(' n t');
disp([(1:length(t))',t'])
disp('1°.-Se obtiene la ecuación de la recta: ');
disp(' ');
Paso numero 10
Para colocar textos pero haciendo que aparezcan resultados ya calculados podemos usar la
función fprintf.
fprintf ('Y =\t'); en este caso el texto que se verá en el Comman Windows es “Y =” y la parte
“ \t” hace que lo escrito debajo de esta función lo lleve en la misma línea.
También se puede escribir “ \n” en este caso estará separado en la siguiente fila de texto.
fprintf('%0.4f\t',Pol(1)); es este caso la función fprintf nos permite hacer que aparezca un
valor ya calculado, para ello se debe tener en cuenta lo que e stá dentro de los paréntesis
“'%0.4f\t',Pol(1)”
%0.4f; el “0” indica la separación, el “4” en número de decimales y “f” nos permite llamar
un dato que queremos que aparezca, para nuestro caso sería Pol(1).
fprintf('Y =\t');
fprintf('%0.4f\t',Pol(1));
fprintf('X\t');
fprintf('%0.4f\t',Pol(2));
disp('(Para Gamma = 0)')
fprintf('Beta (B) =\t');
fprintf('%0.4f\n',beta);
fprintf('Eta (n) =\t');
fprintf('%0.4f\t',eta);
FIGURA (1) GRAFICA DE WEIBULL
%x,yy Representan los puntos de la recta de ajuste
plot(x,yy,'b','LineWidth',2)
hold on
%x,y Representan los puntos que se van a ajustar por la recta
plot(x,y,'bo','LineWidth',1,'MarkerEdgeColor','k','MarkerFaceColor','r','Marker
Size',2)
title; una vez generado la figura podemos ponerle el titulo para ello esta función nos
permite hacer eso y otras cosas más que se detallaran a continuación.
([' A jus te de Weibull con ',' \g amma',' = 0 ']; la primera parte viene a ser el texto del título
“'Ajuste de Weibull con '” la segunda “' \gamma'” nos permite introducir el símbolo y la
última parte “' = 0 '” que también forma parte del texto que estará en el título.
'Color','w'; nos permite colocarle color al texto.
'FontSize',12; tamaño del fondo.
'HorizontalAlignment','center'; alinear y centrar texto.
'B ack g roundColor',[. 55 .100 .0]; no permite elegir el color de fondo, para este caso
podemos colocar el color mediante [.55 .100 . 0] que viene a ser como una paleta
de combinación de los colores primarios.
'EdgeColor','b'; color del marco.
'LineWidth',1; ancho de la línea.
min(yy); escoge el mínimo valor de la matriz seleccionada.
median(x); escoge el valor medio de la matriz seleccionada.
s et(g ca,' yG rid', 'on', 'y C olor', 'k ', 'L ineS tyleOrder', '-'); permite colocar líneas
auxiliares para mejorar la visualización de las gráficas.
Xlabel; permite colocar texto en eje x de la gráfica.
Ylabel; permite colocar texto en eje y de la gráfica.
Text; permite colocar texto adicional o leyenda en el gráfico.
num2str(); permite colocar datos en el texto.
title(['Ajuste de Weibull con ','\gamma',' = 0
'],'Color','w','FontSize',12,'HorizontalAlignment','center','BackgroundColor',[.5
5 .100 .0],'Margin',5,'EdgeColor','b','LineWidth',1)%Crear titulo
posy = min(yy);
posx= median(x);
set(gca,'yGrid','on','yColor','k','LineStyleOrder', '-')
set(gca,'xGrid','on','xColor','k','LineStyleOrder', '-')
xlabel('ln t','Color','k') %crea texto en el eje x de la figura 1
ylabel('ln(ln(1/(1-F(i))))','Color', 'k')%crea texto en el eje y de la figura 1
disp(' ');
disp(' ');
disp('2°.-Evaluamos el ajuste de los datos mediante(r^2)El cuadrado')
disp('del coeficiente de correlación de momento del producto Pearson.' )
R2=(sum((x-(sum(x)/length(x))).*(y-(sum(y)/length(y))))/((sum((x-
(sum(x)/length(x))).^2)*sum((y-(sum(y)/length(y))).^2))^.5))^2;
text(posx,posy,['y = ',num2str(Pol(1)),'x ',num2str(Pol(2)),' \wedge ','r^{2} =
',num2str(R2)],'FontSize',10,'HorizontalAlignment','center','BackgroundColor'
,[.8 .8 .8],'Margin',6,'EdgeColor','k','LineWidth',2)
disp(' ');
fprintf('R2 =\t')
fprintf('%0.4f\n',R2);
Para esta parte se pretende realizar cálculos para gamma mayores que cero, para
ello usamos la función for .
for ; nos permite hacer iteraciones
%Usamos FOR para evaluar el valor de r^2 en j=1:n
n=(t(1,1)-1); %donde n < t(min)
for j=1:n; % "j" representa los posibles valores de gamma
yyy=1:n;
t2=[t-j]; % los nuevos tiempos entre fallas afectado por gamma
x=log(t2);
r2(1,j)=(sum((x-(sum(x)/length(x))).*(y-(sum(y)/length(y))))/((sum((x-
(sum(x)/length(x))).^2)*sum((y-(sum(y)/length(y))).^2))^.5))^2;% R2 para
cada "j"
end
r=r2;
disp(' ');
disp('3°.-Hallando el valor máximo de r^2 que ajuste a la recta' )
disp('con la observación de que se han tomado Gammas positivos y
menores que el valor minimo de t' )
%La siguiente linea calcula el parametro r^2 maximo y Gamma
disp(' ');
%Escribimos "Gama" porque el´programa MATLAB ya posee una funcion
interna Gamma
[R2,Gama] = max(r);
disp('El valor de r^2 que maximiza el ajuste es: ' )
fprintf('R2 =\t')
fprintf('%0.4f\n',R2);
FIGURA (2) GRAFICA DE DISPERCION DE DATOS
plot(r,yyy)
set(gca,'YGrid','on','YColor','k','LineStyleOrder', '-')
title('Sensibilidad y Error color','w','FontSize',12,'HorizontalAlignment','center','Backgroun
dColor',[.55 .100 .0],'Margin',5,'EdgeColor','k','LineWidth',1)
set(gca,'YGrid','on','YColor','k','LineStyleOrder', '-')
posx = min(r);
posy= median(yyy);
text(posx,posy,['Max. aproximación en: ',' (r^{2}) max = ',num2str(R2),'
\wedge ',' \gamma =
',num2str(Gama)],'FontSize',10,'HorizontalAlignment','left','BackgroundColor'
,[.8 .8 .8],'Margin',6,'EdgeColor','k','LineWidth',1)
xlim([]); permite definir el dominio de “x”.
ylim([]) ; permite definir el rango de “y”.
ylim([0 max(Gama)+1]);
xlabel(['Coeficiente de determinacion',' (r^{2})'],'Color','k')
ylabel(['Gamma','(\gamma)',' en [Horas]'],'Color','k')
hold on
plot(R2,Gama,'+','LineWidth',8,'MarkerEdgeColor',[.55 .100 .0],'MarkerSize',2)
disp(' ');
disp('4°.-Ajustamos nuevamente los datos de falla y obtenemos ');
disp('la nueva recta corregida: ');
disp(' ');
x=log(t-Gama);
Pol=polyfit(x,y,1);
yyyy = polyval(Pol,x);
fprintf(' n t t (Gamma=\t');
fprintf('%0.2f\t',Gama);
disp(')');
disp([(1:length(t))',t',(t-Gama)'])
FIGURA (3) GRAFICA DE REPRESENTACION DE AJUSTE DE LA RECTA
title(); se explicó anteriormente pero acotaremos una función más.
int2str(Gama); permite visualizar el resultado de un dato, para nuestro caso
queremos ver el resultado de “Gama”.
title(['Ajuste de Weibull con ','\gamma',' =
',int2str(Gama)],'Color','w','FontSize',12,'HorizontalAlignment','center','Backg
roundColor',[.55 .100 .0],'Margin',5,'EdgeColor','b','LineWidth',1)
posy= min(yyyy);
posx = max(x);
text(posx,posy,['y = ',num2str(Pol(1)),'x ',num2str(Pol(2)),' \wedge ','r^{2} =
',num2str(R2)],'FontSize',10,'HorizontalAlignment','right','BackgroundColor',[.
8 .8 .8],'Margin',6,'EdgeColor','k','LineWidth',2)
set(gca,'YGrid','on','YColor','k','LineStyleOrder', '-')
xlabel('ln (t-\gamma)','Color','k')
ylabel('ln(ln(1/(1-F(i))))','Color','k')
fprintf('y =\t');
fprintf('%0.4f\t',Pol(1));
fprintf('x\t');
fprintf('%0.4f\t',Pol(2));
fprintf('(Ecuación de la recta de Ajuste con Gamma=' );
fprintf('%0.2f\t',Gama);
disp(')');
Beta=Pol(1);
disp(' ');
disp('5°.-Obtenemos los parámetros corregidos de Weibull:')
disp(' ');
fprintf('Beta (B) =\t');
fprintf('%0.4f\n',Beta);
Eta=exp(Pol(2)/(-Pol(1)));
fprintf('Eta (n) =\t');
fprintf('%0.4f\n',Eta);
fprintf('Gamma =\t');
fprintf('%0.2f\n',Gama);
disp(' ');
disp('6°.-Evaluamos el parámetro estadístico de (DAlfa>Dmax)')
disp('Buscar en la tabla K-S el DAlfa para una confianza de Alfa =0.05' )
Ft=1-exp(-((t-Gama)./Eta).^Beta); % infiabilidad
Dni=abs(Ft-F);%La máxima discrepancia obtenida del proceso de comparación entre la
distribución propuesta y la distribución verdadera.
if ; nos evalúa una condición a la cual podemos decir que puede ser verdadera o
falsa.
else; ejecutado si la primera condicional es falsa.
Dmax=max(Dni);
disp(' ');
DAlfa=0.2417;% La máxima discrepancia que se puede aceptar para un
determinado nivel de confianza (saca de la tabla de Kolmogorov Smirnov)
disp(' ');
if (DAlfa>Dmax);
fprintf('El modelo es WEIBULL se ACEPTA (DAlfa>Dmax)\n');
else
fprintf('El modelo no es WEIBULL se RECHAZA (DAlfa<Dmax)\n' );
end
disp(' ');
fprintf('DAlfa =\t');
fprintf('%0.4f\n',DAlfa);
fprintf('Dmax =\t');
fprintf('%0.4f\n',Dmax);
Rt=exp(-((t-Gama)./Eta).^(Beta));
FIGURA (4) GRAFICA DE CONFIABILIDAD
plot(t,Rt,'bo','LineWidth',1,'MarkerEdgeColor','k','MarkerFaceColor','r','MarkerSize',2)
title('Confiabilidad','Color','w','FontSize',12,'HorizontalAlignment','center','BackgroundCol
or',[.55 .100 .0],'Margin',5,'EdgeColor','b','LineWidth',1)
set(gca,'YGrid','on','YColor','k','LineStyleOrder', '-')
posx = 2*median(t);
posy= median(Rt);
text(posx,posy,['R(t) = e ','^{-((\gamma-t)/\eta)}',' ^
','^{\beta}'],'FontSize',10,'HorizontalAlignment','left','BackgroundColor',[.8
.8.8],'Margin',6,'EdgeColor','k','LineWidth',1)
xlabel('Tiempo (t) en [Horas]','Color','k')
ylabel('Confiabilidad R(t)','Color','k')
disp(' ');
disp('Buscar en la tabla para una confianza de Alfa =0.05')
xx=t;%valores ordenados de forma creciente
rot90(rot90(t)); esta función nos permite rotar la matriz en 180, se usó para tener
datos ordenados de forma descendente, es decir de la mayor a menor. Pero antes
se tuvo que ordenar con ayuda de la función s ort(t).
xxinv=rot90(rot90(t));% valores ordenados de forma decreciente
xxdif=(xxinv-xx);% diferencia entre los valores anteriores
ai=[0.4254 0.2944 0.2487 0.2148 0.1870 0.1630 0.1415 0.1219 0.1036
0.0862 0.0697 0.0537 0.0381 0.0227 0.0076 ];
%tomados de tabla
B=ai.*xxdif;
m=median(xx);
D=(xx-m).^2;
%valores requeridos para hallar W
b=sum(B);
d=sum(D);
Wc=b.^2/d;
Wt=0.927;% Distribución del estadístico de Shapiro-Wilk (wt) para alfa=0.05
disp(' ');
if (Wt>Wc);
fprintf('El modelo es WEIBULL se ACEPTA (Wt>Wc)\n' );
else
fprintf('El modelo no es WEIBULL se RECHAZA (Wt<Wc)\n' );
end
disp(' ');
fprintf('Wt =\t');
fprintf('%0.4f\n',Wt);
fprintf('Wc =\t');
fprintf('%0.4f\n',Wc);
disp(' ');
disp('8°.-Evaluamos los parámetro estadístico mediante Test de Anderson
Darling')
disp('comprobar que ADe>ADx para alfa=0.05')
Finv=rot90(rot90(F));
n=length(t);
j=1:n;
coef=(2*j-1)/n;
A=log(F);
B=log(1-Finv);
Res=coef.*(A+B);
S=sum(Res);
ADx =-n-S;
ADe=0.757;% de tabla para un ajuste del 95%
disp(' ');
if (ADe>ADx);
fprintf('El modelo es WEIBULL se ACEPTA (ADe>ADx)\n');
else
fprintf('El modelo no es WEIBULL se RECHAZA (ADe<ADx)\n');
end
disp(' ');
fprintf('ADe =\t');
fprintf('%0.4f\n',ADe);
fprintf('ADx =\t');
fprintf('%0.4f\n',ADx);
disp(' ');
disp(' ');
disp('9°.-Encontramos el Tiempo Promedio Entre Fallas "MTBF" en Horas')
MTBF=Eta*gamma(1+(1/Beta))+Gama
disp(' ');
disp('10°.-Resultados')
TF=(((t-Gama)./Eta).^(Beta-1)).*(Beta/Eta);
disp(' ');
disp(' n t R(t)% F(t)% ');
disp([(1:length(t))',t',Rt'*100,Ft'*100]);
FIGURA (5) GRAFICA DE INFIABILIDAD
plot(t,Ft,'bo','LineWidth',1,'MarkerEdgeColor','k','MarkerFaceColor','r','Marker
Size',2)
title('Infiabilidad Distribución de
Weibull','Color','w','FontSize',12,'HorizontalAlignment','center','BackgroundC
olor',[.55 .100 .0],'Margin',5,'EdgeColor','b','LineWidth',1)
posx = 2*median(t);
posy= median(Ft);
text(posx,posy,['F(t) = 1 - e ','^{-((\gamma-t)/\eta)}',' ^
','^{\beta}'],'FontSize',10,'HorizontalAlignment','left','BackgroundColor',[.8 .8
.8],'Margin',6,'EdgeColor','k','LineWidth',1)
xlim([0 max(t)+30]);
set(gca,'YGrid','on','YColor','k','LineStyleOrder', '-')
xlabel('Tiempo (t) en [Horas]','Color','k')
ylabel('Infiabilidad F(t)','Color','k')
INTRODUCIENDO EN MATLAB:
Se tiene:
% INGENIERIA DE MANTENIMIENTO %%%
%Cálculo de la Fiabilidad de un camion equipo HT044, modelo KOM 930E-
4SE.
%t = tiempo entre fallas (TBF)
clc,clear
t =[0.86 1.14 1.35 1.45 2.96 4.72 6.06 6.79 7.43 7.63 7.74 8.53 8.85
11.75 12.22 13.35 13.9 15.44 18.57 26.95 35.88 37.63 37.68 42.28 43.01
48.33 49.21 54.49 58.88 60.96 72.4 79.57 87.07 92.04 94.71 95.08 97.28
103.51 107.71 117.85 121.37 133.16 135.47 137.26 138.56 154.45 183.23
188.2 218.1 221.35 277.37 284.56 353.55 417.48]*60;
t=sort(t); % t=(tiempo entre fallas de menor a mayor)
n= (1:length(t)); %
% "FORMULA DE AJUSTE MEDIANTE EL MÉTODO DE RANGOS MEDIANOS DE BENARD"
F=([1:length(t)]-0.3)/(length(t)+0.4); % frecuencia acumulada de fallas x
100%
X=log(t);
Y=log(log(1./(1-F)));%Ajustando los datos de falla linealizando
disp('DETERMINACIÓN DE PARAMETROS DE WEIBULL, LA CONFIABILIDAD Y EL MTBF
')
disp(' ')
b=[n ; t ; F*100 ; X ; Y ]';
disp(' n t F*100 X Y ')
disp(b)
P=polyfit(X,Y,1); % muestra del grafico LINEALIZADO y=ax+b ,
% a=beta
beta=P(1); % se calcula beta de la recta
eta=exp((P(2)/(-P(1)))); % % se calcula eta de la recta n
Y2=polyval(P,X); % reemplaza cada valor de X=log(t) en cada funcion p
disp(' ');
disp('A.-Ajustamos los datos de falla al polinomio siguiente: ');
disp(' ');
fprintf('Y =\t');
fprintf('%0.4f\t',P(1));
fprintf('X \t');
fprintf('%0.4f\t',P(2));
disp('(Ecuación de la recta de Ajuste con Gamma = 0)')
disp('Beta (B) =');
fprintf('%0.4f\n',beta);
disp('Eta (n) =');
fprintf('%0.4f\n',eta);
figure(1);
%grafica del ajuste de weibull
%X,Y2 Representan los puntos de la recta de ajuste
plot(X,Y2,'-b','LineWidth',3)
hold on
%X,Y Representan los puntos que se van a ajustar por la recta
plot(X,Y,'o','LineWidth',1,'MarkerEdgeColor','k','MarkerFaceColor','g','M
arkerSize',4)
title(['Ajuste de Weibull con ','\gamma',' = 0
'],'Color','k','FontSize',12,'HorizontalAlignment','center','BackgroundCo
lor','g','Margin',6,'EdgeColor','k','LineWidth',2)
[d11,c11] = min(Y2);
[d12,c12] = max(X);
grid on
xlabel('ln t','Color','b','FontSize',12)
ylabel('ln(ln(1/(1-F(i))))','Color','b','FontSize',12)
disp(' ');
%B
disp('B.-Evaluamos el ajuste de los datos')
disp(' ')
R=(sum((X-(sum(X)/length(X))).*(Y-(sum(Y)/length(Y))))/((sum((X-
(sum(X)/length(X))).^2)*sum((Y-(sum(Y)/length(Y))).^2))^.5));%
coeficiente de correlacion de pearson
Rcuadrado=(sum((X-(sum(X)/length(X))).*(Y-(sum(Y)/length(Y))))/((sum((X-
(sum(X)/length(X))).^2)*sum((Y-(sum(Y)/length(Y))).^2))^.5))^2; %
coeficiente de determinacion
text(d12-2,d11,['y = ',num2str(P(1)),'x ',num2str(P(2)),' \wedge
','r^{2} =
',num2str(Rcuadrado)],'FontSize',10,'HorizontalAlignment','center','Backg
roundColor',[.8 .8 .50],'Margin',6,'EdgeColor','k','LineWidth',1)%texto
en el grafico
disp(' ');
disp('R =')
fprintf('%0.4f\n',R);
disp('Rcuadrado =')
fprintf('%0.4f\n',Rcuadrado);
%C es el máximo valor de r^2
%La sentencia de FOR evalua el r^2 en j=1:n.....NUMERO DE ITERACIONES
n=(t(1,1)-1); % n=(t(1)-1);
for j=1:n;
B=1:n;
%F=[1:length(t)]/(length(t)+1);
X=log(t-j);
e(1,j)=(sum((X-(sum(X)/length(X))).*(Y-(sum(Y)/length(Y))))/((sum((X-
(sum(X)/length(X))).^2)*sum((Y-(sum(Y)/length(Y))).^2))^.5))^2;
end
r2=e;
disp(' ');
disp('C.-Encontramos el valor máximo de r^2 que ajusta la recta')
disp('con la observación de que se ha incluido el parámetro Gamma > 0')
disp('a los datos de falla')
%La siguiente línea calcula el parámetro r^2 máximo y Gamma
disp(' ');
%Escribimos la variable "Gama" de esta forma para salvar un conflicto que
se presentaria
%con la función interna Gamma de MATLAB
[Rcuadrado,Gama] = max(r2);r11=Rcuadrado;g11=Gama;
disp('El valor de r^2 que maximiza el ajuste es el siguiente: ')
disp('Rcuadrado =')
fprintf('%0.4f\n',Rcuadrado);
figure(2);
plot(r2,B,'-o')
title('Sensibilidad y Error
cuadrático','Color','k','FontSize',12,'HorizontalAlignment','center','Bac
kgroundColor',[.20 .80 .20],'Margin',6,'EdgeColor','k','LineWidth',1)
d11 = min(r2);
d12 = median(B);
text(d11,d12,['Max. aproximación en: ',' (r^{2}) max =
',num2str(Rcuadrado),' \wedge ',' \gamma =
',num2str(Gama)],'FontSize',10,'HorizontalAlignment','left','BackgroundCo
lor',[.8 .8 .50],'Margin',6,'EdgeColor','k','LineWidth',1)
grid on
xlabel(['Coeficiente de determinacion ','
(r^{2})'],'Color','b','FontSize',12)
ylabel(['Gamma','(\gamma)',' en [minutos]'],'Color','b','FontSize',12)
hold on
plot(r11,g11,'o','LineWidth',2,'MarkerEdgeColor','r','MarkerFaceColor','r
','MarkerSize',2)
disp(' ');
% D.- AJUSTAMOS NUEVAMENTE LOS DATOS DE FALLA CON POLYFIT
disp('D.-Ajustamos nuevamente los datos de falla y obtenemos ');
disp('la nueva recta corregida: ');
disp(' ');
X=log(t-Gama);
P=polyfit(X,Y,1); % muestra del grafico LINEALIZADO y=ax+b , a i b
fig3 = polyval(P,X);
figure(3);
%X,Y2 Representan los puntos de la recta de ajuste
plot(X,fig3,'-
b','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','r','MarkerSize'
,2)
hold on
%X,Y Representan los puntos que se van a ajustar por la recta
plot(X,Y,'o','LineWidth',1,'MarkerEdgeColor','k','MarkerFaceColor','r','M
arkerSize',4)
title(['Ajuste de Weibull con ','\gamma',' =
',int2str(Gama)],'Color','k','FontSize',12,'HorizontalAlignment','center'
,'BackgroundColor',[.20 .80
.20],'Margin',6,'EdgeColor','k','LineWidth',1)
[d11,c11] = min(Y2);
[d12,c12] = max(X);
text(d12,d11,['y = ',num2str(P(1)),'x ',num2str(P(2)),' \wedge
','r^{2} =
',num2str(Rcuadrado)],'FontSize',10,'HorizontalAlignment','right','Backgr
oundColor',[.8 .8 .50],'Margin',6,'EdgeColor','k','LineWidth',1)
grid on
xlabel('ln (t-\gamma)','Color','b','FontSize',12)
ylabel('ln(ln(1/(1-F(i))))','Color','b','FontSize',12)
fprintf('Y =\t');
fprintf('%0.4f\t',P(1));
fprintf('X=\t');
fprintf('%0.4f\t',P(2));
disp('(Ecuación de la recta de Ajuste con Gamma > 0)')
Beta=P(1);
disp(' ');
% PARAMETROS CORREGIDOS DE WEIBULL
disp('E.-Obtenemos los parametros corregidos de Weibull:')
disp(' ');
disp('Beta (B) =');
fprintf('%0.4f\n',Beta);
Eta=exp(P(2)/(-P(1)));
disp('Eta (n) =');
fprintf('%0.4f\n',Eta);
disp('Gamma =');
fprintf('%0.4f\n',Gama);
disp(' ');
% Evaluamos el parametro estádistico de Kolmogorov Smirnov
disp('F.-Evaluamos el parametro estádistico de Kolmogorov Smirnov
(D_Alfa>Valor_Dmax)')
disp('Buscar en la tabla K-S el D_Alfa para una confianza de Alfa =0.05')
Ft=1-exp(-((t-Gama)./Eta).^Beta);
Dni=((Ft-F).^2).^0.5;
[Valor_Dmax,D_Posicion] = max(Dni);
disp(' ');
q=0.18841;
disp(' ');
if (q>Valor_Dmax);
fprintf('El modelo es WEIBULL se ACEPTA (D_Alfa>Valor_Dmax)\n');
else
fprintf('El modelo no es WEIBULL se RECHAZA (D_Alfa<Valor_Dmax)\n');
end
disp(' ');
disp('Valor_Dmax =');
fprintf('%0.4f\n',Valor_Dmax);
disp('D_Alfa =');
fprintf('%0.4f\n',q);
%Confiabilidadad
Rt=exp(-((t-Gama)./Eta).^(Beta));
figure(4)
plot(t,Rt,'-
bo','LineWidth',1,'MarkerEdgeColor','k','MarkerFaceColor','r','MarkerSize
',2)
title(' Confiabilidad
','Color','k','FontSize',12,'HorizontalAlignment','center','BackgroundCol
or',[.20 .80 .20],'Margin',6,'EdgeColor','k','LineWidth',1)
[d11,c11] = min(t);
d12= median(Rt);
text(d11,d12,['R(t) = e ','^{-((\gamma-t)/\eta)}',' ^
','^{\beta}'],'FontSize',10,'HorizontalAlignment','left','BackgroundColor
',[.8 .8 .50],'Margin',6,'EdgeColor','k','LineWidth',1)
grid on
xlabel('Tiempo (t) en [minutos]','Color','b','FontSize',12)
ylabel('Confiabilidad R(t)','Color','b','FontSize',12)
disp(' ');
% H encontrar en tiempo de falla
disp('H.-Encontramos el Tiempo Promedio Entre Fallas "MTBF" en Horas')
MTBF=Eta*gamma(1+(1/Beta))+Gama
disp(' ');
disp('I.-Encontramos la Tasa de Fallas');
%Tasa de Fallos de Acuerdo a la Distribución de Weibull
TF=((((t-Gama)./Eta).^(Beta-1)).*(Beta/Eta))';
figure(5);
plot(t,TF,'-
bo','LineWidth',1,'MarkerEdgeColor','k','MarkerFaceColor','r','MarkerSize
',2)
title('Tasa de Fallas de Acuerdo a Distribución de
Weibull','Color','k','FontSize',12,'HorizontalAlignment','center','Backgr
oundColor',[.20 .80 .20],'Margin',6,'EdgeColor','k','LineWidth',1)
%Modificación de valor inicial y final de eje X
grid on
xlabel('Tiempo (t) en [minutos]','Color','b','FontSize',12)
ylabel('Tasa de Fallas: \lambda(t)
[Fallas/Hora]','Color','b','FontSize',12)
% Infiabilidad
It=1-exp(-((t-Gama)./Eta).^(Beta));
figure(6)
plot(t,It,'-
bo','LineWidth',1,'MarkerEdgeColor','k','MarkerFaceColor','r','MarkerSize
',2)
title(' Infiabilidad
','Color','k','FontSize',12,'HorizontalAlignment','center','BackgroundCol
or',[.20 .80 .20],'Margin',6,'EdgeColor','k','LineWidth',1)
grid on
[d11,c11] = min(t);
d12= median(Rt);
text(d11,d12,['I(t) = 1 - e ','^{-((\gamma-t)/\eta)}',' ^
','^{\beta}'],'FontSize',10,'HorizontalAlignment','left','BackgroundColor
',[.8 .8 .50],'Margin',6,'EdgeColor','k','LineWidth',1)
xlabel('Tiempo (t) en [minutos]','Color','b','FontSize',12)
ylabel('Infiabilidad I(t)','Color','b','FontSize',12)
fig: 1
Ajuste de Weibull con = 0
2
0
ln(ln(1/(1-F(i))))
-1
-2
-3
y = 0.76381x -6.4358 r2 = 0.97451
-4
-5
3 4 5 6 7 8 9 10 11
ln t
Fig: 02
Sensibilidad y Error cuadrático
50
45
40
Gamma() en [minutos]
35
30
25 Max. aproximación en: (r2) max = 0.98419 = 42
20
15
10
0
0.974 0.976 0.978 0.98 0.982 0.984 0.986 0.988
2
Coeficiente de determinacion (r )
Fig. 03
Ajuste de Weibull con = 42
2
0
ln(ln(1/(1-F(i))))
-1
-2
-3
y = 0.67876x -5.7093 r2 = 0.98419
-4
-5
2 3 4 5 6 7 8 9 10 11
ln (t-)
Fig. 04
Confiabilidad
1
0.9
0.8
0.7
Confiabilidad R(t)
0.6
0.5
R(t) = e -((-t)/)
0.4
0.3
0.2
0.1
0
0 0.5 1 1.5 2 2.5 3
Tiempo (t) en [minutos] 4
x 10
Fig: 05
x 10
-3
Tasa de Fallas de Acuerdo a Distribución de Weibull
1.2
1
Tasa de Fallas: (t) [Fallas/Hora]
0.8
0.6
0.4
0.2
0
0 0.5 1 1.5 2 2.5 3
Tiempo (t) en [minutos] 4
x 10
Fig 06
Infiabilidad
1
0.9
0.8
0.7
Infiabilidad I(t)
0.6
0.5
I(t) = 1 - e -((-t)/)
0.4
0.3
0.2
0.1
0
0 0.5 1 1.5 2 2.5 3
Tiempo (t) en [minutos] 4
x 10
