Offshore Engineering for non-
conventionalRenewable Energies(NCRE)
Introduction to Matlab
Juan Gabriel Rueda Bayona, PhD
Content
1 Relationship of the NCRE with ecosystems and climate.
1. Non-conventional Renewable Energy: NCRE.
2. Oceanography.
3. Hydrology.
2 Energy potential estimation.
1. Introduction to Matlab
2. Time series analysis
3. Image processing
3 Fundamentals in marine structures desing
1. Offshore structures for marine energy
2. Offshore wind turbines.
Matrices en Matlab
APLICACIÓN. Programa: EJEMPLO_2_01.m
Construcción de una tabla de frecuencias y un histogramas usando de
la dirección del viento (fichero [Link])
Paso 1: Lectura de datos y seleccionar la columna donde esta la dirección del viento.
clear all, close all, clc
path = 'c:\CURSOS\EJEMPLO_PARTE_2\datos\';
file = [path '[Link]'];
X = load(file);
dir = X(:,3); % En la columna 3 esta la dirección.
Paso 2: Seleccionar los intervalos
N : 337.5 - 22.5
NE : 22.5 - 67.5
E : 67.5 - 112.5
SE : 112.5 - 157.5 :
S : 157.5 - 202.5
SW : 202.5 - 247.5
W : 247.5 - 292.5
NW : 292.5 - 337.5
rangoNorte = [337.5 360 0 22.5]
rangoRestn = [ 22.5 67.5 112.5 157.5 202.5 247.5 292.5 337.5];
rango_str = char('N','NE','E','SE','S','SW','W','NW');
4
Paso 3: Determinar el número de observaciones en cada intervalo definido.
% Primero se determina la cantidad de datos en la dirección Norte: 337.5 a 22.5
j = find( dir >= rangoNorte(3) & dir < rangoNorte(4) | ...
dir >= rangoNorte(1) & dir <=rangoNorte(2));
% Los resultados obtenidos los colocaremos en la variable Fi, ordenados en
% columnas.
Fi(1,1) = length(j);
% Determina la cantidad de datos en las direcciones faltantes 22.5 a 337.5
for k = 1 : 7;
j = find(dir >= rangoRestn(k) & dir < rangoRestn(k+1));
Fi(k+1,1) = length(j);
end
Paso 4: Determinar las frecuencias relativas y las frecuencias relativas acumuladas en
cada intervalo definido
% Calcula la frecuencia relativa (FiR), frecuencia acumulada (FiA),
% frecuencia relativa acumulada (FrA)
FiR = Fi / sum(Fi);
FiA = cumsum(Fi);
FrA = cumsum(FiR);
Paso 5: Construir la tabla de frecuencias.
% Genera una tabla de salida. Todas las variables calculadas las colocamos en
% una matriz de nombre tabla.
tabla = [Fi FiR FiA FrA];
% Imprime la tabla creada en pantalla.
fprintf('TABLA DE FRECUENCIA: DIRECCION DEL VIENTO \n');
fprintf(' Dir Fi FiR FiA FrA\n')
for k=1:8
fprintf(' %s %6d %7.3f %5d %7.3f\n',rango_str(k,:),tabla(k,:));
end
Paso 6: Crear el histograma.
% Grafica el histograma
figure('Units', 'Normalized','Position',[0.1 0.13 0.5 0.35],...
'PaperPositionMode','auto','color','w')
h = bar(1:8,FiR,'w');
set(h,'BarWidth',1,'linewidth',1.3)
set(gca,'box','off','tickdir','out','ticklength',[0.02 0.3],...
'xlim',[0 9],'xtick',[1:8],'xticklabel',rango_str,...
'fontsize',13,'ylim',[0 0.6],'ytick',[0:0.1:0.6])
ylabel('Porcentaje','fontsize',15)
% Almacena la figura en un archivo
print -dpng ejemplo_2_01
RESULTADOS
TABLA DE FRECUENCIA: DIRECCION DEL VIENTO
Dir Fi FiR FiA FrA
N 487 0.167 487 0.167
NE 138 0.047 625 0.214
E 97 0.033 722 0.247
SE 181 0.062 903 0.309
S 1656 0.567 2559 0.876
SW 206 0.071 2765 0.947
W 67 0.023 2832 0.970
NW 88 0.030 2920 1.000
7
APLICACIÓN. Programa: EJEMPLO_2_02.m
Usa el fichero “[Link]”, para generar una tabla de frecuencias de la
magnitud del viento, y además el histograma.
Intervalo Marca de Clase Fr. Abs
0- 4 2 333
4- 8 6 1038
8-12 10 616
12-16 14 521
16-20 18 209
20-24 22 148
24-28 26 44
28-32 30 7
32-36 34 3
11
Medidas de tendencia central
x1, x2 , x3 , K K , xn
Media Aritmética Media cuadrática
1 n
1 n 2
x=
n
x i x=
n i=1
xi
i=1
y = mean(x)
y = nanmean(x)
Medidas de dispersión
Desviación estándar
N
s= (xi − x ) Varianza = s 2
N −1
2
1
i=1
Covarianza
n s = std(x)
1
c xy =
N −1 i=1
(xi − x )(yi − y ) s
s2
=
=
nanstd(x)
var(x)
s2 = nanvar(x)
c = cov(x,y)
10
APLICACIÓN. Programa: EJEMPLO_2_03.m
Calcula el promedio climatológico cada 6 horas y añade la desviación
estándar a cada lado de la media. El objetivo es observar a que horas
del día se produce la mayor variabilidad de los vientos.
Cuartiles
Los cuartiles q1, q2 , q3 de un conjunto de datos, son números tal que
dividen en cuatro partes iguales dicho conjunto de datos. El rango
intercuartil es la longitud donde está contenido el 50% central de los
datos.
Rintercuartil = q3 − q1
int es la parte entera del valor y n es el número de observaciones
12
q =
quantile(x,[0.25,0.50,0.75]) q1
= q(1);
q2 = q(2);
q3 = q(3);
13
APLICACIÓN. Programa: EJEMPLO_2_04.m
clear all, close all, clc
path = 'c:\CURSOS\EJEMPLO_PARTE_2\datos\';
file = [path '[Link]'];
DATA = load(file);
X = DATA(:,4);
% Determina el cuartil 1 (q1), cuartil 2 (q2) y el cuartil 3 (q3)
q = quantile(X,[0.25 0.50 0.75 ]);
q1 = q(1);
q2 = q(2);
q3 = q(3);
Rango_Intercuartil = q3-q1;
fprintf(' RESULTADOS\n')
fprintf(' q1 (percentil 25) = %2d nudos \n',q1)
fprintf(' q2 (percentil 50) = %2d nudos \n',q2)
fprintf(' q3 (Percentil 75) = %2d nudos \n',q3)
fprintf(' Rango_Intercuartil = %2d nudos \n',Rango_Intercuartil)
RESULTADOS
q1 (percentil 25) = 5 nudos
q2 (percentil 50) = 8 nudos
q3 (Percentil 75) = 13 nudos
Rango_Intercuartil = 8 nudos
14
Diagrama de Cajas: BOXPLOT
X max Es un gráfico que suministra
Datos fuera información sobre los valores
de rango mínimo y máximo, los cuartiles
q1, q2 (Mediana) y q3, y sobre la
Lsup existencia de valores atípicos y la
simetría de la distribución.
q3 Construcción:
X min, X max, q1 , q2 , q3
q2
Linf = q1 − f (q3 − q1 )
q1 Lsup = q3 + f (q3 − q1 )
f = 0.75,ó f = 1.5
Linf
Calcular todos los valores atípicos fuera del rango: Linf ,1L9 sup
h = boxplot(X,'color',color,...
'symbol',marca,'widths',ancho)
X = Matriz de datos.
color = Tipo de color para las líneas y el bigote.
symbol = Tipo de marca para visualizar los valores atípicos
widths = Ancho del boxplot.
Interpretación
Mientras más larga la caja y los bigotes, más dispersa es la
distribución de datos. La línea que representa la mediana indica la
simetría. Si está relativamente en el centro de la caja la distribución
es simétrica. Si por el contrario se acerca al primer o tercer cuartil, la
distribución pudiera ser sesgada a la derecha (asimétrica positiva) o
sesgada a la izquierda (asimétrica negativa) respectivamente.
16
APLICACIÓN. Programa: EJEMPLO_2_05.m
Construcción de un diagrama BOXPLOT.
clear all, close all, clc
path = 'c:\CURSOS\EJEMPLO_PARTE_2\datos\';
file = [path '[Link]'];
DATA = load(file);
X = DATA(:,4);
figure('Units', 'Normalized',...
'Position',[0.1 0.13 0.3 0.4],...
'PaperPositionMode','auto',...
'color','w')
h=boxplot(axes,X,'color','k','symbol','o',...
'widths',0.1)
set(gca,'xtick',[ ],'fontsize',13,...
'ylim',[-1 40],'ytick',[Link],'box','on',...
'ticklength',[0.025 0.3])
ylabel('Magnitud (nudos)','fontsize',15)
title({'Diagrama BoxPlot de la magnitud',...
'del viento'},'fontsize',15)
print -dpng EJEMPLO_2_04
17
APLICACIÓN. Programa: EJEMPLO_2_06.m
Ajuste de Curvas
Principales curvas
y = a +bx y = Acos(t − ) y = aebx
Donde a,b, A, son las
y = axe −bx
constantes a determinar.
20
Método de mínimos cuadrados
Método de mínimos cuadrados es una técnica que nos permite
encontrar los coeficientes del modelo teórico, minimizando el error que
se produce entre las observaciones y el modelo elegido.
e(x j ) = h(x j ) − g( x j )
h(x j ) son nuestras observaciones y g(x j ) es la función que deseamos
ajustar a nuestras observaciones.
m
g( x j ) = c1 f1 (x j ) + c2 f2 (x j ) + L + cm fm (x j ) g(x) = ck f k (x)
k =1
las constantes: c1,c2 , K , cm son los coeficientes del modelo.
21
Transformando esta expresión a un sistema matricial:
F T H = (F T F )C
Despejando C, finalmente nos queda las soluciones de este sistema de
ecuaciones que contiene los coeficientes del modelo teórico.
C = (F F ) F T H
T −1
Esta es la solución fundamental que se utiliza para determinar los
coeficientes de cualquier función que desean ajustar a sus observaciones
Para calcular valores pronosticados usando el modelo:
Ymodelo = FC
24
Error del modelo
residual = Yobs −Ymodelo
residual = 0
Varianza residual
APLICACIÓN. Programa: EJEMPLO_2_09.m
Ajustar modelo cuadrático a nuestras observaciones
% datos.
x = [-3 0 2 4];
y = [3 1 1 3];
Visualizar los datos
Decidir el modelo Determinar los coeficientes
que definen el modelo
C = (F T F ) F T H
Modelo −1
y = c0 + c1 x + c2 x2
Calcular la matriz F
y = c0 f1 (x j ) + c1 f 2 (x j ) + c2 f 2 (x j )
f3 (x j ) = x 2
f1 (x j ) = 1 f 2 (x j ) = x
Reemplazar nuestros datos
Evaluar los valores
del eje-Y en vector H
de nuestros datos
en estas funciones.
Determinar las soluciones del
sistema de ecuaciones
C = F F( T
)
−1
FTH
clear all, close all, clc
% datos.
x = [-3 0 2 4];
y = [3 1 1 3];
n = length(x);
% Definimos la matriz "F"
F = [ones(n,1) x(:) x(:).^2];
% Definimos la matriz H
H = y(:);
% Calculamos los coeficientes
C = (F T F ) F T H
−1 % del modelo.
C = inv(F'*F)*F'*H;
Coeficientes del
modelo c0 = C(1) =
0.851
c1 = C(2) = -0.192
c2 = C(3) =
0.178
Nuestro modelo seleccionado queda:
Y = 0.85 − 0.19x + 0.18x2
modelo
Pronosticar un valor
x = 5; y = 0.85 − 0.19x + 0.18x2
(
X pronos = 1 x x2 )
Ypronos = X pronosC
Ymodelo = XC
29
APLICACIÓN. Programa: EJEMPLO_2_10.m
Ajustar modelo lineal a nuestras observaciones. “datos_ejemplo_2_10.dat”
Modelo
y = c0 + c1 x
Determinar las soluciones del
sistema de ecuaciones
C = (F T F ) F T H
−1
% Definimos nuestras variables
x = DATOS(:,1);
y = DATOS(:,2);
n = length(x);
% Calculamos los coeficientes
H = y;
F = [ones(n,1) x];
C = inv(F'*F)*F'*H;
y = −6.6269 + 8.1185x
Coeficientes del
modelo c0 = C(1) = -
6.6269
c1 = C(2) = 8.1185
Ymodelo = XC
Ymodelo = −6.6269 + 8.1185x
Ajuste Potencial.
Linealización del modelo
y =c 0 x c1 Aplicaremos la función logarítmica para
linealizar esta ecuación
(
log(y )= log c0 x c1 ) renombrando Y = log( y)
estas nuevas variables
log(y ) = log(c0 )+ c1 log(x) a = log(c0 )
X =log(x)
Y a X
Y = a + c1 X Modelo linealizado
APLICACIÓN. Programa: EJEMPLO_2_11.m
Ajustar modelo potencial a nuestras observaciones. “datos_ejemplo_2_11.dat”
Modelo
y = c0 x c1
X =log(x)
Y = log( y)
x = DATOS(:,1);
y = DATOS(:,2); y = 1.8582 + 0.21x
n = length(x);
Y = log(y);
X = log(x);
F = [ones(n,1) X];
c0 = e1.85
C = inv(F'*F)*F'*Y;
a = C(1);
c1 = C(2);
c0 = exp(a);
Coeficientes del modelo
lineal a= C(1) = 1.8c583
c1 = C(2) = 0.2100
Coeficientes del modelo
Potencial c0 = exp(a) = 6.4125
c1 = C(2) = 0.2100
y = 6.4x0.21
y = 6.4125x 0.21
34
Función en MATLAB para el calculo de los coeficientes
de un polinomio de orden “m”
Ajusta un polinomio de orden “m”
c = polyfit(x,y,m)
y = cm x m +c m−1 x m−1 + L + c1 x + c 0
Evalúa el polinomio ajustado
y = polyval(c,x)
m = 1, Regresión Lineal y = c1x + c0
m = 2, Polinomio cuadrático y = c2x2 + c 1x + c 0
m = 3, Polinomio cúbico y = c3x3 + c 2x2 + c 1x + c 0
