Algoritmos de optimización no lineales restringidos

Definición de optimización restringida

La minimización restringida es el problema de encontrar un vector que sea un mínimo local a una función escalar () sujeto a restricciones en el permitido:xfxx

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

tal que una o más de las siguientes retenciones: c(x) ≤ 0, ceq(x) = 0, A·x ≤ b, Aeq·x = beq, l ≤ x ≤ u. Hay aún más restricciones utilizadas en la programación semi-infinita; Ver.fseminf formulación y algoritmo del problema

Algoritmo reflexivo de la región Trust de fmincon

Métodos de la región de confianza para la minimización no lineal

Muchos de los métodos utilizados en los solucionadores se basan en un concepto simple pero potente en la optimización.Optimization Toolbox™trust regions,

Para comprender el enfoque de optimización de la región de confianza, considere el problema de minimización sin restricciones, minimice (), donde la función toma argumentos vectoriales y devuelve escalares.fx Supongamos que está en un punto en el espacio y desea mejorar, es decir, pasar a un punto con un valor de función inferior.xn La idea básica es aproximarse con una función más simple, que refleja razonablemente el comportamiento de la función en un vecindario alrededor del punto.fqfNx Este vecindario es la región de confianza. Un paso de prueba se calcula minimizando (o aproximadamente minimizando) más.sN Este es el subproblema de la región de confianza,

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

(1)

El punto actual se actualiza para ser + sixs f(x + s) < f(x); de lo contrario, el punto actual permanece sin cambios y, la región de confianza, se encogió y se repite el cálculo del paso de prueba.N

Las preguntas clave en la definición de un enfoque específico de la región de confianza para minimizar () son cómo elegir y calcular la aproximación (definida en el punto actual), cómo elegir y modificar la región de confianza y cómo resolver con precisión el subproblema de la región de confianza.fxqxN Esta sección se centra en el problema sin restricciones. Las secciones posteriores discuten complicaciones adicionales debido a la presencia de restricciones en las variables.

En el método de la región de confianza estándar (), la aproximación cuadrática se define por los dos primeros términos de la aproximación de Taylor a en; la vecindad es generalmente esférica o elipsoidal en forma.[48]qFxN Matemáticamente el subproblema de la región de confianza normalmente se indica

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

(2)

donde está el gradiente de en el punto actual, es la matriz de hessian (la matriz simétrica de segundo derivados), es una matriz de escalado diagonal, Δ es un escalar positivo, y ∥. ∥ es la norma 2.gfxHD Existen buenos algoritmos para resolver (ver); Estos algoritmos suelen implicar el cálculo de todos los valores propios y un proceso de Newton aplicado a laEcuación 2[48]H ecuación secular

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

Estos algoritmos proporcionan una solución precisa.Ecuación 2 Sin embargo, requieren tiempo proporcional a varias factorizaciones de.H Por lo tanto, para problemas a gran escala es necesario un enfoque diferente. Varias estrategias de aproximación y heurística, basadas en, se han propuesto en la literatura (y).Ecuación 2[42][50] El enfoque de aproximación seguido en los solucionadores es restringir el subproblema de la región de confianza a un subespacio bidimensional (y).Optimization ToolboxS[39][42] Una vez que se ha calculado el subespacio, el trabajo a resolver es trivial incluso si se necesita información completa de eigenvalue/Eigenvector (ya que en el subespacio, el problema es sólo bidimensional).SEcuación 2 La obra dominante se ha desplazado ahora a la determinación del subespacio.

El subespacio bidimensional se determina con la ayuda de unS proceso de gradiente conjugada precondicionado descrito a continuación. El solucionador define como el espacio lineal extendido porSs₁ Ys₂Dóndes₁ está en la dirección del degradado ygs₂ es un aproximado Newton Dirección, 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 detrás de esta elección es forzar la convergencia global (a través de la dirección de descenso más pronunciada o dirección de curvatura negativa) y lograr una rápida convergencia local (a través del paso de Newton, cuando existe).S

Un boceto de minimización sin restricciones utilizando ideas de la región de confianza es ahora fácil de dar:

1. Formule el subproblema de la región de confianza bidimensional.

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

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 según las reglas estándar. En particular, se disminuye si no se acepta el paso de prueba, es decir, f(x + s) ≥ f(x). Ver y para una discusión de este aspecto.[46][49]

solucionadores tratan algunos casos especiales importantes de funciones especializadas: mínimos cuadrados no lineales, funciones cuadráticas y mínimos cuadrados lineales.Optimization Toolboxf Sin embargo, las ideas algorítmicas subyacentes son las mismas que para el caso general. Estos casos especiales se discuten en secciones posteriores.

Método de gradiente conjugada precondicionado

Una forma popular de resolver los grandes sistemas de ecuaciones lineales definidas de manera simétrica y positiva Hp = –g es el método de los degradados de conjugación preacondicionados (PCG). Este enfoque iterativo requiere la capacidad de calcular los productos de vector de matriz de la forma, donde es un vector arbitrario.H·vv La matriz definida simétrica positiva es unaM preacondicionador Para.H Es decir M = C²Dónde C^–1HC^–1 es una matriz bien acondicionadas o una matriz con valores eigenados agrupados.

En un contexto de minimización, puede suponer que la matriz de Hessian es simétrica.H Sin embargo, está garantizado para ser positivo definitivo sólo en el vecindario de un fuerte minimizador.H El algoritmo PCG sale cuando se encuentra una dirección de curvatura negativa (o cero), es decir, d^THd ≤ 0. La dirección de salida PCG,, o bien una dirección de curvatura negativa o una solución aproximada (controla la aproximación) al sistema Newtonptol Hp = –g. En ambos casos se utiliza para ayudar a definir el subespacio bidimensional utilizado en el enfoque de la región de confianza discutido en.pMétodos de la región de confianza 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 una manera limpia y eficiente. Los métodos de la región de confianza en los solucionadores generan iterados estrictamente factibles.Optimization Toolbox

El problema general de minimización limitada de la igualdad lineal se puede escribir

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

(5)

donde es una-por-matriz (Amnm ≤ n). Algunos solucionadores preprocesan para eliminar las dependencias lineales estrictas mediante una técnica basada en la factorización de LU deOptimization ToolboxA Un^T .[46] Aquí se supone que es de rango.Am

El método utilizado para resolver difiere del enfoque sin restricciones de dos maneras significativas.Ecuación 5 En primer lugar, un punto factible inicialx₀ se calcula, utilizando un paso de mínimos cuadrados dispersos, para que Ax₀ = b. En segundo lugar, el algoritmo PCG se sustituye por degradados de conjugados preacondicionados reducidos (RPCG), ver, para calcular un paso de Newton reducido aproximado (o una dirección de curvatura negativa en el espacio nulo de).[46]A El paso de álgebra lineal clave implica la resolución de 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)

Dónde $$\tilde{A}$$ aproximados (los pequeños ceros de se fijan a cero rango proporcionado no se pierde) y es una aproximación simétrica positiva-definida escasa a, es decir,AACH C = H. Vea para más detalles.[46]

Restricciones de cuadro

El problema de la caja restringida es de la forma

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

(7)

donde es un vector de límites inferiores, y es un vector de los límites superiores.lu Algunos (o todos) de los componentes de puede ser igual a – ∞ y algunos (o todos) de los componentes de puede ser igual a ∞.lu El método genera una secuencia de puntos estrictamente factibles. Se utilizan dos técnicas para mantener la viabilidad mientras se logra un comportamiento de convergencia robusto. En primer lugar, un paso de Newton modificado a escala reemplaza el paso sin restricciones de Newton (para definir el subespacio bidimensional).S Segundo reflexiones se utilizan para aumentar el tamaño del paso.

El paso de Newton modificado a escala surge de examinar las condiciones necesarias de Kuhn-Tucker para,Ecuación 7

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

(8)

Dónde

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

y el vector () se define a continuación, para cadavx 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 no es diferenciable en todas partes.Ecuación 8 La no diferenciabilidad se produce cuando v_i = 0. Puede evitar estos puntos manteniendo una viabilidad estricta, es decir, restringiendo l < x < u.

El paso modificado de Newton escalado s_k para el sistema no lineal de ecuaciones dado por se define como la solución al sistema linealEcuación 8

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

(9)

en la iteración TH, dondek

$$\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 del jacobiano de | |.v Cada componente diagonal de la matriz diagonal J^v es igual a 0, – 1 o 1. Si todos los componentes de y son finitos,lu J^v = diag(sign(g)). En un punto donde g_i = 0, V_i podría no ser diferenciable. $$J_{ii}^v=0$$ se define en ese punto. La no diferenciabilidad de este tipo no es motivo de preocupación porque, para tal componente, no es significativo qué valor V_i Toma. Además, |V_i| seguirá siendo discontinua en este punto, pero la función |v_i|·g_i es continua.

Segundo reflexiones se utilizan para aumentar el tamaño del paso. Un paso de reflexión (único) se define de la siguiente manera. Dado un paso que intersecta una restricción enlazada, considere la primera restricción enlazada cruzada por; suponer que es la restricción de límite de TH (ya sea el límite superior o inferior de la TH).ppiii A continuación, el paso de reflexión p^R = p excepto en el componente TH, dondei p^R_i = –p_i.

fmincon Active Set Algorithm

Introducción

En la optimización restringida, el objetivo general es transformar el problema en un subproblema más fácil que, a continuación, se puede resolver y utilizar como base de un proceso iterativo. Una característica de una clase grande de métodos tempranos es la traducción del problema restringido a un problema básico sin restricciones mediante el uso de una función de penalización para las restricciones que están cerca o más allá del límite de restricción. De este modo, el problema restringido se resuelve utilizando una secuencia de optimizaciones no restringidas parametrizadas, que en el límite (de la secuencia) convergen en el problema restringido. Estos métodos se consideran ahora relativamente ineficaces y se han sustituido por métodos que se han centrado en la solución del Ecuaciones de Karush-Kuhn-Tucker (KKT). Las ecuaciones KKT son condiciones necesarias para la optimalidad de un problema de optimización restringido. Si el problema es un llamado problema de programación convexa, es decir, () yfx G_i(x), i = 1,...,m, son funciones convexas, entonces las ecuaciones KKT son necesarias y suficientes para un punto de solución global.

Refiriéndose a GP (), las ecuaciones Kuhn-Tucker pueden ser indicadas comoEcuación 1

$$\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 degradados entre la función objetiva y las restricciones activas en el punto de solución. Para que se cancelen los degradados, los multiplicadores de Lagrange (Λ_i, = 1,...,) son necesarios para equilibrar las desviaciones en magnitud de la función objetiva y los degradados de restricción.im Dado que solo se incluyen las restricciones activas en esta operación de cancelación, las restricciones que no están activas no deben incluirse en esta operación y, por lo tanto, se les da a los multiplicadores de Lagrange igual a 0. Esto se indica implícitamente en las dos últimas ecuaciones de Kuhn-Tucker.

La solución de las ecuaciones KKT constituye la base de muchos algoritmos de programación no lineales. Estos algoritmos intentan computar los multiplicadores de Lagrange directamente. Los métodos limitados de quasi-Newton garantizan la convergencia superlineal mediante la acumulación de información de segundo orden con respecto a las ecuaciones KKT utilizando un procedimiento de actualización de quasi-Newton. Estos métodos se conocen comúnmente como métodos de programación cuadrática secuencial (SQP), ya que un subproblema de QP se resuelve en cada iteración principal (también conocida como programación cuadrática iterativa, programación cuadrática recursiva y métrica variable restringida métodos).

El algoritmo no es un algoritmo a gran escala; Ver.'active-set'Algoritmos a gran escala frente a mediano escala

Programación cuadrática secuencial (SQP)

Los métodos SQP representan el estado de la técnica en los métodos de programación no lineales. Schittkowski, por ejemplo, ha implementado y probado una versión que supera a todos los demás métodos probados en términos de eficiencia, precisión y porcentaje de soluciones exitosas, en un gran número de problemas de prueba.[36]

Basado en el trabajo de Biggs, han, y Powell (y), el método le permite imitar estrechamente el método de Newton para la optimización restringida al igual que se hace para la optimización sin restricciones.[1][22][32][33] En cada iteración principal, se hace una aproximación de la función de hessian de la Lagrangian utilizando un método de actualización quasi-Newton. A continuación, se utiliza para generar un subproblema de QP cuya solución se utiliza para formar una dirección de búsqueda para un procedimiento de búsqueda de línea. Una visión general de SQP se encuentra en Fletcher, Gill et al.  , Powell, y Schittkowski.[13][19][35][23] El método general, sin embargo, se indica aquí.

Dada la descripción del problema en GP () la idea principal es la formulación de un subproblema de QP basado en una aproximación cuadrática de la función Lagrangia.Ecuación 1

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

(13)

Aquí se simplifica suponiendo que las restricciones enlazadas se han expresado como restricciones de desigualdad.Ecuación 1 Para obtener el subproblema de QP, linealizar las restricciones no lineales.

Programación cuadrática (QP) subproblema

$$\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 se puede resolver utilizando cualquier algoritmo QP (véase, 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 de búsqueda de línea apropiado para que se obtenga una disminución suficiente en una función de mérito (ver).Actualizando la matriz de hessian La matrizH_k es una aproximación definitiva positiva de la matriz Hessiana de la función Lagrangia ().Ecuación 13H_k puede ser actualizado por cualquiera de los métodos quasi-Newton, aunque el método BFGS (ver) parece ser el más popular.Actualizando la matriz de hessian

Un problema no linealmente restringido a menudo se puede resolver en menos iteraciones que un problema sin restricciones usando SQP. Una de las razones de esto es que, debido a los límites en el área factible, el optimizador puede tomar decisiones informadas con respecto a las direcciones de búsqueda y la longitud del paso.

Considerar La función de Rosenbrock con una restricción de desigualdad no lineal adicional, (),gx

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

(15)

Esto se resolvió mediante una implementación de SQP en 96 iteraciones en comparación con 140 para el caso sin restricciones. muestra la ruta de acceso al punto de soluciónFigura 6-3, Método SQP en la función no linealmente restringida Rosenbrock x = [0.9072,0.8228] comenzando en x = [–1.9,2.0].

Figura 6-3, Método SQP en la función no linealmente restringida Rosenbrock

Implementación de SQP

La implementación de SQP consta de tres etapas principales, que se discuten brevemente en las siguientes subsecciones:

• Actualizando la matriz de hessian

• Solución de programación cuadrática

• Inicialización

• Búsqueda de línea y función de mérito

Actualizando la matriz de hessian.  En cada iteración importante, se calcula una aproximación de cuasi-Newton definida de la función de Lagrangia, con el método BFGS, dondeH λ_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)

Dónde

$$\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 recomienda mantener el hessian positivo definitivo a pesar de que podría ser positivo indefinido en el punto de solución.[33] Se mantiene un hessian definitivo positivo que proporciona $$q_k^T{s}_k$$ es positivo en cada actualización y que se inicializa con una matriz definida positiva.H Cuando $$q_k^T{s}_k$$ no es positivo, Q_k se modifica elemento por elemento de forma que $$q_k^T{s}_k>0$$. El objetivo general de esta modificación es distorsionar los elementos de Q_k, que contribuyen a una actualización definitiva positiva, lo más poco posible. Por lo tanto, en la fase inicial de la modificación, el elemento más negativo de q_k*s_k se reduce en repetidas ocasiones. Este procedimiento se continúa hasta $$q_k^T{s}_k$$ es mayor o igual que una pequeña tolerancia negativa. Si, después de este procedimiento, $$q_k^T{s}_k$$ todavía no es positivo, modifique Q_k añadiendo un vector multiplicado por un escalar constante, es decir,vw

$$q_k=q_k+wv,$$

(17)

Dónde

$$\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 aumentar sistemáticamente hastaw $$q_k^T{s}_k$$ se vuelve positivo.

Las funciones, y todos usan SQP.fminconfminimaxfgoalattainfseminf Si se establece en in, a continuación, se da diversa información, como los valores de función y la infracción de restricción máxima.Display'iter'options Cuando el Hessiano tiene que ser modificado utilizando la primera fase del procedimiento anterior para mantenerlo positivo definitivo, entonces Hessian modified se visualiza. Si el Hessiano tiene que ser modificado de nuevo utilizando la segunda fase del enfoque descrito anteriormente, entonces se visualiza.Hessian modified twice Cuando el subproblema de QP es inviable, se visualiza.infeasible Estas pantallas no suelen ser motivo de preocupación, pero indican que el problema es altamente no lineal y que la convergencia puede tardar más de lo habitual. A veces el mensaje no update se muestra, indicando que $$q_k^T{s}_k$$ es casi cero. Esto puede ser una indicación de que la configuración del problema es incorrecta o que está tratando de minimizar una función no continua.

Solución de programación cuadrática.  En cada iteración importante del Método SQP, se resuelve un problema de QP de la siguiente forma, donde Un_i se refiere a la fila TH de la-por-matriz.imnA

$$\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 es una estrategia de conjunto activa (también conocida como método de proyección) similar a la de Gill et al., descrita en y.Optimization Toolbox[18][17] Se ha modificado para los problemas de programación lineal (LP) y de programación cuadrática (QP).

El procedimiento de solución involucra dos fases. La primera fase implica el cálculo de un punto factible (si existe). La segunda fase implica la generación de una secuencia iterativa de puntos factibles que convergen en la solución. En este método, un conjunto activo, $${\overline{A}}_k$$, se mantiene que es una estimación de las restricciones activas (es decir, las que están en los límites de restricción) en el punto de solución. Prácticamente todos los algoritmos de QP son métodos de conjunto activos. Este punto se enfatiza 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, y esto se utiliza para formar una base para una dirección de búsquedak $${\widehat{d}}_k$$. Las restricciones de igualdad siempre permanecen en el conjunto activo $${\overline{A}}_k$$. La notación para la variable $${\widehat{d}}_k$$ se utiliza aquí para distinguirlo de D_k en las iteraciones principales del método SQP. La dirección de búsqueda $${\widehat{d}}_k$$ se calcula y minimiza la función objetiva mientras permanece en los límites de restricción activa. El subespacio factible para $${\widehat{d}}_k$$ se forma a partir de una 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$$). Por lo tanto, 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, se garantiza que permanecerá en los límites de las restricciones activas.

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

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

(19)

Dónde

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

una vez Z_k se encuentra, una nueva dirección de búsqueda $${\widehat{d}}_k$$ se busca que minimiza () dondeqd $${\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 un vector.p

Luego, si ve el cuadrático como una función de, sustituyendo parap $${\widehat{d}}_k$$, usted tiene

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

(20)

Diferenciándolo con respecto a los rendimientosp

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

(21)

∇q(p) se conoce como el degradado proyectado de la función cuadrática porque es el degradado proyectado en el subespacio definido por Z_k. El término $$Z_k^{T}H{Z}_k$$ se llama el hessian proyectado. Suponiendo que la matriz de Hessian es positiva definida (que es el caso en esta implementación de SQP), entonces el mínimo de la función () en el subespacio definido porHqp Z_k se produce 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 toma un paso de la forma

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

(23)

En cada iteración, debido a la naturaleza cuadrática de la función objetiva, sólo hay dos opciones de longitud de paso.α Un paso de unidad junto $${\widehat{d}}_k$$ es el paso exacto al mínimo de la función restringida al espacio nulo de $${\overline{A}}_k$$. Si se puede tomar tal paso, sin infringir las restricciones, esta es la solución a QP ().Ecuación 18 De lo contrario, el paso a lo largo $${\widehat{d}}_k$$ a 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 restricción en cualquier dirección $${\widehat{d}}_k$$ es 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 para las restricciones no en el conjunto activo, y donde la dirección $${\widehat{d}}_k$$ es hacia el límite de restricción, es decir, $$A_i{\widehat{d}}_k>0,i=1,\mathrm{...},m$$.

Cuando se incluyen restricciones independientes en el conjunto activo, sin la ubicación del mínimo, los multiplicadores de Lagrange,n Λ_k, se calculan 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 positivas, X_k es la solución óptima de QP ().Ecuación 18 Sin embargo, si algún componente de Λ_k es negativa y el componente no se corresponde con una restricción de igualdad, el elemento correspondiente se elimina del conjunto activo y se busca una nueva iteración.

Inicialización.  El algoritmo requiere un punto factible para comenzar. Si el punto actual del método SQP no es factible, entonces puede encontrar un punto 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 Un_i indica la fila TH de la matriz.iA Puede encontrar un punto factible (si existe alguno) estableciendo un valor que satisfaga las restricciones de igualdad.Ecuación 26x Puede determinar este valor resolviendo un conjunto de ecuaciones lineales bajo o sobredeterminado formado a partir del conjunto de restricciones de igualdad. Si hay una solución a este problema, la variable de Slack se establece en la restricción de desigualdad máxima en este momento.γ

Puede modificar el algoritmo 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 degradado de la función objetiva (igual a los coeficientes de la función objetiva lineal).

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

(27)

Si se encuentra un punto factible utilizando el método LP anterior, se introduce 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$$ encontrados para resolver el conjunto de ecuaciones lineales

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

(28)

Dónde G_k es el degradado de la función objetiva en la iteración actual X_k (es decir, Hx_k + c).

Si no se encuentra una solución viable para el problema de QP, la dirección de búsqueda de la rutina principal de SQP se toma como una que minimiza.γ

Búsqueda de línea y función de mérito.  La solución al subproblema de QP produce un vector D_k, que se utiliza para formar una nueva iteración

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

(29)

El parámetro de longitud de paso Α_k se determina con el fin de producir una disminución suficiente en un función de mérito. La función de mérito utilizada por han y Powell de la siguiente forma se utiliza en esta implementación.[22][33]

$$\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 establecer 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 han sido recientemente activas. En esta implementación, el parámetro Penalty R_i se establece inicialmente en

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

(32)

Dónde $$\Vert \Vert$$ representa la norma euclidiana.

Esto garantiza mayores contribuciones al parámetro de penalización de las restricciones con degradados más pequeños, que sería el caso de las restricciones activas en el punto de solución.

Algoritmo de fmincon SQP

El algoritmo (y el algoritmo casi idéntico) es similar al algoritmo (para una descripción, consulte).sqpsqp-legacyactive-setfmincon Active Set Algorithm El algoritmo básico se describe en el capítulo 18 de Nocedal y Wright.sqp[31]

El algoritmo es esencialmente el mismo que el algoritmo, pero tiene una implementación diferente.sqpsqp-legacy Normalmente, tiene un tiempo de ejecución más rápido y menos uso de memoria que.sqpsqp-legacy

Las diferencias más importantes entre el y los algoritmos son:sqpactive-set

Viabilidad estricta con respecto a los límites

El algoritmo toma todos los pasos iterativos de la región restringidos por límites.sqp Además, los pasos de diferencia finita también respetan los límites. Los límites no son estrictos; un paso puede estar exactamente en un límite. Esta viabilidad estricta puede ser beneficiosa cuando la función objetiva o las funciones de restricción no lineal no están definidas o son complejas fuera de la región restringida por límites.

Robustez a resultados no dobles

Durante sus iteraciones, el algoritmo puede intentar dar un paso que falla.sqp Esto significa que una función objetiva o una función de restricción no lineal que proporcione devuelve un valor, o un valor complejo.InfNaN En este caso, el algoritmo intenta dar un paso más pequeño.

Las rutinas de álgebra lineal refactorizada

El algoritmo utiliza un conjunto diferente de rutinas de álgebra lineal para resolver el subproblema de programación cuadrática,.sqpEcuación 14 Estas rutinas son más eficientes tanto en el uso de la memoria como en la velocidad que las rutinas.active-set

Las rutinas de viabilidad reformuladas

El algoritmo tiene dos nuevos enfoques para la solución cuando no se satisfacen las restricciones.sqpEcuación 14

• El algoritmo combina las funciones objetivo y restricción en una función de mérito.sqp El algoritmo intenta minimizar la función de mérito sujeto 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 aumenta.Ecuación 14 El aumento del tamaño puede ralentizar la solución del subproblema. Estas rutinas se basan en los artículos de Spellucci y Tone.[60][61] El algoritmo establece el parámetro de penalización para la función de mérito según la sugerencia en.sqpEcuación 30[41]

• Suponga que no se satisfacen las restricciones no lineales y un paso de intento hace que la infracción de restricción crezca. El algoritmo intenta obtener viabilidad utilizando una aproximación de segundo orden a las restricciones.sqp La técnica de segundo orden puede conducir a una solución factible. Sin embargo, esta técnica puede ralentizar la solución al requerir más evaluaciones de las funciones de restricción no lineal.

Algoritmo de punto interior de fmincon

Función de barrera

El enfoque de punto interior para la minimización restringida es resolver una secuencia de problemas aproximados de minimización. El problema original es

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

(33)

μ

$$\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)

variables de Slack s_i ya que hay restricciones de desigualdad.g el s_i se limitan a ser positivos para mantener ln(s_i) Limitada. Como disminuye a cero, el mínimo deμ F_μ debe acercarse al mínimo de.f El término logarítmico añadido se llama a.función de barrera Este método se describe en, y.[40][41][51]

El problema aproximado es una secuencia de problemas de igualdad restringida.Ecuación 34 Estos 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 en (,).Directaxs Este paso intenta resolver las ecuaciones KKT, y, para el problema aproximado a través de una aproximación lineal.Ecuación 2Ecuación 3 Esto también se llama a.Newton paso

• Un paso (degradado conjugada), utilizando una región de confianza.Cg

De forma predeterminada, el algoritmo primero intenta dar un paso directo. Si no puede, intenta un paso CG. Un caso donde no toma un paso directo es cuando el problema aproximado no es convexo localmente cerca de la iteración actual.

En cada iteración, el algoritmo disminuye un Comofunción de mérito

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

El parámetro $$\nu$$ puede aumentar con el número de iteración con el fin de forzar la solución hacia la viabilidad. Si un paso de intento no disminuye la función de mérito, el algoritmo rechaza el paso de intento, e intenta un nuevo paso.

Si la función objetiva o una función de restricción no lineal devuelve un valor complejo, NaN, INF o un error en una iteración X_j, el algoritmo rechaza X_j. El rechazo tiene el mismo efecto que si la función de mérito no disminuye lo suficiente: el algoritmo entonces intenta un paso diferente, más corto. Envolver cualquier código que puede error en-:trycatch

function val = userFcn(x) try     val = ... % code that can error catch     val = NaN; end

El objetivo y las restricciones deben producir valores apropiados () en el punto inicial.double

Paso directo

Las siguientes variables se utilizan para definir el paso directo:

• denota el Hessiano del Lagrangio deH 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)}.$$

  (35)

• J_g denota el jacobiano de la función de restricción.g

• J_h denota el jacobiano de la función de restricción.h

• S = diag(s).

• denota el vector multiplicador de Lagrange asociado con restriccionesλg

• Λ = diag(λ).

• denota el vector multiplicador de Lagrange asociado.yh

• denotan el vector de los del mismo tamaño que.eg

define el paso directoEcuación 36 (Δx, Δs):

$$\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+J_h^{T}y+J_g^T\lambda \\ S\lambda -\mu e\\ h\\ g+s\end{array}\right].$$

(36)

Ecuación 2Ecuación 3

Con el fin de resolver esta ecuación para (Δx, Δs), el algoritmo hace una factorización de LDL de la matriz. (Consulte la página de referencia de la función.)Ejemplo 3 — la estructura de DMATLAB^®ldl Este es el paso más costoso computacionalmente. Un resultado de esta factorización es una determinación de si el hessian proyectado es positivo o no definitivo; Si no, el algoritmo utiliza un paso de degradado conjugada, que se describe en la siguiente sección.

Paso de degradado conjugada

El enfoque de gradiente conjugada para resolver el problema aproximado es similar a otros cálculos de gradiente conjugada.Ecuación 34 En este caso, el algoritmo ajusta ambos y, manteniendo los pantalones positivos.xss El enfoque consiste en minimizar una aproximación cuadrática al problema aproximado en una región de confianza, sujeto a restricciones linearizadas.

En concreto, vamos a denotar el radio de la región de confianza y dejar que otras variables se definan como en.RPaso directo El algoritmo obtiene 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, sujeto a ser positivo.λ Luego da un paso (Δx, Δs) a 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,$$

(37)

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

(38)

Ecuación 38R A continuación, se resuelve con las restricciones que se corresponde con el residuo de la resolución, permanecer dentro de la región de confianza de radio, y mantener estrictamente positivo.Ecuación 37Ecuación 38Rs Para obtener información detallada sobre el algoritmo y la derivación, vea, y.[40][41][51] Para obtener otra descripción de los degradados conjugados, véase.Método de gradiente conjugada precondicionado

Opciones de algoritmo de punto interior

Aquí están los significados y efectos de varias opciones en el algoritmo de punto interior.

• : Cuando se establece en, cada iteración satisface las restricciones enlazadas que ha establecido.HonorBoundstrue Cuando se establece en, el algoritmo puede violar los límites durante las iteraciones intermedias.false

• — Cuando se establece en:HessianApproximation

  • , calcula el Hessiano por una aproximación densa cuasi-Newton.'bfgs'fmincon

  • , calcula el hessian mediante una aproximación de cuasi Newton de memoria limitada y a gran escala.'lbfgs'fmincon

  • , calcula un producto vectorial de tiempos de hessian por diferencias finitas del gradiente (s); otras opciones deben establecerse apropiadamente.'fin-diff-grads'fmincon

• : utiliza el manejador de funciones que especifique para calcular el hessian.HessianFcnfminconHessianFcn Ver.Incluidos los hessianos

• — Dar una función separada para la evaluación de vectores de tiempos de hessian.HessianMultiplyFcn Para obtener más información, consulte y.Incluidos los hessianosFunción de multiplicar de hessian

• : Determina si se intentará o no el paso directo de Newton.SubproblemAlgorithm La configuración predeterminada permite que se intente este tipo de paso.'factorization' El ajuste solo permite los pasos de degradado conjugada.'cg'

Para obtener una lista completa de opciones, consulte en.Interior-Point Algorithmfminconoptions

Algoritmo fminbnd

es un solucionador disponible en cualquier instalación.fminbndMATLAB Resuelve un mínimo local en una dimensión dentro de un intervalo delimitado. No se basa en derivados. En su lugar, utiliza la búsqueda de sección dorada y la interpolación parabólica.

fseminf formulación y algoritmo del problema

fseminf formulación problemática

aborda los problemas de optimización con tipos adicionales de restricciones en comparación con los que se abordan.fseminffmincon La formulación de esfmincon

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

tal que c(x) ≤ 0, ceq(x) = 0, A·x ≤ b, Aeq·x = beq, and l ≤ x ≤ u.

agrega el siguiente conjunto de restricciones semiinfinitas a las que ya se han dado.fseminf Para W_j en un intervalo o rectángulo delimitado de uno o dos dimensiones Me_j, para un vector de funciones continuas K(x, w), las restricciones se

K_j( ,x W_j) ≤ 0 para todos los W_j∈Me_j.

El término "dimensión" de un problema significa la dimensión máxima del conjunto de restricciones:fseminfI 1 si todos los Me_j son intervalos, y 2 si al menos un Me_j es un rectángulo. El tamaño del vector de no entra en este concepto de dimensión.K

La razón por la que esto se llama programación semi-infinita es que hay un número finito de variables (yx W_j), pero un número infinito de restricciones. Esto se debe a que las restricciones están sobre un conjunto de intervalos continuos o rectángulosx Me_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.

Es posible que piense que un problema con un número infinito de restricciones es imposible de resolver. aborda esto reformulando el problema a uno equivalente que tiene dos etapas: una maximización y una minimización.fseminf Las restricciones semi-infinitas 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|,$$

(39)

donde | | es el número de componentes del vector; es decir, el número de funciones de restricción semi-infinitas.KK Para Fixed, se trata de una maximización ordinaria sobre intervalos o rectángulos delimitados.x

simplifica aún más el problema al hacer aproximaciones cuadráticas o cúbicas por tramosfseminf κ_j(x, w_j) a las funciones K_j(x, w_j), para cada uno que el solucionador visita. sólo considera el maxima de la función de interpolaciónxfseminf κ_j(x, w_j)En lugar de K_j(x, w_j)En.Ecuación 39 Esto reduce el problema original, minimizando una función semi-infinitamente restringida, a un problema con un número finito de restricciones.

Puntos de muestreo.  Su función de restricción semi-infinita debe proporcionar un conjunto de puntos de muestreo, puntos utilizados en la realización de las aproximaciones cuadráticas o cúbicas. Para lograrlo, debe contener:

• El espaciado inicial entre los puntos de muestreosw

• Una forma de generar el conjunto de puntos de muestreo dews

El espaciado inicial es una matriz | |-by-2.sK La fila TH representa el espaciado de los puntos de muestreo vecinos para la función de restricciónjs K_j. Si K_j depende de una dimensión unidimensional W_jEstablecer. actualiza la matriz en iteraciones subsiguientes.s(j,2) = 0fseminfs

utiliza la matriz para generar los puntos de muestreo, que luego utiliza para crear la aproximaciónfseminfsw κ_j(x, w_j). El procedimiento para generar desde debe mantener los mismos intervalos o rectángulosws Me_j durante la optimización.

Ejemplo de creación de puntos de muestreo.  Considere un problema con dos restricciones semi-infinitas,K₁ YK₂.K₁ tiene una dimensiónw₁YK₂ tiene dos dimensionesw₂. El código siguiente genera un conjunto de muestreo dew₁ = 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;

especifica como cuando llama por primera vez a la función de restricción.fseminfsNaN La comprobación de esto le permite establecer el intervalo de muestreo inicial.

El código siguiente genera un conjunto de muestreo dew₂ en un cuadrado, con cada componente que va de 1 a 100, inicialmente muestreado con más frecuencia en el primer componente que 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 se pueden simplificar de la siguiente manera:

% 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

reduce esencialmente el problema de la programación semi-infinita a un problema abordado por. toma los siguientes pasos para resolver problemas de programación semi-infinitos:fseminffminconfseminf

1. En el valor actual de, identifica todos losxfseminf W_j,i tal que la interpolación κ_j(x, w_j,i) es un máximo local. (El máximo se refiere a la variación de fijo.)wx

2. toma un paso de iteración en la solución del problema:fseminffmincon

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

   tal que c(x) ≤ 0, ceq(x) = 0, A·x ≤ b, Aeq·x = beq, and l ≤ x ≤ u, donde () se complementa con todos los máximos decx κ_j(x, w_j) tomado todos los w_j∈I_j, que es igual al maxima sobre y deji κ_j(x, w_j,i).

3. comprueba si se cumple algún criterio de detención en el nuevo punto (para detener las iteraciones); Si no, continúa en el paso 4.fseminfx

4. comprueba si la discretización de las restricciones semi-infinitas necesita actualizarse y actualiza los puntos de muestreo apropiadamente.fseminf Esto proporciona una aproximación actualizada κ_j(x, w_j). Luego continúa en el paso 1.