La traducción de esta página aún no se ha actualizado a la versión más reciente. Haga clic aquí para ver la última versión en inglés.

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:

$$\underset{x}{\min }f(x)$$

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,

$$\underset{s}{\min }\left\{q(s),s\in N\right\}.$$

(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

$$\min \left\{\frac{1}{2}{s}^{T}Hs+s^{T}g\text{ such that }\Vert Ds\Vert \le \Delta \right\},$$

(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

$$\frac{1}{\Delta }-\frac{1}{\Vert s\Vert }=0.$$

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 s₁ y s₂, donde s₁ está en la dirección del gradiente g y s₂ es una dirección de Newton aproximada, es decir, una solución para

$$H\cdot s_2=-g,$$

(3)

o una dirección de curvatura negativa,

$$s_2^T\cdot H\cdot s_2<0.$$

(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:

1. Formule el subproblema trust-region bidimensional.

2. Resuelva Ecuación 2 para determinar el paso de prueba s.

3. Si f(x + s) < f(x), entonces x = x + s.

4. 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 = C², 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, d^THd ≤ 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

$$\min \left\{f(x)\text{ such that }Ax=b\right\},$$

(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 A^T. 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 x₀, utilizando un paso de mínimos cuadrados dispersos, de modo que Ax₀ = 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

$$\left[\begin{array}{cc}C& {\tilde{A}}^T\\ \tilde{A}& 0\end{array}\right]\left[\begin{array}{c}s\\ t\end{array}\right]=\left[\begin{array}{c}r\\ 0\end{array}\right],$$

(6)

donde $$\tilde{A}$$ 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

$$\min \left\{f(x)\text{ such that }l\le x\le u\right\},$$

(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,

$${\left(D(x)\right)}^{-2}g=0,$$

(8)

donde

$$D(x)=\text{diag}\left({\left|v_k\right|}^{-1/2}\right),$$

y el vector v(x) se define a continuación, para cada 1 ≤ i ≤ n:

• Si g_i < 0 y u_i < ∞, entonces v_i = x_i – u_i

• Si g_i ≥ 0 y l_i > –∞, entonces v_i = x_i – l_i

• Si g_i < 0 y u_i = ∞, entonces v_i = –1

• Si g_i ≥ 0 y l_i = –∞, entonces v_i = 1

El sistema no lineal Ecuación 8 no puede diferenciarse en todas partes. La no diferenciabilidad se produce cuando v_i = 0. Se pueden evitar estos puntos manteniendo una factibilidad estricta, es decir, restringiendo l < x < u.

El paso de Newton modificado a escala s_k para el sistema de ecuaciones no lineal dado por Ecuación 8 se define como la solución del sistema lineal

$$\widehat{M}D{s}^N=-\widehat{g}$$

(9)

en la k-ésima iteración, donde

$$\widehat{g}=D^{-1}g=\text{diag}\left({\left|v\right|}^{1/2}\right)g,$$

(10)

y

$$\widehat{M}=D^{-1}H{D}^{-1}+\text{diag}(g)J^v.$$

(11)

Aquí J^v desempeña el papel de la matriz jacobiana de |v|. Cada componente diagonal de la matriz diagonal J^v es igual a 0, -1 o 1. Si todos los componentes de l y u son finitos, J^v = diag(sign(g)). En un punto en el que g_i = 0, v_i podría no ser diferenciable. $$J_{ii}^v=0$$ 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 v_i. Además, |v_i| seguirá siendo discontinua en este punto, pero la función |v_i|·g_i 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 p^R = p, excepto en el i-ésimo componente, donde p^R_i = –p_i.

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 G_i(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

$$\begin{array}{c}\nabla f\left(x^{*}\right)+{\displaystyle \sum _{i=1}^m{\lambda }_i\cdot \nabla G_i\left(x^{*}\right)}=0\\ {\lambda }_i\cdot G_i\left(x^{*}\right)=0,i=1,\mathrm{...},m_e\\ {\lambda }_i\ge 0,i=m_e+1,\mathrm{...},m,\end{array}$$

(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.

$$L(x,\lambda )=f(x)+{\displaystyle \sum _{i=1}^m{\lambda }_i\cdot g_i(x)}.$$

(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)

$$\begin{array}{l}\underset{d\in {\Re }^n}{\min }\frac{1}{2}{d}^T{H}_{k}d+\nabla f{\left(x_k\right)}^{T}d\\ \nabla g_i{\left(x_k\right)}^{T}d+g_i\left(x_k\right)=0,i=1,\mathrm{...},m_e\\ \nabla g_i{\left(x_k\right)}^{T}d+g_i\left(x_k\right)\le 0,i=m_e+1,\mathrm{...},m.\end{array}$$

(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

x_{k + 1} = x_k + α_kd_k.

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 H_k es una aproximación definida positiva de la matriz hessiana de la función lagrangiana (Ecuación 13). H_k 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),

$$x_1^2+x_2^2-1.5\le 0.$$

(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

• Solución de programación cuadrática

• Inicialización

• Búsqueda de recta y función de mérito

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.

$$H_{k+1}=H_k+\frac{q_k{q}_k^T}{q_k^T{s}_k}-\frac{H_k{s}_k{s}_k^T{H}_k^T}{s_k^T{H}_k{s}_k},$$

(16)

donde

$$\begin{array}{c}{s}_k=x_{k+1}-x_k\\ q_k=\left(\nabla f\left(x_{k+1}\right)+{\displaystyle \sum _{i=1}^m{\lambda }_i\cdot \nabla g_i\left(x_{k+1}\right)}\right)-\left(\nabla f\left(x_k\right)+{\displaystyle \sum _{i=1}^m{\lambda }_i\cdot \nabla g_i\left(x_k\right)}\right).\end{array}$$

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 $$q_k^T{s}_k$$ sea positivo en cada actualización y que H se inicialice con una matriz definida positiva. Cuando $$q_k^T{s}_k$$ no es positivo, q_k se modifica elemento a elemento para que $$q_k^T{s}_k>0$$. El objetivo general de esta modificación es distorsionar lo menos posible los elementos de q_k, 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 q_k*s_k se reduce repetidamente a la mitad. Este procedimiento se continúa hasta que $$q_k^T{s}_k$$ sea mayor o igual que una pequeña tolerancia negativa. Si, después de este procedimiento, $$q_k^T{s}_k$$ sigue sin ser positivo, modifique q_k añadiendo un vector v multiplicado por un escalar constante w, es decir,

$$q_k=q_k+wv,$$

(17)

donde

$$\begin{array}{l}{v}_i=\nabla g_i\left(x_{k+1}\right)\cdot g_i\left(x_{k+1}\right)-\nabla g_i\left(x_k\right)\cdot g_i\left(x_k\right)\\ \text{ if }{\left(q_k\right)}_i\cdot w<0\text{ and }{\left(q_k\right)}_i\cdot {\left(s_k\right)}_i<0,i=1,\mathrm{...},m\\ v_i=0\text{ otherwise,}\end{array}$$

y aumente w sistemáticamente hasta que $$q_k^T{s}_k$$ 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 $$q_k^T{s}_k$$ 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 A_i se refiere a la i-ésima fila de la matriz A de m por n.

$$\begin{array}{c}\underset{d\in {\Re }^n}{\min }q(d)=\frac{1}{2}{d}^{T}Hd+c^{T}d,\\ A_{i}d=b_i,i=1,\mathrm{...},m_e\\ A_{i}d\le b_i,i=m_e+1,\mathrm{...},m.\end{array}$$

(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, $${\overline{A}}_k$$, 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.

$${\overline{A}}_k$$ se actualiza en cada iteración k y esto se utiliza como base para una dirección de búsqueda $${\widehat{d}}_k$$. Las restricciones de igualdad siempre permanecen en el conjunto activo $${\overline{A}}_k$$. La notación de la variable $${\widehat{d}}_k$$ se utiliza aquí para distinguirla de d_k en las iteraciones principales del método SQP. Se calcula la dirección de búsqueda $${\widehat{d}}_k$$ y minimiza la función objetivo permaneciendo en cualquier límite de restricción activa. El subespacio factible para $${\widehat{d}}_k$$ se forma a partir de la base Z_k cuyas columnas son ortogonales a la estimación del conjunto activo $${\overline{A}}_k$$ (es decir, $${\overline{A}}_k{Z}_k=0$$). 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 Z_k permanezca en los límites de las restricciones activas.

La matriz Z_k se forma a partir de las últimas columnas m – l de la descomposición de QR de la matriz $${\overline{A}}_k^T$$, donde l es el número de restricciones activas y l < m. Es decir, Z_k viene dado por

$$Z_k=Q\left[:,l+1:n\right],$$

(19)

donde

$$Q^T{\overline{A}}_k^T=\left[\begin{array}{c}R\\ 0\end{array}\right].$$

Una vez encontrado Z_k, se busca una nueva dirección de búsqueda $${\widehat{d}}_k$$ que minimice q(d) donde $${\widehat{d}}_k$$ está en el espacio nulo de las restricciones activas. Es decir, $${\widehat{d}}_k$$ es una combinación lineal de las columnas de Z_k: $${\widehat{d}}_k=Z_{k}p$$ para algún vector p.

Entonces, si vemos la cuadrática como función de p, sustituyendo por $${\widehat{d}}_k$$, se obtiene

$$q(p)=\frac{1}{2}{p}^T{Z}_k^{T}H{Z}_{k}p+c^T{Z}_{k}p.$$

(20)

Si se diferencia con respecto a p se obtiene

$$\nabla q(p)=Z_k^{T}H{Z}_{k}p+Z_k^{T}c.$$

(21)

∇q(p) se denomina gradiente proyectado de la función cuadrática porque es el gradiente proyectado en el subespacio definido por Z_k. El término $$Z_k^{T}H{Z}_k$$ 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 Z_k ocurre cuando ∇q(p) = 0, que es la solución del sistema de ecuaciones lineales

$$Z_k^{T}H{Z}_{k}p=-Z_k^{T}c.$$

(22)

A continuación, se da un paso con el formato

$$x_{k+1}=x_k+\alpha {\widehat{d}}_k,\text{ where }{\widehat{d}}_k=Z_{k}p.$$

(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 $${\widehat{d}}_k$$ es el paso exacto para llegar al mínimo de la función restringida al espacio nulo de $${\overline{A}}_k$$. 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 $${\widehat{d}}_k$$ 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 $${\widehat{d}}_k$$ viene dada por

$$\alpha =\underset{i\in \left\{1,\mathrm{...},m\right\}}{\min }\left\{\frac{-\left(A_i{x}_k-b_i\right)}{A_i{\widehat{d}}_k}\right\},$$

(24)

que se define en el caso de las restricciones que no están en el conjunto activo y donde la dirección $${\widehat{d}}_k$$ es hacia el límite de la restricción, es decir, $$A_i{\widehat{d}}_k>0,i=1,\mathrm{...},m$$.

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

$${\overline{A}}_k^T{\lambda }_k=c+H{x}_k.$$

(25)

Si todos los elementos de λ_k son positivos, x_k 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

$$\begin{array}{l}\underset{\gamma \in \Re ,x\in {\Re }^n}{\min }\gamma \text{ such that}\\ A_{i}x=b_i,i=1,\mathrm{...},m_e\\ A_{i}x-\gamma \le b_i,i=m_e+1,\mathrm{...},m.\end{array}$$

(26)

La notación A_i 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 g_k es el gradiente de la función objetivo (igual a los coeficientes de la función objetivo lineal).

$${\widehat{d}}_k=-Z_k{Z}_k^T{g}_k.$$

(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 $${\widehat{d}}_k$$ se inicializa con una dirección de búsqueda $${\widehat{d}}_1$$ hallada a partir de la resolución del conjunto de ecuaciones lineales

$$H{\widehat{d}}_1=-g_k,$$

(28)

donde g_k es el gradiente de la función objetivo en el iterado actual x_k (es decir, Hx_k + 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 $${\widehat{d}}_k$$ 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 d_k, que se utiliza para formar un nuevo iterado

$$x_{k+1}=x_k+\alpha d_k.$$

(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.

$$\Psi (x)=f(x)+{\displaystyle \sum _{i=1}^{m_e}{r}_i\cdot g_i(x)}+{\displaystyle \sum _{i=m_e+1}^m{r}_i\cdot \max [0,g_i(x)]}.$$

(30)

Powell recomienda fijar el parámetro de penalización

$$r_i={\left(r_{k+1}\right)}_i=\underset{i}{\max }\left\{{\lambda }_i,\frac{{\left(r_k\right)}_i+{\lambda }_i}{2}\right\},i=1,\mathrm{...},m.$$

(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 r_i se fija inicialmente en

$$r_i=\frac{\Vert \nabla f(x)\Vert }{\Vert \nabla g_i(x)\Vert },$$

(32)

donde $$\Vert \Vert$$ 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 algoritmo sqp 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

$$\underset{x}{\min }f(x),\text{ subject to }h(x)=0\text{ and }g(x)\le 0.$$

(33)

Para cada μ > 0, el problema aproximado es

$$\begin{array}{l}\underset{x,s}{\min }{f}_{\mu }(x,s)=\underset{x,s}{\min }f(x)-\mu {\displaystyle \sum _i\ln \left(s_i\right)},\text{ subject to }s\ge 0,h(x)=0,\text{ and }g(x)+s=0.\\ \end{array}$$

(34)

Hay tantas variables de holgura s_i como restricciones de desigualdad g. s_i se restringen a positivo para mantener los iterados en el interior de la región factible. A medida que μ disminuye hasta cero, el mínimo de f_μ debe acercarse al mínimo de f. El término logarítmico añadido se denomina función de barrera. Este método se describe en [40], [41] y [51].

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

$$f_{\mu }(x,s)+\nu \Vert \left(h(x),g(x)+s\right)\Vert .$$

(35)

El parámetro $$\nu$$ puede aumentar con el número de iteraciones para forzar a la solución hacia la factibilidad. Si un paso que se ha intentado no disminuye la función de mérito, el algoritmo lo rechaza y prueba un nuevo paso.

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 x_j, el algoritmo rechaza x_j. 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_μ:

  $$H={\nabla }^{2}f(x)+{\displaystyle \sum _i{\lambda }_i{\nabla }^2{g}_i(x)}+{\displaystyle \sum _j{y}_j{\nabla }^2{h}_j(x)}.$$

  (36)

• J_g denota la jacobiana de la función de restricción g.

• J_h 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):

$$\left[\begin{array}{cccc}H& 0& J_h^T& J_g^T\\ 0& \Lambda & 0& S\\ J_h& 0& 0& 0\\ J_g& I& 0& 0\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta s\\ \Delta y\\ \Delta \lambda \end{array}\right]=-\left[\begin{array}{c}\nabla f+J_h^{T}y+J_g^T\lambda \\ S\lambda -\mu e\\ h\\ g+s\end{array}\right].$$

(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:

$$\left[\begin{array}{cccc}H& 0& J_h^T& J_g^T\\ 0& S\Lambda & 0& S\\ J_h& 0& 0& 0\\ J_g& S& 0& 0\end{array}\right]\left[\begin{array}{c}\Delta x\\ S^{-1}\Delta s\\ \Delta y\\ \Delta \lambda \end{array}\right]=-\left[\begin{array}{c}\nabla f(x)+J_h^{T}y+J_g^T\lambda \\ S\lambda -\mu e\\ h\\ g(x)+s\end{array}\right].$$

(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 μ:

$$\max \left(\Vert \nabla f(x)+J_h^{T}y+J_g^T\lambda \Vert ,\Vert S\lambda -\mu e\Vert ,\Vert h\Vert ,\Vert g(x)+s\Vert \right)<\mu .$$

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:

$$\left[\begin{array}{cccc}H& 0& J_h^T& J_g^T\\ 0& \Lambda & 0& S\\ J_h& 0& 0& 0\\ J_g& I& 0& 0\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta s\\ \Delta y\\ \Delta \lambda \end{array}\right]=-\left[\begin{array}{c}\nabla f+J_h^{T}y+J_g^T\lambda \\ S\lambda \\ h\\ g+s\end{array}\right].$$

Defina ɑ_s y ɑ_λ como los tamaños de paso más grandes que no vulneren las restricciones de no negatividad.

$$\begin{array}{l}{\alpha }_s=\max \left(\alpha \in (0,1]:s+\alpha \Delta s\ge 0\right)\\ {\alpha }_{\lambda }=\max \left(\alpha \in (0,1]:\lambda +\alpha \Delta \lambda \ge 0\right).\end{array}$$

Ahora calcule la complementariedad a partir del paso predictor.

$${\mu }_P=\frac{\left(s+{\alpha }_s\Delta s\right)\left(\lambda +{\alpha }_{\lambda }\Delta \lambda \right)}{m},$$

(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

$$\left(s+\Delta s\right)\left(\lambda +\text{}\Delta \lambda \right)=\underset{\text{Linear term set to 0}}{\underbrace{s\lambda +\text{}s\Delta \lambda +\lambda \Delta s}}+\underset{\text{Quadratic}}{\underbrace{\Delta s\Delta \lambda }}.$$

Para corregir el error cuadrático, resuelva el sistema lineal para la dirección del paso corrector.

$$\left[\begin{array}{cccc}H& 0& J_h^T& J_g^T\\ 0& \Lambda & 0& S\\ J_h& 0& 0& 0\\ J_g& I& 0& 0\end{array}\right]\left[\begin{array}{c}\Delta x_{\text{cor}}\\ \Delta s_{\text{cor}}\\ \Delta y_{\text{cor}}\\ \Delta {\lambda }_{\text{cor}}\end{array}\right]=-\left[\begin{array}{c}0\\ \Delta s\Delta \lambda \\ 0\\ 0\end{array}\right].$$

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

$$\left[\begin{array}{cccc}H& 0& J_h^T& J_g^T\\ 0& \Lambda & 0& S\\ J_h& 0& 0& 0\\ J_g& I& 0& 0\end{array}\right]\left[\begin{array}{c}\Delta x_{\text{cen}}\\ \Delta s_{\text{cen}}\\ \Delta y_{\text{cen}}\\ \Delta {\lambda }_{\text{cen}}\end{array}\right]=-\left[\begin{array}{c}\nabla f+J_h^{T}y+J_g^T\lambda \\ S\lambda -\mu e\sigma \\ h\\ g+s\end{array}\right].$$

Aquí, σ se define como

$$\sigma ={\left(\frac{{\mu }_P}{\mu }\right)}^3,$$

donde μ_P se define en la ecuación Ecuación 39 y $$\mu =s^T\lambda /m$$.

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

$${\nabla }_{x}L={\nabla }_{x}f(x)+{\displaystyle \sum _i{\lambda }_i\nabla g_i(x)}+{\displaystyle \sum _j{y}_j\nabla h_j(x)}=0,$$

en el sentido de mínimos cuadrados, siempre que λ sea positivo. A continuación, da un paso (Δx, Δs) para resolver aproximadamente

$$\underset{\Delta x,\Delta s}{\min }\nabla f^T\Delta x+\frac{1}{2}\Delta x^T{\nabla }_{xx}^{2}L\Delta x+\mu e^T{S}^{-1}\Delta s+\frac{1}{2}\Delta s^T{S}^{-1}\Lambda \Delta s,$$

(40)

sujeto a las restricciones linealizadas

$$g(x)+J_g\Delta x+\Delta s=0,h(x)+J_h\Delta x=0.$$

(41)

Para resolver Ecuación 41, el algoritmo intenta minimizar una norma de las restricciones linealizadas dentro de una región con radio escalado por R. Después se resuelve Ecuación 40 con las restricciones de igualar el residual de resolver Ecuación 41, de permanecer dentro de la región de confianza de radio R y de mantener s estrictamente positivo. Para obtener más información sobre el algoritmo y la derivación, consulte [40], [41] y [51]. Para consultar otra descripción de los gradientes conjugados, consulte Método del gradiente conjugado precondicionado.

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 $$r=\left[r_I,r_e^{+},r_e^{-}\right]$$, que corresponden a desigualdades, partes positivas de igualdades y partes negativas de igualdades, respectivamente, el problema es

$$\underset{x,r}{\min }{1}^{T}r=\underset{x,r}{\min }{\displaystyle \sum r}$$

sujeto a las restricciones

$$\begin{array}{c}{r}_I\ge c(x)\\ r_e^{+}-r_e^{-}=ceq(x)\\ r\ge 0.\end{array}$$

Para resolver el problema relajado, el programa utiliza una formulación interior-point con una función de barrera logarítmica y las holguras $$s=\left[s_I,s_R,s_e^{+},s_e^{-}\right]$$ para minimizar

$$\underset{x,r,s}{\min }{1}^{T}r-\mu {\displaystyle \sum \left(\log \left(s_{I,i}\right)+\log \left(s_{R,i}\right)\right)}-{\displaystyle \sum \left(\log \left(s_{e,i}^{+}\right)+\log \left(s_{e,i}^{-}\right)\right)}$$

sujeto a las restricciones

$$\begin{array}{c}ceq(x)=r_e^{+}-r_e^{-}\\ c(x)=r_I-s_I\\ r_e^{+}=s_e^{+}\\ r_e^{-}=s_e^{-}\\ r_I=s_R\\ s\ge 0.\end{array}$$

El proceso de solución del problema relajado comienza con μ inicializado en el valor actual del parámetro de barrera. La variable de holgura s_I se inicializa en el valor de holgura de desigualdad actual, heredado del modo principal. Las variables r se inicializan en

$$\begin{array}{c}{r}_I=\max \left(s_I+c(x),\mu \right)\\ r_e^{+}=\max \left(ceq(x),\mu \right)\\ r_e^{-}=-\min \left(ceq(x),-\mu \right).\end{array}$$

El resto de holguras se inicializan en

$$\begin{array}{l}{s}_R=r_I\\ s_e^{+}=r_e^{+}\\ s_e^{-}=r_e^{-}.\end{array}$$

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 s_I 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 en true, cada iteración satisface las restricciones de límite que ha establecido. Cuando se establece en false, 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 en HessianFcn 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

$$\underset{x}{\min }f(x)$$

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 w_j en un intervalo o rectángulo acotado unidimensional o bidimensional I_j y en el caso de un vector de funciones continuas K(x, w), las restricciones son

K_j(x, w_j) ≤ 0 para todo w_j∈I_j.

Por "dimensión" de un problema fseminf se entiende la dimensión máxima del conjunto de restricciones I: 1 si todas las I_j son intervalos y 2 si al menos una I_j 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 w_j), 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 I_j, que contiene un número infinito de puntos, por lo que hay un número infinito de restricciones: K_j(x, w_j) ≤ 0 para un número infinito de puntos w_j.

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

$$\underset{w_j\in I_j}{\max }{K}_j(x,w_j)\le 0\text{ for all }j=1,\mathrm{...},\left|K\right|,$$

(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, w_j) a las funciones K_j(x, w_j), para cada x que visita el solver. fseminf considera solo los máximos de la función de interpolación κ_j(x, w_j), en lugar de K_j(x, w_j), 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 w

• Una 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 K_j. Si K_j depende de una w_j 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, w_j). Su procedimiento para generar w a partir de s debe mantener los mismos intervalos o rectángulos I_j durante la optimización.

Ejemplo de creación de puntos de muestreo.  Consideremos un problema con dos restricciones semiinfinitas, K₁ y K₂. K₁ tiene una dimensión w₁, y K₂ tiene dos dimensiones w₂. El siguiente código genera un conjunto de muestreo de w₁ = 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 w₂ 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:

1. En el valor actual de x, fseminf identifica todas las w_j,i de tal manera que la interpolación κ_j(x, w_j,i) es un máximo local. (El máximo se refiere a variar w para x fijo).

2. fseminf da un paso de iteración en la solución del problema de fmincon:

   $$\underset{x}{\min }f(x)$$

   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, w_j) tomados sobre todo w_j∈I_j, que es igual a los máximos sobre j e i de κ_j(x, w_j,i).

3. 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.

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, w_j). Luego continúa en el paso 1.