Optimización con restricciones

Traballo Fin de Grao Optimización con Restricciones Manuel Vázquez Mourazos 2018/2019 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao Optimización con Restricciones Manuel Vázquez Mourazos Julio, 2019 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA Trabajo propuesto Área de Coñecemento: Matemática aplicada Título: Optimización con restricciones Title: Constrained optimization Director: Jerónimo Rodríguez García Breve descrición do contido Este trabajo consistirá en el estudio de problemas de optimización con restricciones. Se estudiarán condiciones necesarias y condiciones suficientes para la existencia y unicidad de solución del problema dependiendo de las hipótesis sobre el funcional objetivo y las restricciones que definen el conjunto admisible. Se estudiarán también algunos métodos numéricos básicos para la resolución de este tipo de problemas. Recomendacións Haber cursado la asignatura “Métodos Numéricos en Optimización y Ecuaciones Diferenciales” iii iv v “Puesto que el Universo es perfecto y fue creado por el Creador más sabio, nada ocurre en él sin que esté presente alguna ley de máximo o mínimo”. L. Euler (s. XVIII) Índice general Índice de figuras ix Índice de tablas xi Resumen xiii xiii Notación xv Introducción xvii 1. Optimización sin restricciones 1 1.1. Conceptosgenerales ............................... 1 1.2. Condiciones de optimalidad con derivadas . . . . . . . . . . . . . . . . . . . 2 1.3. Existencia y unicidad de solución . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3.1. Funcionales coercitivos . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3.2. Funcionales convexos . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4. Funcionales cuadráticos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.5. Métodos para la obtención del mínimo . . . . . . . . . . . . . . . . . . . . . 8 1.5.1. Cálculodelpaso ............................. 9 1.5.2. Método del gradiente . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.5.3. Método del gradiente conjugado . . . . . . . . . . . . . . . . . . . . . 11 1.5.4. MétododeNewton............................ 12 1.5.5. Métodos Cuasi-Newton . . . . . . . . . . . . . . . . . . . . . . . . . 13 2. Optimización con restricciones, teoría 15 2.1. Presentación del problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2. Caracterización de óptimos . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3. Multiplicadores de Lagrange . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 vii xiv In addition, these knowledge will be applied to quadratic optimization. This is a special case that usually causes special interest because of its simplicity in many cases, and what is no less important, its use to solve more complex problems. Finally, the optimization problem will be addressed from a practical point of view. In this case, different methods that show different ways of approaching the problem will be studied. Thus, the algorithms that can be implemented in any programming language (such as MATLAB or others) will be explicitly provided to finally obtain the desired solution. Notación Con el motivo de aclarar las posibles dudas que puedan surgir en la lectura del presente trabajo. A continuación se harán algunas aclaraciones con el fin de explicar la notación que se usará a lo largo del mismo: 1. Como se usará con gran frecuencia operaciones que involucran escalares, vectores y matrices, los vectores se escribirán con letra minúscula en negrita, mientras que para denotar a una matriz, se utilizarán letras mayúsculas en negrita. Por ejemplo, vserá un vector (y 0sería el vector idénticamente nulo), Muna matriz, mientras que cualquier escalar se representará con letra normal. 2. Cuando se comparen dos vectores, se entiende que se comparan componente a componente. Es decir, v≥0significaría que todas las componentes del vector vson positivas. 3. Como también se usarán sucesiones de vectores, para referirse al vector k-ésimo de una sucesión, se utilizará el superíndice (k). Es decir, v(k), denota el elemento k-ésimo de la sucesión v(k)k∈N. Por otra parte, la componente i-ésima de un vector v, se denotará por el escalar vi(así mismo, v(k) idenotará el escalar correspondiente a la componente i-ésima del vector k-ésimo). Además, el escalar correspondiente a una posición (i, j)de una matriz cualquiera M, se denotará por mij. 4. Los vectores que se considerarán en el trabajo serán columna, es decir, dado v∈Rn se tiene que v∈ Mn×1. 5. Dada una matriz M∈ Mn×t, se denotará por miel vector correspondiente a la fila i-ésima de la matriz M. Éste será el único caso excepcional donde se considerarán vectores fila, dado que al corresponderse con la fila i-ésima, se entenderá que mies un vector fila. 6. En los capítulos 2 y 3, se considera un problema genérico de optimización con restricciones al que se denotará por (Q). En el caso de que sólo se consideren restricciones xv de igualdad, el problema genérico se denota por (Qig), mientras que, si sólo se consideran las restricciones de desigualdad, se denotarán por (Qdes). Con esta misma filosofía, siempre que se escriba el subíndice ig (resp. des) sólo se estarán a considerar restricciones de igualdad (resp. desigualdad). 7. En las restricciones que se utilizarán en los capítulos 2 y 3, se denotará por ϕ(v)a las restricciones de desigualdad, y por φ(v)a las restricciones de igualdad. Así, también se denotará por Ψ(v) = (ϕ1(v), . . . , ϕm(v))Tal vector que contiene el valor de las restricciones de desigualdad (de forma análoga, Φ(v) = (φ1(v), . . . , φp(v))T). 8. Teniendo en cuenta las funciones Ψ(v)yΦ(v)de las que se ha hablado. Se denotará por Ψ0(v)a la matriz que tiene por filas los gradientes de cada una de las restricciones de desigualdad (respectivamente Φ0(v)con las restricciones de igualdad). De todas formas, a lo largo del trabajo se recordará la notación. Pero se ha considerado pertinente establecer estas bases para tener un punto de referencia en el caso de que surja alguna duda al respecto. Introducción Definida por la RAE como “la acción o efecto de buscar la mejor manera de realizar una actividad”, la optimización es una rama de la matemática aplicada que pretende hallar el elemento del domino de una función (a la que comúnmente se denomina como funcional objetivo) que lo minimice o maximice. Éste, es un problema que ha sido tratado desde hace siglos, pero no ha sido hasta el siglo pasado, cuando más se ha estudiado y profundizado en su resolución. Históricamente, su enfoque ha sido diverso a lo largo de la historia. En tiempos de antaño, 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 de este mismo contexto, como Leonhard Euler (1707-1783), con su famoso problema de los siete puentes de Königsberg. Desde otro punto de vista, cabe resaltar las aportaciones de Isaac Newton (1642-1727) yJohann Carl Friedrich Gauss (1777-1855), los cuales propusieron métodos iterativos que convergiesen a un óptimo. Sin embargo, el máximo desarrollo de la optimización comenzó a producirse en el siglo pasado. Principalmente, a causa de las guerras acaecidas, los problemas de optimización empezaron a tener una gran relevancia para minimizar costes (en el transporte de soldados o víveres) o maximizando beneficios (procurando abordar el mayor objetivo posible que fuese menester conquistar). Además, fue entonces cuando se pasó a conocer la disciplina de la optimización con el nombre de programación, el cuál, se comenzó a utilizar debido al uso de este término por el ejército estadounidense para referirse a los problemas de logística que se deseaban abordar mediante el uso de esta disciplina. Fue precisamente en ese momento, cuando se produjo el mayor desarrollamiento de la matemática aplicada, y en particular, de la optimización. Fundada en 1941, la Oficina xvii xviii INTRODUCCIÓN de Investigación Científica y Desarrollo (IOCD), fue el primer órgano organizador de los matemáticos contratados por las Fuerzas Armadas de los E.E.U.U para el desarrollo de las matemáticas requeridas por la II Guerra Mundial. Sin embargo, fue en 1943 cuando se creó el Comité de Matemática Aplicada (CMA), que siendo un órgano dentro de la IOCD, sirvió para aprovechar el financiamiento del gobierno (a causa del momento) para el desarrollo de la matemática aplicada. En este contexto, en un primer lugar se desarrolló la programación lineal, que comprendía la resolución de problemas de optimización de una función lineal en un conjunto definido por restricciones afines. En este marco, el algoritmo más famoso es el Símplex, que fue desarrollado en 1947 por George Bernard Dantzig (1914-2005). El cual, ya llevaba estudiado el tema al menos desde 1941, fecha en la que fue contratado por las Fuerzas Armadas para llevar a cabo enormes planteamientos logísticos. Además, también cabe destacar la teoría de la dualidad que fue publicada en el mismo año (1947) por John von Neumann (1903-1957). Figura 1: Importantes exponentes de la programación lineal en el siglo XX De forma paralela, y también posterior a la II Guerra Mundial, a parte de continuar estudiando los problemas de programación lineal, también se amplió el estudio a la programación no lineal. La cual, procuraba dar solución a aquellos problemas donde no había linealidad en las funciones a optimizar ni en las restricciones que definían los conjuntos donde se optimizaban dichas funciones. Para ello, fue necesario generalizar conceptos ya conocidos a un caso más genérico. En este marco, una de las aportaciones más importantes fueron las conocidas con- INTRODUCCIÓN xix diciones de Karush-Kuhn-Tucker (conocidas comúnmente como condiciones KKT), que generalizan las condiciones que había dado Lagrange en su día. Éstas, fueron publicadas por primera vez en el verano de 1950 de forma conjunta por Harold William Kuhn (1925- 2014) y Albert William Tucker (1905-1995). Los cuales, utilizaron también en dicho año por primera vez el término de programación no lineal para referirse a esta disciplina. Sin embargo, se descubrió con posterioridad que William Karush (1917-1997) había demostrado el mismo resultado en su trabajo de fin de máster del año 1939. Por consiguiente, hoy en día se atribuye dicho descubrimiento a los tres en conjunto. Además, íntimamente ligadas con las condiciones KKT, el matemático Fritz John (1910- 1994) publicó en un artículo en 1948 un resultado similar al de las condiciones KKT. El cual, las generaliza a un marco más amplio y se conocen en honor a él, como las condiciones de Fritz-John. Figura 2: Importantes exponentes de la programación no lineal en el siglo XX Desde que fueron publicadas, las condiciones KKT adquirieron gran fama y prestigio. Se podría decir, que fueron el primer paso hacia esta nueva rama del saber. Gracias a su reconocimiento, dentro del mundo matemático, se creó oficialmente en 1971 la primera revista que trató a la optimización como tema fundamental, Mathematical Programming. La cuál, precedió a una segunda bautizada con el nombre de Mathematical Programming Society. Finalmente, con el enorme desarrollo de la informática y de las máquinas de cálculo en la última mitad del siglo pasado, se diseñaron muy diversos algoritmos que permitieron dar solución a los problemas de programación. La variedad en dichos algoritmos es muy diversa en función del enfoque que se le da al problema para hallar la solución. En una mayoría, muchos de ellos buscan calcular puntos cumpliendo las condiciones KKT y hacer un estudio xx INTRODUCCIÓN posterior para comprobar si son realmente solución. Por otro lado, otros intentan buscar el mínimo directamente teniendo en cuenta las características del conjunto donde se quiere minimizar, mientras que, también los hay que recurren a la resolución de subproblemas más sencillos que permitan crear un proceso iterativo que llegue a la solución. Hogaño, la optimización esta presente en multitud de facetas de nuestra vida cotidiana. Las compañías aéreas, diseñan los vuelos con la intención de minimizar los costes y maximizar su beneficio. En la medicina, los tratamientos contra el cáncer, pasan por un proceso de optimización para minimizar su impacto nocivo en la salud de los pacientes. Los inversores procuran minimizar los riesgos en sus decisiones a la vez que garantizar una rentabilidad satisfactoria. Pero más allá, incluso la propia naturaleza está optimizando continuamente procurando un estado de mínima energía. Es el caso de los rayos de luz que siguen las trayectorias que minimizan la duración de su recorrido. Todos estos problemas actuales, se visionan en una formulación matemática como problemas de optimización, de hecho, en [5, Cap. 12] se modelizan multitud de problemas en una formulación adaptada a la programación matemática. Pero calcular su solución, puede llegar a ser bastante complejo en ocasiones. Este trabajo, pretende dar respuesta a estos problemas tratando el caso en el que las funciones a minimizar y las restricciones, sean diferenciables. Para empezar, el primer capítulo se dedica a la optimización sin restricciones. Éste es el caso más sencillo dentro de esta disciplina, la cuál desea hallar el valor óptimo del funcional objetivo en todo su dominio. Para ello, se comienza presentando el problema y caracterizando sus óptimos. Como se ha dicho, en este trabajo se trata el caso donde todas las funciones que intervienen en el problema son diferenciables, y éste, será un hecho crucial para poder caracterizar los óptimos. Posteriormente, se aumentarán las hipótesis sobre los problemas a tratar para poder demostrar cuando éstos tienen solución y ésta es única. Para ello, se introducirá la definición de coercitividad y convexidad. Una vez explicados estos conceptos, se introducirán los funcionales cuadráticos, los cuales son interesantes desde el punto de vista teórico, debido a que permiten aplicarles con sencillez las condiciones de optimalidad que se hayan visto con anterioridad. Y para finalizar el capítulo, se presentarán diversos métodos, a través de los cuales, se podrá calcular el óptimo de forma exacta o numérica (aproximada). En este contexto, todos estos métodos están basados en un mismo algoritmo, y lo que se presentará son diferentes formas de abordar cada paso de dicho algoritmo. De esta forma, combinando las diferentes maneras de abordar los pasos, se tendrá una enorme familia de métodos. El interés que subyace detrás de este capítulo inicial, son varios. Por un lado, el de INTRODUCCIÓN xxi comenzar a manejar los conceptos básicos de la optimización, así como, resultados básicos que se utilizarán a lo largo de todo el trabajo. Y por otra parte, el de conocer diferentes formas de resolver los problemas de optimización sin restricciones, ya que será imprescindible para resolver los problemas de optimización con restricciones. Posteriormente, el segundo capítulo ya se introducirá de pleno en el mundo de la optimización con restricciones. En este caso, el objetivo sigue siendo el mismo, hallar el óptimo de un funcional objetivo. Sin embargo, se añade la peculiaridad de que no todos los valores de su dominio son admisibles. Esto quiere decir que sólo se desea optimizar dicho funcional en un subconjunto de su dominio. Este subconjunto, quedará determinado por una serie de funciones que se denominan restricciones, y nuevamente, sólo se tratará en este trabajo el caso en el que dichas restricciones sean diferenciables. El hecho de la diferenciabilidad, vuelve a ser muy interesante en este capítulo, ya que permite dar una variedad de condiciones necesarias y/o suficientes que permiten caracterizar los óptimos del problema. Y precisamente, será ese el primer objetivo del capítulo. Luego, se introducirán los multiplicadores de Lagrange, que volverán a dar condiciones necesarias para óptimos en un caso especial de restricciones. De esta forma, se ampliarán dichas condiciones a un caso genérico incluyendo más tipos de restricciones al introducir los multiplicadores de Karush-Kuhn-Tucker y Fritz-John. Todos esos multiplicadores, permitirán aportar importantes resultados para seguir caracterizando los óptimos de los problemas que se tratan. Estos resultados, serán de vital importancia con posterioridad debido a que muchos métodos que se verán con posterioridad, se basarán en conseguir puntos que cumplan dichas condiciones (los cuales son candidatos a óptimo) y haciendo un estudio posterior, se determinaría si realmente son óptimos o no. En este sentido, también se aportarán condiciones suficientes bajo las cuales, las condiciones vistas aporten exactamente el óptimo del problema. Finalmente y para terminar este segundo capítulo, se introducirán algunas nociones básicas sobre la teoría de la dualidad. Su interés erradica en la creación de un nuevo problema (al que se denominará como problema dual) con el objetivo de que éste sea más sencillo que el problema inicial (al que se conocerá como problema primal). De esta forma, se podrá establecer una relación entre las soluciones de ambos problemas para que en ocasiones fuese conveniente resolver el problema dual en vez del primal. Todo esto, se aplicará también a la programación cuadrática. Este tipo especial de optimización, se basa en un funcional cuadrático (se definirá en los capítulos uno y dos) y un conjunto definido por restricciones afines. La particularidad de este tipo de problemas, permitirá aplicar los conceptos introducidos en el capítulo dos de forma sencilla y simple. Sin embargo, es difícil que en la realidad, los problemas de optimización que se plantean, xxii INTRODUCCIÓN se enmarquen dentro de este marco de problemas. Pero el interés mayoritario de estos, es debido a su utilización en algunos métodos de resolución de problemas genéricos, y es por ello, que se le hace un hueco en el trabajo. Finalmente, todo concluye en el tercer y último capítulo dedicado en exclusividad a los métodos de resolución. Aquí será donde se empleen los conocimientos adquiridos en los capítulos anteriores para poder deducir todos los algoritmos que se planteen. La primera sección se dedicará a los métodos de penalización. Se comenzará por la penalización exterior seguida de la penalización interior. Éstos, son los más intuitivos y abordan directamente el problemas penalizando aquellos puntos que no se encuentren dentro del conjunto admisible donde se desea optimizar el funcional objetivo. Así pues, estos métodos aportan como solución un óptimo calculado de forma numérica resolviendo una secuencia de problemas sin restricciones. Sin embargo, en la práctica a veces resulta problemática la resolución de estos tipos de problemas, por lo que se introduce un nuevo método. Éste es el conocido como Lagrangiano aumentado, y se basa en las condiciones de optimalidad que aportan los multiplicadores de Lagrange y KKT. En este caso, el método es muy interesante debido a su buen funcionamiento en la práctica, pero en vez de converger a un óptimo, lo hace a un punto cumpliendo las condiciones KKT. Es por ello, que es necesario recurrir a los conceptos teóricos explicados en el capítulo dos para estudiar si el punto obtenido es realmente un óptimo. A continuación, se dará otra visión de resolución hablando de métodos basados en direcciones factibles. En este caso, este tipo de métodos, buscan el óptimo mediante un proceso iterativo que produzca una sucesión de puntos del conjunto donde se desea optimizar el funcional. Para ello, un primer método será el desarrollado por Zoutendijk, que debe recurrir a la programación lineal. Nótese, que no es extraño recurrir a la programación lineal, debido a la existencia de métodos que funcionan muy bien en la práctica como el método Símplex. En este trabajo, no se trata la programación lineal (dado que el objetivo es tratar la programación no lineal diferenciable), pero se aportarán referencias que se pueden consultar en las que se explica con detalle el método Símplex, así como, su deducción. Además, también se aportará otro método basado en direcciones factibles utilizado en la resolución de la programación cuadrática. Finalmente, se dedicará la última sección a los métodos de tipo SQP (sequential quadratic programming) que como bien indica su nombre, se basan en la programación cuadrática. Y es que su metodología se basa en la resolución sucesiva de problemas de tipo cuadrático. Por este hecho, se utiliza los conocimientos que se hayan introducido a lo largo del INTRODUCCIÓN xxiii trabajo sobre los problemas cuadráticos. Nuevamente, su convergencia aportará puntos cumpliendo las condiciones que derivan de los multiplicadores de Lagrange y KKT. Es por ello, que se debería hacer una reflexión posterior al igual que sucedía ya con otro tipo de métodos que convergían a estos tipos de puntos. Por otra parte, también se han añadido dos Anexos que contienen pseudocódigos relativos a los algoritmos y métodos introducidos a lo largo del trabajo. Su interés es completar la información o aclarar dudas que podrían surgir en la lectura de los algoritmos. Sin embargo, si se desea consultar los códigos completos de los métodos introducidos, también se puede hacer en el archivo .zip que se adjunta con la entraga del trabajo. En él, se tienen los archivos .m que se hicieron con la programación en Matlab de todos los métodos que se introdujeron en el trabajo, y también, algunas ejecuciones de prueba que muestra su funcionalidad. De esta forma, se pondría punto y final a este trabajo. Lo cierto, es que hay muchos más métodos de los que se podría hablar y basados en muchos más enfoques. Pues bien, el tema es enorme y la diversidad inmensa, por lo que se ha querido centrar el trabajo en aquellos algoritmos y técnicas que se suelen usar con mayor frecuencia dando una visión diversa y útil. Por mi parte, espero que resulte satisfactoria y agradable la lectura del presente trabajo. Con él, deseo transmitir el primor de esta joven rama de las matemáticas que en su corta vida, tanto ha dado de si y para la que se vislumbra un futuro prometedor lleno de investigación y nuevos descubrimientos. Es así, que sólo anhelo transmitir con su lectura, aunque sólo sea una pizca, la enorme pasión y curiosidad que la optimización ha suscitado en mi interior durante la realización de este trabajo de fin de grado. 6CAPÍTULO 1. OPTIMIZACIÓN SIN RESTRICCIONES Demostración. Se puede consultar en [21, pág. 56]. Teorema 1.13. Sea Ω⊆Rnun abierto convexo y J: Ω ⊆Rn→Run funcional diferenciable. Entonces se verifica: (1) Jes convexo, si y sólo si, J(w)≥J(v) + (∇J(v))T(w−v)∀v,w∈Ω. (2) Jes estrictamente convexo, si y sólo si, J(w)> J(v) + (∇J(v))T(w−v)∀v,w∈Ω. Demostración. Se puede ver en [15, págs. 42-43]. Teorema 1.14. Sea Ωun abierto convexo, con J: Ω ⊆Rn→Rclase 2 en Ω, se tiene, (1) Jes convexo, si y sólo si, HJ(v)es semidefinida positiva ∀v∈Ω. (2) Si HJ(v)es definida positiva ∀v∈Ω, entonces, Jes estrictamente convexo. Demostración. Se puede consultar en [15, págs. 43-44]. Los resultados introducidos, son los que se utilizarán a lo largo del trabajo. Pero para mayor información acerca de los conjuntos y funcionales convexos, se puede consultar [4]. Teorema 1.15. Sea Ω⊆Rnun conjunto convexo y J: Ω ⊆Rn→Run funcional convexo. Entonces si u∈Ωes un mínimo local, se verifica: (1) ues un mínimo global. (2) Si Jes estrictamente convexo, ues el único mínimo global. Demostración. Para probar (1), se demostrará por reducción al absurdo, por lo que se supondrá que existe ¯ u∈Ωtal que J(¯ u)< J(u), de modo que uno es un mínimo global. Usando ahora que el funcional es convexo, se obtiene, J(λ¯ u+ (1 −λ)u)≤λJ(¯ u) + (1 −λ)J(u)< λJ(u) + (1 −λ)J(u) = J(u)∀λ∈[0,1] de modo que al tomar λsuficientemente pequeño, se toman puntos tan próximos a ucomo se quiera en los que el funcional Jtoma valores más pequeños que en u, lo que contradice el hecho de que ues mínimo local, de lo que se deduce el resultado. Para probar (2), usando nuevamente la reducción al absurdo, se supondrá que existe ¯ u∈ Ω−{u}de modo que J(¯ u) = J(u). Y usando la hipótesis de convexidad del funcional, J(λ¯ u+ (1 −λ)u)< λJ (¯ u)+(1−λ)J(u) = λJ (u)+ (1−λ)J(u) = J(u)∀λ∈[0,1] lo que contradice el hecho de que usea un mínimo global, luego el mínimo global es único. Corolario 1.16. Sea Ω⊆Rnun conjunto abierto, convexo no vacío y dado un funcional J: Ω ⊆Rn→Rcontinuo, coercitivo y estrictamente convexo, entonces el problema (P) tiene solución y ésta es única. 1.4. FUNCIONALES CUADRÁTICOS 7 1.4. Funcionales cuadráticos A continuación, en esta sección, se estudiarán los funcionales cuadráticos, que debido a sus propiedades, permiten analizar de forma sencilla los conceptos introducidos en la sección anterior. Definición 1.17. Sea J:Rn→Run funcional, se dice que éste es cuadrático si está definido de la siguiente forma: J(v) = 1 2vTQv −cTv donde Q= (qij)∈ Mn(R)es simétrica y c∈Rn. Observación 1.18.Nótese que en la definición anterior, no es necesario que Q∈ Mnsea simétrica. Pues bien, para cualquier v∈Rn, se tiene que vTQv =vTQ+QT 2v, y dado que Q+QTes simétrica independientemente de que Qlo sea, se puede definir el funcional cuadrático para matrices Q∈ Mn(R)generales como: J(v) = 1 4vTQ+QTv−cTv Teorema 1.19. Sea J:Rn→Run funcional cuadrático, entonces se verifica: (1) ∇J(v) = Qv −cyHJ(v) = Q, para todo v∈Rn. (2) Jes convexo, si y sólo si, Qes semidefinida positiva. (3) Jes estrictamente convexo, si y sólo si, Qes definida positiva. Demostración. (1) El funcional J:Rn→Rse escribe como: J(v) = 1 2vTQv −c vT=1 2 n X i,j=1 vjqijvi− n X i=1 civi Derivando ahora respecto de la componente viy teniendo en cuenta que Qes simétrica: dJ dvi =1 2 n X j=1 vjqij +1 2 n X j=1 qjivj−ci=1 2 n X j=1 vjqij +1 2 n X j=1 qijvj−ci= n X j=1 vjqij −ci= (Qv−c)i Por lo que, ∇J(v) = Qv −c. Y si se vuelve a derivar una segunda vez: d2J dvidvj =qij ⇒HJ(v) = Q (2) Por el teorema 1.13, se tiene: Jes convexo ⇔J(w)−J(v)− ∇J(v)T(w−v)≥0∀v,w∈Rn 8CAPÍTULO 1. OPTIMIZACIÓN SIN RESTRICCIONES entonces, aplicando este hecho y haciendo un desarrollo de Taylor (que como J(v)tiene grado dos, entonces el desarrollo de Taylor es exacto si se toma hasta la derivada segunda) se tiene, J(w) = J(v) + ∇J(v)T(w−v) + 1 2(w−v)TQ(w−v)⇒ ⇒J(w)−J(v)− ∇J(v)T=1 2(w−v)TQ(w−v)≥0∀v,w∈Rn lo cual sucede, si y sólo si, Qes semidefinida positiva. (3) Se procede de forma análoga al anterior, llegando a la desigualdad estricta (dado que el funcional es estrictamente convexo), por lo que se obtiene para Qla definición de matriz definida positiva. Observación 1.20.El resultado es análogo para matrices Q∈ Mn(R)teniendo en cuenta la observación 1.18, de la que se deduce que: ∇J(v) = Q+QT 2v−cHJ(v) = Q+QT 2 1.5. Métodos para la obtención del mínimo Hechos los estudios anteriores, es conveniente introducir algunos algoritmos que nos permitan llegar a la obtención del mínimo que se busca. Para ello, la idea fundamental, es tomar un punto inicial u(0) (próximo a la solución si es posible), y a partir de él, construir una sucesión que converja a la solución. Para ello, se diseñará un algoritmo genérico basado en dos etapas. Una de ellas, será calcular una dirección d∈Rn, tal que J(u(k)+αd)< J(u(k))para algún α∈R(a la que se denominará dirección de descenso). La otra etapa, es calcular el valor α(al que se denominará como paso) que garantiza el descenso del valor del funcional. Con este motivo, se introduce el siguiente teorema que será de ayuda. Teorema 1.21 (Condición suficiente de dirección de descenso).Dado un funcional J: Rn→Rdiferenciable, entonces si existe una dirección d∈Rntal que ∇J(u)Td<0, entonces des una dirección de descenso para el punto u. Demostración. Haciendo un desarrollo de Taylor centrado en u, se obtiene que, J(u+αd) = J(u) + α∇J(u)Td+o(α)⇒J(u+αd)−J(u) α=∇J(u)Td+o(α) α(1.3) por lo que teniendo en cuenta que des una dirección de descenso, se tiene que para un α suficientemente pequeño, J(u+αd)< J(u), por lo que haciendo el límite cuando α→0 en la expresión 1.3, se tiene que ∇J(u)Td<0. 1.5. MÉTODOS PARA LA OBTENCIÓN DEL MÍNIMO 9 Algoritmo 1.22 (Algoritmo genérico de descenso).Se procede del siguiente modo: paso 0: Se proporciona el u(0) ∈Rnpara arrancar el algoritmo (si es posible, cerca de la solución que se denotará por u), un valor de tolerancia δyk= 0. paso 1: Se calcula una dirección de descenso d(k), y un paso αk∈(0,∞)que garantiza descenso, es decir, J(u(k)+αkd(k))< J(u(k)). Y se actualiza u(k+1) =u(k)+αkd(k), k=k+ 1. paso 2: Se comprueban los criterios de convergencia, ku(k+1) −u(k)k 1 + ku(k+1)k≤δyk∇J(u(k+1))k< δ que de cumplirse se finaliza, y en caso contrario, se volvería al paso 1. Observación 1.23.Si se desea consultar un pseudocódigo de este algoritmo genérico, se puede ver un ejemplo en el Anexo I. 1.5.1. Cálculo del paso Para calcular el paso αk, es importante definir la función j: [0,∞)→Rcomo sigue: j(α) = J(u(k)+αd(k)) El objetivo, es hallar un αk∈[0,∞)que cumpla que j(αk)< j(0), por lo que se plantea un problema de minimización unidimensional. Para resolverlo, nótese que si Jes suficientemente diferenciable, entonces, j0(α) = d(k)T ∇J(u(k)+αd(k))⇒j0(0) = d(k)T ∇J(u(k))<0 j00(α) = d(k)T HJ(u(k)+αd(k))d(k) y usando estos cálculos, se pueden plantear diferentes situaciones, y a su vez distintas formas de cálculo del paso. Paso óptimo: El paso óptimo, se denota al paso exacto que da el óptimo. Éste vendrá caracterizado por j0(α)=0, y en algunos casos particulares, como los funcionales cuadráticos definidos por una matriz simétrica definida positiva, es posible calcularlo: j0(α) = dT(Q(u+αd)−c) = dTQu+αdTQd−dTc⇒j0(α)=0⇔α=−dT(Qu −c) dTQd 10 CAPÍTULO 1. OPTIMIZACIÓN SIN RESTRICCIONES Dicotomia La situación anterior, no es nada común en la realidad. En su defecto, si se conoce un intervalo (a, b)⊂(0,+∞)en el que se encuentra el paso αk. Una buena opción sería calcular la solución al problema j0(α)=0mediante el método de dicotomia (o de Newton- Raphson si j∈ C2(a, b)). Además, si no se conoce el intervalo (a, b), éste se puede aproximar mediante algún procedimiento numérico. Para conocer detalles sobre su implementación, se puede consultar [7] o en [17, págs. 38-57], en el que se aportan dos versiones, una en la que se utiliza la derivada (j0(α)) y otra en la que no es necesario evaluarla. De dichas versiones, se puede consultar sus pseudocódigos en el Anexo I. Algoritmos de criterio de paso grande y pequeño Uno de los algoritmos más empleados, se puede englobar dentro de una familia de métodos cuyo objetivo será elegir un paso que no sea ni demasiado pequeño (a lo que se conocerá como criterio de paso pequeño (CPP)), ni demasiado grande (a lo que se conocerá como criterio de paso grande (CPG)). Por tanto, se dirá que un paso es admisible (CPA) si no verifica CPP ni CPG. De este modo, toda esta familia de reglas, siguen el siguiente algoritmo genérico: Algoritmo 1.24 (Regla general de CPP y CPG).Se sigue el siguiente procedimiento: paso 0: Se aporta un intervalo inicial (αmin, αmax)de tal manera que αmin cumpla CPP yαmax cumpla CPG. paso 1: Se toma α= (αmin +αmax)/2y se comprueba si αcumple CPP. En caso afirmativo, se actualiza αmin =αy se vuelve al paso 1, y en caso contrario, se pasa al paso paso 2. paso 2: Se comprueba si αcumple el CPG. En caso afirmativo, se actualiza αmax =α, y en caso contrario, αes admisible, por lo que se finaliza el algoritmo y se toma αcomo solución. Observación 1.25.En el paso 1 del algoritmo 1.24, se toma como αde prueba el punto medio del intervalo. Pero en realidad, se podría generar un punto de prueba mediante cualquier otra metodología como la regla de la secante u otras. Por otra parte, en el paso 0, en vez de aportar un intervalo cumpliendo dichas condiciones, se puede hacer un procedimiento previo para inicializar el algoritmo con la búsqueda numérica de un intervalo verificando las condiciones requeridas (se puede ver un pseudocódigo de esta incialización en el Anexo I). El procedimiento anterior, proporciona una gran variedad de reglas de cálculo del paso en función del los diferentes CPP y CPG que se podrían tomar. En este marco, las tres 1.5. MÉTODOS PARA LA OBTENCIÓN DEL MÍNIMO 11 reglas más conocidas y usadas, son las siguientes: 1. Regla del Armijo: No cuenta con CPP y se dice que un paso αes demasiado grande si j(α)> j(0) + ρj0(0)αpara un valor de ρ∈(0,1/2). 2. Regla de Goldstein: Un paso αcumple CPP si j(α)< j(0) + m2j0(0)α, y verificará el CPG si j(α)> j(0) + m1j0(0)α, donde m1∈(0,1/2) ym2∈(1/2,1). Además, en la práctica se suele usar m2= 1 −m1. 3. Regla de Wolfe-Powell: Un paso αcumple el CPP si j0(α)< m2j0(0) ,y verificará el CPG si j(α)> j(0) + m1j0(0)α, donde 0< m1< m2<1. Figura 1.1: Regla de paso con CPG y CPP En la figura, se muestra la idea que subyace detrás de todas estas reglas. Tomando como ejemplo el criterio de Goldstein, aquellos puntos que cumplan los criterios CPP (los que están por debajo de la línea CPP) y CPG (los que se encuentran por encima de la recta CPG), dejan de ser admisibles y se muestran con una línea en rojo. Por el contrario, los puntos que están entre ambas rectas, son puntos admisibles y se representan con en una línea en verde. 1.5.2. Método del gradiente Recibe este nombre el método que toma como dirección de descenso d(k)=−∇J(u(k)), y en base a esta dirección, se seguirá el algoritmo 1.22, tomando como paso αk, el obtenido al usar alguna de las reglas descritas en el apartado 1.5.1. En caso de poder calcular el αkóptimo, el método recibirá el nombre de gradiente con paso óptimo. Método en el cuál, se utiliza la dirección de máximo descenso. 1.5.3. Método del gradiente conjugado Éste método, mejora el visto en el apartado 1.5.2. Esencialmente, está pensado en un primer momento para funcionales cuadráticos. Su convergencia se produce en a lo sumo n iteraciones independientemente del punto inicial dado para arrancar el método, por lo que, no se puede considerar como tal un método iterativo. Su fundamento, se basa en tomar direcciones conjugadas, es decir, tales que, d(i)TQd(j)= 0 con i6=j i, j ∈ {0,1, . . .} 12 CAPÍTULO 1. OPTIMIZACIÓN SIN RESTRICCIONES y explícitamente, el algoritmo vendría dado por: Algoritmo 1.26 (Gradiente conjugado en funcionales cuadráticos).Se sigue el siguiente procedimiento: paso 0: Se toma un u(0) ∈Rnde arranque , con su residuo dado por r(0) =Qu(0) −cy se toma como dirección de descenso d(0) =−r(0),k= 0. paso 1: Se calcula el paso que viene dado explícitamente por, αk=r(k)Tr(k) d(k)T Qd(k) , con el que se actualiza el iterante y se toma u(k+1) =u(k)+αkd(k). Luego, se actualiza r(k+1) = r(k)+αkQd(k). paso 2: Si r(k+1) =0se finaliza, y en caso contrario, se toma βk+1 =r(k+1) Tr(k+1) r(k)Tr(k), con el que se actualiza la nueva dirección d(k+1) =−r(k+1) +βk+1d(k), se toma k=k+ 1 y se vuelve al paso 1. Observación 1.27.Para consultar con detalle los aspectos que hay detrás de este algoritmo (deducción, convergencia, implementación, etc) se puede consultar [13, págs.102-120]. Este caso, es obviamente para el caso en el que el funcional objetivo es cuadrático. Por tanto, una cuestión interesante sería si esta metodología se podría ampliar al caso no cuadrático. La respuesta es positiva y se plantean dos propuestas: 1. Fletcher-Reeves: Partiendo de d(0) =−∇J(u(0)), en cada k-iteración se toma: βk=k∇J(u(k+1))k2 k∇J(u(k))k2yd(k+1) =−∇J(u(k+1)) + βkd(k) 2. Polak-Ribiére: Partiendo de d(0) =−∇J(u(0)), en cada k-iteración se toma: βk=k∇J(u(k+1))k2 k∇J(u(k))k2yd(k+1) =−∇J(u(k+1)) + βkd(k) Obviamente, en los dos casos anteriores, con dichas expresiones de la dirección de descenso, se procede usando el algoritmo 1.22, combinándolo con los vistos en 1.5.1 para calcular el paso. 1.5.4. Método de Newton Otra forma de abordar el problema, sería calculando el u∈Rnque cumpla que ∇J(u) = 0, dado que se vio que esta condición es indispensable para ser mínimo (sin restricciones). Para ello, se utiliza el método de Newton-Raphson aplicado a la función ∇J(v), (u(0) dado u(k+1) =u(k)−HJ(u(k))−1∇J(u(k))k≥0 1.5. MÉTODOS PARA LA OBTENCIÓN DEL MÍNIMO 13 la cual, se enmarca dentro del marco del algoritmo 1.22, tomando αk= 1 yd(k)= −HJ(u(k))−1∇J(u(k)). Por tanto, este método no garantiza descenso en todas sus iteraciones. Es más, la dirección d(k)no tendría porque ser de descenso. Por ello, se puede mejorar el método introduciendo alguna regla de cálculo de paso en vez de tomar paso fijo 1. Por otro lado, el coste computacional es demasiado elevado debido a tener que calcular la inversa de la matriz hessiana en cada iteración. De modo que para solventar esta situación, se puede proceder del siguiente modo,        u(0) dado HJ(u(k))w=−∇J(u(k)) u(k+1) =u(k)+w que se vuelve a enmarcar dentro del algoritmo 1.22 tomando d(k)=wyαk= 1. Y nuevamente, se puede mejorar el método combinándolo con alguna de las reglas de cálculo de paso. 1.5.5. Métodos Cuasi-Newton Sin embargo, en el método de Newton resulta costoso debido a tener que evaluar continuamente la Hessiana. Para resolver este problema, se introducen los métodos Cuasi- Newton que tienen el mismo fundamento que el visto en 1.5.4 con la peculiaridad de usar una aproximación de la hessiana del funcional. Para esta aproximación, se toma una matriz de arranque H0definida positiva y simétrica (la identidad o una aproximación mediante diferencias finitas), para luego calcular en cada iteración una actualización de dicha matriz mediante algunas de las reglas de actualización que que se introducen a continuación,                Hk+1 =Hk+s(k)s(k)T s(k)Ty(k)−Hky(k)y(k)THk y(k)THky(k)Actualización D.F.P Hk+1 = I−s(k)y(k)T s(k)Ty(k)!Hk I−y(k)s(k)T s(k)Ty(k)!+s(k)s(k)T s(k)Ty(k)Actualización B.F.G.S donde yk=∇J(u(k+1))− ∇J(u(k))ysk=u(k+1) −u(k). De esta forma, tomando como aproximación de la hessiana, la dada por la fórmula B.F.P o B.F.G.S en cada k-iteración, se definiría la dirección d(k)como se vio en el apartado 1.5.4. Así pues, tomando dicha dirección y combinándola con una regla de cálculo de paso, se obtendría una gran variedad de métodos de Cuasi-Newton al aplicar el algoritmo 1.22. 14 CAPÍTULO 1. OPTIMIZACIÓN SIN RESTRICCIONES Finalmente, una propiedad interesante que puede resultar útil en ocasiones es que bajo ciertas hipótesis, las actualizaciones de las inversas de la hessiana, pueden mantener el carácter de definidas positivas. El siguiente resultado lo muestra: Teorema 1.28. Si Hkes definida positiva, ∇J(u(k+1))6=0, y el paso αkes elegido de tal forma que el iterante u(k+1) verifica: ∇J(u(k))Td(k)<∇J(u(k+1))Td(k) (lo cual es equivalente a que s(k)Ty(k)>0), entonces la matriz Hk+1 definida por un método B.F.P ó B.F.G.S, es definida positiva. Demostración. Se puede consultar en [3, págs.60-61]. Por otro lado, en casos donde la dimensión del problema es muy elevada, resulta prohibitivo implementar este algoritmo. Esto es debido al elevado coste de operar y guardar matrices con dimensiones muy grandes. Por este motivo, en la práctica se opta por una versión de coste menor al usar la actualización B.F.G.S (ya que es la que mejores resultados aporta). En dicha versión (donde se opta por un algoritmo recursivo), la dirección de descenso usando la actualización B.F.G.S, se obtiene almacenando únicamente un número mde los s(k)ey(k)calculados en las iteraciones anteriores. Así pues, considerando una k-iteración cualquiera, se podría tomar el siguiente algoritmo para obtener d(k)(la dirección de descenso), Algoritmo 1.29 (Cálculo de la dirección de descenso B.F.G.S).Tomando de forma inicial q=∇J(u(k)), se procede como sigue: paso 1: Definiendo ρk=1 y(k)Ts(k), se procede con un bucle en i=k−1, . . . , k−m, donde se toma γi=ρisiTq,q=q−γiy(i). paso 2: Definiendo H0 k=s(k−1) Ty(k−1) y(k−1) Ty(k−1) Iy tomando r=H0 kq, se procede con un bucle en i=k−m, . . . , k−1, donde se toma β=ρiy(i)Tr, con el que se actualiza r=r+s(i)(γi−β). paso 3: Se concluye tomando como dirección de descenso d(k)=−r. Mediante esta forma, se suaviza el coste computacional en dimensiones muy elevadas. De hecho, se ha visto en la práctica que se obtienen resultados muy buenos para un valor mentre 3 y 20. En [13, págs.176-185] se pueden ver más detalles sobre esta forma de implementación de la fórmula B.F.G.S. denominada B.F.G.S de memoria limitada, y se puede consultar un pseudocódigo del método Cuasi-Newton con una fórmula B.F.G.S de coste reducido en el Anexo I. Capítulo 2 Optimización con restricciones: Conceptos teóricos En el anterior capítulo, se ha visto como se abordan los problemas de optimización relativos al problema (P). Pero este tipo de problemas sólo se presentan en situaciones muy idílicas que no suceden en la realidad. Por consiguiente, se pretende a partir de ahora presentar un nuevo tipo de problema que englobe las peculiaridades que surgen en la realidad, y a partir de estos problemas, desarrollar nuevos métodos que nos permitan resolverlos ó en su lugar, reducirlos a problemas equivalentes del tipo (P)para solucionarlos como en el capítulo anterior. 2.1. Presentación del problema Con motivo de lo que se ha explicado, se presenta el siguiente problema de minimización con restricciones, al que se referirá de ahora en adelante como problema (Q): m´ın v∈RnJ(v) Sujeto a: ϕi(v)≤0i= 1, . . . , m φj(v)=0 j= 1, . . . , p (Q) Como se puede observar, en este nuevo problema no tiene porque ser solución el mínimo global de J(v), ya que será solución del problema aquel valor que minimice el funcional cumpliendo las exigencias marcadas por las funciones ϕiyφj. Estas funciones se conocen como restricciones y delimitan el conjunto donde se pretende minimizar el funcional objetivo, al que se denomina como conjunto factible óconjunto 15 22 CAPÍTULO 2. OPTIMIZACIÓN CON RESTRICCIONES, TEORÍA 2.3. Multiplicadores de Lagrange En esta sección, se considerará un caso particular del problema (Q), en el que se tiene el problema de minimización de un funcional objetivo J:Rn→Rsujeto a restricciones de tipo igualdad, por lo que no se tienen restricciones de tipo desigualdad. A este caso, se denotará como problema (Qig), y se verá que existe la posibilidad de su resolución mediante un sistema de ecuaciones no lineales. Para ello, se comienza introduciendo la notación pertinente. En primer lugar, a partir de las funciones que definen las restricciones, se define la función Φ : Rn→Rp, que viene dada por, Φ(v) = (φ1(v), φ2(v), . . . , φp(v))T de modo que el conjunto factible se puede definir de la siguiente manera: S={v∈Rn: Φ(v) = 0} Definición 2.12. Un punto u∈Ses regular para las restricciones {φj(v)}p j=1 si los vectores gradientes ∇φ1(u),...,∇φp(u)son linealmente independientes. Observación 2.13.La definición anterior, no se puede dar si p>n. Teorema 2.14. Dado un punto regular u∈S, el espacio tangente a Sen dicho punto es: T(u) = v∈Rn: Φ0(u)v=0 donde Φ0(u)se corresponde con: Φ0(u) = (∇φ1(u)|. . . |∇φp(u))T Demostración. Se puede consultar en [16, pág. 302]. Lema 2.15. Dado u∈Sun punto regular, supóngase que dicho ues un mínimo de J en S, y que Jes derivable en u. Entonces, para todo v∈Rn, que verifique Φ0(u)v=0 también se tiene que h∇J(u),vi= 0. Demostración. Se puede ver en [16, pág. 304]. Teorema 2.16. Considerando el problema (Qig)donde el funcional y las restricciones son diferenciables, si dado u∈Sun punto regular, éste es un mínimo del funcional J en S. Entonces, existen pnúmeros λi, que se conocerán como multiplicadores de Lagrange asociados a u, definidos de forma única, tales que: ∇J(u) + λ1∇φ1(u) + . . . +λp∇φp(u) = 0 2.3. MULTIPLICADORES DE LAGRANGE 23 Demostración. Si ∇J(u) = 0, entonces la propiedad es obvia tomando λ1=. . . =λp= 0. Entonces, suponiendo que ∇J(u)6=0, haciendo uso del lema anterior, si v∈ker (Φ0(u)), entonces v∈ker ∇J(u)T, por lo que, ker (Φ0(u)) ⊆ker ∇J(u)T. De este modo, tomando el espacio ortogonal, se tiene el contenido contrario, es decir, ker ∇J(u)T⊥⊆ ker (Φ0(u))⊥. Ahora bien, es obvio que ∇J(u)T∈ker ∇J(u)T⊥, por lo que, ∇J(u)T∈ ker (Φ0(u))⊥=Im (Φ0(u)), y entonces, existen coeficientes λ1, . . . , λp∈Rtales que, ∇J(u) = λ1∇φ1(u) + . . . +λp∇φp(u) lo cual es evidente que es equivalente a la tesis del teorema. Ésta, no es la única demostración posible del teorema. Pues existen múltiples versiones para llegar a obtener el mismo resultado. Por ejemplo, en [1, págs. 150-152] se da una demostración haciendo uso del Teorema de la función implícita. Por otra parte, la condición necesaria obtenida, permite obtener un método para la resolución del problema (Qig)mediante la resolución de un sistema de ecuaciones no lineales generalmente. Para ello, se toman como incógnitas (λi1≤i≤p)yu. Ahora bien, como u∈Rnes un vector de ncomponentes, en total se tienen n+p incógnitas. Las cuales, quedarán determinadas por las ecuaciones de las restricciones y la condición de Lagrange, es decir, (∇J(u) + Φ0(u)Tλ=0 Φ(u) = 0(2.3) siendo λel vector que tiene por componentes a los multiplicadores de Lagrange. Escribiendo dicho sistema de una forma más visual, quedaría:                            ∂ ∂v1 J(u) + λ1 ∂ ∂v1 φ1(u) + . . . +λp ∂ ∂v1 φp(u)=0 . . . ∂ ∂vn J(u) + λ1 ∂ ∂vn φ1(u) + . . . +λp ∂ ∂vn φp(u)=0 φ1(u)=0 . . . φp(u) = 0 Observación 2.17.Si se toma el Lagrangiano asociado al problema, que se define de la siguiente forma, L(v,ξ) = J(v) + Φ(v)Tξv∈Rnξ∈Rp 24 CAPÍTULO 2. OPTIMIZACIÓN CON RESTRICCIONES, TEORÍA el sistema 2.3, se puede reescribir en forma compacta como:          L ∂v(u,λ) = ∇J(u) + Φ0(u)Tλ=0 ∂L ∂ξ(u,λ) = Φ(u) = 0 De modo que resolver el sistema de ecuaciones planteado, es equivalente a obtener un punto crítico asociado al problema. Sin embargo, la resolución del sistema no aporta exactamente la solución al problema de minimización. Esto es debido a que la condición aportada por el teorema es necesaria pero no suficiente. Por ello, una vez obtenidos los valores de λyuque verifican las ecuaciones, es necesario evaluar en un entorno de uel funcional para determinar el tipo de punto crítico que se ha obtenido. Teorema 2.18 (Condición suficiente de Lagrange).Se considera el problema (Qig), bajo las hipótesis de derivabilidad del funcional objetivo y de las restricciones. Entonces, suponiendo que u∈Ses un punto verificando la condición de Lagrange con multiplicadores de Lagrange λ, verificando que, dTHLv(u,λ)d>0para todo d/∈0con Φ(u)0d=0 entonces ues un mínimo local estricto del problema (Qig). Demostración. Se puede consultar en [10, pág 272-273]. Este teorema, aporta un caso en el que la condición de Lagrange es suficiente, y por tanto, bajo dichas hipótesis, no habría que comprobar si la solución obtenida es mínimo o no. Por otra parte, los multiplicadores de Lagrange aquí introducidos, tienen diversas aplicaciones en distintas ramas del conocimiento. Por ejemplo, la física o la economía son claros ejemplos de ello, y se puede ver en [9] o [20]. 2.4. Condiciones de Karush-Kuhn-Tucker y Fritz-John En esta sección, se verá una generalización de la condición de Lagrange vista en 2.16, al caso que incluye restricciones de desigualdad. Para ello, se hace uso de las condiciones de Karush-Kuhn-Tucker y posteriormente de las condiciones de Fritz-John que se definen a continuación. 2.4. CONDICIONES DE KARUSH-KUHN-TUCKER Y FRITZ-JOHN 25 Definición 2.19. Considerando el problema (Q), se define su función Lagrangiano asociada, L(v, ν, ξ) : Rn×Rm×Rp→Rcomo sigue: L(v,ν,ξ) = J(v) + m X i=1 νiϕi(v) + p X j=1 ξjφj(v) Definición 2.20. Dado el problema (Q), se supone que el funcional objetivo y las restricciones son derivables. Entonces, se dice que un punto u∈Rnes un punto de Karush- Kuhn-Tucker para el problema (Q), si y sólo si, existen multiplicadores de Lagrange y de Karush-Kuhn-Tucker (a los que se denotará por las letras griegas λyµrespectivamente) verificando las condición!de KKT (KKT) que son las siguientes: 1.- Condición estacionaria: ∇vL(u,µ,λ) = ∇J(u) + m X i=1 µi∇ϕi(u) + p X j=1 λj∇φj(u) = 0(2.4) 2.- Condición de factibilidad (ϕi(u)≤0i= 1, . . . , m φj(u)=0 j= 1, . . . , p (2.5) 3.- Condición de holgura: (µiϕi(u)=0 i= 1, . . . , m µi≥0i= 1, . . . , m (2.6) Para calcular ahora los puntos que cumplen las condiciones de Karush-Kuhn-Tucker (KKT), se puede proceder en dos pasos. El primero, es la resolución de un sistema de ecuaciones (no lineal generalmente), y que se corresponden con la imposición de la condición estacionaria, de factibilidad para las restricciones de igualdad y la de holgura:            ∂ ∂vk J(v) + p X j=1 λj ∂ ∂vk φj(v) + m X i=1 µi ∂ ∂vk ϕi(v)=0 k= 1, . . . , n φj(v)=0 j= 1, . . . , p µiϕi(v)=0 i= 1, . . . , m Nótese que se trata de un sistema de (n+m+p)ecuaciones y (n+m+p)incógnitas. Estas serían, las ncomponentes de u, los pmultiplicadores de Lagrange y los m multiplicadores KKT. Una vez resuelto, el segundo paso sería ver cuáles de las soluciones obtenidas son puntos de KKT. Para ello, hay que comprobar que son puntos factibles, es decir, comprobar que cumplen las restricciones de desigualdad, y finalmente, que todos los multiplicadores KKT son positivos. 26 CAPÍTULO 2. OPTIMIZACIÓN CON RESTRICCIONES, TEORÍA Llegados a este punto, se está en condiciones de probar que todo mínimo del problema (Q), verifica las condiciones KKT. Para ello, se introduce la siguiente versión del Lema de Farkas-Mirkonski que será de utilidad en su demostración. Lema 2.21 (Lema de Farkas-Mirkonski).Considerando el problema (Q)y un punto factible u∈S, el conjunto, Z=       d∈Rn: dT∇J(u)<0 dT∇ϕi(u)≤0i= 1, . . . , m dT∇φj(u) = 0 j= 1, . . . , p        es vacío, si y sólo si, existen multiplicadores de Lagrange y de KKT verificando la condición estacionaria en u. Demostración. Se puede ver en [15, págs. 53-54,391] y en [12, pág. 313-314]. Teorema 2.22 (Karush-Kuhn-Tucker).Sea uun mínimo local del problema (Q). Si se cumple la siguiente condición, a la que se denominará como condición de cualificación, DF S(u, S) = DFL(u, S)(2.7) entonces existen multiplicadores de Lagrange y de KKT verificando las condiciones KKT. Demostración. Dado que ues un mínimo local del problema, éste verifica las condiciones de factibilidad, es decir, 2.5. Ahora bien, sea d∈DFS(u, S), por el teorema 2.9, se tiene que dT∇J(u)≥0. Pero aplicando la condición 2.7, d∈DFL(u, S), y así el sistema,        dT∇J(u)<0 dT∇ϕi(u)≤0i∈I(u) dT∇φj(u)=0 j= 1, . . . , p no tiene solución y por tanto el conjunto Zdefinido en el lema de Farkas–Mirkonski es vacío, y por tanto, se obtiene que existen multiplicadores de Lagrange y de KKT (µi≥0) verificando que, ∇J(u) + X i∈I(u) µi∇ϕi(u) + p X j=1 λj∇φj(u)=0 por lo que tomando µi= 0 i /∈I(u), se obtiene la condición estacionaria y de holgura. De modo que se obtiene el resultado. Esta demostración que se acaba de hacer, no es única. De hecho, existen muy diversas formas de demostrar el Teorema de Karush-Kuhn-Tucker, por ejemplo, se puede consultar 2.4. CONDICIONES DE KARUSH-KUHN-TUCKER Y FRITZ-JOHN 27 [2] para ver otro tipo de demostración donde se utiliza el Teorema de Weierstrass. Por otro lado, la demostración mostrada en este trabajo, requiere de la condición de cualificación 2.7, pero ésta no es única. De hecho, puede ser substituida por otras equivalentes. Con este objetivo, se introduce el siguiente teorema. Teorema 2.23. Sea v∈Sun punto admisible. Si los vectores ∇φj(v)con j= 1, . . . , p y ∇ϕi(v)con i∈I(v)son linealmente independientes, entonces la condición de cualificación 2.7 se cumple, es decir, DF S(v, S) = DFL(v, S). Demostración. Se puede consultar en [15, pág 394-396]. Observación 2.24.Existe una gran variedad de condiciones de cualificación que sirven de utilidad para probar el teorema 2.22, para ello, ver [14, sec. 6.3]. Hasta el momento, se han introducido los multiplicadores KKT, pero también se podría hablar de los multiplicadores de Fritz-John. Éstos, se podrían considerar una generalización de los anteriores, pero que en algunos casos resultan de utilidad. Definición 2.25. Dado el problema (Q), suponiendo que tanto el funcional como las restricciones son derivables. Entonces, se dice que un punto u∈Rnverifica las condiciones de Fritz-John (ó es un punto de Fritz-John), si y sólo si, existen multiplicadores de Lagrange (λ) y de Fritz-John (η0,η)tales que:                  η0∇J(u) + m X i=1 ηi∇ϕi(u) + p X j=1 λj∇φj(u)=0 ϕi(u)≤0yφj(u) = 0 para i= 1, . . . , m j = 1, . . . , p ηiϕi(u)=0 para i= 1, . . . , m η0, ηi≥0para i= 1, . . . , m y además (η0,η,λ)6=0 Como se puede observar, éstas son las mismas condiciones KKT, salvo que en este caso, ∇J(u)está multiplicado por un elemento η0≥0. Por tanto, si η0>0, al dividir la primera ecuación de las condiciones de Fritz-John por η0, se obtienen las condiciones KKT. Por tanto, la única diferencia existirá cuando η0= 0. A continuación, se aporta un resultado similar al visto en el teorema 2.22 para las condiciones de Fritz-John. Teorema 2.26 (Condición necesaria de Fritz-John).Dado el problema (Q), suponiendo que tanto el funcional objetivo como las restricciones sean diferenciables. Si u∈Ses un mínimo local del problema (Q), entonces existen multiplicadores de Lagrange y de Fritz- John, cumpliendo las condiciones de Fritz-John. Demostración. Si los vectores ∇φj(u)con j= 1, . . . , m y∇ϕi(u)con i∈I(u)son linealmente dependientes, entonces existen elementos λ1, . . . , λpyηicon i∈I(u), tales 28 CAPÍTULO 2. OPTIMIZACIÓN CON RESTRICCIONES, TEORÍA que Pi∈I(u)ηi∇ϕi(u) + Pp j=1 λj∇φj(u)=0, siendo alguno de ellos distintos de cero. Por tanto, tomando η0=ηi= 0 con i /∈I(u), se obtienen inmediatamente las condiciones de Fritz-John. Si se supone ahora que los vectores ∇φj(u)con j= 1, . . . , m y∇ϕi(u)con i∈I(u)son linealmente independientes, aplicando el teorema 2.23, se obtiene la condición de cualificación 2.7, y por tanto, se está en condiciones del teorema 2.22. Por consiguiente, se cumplen la condiciones KKT, luego, tomando η0= 1 yη=µ, se obtienen también las condiciones de Fritz-John. Entonces, un punto cumpliendo las condiciones de Fritz-John, es un candidato a mínimo (al igual que sucedía con los puntos KKT). Este hecho, se usará posteriormente en el último capítulo dado que algún método numérico permitirá obtener un punto de Fritz-John. Sin embargo, cuando η0= 0, las condiciones que nos aporta el teorema 2.26, no dependen del ∇J(u), lo que no aporta información sobre el funcional objetivo en el mínimo. Por este motivo, es poco interesante la condición de Fritz-John, y a partir de ahora, se estudiarán los puntos KKT en los que si se ve involucrado ∇J(u). Continuando pues con el estudio de los puntos KKT, en el teorema 2.22, se proporciona una condición necesaria de mínimo. Pues cualquier mínimo cumple las condiciones KKT, sin embargo, ésta no es suficiente en general. Para poder introducir condiciones suficientes, es necesario reforzar las hipótesis sobre el problema. Con este motivo, se introduce la siguiente definición necesaria para introducir una condición suficiente KKT. Definición 2.27. Dado el problema de optimización (Q), se dirá que éste es convexo si el funcional objetivo y el conjunto admisible son convexos. Además, se representará a este tipo especial de optimización como problema (QC). Observación 2.28.En el caso del problema (QC), se supondrá que las restricciones ϕi(v) son convexas y las φj(v)lineales. Estas suposiciones, derivan del hecho de que al ser de estos modos concretos, el conjunto factible que generan es convexo. Teorema 2.29 (Condición suficiente de Karush-Kuhn-Tucker).Dado un problema de optimización convexa (QC), se tiene que dado un punto admisible u∈S, para el cual existen multiplicadores de Lagrange y KKT (λyµrespectivamente). Entonces dicho punto es un mínimo del problema (QC). Demostración. Se considera un punto u∈Sverificando las condiciones KKT, sean λy µsus multiplicadores de Lagrange y KKT asociados. La función L(v,µ,λ)es convexa para cualquiera v. Entonces, usando las propiedades de las funciones convexas 1.13, y las 2.4. CONDICIONES DE KARUSH-KUHN-TUCKER Y FRITZ-JOHN 29 condiciones KKT, se tiene que para cualquier punto factible v, L(v,µ,λ)≥ L(u,µ,λ) + (v−u)T∇vL(u,µ,λ) = L(u,µ,λ) = =J(u) + m X i=1 µiϕi(u) + p X j=1 λjφj(u) = J(u)(2.8) pero como ves un punto factible y µi≥0para i= 1, . . . , m, entonces se deduce que, µiϕi(v)≤0i= 1, . . . , m λjφj(v)=0 j= 1, . . . , p )⇒ L(v,µ,λ)≤J(v) por lo que usando 2.8 se obtiene que J(v)≥J(u), y por consiguiente, ues un mínimo. Vistos estos resultados, la interpretación geométrica que se puede hacer de ellos y de las condiciones KKT, se muestra en la siguiente gráfica. En ella, se presentan dos restricciones de desigualdad que delimitan un conjunto factible. En discontinuo, se presentan los gradientes de las funciones en el punto (1,1) que se supone el mínimo en el conjunto para un cierto funcional. 0 0.5 1 1.5 0 0.5 1 1.5 Figura 2.4: Condiciones KKT De este modo, suponiendo que se está en las hipótesis del teorema 2.22, el gradiente del funcional en el mínimo tiene una dependencia lineal respecto de los gradientes de las restricciones activas. Ahora bien, como las dos restricciones son activas y los multiplicadores KKT son positivos, el conjunto de líneas finas representa el cono de los posibles valores que puede tomar −∇J(1,1). Lo cual, es lógico desde el punto de vista intuitivo, dado que −∇J(1,1) es la dirección de máximo descenso en dicho punto, y si ésta fuese factible, el (1,1) no podría ser mínimo local. De este modo, se tiene una información relativa sobre la dirección de máximo descenso en el mínimo, y a partir de ella, de todas las posibles direcciones de descenso. Por otro lado, hasta el momento sólo se han visto condiciones de primer orden que caracterizan los mínimos. Sin embargo, las condiciones KKT nos van a permitir ampliar la información vista hasta el momento, dando condiciones de segundo orden. 30 CAPÍTULO 2. OPTIMIZACIÓN CON RESTRICCIONES, TEORÍA Definición 2.30. Sea v∈Rnverificando las condiciones de KKT, si existen sucesiones nd(k)ok∈N→dy{δk}k∈N→0tales que v(k)=v+δkd(k)∈S, verificando, ϕi(v(k))≤0i∈I(v)\I+(v) ϕi(v(k))=0 i∈I+(v)yφj(v(k))=0 j= 1, . . . , p donde I+(v) = {i:i∈I(v)tales que µi>0}. Entonces se dice que des una dirección de restricción secuencial nula, y al conjunto de dichas restricciones se denota por S(v,µ,λ)⊆ DF S(v, S). Definición 2.31. Dado un punto v∈Rnverificando las condiciones KKT (con multiplicadores de Lagrange y KKT λyµrespectivamente), si d∈DFL(v, S)yµidT∇ϕi(v) = 0∀i∈I(v), entonces des una dirección de restricción nula linealizada. Al conjunto de dichas restricciones se denota por G(v,µ,λ), y se puede escribir como: G(v,µ,λ) = (d∈Rn:d∈DF L(v, S) dT∇ϕi(v) = 0 i∈I+(v))⊆DFL(v, S) Teorema 2.32 (Condición necesaria de segundo orden).Sea u∈Rnun mínimo local del problema (Q). Si se cumple la condición de cualificación 2.7, se tiene, dTHvL(u,µ,λ)d≥0∀d∈S(u,µ,λ) donde HvL(u,µ,λ)denota la hessiana del Lagrangiano respecto de vevaluada en (u,µ,λ). Demostración. Para cualquier d∈S(u,µ,λ), si d=0el resultado es obvio, por lo que se supone d6=0. Tomando las sucesiones nd(k)ok∈N→dy{δk}k∈N→0en la forma de la definición 2.30, se tiene que Pm i=1 µiϕi(u(k)) + Pp j=1 λjφj(u(k))=0, por lo que, usando este hecho y las condiciones de KKT (pues ues un punto KKT por el teorema 2.22), se llega a: J(u(k)) = J(u+δkd(k)) = L(u+δkd(k),µ,λ) = L(u,µ,λ) + δkd(k)T∇vL(u,µ,λ)+ +1 2δ2 kd(k)THvL(u,µ,λ)d(k)+o(δ2 k) = J(u) + 1 2δ2 kd(k)THvL(u,µ,λ)d(k)+o(δ2 k) (2.9) Pero como ues un mínimo local estricto, para un ksuficientemente grande, se tiene que J(u+δkd(k))≥J(u). De modo que, usando este hecho, 2.9 y aplicando límites, se llega a, dTHvL(u,µ,λ)d≥0 como se quería demostrar. 2.4. CONDICIONES DE KARUSH-KUHN-TUCKER Y FRITZ-JOHN 31 Observación 2.33.Si se tiene que S(u,µ,λ) = G(u,µ,λ), entonces se tendría de forma evidente que dTHvL(u,µ,λ)d≥0para todo d∈G(u,µ,λ). Teorema 2.34 (Condición suficiente de segundo orden).Dado u∈Rnun punto cumpliendo las condiciones de KKT para el problema (Q). Si se verifica que, dTHvL(u,µ,λ)d>0∀d∈G(u,µ,λ) entonces ues un mínimo local estricto del problema. Demostración. Se supone que uno es mínimo, por lo que existe una sucesión de elementos factibles u(k)k∈N→utales que J(u(k))≤J(u). Y sin pérdida de generalidad, se puede asumir que, (u(k)−u ku(k)−uk)k∈N→d por lo que aplicando el razonamiento de la demostración del teorema 2.10, se obtiene que dT∇J(u)≤0,. Usando ahora las condiciones KKT, ∇J(u) = − m X i=1 µi∇ϕi(u)− p X j=1 λj∇φj(u) que al multiplicar por d, si se tiene en cuenta que d∈DFS(u, S)⊆DFL(u, S), entonces, dT∇J(u) = − m X i=1 µidT∇ϕi(u)− p X j=1 λjdT∇φj(u)≥0 de modo que se obtuvieron las dos desigualdades, y por tanto dT∇J(u)=0, lo que implica que Pi∈I(u)µidT∇ϕi(u)=0, y por consiguiente, d∈G(u,µ,λ). Ahora bien, teniendo en cuenta que J(u(k))≤J(u), se puede hacer el siguiente desarrollo, L(u,µ,λ)≥ L(u(k),µ,λ) = L(u,µ,λ) + 1 2δ2 kd(k)THL(u,µ,λ)d(k)+o(δ2 k) por lo que al dividir por δ2 ky tomar límites, se llega a que dTHL(u,µ,λ)d≤0, lo que es una contradicción con la hipótesis, y por consiguiente se obtiene el resultado. Este último resultado, establece una condición suficiente similar a la que se introdujo en el corolario 2.11, con la peculiaridad de que se hace uso de la Hessiana del Lagrangiano. Pues bien, las condiciones KKT casi implican que dT∇J(u)>0salvo en las direcciones de restricción nula linealizadas. Por este motivo, se requiere una condición de segundo orden en ellas. 38 CAPÍTULO 2. OPTIMIZACIÓN CON RESTRICCIONES, TEORÍA De la forma en que se definió la función dual en 2.39, se tiene que ésta viene determinada por, h(ν) = m´ın v∈Rn1 2vTQv −cTv+µT(Av −b)(2.11) pero como se ha supuesto que Qes definida positiva, para un µcualquiera, la función 1/2vTQv −cTv+µT(Av −b)es estrictamente convexa. Además, su mínimo use puede calcular de forma exacta y vendrá determinado por: Qu +ATµ−c=0⇒u=Q−1c−ATµ(2.12) De modo que, substituyendo el valor obtenido en 2.12 en la función dual 2.11, se obtiene que la función dual es, h(ν) = 1 2νTDν+µTd−1 2cTQ−1c donde D=−AQ−1ATyd=AQ−1(b−c). De modo que el problema dual vendría dado por, m´ax ν∈Rm 1 2νTDν−dTν−1 2cTQ−1c Sujeto a: νi≥0i= 1, . . . , m (QPD) cuya resolución, resulta relativamente más sencilla que la del problema original, y para la cuál, hay algoritmos específicos de resolución, que se pueden consultar en [11, págs. 196-207]. Capítulo 3 Optimización con restricciones: Métodos de resolución Hasta ahora, se ha dado respuesta a diferentes cuestiones sobre el problema de optimización con restricciones del tipo (Q). Esto permite conocer las diferentes propiedades que tienen en función de los tipos de funcionales y restricciones. Pero como el objetivo final es poder resolverlos, llegados a este punto, se estudiarán diferentes algoritmos que nos permitan obtener una solución. A continuación, se describirán métodos numéricos que dependiendo del tipo de problema, aportarán distintas visiones para su resolución final. 3.1. Métodos de penalización Un primer método, son los conocidos como métodos de penalización. Hay distintos tipos dentro de esta familia de algoritmos, pero una característica común que los define. Ésta consiste en reducir el problema con restricciones de tipo (Q), a uno sin restricciones del tipo (P)para resolverlo como se vio en el primer capítulo. Para ello, su idea principal consiste en construir un nuevo funcional objetivo de tal manera que si toma un elemento fuera del conjunto admisible, éste penalice esa elección aumentando el valor que toma el funcional. De ahí proviene el nombre de penalización, pues se penaliza la toma de un punto no factible, y dependiendo del tipo de penalización, se pueden diferenciar los distintos tipos de métodos que se presentan a continuación. 39 40 CAPÍTULO 3. OPTIMIZACIÓN CON RESTRICCIONES, MÉTODOS 3.1.1. Penalización exterior Este primer método, reside en la creación de una sucesión de subproblemas del tipo (P), de tal manera que la sucesión formada por las soluciones de los subproblemas, converja a la solución del problema (Q). Con este motivo, se crea una función P:Rn→Rdenominada función de penalización caracterizada por, (P(v)>0si v/∈S P(v)=0 si v∈S y en base a esta función, se define una sucesión de problemas (Pε)del siguiente modo: m´ın v∈RnJε(v) = J(v) + 1 εP(v) (Pε) Por lo que dada una sucesión de elementos {εk}k∈Npositivos convergentes a cero, se obtiene una sucesión u(εk)k∈Nde soluciones de (Pε). La cual, convergerá a un elemento u, siendo éste solución de (Q). Para comenzar, se verán algunas propiedades de las funciones de penalización y las soluciones u(εk). Lema 3.1. Dado 0<1/ε1<1/ε2, se tienen las siguientes propiedades: (1) Jε1u(ε1)≤Jε2u(ε2) (2) Ju(ε1)≤Ju(ε2) (3) Pu(ε1)≥Pu(ε2) Demostración. Aplicando la definición de solución del problema (Pε), se tiene que: Jε1(u(ε1))≤Jε1(u(ε2))≤Jε2(u(ε2))≤Jε2(u(ε1))(3.1) Por tanto, ya se ha obtenido (1). Ahora, de 3.1, se deduce que, 0≤Jε1(u(ε2))−Jε2(u(ε2))−hJε1(u(ε1))−Jε2(u(ε1))i=1 ε1 P(u(ε2))−1 ε2 P(u(ε2))− 1 ε1 P(u(ε1))−1 ε2 P(u(ε1))=1 ε1−1 ε2P(u(ε2))−P(u(ε1)) de lo que se obtiene (3). Y usando ahora (3) y lo visto en 3.1, se tiene, J(u(ε1))≤J(u(ε2)) + 1 ε1P(u(ε2))−P(u(ε1))< J(u(ε2)) lo que demuestra (2) y finaliza la demostración. Lema 3.2. Sea ula solución del problema (Q)considerado. Entonces para cada k∈Nse tiene que: J(u)≥Jεku(εk)≥Ju(εk) 3.1. MÉTODOS DE PENALIZACIÓN 41 Demostración. El resultado es inmediato considerando el siguiente desarrollo: J(u) = J(u) + 1 εk P(u)≥J(u(εk)) + 1 εk P(u(εk))≥J(u(εk)) Lema 3.3. Se considera el problema (Q)yu(ε)la solución al problema (Pε). Entonces, tomando δ=P(u(ε)), se tiene que u(ε)también es solución del problema: m´ın v∈RnJ(v) Sujeto a: P(v)≤δ (3.2) Demostración. Para cualquier vque satisfaga la condición P(v)≤δ, se tiene que, 0≤1 εP(u(ε))−P(v)=Jε(u(ε))−J(u(ε))−Jε(v) + J(v) = =hJε(u(ε))−Jε(v)i+J(v)−J(u(ε))≤J(v)−J(u(ε))⇒J(v)≥J(u(ε)) obteniendo así la definición de mínimo y concluyendo así la demostración. Llegados a este punto, se está en condiciones de deducir el algoritmo del método. Entonces, al considerar el problema (Q), teniendo en cuenta como se define la función de penalización, éste se puede reescribir como: m´ın v∈RnJ(v) Sujeto a: P(v) = 0 (3.3) Ahora bien, tomando δlo suficientemente pequeño, el problema 3.2 será una aproximación del problema original (Q). Así, la idea básica de este método, se basa en que según se reduce el parámetro de penalización εken cada iteración, se reduce el valor de P(u(εk)) como se vio en el lema 3.1. Por este motivo, si se fija un valor de la tolerancia δpara aproximar el problema (Q)mediante 3.2, se puede establecer el siguiente algoritmo. Algoritmo 3.4 (Penalización exterior).Se sigue el siguiente procedimiento: paso 0: Se toma ε0>0,δ > 0yk= 0. paso 1: Se calcula u(εk+1), solución de (Pε). paso 2: Si la penalización P(u(εk+1))< δ, entonces se finaliza. En caso contrario, se toma εk+1 =εk/10,k=k+ 1 y se vuelve al paso 1. 42 CAPÍTULO 3. OPTIMIZACIÓN CON RESTRICCIONES, MÉTODOS Observación 3.5.La actualización del valor ε, no tiene porque ser la marcada en el algoritmo. Ésta puede variar en función de nuestros intereses, haciendo que la sucesión {εk}k∈N decrezca con mayor o menor rapidez. Además, en la resolución de los subproblema (Pε), se tendrá que usar algún método de los vistos en el primer capítulo. Pero en ellos, se debe aportar un iterante inicial, pero en este caso, sería un buen iterante inicial la solución u(εk) calculada en la iteración anterior. Por otro lado, nótese que la función de penalización debe ser elegida de tal forma que el funcional Jε(v)tenga mínimo finito, ya que si no es así, el algoritmo diverge. Esta situación, se puede dar en casos donde el funcional objetivo decrezca según kvk → +∞y la función de penalización no crezca con la misma rapidez que decrece J(v). Ahora bien, una cuestión importante es el estudio de la convergencia del método. Para ello se introduce el siguiente teorema. Teorema 3.6. Se considera el problema (Q)con un funcional continuo. Entonces, cualquier punto límite de la sucesión u(εk)k∈Ngenerada por el algoritmo 3.4 es solución del problema original (Q). Demostración. Se supone que {u(εk0)}es una subsucesión convergente con límite u. Entonces, por la continuidad del funcional J, se tiene que: l´ım k0→+∞J(u(εk0)) = J(u)(3.4) Denotando por J∗el valor del funcional en la solución del problema (Q), y teniendo en cuenta los lemas 3.1 y 3.2, la sucesión {Jεk0(u(εk0))}es creciente, y está limitada superiormente por J∗. Entonces, l´ım k0→+∞Jεk0(u(εk0)) = q≤J∗(3.5) de manera que, extrayendo 3.4 de la ecuación 3.5, se tiene que, l´ım k0→+∞ 1 εk0 P(u(εk0)) = q−J(u)(3.6) y como P(u(εk0))≥0y1/εk0→+∞, por 3.6, implica que: l´ım k0→+∞P(u(εk0))=0 Usando ahora la continuidad de P(v), implica que P(u) = 0, por lo que ues un punto admisible para el problema. Además, por el lema 3.2, J(u(εk0))≤J∗y se tiene: J(u) = l´ım k0→+∞J(u(εk0))≤J∗ Por consiguiente, ues una solución óptima factible del problema original (Q). 3.1. MÉTODOS DE PENALIZACIÓN 43 Observación 3.7.De este teorema, se deduce que si la sucesión generada por el algoritmo 3.4 converge a un punto, éste será solución del problema general (Q). Sin embargo, esta sucesión no tiene por que ser convergente en general. Pues bien, puede darse que el mínimo de algún subproblema no sea finito, y por tanto, que la sucesión comience a divergir. Por este motivo, si se pueden asegurar hipótesis de convexidad y coercitividad sobre los subproblemas (Pεk), se aseguraría la existencia y unicidad de la sucesión {u(εk)}k∈Ngenerada por el algoritmo 3.4. Una cuestión interesante ahora, sería como crear una función de penalización P(v). En realidad, cualquier función que cumpla las propiedades que la definen, se puede usar. Sin embargo, algunas tiene mejores propiedades que otras. Una posible opción, sería usar la siguiente denominada penalización cuadrática: P(v) = m X i=1 (m´ax {0, ϕi(v)})2+ p X j=1 |φj(v)|2(3.7) Tal y como está definida esta penalización, el funcional Jε(v)será derivable y continuo. Pero al no ser continua la derivada segunda, puede conllevar problemas si se usa un método de segundo orden. Para ejemplificar visualmente lo que produce la penalización, tomando la función de penalización 3.7, se va a aplicar ahora al siguiente problema unidimensional, m´ın v∈RJ(v) = ev Sujeto a: ϕ1(v) = −v ϕ2(v) = v−1 (3.8) de lo que se obtienen los siguientes resultados: Figura 3.1: Penalización exterior 44 CAPÍTULO 3. OPTIMIZACIÓN CON RESTRICCIONES, MÉTODOS En la primera figura se muestra el funcional, las restricciones y el conjunto factible. Mientras, en la segunda imagen, se muestra el funcional objetivo J(v)y los diferentes Jεk(v)a medida que se afina el valor de εksegún se indica en la leyenda. Así pues, como se observa, los funcionales Jεk(v)aumentan su valor fuera de la región factible mientras que se mantiene el funcional original dentro del conjunto admisible. 3.1.2. Penalización interior Otro método similar al anterior, es el método de la penalización interior o de la barrera. En este nuevo caso, se vuelve a usar una función de penalización, a partir de la cuál se definirá un nuevo funcional objetivo, y a partir de él, una sucesión de subproblemas cuya sucesión de soluciones converja a la solución del problema inicial. Sin embargo, se introducen diferencias respecto del caso anterior. En primer lugar, sólo se va a considerar el problema (Qdes), y se pretende construir usa sucesión de soluciones factibles, por lo que cada u(εk)∈S. Para introducirse en el método, se comenzará por la definición de la función de penalización, que se define como P:Rn→Rde tal manera que, l´ım d(v,∂S)→0P(v) = +∞yP(v)>0∀v∈S donde d(v, ∂S)denota la distancia del punto va la frontera de S, es decir, la penalización aumenta según se aproxima a la frontera del conjunto admisible. En base a esta penalización, se define el siguiente problema sin restricciones: m´ın v∈RnJε(v) = J(v) + εP(v) (Pε) A continuación se presentan varias propiedades de esta función de penalización que son análogas a la penalización exterior. Lema 3.8. Dado 0< ε2< ε1, se tiene: (1) Jε2(u(ε2))≤Jε1(u(ε1)) (2) J(u(ε2))≤J(u(ε1)) (3) P(u(ε2))≥P(u(ε1)) Demostración. Análoga al lema 3.1 o consultar [15, pág 468] Lema 3.9. Siendo u(ε)solución del problema (Pε)y tomando δ=P(u(ε)). Entonces u(ε) es solución del problema, m´ın v∈RnJ(v) Sujeto a: P(v)≤δ (3.9) 3.1. MÉTODOS DE PENALIZACIÓN 45 Demostración. Para todo v∈Rnverificando que P(v)≤δ, se tiene, 0≤εP(u(ε))−P(v)=Jε(u(ε))−J(u(ε))−Jε(v) + J(v) = hJε(u(ε))−Jε(v)i+J(v)−J(u(ε))≤J(v)−J(u(ε))⇒J(u(ε))≤J(v) De lo que se obtiene la definición de mínimo y por tanto lo que se quería probar. Cuando el valor de δdel lema anterior es lo suficientemente grande, éste sirve para aproximar el siguiente problema: m´ın v∈RnJ(v) Sujeto a: P(v)<+∞ (3.10) Ahora bien, del modo en que se ha definido la función de penalización, la restricción del problema 3.10 es equivalente a que su solución sea factible para el problema inicial (Q). Sin embargo, dicha solución no se encontrará sobre la frontera de S, ya que en dichos puntos la penalización tiende a infinito. Además, si ε > 0es suficientemente pequeño y el δdel lema anterior es lo bastante grande, entonces u(ε)se encuentra dentro de la región factible Sy aproxima el valor de la solución del problema original. Entonces, se puede deducir el siguiente algoritmo. Algoritmo 3.10 (Penalización Interior).Se procede como sigue: paso 0: Se toma ε0>0,δ > 0yk= 0. paso 1: Se calcula u(εk+1), solución de (Pεk). paso 2: Si εkP(u(εk+1))< δ, se finaliza. En caso contrario, se toma εk+1 =εk/10,k=k+1 y se vuelve al paso 1. Observación 3.11.Como sucedía en penalización exterior, el decrecimiento de εken el paso 3, no tiene que ser la indicada, sino que puede ser cualquiera. Además, como se usarán los algoritmos del primer capítulo para resolver los subproblemas (Pεk), en este caso también se puede tomar el iterante anterior como punto de arranque, ya que éste estará próximo a la solución. Una vez visto este algoritmo, se introduce a continuación un teorema que aporta resultados sobre la convergencia de dicho método. Teorema 3.12. Sea J(v)un funcional limitado inferiormente en la región factible S. Entonces, el algoritmo de penalización interior termina de forma finita para δ > 0, y cuando no lo hace, entonces se verifica: 46 CAPÍTULO 3. OPTIMIZACIÓN CON RESTRICCIONES, MÉTODOS (1) l´ımk→∞ εkP(u(εk+1))=0. (2) l´ımk→∞ J(u(εk+1)) = ´ınfv∈Int(S)J(v). Y además, cualquier punto de acumulación de la sucesión {u(εk)}k∈Nes solución del problema (Qdes). Demostración. Sólo será necesario probar (1) y (2) cuando el algoritmo no termina de forma finita. Para ello, se toma η > 0arbitrario, por lo que existirá un elemento (que depende de η) al que se denotará por uη∈Int(S)de tal forma que, J(uη)<´ınf v∈Int(S)J(v) + η 2(3.11) entonces como {εk} → 0yuη∈Int(S), debe existir un ¯ k∈Ntal que: εkP(uη)<η 2∀k≥¯ k(3.12) Usando ahora que u(εk+1)es el mínimo del funcional Jεk(v): Jεku(εk+1)≤Jεk(uη)⇒εkPu(εk+1)≤J(uη) + εkP(uη)−Ju(εk+1)(3.13) Por lo que usando 3.11, 3.12 y 3.13, se obtiene que, εkPu(εk+1)≤J(uη) + εkP(uη)−Ju(εk+1)≤ ≤´ınf v∈Int(S)J(v) + η 2+η 2−J(u(εk+1))≤η∀k≥¯ k y como η > 0es arbitrario, se concluye (1). Ahora, para probar (2), utilizando la desigualdad anterior, se tiene, J(u(εk+1))≤J(uη) + εkP(uη)≤´ınf v∈Int(S)J(v) + η∀k≥¯ k y como η > 0es arbitrario, se concluye (2). Una vez vistos estos resultados, es natural plantear la cuestión sobre como tomar la función de penalización P(v). Con este motivo, se plantean las siguientes dos opciones que son las más generales: P(v) = − m X i=1 1 ϕi(v)yP(v) = − m X i=1 1 log(−ϕi(v)) Entonces, para mostrar el efecto de la penalización interior de un modo visual, se considera el siguiente problema, m´ın v∈RJ(v) = v4 Sujeto a: ϕ1(v) = v−1 ϕ2(v) = −v−1 (3.14) 3.1. MÉTODOS DE PENALIZACIÓN 47 del que se obtiene los siguientes resultados: Figura 3.2: Penalización interior En la primera gráfica, se muestra el funcional objetivo con las restricciones y el conjunto factible. Luego, en la segunda gráfica, se muestra de nuevo el funcional objetivo, y además, los funcionales Jε(v)para distintos valores de ε. De este modo, se observa la convergencia de los funcionales Jε(v)al funcional del problema original dentro de la región factible. Sin embargo, existe una peculiaridad, pues bien, en la frontera del conjunto, los valores de Jε(v)tienden al infinito, lo que muestra que el método sólo convergerá a un punto del interior del conjunto admisible. De este modo, se muestra visualmente como es necesario inicializar el método desde un punto factible. Además, para restricciones de igualdad, el método no serviría, dado que los funcionales Jε(v)no estaría definido en esos puntos. 3.1.3. Método del lagrangiano aumentado Una vez vistos los métodos basados en penalización de esta sección, se puede observar la necesidad de que el parámetro {εk}k∈N→0. Este hecho, es una gran desventaja desde el punto de vista numérico, debido a que según se toma εkmás pequeño, la resolución numérica de los subproblemas generados, puede resultar más compleja de lo previsto. Para solventar esta dificultad, se introduce el concepto de Lagrangiano aumentado. Éste, deriva de la función Lagrangiano introducida en el capítulo dos, y permitirá solventar ese problema que presentan los métodos de penalización. Para ello, se planteará en primer lugar el concepto de Lagrangiano aumentado para el problema con restricciones de igualdad (Qig). Definición 3.13. Se define la función Lagrangiano aumentado para el problema (Qig) 54 CAPÍTULO 3. OPTIMIZACIÓN CON RESTRICCIONES, MÉTODOS del conjunto S, independientemente del punto v∈Sescogido, DF(v, S) = ∅, de modo que no existen direcciones factibles de descenso. Sin embargo, bajo algunas condiciones, la existencia de direcciones factibles de descenso si está garantizada. A continuación, se muestra un resultado que muestra las peculiaridades de cómo debe ser el conjunto admisible. Teorema 3.22. Se considera el problema (Q)y se supone que el conjunto admisible Ses convexo y que el funcional J(v)también lo es. Entonces, dado un punto v∈S, siempre va a existir una dirección factible en v, si y sólo si, vno es un mínimo del problema (Q). Demostración. Obviamente, si ves mínimo del problema, no van a existir direcciones de descenso por el motivo de ser mínimo. Por tanto, se asume que vno es mínimo, y al mínimo, se denota por u∈S. Pero como J(v)es convexo, aplicando 1.13 se deduce, dT∇J(v)<0 siendo d=u−v. Ahora bien, como el conjunto admisible Ses convexo y u,v∈S, entonces d∈DF (v, S). De modo que se verifica que des una dirección factible de descenso. Llegados a este punto, de la forma en la que se ha presentado el método, es fácil presentar un algoritmo genérico para la resolución del problema (Q)a través de las direcciones factibles, el cual es el siguiente: Algoritmo 3.23 (Algoritmo genérico de direcciones factibles).Se siguen los siguientes pasos: paso 0: Se toma un iterante inicial u(0) ∈Syk= 0. paso 1: Si no existe una dirección factible de descenso d(k)para el punto u(k), entonces se finaliza. En caso contrario, se calcula d(k). paso 2: Se calcula el paso αk>0. paso 3: Se actualiza u(k+1) =u(k)+αkd(k),k=k+ 1 y se vuelve al paso 1. Observación 3.24.En el paso 2 del algoritmo anterior, se pide calcular el paso αasociado a una dirección dy un punto v. Sin embargo, este paso no es el mismo del que se habló en el primer capítulo para el caso sin restricciones. Pues bien, dicho valor αes un valor positivo que resuelve el problema unidimensional, m´ın v+αd∈SJ(v+αd)(3.23) por lo que se ve que se exige que el paso no estropee la factibilidad de los puntos de la sucesión creada con el algoritmo. 3.2. MÉTODOS DE DIRECCIONES FACTIBLES 55 Para ello, en la resolución del problema 3.23, se puede optar por hacer una búsqueda numérica del intervalo [0, δ]que asegura la factibilidad del paso α. De este modo, posteriormente se podría aplicar a dicho intervalo un método de búsqueda de paso mediante dicotomia (con derivadas o sin ellas) de la forma que se describieron en 1.5.1. Para ver un pseudocódigo que muestre esta inicialización del algoritmo, se puede consultar el Anexo II al final del trabajo. Sin embargo, no es recomendable aplicar un algoritmo del tipo CPP y CPG. Esto es debido a que si el valor δobtenido no cumple el CPG, el algoritmo no tiene porqué converger a un paso admisible haciendo que pueda tomar un paso excesivamente pequeño. En este sentido, se podría optar por algún tipo de modificación de dichos criterios. Por ejemplo, en [15, pág 493-496] se aporta una propuesta modificada de la regla del Armijo. 3.2.1. Método de Zoutendijk En el contexto de los métodos de direcciones factibles, se considerará en primer lugar el método de Zoutendijk. Éste, puede tener diferentes versiones de su algoritmo según se trate de un problema de optimización con funciones lineales o no lineales. Sin embargo, en este apartado, se considerará el caso genérico no lineal como se ha venido haciendo hasta el momento. Para ello, se comienza considerando el problema de optimización con restricciones de desigualdad, es decir, un problema del tipo (Qdes)para el que se desarrollará el método. En este contexto, el siguiente resultado es la base del método. Teorema 3.25. Sea vun punto factible del problema (Qdes). Suponiendo que el funcional objetivo y las restricciones activas son diferenciables y aquellas que no son activas son continuas. Si ∇J(v)Td<0y∇ϕi(v)Td<0para todo i∈I(v), entonces des una dirección factible de descenso. Demostración. Se supone que dsatisface que ∇J(v)Td<0y∇ϕi(v)Td<0,i∈I(v). Dado que para i /∈I(v),ϕi(v)<0y las funciones ϕi(v)son continuas en v, entonces debe existir un δ > 0tal que ϕi(v+td)≤0para todo t∈[0, δ]. Ahora bien, por la diferenciabilidad de las funciones ϕi(v), para i∈I(v), ϕi(v+td) = ϕi(v) + t∇ϕi(v)Td+o(tkdk) Como ∇ϕi(v)Td<0, entonces ϕi(v+td)< ϕi(v)=0para un t > 0suficientemente pequeño. Ahora bien, ϕi(v+td)≤0para i= 1, . . . , m para un tlo suficientemente pequeño (ya que si i /∈I(v)cualquier dirección es factible), y por tanto existe un t > 0tal que v+tdes factible. Por consiguiente, como por hipótesis ∇J(v)Td<0,dtambién es una dirección de descenso y entonces se concluye el resultado. 56 CAPÍTULO 3. OPTIMIZACIÓN CON RESTRICCIONES, MÉTODOS Tras lo visto, se deduce una posible forma de hallar direcciones factibles de descenso dado un punto admisible v. Pues bien, se podría plantear un problema que minimice el máximo entre ∇J(v)Tdy∇ϕi(v)Tdcon i∈I(v). Por consiguiente, se plantea el siguiente problema: m´ın (z,d)∈Rn+1 z Sujeto a: ∇J(v)Td−z≤0 ∇ϕi(v)Td−z≤0i∈I(v) −1≤dj≤1j= 1, . . . , n (3.24) A partir de este punto, se denotará por ¯z, ¯ da la solución del problema 3.24. Como (z, 0)es factible para cualquier z≥0, entonces ¯z≤0. Por tanto, como en la solución del problema 3.24 ¯z < 0, al aplicar el teorema anterior, se tiene que ¯ des una dirección factible de descenso. Teorema 3.26 (Teorema de Gordan).Dada una matriz Am×n, sólo uno de los siguientes sistemas tiene solución: (1) Ax <0,x∈Rn. (2) ATy=0, para un y∈Rm, con y≥0ey6=0. Demostración. Se puede consultar en [11, pág. 50-51]. Teorema 3.27. Dado un problema del tipo (Qdes), y uun punto admisible. Entonces, u es un punto de Fritz-John, si y sólo si, en el óptimo del problema 3.24, ¯z= 0. Demostración. El valor óptimo ¯zdel problema 3.24 es cero, si y sólo si, el sistema (∇J(u)Td<0 ∇ϕi(u)Td<0∀i∈I(u) no tiene solución. Entonces, haciendo uso del teorema 3.26, éste sistema no tiene solución, si y sólo si, existen escalares η0,ηicon i∈I(u)tales que, η0∇J(u) + X i∈I(u) ηi∇ϕi(u)=0 con η0≥0yηi≥0∀i∈I(u) donde η0oηison estrictamente mayores que 0 para algún i∈I(u). Lo cual es precisamente la condición de punto de Fritz-John. En este caso, no se puede asegurar la existencia de un punto KKT. Además, se ha visto que η0puede no ser estrictamente positivo, lo cuál, sería necesario y suficiente para que el punto de Fritz-John fuese KKT. 3.2. MÉTODOS DE DIRECCIONES FACTIBLES 57 Por otro lado, visto este último resultado, se puede deducir el método de Zoutendijk. Pues bien, dado un punto inicial factible, se puede proceder mediante un procedimiento iterativo en el que se obtenga una dirección factible de descenso mediante la resolución del problema 3.24. De lo que se deduce el siguiente algoritmo: Algoritmo 3.28 (Método de Zoutendijk).Se sigue el siguiente procedimiento iterativo: paso 0: Se toma u(0) ∈Sfactible y te toma k= 0. paso 1: Se calcula I(u(k))y se resuelve el problema 3.24 del que se obtiene una solución (zk,d(k)). En caso de que zk= 0, se finaliza y se obtiene un punto de Fritz-John, y en caso contrario, se pasa al paso 2. paso 2: Se calcula el paso αk, se toma u(k+1) =u(k)+αkd(k)y se vuelve al paso 1. Puede parecer contradictorio que para hallar una dirección factible de descenso, sea necesario resolver un nuevo problema de minimización con restricciones. Sin embargo, el problema planteado en 3.24, tiene una enorme ventaja, pues tanto el funcional objetivo como las restricciones son lineales. Por tanto, se está ante un problema de programación lineal. Así, usando este hecho, se puede obtener la solución de dicho problema mediante el uso del algoritmo del Símplex, que permite su resolución de forma sencilla. Éste algoritmo (el Símplex), no se ha introducido en el trabajo dado que se está a considerar problemas en el contexto no lineal, pero se puede consultar su deducción e implementación en el capítulo 13 de [13, pág 370] o en el capítulo 7 de [19, pág 420] (de éste último, se puede consultar su pseudocódigo en el Anexo II). En lo que se refiere a la convergencia del método, ésta no está asegurada. De hecho, se puede ver aplicando el contraejemplo de Wolfe [11, págs 381-382]. Sin embargo, éste problema se solventa de algún modo haciendo la modificación que Topkis y Veinott propusieron en 1967, y que garantiza la convergencia a un punto de Fritz-John. Esta modificación, consiste en substituir el problema 3.24 en el algoritmo 3.28, por el siguiente: m´ın (z,d)∈Rn+1 z Sujeto a: ∇J(v)Td−z≤0 Ψ0(v)Td−z≤ −Ψ(v) −1≤dj≤1j= 1, . . . , n (3.25) Nótese, que sigue estando dentro del marco de los problemas de programación lineal, por lo que esta modificación no implica una complejidad añadida en la resolución del subproblema. Esta nueva metodología de resolución, se conoce también como el algoritmo 58 CAPÍTULO 3. OPTIMIZACIÓN CON RESTRICCIONES, MÉTODOS de Topkis-Veinott en honor a sus creadores. En lo que se refiere a la convergencia, se puede aportar el siguiente resultado del que se omite su demostración debido a que ésta se basa en conceptos no introducidos en el trabajo, pero que se pueden consultar en la referencia aportada. Teorema 3.29. Se considera un problema del tipo (Qdes)donde el funcional objetivo y las restricciones son suficientemente diferenciables. Entonces, cualquier punto de acumulación de la sucesión {u(k)}k∈Ngenerada por el algoritmo 3.28 con la modificación de Topkis- Veinott, es un punto de Fritz-John. Demostración. Se puede consultar en [11, págs. 386-389]. 3.2.2. Método de los conjuntos activos Otro método de tipo direcciones factibles, es el conocido método de conjuntos activos. En este caso, se desarrollará para problemas de tipo cuadrático (vistos en 2.5), debido a que la resolución de estos problemas se usará en la sección posterior. Primeramente, partiendo de la notación descrita en 2.5, se va a denotar por ai∈Rn con i= 1, . . . , m la fila i-ésima de la matriz A, por lo que A= (a1|. . . |am)T. Y de forma análoga, ej∈Rncon j= 1, . . . , p se corresponde con la fila j-ésima de la matriz E, de modo que E= (e1|. . . |ep)T. Teniendo esto en cuenta, se va a comenzar introduciendo el siguiente resultado básico. Lema 3.30. Sea uun mínimo local del problema (QP), entonces ues un mínimo local del siguiente problema: m´ın v∈Rn 1 2vTQv −cTv Sujeto a: Ev =f aiTv=bii∈I(u) (3.26) Recíprocamente, sea uun punto factible del problema (QP)que verifica las condiciones KKT del problema 3.26, con multiplicadores de Lagrange (λ,µ), de manera que λ∈Rpy µ∈R|I(u)|tal que, µi≥0, i ∈I(u)(3.27) entonces utambién es un punto cumpliendo las condiciones KKT para el problema (QP). Demostración. Como ues solución del problema (QP), éste es factible y por tanto también es un punto admisible para el problema 3.26, pero además, como los funcionales objetivos coinciden para ambos problemas y ues solución de (QP), entonces también lo es de 3.26. 3.2. MÉTODOS DE DIRECCIONES FACTIBLES 59 Veamos ahora el segundo resultado, sea uun punto admisible para el problema (QP) que verifica las condiciones KKT para 3.26 con sus correspondientes multiplicadores de Lagrange cumpliendo 3.27, se tiene que, Qu −c+ETλ+X i∈I(u) µiai T=0 µi(aiu−bi)=0 yµi≥0∀i∈I(u) por lo que definiendo µi= 0 para i∈ {1, . . . , m}\I(u), se obtienen inmediatamente las condiciones KKT para el problema (QP). Una vez introducido este resultado, se puede detallar la idea en la que se basa el algoritmo. Centrándose en ella, partiendo de un punto inicial factible, se desarrolla un proceso iterativo, que en cada iteración, resuelve un subproblema de tipo (QPig). Entonces, sea u(0) un iterante inicial factible. A partir de aquí, en cada k-iteración, se calcula el conjunto de índices de las restricciones de desigualdad activas, es decir, Ik= I(u(k))⊂ {1, . . . , m}. Entonces, se definen los valores, d(k)=v−u(k)c(k)=c−Qu(k)g(k)=1 2(u(k))TQu(k)−cTu(k) que se utilizan para simplificar la notación al hacer el siguiente calculo, JQP (v) = JQP (u(k)+d(k)) = 1 2(d(k))TQd(k)−(c(k))Td(k)+g(k) y teniendo en cuenta que g(k)es una constante, se define el siguiente k-subproblema: m´ın d(k)∈Rn 1 2(d(k))TQd(k)−(c(k))Td(k) Sujeto a: Ed(k)=0 aid(k)= 0 i∈I(u(k)) (3.28) Si en el mínimo de dicho subproblema (al que se denotará por ˜ d(k)) es el vector nulo, entonces dicho valor tiene que cumplir las condición estacionaria de las condiciones KKT. Así, se tiene que −c(k)+ETλ+Pi∈Ikaiµi=0, lo que implica que u(k)cumple la condición estacionaria del problema original (QP). Así, si además verifica 3.27, entonces también será un punto KKT del problema 3.26, y en consecuencia del problema inicial (QP). En caso contrario, se reducirá el conjunto Ikeliminando el índice iktal que µik<0, es decir, Ik+1 =Ik\ {ik}. En este caso, cuando puede que haya varios índices de dicho modo, se podría proceder tomando el ikcorrespondiente a: µik= m´ın i∈Ik µ(k) i<0 µ(k) i(3.29) 60 CAPÍTULO 3. OPTIMIZACIÓN CON RESTRICCIONES, MÉTODOS Por otra parte, si ˜ d(k)6=0, entonces se define un nuevo iterante como u(k+1) =u(k)+ αk˜ d(k)para un αk∈[0,1] de tal manera que u(k+1) sea un punto factible, por lo que el problema en este punto es encontrar el αkque cumple dicha propiedad. Esto es posible dado que ˜ d(k)es una dirección factible. Pues bien, como ˜ d(k)cumple las restricciones del k-problema 3.28, tomando αk= 1 se cumplen las restricciones del problema (QP)que tienen presencia en el subproblema 3.28, es decir: E(u(k)+˜ d(k)) = Eu(k)+E˜ d(k)=0 aiu(k+1) =aiu(k)+ai˜ d(k)=aiu(k)≤bii∈I(u(k)) De modo que, basta constatar la propiedad para las restricciones definidas por las filas i-ésimas de Atales que i /∈I(u(k)). Pues bien, considerando por un lado el caso ai˜ d(k)≤0 para algún i /∈I(u(k)), aiu(k+1) =aiu(k)+αkai˜ d(k)≤aiu(k)≤bi por lo que este caso no influye, pero si ai˜ d(k)>0, entonces, aiu(k+1) =aiu(k)+αkai˜ d(k)≤bi⇔αk=bi−aiu(k) ai˜ d(k) de modo que, para asegurar la factibilidad de u(k+1) se toma: αk= m´ın          1,m´ın i /∈I(u(k)) ai˜ d(k)>0 bi−aiu(k) ai˜ d(k)          Es importante llegado este punto, pensar en el conjunto Ik+1 referente al nuevo iterante u(k+1). Entonces, si αk<1, se tiene que por la forma en la que se ha tomado αkexiste al menos un i /∈Iktal que, aiu(k+1) =aiu(k)+αk˜ d(k)=bi por lo que Ik+1 =Ik∪{i}. Sin embargo, si αk= 1, entonces Ik+1 =Ikya que lo anterior no sucedería. Y llegado a este punto, se está en condiciones de plantear un algoritmo que permita resolver el problema (QP). Algoritmo 3.31 (Método de los conjuntos activos).Se procede a través de los siguientes pasos: paso 0: Se aporta un u(0) ∈Sy se calcula I0,k= 0. paso 1: Se obtiene d(k)como la solución al k-subproblema 3.28. Entonces: Si d(k)6=0, se pasa al paso 2. 3.2. MÉTODOS DE DIRECCIONES FACTIBLES 61 En caso contrario, se comprueba la condición (3.27), que en caso de cumplirse, se para el algoritmo, y en caso contrario, se calcula el índice ikque no cumple la propiedad de la forma descrita en (3.29), se define Ik+1 =Ik\{ik},u(k+1) =u(k)y se pasa al paso 3. paso 2: Se calcula αky se toma u(k+1) =u(k)+αkd(k). Además, si αk= 1, entonces se pasa al paso 3, y en caso contrario, se calcula Ik+1 =I(u(k+1))y se vuelve al paso 1. paso 3: Se toma Ik+1 =Ik,k=k+ 1 y se vuelve al paso 1. Observación 3.32.Para implementar este método, es necesario aportar un punto inicial factible. Este hecho, puede resultar tedioso en algunos casos. Sin embargo, en MATLAB, existen ciertos comandos que al combinarlos aportan ciertos métodos para la obtención de un punto inicial factible. Una propuesta en este sentido, se puede consultar en el Anexo II. Finalmente, para terminar con el estudio del método, se introduce un resultado de convergencia del algoritmo planteado. Teorema 3.33. Si para todo k, los vectores ejcon j= 1, . . . , p yaicon i∈I(u(k))son linealmente independientes, entonces la sucesión generada por el algoritmo (3.31) converge a la solución del problema (QP)en un número finito de iteraciones, o bien, el problema (QP)no tiene mínimo finito. Demostración. Se puede consultar en el teorema 9.4.4 de [15, pág 434]. Sin embargo, para poder aplicar el algoritmo 3.31, es necesario poder resolver los subproblemas cuadráticos sujetos a restricciones afines de igualdad (al que se denotará por (QPig)). Para ello, partiendo de un problema genérico (QPig), se supone que rang(E) = p<n y que las filas de Eson linealmente independientes. Pues bien, si no lo son, o bien el sistema es incompatible y por tanto no hay elementos factibles, o bien hay restricciones redundantes. En el caso p=n, sólo existe un punto admisible, por lo que es trivial. En caso contrario, p<ny se toma Zuna matriz que tiene por columnas una base del ker(E). Ahora bien, dado un elemento factible u(0) cualquiera, la resolución del problema (QPig)es equivalente al problema: m´ınd∈RnJQP (u(0) +d) Sujeto a: Ed =0 (3.30) La restricción sobre d, hace que este elemento pertenezca al ker(E), por lo que debe existir un elemento y∈Rn−ptal que d=Zy. Por consiguiente, el problema 3.30 se 62 CAPÍTULO 3. OPTIMIZACIÓN CON RESTRICCIONES, MÉTODOS reescribe como, m´ın y∈Rn−pJQP (u(0) +Zy) = 1 2yTQZy−cZ Ty+JQP (u(0))(3.31) donde QZ=ZTQZ ycZ=ZTc−Qu(0). Pero dado que JQP (u(0))es una constante, la solución al problema 3.31 es la misma que al del problema cuadrático sin restricciones: m´ın y∈Rn−p 1 2yTQZy−cZ Ty(3.32) Así, una vez obtenida la solución al problema 3.32 (a la que se denotará por ˜ d), la solución al problema original (QPig), será u=u(0) +Z˜ d. Observación 3.34.En este trabajo, se ha optado por esta propuesta para resolver el problema (QPig)dado que se adecua al contexto y los contenidos introducidos con anterioridad. Sin embargo, se podría consultar [23, Cap. 3] para ver otros puntos de vista o [18, págs. 86-88] donde se aplica el Lagrangiano aumentado visto en 3.1.3 al problema (QPig). 3.3. Métodos SQP Uno de los métodos más actuales, son los SQP (del inglés sequential quadratic programming), el cuál consiste en reducir el problema (Q)a una serie de problemas cuadráticos (QP)con restricciones lineales, de tal manera que la sucesión de soluciones de esta serie de problemas converja a la solución del problema inicial. A continuación, se presentarán dos propuestas básicas y su deducción. Sin embargo, se omitirán aspectos de la convergencia debido a que ésta se desmarca del contexto del trabajo, ya que para su demostración se debe usar aspectos muy específicos de dichos métodos, pero se aportarán referencias donde se discuten dichos aspectos. 3.3.1. Método de Lagrange-Newton Un primer método SQP, es el conocido con el nombre de Lagrange-Newton que resuelve problemas del tipo (Qig). Para su deducción, recordando lo visto en la sección 2.3, una condición necesaria para que un punto u∈Rnsea un punto KKT, es la existencia de un vector λ∈Rpde tal forma que: (∇vL(u,λ) = ∇J(u) + Φ0(u)Tλ=0 ∇λL(u,λ) = Φ(u) = 0⇔ ∇L(u,λ) = 0 Para resolver este sistema no lineal, se puede partir de un punto inicial (u(0),λ(0))∈ Rn×Rp(próximo a la solución (u,λ)si es posible) para resolverlo mediante el método 3.3. MÉTODOS SQP 63 de Newton-Raphson. De modo que aplicando dicho procedimiento, para una k-iteración cualquiera, se satisface, HvL(u(k),λ(k)) Φ0(u(k))T Φ0(u(k)) 0 ! ∆u(k) ∆λ(k)!=− ∇J(u(k)) + Φ0(u(k))Tλ(k) Φ(u(k))! donde HvL(u(k),λ(k))denota la matriz hessiana del Lagrangiano respecto de la variable v, ∆u(k)=u(k+1) −u(k)y∆λ(k)=λ(k+1) −λ(k). Ahora bien, la resolución de dicho sistema, aporta un punto KKT para el siguiente problema cuadrático: m´ın d(k)∈Rn 1 2d(k)THvLu(k),λ(k)d(k)+∇vLu(k),λ(k)Td(k) Sujeto a: Φ0u(k)d(k)+ Φ u(k)=0 (3.33) Así, si se denota por ¯ d(k)la solución del problema 3.33 y a ¯ λ(k)sus multiplicadores de Lagrange asociados, se podría actualizar cada iterante como u(k+1) =u(k)+¯ d(k)y λ(k+1) =λ(k)+¯ λ(k). Sin embargo, esto es muy costoso debido al cálculo de la Hessiana en cada iteración. Por ello, se recurre a una aproximación de la Hessiana como en la sección 1.5.5. Pero en este caso, no se desea aproximar la inversa de la Hessiana, sino la Hessiana en sí. Con este objetivo, en [13, págs 536-538], se propone una actualización de la Hessiana partiendo de una matriz definida positiva, que ha resultado muy satisfactoria llevada a la práctica. En ella, denotando por Bkla aproximación de la matriz Hessiana en la k-iteración, viene dada por, Bk+1 =Bk−Bks(k)s(k)TBk s(k)TBks(k)+r(k)r(k)T s(k)Tr(k)(3.34) donde r(k)=θky(k)+ (1 −θk)Bks(k)ys(k),y(k)se definen como en 1.5.5 y θkviene dado por: θk=(1si s(k)Ty(k)≥0,2s(k)TBks(k) (0,8s(k)TBks(k))/(s(k)TBks(k)−s(k)Ty(k))si s(k)Ty(k)<0,2s(k)Bks(k) Sin embargo, la convergencia del método presentado hasta el momento, no esta asegurada. Según el razonamiento planteado, se espera que u(k+1),λ(k+1)sea un mejor aproximador de la solución (u,λ), que u(k),λ(k). Sin embargo, al igual que sucedía en 1.5.4, ésto no tiene que ser cierto a priori. Para asegurar que este echo sea así, se recurre a lo que se conoce como una función de mérito. Ésta, involucra comunmente al funcional objetivo y la cantidad de infactibilidad de las restricciones, ya que si el iterante k+1-ésimo 70 CAPÍTULO 3. OPTIMIZACIÓN CON RESTRICCIONES, MÉTODOS Anexos 71 Anexo I: Pseudocódigos de optimización sin restricciones 1J (v ) , grad_J ( v) ,u , d , a , b , t o l ←% Elementos a aportar 2phi=@( v ) J (u+v . ∗d) ; % D e fi ni ci o n de la funcion phi (v ) 3der_phi=@( v)d ’∗grad_J (u+v∗d) ; % D e fi ni c i o n de la derivada de phi ( v) 4while (abs (b−a ) > t o l ) 5i f (abs( der_phi ( alpha ) ) < 10∗eps) 6return % Se toma como s o l u c i o n alpha 7e l s e i f ( der_phi ( alpha ) > 0) 8b=alpha ; % Se modifica e l extremo i zq ui er do 9e l s e i f ( der_phi ( alpha ) < 0) 10 a=alpha ; % Se modifica e l extremo derecho del intervalo 11 end 12 alpha=(a+b) /2; 13 end Listing 1: Pseudocódigo dicotomia con derivadas 1J ( v) ,u , d , a , b , t o l ←% Elementos a aportar 2phi=@( v ) J (u+v∗dir) ; % D ef in i ci on de l a funcion phi ( v ) 3c=(a+b) /2; % Se toma e l punto medio d el i n t e r v a l o 4while (abs (b−a ) > t o l ) 5d=(a+c ) /2; e=(c+b) /2; 6phic=phi ( c ) ; % Evalucion de phi en c 7i f ( phi (d) < phic ) 8b=c ; c=d ; 9e l s e i f ( phi ( e ) < phic ) 10 a=c ; c=e ; 11 e l s e 12 a=d ; b=e ; 13 end 14 end 73 74 ANEXO I: PSEUDOCÓDIGOS DE OPTIMIZACIÓN SIN RESTRICCIONES 15 alpha=(a+b) /2; % Se toma alpha como s ol uc io n Listing 2: Pseudocódigo dicotomia sin derivadas 1J ( v) , grad_J ( v ) ,u , d , alpha , amplificador , itmax ←% Elementos a aportar 2phi=@( v ) J (u+v∗d) ; % D e fi ni ci o n de l a funcion phi . 3der_phi=@( v)d ’∗grad_J (u+v∗d) ; % D e fi ni ci o n de l a derivada de phi . 4phi0=phi (0) ; der_phi0=der_phi (0) ; % Evaluacion de phi y der_phi en 0. 5alpha_min=0; % Se toma un val o r i n i c i a l para el valor minimo de alpha . 6% Se i n i c i a la busqueda de un alpha_max que cumpla CPG: 7in t=f a l s e ; k=0; 8while (~ in t && (k <= itmax ) ) 9phi_alpha=phi ( alpha ) ; % Evaluacion de phi en alpha 10 der_phi_alpha=der_phi ( alpha ) ; % Evaluacion de l a derivada de phi en alpha 11 CPG=(phi_alpha > phi0+m1∗der_phi0∗alpha ) ; % C r i t e r i o de paso grande 12 CPP=(der_phi_alpha < m2∗der_phi0 ) && (~CPG) ; % C r i t e r i o de paso pequeno 13 k=k+1; 14 i f CPG % Se cumple e l CPG 15 alpha_max=alpha ; 16 in t=true ; 17 e l s e i f CPP % Se cumple e l CPP 18 alpha_min=alpha ; 19 alpha=alpha∗amplificador ; % Se aumenta e l alpha ( a mpl if ic ad or > 1) 20 e l s e % El paso es admisible 21 alpha ←% Se toma alpha como s ol uc io n 22 return ;% Se f i n a l i z a con e xi to 23 end 24 end Listing 3: Pseudocódigo de búsqueda del intervalo de implementación criterios de tipo CPP y CPG 1J ( v) , grad_J ( v ) ,u0 , delta , itmax ←% Elementos a aportar 2k=0; 3% Se r e a l i z a una i t e r a c i o n del metodo . 4d←% Se toma una d i r e c c i o n de descenso . 5alpha ←% Se c a l c u l a e l paso . 6u=u0+alpha ∗d ; k=k+1; 7% Se procede con un proceso i t e r a t i v o : 8while (norm(u−u0 , 2 ) /(1+norm(u , 2 ) ) > d elt a ) | | (norm( grad_J (u ) , I nf ) > d e lta ) & (k < itmax ) 9u0=u ; 10 d←% Se toma una d i r e c c i o n de descenso . 11 alpha ←% Se c a l c u l a e l paso . 12 u=u0+alpha ∗d ; k=k+1; % Se r e a l i z a n l a s a c t u a l i z a c i o n e s p e rt i ne n te s . 75 13 end Listing 4: Pseudocódigo del algoritmo genérico de resolución 1J ( v) , grad_J ( v ) ,u0 , d ,m, tol , itmax ←% Elementos a aportar 2k=0; 3n=length( u0 ) ; % Dimension del problema . 4% R ea liz ac io n de una i t e r a c c i o n del algoritmo : 5grad1=grad_J ( u0 ) ; % Evaluacion del grandiente en u0 . 6H0=eye (n) ; % Matriz i n i c i a l d e f in i d a p o s i t i v a . 7d=−H0∗grad1 ; % Obtencion de l a d i r e c c i o n . 8alpha ←% Se l la ma ri a a una func io n que c a l c u l e e l paso 9u=u0+alpha ∗d ; grad2=grad_J ( u) ; s=u−u0 ; y=grad2−grad1 ; k=k+1; 10 u0=u ; grad1=grad2 ; 11 % Se empieza e l proceso i t e r a t i v o : 12 while (( norm( grad2 , I nf ) > t o l ) | | (norm( s ( : , end ) , I nf ) > t o l ) ) 13 gamma=(s ( : , end ) ’∗y ( : , end ) ) /(y ( : , end ) ’∗y ( : , end ) ) ; 14 H0=gamma∗eye (n) ; 15 q=grad_J ( u0 ) ; a l f a = [ ] ; 16 for j=k−1:−1:k−s i z e ( s , 2 ) 17 i=j+1+s i z e ( s , 2 )−k ; 18 v=1/(y ( : , i ) ’∗s ( : , i ) ) ∗s ( : , i ) ’∗q ; 19 a l f a =[v ; a l f a ] ; 20 q=q−v∗y ( : , i ) ; 21 end 22 r=H0∗q ; 23 for j=k−s i z e ( s , 2 ) : k−1 24 i=j+1+s i z e ( s , 2 )−k ; 25 beta=1/(y ( : , i ) ’∗s ( : , i ) ) ∗y ( : , i ) ’∗r ; 26 r=r+s ( : , i ) ∗( a l f a ( i )−beta) ; 27 end 28 d=−r ; 29 % Comprobacion de la d i r e c c i o n de desc ens o : 30 i f (d ’ ∗grad1 > 0) 31 d=−d ; 32 end 33 alpha ←% Se l la ma ri a a una func io n que c a l c u l e e l paso 34 u=u0+alpha ∗d ; grad2=grad_J ( u) ; 35 i f (k > m) 36 s ( : , 1 ) = [ ] ; s =[ s u−u0 ] ; 37 y ( : , 1 ) = [ ] ; y=[y grad2−grad1 ] ; 38 e l s e 39 s =[u−u0 s ] ; 40 y=[grad2−grad1 y ] ; 41 end 76 ANEXO I: PSEUDOCÓDIGOS DE OPTIMIZACIÓN SIN RESTRICCIONES 42 k=k+1; u0=u ; grad1=grad2 ; % Act ual iza cion de i t e r a n t e s . 43 end 44 u←% Se toma u como so lu ci on obtenida . Listing 5: Pseudocódigo del método Cuasi-Newton con coste reducido B.F.G.S Anexo II: Pseudocódigos de optimización con restricciones 1u , d , amplificador , contraccion , itmax ←% Elementos a aportar 2% Se i n i c i a l i z a l a busqueda de un i n t e r v a l o en e l que a p l i c a r e l alg or itmo . 3a=0; b=1; k=0; 4i f (u+b∗d∈S) 5while (u+amplificador∗b∗d∈S) && (k < itmax ) 6b=amplificador∗b ; % Se a mp li fi ca e l i n t e r v a l o ( am pl if i ca do r > 1) 7k=k+1; 8end 9e l s e 10 b=contraccion∗b ; 11 while ~(u+contraccion∗b∗d∈S) 12 b=contraccion∗b ; % Se con trae e l i n t e r v a l o ( c on tr ac ci on ∈(0 ,1) ) 13 end 14 end Listing 6: Pseudocódigo de búsqueda del intervalo donde calcular el paso 1% Resolucion del problema de programacion l i n e a l : 2% Min fun ’ x , Sujeto a : Ax <= 0 y x>=0 3% Calculo de una s ol uc io n ba s i ca f a c t i b l e de i n i c i o : 4[ f , c]= s i z e (A) ; % Dimensiones de la matriz de r e s t r i c c i o n e s 5fun =[fun ; ze ros ( f , 1 ) ] ; % Re e scrit ura del vector que de sc ri be e l f un ci on al 6A=[A eye ( f ) ] ; % Re e scrit ura de la matriz de r e s t r i c c i o n e s 7B=eye( f ) ; xB=b ; IB=c+1: s i z e (A, 2 ) ; % Se toma una s ol uc io n basica f a c t i b l e 8xN=z ero s ( c , 1 ) ; IN=1:c ; % Se de f i n e la sol uc i o n complementaria 9x=[xN ; xB ] ; I =1: s i z e (A, 2 ) ; 10 11 % Se comienza con e l proceso i t e r a t i v o : 12 for i t e r =1:itmax 13 % Paso 1 : 14 lambda=B’ \ fun ( IB ) ; 15 cr=fun ( IN)−(lambda ’∗A( : , IN) ) ’ ; 77 78 ANEXO II: PSEUDOCÓDIGOS DE OPTIMIZACIÓN CON RESTRICCIONES 16 i f ( cr >= 0) 17 break 18 end 19 % Paso 2 : Determinar l a columna de piv o ta ci on : 20 [ crq , q]=min( cr ) ; 21 y=B\A( : , IN( q) ) ; 22 i f (y <= 0) 23 fprintf(’ El problema es i l i m i t a d o \n ’ ) ; 24 x = [ ] ; z = [ ] ; ←% Se toman s o l u c i o n e s v ac ia s 25 return ←% Se termina 26 end 27 % Paso 3 : Determinar l a f i l a de pi vo ta cio n . A n a l i s i s de Ratios : 28 w=xB./ y ; k=find( y <= 0) ; w(k )=I nf ; 29 [ theta , p]=min(w) ; 30 % Paso 4 : Pivotacion : 31 x(IN( q) )=theta ; 32 x ( IB )=x ( IB )−theta∗y ; 33 ep=z e ros ( f , 1 ) ; ep ( p) =1; 34 B=B+(A( : , IN( q) )−B( : , p ) ) ∗ep ’ ; 35 l=IB ( p) ; IB (p )=IN ( q ) ; xB=x ( IB ) ; 36 IN( q)=l ; 37 end 38 x=x ( 1 : c ) ; ←% Se extraen l a s componentes que son s o lu c io n 39 z=fun ( 1 : c ) ’∗x ; ←% Se c a l c u l a e l va lo r d el f u n c io n a l en la s o l u ci o n Listing 7: Código método Símplex 1K=E; k=f ; q=rank (K) ; % Se i n i c i a la matriz K y e l vector k con l a s r e s t r i c c i o n e s de igualdad 2% Se adieren la s r e s t r i c c i o n e s de desigualdad que son linealmen t e in dependient es con l a s que estan almacenadas en K 3for i =1: length(b) 4i f (rank ( [K;A( i , : ) ] ) ~= q ) 5K=[K;A( i , : ) ] ; 6k=[k ; b( i ) ] ; 7q=q+1; 8end 9end 10 % Se completa K con e l esp ac io ortogon al de K 11 i f (q < length( c ) ) 12 K=[K; n u l l (K) ’ ] ; 13 k=[k ; zeros(length(c)−q , 1 ) ] ; 14 end 15 u=K\k ; % Se toma como punto i n i c i a l f a c t i b l e e l que r es u lv e e l sistema Listing 8: Ejemplo búsqueda de punto inicial factible en el método de conjuntos activos Bibliografía [1] T.M. Apostol, Análisis matemático, Reverté S.A (1960) [2] S.I.Birbil, J.B.G.Frank, G.J.Still, An elemetary proof of the Fritz-John and Karush-Kuhn-Tucker conditions in nonlinear programming, The European Journal of Operational Research, (2006). Se puede consultar en: http://research.sabanciuniv.edu/177/1/3011800000548.pdf [3] D.P.Bertsekas, Constrained optimization and Lagrange multiplier methods,Academic Press (1982) [4] J.V.Tiel, Convex Analysis, An Introduction Text, John Wiley and Sons (1984) [5] E. Castillo, J. Conejo, P. Pedregal, R. García, N. Alguacil, Formulación y resolución de modelos de programación matemática en ingeniería y ciencia, Universidad de Castilla- La Mancha (2012) [6] B.Casas, G. Fiestras Janeiro, I. García Jurado, J. González Díaz, Introducción a la Teoría de Juegos, USC editora, (2012) [7] J. M. Viaño, M. Bruguera, Lecciones de Métodos Numéricos. Volumen 4: Optimización, Andavira (2013) [8] S.G.Nash, A.Sofer, Linear and Nonlinear Programming, McGraw-Hill (1996) [9] K.Sydsaeter, P.Hammond, Matemáticas para el análisis económico, Pretince Hall (1996) [10] D.P.Bertsekas, Nonlinear Programming, Athena Scientific (1995) [11] M.S. Bazaraa, C.M. Shetty, Nonlinear Programming. Theory and algorithms, John Whiley and sons (1979) [12] G. Allaire, A. Craig, Numerical Analysis and Optimization. An introduction to Mathematical Modelling and Numerical Simulation, Oxford University press (2007) 79