En este problema, el objetivo es encontrar el mejor valor (máximo o mínimo) de una función objetivo seleccionando variables reales que satisfagan un conjunto de restricciones de igualdad y desigualdad.
Un problema general de optimización restringida consiste en seleccionar nn variables de decisión reales x0,x1,…,xn−1x_0,x_1,dots,x_{n-1} de una región factible dada de tal manera que se optimice (minimice o maximice) una función objetivo determinada.
[f(x_0,x_1,dots,x_{n-1}).]
Generalmente dejamos que xx denote el vector de nn variables de decisión reales x0,x1,…,xn−1.x_0,x_1,dots,x_{n-1}. Es decir, x=(x0,x1,…,xn−1)x=(x_0,x_1,dots,x_{n-1}) y escribir el programa no lineal general como:
[begin{aligned}text{ Maximizar }&f(x)\text{ Sujeto a }&g_i(x)leq b_i&&qquad(i=0,1,dots,m-1)\&xinmathbb{R}^nend{aligned}]
donde se da cada una de las funciones de restricción g0g_0 hasta gm−1g_{m-1}, y cada bib_i es una constante (Bradley et al., 1977).
Ésta es sólo una forma posible de escribir el problema. Minimizar f(x)f(x) es equivalente a maximizar −f(x).-f(x). Asimismo, una restricción de igualdad h(x)=bh(x)=b se puede expresar como el par de desigualdades h(x)≤bh(x)leq b y −h(x)≤−b.-h(x)leq-b. Además, al agregar una variable de holgura, cualquier desigualdad puede convertirse en una restricción de igualdad (Bradley et al., 1977).
Este tipo de problemas aparecen en muchas áreas de aplicación, por ejemplo, en entornos comerciales donde una empresa apunta a maximizar las ganancias o minimizar los costos mientras opera con recursos o limitaciones de financiamiento (Cherry, 2016).
Si f(x)f(x) es una función lineal y las restricciones son lineales, entonces el problema se llama problema de programación lineal (LP) (Cherry, 2016).
"El problema se denomina problema de programación no lineal (PNL) si la función objetivo no es lineal y/o la región factible está determinada por restricciones no lineales". (Bradley et al., 1977, pág. 410)
Los supuestos y aproximaciones utilizados en la programación lineal a veces pueden producir un modelo adecuado para los rangos de variables de decisión de interés. En otros casos, sin embargo, el comportamiento no lineal en la función objetivo y/o las restricciones es esencial para formular la aplicación como un programa matemático con precisión (Bradley et al., 1977).
La programación no lineal se refiere a métodos para resolver una PNL. Aunque existen muchos solucionadores maduros para LP, los NLP, especialmente aquellos que involucran no linealidades de orden superior, suelen ser más difíciles de resolver (Cherry, 2016).
Surgen problemas desafiantes de programación no lineal en áreas como el diseño de circuitos electrónicos, la optimización de la cartera financiera, la optimización de la red de gas y el diseño de procesos químicos (Cherry, 2016).
Una forma de abordar problemas de programación no lineal es linealizarlos y utilizar un solucionador LP para obtener una buena aproximación. A uno le gustaría hacerlo por varias razones. Por ejemplo, los modelos lineales suelen ser más rápidos de resolver y pueden ser más estables numéricamente.
Para pasar de la teoría a un modelo Python funcional, recorreremos las siguientes secciones:
Este texto supone cierta familiaridad con la programación matemática. Después de definir las funciones separables, presentamos un ejemplo de maximización no lineal separable y describimos un enfoque basado en aproximaciones lineales por partes (PWL) de la función objetivo no lineal.
A continuación, definimos Conjuntos Ordenados Especiales de Tipo 2 y explicamos cómo apoyan la formulación numérica.
Luego presentamos los antecedentes teóricos en Sobre funciones convexas y cóncavas, proporcionando herramientas utilizadas a lo largo del resto de este trabajo sobre programación no lineal (PNL).
Finalmente, presentamos un procedimiento de implementación de Python que utiliza aproximaciones de Gurobipy y PWL de la función objetivo no lineal, lo que permite a los solucionadores de LP/MIP obtener aproximaciones útiles para NLP bastante grandes.
Funciones separables
El procedimiento de solución de este artículo es para programas separables. Los programas separables son problemas de optimización de la forma:
[begin{aligned}text{ Maximizar }&sumlimits_{j=0}^{n-1}f_j(x_j)\text{ Sujeto a }&sumlimits_{j=0}^{n-1} g_{ij}(x_j)leq 0qquad(i=0,1,dots,m-1)\&xinmathbb{R}^nend{aligned}]
donde se conocen cada una de las funciones fjf_j y gijg_{ij}. Estos problemas se denominan separables porque las variables de decisión aparecen por separado: una en cada función de restricción gijg_{ij} y una en cada función objetivo fjf_j (Bradley et al., 1977).
Ejemplo
Considere el siguiente problema de Bradley et al. (1977) que surge en la selección de cartera:
[begin{aligned}text{ Maximizar }&f(x)=20x_0+16x_1-2x_0^2-x_1^2-(x_0+x_1)^2\text{ Sujeto a }&x_0+x_1leq5\&x_0geq0,x_1geq0.end{aligned}]
Como se indicó, el problema no es separable debido al término (x0+x1)2(x_0+x_1)^2 en la función objetivo. Sin embargo, los problemas no separables frecuentemente pueden reducirse a una forma separable utilizando varios trucos de formulación. En particular, al dejar x2=x0+x1x_2=x_0+x_1 podemos volver a expresar el problema anterior en forma separable como:
[begin{aligned}text{ Maximizar }&f(x)=20x_0+16x_1-2x_0^2-x_1^2-x_2^2\text{ Sujeto a }&x_0+x_1leq5\&x_0+x_1-x_2=0\&x_0geq0,x_1geq0,x_2geq 0.end{alineado}]
Claramente, las restricciones lineales son separables y la función objetivo ahora se puede escribir como f(x)=f0(x0)+f1(x1)+f2(x2),f(x)=f_0(x_0)+f_1(x_1)+f_2(x_2), donde
[begin{aligned}f_0(x_0)&=20x_0-2x_0^2\f_1(x_1)&=16x_1-x_1^2end{aligned}]
y
[f_2(x_2)=-x_2^2.]
Para formar el problema de aproximación, aproximamos cada término no lineal fj(xj)f_j(x_j) mediante una función lineal por partes fj^(xj)hat{f_j}(x_j) utilizando un número específico de puntos de interrupción.
Cada función PWL fj^(xj)hat{f_j}(x_j) se define en un intervalo [l,u][l,u] y consta de segmentos de línea cuyos puntos finales se encuentran en la función no lineal fj(xj),f_j(x_j), comenzando en (l,fj(l))(l,f_j(l)) y terminando en (u,fj(u)).(u,f_j(u)). Representamos las funciones PWL utilizando la siguiente parametrización. Se puede escribir cualquier l≤xj≤ulleq x_jleq u como:
[x_j=sumlimits_{i=0}^{r-1}theta_j^i a_i,quad hat f_j(x_j)=sumlimits_{i=0}^{r-1}theta_j^if_j(a_i),quadsumlimits_{i=0}^{r-1}theta_j ^i=1,quadtheta_j=(theta_j^0,theta_j^1,dotstheta_j^{r-1})inmathbb{R}_{geq0}^r]
para cada j=0,1,2,j=0,1,2, donde las funciones PWL están especificadas por los puntos {(ai,fj(ai)}{(a_i,f_j(a_i)} para i=0,1,…,r−1i=0,1,dots,r-1 (Cherry, 2016). Aquí, utilizamos 4 puntos de interrupción distribuidos uniformemente para aproximar cada uno fj(xj)f_j(x_j) Además, dado que x0+x1≤5,x_0+x_1leq5, x0+x1=x2,x_0+x_1=x_2, y x0,x1,x2≥0x_0,x_1,x_2geq0 tenemos que
[0leq x_0leq5,qquad0leq x_1leq 5,quadtext{ y }quad0leq x_2leq5.]
Por tanto, nuestras aproximaciones no necesitan extenderse más allá de estos límites variables. Por lo tanto la función PWL fj^(xj)hat{f_j}(x_j) estará definida por los puntos (0,fj(0)),(0,f_j(0)), (5/3,fj(5/3)),(5/3,f_j(5/3)), (10/3,fj(10/3))(10/3,f_j(10/3)) y (5,fj(5))(5,f_j(5)) para j=0,1,2.j=0,1,2.
Cualquier x0∈[0,5]∩ℝx_0in[0,5]capmathbb{R} se puede escribir como
[begin{aligned}x_0&=theta_0^0cdot0+theta_0^1cdotfrac{5}{3}+theta_0^2cdotfrac{10}{3 }+theta_0^3cdot5\&=theta_0^1cdotfrac{5}{3}+theta_0^2cdotfrac{10}{3}+theta_0^3cdot 5text{ st }sumlimits_{i=0}^3theta_0^i=1end{aligned}]
y f0(x0)f_0(x_0) se puede aproximar como
[begin{alineado}sombrero f_0(x_0)&=theta_0^0f_0(0)+theta_0^1f_0left(frac{5}{3}right)+theta_0^2f_0left(frac{10}{3}right)+th eta_0^3f_0(5)\&=theta_0^1cdotfrac{250}{9}+theta_0^2cdotfrac{400}{9}+theta_0^3cdot50.end{aligned}]
Se sigue el mismo razonamiento para x1x_1 y x2x_2. Por ejemplo, al evaluar la aproximación en x0=1,5x_0=1,5 se obtiene:
[begin{aligned}hat f_0(1.5)&=frac{1}{10}f_0(0)+frac{9}{10}f_0left(frac{5}{3}right)\&=frac{1}{10}cdot 0+frac{9}{10}cdotfrac{250}{9}\&=25end{aligned}]
desde
[1.5=frac{1}{10}cdot0+frac{9}{10}cdotfrac{5}{3}.]
Estas sumas ponderadas se denominan combinaciones convexas, es decir, combinaciones lineales de puntos con coeficientes no negativos que suman 1. Ahora bien, dado que f0(x0),f_0(x_0), f1(x1),f_1(x_1) y f2(x2)f_2(x_2) son funciones cóncavas, se pueden ignorar restricciones adicionales y el problema se puede resolver como un PL (Bradley et al., 1977). Sin embargo, en general, tal formulación no es suficiente. En particular, todavía necesitamos hacer cumplir la condición de adyacencia: como máximo dos pesos θjitheta_j^i son positivos, y si dos pesos son positivos, entonces son adyacentes, es decir, de la forma θjitheta_j^i y θji+1.theta_j^{i+1}. Se aplica una restricción similar a cada aproximación. Por ejemplo, si los pesos θ00=2/5theta_0^0=2/5 y θ02=3/5theta_0^2=3/5 entonces la aproximación da
[begin{aligned}hat f_0(x_0)&=frac{2}{5}cdot 0+frac{3}{5}cdotfrac{400}{9}=frac{80}{3}\x_0&=frac{2}{5}cdot0+frac{3}{5}cdotfrac{10}{3}=2.end{aligned}]
Por el contrario, en x0=2,x_0=2, la curva de aproximación da
[begin{alineado}sombrero f_0(x_0)&=frac{4}{5}cdotfrac{250}{9}+frac{1}{5}cdotfrac{400}{9}=frac{280}{9}\x_0&=frac{4}{5}cdotfrac{5}{3}+frac{1}{5}cdotfrac{10}{3}=2.end{aligned}]
Se puede hacer cumplir la condición de adyacencia introduciendo variables binarias adicionales en el modelo, como lo formula Cherry (2016). Sin embargo, aprovechamos el objeto de restricción SOS (Conjunto ordenado especial) de Gurobipy.
El programa completo es el siguiente
[begin{alignat}{2}text{ Maximizar }&hat f_0(x_0)+hat f_1(x_1)+hat f_2(x_2)\text{ Sujeto a }&x_0=sumlimits_{i=0}^{r-1}theta_0^ia_i\&x_1=sumlimits_{i=0}^{r-1}theta_1^ia_i\&x_2=sumlimits_{i=0}^{r-1}theta_2^ia_i\&hat f_0=sumlimits_{i=0}^{r-1}theta_0^if_0(a_i)\&hat f_1=sumlimits_{i=0}^{r-1}theta_1^if_1(a_i)\&hat f_2=sumlimits_{i=0}^{r-1}theta_2^if_2(a_i)\&sumlimits_{i=0}^{r-1}theta_j^i=1&&qquad j=0,1,2\&{theta_j^0,theta_j^1,dotstheta_j^{r-1}}text{ tipo SOS 2 }&&qquad j=0,1,2\&x_0,x_1,x_2in[0,5]capmathbb{R}\&hat f_0in[0,50]capmathbb{R}\&hat f_1in[0,55]capmathbb{R}\&hat f_2in [-25,0]capmathbb{R}\&theta_j^iin[0,1]capmathbb{R}&&qquad j=0,1,2,;i=0,1,dots,{r-1}.end{alignat}]
Conjuntos de pedido especial de tipo 2
“Un conjunto de variables {x0,x1,…,xn−1}{x_0,x_1,dots,x_{n-1}} es un conjunto ordenado especial de tipo 2 (SOS tipo 2) si xixj=0x_ix_j=0 siempre que |i−j|≥2,|ij|geq2, es decir, como máximo dos variables en el conjunto pueden ser distintas de cero, y si dos variables son distintas de cero deben ser adyacentes en el conjunto." (de Farías et al., 2000, p. 1)
En nuestra formulación, la condición de adyacencia de las variables θjitheta_j^i coincide exactamente con una condición SOS tipo 2: como máximo dos θjitheta_j^i pueden ser positivos y deben corresponder a puntos de interrupción adyacentes. El uso de restricciones SOS tipo 2 nos permite garantizar la adyacencia directamente en Gurobi, que aplica estrategias de ramificación especializadas que capturan automáticamente la estructura SOS, en lugar de introducir variables binarias adicionales manualmente.
Sobre funciones convexas y cóncavas
Para la siguiente discusión, tomamos una instancia de problema de maximización restringida separable general, como se definió previamente en la sección Introducción:
[begin{aligned}text{ Maximizar }&sumlimits_{j=0}^{n-1}f_j(x_j)\text{ Sujeto a }&sumlimits_{j=0}^{n-1} g_{ij}(x_j)leq 0qquad(i=0,1,dots,m-1)\&xinmathbb{R}^n.end{aligned}]
donde se conocen cada una de las funciones fif_i y gijg_{ij}.
En todo momento, el término aproximación lineal por partes (PWL) se refiere a la interpolante secante obtenida conectando pares de puntos de interrupción ordenados (ak,f(ak))(a_k,f(a_k)) con segmentos de línea.
Convexo y cóncavo se encuentran entre las formas funcionales más esenciales en la optimización matemática. Formalmente, una función f(x)f(x) se llama convexa si, para cada yy y zz y cada 0≤λ≤1,0leqlambdaleq1,
[f[lambda y+(1-lambda)z]leqlambda f(y)+(1-lambda)f(z).]
Esta definición implica que la suma de funciones convexas es convexa y los múltiplos no negativos de funciones convexas también son convexos. Las funciones cóncavas son simplemente el negativo de las funciones convexas, para las cuales se sigue la definición anterior excepto por la dirección inversa de la desigualdad. Por supuesto, algunas funciones no son ni convexas ni cóncavas.
Es fácil ver que las funciones lineales son tanto convexas como cóncavas. Esta propiedad es esencial ya que nos permite optimizar (maximizar o minimizar) funciones objetivo lineales utilizando técnicas computacionalmente eficientes, como el método simplex en programación lineal (Bradley et al., 1977).
Además, una función diferenciable de una sola variable f(x)f(x) es convexa en un intervalo exactamente cuando su derivada f′(x)f^prime(x) no es decreciente, lo que significa que la pendiente no disminuye. Del mismo modo, la diferenciable f(x)f(x) es cóncava en un intervalo exactamente cuando f′(x)f^prime(x) no aumenta, lo que significa que la pendiente no aumenta.
De la definición anterior, se puede concluir trivialmente que las aproximaciones PWL sobreestiman las funciones convexas ya que cada segmento de línea que une dos puntos en su gráfica no se encuentra debajo de la gráfica en ningún punto. De manera similar, las aproximaciones PWL subestiman las funciones cóncavas.
Función objetivo no lineal
Esta noción nos ayuda a explicar por qué podemos omitir las restricciones SOS tipo 2 en un problema como el de la sección Ejemplo. Por concavidad, la curva de aproximación PWL se encuentra por encima del segmento que une dos puntos cualesquiera no adyacentes. Por lo tanto, la maximización seleccionará la curva de aproximación solo con pesos adyacentes. Se aplica un argumento similar si tres o más ponderaciones son positivas (Bradley et al., 1977). Formalmente, supongamos que cada fj(xj)f_j(x_j) es una función cóncava y cada gijg_{ij} conocida es una función lineal en xx. Entonces podemos aplicar el mismo procedimiento que en la sección de ejemplo: dejamos que los puntos de interrupción satisfagan
[l= a_0lt a_1lt dotslt a_{r-1}= u,qquad rgeq 2.]
Para números reales l≤xj≤ulleq x_jleq u. Además, sea θj∈ℝ+rtheta_jinmathbb{R}_+^r tal que,
[sumlimits_{i=0}^{r-1}theta_j^i=1]
y asignar
[x_j=sumlimits_{i=0}^{r-1}theta_j^icdot a_i,qquadhat f_j=sumlimits_{i=0}^{r-1}theta_j^if_j(a_i)]
para cada j=0,1,…,n−1.j=0,1,dots,n-1. Y así, el LP transformado queda como sigue.
[begin{alignat}{2}text{ Maximizar }&sumlimits_{j=0}^{n-1}hat f_j\text{ Sujeto a }&sumlimits_{j=0}^{n-1} g_{ij}(x_j)leq 0&&qquad(i=0,1,dots,m-1)\&x_j=sumlimits_{i=0}^{r-1}theta_j^icdot a_i&&qquad(j=0,1,dots,n-1)\&hat f_j=sumlimits_{i=0}^{r-1}theta_j^if_j(a_i)&&qquad(j=0,1,dots,n-1)\&sum_ {i=0}^{r-1}theta_j^i=1&&qquad(j=0,1,dots,n-1)\&theta_jinmathbb{R}_{geq 0}^r&&qquad(j=0,1,dots,n-1)\&xinmathbb{R}^n.end{alignat}]
Ahora, sea pj(xj)p_j(x_j) la curva PWL que pasa por los puntos (ai,fj(ai))(a_i,f_j(a_i)) es decir, para cada xj∈[ak,ak+1],x_jin[a_k,a_{k+1}], j=0,1,…,n−1:j=0,1,dots,n-1:
[p_j(x_j)=frac{a_{k+1}-x_j}{a_{k+1}-a_k}f_j(a_k)+frac{x_j-a_k}{a_{k+1}-a_k}f_j(a_{k+1})]
fj(xj)f_j(x_j) cóncava implica que sus pendientes secantes no son crecientes, y dado que pj(xj)p_j(x_j) se construye uniendo puntos de interrupción ordenados en la gráfica de fj(xj),f_j(x_j), tenemos que pj(xj)p_j(x_j) también es cóncava. Y así, por la desigualdad de Jensen,
[p_j(x_j)=p_jleft(sumlimits_{i=0}^{r-1}theta_j^icdot a_iright)geqsumlimits_{i=0}^{r-1}theta_j^icdot p_j(a_i)=sumlimits_{i=0}^{r-1}theta_j^if_j(a_i)=hat f_j.]
Además, para todo xjx_j en cualquier solución factible ((x0,f0^,θ0),(x1,f1^,θ1),…,(xn−1,f^n−1,θn−1)),((x_0,hat{f_0},theta_0),(x_1,hat{f_1},theta_1),dots,(x_{n-1},hat{f}_{n-1},theta_{n-1})), se puede elegir kk tal que xj∈[ak,ak+1]x_jin[a_k,a_{k+1}] y definir el vector θj^hat{theta_j} por
[hattheta_j^k=frac{a_{k+1}-x_j}{a_{k+1}-a_k},qquad hattheta_j^{k+1}=frac{x_j-a_k}{a_{k+1}-a_k},qquadhattheta_j^i=0,;(inotin {k,k+1})qquad]
con θ^jk,θ^jk+1≥0hattheta_j^k,hattheta_j^{k+1}geq 0 y θ^jk+θ^jk+1=1.hattheta_j^k+hattheta_j^{k+1}=1. Entonces θ^jk⋅ak+θ^jk+1⋅ak+1=xjhattheta_j^kcdot a_k+hattheta_j^{k+1}cdot a_{k+1}=x_j y da
[hat f_j^text{nuevo}=hattheta_j^k f_j(a_k)+hattheta_j^{k+1}f_j(a_{k+1})=p_j(x_j)geqhat f_j.]
Debido a que las restricciones dependen de θjtheta_j solo a través del valor inducido xj,x_j, reemplazar cada θjtheta_j por θ^jhattheta_j preserva la viabilidad (porque deja cada xjx_j sin cambios) y mejora débilmente la función objetivo (porque aumenta cada f^jhat f_j a pj(xj)p_j(x_j)). Por lo tanto, las restricciones SOS tipo 2 que imponen la condición de adyacencia no cambian el valor óptimo en problemas de maximización separables con funciones objetivo cóncavas. Un argumento similar muestra que las restricciones SOS tipo 2 no cambian el valor óptimo en problemas de minimización separables con funciones objetivo convexas.
Sin embargo, tenga en cuenta que esto muestra que existe una solución óptima en la que sólo los pesos correspondientes a los puntos de interrupción adyacentes son positivos. No garantiza que todas las soluciones óptimas del LP satisfagan la adyacencia. Por lo tanto, aún se pueden incluir restricciones SOS tipo 2 para imponer una representación consistente.
En la práctica, resolver modelos con funciones objetivo cóncavas mediante aproximaciones PWL puede producir un valor objetivo que subestima el verdadero valor óptimo. Por el contrario, la resolución de modelos con funciones objetivo convexas mediante aproximaciones PWL produce un valor objetivo que sobreestima el verdadero valor óptimo.
Sin embargo, aunque estas aproximaciones pueden distorsionar el valor objetivo, en modelos donde las restricciones lineales originales permanecen sin cambios, la viabilidad no se ve afectada por las aproximaciones PWL. Cualquier asignación de variables encontrada al resolver el modelo mediante aproximaciones PWL satisface las mismas restricciones lineales y, por lo tanto, se encuentra en la región factible original. En consecuencia, se puede evaluar el objetivo no lineal original f(x)f(x) en la solución devuelta para obtener el valor objetivo, aunque este valor no tiene por qué ser igual al verdadero óptimo no lineal.
Restricciones no lineales
Lo que requiere enfoques diferentes son aproximaciones PWL de restricciones no lineales, ya que estas aproximaciones pueden cambiar la región factible original. Definimos la región/conjunto factible del problema, como el anterior, por:
[F=bigcaplimits_{i=0}^{m-1}{xinmathbb{R}^n:G_i(x)leq0},]
dónde
[G_i(x):=sumlimits_{j=0}^{n-1}g_{ij}(x_j).]
Ahora supongamos (para simplificar la notación) que todas Gi(x)G_i(x) son funciones no lineales. Aproxima cada gij(xj)g_{ij}(x_j) univariado mediante una función PWL g^ij(xj)hat g_{ij}(x_j) y define:
[hat G_i(x):=sumlimits_{j=0}^{n-1}hat g_{ij}(x_j).]
Entonces
[hat F=bigcaplimits_{i=0}^{m-1}{xinmathbb{R}^n:hat G_i(x)leq0},]
es el conjunto factible de LP transformado que utiliza aproximaciones PWL. Para evitar soluciones verdaderamente inviables, queremos
[hat Fsubseteq F.]
Si, en cambio, F⊆F^Fsubseteqhat F entonces puede existir x∈F^∖Fxin hat Fsetminus F, lo que significa que el modelo que usa aproximaciones PWL podría aceptar puntos que violen las restricciones no lineales originales.
Dada la forma en que construimos las aproximaciones PWL, una condición suficiente para F^⊆Fhat Fsubseteq F es que para restricciones de la forma Gi(x)≤0G_i(x)leq0 cada gij(xj)g_{ij}(x_j) sea convexa y para restricciones de la forma Gi(x)≥0G_i(x)geq0 cada una gij(xj)g_{ij}(x_j) es cóncavo.
Para ver por qué esto es válido, tome cualquier x∈F^.xinhat F. Primero, considere las restricciones de la forma Gi(x)≤0G_i(x)leq0 donde cada gij(xj)g_{ij}(x_j) es convexo. Debido a que las aproximaciones PWL sobreestiman las funciones convexas, para cualquier i=0,1,…,m−1,i=0,1,dots,m-1 fijo, tenemos:
[(forall j,;hat g_{ij}(x_j)geq g_{ij}(x_j))rightarrowsumlimits_{j=0}^{n-1}hat g_{ij}(x_j)geqsumlimits_{j=0}^{n-1}g_{ij}(x_j)rightarrowhat G_i(x)geq G_i(x)]
Dado que x∈F^,xinhat F, satisface G^i(x)≤0.hat G_i(x)leq0. Junto con G^i(x)≥Gi(x),hat G_i(x)geq G_i(x), esto implica Gi(x)≤0.G_i(x)leq0.
A continuación, considere restricciones de la forma Gi(x)≥0G_i(x)geq0 donde cada gij(xj)g_{ij}(x_j) es cóncava. Debido a que las aproximaciones PWL subestiman las funciones cóncavas, para cualquier i=0,1,…,m−1,i=0,1,dots,m-1 fijo, tenemos:
[(forall j,;hat g_{ij}(x_j)leq g_{ij}(x_j))rightarrowsumlimits_{j=0}^{n-1}hat g_{ij}(x_j)leqsumlimits_{j=0}^{n-1}g_{ij}(x_j)rightarrowhat G_i(x)leq G_i(x).]
Nuevamente, dado que x∈F^,xinhat F, satisface G^i(x)≥0.hat G_i(x)geq0. Combinado con G^i(x)≤Gi(x),hat G_i(x)leq G_i(x), esto implica Gi(x)≥0.G_i(x)geq0.
Esto motiva la noción de un conjunto convexo: Un conjunto CC es convexo si, para dos puntos cualesquiera x,y∈Cx,yin C, todo el segmento de recta que los conecta se encuentra dentro de CC. Formalmente, CC es convexo si para todos x,y∈Cx,yin C y cualquier λ∈[0,1]lambdain[0,1] el punto λx+(1−λ)ylambda x+(1-lambda)y también está en CC
En particular, si Gi:ℝn→ℝG_i:mathbb{R}^ntomathbb{R} es convexo, entonces
[S_i:={xinmathbb{R}^n:G_i(x)leq b}]
es convexo para cualquier b∈ℝ.binmathbb{R}. Prueba: Tome cualquier x1,x2∈S,x_1,x_2in S, luego Gi(x1)≤bG_i(x_1)leq b y Gi(x2)≤b.G_i(x_2)leq b. Como Gi(x)G_i(x) es convexo, para cualquier λ∈[0,1]lambdain[0,1] tenemos que
[G_i(lambda x_1+(1-lambda)x_2)leqlambda G_i(x_1)+(1-lambda)G_i(x_2)leqlambda b+(1-lambda)b=b.]
Por lo tanto, λx1+(1−λ)x2∈S,lambda x_1+(1-lambda) x_2in S, y por lo tanto, SS es convexo. Se aplica un argumento similar para demostrar que
[T_i:={xinmathbb{R}^n:G_i(x)geq b}]
es convexo para cualquier b∈ℝbinmathbb{R} si Gi:ℝn→ℝG_i:mathbb{R}^ntomathbb{R} es cóncavo.
Y como se podría demostrar fácilmente que la intersección de cualquier colección de conjuntos convexos también es convexa,
"Para un conjunto factible convexo queremos funciones convexas para restricciones menores o iguales que y funciones cóncavas para restricciones mayores o iguales". (Bradley et al., 1977, p.418)
Los problemas de optimización no lineal en los que el objetivo es minimizar un objetivo convexo (o maximizar uno cóncavo) sobre un conjunto convexo de restricciones se denominan problemas no lineales convexos. Estos problemas son generalmente más fáciles de abordar que los no convexos.
De hecho, si F⊆F^Fsubseteqhat F entonces las asignaciones de variables pueden conducir a violaciones de restricciones. En estos casos, sería necesario confiar en otras técnicas, como la simple modificación MAP del algoritmo de Frank-Wolfe proporcionado en Bradley et al. (1977), y/o introducir tolerancias que, hasta cierto punto, permitan violaciones de restricciones.
Implementación de Python
Antes de presentar el código, es útil considerar cómo escalará nuestro enfoque de modelado a medida que aumente el tamaño del problema. Para un programa separable con nn variables de decisión y rr puntos de interrupción por variable, el modelo introduce n⋅rncdot r variables continuas (para las combinaciones convexas) y rr restricciones SOS tipo 2. Esto significa que tanto el número de variables como las restricciones crecen linealmente con nn y r,r, pero su producto puede conducir rápidamente a modelos grandes cuando se incrementa cualquiera de los parámetros.
Como se mencionó anteriormente, usaremos Gurobipy para resolver numéricamente el programa anterior. Generalmente el procedimiento es el siguiente:
Elegir límites [l,u][l,u] Elegir puntos de interrupción aia_i Agregar restricciones de combinación convexa θ,theta Agregar SOS tipo 2 Resolver Evaluar el objetivo no lineal en el xx devuelto Refinar los puntos de interrupción si es necesario
Después de importar las bibliotecas y dependencias necesarias, declaramos y creamos instancias de nuestras funciones no lineales f0,f_0, f1,f_1 y f2:f_2:
# Importaciones import gurobipy as grb import numpy as np # Funciones no lineales f0=lambda x0: 20.0*x0-2.0*x0**2 f1=lambda x1: 16.0*x1-x1**2 f2=lambda x2: -x2**2
A continuación, elegimos el número de puntos de interrupción. Para construir las aproximaciones PWL, necesitamos distribuirlas en el intervalo [0,5]∩ℝ.[0,5]capmathbb{R}. Existen muchas estrategias para hacer esto. Por ejemplo, se podrían agrupar más puntos cerca de los extremos o incluso confiar en procedimientos adaptativos como en Cherry (2016). Por simplicidad, nos atenemos al muestreo uniforme:
lb,ub=0.0,5.0 r=4 # Número de puntos de interrupción a_i=np.linspace(start=lb,stop=ub,num=r) # Distribuir puntos de interrupción uniformemente
Luego, podemos calcular el valor de las funciones fj(xj)f_j(x_j) en los puntos de interrupción elegidos para aproximar las variables de decisión:
# Calcular los valores de la función en los puntos de interrupción f0_hat_r=f0(a_i) f1_hat_r=f1(a_i) f2_hat_r=f2(a_i)
El siguiente paso es declarar y crear una instancia de un modelo de optimización de Gurobipy:
modelo=gp.Modelo()
Después de crear nuestro modelo, necesitamos agregar las variables de decisión. Gurobipy proporciona métodos para agregar múltiples variables de decisión a un modelo a la vez, y aprovechamos esta característica para declarar y crear instancias de algunas de nuestras variables de decisión, es decir, los coeficientes de las aproximaciones PWL. Tenga en cuenta que todas las variables de decisión en nuestro modelo son continuas y pueden asignarse a cualquier valor real dentro de los límites proporcionados. Estos modelos se pueden resolver de manera eficiente utilizando paquetes de software modernos como Gurobi.
# Variables de decisión x0=model.addVar(lb=0.0,ub=5.0,name='x0',vtype=GRB.CONTINUOUS) x1=model.addVar(lb=0.0,ub=5.0,name='x1',vtype=GRB.CONTINUOUS) x2=model.addVar(lb=0.0,ub=5.0,nombre='x2',vtype=GRB.CONTINUOUS) f0_hat=model.addVar(lb=0.0,ub=50.0,name='f0_hat',vtype=GRB.CONTINUOUS) f1_hat=model.addVar(lb=0.0,ub=55.0,nombre='f1_hat',vtype=GRB.CONTINUOUS) f2_hat=model.addVar(lb=-25.0,ub=0.0,name='f2_hat',vtype=GRB.CONTINUOUS) theta=model.addVars(rango(0,3),rango(0,r),lb=0.0,ub=1.0,nombre='theta',vtype=GRB.CONTINUOUS)
Además, necesitamos agregar restricciones que limiten los valores que pueden tomar nuestras variables de decisión. En primer lugar, añadimos las restricciones que aseguran la decisión correcta de las variables que dan lugar a las aproximaciones de las funciones no lineales:
# Restricciones model.addConstr(x0==gp.quicksum(theta[0,i]*a_i[i] for i in range(0,r))) model.addConstr(x1==gp.quicksum(theta[1,i]*a_i[i] for i in range(0,r))) model.addConstr(x2==gp.quicksum(theta[2,i]*a_i[i] para i en el rango(0,r))) model.addConstr(f0_hat==gp.quicksum(theta[0,i]*f0_hat_r[i] para i en el rango(0,r))) model.addConstr(f1_hat==gp.quicksum(theta[1,i]*f1_hat_r[i] para i en el rango(0,r))) model.addConstr(f2_hat==gp.quicksum(theta[2,i]*f2_hat_r[i] para i en el rango(0,r))) model.addConstr(gp.quicksum(theta[0,i] para i en el rango(0,r))==1) model.addConstr(gp.quicksum(theta[1,i] para i en el rango(0,r))==1) model.addConstr(gp.quicksum(theta[2,i] para i en el rango(0,r))==1)
En segundo lugar, agregamos las restricciones que definen el problema original, es decir, aquellas restricciones en la formulación no lineal:
model.addConstr(x0+x1<=5.0) model.addConstr(x0+x1-x2==0.0)
Y finalmente, agregamos las restricciones SOS Tipo 2 para asegurar la condición de adyacencia:
# Conjuntos ordenados especiales model.addSOS(GRB.SOS_TYPE2,[theta[0,i] for i in range(0,r)]) model.addSOS(GRB.SOS_TYPE2,[theta[1,i] for i in range(0,r)]) model.addSOS(GRB.SOS_TYPE2,[theta[2,i] for i in range(0,r)])
Sólo queda definir la función objetivo que nos gustaría maximizar:
model.setObjective(f0_hat+f1_hat+f2_hat,GRB.MAXIMIZE) # Función objetivo
El modelo ahora está completo y podemos usar el solucionador Gurobi para recuperar el valor óptimo y las asignaciones de variables óptimas. El código para hacer esto se muestra a continuación:
model.optimize() # ¡Optimizar! obj_val=model.ObjVal res={} all_vars=model.getVars() valores=model.getAttr('X',all_vars) nombres=model.getAttr('VarName',all_vars) # Recuperar valores de variables después de optimizar para nombre,val en zip(nombres,valores): res[nombre]=val
Después de ejecutar el código anterior e imprimir atributos en la pantalla, el resultado es el siguiente:
– Estado: 2 – Número de puntos de interrupción: 4 – Valor objetivo óptimo: 45,00 – Valor óptimo x0, f0(x0) (aproximación PWL): 1,67, 27,78 – Valor óptimo x1, f1(x1) (aproximación PWL): 3,33, 42,22 – Valor óptimo x2, f2(x2) (aproximación PWL): 5,00, -25.00
La diferencia entre nuestra solución calculada y el valor óptimo de 139/3 es 4/3. Ahora, para aumentar la precisión y obtener una mejor solución, se podrían agregar más puntos de interrupción y, en consecuencia, aumentar el número de segmentos utilizados para aproximar las funciones no lineales. Por ejemplo, el siguiente resultado muestra cómo aumentar el número de puntos de interrupción a 16 conduce a una aproximación prácticamente sin error al verdadero valor óptimo:
– Estado: 2 – Número de puntos de interrupción: 16 – Valor objetivo óptimo: 46,33 – Valor óptimo x0, f0(x0) (aproximación PWL): 2,33, 35,78 – Valor óptimo x1, f1(x1) (aproximación PWL): 2,67, 35,56 – Valor óptimo x2, f2(x2) (aproximación PWL): 5,00, -25.00
Puede encontrar el código Python completo en este repositorio de GitHub: enlace.
Lecturas adicionales
Para un algoritmo más avanzado y robusto que utiliza funciones PWL para aproximar la función objetivo no lineal, Cherry (2016) es una referencia valiosa. Presenta un procedimiento que comienza con aproximaciones gruesas, que se refinan iterativamente utilizando la solución óptima a la relajación LP del paso anterior.
Para una comprensión más profunda de la programación no lineal, incluidos métodos, ejemplos aplicados y una explicación de por qué los problemas no lineales son intrínsecamente más difíciles de resolver, el lector puede consultar Programación no lineal: programación matemática aplicada (Bradley et al., 1977).
Finalmente, de Farías et al. (2000) presentan experimentos computacionales que utilizan restricciones SOS dentro de un esquema de ramificación y corte para resolver un problema de asignación generalizado que surge en la programación de la producción de cables de fibra óptica.
Conclusión
La programación no lineal es un área relevante en la optimización numérica, con aplicaciones que van desde la optimización de carteras financieras hasta el control de procesos químicos. Este artículo presentó un procedimiento para abordar programas no lineales separables reemplazando términos no lineales con aproximaciones PWL y resolviendo el modelo resultante usando solucionadores LP/MIP. También proporcionamos antecedentes teóricos sobre funciones convexas y cóncavas y explicamos cómo estas nociones se relacionan con la programación no lineal. Con estos elementos, un lector debería poder modelar y resolver aproximaciones compatibles con LP/MIP de muchas instancias de programación no lineal sin depender de solucionadores no lineales dedicados.
¡Gracias por leer!
Referencias
Cherry, A., 2016. Aproximación lineal por partes para problemas de programación no lineal. Una disertación. Universidad Tecnológica de Texas. Disponible en: https://ttu-ir.tdl.org/server/api/core/bitstreams/2ffce0e6-87e9-4e4c-b41a-60ea670c5a13/content (Consulta: 27 de enero de 2026).
Bradley, SP, Hax, AC y Magnanti, TL, 1977. Programación no lineal. En: Programación Matemática Aplicada. Reading, MA: Addison-Wesley Publishing Company, págs. 410–464. Disponible en: https://web.mit.edu/15.053/www/AMP-Chapter-13.pdf (Consulta: 27 de enero de 2026).
de Farias, IR, Johnson, EL y Nemhauser, GL, 2000. Un problema de asignación generalizada con conjuntos ordenados especiales: un enfoque poliédrico. Programación matemática, 89 (1), págs. 187-203. Disponible en: https://doi.org/10.1007/PL00011392 (Consulta: 27 de enero de 2026).