Optimización en el espacio de parámetros de un modelo
(El modelo de NetLogo asociado a esta entrada lo puedes encontrar aquí.)
Uno de los problemas más habituales cuando se construye un modelo para simular un proceso es que, tras haber definido un buen número de parámetros para darle más flexibilidad y generalidad con el fin de abarcar las situaciones más variopintas del proceso, no sabemos qué valores de esos parámetros serán los que puedan generar un comportamiento que resulte interesante. Para resolverlo, la primera tarea que abordamos es recorrer el espacio de parámetros como si fuera un bizcocho, pinchando aquí y allá para ver cómo se comporta el modelo en esos puntos aislados y esperando reconocer qué relación hay entre las zonas exploradas y los comportamientos observados. Éste problema de búsqueda en el espacio paramétrico, aunque bastante común, no está resuelto en general, ya que la dimensión del espacio suele ser alta (tanta como parámetros hayamos introducido), muchos de los parámetros pueden tomar valores continuos (habitualmente, muchos parámetros suelen tomar valores en un intervalo de \({\bf R}\)), y con una estructura de interrelaciones entre los diversos parámetros que muchas veces es compleja y alejada de un comportamiento lineal.
En esta ocasión vamos a ver un ejemplo muy particular de cómo ajustar estos parámetros para conseguir un fin concreto, que es el que se deriva de la optimización de alguna medición (en forma de variable medible) que podemos realizar sobre nuestro modelo. Aunque se presentará en forma de ejemplo concreto, esperamos que se pueda extraer una metodología relativamente general que puede ser aplicada a muchos otros casos.
Nuestro objetivo no es presentar una herramienta o metodología de carácter general, sino solo mostrar una idea de cómo combinar adecuadamente diversos modelos para facilitar la exploración de los comportamientos de uno de ellos. Una "composición" de modelos que, en este caso particular, se hará adaptando implementaciones de ellos realizadas previamente en NetLogo. Sin duda, no es el lenguaje más apropiado para realizar este tipo de composiciones de forma general, ya que no permite nada parecido al encapsulado de procedimientos ni a la utilización de espacios de trabajo independientes, pero a cambio NetLogo requiere muy poco esfuerzo para crear, adaptar y visualizar los resultados de nuestros modelos.
Usando la librería PSO que puedes encontrar en los recursos del curso de IA es bastante fácil aplicar el algoritmo PSO para optimizar los parámetros de cualquier modelo. En cualquier caso, a continuación explicamos cómo puede realizarse la misma optimización adaptando directamente el modelo que queremos optimizar y las funciones de PSO que necesitaremos.
Aplicación manual de PSO sobre modelos existentes
En particular, y solo por el hecho de seleccionar un sistema de optimización que es fácilmente implementable en NetLogo, vamos a mostrar cómo podemos aplicar una optimización basada en enjambres de partículas (PSO) sobre un segundo modelo (M). Para fijar ideas vamos a hacer las siguientes suposiciones:
- M dispone de varios parámetros que se pueden ajustar, y que determinan la evolución del modelo.
- En el mismo modelo podemos tomar varias medidas como resultado de su ejecución.
- M dispone de un procedimiento (Setup) para inciar el modelo, y de un procedimiento (Go) para ejecutar cada paso de evolución del mismo.
- En M fijamos algunos de sus parámetros y dejamos libres algunos otros: \(x_1,...,x_n\);
- En M fijamos uno de sus resultados medibles: \(z_1\).
- Como resultado, podemos interpretar el modelo como una función que recibe como dato de entrada \((x_1,...,x_n)\) y devuelve \(z_1\): notaremos la función por la misma letra \(M(x_1,...,x_n)=z_1\).
La idea de optimizar el modelo M (respecto a \(z_1\)) es ver para qué valores de los parámetros obtenemos el máximo de esa medida... o lo que es lo mismo, optimizar la función anterior asociada a M.
Cualquier curso básico de optimización (o una búsqueda rápida acerca de recursos de optimización en internet) muestra muchos y diversos métodos de optimización de funciones. Como hemos comentado, utilizaremos para este ejemplo un algoritmo PSO para optimizar la función que calcula M debido a que tiene un comportamiento muy sencillo y lo tenemos ya implementado en NetLogo. Como en general la función que calcula M no la podemos calcular directamente por medio de una expresión matemática, tendremos que ejecutar el modelo asociado y medir el valor de \(z_1\) de forma directa sobre el modelo tras haber permitido su evolución. Para poder hacer esto en NetLogo hemos de unir los modelos PSO y M para poder ejecutarlos en un espacio de programación común; es más, debemos poder ejecutar M desde el modelo PSO, y para conseguirlo se debe tener cuidado con unos pocos detalles:
Deben separarse ambos modelos desde el punto de vista de programación. Es decir, no deben tener objetos en común: ni razas de tortugas, ni variables globales, ni variables de patches (que es el único conjunto de agentes que no podemos personalizar completamente), ni nombres de procedimientos.- Como la evaluación de la función se debe llevar a cabo desde cada partícula, se deben añadir los parámetros que hemos dejado libres del modelo M como propiedades de las partículas. De hecho, en el caso concreto de PSO, estas variables serán las que se usarán como coordenadas de las partículas, ya que representan el espacio de movimiento de las mismas. Si solo son 2, podemos reutilizar adecuadamente xcor/ycor para que el movimiento de las partículas en el espacio de parámetros se corresponda con el movimiento de la raza de tortugas asociada al PSO en el mundo de NetLogo.
- Se debe definir una función de evaluación(que será ejecutada por las partículas) que haga lo siguiente:
- Ejecuta el Setup de M de acuerdo a los parámetros almacenados en la partícula, de forma que prepara un estado inicial de M acorde a la información que tiene la partícula.
- Ejecuta el modelo M a partir de ese estado inicial. Para ello, normalmente se ejecutará el Go de M un número de pasos preestablecido, o hasta que se cumpla una cierta condición.
- Una vez que M se ha detenido (o lo hemos detenido) se mide y devuelve el valor de \(z_1\). En algunos casos podríamos estar interesados en devolver como medida algún agregado considerado a lo largo de la ejecución del modelo, por ejemplo, la media de valores de una cierta variable a lo largo de la ejecución.
- Si el modelo incluye algún tipo de aleatoriedad, quizás es deseable hacer que este paso se repita algunas veces y la función devuelva el valor medio de los \(z_1\) medidos (y quizás también desviación estándar o cualquier medida adicional).
Al igual que ocurre con los Algoritmos Genéticos, a la hora de aplicar PSO a un problema de optimización hay dos pasos clave: la representación del espacio de soluciones, y la función de bondad (fitness). Una de las ventajas que presenta PSO es que puede trabajar de forma natural con números reales, y las partículas se mueven por tanto en el espacio paramétrico sin necesidad de codificaciones intermedias. No ocurre así en el caso de los Algoritmos Genéticos, donde es necesario usar codificaciones de los problemas en cadenas (normalmente binarias) que puedan ser manipuladas por operadores genéticos. Por ejemplo, si intentamos encontrar una optimización de la función \(f(x) = x_1^2 + x_2^2+x_3^2\), las partículas se definen inmediatamente como ternas \((x_1, x_2, x_3)\), y la función de fitness es directamente la propia \(f(x)\). En consecuencia, el proceso de búsqueda se convierte en un proceso iterado muy cercano al problema real que se pretende resolver, donde la condición de parada vendrá dada por haber alcanzado el número máximo de iteraciones, o hemos satisfecho una condición de error mínimo.
Esta proyección natural de muchos problemas a las partículas que se usan en PSO hace que no haya muchos parámetros que deban ser ajustados para conseguir buenos resultados y, como veremos a continuación, la mayoría de ellos o bien quedan determinados por el propio problema que se está optimizando, o bien pueden prefijarse de forma independiente al problema en cuestión y se conocen (experimentalmente) valores que dan buenos resultados. Los más típicos son:
- El número de partículas: el rango típico es entre 20 y 40. Realmente, para la mayoría de los problemas basta con 10 partículas para tener buenos resultados, y para los casos más difíciles se puede intentar con rangos entre 100 y 200 (muy por debajo de la "población" necesaria para otros métodos de optimización basados en poblaciones, como los algortimos genéticos).
- Dimensión de las partículas: viene determinado por la dimensión del espacio paramétrico del problema a optimizar (número de parámetros).
- Rango de las partículas (valores que toman): también viene determinado directamente por el problema, y tenemos la ventaja de poder especificar diferentes rangos para cada una de las dimensiones (o parámetros) de las partículas.
- Vmax: determina la máxima variación que puede sufrir la posición de una partícula (sus parámetros) en cada iteración.
- Factores de aprendizaje/atracción: Normalmente, y por razones experimentales, suelen tomarse como \(c_1=c_2=2\) (recordemos que \(c_1\) es la atracción al mejor personal de la historia de la partícula, y \(c_2\) la atracción al mejor global encontrado por todas las partículas hasta el momento). Sin embargo, se tiene libertad para acomodarlos a las necesidades del problema. En la mayoría de los trabajos se suelen tomar \(c_1=c_2\) y ambos en \([0, 4]\).
- Condición de parada: habitualmente suele venir dada por la combinación de un máximo número de iteraciones y un error mínimo que se debe alcanzar, que depende del problema a resolver.
- Comunicación Global / Comunicación local: la comunicación global (que el máximo global se considere entre todas las partículas) es mucho más rápida, pero tiene el problema de que cae fácilmente en óptimos locales en algunos problemas. La comunicación local (que el máximo global para cada partícula se calcula únicamente entre unas cuantas partículas, el "entorno" de la partícula) es más lenta, pero no queda atrapada con tanta facilidad en esos óptmos locales. Puede usarse una combinación de ambas, usando la global para encontrar resultados rápidamente, y la local para refinar la búsqueda.
- Factor de inercia: se comporta como un parámetro similar a los factores de aprendizaje, varía en \([0,1]\) y si es próximo a 1 da mayor independencia a las partículas, permitiendo explorar más secciones del espacio paramétrico, pero a cambio de obtener las soluciones de forma más lenta.
Como demostración de lo anterior sobre un modelo sencillo, vamos a usar Traffic Basic (que se puede encontrar en la biblioteca de modelos de NetLogo bajo el nombre /Social Science/Traffic Basic).
Este modelo muestra el movimiento de coches en una carretera. Cada coche sigue un conjunto de reglas muy simple: frena si ve un coche cerca delante suya y acelera en caso contrario (tienen una visión de 1 patch). El modelo demuestra que pueden producirse atascos sin necesidad de accidentes en la carretera, cortes o vehículos excesivamente lentos... decir, no es necesaria una causa centralizada para la existencia de atascos, sino que surje por los efectos de la aceleración y deceleración de los propios coches. El modelo depende de 3 parámetros:
- Número de coches: como el espacio es limitado, a mayor número de coches mayor probabilidad de tener que frenar cada coche.
- Aceleración: indica la aceleración que hace cada coche si no tiene otro delante.
- Deceleración: indica la frenada que hace cada coche si encuentra otro delante.
Vamos a trabajar solo con 2 parámetros en el PSO: aceleración/deceleración con unos rangos permitidos de \([0,0.01]\times [0,0.1]\), y como resultado mediremos la velocidad media de todos coches al cabo de 500 pasos. Por tanto, la función que vamos a optimizar es:
M(a,d) = mean [speed] of cars
Fijamos previamente el número de coches, y el número de iteraciones que permitimos hacer al modelo que estamos optimizando. El hecho de quedarnos solo con 2 parámetros es simplemente por facilitar la visualización del proceso de optimización, ya que así podemos asociar las coordenadas de las partículas (xcor, ycor) con los parámetros que definen cada ejecución del modelo de tráfico. Para ello, habremos de tener en cuenta que hemos de transformar las coordenadas que la partícula tiene en el mundo en los valores adecuados del espacio paramétrico. Concretamente, nuestro mundo tiene unas dimensiones de \([-25,25]\times [-25,25]\), por tanto, para las coordenadas \((x,y)\) de las partículas en el mundo hemos de aplicar la transformación:
\[(x,y) \mapsto (\frac{25 + x}{5000}, \frac{25 + y}{500})\]
Tras haber unido en un solo fichero ambos modelos (el de tráfico y el de PSO) teniendo en cuenta los detalles a los que haciamos mención antes, creamos una función (que pueden ejecutar las partículas) que les permite simular el modelo de tráfico y devolver la media de velocidad de los coches:
to-report funcion-caja-negra [x y]
; Para que vaya más rápido, se impide la actualización gráfica
no-display
; Generación de los coches para este experimento
setup-cars
; Se fijan los valores de aceleración y deceleración asociados a esta partícula
set acceleration x
set deceleration y
; Se deja correr el modelo de coches 500 pasos
repeat 500 [go-modelo]
display
; Se devuelve la función que se quiere optimizar: la velocidad media de los coches
report mean [speed] of cars
end
donde setup-cars es el nombre del procedimiento que inicia el modelo de tráfico, y go-modelo es el procedimiento que permite dar un paso en ese modelo. Debido a que el modelo de tráfico no se muestra (se ejecuta en segundo plano, pero sin visualizarse) su código ha sido limpiado de todo lo que afecta a la visualización y reducido a aquellos procedimientos que son imprescindibles para el cálculo del resultado.
Posteriormente, en los procedimientos que ejecutan el PSO, cada vez que una partícula quiere calcular el valor asociado a sus coordenadas/parámetros debe ejecutar:
let val funcion-caja-negra ((xcor + 25) / 5000) ((ycor + 25) / 500)
En el código que se adjunta a esta entrada también se puede encontrar un reporte que permite calcular la media de \(n\) repeticiones de la funcion-caja-negra. En la llamada anterior se puede sustituir funcion-caja-negra por f-repetida:
to-report f-repetida [n x y]
report mean (n-values n [funcion-caja-negra x y])
end
Observa que, por lo demás, ambos modelos, el de tráfico y el de PSO, son prácticamente iguales a los de los modelos originales.
El modelo de NetLogo asociado a esta entrada ejemplo lo puedes encontrar aquí.