Proyecto Fin de Carrera Implementación de algoritmos MPC con restricciones en mbed NXP LPC1768

Capítulo 4

Solver: Algoritmo de Lemke

Dpto. Ing. de Sistemas y Automática Universidad de Sevilla

Página 45

Ramón Jesús Cuesta Pérez Curso 2011/2012

Proyecto Fin de Carrera Implementación de algoritmos MPC con restricciones en mbed NXP LPC1768

4.1 Formulación teórica En este capítulo vamos a mostrar y argumentar la implementación del algoritmo de Lemke como solver para nuestro algoritmo de control predictivo. Para ello vamos a realizar un enfoque teórico previo en el que comenzaremos introduciendo el concepto de Problema Lineal Complementario, plementario, el cual nos permitirá aplicar el algoritmo de Lemke para resolver el problema blema QP que nos plantea el algoritmo de control predictivo.

Problema Lineal Complementario Sean q y M un vector de dimensión m dado y una matriz NxN dada respectivamente. El Problema Lineal Complementario (LCP) consiste en encontrar dos vectores s y z tales que: s − Mz = q

s, z ≥ 0

,

,

s, z

=0

Se dice que una na solución (s, ( z) al sistema anterior es una solución básica factible complementaria si para cada par de las variables complementarias (si, zi) una de ellas es básica para i = 1,…,m m , donde si y zi son las componentes de los vectores s y z, respectivamente.

Si q es no negativo, se puede hallar una solución básica factible complementaria haciendo s = q y z = 0.. Si este no es el caso podemos utilizar el algoritmo de Lemke, en el cual se introduce una variable artificial z0, lo cual nos lleva a: s − Mz − 1z0 = q

,

s, z, z0 ≥ 0

,

s, z

=0

Obtenemos una solución inicial para el sistema anterior haciendo z0 = max(−q max( i), z = 0 y s = q + 1z0. Si mediante una secuencia de pivoteo compatible con el sistema anterior llevamos la variable artificial z0 a cero habremos hallado una solución al Problema Lineal Complementario.

Dpto. Ing. de Sistemas y Automática Universidad de Sevilla

Página 46

Ramón Jesús Cuesta Pérez Curso 2011/2012

Proyecto Fin de Carrera Implementación de algoritmos MPC con restricciones en mbed NXP LPC1768

Una manera eficiente de encontrar una secuencia de pivoteo que converja a una

solución en un número finito de iteraciones es usar el algoritmo de Lemke partiendo de la siguiente tabla: s Imxm

s

z −M

z0 −1

q

Aplicación del LPC al problema QP del algoritmo de control predictivo

Una vez asimilado el concepto del LPC volvamos al problema QP que nos ocupa, que

como vimos en el capítulo 2 no es otro que:

u** = arg min u

1 T u Hu + bu + f0 2

sujeto a Ru ≤ c

Donde H = 2·(GTG + λI),

b = 2·(f − w)TG,

f0 = (f − w)T(f − w)

  1u   − 1u    1U − 1u (t − 1)   c=  − 1U + 1u (t − 1)    1y − f    − 1 y + f 

 I NxN  − I   NxN   T  R=    −T   G     − G 

Podemos reducir el número de restricciones (las segundas N restricciones de R) haciendo el cambio de variable u = 1 u + x, con lo que el problema QP quedaría:

x* = arg min x

1 T x Hx + ax + f1 2

sujeto a Rx ≤ c

x≥0

Dpto. Ing. de Sistemas y Automática Universidad de Sevilla

Página 47

Ramón Jesús Cuesta Pérez Curso 2011/2012

Proyecto Fin de Carrera Implementación de algoritmos MPC con restricciones en mbed NXP LPC1768

Con

a = b + u 1TH,,

f1 = f0 + u 1TH1 + b u 2

(

)

  1u −u    1U − T 1u − 1u (t − 1)  c = − 1U + T 1u + 1u (t − 1)   1y − f − Gu    − 1y + f + Gu   

 I NxN   T    R =  −T     G   − G 

Este es el problema QP que vamos a resolver aplicando el algoritmo de Lemke tras usar

el concepto de Problema Lineal Complementario. Definimos en primer lugar un vector de variables de holgura para convertir la restricción Rx ≤ c en una restricción de igualdad. Denominando v2 a este vector el problema QP queda:

x* = arg min x

1 T x Hx + ax + f1 2

sujeto a Rx + v2 = c

x≥0 Para poder usar el concepto de Problema Lineal Complementario para hallar x* imponemos las condiciones de Kuhn Kuhn-Tucker Tucker de este problema. Matemáticamente, dichas condiciones para un problema de optimización de la forma min f(x) s.a. g(x) ≤ 0

h(x) = 0

son ∇f ( x*) + λ∇g ( x*) + µ∇h( x*) = 0 , donde λ y µ son los vectores de multiplicadores de Lagrange de las respectivas restricciones.

Si denominamos v1 al vector de multiplicadores de Lagrange de la restricción −x ≤ 0 y v al de la restricción Rx + v2 − c = 0, imponiendo las condiciones de Kuhn-Tucker Kuhn tenemos: Hx + a + RTv − v1 = 0 Dpto. Ing. de Sistemas y Automática Universidad de Sevilla

Página 48

Ramón Jesús Cuesta Pérez Curso 2011/2012

Proyecto Fin de Carrera Implementación de algoritmos MPC con restricciones en mbed NXP LPC1768

xTv1 = 0 , vTv2 = 0 (han de ser vectores l.i.) x, v, v1, v2 ≥ 0 La primera ecuación puede reescribirse como −Hx − RTv + v1 = a En definitiva, para hallar x** se ha de verificar Rx + v2 = c −Hx − RTv + v1 = a xTv1 = 0 , vTv2 = 0 x, v, v1, v2 ≥ 0 Estas expresiones se pueden expresar de forma matricial como

 I mxm O  Nxm

OmxN I NxN

v2  R   v1   c  = − H   v  a    x

Omxm − RT

Lo cual es finalmente susceptible de ser expresado como un Problema Lineal Complementario s − Mz = q , sTz = 0 , s, z ≥ 0 O M=  T R

− R , H 

c  q=   , a 

s

z −M

con

v2  s=    v1 

v  z=    x

,

Así pues, la tabla s

Imxm

z0 −1

q

quedaría: v2 v2 v1

Imxm ONxm

v1 OmxN INxN

v Omxm −RT

x R −H

z0 − −1 − −1

c a

Sobre esta tabla será sobre la que empecemos a iterar al aplicar el algoritmo de Lemke.

Dpto. Ing. de Sistemas y Automática Universidad de Sevilla

Página 49

Ramón Jesús Cuesta Pérez Curso 2011/2012

Proyecto Fin de Carrera Implementación de algoritmos MPC con restricciones en mbed NXP LPC1768

Obtención de la solución mediante pivoteo en la tabla

Una vez hemos construido la tabla obtendremos la solución pivotando en la tabla durante un numero finito de iteraciones (siempre y cuando H sea definida positiva). A continuación expondremos cómo realizar dichas operaciones de pivoteo.

En primer lugar, como vimos en la explicación del LPC, tenemos que adaptar la tabla de modo que comencemos con la solución inicial inic z0 = max(−qi), z = 0 y s = q + 1z0. En la tabla que hemos construido (ver página anterior) todas las variables básicas son las correspondientes a s,, por lo tanto tenemos s = q, z = 0, z = 0 y z0 = 0,, lo cual sería una c  solución básica factible complementaria si s = q =   ≥ 0, lo cual en general no será a  cierto (tanto c como a tendrán alguna componente negativa).

Así pues, para hacer z0 = max(−qi) tendremos que introducir la columna de z0 en la base. Para ello tomaremos como pivote la fila correspondiente con el valor más negativo de q.. De este modo al hacer este primer pivoteo nos quedará q ≥ 0 y una tabla tal que así: v2 v2

v1

v1

v

1

0

...

*

...

0

0 ... ...

1 ... ...

... ... ...

* * *

... ... ...

0 ... ...

... 0

... 0

... ...

* *

... ...

... 1

x

z0 c

Omxm

R

−RT

−H

ei a

Ahora tenemos ya la solución inicial y empezamos a iterar hasta que z0 salga de la base. El criterio para cambiar de pivote tras cada iteración es el siguiente: la columna del nuevo pivote será igual a la suma de la fila del pivote anterior más (menos) el número total de filas si la columna del pivote superior es mayor (menor) (menor) el número total de filas. De esto modo nunca nos saldremos de la tabla. La fila del nuevo pivote será aquella de la columna anteriormente seleccionada en la que, si el elemento es positivo, el cociente qi/elemento es mínima para dicha columna. Dpto. Ing. de Sistemas y Automática Universidad de Sevilla

Página 50

Ramón Jesús Cuesta Pérez Curso 2011/2012

Proyecto Fin de Carrera Implementación de algoritmos MPC con restricciones en mbed NXP LPC1768

Nótese quee el proceso de pivoteo es análogo al que se realiza en otros algoritmos como el Simplex, es decir, consiste en operar por filas en la matriz de modo que la columna en la que se encuentre el pivote se convierta en el vector ei, donde i es la fila en la que se halla el pivote.

Cuando, tras una iteración, z0 haya salido de la base habremos obtenido la solución al LPC y, por tanto, x*. *. Para identificar el valor de las componentes de x* lo único que tenemos que hacer es observar qué valor de qi corresponde a cada xi básico. Si alguna columna xi no pertenece a la base el valor de dicha componente será cero. Una vez realizada esta sencilla identificación tendremos nuestro vector x*, al cual deberemos en el algoritmo MPC deshacer el cambio de variable vari u = 1 u + x para obtener el vector u** de incrementos de acción de control óptimos, que es el objetivo del solver.

4.2 Programación del algoritmo en C++ Una vez explicado el fundamento teórico del algoritmo, así como la manera de proceder una vez construida la tabla para hallar la solución, incluiremos la programación en lenguaje C++ del algoritmo. El código ha sido debidamente comentado para que pueda ser fácilmente cilmente seguido en conjunción con las explicaciones anteriores.

Sin más dilación presentamos la función qp_lemke,, a la cual llamamos en el algoritmo MPC cuya estructura se explicó en el capítulo anterior. Junto con ella se presentan el resto de funciones a la que ésta, a su vez, llama.

Dpto. Ing. de Sistemas y Automática Universidad de Sevilla

Página 51

Ramón Jesús Cuesta Pérez Curso 2011/2012

Proyecto Fin de Carrera Implementación de algoritmos MPC con restricciones en mbed NXP LPC1768

Función qp_lemke double *qp_lemke(double double *u, const double H[][N], double *a, const double R[][N], double *c, int dim, int a_rows) { int i,j; /* contadores */ double **tab; /* tabla de Lemke */ int m_dim; /* dimensión de m */ int *act_base; /* vector donde almacenar las variables basicas */ /* Inicializamos a cero el vector que devolveremos como solución */ vec_null(u,dim); /* Si q >= 0 entonces s = q y hemos terminado */ if (vec_pos_p(a,dim) (a,dim) && vec_pos_p(c,a_rows)) return(u); else { /* Número de filas que tendrá la tabla */ m_dim=a_rows+dim; /* Creamos la matriz que representa la tabla dinámicamente */ tab = mat_din(m_dim,2*m_dim+2); (m_dim,2*m_dim+2); /* Creamos el vector donde guardamos las variables básicas */ act_base = vec_din_int(m_dim); vec_din_int /* Construimos la tabla */ for (i=0;i