Optimización con restricciones y control óptimo de sistemas discretos

Traballo Fin de Grao Optimización con restricciones y control óptimo de sistemas discretos Sandra Díaz Seoane 2019/2020 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao Optimización con restricciones y control óptimo de sistemas discretos Sandra Díaz Seoane Septiembre,2020 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA Trabajo propuesto Área de Coñecemento: Matemática aplicada Título: Optimización con restricciones y control óptimo de sistemas discretos Breve descrición do contido Nos interesamos por problemas de optimización del tipo Minimizar ˜ J(u)u∈U donde u , la variable de minimización o variable de control, pertenece a un conjunto de controles llamados admisibles. Con frecuencia la expresión del funcional ˜ J(u) no admite una expresión analítica sencilla pero sin embargo verica ˜ J(u) = J(y(u), u) , donde y(u) , llamada la variable de estado, es la única solución de una ecuación que representamos por F(y, u) = 0 , y que recibe el nombre de sistema de estado. Los problemas de optimización con estas características reciben el nombre de problemas de control óptimo y debido a su estructura particular, admiten técnicas de análisis y resolución especiales. Nos centraremos en el caso en el que tanto la variable de estado como la de control son vectores y el sistema de estado vendrá dado por un sistema de ecuaciones algebraicas. En el trabajo se particularizarán los conceptos propios de problemas de optimización (tales como condiciones de optimalidad o el cálculo del gradiente del funcional coste) a problemas con esta estructura. Por otra parte, se propondrán e implementarán métodos numéricos sencillos para su resolución numérica. iii Índice general Resumen viii Introducción xi 1. Optimización sin restricciones 1 1.1. Preliminares.................................... 1 1.2. Existenciademínimos.............................. 6 1.3. Métodosderesolución .............................. 8 1.3.1. Seleccióndepaso............................. 11 1.3.2. Selección de dirección . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2. Optimización con restricciones 23 2.1. Restricciones de igualdad . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.1.1. Condiciones basadas en los multiplicadores de Lagrange . . . . . . . 24 2.1.2. Método de penalización, caso de restricciones de igualdad . . . . . . 25 2.2. Restricciones de igualdad y desigualdad . . . . . . . . . . . . . . . . . . . . 26 2.2.1. Condiciones de optimalidad de primer orden . . . . . . . . . . . . . . 27 2.2.2. Condiciones de optimalidad de segundo orden . . . . . . . . . . . . . 31 2.2.3. Método de penalización . . . . . . . . . . . . . . . . . . . . . . . . . 35 3. Control óptimo de sistemas discretos 39 3.1. Planteamiento del problema . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.1.1. Planteamiento como problema de optimización con restriciones . . . 40 3.1.2. Planteamiento como problema de optimización sin restricciones . . . 42 3.2. Métododeresolución............................... 43 3.3. Aplicaciones.................................... 49 3.3.1. Sistemas dinámicos discretos . . . . . . . . . . . . . . . . . . . . . . 49 3.3.2. Aprendizaje de una red neuronal articial . . . . . . . . . . . . . . . 51 v vi ÍNDICE GENERAL 3.3.3. Discretización de un sistema con EDO . . . . . . . . . . . . . . . . . 60 Bibliografía 67 2 CAPÍTULO 1. OPTIMIZACIÓN SIN RESTRICCIONES Demostración. 1. Véase la proposición 6.2.5, capítulo 6, página 216 del libro [1]. 2. Véase el resultado 173, capítulo 9, página 347 del libro [2]. Vemos ahora resultados relacionados con matrices simétricas ya que habrá menciones a ellas a lo largo de este texto. Vemos posteriormente la denición, así como relaciones con autovalores y autovectores, de matriz denida positiva y demás. Estos conceptos se usarán en resultados en los cuales se habla de condiciones necesarias y sucientes para la existencia de mínimos (Teoremas 1.23 y 1.24). Proposición 1.3. (Propiedades de las matrices simétricas) Sea A∈Rn×n simétrica. Entonces 1. Todos los autovalores de A son reales. 2. A es ortogonalmente similar a una matriz diagonal, i.e., existe una matriz ortogonal Q∈Rn×n tal que Q−1AQ =QTAQ =    λ1 ... λn     , con λ1, . . . , λn autovalores de A . Demostración. Véase el capítulo 8, teorema 8.4.7, página 380 del libro [1]. Véase el Capítulo 8, teorema 8.4.8, página 380 del libro [1]. Denición 1.4. Sea A∈Rn×n simétrica. A se dice denida positiva si vTAv > 0 , ∀v∈Rn , v6= 0 . Se dice que es semidenida positiva si vTAv ≥0,∀v∈Rn . A es denida negativa o semidenida negativa si −A es denida positiva o semidenida positiva respectivamente. A se dice indenida si no es ni semidenida positiva ni semidenida negativa. Dada la denición 1.4 tenemos a partir de la proposición 1.3 el siguiente corolario. Corolario 1.5. (Matrices denidas y semidenidas positivas y negativas e indenidas) Sea A∈Rn×n simétrica: 1. A es denida positiva ⇔ todos sus autovalores son positivos. 1.1. PRELIMINARES 3 2. A es semidenida positiva ⇔ todos sus autovalores son no negativos. 3. A es denida negativa o semidenida negativa si todos sus autovalores son negativos o no positivos respectivamente. 4. A es entonces indenida si tiene autovalores tanto positivos como negativos. Veremos ahora las deniciones y resultados relacionados con el análisis matemático. Comenzamos con deniciones de función continua, continuamente diferencianble y dos veces continuamente diferenciable. Denición 1.6. Una función f:Rn→R se dice continua en ˆx∈Rn si, para cada  > 0 , existe un δ > 0 tal que kx−ˆxk< δ implica kf(x)−f(ˆx)k<  . Si f es continua en cada punto de un conjunto abierto B⊂Rn , entonces de dice que f es continua en B , f∈C(B) . Denición 1.7. Una función continua f se dice continuamente diferenciable en x∈Rn , si ∂f ∂xi(x) existe y es continua para i= 1, . . . , n . El gradiente de f en x se dene como ∇f(x) = ∂f ∂x1 (x),..., ∂f ∂xn (x)T . Si f es continuamente diferenciable en todo punto de un conjunto B⊂Rn , entonces se dice que f es continuamente diferenciable en B , f∈C1(B) . Denición 1.8. Una función continuamente diferenciable f se dice dos veces continuamente diferenciable en x∈Rn , si ∂2f ∂xi∂xj(x) existen y son continuas para 1≤i, j ≤n . La Hessiana de f en x se dene como la matriz n×n cuyos elementos son ∇2f(x)ij =∂2f ∂xi∂xj(x),1≤i, j ≤n. Denotaremos también a la Hessiana de f en x como Hf(x) . Si f es dos veces continuamente diferenciable en todo punto de un conjunto B⊂Rn , entonces se dice que f es dos veces continuamente diferenciable en B , f∈C2(B) . Observación 1.9 . El teorema de Schwarz que podemos encontrar en el capítulo 5, teorema 5.6, página 194 de [3] arma que la matriz hessiana de una función dos veces continuamente diferenciable es simétrica. Finalizando esta sección tenemos deniciones de convexidad tanto para conjuntos como para funciones, así como resultados relacionados con estás. En primer lugar la denición y resultados de conjuntos convexos. Denición 1.10. Sea S⊂Rn . S es convexo si ∀x, y ∈S, λx+(1−λ)y∈S, ∀λ∈[0,1] . 4 CAPÍTULO 1. OPTIMIZACIÓN SIN RESTRICCIONES Observación 1.11 . Geométricamente la denición anterior nos indica que para cada dos puntos del conjunto el segmento que los une está contenido en dicho conjunto. Por lo cual vemos que todo conjunto convexo es conexo. Proposición 1.12. Sean S, S1, S2∈Rn conjuntos convexos. Entonces: 1. S1TS2 es convexo. 2. S1±S2={x1±x2|x1∈S1, x2∈S2} es convexo. 3. El interior de S es convexo. 4. La clausura de S , ¯ S , es un conjunto convexo. Demostración. 1) Sean x, y ∈S1TS2 . Tenemos que x, y ∈S1 y por ser S1 convexo tenemos que λx + (1 −λ)y∈S1,∀λ∈[0,1] . Análogamente para S2 tenemos λx + (1 −λ)y∈S2,∀λ∈[0,1] . Entonces λx+(1−λ)y∈S1TS2,∀λ∈[0,1] y hemos demostrado que S1TS2 es convexo. 2) Sean x, y ∈S1±S2 . Tenemos x=x1±x2 para ciertos x1∈S1 , x2∈S2 e y=y1±y2 para ciertos y1∈S1 , y2∈S2 . Por ser S1 convexo tenemos λx1+ (1 −λ)y1∈S1 , ∀λ∈ [0,1] ; análogamente para S2 ; λx2+ (1 −λ)y2∈S2 , ∀λ∈[0,1] . Entonces tenemos λx1+(1−λ)y1±λx2+(1−λ)y2=λ(x1±x2)+(1−λ)(y1±y2) = λx+(1−λ)y∈S1±S2. Con lo cual queda demostrado que S1±S2 es convexo. 3) Sean x e y dos puntos en el interior de S y z=λx+(1−λ)y con λ∈[0,1] . Escogemos un δ > 0 tal que B(y, δ)⊂S , siendo B(y, δ) el entorno de y con radio δ . Tenemos ||z−x|| ||y−x|| = 1 −λ . Entonces B(z, (1 −λ)δ) es el conjunto λx + (1 −λ)B(y, δ) que está contenido en S . Finalmente B(z, (1 −λ)δ)⊂S prueba que z pertenece al interior de S y que este es convexo. 4) Sean x, y ∈¯ S escogemos en S dos sucesiones xk e yk convergente a x e y respectivamente. Entonces sea λ∈[0,1] tenemos ||[λxk+ (1 −λ)yk]−[λx + (1 −λ)y]|| =||λ(xk−x) + (1 −λ)(yk−y)|| ≤λ||xk−x||+ (1 −λ)||yk−y||. Haciendo el límite tenemos l´ım k→∞ ||[λxk+ (1 −λ)yk]−[λx + (1 −λ)y]|| = 0, lo cual demuestra que λx + (1 −λ)y∈¯ S , es decir, ¯ S es convexo. 1.1. PRELIMINARES 5 Seguimos ahora con funciones convexas. Denición 1.13. Sea S⊂Rn un conjunto no vacío y convexo. La función f:S→R es convexa si: ∀x, y ∈S, f (λx + (1 −λ)y)≤λf(x) + (1 −λ)f(y),∀λ∈[0,1], y es estrictamente convexa si: ∀x, y ∈S tales que x6=y, f (λx + (1 −λ)y)< λf(x) + (1 −λ)f(y),∀λ∈(0,1). Proposición 1.14. Dadas J1, . . . , Jn funciones convexas y los números reales α1≥0, . . . , αn≥0 , entonces n P i=1 αiJi es una función convexa. Demostración. Sean x, y ∈S y λ∈[0,1] , entonces n X i=1 αiJi(λx + (1 −λ)y)≤ n X i=1 αiλJi(x) + αi(1 −λ)Ji(y) =λ n X i=1 αiJi(x) + (1 −λ) n X i=1 αiJi(y). Se tienen los resultados siguientes cuyas demostraciones las encontraremos en el libro [4]. Proposición 1.15. Sean S⊂Rn abierto y convexo y f:S⊂Rn→R convexa. Entonces f es continua en S . Demostración. Véase el capítulo 1, teorema 1.3.12, página 40 de [4]. Proposición 1.16. Sea S abierto y convexo, f:S→R diferenciable. Entonces f es convexa si, y solo si, f(y)≥f(x)+∇f(x)t(y−x),∀x, y ∈S . Además f es estrictamente convexa si, y solo si, f(y)> f(x) + ∇f(x)t(y−x),∀x, y ∈S con y6=x . Demostración. Véase el capítulo 1, teorema 1.3.13, página 42 de [4]. Proposición 1.17. Sea S abierto y convexo, f:S→R , f∈C2(S) . Entonces f es convexa si, y solo si, Hf(x) es semidenida positiva ∀x∈S . Además f es estrictamente convexa si Hf(x) es denida positiva ∀x∈S . Demostración. Véase el capítulo 1, teorema 1.3.14, página 44 de [4]. 6 CAPÍTULO 1. OPTIMIZACIÓN SIN RESTRICCIONES Observación 1.18 . Diremos que f:Rn→R es cuadrática si f(x) = 1 2xTAx +xTb+c con A∈Rn×n simétrica, b∈Rn y c un número real. Para este tipo de funciones, tenemos: 1. ∇f(x) = Ax +b, Hf(x) = A∀x∈Rn . 2. f en convexa ⇔A es semidenida positiva. 3. f es estrictamente convexa ⇔A es denida positiva. 1.2. Condiciones necesarias y sucientes para la existencia de mínimos La forma general de un problema de optimización sin restricciones es m´ın x∈XJ(x) donde x∈Rn es la variable de optimización, J(x) la función objetivo y X≡Rn su conjunto o región factible. Vamos a realizar un estudio de las condiciones necesarias y sucientes para la existencia de mínimos locales y globales. Comenzamos deniendo estos mínimos. Denición 1.19. Un punto ˆx se llama mínimo local si ∃δ > 0 tal que J(ˆx)≤J(x)∀x∈Rn que satisface ||x−ˆx|| < δ . Un punto ˆx se llama mínimo local estricto si ∃δ > 0 tal que J(ˆx)< J(x)∀x∈Rn con x6= ˆx y ||x−ˆx|| < δ . Denición 1.20. Sea ˆx∈X , un mínimo local tal que existe un entorno de ˆx , B(ˆx, δ) con δ > 0 , en el cual ˆx es el único mínimo local en XTB(ˆx, δ) ; entonces diremos que ˆx es un mínimo local aislado . Denición 1.21. Un punto ˆx se llama mínimo global si J(ˆx)≤J(x)∀x∈Rn . Un punto ˆx se llama mínimo global estricto si J(ˆx)< J(x)∀x∈Rn con x6= ˆx . Teorema 1.22 (Condición necesaria de primer orden) . Sea J:D⊂Rn→R continuamente diferenciable en el abierto D , J∈C1(D) . Si ˆx∈D es mínimo local de J , entonces ∇J(ˆx)=0 . Demostración. Sea ˆx∈D un mínimo local y consideramos la sucesión xk= ˆx−αk5J(ˆx), αk>0 , con {αk} tendiendo a cero. Usando el desarrollo de Taylor, para k sucientemente grande, 0≤J(xk)−J(ˆx) = −αk∇J(ηk)T∇J(ˆx), 1.2. EXISTENCIA DE MÍNIMOS 7 donde ηk es una combinación convexa de xk y ˆx . Ahora bien pasando al límite cuando xk tiende a ˆx después de haber dividido por αk , se sigue de J∈C1(D) que 0≤ −||∇J(ˆx)||2 lo cual quiere decir que ∇J(ˆx) = 0 . Seguimos ahora con condiciones necesarias de primer y segundo orden para funciones generales. Teorema 1.23 (Condición necesaria de segundo orden) . Sea J:D⊂Rn→R , D un abierto y J∈C2(D) . Si ˆx∈D es mínimo local de J , entonces ∇J(ˆx)=0 y Hf(ˆx) es semidenida positiva. Demostración. Por el teorema anterior tenemos que ∇J(ˆx)=0 por tanto, solo tenemos que demostrar que HJ(ˆx) es semidenida positiva. Consideramos la sucesión xk= ˆx−αkd, αk>0 , con {αk} tendiendo a cero, donde d es arbitrario. Como J∈C2(D) y ∇J(ˆx)=0 , entonces el desarrollo de Taylor que tenemos para k sucientemente grande será 0≤J(xk)−J(ˆx) = 1 2α2 kdTHJ(ηk)d, donde ηk es una combinación convexa de xk y ˆx . Dividiendo ahora entre 1 2α2 k y pasando al límite, tenemos dTHJ(ˆx)d≥0,∀d∈Rn . Tenemos a continuación una condición suciente de segundo orden para funciones generales. Teorema 1.24 (Condición suciente de segundo orden) . Sea J:D⊂Rn→R , D un abierto y J∈C2(D) . Si ∇J(ˆx)=0 y HJ(ˆx) es denida positiva, entonces ˆx∈D es mínimo local estricto de J . Demostración. Supongamos que ∇J(ˆx)=0 y HJ(ˆx) es denida positiva. Usando el desarrollo de Taylor, para cualquier vector d∈Rn tal que ˆx+d esté en un entorno de ˆx en el cual HJ(ˆx+d) es denida positiva, tenemos J(ˆx+d) = J(ˆx) + 1 2dTHJ(ˆx+θd)d , donde θ∈(0,1) . Entonces podemos escoger un δ > 0 tal que ˆx+d∈B(ˆx, δ) y dTHJ(ˆx+θd)d > 0 . Entonces, J(ˆx+d)> J(ˆx) y tenemos demostrado el teorema. Terminamos esta sección con condiciones necesarias y sucientes pero en este caso en relación a funciones convexas. La primera es una condición suciente para la existencia de un mínimo global y la segunda es necesaria y suciente. Proposición 1.25. Sea S⊂Rn un conjunto convexo y J:S→R una función convexa. Supongamos que existe ˆx∈S mínimo local de J . Entonces: 1. ˆx es mínimo global. 8 CAPÍTULO 1. OPTIMIZACIÓN SIN RESTRICCIONES 2. Si J es además estrictamente convexa, ˆx es el único mínimo global. Demostración. 1. Supongamos que existe x∈S tal que J(x)< J(ˆx) . Para todo λ∈(0,1) se tiene: J(λx + (1 −λ)ˆx)≤λJ(x)−(1 −λ)J(ˆx)< λJ(ˆx) + (1 −λ)J(ˆx) = J(ˆx). Luego existen puntos tan próximos a ˆx como queramos en los que J toma un valor menor que J(ˆx) , lo que contradice el hecho de que ˆx sea mínimo local. 2. Supongamos existe x∈S\{ˆx}:J(x) = J(ˆx) . Para todo λ∈(0,1) se tiene: J(λx + (1 −λ)ˆx)< λJ(ˆx) + (1 −λ)J(ˆx) = J(ˆx). Llegamos de nuevo a una contradicción. Teorema 1.26. Sea S⊂Rn un conjunto abierto y convexo y J:S→R una función diferenciable y convexa. Entonces: ˆx es mínimo global ⇔ ∇J(ˆx)=0. Demostración. " ⇒ " ˆx mínimo global ⇒ˆx mínimo local ⇒ ∇J(ˆx)=0 . " ⇐ " Por la Proposición 1.16 sabemos que J(x)≥J(ˆx) + ∇J(ˆx)T(x−ˆx) = J(ˆx)∀x∈S , pues ∇J(ˆx)=0 . 1.3. Métodos de resolución Recordemos que nuestro problema consiste en hallar un ˆx∈S tal que J(ˆx) = m´ın x∈SJ(x) . Nuestros métodos de resolución buscarán disminuir el funcional coste. Resulta interesante destacar que los métodos que presentamos supondrán una regularidad mínima sobre J , J∈C1 , o incluso J∈C2 . Bajo estás condiciones si existe ese ˆx tenemos que ∇J(ˆx) = 0 . Así ∇J(ˆx)=0 se usará como test de parada. Las siguientes deniciones son interesantes para la implementación de los métodos. Denición 1.27. Sea S abierto, J:S→R diferenciable, x∈S , u∈Rn y α el ángulo que forman ∇J(x) y u . ∂J ∂u(x) := h∇J(x), ui=||∇J(x)|| ||u||cos(α) es la derivada direccional de J en el punto x en la dirección u . 1.3. MÉTODOS DE RESOLUCIÓN 9 Observación 1.28 . Nótese que: ∂J ∂u (x) es máximo cuando α= 0 , es decir, ∇J(x) indica la dirección de máximo crecimiento. ∂J ∂u (x) es mínimo cuando α=π , es decir, −∇J(x) indica la dirección de mínimo crecimiento o máximo decrecimiento. Denición 1.29. Sea S abierto, J:S→R diferenciable y x∈S . Dado d∈Rn diremos que J decrece desde x en la dirección d si h∇J(x), di<0 . En tal caso, d se llama dirección de descenso para J desde el punto x . Vamos a ver ahora el uso de métodos iterativos para encontrar el mínimo en un problema de optimización. Nos centraremos en los llamados métodos de descenso. La idea es, dado un punto inicial, crear una sucesión haciendo uso de un método iterativo. Queremos con esto llegar a la solución del problema, construyendo así una sucesión que converge al mínimo. Para nalizar con las iteraciones tendrá que cumplirse un cierto criterio de convergencia. También se dará un número máximo de iteraciones el cual si es alcanzado se dirá que el método usado no converge. Uno de los criterios que vamos a usar y quizás el más intuitivo será ||∇J(xk)|| ≤ δ , donde δ será lo que denominaremos parámetro de tolerancia . Que este criterio se cumpla nos indica que el gradiente tiende a cero y que la sucesión iterativa converge a un punto estacionario. Existen distintos métodos iterativos de optimización. Muchos de ellos se engloban dentro de los métodos de descenso denominados así debido a que el valor de la función objetivo en cada iterante es menor que su valor en el anterior. Veamos el algoritmo básico de optimización esquematizado Algoritmo 1.30. En primer lugar se jan un número máximo de iteraciones nitmax y los parámetros de tolerancia  > 0 y δ > 0 (pequeños) que se usarán como criterios de parada. Paso 0 (Inicialización) Se proporciona una semilla x0∈Rn , si es posible cerca de ˆx , k:= 0 . Paso 1 (Dirección de descenso) Se elige una dirección de descenso dk para J desde xk , se debe cumplir h∇J(xk), dki<0 . Paso 2 (Cálculo de paso) Se calcula αk∈(0,+∞) paso tal que J(xk+αkdk)< J(xk) . Paso 3 (Bucle) Se actualiza xk+1 =xk+αkdk , k:= k+ 1 . 10 CAPÍTULO 1. OPTIMIZACIÓN SIN RESTRICCIONES Paso 4 (Criterios de convergencia) Se comprueban que ||xk+1−xk|| 1+||xk+1|| ≤ y ||∇J(xk+1)|| ≤ δ si se cumple se terminan las iteraciones y xk+1 se considera una buena aproximación de ˆx . En otro caso se vuelve al Paso 1. Si el número de iteraciones alcanza nitmax sin que se cumplan los criterios de convergencia concluimos que el método no converge. Veremos ahora dos resultados que nos resultará de utilidad para demostraciones posteriores (Teorema 1.38 y Teorema 1.47 respectivamente) en relación con el algoritmo anterior. Teorema 1.31. Sea αk>0 la solución de J(xk+αkdk) = m´ın α∈[0,+∞)J(xk+αdk) , donde dk es una dirección de descenso. Sea M > 0 tal que ||HJ(xk+αdk)|| ≤ M∀α > 0 . Entonces J(xk)−J(xk+αkdk)≥1 2M||∇J(xk)||2cos2hdk,−∇J(xk)i. Demostración. Por hipótesis tenemos que J(xk+αdk)≤J(xk) + αdT k∇J(xk) + α2 2M||dk||2,∀α > 0. Dado ˆα=−dT k∇J(xk)/(M||dk||2) ; tenemos que J(xk)−J(xk+αkdk)≥J(xk)−J(xk+ ˆαdk) ≥ −ˆαdT k∇J(xk)−ˆα2 2M||dk||2 =1 2 (dT k∇J(xk))2 M||dk||2 =1 2M||∇J(xk)||2(dT k∇J(xk))2 ||dk||2||∇J(xk)||2 =1 2M||∇J(xk)||2cos2hdk,−∇J(xk)i. Teorema 1.32. Sea ∇J(x) uniformemente continua en el conjunto de nivel L(x0) = {x∈ Rn|J(x)< J(x0)} . Sea θk el ángulo que forman ∇J(xk) y la dirección dk elegida en el algoritmo 1.30 cumpliendo θk≤π 2−µ , para algún µ > 0 . Entonces ∇J(xk)=0 para un cierto k ; o bien J(xk)→ −∞ ; o ∇J(xk)→0 . Demostración. Supongamos que para todo k , ∇J(xk)6= 0 y J(xk) está acotada inferiormente. Como {J(xk)} es monótona decreciente, existe su límite. Además J(xk)−J(xk+1)→0. (1.1) 1.3. MÉTODOS DE RESOLUCIÓN 11 Supongamos que no se cumple ∇J(xk)→0 . Entonces existe  > 0 y un subconjunto K , tal que ||∇J(xk)|| ≥ ∀k∈K . Además −∇J(xk)Tdk/||∇J(xk)||cos θk≥sin µ=1 (1.2) Notese además que J(xk+αdk) = J(xk) + α∇J(ξk)Tdk =J(xk) + α∇J(xk)Tdk+α[∇J(ξk)−∇J(xk)]Tdk ≤J(xk) + α||dk||∇J(xk)Tdk ||dk|| +||∇J(ξk)−∇J(xk)||, (1.3) donde ξk se encuentra entre xk y xk+αdK . Como ∇J(x) es uniformemente continuo en el conjunto de nivel L(x0) , existe ¯α tal que cuando 0≤α||dk|| ≤ ¯α , tenemos ||∇J(ξk)−∇J(xk)|| ≤ 1 21 (1.4) Por (1.2)-(1.4) tenemos Jxk+ ¯αdk ||dk||≤J(xk) + ¯α∇J(xk)Tdk ||dk|| +1 21≤J(xk)−1 2¯α1. Entonces J(xk+1)≤Jxk+ ¯αdk ||dk|| ≤J(xk)−1 2¯α1 , lo cual contradice (1.1). Por lo tanto ∇J(xk)→0 . 1.3.1. Selección de paso En esta sección vamos a ver distintas reglas que nos permitirán calcular el paso a usar en los algoritmos de optimización. Nos va a resultar útil considerar la función j: (0,+∞)→R denida por j(α) = J(xk+αdk). Para que α sea adecuado tendrá que garantizar descenso, es decir, debería cumplirse j(α)< j(0) . Observación 1.33 . Si J es diferenciable también lo será j(.) y j0(α) = dt k∇J(xk+αdk) , entonces j0(0) = dt k∇J(xk)<0 . Si J es dos veces diferenciable también lo será j(.) y j00(α) = dt kHJ(xk+αdk)dk . Omitiremos los subíndices k para simplicar la notación. 18 CAPÍTULO 1. OPTIMIZACIÓN SIN RESTRICCIONES Gradiente conjugado Veremos ahora el método del gradiente conjugado que presenta una convergencia más rápida que el gradiente y evita el calculo de hessianas que puede ser costoso ya que este método solo usa las derivadas de primer orden. Consideraremos para la presentación de estos métodos el caso cuadrático asumiendo que J(x) = 1 2xTAx +bTx+c . Introducimos en primer lugar las direcciones conjugadas y los métodos que usan estas direcciones pues el método del gradiente conjugado es uno de ellos. Denición 1.48. Sean A una matriz simétrica y denida positiva n×n , d1, . . . , dm∈Rn vectores distintos de cero, m≤n . Si dT iAdj= 0,∀i6=j , los vectores d1, . . . , dm se dicen A− conjugadas o simplemente conjugadas . Observación 1.49 . Si A=mI ser A− conjugado equivale a ser ortogonal. Algoritmo 1.50. (Método de dirección conjugada general) Paso 0 Damos x0 ,  > 0 , k:= 0 . Paso 1 Calculamos g0=g(x0) , d0 tal que dT 0g0<0 . Paso 2 Si ||gk|| ≤  , ˆx=xk y terminamos. Paso 3 Calculamos αk tal que J(xk+αkdk) = m´ın α∈[0,+∞)J(xk+αdk) . Paso 4 Construimos xk+1 =xk+αkdk . Paso 5 Calculamos dk+1 con algún método de dirección conjugada, tal que dT k+1Adj= 0 , j= 0,1, . . . , k . Paso 6 Actualizamos k:= k+ 1 y volvemos las Paso 2. Teorema 1.51. Dada una función cuadrática con matriz Hessiana denida positiva el método de dirección conjugada converge en a lo sumo n iteraciones de búsqueda lineal exacta. Demostración. Demostración en el capítulo 4, teorema 4.1.3, página 176 de [4]. Veremos ahora el método del gradiente conjugado con una explicación de como calcular las direcciones conjugadas. En primer lugar nos centraremos en el caso de funcionales cuadráticos. Tenemos J(x) = 1 2xTAx +bTx+c , donde A es una matriz simétrica denida positiva n×n , b∈Rn y c un número real. Denimos g(x) = Ax+b el gradiente de J en x . Entonces 1.3. MÉTODOS DE RESOLUCIÓN 19 con d0=−g0 tenemos x1=x0+α0d0 , donde α0 es óptimo. Así tenemos gT 1d0= 0 . Ahora d1=−g1+β0d0 y escogemos β0 tal que dT 1Ad0= 0 . Podemos calcular β0 multiplicando la formula de d1 por dT 0A y obtenemos β0=gT 1Ad0 dT 0Ad0=gT 1(g1−g0) dT 0(g1−g0)=gT 1g1 gT 0g0 nos interesa esta equivalencia ya que tenemos una formulación solo involucrando gradientes. En general tendremos βk=gk+1Tgk+1 gT kgk (fórmula de Fletcher-Reeves) y siendo αk paso óptimo αk= −gT kdk dT kAdK . Algoritmo 1.52. (Gradiente conjugado caso cuadrático) Paso 0 Damos x0∈Rn . Paso 1 Inicializamos g0=Ax0+b , d0=−g0 . ( d0=−5J(x0) ) Paso 2 Si d0= 0 , ˆx=x0 y terminamos. Paso 3 Calculamos el paso óptimo αk=−gT kdk dT kAdk . Paso 4 Actualizamos xk+1 =xk+αkdk y gk+1 =Axk+1 +b . Paso 5 Si gk+1 = 0 , ˆx=xk+1 y terminamos. Paso 6 Calculamos βk=gT k+1gk+1 gT kgk . Paso 7 Actualizamos dk+1 =−gk+1 +βkdk y volvemos al Paso 3. Algoritmo 1.53. (Gradiente conjugado caso general) Paso 0 Damos x0∈Rn y  < 0 . Paso 1 g0=∇J(x0) , d0=−g0 . Paso 2 Si ||g0|| ≤ ˆx=x0 y terminamos. Paso 3 Calculamos el paso óptimo o aproximado. Paso 4 Actualizamos xk+1 =xk+αkdk . Paso 5 gk+1 =∇J(xk+1) si ||gk+1|| ≤ ˆx=xk+1 y terminamos. Paso 6 Calculamos βk=||∇J(xk+1)||2 ||∇J(xk)||2 . Paso 7 Actualizamos dk+1 =−gk+1 +βkdk y volvemos al Paso 3. 20 CAPÍTULO 1. OPTIMIZACIÓN SIN RESTRICCIONES Métodos cuasi-Newton El coste computacional de calcular las Hessianas puede ser muy alto al igual que puede ser difícil evaluarlas e incluso que no sean calculables analíticamente. Para solucionar estos problemas que se pueden presentar en el método de Newton usamos los cuasi-Newton. Esta familia de métodos solo usa la evaluación de las funciones y de los gradientes y en relación con la Hessianas lo que se hace es generar una sucesión de aproximaciones. Cabe resaltar que a pesar de realizar estas aproximaciones la velocidad de convergencia comparándola con la del método del gradiente se mejora notablemente. Vamos a aproximar las Hessianas procurando que las aproximaciones Bk sean simétricas, denidas positivas y que la dirección dk=−B−1 kgk sea de descenso. Veremos ahora qué condición debe satisfacer la sucesión {Bk} . Sea J:S→R , S⊂Rn abierto, J∈C2(S) . En virtud del teorema de Taylor sobre J(x) J(x)≈J(xk+1)+(∇J(xk+1))T(x−xk+1) + 1 2(x−xk+1)THJ(xk+1)(x−xk+1). y sobre F(x) = ∇J(x) ∇J(x)≈ ∇J(xk+1) + HJ(xk+1)(x−xk+1) y evaluando en x=xk , (HJ(xk+1)−1(∇J(xk+1)−∇J(xk)) ≈xk+1 −xk. Para simplicar la notación tomamos sk=xk+1 −xk y yk=∇J(xk+1)−∇J(xk) . Ahora tenemos HJ(xk+1)−1yk=sk y le pedimos a la aproximación de la inversa de la Hessiana Hk+1 que cumpla Hk+1yk=sk que se llama ecuación o condición cuasi-Newton. Considerando la función mk+1(x) = J(xk+1) + gT k+1(x−xk+1) + 1 2(x−xk+1)TBk+1(x−xk+1) que satisface mk+1(xk+1) = J(xk+1) y ∇mk+1(xk+1) = gk+1 . Pedimos que satisfaga ∇mk+1(xk) = gk , así gk=gk+1 +Bk+1(xk−xk+1). Tenemos entonces que se cumple la condición cuasi-Newton. Si multiplicamos por sT k tenemos sT kBk+1sk=sT kyk . Lo cual signica que si sT kyk>0 , la matriz Bk+1 es denida positiva. Este método usa la aproximación de la inversa de la hessiana, Hk , para calcular el descenso, dk , como −Hkgk . Se da un valor inicial de la aproximación que se actualiza 1.3. MÉTODOS DE RESOLUCIÓN 21 en cada iteración de manera que cumpla la condición cuasi-Newton. Tenemos entonces el siguiente algoritmo del método. Será un algoritmo con la estructura del de un método de descenso, con el descenso calculado como hemos visto y que ademas debe añadir el paso de actualización de la inversa de la hessiana. Algoritmo 1.54. Paso 0 Damos x0∈Rn , H0∈Rn×n (simétrica y denida positiva, por ejemplo I ), 0≤ < 1 e inicializamos k:= 0 . Paso 1 Si ||gk|| ≤  y terminamos. Paso 2 Calculamos dk=−Hkgk . Paso 3 Buscamos el paso óptimo αk>0 y actualizamos xk+1 =xk+αkdk . Paso 4 Actualizamos Hk a Hk+1 de manera que se cumple la condición cuasi-Newton. Paso 5 k:= k+ 1 y volvemos al Paso 1. También se puede usar este algoritmo para encontrar Bk en lugar de Hk . En el Paso 2 se resuelve Bkdk=−gk y en el resto del algoritmo cambiamos Hk por Bk . Veremos ahora dos formas distintas de actualizar Hk . Actualización DFP (Davidon-Fletcher-Powell) Hk+1 se forma añadiendo a Hk dos matrices simétricas, cada una de rango uno. Consideramos la actualización simétrica de rango dos Hk+1 =Hk+auuT+bvvT , donde u, v ∈Rn , a y b escalares. Haciendo que se cumpla la condición cuasi-Newton Hkyk+auuTyk+bvvTyk= sk . u y v no están unicamente determinadas, pero escogeremos u=sk , v=Hkyk por ser las opciones más intuitivas. Entonces tenemos a=1 uTyk=1 sT kyk , b=−1 vTyk=−1 yT kHkyk . Por tanto, Hk+1 =Hk+sksT k sT kyk−HkykyT kHk yT kHkyk (1.7) A la hora de la implementación del algoritmo la actualización en el Paso 4 se hace con la fórmula anterior teniendo en cuenta que sk=αkdk , xk+1 =xk+sk , yk=gk+1 −gk . Teorema 1.55. La actualización 1.7 mantiene las matrices denidas positivas si ,y solo si, sT kyk>0 . 22 CAPÍTULO 1. OPTIMIZACIÓN SIN RESTRICCIONES Actualización BFGS(Broyden-Fletcher-Goldfard-Shanno) Inicialmente haciendo el cambio en la actualización anterior de Hk por Bk y de yk y sk entre sí obtenemos Bk+1 =Bk+ykyT k yT ksk−BksksT kBk sT kBksk . Sean A∈Mn×n(R) invertible, x, y ∈Rn y yTA−1x6=−1 . Aplicando dos veces la fórmula de Sherman-Morrison (A+xyt)−1=A−1−(1 + yTAx)−1A−1xyTA−1 tenemos Hk+1 =Hk+1 + yT kHkyk sT kyksksT k sT kyk−skyT kHk+HkyksT k sT kyk =Hk+(sk−Hkyk)sT k+sk(sk−Hkyk)T sT kyk−(sk−Hkyk)Tyk (sT kyk)2sksT k =I−skyt k sT kykHkI−yksT k st kyk+sksT k sT kyk Capítulo 2 Optimización con restricciones En este capitulo haremos una revisión de los problemas de optimización con restricciones, su formulación, condiciones de optimalidad y hablaremos de uno de sus métodos de resolución que lo reduce a un problema de optimización sin restricciones. Trataremos además el método de los multiplicadores de Lagrange que nos dará una nueva condición necesaria de optimalidad. En primer lugar vamos a centrarnos en el caso en el que solo nos encontramos con restriciones de igualdad para a continuación abordar el caso en el que tenemos también restricciones de desigualdad. 2.1. Restricciones de igualdad Sean J:Rn→R la función objetivo y ci:Rn→R, i = 1; ...,m las funciones de restricción. La formulación general de este tipo de problemas de optimización con restricciones será:    m´ın x∈RnJ(x) sujeto a ci(x)=0, i = 1, . . . , m (2.1) Denición 2.1. Un punto x∈Rn se dice factible si, y solo si, ci(x)=0 para i= 1, . . . , m . El conjunto de todos los puntos factibles denotados por S={x∈Rn|ci(x) = 0 ∀i∈I} siendo I={1, . . . , m} , se denomina conjunto factible . Vamos a suponer S no vacío y reescribiremos el problema como m´ın x∈SJ(x) . Daremos a continuación un par de teoremas que nos serán de utilidad más adelante. Teorema 2.2. Supongamos que ci(x)≥0∀x∈Rn∀i∈I , ci convexa ∀i∈I entonces S es convexo. 23 24 CAPÍTULO 2. OPTIMIZACIÓN CON RESTRICCIONES Demostración. Sean α1≥0, α2≥0 tal que α1+α2= 1 [x∈S, y ∈S]⇒ci(x) = 0 = ci(y)⇒0≤ci(α1x+α2y)≤α1ci(x) + α2ci(y) = 0 ⇒ci(α1x+α2y)=0⇒α1x+α2y∈S Teorema 2.3. Sea c función no nula convexa, entonces c2 es convexa. Demostración. Sean α1≥0, α2≥0 tal que α1+α2= 1 c(α1x+α2y)≤α1c(x) + α2c(y) ⇒c2(α1x+α2y)≤(α1c(x) + α2c(y))2=α2 1c2(x)+2α1α2c(x)c(y) + α2 2c2(y) Usando (a−b)2≥0⇔2ab ≤a2+b2 siendo a=c(x) y b=c(y) tenemos α2 1c2(x)+2α1α2c(x)c(y) + α2 2c2(y)≤α2 1c2(x) + α1α2(c2(x) + c2(y)) + α2 2c2(y) = (α2 1+α1α2)c2(x)+(α2 2+α1α2)c2(y) = α1c2(x) + α2c2(y) Observación 2.4 . En el teorema 2.2 tenemos condiciones sobre las restricciones que de cumplirse el conjunto factible sería convexo. En el teorema 2.3 se habla de la convexidad de las restricciones al cuadrado. Resulta interesante ver que cumpliendo las restricciones las condiciones del teorema 2.3 entonces las restricciones al cuadrado cumplirían las condiciones del teorema 2.2. Teniendo en cuenta además la proposición 1.14 estos dos teoremas serán aplicables después con la denición de la función ψ(.) asociada al método de penalización (ver (2.2)). 2.1.1. Condiciones basadas en los multiplicadores de Lagrange Proporcionamos ahora condiciones necesarias que deben cumplirse en el óptimo para los problemas de optimización con restricciones de igualdad. Denición 2.5. Consideraremos λ∈Rm , la función L(x, λ) = J(x) + m X i=1 λici(x). se denomina el Lagrangiano asociado al problema (2.1) y llamaremos a los λi multiplicadores de Lagrange . 2.1. RESTRICCIONES DE IGUALDAD 25 Observación 2.6 . Vista la denición anterior está claro que ∀x∈Rn tal que ci(x)=0 ∀i∈I tenemos que L(x, λ) = J(x) . Teorema 2.7. Será condición necesaria para que ˆx sea mínimo del problema (2.1) que las derivadas parciales con respecto a todas las variables del Lagrangiano sean cero, ∇L(ˆx, ˆ λ)=0 . Observación 2.8 . Podemos observar que esto aplicado solo a la parciales con respecto a los multiplicadores de Lagrange es equivalente a que se cumplan las restricciones. Cumpliéndose entonces que ci(ˆx)=0 y ∂xJ(ˆx) + Piˆ λi∂xci(ˆx)=0 tendremos un sistema con el mismo número de ecuaciones que de incógnitas tal que si los ci son lineales y los 5ci son linealmente independientes tendremos una única solución. 2.1.2. Método de penalización, caso de restricciones de igualdad En este apartado daremos solo una breve descripción de este método. Nos centraremos más en él cuando a continuación estudiemos el caso más general. Tomando la función ψ(x) = m X i=1 c2 i(x) (2.2) podemos reescribir el conjunto S como S={x∈Rn|ψ(x)=0} . Podemos también reescribir el problema de la siguiente manera m´ın x∈RnJ(x) sujeto a ψ(x) = 0 . Denición 2.9. Sea J:S⊂Rn→R con S no acotado, se dice que es coercitiva sobre S si l´ım |x|→∞ J(x) = +∞ . Teorema 2.10. Consideramos ahora J:Rn→R coercitiva y estrictamente convexa. Sea S⊂Rn no vacío, cerrado y convexo, ψ:Rn→R convexa tal que ψ(x)≥0∀x∈Rn . Entonces ∀η > 0∃ˆxη∈Rn tal que Jη(ˆxη) = m´ın x∈RnJη(x) donde Jη(x) = J(x) + ηψ(x), y además, l´ım η→∞ ˆxη= ˆx , donde ˆx∈S es el único tal que J(ˆx) = m´ın x∈SJ(x) . En este teorema 2.10 vemos que dado el problema de optimización con restricciones de igualdad (2.1) podemos aproximar su solución mediante la solución de un problema de optimización sin restricciones m´ın x∈RnJ(x) + ηψ(x). (2.3) 26 CAPÍTULO 2. OPTIMIZACIÓN CON RESTRICCIONES Si nos jamos bien la manera de construir el funcional objetivo en (2.3), cuya solución aproximará la de (2.1), es añadirle al funcional objetivo de (2.1) un sumatorio de las funciones de restricción al cuadrado multiplicado por un número grande. Vemos entonces que lo que buscamos es obtener un funcional objetivo que sea igual al de (2.1) si se cumplen las condiciones de restricción y considerablemente mayor si no se cumplen. A una función construida de esta manera la llamaremos función de penalización. Veremos ahora el algoritmo de este método para resolver (2.1) Algoritmo 2.11. Paso 0 Damos x0∈Rn , η0>0 , ρ≥0 e inicializamos k:= 0 . Paso 1 Empezando en xk buscamos la solución del problema m´ın x∈RnJ(x) + ηkψ(x) usando alguno de los métodos vistos en el capítulo anterior. Obtenemos xk+1 . Paso 2 Si ||ψ(xk+1)|| ≤ ρ terminamos. Paso 3 Actualizamos ηk+1 = 10ηk , k:= k+ 1 y volvemos al Paso 1. 2.2. Restricciones de igualdad y desigualdad En este caso trabajamos con J(x) la función objetivo y las funciones de restricción ci(x) (i= 1, . . . , m) regulares, de valor real evaluadas en Rn . Sean me y m enteros no negativos tales que 0≤me≤m tenemos la formulación general          m´ın x∈RnJ(x) sujeto a ci(x) = 0, i = 1, . . . , me; ci(x)≥0, i =me+ 1, . . . , m. Simplicamos la notación considerando los conjuntos E={1, . . . , me} I={me+ 1, . . . , m} conjuntos de indices de las restricciones de igualdad y de desigualdad respectivamente. El conjunto factible será así S={x|ci(x)=0, i ∈E;ci(x)≥0, i ∈I} y podemos reescribir el problema como m´ın x∈SJ(x) . Denición 2.12. Sea x∈S , si ci(x) = 0 , i∈I , diremos que ci es una restricción activa en x en caso contrario diremos que ci en una restricción inactiva en x . Siendo I(x) = {i∈I|ci(x)=0} , A(x) = ESI(x) es el conjunto de los índices de las restricciones activas en x . 2.2. RESTRICCIONES DE IGUALDAD Y DESIGUALDAD 27 Supongamos que ˆx es un mínimo local, si hay un índice i0∈I tal que ci0(ˆx)>0 , entonces, si eliminamos la restricción i0 -ésima, ˆx sigue siendo el mínimo local del problema asociado. Vemos entonces que las i0 -ésima restricción es inactiva en ˆx y una motivación para denominarla de esta manera. Sea A(ˆx) el conjunto de indices de las restricciones activas en ˆx , con ˆx mínimo del problema. Entonces será suciente con resolver    m´ın x∈RnJ(x) sujeto a ci(x)=0, i ∈A(ˆx). En general será más fácil resolver el problema de optimización con restricciones de igualdad que el original. El problema es que ˆx no se conoce. 2.2.1. Condiciones de optimalidad de primer orden Empezamos esta sección dando una serie de deniciones que nos resultarán útiles para estudiar después las condiciones de optimalidad. Nos centraremos en dar los enunciados de estas condiciones pasando por alto las demostraciones de las mismas haciendo referencia a dónde se pueden encontrar. Denición 2.13. Sean x∈S , 06=d∈Rn . Si existe δ > 0 tal que x+td ∈S , ∀t∈[0, δ] , entonces d se dice dirección factible de S en x . El conjunto de todas las direcciones factibles de S en x es DF(x, S) = {d|∃δ > 0 tal que x+td ∈S, ∀t∈[0, δ]}. Observación 2.14 . Con restricciones de igualdad en muchos casos no hay direcciones factibles como por ejemplo en optimmización sobre una circunferencia. Denición 2.15. Sea x∈S y d∈Rn . Si dT∇ci(x)=0, i ∈E dT∇ci(x)≥0, i ∈I(x) entonces d se dice dirección factible linealizada de S en x . El conjunto de todas las direcciones factibles linealizadas de S en x es DF L(x, S) = {d|dT∇ci(x) = 0, i ∈E;dT∇ci(x)≥0, i ∈I(x)}. Denición 2.16. Sea x∈S y d∈Rn . Si existen sucesiones {dk} y {δk} , δk>0 tal que x+δkdk∈S, ∀k y dk→d , δk→0 , entonces la dirección límite d de llama dirección 34 CAPÍTULO 2. OPTIMIZACIÓN CON RESTRICCIONES Probaremos el resultado demostrando que para toda sucesión de puntos factibles {xk} tendiendo a ˆx se tiene J(xk)≥J(ˆx) + σ 4||xk−ˆx||2, para todo k sucientemente grande. (2.10) Supongamos que lo anterior no se cumple y existe una sucesión {xk} tendiendo a ˆx tal que J(xk)< J(ˆx) + σ 4||xk−ˆx||2, para todo k sucientemente grande. (2.11) Tomando, si fuese necesario, una subsucesión podemos denir una dirección límite tal que l´ım k→∞ xk−ˆx ||xk−ˆx|| =d . Tenemos entonces xk= ˆx+tkd+o(tk) donde tk=||xk−ˆx|| . Tomando i∈E y usando el teorema de Taylor tenemos 0 = 1 tk ci(xk) = 1 tkci(ˆx) + tk∇ci(ˆx)Td+o(tk)=∇ci(ˆx)Td+o(tk) tk . Haciendo el límite cuando k→ ∞ , el último término de esta expresión tiende a cero y tenemos ∇ci(ˆx)Td= 0 . Para i∈I(x) , tenemos 0≤1 tk ci(xk) = 1 tkci(ˆx) + tk∇ci(ˆx)Td+o(tk)=∇ci(ˆx)Td+o(tk) tk y por argumento análogos ∇ci(ˆx)Td≥0 . Con lo cual d∈DF L(ˆx, S) . Por como denimos el Lagrangiano y dado que ˆ λi≥0 , ci(xk)≥0 para i∈I y ci(xk)=0 para i∈E tenemos L(xk,ˆ λ) = J(xk)−X i∈A(xk) ˆ λici(xk)≤J(xk), (2.12) y al igual que en la demostración anterior tenemos (2.9). Si d no estuviese en G(ˆx, ˆ λ) , podemos encontrar un j∈I(ˆx) tal que se cumple ˆ λj5cj(ˆx)Td > 0, (2.13) mientras que para el resto de indices i∈A(ˆx) tenemos ˆ λi5ci(ˆx)Td≥0. (2.14) Por el teorema de Taylor y l´ım k→∞ xk−ˆx ||xk−ˆx|| =d , para este j tenemos ˆ λjcj(xk) = ˆ λjcj(ˆx)+ˆ λj∇cj(ˆx)T(xk−ˆx)+o(||xk−ˆx||) = ||xk−ˆx||ˆ λj∇cj(ˆx)Td+o(||xk−ˆx||). Entonces por (2.12) L(xk,ˆ λ) =J(xk)−X i∈A(ˆx) ˆ λici(xk)≤J(xk)−ˆ λjcj(xk) ≤J(xk)−||xk−ˆx||ˆ λj5cj(ˆx)Td+o(||xk−ˆx||) (2.15) 2.2. RESTRICCIONES DE IGUALDAD Y DESIGUALDAD 35 Por la estimación basada en el teorema de Taylor (2.9) tenemos L(xk,ˆ λ) = J(ˆx)+O(||xk− ˆx||2) y combinando con (2.15) obtenemos J(xk)≥J(ˆx) + ||xk−ˆx||ˆ λj5cj(ˆx)Td+o(||xk−ˆx||). Esto sería incompatible con (2.11) si se cumpliese (2.13) por tanto tenemos que para toda sucesión de puntos factibles {xk} tendiendo a ˆx debe satisfacer J(xk)≥J(ˆx)+ σ 4||xk−ˆx||2 , para todo k sucientemente grande. Por tanto ˆx es un mínimo local estricto. En el libro [5] encontramos varios ejemplos en los cuales se comprueban las condiciones de segundo orden que resultan interesantes. 2.2.3. Método de penalización Veremos ahora un estudio más amplio del método de penalización centrado en el caso que estamos viendo actualmente. En esta clase de métodos resolveremos en lugar del problema original (2.2) una serie de problemas sin restricciones que minimizan las funciones de penalización. Función de penalización La función de penalización, P(x) = P(J(x), c(x)) , se construye a partir de la función objetivo y de las restricciones y es una función con propiedades de penalización. La propiedad de penalización pide que P(x) = J(x) para todos los puntos factibles y P(x) mucho mayor que J(x) cuando se incumplan las restricciones. Denimos la siguiente función para de alguna manera cuanticar los incumplimientos c∗(x)=(c∗ 1(x), . . . , c∗ m(x))T c∗ i(x) = ci(x), i ∈E; c∗ i(x) = m´ın{ci(x),0}, i ∈I. Para cada restricción la correspondiente función c∗ i es no nula cuando no se cumple y cero cuando sí. La función de penalización consiste en una suma de la función objetivo original y un término de penalización, i. e., P(x) = J(x) + h(c∗(x)) tenemos que h(c∗(x)) será una función denida en Rm y satisface h(0) = 0 , l´ım ||c||→∞ h(c) = ∞ . Un ejemplo será P2(x) = J(x) + η||c∗(x)||2 2 con η > 0 constante, parámetro de penalización. También tendremos P1(x) = J(x) + η||c∗(x)||1 y P∞(x) = J(x) + η||c∗(x)||∞ . 36 CAPÍTULO 2. OPTIMIZACIÓN CON RESTRICCIONES Método de penalización simple Este método permite obtener el mínimo del problema de optimización original (2.2) minimizando una serie de funciones de penalización. Consideremos la función de penalización simple Pη(x) = J(x) + η||c∗(x)||α , donde η > 0 parámetro de penalización destinado a tender a innito (o hacerse grande), α > 0 constante y ||.|| una norma dada en Rm . Escribiremos x(η) solución de m´ın x∈RnPη(x) . Lema 2.37. Sea 0< σ1< σ2 . Entonces Pσ1(x(σ1)) ≤Pσ2(x(σ2)), J(x(σ1)) ≤J(x(σ2)), ||c∗(x(σ1))|| ≥||c∗(x(σ2))||. Demostración. Demostración en el capítulo 10, lema 10.2.1, página 461 de [4]. Lema 2.38. Sea δ=||c∗(x(η))|| . Entonces x(η) es también solución del problema    m´ın x∈RnJ(x) sujeto a ||c∗(x)|| ≤ δ. Demostración. Demostración en el capítulo 10, lema 10.2.2, página 461 de [4]. Por la denición de c∗(x) , el problema original (2.2) se puede reescribir    m´ın x∈RnJ(x) sujeto a ||c∗(x)|| = 0. De manera que si δ es lo sucientemente pequeño obtenemos una aproximación y x(η) se podrá considerar una aproximación de la solución del problema original (2.2). De hecho si c∗(x(η)) = 0 , x(η) será la solución del problema original. La idea básica del método es que el parámetro de penalización η se incremente en cada iteración hasta que ||c∗(x(η))|| sea menor que una tolerancia dada. Algoritmo 2.39. Paso 0 Damos x0∈Rn , η1>0 , ≥0 y inicializamos k:= 1 . Paso 1 Buscamos la solución x(ηk) de m´ın x∈RnPηk(x) empezando en xk . Paso 2 Si ||c∗(x(ηk))|| ≤  terminamos. Paso 3 Actualizamos xk+1 =x(ηk) , ηk+1 = 10ηk ; k:= k+ 1 y volvemos al Paso 1. 2.2. RESTRICCIONES DE IGUALDAD Y DESIGUALDAD 37 Observación 2.40 . Podemos escoger ηk dependiendo de la dicultad de la minimización de la función de penalización en la iteración k-ésima ( ηk+1 =Cηk , C > 1 ). Teorema 2.41. Supongamos que la tolerancia,  del algoritmo satisface  > m´ın x∈Rn||c∗(x)|| entonces el algoritmo debe terminar de manera nita. Demostración. Véase el capítulo 10, teorema 10.2.4, página 463 de [4]. Teorema 2.42. Si el algoritmo no termina de manera nita, entonces m´ın x∈Rn||c∗(x)|| ≥  y l´ım k→∞ ||c∗(x(ηk))|| = m´ın x∈Rn||c∗(x)|| , y cualquier punto de acumulación ˆx de {x(ηk)} es solución del problema      m´ın x∈RnJ(x) sujeto a ||c∗(x)|| = m´ın y∈Rn||c∗(y)||. Demostración. Véase el capítulo 10, teorema 10.2.5, página 463 de [4]. De los dos teoremas anteriores obtenemos Corolario 2.43. Si el problema tiene puntos factibles ( m´ın x∈Rn||c∗(x)|| = 0 ). Entonces o el algoritmo encuentra de manera nita la solución o cualquier punto de acumulación de la secuencia generada será solución del problema original (2.2). 38 CAPÍTULO 2. OPTIMIZACIÓN CON RESTRICCIONES Capítulo 3 Control óptimo de sistemas discretos En este capítulo estudiaremos el control óptimo de sistemas. Empezaremos dando la denición del problema de control para pasar a dar su planteamiento como problemas de optimización con y sin restricciones. Obtendremos las ecuaciones que deben satisfacer los controles óptimos y las relacionaremos con el sistema proporcionado por los multiplicadores de Lagrange en el capítulo anterior. Veremos algún método de resolución de problemas de control que será implementado y usado para resolver problemas académicos. Por último se aplicarán estas técnicas a problemas de interés tales como la resolución de un problema de control en EDO tras haber sido discretizado o el entrenamiento de una red neuronal de clasicación. 3.1. Planteamiento del problema Con frecuencia nos enfrentamos a problemas de minimización m´ın u∈S ˜ J(u) donde la expresión del funcional ˜ J es compleja pero que admite una reescritura ˜ J(u) = J(y(u), u) donde y(u) es la solución de F(y(u), u)=0. (3.1) En estos casos u∈Rn recibe el nombre de variable de control mientras que y∈Rm es la llamada variable de estado denida por la ecuación de estado (3.1). Tendremos F:Rm×Rn→Rm diferenciable, denotaremos por yu≡y(u) a la única (por hipótesis) 39 40 CAPÍTULO 3. CONTROL ÓPTIMO DE SISTEMAS DISCRETOS solución de la ecuación. Así el problema será    m´ın u∈Rn ˜ J(u) := J(yu, u) donde yu es la única solución de F(yu, u)=0, (3.2) con ˜ J también diferenciable. Vectorialmente tenemos y= (y1, y2,...ym), F(y, u)=(f1(y, u), f2(y, u), . . . , fm(y, u))T identicaremos fi como una función que dependerá de la variable de control y de la variable de estado. 3.1.1. Planteamiento como problema de optimización con restriciones El problema (3.2) se puede ver como el siguiente problema de optimización con restricciones de igualdad,      m´ın (y,u)J(y, u) sujeto a F(y, u)=0 (3.3) Con la nueva formulación (3.3) será condición necesaria que las parciales del Lagrangiano con respecto a todas la variables se anulen (Teorema 2.7). λ∈Rm será el vector que contiene a los multiplicadores de Lagrange, en el marco de la teoría de control óptimo recibe el nombre de estado adjunto, y el Lagrangiano será L(y, u, λ) = J(y, u)−λTF(y, u). Así tendremos        dλL(y, u, λ) = −F(y, u)=0 Ecuación de estado dyL(y, u, λ) = dyJ(y, u)−λTdyF(y, u)=0 Sistema adjunto duL(y, u, λ) = duJ(y, u)−λTduF(y, u) = 0 Condición de optimalidad (3.4) La terminología usada aquí se entenderá más tarde. Una manera más eciente de ver esto es de manera vectorial y matricial, que nos ayudará sobre todo después a entender el método de resolución y a aplicarlo a casos más concretos. La ecuación de estado la podremos ver como F(y, u) = 0 y teniendo en cuenta su formulación vectorial se tendrá que cumplir fi(y, u)=0 ∀i∈ {1, . . . , m}. (3.5) 3.1. PLANTEAMIENTO DEL PROBLEMA 41 El sistema adjunto lo podemos ver como dyFT(y, u)λ=dyJT(y, u) . Obviamos la dependencia de y y de u para simplicar la notación y tenemos dyF(y, u) =         ∂y1f1∂y2f1. . . ∂ymf1 ∂y1f2∂y2f2. . . ∂ymf2 . . .. . ..... . . ∂y1fm∂y2fm. . . ∂ymfm         dyJ(y, u)=(∂y1J, ∂y2J, . . . , ∂ymJ) Entonces podemos expresar dyFT(y, u)λ=dyJT(y, u) como         ∂y1f1∂y2f1. . . ∂ymf1 ∂y1f2∂y2f2. . . ∂ymf2 . . .. . ..... . . ∂y1fm∂y2fm. . . ∂ymfm                 λ1 λ2 . . . λm         =         ∂y1J ∂y2J . . . ∂ymJ         (3.6) En cuanto a la condición de optimalidad tenemos que la podemos expresar como duJT(y, u)−duFT(y, u)λ= 0 vectorialmente y omitendo como antes la dependencia tenemos duJ(y, u) = (∂u1J, ∂u2J, . . . , ∂unJ) duF(y, u) =         ∂u1f1∂u2f1. . . ∂unf1 ∂u1f2∂u2f2. . . ∂unf2 . . .. . ..... . . ∂u1fm∂u2fm. . . ∂unfm         Entonces tendremos         ∂u1J ∂u2J . . . ∂unJ         −         ∂u1f1∂u1f2. . . ∂u1fm ∂u2f1∂u2f2. . . ∂u2fm . . .. . ..... . . ∂unf1∂unf2. . . ∂unfm                 λ1 λ2 . . . λm         =         0 0 . . . 0         (3.7) Observación 3.1 . Tomamos todas las fórmulas traspuestas con respecto a lo que aparece en (3.4) para que nos coincidan con como vamos a calcular el gradiente de ˜ J en la sección 3.1.2 donde se reformulará el problema como un problemas de optimización sin restricciones. 42 CAPÍTULO 3. CONTROL ÓPTIMO DE SISTEMAS DISCRETOS 3.1.2. Planteamiento como problema de optimización sin restricciones También podemos ver el problema (3.2) como un problema de optimización sin restricciones. En este caso habría que obtener en primer lugar yu resolviendo F(yu, u) = 0 y tomando luego el problema m´ın u ˜ J(u) = J(y, u) (3.8) En este caso tenemos ahora que una condición necesaria será ∇˜ J= 0 . A parte de por esto nos resultará de especial interés el cálculo del gradiente ya que para la mayoría de métodos de resolución para problemas de optimización sin restricciones vistos en el capitulo 1 es necesario. Veremos dos formas distintas de calcular ∇˜ J . Cálculo del gradiente con el linelizado Ésta se podría considerar la forma más directa de hacer el cálculo. Tenemos que ∇˜ J(u) = du˜ JT(u) calcularemos entonces du˜ J(u) = dyJ(yu, u)duyu+duJ(yu, u) (3.9) en donde todo se obtendría de manera relativamente fácil a excepción de duyu . Para calcular duyu se empieza derivando implícitamete la ecuación de estado (3.1) con lo cual obtenemos dyF(yu, u)duyu+duF(yu, u)=0 , que nos proporciona el sistema linealizado dyF(yu, u)duyu=−duF(yu, u). (3.10) En este sistema podemos obtener dyF(yu, u) y duF(yu, u) , por tanto resolviendo la ecuación conseguiríamos duyu . Tendríamos entonces todo lo necesario para el cálculo de (3.9) y por tanto también para obtener el gradiente. En ocasiones puede resultar muy costosa la resolución de (3.10) en estas ocasiones resultará más adecuado calcular el gradiente como explicamos a continuación. Cálculo del gradiente con el adjunto Tenemos al igual que antes que ∇˜ J(u) = du˜ JT(u) por tanto querremos de nuevo usar (3.9) para calcular du˜ J(x) pero en este caso abordaremos el problema con duyu de manera diferente. Buscaremos trasformar dyJ(yu, u)duyu en algo que resulte más fácil calcular ( duJ(yu, u) ya aparece en la fórmula (3.7) y por tanto querremos substituir el otro sumando por −λTduF(yu, u) para que coincida). 3.2. MÉTODO DE RESOLUCIÓN 43 Usamos en primer lugar el sistema linealizado 3.10. Multiplicamos escalarmente este sistema por λT∈Rm arbitrario y obtenemos λTdyF(yu, u)duyu=−λTduF(yu, u). La idea es considerar un λ tal que coincidan dyJ(yu, u)duyu y λTdyF(yu, u)duyu para reemplazarlos por −λTduF(yu, u) . Tomamos λ denido por dyFT(yu, u)λ=dyJT(yu, u) que coincide con el sistema adjunto y es un sistema lineal. En este caso, llamando λu≡λ(u) a la solución del sistema obtenemos du˜ J(u) = −λT uduF(yu, u) + duJ(yu, u) y por tanto ∇˜ J(u) = du˜ JT(u) = −duFT(yu, u)λu+duJT(yu, u). (3.11) Tenemos entonces que podemos calcular el gradiente resolviendo la ecuación de estado (3.1) para calcular y , el sistema (3.6) para calcula λ y el gradiente con (3.11). 3.2. Método de resolución Para la resolución del problema resultará interesante considerar el planteamiento como problema de optimización sin restricciones. Usaremos entonces los métodos vistos en el capítulo 1 para este tipo de problemas calculando el gradiente de alguna de las formas vistas en la sección anterior. De forma general tendremos que calcular el gradiente con uno de los dos métodos anteriores para después pasar a usar uno de los métodos de resolución del capítulo 1. Además dentro de los métodos tendremos que calcular el paso en cada iteración para lo cual también tendremos que escoger una de las reglas vistas en la correspondiente sección. Vamos entonces a crear un programa que nos permita solucionar el problema. Buscaremos que sea aplicable a casos académicos con los que comprobaremos su funcionamiento. Nos permitirá escoger como calcular el gradiente (linealizado o adjunto), que método usar (gradiente o gradiente conjugado) y como calcular el paso (Armijo, Wolfe o paso jo). Debido a la longitud de este programa a continuación se incluye solo el código de uno de los casos que hemos hablado. El código completo se tendrá como material adicional. Por tanto vemos ahora un programa que resuelve el problema de control calculando el gradiente con el linealizado, usando el método del gradiente conjugado y la regla de Armijo. Código control óptimo (linelizado, gradiente conjugado y Armijo) 1 function [u ,y , j , conv , ng ,d]=conopt_lin_gradconj_arm(n ,m, J ,duJ ,dyJ , . . . 50 CAPÍTULO 3. CONTROL ÓPTIMO DE SISTEMAS DISCRETOS adjunto (3.7). La dependencia de las funciones fi nos dará una estructura muy particular de la matriz que interviene en (3.6) y de las funciones de (3.5). Ahora bien, lo que vamos a usar es el método del gradiente y la variable en la cual vamos a buscar el mínimo sera la variable de control u . Por tanto lo primero que tendremos que hacer será dar u0 inicial, el cual usaremos para construir la sucesión que convergerá al mínimo. No usaremos uk para nombrar a los elemento de la sucesión pues esta es la notación que estamos usando para lo elementos del vector u . En su lugar usaremos uk para indicar que estamos hablando del valor de u en concreto en la iteración k y actualizaremos sobre u entendiendo que esto creará una sucesión convergente a la solución. Aclarado esto podemos pasar a la parte más especica del método. En cada iteración k , por como tenemos denida la función F(y, u) , sabemos que f1 solo dependerá de y1 y u . Al tener dado el valor de u como uk tendremos que F(y, uk) = 0 en su primera coordenada, f1(y1, uk) = 0 , nos dará el valor de y1 , yk 1 . Esto tiene sentido puesto que en el instante inicial tendremos que tener el valor de la variable de estado. No lo podremos dar en función de los anteriores pues no hay ninguno. A partir de esto, por como tenemos denidos los fi podemos obtener de manera recurrente los yk i pues tenemos el valor de la variable de control y los valores de la variable de estado el los instantes anteriores. Obtendremos así todos los valores de la variable de estado yk . A continuación tendremos que calcular λk y lo calcularemos como la solución de (3.6). Ahora bien teniendo en cuenta de nuevo como asumimos que estaban construidas las funciones fi tendremos ∂ytfl= 0 para t>l . Esto nos dará una construcción de la ecuación (3.6) que nos permitirá calcular los λi de manera recurrente retrógrada. Tendremos la representación matricial         ∂y1f1∂y1f2. . . ∂y1fm 0∂y2f2. . . ∂y2fm . . ........ . . 0. . . 0∂ymfm                 λ1 λ2 . . . λm         =         ∂y1J ∂y2J . . . ∂ymJ         de aquí obtendremos las ecuaciones ∂ymfmλm=∂ymJ . . . ∂y2f2λ2+···+∂y2fmλm=∂y2J ∂y1f1λ1+∂y1f2λ2+···+∂y1fmλm=∂y1J (3.12) en las cuales se observa de manera clara la recurrencia retrógrada de la que habíamos hablado. Teniendo en cuenta que ∂ytfl y ∂ytJ estarán evaluadas en uk y yk podremos 3.3. APLICACIONES 51 calcular λk y tendremos todo lo necesario para calcular el gradiente 5˜ J(uk) mediante la fórmula (3.7). Una vez hecho esto podremos comprobar si se cumple el criterio de convergencia ||5 ˜ J(uk)|| <  para un cierto  > 0 , tolerancia dada. Si el criterio de convergencia no se cumple podremos usar alguno de los métodos vistos en el primer capítulo para calcular el paso, αk y actualizaremos uk+1 =uk−α5˜ J(uk) ; k=k+1 . Después de esto volveremos a empezar el proceso calculando los yk i . Veremos ahora el algoritmo de manera más esquematizada para que nos quede más claro. Algoritmo 3.2. Paso 0 Damos u0 ,  > 0 la tolerancia e inicializamos k:= 0 . Paso 1 Calculamos los yk i de manera recurrente con las funciones fi(yi, ..., y1, uk) igualadas a cero. Paso 2 Calculamos λk de manera recurrente retrograda con las ecuaciones 3.12. Paso 3 Calculamos 5˜ J(uk) con la fórmula 3.7. Paso 4 Si se cumple ||5 ˜ J(uk)|| <  terminamos y uk es el mínimo. Paso 5 Calculamos el paso, αk , de alguna de las formas vistas en el capítulo 1. Paso 6 Actualizamos uk+1 =uk−αk5˜ J(uk) , k:= k+ 1 y volvemos al Paso 1. 3.3.2. Aprendizaje de una red neuronal articial Vamos a dar ahora un ejemplo de uso del control óptimo de sistemas en el aprendizaje de redes neuronales articiales. Este ejemplo está sacado de el artículo [9]. Buscamos en primer lugar dar una descripción del problema general del aprendizaje de una red neuronal articial relacionándolo con lo visto anteriormente. Finamente damos un código basado en el que aparece en el artículo pero generalizado. Consideramos un conjunto de puntos divididos en dos categorías, A y B . Lo que queremos es construir una función que dado un punto cualquiera nos lo clasique dentro de una de estas categorías. Basaremos la red en la función σ(x) = 1 1+e−x . Podemos ver la función σ(x) como una versión regular de una función salto, que replica el comportamiento de una neurona en el cerebro, activa (1) si la entrada es suciente e inactiva (0) en otro caso. Además es conveniente que su derivada sea de la forma σ0(x) = σ(x)(1 −σ(x)) . El escalón y su posición en la función sigma se puede alterar escalando y cambiando el argumento o, en lenguaje de redes neuronales añadiendo un peso y desplazamiento a la entrada. Para 52 CAPÍTULO 3. CONTROL ÓPTIMO DE SISTEMAS DISCRETOS mantener una notación sencilla se dene σ:Rm→Rm en z∈Rm de manera vectorial siendo σi(z) = σ(zi) . Tenemos ahora una serie de capas de neuronas. En cada capa, cada neurona devuelve un único numero real, que se pasa a cada neurona de la capa siguiente. Ahora cada neurona en esta capa crea su propia combinación de estos valores (operación afín) con un peso y desplazamiento propios y aplica la función sigma. Matemáticamente tenemos que si a es el vector de los números reales producidos por una capa el vector de salida de la siguiente tendrá la forma σ(Wa +b) . W será una matriz que contiene los pesos y b un vector con los desplazamientos. Sean ni el número de neuronas de una capa y ni+1 el de la capa siguiente. Tendremos a∈Rni pues recoge los valores en las neuronas de la primera capa. W∈Rni+1×ni pues multiplicará a a y la dimensión del vector resultado de la multiplicación deberá ser el número de neuronas de la segunda capa. Además la dimensión de b también debe coincidir con el número de neuronas de la segunda capa, por tanto b∈Rn+1 . Aplicando esto ahora a un caso en el que tenemos L capas hay que distinguir dos capas que se comportarán de manera diferente. La primera capa, la de entrada no tendrá una capa que le pase sus valores, contendrá los datos de entrada y por tanto el número de neuronas de esta capa coincidirá con la dimensión de estos. En la capa de salida, la L - 3.3. APLICACIONES 53 ésima no hay capa siguiente y estas neuronas proporcionan la salida nal. Suponemos que la capa l , para l= 1,2,3, . . . , L , contiene nl neuronas. Usamos Wl∈Rnl×nl−1 para denotar la matriz de pesos en la capa l . Especícamente Wl jk es el peso que la neurona j en la capa l aplica a la salida de la neurona k en la capa l−1 . Analogamente, bl∈Rn l es el vector de desplazamientos para la capa l , así la neurona j de la capa l usa el desplazamiento bl j . Dada una entrada x∈Rn1 , podemos resumir la acción de la red siendo al j la salida, o activación, de la neurona j en la capa l . Así tenemos a1=x∈Rn1 (3.13) al=σ(Wlal−1+bl)∈Rnl para l= 2,3, . . . , L. (3.14) Estás ecuaciones serán la relacionadas con el sistema de estado y son las que nos permitirán calcular las variables de estado de manera recurrente. Queda claro que con ellas y un algoritmo para llevar la entrada a través de la red obtendremos la salida aL∈RnL . Lo que nos interesará será tener unos pesos y desplazamientos que hagan que las salidas nos clasiquen correctamente los puntos que introducimos, es decir, tener salida (1,0) si el punto está en la categoría A y (0,1) si está en la B. Ahora suponemos que tenemos N datos, o puntos dato, en RnL , {xi}N i=1 , para los que nos dan objetivos de salida {y(xi)}N i=1 en RnL , donde y(xi) =                  1 0  si xi está en la categoría A   0 1  si xi está en la categoría B. Nuestra función de coste toma entonces la forma Cost =1 NPN i=1 1 2||y(xi)−aL(xi)||2 2 en donde para no alargarnos con la notación no se especica que el coste es función de los pesos y desplazamientos. Aquí se introduce el factor 1 2 por conveniencia. Así tendríamos un problema de control óptimo. Las variables de control serán los pesos y desplazamientos, la variable de estado la salida aL , el funcional objetivo será la función coste y aL la solución del sistema compuesto por (3.13) y (3.14), que nos dará el sistema de estado. Para optimizar querremos usar 5Cost pero esto puede suponer un gran coste computacional al contener un sumatorio a lo largo de los puntos dato. Usamos entonces el gradiente estocástico y en cada iteración escogeremos un punto xi aleatoriamente y substituiremos la función Cost por uno de sus sumandos Cxi=1 2||y(xi)−aL(xi)||2 2 para simplicar tomaremos C=1 2||y−aL||2 2 evitando la dependencia de xi . Vamos a calcular a continuación derivadas parciales con respecto a la variables de control, pero para simplicar la notación es interesante introducir dos nuevos conjuntos 54 CAPÍTULO 3. CONTROL ÓPTIMO DE SISTEMAS DISCRETOS de variables. zl=Wlal−1+bl∈Rnl para l= 2,3, . . . , L , donde zk j es la entrada escalada de la neurona j en la capa l . Así, al=σ(zl) con l= 2,3, . . . , L y δl∈Rnl denido como δl j=∂C ∂xl j ,1≤j≤nl,2≤l≤L . Esta expresión que se acostumbra a llamar error de la j -ésima neurona de la capa l es una cantidad intermedia que es útil en el análisis y en la computación. El termino error en este caso no es el más acertado pero cobra sentido porque la función coste solo purde minimizarse si todos los δl j son cero (podremos observar que las variables δ son el estado adjunto en la terminología de control óptimo). Este paso se correspondería al del cálculo del estado adjunto usando el sistema adjunto, debido a las particularidades de este ejemplo no se corresponde exactamente con lo visto en el caso teórico. Veamos entonces como calcular estos δl , que se corresponderán con los multiplicadores de Lagrange. Tenemos que aL=σ(zL) entonces ∂aL j ∂zL j =σ0(zL j) y por la denición de C , ∂C ∂aL j =∂ ∂aL j 1 2 nL X k=1 (yk−aL k)2=−(yj−aL j). Usando la regla de la cadena tenemos δL j=∂C ∂zL j =∂C ∂aL j ∂aL j ∂zL j = (aL j−yj)σ0(zL j) . Usamos de nuevo la regla de la cadena, esta vez para pasar de la derivada con respecto a zj j a las parciales con respecto a {zl+1 k}k= 1nl+1 . δl j=∂C ∂zl j = nl+1 X k=1 ∂C ∂zl+1 k ∂zl+1 k ∂zl j = nl+1 X k=1 δl+1 j ∂zl+1 k ∂zl j . Ahora como zl=Wlal−1+bl sabemos que zl+1 k y zl j están conectados mediante zl+1 k= Pnl s=1 Wl+1 ks σ(zl s) + bl+1 k . Por tanto ∂zl+1 k ∂zl j =Wl+1 kj σ0(zl j) y tenemos δl j= nl+1 X k=1 δl+1 jWl+1 kj σ0(zl j). La fórmula anterior la podemos reajustar como δl j=σ0(zl j)(Wl+1)Tδl+1. 3.3. APLICACIONES 55 Podemos ver entonces el sistema adjunto como δL j= (aL j−yj)σ0(zL j), (3.15) δl j=σ0(zl j)(Wl+1)Tδl+1. (3.16) y observamos como efectivamente nos permite calcular los δl de manera recurrente retrograda, como hacíamos con los multiplicadores de Lagrange. Una vez calculados los δl j pasamos al cálculo de las derivadas parciales de la función coste con respecto a las componentes de las variables de control. Con estas derivadas construiremos los gradientes con respecto a cada una de las variables de control que después usaremos en el método del gradiente estocástico con paso jo. Sabemos que zl j está conectado a bl j por zl j= (Wlσ(zl−1))j+bl j. Tenemos que ∂zl j ∂bl j = 1 . Entonces usando la regla de la cadena ∂C ∂bl j =∂C ∂zl j ∂zl j bl j =∂C ∂zl j =δl j. Finalmente tenemos zl j=Pnl−1 k=1 Wl jkal−1 k+bl j , que nos da ∂zl j ∂W l jk =al−1 x , independientemente de j y ∂zl s ∂W l jk = 0 para s6=j . Entonces aplicando la regla de la cadena de nuevo tenemos ∂C ∂Wl jk = nl X s=1 ∂C ∂zl s ∂zl s ∂Wl jk =∂C ∂zl j ∂zl j ∂Wl jk =δl jal−1 k. Resumiendo tenemos: δL j= (aL j−yj)σ0(zL j), δl j=σ0(zl j)(Wl+1)Tδl+1j, ∂C ∂bl j =δl j, ∂C ∂Wl jk =δl jal−1 k. Introduciendo ahora la matriz Dl∈Rnl×nl , matriz diagonal con σ0(zl i) en la entrada (i, i) podemos reescribir δL=DL(aL−y), δl=DlWl+1Tδl+1. Es interesante recordar además que cuando denimos la función σ vimos que su derivada era fácilmente calculable. 56 CAPÍTULO 3. CONTROL ÓPTIMO DE SISTEMAS DISCRETOS Con lo visto estamos ya en posición de dar un pseudocódigo usando el gradiente estocástico (Observación 1.39) con paso jo η . For counter = 1 a Nit k escogido aleatoriamente en {1, . . . , N} xk dato de entrenamiento a1=xk For l =2:L zl=Wlal−1+bl al=σ(zl) Dl=diag(σ0(zl)) end For l =L−1 : 2 δl=Dl(Wl+1)Tδl+1 end For l =L: 2 Wl→Wl−ηδl(al−1)T bl→bl−ηδl end end A continuación tenemos el código en Matlab, que como dijimos antes es una generalización del que se da el el artículo citado. Aprendizaje de una red neuronal articial 1 function [W,b]=rna_adjunto (nxc , clas , x1 , x2 , N_it , eta ) 2 % rna_adjunto (nxc , clas , x1 , x2 , N_iter , eta ) 3 % Usa el metodo del adjunto para entrenar una red neuronal a r t i f i c i a l 4 % 5 % Vamos a introducir el numero de capas y el numero de neuronas en cada 6 % capa para que sea variable en el programa , las coordenadas de los 7 % elementos distinguiendo a que categoria pertenecen , el numero maximo de 8 % iteraciones que queremos que ejecute el programa y la tasa de 9 % aprendizaje. 3.3. APLICACIONES 57 10 % 11 % nxc : Vector de longitud el numero de capas de nuestra RNA y que 12 % contiene en cada elemento i el numero de neuronas de la capa 13 % i . 14 % clas : Vector de longitud el numero distinto de clases a las que pueden 15 % pertenecer los puntos dato . Cada elemento nos da el numero de 16 % puntos dato en cada clase . 17 % x1 : Vector con las primeras componentes de los puntos dato . Su 18 % longitud debe coincidir con la suma de los elemntos en clas , los 19 % puntos deben estar ordenados por clase y las clases en el mismo 20 % orden que el el vector clas . De manera que los clas (1) primeros 21 % elementos de x1 correspondan con los puntos dato de la clase "1" 22 % los clas (2) siguentes con los de la "2" y asi sucesivamente . 23 % x2 : Vector con las segundas coordenadas de los puntos dato . Debe estar 24 % construido igual que x1 de forma que elemento i de x2 sea la 25 % segunda coordenada del punto al que corresponde la primera 26 % coordenada almacenada en x1( i ). 27 % N_it : Numero de iteraciones que queremos llevar a cabo . 28 % eta : Tasa de aprendizaje (tamanho de paso ). 29 30 31 %Matriz de cla si ficacio n de los elementos 32 y = zeros ( length ( clas ) , length (x1 ) ); 33 u=0; 34 35 for j =1: length ( clas ) 36 l=u+1; 37 u=u+clas ( j ); 38 y( j , l :u)=1; 39 end 40 41 % Inicializacion de los pesos y desplazamientos 42 rng (5000); 43 capas= length (nxc ); 44 45 w=c el l ( capas = 1 ,1); 58 CAPÍTULO 3. CONTROL ÓPTIMO DE SISTEMAS DISCRETOS 46 for j =1:capas = 1 47 w( j )={0.5 * randn (nxc( j +1),nxc( j ))}; 48 end 49 b=c ell ( capas = 1,1); 50 for j =1:capas = 1 51 b( j )={0.5 * randn (nxc( j +1) ,1)}; 52 end 53 54 %Propagacion Forward and Backward 55 savecost = zeros (N_it , 1); %Valor de la funcion coste en cada iteracion 56 N= length (x1 ); 57 for counter = 1:N_it 58 k = randi (N); %Escoger un punto aleatorio para la evaluacion 59 x = [ x1(k ); x2(k ) ] ; 60 %Paso Forward 61 a=cell ( capas = 1,1); 62 x_act=x ; 63 for j =1:capas = 1 64 a( j)={activacion (x_act , cell2mat (w( j )) , cell2mat (b( j )))}; 65 x_act=cell2mat (a( j ) ); 66 end 67 68 %Paso Backward 69 delta=c ell (capas = 1,1); 70 delta ( capas = 1)={cell2mat (a(capas = 1)). * (1 = cell2mat (a( capas = 1))) . . . 71 . * ( cell2mat (a(capas = 1)) = y (: , k ))}; 72 for j=capas = 2: = 1:1 73 delta ( j ) = {cell2mat (a( j )). * (1 = cell2mat (a( j ))) . . . 74 . * ( cell2mat (w( j +1)) ' * cell2mat ( delta ( j +1)))}; 75 end 76 77 %Actualizacon metodo del gradiente 78 w(1) = {cell2mat (w(1)) = eta * cell2mat ( delta (1)) * x '} ; 79 b(1) = {cell2mat (b(1)) = eta * cell2mat ( delta (1))}; 80 for j =2:capas = 1 81 w( j ) = {cell2mat (w( j )) = eta * cell2mat ( delta ( j )) * cell2mat (a( j = 1)) '}; 3.3. APLICACIONES 59 82 b( j ) = {cell2mat (b( j )) = eta * cell2mat ( delta ( j ))}; 83 end 84 85 %Visualizacion y gardado del progreso 86 newcost = cost (w,b) 87 savecost ( counter ) = newcost ; 88 end 89 %Muestra la progresion del vector coste 90 save costvec 91 semilogy (1:1 e4 : N_it , savecost (1:1 e4 : N_it)) 92 function costval = cost (w,b) 93 costvec = zeros (10,1); 94 for i = 1:10 95 x =[x1( i ); x2( i ) ] ; 96 x_act=x ; 97 for jk=1:capas = 1 98 a( jk)={activacion (x_act , cell2mat (w( jk )) , cell2mat (b( jk )))}; 99 x_act=cell2mat (a( jk )) ; 100 end 101 costvec ( i ) = norm (y (: , i ) = cell2mat (a( capas = 1)),2); 102 end 103 costval = norm ( costvec ,2)^2; 104 end 105 end En está función se usa además la función de activación que será. Función de activación 1 function y = activacion (x ,W, b) 2 %y=activacion (x ,W, b) 3 %Evalua la funcion sigma . 4 % 5 % x es el vector de entrada , y es el vector de salida 6 %W contiene los pesos , b contiene los shifts 7 % 8 % La componente i de y es activacion ((Wx+b)_i) 9 % donde activacion (z) = 1/(1+exp( = z )) 66 CAPÍTULO 3. CONTROL ÓPTIMO DE SISTEMAS DISCRETOS Bibliografía [1] Álgebra lineal y geometría , 3 ª edición, Eugenio Hernández Rodríguez, María Jesús Vázquez Gallo, María Ángeles Zurro Moro, PEARSON EDUCATION, S.A., Madrid, 2012. [2] Álgebra lineal y geometría cartesiana , 3 ª edición, Juan de Burgos Román, McGranw- Hill, 2006. [3] Apuntes de Análisis Matemático I , María D. Acosta, Camilo Aparicio, Antonio Moreno, Universidad de Granada. [4] Optimization theory and methods. Nonlinear programming , volumen 1, Wenyu Sun, Ya-Xiang Yuan, Springer, 2006. [5] Numerical Optimization , Nocedal, J., Wright, S. J., Springer-Verlag, 1999. [6] Optimización dinámica , E. Cerdá Tena, Prentice Hall, 2001. [7] Ingeniería de control moderna , K. Ogata, Pearson-Prentice-Hall, 2010. [8] Commande optimale. Notes du cours A08 , E. Trélat, 2007/2008. [9] Deep Learning: An Introduction for Applied Mathematicians de Catherine F. Higham y Desmond J. Higham publicado en Siam Review Vol. 61, No. 4, pp. 860891 67