Del cálculo variacional a la teoría del control óptimo: Aplicación a sistemas lejos del equilibrio

DEL CÁLCULO VARIACIONAL A LA TEORÍA DEL CONTROL ÓPTIMO: APLICACIÓN A SISTEMAS LEJOS DEL EQUILIBRIO Natalia Ruiz Pino DEL CÁLCULO VARIACIONAL A LA TEORÍA DEL CONTROL ÓPTIMO: APLICACIÓN A SISTEMAS LEJOS DEL EQUILIBRIO Natalia Ruiz Pino Memoria presentada como parte de los requisitos para la obtención del título de Doble Grado en Física y Matemáticas por la Universidad de Sevilla. Tutorizada por Prof. Antonio Prados Montaño Dpto. Física Atómica, Molecular y Nuclear Índice general English Abstract 1 Resumen 3 I Del cálculo variacional a la teoría del control óptimo 5 1. Introducción. Cálculo de variaciones 7 2. Teoría del control óptimo 15 2.1. Problema de control óptimo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2. Tratamiento variacional del problema de control óptimo .......................................... 18 2.3. Elprincipiodelmáximo ................................ 23 2.3.1. Enunciado del principio del máximo . . . . . . . . . . . . . . . . . . . . . 24 2.3.2. El principio del máximo en acción: aparcando un coche con aceleración acotada ..................................... 25 II Aplicación a sistemas lejos del equilibrio 31 3. Minimización del tiempo de conexión 33 3.1. Mediosgranulares.................................... 33 3.2. Descripcióndelsistema................................. 35 3.3. Tratamiento matemático del sistema . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.3.1. Aplicación del principio del máximo a sistemas lineales: comprobación de hipótesis . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.3.2. Aplicación del principio del máximo a sistemas lineales: resolución del sistema y cálculo de trayectorias . . . . . . . . . . . . . . . 41 ii del cálculo variacional a la teoría del control óptimo: aplicación a sistemas lejos del eqilibrio 4. Análisis de los resultados 49 4.1. Tiempo de cambio y final como función de la intensidad del termostato . . . . . . 49 4.2. Trayectorias de la temperatura y la kurtosis . . . . . . . . . . . . . . . . . . . . . 52 4.3. Rango de validez de la aproximación lineal . . . . . . . . . . . . . . . . . . . . . . 54 5. Conclusiones 63 A. Apéndice. Valores límite de la intensidad del termostato 65 Anexo 67 English Abstract Throughout history, time has marked a determining factor in multiple processes, from optimization of tasks in assembly lines, search for routes of minimum time, etc. In this fierce fight against time, mathematics has developed a series of powerful tools to tackle all the problems that require minimizing time. In our case, the problem we deal with comes from a granular system that will be led from an initial non-equilibrium stationary state (NESS) to another NESS in the shortest possible time, by controlling the intensity of an external thermostat. In order to face this problem, it is necessary to become familiar with mathematical optimization methods that fit the theoretical framework imposed by our problem. This requirement lead us to study the so-called maximum principle, the fundamental pillar of the branch of mathematics known as optimal control theory. We also study its linear variants, which allow us to solve the problem when the final and initial states are close to each other. Although we do not present a rigorous demonstration of the principle of the maximum, we reach it in an intuitive way using the tools of variational calculus. Along the way, we also illustrate the limitations of variational calculus when dealing with non-holonomic constraints, like the bounded interval of possible values for the intensity of the thermostat in our problem. Once we are equipped with the necessary mathematical tools, we address the study of the specific physical problem we are interested in. First, we briefly review the dynamics of granular systems and the equations that govern it, in order to be able to apply the mathematical tools presented above. Second, applying the maximum principle, we solve the optimization problem and thus engineer a protocol for the intensity of the thermostat that minimizes the connection time between the two NESS. We analyze the behavior of the system when applying the optimal protocol and solve its dynamics analytically. Moreover, the dependence of these results with the available range for the intensities of the thermostat is investigated. Finally, we explore the limit of applicability of the linear approximation. Resumen El tiempo, a lo largo de la historia, ha marcado un factor determinante en múltiples procesos y ha impuesto límites a batir en otros tantos, desde optimización de tareas en cadenas de montaje, búsqueda de rutas de mínimo tiempo, etc. En esa lucha encarnizada a contrarreloj del hombre contra el tiempo, la matemática ha puesto a nuestro servicio toda una serie de potentes herramientas para abordar, desde ingeniosos puntos de vista, todos los problemas que requieren minimizar el tiempo. En nuestro caso, el problema que tratamos viene de la mano de un sistema granular que, mediante un termostato externo, conducimos desde un estado estacionario de no-equilibrio inicial a otro final en el menor tiempo posible. Para poder abordar esta cuestión es necesario familiarizarse con los métodos matemáticos de optimización que se ajustan mejor al marco teórico impuesto por nuestro problema. Este requerimiento nos lleva a estudiar el principio del máximo, pilar fundamental de la conocida como teoría del control óptimo, y sus variantes lineales, que nos permiten resolver el problema cuando los estados finales e iniciales son próximos entre sí. Aunque no presentamos una demostración rigurosa del principio del máximo, llegamos a él de forma intuitiva empleando las herramientas del cálculo variacional. En este camino mostramos también las limitaciones del cálculo variacional para tratar problemas con ligaduras “no holónomas”, como las presentes en nuestro sistema debido a la acotación en la intensidad del termostato. Una vez equipados con las herramientas matemáticas necesarias, abordaremos el estudio del problema físico en cuestión. Exponemos, sin entrar en demasiado formalismo, la dinámica de los gases granulares y las ecuaciones que la gobiernan, a fin de poder encuadrarlas en el marco matemático expuesto anteriormente. Aplicando el principio del máximo, resolvemos el problema de optimización, diseñando un protocolo para la intensidad del termostato que nos permite minimizar el tiempo de conexión entre dos estados estacionarios y analizamos el comportamien- to del sistema al aplicar dicho protocolo, resolviendo analíticamente su dinámica. Así mismo, se investiga la dependencia de los resultados obtenidos con el rango de valores disponibles para la intensidad del termostato. Por último, exploramos el límite de validez de la aproximación lineal. 10 del cálculo variacional a la teoría del control óptimo: aplicación a sistemas lejos del eqilibrio En esta ocasión 𝜖2es una cantidad positiva de modo que la condición necesaria para que la desigualdad anterior se mantenga es que 𝛿2𝐽𝑦∗(𝜂)𝜖2≥0,∀𝜂∈𝑃𝑎𝑑𝑚,(1.10) donde hemos usado la definición de 𝑜(𝜖2)para afirmar que existirá un valor de 𝜖límite tal que para todo 𝜖con valor absoluto menor que dicho límite el término cuadrático en 𝜖(esto es, la segunda variación), dominará a 𝑜(𝜖2), imponiendo su signo. De este modo, (1.10) conforma la condición necesaria de segundo orden. Estas condiciones generales pueden aplicarse a una gran variedad de problemas. Si fijamos como marco general el espacio normado 𝑉= (𝐶1[𝑎, 𝑏],||⋅||1), esto es, el espacio de funciones 𝑦∶ [𝑎, 𝑏]→ℝ𝑛continuas, con valor ||𝑦||1∶= 𝑚𝑎𝑥[𝑎,𝑏]|𝑦|+𝑚𝑎𝑥[𝑎,𝑏]|𝑦′|bien definido, y como funcional uno de la forma 𝐽(𝑦) = ∫𝑏 𝑎 𝐿(𝑥, 𝑦(𝑥), 𝑦′(𝑥))𝑑𝑥, (1.11) donde 𝐿(𝑥1, 𝑥2, 𝑥3)será una función que tomará valores sobre ℝ×ℝ𝑛×ℝ𝑛,1podemos particularizar las condiciones anteriores a los problemas que presentamos a continuación. Condiciones de primer orden: Ecuaciones de Euler-Lagrange.- Llamamos problema básico del cálculo variacional al problema de minimizar 𝐽en el espacio 𝑉sujeto a las condiciones de contorno 𝑦(𝑎) = 𝑦0y, 𝑦(𝑏) = 𝑦𝑓. Al particularizar la condición necesaria de primer orden obtenemos las ecuaciones de Euler-Lagrange. Para deducirlas bastará reconocer el conjunto 𝐴⊂𝑉 sobre el que deseamos minimizar, que en este caso será 𝐴= {𝑦∈𝐶1[𝑎, 𝑏]|𝑦(𝑎) = 𝑦0, 𝑦(𝑏) = 𝑦𝑓}(1.12) y, por ende, el conjunto de perturbaciones admisibles es 𝑃𝑎𝑑𝑚 = {𝜂∈𝐶1[𝑎, 𝑏]|𝜂(𝑎) = 𝜂(𝑏) = 0}.(1.13) Por otro lado, la forma concreta que adopta la primera variación para un funcional del tipo (1.11) (suponiendo que 𝐿es derivable respecto de 𝑦e𝑦′y llamando 𝑦al mínimo a fin de no recargar la notación) será 𝐽(𝑦+𝜖𝜂) = ∫𝑏 𝑎 𝐿(𝑥, 𝑦(𝑥) + 𝜖𝜂(𝑥), 𝑦′(𝑥) + 𝜖𝜂′(𝑥))𝑑𝑥 =𝐽(𝑦) + ∫𝑏 𝑎(𝜕𝐿 𝜕𝑦 𝜖𝜂 +𝜕𝐿 𝜕𝑦′𝜖𝜂′)𝑑𝑥. Por tanto, 𝛿𝐽𝑦(𝜂) = ∫𝑏 𝑎(𝜕𝐿 𝜕𝑦 𝜂+𝜕𝐿 𝜕𝑦′𝜂′)𝑑𝑥. (1.14) 1Esto es, en general 𝑦tiene 𝑛componentes y nos estamos enfrentando a un problema “multidimensional”. 1. introducción. cálculo de variaciones 11 Basta aplicar integración por partes y nuestra definición de perturbación admisible (con el fin de que desaparezcan los términos de contorno) para obtener, imponiendo que dicha variación sea nula, las conocidas como ecuaciones de Euler-Lagrange ∫𝑏 𝑎[𝜕𝐿 𝜕𝑦 −𝑑 𝑑𝑥 (𝜕𝐿 𝜕𝑦′)]𝜂 𝑑𝑥 = 0.(1.15) Como esto ocurre ∀𝜂∈𝑃𝑎𝑑𝑚 podemos concluir que 𝜕𝐿 𝜕𝑦 −𝑑 𝑑𝑥 (𝜕𝐿 𝜕𝑦′)= 0,(1.16) usando la arbitrariedad de 𝜂como función diferenciable. A aquellas funciones diferenciables que verifiquen las ecuaciones de Euler-Lagrange las denominaremos extremos del funcional J.2 Condiciones de segundo orden: condición de Legendre.- Para obtener esta condición comencemos computando la segunda variación de nuestro funcional, 𝛿2𝐽𝑦(𝜂) = 1 2∫𝑏 𝑎[𝐿𝑦𝑦𝜂2+ 2𝐿𝑦𝑦′𝜂𝜂′+𝐿𝑦′𝑦′𝜂′2]𝑑𝑥, (1.17) donde hemos empleado por primera vez la notación 𝐿𝑦𝑦 para referirnos a la segunda derivada parcial respecto de 𝑦. Usaremos la notación de subíndices para las derivadas cuando esto nos evite sobrecargar innecesariamente la notación. Para realizar este cálculo basta evaluar el funcional, fijada una perturbación arbitraria 𝜂∈𝑉, en la función 𝑦+𝜖𝜂. Esto convierte a 𝐽en una función escalar de 𝜖a la que es legítimo calcularle su desarrollo de Taylor en torno a 𝜖= 0 hasta segundo orden en el parámetro. Si integramos por partes ∫𝑏 𝑎2𝐿𝑦𝑦′𝜂𝜂′𝑑𝑥 , teniendo en cuentas las condiciones de contorno (1.13), con el fin de eliminar el término de la segunda variación que depende del producto de la perturbación y su derivada, nos queda una expresión para 𝛿2𝐽𝑦de la forma 𝛿2𝐽𝑦=∫𝑏 𝑎[𝑃(𝑥)𝜂′2+𝑄(𝑥)𝜂2]𝑑𝑥, (1.18) donde 𝑃(𝑥) = 1 2𝐿𝑦′𝑦′(𝑥, 𝑦(𝑥), 𝑦′(𝑥)), 𝑄(𝑥) = 1 2[𝐿𝑦𝑦(𝑥, 𝑦(𝑥), 𝑦′(𝑥)) − 𝑑 𝑑𝑥𝐿𝑦𝑦′(𝑥, 𝑦(𝑥), 𝑦′(𝑥))]. (1.19) El hecho de que aparezcan derivadas de tercer orden del Lagrangiano nos obligará a exigir unas hipótesis extra de regularidad que pueden verse ampliamente satisfechas si imponemos que 𝐿∈𝐶3. 2Notemos que hay 𝑛ecuaciones de Lagrange en general, ya que 𝑦tiene 𝑛componentes. 12 del cálculo variacional a la teoría del control óptimo: aplicación a sistemas lejos del eqilibrio Figura 1.1: Familia de perturbaciones adecuadas para la deducción de la condición de Legendre. En el panel superior aparece recogida una familia de perturbaciones de valor nulo fuera de un intervalo [𝑐, 𝑑]y valor acotado superiormente por uno en dicho intervalo, mientras que en el panel inferior se presentan las derivadas de dichas perturbaciones, que no se hallan acotadas, pues dependiendo de la pendiente de las funciones 𝜂𝜖, estas podrían disparar su valor. Esta asunción nos da una cierta regularidad sobre 𝑃y𝑄, que deben ser funciones necesariamente continuas. Habiendo calculado ya la expresión de la segunda derivada podemos aplicar la condición general necesaria de segundo orden. Para ello consideremos una función 𝑦∗∈ (𝐶1([𝑎, 𝑏]),||.||1) mínimo local, debe entonces cumplirse que 𝛿2𝐽∗ 𝑦=∫𝑏 𝑎[𝑃(𝑥)𝜂′2+𝑄(𝑥)𝜂2]𝑑𝑥 ≥0.(1.20) La idea reside, aunque no nos detendremos en la construcción pues no es objeto de este trabajo (para una deducción más rigurosa consultar [2, 3]), en diseñar una familia de perturbaciones admisibles dependientes de un parámetro, pongamos 𝜖, que deben cumplir (al ser admisibles) la condición necesaria de segundo orden 𝛿2𝐽𝑦∗(𝜂𝜖)≥0,(1.21) cuyo rango se encuentre acotado entre 0y1pero cuya derivada crezca al disminuir 𝜖, como muestra la figura 1.1. Este tipo de perturbación nos sugiere que el término dominante y que por tanto impondrá su signo dentro de la integral será el termino 𝑃(𝑥)𝜂′ 𝜖 2que dominará frente a 𝑄(𝑥)𝜂𝜖, pues mientras que para valores de 𝜖cada vez mas pequeños el valor de 𝜂𝜖permanece 1. introducción. cálculo de variaciones 13 acotado por la unidad, su derivada se dispara. De este modo acabamos de deducir la conocida como condición necesaria de Legendre: Si 𝑦∗(𝑥)es un mínimo, entonces debe cumplirse que 𝐿𝑦′𝑦′(𝑥, 𝑦∗(𝑥), 𝑦∗(𝑥)′)≥0.(1.22) Para entender la información que nos ofrece esta condición debemos definir lo que llamaremos el hamiltoniano de una función, así como el momento: Definición 1.5 (Hamiltoniano). Definimos la función hamiltoniano, 𝐻, como sigue: 𝐻∶ℝ×ℝ𝑛×ℝ𝑛×ℝ𝑛→ℝ, 𝐻(𝑥1, 𝑥2, 𝑥3, 𝑥4) = 𝑥4𝑥3−𝐿(𝑥1, 𝑥2, 𝑥3),con 𝑛≥1.(1.23) Definición 1.6 (Momento). Definimos la función momento, 𝑝, como sigue: 𝑝∶ℝ×ℝ𝑛×ℝ𝑛→ℝ𝑛𝑝=𝐿𝑥3(𝑥1, 𝑥2, 𝑥3),con 𝑛≥1.(1.24) Evaluando las definiciones anteriores para el caso en que las variables describen una trayectoria parametrizada en 𝑥se tiene: 𝑝(𝑥, 𝑦, 𝑦′) = 𝐿𝑦′(𝑥, 𝑦, 𝑦′), 𝐻(𝑥, 𝑦(𝑥), 𝑦′(𝑥), 𝑝(𝑥, 𝑦, 𝑦′)) = 𝑝(𝑥, 𝑦, 𝑦′)𝑦′−𝐿(𝑥, 𝑦, 𝑦′).(1.25) Si consideramos (𝑥, 𝑦, 𝑦′, 𝑝)como variables independientes (esto sería, antes de ser evaluadas sobre la trayectoria) y fijamos (𝑥, 𝑦, 𝑝)a unos valores determinados, es legítimo estudiar el comportamiento de 𝐻respecto de la variable 𝑦′(como variable independiente respecto de las otras tres). Lo que nos dice la condición de Legendre es que si 𝑦∗es mínimo local, para cualquier 𝑥fijo del intervalo [𝑎, 𝑏]se debe verificar que si (𝑦∗(𝑥))′es un extremo del hamiltoniano, entonces debe necesariamente ser un máximo, puesto que 𝜕𝐻(𝑥, 𝑦, 𝑦′, 𝑝) 𝜕𝑦′|𝑥,𝑦,𝑝 =𝑝−𝜕𝐿(𝑥, 𝑦, 𝑦′) 𝜕𝑦′|𝑥,𝑦,𝑝 =𝑝−𝐿𝑦′= 0, 𝜕2𝐻(𝑥, 𝑦, 𝑦′, 𝑝) 𝜕2𝑦′|𝑥,𝑦,𝑝 = −𝐿𝑦′𝑦′(𝑥, 𝑦, 𝑦′). (1.26) La primera ecuación es nula por la propia definición de momento, y la segunda ecuación es menor o igual a 0 al ser evaluada sobre 𝑦∗, por la condición de Legendre. 2 Teoría del control óptimo En esta sección reformularemos el problema de minimización de una forma más conveniente de cara a nuestros futuros cálculos, estableciendo así el marco general del conocido como problema de control óptimo. 2.1 Problema de control óptimo La idea de partida en la resolución de cualquier problema en cálculo de variaciones era fijar un funcional a minimizar, definir un conjunto de funciones con una cierta regularidad y buscar, de entre todas ellas, aquella que generase el valor mínimo. Sin embargo, podríamos cambiar este enfoque estático, en el que simplemente buscamos entre la colección de curvas permitidas, por uno en el que vamos progresivamente construyendo la curva en base de pequeñas decisiones óptimas que tomamos en cada instante. Puede entenderse de la siguiente forma: supogamos que partimos de un cierto punto 𝑥0en un instante de tiempo 𝑡0y deseamos mínimizar un cierto funcional 𝐽. La pregunta que surge de manera natural dado un punto de partida es: ¿cuál es la dirección de avance óptima? Es decir, ¿qué 𝑥′(de aquí en adelante notaremos 𝑥′como la derivada temporal de 𝑥) debería escoger en el tiempo 𝑡0para avanzar en la dirección de minimización? Escogida esta velocidad nos hallamos en un punto 𝑥0+𝑥′(𝑡0)𝑑𝑡 en un tiempo 𝑡0+𝑑𝑡. Se plantea por tanto el mismo dilema, ¿ qué rumbo fijar ahora desde las nuevas coordenadas?. Repitiendo este proceso de forma iterativa, acabaríamos contruyendo una trayectoria óptima basada en la toma de decisiones infinitesimalmente óptimas. Esta forma dinámica de abordar el problema nos va a permitir acomodar restricciones de manera muy sencilla e introducir de forma natural un sistema controlado mediante un parámetro, que en ocasiones modeliza de manera más natural ciertos problemas que se nos pueden presentar. Es precisamente esta forma de proceder la que evidencia el uso de la variable temporal como variable natural independiente para parametrizar nuestra trayectoria, que notaremos, en lugar de 𝑦(𝑥)como hicimos anteriormente, como 𝑥(𝑡). Antes de plantear más formalmente el problema de control óptimo resulta interesante hacer un breve comentario: como ya hemos notado antes, nuestras decisiones se basarán en fijar una 16 del cálculo variacional a la teoría del control óptimo: aplicación a sistemas lejos del eqilibrio determinada velocidad o, de forma más general, de fijar un valor del parámetro externo, 𝑢, que controle mediante alguna función propia del sistema considerado a la velocidad en cuestión. La clave en todo esto reside, por tanto, en determinar qué valor de 𝑢tomamos una vez alcanzado un estado (𝑡, 𝑥(𝑡), 𝑥′(𝑡)) para que la trayectoria que tracemos durante el siguiente intervalo infinitesimal sea la adecuada. Vamos a ver que la herramienta matemática que nos acotará el número de valores candidatos para 𝑢en ese estado es precisamente el hamiltoniano. El cómo genera esta información el hamiltoniano puede resultarnos familiar si recordamos la condición de Legendre y suponemos un sistema de control sencillo donde el parámetro de control sea directamente la velocidad 𝑥′=𝑢:si 𝑢∗(𝑡)es la decisión óptima a tomar en el instante t, esto es, si x*(t) es óptima, 𝑢∗(𝑡)será un máximo del Hamiltoniano. Ahora sí, formalicemos el marco general de un problema de control óptimo. Existen tres partes que definen el problema: Sistema de control óptimo.- Nos indicará la relación entre el parámetro de control y la velocidad de la trayectoria. 𝑥′=𝑓(𝑡, 𝑥, 𝑢), 𝑥(𝑡0) = 𝑥0.(2.1) Evidentemente desearemos que, una vez escogido un control 𝑢(𝑡)(ya veremos qué condiciones de regularidad le impondremos), el sistema presente una solución única. Hay pues, tres cuestiones principales que debemos atender llegados a este punto: qué regularidad imponemos sobre 𝑢como función del tiempo, qué hipótesis de regularidad imponemos sobre 𝑓para que el sistema tenga solución y, más importante, qué entenderemos como solución del sistema. Con respecto a la regularidad de 𝑢(𝑡)trabajaremos en el contexto de las funciones medibles en 𝑡, aunque en el problema que a nosotros nos compete nos bastará exigir que sea continua a trozos, y esa constituirá nuestra hipótesis de regularidad a la hora de deducir el principio del máximo. La cuestión sobre la existencia de solución del sistema quedará ampliamente satisfecha si 𝑓(𝑡, 𝑥, 𝑢)es continua en 𝑡,𝐶1en la variable 𝑥y continua en 𝑢. Bajo estas suposiciones podemos garantizar la existencia de una única función 𝑥(𝑡)absolutamente continua,1es decir, definida sobre un cierto intervalo [𝑡0, 𝑡𝑓], que es continua en todo tiempo 𝑡y diferenciable para casi todo 𝑡, y que cumple la ecuación integral 𝑥(𝑡) = 𝑥0+∫𝑡 𝑡0 𝑓(𝑠, 𝑥(𝑠), 𝑢(𝑠))𝑑𝑠. (2.2) Esta función 𝑥(𝑡), cuya existencia y unicidad queda garantizada en el marco anterior, es la solución de nuestro sistema diferencial. Como puede comprobarse, en el planteamiento del problema de control óptimo no se exige ningún tipo de regularidad sobre la trayectoria, sino que su regularidad surge de forma natural al introducir el parámetro de control y resolver la ecuación. De esta forma, 1Intuitivamente, diferenciable salvo en un conjunto de puntos de medida nula. 2. teoría del control óptimo 17 si exigimos que 𝑢(𝑡)sea continuo a trozos, la mínima regularidad que obtenemos para 𝑥(𝑡)es ser absolutamente continua, lo cual es inherente a su condición de solución del sistema. De esta forma, desaparece toda posible discusión sobre espacios normados donde plantear el problema: habrá un único espacio, el de los controles continuos a trozos (o si queremos generalizar, el de los controles medibles) con imagen un conjunto cualquiera, 𝑈, sobre el que definiremos un problema de mínimización global que nos liberará del “yugo” de la norma. Coste.- El funcional 𝐽con el que trabajaremos conforma el conocido como problema de Bolza, 𝐽(𝑢;𝑥0, 𝑡0, 𝑥𝑓, 𝑡𝑓) = ∫𝑡𝑓 𝑡0 𝐿(𝑡, 𝑥(𝑡), 𝑢(𝑡))𝑑𝑡 +𝐾(𝑥𝑓, 𝑡𝑓).(2.3) Destaquemos dos aspectos referentes a la ecuación anterior. Primero, la similitud entre el término integral y el funcional tratado en el contexto del cálculo variacional es evidente, basta recordar que 𝑥′(𝑡) = 𝑓(𝑡, 𝑥(𝑡), 𝑢(𝑡)) y de este modo, el lagrangiano que aparece en calculo de variaciones quedaría como: 𝐿(𝑡, 𝑥(𝑡), 𝑥′(𝑡)) = 𝐿(𝑡, 𝑥, 𝑓(𝑥, 𝑥(𝑡), 𝑢(𝑡)) =  𝐿(𝑡, 𝑥(𝑡), 𝑢(𝑡)).(2.4) Nótese que el camino inverso de paso de un lagrangiano cualquiera con dependencia en 𝑢(𝑡)a otro con dependencia en 𝑥′(𝑡)no tiene porque ser siempre posible. Segundo, aparece el término 𝐾que se conoce como coste final. Un problema de Bolza donde solo tenemos el término integral se conoce como problema de Lagrange, mientras que el problema que solo presente un coste final se como como problema de Mayer. Vamos a poder pasar de una formulación a otra sin problema, puesto que ambos planteamientos resultan equivalentes. Conjunto objetivo.- Este último ingrediente es de vital importancia, pues nos permitirá acomodar una gran variedad de problemas en el contexto del control óptimo. El conjunto objetivo o conjunto blanco (target set en inglés) nos brinda una forma muy cómoda de recoger las condiciones de contorno del problema, englobándolas en un conjunto 𝑆 ⊂ [𝑡0,∞) × 𝑅𝑛. De este modo un problema de estado y tiempos finales fijos presenta un conjunto objetivo 𝑆=𝑡𝑓×𝑥𝑓, en un problema de tiempo final fijo y extremo libre tenemos 𝑆=𝑡𝑓×ℝ𝑛,y un problema donde a tiempo final fijo el estado final deba reposar sobre un determinado subconjunto 𝑀 ⊂ ℝ𝑛tendríamos 𝑆=𝑡𝑓×𝑀. También podemos imaginar una situación donde 𝑡𝑓no sea fijo, de modo que se necesitará una definición concreta para este tiempo final. Para ello bastará exigir que la parte espacial de S sea un conjunto cerrado, de modo que definiremos 𝑡𝑓como el menor de los tiempos tal 𝑥(𝑡)interseca dicha parte espacial (que está bien definido por ser S cerrado).2Notemos que esta definición de 2Podría ocurrir que si el conjunto es abierto el tiempo 𝑡𝑓no estuviese bien definido pues la trayectoria podría acercarse indefinidamente a 𝑆sin llegar a tocarlo, por ejemplo 𝑥(𝑡) = √𝑡y𝑆= (−∞,0). 18 del cálculo variacional a la teoría del control óptimo: aplicación a sistemas lejos del eqilibrio tiempo final presenta como caso particular aquel en el que 𝑡𝑓es fijo. Habiendo tratado los tres ingredientes principales punto por punto estamos en condiciones de formular el problema de control óptimo completo: Definición 2.1 (Problema de control óptimo). Buscamos 𝑢:ℝ→𝑈 ⊂ ℝ𝑚una función continua a trozos (o, más generalmente, medible) con rango contenido en U tal que, al ser introducido en el sistema 𝑥′=𝑓(𝑡, 𝑥, 𝑢), 𝑥(𝑡0) = 𝑡0,(2.5) genere una trayectoria 𝑥(𝑡)que, por un lado, interseque al conjunto objetivo y por otro minimice, conjuntamente con 𝑢(𝑡), el funcional 𝐽. Como puede comprobarse, la función clave que minimiza el funcional es el control y no la trayectoria 𝑥(𝑡), pues esta última queda unívocamente determinada una vez se ha fijado el primero. 2.2 Tratamiento variacional del problema de control óptimo Para introducir el principio del máximo, que se trata de un teorema que ofrece unas condiciones necesarias muy útiles para la resolución de una gran cantidad de problemas de control óptimo, desde la óptica del cálculo de variaciones hemos de restringirnos a un problema de control que encaje dentro de esta teoría. Consideremos pues el siguiente problema: Tomemos S como 𝑆= {𝑡𝑓} × ℝ𝑛, esto es, un problema de tiempo final fijo y estado final libre. Tomemos el rango de 𝑈 como ℝ𝑚, esto es, tenemos un problema sin restricciones sobre 𝑢(no hay restricción en su imagen ni tampoco sobre el conjunto blanco que debe golpear la trayectoria que este genere). Consideremos un funcional de la forma 𝐽(𝑢) = ∫𝑡𝑓 𝑡0 𝐿(𝑡, 𝑥(𝑡), 𝑢(𝑡))𝑑𝑡 +𝐾(𝑥(𝑡𝑓)),(2.6) donde el término de coste final, 𝐾(𝑥𝑓(𝑡𝑓)) no depende explícitamente del tiempo 𝑡𝑓. Nuestro objetivo es deducir unas condiciones necesarias de mínimo global para 𝑢(esto es, el control óptimo 𝑢∗debe cumplir 𝐽(𝑢∗)≤𝐽(𝑢)para cualquier función 𝑢que sea continua a trozos), empleando las herramientas que hasta ahora hemos descrito. Es claro que si 𝑢∗es un mínimo en sentido global lo será en sentido local si dotamos al espacio de las funciones continuas a trozos (espacio de controles aceptables) de una norma cualquiera, por ejemplo, la norma 0 (esto es ||𝑥||0∶= 𝑚𝑎𝑥[𝑎,𝑏]|𝑥|). Esto nos permite llevarnos el problema a un terreno conocido donde podemos plantearnos extraer condiciones necesarias de primer y segundo orden de la mera particularización de las condiciones necesarias ya deducidas. Esto es precisamente lo que vamos a hacer. 2. teoría del control óptimo 19 Nuestro esquema de razonamiento será el siguiente: consideraremos 𝑢∗el óptimo global del problema, perturbaremos dicho 𝑢∗y estudiaremos como esa perturbación altera la trayectoria óptima (a primer orden en un cierto parámetro). Reescribiremos el funcional incorporando el sistema de control en forma de restricción en el integrando, multiplicada por unas funciones 𝑝(𝑡) cualesquiera (que pueden dar una falsa impresión de ser los multiplicadores de Lagrange, pero su verdadera naturaleza será esclarecida más adelante). Sustituyendo el control y trayectorias perturbadas, calcularemos la primera variación e impondremos que esta se anule. Haciendo uso de la libertad con la que los 𝑝(𝑡)fueron escogidos, seremos capaces de deducir un sistema diferencial de primer orden y dos ecuaciones que involucren a unos 𝑝∗(𝑡)muy particulares y a la trayectoria óptima 𝑥∗(𝑡), junto con unas condiciones iniciales que derivaremos del término de coste final. Las ecuaciones diferenciales obtenidas se corresponderán con las ecuaciones canónicas ampliamente conocidas en el marco de la mecánica teórica. Posteriormente, el estudio del término dominante en la segunda variación nos permitirá entender la dependencia del hamiltoniano con el control, una vez este primero ha sido evaluado sobre la trayectoria óptima. Perturbación del control y la trayectoria.- Comenzamos perturbando 𝑢en lugar de 𝑥pues ahora minimizamos el funcional respecto de 𝑢y no respecto de la trayectoria. De este modo, definimos el conjunto de perturbaciones admisibles para el control como 𝑃𝑎𝑑𝑚(𝑢)= {𝜒continuas a trozos sobre [𝑡0, 𝑡𝑓]}.(2.7) Así los controles que se obtienen aplicando este tipo de perturbaciones serán de la forma: 𝑢=𝑢∗+𝜖𝜒, (2.8) con 𝜖un parámetro real arbitrario. Podemos, fijada una perturbación cualquiera 𝜒, definir una familia de trayectorias parametrizadas en 𝜖,𝑥(𝑡, 𝜖), resultado de resolver el sistema de control aplicando los distintos controles. Dicha familia cumple que 𝑥(𝑡, 𝜖 = 0) = 𝑥∗(𝑡). Si tomamos 𝜖suficientemente pequeño y asumimos que la dependencia de 𝑥(𝑡, 𝜖)con dicho parámetro es diferenciable, podemos hacer un desarrollo de Taylor a orden uno en 𝜖, 𝑥(𝑡, 𝜖) = 𝑥∗(𝑡) + 𝑥𝜖(𝑡, 0)𝜖+𝑜(𝜖).(2.9) Por ende, llamando 𝜂a la perturbación de la trayectoria óptima a primer orden en 𝜖, 𝜂=𝑥𝜖(𝑡, 0),(2.10) podemos plantearnos encontrar la relación entre la perturbación del control 𝜒y la de la trayectoria 𝜂, a 𝑜(𝜖). Para ello disponemos del sistema de ecuaciones diferenciales 𝑥′(𝑡) = 𝑓(𝑡, 𝑥(𝑡), 𝑢(𝑡)), 𝑥(𝑡0) = 𝑥0,(2.11) 26 del cálculo variacional a la teoría del control óptimo: aplicación a sistemas lejos del eqilibrio tiempo 𝑡0= 0 en una posición 𝑥10 y con una velocidad 𝑥20. Supongamos también que deseamos aparcar este coche en el origen en el menor tiempo posible, teniendo en cuenta que la capacidad de aceleración de nuestro coche no es infinita y que se encuentra acotada en el intervalo 𝑈= [−1,1]. Si deseásemos abordar el problema empleando las herramientas del cálculo variacional nos encontramos ante distintas dificultades. Notemos por 𝑥1(𝑡)y𝑥2(𝑡)la posición y velocidad del coche en el instante de tiempo 𝑡, respectivamente. Para minimizar el tiempo, podríamos definir el funcional como 𝐽(𝑥1, 𝑥2, 𝑢) = ∫𝑡∗ 𝑡01𝑑𝑡, pero no aplicar ninguna de las dos condiciones (Euler-Lagrange y Legendre) obtenidas anteriormente, pues 𝛿𝐽 y𝛿𝐽2son nulos. Esta dificultad es salvable introduciendo las ecuaciones del movimiento como “restricciones holónomas” con una extensión del método de los multiplicadores de Lagrange.4Pero tenemos otra dificultad que no es salvable: no podemos incluir las restricciones de desigualdad, “no holónomas”, para la aceleración, −1 ≤𝑥′′ ≤1. En el marco de la teoría del control óptimo podemos plantear el problema anterior de manera natural. La posición y velocidad del coche en el instante de tiempo 𝑡evolucionan en el tiempo de acuerdo con el sistema 𝑥′ 1=𝑥2, 𝑥′ 2=𝑢, (2.36) donde el control presenta un rango 𝑢∈ [−1,1]. Se desea llevar al sistema desde un punto (𝑥1(𝑡0), 𝑥2(𝑡0)) hasta el punto (𝑥1(𝑡∗)=0, 𝑥2(𝑡∗) = 0) en el menor tiempo posible (donde 𝑡∗se define, como ya se comento en la sección 3.1, como el menor de los tiempos en los que la trayectoria interseca el punto (𝑥1= 0, 𝑥2= 0) del espacio de las fases). En vista del sistema de control, nos encontramos ante un problema lineal de la forma 𝑥′=𝐴𝑥 +𝑏𝑢, con 𝐴=(0 1 0 0), 𝑏 =(0 1),(2.37) que verifica las condiciones para la aplicación del teorema 2.2 puesto que el autovalor de A es 0, el conjunto {𝑏, 𝐴𝑏}conforma una base de ℝ2y el rango del control es un intervalo de ℝ. Por tanto, sabemos que todos los controles óptimos, en caso de existir, adoptarán los valores 1y−1 con, a lo sumo, un cambio entre estos dos valores. Aunque no nos extenderemos en la discusión de aspectos sobre existencia y unicidad (para más información consultar [4]), sí que nos gustaría indicar que la existencia de un control cualquiera que conecte los estados inicial y final garantiza la existencia del control 𝑢∗óptimo. En este caso concreto la existencia de este control queda asegurada por el teorema de Filippov [2]. De este modo, la forma más genérica de estos controles será, dependiendo de las condiciones iniciales (y pudiendo 𝑡𝑗ser nulo) 4La idea intuitiva es que el tiempo se hace mínimo sobre el conjunto de soluciones de las ecuaciones del movimiento. Una discusión general del uso de multiplicadores de Lagrange en el contexto del cálculo variacional se puede encontrar en la Ref. [1], y una aplicación a un problema físico para la minimización del tiempo de conexión en la Ref. [5]. 2. teoría del control óptimo 27 𝑢∗(𝑡) = {−1, 𝑠𝑖 𝑡 ≤𝑡𝑗, 1, 𝑠𝑖 𝑡 > 𝑡𝑗,(2.38) 𝑢∗(𝑡) = {1, 𝑠𝑖 𝑡 ≤𝑡𝑗, −1, 𝑠𝑖 𝑡 > 𝑡𝑗,(2.39) anulándose ambos en 𝑡∗para que el coche quede inmóvil en el origen. Establecido lo anterior, el único parámetro que debemos determinar para poder definir de manera unívoca el control óptimo es el tiempo de cambio 𝑡𝑗. Para ello podemos aplicar la condición de máximo sobre el hamiltoniano, que en este problema adopta la forma 𝐻(𝑥0, 𝑥1.𝑥2, 𝑝0, 𝑝1, 𝑝2, 𝑢) = 𝑝0+𝑝1𝑥2+𝑝2𝑢. Claramente el máximo de la función 𝑢→𝐻(𝑢)se alcanza en los extremos del intervalo 𝑈, dependiendo del signo de 𝑝2. Por tanto la expresión de 𝑢será 𝑢∗(𝑡) = {1, 𝑠𝑖 sgn(𝑝∗ 2(𝑡))>0, −1, 𝑠𝑖 sgn(𝑝∗ 2(𝑡))<0,(2.40) Podemos hallar la expresión de 𝑝2(𝑡)empleando las ecuaciones canónicas 𝑝∗ 1 ′= −𝐻𝑥1= 0, 𝑝∗ 2 ′= −𝐻𝑥2=𝑝∗ 1.(2.41) Por tanto, 𝑝∗ 1(𝑡) = 𝑐1y𝑝∗ 2(𝑡) = 𝑐1𝑡+𝑐2y el tiempo de cambio será 𝑡𝑗= −𝑐2 𝑐1 . Del mismo modo podemos hallar las trayectorias (𝑥1(𝑡), 𝑥2(𝑡)), teniéndose 𝑥2(𝑡) = 𝑥20 +𝑢𝑡 y𝑥1(𝑡) = 𝑥10 +𝑥20𝑡+1 2𝑢𝑡2. Por otro lado, sabemos que 𝐻∗≡0para todo tiempo, particularmente para 𝑡𝑗 𝐻∗(𝑡𝑗= −𝑐2 𝑐1)= −1 + 𝑐1(𝑥10 +𝑥20𝑡𝑗+1 2𝑢 𝑡𝑗2)= 0,(2.42) donde hemos tomado 𝑝0= −1 pues 𝐻está indeterminado en una constante multiplicativa5y donde 𝑢= ±1 dependiendo de las condiciones iniciales. Llegados a este punto podemos hacer uso de la intuición física para saber qué valor comenzará adoptando 𝑢durante la primera ventana temporal: para 𝑥10 <0se tiene que 𝑢= +1 durante el primer intervalo de tiempo y para 𝑥10 >0se tendrá que 𝑢= −1. Por tanto tenemos dos ecuaciones a resolver según la condición inicial para 𝑥1, una en la que se fija 𝑢= 1 en (2.42) y otra en la que se fija 𝑢= −1. Sin embargo, en cada caso tenemos dos incógnitas, 𝑐1y𝑐2, por lo que necesitaremos una segunda ecuación para poder despejar la cual puede obtenerse imponiendo la condición 𝐻∗≡0en el tiempo inicial 𝑡0= 0, 𝐻∗(𝑡= 0) = −1 + 𝑐1𝑥20 +𝑐2𝑢= 0,(2.43) donde 𝑢adquiere nuevamente los valores discutidos en el párrafo anterior. De este modo los sistemas a resolver son, según las condiciones iniciales: 5Basta comprobar que, si multiplicamos 𝐻y𝑝por la misma constante, las ecuaciones canónicas quedan invariantes, por lo que existe libertad para fijar la componente 𝑝0del vector de momentos. 28 del cálculo variacional a la teoría del control óptimo: aplicación a sistemas lejos del eqilibrio Para 𝑥10 <0 −1+𝑐1[𝑥10 +𝑥20 (−𝑐2 𝑐1)+1 2(−𝑐2 𝑐1)2]= 0, −1+𝑐1𝑥20 +𝑐2= 0, (2.44) Para 𝑥10 >0 −1+𝑐1[𝑥10 +𝑥20 (−𝑐2 𝑐1)−1 2(−𝑐2 𝑐1)2]= 0, −1+𝑐1𝑥20 −𝑐2= 0. (2.45) Podemos despejar 𝑐1y𝑐2y el tiempo 𝑡𝑗en función de las condiciones iniciales: Para 𝑥10 <0. Si se cumple que 3𝑥2 20 + 2𝑥10 ≠0, entonces se tiene 𝑐1=(1+2𝑥20 ±√1−4𝑥20 +𝑥2 20 + 2𝑥10) 3𝑥2 20 − 2𝑥10 , 𝑐2= 1 − 𝑐1𝑥20.(2.46) y por tanto 𝑡𝑗puede escribirse como el menor positivo del siguiente conjunto de dos elementos, 𝑡𝑗=min𝑡𝑗>0⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝑥20 −3𝑥2 20 − 2𝑥10 (1+2𝑥20 ±√1−4𝑥20 +𝑥2 20 + 2𝑥10)⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ .(2.47) Si, por lo contrario, se verifica la relación 3𝑥2 20 + 2𝑥10 = 0, se tiene que no hay solución. En consecuencia, la curva descrita por 3𝑥2 20 + 2𝑥10 = 0 en el espacio de las fases conforma la curva de cambio para 𝑥10 <0. Para 𝑥10 >0. Si se cumple que 3𝑥2 20 − 2𝑥10 ≠0, entonces se tiene 𝑐1= ±(1−2𝑥20 +√1−4𝑥20 +𝑥2 20 + 2𝑥10) 3𝑥2 20 − 2𝑥10 , 𝑐2= −1 + 𝑐1𝑥20.(2.48) Por tanto, 𝑡𝑗puede escribirse como el menor positivo del siguiente conjunto de dos elementos, 𝑡𝑗=min𝑡𝑗>0⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ −𝑥20 ±3𝑥2 20 − 2𝑥10 (1−2𝑥20 +√1−4𝑥20 +𝑥2 20 + 2𝑥10)⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ .(2.49) 2. teoría del control óptimo 29 Figura 2.1: Trayectorias en el espacio de las fases. Curva de cambio (linea roja punteada) y trayectorias (línea sólida azul y roja) para dos condiciones iniciales diferentes, (𝑥10, 𝑥20)y(𝑥10, 𝑥20). Las flechas indican el sentido en que se recorren las curvas. Si, por lo contrario, se verifica la relación 3𝑥2 20 − 2𝑥10 = 0 se tiene que no hay solución: la curva descrita por 3𝑥2 20 − 2𝑥10 = 0 en el espacio de las fases conforma la curva de cambio para 𝑥10 >0. Podemos entender de manera más gráfica lo que ocurre con las trayectorias pasando al espacio de las fases (𝑥1, 𝑥2), donde se verifica la relación 𝑥1= ±1 2𝑥22+𝑐para 𝑢= ±1 respectivamente y 𝑐es una contante función de las condiciones de contorno. Sabemos que bajo la acción, o bien de 𝑢= +1 o de 𝑢= −1, el sistema debe llegar al punto (0,0). Bastará por tanto ajustar la constante 𝑐 para obtener las curvas en el espacio de las fases que para ambos valores del control intersequen el (0,0). La unión de estas dos curvas (esto es, 3𝑥2 20 +2𝑥10 = 0 si 𝑥10 <0y3𝑥2 20 −2𝑥10 = 0 si 𝑥10 >0) conforma la curva de cambio. Además sabemos que cuando 𝑥10 >0debe actuar primero 𝑢= −1, por ende la trayectoria fásica será en este caso 𝑥1= −1 2𝑥22+𝑐para la c correspondiente hasta alcanzar la curva de cambio, donde se invierte el valor del control a 𝑢= +1 y el sistema se mueve bajo la acción del nuevo control hasta el origen. Análogamente, si 𝑥10 <0el sistema comienza moviéndose bajo el efecto de 𝑢= +1 describiendo la parábola 𝑥1=1 2𝑥22+𝑐hasta intersecar la curva de cambio, a partir de ese momento 𝑢= −1 y con este valor del control alcanza el origen, tal y como se muestra en la figura 2.1. Nos gustaría destacar que una vez alcanzado el punto (0,0) el control deberá cambiar instantáneamente al valor 𝑢(𝑡) = 0 para 𝑡≥𝑡∗, a fin de mantener el coche en el origen (velocidad nula) para tiempos superiores a 𝑡∗. Este salto del control no afecta al carácter óptimo del mismo. Parte II Aplicación a sistemas lejos del equilibrio 3 Minimización del tiempo de conexión 3.1 Medios granulares El estudio de medios granulares, más allá del propio interés teórico que suscita, es de vital importancia en el ámbito de la industria y en la mejora del transporte y almacenamiento de multitud de productos, desde agrícolas hasta farmacéuticos. Los materiales granulares son conglomerados discretos de partículas macroscópicas que exhiben dos características fundamentales: La temperatura ordinaria no desempeña ningún papel relevante: a diferencia de los gases ordinarios, donde el término 𝑘𝐵𝑇establece una escala energética, en los medios granulares esta cantidad resulta irrelevante frente a, por ejemplo, la energía potencial gravitatoria, que es típicamente 1012 veces superior en magnitud a la térmica. En vista de este hecho es posible definir una temperatura granular que mide de las fluctuaciones de la velocidad respecto de su media, permitiendo así recuperar ciertas herramientas “termodinámicas” para abordar el estudio de estos medios. La interacción entre los granos es disipativa: a causa de las colisiones inelásticas y la fricción estática. La no conservación de la energía tiene importantes consecuencias conceptuales, ya que su conservación en los sistema moleculares permite la separación de escalas de tiempo necesaria para, por ejemplo, establecer una descripción hidrodinámica. Dependiendo del estado de movimiento de los granos vamos a poder diferenciar entre estado sólido, líquido o gaseoso de un medio granular. Examinemos a continuación algunas de sus características más peculiares, que hacen de los medios granulares un interesante campo de estudio. Sólidos granulares.- Consideremos un medio granular en reposo (por ejemplo, arena) contenido en un cilindro. Puede comprobarse como en ese estado la presión en la base del cilindro se vuelve, a partir de una cierta altura, independiente de esta (a diferencia de lo que ocurre con los líquidos convencionales) [6]. Este comportamiento particular implica que el exceso de presión que no recae en la base sea soportado por las paredes del recipiente. Otras características propias de los sólidos granulares son, por una parte, la gran sensibilidad que presentan las redes de fuerzas generadas en el medio a la disposición de los granos del mismo, cuando este es sometido a alguna 34 del cálculo variacional a la teoría del control óptimo: aplicación a sistemas lejos del eqilibrio tensión. Por otra parte, estos medios presentan una gran aleatoriedad en el empaquetamiento: cambia dependiendo del procedimiento seguido a la hora de verter los granos en el cilindro. Se ha comprobado que la fracción de empaquetamiento, 𝜂, varía entre dos valores extremos: 𝜂= 0.55 y 𝜂= 0.64, pudiendo pasar de una fracción a otra a través de vibraciones externas que permitan al sistema explorar el “espacio de las fases” [7]. Líquidos granulares.- Los medios granulares pueden fluir, lo que abre una puerta a la creación de modelos continuos en el seno de lo que se conoce como hidrodinámica granular. Estas teorías deben enfrentar el problema de la escala espacial y temporal adecuada para el tratamiento de los flujos granulares, que no se haya tan bien diferenciada de la escala macroscópica y microscópica, como sí ocurre con los fluidos “ordinarios” o moleculares. Del mismo modo, debe lidiar con el hecho de que los flujos lentos y densos son no ergódicos y los flujos rápidos experimentan fenómenos de aglutinamiento, que deben tener cabida en esta teoría [6]. Así mismo, debe ser capaz de explicar la transición entre ambos estados (flujos rápidos y lentos), pues en sistemas reales, debido a la disipación de energía, ambos pueden convivir de manera simultánea. Los fluidos granulares muestran comportamientos de difícil explicación y que son objeto de debate e investigación: la causa de los procesos convectivos y de aglutinamiento inducidos por vibraciones externas y el papel que juega en estos procesos el gas intersticial, así como la causa de la separación por tamaño de grano, lo que se conoce como fenómeno de segregación o “Brazil nut problem” [8, 9]. Gases granulares.- La principal diferencia entre gases granulares y gases convencionales reside en el carácter no conservativo de las interacciones grano-grano. En las colisiones grano-grano parte de la energía se “pierde”, almacenándose en grados internos de libertad o disipándose al entorno. Esto impide, en principio, aplicar teorías fuertemente asentadas en los gases convencionales, como por ejemplo la teoría de gases ideales. Hay grandes interrogantes abiertos en el campo de los gases granulares, como por ejemplo el papel de los grados de libertad de rotación durante una colisión o la dependencia respecto a la velocidad de los granos de las fuerzas de fricción. El modelo que emplearemos para describir el gas granular con el que vamos a tratar en las siguientes secciones está basado en un parámetro, 𝛼, conocido como coeficiente de restitución. Los granos se modelan como esferas duras lisas: en una colisión binaria, la velocidad relativa tangencial permanece inalterada mientras que el signo de la velocidad relativa normal se invierte y su módulo se reduce por un factor 𝛼.1La pérdida de energía en una colisión es entonces proporcional a1 − 𝛼2,𝛼= 1 corresponde al límite elástico y 𝛼= 0 al límite completamente inelástico. Dicho coeficiente será considerado como una constante propia del gas en nuestro modelo, aunque esta simplificación resulta en ocasiones insuficiente para explicar ciertos fenómenos que presentan estos gases [10]. 1Aquí tangencial (normal) se refiere a perpendicular (paralelo) a la recta que une los centros de las dos esferas en la colisión. 3. minimización del tiempo de conexión 35 3.2 Descripción del sistema Consideremos un gas granular de densidad 𝑛, coeficiente de restitución 𝛼, bi o tridimensional, conformado por discos o esferas duras respectivamente, de masa 𝑚y diámetro 𝜎. Para un gas de este tipo la función de distribución de velocidades de una partícula no se corresponde con la de Maxwell, debido a la acción de las colisiones inelásticas. Sin embargo, es posible presentar una solución aproximada en términos de los denominados polinomios de Sonine, como veremos a continuación [10]. A la hora de estudiar el gas granular consideraremos dos situaciones. Una primera, a modo ilustrativo, donde el sistema estará libre de fuerza alguna y la evolución de la función de distribución de velocidades monoparticular 𝑓(𝑟, 𝑣, 𝑡)vendrá descrita por la ecuación de Enskog-Boltzmann inelástica. Como segundo caso, sobre el que trabajaremos con más profundidad, abordaremos la situación en la que se inyecta energía al sistema mediante la aplicación de una fuerza de tipo “ruido blanco” (lo que se conoce en la bibliografía como “termostato estocástico” [11]). Esto introducirá en la ecuación de evolución de la distribución de velocidades un termino de difusión de Fokker-Planck, pero cuya solución puede desarrollarse de manera análoga al caso sin termostato. Estado de enfriamiento homogéneo (HCS).- Consideremos un gas granular homogéneo de densidad 𝑛y coeficiente de restitución 𝛼1sobre el cual no se aplica ninguna fuerza y que se deja evolucionar libremente (estado de enfriamiento homogéneo, HCS por sus siglas en inglés). Tras un transitorio, el sistema llega a un estado en que toda la dependencia temporal en la función de distribución de velocidades se produce a través de la temperatura (granular) [10]. En este régimen, podemos reescribir dicha función como 𝑓(𝑟, 𝑣, 𝑡) = 𝑛 (𝑣0(𝑡))𝑑 𝑓(𝑣 𝑣0(𝑡)),(3.1) donde 𝑑= 2,3es la dimensión donde estamos tratando el problema y 𝑇(𝑡) = 1 2𝑚𝑣0(𝑡)2. La temperatura granular se relaciona con el momento de segundo orden de la distribución de velocidades a través de la expresión 𝑑 2𝑛𝑇 (𝑡) = ∫𝑑 𝑣 1 2𝑚𝑣2𝑓(𝑣, 𝑡).(3.2) Es esperable que, bajo la hipótesis de baja inelasticidad, la solución dada por (3.1) sea próxima a la distribución de Maxwell. Para ello resulta muy conveniente desarrollar la función escalada  𝑓(𝑐), siendo 𝑐 =𝑣∕𝑣0(𝑡), en términos de los polinomios de Sonine 𝑆𝑝[11],  𝑓(𝑐) = 𝜙(𝑐)(1 + ∞ ∑ 𝑝=1 𝑎𝑝𝑆𝑝(𝑐2)).(3.3) 42 del cálculo variacional a la teoría del control óptimo: aplicación a sistemas lejos del eqilibrio puede ser 𝜒max −𝜒min o𝜒min −𝜒max (entendiendo el orden de aparición como el orden en el que toma los valores el control) puede intuirse, a nivel físico, la unicidad de dicho control óptimo una vez han sido fijadas unas condiciones inciales (para un argumento matemático más riguroso, el cual garantiza la unicidad del control, consultar el capítulo 3de la Ref. [4]) . La forma particular que adopte el control entre esas dos variantes dependerá de las condiciones iniciales que presente el sistema, de modo análogo a lo que vimos en nuestra discusión del problema de aparcamiento con aceleración acotada de la sección 2.3.2. Puntualizado este aspecto nos disponemos a calcular el tiempo 𝑡𝑗de cambio de los controles, que es el parámetro que falta para poder definir de manera unívoca 𝜒(𝑡). Dependiendo de las condiciones iniciales se presentan dos protocolos distintos:5 Para 𝑇𝑖>1: 𝜒(𝑡) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝑇3∕2 𝑖, 𝑠𝑖 𝑡 < 0, 𝜒min, 𝑠𝑖 0≤𝑡<𝑡𝑗, 𝜒max, 𝑠𝑖 𝑡𝑗≤𝑡<𝑡𝑓, 1, 𝑠𝑖 𝑡 ≤𝑡𝑓, (3.16) Para 𝑇𝑖<1: 𝜒(𝑡) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝑇3∕2 𝑖, 𝑠𝑖 𝑡 < 0, 𝜒max, 𝑠𝑖 0≤𝑡<𝑡𝑗, 𝜒min, 𝑠𝑖 𝑡𝑗≤𝑡<𝑡𝑓, 1, 𝑠𝑖 𝑡 ≤𝑡𝑓. (3.17) Los valores para 𝑡 < 0o𝑡≥𝑡𝑓garantizarán que el sistema parta de un estado estacionario 𝑇𝑖y, tras la aplicación del bang-bang, alcance otro estado estacionario con 𝑇𝑓= 1. Para los cálculos que siguen expresaremos todos los resultados finales en términos de las variables que, a nivel práctico, podemos fijar: 𝛿𝑇𝑖= (1 − 𝑇𝑖),𝜒min (o 𝛿𝜒min =𝜒min − 1) y 𝜒max (o 𝛿𝜒max =𝜒max − 1). 𝑇𝑖>1. Procesos de enfriamiento C-H (“cooling-heating”) En este apartado vamos a integrar la solución en dos ventanas de tiempo: una primera ventana temporal [0, 𝑡𝑗), en la cual actúa 𝜒min, y una segunda ventana [𝑡𝑗, 𝑡𝑓), donde actúa 𝜒max. Determinaremos el tiempo 𝑡𝑗de cambio así como el valor de las variables 𝛿𝑇𝑗y𝛿𝐴2𝑗en ese instante de tiempo, valores que conformarán el conjunto de condiciones iniciales para la resolución del sistema en la segunda ventana. Comencemos homogeneizando el sistema y hallando sus autovalores y autovectores. Para ello redefinimos las variables de forma que el sistema diferencial sea homogéneo. Este proceso de homogeneizado de las variables puede entenderse de manera gráfica a partir de la figura 3.4, donde se presenta el valor del control y de la variable temperatura a lo largo del tiempo. Las nuevas variables serán 𝛿𝑇 ∗=𝛿𝑇 −𝛿𝑇0, 𝛿𝐴∗ 2=𝛿𝐴2−𝛿𝐴20,(3.18) 5Hemos tomado como convenio que los controles sean continuos por la derecha. 3. minimización del tiempo de conexión 43 Figura 3.4: Control bang-bang y temperatura en función del tiempo para 𝑇𝑖>1(enfriamiento). En el panel izquierdo se presenta un control de tipo bang-bang, que pasa de un valor mínimo, 𝜒min = 1 + 𝛿𝜒min (línea azul) para 0≤𝑡<𝑡𝑗a uno máximo 𝜒max = 1 + 𝛿𝜒max (línea roja) para 𝑡𝑗≤0< 𝑡𝑓. Posteriormente, pasa a ser igual a la unidad (línea verde) una vez alcanzado el estado objetivo (𝑇𝑓= 1 y𝐴2𝑓= 1), a fin de mantener el sistema en el estado estacionario. En el panel derecho se recoge la evolución de la temperatura del sistema al ser sometida al control bang-bang. En primer lugar, relaja hacia la temperatura 𝑇min =𝜒min2∕3 correspondiente a 𝜒min (línea azul), relajación se ve interrumpida en 𝑡=𝑡𝑗. A continuación, comienza a relajar hacia la temperatura 𝑇max =𝜒max2∕3 correspondiente a 𝜒max (línea roja). Esta relajación también se ve interrumpida cuando 𝑇=𝑇𝑓= 1 (simultáneamente 𝐴2=𝐴2𝑓= 1), lo que determina el tiempo de conexión 𝑡𝑓. El sistema permanece estacionario para 𝑡>𝑡𝑓(línea verde). 44 del cálculo variacional a la teoría del control óptimo: aplicación a sistemas lejos del eqilibrio donde 𝛿𝑇0y𝛿𝐴20 se obtienen desarrollando a primer orden, 𝛿𝑇0= 1 − (𝜒𝑒𝑥𝑡𝑟𝑒𝑚𝑜)2∕3 =2 3𝛿𝜒𝑒𝑥𝑡𝑟𝑒𝑚𝑜, 𝛿𝐴20 = 0.(3.19) Con el subíndice extremo pretendemos engoblar tanto 𝛿𝑇0como 𝛿𝑇 ′ 0que aparecen en el dibujo, esto es, el caso en que el control se hace mínimo y el caso en que se hace máximo. De este modo el sistema diferencial a resolver quedaría como (𝛿𝑇 ∗′ 𝛿𝐴∗ 2 ′)=(−3 2𝛽1 − 𝛽 3 −2𝐵)(𝛿𝑇 ∗ 𝛿𝐴∗ 2),(3.20) y sus autovalores y autovectores serían: Para −𝜆1=1 2(−𝑘−3 2𝛽− 2𝐵)<0 𝑣1=(− 3 2𝛽−2𝐵+𝑘 6 1),(3.21) Para −𝜆2=1 2(𝑘−3 2𝛽− 2𝐵)<0 𝑣2=(− 3 2𝛽−2𝐵−𝑘 6 1),(3.22) donde hemos definido 𝑘≡𝜆1−𝜆2=√9 4𝛽2− (12 + 6𝐵)𝛽+(4𝐵2+ 12)>0.(3.23) Nótese que, en consecuencia, 𝜆1> 𝜆2. La solución general del sistema anterior es (𝛿𝑇 ∗ 𝛿𝐴∗ 2)=𝐶1𝑣1𝑒−𝜆1𝑡+𝐶2𝑣2𝑒−𝜆2𝑡.(3.24) A la hora de plantear el problema en cada ventana temporal hemos de ser cuidadosos con las variables homogeneizadoras, pues el término homogeneizador (𝛿𝑇0y𝛿𝑇 ′ 0) cambia en cada ventana al ser este función de 𝛿𝜒 y por tanto las variables homogeneizadas cambian en cada intervalo de tiempo. A fin de evitar confusiones rebautizamos las variables definidas en (3.18) en cada ventana, como se detalla en el cuadro 3.1. De este modo, tendremos los sistemas: En la primera ventana 𝑡∈ [0, 𝑡𝑗) (𝛿𝑇 − 𝛿𝐴− 2)=𝐶1𝑣1𝑒−𝜆1𝑡+𝐶2𝑣2𝑒−𝜆2𝑡. (3.25) En la segunda ventana 𝑡∈ [𝑡𝑗, 𝑡𝑓) (𝛿𝑇 + 𝛿𝐴+ 2)=𝐶3𝑣1𝑒−𝜆1𝑡+𝐶4𝑣2𝑒−𝜆2𝑡. (3.26) Primera ventana de tiempo [0, 𝑡𝑗).- Sustituyendo en el sistema (3.25) las condiciones iniciales (𝛿𝑇𝑖=𝑇𝑖− 1 >0, 𝛿𝐴2𝑖= 0) determinamos 𝐶1y𝐶2, 𝐶1= 3𝛿𝑇0−𝛿𝑇𝑖 𝑘= −3𝑇𝑖− 2𝜒min − 1 𝑘, 𝐶2=3𝑇𝑖− 2𝜒min − 1 𝑘.(3.27) 3. minimización del tiempo de conexión 45 Primera ventana Segunda ventana Variable de control 𝛿𝜒 =𝛿𝜒min =𝜒min − 1 𝛿𝜒 =𝛿𝜒max =𝜒max − 1 Variable homogeneizadora (T) 𝛿𝑇0=2 3𝛿𝜒min 𝛿𝑇 ′ 0=2 3𝛿𝜒max Variable homogeneizadora (𝐴2)𝛿𝐴20 = 0 𝛿𝐴′ 20 = 0 Variable homogeneizada (T) 𝛿𝑇 −=𝛿𝑇 −𝛿𝑇0𝛿𝑇 +=𝛿𝑇 −𝛿𝑇 ′ 0 Variable homogeneizada (𝐴2)𝛿𝐴− 2=𝛿𝐴2𝛿𝐴+ 2=𝛿𝐴2 Cuadro 3.1: Variables en cada ventana temporal del protocolo C-H. La primera ventana corresponde al intervalo [0, 𝑡𝑗)(cooling) y la segunda ventana al intervalo [𝑡𝑗, 𝑡𝑓)(heating). El punto sobre la trayectoria en el instante 𝑡𝑗será 𝛿𝑇𝑗=𝛿𝑇 − 𝑗+𝛿𝑇0=𝐶1(𝑣1(1)𝑒−𝜆1𝑡𝑗−𝑣2(1)𝑒−𝜆2𝑡𝑗)+𝛿𝑇0, 𝛿𝐴2𝑗=𝛿𝐴− 2𝑗+𝛿𝐴20 =𝐶1(𝑒−𝜆1𝑡𝑗−𝑒−𝜆2𝑡𝑗).(3.28) Segunda ventana de tiempo [𝑡𝑗, 𝑡𝑓).- Las condiciones iniciales para el sistema (3.26) vienen suministradas por el apartado anterior (𝛿𝑇𝑗, 𝛿𝐴2𝑗). Teniendo en cuenta que conocemos los valores que deben tomar las variables en el tiempo 𝑡𝑓(𝛿𝑇𝑓= 0, 𝛿𝐴2𝑓= 0) (o equivalentemente para las variables homogeneizadas (𝛿𝑇 + 𝑓= −𝛿𝑇 ′ 0, 𝛿𝐴+ 2𝑓= 0)), podemos plantear un sistema de cuatro ecuaciones con cuatro variables (𝐶3, 𝐶4, 𝑡𝑗, 𝑡𝑓), (𝛿𝑇𝑗−𝛿𝑇 ′ 0 𝛿𝐴2𝑗)=𝐶3𝑣1𝑒−𝜆1𝑡𝑗+𝐶4𝑣2𝑒−𝜆2𝑡𝑗,(−𝛿𝑇 ′ 0 0)=𝐶3𝑣1𝑒−𝜆1𝑡𝑓+𝐶4𝑣2𝑒−𝜆2𝑡𝑓.(3.29) Resolvemos para las variables 𝐶3,𝐶4en términos del tiempo de cambio 𝑡𝑗, 𝐶3=2(𝜒max − 1) 𝑘⎛⎜⎜⎝ 2(𝜒max−𝜒min)𝑒𝜆1𝑡𝑗 3𝑇𝑖−2𝜒min−1 − 1 2(𝜒max−𝜒min)𝑒𝜆2𝑡𝑗 3𝑇𝑖−2𝜒min−1 − 1⎞⎟⎟⎠ 𝜆1 𝑘 , 𝐶4= −𝐶3𝑒−𝑘𝑡𝑓,(3.30) así como para el tiempo final de conexión 𝑡𝑓=1 𝑘log ⎛⎜⎜⎝ 2(𝜒max−𝜒min)𝑒𝜆1𝑡𝑗 3𝑇𝑖−2𝜒min−1 − 1 2(𝜒max−𝜒min)𝑒𝜆2𝑡𝑗 3𝑇𝑖−2𝜒min−1 − 1⎞⎟⎟⎠ ,(3.31) mientras que 𝑡𝑗está determinado de modo implícito por la solución de la ecuación (3𝑇𝑖− 2𝜒min − 1) (2(𝜒max −𝜒min)𝑒𝜆2𝑡𝑗 3𝑇𝑖− 2𝜒min − 1 − 1)𝜆1 𝑘 = 2(𝜒max − 1) (2(𝜒max −𝜒min)𝑒𝜆1𝑡𝑗 3𝑇𝑖− 2𝜒min − 1 − 1)𝜆2 𝑘 .(3.32) 46 del cálculo variacional a la teoría del control óptimo: aplicación a sistemas lejos del eqilibrio Figura 3.5: Control bang-bang y temperatura en función del tiempo para 𝑇𝑖<1(calentamiento). En el panel izquierdo se presenta un control de tipo bang-bang, que pasa de un valor máximo, 𝜒max = 1 + 𝛿𝜒max (línea roja) para 0≤𝑡<𝑡𝑗a uno mínimo 𝜒min = 1 + 𝛿𝜒min (línea azul) para 𝑡𝑗≤𝑡<𝑡𝑓. Posteriormente, pasa a ser igual a la unidad (línea verde) una vez alcanzado el estado objetivo (𝑇𝑓= 1 y𝐴2𝑓= 1), a fin de mantener el sistema en el estado estacionario. En el panel derecho se recoge la evolución de la temperatura del sistema al ser sometida al control bang-bang. En primer lugar, relaja hacia la temperatura 𝑇max =𝜒max2∕3 correspondiente a 𝜒max (línea roja), relajación se ve interrumpida en 𝑡=𝑡𝑗. A continuación, comienza a relajar hacia la temperatura 𝑇min =𝜒min2∕3 correspondiente a 𝜒min (línea azul). Esta relajación también se ve interrumpida cuando 𝑇=𝑇𝑓= 1 (simultáneamente 𝐴2=𝐴2𝑓= 1), lo que determina el tiempo de conexión 𝑡𝑓. El sistema permanece estacionario para 𝑡>𝑡𝑓(línea verde). Nótese la analogía con la Figura 3.4, con la diferencia fundamental de la inversión del orden de los bangs; en ambas figuras, las curvas rojas (azules) corresponden al bang con 𝜒max (𝜒min). Recordemos que 𝑘está definido en la ecuación (3.23). El conjunto de las ecuaciones (3.31) y (3.32) nos proporciona una solución analítica para el tiempo mínimo de conexión 𝑡𝑓en el protocolo C-H, que es válida en la aproximación lineal que estamos considerando en esta memoria. 𝑇𝑖<1. Procesos de calentamiento H-C (“heating-cooling”) Partamos ahora de un estado con 𝑇𝑖<1, de modo que el protocolo a seguir será calentar el sistema, sometiéndolo a 𝜒max, para después enfriarlo a 𝜒min como se muestra en la figura 3.5. Para estas condiciones iniciales se invierten, respecto del caso anterior, las variables homogeneizadoras y homogeneizadas en cada ventana de tiempo, tal y como se recoge en el cuadro 3.2, teniéndose los sistemas: En la primera ventana 𝑡∈ [0, 𝑡𝑗) (𝛿𝑇 + 𝛿𝐴+ 2)=𝐶5𝑣1𝑒−𝜆1𝑡+𝐶6𝑣2𝑒−𝜆2𝑡. (3.33) En la segunda ventana 𝑡∈ [𝑡𝑗, 𝑡𝑓) (𝛿𝑇 − 𝛿𝐴− 2)=𝐶7𝑣1𝑒−𝜆1𝑡+𝐶8𝑣2𝑒−𝜆2𝑡. (3.34) 3. minimización del tiempo de conexión 47 Primera ventana Segunda ventana Variable de control 𝛿𝜒 =𝛿𝜒max =𝜒max − 1 𝛿𝜒 =𝛿𝜒min =𝜒min − 1 Variable homogeneizadora (T) 𝛿𝑇 ′ 0=2 3𝛿𝜒max 𝛿𝑇0=2 3𝛿𝜒min Variable homogeneizadora (𝐴2)𝛿𝐴′ 20 = 0 𝛿𝐴20 = 0 Variable homogeneizada (T) 𝛿𝑇 +=𝛿𝑇 −𝛿𝑇 ′ 0𝛿𝑇 −=𝛿𝑇 −𝛿𝑇0 Variable homogeneizada (𝐴2)𝛿𝐴+ 2=𝛿𝐴2𝛿𝐴− 2=𝛿𝐴2 Cuadro 3.2: Variables en cada ventana temporal para el protocolo H-C. La primera ventana corresponde al intervalo [0, 𝑡𝑗)(heating) y la segunda ventana al intervalo [𝑡𝑗, 𝑡𝑓)(cooling). Primera ventana de tiempo [0, 𝑡𝑗)-. En esta ventana de tiempo actúa 𝜒max y las condiciones iniciales con las que completar el sistema lineal serán (𝑇𝑖, 𝐴2𝑖= 1), lo que en términos de las variables infinitesimales es (𝛿𝑇𝑖=𝑇𝑖− 1 <0, 𝛿𝐴2𝑖= 0). Sustituyendo la condición inicial en (3.33) despejamos 𝐶5y𝐶6, 𝐶5=2𝜒max + 1 − 3𝑇𝑖 𝑘, 𝐶6= −2𝜒max + 1 − 3𝑇𝑖 𝑘.(3.35) Evaluando el sistema en el tiempo 𝑡𝑗(a determinar aún) tendríamos la condición inicial para el siguiente intervalo temporal. Segunda ventana de tiempo [𝑡𝑗, 𝑡𝑓).- En esta ventana de tiempo actúa 𝜒min y las condiciones iniciales con las que completar el sistema lineal serán (𝛿𝑇𝑗, 𝛿𝐴2𝑗). Teniendo en cuenta que conocemos los valores que deben tomar las variables en el tiempo 𝑡𝑓(𝛿𝑇𝑓= 0, 𝛿𝐴2𝑓= 0) podemos evaluar el sistema en el tiempo 𝑡𝑗y𝑡𝑓, de manera análoga a como procedimos para 𝑇𝑖>1, y obtener un sistema de 4ecuaciones con 4incógnitas (𝐶7, 𝐶8, 𝑡𝑓y𝑡𝑗), (𝛿𝑇𝑗−𝛿𝑇0 𝛿𝐴2𝑗)=𝐶7𝑣1𝑒−𝜆1𝑡𝑗+𝐶8𝑣2𝑒−𝜆2𝑡𝑗,(−𝛿𝑇0 0)=𝐶7𝑣1𝑒−𝜆1𝑡𝑓+𝐶8𝑣2𝑒−𝜆2𝑡𝑓.(3.36) Despejando obtenemos las constantes 𝐶7, 𝐶8, 𝐶7= −2(1 − 𝜒min) 𝑘⎛⎜⎜⎝ 2(𝜒max−𝜒min)𝑒𝜆1𝑡𝑗 2𝜒max−3𝑇𝑖+1 − 1 2(𝜒max−𝜒min)𝑒𝜆2𝑡𝑗 2𝜒max−3𝑇𝑖+1 − 1⎞⎟⎟⎠ 𝜆1 𝑘 , 𝐶8= −𝐶7𝑒−𝑘𝑡𝑓,(3.37) así como el tiempo de conexión 𝑡𝑓en términos del tiempo de cambio 𝑡𝑗, 𝑡𝑓=1 𝑘log ⎛⎜⎜⎝ 2(𝜒max−𝜒min)𝑒𝜆1𝑡𝑗 2𝜒max−3𝑇𝑖+1 − 1 2(𝜒max−𝜒min)𝑒𝜆2𝑡𝑗 2𝜒max−3𝑇𝑖+1 − 1⎞⎟⎟⎠ ,(3.38) 48 del cálculo variacional a la teoría del control óptimo: aplicación a sistemas lejos del eqilibrio para el cual obtenemos de nuevo una expresión implícita de la forma 𝑔(𝑡𝑗)=0, (2𝜒max − 3𝑇𝑖+ 1) (2(𝜒max −𝜒min)𝑒𝜆2𝑡𝑗 2𝜒max − 3𝑇𝑖+ 1 − 1)𝜆1 𝑘 = 2(1 − 𝜒min)(2(𝜒max −𝜒min)𝑒𝜆1𝑡𝑗 2𝜒max − 3𝑇𝑖+ 1 − 1)𝜆2 𝑘 (3.39) De nuevo, el conjunto de las ecuaciones (3.38) y (3.39) nos da el tiempo mínimo de conexión, en este caso para el protocolo H-C, dentro de la aproximación lineal. 4 Análisis de los resultados Una vez obtenidas predicciones analíticas para las variables físicas de interés (fundamentalmente las trayectorias de la temperatura 𝑇y de 𝐴2, y los tiempos de cambio 𝑡𝑗y final 𝑡𝑓), nos disponemos a inspeccionar el comportamiento de las mismas para entender la respuesta del sistema al control óptimo que hemos diseñado.1A nivel intuitivo, esperamos encontrar una tendencia decreciente en el tiempo final conforme aumentamos la capacidad de calentamiento o enfriamien- to del termostato. Del mismo modo, esperamos comprobar la necesidad de un protocolo en dos partes al estudiar la evolución temporal de las variables 𝑇y𝐴2y poder examinar los límites del cálculo anteriormente expuesto, comparándolos con el caso límite 𝜒max →∞y𝜒min = 0 investigado en la Ref. [14]. 4.1 Tiempo de cambio y final como función de la intensidad del termostato Comenzamos estudiando el comportamiento de los tiempos finales y de cambio frente a la variación de la intensidad del termostato en la figura 4.1 para ambos protocolos, C-H y H-C. Una primera cuestión relevante es el rango de valores de 𝜒min y𝜒max que permiten conectar los estados inicial y final. En el régimen no lineal y en el caso límite (𝜒min = 0, 𝜒max = ∞) siempre se pueden conectar los estados estacionarios correspondientes a dos temperaturas (granulares) 𝑇𝑖y𝑇𝑓[14]. Sin embargo, no es trivial que esto suceda cuando no se dispone de toda la potencia del termostato, esto es, cuando (𝜒min >0, 𝜒max <∞). En consecuencia, nos proponemos estudiar este problema dentro de la aproximación lineal realizada en este trabajo. En el caso C-H, esperamos físicamente que el tiempo de conexión diverja si 𝜒min ≥1. Esto se entiende intuitivamente: si la menor intensidad del termostato es mayor que la correspondiente a la temperatura final, esto es, la unidad, no podremos realizar la conexión. Si la intensidad mínima es exactamente igual a la unidad, podemos llegar al estado estacionario objetivo pero tardaremos un tiempo infinito, esto es, 𝑡𝑓→∞para 𝜒min →1−. Demostramos rigurosamente que este es 1En el Anexo a esta memoria se pueden encontrar algunos ejemplos relevantes de los códigos de Matlab desarrollados y utilizados para realizar este análisis. 50 del cálculo variacional a la teoría del control óptimo: aplicación a sistemas lejos del eqilibrio efectivamente el caso en el Apéndice, donde además vemos que el tiempo de conexión diverge en esta situación debido a que 𝑡𝑗→∞. Es la primera parte del bang-bang, en que el sistema se enfría con 𝜒min, la responsable de la divergencia del tiempo de conexión. Análogamente, 𝑡𝑓→∞cuando 𝜒max →1+, como se demuestra en el Apéndice. Sin embargo, es la segunda parte del bang-bang la que produce la divergencia del tiempo de relajación, ya que el tiempo de cambio 𝑡𝑗tiende en este caso a una constante. De manera análoga, en el protocolo H-C el tiempo de conexión diverge cuando (i) la cota superior 𝜒max →1+o (ii) la cota inferior 𝜒min →1−. En el caso (i), la incapacidad del sistema para calentar más allá de 𝑇𝑓= 1 en la primera parte del bang-bang es la responsable de la divergencia, y consistentemente 𝑡𝑗→∞también. En el caso (ii), la divergencia está marcada por la segunda parte del bang-bang, ya que la incapacidad del sistema para enfriar implica que esta debe durar un tiempo infinito. Fijados dichos límites procedemos a estudiar el comportamiento de los tiempos 𝑡𝑗y𝑡𝑓como función de la intensidad del termostato. A la luz de los paneles a la derecha de la figura 4.1 podemos comprobar como, conforme 𝜒max tiende a 1+, los tiempos 𝑡𝑓divergen en ambos protocolos. Por otro lado, en lo que a los tiempos de cambio 𝑡𝑗respecta, el comportamiento para los procesos de tipo C-H y H-C es distinto. En el caso C-H, para alcanzar el estado estacionario cuando 𝜒max tiende a 1+, el tiempo 𝑡𝑗de aplicación de 𝜒min debe hacerse mínimo al carecer de capacidad de calentamiento y el tiempo final diverge debido a la segunda parte del bang-bang. En el caso H-C, esta tendencia para 𝑡𝑗se invierte, al cambiar el orden del protocolo, y 𝑡𝑗tiende a 𝑡𝑓: el tiempo de conexión diverge debido a la primera parte del bang-bang. Esto es consistente con la discusión general que hemos realizado anteriormente y con la intuición física: 𝑡𝑓diverge cuando 𝜒max tiende a1+debido a que la etapa del bang-bang del calentamiento debe durar un tiempo infinito. Por otra parte, podemos comprobar en ambas gráficas como un aumento de 𝜒max genera una disminución de 𝑡𝑓. Esta tendencia decreciente puede explicarse teniendo en cuenta que, al incrementar el valor máximo del termostato, aumentamos el conjunto de controles sobre el cual optimizamos el tiempo y, por tanto, obtenemos tiempos finales menores que los obtenidos para rangos menores de 𝜒max. El tiempo de actuación de 𝜒max debe ser menor, para no sobrepasar 𝑇𝑓= 1, lo que explicaría el crecimiento de 𝑡𝑗para el protocolo C-H y el decrecimiento del mismo en el protocolo H-C. Para los paneles de la izquierda, donde se recoge la evolución de los tiempos respecto a 𝜒min, el análisis es similar. Para ambos protocolos tenemos que 𝑡𝑓diverge cuando 𝜒min →1−. La ausencia de capacidad para enfriar provoca nuevamente que el sistema tarde un tiempo infinito en conectar los estados inicial y final; en el caso H-C el segundo paso del protocolo deberá aplicarse un tiempo infinito y es el responsable de que 𝑡𝑓→∞, mientras que en el caso C-H es el primer paso del protocolo el que consume un tiempo 𝑡𝑗→∞. En ambos casos es la parte C (“cooling”) del protocolo la responsable de la divergencia. Por otra parte, el tiempo de conexión 𝑡𝑓disminuye conforme 𝜒min lo hace, ya que estaríamos nuevamente ampliando el conjunto de controles sobre el que minimizamos 𝑡𝑓. Con respecto a los tiempos de cambio podemos verificar como, análoga- 4. análisis de los resultados 51 0 0.2 0.4 0.6 0.8 0.5 1 1.5 2 2.5 3 3.5 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1 1.2 1.4 1.6 1.8 2 0 0.5 1 1.5 2 2.5 Figura 4.1: Tiempo de cambio 𝑡𝑗(azul) y tiempo de conexión 𝑡𝑓(rojo) como función de los valores límite del termostato. Se representan cuatro paneles: los paneles de la izquierda representan los tiempos como función de 𝜒min, para un valor fijo de 𝜒max = 1.1, mientras que los paneles de la derecha muestran los tiempos como función de 𝜒max, para un valor fijo de 𝜒min = 0.9; los paneles de arriba corresponden al caso C-H, con 𝑇𝑖= 1.01, mientras que los de abajo corresponden al caso H-C, con 𝑇𝑖= 0.99. En todos los casos, se ha tomado 𝛼= 0.9y𝑑= 2. No existe un cambio cualitativo al variar el valor de 𝛼, creciendo el tiempo final al disminuir 𝛼. El cambio a 𝑑= 3 apenas modifica las figuras. 58 del cálculo variacional a la teoría del control óptimo: aplicación a sistemas lejos del eqilibrio 0 0.5 1 0.2 0.4 0.6 0.8 1 0 0.5 1 0.5 0.505 0.51 0.515 0.52 0 0.5 1 0.2 0.4 0.6 0.8 1 0 0.5 1 0.5 0.505 0.51 0.515 0.52 0.2 0.4 0.6 0.8 0.01 0.014 0.018 0.022 0.2 0.4 0.6 0.8 0.01 0.014 0.018 0.022 Figura 4.6: Convergencia del ajuste (4.2) para los tiempos finales de conexión hacia la expresión (4.1), para el caso C-H. De izquierda a derecha, consideramos valores crecientes de 𝜒max y decrecientes de 𝜒min:𝜒max = 5 y𝜒min = 0.05 (izquierda) y 𝜒max = 100 y𝜒min = 0.00 (derecha). Arriba: en cada uno de los paneles izquierdos se representa 𝑝1(rojo) y 1 √𝐵(azul) frente a 𝛼mientras que en los derechos se representa 𝑝2(rojo) y1∕2 (azul) frente a 𝛼. Abajo: Tiempo final 𝑡𝑓frente a 𝑇𝑖para 𝛼= 0.6, según la ecuación (4.1) (línea azul), según (3.31) (estrellas negras), y ajuste según la ecuación (4.2) (línea roja). La convergencia hacia la ecuación (4.1) es “desde arriba”, ya que el tiempo óptimo de conexión debe ser mayor si el conjunto de controles admisibles es más restrictivo. En todos los casos 𝑑= 2. 4. análisis de los resultados 59 0 0.5 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0 0.5 1 0.45 0.46 0.47 0.48 0.49 0.5 0 0.5 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0 0.5 1 0.45 0.46 0.47 0.48 0.49 0.5 2 4 6 8 10 0.03 0.04 0.05 0.06 0.07 2 4 6 8 10 0.03 0.04 0.05 0.06 0.07 Figura 4.7: Convergencia del ajuste (4.2) para los tiempos finales de conexión hacia la expresión (4.1), para el caso H-C. De izquierda a derecha, consideramos valores crecientes de 𝜒max y decrecientes de 𝜒min:𝜒max = 5 y𝜒min = 0.05 (izquierda) y 𝜒max = 100 y𝜒min = 0.00 (derecha). Arriba: en cada uno de los paneles izquierdos se representa 𝑝1(rojo) y 1 √𝐵(azul) frente a 𝛼mientras que en los derechos se representa 𝑝2(rojo) y1∕2 (azul) frente a 𝛼. Abajo: Tiempo final 𝑡𝑓frente a 𝑇𝑖para 𝛼= 0.6, según la ecuación (4.1) (línea azul), según (3.38) (estrellas negras), y ajuste según la ecuación (4.2) (línea roja). Observamos cómo la tendencia natural de convergencia se invierte para valores muy altos de 𝜒max. En todos los casos 𝑑= 2. 60 del cálculo variacional a la teoría del control óptimo: aplicación a sistemas lejos del eqilibrio el caso C-H, como se ha mostrado en las figuras 4.6 y 4.7. Si eliminásemos la restricción física impuesta sobre 𝜒min y estudiásemos a nivel puramente matemático el sistema cuando el valor de la intensidad puede variar entre −∞ y+∞ comprobaríamos cómo desaparece este comportamiento asimétrico del sistema ante ambos protocolos. También hemos analizado la dependencia respecto de 𝛼de la tendencia de los tiempos finales, dados por (3.31) y (3.38), hacia a la predicción no lineal (4.1). Se encuentra que, para ambos protocolos, al disminuir 𝛼(aumentar la inelasticidad) disminuye el valor de 𝜒max y aumenta el de 𝜒min necesarios para invertir el orden natural de convergencia, es decir, la aproximación lineal realizada deja de ser válida para valores más tempranos de intensidad. Sin embargo, también se aprecian diferencias entre ambos protocolos. Si fijamos en ambos casos 𝛼= 0.1(que puede considerarse como un régimen de alta inelasticidad) podemos verificar cómo, para el protocolo C-H, los tiempos finales de la aproximación lineal son superiores a los de la fórmula (4.1) hasta valores 𝜒max = 85 y𝜒min = 0.3Por otra parte, en el caso H-C los tiempos finales de la aproximación lineal exceden los de la fórmula (4.1) hasta valores 𝜒max = 9 y𝜒min = 0, que son menores que los del protocolo C-H.4La asimetría entre ambos protocolos vuelve a ser reflejo de la asimetría existente entre 𝜒min y𝜒max anteriormente expuesta. Resumiendo, como tendencia general los tiempos finales del caso C-H presentan un mejor ajuste a la fórmula (4.1) para valores más extremos del termostato (en todo el rango de 𝛼) mientras que los tiempos finales del protocolo H-C, aunque también presentan un buen ajuste a la predicción (4.1), se alejan del comportamiento marcado por esta para valores menos extremos de los controles, como puede verificarse en los paneles inferiores de la figura 4.7 para 𝛼= 0.6. Otro rasgo general en ambos protocolos es que los controles a partir de los cuales deja de ser válida la aproximación lineal (esto es, se invierte la convergencia natural “desde arriba”) se vuelven menos extremos al disminuir 𝛼. Este comportamiento con 𝛼puede entenderse teniendo en cuenta que 𝑎𝑠 2aumenta al disminuir 𝛼y, en consecuencia, la importancia del término proporcional a la intensidad 𝜒en la ecuación (3.5) de evolución de la temperatura también aumenta: con un valor más pequeño de 𝜒tenemos el mismo valor de la derivada de la temperatura. 3Para 𝛼= 0.6no hemos encontrado valores para los controles que invirtiesen la tendencia “desde arriba”. El primer valor de 𝛼para el que encontramos la inversión, al fijar 𝜒min = 0 y𝜒max →∞, es 𝛼= 0.45. 4Por otra parte, también son menores que los encontrados para 𝛼= 0.6en el caso H-C. 4. análisis de los resultados 61 9.975 9.98 9.985 9.99 9.995 10 0.07085 0.0709 0.07095 0.071 Figura 4.8: Inversión de la convergencia del tiempo óptimo de conexión. Representamos los tiempos finales 𝑡𝑓frente a 𝑇𝑖según las ecuaciones (3.38) (estrellas negras) y según (4.1) (línea azul), para 𝛼= 0.6en el caso H-C. Encontramos que para 𝜒max = 10.1y𝜒min = 0.0el tiempo final es menor que el proporcionado por la ecuación (4.1), lo que no tiene sentido físico ya que esta corresponde al caso 𝜒max = ∞ y𝜒min = 0. Esto marca un límite de validez para la teoría lineal desarrollada en este trabajo. 5 Conclusiones El primer objetivo de este trabajo ha sido introducir la teoría del control óptimo, que culmina en el principio del máximo, desde la óptica que ofrece el cálculo variacional. Esta vía, sin ser la forma más rigurosa de exponer el formalismo matemático, ofrece una perspectiva más intuitiva desde el punto de vista físico. El segundo objetivo de este trabajo ha sido utilizar la potencia del método matemático expuesto para resolver un problema físico relevante: la conexión óptima (en el sentido de minimizar el tiempo de conexión) entre dos estados estacionarios de un sistema fuera del equilibrio. Concretamente, hemos considerado la conexión para un gas granular sometido a un termostato externo cuando los dos estados estacionarios son muy próximos, lo que nos ha permitido introducir una aproximación lineal y obtener soluciones analíticas exactas. A continuación, presentamos los principales resultados recogidos en este trabajo: Hemos deducido las condiciones de primer y segundo orden para que un funcional alcance su valor mínimo, a partir del concepto de primera y segunda variación del mismo. Estas condiciones conforman la columna vertebral del cálculo variacional y dan lugar a las famosas ecuaciones de Euler-Lagrange y a las condiciones de Legendre. Hemos empleado las herramientas que nos ofrece el cálculo en variaciones para esbozar un primer esquema del principio del máximo, desarrollar una intuición del mismo a partir de los conocimientos adquiridos y mostrar la conexión entre ambas teorías que pudiera diluirse tras el formalismo matemático. Hemos formulado el principio del máximo en dos casos particulares, cuyas condiciones son fácilmente generalizables a otros problemas, y una variante del mismo aplicable a sistemas lineales. Así mismo, hemos mostrado las ventajas que la teoría del control ofrece frente al cálculo variacional para un ejemplo sencillo, como puede ser aparcar un coche en el menor tiempo posible y con una aceleración restringida. Esto es importante porque este tipo de restricciones no holónomas aparecen en muchos problemas físicos relevantes. Hemos introducido una descripción de sistemas granulares en el régimen de flujo rápido (gases granulares), exponiendo sus rasgos característicos y las ecuaciones que gobiernan su evolución temporal, tanto en el caso en el que dicho sistema se encuentra aislado (enfriamiento homogéneo) como en el caso en el cual está sometido a un forzamiento externo, debido a la presencia de un termostato. 64 del cálculo variacional a la teoría del control óptimo: aplicación a sistemas lejos del eqilibrio Hemos linealizado el sistema de ecuaciones para la temperatura y la kurtosis de un sistema granular sometido a un termostato de intensidad próxima a la unidad, hemos verificado las hipótesis del teorema lineal enunciado anteriormente y lo hemos aplicado para resolver el sistema. Hemos resuelto analíticamente el sistema lineal planteado, obteniendo unas soluciones cerradas para el protocolo a seguir y para la evolución completa del sistema. Hemos analizado el comportamiento de dichas soluciones para contrastarlas con la intuición que pudiéramos tener previa a la resolución. Hemos explorado los límites de la aproximación lineal, forzando la intensidad del termosta- to a valores lejos de la unidad y comparando los tiempos finales obtenidos bajo la hipótesis de linealidad con los resultados obtenidos en la Ref. [14], encontrando un muy buen compromiso entre ambos enfoques. Los resultados obtenidos en este trabajo son aplicables cuando tratamos gases granulares con temperaturas iniciales y finales próximas entre sí y un rango muy estrecho para la intensidad del termostato. Por otra parte, en la Ref. [14] se cubre el caso en el que, siendo las temperaturas final e inicial cualesquiera, se dispone de toda lo potencia posible del termostato, lo cual resulta esencialmente diferente del planteamiento considerado en esta memoria. Esto abre, a su vez, nuevas posibilidades de estudio, en particular la posibilidad de resolver el problema de conexión óptima (en el sentido de minimizar el tiempo de conexión) en una situación completamente general. Esto es, encontrar el tiempo mínimo de conexión 𝑡𝑓para temperaturas inicial y final arbitrarias, así como valores máximo y mínimo de la intensidad del termostato arbitrarias, lo que nos permitiría caracterizar de manera completa el comportamiento del sistema. En este contexto, es relevante verificar si la solución del “problema general” converge, y (de hacerlo) cómo lo hace, al régimen de respuesta lineal considerado aquí y al régimen de máxima potencia del termostato estudiado en la Ref. [14]. A Apéndice. Valores límite de la intensidad del termostato En este apéndice, estudiamos el rango de valores de la intensidad del termostato en ambos protocolos bang-bang, C-H y H-C, que permiten la conexión entre los dos estados estacionarios de no equilibrio. Esto es, queremos analizar si es posible realizar la conexión para cualquier pareja de valores (𝜒min, 𝜒max), o si aparecen zonas en el espacio de parámetros que hacen imposible la conexión. Analizamos el caso C-H (enfriamiento, 𝑇𝑖>1) en detalle, ya que el caso H-C es completamente análogo. Físicamente es claro que para poder alcanzar la temperatura final 𝑇𝑓= 1 objetivo, es imprescindible que 𝜒min <1. Esto es, necesitamos que el valor mínimo de la intensidad del termostato sea menor que el correspondiente a la temperatura final deseada. Pero, ¿qué sucede con 𝜒m ax? ¿Es posible la conexión siempre que 𝜒max >1o hay una cota inferior, bajo la cual no es posible conectar los dos estados? Los valores límite (𝜒min, 𝜒max)son aquellos que conduzcan a la divergencia del tiempo de conexión óptimo 𝑡𝑓. A partir de la ecuación (3.31), que define el tiempo 𝑡𝑓como función de 𝑇𝑖,𝜒min y𝜒max, 𝑡𝑓=1 𝑘log ⎛⎜⎜⎝ 2(𝜒max−𝜒min)𝑒𝜆1𝑡𝑗 3𝑇𝑖−2𝜒min−1 − 1 2(𝜒max−𝜒min)𝑒𝜆2𝑡𝑗 3𝑇𝑖−2𝜒min−1 − 1⎞⎟⎟⎠ ,(A.1) hallamos los valores de los controles que hacen que 𝑡𝑓→∞. La única manera de conseguir que 𝑡𝑓diverja es hacer que, o bien el numerador de la fracción de logaritmo tienda a ∞o bien el denominador de la misma tienda a 0. A continuación analizamos ambas posibilidades. Para que el numerador del logaritmo en (A.1) tienda a ∞es necesario que, o bien 𝜒m ın → (3𝑇𝑖− 1)∕2 >1,1o que 𝑡𝑗tienda a ∞. En este último caso es evidente que 𝑡𝑓divergirá siempre. Físicamente, es claro que el tiempo de conexión no puede ser menor que el tiempo de cambio, 𝑡𝑓≥𝑡𝑗. Matemáticamente, como 𝜆1> 𝜆2, se cumple que la divergencia del numerador siempre gana, de hecho 𝑡𝑓∼𝑡𝑗(el tiempo de la segunda parte del bang-bang se vuelve despreciable frente 1De acuerdo con la discusión del párrafo anterior, esta solución no tiene sentido físico. 66 del cálculo variacional a la teoría del control óptimo: aplicación a sistemas lejos del eqilibrio a𝑡𝑗). La ecuación (3.32) que define 𝑡𝑗, esto es, (3𝑇𝑖− 2𝜒min − 1) (2(𝜒max −𝜒min)𝑒𝜆2𝑡𝑗 3𝑇𝑖− 2𝜒min − 1 − 1)𝜆1 𝑘 = 2(𝜒max − 1) (2(𝜒max −𝜒min)𝑒𝜆1𝑡𝑗 3𝑇𝑖− 2𝜒min − 1 − 1)𝜆2 𝑘 ,(A.2) nos dice que cuando 𝑡𝑗→∞tenemos (𝜒max −𝜒min)∼(𝜒max − 1).(A.3) Por tanto, 𝑡𝑗→∞cuando 𝜒min →1−, de acuerdo con la intuición física antes descrita.2 Exploremos la segunda posibilidad, esto es, que el denominador del logaritmo en (A.1) se anule. Vemos que esto ocurre cuando 𝑡𝑗→1 𝜆2 log (3𝑇𝑖− 2𝜒min − 1 2(𝜒max −𝜒min)),(A.4) lo que mantiene el numerador del logaritmo en (A.1) no nulo. La anulación del denominador implica que el primer miembro de la ecuación (A.2) es nulo y, en consecuencia 2(𝜒max − 1) = 0, ya que el factor que lo acompaña en el segundo miembro es básicamente el numerador del logaritmo en (A.1). Por tanto, 𝑡𝑓→∞cuando 𝜒max →1+. Nuestro análisis implica que es posible conectar los estados estacionarios de no equilibrio inicial y final, siempre que la cota inferior (superior) de la intensidad del termostato, 𝜒min (𝜒max) esté por debajo (encima) de la correspondiente al estado final (la unidad con nuestra elección de variables adimensionales). Esto es, no aparecen zonas de conexión prohibidas más allá de la imposibilidad de enfriar el sistema a una temperatura final que sea menor (mayor) que la correspondiente a la máxima (mínima) intensidad del termostato disponible. Recordamos que el caso H-C (calentamiento, 𝑇𝑖<1) se trataría de un modo completamente análogo, y obtendríamos como conclusión también que 𝜒min <1y𝜒max >1, con sus papeles intercambiados respecto al caso C-H. 2Esto hace innecesario, matemáticamente, estudiar la otra posibilidad de divergencia de 𝑡𝑗,𝜒m ın →(3𝑇𝑖−1)∕2 >1. Anexo Adjuntamos a continuación los códigos de Matlab empleados para la resolución de la dinámica del sistema y la generación de las gráficas expuestas a lo largo de esta memoria. Nos ceñiremos al caso C-H, pues el caso H-C es completamente análogo. Hemos de indicar que la notación empleada en los códigos es ligeramente distinta a la utilizada en la memoria, incluyendo un parámetro 𝜖que se fija igual a 1en los desarrollos expuestos con anterioridad (puesto que su valor, realmente, no influye de modo alguno en los resultados). Tiempo 𝑡𝑓y𝑡𝑗frente a 𝜒max (el caso de 𝜒min es completamente análogo): 1%% Valores fijados 2Tfin =1; 3Ti=1.01; 4alfa =0.9; d =2; 5XIMAX =2; 6chimin =0.9; 7%% Par á metros 8epsilon=abs (1-Ti); 9deltaTi =(Ti -1) / epsilon ; % condici ón inical para delta T 10 deltaximin =( chimin -1) / epsilon ; % Valor mí nimo de delta xi 11 deltaximax =@( chimax ) ( chimax -1) ./ epsilon ; % Valor maximo de delta xi 12 deltaT0 =2/3* deltaximin ; 13 deltaT0_p =@( chimax ) 2/3* deltaximax ( chimax ); 14 a_s =(16*(1 - alfa )*(1 -2* alfa ^2) ) /(73+56*d -24* d*alfa -105* alfa +30*(1 - alfa )* alfa ^2); 15 a_h =(16*(1 - alfa )*(1 -2* alfa ^2)) /(25+2* alfa ^2*( alfa -1) +24*d+ alfa *(8*d -57)); 16 B= a_h /( a_h - a_s); 17 beta=1+3/16* a_s ; 18 b1 = -3/2* beta; 19 b2 =-( beta -1); 20 b3 =3; 21 b4 = -2*B; 22 k=( sqrt(+ b1 ^2 -2* b1* b4 +4* b2 *b3 +b4 ^2) ); 23 C1 = -3/k *(( Ti -1) / epsilon - deltaT0 ); 24 P=@( chimax )( b3 * deltaT0_p ( chimax ))/k; 74 del cálculo variacional a la teoría del control óptimo: aplicación a sistemas lejos del eqilibrio 62 axis([0 1 0.2 1]) 63 xlabel(’$\ alpha$ ’,’interpreter’,’latex ’,’Fontsize ’ ,16) 64 ylabel(’$p_1$ ’,’interpreter’,’latex ’,’Fontsize ’ ,16) 65 subplot(1 ,2 ,2) 66 box on 67 plot(alfav ,param2 ,’r’,alfav ,0.5* ones (1,n),’b’,’LineWidth ’ ,2) 68 axis([0 1 0.5 0.52]) 69 xlabel(’$\ alpha$ ’,’interpreter’,’latex ’,’Fontsize ’ ,16) 70 ylabel(’$p_2$ ’,’interpreter’,’latex ’,’Fontsize ’ ,16) 71 figure(2) 72 nume=23; %nos permite seleccionar el valor de alfa deseado 73 unos=ones(1 , length(Tiv)); 74 hold on 75 box on 76 plot(( Tiv - unos )*10^3 , tf(nume ,:) ,’k*’,’MarkerSize ’ ,10, ’LineWidth ’ ,2) 77 plot(( Tiv - unos ) *10^3 , (B( alfav ( nume ))) ^( -1/2) *( Tiv -1) .^(1/2) ,’b -’ ,(Tiv - unos ) *10^3 , param1 ( nume )*(Tiv -1) .^( param2 ( nume )) ,’r-’,’ LineWidth ’ ,2) 78 axis([- inf inf -inf inf ]) 79 xlabel(’$(T_i -1) 10^3$’,’interpreter’,’latex ’,’ Fontsize ’ ,16) 80 ylabel(’$t_f$ ’,’interpreter’,’latex ’,’Fontsize ’ ,16) 1function [f] = funcion_error_CH (Ti ,tf ,p) 2p1=p (1); 3p2=p (2); 4%% Aproximamos por funci ón teó rica 5tpf =p1 .* abs (1 -Ti ).^ p2; 6%% funci ón a optimizar 7f=tpf -tf; 8end Bibliografía [1] C. Lanczos, The variational principles of mechanics (University of Toronto Press, 1949). [2] D. Liberzon Calculus of variation and optimal control theory (Princeton University Press,2012). [3] I. M. Gelfand and S. V. Fomin, Calculus of variations (Prentice-Hall INC, 1963). [4] L. S Pontryagin, Mathematical theory of optimal processes (CRC Press, 1987). [5] C. A. Plata, D. Guéry-Odelin, E. Trizac and A. Prados, Finite time adiabatic processes: derivation and speed limit, Physical Review E 101, 032129 (2020). [6] H. M. Jaeger, Sidney R. Nagel and Robert O. Behringer, Granular solids, liquids and gases, Reviews of Modern Physics 68, 1259 (1996). [7] G. Y. Onoda and E. G. Liniger, Random loose packings of uniform spheres and the dilatancy onset, Physical Review Letters 64, 2727 (1990). [8] A. Rosatp, K. J. Strandburg, F. Prinz and R. H. Swendsen, Why the Brazil Nut are on the top: size segragation of particulate matter by shaking, Physical Review Letters 8, 1038 (1987). [9] Matthias E. Möbius, Benjamin E. Laudarle, Sidney R. Nagel, Heinrich M. Jaeger, Size separation of granular particles, Nature 414, 270 (2001). [10] N. V. Brilliantov and T. Pöschel, Kinetic theory of granular gases (Oxford University Press, 2004). [11] T.P.C Van Noije and M. H. Ernst, Velocity distributions in homogeneous granular fluids: the free and heated case, Granular Matter 1, 57 (1998). [12] J. M. Montanero and A. Santos, Computer simulation of uniformly heated granular fluids, Granular Matter 2, 53 (2000). [13] E. Trizac and A. Prados, Memory effect in uniformly heated granular gases, Physical Review Letters 112, 198001 (2014). 76 del cálculo variacional a la teoría del control óptimo: aplicación a sistemas lejos del eqilibrio [14] A. Prados, Optimizing the relaxation route with optimal control, Physical Review Research 3, 023128 (2021).