Aplicación de redes neuronales al ajuste de parámetros en modelos diferenciales

TRABAJO DE FIN DE GRADO Doble grado en f´ ısica y matem´ aticas Aplicaci´ on de redes neuronales al ajuste de par´ ametros en modelos diferenciales Autor: Sergio Murillo Garc´ıa Tutores: Francisco Manuel Guill´en Gonz´alez Anna Doubova Krasotchenko Junio 2023 Por todos las personas que han hecho que estos cinco a˜nos sean tan bonitos. Resumen El ´exito del desarrollo de la inteligencia artificial durante los ´ultimos 15 a˜nos ha revolucionado todas las ´areas del conocimiento. Debido a la ((maldici´on de la dimensi´on)), que hace referencia al desaf´ıo que surge al tratar de resolver sistemas de ecuaciones diferenciales con un gran n´umero de variables, su aplicaci´on a la aproximaci´on num´erica de ecuaciones diferenciales juega un papel muy importante en problemas de alta dimensi´on. En este trabajo se realizar´a una introducci´on a los conceptos b´asicos de las redes neuronales generales. Se estudiar´a el algoritmo de diferenciaci´on autom´atica (backpropagation), clave de su eficiencia, y se demostrar´a el Teorema de aproximaci´on universal, base matem´atica que justifica la utilidad de las redes neuronales en la aproximaci´on de funciones. Se utilizar´an las redes neuronales denominadas como Physical Informed Neural Networks (PINN) para la resoluci´on de ecuaciones y sistemas diferenciales y el ajuste de par´ametros en ellos. Adem´as, se implementar´an utilizando Python y sus librer´ıas en problemas concretos. Palabras clave: redes neuronales, entrenamiento, hiperpar´ametros, derivaci´on autom´atica, ajuste de par´ametros, sistema diferencial, python. Abstract The successful development of artificial intelligence over the last 15 years has revolutionised all areas of knowledge. Due to the ((curse of dimensionality)), which refers to the challenge that arises when trying to solve systems of differential equations with a large number of variables, its application to the numerical approximation of differential equations plays a very important role in high-dimensional problems. In this work, an introduction to the basic concepts of general neural networks will be given. The automatic differentiation algorithm (backpropagation), the key to its efficiency, will be studied and the universal approximation theorem, the mathematical basis that justifies the usefulness of neural networks in the approximation of functions, will be demonstrated. Neural networks known as Physically Informed Neural Networks (PINN) will be used to solve differential equations and systems and to adjust parameters in them. Furthermore, they will be implemented using Python and its libraries in concrete problems. Keywords: neural networks, training, hyperparameters, automatic derivation, parameter tuning, differential system, python. ´ Indice general 1. Introducci´on 8 2. Redes neuronales 10 2.1. Entrenamiento .............................. 12 2.2. Hiperpar´ametros ............................. 13 2.2.1. Funci´on de activaci´on . . . . . . . . . . . . . . . . . . . . . . . 13 2.2.2. Inicializaci´on de par´ametros . . . . . . . . . . . . . . . . . . . 16 2.2.3. Funci´on de p´erdida . . . . . . . . . . . . . . . . . . . . . . . . 16 2.2.4. Datos de entrenamiento . . . . . . . . . . . . . . . . . . . . . 17 2.2.5. N´umero de neuronas y capas . . . . . . . . . . . . . . . . . . . 18 2.2.6. N´umero de ´epocas . . . . . . . . . . . . . . . . . . . . . . . . 19 2.2.7. Optimizador............................ 20 2.2.8. Factor de entrenamiento (learning rate)............. 23 3. Derivaci´on autom´atica: Backpropagation 26 4. Teorema de aproximaci´on universal 32 5. Resoluci´on de problemas diferenciales 36 5.1. PINN ................................... 36 5.2. Problema directo para una EDO . . . . . . . . . . . . . . . . . . . . . 39 5.2.1. Estructura de la red . . . . . . . . . . . . . . . . . . . . . . . 41 5.2.2. N´umero de ´epocas . . . . . . . . . . . . . . . . . . . . . . . . 43 5.2.3. Factor de entrenamiento (learning rate)............. 44 5.2.4. Optimizador............................ 45 Sergio Murillo Garc´ıa ´ INDICE GENERAL 5.2.5. Funci´on de activaci´on . . . . . . . . . . . . . . . . . . . . . . . 46 5.3. Sistema diferencial . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 6. PINN para el ajuste de par´ametros 50 6.1. Obtenci´on de par´ametros a posteriori . . . . . . . . . . . . . . . . . . 50 6.2. Ajuste de par´ametros discretos . . . . . . . . . . . . . . . . . . . . . . 51 6.2.1. Ejemplo de ajuste de par´ametros discretos para una EDO . . 52 6.2.2. Ejemplo de ajuste de par´ametros discretos para un sistema de EDO’s............................... 55 6.3. Ajuste de par´ametros continuos . . . . . . . . . . . . . . . . . . . . . 59 6.3.1. Ejemplo ((educativo)) de ajuste de par´ametros continuos . . . . 59 6.3.2. Ejemplo ((real)) de ajuste de par´ametros continuos . . . . . . . 63 7. Ecuaciones en derivadas parciales 68 7.1. Problemadirecto ............................. 68 7.2. Ajuste de par´ametros . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 8. Herramientas inform´aticas y c´odigos 74 8.1. Problemadirecto ............................. 75 8.2. Ajuste de par´ametros discretos . . . . . . . . . . . . . . . . . . . . . . 80 8.3. Ajuste de par´ametros continuos . . . . . . . . . . . . . . . . . . . . . 83 7 Cap´ıtulo 1 Introducci´on El auge de la inteligencia artificial en los ´ultimos tiempos es algo irrefutable. Con la creaci´on de la plataforma ((ChatGPT)) a finales de 2022 (ver [1]), cualquier persona puede acceder a la potencia de la inteligencia artificial (IA) con un simple clic. Sin embargo, la IA no es algo que se haya desarrollado en cuesti´on de un a˜no. La expresi´on ((inteligencia artificial)) fue acu˜nada en 1956 en la Conferencia de Dartmouth por el inform´atico estadounidense John McCarthy. Se entiende por IA un conjunto de capacidades cognoscitivas e intelectuales expresados por sistemas inform´aticos cuyo prop´osito es imitar la inteligencia humana para realizar tareas, y que pueden mejorar a medida que recopilan informaci´on. Anteriormente, los estudios publicados en 1936 por parte del ingl´es Alan Turing formalizaron el concepto de ((algoritmo)) y establecieron las bases te´oricas de la computaci´on, es por ello que se le considera ((el padre)) de la computaci´on moderna (ver [2]). El desarrollo de la inteligencia artificial ha ido de la mano de cuestiones filos´oficas y ´eticas. En 1950, Alan Touring propone el conocido como ((Test de Touring)) cuyo objetivo era diferenciar si un interlocutor era una inteligencia artificial o un humano (ver [3]). El programa ((Eliza)) es considerado como el primer programa inform´atico que pasa el test, fue desarrollado en 1966 por el Alem´an Joseph Weizenbaum (ver [4]). En la historia de la inteligencia artificial se distinguen 3 periodos principales: la ((cibern´etica)) en 1940-1960, ((conexionismo)) en 1980-1990, y la corriente actual, bajo el nombre de ((aprendizaje profundo)) (Deep Learning). Los algoritmos de inteligencia artificial surgen con la idea de imitar c´omo el cerebro es capaz de aprender y procesar la informaci´on, es por ello que uno de los nombres que recibe la IA es el de ((Redes neuronales artificiales)) (Artificial Neuronal Networks). Desde el desarrollo por parte de Frank Rosenblatt de la primera red neuronal en 1959, conocida como ((perceptr´on)), la complejidad de las mismas ha ido en un aumento sin l´ımites. Desde el ((perceptr´on)), que consta de una sola neurona intermedia, hasta llegar a tener m´as de un mill´on de neuronas (para hacerse una idea del orden de magnitud, el ser humano tiene del orden de 1011 neuronas). Adem´as, se estima que el n´umero de neuronas intermedias se duplica cada dos a˜nos y medio (ver [5]). Sergio Murillo Garc´ıa En el ´ambito de la resoluci´on num´erica de las ecuaciones diferenciales, los m´etodos cl´asicos m´as importantes son los m´etodos de diferencias finitas y elementos finitos. Estos m´etodos discretizan las ecuaciones para aproximar la soluci´on, cuando el n´umero de variables aumenta, el n´umero de puntos de la cuadr´ıcula de discretizaci´on aumenta de forma exponencial. Esto se conoce como la ((maldici´on de la dimensi´on)) (Curse of dimensionality) y hace que la resoluci´on de ecuaciones diferenciales en altas dimensiones sea ineficiente. La aplicaci´on de la inteligencia artificial a la aproximaci´on num´erica de ecuaciones diferenciales no consiste, por tanto, en tratar de encontrar m´etodos que superen los m´etodos cl´asicos en dimensiones bajas, es decir, con pocas variables, sino tratar de resolver ecuaciones diferenciales de alta dimensi´on, donde los m´etodos cl´asicos presentan dificultades. En el Cap´ıtulo 2 se plantean las bases del funcionamiento de las redes neuronales artificiales y sus principales elementos: entrenamiento, hiperpar´ametros. . . En el Cap´ıtulo 3 se introduce el proceso de derivaci´on autom´atica utilizado en las redes neuronales, que se denomina algoritmo de backpropagation. El Cap´ıtulo 4 est´a dedicado al Teorema de aproximaci´on universal, que sirve como fundamento matem´atico para asegurar que una red neuronal puede aproximar funciones bajo determinadas condiciones. Una vez realizado el desarrollo te´orico, en los siguientes cap´ıtulos se implementar´a en Python algoritmos de redes neuronales aplicados a problemas concretos. En el Cap´ıtulo 5 se introducen las redes neuronales utilizadas para la resoluci´on de ecuaciones y sistemas diferenciales (PINN), y se realiza un estudio de c´omo afectan los diferentes hiperpar´ametros en la aproximaci´on obtenida. En el Cap´ıtulo 6 se utilizan las redes neuronales para el ajuste de par´ametros en ecuaciones y sistemas diferenciales ordinarios. Las redes neuronales pueden adaptarse de forma sencilla a la resoluci´on de ecuaciones en derivadas parciales, estos algoritmos se implementan en el Cap´ıtulo 7. Por ´ultimo, el Cap´ıtulo 8 est´a dedicado a las herramientas inform´aticas y los c´odigos utilizados en el proyecto. 9 CAP´ ITULO 2. REDES NEURONALES Figura 2.6: Funci´on ReLU. La utilizaci´on de una funci´on de activaci´on u otra depende del problema concreto. En redes neuronales feed forward (aquellas en las que las conexiones de las neuronas no forman ciclos) se utilizan, en general, la sigmoide, tanh o ReLU. 2.2.2. Inicializaci´on de par´ametros La inicializaci´on de los pesos y sesgos es importante en el aprendizaje de la red. En general, se realiza de forma aleatoria, ya que en caso de que dos neuronas tuviesen iguales par´ametros, ambas se optimizar´ıan siempre de igual forma (tendr´ıan el mismo gradiente asociado) y, por tanto, se estar´ıa perdiendo poder de mejora. Deben ser n´umeros no excesivamente grandes, ni excesivamente peque˜nos, ya que aparecer´ıa lo que se conoce como desvanecimiento del gradiente (vanishing gradient) y gradiente explosivo (exploding gradients), respectivamente (ver [6]). La inicializaci´on m´as extendida es la denominada Xavier normal oGlorot normal. Consideremos las neuronas de una determinada capa, y llamemos fin yfout el n´umero de entradas y salidas de dicha capa. Se toma entonces, para los valores iniciales de los pesos y bias de las neuronas de esta capa, n´umeros aleatorios siguiendo una distribuci´on normal, acotados entre ±q6 fin+fout (ver [13]). 2.2.3. Funci´on de p´erdida La funci´on de p´erdida es la que se utiliza en el proceso de optimizaci´on y entrenamiento una vez se ha realizado una etapa del entrenamiento de la red. En la literatura tambi´en se denomina funci´on de coste, objetivo o de error. La determinaci´on de esta funci´on depende mucho del problema concreto que se est´e tratando y de la funci´on de activaci´on que se utilice. A continuaci´on pasamos a describir algunas de las m´as comunes. Error cuadr´atico medio El mean squared error o((MSE)) es el m´as usado en los modelos de regresi´on. Su expresi´on es 16 Sergio Murillo Garc´ıa 2.2. HIPERPAR´ AMETROS MSE =1 n n X i=1 (y(xi)−yi)2, donde los (xi, yi) son los ndatos utilizados en el entrenamiento (yies la salida asociada a la entrada xi) e y(x) representa la funci´on de la red neuronal con una entrada y una salida. En el caso de que la red tuviese jentradas, xipasar´ıa a ser un vector con jcomponentes e igual para yien una red neuronal con m´as de una neurona en la capa de salida. Error absoluto medio Tambi´en puede usarse el mean absolut error o((MAE)) cuya expresi´on es MAE =1 n n X i=1 |y(xi)−yi|, donde los (xi, yi) son los ndatos utilizados en el entrenamiento e y(x) representa la funci´on de la red neuronal con una entrada y una salida. De nuevo, en el caso de que la red tuviese jentradas, xipasar´ıa a ser un vector con j componentes y an´alogamente para yien una red con m´as de una salida. Existen otras muchas funciones de p´erdida dependiendo del objetivo de la red, por ejemplo, categorical cross-entropy osparse categorical cross-entropy usadas en el contexto de la clasificaci´on de im´agenes [12]. 2.2.4. Datos de entrenamiento Ya se ha mencionado que para el proceso de entrenamiento se necesitan datos de nuestro problema espec´ıfico, estos son el carburante de la red. Por ejemplo, si queremos un sistema que distinga entre fotos de gatos y perros, necesitaremos una gran cantidad de im´agenes de ambos animales para que el sistema ((aprenda)) a diferenciarlos durante el entrenamiento. El desaf´ıo del Machine Learning es predecir de forma acertada el comportamiento del sistema para datos que la red no ha visto anteriormente y no solo para los datos usados en el entrenamiento, para los que, tras dicho proceso, el error habr´a sido minimizado. Esta habilidad de la red neuronal para predecir de forma correcta para datos ((desconocidos)) para la red se denomina generalization o((generalizaci´on)). Normalmente se dividen los datos disponibles en 2 grupos: ((datos de entrenamiento)) y((datos de prueba)). Los datos de entrenamiento, que suponen entorno al 80 %, son los usados en el proceso de entrenamiento de la red para ajustar par´ametros y elegir adecuadamente los hiperpar´ametros, mientras que los datos de prueba se utilizan para comprobar la capacidad de generalizaci´on de la red (ver [5]). Se denominan training error ytest error el error de los datos de entrenamiento y los datos de test, respectivamente. En general, se espera que el error del conjunto de test sea superior al de entrenamiento. Se tienen entonces dos objetivos principales: 17 CAP´ ITULO 2. REDES NEURONALES 1. Minimizar el training error. 2. Minimizar la diferencia entre el training error y el test error. Estos dos objetivos se corresponden con los dos problemas centrales del Machine Learning, conocidos como ((subajuste)) y((sobreajuste)) (underfitting yoverfitting). El underfitting es no poder obtener un error de entrenamiento suficientemente bajo y el overfitting ocurre cuando la diferencia entre ambos errores es grande, es decir, cuando la red no es capaz de predecir adecuadamente valores no vistos. Se define en este contexto la ((capacidad)) de la red como la habilidad de la misma para ajustar m´as o menos variedad de funciones. Por ejemplo, una red con una alta capacidad podr´ıa ajustar polinomios de grado 20, mientras que una red con poca capacidad ajustar´ıa hasta polinomios lineales. Naturalmente, esto depender´a del n´umero de par´ametros entrenables de la red, aspecto que se tratar´a con mas detalle en el pr´oximo apartado. Por ´ultimo, observamos que es interesante saber que, en ocasiones, se trabaja con cantidades enormes de datos, es por ello que el ordenador puede tener problemas para trabajar con todos ellos al mismo tiempo en la memoria. Se realiza entonces el entrenamiento tomando ((lotes)) (batches) del total de datos de entrenamiento. El tama˜no ´optimo depende de muchos factores, entre ellos, de la memoria del ordenador, como ya hemos comentado. 2.2.5. N´umero de neuronas y capas El n´umero de par´ametros internos de la red neuronal se ve reflejado en el n´umero de neuronas y capas. En el apartado anterior se han introducido los conceptos de underfitting yoverfitting. Veamos un ejemplo para poner de manifiesto qu´e puede ocurrir con la red neuronal cuando se tienen estos fen´omenos. Supongamos que tenemos una red cuya salida son polinomios de grado 9 y otra red cuya salida es una funci´on lineal. Supongamos tambi´en que la soluci´on exacta del problema es una funci´on cuadr´atica (el ajuste buscado) y que contamos con 6 datos de entrenamiento (valores que satisfacen la soluci´on cuadr´atica). Es claro que la red lineal no podr´a ajustar bien los datos. Mientras que la red de grado 9 s´ı que es capaz de representarlos de forma adecuada, pero adem´as, es capaz de representar otras muchas funciones que tambi´en pasan de forma exacta por los puntos, ya que hay m´as par´ametros que datos. Esto se puede observar mejor en la Figura 2.7 [5]. Como se observa en dicha figura, para la red lineal el error de entrenamiento es grande, pues la red no es capaz de ajustar los datos, se tiene el underfitting. En la imagen de la derecha de la Figura 2.7, se puede ver que la red de grado 9 ajusta de forma exacta los datos de entrenamiento, es decir, se tiene error de entrenamiento 0. Pero es claro que el error de test ser´a muy grande, pues la curva no sigue una ecuaci´on cuadr´atica. Luego, no tendr´a buena capacidad de generalizar, se da el overfitting. En este caso la capacidad apropiada ser´ıa la del centro, donde no hay un n´umero de par´ametros excesivos, ni demasiado pocos. 18 Sergio Murillo Garc´ıa 2.2. HIPERPAR´ AMETROS Figura 2.7: Ejemplo de underfitting (izquierda, red lineal); aproximaci´on adecuada (centro, red de grado 2) y overfitting (derecha, red de grado 9). Los puntos son los datos de entrenamiento y la curva la salida de la red [5]. Figura 2.8: Relaci´on t´ıpica entre los errores de entrenamiento y de test (generalization error) y la capacidad [5]. En conclusi´on, el n´umero de par´ametros no puede ser demasiado elevado, ni demasiado bajo. Se debe ajustar en relaci´on a la cantidad de datos de los que se disponga. En la Figura 2.8 se puede ver la relaci´on habitual existente en las redes para la capacidad ´optima de un problema dado. 2.2.6. N´umero de ´epocas Como ya se ha visto, el n´umero de ´epocas o epochs es el n´umero de veces que los datos de entrenamiento pasan por la red y se actualizan los par´ametros. En otras palabras, es el n´umero de veces que se resuelve el problema de optimizaci´on asociada a la funci´on de p´erdida. En general, el error de los datos de entrenamiento disminuye con las ´epocas, pero llegados a un punto, la mejora puede no ser rentable comparado con el tiempo necesario. Por otro lado, un n´umero de ´epocas por debajo del ´optimo puede limitar el potencial del modelo. En general, el n´umero de ´epocas adecuado se decide haciendo diferentes pruebas, en funci´on del tiempo y de los recursos computacionales que se tenga. 19 CAP´ ITULO 2. REDES NEURONALES 2.2.7. Optimizador Se ha introducido anteriormente que el optimizador ser´a el encargado de minimizar la funci´on de coste, encontrando los valores adecuados para los par´ametros internos del modelo (pesos y sesgos). Tambi´en, se ha mencionado que estos optimizadores se basan en el algoritmo del descenso del gradiente: cada par´ametro interno se actualiza rest´andole la derivada de la funci´on de error respecto de dicho par´ametro multiplicado por un factor que se denomina ((factor de entrenamiento)) (learning rate) del que hablaremos en el apartado siguiente. Existen variantes del algoritmo de descenso que tratan de mejorar el rendimiento, pues el descenso de gradiente b´asico, tiene una convergencia lenta, en general, y costosa num´ericamente. Antes de ver algunos ejemplos de los optimizadores m´as usuales, veamos una de las variantes m´as importantes al m´etodo b´asico: el m´etodo del momento (momentum). Figura 2.9: Esquema de la mejora del m´etodo del gradiente con el uso del momentum en el descenso al m´ınimo global, evitando el m´ınimo local [14]. El momentum o((cantidad de movimiento)) hace referencia al momento lineal de la F´ısica, cuanto m´as velocidad y m´as masa tenga un cuerpo, m´as ((movimiento)) tendr´a, es decir, m´as dif´ıcil ser´a frenarlo. En el contexto de la optimizaci´on de redes neuronales, el momentum consiste en considerar no solo el gradiente obtenido para una observaci´on, sino un promedio de los gradientes anteriores, de forma que se tiene esta noci´on de ´ımpetu y de ((cantidad de movimiento)) [12]. De esta manera, para actualizar uno de los par´ametros internos, se tiene en cuenta la ((historia)) del gradiente, el promedio de los gradientes calculados anteriormente, adem´as del gradiente actual: ω=ωanterior −α∂C ∂ω +mom ∂C ∂ω Donde ωes uno de los pesos, αes el paso del m´etodo de descenso (learning rate), C hace referencia a la funci´on de coste, ∂C ∂ω hace referencia al promedio de los gradientes anteriores y mom ∈[0,1] es el coeficiente que mide la ((importancia)) que se le da al momento (si se toma como 0 ser´ıa el gradiente habitual). 20 Sergio Murillo Garc´ıa 2.2. HIPERPAR´ AMETROS Hasta ahora hemos hablado del descenso del gradiente pensando ´unicamente en un m´ınimo global, es decir, sin tener en cuenta posibles m´ınimos locales que ((atasquen)) la optimizaci´on. Los problemas reales son, en general, complejos y existen m´ınimos locales, puesto que la funci´on de p´erdida no es convexa en general. La introducci´on del momentum puede solventarlos, como puede observarse en la Figura 2.9. Veamos a continuaci´on las caracter´ısticas principales de algunos de los optimizadores m´as usuales (ver [14, 15, 31]). Stochastic Gradient Descent (SGD): Dada la gran cantidad de par´ametros internos, calcular todas y cada una de las derivadas del error respecto a cada par´ametro y evaluarlas para todas las observaciones es inviable. Una primera aproximaci´on consiste en evaluar las derivadas ´unicamente para una observaci´on en cada ((lote)) (batch), haciendo el gradiente del funcional objetivo reducido a dicho batch. Este m´etodo se utiliza en ocasiones implementando tambi´en el m´etodo del momentum. Adaptive Gradient Algorithm (AdaGrad): La innovaci´on que introduce este optimizador es considerar que el factor de entrenamiento (learning rate), no es constanste para cada par´ametro interno, se actualiza de la siguiente forma: α=α 1+(epoch −1)µ, donde µse denomina learning rate decay (ver Secci´on 2.2.8), αes el paso del m´etodo de descenso constante y epoch es la ´epoca. Adem´as se multiplica el learning rate por un factor que escala seg´un la magnitud del gradiente. AdaDelta: Este m´etodo se basa en el AdaGrad, pero considerando ´unicamente los ngradientes anteriores, y no todos. Root Mean Square Propagation (RMSProp): Este optimizador es similar a AdaDelta. La diferencia est´a en la forma de asignar el factor de entrenamiento a cada par´ametro interno, de forma que se le da mayor prioridad a los ´ultimos valores del gradiente a˜nadidos al promedio que a los m´as antiguos. Adaptive moment estimation (Adam): Este conocido optimizador combina las bondades de AdaGrad y RMSProp y adem´as hace uso del momentum explicado anteriormente. AdaMax (Adamax): Este optimizador es una variante del m´etodo Adam, pero en vez de usar la norma 2 se usa la norma infinito en el c´alculo del momento. 21 CAP´ ITULO 2. REDES NEURONALES Como se ha visto, en general, los optimizadores se basan en optimizadores ya creados anteriormente introduciendo algunas variantes, para que su rendimiento sea mayor. En la pr´actica, el mejor optimizador depende del problema concreto con el que se est´e tratando. Adem´as, se hace notar que, en todos los casos, son optimizadores para calcular m´ınimos sin restricciones. (a) Clasificador binario de textos. (b) Clasificador multi-clase de im´agenes. (c) Pron´ostico de series temporales (Forecasting). (d) Regresi´on num´erica simple. Figura 2.10: Error de entrenamiento (loss) y error de test (val loss) con varios optimizadores para 4 tipos de problemas cl´asicos en redes neuronales [14]. En la Figura 2.10 se adjuntan los resultados de los errores (de entrenamiento y de test) obtenidos con diferentes optimizadores (los disponibles en la libreria ((Keras)) de Python) aplicados a diferentes problemas. Se hace notar que en todos ellos se usa un learning rate de 0.001 y 50 ´epocas. Se puede observar como, dependiendo del experimento concreto, el optimizador m´as adecuado var´ıa. Luego, la conclusi´on general es que la decisi´on de qu´e optimizador elegir depende del objetivo de la red. Como un primer intento, se puede utilizar el optimizador Adam, pues se observa que tiene una actuaci´on adecuada en los 4 problemas de la Figura 2.10. 22 Sergio Murillo Garc´ıa 2.2. HIPERPAR´ AMETROS 2.2.8. Factor de entrenamiento (learning rate) Seg´un lo visto anteriormente, el learning rate,((factor de entrenamiento)) o((paso)) es el factor por el que se multiplica la magnitud del gradiente para actualizar los par´ametros internos del modelo. El valor adecuado del learning rate depende mucho del problema en cuesti´on. Si los pasos son demasiado peque˜nos, el entrenamiento puede ser demasiado lento, adem´as de aumentar la probabilidad de ((atascarse)) en un m´ınimo local. Por el contrario, si el ((factor de entrenamiento)) es demasiado grande, al comienzo puede ser positivo, pues la red aprender´a r´apidamente, pero en ubicaciones cercanas al m´ınimo, puede ser que las actualizaciones lleven a los par´ametros a ubicaciones aleatorias sin llegar al m´ınimo global, pudiendo incluso divergir (ver [12]). En la Figura 2.11 se observan los comportamientos descritos. Figura 2.11: Esquema de la importancia del learning rate en el proceso de entrenamiento. L(ω) hace referencia a la funci´on de p´erdida y ωa los par´ametros internos. En la parte superior, distintos tipos de descensos seg´un el tama˜no del paso. En la gr´afica central se observa el comportamiento de convergencia de L(ω) en funci´on de las ´epocas, para cada uno de los 4 tama˜nos de las gr´aficas superiores [16]. Se puede ver en la Figura 2.11 (gr´afica central inferior) como, un factor de entrenamiento alto, puede ser ´util al comienzo, ya que permite acercarse r´apido a entornos del m´ınimo, una vez all´ı, es mejor un factor m´as peque˜no, pues nos permite alcanzar la convergencia a valores ´optimos. Por ello, un buen m´etodo es utilizar un factor de entrenamiento variable, de forma que disminuya a medida que aumentan las ´epocas. Esto es lo que se conoce como ((decaimiento del factor de entrenamiento)) (learning rate decay). Existen muchos m´etodos para implementarlo, algunos de los m´as comunes son los siguientes (ver [17]). 23 CAP´ ITULO 2. REDES NEURONALES M´etodo b´asico: El factor de entrenamiento αvar´ıa seg´un la ´epoca de la siguiente manera: α(epoch) = α0 1 + decayRate ·epoch, donde α0∈Res el factor inicial y decayRate ∈Res el learning rate decay (suele tomarse entre 0 y 1). De esta forma el factor disminuye a medida que se aumentan las ´epocas. Por ejemplo, para α0= 0.2 y tomando decayRate = 0.2, se tiene un factor de entrenamiento de 0.001 en la ´epoca 1000 y de 0.0003 en la ´epoca 3000. Decaimiento exponencial: En este caso el decaimiento es exponencial dependiendo de las ´epocas siguiendo la expresi´on: α(epoch) = decayRateepoch ·α0. Decaimiento manual: En este caso se decide de antemano el decaimiento seg´un las ´epocas. Un ejemplo podr´ıa ser: α(epoch) =    0.01 si epoch ≤1000, 0.001 si 1000 <epoch ≤3000, 0.0005 si epoch >3000, pudi´endose ajustar seg´un las ´epocas que vayan a realizarse y la velocidad de convergencia del problema concreto. 24 Sergio Murillo Garc´ıa 2.2. HIPERPAR´ AMETROS 25 Cap´ıtulo 4 Teorema de aproximaci´on universal El objetivo del trabajo es la utilizaci´on de las redes neuronales en la aproximaci´on de funciones, es por ello que se necesita saber de antemano que es posible acercarse tanto como se quiera a una funci´on dada a trav´es de la funci´on de una red neuronal. En este cap´ıtulo se dar´a una prueba de un teorema que lo asegura. Se presentan algunas definiciones previas (ver [7]). Se denominar´a Inal cubo unidad de dimensi´on n: [0,1]n. El espacio de funciones continuas en Inse denota por C(In). Se usar´a la norma del m´aximo y se denotar´a por ∥·∥ dicha norma de una funci´on en su dominio. Se denotar´a M(In) al conjunto de medidas Borel finitas. Definici´on 4.1. Se dice que una funci´on σes discriminadora si para una medida µ∈M(In)se tiene que para todo y∈Rnyθ∈R, la igualdad ZIn σ(yTx+θ)dµ(x) = 0 implica que µ= 0. Definici´on 4.2. Se dice que una funci´on σ:R→Res sigmoide si cumple σ(t)→1si t →+∞ 0si t → −∞ A continuaci´on, se presentan dos resultados que permitir´an demostrar el Teorema de Aproximaci´on Universal. Teorema 4.1. Dada una funci´on σ:R→Rcontinua y discriminadora, las siguientes funciones G(x) = N X i=1 αiσ(yT ix+θi) (4.1) Sergio Murillo Garc´ıa son densas en C(In)respecto de la norma del m´aximo (con N > 1natural). En otras palabras, para cualquier f∈C(In)yε > 0∃α, y, θ ∈RNtales que G(x) = G(x;α, y, θ)con la forma anterior satisface |G(x)−f(x)|< ε ∀x∈In. Demostraci´on. Sea S⊂C(In) el conjunto de funciones G(x) de la forma (4.1). Se tiene que Ses un subespacio vectorial de C(In). Por definici´on, Ses denso en C(In) si S=C(In), donde Ses la clausura de S. Se har´a la demostraci´on por reducci´on al absurdo. Sea R=Stal que R=C(In). Es claro que, como Res la clausura de S,Res un subespacio cerrado de C(In). Por el Teorema de Hahn Banach (ver [8], Corolario 1.8) se deduce la existencia de un funcional lineal acotado L∈L1(Rn) no nulo que cumple L(S) = L(R) = 0. Por el Teorema de representaci´on de Riesz se puede escribir el funcional anterior como L(h) = ZIn h(x)dµ(x)∀h∈C(In) para alg´un µ∈M(In). En particular, como σ(yTx+θ)∈R∀y, θ, se tiene ZIn σ(yTx+θ)dµ(x) = 0 ∀y, θ. Como σes discriminadora, de la igualdad anterior se obtiene µ= 0, que implica L id´enticamente nulo, luego se llega a una contradicci´on. Se tiene entonces S=C(In) y Sdenso en C(In). ■ A continuaci´on, se presenta un resultado que muestra que cualquier funci´on sigmoide continua es discriminadora. Cabe se˜nalar que, en las aplicaciones de redes neuronales, las funciones de activaci´on sigmoide continuas se toman mon´otonas crecientes, pero no se requiere la monoton´ıa en los resultados presentados aqu´ı. Lema 4.1. Toda funci´on sigmoide σacotada y medible es discriminadora. En particular, toda funci´on continua sigmoide es discriminadora. Demostraci´on. Sea σfunci´on sigmoide. Se observa que ∀x, y ∈RNyθ, φ ∈Rse tiene: σ(λ(yTx+θ) + φ)   →1 si yTx+θ > 0 cuando λ→+∞, →0 si yTx+θ < 0 cuando λ→+∞, =σ(φ) si yTx+θ= 0 para todo λ. 33 CAP´ ITULO 4. TEOREMA DE APROXIMACI´ ON UNIVERSAL Entonces, la funci´on σλ(x) = σ(λ(yTx+θ)+φ) converge puntualmente a la funci´on γcuando λ→+∞, donde γ:R→Rse define como sigue: γ(x) =    1 si yTx+θ > 0, 0 si yTx+θ > 0, σ(φ) si yTx+θ= 0. Se definen entonces el hiperplano y los semiespacios siguientes: Πy,θ =x∈RN|yTx+θ= 0, H′ y,θ =x∈RN|yTx+θ < 0, Hy,θ =x∈RN|yTx+θ > 0. Se tiene entonces que l´ımλ→∞ σλ(x) = γ(x) y adem´as la funci´on |σλ(x)|est´a acotada, pues σes sigmoide, luego por el Teorema de convergencia dominada de Lebesgue (ver [10], apartado 1.34) se tiene 0hip = l´ım λ→∞ ZIn σλ(x)dµ(x) = ZIn γ(x)dµ(x) = 1 ZHy,θ dµ(x)+0ZH′ y,θ dµ(x) + σ(φ)ZΠy,θ dµ(x) = µ(Hy,θ) + σ(φ)µ(Πy,θ), para todo φ, θ, y. Se ver´a ahora que, como la medida de todos los semiplanos es nula, la medida µdeber´a ser id´enticamente nula. Fijado y, para una funci´on medible y acotada hse define el funcional lineal Fcomo sigue: F(h) = ZIn h(yTx)dµ(x). Como µes una medida finita, Fes un funcional acotado en L∞(R). Tomemos hla funci´on indicatriz del intervalo [θ, ∞) (esto es h(u) = 1 si u≥θyh(u) = 0 si u < θ). En ese caso, se tiene: h(yTx) = 1 si yTx−θ≥0, 0 si yTx−θ < 0. En particular hes funci´on sigmoide. Entonces se tiene 0hip =F(h) = ZIn h(yTx)dµ(x) = ZHy,−θ h(yTx)dµ +ZH′ y,−θ h(yTx)dµ +ZΠy,−θ h(yTx)dµ =µ(Hy,θ) + µ(Πy,θ). Si se hace lo mismo esta vez con hla funci´on indicatriz del intervalo abierto (θ, ∞), se tiene 0hip =F(h) = ZIn h(yTx)dµ(x) = µ(Hy,θ). Luego, se deduce que F(h) = 0 para toda funci´on indicatriz de cualquier intervalo. Por linealidad, se deduce que tambi´en es cero para cualquier combinaci´on de ellas, 34 Sergio Murillo Garc´ıa es decir, para cualquier funci´on simple. Como las funciones simples son densas en L∞(R) (ver [9], apartado (2.4.13)) pues L∞⊂Lp, se deduce que Fes id´enticamente nulo. En particular, para las funciones acotadas y medibles s(u) = sen (m·u) y c(u) = cos(m·u), se puede escribir F(s+ic) = F(s)+iF(c) = ZIn cos(mTx)+isen(mTx)dµ(x) = ZIn eimTxdµ(x)=0 ∀m Es decir, la transformada de Fourier de µes 0 (ver [9], apartado 8.1.1). Luego µ debe ser 0. Luego σes discriminadora. ■ Se puede entonces probar ahora el teorema principal de esta secci´on, que afirma que cualquier funci´on continua fpuede ser aproximada por una combinaci´on lineal de funciones sigmoide. Teorema 4.2 (de aproximaci´on universal).Sea σuna funci´on sigmoide continua. Entonces las siguientes funciones G(x) = N X i=1 αiσ(yT ix+θi) son densas en C(In). Es decir, dada una funci´on f∈C(In)y dado ϵ > 0, existe G(x)de la forma anterior tal que, |G(x)−f(x)|< ϵ x ∈In. Demostraci´on. Por el Lema 4.1 se tiene que toda funci´on sigmoide continua es discriminadora y por el Teorema 4.1 se tiene lo que se quiere. ■ Existen muchas generalizaciones de este teorema utilizando otro tipo de funciones de activaci´on, ya que las funciones sigmoides no son las ´unicas que se implementan en las redes neuronales (ver [7]). 35 Cap´ıtulo 5 Redes Neuronales para la resoluci´on de problemas diferenciales: PINN Hasta ahora se han introducido los fundamentos de las redes neuronales y su funcionamiento de forma general, en este cap´ıtulo veremos c´omo se pueden aplicar a la aproximaci´on num´erica de ecuaciones diferenciales. Una de las ventajas m´as importantes de utilizar este m´etodo es que es f´acilmente adaptable a problemas diversos, lineales y no lineales, ecuaciones y sistemas diferenciales, ecuaciones y sistemas de ecuaciones en derivadas parciales, problemas inversos y problemas de control. . . Adem´as, veremos ejemplos de aplicaci´on a la resoluci´on de problemas concretos. 5.1. PINN Las redes neuronales que se usan para resolver este tipo de problemas se denominan ((Redes neuronales informadas por la f´ısica)) o PINN por sus siglas en ingl´es Physics-Informed Neural Networks. Para poder comprender mejor c´omo se adaptan los hiperpar´ametros que vimos en el Cap´ıtulo 2 a nuestro objetivo, comenzaremos aplicando las redes neuronales a un ejemplo sencillo de una ecuaci´on diferencial. Tomaremos como ejemplo la resoluci´on de una ecuaci´on diferencial sencilla de la que se conoce la soluci´on expl´ıcita, con una condici´on inicial dada. M´as concretamente, consideramos el problema de Cauchy: u′(t) = λt, t ∈[tmin, tmax], u(tmin) = u0,(5.1) donde λ, u0, tmin, tmax ∈Rson par´ametros conocidos. La soluci´on anal´ıtica es u(t) = λ(t2−t2 min) 2+u0, t ∈[tmin, tmax].(5.2) Sergio Murillo Garc´ıa 5.1. PINN La pregunta ahora es: ¿cu´ales son la funci´on de p´erdida, las entradas y salidas de la red neuronal y los datos de entrenamiento? Como ya dijimos en el Cap´ıtulo 2, la red neuronal es, en realidad, una funci´on que toma la entrada (o entradas), opera con un conjunto de par´ametros internos y nos devuelve la salida (o salidas) correspondiente. Por tanto, el n´umero de entradas y salidas de la red, es decir, el n´umero de neuronas en la capa de entrada y salida, debe ser igual al n´umero de variables independientes y dependientes de nuestro problema, respectivamente. En nuestro caso, 1 neurona de entrada (t) y 1 neurona de salida (u(t)). Veamos ahora cu´al deber´ıa ser la funci´on de coste adecuada. Denominemos ˆu(t) la salida de la red, que ser´ıa la aproximaci´on de la soluci´on exacta, u(t). El objetivo de la red es hacer que se cumpla la ecuaci´on diferencial y la condici´on inicial del problema de Cauchy (5.1). Por tanto, el funcional de coste podr´ıa ser el siguiente: C=1 Nc Nc X i=1 ∥ˆu′(ti)−λti∥2, donde ˆu′(t) es la derivada, obtenida mediante el algoritmo de backpropagation, de la funci´on de la red neuronal ˆu(t), Nces el n´umero total de tiempos de entrada, ti, usados en el entrenamiento y ∥·∥ hace referencia a la norma 2. Es decir, se toma la funci´on de coste como el error cuadr´atico medio entre ambos miembros de la ecuaci´on diferencial, que es lo que se denomina ((residuo)) de la ecuaci´on diferencial. Adem´as, como debe cumplirse la condici´on inicial, le sumamos a la funci´on de coste anterior el t´ermino correspondiente, obteniendo el funcional de coste asociado a este problema C=1 Nc Nc X i=1 ∥ˆu′(ti)−λti∥2+∥ˆu(tmin)−u0∥2.(5.3) Durante el entrenamiento de la red, se minimizar´a esta funci´on de coste de forma que se satisfaga, ((idealmente)), la ecuaci´on diferencial y la condici´on inicial. M´as concretamente, tras el proceso de entrenamiento se tendr´ıa C−→ 0⇒∥ˆu′(ti)−λti∥ −→ 0 ∥ˆu(tmin)−u0∥ −→ 0⇒ˆu′(ti)−→ λti ˆu(tmin)−→ u0para i= 1, . . . , Nc. Que es justo lo que queremos que cumpla la soluci´on. Vemos que, en principio, la red neuronal cumple el problema de Cauchy (5.1), al menos de forma aproximada, en los tiempos de entrenamiento ti, pero, si los hiperpar´ametros son adecuados, cabe esperar que se extrapole bien la soluci´on al resto de puntos del dominio. Podemos preguntarnos qui´enes son los datos de entrenamiento tien este caso. La respuesta es que, a priori, no conocemos la soluci´on exacta, luego no tenemos datos que imponer, solo sabemos que la soluci´on debe cumplir la ecuaci´on diferencial y debe valer u0en el tiempo inicial tmin. Lo ´unico en lo que debemos pensar es en qui´enes van a ser los datos de entrada con los que vamos a realizar el entrenamiento y en los que se aplicar´a la funci´on de p´erdida especificada anteriormente. En general, 37 CAP´ ITULO 5. RESOLUCI´ ON DE PROBLEMAS DIFERENCIALES tenemos dos opciones principales: la primera es tomar una partici´on del intervalo [tmin, tmax], normalmente uniforme; la segunda es tomar puntos aleatorios en todo el intervalo. Para tomar los puntos de forma aleatoria hay distintos m´etodos, algunos de los usuales son los que se describen a continuaci´on. Puntos de Sobol Es una sucesi´on casi aleatoria de n´umeros (aparentemente aleatorios, pero que provienen de un algoritmo no aleatorio). Esta sucesi´on surge con el objetivo de rellenar [0,1]sminimizando los huecos, donde ses la dimensi´on del hipercubo (ver [21]). Se puede ver una representaci´on de estos puntos para el caso de una dimensi´on en la Figura 5.1. (a) 20 primeros puntos. (b) 50 primeros puntos. Figura 5.1: Puntos de Sobol en 1 dimensi´on. Muestreo de hipercubo latino (Latin Hypercube sampling)Es un m´etodo estad´ıstico para generar una muestra de n´umeros casi aleatorios. Si tenemos una cuadr´ıcula que contiene las posiciones de una muestra dada formada por n s´ımbolos, se dice que es ((cuadrado latino)) si es una matriz n×nen la que que, en cada casilla, hay una muestra de los ns´ımbolos que aparece exactamente una vez en cada fila y columna (en el caso de dos dimensiones) (ver [22]). Naturalmente, ((hipercubo latino)) es la generalizaci´on de este concepto a un n´umero arbitrario de dimensiones. Podemos ver una representaci´on gr´afica de un muestreo de hipercubo latino en el caso de una dimensi´on en la Figura 5.2. El ejemplo considerado muestra la resoluci´on del problema directo, por tanto no necesitamos datos ((concretos)) (adem´as de los puntos de entrenamiento ti) para realizar el entrenamiento. En el caso de que se buscase resolver un problema inverso, s´ı que necesitaremos datos experimentales para determinar los par´ametros desconocidos, como veremos m´as adelante. En conclusi´on, para la aplicaci´on de las redes neuronales a la resoluci´on de ecuaciones diferenciales debemos tener en cuenta los siguientes aspectos: 38 Sergio Murillo Garc´ıa 5.2. PROBLEMA DIRECTO PARA UNA EDO (a) 20 puntos. (b) 50 puntos. Figura 5.2: Muestreo de hipercubo latino (LHS) en 1 dimensi´on. 1. El n´umero de neuronas en la capa de entrada y en la capa de salida deben ser el n´umero de variables independientes y dependientes, respectivamente. 2. La funci´on de coste debe contener el residuo de la ecuaci´on diferencial, generalmente, el error cuadr´atico medio del mismo. Adem´as, las condiciones adicionales (iniciales y/o de contorno) deben a˜nadirse con un t´ermino m´as a la funci´on de coste, como el error cuadr´atico medio del residuo de estas condiciones 3. Los datos de entrenamiento son, ahora, los datos de entrada que se usar´an en dicho proceso (en el caso del problema directo). Para generarlos, debemos tomar, bien puntos aleatorios en el dominio en el que se busque solucionar el problema, bien una partici´on uniforme del mismo. 5.2. Problema directo para una EDO Veamos un ejemplo de utilizaci´on de las redes neuronales para el caso de resoluci´on de un problema directo similar al que se explica en la secci´on precedente. Adem´as, se comparar´a el impacto que tiene la variaci´on de distintos hiperpar´ametros que hemos visto anteriormente en la aproximaci´on obtenida. Consideramos el siguiente problema de Cauchy: u′(t)=3t, t ∈[0,5], u(0) = 5.(5.4) Veamos primero la aproximaci´on obtenida con la elecci´on de par´ametros que se ha considerado adecuada. Se tienen los datos y resultados de la simulaci´on en la Tabla 5.1 y en la Figura 5.3. La inicializaci´on de los par´ametros internos de la red se realizar´a en todos los casos con la inicializaci´on ((Xavier Normal)) (ver la Secci´on 2.2.2). Adem´as, todas las gr´aficas en las que se represente la funci´on de p´erdida frente a las ´epocas, se realizar´an en escala semi-logar´ıtmica para poder apreciar bien la variaci´on de la funci´on. 39 CAP´ ITULO 5. RESOLUCI´ ON DE PROBLEMAS DIFERENCIALES Hiperpar´ametros Resultados Activaci´on Sigmoide Estructura [1,10,10,1] F. coste final 3·10−6 Optimizador Adamax ´ Epocas 20000 |ˆu(0) −u0|7·10−6 Learning rate lr= 0.01, ×0.5/(2000 ´epocas) Nºpuntos 1000 Tiempo 64s Tabla 5.1: Hiperpar´ametros para la simulaci´on con los valores adecuados. Resoluci´on del problema (5.4). Los hiperpar´ametrios recogidos en la Tabla 5.1 fueron explicados en el Cap´ıtulo 2, aunque se hacen las siguientes aclaraciones: la Estructura hace referencia a la estructura de capas y neuronas de la red, as´ı [1,10,10,1] nos indica que tiene una capa de entrada con 1 neurona, luego dos capas intermedias de 10 neuronas y una ´ultima capa de salida de 1 neurona; el Nºpuntos es el n´umero de datos de entrada tique se hacen pasar por la red en cada ´epoca del proceso de entrenamiento, ser´ıa Nccon la notaci´on de la funci´on de p´erdida (5.3); la F. coste final es la evaluaci´on de la funci´on de coste tras el proceso de entrenamiento de la red neuronal; |ˆu(0)−u0| es el error absoluto cometido por la red neuronal en la condici´on inicial; el Tiempo es el tiempo total empleado por el ordenador en el proceso de entrenamiento. En la Figura 5.3 podemos ver como el learning rate variable permite al optimizador ajustar la funci´on de coste con mucha precisi´on en un tiempo adecuado. En la gr´afica c) de dicha figura se observa como la soluci´on aproximada por la red neuronal coincide pr´acticamente con la soluci´on exacta. Adem´as, se puede comprobar este hecho teniendo en cuenta el valor de la funci´on de coste tras el entrenamiento y el error en la condici´on inicial que se ve en la Tabla 5.1. 40 Sergio Murillo Garc´ıa 5.2. PROBLEMA DIRECTO PARA UNA EDO (a) Learning rate frente a las ´epocas. (b) Funci´on de p´erdida frente a las ´epocas. (c) Soluci´on exacta (((exacta))) y aproximaci´on obtenida del problema (((PINN))). Figura 5.3: Simulaci´on con los valores adecuados e hiperpar´ametros de la Tabla 5.1. Resoluci´on del problema (5.4). Estudiemos a continuaci´on c´omo afecta la variaci´on de algunos hiperpar´ametros a la aproximaci´on obtenida por la red neuronal. 5.2.1. Estructura de la red Veamos qu´e ocurre si se var´ıa el n´umero de par´ametros internos de la red neuronal. M´as concretamente, si se cambia la estructura de la misma, mientras que el resto de hiperpar´ametros se dejan sin variaciones. En primer lugar, se utilizar´a una estructura con pocas neuronas, para poder apreciar la falta de ((capacidad)) de la red (no ser´a capaz de ajustar bien con tan pocos par´ametros), se usar´a una muy sencilla: [1,1,1], es decir, una neurona de entrada, una neurona en la capa intermedia y una neurona de salida. Vemos la soluci´on obtenida en la Figura 5.4. Se puede observar que la red ya no ajusta tan bien como lo hac´ıa antes. Ahora la funci´on de la red neuronal puede obtenerse f´acilmente: ω2tanh(ω1t+b1) + b2, y con esta expresi´on no podemos ajustar tan bien la soluci´on exacta. Al haber reducido 41 CAP´ ITULO 5. RESOLUCI´ ON DE PROBLEMAS DIFERENCIALES adaptado al sistema: C=1 Nc Nc X i=1    ˆ S′(ti) + βˆ S(ti)ˆ I(ti)   2+1 Nc Nc X i=1    ˆ I′(ti)−βˆ S(ti)ˆ I(ti) + γˆ I(ti)   2 +1 Nc Nc X i=1    ˆ R′(ti)−γˆ I(ti)   2+   ˆ S(tmin)−S0   2+   ˆ I(tmin)−I0   2+   ˆ R(tmin)−R0   2, donde ˆ S(t),ˆ I(t),ˆ R(t) son las salidas de la red en el tiempo t, soluciones aproximadas del modelo SIR (5.5) y los tiempos tison los Ncvalores utilizados en el entrenamiento de la red neuronal. Este sistema diferencial no se puede resolver de manera exacta, por tanto para comparar y comprobar el funcionamiento de la red se utilizar´a un solucionador externo ya implementado en Python, como puede ser odeint, del que hablaremos en el cap´ıtulo final dedicado a las herramientas inform´aticas. El sistema que vamos a resolver es el siguiente:            S′(t) = −0.1S(t)I(t), I′(t)=0.1S(t)I(t)−I(t), t ∈[0,3], R′(t) = I(t), S(0) = 100, I(0) = 50, R(0) = 0. (5.6) Los par´ametros de la red escogidos para la simulaci´on se recogen en la Tabla 5.2 y los resultados de la simulaci´on en la Figura 5.14. Hiperpar´ametros Resultados Activaci´on Tanh Estructura [1,10,10,3] F. coste final 2.3·10−2 Optimizador Adamax ´ Epocas 20000 ||ˆu(0) −u0|| 1.5·10−2 Learning rate lr= 0.05, ×0.99/(20 ´epocas) Nºpuntos 1000 Tiempo 138s Tabla 5.2: Hiperpar´ametros para la simulaci´on con los valores adecuados. Resoluci´on del problema (5.6). Podemos observar como la aproximaci´on obtenida por la red neuronal se ajusta bien a la esperada, obtenida con odeint. La funci´on de p´erdida es minimizada de forma adecuada. Es comprensible que el tiempo de computaci´on y el m´ınimo alcanzado del funcional de p´erdida sean mayores, ya que ahora hay 3 residuos (tenemos 3 ecuaciones diferenciales) en la funci´on de p´erdida y 3 condiciones iniciales impuestas. 48 Sergio Murillo Garc´ıa 5.3. SISTEMA DIFERENCIAL (a) Learning rate frente a las ´epocas. (b) Funci´on de p´erdida frente a las ´epocas. (c) Soluci´on aproximada por la red neuronal (colores), soluci´on aproximada dada por un solucionador externo a la red (lineas continuas negras) y poblaci´on total (discontinua). Figura 5.14: Simulaci´on con los valores adecuados e hiperpar´ametros de la Tabla 5.2. Resoluci´on del problema (5.6). 49 Cap´ıtulo 6 PINN para el ajuste de par´ametros Una vez se ha visto c´omo se pueden resolver problemas directos, tanto para ecuaciones como para sistemas diferenciales, estamos en las condiciones de poder resolver problemas de ajuste para la determinaci´on de par´ametros desconocidos. Vamos a diferenciar tres tipos de problemas de ajuste de par´ametros: Tipo 1 En estos problemas el par´ametro desconocido puede obtenerse a posteriori a partir de la soluci´on aproximada. Por ejemplo, si el par´ametro es la condici´on inicial, se puede obtener tras el entrenamiento evaluando la funci´on obtenida por la red neuronal en el tiempo inicial (ver Secci´on 6.1). Tipo 2 El par´ametro desconocido es un par´ametro discreto (un valor real) que no puede obtenerse a posteriori. Por ejemplo, los par´ametros β, γ del modelo SIR (5.5) (ver Secci´on 6.2). Tipo 3 En este caso el par´ametro es una funci´on continua. Por ejemplo, en un problema de control (ver Secci´on 6.3). Se ver´an a continuaci´on cada uno de ellos con m´as detalles. 6.1. Obtenci´on de par´ametros a posteriori En estos problemas tipo 1, podemos obtener el par´ametro en cuesti´on a partir de la soluci´on aproximada. Desde el punto de vista computacional, es un problema an´alogo al problema directo, ya que para obtener dicho par´ametro s´olo hay que evaluar la aproximaci´on obtenida por la red neuronal donde sea necesario. Por ejemplo, podemos resolver el siguiente problema, donde el par´ametro desconocido es la condici´on inicial e imponemos la condici´on final, es decir, suponemos conocida la funci´on en el tiempo final: Sergio Murillo Garc´ıa 6.2. AJUSTE DE PAR´ AMETROS DISCRETOS      u′(t) = 3t, t ∈[0,5], u(5) = 20, u(0) = λ, (6.1) siendo λ∈Rdesconocido. De acuerdo con la soluci´on exacta (5.2), el valor exacto del par´ametro que buscamos es λexa =−17.5. Los hiperpar´ametros y resultados de la simalaci´on pueden verse en la Tabla 6.1 y en la Figura 6.1. Hiperpar´ametros Resultados Activaci´on Sigmoid Estructura [1,10,10,1] F. coste final 1·10−4 Optimizador Adamax ´ Epocas 20000 |ˆu(0) −λexa|1.9·10−3 Learning rate lr= 0.1, ×0.8/(500 ´epocas) Nºpuntos 1000 Tiempo 62s Tabla 6.1: Hiperpar´ametros para la simulaci´on con los valores adecuados. Resoluci´on del problema (6.1). (a) Funci´on de p´erdida frente a las ´epocas. (b) Soluci´on aproximada por la red neuronal y soluci´on exacta. Figura 6.1: Simulaci´on con los hiperpar´ametros de la Tabla 6.1. Resoluci´on del problema (6.1). Se puede observar como la red neuronal aproxima de forma eficaz el problema con la condici´on final conocida, alcanzando valores muy peque˜nos para la funci´on de p´erdida, con una actuaci´on similar a la del problema directo. 6.2. Ajuste de par´ametros discretos En este tipo de problemas tipo 2 buscamos aproximar uno o varios par´ametros desconocidos. Hasta ahora no hemos necesitado observaciones de las variables de la ecuaci´on para obtener una aproximaci´on de la soluci´on, bastaba con imponer a la funci´on de la red neuronal que debe cumplir la ecuaci´on diferencial y las condiciones 51 CAP´ ITULO 6. PINN PARA EL AJUSTE DE PAR´ AMETROS iniciales. Ahora el problema cambia, la funci´on de p´erdida no s´olo debe incluir el residuo de las ecuaciones diferenciales y las condiciones iniciales, para poder determinar par´ametros desconocidos, necesitamos darle observaciones adicionales de la soluci´on en ciertos tiempos, de esta forma podremos obtener los par´ametros que mejor ajusten la soluci´on a las observaciones. La funci´on de p´erdida para el caso de una ecuaci´on diferencial ordinaria u′(t) = f(u, t), junto con la condici´on inicial u(tmin) = u0ser´a, por tanto, C=1 Nc Nc X i=1 ∥ˆu′(ti)−f(ˆu, ti)∥2+∥ˆu(tmin)−u0∥2+1 No No X k=1 ∥ˆu(tk)−uk∥2, donde Ncdenota el n´umero total de ((puntos de colocaci´on)) (collocations points), que son los valores de entrada utilizados para evaluar el residuo de la ecuaci´on diferencial; Noes el n´umero total de observaciones uken el intervalo de resoluci´on, es decir, los valores experimentales que sabemos (o intuimos) que sigue la soluci´on [23]. Las preguntas ahora son: ¿c´omo introducimos los par´ametros desconocidos en la red y c´omo se ajustar´an sus valores? Los par´ametros desconocidos se introducir´an como par´ametros entrenables de la red, de manera que podr´an ser ajustados, de acuerdo con el descenso de gradiente, junto con el resto de par´ametros internos. Al igual que el resto de par´ametros de la red neuronal, debemos darle un valor inicial. Veamos a continuaci´on la resoluci´on de dos ejemplos para este tipo de ajuste de par´ametros, uno aplic´andolo a una ecuaci´on diferencial ordinaria y otro a un sistema. 6.2.1. Ejemplo de ajuste de par´ametros discretos para una EDO Resolveremos el problema de Cauchy (5.1), donde λser´a un par´ametro desconocido esta vez, concretamente se resolver´a el siguiente problema: dados uk, hallar λ∈R tal que u(tk) = uk, con kdesde 1 hasta el n´umero de observaciones que se utilice, siendo ula soluci´on de u′(t) = λt, t ∈[0,5], u(0) = 2.(6.2) El objetivo ahora es: dado un conjunto de observaciones uk, que se supone que sigue este modelo, determinar el par´ametro λque hace que la soluci´on los ajuste de la mejor forma posible. En la ((vida real)), cuando trabajamos con modelos biol´ogicos como el modelo SIR introducido anteriormente en (5.5), las observaciones las podemos obtener a partir de los datos reales que, a priori, sigue el modelo. En nuestro caso, no contamos con observaciones reales, pero, como tenemos la soluci´on exacta (5.2), lo que haremos ser´a utilizar como datos los valores de la soluci´on exacta para un par´ametro dado, y veremos si la red consigue aproximar dicho valor. Los datos los obtendremos con λ=−2, que ser´a el valor que debe aproximar la red neuronal. Adem´as, para los tiempos tise toma una partici´on uniforme del 52 Sergio Murillo Garc´ıa 6.2. AJUSTE DE PAR´ AMETROS DISCRETOS intervalo con un n´umero determinado de puntos (especificados en la Tabla 6.2 de hiperpar´ametros de la simulaci´on). Los hiperpar´ametros escogidos para la simulaci´on se recogen en la Tabla 6.2, y los resultados obtenidos en la Figura 6.2. Hiperpar´ametros Resultados Activaci´on Sigmoid Estructura [1,10,10,10,1] F. coste final 4·10−4 Optimizador Adamax ´ Epocas 15000 |ˆ λ−λexa|3·10−4 Learning rate lr= 0.06, ×0.9/(200 ´epocas) Nºpuntos 1000 Tiempo 65s Inicializaci´on λ λ0= 3 Nºobs. 50 Tabla 6.2: Hiperpar´ametros para la simulaci´on con los valores adecuados. Resoluci´on del problema (6.2). Podemos observar en la Figura 6.2 como el par´ametro desconocido λllega r´apidamente al valor real −2, con un error final de menos del 0.02 % que es bastante aceptable. Como podemos ver, la aproximaci´on de la red neuronal se ajusta satisfactoriamente a la soluci´on exacta. 53 CAP´ ITULO 6. PINN PARA EL AJUSTE DE PAR´ AMETROS (a) Observaciones usadas para el entrenamiento, a partir de la soluci´on real con λ=−2. (b) Learning rate frente a las ´epocas. (c) Funci´on de p´erdida frente a las ´epocas. (d) Evoluci´on del par´ametro desconocido λen funci´on de las ´epocas. (e) Soluci´on aproximada por la red neuronal y soluci´on exacta. Figura 6.2: Simulaci´on con los hiperpar´ametros de la Tabla 6.2. Resoluci´on del problema (6.2). 54 Sergio Murillo Garc´ıa 6.2. AJUSTE DE PAR´ AMETROS DISCRETOS 6.2.2. Ejemplo de ajuste de par´ametros discretos para un sistema de EDO’s Veamos en esta secci´on un problema inverso de tipo 2 aplicado a un sistema diferencial. Adem´as, en este caso tendremos 2 par´ametros desconocidos, en lugar de 1 y s´olo contaremos con informaci´on de algunas variables. Se realizar´an tambi´en algunas simulaciones para ver el funcionamiento de la red neuronal en el caso de que las observaciones introducidas tengan un porcentaje de ruido determinado, y en el caso de que solo se tenga informaci´on en una parte del dominio de resoluci´on. El caso de sistema que vamos a resolver es el modelo SIR que ya ha sido introducido en (5.5), pero ahora tanto β, γ, como S0ser´an desconocidos. Concretamente:        S′(t) = −βS(t)I(t), I′(t) = βS(t)I(t)−γI(t), t ∈[0,3], R′(t) = γI(t), S(0) = S0, I(0) = 50, R(0) = 0. (6.3) Adem´as, no tendremos observaciones de todas las variables S, I, R, s´olo dispondremos de observaciones de los infectados Iy de los recuperados R. Esta situaci´on puede darse en la vida real, ya que podr´ıamos no poder medir f´acilmente todas las variables de nuestro modelo. ¿De d´onde obtenemos las observaciones? El modelo SIR no tiene una soluci´on anal´ıtica y, como no tenemos acceso a datos reales, usaremos un solucionador externo a la red para obtener los datos, como puede ser la funci´on odeint de Python, que veremos en el Cap´ıtulo 8 de ´utiles inform´aticos. Estas observaciones se obtendr´an para unos valores de β,γyS0concretos, que son los valores reales que la red neuronal deber´a aproximar. Adem´as, como no tenemos datos de S(t), el valor de la condici´on inicial S(0) tambi´en es desconocido, aunque, como ya vimos en los problemas de ajuste tipo 1, no es necesario que forme parte de los par´ametros entrenables de la red neuronal, pues lo podremos obtener a posteriori, una vez la red neuronal est´e entrenada. Los hiperpar´ametros de la simulaci´on est´an recogidos en la Tabla 6.3 y los resultados de la simulaci´on y las observaciones utilizadas se encuentran en la Figura 6.3. Hiperpar´ametros Resultados Activaci´on Sigmoid Estructura [1,15,15,15,3] F. coste final 4.4·10−3 Optimizador Adamax ´ Epocas 15000 ||ˆ λ−λexa|| 4·10−5 Learning rate lr= 0.1, ×0.83/(230´epocas) Nºpuntos 1000 |ˆ S(0) −S0|3·10−2 Inicializaci´on λβ= 2 γ=−3Nºobs. 50 Tiempo 114s Tabla 6.3: Hiperpar´ametros para la simulaci´on con los valores adecuados (λ= (β, γ)). Resoluci´on del problema (6.3). 55 CAP´ ITULO 6. PINN PARA EL AJUSTE DE PAR´ AMETROS (a) Observaciones obtenidas con solucionador externo para β= 0.05, γ = 2, S0= 100. (b) Learning rate frente a las ´epocas. (c) Par´ametro desconocido βfrente a ´epocas. (d) Par´ametro desconocido γfrente a ´epocas. (e) Funci´on de p´erdida frente a ´epocas. (f) Soluci´on aproximada por la red neuronal y soluci´on obtenida con una funci´on externa a la red. Figura 6.3: Simulaci´on con los hiperpar´ametros de la Tabla 6.3. Resoluci´on del problema (6.3). Se observa, al igual que en el caso de una EDO, que ambos par´ametros son ajustados de forma adecuada y en pocas ´epocas, ya que podemos apreciar que en torno a la ´epoca 2000 se estabilizan en el valor exacto (β= 0.05, γ = 2). Adem´as, la red neuronal es capaz de aproximar de forma satisfactoria la poblaci´on de susceptibles 56 Sergio Murillo Garc´ıa 6.2. AJUSTE DE PAR´ AMETROS DISCRETOS S(t) aunque no hayamos introducido ninguna observaci´on sobre ella. La funci´on de p´erdida llega a valores muy peque˜nos, lo que nos indica que se minimiza correctamente. Otro de los aspectos en los que podemos pensar, es en que en la vida real, las observaciones que tengamos de medidas no seguir´an de forma tan precisa las ecuaciones diferenciales de nuestro problema, adem´as no tenemos por qu´e tener observaciones de las variables en todo el dominio. Por esa raz´on, se va a realizar otra simulaci´on introduciendo ruido en los datos obtenidos con el solucionador externo, y, adem´as, se van a tomar estas observaciones solo en la mitad del intervalo de resoluci´on, esto es, en [0,1.5]. Los hiperpar´ametros escogidos para la simulaci´on se recogen en la Tabla 6.4, y las gr´aficas obtenidas en la Figura 6.4. Hiperpar´ametros Resultados Activaci´on Sigmoid Estructura [1,15,15,15,3] F. coste final 5.22 Optimizador Adamax ´ Epocas 20000 ||ˆ λ−λexa|| 2·10−2 Learning rate lr= 0.1, ×0.7/(500 ´epocas) Nºpuntos 1000 |ˆ S(0) −S0|1.6 Inicializaci´on λβ= 2, γ=−3Nºobs. 50, en [0,1.5] Tiempo 160s Ruido 3 % Tabla 6.4: Hiperpar´ametros para la simulaci´on con los valores adecuados (λ= (β, γ)). Resoluci´on del problema (6.3). En este caso, no es de extra˜nar que el valor final de la funci´on de p´erdida sea mayor, ya que el error que se comete al aproximar los puntos con ruido es mayor que sin ´el. Adem´as, al tener ruido, puede haber varias soluciones muy pr´oximas que ((disten)) poco de la soluci´on concreta que buscamos, luego se entiende que el error sea mayor. Sin embargo, obtenemos los par´ametros buscados (β, γ yS0) con un error relativo menor al 3 % en los 3 casos. Se puede observar como la funci´on de p´erdida y los par´ametros alcanzan valores ((estables)) en pocas ´epocas, aun as´ı, se han dado m´as ´epocas para intentar que la red nos de valores m´as cercanos al exacto, ya que este problema es, necesariamente, m´as complejo. 57 CAP´ ITULO 6. PINN PARA EL AJUSTE DE PAR´ AMETROS Figura 6.9: Curva deseada para los infectados. de los susceptibles (que pasan a ser infectados), cuanto mayor sea βmayor ser´a el decrecimiento de los susceptibles, es decir, m´as infecciones se producen en (6.4). Luego, podemos asociar un βm´as grande, a un menor n´umero de restricciones y viceversa. Podemos a˜nadir, por tanto, a la funci´on de p´erdida un t´ermino donde reflejemos el coste que nos supone un beta mayor. La funci´on de p´erdida es entonces: C=1 Nc Nc X i=1    ˆ S′(ti) + ˆ β(t)ˆ S(ti)ˆ I(ti)   2+1 Nc Nc X i=1    ˆ I′(ti) + ˆ β(t)ˆ S(ti)ˆ I(ti)−γˆ I(ti)   2 +1 Nc Nc X i=1    ˆ R′(ti)−γˆ I(ti)   2+   ˆ S(tmin)−S0   2+   ˆ I(tmin)−I0   2+   ˆ R(tmin)−R0   2+ W 1 Nd Nd X j=1    ˆ I(tj)−Ideseado(tj)   2−   ˆ β(t)   2!. Podemos realizar un estudio similar al realizado en la secci´on anterior para encontrar un valor ´optimo del factor W. Realizando simulaciones con distintos valores con 10000 ´epocas de simulaci´on, obtenemos la Figura 6.10. El valor ´optimo es W= 10−5. Las simulaciones se realizan con los hiperpar´ametros de la Tabla 6.6, que se han considerado adecuados. Las gr´aficas obtenidas est´an recogidas en la Figura 6.11. En la Figura 6.11 podemos apreciar que el control calculado difiere bastante del obtenido en la secci´on anterior. Para ver si la aproximaci´on obtenida satisface el sistema diferencial, podemos resolver el sistema con el control aproximado que se obtiene, utilizando un algoritmo externo (l´ınea negra). Gracias a estas curvas, observamos que la aproximaci´on obtenida satisface de forma m´as o menos adecuada el sistema SIR. Por otra parte, los infectados evolucionan de forma que se ajustan lo m´aximo posible al estado final deseado (en azul claro), como tiene que cumplir m´as condiciones, puede no ser posible un mejor ajuste. En cuanto al control, vemos que 64 Sergio Murillo Garc´ıa 6.3. AJUSTE DE PAR´ AMETROS CONTINUOS Figura 6.10: Variaci´on del error en funci´on del factor W. Hiperpar´ametros Resultados Activaci´on Sigmoide Estructura [1,20,20,20,4] F. coste final 1.9·10−3 Optimizador Adamax ´ Epocas 15000 Error Pb. directo 1.9·10−3 Learning rate lr= 0.01, ×0.5/(800´epocas) Nºpuntos dom 1000 Error Pb. inverso 2.4·10−3 W10−5Nºpuntos cd. deseada 1000 Tiempo 222s Tabla 6.6: Hiperpar´ametros para la simulaci´on con los valores adecuados. Resoluci´on del problema (6.4). aumenta con el tiempo, lo que supondr´ıa mayor n´umero de restricciones al comienzo de la epidemia, y menos restricciones al final, lo cual tambi´en es coherente. 65 CAP´ ITULO 6. PINN PARA EL AJUSTE DE PAR´ AMETROS (a) Evoluci´on del ((learning rate)) frente a las ´epocas. (b) Evoluci´on del error frente a las ´epocas. (c) Control aproximado. (d) Soluci´on aproximada por la red neuronal, soluci´on con aproximador externo (para el control aproximado) y estado deseado. Figura 6.11: Gr´aficas para la simulaci´on con los hiperpar´ametros de la Tabla 6.6. Resoluci´on del problema (6.4). 66 Sergio Murillo Garc´ıa 6.3. AJUSTE DE PAR´ AMETROS CONTINUOS 67 Cap´ıtulo 7 Ecuaciones en derivadas parciales Como ya se ha dicho anteriormente, las ((PINN)) son f´acilmente adaptables a las ecuaciones en derivadas parciales. Basta adaptar la funci´on de p´erdida a la ecuaci´on correspondiente. De forma pr´actica, la red neuronal no distingue entre ecuaciones diferenciales ordinarias y ecuaciones en derivadas parciales, trata ambos problemas de la misma forma, a diferencia de los m´etodos cl´asicos de resoluci´on y aproximaci´on de problemas diferenciales. En este cap´ıtulo se va a resolver un problema directo relativo a una ecuaci´on en derivadas parciales y, posteriormente, se utilizar´a dicha aproximaci´on para obtener observaciones para un problema de ajuste de par´ametros. 7.1. Problema directo Consideramos una aproximaci´on usando ((PINN)) de la soluci´on de la ecuaci´on en derivadas parciales siguiente:          −∂2u ∂x2−∂2u ∂y2+ sen (x+y)u= 0, x, y ∈[0,1], u(x, 1) = sen (πx), u(x, 0) = u(0, y) = u(1, y)=0. (7.1) Para poder obtener una aproximaci´on externa con la que poder comparar la soluci´on obtenida con la red neuronal, se utiliza el m´etodo de las diferencias finitas (MDF). La funci´on de p´erdida de este problema ser´ıa: C=1 Nc Nc X i=1 ∥−ˆuxx(xi, yi)−ˆuyy(xi, yi) + sen (xi+yi)ˆu(xi, yi)∥2+ W·1 No No X k=1 ∥ˆu(xk,1) −sen (πxk)∥2+∥ˆu(xk,0)∥2+∥ˆu(0, yk)∥2+∥ˆu(1, yk)∥2, Sergio Murillo Garc´ıa 7.1. PROBLEMA DIRECTO donde W es el factor que ya se introdujo en el cap´ıtulo anterior. Podemos ver en la Tabla 7.1 los hiperpar´ametros que se utilizar´an para realizar la simulaci´on. Hiperpar´ametros Resultados Activaci´on Sigmoide Estructura [2,25,25,25,1] F. coste final 0.0051 Optimizador Adamax ´ Epocas 15000 Error Residuo 0.0046 Learning rate lr= 0.15, ×0.99/(20 ´epocas) Nºpuntos dom 1000 Error cd inicial 5·10−5 W10 Nºcd inicial 1000 Tiempo 488s Tabla 7.1: Hiperpar´ametros para la simulaci´on con los valores adecuados. Resoluci´on del problema (7.1). En las Figuras 7.1 y 7.2 se recogen los resultados de la simulaci´on realizada. (a) Evoluci´on del ((learning rate)) por ´epocas. (b) Evoluci´on del error por ´epocas. (c) Soluci´on obtenida con el MDF. (d) Soluci´on obtenida con la red neuronal. Figura 7.1: Gr´aficas para la simulaci´on con los hiperpar´ametros de la Tabla 7.1. Resoluci´on del problema (7.1). Podemos observar una ligera variaci´on en la soluci´on obtenida, principalmente en la condici´on de contorno, en la Figura 7.1 (d). 69 CAP´ ITULO 7. ECUACIONES EN DERIVADAS PARCIALES (a) Soluci´on obtenida con el MDF. (b) Soluci´on obtenida con la red neuronal. Figura 7.2: Gr´aficas para la simulaci´on con los hiperpar´ametros de la Tabla 7.1. Resoluci´on del problema (7.1). 7.2. Ajuste de par´ametros En esta secci´on vamos a ajustar un par´ametro continuo (funci´on ((control))) en una ecuaci´on en derivadas parciales. El par´ametro que ajustaremos ser´a la funci´on a(x, t) del siguiente problema:          −∂2u ∂x2−∂2u ∂y2+a(x, y)u= 0, x, y ∈[0,1], u(x, 1) = sen (πx), u(x, 0) = u(0, y) = u(1, y) = 0. (7.2) Se hace notar que el problema (7.1) que se resolvi´o en la secci´on anterior, es un caso particular para a(x, y) = sen (x+y). Al igual que hicimos en el cap´ıtulo anterior en la Secci´on 6.3, introduciremos una nueva neurona de salida a la red neuronal, de forma que una salida se corresponder´a con la propia aproximaci´on a la red u(x, y) y la segunda ser´a la aproximaci´on de la funci´on a(x, y). Como se trata de un problema inverso, daremos a la red neuronal observaciones (ul) de la soluci´on del problema para una funci´on a(x, y) concreta, y entrenaremos la red con el objetivo de que el mejor ajuste sea el que buscamos. Los valores de observaciones tomados ser´an los de la aproximaci´on obtenida con la red neuronal obtenida en la Secci´on 7.1. 70 Sergio Murillo Garc´ıa 7.2. AJUSTE DE PAR´ AMETROS La funci´on de coste pasa a ser ahora: C=1 Nc Nc X i=1 ∥−ˆuxx(xi, yi)−ˆuyy(xi, yi) + ˆa(xi, yi)ˆu(xi, yi)∥2+ W·1 No No X k=1 ∥ˆu(xk,1) −sen (πxk)∥2∥ˆu(xk,0)∥2+∥ˆu(0, yk)∥2+∥ˆu(1, yk)∥2 +1 Nd Nd X l=1 ∥ˆu(xl, yl)−ul∥2. La diferencia respecto al problema directo es el ´ultimo t´ermino que se ha a˜nadido para imponer que la red neuronal se ajuste a las Ndobservaciones que se van a usar en el entrenamiento. Los hiperpar´ametros con los que se realiza la simulaci´on se han recogido en la Tabla 7.2. Hiperpar´ametros Resultados Activaci´on Sigmoide Estructura [2,30,30,30,2] F. coste final 0.11 Optimizador Adamax ´ Epocas 5000 Error Obs 0.10 Learning rate lr= 0.05, ×0.99/(20 ´epocas) Nºpuntos dom 1000 Tiempo 133s W10 Nºpuntos cd inicial 1000 Nºpuntos observaciones 50 ×50 Tabla 7.2: Hiperpar´ametros para la simulaci´on con los valores adecuados. Resoluci´on del problema (7.2). Los puntos tomados en las observaciones, son el resultado de hacer un mallado para 50 puntos en cada uno de los ejes, es decir, un total de 2500 observaciones. Los resultados de la simulaci´on se pueden observar en la Figura 7.3. Se puede apreciar que la aproximaci´on obtenida para la funci´on a(x, y) no es igual que la funci´on sen (x+y) que esper´abamos. Esto puede deberse a varios factores. En primer lugar, no podemos descartar que la red neuronal se haya ((atascado)) en un m´ınimo local de la funci´on de coste, aunque teniendo en cuenta el error tan peque˜no que se obtiene, el m´ınimo local deber´ıa estar pr´oximo a uno global. Otro factor, quiz´as el m´as importante, es que el problema no tiene por qu´e tener soluci´on ´unica, ya que vemos que aunque la funci´on no es la misma, la soluci´on obtenida se ajuste mucho a la esperada. 71 CAP´ ITULO 7. ECUACIONES EN DERIVADAS PARCIALES (a) Evoluci´on del ((learning rate)) por ´epocas. (b) Evoluci´on del error por ´epocas. (c) Soluci´on obtenida con la RN. (d) Soluci´on obtenida con la RN. (e) Ajuste exacto sen (x+y). (f) Ajuste obtenido con la red neuronal. Figura 7.3: Gr´aficas para la simulaci´on con los hiperpar´ametros de la Tabla 7.2. Resoluci´on del problema (7.2). 72 Sergio Murillo Garc´ıa 7.2. AJUSTE DE PAR´ AMETROS 73 CAP´ ITULO 8. HERRAMIENTAS INFORM´ ATICAS Y C´ ODIGOS 23 24 U=integrate.odeint(aux, y0, tspace) 25 # Graficamos la aproximaci´on y la soluci´on exacta 26 fig, ax1 =plt.subplots() 27 ax1.plot(x,yh_plot,color='blue',label='PINN',linewidth=6.0) 28 ax1.plot (tspace,U, color='orange', label='exacta',linewidth=3.0) 29 ax1.set_xlabel('t',color='black') 30 ax1.set_ylabel('u(t)',color='black') 31 ax1.legend(loc ='upper left') Y con estos c´odigos obtendr´ıamos gr´aficas para la evoluci´on del learning rate y de la funci´on de p´erdida en funci´on de la ´epoca. Adem´as, tambi´en tenemos la representaci´on de la aproximaci´on obtenida por la red neuronal, el aproximador externo en este caso es simplemente integrar la funci´on, dada la simplicidad de la ecuaci´on diferencial. 8.2. Ajuste de par´ametros discretos Para estos problemas necesitaremos a˜nadir de forma manual un par´ametro nuevo a los ya existentes en la red neuronal, que se ajustar´a junto con el resto. La mayor parte del c´odigo es an´aloga al de la secci´on anterior, se puntualizar´an, por tanto, las diferencias. 1# Definimos los par´ametros necesarios para la resoluci´on 2 3steps=5000 # n´umero de iteraciones 4ruido=0 # % de ruido para las observaciones 5 6A=-2.0 # valor exacto para el parametro lambda 7y0=2 #valor de la condicion inicial 8min=0 # Dominio de resolucion 9max=5 10 11 lambda1= 3.0 # inicializacion del par´ametro 12 13 lr=0.06 # learning rate inicial 14 layers =np.array([1,10,10,10,1]) # estructura de la red 15 16 Ndom= 1000 # puntos para evaluar las p´erdidas 17 Nobs=50 # puntos de observaciones Se han definido los valores concretos del problema. Se define adem´as el valor inicial para el par´ametro que se quiere ajustar lambda1. Debemos ahora definir el funcionamiento de la red neuronal, al igual que en la secci´on anterior. En este momento es cuando a˜nadimos de forma manual un par´ametro extra a la red. 80 Sergio Murillo Garc´ıa 8.2. AJUSTE DE PAR´ AMETROS DISCRETOS 1class FCN(nn.Module): 2##Neural Network 3def __init__(self,layers): 4super().__init__() #call __init__ from parent class 5 6'Funci´on de activaci´on' 7self.activation =nn.Sigmoid() 8 9'loss function'# tipo de funci´on de p´erdida 10 self.loss_function =nn.MSELoss(reduction ='mean') 11 12 'Initialise neural network as a list using nn.Modulelist' 13 self.linears =nn.ModuleList([nn.Linear(layers[i], layers[i+1]) for iin range(len(layers)-1)]),→ 14 self.iter = 0 15 16 'Inicializamos el par´ametro (Inverse problem)' 17 self.lambda1 =torch.tensor([lambda1], requires_grad=True).float().to(device),→ 18 19 'Registramos para el entrenamiento' 20 self.lambda1 =nn.Parameter(self.lambda1) 21 22 'Xavier Normal Initialization' 23 for iin range(len(layers)-1): 24 nn.init.xavier_normal_(self.linears[i].weight.data, gain=1.0),→ 25 nn.init.zeros_(self.linears[i].bias.data) 26 27 'foward pass' # programamos el paso hacia delante 28 def forward(self,x): 29 if torch.is_tensor(x) != True: 30 x=torch.from_numpy(x) 31 a=x.float() 32 for iin range(len(layers)-2): 33 z=self.linears[i](a) 34 a=self.activation(z) 35 a=self.linears[-1](a) 36 return a 37 38 'Loss Functions' # se programa la funci´on de p´erdida 39 def loss_BC(self): # p´erdida en la cd inicial 40 loss_i=self.loss_function(self.forward(min*torch.ones(1,1)),\ 41 y0*torch.ones(1,1)) 42 return loss_i 43 44 def loss_PDE(self,x_PDE): # p´erdida en la ecuacion diferencial 45 lambda1=self.lambda1 # nuestro parametro a optimizar 81 CAP´ ITULO 8. HERRAMIENTAS INFORM´ ATICAS Y C´ ODIGOS 46 g=x_PDE.clone() # la g es la var. independiente (x o t) 47 g.requires_grad=True #Enable differentiation 48 f=self.forward(g) # f es la salida de la red 49 f_x=autograd.grad(f,g,torch.ones([x_PDE.shape[0],1]).\ 50 to(device),retain_graph=True, create_graph=True)[0] 51 loss_f=self.loss_function(f_x,(lambda1)*g) 52 return loss_f 53 54 def loss_data(self,x_obs,y_obs): # p´erdida en las observaciones 55 loss_obs =self.loss_function(self.forward(x_obs), y_obs) 56 return loss_obs 57 58 def loss(self,x_PDE,x_obs,y_obs): # p´erdida total 59 loss_i=self.loss_BC() 60 loss_f=self.loss_PDE(x_PDE) 61 loss_obs=self.loss_data(x_obs,y_obs) 62 return loss_i +loss_f +loss_obs La diferencia con el problema directo es el par´ametro extra que hemos a˜nadido, el resto de partes, como la definici´on del paso hacia delante y la funci´on de p´erdida, se hacen de la misma forma. Como para el problema de ajuste de par´ametros se necesitan observaciones del problema, tambi´en debemos introducirlas: x obs ey obs. Solo queda definir las observaciones y tiempos de entrenamiento que introduciremos en la red. 1# Generamos los datos 2x_obs=torch.linspace(min,max,Nobs).view(-1,1) 3y_obs=(A/2)*(x_obs**2)+y0 #observaciones 4 5# Introducimos ruido 6for iin range(Nobs): 7y_obs[i]=y_obs[i]+random.randint(-1,1)*ruido*y_obs[i] 8y_obs=y_obs.view(-1,1) 9 10 # Cargamos los datos 11 x_PDE=torch.linspace(min,max,Ndom).view(-1,1) 12 x_PDE=x_PDE.float().to(device) 13 x_obs=x_obs.float().to(device) 14 y_obs=y_obs.float().to(device) 15 16 # Cargamos el modelo y elegimos optimizador 17 model =FCN(layers) 18 model.to(device) 19 params =list(model.parameters()) 20 optimizer =torch.optim.Adam(model.parameters(),lr=lr,amsgrad=False) 21 scheduler =optim.lr_scheduler.StepLR(optimizer, step_size=200, gamma=0.9),→ 82 Sergio Murillo Garc´ıa 8.3. AJUSTE DE PAR´ AMETROS CONTINUOS Las observaciones introducidas se han generado con el valor real del par´ametro que buscamos, ya que en este ejemplo ((educativo)) sabemos con anterioridad lo que debemos obtener, en una situaci´on real estos ser´an datos experimentales. Adem´as, se programa la introducci´on de ruido en los datos. Al igual que en el apartado anterior, programamos el entrenamiento y almacenamos las variables para las que queramos ver la evoluci´on. 1# Programamos el entrenamiento 2from time import time 3start_time =time() 4t0=time() 5hist=torch.Tensor([]) 6par=torch.Tensor([]) 7lrate=torch.Tensor([]) 8 9for iin range(steps): 10 loss =model.loss(x_PDE,x_obs,y_obs) 11 nuevovalor=model.lambda1 12 optimizer.zero_grad() 13 14 loss.backward() 15 optimizer.step() 16 scheduler.step() 17 18 par=torch.cat((par,torch.Tensor([nuevovalor])),0) 19 hist=torch.cat((hist,torch.Tensor([loss])),0) 20 lrate=torch.cat((lrate,torch.\ 21 Tensor([optimizer.param_groups[0]['lr']])),0) 22 if i%(steps/10)==0: 23 print('It {:05d}: loss = {:10.8e}'.format(i,loss)) 24 print('\nComputation time: {} seconds'.format(time()-t0)) En este caso se va a almacenar tambi´en la evoluci´on del par´ametro a ajustar. A partir de este punto, se deben representar la aproximaci´on obtenida con la red neuronal y la evoluci´on de los par´ametros, de igual forma que en la secci´on anterior. 8.3. Ajuste de par´ametros continuos Estos problemas se pueden programar partiendo del problema directo, a˜nadiendo simplemente una nueva neurona en la capa de salida. Esta nueva salida ser´ıa la funci´on de ((control)) que buscamos obtener. Las condiciones extra que van a determinar esta funci´on se a˜naden en la funci´on de p´erdida. 83 Bibliograf´ıa [1] OpenAI,ChatGPT, 2022, https://openai.com/blog/chatgpt. [2] Claudio Guti´ errez, Andr´ es Abeliuk,Historia y evoluci´on de la inteligencia artificial,Universidad de Chile, Revista BITS, 2021, https://revist asdex.uchile.cl/index.php/bits/article/download/2767/2700/10150. [3] Varol Akman, Patrick Blackburn,Editorial: Alan Turing and Artificial Intelligence,Journal of Logic, Language, and Information, 2000, https://ww w.researchgate.net/figure/Turings-1950-paper-is-one-of-the-mos t-cited-in-philosophical-AI-literature-Reproduced_fig2_28762974. [4] Joseph Weizenbaum,ELIZA-A Computer program for the study of natural lenguage communication between man and machine,Communications of the ACM, 1966, https://dl.acm.org/doi/10.1145/365153.365168 [5] Ian Goodfellow and Yoshua Bengio and Aaron Courville,Deep Learning,MIT Press, 2016, http://www.deeplearningbook.org. [6] Anthony L. Caterini and Dong Eui Chang,Deep Neural Networks in a Mathematical Framework,Springer International Publishing, 2018, http://li nk.springer.com/10.1007/978-3-319-75304-1. [7] G Cybenkot,Approximation by Superpositions of a Sigmoidal Function, Math. Control Signals Systems, 1989. [8] Haim Br´ ezis,An´alisis funcional. Teor´ıa y aplicaciones, 84-206-8088-5, Alianza Editorial, 1983. [9] Robert B Ash,Real Analysis and Probability. [10] Walter Rudin,Real and complex analysis, 0-07-100276-6, McGRAW-HILL internacional editions, 1987. [11] Jos´ e Luis Sarmiento-Ramos,Aplicaciones de las redes neuronales y el deep learning a la ingenier´ıa biom´edica, 10.18273/revuin.v19n4-2020001, 16574583, Revista UIS Ingenier´ıas, 2020. [12] Jordi Torres,PYTHON DEEP LEARNING. Introducci´on pr´actica con Keras y TensorLow 2, 978-607-538-614-0, Alfaomega, 2020. Sergio Murillo Garc´ıa BIBLIOGRAF´ IA [13] James Dellinger,Weight Initialization in Neural Networks: A Journey From the Basics to Kaiming,Towards Data Science, 2019, https://towardsdatas cience.com/weight-initialization-in-neural-networks-a-journey-f rom-the-basics-to-kaiming-954fb9b47c79. [14] Luis Velasco,Optimizadores en redes neuronales profundas: un enfoque pr´actico, 2020, https://velascoluis.medium.com/optimizadores-en-r edes-neuronales-profundas-un-enfoque-pr%C3%A1ctico-819b39a3eb5. [15] Jes´ us Martinez,¿Qu´e es un optimizador y para qu´e se usa en Deep Learning?,DATASMARTS, 2020, https://datasmarts.net/es/que-es-un-opt imizador-y-para-que-se-usa-en-deep-learning/. [16] B.B. Hammel,What learning rate should I use?, 2019, http://www.bdhamm el.com/learning-rates/. [17] Vaibhav Haswani,Learning Rate Decay and methods in Deep Learning, 2020, https://medium.com/analytics-vidhya/learning-rate-decay-and-met hods-in-deep-learning-2cee564f910b#:~:text=Learning%20rate%20dec ay%20is%20a,help%20both%20optimization%20and%20generalization.. [18] Shree Nayar,Lecture series ’First Principles of Computer Vision’,Computer Science Department, School of Engineering and Applied Sciences, Columbia University, 2021, https://www.youtube.com/watch?v=sIX_9n-1UbM. [19] Michael A. Nielsen,Neural Networks and Deep Learning,Determination Press, 2015, http://neuralnetworksanddeeplearning.com/. [20] Davis E. Rumelhart, Geoffrey E. Hinton, Ronald J. Williams, Learning representations by back-propagating errors Letters to nature, 1986, https://www.iro.umontreal.ca/~vincentp/ift3395/lectures/backprop_ old.pdf. [21] Yu.L. Levitan, N.I. Markovich, S.G. Rozin, I.M. Sobol,On quasirandom sequences for numerical computations,Pergamon Press, 1988. [22] Wei-Liem Loh,ON LATIN HYPERCUBE SAMPLING The Annals of Statistics, 1996. [23] M. Raissi, P. Perdikaris, G.E. Karniadakis,Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, ISSN:0021-9991, ELSEVIER: Journal of Computational Physics, 2091, https://www.sciencedirect.com/ science/article/pii/S0021999118307125. [24] Howard Weiss,The SIR model and the Foundations of Public Health, ISSN: 1887-1097, Materials Matem`atics, 2013, https://mat.uab.cat/web/matmat/e s/v2013n03/. 85 BIBLIOGRAF´ IA [25] Saviz Mowlavi,Saleh Nabi,Optimal control of PDEs using physicsinformed neural networks, ISSN:0021-9991, Journal of Computational Physics, 2023, https://www.sciencedirect.com/science/article/pii/S002199912 200794X. [26] Python Software Foundation,PYTHON,https://www.python.org/. [27] Eleonora Fanouraki,State of the Developer Nation 23rd edition: the fall of web frameworks, coding languages, blockchain, and more!,Towards Data Science, 2022, https://www.slashdata.co/blog/state-of-the-developer-nat ion-23rd-edition-the-fall-of-web-frameworks-coding-languages-blo ckchain-and-more?. [28] John Terra,Keras vs Tensorflow vs Pytorch: Key Differences Among Deep Learning, 2023, https://www.simplilearn.com/keras-vs-tensorflow-v s-pytorch-article#what_is_pytorch. [29] Juan Diego Toscano,Repositorio GitHub, 2023, https://github.com/jdt oscano94/Learning-Python-Physics-Informed-Machine-Learning-PINNs -DeepONets. [30] Raissi, M., Perdikaris, P., Karniadakis, G. E.,Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, ISSN=0021-9991, Journal of Computational Physics, 2019, https://www.sciencedirect.com/science/ar ticle/pii/S0021999118307125. [31] PyTorch Contributors,TORCH.OPTIM,https://pytorch.org/docs/s table/optim.html#torch.optim.Optimizer. 86