Algoritmo Genético Bio-Inspirado en Cáncer (CancerGA)

Planteamiento matemático (implementación v0.1.0)

1 Visión general

CancerGA es un algoritmo genético co-evolutivo bio-inspirado en la dinámica de un tumor dentro de su microambiente. Dos poblaciones coexisten en el mismo espacio de búsqueda:

• Células cancerosas \(\mathcal{C}\): portan las soluciones candidatas al problema.
• Linfocitos \(\mathcal{L}\): agentes de control que ejercen presión selectiva sobre el cáncer.

Un microambiente de fibroblastos (campo derivado de los nichos, activo por defecto) refuerza localmente a las regiones que van mejor que el promedio.

El algoritmo resuelve problemas mono-objetivo y multiobjetivo, con variables continuas, enteras o mixtas, y restricciones lineales o no lineales, de igualdad o desigualdad, declaradas como funciones escalares (sección 2.3). Cada módulo biológico tiene una contraparte algorítmica y un lugar en el repositorio:

Módulo biológico	Operador algorítmico	Sección	Implementación (cancer_ga/)
Evaluación y ordenamiento	Dominancia restringida (Deb) + NSGA-II	4	utils/nsga2.py, algorithm.py::CancerGA.evaluate
Mutación adaptativa	Estrategia evolutiva auto-adaptativa	5.1	operators/mutation.py::adaptive_mutation
Quimiotaxis	Movimiento guiado por gradiente (diferencias finitas)	5.2	operators/chemotaxis.py::chemotaxis
Duelo biológico	Interacción probabilística cáncer–linfocito	5.3	operators/duel.py::biological_duel
Cruce	BLX-\(\alpha\) intra-población	5.4	operators/crossover.py::crossover
Proliferación por nichos	Crecimiento logístico local	5.5	operators/proliferation.py::niche_proliferation, assign_niches
Metástasis	Saltos con vuelos de Lévy	5.6	operators/metastasis.py::metastasis
Muerte por inanición	Balance de energía	5.7	operators/starvation.py::starvation
Fibroblastos	Campo de soporte por nicho	5.8	operators/proliferation.py::compute_fibroblast_field
Ciclo evolutivo	Bucle principal	7	algorithm.py::CancerGA.run

Alcance y versión. Este documento describe la implementación de cancer_ga versión 0.1.0 (https://github.com/nanodom/Cancer-GA). Cuando documento y código discrepen, prevalece el código. Las limitaciones conocidas de la implementación se listan en la sección 12.

2 Planteamiento del problema

2.1 Problema de optimización

Se busca resolver

\[ \begin{aligned} \min_{x}\quad & F(x)=\big(f_1(x),\dots,f_M(x)\big) \\ \text{s.a.}\quad & g_j(x)\le 0, \quad j=1,\dots,J \\ & h_k(x)=0, \quad k=1,\dots,K \\ & l_i\le x_i\le u_i, \quad i=1,\dots,n \\ & x_i\in\mathbb{Z}, \quad i\in\mathcal{I}\subseteq\{1,\dots,n\} \end{aligned} \]

con \(x\in\mathbb{R}^n\), \(M\ge 1\) objetivos y \(\mathcal{I}\) el conjunto de índices de variables enteras. Para \(M=1\) el problema es mono-objetivo.

2.2 Convención de optimización

El algoritmo siempre minimiza y no hay un indicador de signo por objetivo: para maximizar \(f_m\) basta entregar \(-f_m\) a la función objetivo. En lo que sigue, \(F(x)\) es el vector de objetivos tal como lo devuelve objective_function (un escalar o un arreglo de \(M\) valores).

Si la función objetivo devuelve una tupla (objs, g) (como los benchmarks CTP, CF y MW del repositorio), solo se usa objs: las restricciones g se ignoran y deben declararse por separado mediante constraints y eq_constraints (ver sección 12 y Apéndice B).

2.3 Violación de restricciones

Las restricciones se declaran como funciones escalares: constraints para \(g_j(x)\le 0\) y eq_constraints para \(h_k(x)=0\). Para una tolerancia de igualdad \(\varepsilon_h\) (eq_tol), la violación total es

\[ cv(x)=\sum_{j=1}^{J}\max\{0,\,g_j(x)\}+\sum_{k=1}^{K}|h_k(x)|\;\mathbb{1}\{|h_k(x)|>\varepsilon_h\}\;\ge 0 . \]

Una igualdad con \(|h_k|\le\varepsilon_h\) no aporta violación, pero una que excede la tolerancia aporta \(|h_k|\) completo (no \(|h_k|-\varepsilon_h\)), de modo que \(cv\) es discontinua en \(\varepsilon_h\). Un punto es factible si \(cv(x)\le\varepsilon_f\) (feas_eps, por defecto \(10^{-6}\)). Denotamos \(\mathcal{F}\) al conjunto de individuos factibles.

2.4 Operador de reparación

Tras cada operador que modifica \(x\) se aplica

\[ \mathcal{R}(x)=\operatorname{round}_{\mathcal{I}}\big(\operatorname{clip}(x,\,l,\,u)\big), \]

donde \(\operatorname{clip}\) proyecta cada coordenada a \([l_i,u_i]\) y \(\operatorname{round}_{\mathcal{I}}\) redondea al entero más cercano únicamente las coordenadas \(i\in\mathcal{I}\). El redondeo no vuelve a proyectar, por lo que conviene que los límites de las variables enteras sean enteros.

3 Representación

3.1 Individuos

Cada célula es una tupla

\[ \text{ind}=(x,\ \sigma,\ e,\ w), \]

con posición \(x\in\mathbb{R}^n\) (genotipo), paso de mutación auto-adaptativo \(\sigma>0\), reserva de energía \(e\) y vector de pesos de escalarización \(w\in\Delta^{M-1}\) (símplex; \(w=1\) si \(M=1\)). En cada evaluación se calculan además: \(F(x)\), \(cv(x)\), el rango \(r\), la distancia de aglomeración \(c\), el fitness compuesto \(\Phi\) y el nicho asignado.

3.2 Poblaciones e inicialización

En la generación \(t\) existen \(\mathcal{C}_t\) y \(\mathcal{L}_t\), con tamaños iniciales \(N_C\) (n_cancer) y \(N_L\) (n_lymph). La inicialización es uniforme sobre la caja y luego se repara:

\[ x_i\sim\mathcal{U}(l_i,u_i)\ \ (i=1,\dots,n),\qquad x\leftarrow\mathcal{R}(x), \]

de modo que las coordenadas enteras resultan de redondear una variable uniforme continua. Además \(\sigma_0=\sigma_{\text{init}}\), \(e_0=e_{\text{init}}\) y, si \(M>1\), \(w\sim\operatorname{Dirichlet}(\mathbf{1}_M)\).

4 Evaluación y ordenamiento

Esta sección unifica el caso mono y multiobjetivo, con y sin restricciones. Todo operador de selección posterior depende únicamente del fitness compuesto \(\Phi\) definido aquí.

4.1 Dominancia restringida (regla de Deb)

Con tolerancia \(\delta=10^{-12}\), el individuo \(a\) domina al individuo \(b\) (\(a\succ b\)) si:

\[ a\succ b\iff \begin{cases} a,b\in\mathcal{F}\ \wedge\ \tilde f_m(a)\le \tilde f_m(b)+\delta\ \forall m\ \wedge\ \exists m:\ \tilde f_m(a)< \tilde f_m(b)-\delta & \text{(Pareto, ambos factibles)}\\[4pt] a\in\mathcal{F}\ \wedge\ b\notin\mathcal{F} & \text{(factible domina a infactible)}\\[4pt] a,b\notin\mathcal{F}\ \wedge\ cv(a)<cv(b)-\delta & \text{(entre infactibles, menor violación)} \end{cases} \]

4.2 Frentes de no-dominancia

El ordenamiento rápido (tipo NSGA-II) asigna el rango \(r=0\) a los individuos no dominados; el frente \(r+1\) está formado por los individuos cuyos dominadores pertenecen todos a frentes \(\le r\). Si existe al menos un individuo factible, el frente \(r=0\) contiene únicamente individuos factibles.

4.3 Distancia de aglomeración

Dentro de un frente \(\mathcal{F}_r\) con \(|\mathcal{F}_r|>2\), para cada objetivo \(m\) se ordenan los individuos por \(f_m\); los extremos reciben \(c=\infty\) y los interiores acumulan

\[ c_i\;{+}{=}\;\frac{\tilde f_m(i^{+})- \tilde f_m(i^{-})}{ \tilde f_m^{\max}- \tilde f_m^{\min}}, \]

donde \(i^{\pm}\) son los vecinos inmediatos en el orden (se omite el objetivo si el denominador es menor que \(10^{-12}\)). Si \(|\mathcal{F}_r|\le 2\), todos reciben \(c=\infty\).

4.4 Fitness compuesto \(\Phi\)

Sea \(R_{\max}=\max\{(\text{nº de frentes})-1,\ 1\}\) y \(\tilde c_i=c_i/c_{\max}^{(r)}\) la distancia normalizada dentro de su frente (\(\tilde c_i=1\) si \(c_i=\infty\); \(c_{\max}^{(r)}\) es el máximo finito del frente). Para \(M\ge2\):

\[ \Phi_i \;=\; 1-\frac{r_i}{R_{\max}}\;+\;0.05\,\tilde c_i . \]

El rango domina la comparación; la aglomeración actúa solo como desempate que favorece la diversidad.

Caso mono-objetivo (\(M=1\)). Se evita el costo \(O(N^2)\) ordenando por la clave lexicográfica

\[ \kappa_i=\begin{cases}(0,\ \tilde f_1(x_i)) & i\in\mathcal{F}\\ (1,\ cv(x_i)) & i\notin\mathcal{F}\end{cases} \]

y, si \(\rho_i\in\{0,\dots,N-1\}\) es la posición de \(i\) en dicho orden,

\[ \Phi_i=1-\frac{\rho_i}{\max\{N-1,\,1\}} . \]

Esto reproduce exactamente la regla de Deb para un objetivo.

4.5 Fitness normalizado por población

Para el duelo y la inanición se usa la normalización min–max dentro de cada población:

\[ \hat\phi_i=\frac{\Phi_i-\Phi_{\min}}{\Phi_{\max}-\Phi_{\min}}\in[0,1], \]

con \(\hat\phi_i=0.5\) si \(\Phi_{\max}-\Phi_{\min}<10^{-12}\).

5 Operadores

5.1 Mutación adaptativa

Se calculan \(\tau=1/\sqrt{2\sqrt{n}}\) y \(\tau'=1/\sqrt{2n}\). Para cada individuo, con \(\xi,\xi'\sim\mathcal{N}(0,1)\) independientes:

\[ \sigma\leftarrow\max\Big\{\sigma_{\min},\ \sigma\,\exp\big(\tau'\xi+\tau\xi'\big)\Big\}. \]

• Variables continuas (\(i\notin\mathcal{I}\)): \(\;x_i\leftarrow x_i+\sigma\,\eta_i\), con \(\eta_i\sim\mathcal{N}(0,1)\).
• Variables enteras (\(i\in\mathcal{I}\)): \(\;x_i\leftarrow x_i+\zeta_i\), con \(\zeta_i\sim\mathcal{U}\{-m,\dots,m\}\) y \(m=\max\{1,\operatorname{round}(\sigma)\}\).

Finalmente \(x\leftarrow\mathcal{R}(x)\). Así el tamaño del paso evoluciona junto con el individuo y, en variables enteras, la magnitud del salto escala con \(\sigma\).

5.2 Quimiotaxis (movimiento guiado)

Escalar de guía. Para obtener un gradiente se define un escalar “mayor es mejor” que prioriza la factibilidad:

\[ S(x;w)=\begin{cases}-\,cv(x) & cv(x)>\varepsilon_f\quad\text{(infactible: reducir violación)}\\[2pt] -\,w^{\top} \tilde F(x) & cv(x)\le\varepsilon_f\quad\text{(factible: mejorar objetivos)}\end{cases} \]

Para \(M>1\), cada individuo usa sus propios pesos \(w\) (escalarización tipo MOEA/D), de modo que distintas células persiguen distintas zonas del frente de Pareto.

Gradiente por diferencias finitas (hacia adelante, paso \(\varepsilon_{fd}\)), anulado en variables enteras:

\[ \hat g_i=\begin{cases}\dfrac{S(x+\varepsilon_{fd}e_i;w)-S(x;w)}{\varepsilon_{fd}} & i\notin\mathcal{I}\\[8pt] 0 & i\in\mathcal{I}\end{cases} \qquad \bar g=\frac{\hat g}{\lVert\hat g\rVert}. \]

Movimiento. Con \(d\sim\mathcal{N}(0,I)/\lVert\mathcal{N}(0,I)\rVert\) (dirección aleatoria unitaria):

\[ \beta=\min\{1,\ \beta_0+\beta_f\,\varphi_f\},\qquad x\leftarrow\mathcal{R}\Big(x+s_0\,\frac{\sigma}{\sigma_{\text{init}}}\big[(1-\beta)\,d+\beta\,\bar g\big]\Big), \]

donde \(\beta_0\) es el peso guiado base, \(s_0\) el paso base y \(\varphi_f\in[0,1]\) la activación de fibroblastos del nicho del individuo (sección 5.8; \(\varphi_f=0\) si están desactivados).

5.3 Duelo biológico

Para cada célula cancerosa \(c\) (en el orden de la lista), se busca el linfocito no emparejado más cercano dentro del radio de interacción \(r_d\):

\[ l^*(c)=\operatorname*{arg\,min}_{l\in\mathcal{L}\ \text{libre},\ d(c,l)<r_d} d(c,l),\qquad d(c,l)=\lVert x_c-x_l\rVert . \]

Si no existe, \(c\) sobrevive sin duelo. En caso contrario, se calcula la afinidad gaussiana y la probabilidad de muerte del cáncer:

\[ A(c,l)=\exp\!\Big(-\frac{d(c,l)^2}{2\,(r_d/2)^2}\Big),\qquad z=k\,\big(A\,\hat\phi_l-\hat\phi_c\big), \]

\[ P(\text{muere } c)=\operatorname{sigm}(z)\,\big(1-\kappa\,\varphi_f(c)\big),\qquad \operatorname{sigm}(z)=\frac{1}{1+e^{-z}} . \]

El factor \((1-\kappa\varphi_f)\) es el escudo inmune parcial de los fibroblastos (\(\kappa=\) fibro_shield). Con un sorteo \(U\sim\mathcal{U}(0,1)\):

• Si \(U<P(\text{muere }c)\): \(c\) es eliminada y \(e_l\leftarrow\min\{e_l+0.3,\ 2\}\).
• En caso contrario: \(e_c\leftarrow\min\{e_c+0.2,\ 2\}\), \(e_l\leftarrow e_l-0.4\); \(c\) sobrevive y \(l\) sobrevive solo si \(e_l>0\).

Los linfocitos sin pareja sobreviven. Así, el cáncer fuerte (alto \(\hat\phi_c\)) resiste, y el linfocito fuerte y cercano lo elimina.

k <- 1.5; kappa <- 0.5
phi_c <- seq(0, 1, length.out = 200)
p_sin <- sigm(k * (1 * 0.5 - phi_c))
p_con <- p_sin * (1 - kappa * 1)
plot(phi_c, p_sin, type = "l", lwd = 2, col = "#E63946", ylim = c(0, 1),
     xlab = "Fitness normalizado de la célula cancerosa",
     ylab = "Probabilidad de muerte")
lines(phi_c, p_con, lwd = 2, lty = 2, col = "#2E86AB")
legend("topright", bty = "n", lwd = 2, lty = c(1, 2), col = c("#E63946", "#2E86AB"),
       legend = c("Sin fibroblastos", "Con fibroblastos (phi_f = 1)"))

Probabilidad de muerte de una célula cancerosa en el duelo (afinidad A = 1, fitness normalizado del linfocito = 0.5, k = 1.5), con y sin escudo de fibroblastos (kappa = 0.5, phi_f = 1).

5.4 Cruce (BLX-\(\alpha\))

Se realizan \(\lfloor|\mathcal{P}|/2\rfloor\) intentos por población \(\mathcal{P}\); cada uno ocurre con probabilidad \(p_c\). Se eligen dos padres \(p_1\neq p_2\) de forma uniforme y, con \(\lambda_i\sim\mathcal{U}(-\alpha,\,1+\alpha)\) independiente por coordenada:

\[ x^{\text{hijo}}=\mathcal{R}\big(\lambda\odot x^{p_1}+(1-\lambda)\odot x^{p_2}\big),\quad \sigma^{\text{hijo}}=\tfrac12(\sigma^{p_1}+\sigma^{p_2}),\quad e^{\text{hijo}}=\tfrac12(e^{p_1}+e^{p_2}),\quad w^{\text{hijo}}=\frac{w^{p_1}+w^{p_2}}{\mathbf{1}^{\top}(w^{p_1}+w^{p_2})}. \]

El cruce se restringe a individuos de la misma población. En el ciclo principal solo se aplica a la población cancerosa (los linfocitos no se cruzan). Los hijos se añaden a la población y no tienen fitness asignado hasta la siguiente evaluación.

5.5 Nichos y proliferación

Asignación de nichos. Se recorre la población en orden; cada individuo aún sin nicho funda uno nuevo y adopta en él a todos los individuos posteriores sin nicho a distancia menor que \(r_n\) de él (agrupamiento codicioso por semilla).

Crecimiento logístico. Para un nicho con \(n_\ell\) miembros y capacidad de carga \(K\), el número de descendientes es

\[ n_{\text{new}}=\max\Big\{0,\ \operatorname{round}\Big(r\,n_\ell\,\Big(1-\frac{n_\ell}{K}\Big)\Big)\Big\}. \]

Los padres se toman en orden de fitness decreciente (cíclicamente) y cada hijo es

\[ x^{\text{hijo}}=\mathcal{R}\big(x^{\text{padre}}+\mathcal{N}(0,\ (0.5\,\sigma^{\text{padre}})^2 I)\big),\quad \sigma^{\text{hijo}}=\sigma^{\text{padre}},\quad e^{\text{hijo}}=0.7\,e_{\text{init}},\quad w^{\text{hijo}}=w^{\text{padre}}. \]

Si la población resultante excede \(N_C\) (su tamaño inicial, n_cancer), se reevalúa y se conservan los \(N_C\) individuos con mayor \(\Phi\) (mejores frentes primero, y a igual frente mayor aglomeración). Solo la población cancerosa prolifera; los linfocitos no se reproducen.

r <- 0.3; K <- 10
n <- 0:20
nuevos <- pmax(0, round(r * n * (1 - n / K)))
plot(n, nuevos, type = "h", lwd = 6, lend = 1, col = "#2E86AB",
     xlab = "Tamaño del nicho", ylab = "Descendientes por generación")
abline(v = K, lty = 2, col = "grey40")

Descendientes por generación en un nicho según su tamaño (r = 0.3, K = 10). El crecimiento se anula al alcanzar la capacidad de carga.

5.6 Metástasis (vuelos de Lévy)

Sea \(n_\ell(i)\) el tamaño del nicho del individuo \(i\). Cada individuo de la población cancerosa “metastatiza” con probabilidad dependiente del hacinamiento:

\[ p_{\text{meta}}(i)=p_0\cdot\min\Big\{1,\ \frac{n_\ell(i)}{K}\Big\}. \]

Si ocurre, \(x\leftarrow\mathcal{R}(x+L)\), \(\sigma\leftarrow\sigma_{\text{init}}\) y, si \(M>1\), \(w\sim\operatorname{Dirichlet}(\mathbf{1}_M)\) (una nueva “estrategia” en el nuevo sitio). El salto \(L\in\mathbb{R}^n\) sigue el algoritmo de Mantegna con exponente \(\beta_L\in(1,2]\) y escala \(s_L\):

\[ \sigma_u=\left[\frac{\Gamma(1+\beta_L)\,\sin(\pi\beta_L/2)}{\Gamma\!\big(\tfrac{1+\beta_L}{2}\big)\,\beta_L\,2^{(\beta_L-1)/2}}\right]^{1/\beta_L},\qquad u_i\sim\mathcal{N}(0,\sigma_u^2),\quad v_i\sim\mathcal{N}(0,1),\qquad L_i=s_L\,\frac{u_i}{|v_i|^{1/\beta_L}} . \]

Las colas pesadas permiten tanto saltos cortos como colonizaciones lejanas, lo que ayuda a escapar de óptimos locales.

set.seed(1)
beta <- 1.5
sigma_u <- (gamma(1 + beta) * sin(pi * beta / 2) /
              (gamma((1 + beta) / 2) * beta * 2^((beta - 1) / 2)))^(1 / beta)
u <- rnorm(20000, 0, sigma_u); v <- rnorm(20000)
salto <- abs(u / abs(v)^(1 / beta))
hist(pmin(salto, 8), breaks = 60, col = "#2E86AB", border = "white", main = "",
     xlab = "Magnitud del salto", ylab = "Frecuencia")

Distribución de la magnitud del salto de Lévy (beta = 1.5, escala 1, truncada en 8 para visualización).

5.7 Muerte por inanición

Cada generación, la energía de cada individuo de una población \(\mathcal{P}\) se actualiza con un balance de mantenimiento, ganancia por calidad, bonificación de fibroblastos y competencia local:

\[ e_i\leftarrow e_i-e_m+e_g\,\hat\phi_i+e_f\,\varphi_f(i)-e_c\,\rho_i,\qquad \rho_i=\#\{\,j\in\mathcal{P},\ j\ne i:\ \lVert x_i-x_j\rVert<r_n\,\}. \]

El individuo muere si \(e_i\le 0\). Aquí \(e_m\) es el costo de mantenimiento, \(e_g\) la ganancia máxima por fitness, \(e_f\) la bonificación por fibroblastos (\(\texttt{fibro\_energy\_boost}\)) y \(e_c\) el costo por vecino cercano (competencia por recursos).

5.8 Fibroblastos (microambiente de soporte)

El campo de fibroblastos se calcula una vez por generación sobre la población cancerosa y se aplica solo a ella (\(\varphi_f=0\) para los linfocitos). Sea \(\{\mathcal{N}_\ell\}\) la partición en nichos de \(\mathcal{C}_t\) (sección 5.5). La calidad del nicho y su promedio global son

\[ Q_\ell=\frac{1}{|\mathcal{N}_\ell|}\sum_{i\in\mathcal{N}_\ell}\Phi_i,\qquad \bar Q=\frac{1}{L}\sum_{\ell=1}^{L}Q_\ell, \]

y la activación del nicho es una sigmoide de su calidad relativa:

\[ \varphi_f(\ell)=\operatorname{sigm}\!\big(k_f\,(Q_\ell-\bar Q)\big)\in(0,1). \]

Los nichos mejores que el promedio (\(Q_\ell>\bar Q\)) reciben \(\varphi_f>0.5\) y los peores \(\varphi_f<0.5\). La activación produce tres efectos biológicamente motivados:

1. Nutrientes (inanición): bonificación de energía \(+e_f\,\varphi_f\) (sección 5.7).
2. ECM remodelada (quimiotaxis): mayor peso del movimiento guiado, \(\beta=\min\{1,\beta_0+\beta_f\varphi_f\}\) (sección 5.2).
3. Inmunosupresión (duelo): reducción de la probabilidad de muerte por el factor \((1-\kappa\varphi_f)\) (sección 5.3).

Como \(\varphi_f\) proviene del clustering de nichos que ya se calcula para la proliferación y la metástasis, el costo adicional es \(O(N)\) por generación.

dQ <- seq(-1, 1, length.out = 200)
cols <- c("#999999", "#2E86AB", "#E63946"); kfs <- c(1, 3, 8)
plot(NA, xlim = c(-1, 1), ylim = c(0, 1),
     xlab = "Calidad del nicho menos calidad media",
     ylab = "Activación de fibroblastos")
for (i in seq_along(kfs)) lines(dQ, sigm(kfs[i] * dQ), lwd = 2, col = cols[i])
abline(h = 0.5, v = 0, lty = 3, col = "grey60")
legend("topleft", bty = "n", lwd = 2, col = cols, legend = paste("k_f =", kfs))

Activación de fibroblastos según la calidad relativa del nicho, para distintas ganancias k_f.

6 Selección de la solución final y registro

El algoritmo no mantiene un archivo élite: la solución devuelta se elige entre los individuos vivos de la población cancerosa al terminar la última generación. Un buen individuo que muera en el duelo o por inanición puede, por tanto, perderse, y la curva de mejor valor no es necesariamente monótona.

Sea \(\mathcal{F}_C\) el conjunto de individuos factibles de la población cancerosa final.

• \(M=1\). Si \(\mathcal{F}_C\neq\emptyset\) se devuelve \(\arg\min_{i\in\mathcal{F}_C} f_1(x_i)\), con indicador de factibilidad verdadero. En otro caso, el individuo de menor \(cv\), con indicador falso.
• \(M>1\). Si \(\mathcal{F}_C\neq\emptyset\) se devuelve el individuo factible de mayor fitness compuesto \(\Phi\) (mejor frente y, dentro de él, mayor aglomeración). En otro caso, el de menor \(cv\), con indicador falso.

La salida es la terna (x, F(x), factible). Para \(M>1\) es un único punto, no una aproximación del frente de Pareto.

En cada generación se calcula, además, un valor de seguimiento: el valor \(f_1\) del mejor factible según \(\Phi\) para \(M=1\) (o su \(\Phi\) para \(M>1\)), o \(-cv\) del menos infactible si no hay factibles; también los tamaños de ambas poblaciones. Ver la sección 12 sobre este registro.

7 Ciclo evolutivo completo

Se crean \(\mathcal{C}_0\) y \(\mathcal{L}_0\) (sección 3.2) y se evalúan. Dado el estado \((\mathcal{C}_t,\mathcal{L}_t)\), una generación realiza, en este orden (CancerGA.run):

1. Campo de fibroblastos \(\varphi_f\) sobre \(\mathcal{C}_t\) (sección 5.8), si está activado; usa el \(\Phi\) de la última evaluación.
2. Mutación adaptativa en ambas poblaciones (sección 5.1).
3. Quimiotaxis en ambas poblaciones (sección 5.2).
4. Evaluación de ambas poblaciones: \(F\), \(cv\), rangos, aglomeración y \(\Phi\) (sección 4).
5. Duelo biológico (sección 5.3).
6. Cruce BLX-\(\alpha\), solo en la población cancerosa (sección 5.4).
7. Proliferación por nichos con tope \(N_C\), solo cancerosa (sección 5.5).
8. Metástasis, solo cancerosa (sección 5.6).
9. Evaluación de ambas poblaciones.
10. Inanición en ambas poblaciones (sección 5.7).
11. Reabastecimiento: si \(|\mathcal{C}|<5\) se completa hasta \(N_C\) con individuos aleatorios (como en 3.2); si \(|\mathcal{L}|<2\), hasta \(N_L\).
12. Evaluación final y registro de la generación.

Se repite durante \(T\) generaciones (n_generations). No hay criterio de paro por convergencia y la extinción de una población no termina la ejecución (se repone en el paso 11).

8 Ejemplo de comportamiento sobre la función Rastrigin

8.1 Paisaje multimodal

La función de Rastrigin en dos dimensiones,

\[ f(x,y)=20+x^2-10\cos(2\pi x)+y^2-10\cos(2\pi y),\qquad (x,y)\in[-5.12,\,5.12]^2, \]

tiene mínimo global \(f(0,0)=0\) y del orden de cien mínimos locales distribuidos en una rejilla casi regular. Es el prototipo de paisaje multimodal: un método puramente explotador queda atrapado en la primera cuenca que encuentra; CancerGA debe combinarlo con exploración (mutación adaptativa, metástasis de Lévy) para saltar entre cuencas.

8.2 Simulación ilustrativa

Para que este documento se compile sin depender del repositorio, se incluye aquí una simulación en R que reproduce los ingredientes esenciales de CancerGA sobre este paisaje: mutación gaussiana, cruce BLX-\(\alpha\), selección por torneo y metástasis con vuelos de Lévy. No modela el duelo ni los fibroblastos —eso requeriría el algoritmo completo—, pero sí comparte la mecánica de exploración-explotación que interesa visualizar. La figura muestra únicamente la población cancerosa (soluciones candidatas), ya que es la que converge. Para ejecutar el algoritmo real sobre Rastrigin véase main.py y el Apéndice B.

set.seed(2026)

rastrigin <- function(X) {
  X <- matrix(X, ncol = 2)
  20 + (X[,1]^2 - 10*cos(2*pi*X[,1])) + (X[,2]^2 - 10*cos(2*pi*X[,2]))
}

L <- -5.12; U <- 5.12
n_pop <- 80; n_elite <- 20; n_children <- n_pop - n_elite

x <- cbind(runif(n_pop, L, U), runif(n_pop, L, U))
snap_gens <- c(0, 10, 30, 100, 300, 1000)
snaps <- list(`0` = x)
best_fit <- numeric(max(snap_gens) + 1)
best_fit[1] <- min(rastrigin(x))

gen <- 0
for (target in snap_gens[-1]) {
  while (gen < target) {
    gen <- gen + 1
    # Mutación gaussiana (simplificación de la mutación adaptativa)
    x <- pmin(pmax(x + cbind(0.3*rnorm(n_pop), 0.3*rnorm(n_pop)), L), U)
    fit <- rastrigin(x)
    ord <- order(fit)
    elite <- x[ord[1:n_elite], , drop = FALSE]
    # Cruce BLX-alpha (alpha = 0.3): cada hijo muestrea dos padres independientes
    i1 <- sample(n_pop, n_children, replace = TRUE)
    i2 <- sample(n_pop, n_children, replace = TRUE)
    p1 <- x[i1, , drop = FALSE]
    p2 <- x[i2, , drop = FALSE]
    lam <- matrix(runif(2*n_children, -0.3, 1.3), ncol = 2)
    child <- lam*p1 + (1 - lam)*p2
    child <- pmin(pmax(
      child + cbind(0.15*rnorm(n_children), 0.15*rnorm(n_children)), L), U)
    x <- rbind(elite, child)            # 20 + 60 = 80 filas, como n_pop
    # Metástasis con vuelos de Lévy
    i <- which(runif(n_pop) < 0.03)
    if (length(i) > 0) {
      beta <- 1.5
      su <- (gamma(1+beta)*sin(pi*beta/2) /
             (gamma((1+beta)/2)*beta*2^((beta-1)/2)))^(1/beta)
      u <- rnorm(length(i)*2, 0, su); v <- rnorm(length(i)*2)
      jump <- matrix(u/abs(v)^(1/beta), ncol = 2) * 0.5
      x[i, ] <- pmin(pmax(x[i, , drop = FALSE] + jump, L), U)
    }
    best_fit[gen + 1] <- min(rastrigin(x))
  }
  snaps[[as.character(target)]] <- x
}

gx <- seq(L, U, length.out = 250)
Z  <- outer(gx, gx, function(a, b)
  20 + a^2 - 10*cos(2*pi*a) + b^2 - 10*cos(2*pi*b))

pal <- hcl.colors(64, "YlOrRd", rev = TRUE)
op  <- par(mfrow = c(2, 3), mar = c(1.8, 1.8, 2, 0.8), oma = c(0.5, 0.5, 2, 0.5))

for (g in snap_gens) {
  image(gx, gx, Z, col = pal, xlab = "", ylab = "", axes = FALSE,
        main = paste0("Generación ", g), useRaster = TRUE)
  box()
  contour(gx, gx, Z, add = TRUE, col = "grey55",
          drawlabels = FALSE, nlevels = 15)
  points(snaps[[as.character(g)]], pch = 19, col = "#7B2CBF", cex = 0.7)
  points(0, 0, pch = 3, col = "black", cex = 1.4, lwd = 2)
}
mtext("Rastrigin 2D — evolución de la población cancerosa",
      side = 3, outer = TRUE, line = 0.4, cex = 1.05, font = 2)

Evolución de la población cancerosa sobre el paisaje de Rastrigin. Cada panel muestra las 80 células (puntos rojos) tras la generación indicada; la cruz negra marca el óptimo global (0,0). Los colores oscuros corresponden a valores bajos de f.

par(op)

plot(0:max(snap_gens), pmax(best_fit[seq_len(max(snap_gens) + 1)], 1e-8),
     type = "l", lwd = 1.6, col = "#2E86AB", log = "y",
     xlab = "Generación", ylab = "Mejor f (escala log)",
     main = "Convergencia del mejor individuo")
grid(col = "grey85")
points(snap_gens, pmax(best_fit[snap_gens + 1], 1e-8),
       pch = 19, col = "#E63946", cex = 1.2)

Convergencia del mejor individuo (escala logarítmica, valores acotados inferiormente por 1e-8 solo para poder usar escala logarítmica). Los puntos rojos corresponden a las generaciones mostradas en la figura anterior.

La primera figura permite seguir cómo una población que parte uniforme sobre el dominio se concentra en cuencas del paisaje: la mutación y el cruce explotan la cuenca en curso, y los saltos de Lévy (metástasis) son los que permiten abandonar cuencas subóptimas. La segunda figura registra el mejor valor por generación; se espera una caída rápida al inicio (convergencia dentro de una cuenca) seguida de escalones espaciados cuando un salto alcanza una cuenca mejor. Los valores concretos dependen de la semilla (set.seed(2026)) y se obtienen al compilar el documento. Esta simulación es solo ilustrativa y no constituye un resultado de desempeño de CancerGA.

9 Dinámica poblacional agregada (interpretación heurística)

El algoritmo es basado en individuos; no obstante, su comportamiento medio se interpreta cualitativamente con un modelo tipo Lotka–Volterra con capacidad de carga:

\[ \dot N_C=r\,N_C\Big(1-\frac{N_C}{K}\Big)-\alpha\,N_C N_L,\qquad \dot N_L=\beta\,N_C N_L-\delta\,N_L , \]

donde \(r\) es la tasa de proliferación, \(\alpha\) la letalidad del duelo, \(\beta\) la ganancia de linfocitos por duelos ganados y \(\delta\) su mortalidad. Es un marco conceptual para razonar sobre el equilibrio entre exploración (proliferación, metástasis) y presión selectiva (duelo, inanición): no se deriva rigurosamente del modelo individual ni garantiza convergencia.

10 Complejidad computacional

Sean \(N\) el tamaño de población, \(n\) el número de variables y \(n_c\) el de variables continuas.

Componente	Costo por generación
Quimiotaxis (gradiente)	\(O\big((n_c+1)\,N\big)\) evaluaciones de \(S\)
Evaluaciones de objetivo/restricciones	\(O(N)\) por evaluación; 3 evaluaciones completas por ciclo (más una si hay truncamiento en la proliferación)
Ordenamiento, \(M=1\)	\(O(N\log N)\)
Ordenamiento por dominancia, \(M\ge2\)	\(O(M\,N^2)\) (vectorizado)
Duelo	\(O(\lvert\mathcal{C}\rvert\,\lvert\mathcal{L}\rvert\,n)\)
Nichos / metástasis	\(O(N^2\,n)\) en el peor caso
Inanición (densidad local)	\(O(N^2\,n)\)
Fibroblastos	\(O(N)\) adicional (reutiliza los nichos)

En problemas con función objetivo costosa, el término dominante es la quimiotaxis, proporcional a \((n_c+1)N\) evaluaciones por generación.

11 Parámetros y valores por defecto

Valores del constructor de CancerGA (algorithm.py).

Símbolo	Parámetro en el código	Valor	Rol
\(N_C,\ N_L\)	n_cancer, n_lymph	60, 20	Tamaños iniciales
\(\sigma_{\text{init}},\ \sigma_{\min}\)	sigma_init, sigma_min	0.5, \(10^{-5}\)	Paso de mutación
\(\beta_0\)	chemo_beta	0.3	Peso guiado base
\(s_0\)	chemo_step	0.2	Paso base de quimiotaxis
\(\varepsilon_{fd}\)	chemo_fd_eps	\(10^{-5}\)	Paso de diferencias finitas
\(r_d\)	duel_radius	1.0	Radio de interacción del duelo
\(k\)	duel_k	2.0	Ganancia de la sigmoide del duelo
\(p_c\)	crossover_rate	0.7	Probabilidad de cruce
\(\alpha\)	blx_alpha	0.5	Extensión de BLX-\(\alpha\)
\(r_n\)	niche_radius	1.5	Radio de nicho
\(K\)	niche_capacity	15	Capacidad de carga del nicho
\(r\)	growth_rate	0.4	Tasa de proliferación
\(p_0\)	p_meta0	0.08	Probabilidad base de metástasis
\(\beta_L,\ s_L\)	levy_beta, levy_scale	1.5, 0.3	Vuelo de Lévy
\(e_m,\ e_g,\ e_c\)	e_maintenance, e_gain, e_comp	0.1, 0.3, 0.05	Balance de energía
\(e_{\text{init}}\)	e_init	1.0	Energía inicial
\(\varepsilon_h\)	eq_tol	\(10^{-4}\)	Tolerancia de igualdades
\(\varepsilon_f\)	feas_eps	\(10^{-6}\)	Tolerancia de factibilidad
\(k_f\)	fibro_k	3.0	Ganancia de reclutamiento
\(e_f\)	fibro_energy_boost	0.1	Nutrientes extra
\(\kappa\)	fibro_shield	0.3	Escudo inmune máximo
\(\beta_f\)	fibro_beta_boost	0.4	Aumento máximo del peso guiado
—	use_fibroblasts	True	Activa el microambiente
—	seed	None	Semilla del generador aleatorio

12 Limitaciones

12.1 Del método

• Es una metaheurística: no hay garantía de convergencia al óptimo global y el comportamiento depende de los parámetros de la sección 11.
• Con restricciones de igualdad y tolerancia pequeña, la región factible es muy delgada y encontrarla puede requerir más generaciones o población.
• En \(M=1\) el fitness \(\Phi\) es ordinal: depende solo de la posición en el orden, no de la magnitud de las diferencias entre objetivos, y el duelo y la inanición heredan esa propiedad.
• Algunos problemas de la familia MW tienen regiones factibles casi inexistentes por diseño; pueden requerir muchas más generaciones.

12.2 De la implementación actual (v0.1.0)

• Restricciones de benchmarks. Las funciones CTP, CF y MW devuelven (objs, g), pero el algoritmo ignora g. Deben pasarse como constraints (Apéndice B).
• Sin maximización nativa. Para maximizar, se entrega \(-f\).
• Salida de un solo punto. No hay archivo élite, frente de Pareto ni punto codo. plot_pareto_front espera un archivo y un punto codo calculados externamente.
• Registro no devuelto. run() acumula por generación el mejor valor y los tamaños poblacionales, pero retorna solo la terna final, de modo que plot_convergence requiere acceder a ese historial por otra vía.
• Nichos y fibroblastos. La tabla \(\varphi_f\) se calcula al inicio de la generación; niche_proliferation y metastasis reasignan después los identificadores de nicho. Cuando la inanición consulta \(\varphi_f\) (efecto de energía), usa identificadores ya reasignados, por lo que ese término corresponde solo aproximadamente al nicho original. Los efectos sobre la quimiotaxis y el duelo no se ven afectados.
• Linfocitos. Mutan, se mueven, duelan, mueren por inanición y se reponen, pero no se cruzan, no proliferan ni metastatizan.
• Extinción. Una población colapsada se repone con individuos aleatorios (paso 11 del ciclo) en lugar de terminar la ejecución.
• Quimiotaxis. Los puntos de la diferencia finita no se reparan, por lo que pueden quedar fuera de los límites.
• Benchmarks. Las implementaciones de CTP y CF difieren en parte de las definiciones originales y presentan los problemas indicados en el Apéndice A.

Apéndice A: familias de benchmark

Todas se minimizan. Las funciones CTP y CF devuelven (objs, g) con las restricciones en la forma \(g(x)\le 0\); como se indicó en la sección 12, el algoritmo no las aplica automáticamente (ver Apéndice B). Las fórmulas de este apéndice son las implementadas en el repositorio (benchmarks/) y no siempre coinciden con las definiciones originales (Deb et al., 2001; Zhang et al., 2008). Antes de comparar resultados con la literatura, conviene contrastarlas con una implementación de referencia. En CF los índices son de base 0, como en el código (\(x_0\) es la primera variable).

A.0 Mono-objetivo

Función (single_objective.py)	Definición
sphere	\(\sum_i x_i^2\)
rastrigin	\(10n+\sum_i\big(x_i^2-10\cos(2\pi x_i)\big)\)
rosenbrock	\(\sum_{i}\big[100\,(x_{i+1}-x_i^2)^2+(1-x_i)^2\big]\)
ackley	\(-20\,e^{-0.2\sqrt{\sum x_i^2/n}}-e^{\sum\cos(2\pi x_i)/n}+20+e\)
constrained_mixed_example	\(x_0^2+x_1^2+3x_2\)

A.1 Familia CTP

Con \(x\in[0,1]^n\), todas las funciones CTP del repositorio comparten

\[ f_1=x_0,\qquad \gamma=1+\frac{9}{n-1}\sum_{i=1}^{n-1}x_i,\qquad f_2=\gamma\Big(1-\sqrt{f_1/\gamma}\Big). \]

• CTP1: \(g_1=f_2-0.858+0.541\,f_1\le 0\) y \(g_2=f_2-0.728+0.295\,f_1\le 0\) (restricciones lineales; la definición original usa \(f_2\ge a_j e^{-b_jf_1}\) con \(f_2=\gamma e^{-f_1/\gamma}\)).
• CTP2–CTP7: con \(\theta=-0.2\pi\),

\[ c(x)=(f_2-f_1)\sin\theta+(f_2+f_1)\cos\theta-a\,\big|\cos\big(e\,(f_2-f_1)\big)\big|^{b}-d, \]

y se devuelve la restricción \(g=-c\le 0\), con

Problema	\(a\)	\(b\)	\(d\)	\(e\)
CTP2	0.2	10	1	6
CTP3	0.1	10	1	0.5
CTP4	0.75	0.5	0.75	0.5
CTP5	0.1	0.5	0.75	2
CTP6	40	0.5	1	6
CTP7	20	2	1	6

• CTP8: dos restricciones \(g_1=-c_1\le0\), \(g_2=-c_2\le0\), con

\[ c_1=(f_2-f_1)\sin(-0.2\pi)+(f_2+f_1)\cos(-0.2\pi)-0.1,\quad c_2=(f_2-f_1)\sin(0.1\pi)+(f_2+f_1)\cos(0.1\pi)-0.1 . \]

A.2 Familia CF

Por defecto \(n=10\). Conjuntos de índices (base 0):

• Bi-objetivo (CF1–CF7): \(J_1=\{2,4,6,\dots\}\) y \(J_2=\{1,3,5,\dots\}\), con índices \(\le n-1\).
• Tri-objetivo (CF8–CF10): \(J_1=\{2,5,8,\dots\}\), \(J_2=\{3,6,9,\dots\}\) y \(J_3=\{1,4,7,\dots\}\), con índices \(\le n-1\).

Notación: \(\Sigma_k[\phi]=\dfrac{2}{|J_k|}\sum_{i\in J_k}\phi(y_i)\) y \(\psi_c(y)=2y^2-\cos(4\pi y)+1\).

Bi-objetivo (CF1–CF7)

Para CF2, CF4, CF5, CF6 y CF7: \(y_i=x_i-\sin\!\big(6\pi x_0+\tfrac{(i+1)\pi}{n}\big)\) para \(i\ge1\). Para CF3: \(y_i=x_i-x_0^{\,0.5\,(1+3(i-1)/(n-2))}\) para \(i\ge1\). En CF3 se usa además \(\Xi_k=\dfrac{2}{|J_k|}\Big(4\sum_{J_k}y_i^2-2\prod_{J_k}\cos\!\tfrac{20\pi y_i}{\sqrt{i+1}}+2\Big)\).

Problema	\(f_1\)	\(f_2\)	Restricción (se exige \(\ge 0\))
CF1	\(x_0+\Sigma_1[y^2]\)	\(1-x_0+\Sigma_2[y^2]\)	\(f_2+f_1^2-1\)
CF2	\(x_0+\Sigma_1[y^2]\)	\(1-\sqrt{x_0}+\Sigma_2[y^2]\)	\(f_2+\sqrt{f_1}-\sin\!\big(2\pi(f_2-\sqrt{f_1}+1)\big)-1\)
CF3	\(x_0+\Xi_1\)	\(1-x_0^2+\Xi_2\)	\(f_2+f_1^2-1\)
CF4	\(x_0+\Sigma_1[y^2]\)	\(1-x_0^2+\Sigma_2[y^2]\)	\(f_2+f_1^2-\sin\!\big(2\pi(f_2-f_1^2+1)\big)-1\)
CF5	\(x_0+\Sigma_1[\psi_c]\)	\(1-x_0+\Sigma_2[\psi_c]\)	\(f_2+f_1-1\)
CF6	\(x_0+\Sigma_1[y^2]\)	\((1-x_0)^2+\Sigma_2[y^2]\)	\(f_2+f_1^2-1\)
CF7	idéntica a CF6	idéntica a CF6	idéntica a CF6

En CF1 el código define \(y=x-\sin\!\big(6\pi x_0+\tfrac{k\pi}{n}\big)\) con \(k=1,\dots,n-1\), es decir, resta un vector de \(n-1\) componentes a uno de \(n\): para \(n>2\) la evaluación falla con un error de dimensiones.

Tri-objetivo (CF8–CF10)

Con \(y_i=x_i-2x_1\sin\!\big(2\pi x_0+\tfrac{(i+1)\pi}{n}\big)\) para \(i\ge2\) y \(y_i=x_i\) para \(i<2\):

\[ \begin{aligned} f_1&=\cos\tfrac{\pi x_0}{2}\,\cos\tfrac{\pi x_1}{2}+\Sigma_1[\psi],\\ f_2&=\cos\tfrac{\pi x_0}{2}\,\sin\tfrac{\pi x_1}{2}+\Sigma_2[\psi],\\ f_3&=\sin\tfrac{\pi x_0}{2}+\Sigma_3[\psi], \end{aligned} \]

con \(\psi(y)=y^2\) en CF8 y CF9, y \(\psi=\psi_c\) en CF10. La restricción, igual en los tres, exige

\[ \frac{f_1^2+f_2^2}{1-f_3^2}-1\;\ge\;0 , \]

y se devuelve con signo cambiado (\(g\le0\)). En el código CF9 es idéntica a CF8; el denominador se anula si \(|f_3|=1\).

A.3 Familia MW

get_mw_problem(nombre, n_var=10) envuelve el problema homónimo de pymoo y devuelve una función que entrega (objs, g) con \(g\le0\). Requiere pip install pymoo (extra benchmarks). Las fórmulas no se reproducen aquí.

Apéndice B: instalación y uso desde Python

Instalación desde GitHub

git clone https://github.com/nanodom/Cancer-GA.git
cd Cancer-GA
pip install -e .                  # núcleo (numpy, matplotlib)
pip install -e ".[benchmarks]"    # opcional: pymoo, para la familia MW

Mono-objetivo

from cancer_ga import CancerGA
from cancer_ga.benchmarks import rastrigin

ga = CancerGA(
    objective_function=rastrigin,
    n_dim=10,
    bounds=(-5.12, 5.12),
    n_cancer=80,
    n_lymph=25,
    seed=42,
)
best_x, best_f, feasible = ga.run(n_generations=100, verbose=True)

Restricciones y variables enteras

from cancer_ga import CancerGA
from cancer_ga.benchmarks import constrained_mixed_example
from cancer_ga.utils import linear_leq_constraints

ga = CancerGA(
    objective_function=constrained_mixed_example,
    n_dim=3,
    bounds=[(-5, 5), (-5, 5), (0, 10)],
    is_integer=[False, False, True],
    constraints=linear_leq_constraints([[-1, -1, 0]], [-1]),  # x0 + x1 >= 1
    seed=1,
)
best_x, best_f, feasible = ga.run(n_generations=100, verbose=False)

Las igualdades se declaran con eq_constraints=linear_eq_constraints(A, b) o con funciones \(h(x)\) propias.

Multiobjetivo con restricciones de un benchmark

Como el algoritmo ignora el segundo elemento que devuelven los benchmarks, las restricciones se declaran explícitamente:

from cancer_ga import CancerGA
from cancer_ga.benchmarks import CTP2

def f(x):
    return CTP2(x)[0]            # vector de objetivos

def c(x):
    return CTP2(x)[1][0]         # g(x) <= 0

ga = CancerGA(
    objective_function=f,
    constraints=[c],
    n_dim=10,
    bounds=(0.0, 1.0),
    n_obj=2,
    seed=1,
)
best_x, best_F, feasible = ga.run(n_generations=200, verbose=False)

Para benchmarks con varias restricciones (p. ej. CTP8) se usa una función por restricción: [lambda x, k=k: CTP8(x)[1][k] for k in range(2)].

Pruebas

pip install -e ".[dev]"
pytest

Referencias

• Bäck, T. y Schwefel, H.-P. (1993). An overview of evolutionary algorithms for parameter optimization. Evolutionary Computation, 1(1), 1–23.
• Deb, K. (2000). An efficient constraint handling method for genetic algorithms. Computer Methods in Applied Mechanics and Engineering, 186(2–4), 311–338.
• Deb, K., Pratap, A., Agarwal, S. y Meyarivan, T. (2002). A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Transactions on Evolutionary Computation, 6(2), 182–197.
• Deb, K., Pratap, A. y Meyarivan, T. (2001). Constrained test problems for multi-objective evolutionary optimization. EMO 2001, LNCS 1993, 284–298.
• Eshelman, L. J. y Schaffer, J. D. (1993). Real-coded genetic algorithms and interval-schemata. Foundations of Genetic Algorithms 2, 187–202.
• Goldberg, D. E. y Richardson, J. (1987). Genetic algorithms with sharing for multimodal function optimization. ICGA, 41–49.
• Ma, Z. y Wang, Y. (2019). Evolutionary constrained multiobjective optimization: Test suite construction and performance comparisons. IEEE Transactions on Evolutionary Computation, 23(6), 972–986.
• Mantegna, R. N. (1994). Fast, accurate algorithm for numerical simulation of Lévy stable stochastic processes. Physical Review E, 49(5), 4677–4683.
• Passino, K. M. (2002). Biomimicry of bacterial foraging for distributed optimization and control. IEEE Control Systems Magazine, 22(3), 52–67.
• Zhang, Q. y Li, H. (2007). MOEA/D: A multiobjective evolutionary algorithm based on decomposition. IEEE Transactions on Evolutionary Computation, 11(6), 712–731.
• Zhang, Q., Zhou, A., Zhao, S., Suganthan, P. N., Liu, W. y Tiwari, S. (2008). Multiobjective optimization test instances for the CEC 2009 special session and competition. Informe técnico CES-487, Universidad de Essex.