Simulacion planta solar atraves de matlab

Descripción del sistema
Una empresa dedicada a la instalación de plantas solares está desarrollando un sistema de control de temperatura para las plantas de colectores solares cilindro-parabólicos. Este tipo de plantas tiene como objetivo calentar un líquido, típicamente un aceite, que circula por una tubería concentrando sobre ella los rayos del sol. Esto se hace gracias al uso de espejos cilindro-parabólicos situando la tubería en el foco de la sección parabólica.
En la figura 1 puede observar una de estas plantas.

Fig 1. - Planta de colectores cilindroparabólicos
Este tipo de plantas solares consta de varias líneas de tuberías y espejos dispuestas en paralelo, de forma que todas las líneas toman como entrada el líquido frío que se extrae un depósito acumulador y tras calentarse en la línea, se almacena en el depósito acumulador. Con el fin de controlar la temperatura de salida de cada línea
T
s
Ts
, se plantea regular el ángulo de elevación de los colectores,
θ
θ
, para modificar la cantidad de radiación efectiva absorbida por cada tramo de colector.
De este modo, cuando la elevación de los colectores,
θ
θ
, coincide con la del sol en una instante determinado,
θ
S
O
L
θSOL
, la concentración de los rayos solares es máxima y así lo es también la cantidad de energía absorbida por el fluido circulante en los colectores. Desviaciones de
θ
θ
respecto a
θ
S
O
L
θSOL
permiten modular la cantidad de energía absorbida y, por tanto, la temperatura de salida,
T
s
Ts
de los colectores. El sistema que se propone controlar consiste en una línea compuesta de dos segmentos cilindroparabólicos en serie.
La temperatura de entrada del primer segmento,
T
e
Te
, es conocida. Al estar dispuestos en serie, la temperatura de salida del primer segmento,
T
i
Ti
, es la misma que de la de entrada del segundo segmento. El objetivo de control es regular la temperatura de salida del segundo segmento,
T
s
Ts
, modificando el ángulo de inclinación de los colectores,
θ
θ
.
Este ángulo de inclinación se regula a través de la tensión de alimentación,
u
u
, de un motor eléctrico que actúa sobre la estructura de soporte de los segmentos. Es por tanto esta tensión,
u
u
, la señal de control que emplearemos para regular el sistema. La figura 2 a continuación, ilustra la disposición de los elementos de la planta.


Fig 2. - Sistema a controlar
Modelo de la planta
Las ecuaciones que modelan el comportamiento de la planta vienen descritas por:
ρ
C
p
V
d
T
i
(
t
)
d
t

=
Q
C
p
(
T
e
(
t
)
−
T
i
(
t
)
)
+
I
(
θ
(
t
)
)
η
S
e
−
H
(
T
i
(
t
)
+
T
e
(
t
)
2

−
T
a
)
ρCpVdTi(t)dt=QCp(Te(t)−Ti(t))+I(θ(t))ηSe−H(Ti(t)+Te(t)2−Ta)
ρ
C
p
V
d
T
s
(
t
)
d
t

=
Q
C
p
(
T
i
(
t
)
−
T
s
(
t
)
)
+
I
(
θ
(
t
)
)
η
S
e
−
H
(
T
i
(
t
)
+
T
s
(
t
)
2

−
T
a
)
ρCpVdTs(t)dt=QCp(Ti(t)−Ts(t))+I(θ(t))ηSe−H(Ti(t)+Ts(t)2−Ta)
J
d
2
θ
(
t
)
d
t
2

+
B
d
θ
(
t
)
d
t

+
g
m
e
l
c
⋅
sin

(
θ
)
=
K
t
u
(
t
)
Jd2θ(t)dt2+Bdθ(t)dt+gmelc⋅sin⁡(θ)=Ktu(t)
En estas ecuaciones, el término
I
(
θ
(
t
)
)
I(θ(t))
representa el valor de la radiación solar efectiva que recibe cada segmento y que viene dada por:
I
(
θ
(
t
)
)
=
I
N
cos

(
θ
(
t
)
−
θ
S
O
L
)
I(θ(t))=INcos⁡(θ(t)−θSOL)
De este modo, las dos primeras ecuaciones describen la dinámica asociada a la evolución de las temperaturas de los segmentos y la última, la dinámica del actuador eléctrico. El valor y explicación de los parámetros del sistema dependen el DNI del alumno y se muestra en la tabla siguiente.

Tabla de parámetros




Parámetro
Valor
Unidad
Descripción
g
g
9.81
m
/
s
2
m/s2
Aceleración de la gravedad
C
p
Cp
2200
(
J
/
(
K
k
g
)
(J/(Kkg)
Calor específico del líquido
ρ
ρ
750
k
g
/
m
3
kg/m3
Densidad del líquido
H
H
350
W
/
K
W/K
Coeficiente global de pérdidas
I
n
In
850
W
/
m
2
W/m2
Radiación solar efectiva
T
e
Te
450
K
K
Temperatura de entrada a la línea
Q
Q
12
K
g
/
s
Kg/s
Caudal másico de líquido en la línea
η
η
0.85
−
−
Eficiencia óptica del espejo
S
e
Se
250
m
2
m2
Superficie de los espejos de cada segmento
θ
S
O
L
θSOL
0.78
r
a
d
rad
Elevación del sol
m
e
me
220
K
g
Kg
Masa de cada segmento
l
c
lc
0.65
m
m
Distancia del cdg del colector al eje rotación espejos
k
t
kt
100
N
m
/
A
Nm/A
Cte de par del motor actuador
T
a
Ta
310
K
K
Temperatura ambiente
V
V
1.05
m
3
m3
Volumen interior de la tubería
u
0
u0
7.014150000000001
V
V
Tensión de alimentación en el equlibrio
J
J
1.5
N
m
s
2
Nms2
Inercia del conjunto motor-segmentos
B
B
3
N
m
s
Nms
Fricción viscosa del conjunto motor-segmentos
Simulación del sistema
Con el fin de simular el comportamiento de la planta a controlar, se dispondrá de un modelo de Matlab/Simulink que puede descargar en la siguiente dirección

Fichero .zip con modelo Simulink del Sistema (Nueva versión corregida de 10/05/2018)

Fichero entregaMatlab.m

Las simulaciones de la planta así como los ensayos de los controladores se deben realizar utilizando este modelo, que a efectos prácticos se considera como si fuese la propia planta. Recuerde introducir su DNI en el bloque de simulink que simboliza la planta. Para ello haga doble click sobre el bloque del modelo.

Trabajo a realizar
Dado el sistema anteriormente descrito y facilitado mediante un bloque de Simulink, se pide:
Determinar analíticamente las ecuaciones que describen la característica estática del sistema. Para ello, siga los siguientes pasos
Calcule en primer lugar la expresión que relaciona el ángulo de elevación de los colectores en el equlibrio,
θ
0
θ0
, con la tensión de entrada al motor en el equilibrio,
u
0
u0
. Razone a partir de la expresión obtenida qué rango de valores de
u
0
u0
son admisibles para que el sistema alcance un ángulo de elevación
θ
0
θ0
estable.

Calcule en segundo lugar la ecuación que relaciona la temperatura de salida de los colectores en el equilibrio,
T
s
0
Ts0
, con la elevación de los mismos,
θ
0
θ0
.

Use las fórmulas anteriores para calcular algunas parejas de valores (
u
0
u0
,
T
s
0
Ts0
) en el equilibrio. Calcule cinco parejas de valores apropiados. Compruebe experimentalmente (con el modelo Simulink proporcionado) que los valores obtenidos con la fórmula analítica y el modelo coinciden.
Representar gráficamente la característica estática del sistema.

Nota: Use tiempos de simulación de al menos 1000 seg para garantizar que el sistema alcanza el régimen permanente. Los resultados obtenidos de los 5 puntos elegidos, deben transcribirse en el fichero entregaMatlab.m como corresponda.

Linealizar el sistema en torno al punto de equilibrio obtenido para
u
0
u0
(ver tabla de parámetros) y obtener las siguientes funciones de transferencia (observe el esquema del apartado siguiente):
G
11
(
s
)
=
T
s
(
s
)
T
e
(
s
)

G11(s)=Ts(s)Te(s)
G
12
(
s
)
=
T
s
(
s
)
θ
(
s
)

G12(s)=Ts(s)θ(s)
G
2
(
s
)
=
θ
(
s
)
U
(
s
)

G2(s)=θ(s)U(s)
Nota: Escriba la funciones de transferencia en el fichero entregaMatlab.m

Compare en este apartado el comportamiento del sistema no lineal y el linealizado. Partiendo del esquema Simulink que se le proporciona (el sistema no lineal), añada los bloques necesarios para simular el comportamiento del sistema linealizado. Observe que puede obtener el sistema lineal a partir de la funciones de transferencia del apartado anterior considerando el esquema siguiente

Suponiendo que la temperatura de entrada a los colectores,
T
e
(
t
)
Te(t)
, es constante, y que el sistema opera en el punto de funcionamiento proporcionado, aplique un salto de un
5
%
5%
y otro de
−
5
%
−5%
del valor
u
(
t
)
u(t)
correspondiente al punto de funcionamiento,
u
0
u0
.
Para estos dos experimentos calcule a partir de las simulaciones:
Ganancias estáticas de los sistemas no lineal y linealizado
Sobreoscilación (si la hubiese) para ambos sistemas
Tiempo de subida para ambos sistemas (Definido como el tiempo que se tarda en alcanzar del
5
%
5%
al
95
%
95%
del valor final de la respuesta).

Nota: Transcriba estos valores al fichero entregaMatlab.m

Obtenga la función de transferencia del sistema,
G
(
s
)
=
G
2
(
s
)
⋅
G
12
(
s
)
G(s)=G2(s)⋅G12(s)
, que representa la evolución de la temperatura de salida,
T
s
(
t
)
Ts(t)
, frente a la tensión de alimentación del motor actuador
u
(
t
)
u(t)
, suponiendo
T
e
(
t
)
Te(t)
constante (
T
e
(
s
)
=
0
Te(s)=0
).
Analice los polos/ceros de
G
(
s
)
G(s)
y determine si es posible encontrar una función de transferencia,
G
s
(
s
)
Gs(s)
, simplificada despreciando las dinámicas rápidas y/o los ceros que no afecten significativamente a la dinámica. Razone la respuesta.

Para la función de transferencia,
G
s
(
s
)
Gs(s)
, simplificada en el apartado anterior diseñe un controlador PID por métodos analíticos que satisfaga las siguientes especificaciones
Error en régimen permanente ante escalón cero
Sobreoscilación inferior al
10
%
10%
Tiempo de subida menor o igual a 150 seg

Nota: Puede despreciar a efectos de diseño los posibles ceros que aparezcan en la función
G
s
(
s
)
Gs(s)
.

Verifique el comportamiento del controlador diseñado sobre el modelo no lineal. Emplee para las simulaciones el esquema Simulink del apartado 3, cerrando el lazo de control de manera conveniente.
Suponga de nuevo
T
e
(
t
)
Te(t)
constante (
T
e
(
s
)
=
0
Te(s)=0
), y aplique como referencia un escalón de un
5
%
5%
sobre el valor de equilibrio
T
s
0
Ts0
calculado en el punto 1.
Compruebe sobre las simulaciones el error en régimen permanente observado, la sobreoscilación y el tiempo de subida.

Nota: Transcriba al fichero entregaMatlab.m la expresión del controlador diseñado, y los valores de error en régimen permanente observado, la sobreoscilación y el tiempo de subida.


Diseñe ahora un controlador empleando el método del lugar de las raices que satisfaga las mismas especificaciones que el apartado anterior. Para ello siga los siguiente pasos:
Dibuje con MATLAB el lugar de las raices del sistema
Convierta las especificaciones a restricciones geométricas
Diseñe el controlador más simple posible (de menor a mayor complejidad: P, pd, pi, pid) que satisfaga las especificaciones del apartado 5. Si no fuese posible satisfacer todas las especificaciones, explique por qué y diseñe en su lugar el controlador más rápido posible que satisfaga el error ante escalón nulo y sobreoscilación inferior al
10
%
10%
.



Al igual que en el apartado anterior verifique el comportamiento del controlador diseñado sobre el modelo no lineal. Suponga de nuevo
T
e
(
t
)
Te(t)
constante (
T
e
(
s
)
=
0
Te(s)=0
), y aplique un cambio de referencia en escalón de un
5
%
5%
sobre el valor de equilibrio
T
s
0
Ts0
calculado en el punto 1.
Nota: Consigne estos valores al fichero entregaMatlab.m


Analice ahora el problema de rechazo de perturbaciones. Suponga que la referencia,
R
(
s
)
R(s)
es cero y que se aplica un incremento en escalón para la temperatura de entrada,
T
e
Te
, de
10
o
10o
. Tome el esquema de Simulink del apartado anterior con el mismo controlador diseñado y determine sobre la respuesta:
Error en régimen permanente
Sobreoscilación (Si la hubiese)
Tiempo de subida

Nota: Consigne estos valores al fichero entregaMatlab.m

Memoria del proyecto
Los resultados del proyecto se explicarán en una memoria justificativa a realizar de forma individual.
En ésta se deben presentar los procedimientos de cálculo y resultados obtenidos de cada uno de los puntos del trabajo a realizar. Además se deben ilustrar los resultados obtenidos gráficamente; por ejemplo, las gráficas de la respuesta temporal del sistema, diagrama del Lugar de las Raices obtenido, etc.

Los datos numéricos del proyecto en cada apartado deberán rellenarse en el fichero adjunto de nombre
entregaMatlab.m
.
Entrega del trabajo
Tanto la memoria del proyecto, como los ficheros de MATLAB empleados en la realización del mismo (especialmente el fichero
entregaMatlab.m
arriba mencionado), se entregarán a través de Enseñanza Virtual en el enlace siguiente (debe identificarse antes):

Actividad de entrega del trabajo en Enseñanza Virtual

Los resultados se entregarán en un fichero comprimido .zip que contenga el documento de la memoria en pdf y el fichero entregaMatlab.m.

La fecha límite de entrega es el 31 de mayo de 2018 a las 23:59