Modelado y control de un reactor nuclear tipo Stellarator

Debido al continuo crecimiento de la demanda energética, se augura que entre los años 2030 y 2050 la humanidad entrará en un estado de déficit de energía que irá aumentando con los años. Esto hace que la búsqueda de nuevas fuentes de energía sostenibles, limpias y rentables sea una prioridad.

A día de hoy, existen muchas fuentes renovables de las que se ha conseguido extraer energía de manera más o menos exitosa, siendo la eólica y la solar las más explotadas en España.  Aunque el crecimiento de la explotación de este tipo de energías haya crecido en un 174% entre 1990 y 2014 todavía queda un largo camino que recorrer en este ámbito. Esto se debe primordialmente a que este tipo de fuentes de energía no son capaces de proporcionar la suficiente demanda de energía con regularidad ya que depende de las condiciones climatológicas, así como de la localización geográfica.

Aunque se esté avanzando por el buen camino, a lo que a energías renovables se refiere, todavía sigue siendo necesario encontrar una fuente de energía primaria que sea capaz de responder a la demanda de una forma constante y sin sobresaltos. Esto se debe a la característica de irregularidad que presentan las energías renovables en la generación eléctrica. Por ello, en la actualidad, están naciendo y consolidándose líneas de investigación que se centran en obtener energía mediante una fuente energética alternativa: la fusión nuclear.

La fusión nuclear es una reacción en la que dos núcleos de átomos ligeros, en general dos isótopos del hidrógeno (deuterio y tritio), se unen para formar otro núcleo más pesado, liberando partículas altamente energéticas en el proceso. Un ejemplo de este tipo de reacciones se encuentra en el sol donde constantemente se están produciendo reacciones entre átomos de hidrógeno para formar helio.

Para conseguir que las reacciones de fusión se den con éxito es necesario garantizar una serie de condiciones. Para empezar, hay que conseguir vencer la repulsión electroestática de los núcleos que se dese fusionar, para lo que hay que proporcionales mucha energía calentándolos a altas temperaturas. Al calentar tanto el gas que se usa como combustible en estas reacciones, éste pasa a un nuevo estado, el plasma. La necesidad del calentamiento acarrea el problema primordial que tiene a día de hoy la energía por fusión de núcleos, el rendimiento energético.

El rendimiento es un problema debido a que se gasta más energía calentando el gas de la que se produce mediante la fusión, lo que convierte a esta estrategia de generación energética en no-rentable. Cabe destacar que el obtener un beneficio energético de la fusión no sólo depende de la temperatura del plasma, sino también de su densidad y el periodo de tiempo durante el que se produce la reacción. Por ello, con el fin de reducir la cantidad de energía necesaria para el calentamiento de plasma es inevitable hacer un buen confinamiento del plasma.

En la actualidad, existen varios métodos para el confinamiento del plasma en un reactor de fusión, siendo el confinamiento magnético el más desarrollado. Este tipo de confinamiento consiste en aprovechar la polaridad y descompensación eléctrica del plasma, para conseguir retenerlo en un espacio concreto predeterminado (dentro del reactor pero sin entrar en contacto con las paredes del mismo) mediante campos magnéticos.

Este método de confinamiento es el que se está usando en el reactor de fusión tipo Stellarator que se encuentra en la UPV/EHU, y cuyo modelo ya está desarrollado y validado en el artículo Modelling and control of the UPV/EHU Stellarator. Debido a la necesidad de controlar las magnitudes de tensión e intensidad de las bobinas generadoras del campo magnético del Stellarator en este post se propone una forma de llevar a cabo ese contol.

2.  Validación de los modelos

Se dispone de 3 modelos distintos que caracterizan la dinámica del mismo reactor del tipo Stellarator. En este apartado se va a tratar la validación de estos modelos. Para ello, se utilizarán los tres modelos propuestos y el controlador tipo PID diseñado en el artículo Modelling and control of the UPV/EHU Stellarator.

Modelo caja blanca:

$$\left( \begin{matrix} \dot { { x }_{ 1 } }  \\ \dot { { x }_{ 2 } }  \end{matrix} \right) =\begin{pmatrix} -11.1 & 0 \\ 0 & -11.1 \end{pmatrix}·\left( \begin{matrix} { x }_{ 1 } \\ { x }_{ 2 } \end{matrix} \right) +\left( \begin{matrix} 47.62 \\ 0 \end{matrix} \right) { V }_{ in }$$

$$\left( \begin{matrix} { I }_{ OUT } \\ { V }_{ OUT } \end{matrix} \right) =\begin{pmatrix} 1 & 0 \\ 0.11 & 0.00728 \end{pmatrix}·\left( \begin{matrix} { x }_{ 1 } \\ { x }_{ 2 } \end{matrix} \right) $$

Modelo caja gris:

$$\left( \begin{matrix} \dot { { x }_{ 1 } }  \\ \dot { { x }_{ 2 } }  \end{matrix} \right) =\begin{pmatrix} -15.5 & 0 \\ 0 & -15.5 \end{pmatrix}·\left( \begin{matrix} { x }_{ 1 } \\ { x }_{ 2 } \end{matrix} \right) +\left( \begin{matrix} 69.3 \\ 0 \end{matrix} \right) { V }_{ in }$$

$$\left( \begin{matrix} { I }_{ OUT } \\ { V }_{ OUT } \end{matrix} \right) =\begin{pmatrix} 1 & 0 \\ 0.11 & 0.00632 \end{pmatrix}·\left( \begin{matrix} { x }_{ 1 } \\ { x }_{ 2 } \end{matrix} \right) $$

Modelo de caja negra:

$$\left( \begin{matrix} \dot { { x }_{ 1 } }  \\ \dot { { x }_{ 2 } }  \end{matrix} \right) =\begin{pmatrix} -29.7 & 12.71 \\ -74.41 & -29.3 \end{pmatrix}·\left( \begin{matrix} { x }_{ 1 } \\ { x }_{ 2 } \end{matrix} \right) +\left( \begin{matrix} 0.4717 \\ 2.815 \end{matrix} \right) { V }_{ in }$$

$$\left( \begin{matrix} { I }_{ OUT } \\ { V }_{ OUT } \end{matrix} \right) =\begin{pmatrix} 169 & -10.52 \\ 11.7 & 6.423 \end{pmatrix}·\left( \begin{matrix} { x }_{ 1 } \\ { x }_{ 2 } \end{matrix} \right) $$

Lo primero que se hará será comparar la respuesta de los modelos ante un escalón unitario de tensión a la entrada, para así comprobar que, efectivamente, los tres modelos están correctamente implementados y caracterizan la dinámica del mismo sistema:

El controlador PID propuesto presenta los siguientes valores de parámetros:

P = 2.9772 I = 58.228 D = 0.018578

El resultado obtenido es el siguiente:

Si se repara a los resultados obtenidos se puede ver cómo las salidas de los tres modelos son muy similares, pero no idénticas. En cuanto al rebose, el modelo de caja negra no presenta, el de caja gris tiene alrededor de un 2% y el de caja blanca un 20%. En lo referido al tiempo de

establecimiento, el del modelo de caja blanca es de 0.3s, el de caja gris de 0.2s y el de caja negra de 0.35s.

Si se analiza lo que ocurre con las corrientes se puede apreciar que existe diferencia incluso en el valor estacionario de la corriente. En el caso de modelo de caja blanca se ve como en estado estacionario alcanza 8.6A aproximadamente, mientras que el de caja gris alcanza los 9A y el modelo de caja negra llega a 8.8A. Además, como es lógico, la dinámica cambia dependiendo del modelo al igual que sucede con la tensión.

3.  Diseño de controladores tipo PID para cada modelo

Una vez comprobado que los modelos funcionan como deben, se pasará a diseñar un controlador tipo PID. Para ello, se pasarán los sistemas en variables de estado a funciones de transferencia. Una vez conseguidas dichas funciones de transferencia se diseñarán los controladores analíticamente utilizando el principio de la cancelación polo-cero.

Para conseguir las funciones correspondientes a cada uno de los modelos se ha utilizado el comando ss2tf de Matlab. 

Cabe destacar que dicha función calculará dos funciones de transferencia por cada uno de los modelos, esto se debe a que los modelos tienen dos salidas, intensidad y tensión, tal y como se ha comentado anteriormente. Por lo tanto, debido a que se controlará directamente el voltaje para controlar indirectamente la intensidad, sólo se utilizarán las funciones de transferencia referidas a la tensión como salida. En este caso dichas funciones de transferencia serán G_w, G_g y G_b. 

$${ G }_{ w }(s)\quad =\quad \frac { 5.552 }{ s\quad +\quad 11.1 } \quad \quad { { G }_{ g }(s)\quad =\quad \frac { 7.692 }{ s\quad +\quad 15.5 } \quad { G }_{ b }(s)\quad =\quad \frac { 23.6s\quad +\quad 892 }{ { s }^{ 2 }\quad +\quad 59.02s\quad +\quad 1817 } \quad \quad  }$$

Para comenzar el diseño de los controladores utilizando el método de asignación de polos y ceros, en primer lugar se seleccionará la estructura del controlador tipo PID dependiendo del número de polos y ceros que presenten las funciones de transferencia, y teniendo en cuenta que se desea eliminar el error del sistema en estado estacionario. Como se ha mencionado con anterioridad, este proceso se llevará a acabo analíticamente. A continuación se muestran los requisitos de la respuesta para los 3 modelos: 

t_s = 0.15s   M_p =  2%   t_subida = 0.086s

Modelo caja blanca

En primer lugar se diseñara el controlador para el modelo de caja blanca. Como se puede apreciar en la figura se trata de una función de transferencia de primer orden, por lo que se diseñará un controlador del tipo PI, con el fin de cancelar el polo e intentar cumplir las especificaciones de diseño.

La función de transferencia de la planta con el controlador tipo PI tendrá la siguiente forma en lazo abierto: 

$$G_{ w }(s)\quad =\quad \frac { { K }_{ p }·0.5\left( { T }_{ i }s\quad +\quad 1 \right)  }{ { T }_{ i }s(0.09s\quad +\quad 1) } $$

Por un lado, para cancelar el polo la T_i deberá ser 0.09. Por otro lado para conseguir el tiempo de establecimiento deseado en primer lugar se calculará el lazo cerrado de la función de transferencia:

$$G_{ wlc }(s)\quad =\quad \frac { { 0.5·K }_{ p }/{ T }_{ i } }{ s\quad +\quad { 0.5·K }_{ p }/{ T }_{ i } }$$

Teniendo en cuenta que en el caso de una función de transferencia de primer orden el tiempo de establecimiento, con criterio del 2%, es cuatro veces la constante de tiempo de la función de transferencia, la expresión para calcular la constante proporcional del controlador quedará así: 

$${ K }_{ P }\quad =\quad \frac { 4·{ T }_{ i } }{ { t }_{ s }·0.5 } \quad =\quad 4.79$$

Entonces, el PI para controlar la intensidad en el modelo de caja blanca tendrá esta forma:

$${ P }_{ Iw }(s)\quad =\quad 4.79\left( 1\quad +\quad \frac { 1 }{ 0.09s }  \right) $$

Modelo caja gris

Una vez diseñado el controlador para el modelo de caja blanca, se pasará a diseñar el controlador  para el modelo de caja gris. Si se analiza la función e transferencia de este modelo se ve claramente que al igual que con el modelo de caja blanca, este también es una función de transferencia de primer orden. Por lo tanto, se seguirá exactamente el mismo procedimiento para el diseño del controlador. Por simplicidad, solamente se mostrará la expresión final del controlador obtenido:

$${ P }_{ Ig }(s)\quad =\quad 3.49\left( 1\quad +\quad \frac { 1 }{ 0.065s }  \right) $$

Modelo caja negra

Para el diseño del controlador del modelo de caja negra, se debe tener en cuenta que el modelo en este caso no es de primer orden, por lo que en lugar de diseñar un PI se diseñara un PID, para poder cancelar así, los dos polos que presenta la función de transferencia del modelo de caja negra. El procedimiento a seguir será el mismo que se ha utilizado con los dos modelos anteriores, con la diferencia que los cálculos serán más complejos. Entonces, la función de transferencia en lazo abierto con el controlador tendrá la siguiente forma: 

$$PID(s)\quad =\quad \frac { { K }_{ p }{ T }_{ d }\left( { s }^{ 2 }+\frac { 1 }{ { T }_{ d } }s +\frac { 1 }{ { T }_{ i }{ T }_{ d } }  \right) \left( 23.6s+892 \right)  }{ s\left( { s }^{ 2 }+59.02s+1817 \right)  } $$

Para calcular los valores de las constantes de tiempo integral y derivativa, se igualarán los coeficientes de los ceros introducidos por el controlador y los polos del sistema.

$${ s }^{ 2 }+\frac { 1 }{ { T }_{ d } } s+\frac { 1 }{ { T }_{ i }{ T }_{ d } } ={ s }^{ 2 }+59.02s+1817$$

$${ T }_{ d }=\frac { 1 }{ 59.02 } \quad y\quad { T }_{ i }=\frac { 59.02 }{ 1817 } $$

Una vez obtenidas las constantes de tiempo integral y derivativa, se calculará la función e transferencia del sistema en lazo cerrado: 

$${ G }_{ blc }(s)\quad =\quad \frac { 1 }{ \left( \frac { 1+{ K }_{ p }{ T }_{ d }23.6 }{ { K }_{ p }{ T }_{ d }892 }  \right) s+1 } $$

Cabe destacar que para este último modelo, debido al cero que presenta, la ganancia proporcional se deberá disminuir un poco, ya que el efecto de este cero hará que el sistema sea más rápido y que presente más rebose. Por lo tanto, para conseguir que la salida se ajuste a las especificaciones la constante proporcional final la constante proporcional se disminuirá a la mitad: 

$${ K }_{ p }\quad =\quad \frac { 4/2 }{ { 892 }·t_{ s }{ T }_{ d }-4·23.6{ ·T }_{ d } } =2.9959$$

Para comprobar que los controladores funcionan correctamente, los tres modelos serán sometidos a una entrada escalón y el resultado se comparará con los resultados obtenidos en el artículo Modelling and control of the UPV/EHUStellarator.

Si se analizan los resultados obtenidos, se pude ver cómo los controladores propuestos en este trabajo son más rápidos, ya que los tres se ajustan perfectamente al tiempo de establecimiento especificado, 0.15s. Además, en los modelos de caja blanca y gris se ha conseguido eliminar por completo el rebose y en el modelo de caja negra reducirlo hasta el 1%.

Si se analiza lo que sucede con las corriente a través de las bobinas, como es lógico, la dinámica de las misma, con los controladores propuestos en este trabajo son más rápidas que en el artículo Modelling and control of the UPV/EHU Stellarator. Además al igual que las tensiones, los modelos de caja blanca y gris no presentan rebose mientras que el modelo de caja negra presenta un rebose del 1%. Aun así los valores de las corrientes en estado estacionario no varían respecto de los obtenidos con el controlador del artículo.