1. Introducción

El Problema de Complementariedad No Lineal, asociado a F: ℝⁿ → ℝⁿ, F(x) = (F₁(x),...,F_n (x))^T , no lineal y continuamente diferenciable, (PCNL(F)), consiste en encontrar un vector x Є ℝⁿ tal que

La tercera condición en (l) exige que x y F(x) sean ortogonales, y es la que le da el calificativo de "complementariedad" al problema. Son numerosas sus aplicaciones en ingeniería, economía y física [1], [2], [3], [4]. Su importancia se debe a que el concepto de complementariedad es sinónimo de sistema en equilibrio.

La estrategia más popular para resolver el PCNL(F) consiste en reformularlo como un sistema de ecuaciones no lineales mediante un operador Ф: ℝⁿ → ℝⁿ, definido por Ф(х) = (φ(x₁,F₁(x)), ..., φ (x_n,F_n(x)))^T , con φ: ℝ² → ℝ, llamada función de complementariedad, que satisface la equivalencia: φ(a, b) = 0 ⇔ α, ≥ 0,b ≥ 0,ab=0 la cual permite demostrar que resolver el PCNL(F) es equivalente a resolver el sistema de ecuaciones no lineales [5] :

que resulta ser no diferenciable debido a la no suavidad de la función de complementariedad φ. El sistema (2), se conoce como reformulación del PCNL(F) como un sistema de ecuaciones no lineales.

Existen numerosas funciones de complementariedad [6], [7], [8] pero, quizá, las más populares debido a su amplio uso en problemas de complementariedad son las funciones Mínimo [9] y Fischer [10] definidas, respectivamente, por: φ(a,b) = min{a,b} y φ(a,b) =

Para resolver el sistema (2) e indirectamente el PCNL(F), se han propuesto algoritmos tipo Newton generalizados [11], cuasi-Newton generalizados [12], [13] y recientemente, algoritmos suavizados y con jacobiano suavizado [14], [15]. Estos últimos como alternativa a la no diferenciabilidad de la reformulación del PCNL(F).

La idea básica de los métodos suavizados es aproximar Ф por Ф_μ ℝⁿ → ℝⁿ ,μ >0,definido por Ф_μ(x)=( φ_μ (x₁, F₁(x)),… φ_μ (x_n, F_n(x))^T donde φ_μ es una suavización de la función de complementariedad (p, y así, aproximar el sistema de ecuaciones no lineales(2) por una sucesión de sistemas suavizados Ф_μ (𝐱)=0 y hacer tender el parámetro 𝜇 a cero. Estos métodos han resultado muy eficientes, en particular los llamados en [14] de jacobiano suavizado que se caracterizan por resolver, en cada iteración, un sistema de ecuaciones lineales suavizado de la forma:

En lugar del sistema Newtoniano,

Con 𝐻_𝑘∈𝜕𝛷(𝐱_𝑘) o 𝐻_𝑘∈𝜕_𝐶𝛷(𝐱_𝑘), donde 𝜕𝛷(𝐱) es el jacobiano generalizado de 𝛷 en 𝐱, definido para 𝛷 Lipschitz continua como el siguiente conjunto [16], ∂Φ(𝒙)=𝑐𝑜𝑛𝑣{𝐻=lim𝑘→∞Φ′(𝒙_𝑘)∈ℝ^𝑛×𝑛: 𝒙_𝑘→𝒙, 𝒙_𝑘∈𝐷_Φ}, donde 𝐷𝛷 es el conjunto de puntos de ℝ𝑛 en los que 𝛷 es diferenciable y 𝑐𝑜𝑛𝑣{𝐴} es la envolvente convexa del conjunto A. En general, 𝜕𝛷(𝐱) es difícil de calcular. Una alternativa basada en el hecho de que 𝜕𝛷(𝐱)𝑇⊆𝜕𝛷1(𝐱)×…×𝜕𝛷𝑛(𝐱) [16], es usar el C-subdiferencial [9] de 𝛷 en 𝐱 definido por, 𝜕𝐶𝛷(𝐱)𝑇=𝜕𝛷1(𝐱)×…×𝜕𝛷𝑛(𝐱). Es decir, el conjunto de matrices en ℝ𝑛×𝑛 cuya 𝑖-ésima columna es el gradiente generalizado de la 𝑖-ésima función componente 𝛷𝑖.

A partir de (3) y (4), se infiere que los métodos con jacobiano suavizado son cuasi-Newton suavizados. Por ello, en este documento se dirá métodos con jacobiano suavizado o cuasi-Newton suavizados.

Debido a la importancia y efectividad de los métodos cuasi-Newton suavizados para resolver indirectamente el PCNL(F), y que la función Mínimo no ha sido usada en conexión con dichos métodos, en el presente trabajo se propone un algoritmo de este tipo, en el cual se usa la función Mínimo y una suavización de la misma propuesta en [17]. Se demuestra, bajo ciertas hipótesis, que el algoritmo propuesto converge local y q-cuadráticamente. Además, se presentan pruebas numéricas que muestran el buen desempeño del algoritmo propuesto.

Este documento se organiza de la siguiente forma. En la Sección 2 se presentan algunos resultados teóricos importantes acerca de la función Mínimo y su suavización; también, se deducen cotas útiles para el desarrollo teórico presentado en la Sección 4. En la Sección 3 se reformula el PCNL(F) como un sistema de ecuaciones no lineales usando la función Mínimo y se c aracteriza la matriz jacobiana de Ф_μ en un punto x. En la primera parte de la Sección 4, se propone un algoritmo cuasi-Newton suavizado para resolver el PCNL(F) y luego, bajo ciertas hipótesis se desarrolla su teoría de convergencia local. En la Sección 5 se analiza numéricamente el desempeño local del algoritmo propuesto. Finalmente, en la Sección 6, se presentan algunos comentarios finales y propuestas de trabajos futuros.

2. La función Mínimo y su suavización

En esta sección se analizan algunas propiedades de la función Mínimo y de una suavización de la misma propuesta en [17]. Además, se deducen algunas cotas útiles en el desarrollo teórico posterior. Para mayor claridad en la lectura del capítulo se incluyen a continuación, a manera de definición, la función Mínimo y el concepto de suavización de una función de complementariedad.

Definición 1. La función 𝜑:ℝ2→ℝ, definida por

se conoce como función Mínimo.

La función Mínimo se puede escribir, equivalentemente, como:

que permite observar su no diferenciabilidad en los puntos de la forma (α, α) como lo ilustra la Figura 1.

Definición 2. Una suavización de una función de complementariedad 𝜑, es una función 𝜑_𝜇:ℝ²→ℝ tal que

A continuación se presenta la suavización 𝜑_𝜇 propuesta en [17] y conocida como función suavizada CHKS (Chen-Harker-Kanzow-Smale).

Definición 3. La función 𝜑𝜇:ℝ²→ℝ, definida por

es una suavización de la función Mínimo.

Un caso particular de esta suavización, se ilustra en la Figura 2, para el parámetro μ = 0.5, en ella se observa la suavidad de dicha superficie, en contraste con la no suavidad de la Figura 1.

El siguiente lema garantiza que 𝜑_𝜇, es en efecto una suavización de la función mínimo.

Lema 1. La función 𝜑𝜇, es una suavización de la función Mínimo.

Demostración. Se debe verificar que la suavización 𝜑𝜇 satisface la equivalencia (7).

de donde

Elevando al cuadrado ambos lados de (9) y organizando términos, se tiene que

con lo cual ab = μ, y como μ > 0, α y b deben tener el mismo signo. De la desigualdad (9) se sigue que a ≥ 0,Ъ≥0.

(ii) Suponga que a ≥ 0,b ≥ 0,ab = μ.

con lo cual, queda demostrado el lema.

Con el fin de simplificar notación en el análisis de resultados posteriores se definen dos nuevas funciones G y G_μ, que hacen parte de las funciones φ y φ_μ, respectivamente.

Definición 4. La función G: ℝ²→ℝ está definida por,

para todo (a,b) ∈ℝ².

Definición 5. La función G_μ: ℝ²→ℝ está definida por

Observación 1. Note que (a - b)²≤ (a - b)²+ 4μ, para todo (a,b) ∈ℝ², de ahí que

El siguiente lema proporciona una cota superior, en términos de μ, para la distancia entre la función <p y su suavización φ_μ.

Lema 2. Las funciones φ y φ_μ satisfacen la siguiente desigualdad

para todo (a,b) ∈ℝ²,

Usando las Definiciones 4 y 5, se tiene que

la tercera igualdad sigue de la Observación 1.

Por lo tanto, |𝜑𝜇(𝑎,𝑏)−𝜑(𝑎,𝑏)|≤√𝜇.

Corolario 1. Las funciones G y G_μ satisfacen la siguiente desigualdad para todo (a,b) ∈ℝ².

Corolario 2. Cuando μ tiende a 0, φμ tiende a φ.

Como se mencionó anteriormente, el propósito de las suavizaciones es lidiar con la no diferenciabilidad de las funciones de complementariedad. Recuerde que, en el caso particular de la función Mínimo su no diferenciabilidad se presenta en los puntos de la forma (α, α), lo cual se soluciona con la suavización φμ definida en (8). Esto se formaliza en el siguiente lema.

Lema 3. La función φμ definida por (8), para todo (a,b) ∈ℝ², es continuamente diferenciable en los puntos de la forma (α, α).

Demostración. Cálculos sencillos permiten verificar las derivadas parciales de φμ en (a,b) existen y son continuas lo que implica el resultado del lema.

Dado que, la función φμ es una función continuamente diferenciable en ℝ², se tiene que para todo (a,b) ∈ℝ², su vector gradiente está definido por:

Notemos que para G_μ definida en (11),

Observación 2. Para todo (a, b) ∈ℝ², la norma euclidiana del gradiente de G_μ en (a,b) satisface la siguiente desigualdad

Lo anterior se tiene usando la Observación 1 en la siguiente igualdad

De (12) y (13)

A partir de (15), usando la desigualdad triangular y (14), se puede calcular una cota superior para la norma euclidiana del gradiente de φμ en (a,b), lo cual se garantiza en la siguiente observación.

Observación 3. Para todo (a,b) ∈ℝ² la norma euclidiana del gradiente de φμ en (a, b) satisface la siguiente desigualdad

Para su uso posterior, las derivadas parciales de G y de G_μ, se denotarán por:

El siguiente lema da una cota superior para 𝐺_𝜇(𝑎,𝑏) la cual será muy útil en la demostración del Lema 6.

Lema 4. Para todo (𝑎,𝑏)∈ℝ², se satisface que 𝐺𝜇(𝑎,𝑏)≥2√𝜇.

Demostración. Para todo (𝑎,𝑏)∈ℝ², se tiene que

por lo tanto 𝐺𝜇(𝑎,𝑏)≥2√𝜇.

Por el Lema 4 se deduce que, para todo (𝑎,𝑏)∈ℝ²

Por otra parte,

Donde

es una matriz simétrica y semidefinida positiva con valores propios 𝜆₁=0 y 𝜆₂=2. Luego,

Lema 5. La función φ_μ es Lipschitz continua con

constante √2. Es decir, para todo 𝐱,𝐲 ∈ℝ², se satisface la desigualdad

Demostración. Sean x,y ∈ℝ². Por el Teorema del Valor Medio existe un vector z en el segmento [x, y] tal que

Usando la desigualdad de Cauchy-Schwartz y la cota dada en (16), se tiene que

con lo cual queda demostrado el lema.

Corolario 3. La función G_μ es Lipschitz continua con constante √2. Es decir, para todo x,y ∈ℝ², se satisface

Demostración. La prueba es análoga a la del lema anterior.

Usando (19), para (𝑎,𝑏) con 𝑎≠𝑏, se tiene que

Lema 6. Sea φ𝜇 la función definida en (7), 𝒘∈ℝ², con 0<‖𝒘‖2≤𝑟, y ℬ(𝒘;𝑟) una bola con centro en 𝐰 y radio 𝑟. Entonces existe 𝜅>0 tal que para todo 𝒖,𝒗∈ℬ(𝒘;𝑟), se satisface que

Demostración. Sean 𝒖 y 𝒗∈ℝ². De (15), (20) y (21), se tiene

Sumando y restando la expresión se deduce que

Para el primer sumando de la anterior desigualdad, y puesto que, por Corolario 3, G_μ es Lipschitz continua con constante √2, se tiene que

Por lo anterior y el Lema 4,

Por otro lado, 𝒗∈ℬ(𝒘;𝑟), ‖𝒗−𝒘‖2<𝜀, de donde ‖𝒗‖2<𝑟+‖𝒘‖2. Dado que ‖𝒘‖2<𝑟, se tiene que ‖𝒗‖2<2𝑟. Finalmente, para 𝒖,𝒗∈ℬ(𝒘;𝑟), se concluye que

luego

3. Reformulación usando la función Mínimo

Si en (2), φ es la función Mínimo (5), se tiene la siguiente reformulación del PCNL(F),

la cual se usa en el desarrollo de este trabajo. Por ello, a partir de ahora, cada vez que se haga referencia a la reformulación Ф = 0, se estará pensando en (22).

Observe que, si existe i Є {1,..., n) tal que x_i = F_i (x)la función es no diferenciable, por lo cual el sistema de ecuaciones no lineales (22) será no diferenciable. El conjunto de estos índices se denotará β. Es decir, β = {i:x_i=F_i(x)}.

Dado que por el Lema 5 la función Mínimo es Lipschitz continua, Ф también lo es, y por tanto, su jacobiano generalizado existe. Formalizando los cálculos presentados en (17), se tiene el siguiente Lema.

Lema 7. Para todo (a, b) ∈ℝ² con а ≠ b,

Demostración. El resultado sigue inmediatamente de

(17).

Corolario 4. Sea G la función definida por (10). Para todo (a, b) ∈ℝ² con а≠Ь, se satisface que

Demostración, Sigue directamente de (23) y de la definición de norma 2.

El siguiente lema permite caracterizar las matrices en 𝜕_𝐶𝛷(x),

Lema 8. Dado x ∈ℝ², H e 𝜕_𝐶𝛷(x), es de la forma

donde la fila i-ésima está dada por

donde c₁ y c₂ son constantes no negativas tales que C₁ + C₂ = 1.

Demostración. A partir de la Proposición 2.6.2 (e) en [16], tenemos que

donde Ф_i denota la i-ésima función componente de Ф.

Si (x_i,F_i(x)) es tal que 𝑥_𝑖≠𝐹_𝑖(x), 𝛷 continuamente diferenciable, por lo tanto,

donde e_i es el i-ésimo vector canónico de ℝⁿ. Para x_i≠ F_i (x), tenemos que

• Si x_i > F_i(x) entonces ∇ Ф_i (x) = ∇F_i (x).

• Si, x_i < F_i(x) entonces ∇ Ф_i (x) = e_i

Por otra parte, si (x_i ,F_i (x)) es tal que x_i = Fi(x), cada H_i∈∂Φ_i(x) es de la forma ∂Φ𝑖(𝑥) = conv[H]i =, cuando к → ∞.

Ahora, denotando la componente i-ésima del vector y_k, por y^k_i, se tiene que

Por el Lema 7, (𝜉𝑖,𝜏𝑖)=(1,−1) o (𝜉𝑖,𝜏𝑖)=(−1,1). En el primer caso, [𝐻]𝑖=∇𝐹𝑖(x)y para el segundo, [𝐻]𝑖=𝒆𝑖. Luego, para (𝑥𝑖,𝐹𝑖(x)), con 𝑥𝑖=𝐹𝑖(x), se tiene 𝜕𝛷𝑖(x),𝑐𝑜𝑛𝑣 {𝒆𝑖∇𝐹𝑖(x)} definido por:

Por lo tanto, una matriz 𝐻∈∂_𝐶Φ(x)es de la forma:

donde

donde 𝑐₁ y 𝑐₂ son constantes no negativas tales que 𝑐₁+𝑐₂=1.

Sin perder generalidad, en lo que sigue, se considera

Es decir, se escoge 𝑐1=1 y 𝑐₂=0 o 𝑐₁=0 y 𝑐₂=1. Esta elección es natural ya que en el caso 𝑥_𝑖=𝐹_𝑖(x),se tiene que 𝑚í𝑛 {𝑥_𝑖,𝐹_𝑖(x)}=𝑥_𝑖=𝐹_𝑖(x). Por ello, es razonable escoger [𝐻]𝑖=𝐞_𝑖^𝑇 o [𝐻]𝑖=∇𝐹_𝑖(x)^𝑇.

Por otra parte, asociada a una solución x∗definimos 𝐻∗∈∂𝐶Φ(x∗) por:

donde

Para finalizar la sección, se deduce una forma explícita de la matriz jacobiana de la suavización 𝛷𝜇 en x. A partir de (8), se define 𝛷𝜇:ℝ^𝑛→ℝ^𝑛 de la siguiente forma

y con ello, el sistema ф_𝜇(𝒙)=𝟎 es una suavización de (22). En general, la matriz jacobiana de ф_𝜇 está dada por,

Para cualquier (𝑎,𝑏)∈ℝ² de (12) y (18) se tiene que

con 𝜒_𝜇(𝑎,𝑏) y 𝜓_𝜇(𝑎,𝑏) dadas por (18).

Para obtener una forma explícita de la matriz (28), se definen las variables auxiliares 𝑢 y 𝑣 por 𝑢=𝑥₁ y 𝑣=𝐹₁(𝒙). Usando la regla de la cadena, se tiene que

De ahí que,

Luego, la primera fila de matriz Φ_𝜇′(𝒙) está dada por:

y de forma análoga, la fila i -ésima de la matriz Φ_𝜇′(x)está dada por:

Por lo tanto, una expresión explícita de la matriz jacobiana Φ_𝜇′(x) es la siguiente

donde 𝜒_𝜇, 𝜓_𝜇 y 𝐺_𝜇 están definidas en (18) y (11),respectivamente.

3.1. Algoritmo y resultados de convergencia

En esta sección, se presenta un nuevo algoritmo tipo cuasi-Newton suavizado para resolver el problema de complementariedad no lineal, indirectamente a través de su reformulación como un sistema de ecuaciones no lineales que usa la función de complementariedad Mínimo, el cual, a su vez, se aproxima por una sucesión de sistemas suavizados que usan una suavización de dicha función. Siguiendo la filosofía de los métodos con jacobiano suavizado [14], una iteración básica del método que se propone es de la forma:

donde de Φ_𝜇′ (х_к) es la matriz jacobiana de Φ_𝜇 en x_k , dada en (30).

A diferencia del método de Newton generalizado, se usa la matriz 𝛷𝜇′(𝒙𝑘) en lugar de una matriz 𝐻_𝑘∈∂_𝐶Φ(𝒙_𝑘). Desde este punto de vista, la iteración (31), puede verse como una iteración cuasi-Newton, donde la matriz Φ′ μ(𝒙_𝑘) es una aproximación de 𝐻_𝑘∈∂_𝐶Φ(𝒙_𝑘), o 𝐻_𝑘∈∂Φ(𝒙_𝑘).

El siguiente lema garantiza que cuando el parámetro μ se acerca a cero, la distancia entre la matriz 𝛷𝜇′(𝒙_𝒌)y el conjunto ∂𝐶Φ(𝒙_𝒌), también tiende a cero, lo cual justifica reemplazar la iteración del método de Newton generalizado por (31).

Lema 9. Sean x ∈ ℝ^𝑛 arbitrario pero fijo 𝑦 μ>0 entonces

Demostración, Puesto que

para que (32) se satisfaga, basta demostrar que existe 𝐻∈∂𝐶Φ(𝒙), la cual es el límite de Φ_𝜇′(𝒙) cuando 𝜇→0. Para ello, se considera el conjunto de índices 𝛽=𝛽(𝒙)={𝑖 :𝑥_𝑖=𝐹_𝑖(𝒙)=0},y usando (29), se tiene que[Φ𝜇′(𝒙)]𝑖está definido por:

equivalentemente

Note que la matriz 𝐻 tiene la estructura que describe el Lema 8, donde los escalares de la combinación lineal son iguales a para 𝑖∈𝛽(𝒙).Por lo tanto, existe H∈𝜕_𝐶𝛷(𝒙),la cual es el límite de 𝛷𝜇′(𝒙)cuando 𝜇→0, de manera que el ínfimo buscado es cero y satisface (32).

Las hipótesis para la teoría de convergencia del algoritmo propuesto son las siguientes.

H1. Existe 𝒙^∗∈ℝ^𝑛 tal que 𝛷(𝒙^∗)=𝟎.

H2. La matriz jacobiana de 𝐹 es Lipschitz continua en una vecindad de 𝒙^∗∈ℝ^𝑛.Es decir, existen constantes 𝛾>0 y 𝛿>0 tales que ∥𝐹′(𝒙)−𝐹′(𝒙∗)∥≤𝛾 ∥𝒙−𝒙^∗∥, para todo 𝒙∈ℬ(𝒙^∗;𝛿) donde ‖ ⋅ ‖ denota una norma en ℝ𝑛 y también su respectiva norma matricial inducida.

H3. Las matrices de 𝜕_𝐶𝛷(𝒙^∗).son no singulares.

Por H3 y la compacidad de ∂Φ(x^∗), existe una constante positiva σ tal que, para toda 𝐻∗∈𝜕𝐶𝛷(𝒙∗) se satisface que ‖𝐻∗−1‖≤𝜎.

Observación 4. Las matrices en ∂_𝐶Φ(𝒙^∗) forman un conjunto finito de 2^𝑟 matrices.

En efecto, suponga que 𝒙∗ es una solución degenerada del PCNL(𝐹), la cual tiene 𝑟 degeneraciones, que sinpérdida de generalidad, se presentan en las primeras 𝑟 componentes de 𝒙^∗ es decir,

Note que por cada degeneración en 𝒙∗ se obtienen dos matrices en ∂_𝐶Φ(𝑥^∗)entonces existirán 2^𝑟 matrices en dicho conjunto.

El siguiente lema da una cota superior para la distancia, en norma infinito, entre las matrices 𝐻 y 𝐻_∗.

Lema10. Sean 𝐹:ℝ^𝑛→ℝ^𝑛 una función continuamente diferenciable cuya matriz jacobiana satisface H2, 𝜖>0 dado, 𝐻∈∂Φ(x )definida por (24) y𝐻∗∈𝜕_𝐶𝛷(𝒙^∗), definida por (26). Entonces para cada 𝒙∈𝐵(𝒙^∗;𝜖), existe una constante 𝛼 tal que

Por definición de norma infinito:

para algún 𝑗∈{1,…,𝑛}. Por (24) y ([26]), se tiene que

Por tanto, existe 𝛼=𝛾𝜖, tal que ‖𝐻∗−𝐻‖_∞<𝛼.

El siguiente es un lema técnico que será útil en la demostración del Teorema 1.

Lema 11. Sea T(x)=Φ(x)−Φ(x^∗)−H(x−x^∗),Existe una constante positiva 𝜌 tal que, si , ∥x−x∗∥_∞<Є entonces T(x) satisface

𝐴demás, si 𝐹′ es Lipschitz continua en una vecindad de x^∗ existe 𝛾>0 tal que

Demostración. Por la definición de 𝛷, se tiene que

luego, la fila 𝑖-ésima de este vector está dada por:

Por otra parte, la fila 𝑖-ésima de H(x-x*) es dada por:

Además, cuando 𝑥_𝑖^∗=𝐹𝑖(𝑥^∗), se puede escoger 𝒆_𝑖^𝑇 o ∇F_i(x)T (o una combinación lineal de estos). Ahora, por la definición de la función T, (33) y (34), se tiene que

luego

para algún 𝑘∈{1,…,𝑛}. Dado que 𝐹 es continuamente diferenciable, para cualquier 𝜌>0 existe 𝜀>0 tal que si ∥x−x∗∥_∞<𝜖 entonces∥𝑇(x)∥_∞≤𝜌 ∥x−x^∗∥_∞.

Por otra parte, como 𝐹′ es Lipschitz continua, para cualquier 𝜌₁ existe >0 tal que, si ∥x−x^∗∥_∞< entonces ∥𝑇(x)∥_∞≤𝜌1 ∥x−x^∗∥_∞².

Teorema 1. Existe una constante positiva 𝜀 tal que, si ∥x−x∗∥_∞<𝜀 la función 𝒬 definida por:

está bien definida, y dado 𝑟∈(0,1)

Además, si la matriz jacobiana de 𝐹 es Lipschitz continua en una vecindad de x^∗∈ℝ^𝑛,

Demostración. En primer lugar, se demuestra que (35) está bien definida, es decir, existe ??𝜇′(x)−1. Por la definición del C-subdiferencial, se tiene que la matriz 𝐻^∗∈∂_𝐶Φ(x) es de la forma

donde {y^k } es una sucesión de vectores de ℝ^𝑛, que converge a x^∗ y tal que Φ′(y^k ) existe. Se analiza a continuación ∥[Φ𝜇′(x)]−𝐻∥_∞.

Para algún 𝑗∈{1,…,𝑛}, se tiene que

Por el Lema 6,

Dado que 𝐹 es continua, para todo existe >0 tal que, si ∥x−x^∗∥_∞< entonces ∣𝐹𝑗(x )−𝐹𝑗(x^∗)∣<. Sea =mín{,}. Se consideran las dos posibilidades para el máximo en (38).

Por las consideraciones anteriores y por (38) para cualquier , se tiene que

Ahora

De donde ∥𝐻∗⁻¹𝛷𝜇′(x)−𝐼∥_∞<1.. Por el Lema de Banach [18] , existe 𝛷𝜇′(x)⁻¹ y por lo tanto, la función 𝒬 está bien definida. Además,

Para la segunda parte de la prueba, se demuestra (36) usando (35). Mediante algunas manipulaciones algebraicas, se tiene que

Por el Lema 9, para todo 𝛿>0, existe 𝜀₁>0 tal que si ‖x−x∗‖∥_∞<𝜀₁ entonces ∥𝛷𝜇′(x^∗)−𝐻∥<𝛿, para 𝐻∈𝜕_𝐶𝛷(x). En particular, para 𝛿= existe 𝜀_𝑟 tal que, si ∥x−x∗∥_∞<𝜀_𝑟,

Además, por el Lema 11, para todo 𝜌>0, existe 𝜀₂>0 tal que si ∥x−x∗∥_∞<𝜀₁ se tiene que

En particular, para 𝜌< existe >0 tal que si ∥x−x∗∥∞< se tiene que

Luego para 𝜀=mín{𝜀′,𝜀𝑟,}, por (40), (41) y (42) se tiene,

La tercera parte de la prueba consiste en demostrar (37) bajo la hipótesis de que la matriz jacobiana de 𝐹 es Lipschitz continua en x^∗∈ℝ^𝑛.Por el Lema 11 existe una constante positiva 𝛾 tal que,

Por otra parte, sumando y restando la matriz 𝐻^∗ en ∥Φ𝜇′(x^∗)−𝐻∥_∞ y aplicando la desigualdad triangular, se tiene que

Por (38), el primer término del lado derecho de (44) está acotado superiormente por √2𝑛𝜅𝑀, donde 𝑀=máx {∣𝑥_𝑗−𝑥_𝑗∗∣,∣𝐹𝑗(x)−𝐹𝑗(x^∗)∣} y por el Lema 10, el segundo término está acotado por 𝛾𝑛∥x−x∗∥∞. Observe que,

La última desigualdad se tiene porque 𝐹 es continuamente diferenciable. Sea =max{1,𝜁},

con 𝜃=𝑛(√2𝜅+𝛾), luego de (40), (43) y (45),

El siguiente teorema garantiza que el algoritmo propuesto está bien definido, converge linealmente y además da una condición suficiente para la convergencia cuadrática.

Teorema 2. Existe una constante ε₀ > 0 tal que si ‖ x₀- x* ‖_∞ ≤ ε₀, entonces la sucesión {x_k), generada por

está bien definida, converge a x∗ y dado 𝑟 ∈(0,1), satisface

Además, si la matriz jacobiana de 𝐹 es Lipschitz continua en una vecindad de x^∗ entonces

donde с es la constante del Teorema 1.

Demostración. Sea Q definida en (35). Para k = 0,1,...

Sean r Є ( 0,1) y ε₀ Є (0, ε), donde e es la constante del Teorema 1. La demostración es inducción sobre k.

•Para 𝑘=0. Si ∥x₀−x^∗∥_∞≤𝜀₀<𝜀, por el Lema 11, x1=𝒬(x0) está bien definida y satisface ∥x1−x^∗∥_∞≤𝑟∥x₀−x∗∥_∞. Además, si la matriz jacobiana de 𝐹 es Lipschitz continua en una vecindad de x^∗,∥x1−x^∗∥_∞≤𝑐∥x₀−x∗∥_∞².

•Hipótesis inductiva. Supongamos que ∥x_𝑘−x∗∥_∞<𝜀₀, para todo 0≤𝑘≤𝑚−1. Entonces por el Lema 11, x𝑚=𝒬(x𝑚−1) está bien definida y satisface

Además, ∥x𝑚−x^∗∥_∞≤𝑐∥x𝑚−1−x^∗∥_∞2, si la matriz jacobiana de 𝐹 es Lipschitz continua en una vecindad de x^∗.

A partir de (49), se tiene

así, el Lema 11 garantiza que x_m+1 ,está bien definido y satisface que

y si la matriz jacobiana de F es Lipschitz continua en una vecindad de x^*,entonces x_m+1, satisface que

Por lo tanto, se concluye que (47) y (48) se cumplen para todo k = 0,1,....

3.2. Experimentación numérica

En esta sección se analiza el desempeño numérico del algoritmo cuasi-Newton propuesto en la sección anterior que se llamará Algoritmo CNS. Para ello, se compara con dos algoritmos, el primero usa una suavización de la función de Fischer-Burmeister en su reformulación [15], el cual se llamará Algoritmo SF, y el segundo es un algoritmo tipo Newton generalizado que usa la función Mínimo, Algoritmo NM.

Se resolvieron diez problemas de complementariedad asociados con las funciones Billups [19] , Nash-Cornout (Nash-Co) [20], Kojima-Shindo (Koj-Sh), Kojima-Josephy (Koj-Jo), Mathiesen modificado (Math mod) [9], Hock-Schittkowski (HH-66) [21], Geiger-Kanzow (Gei-Kan) [22], Ahn [23], la función del Ejemplo A en [24] y finalmente, un problema que se llamó Ejemplo 1, asociado con la función F: ℝ² → ℝ², definida por F(x) = ((x₁- 2)² + x₂² - 4, (x₁- 3)² + 4x₂² - 4)^T el cual se construyó para hacer ilustraciones geométricas.

Los algoritmos y las funciones de prueba fueron codificados utilizando el software MATLAB R2022b, y las pruebas numéricas fueron realizadas en un computador con procesador Intel(R) Core(TM) i5-1035G1 CPU @ 1.00GHz 1.19 GHz. Se dirá que el algoritmo converge si ‖ Ф(х_к) ‖< 10^-6y diverge, si el número de iteraciones es mayor que 200, lo cual se indica con el símbolo (-). Para la variación de μ_k en los Algoritmos ACNS y ASF se usaron las sucesiones {2^-k} y {100^-k}.

3.2.1. Experimento 1

En primer experimento numérico se analiza el comportamiento local de los Algoritmos CNS, SF y NM, en términos de número de iteraciones usando los 10 problemas mencionados anteriormente y los siguientes puntos iniciales.

Los resultados obtenidos se presentan en la Tabla 1 que contiene la siguiente información: nombre del problema (Problema), dimensión del problema (Dim), punto inicial (x₀) , número de iteraciones con el Algoritmo CNS (K_CNS), número de iteraciones con el Algoritmo SF(K_SF), número de iteraciones con el Algoritmo NM(K_Nm). La sucesión {μ_k} usada se indica en la parte superior de los resultados respectivos.

En la Tabla 1, se observa que el Algoritmo CNS que se propone en este trabajo, tiene un buen comportamiento respecto al número de iteraciones, comparado con los otros dos. En general, dicho número es menor o igual que el del Algoritmo FS. No obstante, este último converge en casi todos los casos. En contraste el Algoritmo NM, cuando converge lo hace en menos iteraciones que los suavizados, lo cual no es sorprendente puesto que dicho algoritmo es tipo Newton, y los otros son cuasi-Newton. También se observa que el algoritmo propuesto converge en un número mayor de casos que el Algoritmo NM.

Un aspecto interesante que ilustra la Tabla 1 es que la elección: μ_k= {100-^k} influye en el número de iteraciones de los algoritmos suavizados, y el tiempo de ejecución dado en segundos, en ambos casos este número es menor o igual que los reportados con la sucesión μ_k = {2-^k}.

Para complementar este experimento, a partir de los datos reportados en la Tabla 1, se calcularon los índices de robustez (porcentaje de éxito del método), eficiencia (velocidad del método en sentido de número de iteraciones) y eficiencia combinada (combinación de los dos previos) [25] de los tres algoritmos. Para todos los índices, entre más cercano a uno esté su valor, el método será mejor.

Los resultados se ilustran en la Figura 3, en la cual se observa que el Algoritmo CNS tiene el mayor índice de eficiencia y eficiencia combinada mientras que, el Algoritmo SF es el más robusto.

3.2.2. Experimento 2

Un aspecto importante en el trabajo con algoritmos locales es determinar su radio de convergencia que corresponde al radio de una vecindad de la solución, en la cual cualquiera de sus puntos puede tormarse como inicial y el algoritmo converge. El propósito de este segundo experimento es determinar y comparar los radios de convergencia de los tres algoritmos usando los diez problemas mencionados anteriormente. Se considera la solución de cada problema x^* indicada en la Tabla 2 y se procede de la siguiente manera. Se generan 1000 puntos iniciales arbitrarios en una vecindad de esa solución, B_r(x^*), con un radio inicial r = 0.1, y se evalúa cada uno de los algoritmos; si para la totalidad de los puntos iniciales se logra convergencia, el radio se incrementa y se repite el proceso.

En caso contrario, se detiene y se registra el valor del radio antes del incremento. Para la métrica de las vecindades se usa la norma infinita y cada uno de los 1000 puntos aleatorios x se genera teniendo en cuenta que x = (x₁,x₂,x₃,...,x_n)^T pertenece a B_r(x^*) si x_i pertenece al intervalo [x^* -r,x^* - r], para todo i = 1,..,n. Estas componentes se calcularon mediante la parametrización, x_i = x* - r + 2t_i, con t_i Є [0,1]. Es de mencionar que en este experimento, se utilizó la sucesión {μ_k} = (100-^k).

Los resultados de este experimento se presentan en la Tabla 3, donde r_CNS, r_SF y r_NM indican el valor del radio en el que hubo convergencia para todos los puntos iniciales, el símbolo 0.01 - indica que el radio de convergencia es menor que 0.01 y 3 + indica que el radio de convergencia es mayor que 3.

La Tabla 3 permite observar que en general, los radios de convergencia de los tres algoritmos son pequeños, lo cual está de acuerdo con el carácter local de dichos algoritmos. Además, para la mitad de los problemas, el radio de convergencia del Algoritmo SF supera el de los otros dos, y se presenta un radio de convergencia grande en los dos problemas asociados a funciones lineales.

3.2.3. Experimento 3

Para ilustrar la región de convergencia de los tres algoritmos se considera el problema llamado Ejemplo1, descrito al inicio de esta Sección, el cual tiene dos soluciones: x*₁= (3.7748, 0.9219)^T y x₂* = (0, 0)^T . Se construye una malla rectangular de puntos iniciales igualmente espaciados a una distancia de 0.01 unidades, contenida en [1.77, 5.77] [- 1.07, 2.92] para un total de 160.801 puntos.

La Figuras 4, 5 y 6 ilustran la región de convergencia los Algoritmo CNS, FS y NM, respectivamente. En color verde limón aparecen los puntos iniciales para los cuales el algoritmo converge a x^*₁; en aguamarina los que convergen a x^*₂ y en azul los que divergen. Además, las curvas (x₁ - 2)² + x²₂ = 4 y (x₁ - 3)² + 4x²₂ = 4 aparecen dibujadas en color negro.

La Figura 4 ilustra la región de convergencia del algoritmo propuesto. En este caso, hay convergencia para el 99 % de los puntos iniciales generados y divergencia para el 1 % restante. Además, los resulta- dos muestran que, en caso de convergencia, solo con el 5.2 % de los puntos iniciales el algoritmo converge a la solución x^*₁.

En la Figura 5ab, correspondiente al Algoritmo FS, se observa una amplia región de convergencia formada por el 99.95 % de los puntos iniciales generados, a partir de los cuales con 67.8 % el algoritmo converge a x^*₁. Por otra parte, el algoritmo diverge para el 0.05 % de dichos puntos.

La Figura 6 muestra que el Algoritmo NM tiene una región de convergencia más pequeña que los otros dos algoritmos, y una región de divergencia mayor que los mismos. En efecto, converge para el 75 % de los puntos iniciales (con el 5.48% de estos puntos converge a x1), y diverge para el 25% restante. Este experimento permite observar dos aspectos interesantes: la región de convergencia de los algoritmos no necesariamente es un conjunto convexo y su tamaño puede aumentar usando una suavización, lo que hace de esta estrategia una buena alternativa para tener en cuenta en las propuestas algorítmicas que resuelven el problema de complementariedad.

3.2.4. Experimento 4

Finalmente, con el fin de analizar numéricamente la rapidez de convergencia del algoritmo propuesto, se calcula para cada problema usando la sucesión μ₀^2-k, el cociente R_k definido por

Los resultados para cuatro problemas (para los otros, ver [26]), se presentan en la Tabla 4, donde se reporta el número de iteraciones (k) y el valor de R_k. A pesar de que el algoritmo propuesto converge en un número menor de iteraciones que las reportadas en las siguientes tablas, se ha decidido mostrar más iteraciones para tener una mejor idea de la posible acotación del cociente R_k.

Los resultados de la Tabla 4 evidencian que el cociente Rk, para los cuatro problemas considerados, parece estar acotado, lo cual significa que el algo- ritmo propuesto puede ser cuadráticamente convergente. En efecto, dicho algoritmo converge cuadráticamente, bajo ciertas hipótesis, como se demostró en el Sección 4.

4. Conclusiones

Recientemente métodos con jacobiano suavizado han resultado muy eficientes para resolver el problema de complementariedad no lineal, vía su reformulación como un sistema de ecuaciones no lineales, el cual resulta ser no diferenciable.

En este trabajo se propone un nuevo algoritmo lo- cal cuasi-Newton suavizado para complementariedad no lineal usando la función Mínimo y una la suavización de la misma. Se demuestran algunas propiedades de estas dos funciones junto con la estrategia del jacobiano suavizado para obtener resultados de convergencia.

Pruebas numéricas preliminares muestran un buen desempeño local del algoritmo propuesto; además, se determina su región de convergencia, y se analiza numéricamente su rapidez de convergencia, la cual parece ser cuadrática, coincidiendo con lo demostrado teóricamente. En las pruebas numéricas se evidencia, que en la mayoría de casos, el método con jacobiano suavizado que utiliza una suavización de la función de complementariedad de Fischer es más exitoso que el propuesto, pero en cuanto a los promedios de iteraciones y tiempo, el algoritmo que proponemos lo supera.

Finalmente, sería interesante usar una estrategia híbrida, quizá usando la función de Fischer, para globalizar el algoritmo y analizar su desempeño numérico. Esto permitiría usar las ventajas ya conocidas mencionadas en la introducción de cada una de esas funciones.

Referencias

M. M. Kostreva, "Elasto-hydrodynamic lubrication: A non-linear complementarity problem," International Journal for Numerical Methods in Fluids, vol. 4, no. 4, pp. 377-397, 1984, doi: https://doi.org/10.1002/fld. 1650040407

M. Anitescu, J. F. Cremer, F. A. Potra, "On the existence of solutions to complementarity formulations of contact problems with friction," Complementarity and Variational Problems: State of the art. SIAM Publications, 1997, pp. 12-21.

M. C. Ferris and J. S. Pang, "Engineering and economic applications of complementarity problems," SIAM Review, vol. 39, pp. 669-713, 1997, doi: https://doi.org/10.1137/S0036144595285963

A. Chen, J.-S. Oh, D. Park, W. Recker, "Solving the bicriteria traffic equilibrium pro- blem with variable demand and nonlinear path costs," Applied Mathematics and Computation, vol. 217, no. 7, pp. 3020-3031, 2010, doi: https://doi.org/10.1016/j.amc.2010.08.035

R. Pérez, F. Arenas, H. J. Martínez, C. Arias, El problema de complementariedad no lineal: Teoría, aplicaciones y nuevos algoritmos para su solución. Programa Editorial de la Universidad del Valle, 2019.

C. Kanzow and H. Kleinmichel, "A new class of semismooth newton-type methods for nonlinear complementarity problems," Computational Optimization and Applications, vol. 11, no. 3, pp. 227251, 1998, doi: https://doi.org/10.1023/A:1026424918464

S. L. Hu, Z. H. Huang, and J. S. Chen, "Properties of a family of generalized ncp functions and a derivative free algorithm for complementarity problems," Journal of Computational and Applied Mathematics, vol. 230, no. 1, pp. 69-82, 2009, doi: https://doi.org/10.1016/j.cam.2008.10.056

L. Yong, "Nonlinear complementarity problem and solution methods," in Artificial Intelligence and Computational Intelligence, F. L. Wang, H. Deng, Y. Gao, and J. Lei, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2010, pp. 461 -469. doi: https://doi.org/10.1007/978-3-642-16530-6_55

J. S. Pang and L. Qi, "Nonsmooth equations: Motivation and algorithms," SIAM Journal on Optimization, vol. 3, no. 3, pp. 443-465, 1993, doi: https://doi.org/10.1137/0803021

A. Fischer and C. Kanzow, "On finite termina- tion of an iterative method for linear complementarity problems," Math. Program., vol. 74, no. 3, pp. 279-292, 1996, doi: https://doi.org/10.1007/BF02592200

A. Sherman, "On newton-iterative methods for the solution of systems of nonlinear equations," SIAM Journal on Numerical Analysis, vol. 15, no. 4, pp. 755711, 1978, doi: https://doi.org/10.1137/0715050

C. G. Broyden, J. E. Dennis, and J. J. Moré, "On the Local and Superlinear Convergence of Quasi-Newton Methods," IMA J. Appl. Math., vol. 12, no. 3, pp. 223245, 1973, doi: https://doi.org/10.1093/imamat/12.3.223

V. L. R. Lopes, J. M. Martínez, and R. Pérez, "On the local convergence of quasi-newton methods for nonlinear complementary problems," Applied Numerical Mathematics, vol. 30, pp. 3-22, 1999, doi: https://doi.org/10.1016/S0168-9274(98)00080-4

C. Kanzow and H. Pieper, "Jacobian smoothing methods for nonlinear complementarity problems," SIAM Journal on Optimization, vol. 9, no. 2, pp. 342373, 1999, doi: https://doi.org/10.1137/S1052623497328781

F. E. Arenas, H. J. Martínez, and R. Pérez, "A local jacobian smoothing method for solving nonlinear complementarity problems," Universitas Scientiarum, vol. 25, no. 1, p. 149-174, May 2020, doi: https://doi.org/10.11144/JAVERIANA.SC25-1.ALJS

F. H. Clarke, Optimization and Nonsmooth Analysis. Society for Industrial and Applied Mathematics, 1990.

X. Chen, "Smoothing methods for complementarity problems and their applications: a survey," J. Oper. Res. Soc. Japan, vol. 43, no. 1, pp. 32-47, 2000.

J. E. Dennis and R. B. Schnabel, Numerical methods for unconstrained optimization and non-linear equations. Society for Industrial and Applied Mathematics, 1996.

S. C. Billups, "Algorithms for complementarity problems and Generalized Equations," Ph.D. dissertation, University of Wisconsin, 1995.

P. T. Harker, "Accelerating the convergence of the diagonalization and projection algorithms for finite-dimensional variational inequalities," Mathematical Programming, vol. 41, pp. 29-59, 1988, doi: https://doi.org/10.1007/BF01580752

J. Zhu and B. Hao, "A new smoothing method for solving nonlinear complementarity problems," Open Mathematics, vol. 17, no. 1, pp. 104-119, 2019, doi: https://doi.org/10.1515/math-2019-0011

C. Geiger and C. Kanzow, "On the resolution of monotone complementarity problems," Computational Optimization and Applications, vol. 5, no. 2, pp. 155173, 1996, doi: https://doi.org/10.1515/math-2019-0011

A. Byong-Hun, "Iterative methods for linear complementarity problems with upperbounds on primary variables," Mathematical Programming, vol. 26, no. 3, pp. 295-315, oct 1983.

Y. Xia, H. Leung, and J. Wang, "A projection neural network and its application to constrained optimization problems," IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, vol. 49, no. 4, pp. 447-458, 2002, doi: https://doi.org/10.1109/81.995659

S. Buhmiler and N. Krejic, "A new smoothing quasi-newton method for nonlinear complementarity problems," Journal of Computational and Applied Mathematics, vol. 211, no. 2, pp. 141-155, feb 2008, doi: https://doi.org/10.1016/j.cam.2006.11.007

C. Quirá, "Un algoritmo local cuasi-Newton suavizado para resolver el problema de complementariedad no lineal," Ph.D. dissertation, Universidad del Cauca, Popayán, Colombia, Agosto 2021.

Notes

Author notes

^arosana@unicauca.edu.co^bfarenas@unicauca.edu.co^cdicorrea@unicauca.edu.co, quira@uaiinpebi-cric.edu.co

Conflict of interest declaration

IR