Algoritmos de optimización no lineal con restricciones
Definición de la optimización con restricciones
La minimización con restricciones es el problema de encontrar un vector x que es un mínimo local de una función escalar f(x) sujeta a restricciones sobre la x admisible:
de tal manera que se cumpla una o más de las siguientes condiciones: c(x) ≤ 0, ceq(x) = 0, A·x ≤ b, Aeq·x = beq, l ≤ x ≤ u. En la programación semiinfinita se utilizan aún más restricciones; consulte Algoritmo y formulación de problemas fseminf.
Algoritmo fmincon trust-region-reflective
Métodos trust-region para la minimización no lineal
Muchos de los métodos utilizados en los solvers de Optimization Toolbox™ se basan en regiones de confianza, un concepto sencillo, pero potente de optimización.
Para comprender el enfoque trust-region de la optimización, considere el problema de minimización no restringida, minimice f(x), donde la función toma argumentos de vector y devuelve escalares. Suponga que está en un punto x en el espacio n y desea mejorar, es decir, moverse a un punto con un valor de función inferior. La idea básica es aproximar f con una función más sencilla q, que refleja razonablemente el comportamiento de la función f en un entorno N alrededor del punto x. Este entorno es la región de confianza. Se calcula un paso de prueba s minimizando (o minimizando aproximadamente) en N. Este es el subproblema trust-region,
(1) |
El punto actual se actualiza para ser x + s si f(x + s) < f(x); de lo contrario, el punto actual permanece sin cambiar y N, la región de confianza, disminuye y se repite el cálculo del paso de prueba.
Las principales cuestiones al definir un enfoque específico trust-region para minimizar f(x) son cómo elegir y calcular la aproximación q (definida en el punto actual x), cómo elegir y modificar la región de confianza N y cómo resolver de forma precisa el subproblema trust-region. Esta sección se centra en el problema no restringido. En secciones posteriores, se abordan otras complicaciones causadas por la presencia de restricciones en las variables.
En el método estándar trust-region ([48]), la aproximación cuadrática q se define por los dos primeros términos de la aproximación de Taylor a F en x; el entorno N suele tener forma de esfera o de elipsoide. En términos matemáticos, el subproblema trust-region se indica, por lo general, como
(2) |
donde g es el gradiente de f en el punto actual x, H es la matriz hessiana (la matriz simétrica de segundas derivadas), D es una matriz de escalado diagonal, Δ es un escalar positivo y ‖ . ‖ es la norma euclídea. Existen buenos algoritmos para resolver Ecuación 2 (consulte [48]); estos algoritmos incluyen, por lo general, el cálculo de todos los valores propios de H y un proceso de Newton aplicado a la ecuación secular
Estos algoritmos proporcionan una solución precisa para Ecuación 2. Sin embargo, requieren tiempo proporcional a varias factorizaciones de H. Por lo tanto, para los problemas a gran escala se necesita un enfoque diferente. Varias estrategias heurísticas y de aproximación, basadas en Ecuación 2, se han propuesto en la literatura ([42] y [50]). El enfoque de aproximación seguido en los solvers de Optimization Toolbox es restringir el subproblema trust-region a un subespacio bidimensional S ([39] y [42]). Una vez calculado el subespacio S, el trabajo para resolver Ecuación 2 es trivial incluso si se necesita información completa del valor propio/vector propio (ya que, en el subespacio, el problema solo es bidimensional). El trabajo dominante se ha desplazado ahora a la determinación del subespacio.
El subespacio bidimensional S se determina con la ayuda de un proceso de gradiente conjugado precondicionado descrito a continuación. El solver define S como el espacio lineal que abarcan s1 y s2, donde s1 está en la dirección del gradiente g y s2 es una dirección de Newton aproximada, es decir, una solución para
(3) |
o una dirección de curvatura negativa,
(4) |
La filosofía sobre la que se basa esta elección de S es forzar la convergencia global (mediante la dirección de descenso más pronunciado o la dirección de curvatura negativa) y lograr una convergencia local rápida (mediante el paso de Newton, cuando exista).
Ahora es fácil dar una aproximación de minimización no restringida utilizando ideas trust-region:
Formule el subproblema trust-region bidimensional.
Resuelva Ecuación 2 para determinar el paso de prueba s.
Si f(x + s) < f(x), entonces x = x + s.
Ajuste Δ.
Estos cuatro pasos se repiten hasta la convergencia. La dimensión Δ de la región de confianza se ajusta de acuerdo con las reglas estándar. En particular, disminuye si el paso de prueba no se acepta, es decir, f(x + s) ≥ f(x). Consulte [46] y [49] para ver una exposición sobre este aspecto.
Los solvers de Optimization Toolbox tratan unos pocos casos importantes especiales de f con funciones especializadas: mínimos cuadrados no lineales, funciones cuadráticas y mínimos cuadrados lineales. Sin embargo, las ideas algorítmicas subyacentes son las mismas que para el caso general. Estos casos especiales se abordarán en secciones posteriores.
Método del gradiente conjugado precondicionado
Una forma popular de resolver sistemas definidos grandes, simétricos y positivos de ecuaciones lineales Hp = –g es el método de los gradientes conjugados precondicionados (PCG). Este enfoque iterativo requiere la capacidad de calcular productos de matriz-vector de la forma H·v, donde v es un vector arbitrario. La matriz definida simétrica positiva M es un precondicionador para H. Es decir, M = C2, donde C–1HC–1 es una matriz bien condicionada o una matriz con valores propios agrupados.
En un contexto de minimización, puede dar por supuesto que la matriz hessiana H es simétrica. Sin embargo, está garantizado que H solo sea definida positiva en el entorno de un minimizador fuerte. El algoritmo PCG existe cuando encuentra una dirección de curvatura negativa (o cero), es decir, dTHd ≤ 0. La dirección de salida PCG p es una dirección de curvatura negativa o una solución aproximada al sistema de Newton Hp = –g. En cualquier caso, p ayuda a definir el subespacio bidimensional utilizado en el enfoque trust-region abordado en Métodos trust-region para la minimización no lineal.
Restricciones de igualdad lineales
Las restricciones lineales complican la situación descrita para la minimización sin restricciones. Sin embargo, las ideas subyacentes descritas anteriormente pueden llevarse a cabo de forma limpia y eficaz. Los métodos trust-region en solvers de Optimization Toolbox generan iterados estrictamente factibles.
El problema general de minimización lineal con restricciones de igualdad se puede escribir como
(5) |
donde A es una matriz de m por n (m ≤ n). Algunos solvers de Optimization Toolbox preprocesan A para eliminar las dependencias lineales estrictas mediante una técnica basada en la factorización LU de AT. Aquí A se supone que es de rango m.
El método usado para resolver Ecuación 5 difiere del enfoque sin restricciones en dos aspectos significativos. En primer lugar, se calcula un punto factible inicial x0, utilizando un paso de mínimos cuadrados dispersos, de modo que Ax0 = b. En segundo lugar, el algoritmo PCG se sustituye por gradientes conjugados precondicionados reducidos (RPCG, por sus siglas en inglés), consulte [46], para calcular un paso de Newton reducido aproximado (o una dirección de curvatura negativa en el espacio nulo de A). El paso clave del álgebra lineal consiste en resolver sistemas de la forma
(6) |
donde se aproxima a A (los pequeños elementos distintos de cero de A se establecen en cero siempre que no se pierda rango) y C es una aproximación positiva, definida y simétrica dispersa a H, es decir, C = H. Consulte [46] para ver más detalles.
Restricciones de caja
El problema de caja restringida tiene la forma
(7) |
donde l es un vector de límites inferiores y u es un vector de límites superiores. Algunos (o todos) los componentes de l pueden ser iguales a -∞ y algunos (o todos) los componentes de u pueden ser iguales a ∞. El método genera una secuencia de puntos estrictamente factibles. Se emplean dos técnicas para mantener la factibilidad a la vez que se consigue un comportamiento de convergencia robusto. En primer lugar, un paso de Newton modificado a escala sustituye al paso de Newton sin restricciones (para definir el subespacio bidimensional S). En segundo lugar, se utilizan reflexiones para aumentar el tamaño del paso.
El paso de Newton modificado a escala surge del examen de las condiciones necesarias de Kuhn-Tucker para Ecuación 7,
(8) |
donde
y el vector v(x) se define a continuación, para cada 1 ≤ i ≤ n:
Si gi < 0 y ui < ∞, entonces vi = xi – ui
Si gi ≥ 0 y li > –∞, entonces vi = xi – li
Si gi < 0 y ui = ∞, entonces vi = –1
Si gi ≥ 0 y li = –∞, entonces vi = 1
El sistema no lineal Ecuación 8 no puede diferenciarse en todas partes. La no diferenciabilidad se produce cuando vi = 0. Se pueden evitar estos puntos manteniendo una factibilidad estricta, es decir, restringiendo l < x < u.
El paso de Newton modificado a escala sk para el sistema de ecuaciones no lineal dado por Ecuación 8 se define como la solución del sistema lineal
(9) |
en la k-ésima iteración, donde
(10) |
y
(11) |
Aquí Jv desempeña el papel de la matriz jacobiana de |v|. Cada componente diagonal de la matriz diagonal Jv es igual a 0, -1 o 1. Si todos los componentes de l y u son finitos, Jv = diag(sign(g)). En un punto en el que gi = 0, vi podría no ser diferenciable. se define en un punto de este tipo. La no diferenciabilidad de este tipo no es motivo de preocupación, porque, en el caso de un componente de este tipo, no es significativo el valor que tome vi. Además, |vi| seguirá siendo discontinua en este punto, pero la función |vi|·gi es continua.
En segundo lugar, se utilizan reflexiones para aumentar el tamaño del paso. Se define un (único) paso de reflexión del siguiente modo. Dado un paso p que cruza una restricción límite, considere la primera restricción límite cruzada por p; suponga que es la i-ésima restricción límite (ya sea la i-ésima superior o la i-ésima inferior). Entonces, el paso de reflexión pR = p, excepto en el i-ésimo componente, donde pRi = –pi.
Algoritmo fmincon Active-Set
Introducción
En la optimización con restricciones, el objetivo general es transformar el problema en un subproblema más sencillo que pueda resolverse y utilizarse como base de un proceso iterativo. Una característica de los primeros métodos es trasladar el problema con restricciones a un problema básico sin restricciones mediante el uso de una función de penalización con las restricciones que están cerca o más allá del límite de la restricción. De este modo, el problema con restricciones se resuelve mediante una secuencia de optimizaciones sin restricciones parametrizadas, que en el límite (de la secuencia) convergen al problema con restricciones. Actualmente, estos métodos se consideran relativamente ineficaces y se han sustituido por otros centrados en la solución de las ecuaciones de Karush-Kuhn-Tucker (KKT). Las ecuaciones de KKT son condiciones necesarias para la optimalidad de un problema de optimización con restricciones. Si el problema es de los llamados de programación convexa, es decir, si f(x) y Gi(x), i = 1,...,m son funciones convexas, las ecuaciones de KKT son necesarias y suficientes para un punto de solución global.
En lo que respecta a GP (Ecuación 1), las ecuaciones de Kuhn-Tucker pueden enunciarse como
(12) |
además de las restricciones originales en Ecuación 1.
La primera ecuación describe una cancelación de los gradientes entre la función objetivo y las restricciones activas en el punto de solución. Para que los gradientes se cancelen, son necesarios multiplicadores de Lagrange (λi, i = 1,...,m) para equilibrar las desviaciones en magnitud de los gradientes de las restricciones y de la función objetivo. Dado que solo las restricciones activas se incluyen en esta operación de cancelación, las restricciones que no están activas no deben incluirse en esta operación, por lo que reciben multiplicadores de Lagrange iguales a 0. Esto está implícito en las dos últimas ecuaciones de Kuhn-Tucker.
La solución de las ecuaciones de KKT constituye la base de muchos algoritmos de programación no lineal. Estos algoritmos intentan calcular directamente los multiplicadores de Lagrange. Los métodos cuasi-Newton con restricciones garantizan una convergencia superlineal mediante la acumulación de información de segundo orden relativa a las ecuaciones de KKT utilizando un procedimiento de actualización cuasi-Newton. Estos métodos suelen denominarse métodos de programación cuadrática secuencial (SQP, por sus siglas en inglés), ya que en cada iteración principal se resuelve un subproblema de QP (también conocidos como métodos de programación cuadrática iterativa, programación cuadrática recursiva y métrica variable con restricciones).
El algoritmo 'active-set'
no es un algoritmo a gran escala; consulte Algoritmos a gran escala frente a algoritmos a media escala.
Programación cuadrática secuencial (SQP)
Los métodos SQP son los más avanzados en programación no lineal. Schittkowski [36], por ejemplo, ha implementado y probado una versión que supera a todos los demás métodos probados en términos de eficacia, precisión y porcentaje de soluciones acertadas en un gran número de problemas de prueba.
Este método, basado en los trabajos de Biggs [1], Han [22] y Powell ([32] y [33]), permite imitar fielmente el método de Newton en la optimización con restricciones, al igual que se hace en la optimización sin restricciones. En cada iteración principal, se realiza una aproximación de la matriz hessiana de la función lagrangiana mediante un método de actualización cuasi-Newton. Luego esto se utiliza para generar un subproblema de QP cuya solución se emplea para formar una dirección de búsqueda en un procedimiento de búsqueda de recta. Se puede encontrar una visión general de SQP en los trabajos de Fletcher [13] Gill et al. [19], Powell [35] y Schittkowski [23]. No obstante, a continuación se expone el método general.
Dada la descripción del problema en GP (Ecuación 1), la idea principal es la formulación de un subproblema de QP basado en una aproximación cuadrática de la función lagrangiana.
(13) |
Aquí se simplifica Ecuación 1 suponiendo que las restricciones de límite se han expresado como restricciones de desigualdad. Se obtiene el subproblema de QP linealizando las restricciones no lineales.
Subproblema de programación cuadrática (QP)
(14) |
Este subproblema puede resolverse mediante cualquier algoritmo de QP (consulte, por ejemplo, Solución de programación cuadrática). La solución se utiliza para formar una nueva iteración
xk + 1 = xk + αkdk.
El parámetro de longitud de paso αk se determina mediante un procedimiento adecuado de búsqueda de recta, de modo que se obtenga una disminución suficiente de una función de mérito (véase Actualización de la matriz hessiana). La matriz Hk es una aproximación definida positiva de la matriz hessiana de la función lagrangiana (Ecuación 13). Hk puede actualizarse mediante cualquiera de los métodos cuasi-Newton, aunque el método BFGS (consulte Actualización de la matriz hessiana) parece ser el más popular.
Un problema con restricciones no lineales a menudo puede resolverse en menos iteraciones que un problema sin restricciones mediante SQP. Una de las razones de ello es que, gracias a los límites del área factible, el optimizador puede tomar decisiones informadas sobre las direcciones de búsqueda y la longitud de los pasos.
Consideremos la función de Rosenbrock con una restricción de desigualdad no lineal adicional, g(x),
(15) |
Se resolvió mediante una implementación de SQP en 96 iteraciones, frente a las 140 del caso sin restricciones. Figura 5-3, Método SQP para la función de Rosenbrock con restricciones no lineales muestra el camino hacia el punto de solución x = [0.9072,0.8228] a partir de x = [–1.9,2.0].
Figura 5-3, Método SQP para la función de Rosenbrock con restricciones no lineales
Implementación de SQP
La implementación de SQP consta de tres etapas principales, las cuales se comentan brevemente en los siguientes subapartados:
Actualización de la matriz hessiana. En cada iteración principal se calcula una aproximación cuasi-Newton definida y positiva de la hessiana de la función lagrangiana, H, utilizando el método BFGS, donde λi, i = 1,...,m es una estimación de los multiplicadores de Lagrange.
(16) |
donde
Powell [33] recomienda mantener la hessiana definida positiva aunque pueda ser indefinida positiva en el punto de solución. Se mantiene una hessiana definida positiva siempre que sea positivo en cada actualización y que H se inicialice con una matriz definida positiva. Cuando no es positivo, qk se modifica elemento a elemento para que . El objetivo general de esta modificación es distorsionar lo menos posible los elementos de qk, lo que contribuye a que la actualización sea definida positiva. Por tanto, en la fase inicial de la modificación, el elemento más negativo de qk*sk se reduce repetidamente a la mitad. Este procedimiento se continúa hasta que sea mayor o igual que una pequeña tolerancia negativa. Si, después de este procedimiento, sigue sin ser positivo, modifique qk añadiendo un vector v multiplicado por un escalar constante w, es decir,
(17) |
donde
y aumente w sistemáticamente hasta que se vuelva positivo.
Las funciones fmincon
, fminimax
, fgoalattain
y fseminf
usan todas SQP. Si Display
se establece en 'iter'
en options
, se proporciona información diversa, como los valores de la función y la vulneración máxima de restricciones. Cuando hay que modificar la hessiana utilizando la primera fase del procedimiento anterior para mantenerla definida y positiva, se muestra Hessian modified
. Si hay que volver a modificar la hessiana mediante la segunda fase del enfoque descrito anteriormente, se muestra Hessian modified twice
. Cuando el subproblema de QP no es factible, se muestra infeasible
. Este tipo de visualizaciones no son preocupantes, pero indican que el problema es muy poco lineal y que la convergencia puede tardar más de lo normal. A veces aparece el mensaje no update
, que indica que es casi cero. Esto puede indicar que la configuración del problema es incorrecta o que se está intentando minimizar una función no continua.
Solución de programación cuadrática. En cada iteración principal del método SQP, se resuelve un problema de QP de la forma descrita más abajo, donde Ai se refiere a la i
-ésima fila de la matriz A de m por n.
(18) |
El método utilizado en las funciones de Optimization Toolbox es una estrategia de conjuntos activos (también conocida como método de proyección) similar a la de Gill et al., descrita en [18] y [17]. Se ha modificado tanto para problemas de programación lineal (LP) como de programación cuadrática (QP).
El procedimiento de solución consta de dos fases. La primera fase consiste en calcular un punto factible (si existe). La segunda fase consiste en generar una secuencia iterativa de puntos factibles que convergen a la solución. En este método se mantiene un conjunto activo, , que es una estimación de las restricciones activas (es decir, las que están en los límites de las restricciones) en el punto de solución. Prácticamente todos los algoritmos de QP son métodos de conjuntos activos. Se hace hincapié en este punto porque existen muchos métodos diferentes que son muy similares en estructura pero que se describen en términos muy diferentes.
se actualiza en cada iteración k y esto se utiliza como base para una dirección de búsqueda . Las restricciones de igualdad siempre permanecen en el conjunto activo . La notación de la variable se utiliza aquí para distinguirla de dk en las iteraciones principales del método SQP. Se calcula la dirección de búsqueda y minimiza la función objetivo permaneciendo en cualquier límite de restricción activa. El subespacio factible para se forma a partir de la base Zk cuyas columnas son ortogonales a la estimación del conjunto activo (es decir, ). Así, se garantiza que una dirección de búsqueda que se forma a partir de una suma lineal de cualquier combinación de las columnas de Zk permanezca en los límites de las restricciones activas.
La matriz Zk se forma a partir de las últimas columnas m – l de la descomposición de QR de la matriz , donde l es el número de restricciones activas y l < m. Es decir, Zk viene dado por
(19) |
donde
Una vez encontrado Zk, se busca una nueva dirección de búsqueda que minimice q(d) donde está en el espacio nulo de las restricciones activas. Es decir, es una combinación lineal de las columnas de Zk: para algún vector p.
Entonces, si vemos la cuadrática como función de p, sustituyendo por , se obtiene
(20) |
Si se diferencia con respecto a p se obtiene
(21) |
∇q(p) se denomina gradiente proyectado de la función cuadrática porque es el gradiente proyectado en el subespacio definido por Zk. El término se denomina hessiana proyectada. Si la matriz hessiana H es definida y positiva (como es el caso en esta implementación de SQP), el mínimo de la función q(p) en el subespacio definido por Zk ocurre cuando ∇q(p) = 0, que es la solución del sistema de ecuaciones lineales
(22) |
A continuación, se da un paso con el formato
(23) |
Debido a la naturaleza cuadrática de la función objetivo, en cada iteración solo hay dos opciones de longitud de paso α. Un paso de unidad a lo largo de es el paso exacto para llegar al mínimo de la función restringida al espacio nulo de . Si ese paso se puede dar sin vulnerar las restricciones, esa es la solución a QP (Ecuación 18). En caso contrario, el paso que va desde hasta la restricción más cercana es menor que la unidad y se incluye una nueva restricción en el conjunto activo en la siguiente iteración. La distancia a los límites de la restricción en cualquier dirección viene dada por
(24) |
que se define en el caso de las restricciones que no están en el conjunto activo y donde la dirección es hacia el límite de la restricción, es decir, .
Cuando se incluyen n restricciones independientes en el conjunto activo, sin localización del mínimo, se calculan los multiplicadores de Lagrange, λk, que satisfacen el conjunto no singular de ecuaciones lineales
(25) |
Si todos los elementos de λk son positivos, xk es la solución óptima de QP (Ecuación 18). Sin embargo, si algún componente de λk es negativo y dicho componente no corresponde a una restricción de igualdad, el elemento correspondiente se elimina del conjunto activo y se busca un nuevo iterado.
Inicialización. El algoritmo requiere un punto factible para empezar. Si el punto del método SQP no es factible, puede encontrar uno resolviendo el problema de programación lineal
(26) |
La notación Ai indica la i-ésima fila de la matriz A. Puede encontrar un punto factible (si existe) para Ecuación 26 estableciendo x en un valor que satisfaga las restricciones de igualdad. Puede determinar este valor resolviendo un conjunto subdeterminado o sobredeterminado de ecuaciones lineales formadas a partir del conjunto de restricciones de igualdad. Si hay una solución a este problema, la variable de holgura γ se establece en la restricción de desigualdad máxima en este punto.
Puede modificar el algoritmo de QP anterior para problemas de LP estableciendo la dirección de búsqueda en la dirección de descenso más pronunciada en cada iteración, donde gk es el gradiente de la función objetivo (igual a los coeficientes de la función objetivo lineal).
(27) |
Si se encuentra un punto factible utilizando el método de LP anterior, se entra en la fase principal de QP. La dirección de búsqueda se inicializa con una dirección de búsqueda hallada a partir de la resolución del conjunto de ecuaciones lineales
(28) |
donde gk es el gradiente de la función objetivo en el iterado actual xk (es decir, Hxk + c).
Si no se encuentra una solución factible para el problema de QP, la dirección de búsqueda de la rutina principal SQP se toma como aquella que minimiza γ.
Búsqueda de recta y función de mérito. La solución del subproblema de QP produce un vector dk, que se utiliza para formar un nuevo iterado
(29) |
El parámetro de longitud de paso αk se determina para producir un descenso suficiente en una función de mérito. En esta implementación se utiliza la función de mérito empleada por Han [22] y Powell [33] de la siguiente forma.
(30) |
Powell recomienda fijar el parámetro de penalización
(31) |
Esto permite una contribución positiva de las restricciones que están inactivas en la solución de QP pero que estaban activas hace poco. En esta implementación, el parámetro de penalización ri se fija inicialmente en
(32) |
donde representa la norma euclidiana.
De este modo, se garantiza una mayor contribución al parámetro de penalización de las restricciones con gradientes más pequeños, como sería el caso de las restricciones activas en el punto de solución.
Algoritmo SQP de fmincon
El algoritmo sqp
(y el algoritmo casi idéntico sqp-legacy
) es similar al algoritmo active-set
(para obtener una descripción, consulte Algoritmo fmincon Active-Set). El algoritmo básico sqp
se describe en el capítulo 18 de Nocedal y Wright [31].
El algoritmo sqp
es esencialmente el mismo que el algoritmo sqp-legacy
, pero tiene una implementación diferente. Normalmente, sqp
tiene un tiempo de ejecución más rápido y un menor uso de memoria que sqp-legacy
.
Las diferencias más importantes entre los algoritmos sqp
y active-set
son:
Factibilidad estricta con respecto a los límites
El algoritmo sqp
realiza cada paso iterativo en la región restringida por límites. Además, los pasos de diferencias finitas también respetan los límites. Los límites no son estrictos; un paso puede estar exactamente en un límite. Esta factibilidad estricta puede ser beneficiosa cuando la función objetivo o las funciones con restricciones no lineales no están definidas o son complejas fuera de la región restringida por los límites.
Solidez de los resultados no dobles
Durante sus iteraciones, el algoritmo sqp
puede intentar dar un paso que falle. Esto significa que una función objetivo o una función de restricción no lineal que usted proporcione devuelve un valor de Inf
, NaN
o un valor complejo. En este caso, el algoritmo intenta dar un paso más pequeño.
Rutinas de álgebra lineal refactorizadas
El algoritmo sqp
utiliza un conjunto diferente de rutinas de álgebra lineal para resolver el subproblema de programación cuadrática, Ecuación 14. Estas rutinas son más eficientes tanto en uso de memoria como en velocidad que las rutinas active-set
.
Rutinas de factibilidad reformuladas
El algoritmo sqp
tiene dos nuevos enfoques para la solución de Ecuación 14 cuando no se satisfacen las restricciones.
El algoritmo
sqp
combina las funciones objetivo y de restricción en una función de mérito. El algoritmo intenta minimizar la función de mérito sujeta a restricciones relajadas. Este problema modificado puede conducir a una solución factible. Sin embargo, este enfoque tiene más variables que el problema original, por lo que el tamaño del problema en Ecuación 14 aumenta. El aumento de tamaño puede ralentizar la solución del subproblema. Estas rutinas se basan en los artículos de Spellucci [60] y Tone [61]. El algoritmosqp
establece el parámetro de penalización de la función de mérito Ecuación 30 según la sugerencia de [41].Supongamos que no se satisfacen las restricciones no lineales y que un intento de paso hace que la vulneración de la restricción aumente. El algoritmo
sqp
intenta obtener la factibilidad utilizando una aproximación de segundo orden a las restricciones. La técnica de segundo orden puede conducir a una solución factible. Sin embargo, esta técnica puede ralentizar la solución, ya que requiere más evaluaciones de las funciones de restricción no lineales.
Algoritmo interior-point de fmincon
Función de barrera
El enfoque de interior-point para la minimización con restricciones consiste en resolver una secuencia de problemas de minimización aproximada. El problema original es
(33) |
(34) |
El problema aproximado Ecuación 34 es una secuencia de problemas con restricciones de igualdad. Son más fáciles de resolver que el problema original con restricciones de desigualdad Ecuación 33.
Para resolver el problema aproximado, el algoritmo utiliza uno de los dos tipos principales de pasos en cada iteración:
Un paso directo en (x, s). Este paso intenta resolver las ecuaciones de KKT, Ecuación 2 y Ecuación 3 para el problema aproximado a través de una aproximación lineal. También se denomina paso de Newton.
Un paso de GC (gradiente conjugado), utilizando una región de confianza.
De forma predeterminada, el algoritmo intenta primero dar un paso directo. Si no puede, intenta dar un paso de GC. Un caso en el que no da un paso directo es cuando el problema aproximado no es localmente convexo cerca del iterado actual.
En cada iteración, el algoritmo disminuye una función de mérito, como
(35) |
Si la función objetivo o una función de restricción no lineal devuelve un valor complejo, NaN, Inf o un error en un iterado xj, el algoritmo rechaza xj. El rechazo tiene el mismo efecto que si la función de mérito no disminuyera lo suficiente: el algoritmo intenta entonces un paso diferente, más corto. Envuelva cualquier código que pueda dar error en try
-catch
:
function val = userFcn(x) try val = ... % code that can error catch val = NaN; end
El objetivo y las restricciones deben dar valores adecuados (double
) en el punto inicial.
Paso directo
Para definir el paso directo se utilizan las siguientes variables:
H denota la hessiana del lagrangiano de fμ:
(36) Jg denota la jacobiana de la función de restricción g.
Jh denota la jacobiana de la función de restricción h.
S = diag(s).
λ denota el vector multiplicador de Lagrange asociado a las restricciones de g
Λ = diag(λ).
y denota el vector multiplicador de Lagrange asociado a h.
e denota el vector que tiene el mismo tamaño que g.
Ecuación 38 define el paso directo (Δx, Δs):
(37) |
Esta ecuación proviene directamente de intentar resolver Ecuación 2 y Ecuación 3 utilizando una lagrangiana linealizada.
Puede simetrizar la ecuación premultiplicando la segunda variable Δs por S-1:
(38) |
Para resolver esta ecuación para (Δx, Δs), el algoritmo realiza una factorización LDL de la matriz. (Consulte Ejemplo 3: la estructura de D en la página de referencia de la función ldl
de MATLAB®). Este es el paso más exigente desde el punto de vista computacional. Uno de los resultados de esta factorización es la determinación de si la hessiana proyectada es definida positiva o no; si no lo es, el algoritmo utiliza un paso de gradiente conjugado, descrito en Paso de gradiente conjugado.
Actualizar un parámetro de barrera
Para que el problema aproximado Ecuación 34 se aproxime al problema original, el parámetro de barrera μ debe disminuir hacia 0 a medida que avanzan las iteraciones. El algoritmo dispone de dos opciones de actualización de los parámetros de barrera, que se especifican mediante la opción BarrierParamUpdate
: 'monotone'
(de forma predeterminada) y 'predictor-corrector'
.
La opción 'monotone'
disminuye el parámetro μ mediante un factor de 1/100 o 1/5 cuando el problema aproximado se resuelve con suficiente precisión en la iteración anterior. La opción utiliza un factor de 1/100 cuando el algoritmo tarda solo una o dos iteraciones en alcanzar una precisión suficiente, y usa 1/5 en caso contrario. La medida de la precisión es la siguiente prueba, que determina si el tamaño de todos los términos del lado derecho de Ecuación 38 es menor que μ:
Nota
fmincon
sustituye el ajuste BarrierParamUpdate
por 'monotone'
en cualquiera de estos casos:
El problema no tiene restricciones de desigualdad, incluidas las restricciones de límite.
La opción
SubproblemAlgorithm
es'cg'
.
El algoritmo 'predictor-corrector'
para actualizar el parámetro de barrera μ es similar al algoritmo de programación lineal Predictor-Corrector.
Los pasos predictores-correctores pueden acelerar el enfoque existente de Fiacco-McCormick (monótono) ajustando el error de linealización en los pasos de Newton. Los efectos del algoritmo predictor-corrector son dobles: a menudo mejora las direcciones de los pasos y, al mismo tiempo, actualiza el parámetro de barrera de forma adaptativa con el parámetro de centrado σ para hacer que los iterados sigan el camino central. Consulte el análisis de Nocedal y Wright [31] sobre los pasos predictores-correctores para programas lineales para entender por qué el camino central permite tamaños de paso mayores y, en consecuencia, una convergencia más rápida.
El paso predictor utiliza el paso linealizado con μ = 0, es decir, sin función de barrera:
Defina ɑs y ɑλ como los tamaños de paso más grandes que no vulneren las restricciones de no negatividad.
Ahora calcule la complementariedad a partir del paso predictor.
(39) |
donde m es el número de restricciones.
El primer paso corrector ajusta el término cuadrático omitido en la linealización de búsqueda de raíces de Newton
Para corregir el error cuadrático, resuelva el sistema lineal para la dirección del paso corrector.
El segundo paso corrector es un paso de centrado. La corrección de centrado se basa en la variable σ de la parte derecha de la ecuación
Aquí, σ se define como
donde μP se define en la ecuación Ecuación 39 y .
Para evitar que el parámetro de barrera disminuya demasiado rápido, lo que podría desestabilizar el algoritmo, este mantiene el parámetro de centrado σ por encima de 1/100. Esta acción hace que el parámetro de barrera μ disminuya como máximo en un factor de 1/100.
Desde el punto de vista algorítmico, los primeros pasos de corrección y centrado son independientes entre sí, por lo que se calculan juntos. Además, la matriz de la izquierda del predictor y de los dos pasos de los correctores es la misma. Por lo tanto, desde el punto de vista algorítmico, la matriz se factoriza una vez, y esta factorización se utiliza en todos estos pasos.
El algoritmo puede rechazar el paso predictor-corrector propuesto cuando el paso aumenta el valor de la función de mérito Ecuación 35, cuando aumenta la complementariedad en al menos un factor de dos o cuando la inercia calculada es incorrecta (el problema parece no convexo). En estos casos, el algoritmo intenta dar un paso diferente o un paso de gradiente conjugado.
Paso de gradiente conjugado
El enfoque de gradiente conjugado empleado para resolver el problema aproximado Ecuación 34 es similar a otros cálculos de gradiente conjugado. En este caso, el algoritmo ajusta tanto x como s, manteniendo positivas las holguras de s. El enfoque consiste en minimizar una aproximación cuadrática al problema aproximado en una región de confianza, sujeta a restricciones linealizadas.
En concreto, R denota el radio de la región de confianza y las demás variables se definen como en Paso directo. El algoritmo obtiene los multiplicadores de Lagrange resolviendo aproximadamente las ecuaciones de KKT
en el sentido de mínimos cuadrados, siempre que λ sea positivo. A continuación, da un paso (Δx, Δs) para resolver aproximadamente
(40) |
(41) |
Modo de factibilidad
Cuando la opción EnableFeasibilityMode
es true
y las iteraciones no disminuyen la infactibilidad con suficiente rapidez, el algoritmo pasa al modo de factibilidad. Este cambio se produce después de que el algoritmo no consiga disminuir la infactibilidad en modo normal y vuelva a fallar después de cambiar al modo de gradiente conjugado. Por lo tanto, para obtener el mejor rendimiento cuando el solver no encuentra una solución factible sin el modo de factibilidad, establezca el SubproblemAlgorithm
en 'cg'
cuando utilice el modo de factibilidad. De esta forma se evitan búsquedas sin resultado en el modo normal.
El algoritmo del modo de factibilidad se basa en Nocedal, Öztoprak y Waltz [1]. El algoritmo ignora la función objetivo y en su lugar intenta minimizar la infactibilidad, definida como la suma de las partes positivas de las funciones de restricción de desigualdad y el valor absoluto de las funciones de restricción de igualdad. En términos de las variables de relajación , que corresponden a desigualdades, partes positivas de igualdades y partes negativas de igualdades, respectivamente, el problema es
sujeto a las restricciones
Para resolver el problema relajado, el programa utiliza una formulación interior-point con una función de barrera logarítmica y las holguras para minimizar
sujeto a las restricciones
El proceso de solución del problema relajado comienza con μ inicializado en el valor actual del parámetro de barrera. La variable de holgura sI se inicializa en el valor de holgura de desigualdad actual, heredado del modo principal. Las variables r se inicializan en
El resto de holguras se inicializan en
A partir de este punto inicial, el algoritmo del modo de factibilidad reutiliza el código del algoritmo interior-point normal. Este proceso requiere cálculos de pasos especiales porque las variables r son lineales y, por lo tanto, sus segundas derivadas asociadas son cero. En otras palabras, la función objetivo hessiana del problema de factibilidad es de rango deficiente. Por lo tanto, el algoritmo no puede dar un paso Newton. En su lugar, el algoritmo da un paso en la dirección de descenso más pronunciada. El algoritmo parte del gradiente del objetivo con respecto a las variables, proyecta el gradiente sobre el espacio nulo de la jacobiana de las restricciones y reescala el vector resultante para que tenga una longitud de paso adecuada. Este paso puede ser eficaz para reducir la infactibilidad.
El algoritmo del modo de factibilidad finaliza cuando reduce la infactibilidad por un factor de 10. Cuando finaliza el modo de factibilidad, el algoritmo pasa las variables x y sI al algoritmo principal y descarta las demás variables de holgura y las variables de relajación r.
Referencias
[1] Nocedal, Jorge, Figen Öztoprak, and Richard A. Waltz. An Interior Point Method for Nonlinear Programming with Infeasibility Detection Capabilities. Optimization Methods & Software 29(4), July 2014, pp. 837–854.
Opciones del algoritmo interior-point
Estos son los significados y efectos de varias opciones del algoritmo interior-point.
HonorBounds
: cuando se establece entrue
, cada iteración satisface las restricciones de límite que ha establecido. Cuando se establece enfalse
, el algoritmo puede vulnerar los límites durante las iteraciones intermedias.HessianApproximation
: cuando se establece en:'bfgs'
,fmincon
calcula la matriz hessiana mediante una aproximación cuasi-Newton densa.'lbfgs'
,fmincon
calcula la matriz hessiana mediante una aproximación cuasi-Newton a gran escala con memoria limitada.'fin-diff-grads'
,fmincon
calcula un producto del vector por la matriz hessiana mediante diferencias finitas del gradiente o los gradientes; las demás opciones deben configurarse adecuadamente.
HessianFcn
:fmincon
usa el identificador de función que especifique enHessianFcn
para calcular la hessiana. Consulte Incluir matrices hessianas.HessianMultiplyFcn
: da una función separada para evaluar el vector por la matriz hessiana. Para obtener más detalles, consulte Incluir matrices hessianas y Función de multiplicación de matriz hessiana.SubproblemAlgorithm
: determina si se intenta o no el paso directo de Newton. El ajuste por defecto'factorization'
permite intentar este tipo de paso. El ajuste'cg'
solo permite pasos de gradiente conjugado.
Para obtener una lista completa de opciones, consulte Algoritmo interior-point en options
de fmincon
.
Algoritmo fminbnd
fminbnd
es un solver disponible en cualquier instalación de MATLAB. Resuelve un mínimo local en una dimensión dentro de un intervalo acotado. No se basa en derivadas. En su lugar, utiliza la búsqueda por sección áurea y la interpolación parabólica.
Algoritmo y formulación de problemas fseminf
Formulación de problemas fseminf
fseminf
aborda problemas de optimización con tipos de restricciones adicionales en comparación con los que aborda fmincon
. La formulación de fmincon
es
de forma que c(x) ≤ 0, ceq(x) = 0, A·x ≤ b, Aeq·x = beq, and l ≤ x ≤ u.
fseminf
añade el siguiente conjunto de restricciones semiinfinitas a las ya dadas. En el caso de wj en un intervalo o rectángulo acotado unidimensional o bidimensional Ij y en el caso de un vector de funciones continuas K(x, w), las restricciones son
Kj(x, wj) ≤ 0 para todo wj∈Ij.
Por "dimensión" de un problema fseminf
se entiende la dimensión máxima del conjunto de restricciones I: 1 si todas las Ij son intervalos y 2 si al menos una Ij es un rectángulo. El tamaño del vector de K no entra en este concepto de dimensión.
La razón por la que se llama programación semiinfinita es que hay un número finito de variables (x y wj), pero un número infinito de restricciones. Esto se debe a que las restricciones sobre x son sobre un conjunto de intervalos continuos o rectángulos Ij, que contiene un número infinito de puntos, por lo que hay un número infinito de restricciones: Kj(x, wj) ≤ 0 para un número infinito de puntos wj.
Se podría pensar que un problema con un número infinito de restricciones es imposible de resolver. fseminf
lo aborda reformulando el problema a uno equivalente que tiene dos etapas: una maximización y una minimización. Las restricciones semiinfinitas se reformulan como
(42) |
donde |K| es el número de componentes del vector K; es decir, el número de funciones de restricción semiinfinitas. En el caso de x fijo, se trata de una maximización ordinaria sobre intervalos o rectángulos acotados.
fseminf
simplifica aún más el problema haciendo aproximaciones cuadráticas o cúbicas por tramos κj(x, wj) a las funciones Kj(x, wj), para cada x que visita el solver. fseminf
considera solo los máximos de la función de interpolación κj(x, wj), en lugar de Kj(x, wj), en Ecuación 42. De este modo, el problema original, minimizar una función con restricciones semiinfinitas, se reduce a un problema con un número finito de restricciones.
Puntos de muestreo. Su función de restricción semiinfinita debe proporcionar un conjunto de puntos de muestreo, que se utilizan para realizar las aproximaciones cuadráticas o cúbicas. Para ello, debe contener lo siguiente:
El espaciado inicial
s
entre puntos de muestreo wUna forma de generar el conjunto de puntos de muestreo w a partir de
s
El espaciado inicial s
es una matriz de |K| por 2. La j-ésima fila de s
representa el espaciado de los puntos de muestreo vecinos para la función de restricción Kj. Si Kj depende de una wj unidimensional, establezca s(j,2) = 0
. fseminf
actualiza la matriz s
en iteraciones posteriores.
fseminf
utiliza la matriz s
para generar los puntos de muestreo w, que luego emplea para crear la aproximación κj(x, wj). Su procedimiento para generar w a partir de s
debe mantener los mismos intervalos o rectángulos Ij durante la optimización.
Ejemplo de creación de puntos de muestreo. Consideremos un problema con dos restricciones semiinfinitas, K1 y K2. K1 tiene una dimensión w1, y K2 tiene dos dimensiones w2. El siguiente código genera un conjunto de muestreo de w1 = 2 a 12:
% Initial sampling interval if isnan(s(1,1)) s(1,1) = .2; s(1,2) = 0; end % Sampling set w1 = 2:s(1,1):12;
fseminf
especifica s
como NaN
cuando llama por primera vez a su función de restricción. Esta comprobación permite establecer el intervalo de muestreo inicial.
El siguiente código genera un conjunto de muestreo de w2 en un cuadrado, con cada componente que va de 1 a 100, que inicialmente se muestrea más a menudo en el primer componente que en el segundo:
% Initial sampling interval if isnan(s(1,1)) s(2,1) = 0.2; s(2,2) = 0.5; end % Sampling set w2x = 1:s(2,1):100; w2y = 1:s(2,2):100; [wx,wy] = meshgrid(w2x,w2y);
Los fragmentos de código anteriores pueden simplificarse de la siguiente forma:
% Initial sampling interval if isnan(s(1,1)) s = [0.2 0;0.2 0.5]; end % Sampling set w1 = 2:s(1,1):12; w2x = 1:s(2,1):100; w2y = 1:s(2,2):100; [wx,wy] = meshgrid(w2x,w2y);
Algoritmo fseminf
fseminf
reduce el problema de la programación semiinfinita a un problema abordado por fmincon
. fseminf
sigue estos pasos para resolver problemas de programación semiinfinita:
En el valor actual de x,
fseminf
identifica todas las wj,i de tal manera que la interpolación κj(x, wj,i) es un máximo local. (El máximo se refiere a variar w para x fijo).fseminf
da un paso de iteración en la solución del problema defmincon
:tal que c(x) ≤ 0, ceq(x) = 0, A·x ≤ b, Aeq·x = beq, and l ≤ x ≤ u, donde c(x) se aumenta con todos los máximos de κj(x, wj) tomados sobre todo wj∈Ij, que es igual a los máximos sobre j e i de κj(x, wj,i).
fseminf
comprueba si se cumple algún criterio de detención en el nuevo punto x (para detener las iteraciones); en caso contrario, continúa con el paso 4.fseminf
comprueba si es necesario actualizar la discretización de las restricciones semiinfinitas y actualiza los puntos de muestreo según corresponda. Esto proporciona una aproximación actualizada κj(x, wj). Luego continúa en el paso 1.