Algoritmos de primer orden para implementación del Control Predictivo para Seguimiento

Trabajo de Fin de Grado Ingeniería de Tecnologías Industriales Algoritmos de primer orden para implementación del Control Predictivo para Seguimiento Autor: Javier Maestro Valenzuela Tutor: Teodoro Álamo Cantarero Sevilla, 2021 iii Trabajo de Fin de Grado Ingeniería de Tecnologías Industriales Algoritmos de primer orden para implementación del Control Predictivo para Seguimiento Autor: Javier Maestro Valenzuela Tutor: Teodoro Álamo Cantarero Profesor Catedrático Dpto. de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2021 v Trabajo de Fin de Grado: Algoritmos de primer orden para implementación del Control Predictivo para Seguimiento Autor: Javier Maestro Valenzuela Tutor: Teodoro Álamo Cantarero El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2021 El Secretario del Tribunal vii A mi familia A mis maestros ix Agradecimientos En primer lugar, me gustaría dedicar este apartado a todos y cada uno de los profesores que he tenido en esta universidad. Todos ellos ofrecieron su esfuerzo, dedicación y sabiduría para hacerme uno de los mayores regalos que una persona puede hacerle a otra, dar conocimiento. En particular de entre la larga lista de profesores, quiero hacer un mención muy especial al tutor de este proyecto, Teodoro me ha ofrecido más esfuerzo, dedicación y sabiduría que ningún otro. Y mucha paciencia. Por último, a mi familia, amigos y compañeros, los que me han acompañado en el día a día y sin los que este trabajo no habría sido redactado. ÍNDICE DE FIGURAS Ilustración 1-1: Conjunto convexo. Ilustración 1-2: Conjunto no convexo. Ilustración 3-1: FISTA. Ilustración 4-1: Dinámica MPC. xvii Notación ()dom h Dominio de h ()epi h Epígrafo de h n Espacio real de n dimensiones min max argmin Mínimo de una función Máximo de una función Argumento mínimo h Norma de h  s.a. Perteneciente a Sujeto a : x y    Tal que Derivada de x respecto de y Para todo  Menor o igual  Mayor o igual ⇔ Si y sólo si 1 1 PROBLEMA DE OPTIMIZACIÓN CONVEXA omienzo definiendo el objeto de estudio de este primer capítulo: la función convexa. Para ello emplearé una serie de descripciones sobre los diferentes conceptos que han de estar implicados en la definición. Por una parte, está el concepto de convexidad, qué es y qué conlleva que un conjunto sea convexo. Por otro lado, la presentación del problema general, que será un problema de minimizar una función convexa. Todos estos conceptos son bien conocidos, ya que se profundizado mucho en ellos en trabajos como [4],[7],[9]. 1.1 Definición de convexidad La convexidad es una propiedad extremadamente utilizada y desarrollada en el campo de la optimización. Antes de pasar a la función es necesario explicar un concepto fundamental, el conjunto convexo. Definición 1-1. Un conjunto n S será convexo si para cualquier par de elementos de S (x, y) : (1 ) , [0,1]x y S    − +    . Ilustración 1-1 1 Ilustración 1-2 2 Se puede demostrar, esto será muy útil más adelante, que un conjunto definido por   :S x Ax b   , es convexo. Sean dos elementos de S, ,x y S por lo que cumplen la definición del conjunto y [0,1]Ax b     1 Representación gráfica de un conjunto que satisface las condiciones de convexidad. 2 Representación gráfica de un conjunto que no satisface las condiciones de convexidad. C Problema de optimización convexa 2 (1 ) (1 )Ay b  −  − . Sumando ambas, tenemos (1 ) ( (1 ) ) . Ax Ay b A x y b   + −  + −  La expresión anterior es la definición de conjunto convexo. Ahora paso a definir conceptos relativos a la convexidad de una función, la que será objeto de este capítulo. Definición 1-2. El dominio de una función lo definimos como el conjunto de puntos para los que la función toma valores finitos: ( ) { :| ( )| } n dom h x h x=  + . Una función será convexa si su dominio es convexo y para todo [0,1]   se cumple que: ( (1 ) ) ( ) (1 ) ( ), ( ), ( )h x y h x h y x dom h y dom h     + −  + −     . El epígrafo de una función se define como el conjunto: ( ) {(x,t) dom(h) :h(x) t}epi h =    . Se concluye que una función será convexa si y sólo si su epígrafo es un conjunto convexo. 1.2 Definición del problema de optimización convexa En segundo lugar, presento aquí el problema de optimización de una función convexa que vamos a abordar: * min ( ) min ( ) x X x X f f x h x  == . 3 3 Algoritmos de primer orden para implementación del Control Predictivo para Seguimiento Trataré de encontrar el mínimo de la función objetivo ()fx , teniendo en cuenta una serie de restricciones que iré definiendo más adelante. Sólo se tratarán problemas de minimización y con restricciones de igualdad y restricciones en caja, debido a que cualquier problema de maximización puede convertirse en un problema de minimización, no se pierde la generalidad. Asumiendo lo siguiente: 1. ()hx es una función convexa diferenciable suave. Esto es, hay una 0R tal que se cumple que: 2 1 ( ) ( ) ( ), , , 2nn R h x h y h y x y x y x y + − + −     . 2. El problema de minimización: min ( ) xXfx  . es resoluble. Esto es, hay una * ( )x X dom   tal que * ( *) inf ( ) xX f f x f x  == . Para que el problema de optimización sea convexo es necesario, por una parte, que la función ()fx sea convexa, cosa que se afirma en la asunción 1. Por otra parte, que la región factible sea un conjunto convexo, esto se cumplirá cuando las restricciones de igualdad sean funciones afines y las restricciones de desigualdad sean convexas. 1.3 Reformulación del problema y resolución Una vez formulado el problema, procedo a su resolución. Teniendo en cuenta lo anterior podemos concluir que ( ) ( )f x h x= . También propongo que ()hx es de la forma: 1 () 2 h x x Rx q x=+ ΤΤ . Primero, tratamos de poner esta expresión en la forma original para identificar términos: Problema de optimización convexa 4 1 ( ) ( ) ( ) ( ) 2 h x x y y R x y y q x y y= − + − + + − + = ΤΤ () () 1 1 1 1( ) ( ) ( ) ( ) ( ) 2 2 2 2 hy y R x y y Ry q y x y R x y q x y x y Ry y R x y − = + + − − + − + − + − = Τ Τ Τ Τ Τ ΤΤ 1 ( ) ( )( ) ( ) ( ) 2 h y y R q x y x y R x y= + + − + − − = Τ Τ Τ 2 1 ( ) , 2R h y q Ry x y x y= + + − + − . Entonces podemos concluir que ()h y Ry q = + . Segundo, dado y , el opt x sin restricciones se obtendrá derivando la función ()hx e igualando a 0. Denotando z x y=− : 2 ( ) 1 ( ( ) , ) 2R hz h y q Ry z z zz  = + + + =  1 0 ( ) opt opt opt q Ry Rz z R q Ry x z y − = + + = → =− + → = + . Ahora introducimos las restricciones en caja, acotando x : x x x −+  x y z x y −+ −   − z z z −+  . 5 5 Algoritmos de primer orden para implementación del Control Predictivo para Seguimiento Reformulamos el problema de optimización convexa en función de z : 2 1 * min ( ) min ( ) , 2R x X z Z f h z h y q Ry z z  = = + + + s.a. z z z −+  . Definiendo la matriz D diagonal, definida positiva e igual al autovalor máximo de R por la matriz identidad. max 0 0 0 0 1 max 0 0 0 0 0 0 max 00 dR RdnR DI         =         == . Y defino:   max 1[1...1] n nR d d d  = =  . El punto óptimo de nuestra función será: 2 1 * argmin ( ) , 2D xX x h y q Ry x y x y  = + + − + − s.a. x x x −+  . Este problema es separable, es decir, cada componente del óptimo puede calcularse por separado resolviendo n problemas escalares, como veremos ahora. En el primer programa adjunto en el Anexo I expongo el programa en MATLAB que implementa la resolución del problema de encontrar el punto óptimo de una función convexa como la que hemos estado tratando. Para ello he usado el método que he expuesto anteriormente, usar la derivada de la función e igualarla a 0 para hallar * i x y comparando con el resultado del comando quadprog. Tal como concluimos antes, para el caso sin restricciones: Problema de optimización convexa 6 1() opt opt opt z R q Ry x z y − =− + → = + . Por lo tanto, en nuestra nueva formulación: *1 ˆ() i i i i x y q Ry d = − + . En cada iteración se define * i x : * * * * * ˆ ˆˆ ˆ i i i i i i i i i i i x x x x x x x x x x x −−   −+ = ++      . 7 2 DUALIDAD EN OPTIMIZACIÓN CONVEXA n la optimización convexa con restricciones hay un concepto que es de gran importancia y que se recoge en la teoría de la dualidad. Esta consiste en conectar una formulación original, también denominada primal, a través de dualidades con otra formulación llamada dual. Conceptualmente se basa en transformar las restricciones del problema primal en variables del problema dual, así es como quedarían conectadas ambas formulaciones. Para hacer esta transformación se recurre al lagrangiano del problema original y usar una serie de propiedades que describiremos ahora. [4],[6],[7],[8]. Nuestro problema canónico incluye la minimización de una función cuadrática que es convexa y suave, con restricciones de igualdad y restricciones en caja. De esta forma: *1 min2 z J z Hz q z  =+ ΤΤ s.a.   : Az b z z z z −+ = =   Siendo la matriz H diagonal y definida positiva. Así pues, este problema puede dualizarse de manera que las restricciones de igualdad en el problema primal se conviertan en una variable del problema dual, relajando las restricciones. Definiendo los teoremas de dualidad débil y dualidad fuerte sentaré las bases de lo que va a ser la dualización del problema canónico, que será un pilar fundamental en su resolución. Este problema dual queda definido a partir de la lagrangiana del primal, basandose en la propiedad de la lagrangiana que sigue, sabiendo que  y  son las variables del dual: Propiedad 2-1: Sea un problema de optimización como nuestro problema canónico, su función lagrangiana asociada ( , , )Lz  , 0, z : ( , , ) ( ) i i Z L z f z        . Se define la función dual: La función dual es una cota inferior de la función objetivo evaluada en el óptimo E 0( , ) min ( , , ). zZ g L z      = Algoritmo FISTA 14 15 4 INTRODUCCIÓN A LA FORMULACIÓN MPC l control predictivo es un conjunto de métodos de control que fue conceptualmente concebido en la industria petroquímica a finales de los 70. Está desarrollado sobre ciertos principios comunes tales como el uso de un modelo de proceso para pronosticar la salida del sistema en momentos futuros de la dinámica y el cálculo de una acción de control óptima basado en la minimización de una función de coste y con la posibilidad de añadir restricciones sobre las variables de ese proceso. Este modelo de la planta se denomina modelo de predicción, y su comportamiento futuro se predice dentro de un intervalo de tiempo llamado horizonte de predición. Una vez aplicada la estrategia para encontrar la acción de control óptima se recalculará ésta desplazando el horizonte de predicción ajustándolo a las condiciones iniciales. Esto se llama horizonte deslizante.[1],[2],[3],[4],[5]. El control predictivo presenta una serie de ventajas respecto de otros métodos de control, entre ellas están. • Permite controlar procesos con comportamientos dinámicos tales como procesos de fase no-mínima o procesos altamente oscilatorios. • Aunque se basa en principios básicos, es una metodología abierta que permite contribuciones y mejoras en el futuro. • Su carácter predictivo hace que compense intrínsecamente los tiempos muertos. • Posibilidad de incorporar restriccionesa a las variables de control. Por otra parte, también tiene sus desventajas, como el alto coste computacional que requiere la complejidad del algoritmo. Aunque lo cierto es que esta desventaja se está paliando con los avances tecnológicos en los últimos años, por lo que podría no considerarse una gran desventaja actualmente. 4.1 Fundamento del MPC El objetivo del MPC consiste en el hallazgo de una trayectoria futura de la variable manipulada u . Para lograr este objetivo, se sigue la siguiente metodología: Se fija un horizonte de predicción N , para el cual se van a predecir las salidas futuras en cada instante. Esto se logra utilizando un modelo de predicción que describe el comportamiento de las variables del proceso a controlar y el cual reside en el controlador. Las predicciones dependen de los valores conocidos hasta el instante i k y de las señales de control futuras. Para calcular las señales de control futuras se utiliza una función de coste que en el caso estándar es cuadrática y penaliza los errores entre la salida predicha y la trayectoria de referencia. En el caso de este trabajo, iré exponiendo gradualmente diferentes funciones de coste y sus resultados, concluyendo en una función más rica y propia de este tipo de estrategias. Este problema, en ausencia de restricciones, puede ser resuelto analíticamente. En caso contrario es necesario un algoritmo iterativo de optimización. Para ambas soluciones, las ecuaciones devuelven un vector que E Introducción a la formulación MPC 16 16 contiene las acciones futuras, cuya dimension depende del horizonte de control. Sin embargo, solamente el primer elemento de este vector tiene interés para nuestra estrategia y se calculará uno nuevo en cada instante de la dinámica. Esto se refleja en la siguiente figura. Ilustración 4-1 4 4.2 Formulaciones del problema MPC Vamos a considerar 3 formulaciones diferentes para plantear la función de coste que queremos minimizar para nuestro problema MPC, desde una más básica hasta otra más elaborada. El procedimiento consistirá en transformar la formulación inicial en un problema cuadrático QP que dualizaremos para poder resolver. Para la resolución dual utilizaré el algoritmo FISTA en las tres formulaciones, comparándolo con el resultado del commando quadprog . En el caso de la tercera formulación, utilizaré también el algoritmo ADMM (Alternating Direction Method of Multipliers), que será objeto del siguiente capítulo. 4 Esquema de explicación de la dinámica típica de un sistema controlado con MPC. 17 Algoritmos de primer orden para implementación del Control Predictivo para Seguimiento 4.2.1 Formulación 1 Vamos a partir de una formulación simplificada de un problema MPC, en la que tendremos una serie de restricciones y una simplificación inicial que dice que nuestro punto de equilibrio estará en cero. El objetivo será no circunscribirnos a la restricción de esta formulación y poder hacer algo con condiciones más libres, lo que veremos en la siguiente formulación (2). 0 0 ,..., 0 ,..., min ( ) ( ) ( ) ( ) N N N i e i e i e i e uu i xx J x x Q x x u u R u u = = − − + − − = ΤΤ 0 0 11 ,..., ,..., 1 min 2 N N uu xx z H z q z=+ ΤΤ . Sabiendo que 0, 0QR que son diagonales y denominando 0 0 1 1 2 2 1 1 , , , , , , , , , , N N N N z x u x u x u x u x u −− =  1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 20 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Q R Q HR Q R      =      1[ 2 2 2 2 2 2 ] e e e e e e q x Q u R x Q u R x Q u R= − − − − − − Τ Τ Τ Τ Τ Τ Τ . El problema está sujeto a las siguientes restricciones de igualdad y en caja. 0 xx= , condición inicial que es dato. Introducción a la formulación MPC 18 18 1, 0,..., 1 i i i x Ax Bu i N += + = − Ne Ne xx uu = = , 0,..., , 0,..., . i i x x x i N u u u i N   =   = Que podemos reescribir definiendo las matrices 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 A B I A B I AB A D A B I −   −     =  −      0 0 0 0 0 e e b x u       =       . Como simplificación inicial, tenemos que 0, 0 ee xu== . Este es el punto de equilibrio y ha de cumplir que 19 Algoritmos de primer orden para implementación del Control Predictivo para Seguimiento e e e Ax Bu x+= . Reescribiendo el problema, tenemos: 0 0 11 ,..., ,..., 1 min 2 N N uu xx z H z q z+ ΤΤ s.a. Dz b= z z z . Sabiendo que , xx uu xx z u z u xx uu                     ==                         . Tal como vimos en capítulos anteriores, vamos a aplicar el algoritmo FISTA para resolver el problema de optimización del dual a este, pero de la forma que nos explica el método MPC descrito anteriormente. Para ello voy a comenzar por dualizar el problema y formularlo de manera que sea sencillo de resolver en función de las variables duales. El funcional dual sería 11 0 1 ( ) min ( ) 2 N ziz H z q z Dz b     = = + + − Τ Τ Τ 2 () 1( ) 1( ) ( ) ( ) 0 argmin 2 ni ii i i zZi z h q c z = = + +  . Ya vimos que este problema es separable, por lo que hallamos el óptimo resolviendo un problema escalar por Introducción a la formulación MPC 20 20 cada componente. Teniendo en cuenta que h→ diagonal de H y que cD  =Τ . Por lo tanto, para calcular el óptimo del problema dual: ( ) ( ) ( ) ( ) 1( ) 1( ) () ( ) ( ) () 1( ) () ( ) ( ) 1( ) () ii ii ii i ii ii i ii i cc zz hh c y z z h c zz h  − −      = −     −   Como ya se explica en el capítulo 2. Este proceso se hará para todas las formulaciones que presentaré aquí, para así poder proceder a su resolución con el agortimo FISTA, que trabaja con la formulación dualizada del problema de optimización. Este algoritmo va a ser utilizado de cara a resolver el problema basado en MPC, que nos devolverá las acciones de control que minimizan la función de coste en función de predicciones futuras del sistema. Para ello seguimos el siguiente procedimiento: 1) Elegir el horizonte de predicción y el número de componentes del vector de control: N y n respectivamente. 2) Definir el punto de equilibrio: 0, 0 ee xu== en este primer caso. 3) Inicializar vectores de salida ( , )x n N y de actuación ( , )u n N . Así como sus respectivos límites para las restricciones en caja: , , ,x x u u . 4) Definir la condición inicial, que será dato para el problema, por ejemplo, elijo el valor medio entre los límites inferior y superior de la salida. 02 xx x+ = . 5) Construimos las matrices 1 ,,,, ,Q R A B H D . 6) Elegimos el parámetro de tolerancia  , 5 10− por ejemplo. 7) REPETIR mientras la norma de la diferencia entre la primera componente del vector que contiene la salida del sistema y la acción de control ( z ), y el vector de referencias l sea mayor que  . 8) FOR 1: 1iN=− 21 Algoritmos de primer orden para implementación del Control Predictivo para Seguimiento 9) Ejecutar algoritmo FISTA con condición que cumple las restricciones de igualdad, ( (:, ) )norm Dz i b  − . 4.2.2 Formulación 2 Abordamos ahora una formulación más propia de MPC, aún con sus limitaciones. Siendo , ss xu un punto de equilibrio que cumple estrictamente las restricciones. Esta formulación garantiza que el problema va a ser factible, asumiendo que existe un punto de equilibrio que cumple las restricciones,cosa que tenemos en , ss xu . Usaremos una estrategia de control que se fundamenta en el uso de un objetivo artificial de control s x . Dicha estrategia se basa en ponderar tanto una distancia del estado k-ésimo y de la actuación k-ésima al estado y actuación objetivo artificial con respecto a la referencia. Si esta última ponderación es elevada, la salida final será cercana a la referencia, mientras que si es baja la salida no tiene por qué acercarse a la referencia requerida. [4] 5 . 01 0 1 ,..., 0 ,..., , min ( ) ( ) ( ) ( ) N N N i s i s i s i s uui xx xu ss J x x Q x x u u R u u − − = = − − + − − + ΤΤ 01 ( ) ( ) ( ) ( ) ( ) ( ) e s e s s e s e N e N e x x T x x u u S u u x x x x T x x  + − − + − − + − + − − = Τ Τ Τ 01 0 22 ,..., ,..., , 1 min 2 N N uu xx xu ss z H z q z − =+ ΤΤ . Sabiendo que: • x (condición inicial) es dato. • ( , ee xu ), punto de equilibrio que también viene dado. •  es un escalar 0 . 5 Esta estrategia fue desarrollada por el Departamento de Sistemas y Automática de la Escuela Técnica Superior de Ingeniería de la Universidad de Sevilla. Esta publicada en el artículo ‘’MPC for tracking piecewise constant references for constrained linear systems’’, escrito por D. Limón, I. Alvarado, T. Álamo y E. F. Camacho. Introducción a la formulación MPC 22 22 La función está sujeta a una serie de restricciones 1, 0,..., 1 i i i x Ax Bu i N += + = − s s s x Ax Bu=+ Ne xx= , 0,..., , 0,..., 1 i i x x x i N u u u i N   =   = − s s x x x u u u   +   − +   − La formulación número 2 tan sólo tiene una intención de acercamiento a un función más propia de MPC, el cálculo del óptimo en cada iteración se realizará tanto con el algoritmo ISTA como con el FISTA. 4.2.3 Formulación 3 Para este bloque se va a trabajar con una formulación más compleja y adaptada a lo que es una expresión típica basada en MPC. Tras definirla, se calculará su expresión dual para proceder a su resolución en el siguiente capítulo. El capítulo 5 estará basado en el uso de método ADMM a través de la expresión dual de esta formulación. • Definimos 2 vectores que incluirán las variables de decisión que aparecen en la formulación 2 ( v ) y las variables que están restringidas ( y ). Estos vectores pueden colocarse según 2 disposiciones diferentes: separada e intercalada. Aquí usaremos la ordenación intercalada de manera que H nos quede como una matriz en banda, esto lo veremos luego. 0 0 1 1 2 2 1 1 , , , , , , , , , , , N N N s s v x u x u x u x u x x u −− =  . Teniendo que 1, 0,..., 1 s s s i i i Ns x Ax Bu x Ax Bu i N xx + =+ = + = − = 23 Algoritmos de primer orden para implementación del Control Predictivo para Seguimiento Estas restricciones pueden reescribirse de modo que queden de la forma 0Dv= . Donde D es una matriz de la forma 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0000000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ( ) A B I A B I A B I A B I II A I B D −   −   −    −  −   −  = . El segundo vector es 0 0 1 1 2 2 1 1 , , , , , , , , , , , N N N s s y x u x u x u x u x x u −− =  . Cuyas restricciones son las siguientes , 0,..., , 0,..., 1 i i s s x x x i N u u u i N x x x u u u     =   = − +   − +   − Al igual que antes, reescribimos las restricciones de manera que quedan de la forma. Introducción a la formulación MPC 30 30 Ahora debemos analizar una casuística, ya que la incógnita a puede ser igual o distinto de cero. Sólo será distinto de cero para 0 x . Caso 1: 0a= 0 2* ˆ min ( ) r( ) 0 2z z rz d gz z d g =  − + ⎯⎯⎯→ − + = * ... [ , ] ˆ ... ... gg d d z z rr g g z d z d z r rg z d z r − −     = − = −    −   . Dicho de otra forma * ˆmax , min , g z z d z r   =−     . Caso 2: 0a Se puede resolver el problema 2 veces en función del signo de * ˆ zs− . 1. Asumimos que * ˆ0zs− 0 2* ˆ min ( ) r( ) 0 2z z rz d az as gz z d a g =  − + − + ⎯⎯⎯→ − + + = . Despejando y denotando g g a=+ 31 Algoritmos de primer orden para implementación del Control Predictivo para Seguimiento * ... [ , ] ˆ ... max , min , ... gg d d z z rr gg g z d z d z z d z rr rg z d z r  − −       = − = −  = −          −   . 1. Asumimos que * ˆ0zs− 0 2* ˆ min ( ) r( ) 0 2z z rz d az as gz z d a g =  − − + + ⎯⎯⎯→ − − + = . Despejando y denotando g a g=−     * ... [ , ] ˆ ... max , min , ... gg d d z z rr g g z d z d z z d g z r rg z d z r  + +     = + = +  = +    +   . Para las siguientes resoluciones hará falta poner el funcional ()Jy  de la forma QP. Las matrices y H y y q quedarían así 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 Q R Q R Q R T T S I I I I yI I I I I H                          = Introducción a la formulación MPC 32 32 0 0 1 1 11 2 2 2 2 2 2 2 ( ) 2 ( ) 2 ( )[] N N N S S e Q e R e Q e R e R e R e T e T e S x u x u x x u xu x I u I x I u I u I u I x I x I u I y q           −− − + − − − − − − − −=Τ Τ Τ Τ Τ Τ Τ Τ ΤΤ Siendo  un vector cuyas componentes pueden ser 1 o -1, en función del valor de 0,1 1 xx− . 1 2nx j Q e e R e e T e e S e e j c Nx x Nu u x x u u x      = =− − + + −  Τ Τ Τ Τ . Una vez determinada la componente óptima en ambas suposiciones, vemos si * ˆ zs− es mayor o menor que cero en los 2 casos. Cuando tengamos el óptimo que cumple su hipótesis, descartamos la otra opción. Hecho esto con todas las componentes, ya tendremos determinada la y  . Como voy a resolver el problema con el algoritmo FISTA, el próximo paso será calcular la   utilizando la siguiente propiedad: Propiedad 4-1. Denotanto 1 argmin ( ) 2 zZ z z Hz q z Az b   = + + − Τ Τ Τ . Entonces, si H es definida positiva, según lo expresado en [7],[8]. 1 1 ( ) ( ) ( ) 2 Az b AH A          − +  + − −   Τ Τ Τ . Al ser nuestro problema de la forma , 11 min ( ) ( ) 22 v y v y v y vy J J J v y v H v y H y q v q y v y   = + + − = + + + + − Τ Τ Τ Τ Τ Τ Definiendo las matrices []A I I=− 0b= . 33 Algoritmos de primer orden para implementación del Control Predictivo para Seguimiento De esta forma [ ] 0 0 v Az b I I v y y  = → − = → − =   . Tal como ya hemos comentado, dado un  podemos calcular tanto v  como y  . Ahora definimos ( ) ( ) ( ) ( ) vy J v J y v y        = + + − Τ . Utilizando la propiedad anterior y las matrices de la dinámica 1 1 ( ) ( ) ( ) 2 v y AH A          − +  + − −   Τ Τ Τ . Con 0 0 v y H HH  =  , tenemos que 1 10 0 v vy y HI AH A I I H H HI − −  = − = +     −    Τ . Reformulando 1 ( ) ( ) ( ) ( ) 2vy v y H H          +  + − −  +  ΤΤ . Derivando respecto de   podemos concluir que el valor que maximiza la cota inferior de ()    + es 1 ( ( ) ( ) ( ) ) 0 2vy v y H H        + − −  +  = →  ΤΤ ( ) 0 vy v y H H   → − − −  = . Introducción a la formulación MPC 34 34 Despejando obtenemos Planteo aquí el algoritmo a implementar para resolver esta última formulación con FISTA -------------------------------------------------------------------------------------------------------------------------------------- 1) ELEGIMOS PARÁMETRO DE TOLERANCIA  . 0 0 1 0 2) 0, 1, 1tk    = = = = = 3) REPETIR 21 1 4) (1 1 4 ) 2k k tt − = + + 11 1 5) ( ) k k k k k k t t     −− − = + − 1 6) argmin 2 k v v k v v H v q v v  = + + Τ Τ Τ 7) 1 argmin 2 k y y k y y H y q y y  = + − Τ Τ Τ 8) 1 1( ) ( ) k k v y k k H H v y  − += + + − 9) SI 2 kk vy  − ENTONCES EXIT 10) 1kk=+ -------------------------------------------------------------------------------------------------------- Cuya implementación se encuentra en el anexo de códigos al final del documento. Con esta formulación MPC que se expone, se garantiza que el problema es siempre factible, ya que si la condición inicial no puede llevar el sistema al punto de equilibrio sin violar las restricciones, se va a calcular la condición inicial más próxima que sí es capaz de hacerlo. Para trabajos futuros, se podría aprovechar la estructura en semi-banda de problemas intermedios de optimización para poder implementar algoritmos más eficientes, cuyo coste sea menor. *1 ( ) ( ). vy H H v y   −  = + − 35 5 INTRODUCCIÓN AL MÉTODO ADMM l método de multiplicadores de dirección alterna usa actualizaciones parciales para las variables duales. Este método se aplica a menudo para resolver problemas como el nuestro, es decir, del tipo , min ( ) ( ) , s.a. xy f x g y x y+= . Este problema puede atacarse con métodos de optimización restringida, cuya función objetivo es separable en x e y . [2],[3]. La actualización dual requiere resolver una función de proximidad en x e y al mismo tiempo. La técnica ADMM permite resolver este problema aproximadamente resolviendo primero x con y fijo, y luego resolviendo y con x fijo. En lugar de iterar hasta la convergencia, el algoritmo procede a actualizar la variable dual y luego repite el proceso. Esto no es equivalente a la minimización exacta, pero se puede demostrar que este método converge a la respuesta correcta. [2]. Si particularizo a nuestro caso, tendríamos , 0 min + s.a. vy vy Dv J J J y y y vy =   =    =  . Sabiendo que v J y y J se corresponden con las expresiones del capítulo anterior. Ahora, dualizando el problema nos quedará la expresión que pasaremos a la formulación ADMM. , 0 min ( ) s.a. vy vy Dv J J v y yyy  =   + + −    Τ . El problema ADMM será de la forma 2 2 , ( ) min ( ) 2 vy vy J J v y v y      = + + − + − Τ . E Introducción al método ADMM 36 36 El algoritmo se basa en dejar fija una de las variables y calcular la otra ( fijo y new v→ ), para después volver a calcular la anterior ( new new vy→ ). A continuación, actualizamos  . Tendrá la siguiente estructura -------------------------------------------------------------------------------------------------------------------------------------- 1) ELEGIMOS 0   ( parámetro de tolerancia 5 10− ). 2) ELEGIMOS y 0 0yY   3) 0 0, 0k  == 4) 2 12 argmin 2k k v k v v J v v y   += + + − Τ 5) 2 1 12 argmin 2k k y k y y J y v y   + += − + − Τ 6) 1 1 1 () k k k k vy    + + + = + − 7) SI   1 1 1 22 max , k k k k v y y y  + + + − −  → EXIT 8) 1 y (4)k k GOTO=+ -------------------------------------------------------------------------------------------------------------------------------------- A continuación, desarrollo el procedimiento para resolver los puntos 4 y 5 del algoritmo.[2],[3]. 4) 1 2 v v v J v H v q v c= + + ΤΤ 2 2 () 2k k v k v J v v y    = + + − = Τ 1 2 2 2 v v k k k k v H v q v c v v v y v y y   = + + + + − + → Τ Τ Τ Τ Τ Τ 37 Algoritmos de primer orden para implementación del Control Predictivo para Seguimiento 0( ) ( ) 0 0 vv v k k H I v q y Dv    =   ⎯⎯⎯→ + + + − = =   = 1 0 0 k k v k vyq v HD D   +   − + −   ==        Τ . Desarrollando el sistema de ecuaciones 1 () v k k k v H I v D y q     + + + =− + − Τ 10 k Dv += . Despejando 1 11 ( ) ( ) 0 k v k k v k v H I y q D Dv     − ++ = + − + − − → = → Τ 1 ( ) ( ) v k k v D H I y q D     − → + − + − − = Τ 11 ( ) ( ) ( ) 0 v v k k v D H I D D H I y q      −− =− + + + − + − = → Τ 11 ( ) ( ) ( ) v v k k v D H I D D H I y q      −− → + = + − + − Τ . Simplificando 1 1 1 ( ( ) ) ( ) ( ) v v k k v D H I D D H I y q      − − − = + + − + − Τ 1 1 1 1 1( ) ( ( ( ) ) ( ) ( )) k v k k v v v k k v v H I y q D D H I D D H I y q        − − − − += + − + − − + + − + − ΤΤ Introducción al método ADMM 38 38 Ahora, vamos a ver como desarrollar el punto 5. 5) 1 2 y y y J y H y q y c= + + ΤΤ 2 12 () 2k k y k y J y v y   + = − + − = Τ 1 1 1 1 2 2 2 y y k k k k y H y q y c y y y y v v v   + + + = + + − + − + Τ Τ Τ Τ Τ Τ Denotando 1kk fv  + =− 2 1argmin 2 ify yii kyY y y y f y    + =− ΤΤ La componente i-ésima de 1k y+ será 1 [ , ] ii i k i ff yy f y y y f yy    +     =      En MATLAB 1max ,min , i k f y y y  +    =     39 Algoritmos de primer orden para implementación del Control Predictivo para Seguimiento Para mostrar resultados del algoritmo anterior, expongo una serie de representaciones gráficas que muestran el comportamiento del algoritmo en función de una serie de parámetros. Para caracterizar el problema que se ha solucionado para obtenerlas, defino aquí los siguientes parámetros: Horizonte de predicción: N=50 Dimensión del problema: n=25 Para la construcción de las matrices , , y Q R T S : () 2 max( ) min( ) 2 Q Q Q Q randn n Q QQ I QI    = =+ = = Τ () 2 max( ) min( ) 2 R R R R randn n R RR I RI    = =+ = = Τ 2 max( ) min( ) 2 T T T T cI T TT I TI    = =+ = = Τ 2 max( ) min( ) 2 S S S S cI S SS I SI    = =+ = = Τ 100  = 0.1  = COND. INICIAL TIEMPO DE EJECUCIÓN (sg) -0.333 92.4552 0.3273 99.2307 2.5777 113.2307 0.7436 113.2016 -1.3344 94.2417 -4.4671 99.7304 COND. INICIAL TIEMPO DE EJECUCIÓN (sg) -1.0265 15.7848 2.3128 16.6988 4.8788 20.9689 -7.1629 16.2011 14.5903 16.0353 -21.4026 15.8792 ANEXO: Códigos de MATLAB 46 46 uplus=7*ones(n,1); umin=-8*ones(n,1); u=zeros(n,N); %========================================================================= %Construcción de Q: Q=randn(n); Q=Q*Q'+2*eye(n); Q=max(eig(Q))*eye(n); %Construcción de R: R=randn(n); R=R*R'+2*eye(n); R=max(eig(R))*eye(n); %========================================================================= A=randn(n); B=randn(n); C=randn(n); K=randn(n); %Punto de equilibrio(xe,ue): ye=randn(n,1); s=[(eye(n)-A) -B; C K]\[zeros(n,1);ye]; xe=s(1:n); ue=s(n+1:2*n); %Matrices H y D: I=zeros(n); H(:,:)=[Q I; I R]; q=zeros(n*2,1); D(:,:)=[A B]; sum_qp=0; sum=0; %Parámetros de tolerancia y contador de iteraciones: epsilon1=10^(-5); epsilon2=10^(-5); k=0; z(:,:)=[x;u]; l(:,:)=[xe;ue]; %========================================================================= %Bucle de control: while norm(z(:,1)-l)>epsilon2 for i=1:N-1 zplus=[xplus-xe ; uplus-ue]; zmin=[xmin-xe ; umin-ue]; b=x(:,i+1)-xe; zopt_qp(:,i)=quadprog(H,q,[],[],D,b,zmin,zplus); f=zeros(n,1); fpost=f; lambda=zeros(n,1); tant=0; t=tant; 47 Algoritmos de primer orden para implementación del Control Predictivo para Seguimiento z(:,i)=[x(:,i)-xe;u(:,i)-ue]; lambda_1=0; g=0; %Algoritmo ISTA: while norm(D*z(:,i)-b)>epsilon1 z(:,i)=ARGMIN2(D,lambda,H,zplus,zmin,n); lambda_1=lambda_1+inv(D*inv(H)*D')*(D*z(:,i)-b); g=g+1; end %Algoritmo FISTA: w=0; while norm(D*z(:,i)-b)>epsilon1 tant=t; fant=f; f=fpost; t=0.5*(1+sqrt(1+4*(tant^2))); lambda=f+((tant-1)/t)*(f-fant); z(:,i)=ARGMIN2(D,lambda,H,zplus,zmin,n); fpost=lambda+inv(D*inv(H)*D')*(D*z(:,i)-b); w=w+1; end zopt(:,i)=z(:,i); xopt_qp(:,i)=zopt_qp(1:n,i)+xe; uopt_qp(:,i)=zopt_qp(n+1:2*n,i)+ue; xopt(:,i)=zopt(1:n,i)+xe; x(:,i)=xopt(:,i); uopt(:,i)=zopt(n+1:n*2,i)+ue; u(:,i)=uopt(:,i); sum_qp=sum_qp+zopt_qp(:,i)'*H*zopt_qp(:,i); sum=sum+zopt(:,i)'*H*zopt(:,i); end x(:,1)=x(:,2); k=k+1; end %Comparación algoritmo FISTA implementado con quadprog: S(:,:)=[uopt_qp(:,1) uopt(:,1)]; F(:,:)=[xopt_qp(:,1) xopt(:,1)]; SUM(:,:)=[sum_qp sum]; ARGMIN2 function [ z ] = ARGMIN2( D,lambda,H,zplus,zmin,n ) c=D'*lambda; h=diag(H); for i=1:(n*2) z(i)=max(zmin(i),min(-c(i)/h(i),zplus(i))); end z=z'; end ANEXO: Códigos de MATLAB 48 48 Formulación 3 Main program N=50; n=25; tau=10000000; epsilon=10^(-3); %A,B,C y K: A=randn(n); B=randn(n); C=randn(n); K=randn(n); %Punto de equilibrio(xe,ue): ye=randn(n,1); s=[(eye(n)-A) -B; C K]\[zeros(n,1);ye]; xe=s(1:n); ue=s(n+1:2*n); %Valores para restricciones en caja: xmax=xe+epsilon+3*ones(n,1); xmin=xe-epsilon-3*ones(n,1); x_ini=xmax-1; %x dato inicial. umax=ue+epsilon+3*ones(n,1); umin=ue-epsilon-3*ones(n,1); xsmin=xmin+epsilon; xsmax=xmax-epsilon; usmin=umin+epsilon; usmax=umax-epsilon; 49 Algoritmos de primer orden para implementación del Control Predictivo para Seguimiento %Restricciones en caja de y: ymin=zeros((2*N+3)*n,1); for i=1:2*n:(2*N+2)*n+1 if i<(2*N+1)*n-n+1 ymin(i:i+n-1)=xmin; ymin(i+n:i+2*n-1)=umin; end if i>(2*N+1)*n-n && i<(2*N+1)*n+1 ymin(i:i+n-1)=xmin; end end for i=(2*N+1)*n+1:n:(2*N+3)*n if i>(2*N+1)*n && i<(2*N+2)*n+1 ymin(i:i+n-1)=xsmin; end if i>(2*N+2)*n ymin(i:i+n-1)=usmin; end end ymax=zeros((2*N+3)*n,1); for i=1:2*n:(2*N+2)*n+1 if i<(2*N+1)*n-n+1 ymax(i:i+n-1)=xmax; ymax(i+n:i+2*n-1)=umax; end if i>(2*N+1)*n-n && i<(2*N+1)*n+1 ymax(i:i+n-1)=xmax; end end for i=(2*N+1)*n+1:n:(2*N+3)*n if i>(2*N+1)*n && i<(2*N+2)*n+1 ymax(i:i+n-1)=xsmax; end if i>(2*N+2)*n ymax(i:i+n-1)=usmax; ANEXO: Códigos de MATLAB 50 50 end end %Construcción de Q y de tauQ_eye: Q=randn(n); Q=Q*Q'+2*eye(n); Q=max(eig(Q))*eye(n); tauQ=min(eig(Q))/2; tauQ_eye=tauQ*eye(n); %Construcción de R y de tauR_eye: R=randn(n); R=R*R'+2*eye(n); R=max(eig(R))*eye(n); tauR=min(eig(R))/2; tauR_eye=tauR*eye(n); %Construcción de T y de tauT_eye: T=350*eye(n); %T=cte*I, la constante la elegí hasta que Hv fuera positiva. T=T*T'+2*eye(n); T=max(eig(T))*eye(n); tauT=min(eig(T))/2; tauT_eye=tauT*eye(n); %Construcción de S y de tauS_eye: S=350*eye(n); %S=cte*I, la hice igual que T. S=S*S'+2*eye(n); S=max(eig(S))*eye(n); tauS=min(eig(S))/2; tauS_eye=tauS*eye(n); %===================================================================== %Construcción de las matrices H_v, q_v y D: %H_v: H_v=zeros((2*N+3)*n); %H_v: for i=1:2*n:(2*N+1)*n+1 51 Algoritmos de primer orden para implementación del Control Predictivo para Seguimiento if i<(2*N)*n-n+1 H_v(i:i+n-1,i:i+n-1)=Q-tauQ_eye; H_v((2*N+1)*n+1:(2*N+2)*n,i:i+n-1)=-Q; H_v(i+n:i+2*n-1,i+n:i+2*n-1)=R-tauR_eye; H_v((2*N+2)*n+1:(2*N+3)*n,i+n:i+2*n-1)=-R; end if i>(2*N)*n-n && i<(2*N+1)*n+1 H_v(i:i+n-1,i:i+n-1)=T-tauT_eye; end end for i=(2*N+1)*n+1:n:(2*N+3)*n if i>(2*N+1)*n && i<(2*N+2)*n+1 for j=1:2*n:(2*N)*n-n H_v(j:j+n-1,(2*N+1)*n+1:(2*N+2)*n)=-Q; end H_v((2*N+1)*n+1:(2*N+2)*n,i:i+n-1)=N*Q+(T-tauT_eye); end if i>(2*N+2)*n for j=1:2*n:(2*N)*n-n H_v(j+n:j+2*n-1,(2*N+2)*n+1:(2*N+3)*n)=-R; end H_v((2*N+2)*n+1:(2*N+3)*n,i:i+n-1)=N*R+(S-tauS_eye); end end H_v=2*H_v; %q_v: q_v=zeros(1,(2*N+3)*n); for i=1:2*n:(2*N+1)*n+1 if i<(2*N)*n-n+1 q_v(i:i+n-1)=2*xe'*tauQ_eye; q_v(i+n:i+2*n-1)=2*ue'*tauR_eye; end if i>(2*N)*n-n && i<(2*N+1)*n+1 q_v(i:i+n-1)=-2*xe'*(T-tauT_eye); end end for i=(2*N+1)*n+1:n:(2*N+3)*n ANEXO: Códigos de MATLAB 52 52 if i>(2*N+1)*n && i<(2*N+2)*n+1 q_v(i:i+n-1)=-2*xe'*(T-tauT_eye); end if i>(2*N+2)*n q_v(i:i+n-1)=-2*ue'*(S-tauS_eye); end end q_v=q_v'; %D: gg=-eye(n); j=0; D=zeros((N+2)*n,(2*N+3)*n); for i=1:n:(N+1)*n+1 if i<N*n+1 D(i:i+n-1,i+j*n:i+j*n+n-1)=A; D(i:i+n-1,i+j*n+n:i+j*n+2*n-1)=B; D(i:i+n-1,i+j*n+2*n:i+j*n+3*n-1)=gg; j=j+1; end if i>N*n && i<(N+1)*n+1 D(i:i+n-1,i+N*n:i+(N+1)*n-1)=eye(n); D(i:i+n-1,i+(N+1)*n:i+(N+2)*n-1)=gg; end if i>(N+1)*n D(i:i+n-1,(2*N+1)*n+1:(2*N+2)*n)=A-eye(n); D(i:i+n-1,(2*N+2)*n+1:(2*N+3)*n)=B; end end %Construcción de H_y, q_y: %H_y: H_y=zeros((2*N+3)*n); for i=1:2*n:(2*N+1)*n+1 if i<(2*N+1)*n-n+1 H_y(i:i+n-1,i:i+n-1)=tauQ_eye; 53 Algoritmos de primer orden para implementación del Control Predictivo para Seguimiento H_y(i+n:i+2*n-1,i+n:i+2*n-1)=tauR_eye; end if i>(2*N+1)*n-n && i<(2*N+1)*n+1 H_y(i:i+n-1,i:i+n-1)=tauT_eye; end end for i=(2*N+1)*n+1:n:(2*N+3)*n if i>(2*N+1)*n && i<(2*N+2)*n+1 H_y((2*N+1)*n+1:(2*N+2)*n,i:i+n-1)=tauT_eye; end if i>(2*N+2)*n H_y((2*N+2)*n+1:(2*N+3)*n,i:i+n-1)=tauS_eye; end end H_y=2*H_y; %q_y: q_y=zeros(1,(2*N+3)*n); for i=1:2*n:(2*N+1)*n+1 if i<(2*N+1)*n-n+1 if i==1 q_y(i:i+n-1)=-2*xe'*tauQ_eye+tau; else q_y(i:i+n-1)=-2*xe'*tauQ_eye; end q_y(i+n:i+2*n-1)=-2*ue'*tauR_eye; end if i>(2*N+1)*n-n && i<(2*N+1)*n+1 q_y(i:i+n-1)=-2*xe'*tauT_eye; end end for i=(2*N+1)*n+1:n:(2*N+3)*n if i>(2*N+1)*n && i<(2*N+2)*n+1 q_y(i:i+n-1)=-2*xe'*tauT_eye; end if i>(2*N+2)*n q_y(i:i+n-1)=-2*ue'*tauS_eye; ANEXO: Códigos de MATLAB 54 54 end end q_y=q_y'; FISTA % FISTA: %Inicialización de variables algoritmo: lambda=0*ones((2*N+3)*n,1); phi2=0*ones((2*N+3)*n,1); phi3=0*ones((2*N+3)*n,1); t1=0; t2=0; k=1; v_opt=7*ones((2*N+3)*n,1); y_opt=5*ones((2*N+3)*n,1); b=1; while (norm(v_opt-y_opt)>epsilon) && (b==1) t1=t2; phi1=phi2; phi2=phi3; t2=0.5*(1+sqrt(1+4*t1^2)); lambda=phi2+((t1-1)/t2)*(phi2-phi1); v_opt=calculo_v(N,n,lambda,q_v,H_v,D); [y_opt,q_y]=calculo_y(N,n,q_y,lambda,xe,ue,x_ini,ymin,ymax,tauQ,tauR,tauT,tauS,tau); phi3=lambda+(inv(inv(H_v)+inv(H_y))*(v_opt-y_opt)); [Jv,b]=varphy(H_v,q_v,H_y,q_y,lambda,v_opt,y_opt); k=k+1; end L=[v_opt y_opt]; Calculo_v function [ v_opt ] = calculo_v( N,n,lambda,q_v,H_v,D ) %Matrices de ceros para el sistema de ecuaciones: m=zeros((N+2)*n); w=zeros((N+2)*n,1); 55 Algoritmos de primer orden para implementación del Control Predictivo para Seguimiento %Matriz de coeficientes: F=[H_v D';D m]; %Matriz solución: e=[-q_v-lambda;w]; %Resultado: o=F\e; v_opt=o(1:(2*N+3)*n); end Calculo_y function [ y_opt,q_y ] = calculo_y( N,n,q_y,lambda,xe,ue,x,ymin,ymax,tauQ,tauR,tauT,tauS,tau ) %x -> dato inicial. %ymin,ymax -> las restricciones en caja, están construidas en el main. for i=1:(2*N+3)*n %bucle que recorre el vector y. %El siguiente if es para ir haciendo que j vaya de 1 a n cada vez que %llegamos a una nueva variable de decisión: if (i<n+1) %Si el índice es menor que n+1, entonces estoy en x(0). j=i; else j=mod(i,n); %Si el índice es superior, miro el resto de dividir i entre n %y la j será igual a ese resto, excepto si es cero. if (j==0) %Si el resto es 0, j=n. j=n; end end if (mod(i,2*n)<=n) && (mod(i,2*n)~=0) %Para ver la variable de decisión que estoy optimizando es x(1,2,...,N-1). %Los parámetros son los que acompañan a los términos de x en el %funcional. d=xe(j); r=2*tauQ; if (i>2*N*n) && (i<(2*N+1)*n+1) %Si es x(N), r es otra. r=2*tauT; end else %Si estoy optimizando una de las u(1,2,...,N-1). d=ue(j); r=2*tauR; end %Lo que sigue es para las 2 últimas variables, xs e us, que por el orden en %el que están puestas en y, no corresponden al lugar que se le asignaría de %la forma anterior. if (i>(2*N+1)*n) && (i<(2*N+2)*n+1)