1. Introducción

La optimización es una disciplina matemática estudiada desde hace tiempo. Importantes matemáticos como Pierre de Fermat (1601-1665) en su memoria Methodus ad disquirendam maximan et minimam (escrita sobre 1629 pero publicada en 1679 después de su muerte por su hijo Samuel), ya estableció reglas para el cálculo de máximos y mínimos. También Joseph Louis Lagrange (1736-1813) en su publicación Mécanique Analytique (1788-89), fue capaz de hallar fórmulas basadas en el cálculo que permitieron caracterizar los valores óptimos que optimizaban ciertos problemas de optimización. Otros, dieron solución a problemas clásicos que ya se enmarcaban dentro del contexto de la optimización, como Leonhard Euler (1707-1783), con su famoso problema de los siete puentes de Königsberg. Además, cabe resaltar las aportaciones de Isaac Newton (1642-1727) y Johann Carl Friedrich Gauss (1777-1855) los cuales propusieron métodos iterativos que convergiesen a un óptimo.

Posteriormente, a mediados del siglo XX se resolvieron problemas más complejos dentro de esta disciplina. En esta época, se produjo un impulso del estudio de estos problemas; así, se puede consultar [3, pags. 326 - 332] para comprobar la aplicabilidad de los mismos. Además, desde el punto de vista matemático, la resolución de estos problemas resulta especialmente interesante al requerir de distintos conocimientos matemáticos como la topología, el análisis o la matemática computacional. Este hecho, se puede apreciar en [5], ya que en esta obra se puede observar la optimización desde un punto de vista teórico, estudiando casos especiales que ofrecen situaciones idílicas en las que es fácil calcular un óptimo.

Hogaño, debido a la amplia (y creciente) uso de la optimización en ciencia, ingeniería, economía e industria, tal y como se comenta en el prefacio de [12], se hace esencial el estudio de nuevos algoritmos que resuelvan este tipo de problemas. De hecho, hoy en día todas las empresas desean utilizar la optimización en sus procesos de fabricación, desde ahorrar los costes hasta maximizar los beneficios obtenidos. El estudio e implicación de esta disciplina en otras ramas del conocimiento, se puede ver en estudios como los que se incluyen en [7], [6], [9] o [15].

En este artículo se presentará la deducción de un nuevo algoritmo de resolución en los problemas de optimización en varias variables con restricciones. Para ello se realizará una reseña de los resultados más importantes que se utilizarán de forma reiterada. En concreto, se introducirá la definición genérica del problema de estudio, junto con definiciones topológicas que permitan caracterizar el conjunto admisible y las direcciones de descenso, las cuales cobrarán gran importancia en el algoritmo que se planteará. A continuación, se introducirá la deducción del método que presenta este artículo. En esta sección se propondrá la construcción de una sucesión que converja al valor óptimo del problema.

Para ello, se propondrá una idea intuitiva, cuya eficacia quedará entredicho al plantear un contraejemplo, lo que permitirá visualizar una carencia en este sentido. Luego, a partir de la visión del contraejemplo se introducirá una definición que permita solventar la dificultad observada. A partir de este punto, se incluirán una serie conceptos como las coordenadas esféricas en un espacio n-dimensional o el método de Gram-Schmidt para poder exponer un algoritmo aplicable. Finalmente, se estudiarán los aspectos computacionales del algoritmo para poder implementarlo en un lenguaje de programación numérica, concretamente se utilizará MatLab. De este modo, se hará una serie de ejecuciones para comprobar el funcionamiento del método planteado y se comparará con las soluciones exactas a los problemas resueltos y con los resultados obtenidos por otros métodos ya conocidos.

Con este artículo se pretende aportar una nueva visión en la resolución de algunos problemas de optimización. A partir del método propuesto, este queda a objeto de estudio del lector con vistas a una posible mejora en sus aspectos computacionales u otros que se consideren.

Notación 1

Antes de comenzar con la redacción, se hacen las siguientes aclaraciones correspondientes a lanotación que se empleará en lo sucesivo:

1. A lo largo del artículo se emplearán escalares, vectores y matrices de forma frecuente. Para diferenciarlos, se escribirá a los vectores en letra negrita minúscula (por ejemplo v), a las matrices por letra negrita mayúscula (por ejemplo M) y en letra normal minúscula a los escalares (por ejemplo k).

2. El empleo de los vectores será muy reiterado. Dado un vector$v$, si se desea referirse a su componente i-ésima, esta se denotará en letra normal minúscula con el subíndice, es decir,$v_i$.

3. En una sucesión de vectores, se denotará al vector k-ésimo de la sucesión con la siguiente escritura:$v^k$.

4. A los conjuntos se les denotará en letra normal mayúscula, por ejemplo S. El término Int(S) se refiere al interior del conjunto y $\partial \text{S}$ para denotar a la frontera del conjunto.

5. Dada una función$f:{ℝ}^n\to ℝ$, se denotará por $\nabla f\left(x\right)$ al valor del gradiente de la función $f$ evaluado en el vector (o también se denominará por punto)$x$.

6. Si $f\left(x\right)$ y $g\left(x\right)$ son funciones definidas en un entorno de un punto$x_0$, entonces $f=o\left(g\right)$ cuando$x\to x_0$, si y sólo si, para todo $\epsilon >0$ se tiene que $\left|f\right(x\left)\right|<\epsilon \left|g\right(x\left)\right|$ para todo $x$ en un entorno de $x_0$.

7. Se utilizará el concepto de la esfera unitaria $n$-dimensional, que se denotará por ${\text{S}}^n$ y se definirá del siguiente modo: ${\text{S}}^n=\{x\in {ℝ}^n:|\left|x\right||=1\}$

2. Preliminares

Antes de comenzar con la deducción del método que se expondrá en este artículo, se hará una recopilación de preliminares que se utilizarán con cierta frecuencia. Para la realización de estos preliminares se consultaron las referencias: [2], [12], [16], [17], [13], las cuales se pueden aportar una mayor ampliación de los contenidos aquí mostrados.

En este artículo se considerará el problema (1),

Este problema pretende minimizar la función$\text{J}\left(v\right)$, denominada funcional objetivo, en un subconjunto de su dominio, el cual se denominará conjunto admisible o factible. Según las restricciones,${\phi }_i\left(v\right)$, el conjunto admisible quedaría definido del siguiente modo:

$\text{S}:=\{v\in {ℝ}^n:{\phi }_i(v)\le 0,i=1,...,m\}$

De esta manera, el problema inicial (1) se podría reescribir considerando la notación mostrada en (2),

Para resolver el problema (2) se cuentan con múltiples opciones de algoritmos. Algunos de ellos se basan en la creación de una función auxiliar cuyo mínimo en todo su dominio coincida con el mínimo del problema de estudio, otros se caracterizan por la transformación del problema (2) en un problema cuya solución se pueda calcular de un modo más sencillo (como los métodos SQP). Sin embargo, en este artículo se presentará un método que tiene en consideración el concepto de dirección factible, se podría consultar [17, Cap. XI] para estudiar métodos inscritos en esta disciplina. Por este motivo, a continuación se introducirán una serie de definiciones y resultados que serán de enorme importancia para la deducción del algoritmo final.

Definición 1

Dado un funcional $J:{ℝ}^n\to ℝ$ y un conjunto admisible$\text{S}$, se dice que $u$ es solución al problema (1) si se verifica que existe un $\epsilon >0$ tal que,

$J\left(u\right)\le J\left(v\right),\forall v\in S\cap B(u,\epsilon ),$

donde $B(u,\epsilon )$ denota a la bola de centro $u$ y radio $\epsilon$

Para resolver el problema (2) se cuentan con múltiples opciones de algoritmos. Algunos de ellos se basan en la creación de una función auxiliar cuyo mínimo en todo su dominio coincida con el mínimo del problema de estudio, otros se caracterizan por la transformación del problema (2) en un problema cuya solución se pueda calcular de un modo más sencillo (como los métodos SQP). Sin embargo, en este artículo se presentará un método que tiene en consideración el concepto de dirección factible, se podría consultar [17, Cap. XI] para estudiar métodos inscritos en esta disciplina. Por este motivo, a continuación se introducirán una serie de definiciones y resultados que serán de enorme importancia para la deducción del algoritmo final.

Definición 2

Dado un punto $v\in S$ y $d\in {ℝ}^n$ no nulo, si existe un $\delta >0$ de tal forma que$v+td\in S$, para todo$t\in [0,\delta ]$, entonces se dice que $d$ es una dirección factible de $S$ en$v$. Además, el conjunto de dichas direcciones factibles se denotará como:

tal que $v+td\in S;\forall t\in [0,\delta ]\}$.

Observación 1

Para profundizar en la visión geométrica del concepto de dirección factible, se puede consultar [8, pags. 387-393].

Definición 3

Dado un punto$v\in S$, se define el conjunto de restricciones activas de $v$ como sigue:

$I\left(v\right):=\{i\in \{1,...,m\}:{\phi }_i(v)=0\}$

La Definición 2 hace referencia a aquellas direcciones, a partir de las cuales, se puede desplazar una cierta distancia, desde un punto admisible, permaneciendo en el conjunto factible. En la Figura 1 se muestran algunos ejemplos. En ella se vislumbra un conjunto admisible (el conjunto rayado) y desde algunos puntos se dispersan una serie de direcciones que permiten desplazarse a través de ellas permaneciendo en el conjunto.

En este artículo, el método que se planteará utilizará las direcciones factibles como punto de partida. Esta idea se inscribe en aquellos métodos que se conocen con el nombre de métodos de direcciones factibles. La idea consiste en partir de un punto cualquiera del conjunto admisible y, a partir de él,

moverse a través de direcciones factibles que permitan descender el valor del funcional. Para ello, es necesario establecer la Definición 4.

Definición 4

Dado un punto v 2 S, se define el conjunto de restricciones activas de v como sigue:

$J(v+td)<J\left(v\right)$ para algún$t\in ℝ$,

entonces $d$ es una dirección de descenso de $J$ en el punto$v$.

La Definición 4 caracteriza a las direcciones correspondientes a un punto dado que permiten descender el valor del funcional. Este hecho se puede comprobar utilizando la fórmula de Taylor. Pues bien, partiendo de un punto arbitrario $v\in S$ se puede desarrollar la expresión (3),

Así, es fácil comprobar que,

si y sólo si, $d$ es una dirección de descenso de $J$ correspondiente al punto$v$.

Por consiguiente, dado un punto factible cualquiera, la idea que se enmarca dentro de los métodos de direcciones factibles, es la de moverse dentro del conjunto admisible a través de direcciones de descenso que sean factibles. De este modo, una vez que se obtenga un punto para el que no existen direcciones de descenso que sean admisibles, se habrá conseguido un mínimo al problema de partida.

En la Figura 2 se muestra un ejemplo del hecho comentado. En la gráfica se aprecian las curvas de nivel del funcional objetivo $J\left(v\right)={1.8}^{-1.5\sqrt{v_1^2+v_2^2}}·sin\left(v_1\right)cos\left(v_2\right)$ y una elipse que delimita el conjunto admisible. El mínimo del problema se sitúa en el origen, el punto marcado en la gráfica, y de él salen distintas direcciones que delimitan el cono de las direcciones de descenso en el punto. Sin embargo, ninguna de esas direcciones permite permanecer en el conjunto admisible, por lo que no existen direcciones factibles de descenso.

3. Deducción del método

En esta sección se deducirá el método a estudiar en este artículo. Para ello se comenzará aproximando el problema original (1) por el problema (5),

para un valor suficientemente pequeño de $\epsilon >0$.

La ventaja de esta notación radica en que el conjunto admisible $S$ es abierto. Por tanto, dado un punto $v\in S$ existe un conjunto abierto $W\subset S$ tal que$v\in W\subset S$. Entonces, tomando $\delta =\underset{y\in \partial \mathrm{W}}{\mathrm{mín}}\left|\right|y-v\left|\right|$ se tiene que cualquier dirección$d\in {ℝ}^n$, con$\left|\right|d\left|\right|=1$, es admisible, pues $v+td\in S$ para todo$t\in [0,\delta ]$.

Hecho este análisis, el objetivo es conseguir, a partir de un punto factible inicial $u^{\left(0\right)}\in S$ cualquiera, una sucesión${\left\{u^{\left(k\right)}\right\}}_{k\in \mathbb{N}}\to u$, donde $u$ es un mínimo local del problema (5). En este sentido, se optará por construir la sucesión de tal forma que se reduzca el valor del funcional en cada iteración del método, es decir, la sucesión ${\left\{u^{\left(k\right)}\right\}}_{k\in \mathbb{N}}$ se construirá de tal manera que:

$J\left(u^{(k+1)}\right)<J\left(u^{\left(k\right)}\right)<...<J\left(u^{\left(1\right)}\right)<J\left(u^{\left(0\right)}\right)$

Con este objetivo, se propone descender por el conjunto admisible a través de direcciones de descenso. Generalmente no todas las direcciones de descenso son admisibles, pero debido a que $\mathrm{DF}(v,S)={ℝ}^n$ para todo$v\in S$, entonces cualquier dirección de descenso es factible.

Intuitivamente se podría considerar, en un primer momento, como dirección de descenso en la iteración k-ésima el valor de$-\nabla J\left(u^{\left(k\right)}\right)$, pues esta es la dirección de máximo descenso en un punto. Sin embargo, este razonamiento intuitivo (que parecería funcionar) no resulta de utilidad debido a las posibles formas que puede tomar el conjunto factible. A continuación se presenta el siguiente Contraejemplo 1 propio, propuesto en este artículo, para visualizar un caso.

Contraejemplo 1

Se considera el problema de minimización

En la Figura 3 se representa el problema 6. Tomando como punto inicial $u^{\left(0\right)}={(1,5/2)}^{\text{T}}$ (el punto verde) y como dirección de descenso$d^{\left(0\right)}=-\nabla J\left(u^{\left(0\right)}\right)$, lo que se obtendría es un punto contenido en la recta roja de la imagen (dentro del conjunto admisible). Procediendo sucesivamente, sin salir del conjunto factible, se obtendría como punto de acumulación el ${(1/2,5/4)}^{\text{T}}$ (el punto morado de la representación) que dista una distancia de 0.55902 unidades del mínimo del problema (que es el${(0,1)}^{\text{T}}$). Por consiguiente, no se está convergiendo a la solución que se busca.

Entonces, por el motivo visto en el Contraejemplo 1, es necesario ampliar el concepto de dirección de descenso a un marco mayor. Con este objetivo, se plantea la Definición 5.

Definición 5

Se denominará dirección de descensoen el conjunto $S$ respecto de un punto $v$ a una dirección de descenso $d\in {ℝ}^n$ con$\left|\right|d\left|\right|=1$, tal que, siendo $\alpha \in ℝ$ el máximo valor verificando que $v+\alpha d\in S$ y$J(v+\alpha d)<J\left(v\right)$, esta cumple que:

$J(v+\alpha d)<J(v+\alpha \text{'}d\text{'})$, $\forall d\text{'}$ dirección de descenso.

Observación 2

La Definición 5 toma como dirección de descenso aquella a través de la que más se desciende el valor del funcional. Como el conjunto admisible no es el total, puede que dicha dirección no sea la de máximo descenso (en el sentido habitual). Dicho de otro modo, puede que otra dirección de descenso recorra una mayor distancia permaneciendo en el conjunto admisible, y esto haga que el valor del funcional descienda más que usando simplemente el valor negativo del gradiente del funcional en el punto.

Teniendo en cuenta la Definición 5, sin pérdida de generalidad, se van a considerar a partir de ahora las direcciones de descenso normalizadas. Como el objetivo es crear una sucesión ${\left\{u^{\left(k\right)}\right\}}_{k\in \mathbb{N}}\to u$ donde $u$ es solución al problema (5), en cada iteración se va a procurar aproximar la dirección $d^{\left(k\right)}=\frac{u^{\left(k\right)}-u}{\left|\right|u^{\left(k\right)}-u\left|\right|}$ y el valor${\alpha }_k=\left|\right|u^{\left(k\right)}-u\left|\right|$, para tomar$u^{(k+1)}=u^{\left(k\right)}+{\alpha }_k{u}^{\left(k\right)}$.

El mayor inconveniente de este procedimiento erradica en la aproximación de las direcciones$d^{\left(k\right)}$. Para ello se usará el concepto de discretización en el sentido de ángulos. Pues bien, $d^{\left(k\right)}=\frac{u^{\left(k\right)}-u}{\left|\right|u^{\left(k\right)}-u\left|\right|}$ es una dirección de descenso, luego${d^{\left(k\right)}}^{\text{T}}\nabla J\left(u^{\left(k\right)}\right)<0$. Ahora bien, como

${d^{\left(k\right)}}^{\text{T}}\nabla J\left(u^{\left(k\right)}\right)=\left|\right|d^{\left(k\right)}\left|\right|\left|\right|\nabla J\left(u^{\left(k\right)}\right)\left|\right|cos(\measuredangle (d^{\left(k\right)},-\nabla J\left(u^{\left(k\right)}\right)\left)\right)<0$

se deduce que$\measuredangle (d^{\left(k\right)},-\nabla J(u^{\left(k\right)}\left)\right)\in (-\pi /2,\pi /2)$.

Considerando entonces$-\nabla J\left(u^{\left(k\right)}\right)$, se pretende discretizar los valores de los ángulos que forman con$-\nabla J\left(u^{\left(k\right)}\right)$ un ángulo comprendido en el intervalo$(-\pi /2,\pi /2)$. Es decir, se calculan direcciones separadas por un mismo ángulo entre ellas cumpliendo que su producto escalar con la dirección demáximo descenso sea negativo.

Por otra parte, como las direcciones calculadas deben tener norma 1, se pueden interpretar como puntos de la esfera$S^n$. Así, si se tiene en cuenta lo que se pretende hacer con la discretización, lo que se busca es calcular una discretización de puntos del hemisferio de la esfera $S^n$ que tiene como polo el punto${\lceil }^{\left(k\right)}=-\nabla J\left(u^{\left(k\right)}\right)/\left|\right|-\nabla J\left(u^{\left(k\right)}\right)\left|\right|$. Con este objetivo, se recurre a las coordenadas polares, las cuales resultan cómodas desde el punto de vista práctico.

Ejemplo 1

Si se tiene un problema de dimensión 2, se consideraría la esfera $S^2$ (la circunferencia usual).Entonces, se parte de ${\theta }_k\in [0,2\pi ]$ el valor que define el punto ${\lceil }^{\left(k\right)}$ en coordenadas polares, es decir,${\left(cos\right({\theta }_k),sen({\theta }_k\left)\right)}^{\text{T}}={\lceil }^k$. Luego, se divide el intervalo $(0,\pi /2)$ en una serie de puntos igualmente espaciados$0<t_1<...<t_m<\pi /2$, por lo que tomando

se obtienen $2m+1$ puntos que discretizan el hemisferio de la esfera$S^2$, cuyo polo es el punto${\lceil }^{\left(k\right)}$.

En la Figura 4 se muestra un ejemplo de discretización de la esfera $S^2$ entorno al punto ${(\sqrt{2}/2,\sqrt{2}/2)}^{\text{T}}$ (punto rojo). En torno a él, los puntos morados son los que generan la discretización buscada (con 4 puntos en este caso, 8 en total) y los puntos en amarillo serán los que tendrían producto escalar nulo al operarlos con el punto de estudio en este caso.

Para un espacio de dimensión$n>2$, es necesario recurrir a la expresión de las coordenadas esféricas. En el caso general, para una esfera de dimensión$n$, se tiene que las coordenadas de un punto arbitrario

$v$ vienen dadas por,

donde $r$ es el radio de la esfera (en este caso 1, por lo que se puede omitir), y los valores de ${\tau }_1,...,{\tau }_{n-1}$ vienen dados por:

donde ${\tau }_1,...,{\tau }_{n-1}\in [0,\pi ]$ y ${\tau }_{n-1}\in [0,2\pi )$.

Con estas fórmulas es posible crear una discretización de un hemisferio de la esfera$S^2$, con$n>2$. Para ello se procede como en el caso de la esfera$S^3$, salvo que en este caso existen más variables. Para los ángulos${\tau }_1,...,{\tau }_{n-2}$, se divide el intervalo $(0,\pi /2)$ en $m$ partes iguales $(0<t_1<...,t_m<\pi /2)$

y se toman los ángulos ${\tau }_j+t_i$ con $j=1,...,n-2$ e$e=1,...,m$. Por otra parte, para ${\tau }_{n-1}$ se divide el intervalo $[0,2\pi )$ de tal modo que $0<q_1<...<q_{4m}<2\pi$ y se toman los ángulos ${\tau }_{n-1}+q_1$ con$l=1,...,4m$. Entonces, haciendo todas las posibles combinaciones entre los ángulos obtenidos, se obtiene toda la discretización buscada de un hemisferio de la esfera$S^n$. En la Figura 5 se muestra un ejemplo visual de una discretización de $S^3$ con $m=4$ para el hemisferio que tiene por polo el punto$e^{\left(3\right)}={(0,0,1)}^T$.

Sin embargo, el problema que resulta de esta construcción de la esfera, es que esta se entiende como una superficie de revolución. En este caso el eje bajo el que se gira es entorno al vector$e^{\left(1\right)}={(1,0,...,0)}^T$, pero el interés es que sea sobre el vector que determina${\lceil }^k$. Para ello, bajo la discretización obtenida del hemisferio cuyo polo es el vector$e^{\left(1\right)}$, lo que se hará es una isometría $f:{ℝ}^n\to {ℝ}^n$ de tal manera que$f\left(e^{\left(1\right)}\right)={\lceil }^{\left(k\right)}$. Para construir dicha isometría se comenzará tomando el vector ${\lceil }^{\left(k\right)}$ y se completará a una base del espacio vectorial${ℝ}^n$. Luego, para conseguir que esta nueva base sea ortonormal, se puede usar algún procedimiento conocido como el método de ortogonalización de Gram-Schmidt que se muestra en el Algoritmo 1, y cuya demostración y deducción se puede consultaren [4, pag. 273] o [10, pág. 365-366].

Algoritmo 1 Método de ortogonalización de Gram-Schmidt

Si $B=\{u^{\left(1\right)},...,u^{\left(n\right)}\}$ es una base del espacio euclideano$\text{V}$, entonces se puede construir una ortogonal $B\text{'}=\{v^{\left(1\right)},...,v^{\left(n\right)}\}$ tomando

donde ${\lambda }_{ij}=-\frac{{u^{\left(i\right)}}^{\text{T}}{v}^{\left(j\right)}}{{v^{\left(j\right)}}^{\text{T}}{v}^{\left(j\right)}}$

Entonces se completa ${\lceil }^{\left(k\right)}$ a una base y, con el método de Gram-Schimtd, a una base ortonormal. De este modo, sea $A$ la matriz que tiene por columnas los vectores de esta base ortonormal, si el determinante$\left|A\right|=-1$, se intercambian dos columnas (sin incluir la primera). Así se obtiene que$\left|A\right|=1$, por lo que esta matriz determina una aplicación lineal que mantiene la orientación del espacio y, en consecuencia, define una isometría $f:{ℝ}^n\to {ℝ}^n$ tal que$f\left(e^{\left(1\right)}\right)={\lceil }^{\left(k\right)}$. De esta manera, si se emplea la aplicación definida por $A$ a cada punto de la discretización obtenida por el procedimiento que tiene por polo el punto que se buscaba.

Una vez vistos todos los aspectos relativos a la discretización, sólo queda comparar todas las direcciones calculadas por el método de discretización para escoger la mejor. Es decir, se comprobará (numéricamente) cuál, de entre todas las direcciones obtenidas, cumple la Definición 5. Una vez hallada dicha dirección,$d^{\left(k\right)}$, se toma$u^{(k+1)}=u^{\left(k\right)}+{\alpha }_k{d}^{\left(k\right)}$, de tal manera que$J\left(u^{(k+1)}\right)<J\left(u^{\left(k\right)}\right)$, por lo que se crea una sucesión ${\left\{u^{\left(k\right)}\right\}}_{k\in \mathbb{N}}$ de tal forma que

$J\left(u^{(k+1)}\right)<J\left(u^{\left(k\right)}\right)<...<J\left(u^{\left(2\right)}\right)<J\left(u^{\left(1\right)}\right)$

Cuando el valor de ${\alpha }_k$ es muy pequeño, indica que, de entre todas la direcciones, en la que mayor descenso se produce y se avanza, el avance es mínimo y, por tanto, se alcanza una convergencia. Por consiguiente, un criterio de parada que se podría tomar, sería que ${\alpha }_k$ fuese menor que una cierta tolerancia.

Con esto, se concluiría la deducción del método que se plantea. Haciendo una recopilación de lo visto, se puede presentar el siguiente algoritmo:

Algoritmo 2

Se procede con el siguiente procedimiento iterativo:

PASO 0: Se parte de un$u^{\left(0\right)}\in S$, un valor de tolerancia $tol>0$ pequeño y $k=0$.

PASO 1: Se calcula ${\lceil }^{\left(k\right)}$ y, a partir de él, se construye una discretización.

PASO 2: A partir de la discretización obtenida, se comparan todas las direcciones para encontrar aquella que cumple la Definición 5, de lo que se obtiene como solución la dirección $d^{\left(k\right)}$ y el paso ${\alpha }_k$.

PASO 3: Se toma$u^{(k+1)}=u^{\left(k\right)}+{\alpha }_k{d}^{\left(k\right)}$. Si ${\alpha }_k<tol$ se finaliza y, en caso contrario, se vuelve al PASO 1.

4. Aspectos computacionales

Para poder aplicar en la práctica el Algoritmo 2, este se debe poder programar en algún lenguaje de programación para su utilización. Por ello, en esta sección se presentará una manera de programar el método planteado en el lenguaje de MatLab.

En Código 1 se muestra un ejemplo de programa principal. En este programa se muestra una forma de programar el Algoritmo 2. Sin embargo, en él se utilizan diferentes funciones auxiliares. En concreto se utiliza la función “discretizacion” para crear la discretización de las direcciones y la función “conjunto” para determinar si el punto obtenido al desplazarse en una dirección un determinado valor, a partir de un punto dado, pertenece al conjunto admisible.

En Código 2 se expone la función empleada para calcular la discretización del hemisferio de la esfera $S^n$ cuyo polo es un vector determinado. Los datos de entrada son$d$, el vector que determina el polo en torno al que se desea discretizar la esfera, y el valor de $k$ que se corresponde con el número de puntos de la discretización. De este modo, para un $k$ mayor se obtienen más puntos de la discretización, mientras que si se reduce su valor, también lo hace el número de direcciones a comprobar. Por otra parte, se utiliza una función con el nombre “gramscb” que lo que hace es aplicar el método de ortogonalización de Gram-Schmidt explicado en el Algoritmo 1 (dado que el método de Gram-Schmidt es de fácil aplicación, se deja su programación al lector).

En Código 3 se expone el programa utilizado para comprobar la factibilidad de un nuevo elemento obtenido a partir de un punto arbitrario, una dirección y un valor “alpha” que se corresponde con el desplazamiento desde el punto arbitrario en la dirección dada. Este programa verifica si el siguiente iterante del Algoritmo 2 es admisible para el problema (2) tomando como valor de $\epsilon$ el propio épsilon de la máquina.

De esta forma se concluiría con los aspectos computacionales referentes al Algoritmo 2. Es importante resaltar que la programación aquí mostrada es un ejemplo de cómo se podría llevar a la práctica, pero en realidad podría haber otras versiones. Un ejercicio interesante resultaría de reflexionar la forma en la que se podría abaratar el coste computacional de esta programación. Estos aspectos se dejan a reflexión de los lectores y a juicio del público para mejorar el método y la programación planteados.

5. Pruebas y ejecuciones

Con el objetivo de testear el buen funcionamiento del Algoritmo 2, se ha probado su ejecución en una serie de problemas de sencillo cálculo empleando la programación expuesta en la Sección 4. Para ello se seleccionaron una serie de problemas, extraídos de [11], con los que corregir el método creado y comprobar su funcionamiento, además de incluir el Contra-Ejemplo 1. Así mismo, esos mismos problemas se resolverán con otros algoritmos ya conocidos para comparar los resultados obtenidos. En este sentido, en [1] se pueden apreciar una recopilación de métodos ya conocidos que el lector puede consultar como ampliación de esta sección. Sin embargo, en este artículo se mostrarán los resultados obtenidos con la función fmincon de MatLab y el método del Lagrangiano aumentado presentado en el artículo [14].

- Solución al problema del Contraejemplo 1

Se considera el problema (6), que se tomó en el Contraejemplo 1, para ver el funcionamiento del algoritmo en este caso:

El iterante inicial tomado es $u^{\left(0\right)}={\left(\begin{array}{cc}1& 5/2\end{array}\right)}^T$ de tal forma que$J\left(u^{\left(0\right)}\right)=29/4$. La solución exacta al problema es$u={\left(\begin{array}{cc}1& 5/2\end{array}\right)}^T$, lo que implica que el valor del funcional objetivo en el mínimo es$J\left(u\right)=1$.

- Problema de 2 dimensiones:

Se considera el problema (8), que se corresponde con el ejemplo número 22 en [11],

El iterante inicial tomado es $u^{\left(0\right)}={\left(\begin{array}{cc}0& 0\end{array}\right)}^T$ de tal forma que$J\left(u^{\left(0\right)}\right)=0$. La solución exacta al problema es $u={\left(\begin{array}{cc}1& 1\end{array}\right)}^T$ lo que implica que el valor del funcional en el mínimo es$J\left(u\right)=1$.

- Problema de 3 dimensiones:

Se considera el problema (9), que se corresponde con el ejemplo número 35 en [11],

El iterante inicial tomado es$u^{\left(0\right)}={\left(\begin{array}{ccc}1/2& 1/2& 1/2\end{array}\right)}^T$, de tal forma que$J\left(u^{\left(0\right)}\right)=2.25$. La solución exacta al problema es$u={\left(\begin{array}{ccc}4/3& 7/9& 4/9\end{array}\right)}^T$, lo que implica que el valor del funcional en el mínimo es $J\left(u\right)=1/9$.

- Problema de 4 dimensiones:

Se considera el problema (10), que se corresponde con el ejemplo número 43 en [11],

El iterante inicial tomado es$u^{\left(0\right)}={\left(\begin{array}{cccc}0& 0& 0& 0\end{array}\right)}^T$, de tal forma que$J\left(u^{\left(0\right)}\right)=0$. La solución exacta al problema es$u={\left(\begin{array}{cccc}0& 1& 2& -1\end{array}\right)}^T$, lo que implica que el valor del funcional en el mínimo es$J\left(u\right)=-44$.

A la vista de los resultados obtenidos se puede comprobar como el testeo del Algoritmo 2 ha respondido de forma positiva. Se podrían hacer más testeos, e incluso más complejos, pero se deja como ejercicio al lector. Lo que se ha querido mostrar con estas ejecuciones es el correcto funcionamiento del método planteado que se aproxima al valor del mínimo del problema original.

Por otra parte, un aspecto importante son los tiempos de ejecución, que no se incluyen en los resultados debido a que estos varían en función de cada computadora. Sin embargo, hablando de los tiempos obtenidos mediante el Algoritmo 2, es conveniente decir que en el problema (8) se obtiene el resultado en un tiempo razonable, algo que incluso sucede con el problema (9). Sin embargo, en la resolución del problema (10) el tiempo de ejecución se dispara a más de un minuto, algo extraño en un problema de sólo cuatro dimensiones. Este hecho se debe al número de direcciones que el algoritmo debe de evaluar en cada iteración. Con dimensiones pequeñas el número de direcciones de descenso, creadas con la discretización, no es demasiado elevado, pero este número crece de forma exponencial según se aumenta la dimensión del problema.

Esta carencia es importante tenerla en cuenta. De hecho, para que el Algoritmo 2 aporte resultados concluyentes, es importante construir una discretización amplia para poder comparar el mayor número de direcciones de descenso posibles. Por este motivo, el algoritmo se ve limitado en el tiempo de ejecución y, al mismo tiempo, en la memoria requerida en cada iteración para almacenar las direcciones que se deben de comprobar.

Por consiguiente, una futura línea de investigación para el método planteado en este artículo es la forma en la que se podría reducir la memoria y el tiempo de ejecución empleado por el algoritmo. En este sentido se podría hacer un estudio acerca de la forma en la que se podría reducir el número de direcciones de descenso en la discretización sin que esto implique menor precisión del método.

6. Conclusiones y futuras líneas de investigación

En base al trabajo desarrollado en el presente artículo, se obtienen las siguientes conclusiones:

1. 1. El método expuesto se enmarca dentro del conjunto de métodos de direcciones factibles, permitiendo resolver una serie de problemas de optimización con restricciones cuyos conjuntos admisibles sean cerrados y conexos.
2. 2. El interés prioritario del algoritmo planteado se basa en el aporte teórico de su deducción. A lo largo del artículo se muestran diferentes interpretaciones geométricas y se relacionan distintas disciplinas matemáticas para obtener el método final.
3. 3. Una de las ventajas del algoritmo es que permite crear una sucesión de puntos factibles que descienden el valor del funcional objetivo en cada iteración. Así pues, siempre se está descendiendo el valor del funcional en cada iteración.
4. 4. La principal desventaja que presenta el método es su coste computacional, lo que implica altos tiempos de ejecución en la práctica según aumenta la dimensión de los problemas a resolver. Por otra parte, el coste de memoria de la programación expuesta es también una desventaja del método según se aumenta el de la dimensión de los problemas a resolver.

Una vez expuestas las principales conclusiones, se muestran las principales líneas de investigación para mejorar el método expuesto en el artículo:

1. 1. Resulta primordial poder reducir la memoria requerida por la programación del método. Sería muy interesante descartar, en cada iteración, una serie de direcciones en lugar de tener que considerarlas.
2. 2. La reducción de direcciones de descenso debe de hacerse de tal forma que la discretización no se vea reducida, para no menguar el nivel de precisión del método.
3. 3. En lugar de reducir el valor de α Importar imagen a la mitad cómo se hace en el Código 1, sería conveniente desarrollar un método más efectivo de calcular el valor de α. Para ello se podría investigar acerca de cómo aplicar la regla del Armijo, de Goldestein o de Wolfe-Powell a este algoritmo, las cuales se pueden estudiar en [17, Sec. 2.5], que podrían aportar un mejor método de reducir el valor inicial de α.

El método expuesto se enmarca dentro del conjunto de métodos de direcciones factibles, permitiendo resolver una serie de problemas de optimización con restricciones cuyos conjuntos admisibles sean cerrados y conexos.

El interés prioritario del algoritmo planteado se basa en el aporte teórico de su deducción. A lo largo del artículo se muestran diferentes interpretaciones geométricas y se relacionan distintas disciplinas matemáticas para obtener el método final.

Una de las ventajas del algoritmo es que permite crear una sucesión de puntos factibles que descienden el valor del funcional objetivo en cada iteración. Así pues, siempre se está descendiendo el valor del funcional en cada iteración.

La principal desventaja que presenta el método es su coste computacional, lo que implica altos tiempos de ejecución en la práctica según aumenta la dimensión de los problemas a resolver. Por otra parte, el coste de memoria de la programación expuesta es también una desventaja del método según se aumenta el de la dimensión de los problemas a resolver.

Bibliografía

[1] Ancau, M. Practical Optimization with MatLab. Cambridge Scholars Publishing. 2019.

[2] Boglárka G.-Tóth and Eligius M.T. Hendrix. Introduction to Nonlinear and Global Optimization. Springer. 2010.

[3] Crilly, A.J, and Pérez, E.H. 50 cosas que hay que saber sobre matemáticas, Ariel, 2014.

[4] De Burgos J. álgebra lineal y geometría cartesiana. Ed. MacGraw-Hill, Madrid, 1999.

[5] Eiselt, H.A. and Sandblom, C..L. Nonlinear Optimization. Methods and Applications. Springer. 2019.

[6] Eremeev, A. Khachay, M., Kochetov, Y. and Pardalos, P. Optimization Problems and Their Applications. Springer. 2018.

[7] Escudero Andino, F. F.Diseño de un modelo matemático para optimizar las rutas de recorrido del proceso de recolección de desechos sólidos para el Cantón Valencia. Master's thesis, Universidad Técnica de Ambato. Facultad de Ingeniería en Sistemas, Electrónica e Industrial. Maestría en Matemática Aplicada. 2021

[8] Gallier, J. and Quaintance, J. Fundamentals of Optimization TheoryWith Applications to Machine Learning. University of Pennsylvania Philadelphia. 2019.

[9] González Castaño, F. A. Análisis de ciclo de vida y optimización matemática como herramientas para la producción sustentable de energía y de productos químicos. 2019

[10] Hernández, E. álgebra y geometría. Ed. Addison Wesley, Madrid, 1994.

[11] Hock W. and Schittkowski K. Test Examples for Nonlinear Programming Codes, Lecture Notes in Economical and Mathematical Systems 187, Springer-Verlag, Berlin, 1981.

[12] Nocedal, Jorge and Wright, Stephen J. Numerical Optimization. Springer. 1999.

[13] Peressini, A.L; Sullivan, F.E and Uhl, J.J. The mathematics of nonlinear programming. Springer-Verlag. 1988

[14] Vázquez, M. Resolución de los problemas de optimización con restricciones mediante los métodos de penalización y del Lagrangiano aumentado. Revista digital Matemática, Educación e Internet, 21(2). 2021. Recuperado de: https://tecdigital.tec.ac.cr/revistamatematica/

[15] Silva Melgarejo, M. A., and Zevallos Diaz, A. J. E.Programación lineal para optimizar separación de grasas del proceso de fabricación de harina de pescado. Corporación Hayduk SA Coishco. 2019.

[16] Yang, X-S. Optimization Techniques and Applications with Examples. John Wiley and Sons. 2018.

[17] Sun W. and Y-X Yuan. Optimization Theory and Methods. Nonlinear Programming, Springer, 2006.