Función de distribución de Weibull

Media aritmética y desviación estándar

Si tenemos un conjunto de n números [x1 x2 x3....xn], se define la media aritmética <x>

<x>= 1 n i=1 n x i

Si queremos conocer la desviación del conjunto de datos respecto de la media, definimos la desviación estándar σ

σ 2 = 1 n1 i=1 n ( x i <x>) 2

MATLAB calcula la media y desviación estándar de un conjunto de datos llamando a las funciones mean y std.

>> x=[1.65 1.82 1.72 1.75 1.73 1.85 1.90 1.74 1.76 1.77];
>> mean(x)
ans=1.7690 
>> std(x)
ans=0.0713 

En ocasiones las velocidades del viento vienen dadas por números enteros, de modo que cada valor entero se mide varias veces durante el periodo de observación. Así, si el dato xi se ha medido wi veces. Tendremos dos vectores x de las medidas y w de las observaciones de dicha medida (pesos) de la misma dimensión m.

x=[x1, x2, x3,....xm]

w=[w1, w2, w3,....wm]

Las fórmulas de la media aritmética y desviación estándar se expresan:

<x>= 1 n i=1 n w i x i σ 2 = 1 n1 i=1 n w i ( x i <x>) 2

orden, i medida, x observaciones, w
1 6 19
2 7 54
3 8 42
total, n 115

Definimos dos funciones w_media y w_estandar que nos devuelvan la media aritmética y la desviación estándar cuando se nos proporcionan los datos de las medidas xi y su correspondiente número de observaciones wi.

Para calcular el valor medio definimos la función w_media

function res=w_media(x,w)
    res=sum(w.*x)/sum(w);
end

Para calcular la desviación estándar definimos la función w_estandar

function res=w_estandar(x,w)
    n=sum(w);    
    media=sum(w.*x)/n;    
    suma2=sum(w.*(x-media).^2);
    res=sqrt(suma2/(n-1));
end

En la ventana de comandos probamos las dos funciones que hemos definido

>> x=[6 7 8];
>> w=[19 54 42];
>> w_media(x,w)
ans =    7.2000
>> w_estandar(x,w)
ans =    0.7034

Frecuencias

Se define frecuencia f(xi) de una medida xi como el cociente del número wi de observaciones de dicha medida dividida por el número total n de medidas. La frecuencia acumulada F(xi) se define como la suma de las frecuencias desde k=1 a i

f( x i )= w i n i=1 m f( x i )=1 F( x i )= k=1 i f( x k )

La tabla anterior de medidas xi y número de observaciones de cada medida wi la convertimos en tabla de frecuencias f(xi) y frecuencias acumuladas F(xi).

orden, i medida, x observaciones, w frecuencia, f frecuencia acumulada, F
1 6 19 0.165 0.165
2 7 54 0.470 0.635
3 8 42 0.365 1.0
total, n 115 1.0  

Modelo estadístico

Es conveniente establecer un modelo de las frecuencias de las velocidades del viento (última figura de la página 'Análisis de los datos del viento') que venga descrito por una función matemática continua en vez de por una tabla de valores discretos.

La función f(x) representa la probabilidad de que la velocidad del viento esté en un intervalo entre x y x+dx. El área bajo f(x) es la unidad.

0 f(x)·dx=1

La probabilidad de que la velocidad del viento esté comprendida entre x0 y x1 viene dada por la integral

P( x 0 x x 1 )= x 0 x 1 f(x)·dx

Si el número N de medidas de la velocidad del viento es grande y se ha hecho una medida en el intervalo de tiempo Δt (10 minutos, una hora). El tiempo en el que el viento está soplando con una velocidad comprendida entre entre x0 y x1 es N·Δt·P(x0x< x1)

El valor medio de la velocidad <x> y la desviación estándar σ tienen ahora la siguiente definición:

<x>= 0 x·f(x)·dx σ 2 = 0 (x<x>) 2 ·f(x)·dx

La distribución de probabilidad acumulada F(x) se define

F(x)= 0 x f(x)·dx

Hay varias funciones f(x) que se pueden utilizar para describir la frecuencia de la distribución de velocidades del viento (última figura de la página 'Análisis de los datos del viento'). Las más utilizadas son las funciones de Weibull y Rayleigh.

Función de distribución de Weibull

La función de distribución Weibull depende de dos parámetros denominados c y k y la función de distribución de Rayleigh de un sólo parámetro. Esto hace que la primera sea más versátil y preferida que la segunda por lo que la estableceremos como modelo.

f(x)= k c ( x c ) k1 exp[ ( x c ) k ](k>0,x>0,c>1)

Para dibujar las gráficas de esta función con c=1 y variando el parámetro k, escribimos el script

c=1;
K=[1.2 1.6 2.0 2.4 2.8];
f=@(k,x) (k/c)*((x/c).^(k-1)).*exp(-(x/c).^k);
x=linspace(0,2.5,100);
hold on
for i=1:length(K)
    plot(x,f(K(i),x),'displayName',num2str(K(i)))
end
ylim([0 1.2])
xlabel('x')
ylabel('f(x)')
title('Función de Weibull')
legend('-DynamicLegend','location','NorthEast')
hold off

La velocidad para la cual la función de distribución de Weibull alcanza un máximo se obtiene derivando f(x) e igualando a cero.

k c k exp[ ( x c ) k ] x k2 { (k1) k c k x k }=0 x max =c ( k1 k ) 1/k

El valor medio <x>

<x>= 0 x·f(x)·dx= 0 x· k c ( x c ) k1 exp[ ( x c ) k ]·dx= c 0 y 1/k e y dy=cΓ( 1+ 1 k )

Donde se ha hecho el cambio de variable y=(x/c)k,

La función gamma Γ(x) se define

Γ(x)= 0 e t t x1 dt

Si n es un número entero Γ(n+1) =n! (factorial de n). Para otros valores x la función gamma de MATLAB nos devuelve el valor de la integral.

Comprobamos en la ventana de comandos que el área bajo la curva de cualquiera de las gráficas es la unidad, y que el valor medio <x>=c·Γ(1+1/k)

>> k=2.8;c=1;
>> f=@(x) (k/c)*((x/c).^(k-1)).*exp(-(x/c).^k);
>> quad(f,0,3) %area bajo la curva f(x)
ans =    1.0000
>> g=@(x) x.*f(x);
>> quad(g,0,3) %valor medio
ans =    0.8905
>> c*gamma(1+1/k)  %valor medio
ans =    0.8905

La desvación estándar σ

σ 2 = 0 (x<x>) 2 ·f(x)·dx= c 2 [ Γ( 1+ 2 k ) Γ 2 ( 1+ 1 k ) ]

Otros resultados interesantes son los siguientes:

Representamos f(x) y calculamos la probabilidad P(x) en el intervalo (0.75, 1.25), es decir, centrado en x=1 y de anchura 0.5

k=2.8;c=1;
f=@(x) (k/c)*((x/c).^(k-1)).*exp(-(x/c).^k);
x=linspace(0,3,100);
y=f(x);
hold on
plot(x,y,'r')
x0=0.75;x1=1.25;
xx=[x0 x0 x(x>x0 & x<x1) x1 x1];
yy=[0 f(x0) y(x>x0 & x<x1) f(x1) 0];
fill(xx,yy,'y');
res=quad(f,x0,x1)
prob=exp(-(x0/c)^k)-exp(-(x1/c)^k)
text(1.5, max(y)-0.1,num2str(res));
title('Probabilidad')
xlabel('x')
ylabel('f(x)')
hold off

Obtenemos el mismo resultado efectuando la integración numérica de la función de Weibull f(x) mediante la función quad, que a partir de la expresión de la probabilidad P(x0xx1) deducida anteriormente.

Finalmente, la probabilidad acumulada F(x) vale

F(x)= 0 x f(x)·dx=1exp[ ( x c ) k ]

Las funciones f(x) y F(x) tienen las siguientes características:

Representamos esta función con c=1 y variando el parámetro k, creando script similar al anterior, cambiando la definición de la función de f(x) a F(x).

c=1;
K=[1.2 1.6 2.0 2.4 2.8];
f=@(k,x) 1-exp(-(x/c).^k);
x=linspace(0,2.5,100);
hold on
for i=1:length(K)
    plot(x,f(K(i),x),'displayName',num2str(K(i)))
end
ylim([0 1.2])
xlabel('x')
ylabel('f(x)')
title('Probabilidad acumulada')
legend('-DynamicLegend','location','SouthEast')
hold off

Ajuste de los datos a la función de Weibull

Mediante el script viento_1, determinamos el número de horas (número de datos) en las que la velocidad del viento tenía un valor medio comprendido entre 0 y 1 m/s, entre 1 y 2, ... entre 22 y 23. Los datos guardados en el vector horas los utilizamos para crear la siguiente tabla.

orden, i velocidad, x observaciones, w frecuencia , f frecuencia acumulada, F
1 0.5 6 0.0081 0.0081
2 1.5 31 0.0417 0.0479
3 2.5 30 0.0403 0.0901
4 3.5 50 0.0672 0.1573
5 4.5 69 0.0927 0.2500
6 5.5 75 0.1008 0.3508
7 6.5 60 0.0806 0.4315
8 7.5 59 0.0793 0.5108
9 8.5 69 0.0927 0.6035
10 9.5 61 0.0820 0.6855
11 10.5 43 0.0578 0.7443
12 11.5 44 0.0591 0.8024
13 12.5 54 0.0726 0.8750
14 13.5 43 0.0578 0.9328
15 14.5 21 0.0282 0.9610
16 15.5 11 0.0148 0.9758
17 16.5 7 0.0094 0.9852
18 17.5 4 0.0054 0.9906
19 18.5 2 0.0027 0.9933
20 19.5 1 0.0013 0.9946
21 20.5 1 0.0013 0.9960
22 21.5 3 0.0040 1.0000

Para calcular la frecuencia acumulada se ha llamado a la función cumsum de MATLAB.

>> horas/sum(horas)  %frecuencia
>> cumsum(horas)/sum(horas)   %frecuencia acumulada

Ajuste no lineal

Para ajustar los datos de la columna frecuencia a la función de Weibull tendremos que determinar los valores de los parámetros c y k. Haremos un ajuste no lineal con la función nlinfit de MATLAB, para lo cual llamaremos a(1) al parámetro k y a(2) al parámetro c.

En color azul claro vemos el diagrama de frecuencias y en color rojo la función de Weibull que mejor ajusta.

clear,clc
velocidad=xlsread('WhiteDeer2013','Mar','F2:F745');
%interpolar si es necesario
if any(isnan(velocidad)) %si hay algún NaN
    x=1:length(velocidad);
    i=find(~isnan(velocidad));
    velocidad=interp1(x(i),velocidad(i),x);
end
%histograma
x=0.5:1:max(velocidad);
horas=hist(velocidad,x);

%convierte a frecuencias y ajusta a la función de Weibull
frec=horas/sum(horas);
f=@(a,x) (a(1)/a(2))*((x/a(2)).^(a(1)-1)).*exp(-(x/a(2)).^a(1));
a0=[2 8];  %valor inicial de los parámetros
af=nlinfit(x,frec,f,a0)
hold on 

%diagrama de frecuencias
bar(x,frec,'c');

%representa la curva de ajuste
x=linspace(0,max(velocidad),100);
y=f(af,x);
plot(x,y,'r')

title('Ajuste a la función Weibull')
xlabel('Velocidad')
ylabel('Frecuencia')
hold off

En la ventana de comandos obtenemos los valores de los parámetros k y c de ajuste

    af =    2.0486    9.4165

El dato de k=af(1)=2.0486 y c=af(2)=9.4165

Ajuste lineal

En la expresión de las frecuencias acumuladas F(x), despejamos 1-F(x) y tomamos dos veces logaritmos neperianos.

F(x)=1exp[ ( x c ) k ] ln[ ln( 1F(x) ) ]=kln(x)kln(c)

Esta ecuación es la de una línea recta de pendiente a y de ordenda en el origen b.

z=au+b

donde u y z son variables

z=ln[-ln(1-F(x))]
a=k
u=lnx
b=-klnc

Hay que tener en cuenta que se produce un error cuando la frecuencia acumulada F(xn)=1, ya que MATLAB tiene que calcular el logaritmo neperiano de 0 que da -Inf. Antes del realizar el ajuste tenemos que eliminar el último dato del vector de frecuencias acumuladas y su correspondiente valor en el vector de velocidades. Véase el código marcado en letra negrita. La función MATLAB polyfit ajusta los pares de datos a un polinomio de grado uno, es decir, a una línea recta.

clear,clc
velocidad=xlsread('WhiteDeer2013','Mar','F2:F745');
%interpolar si es necesario
if any(isnan(velocidad)) %si hay algún NaN
    x=1:length(velocidad);
    i=find(~isnan(velocidad));
    velocidad=interp1(x(i),velocidad(i),x);
end
%histograma
x=0.5:1:max(velocidad);
horas=hist(velocidad,x);
frec=horas/sum(horas); %frecuencias
%convierte a frecuencias acumuladas y ajusta a una línea recta
fcum=cumsum(frec);
z=log(-log(1-fcum));
z(length(z))=[];
u=log(x);
u(length(z))=[];
p=polyfit(u,z,1) %p(1) es la pendiente y p(2) es la ordenada en el origen
hold on 
%representa la juste
plot(u,z,'ro','markersize',4,'markerfacecolor','r')
z=p(1)*u+p(2);
plot(u,z,'b')
title('Ajuste lineal')
xlabel('x')
ylabel('y')
hold off

    p =    1.8253   -3.8557

La pendiente de la recta de ajuste p(1) es el parámetro k. El parámetro c se calcula a partir de k y la ordenada en el origen p(2).

k=p(1)
c=exp(-p(2)/k)

k=1.8253 y c=8.2680

que difieren ligeramente de los calculados mediante el ajuste no lineal

Nos queda añadir al script la representación de las frecuencias y la función de Weibull con los parámetros calculados por el procedimiento de ajuste lineal para producir una figura similar a la del apartado anterior.

El código que añadimos al script para producir esta figura es el siguiente

%representación gráfica de frecuencias y función de Weibull
figure
k=p(1)
c=exp(-p(2)/k)
f=@(x) (k/c)*((x/c).^(k-1)).*exp(-(x/c).^k); %función de Weibull
hold on
%diagrama de barras
bar(x,frec,'c');
%representa la curva de juste
x=linspace(0,max(velocidad),100);
plot(x,f(x),'r')
title('Ajuste a la función Weibull')
xlabel('Velocidad')
ylabel('Frecuencia')
hold off

Se producen dos figuras: la recta de ajuste y la representación del diagrama de frecuencias junto a la función Weibull que mejor ajusta. Nos devuelve los valores de los parámetros k y c.

k =    1.8253
c =    8.2680

Comparación de las medidas y el modelo estadístico

En este sección vamos a comparar los resultados estadísticos obtenidos con las medidas de la velocidad del viento con los proporcionados por el modelo estadístico de Weibull que hemos ajustado a las frecuencias por dos procedimientos distintos: no lineal y lineal.

Valor medio y desviación estándar

Calculamos el valor medio de las 744 medidas de la velocidad del viento a lo largo del mes de marzo con las funciones mean y std de MATLAB.

Utilizando el modelo estadístico de Weibull

El valor medio de la velocidad del viento se obtiene:

<x>=c·Γ(1+1/k)

La desvación estándar σ

σ 2 = 0 ( x < x > ) 2 · f ( x ) · d x = c 2 [ Γ ( 1 + 2 k ) Γ 2 ( 1 + 1 k ) ]

En la sección anterior hemos obtenido los valores de los parámetros k y c:

Creamos el script para comparar los resultados

clear,clc
velocidad=xlsread('WhiteDeer2013','Mar','F2:F745');
%interpolar si es necesario
if any(isnan(velocidad)) %si hay algún NaN
    x=1:length(velocidad);
    i=find(~isnan(velocidad));
    velocidad=interp1(x(i),velocidad(i),x);
end
%media y desviación estándar de las medidas de la velocidad del viento
media=mean(velocidad)
estandar=std(velocidad)
%Modelo estadístico de Weibull
%ajuste no lineal
k=2.0486,c=9.4165
media=c*gamma(1+1/k)
estandar=c*sqrt(gamma(1+2/k)-gamma(1+1/k)^2)
%ajuste lineal
k=1.8253,c=8.2680
media=c*gamma(1+1/k)
estandar=c*sqrt(gamma(1+2/k)-gamma(1+1/k)^2)
media =    8.1741
estandar =    3.9761

k =    2.0486  c =    9.4165
media =     8.3421
estandar =    4.2674

k =    1.8253  c =  8.2680
media =     7.3479
estandar =    4.1711

Tiempo de funcionamiento de un aerogerador

Un aerogerador funciona cuando la velocidad del viento es superior a un valor mínimo (cut-in) e inferior a un valor máximo (cut-out). Para un hipotético aerogenerador estas velocidades son 4 m/s y 18 m/s, respectivamente. Calcular el número de horas a lo largo del mes de marzo que el aerogenerador ha estado generando energía eléctrica.

A partir de las medidas de la velocidad del viento guardadas en el vector velocidad calculamos el número de valores (equivalente a horas) que son mayores que 4 y menores que 18.

En el modelo estadístico de Weibull, la probabilidad de que la velocidad del viento x esté en el intervalo comprendido entre x0 y x1 se calcula mediante la siguiente expresión

P ( x 0 x x 1 ) = x 0 x 1 f ( x ) · d x = exp [ ( x 0 c ) k ] exp [ ( x 1 c ) k ]

Hacemos estos cálculos con MATLAB

clear,clc
velocidad=xlsread('WhiteDeer2013','Mar','F2:F745');
%interpolar si es necesario
if any(isnan(velocidad)) %si hay algún NaN
    x=1:length(velocidad);
    i=find(~isnan(velocidad));
    velocidad=interp1(x(i),velocidad(i),x);
end
%horas a partir de las medidas de la velocidad del viento
x0=4;x1=18;
horas=sum(velocidad>=x0 & velocidad<=x1)
%Modelo estadístico de Weibull
%ajuste no lineal
k=2.0486,c=9.4165
prob=exp(-(x0/c)^k)-exp(-(x1/c)^k);
horas=length(velocidad)*prob
%ajuste lineal
k=1.8253,c=8.2680
prob=exp(-(x0/c)^k)-exp(-(x1/c)^k);
horas=length(velocidad)*prob
>> weibull_6
horas =    622

k =    2.0486  c =    9.4165
horas =  608.6147

k =    1.8253  c =     8.2680
horas =  558.5172

El ajuste no lineal parece que proporciona unos resultados más cercanos a los que se obtienen directamente a partir de las medidas de las velocidades del viento.