Teoría de Estabilidad y Control Isaac A. García

eines 52

Teoría de estabilidad y control

Isaac A. García Seminari de Sistemes Dinàmics Departament de Matemàtica Universitat de Lleida

ISBN: 978-84-8409-422-7 © Edicions de la Universitat de Lleida, 2005 © El autor

Maquetación: Servei de Publicacions (UdL) Diseño cubierta: cat & cas Impresión: Cargraphics

La reproducción total o parcial de esta obra por cualquier procedimiento, incluidos la reprografía y el tratamiento informático, y la distribución de ejemplares mediante alquiler o préstamo público, queda rigurosamente prohibida sin la autorización de los titulares del copyright, y será sometida a las sanciones establecidas por la ley.

EINES es una col·lección del Institut de Ciències de l’Educació de la Universitat de Lleida.

A mi mujer Susanna y a mi hija Nadia

Índice general IX

Pr´ ologo 1. Introducci´ on y Ejemplos 1.1. Breve historia de control autom´ atico . . . . . . 1.2. Algunos ejemplos f´ısicos . . . . . . . . . . . . . 1.2.1. Part´ıcula en movimiento unidimensional 1.2.2. Termostato y transferencia de calor . . .

. . . .

1 1 3 3 3

2. Control Cl´ asico y Transformada de Laplace 2.1. El problema clásico de control . . . . . . . . . . . . . . . . . . . . 2.2. Transformada de Laplace y funci´ on de transferencia . . . . . . .

5 5 5

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

3. Soluci´ on de Sistemas Lineales 3.1. Solución espectral de sistemas lineales . . . . . . . . . . . . . . 3.1.1. La matriz exponencial . . . . . . . . . . . . . . . . . . . 3.1.2. El teorema de Cayley-Hamilton y la matriz exponencial 3.2. Solución de sistemas controlados . . . . . . . . . . . . . . . . . 3.2.1. Existencia y unicidad de la soluci´ on . . . . . . . . . . . 3.3. Relación entre espacio de estados y el control clásico . . . . . . 3.4. Problemas resueltos . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . .

7 7 8 11 12 13 13 16

4. Sistemas de Control Lineal 4.1. Introducción . . . . . . . . 4.2. Controlabilidad . . . . . . 4.3. Equivalencia Algebraica . 4.4. Problemas resueltos . . .

. . . .

23 23 23 31 32

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

5. Observabilidad 37 5.1. Observabilidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 5.2. Problemas resueltos . . . . . . . . . . . . . . . . . . . . . . . . . 40 6. Realimentaci´ on Lineal 6.1. Definición y preliminares . . . . . . . . . . . . . . 6.1.1. Realimentación y funci´ on de transferencia 6.2. Realimentación frente a control precalculado . . . 6.2.1. Sensibilidad a las condiciones iniciales . . VII

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

43 43 46 46 48

´ INDICE GENERAL

VIII

6.2.2. Sensibilidad a perturbaciones externas . . . . . . . . . . . 6.3. Problemas resueltos . . . . . . . . . . . . . . . . . . . . . . . . .

48 49

7. Estabilidad 7.1. Introducción y definiciones . . . . . . . . . . . . . . . . 7.2. Estabilidad en sistemas lineales . . . . . . . . . . . . . 7.3. Teor´ıa de Liapunov de la estabilidad . . . . . . . . . . 7.3.1. Aplicación a sistemas lineales . . . . . . . . . . 7.3.2. Estabilidad mediante linearizaci´ on . . . . . . . 7.4. Estabilidad y control . . . . . . . . . . . . . . . . . . . 7.4.1. Estabilidad entrada–salida . . . . . . . . . . . . 7.4.2. Estabilización por realimentaci´ on lineal . . . . 7.4.3. Linealización de sistemas de control no lineales 7.5. Problemas resueltos . . . . . . . . . . . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

53 53 54 56 58 59 61 61 63 63 64

8. C´ alculo de Variaciones 8.1. Un ejemplo: braquistocrona . . . . . . . . . . . . . . 8.2. Introducción . . . . . . . . . . . . . . . . . . . . . . . 8.3. Ecuaciones de Euler–Lagrange . . . . . . . . . . . . 8.3.1. Lema fundamental del c´ alculo de variaciones 8.3.2. Ecuaciones de Euler–Lagrange . . . . . . . . 8.3.3. Integrales primeras en casos simples . . . . . 8.4. Extremos con restricciones . . . . . . . . . . . . . . . 8.5. Apéndice: Multiplicadores de Lagrange . . . . . . . . 8.5.1. Problemas isoperimétricos . . . . . . . . . . . 8.6. Problemas resueltos . . . . . . . . . . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

75 75 77 77 77 78 80 81 82 82 85

´ 9. Control Optimo 9.1. Introducción . . . . . . . . . . . . . . . 9.2. Control óptimo: método hamiltoniano 9.3. El regulador lineal . . . . . . . . . . . 9.4. Teor´ıa de Pontryagin . . . . . . . . . . 9.5. Problemas resueltos . . . . . . . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

95 . 95 . 97 . 99 . 101 . 103

. . . . . . . . . . . . .

115 115 116 117 117 119 119 120 121 122 123 123 125 126

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

10.Pr´ acticas 10.1. Control de la temperatura de una c´ amara . . . . . . . . . . 10.1.1. Realización de la pr´ actica . . . . . . . . . . . . . . . 10.2. Dinámica de satélites de comunicación . . . . . . . . . . . . 10.2.1. Un modelo matemático para la din´ amica de satélites 10.2.2. Realización de la pr´ actica . . . . . . . . . . . . . . . 10.3. El péndulo doble invertido . . . . . . . . . . . . . . . . . . . 10.3.1. Realización de la pr´ actica . . . . . . . . . . . . . . . onico . . . . . . . . . . . . . . . 10.4. Estabilidad en el péndulo c´ 10.4.1. Realización de la pr´ actica . . . . . . . . . . . . . . . 10.5. Un modelo de gr´ ua con servomecanismo y regulador . . . . 10.5.1. Realización de la pr´ actica . . . . . . . . . . . . . . . 10.6. El regulador centr´ıfugo de Watt . . . . . . . . . . . . . . . . 10.6.1. Realización de la pr´ actica . . . . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

IX

Pr´ ologo El presente libro de teor´ıa y problemas corresponde a los temas básicos de un primer curso de introducci´ on al la teor´ıa de control y estabilidad en variables de estado. El autor imparte dicha asignatura en la titulaci´ on de Ingenier´ıa Técnica Industrial aunque el libro es igualmente recomendable para estudiantes de otras titulaciones técnicas. El objetivo principal es que el alumno disponga de un material preliminar para el curso, con los resultados principales y algunas demostraciones de estos. Además es interesante que el alumno pueda seguir paso a paso la resolución de numerosos problemas de los temas mencionados como ejemplificación de los conceptos y resultados teóricos, complementando as´ı los realizados en clase. Se ha procurado presentar las soluciones en la forma m´ as práctica y directa. Me gustar´ıa que este libro facilitase el aprendizaje de la asignatura y, agradecer´ıa cualquier sugerencia o comentario que pueda mejorarlo dirigiéndose a la siguiente dirección electrónica: [email protected].

Dr. Isaac A. Garc´ıa, Septiembre de 2004

Cap´ıtulo 1

Introducci´ on y Ejemplos 1.1.

Breve historia de control autom´ atico

En este cap´ıtulo introductorio revisaremos, en primer lugar, una breve introducci´ on hist´ orica de la teor´ıa del control autom´ atico. El llamado control por retroalimentaci´ on es un mecanismo básico a través del cual sistemas de naturaleza muy diferente (mecánicos, eléctricos, biológicos, etc...) mantienen su equilibrio. Por ejemplo, un cambio en la temperatura corporal de 1 grado es, en general, una se¯ nal de alg´ un tipo de enfermedad. De hecho, de la teor´ıa de la evoluci´ on de las especies de Darwin se desprende que la retroalimentaci´ on sobre largos periodos temporales es responsable de tal evoluci´ on. Se puede definir el control por retroalimentaci´ on como el uso de diferentes se¯ nales, determinadas a partir de la comparaci´ on del estado actual del sistema y del estado deseado, para controlar el sistema. Un ejemplo cotidiano de sistema controlado por retroalimentaci´ on es el control de la velocidad de un coche mediante el uso de la diferencia entre la velocidad actual y la deseada para variar el flujo de gasolina. El hecho de que la salida del sistema sea usada para regular la entrada del sistema da lugar a lo que se llama sistema de control en lazo cerrado. El control por retroalimentaci´ on es una disciplina de la ingenier´ıa y, como tal, está estrechamente relacionada con diversos problemas aplicados que la humanidad ha querido resolver a lo largo de su historia. M´ as concretamente, existi´ o una época muy importante para el desarrollo de la teor´ıa de control que se comprende entre la Revoluci´ on Industrial y la primera y segunda Guerra Mundial. El control necesit´ o adquirir el lenguaje de las matem´ aticas para poder expresarse correctamente. De este modo J.C. Maxwell introdujo el análisis riguroso de la teor´ıa de control en 1868. Existen diversos periodos remarcables: época primitiva desde 1868–1900; época cl´ asica desde 1900–1960; época moderna desde 1960 hasta la actualidad.

1

2

Introducci´ on y Ejemplos

La Revoluci´ on Industrial en Europa introdujo nuevas m´ aquinas que no se pod´ıan regular con la mano. De este modo se iniciaron los primeros dispositivos de control autom´ atico. Algunos de ellos son los reguladores de temperatura, de presión, de velocidad, etc...

En 1840, el astr´ onomo G.B. Airy desarroll´ o un sistema de retroalimentaci´ on para colocar su telescopio mediante un control de la velocidad con el fin de compensar la rotación terrestre. Airy descubri´ o que un dise¯ no poco adecuado del sistema de lazo cerrado de retroalimentación produc´ıa oscilaciones indeseables en el sistema. De este modo introdujo el concepto de inestabilidad y la herramienta de las ecuaciones diferenciales en su análisis. J.C. Maxwell (1868) fue el primero en linearizar las ecuaciones diferenciales del movimiento para hallar la ecuaci´ on caracter´ıstica del sistema y analizar los efectos de los parámetros del sistema en su estabilidad. En 1877, E.J. Routh muestra una técnica numérica para determinar cuando una ecuación caracter´ıstica tiene ra´ıces estables. Ese mismo a¯ no, I.I. Vishnegradsky analizó la estabilidad de los reguladores usando ecuaciones diferenciales. En 1892, A.M. Liapunov introdujo un trabajo b´ asico en teor´ıa de control, estudiando la estabilidad de sistemas no lineales usando una generalizaci´ on del concepto de energ´ıa. Entre los a¯ nos 1892–1898, O. Heaviside invent´ o el cálculo operacional e introdujo la llamada actualmente funci´ on de transferencia.

La teor´ıa de control cl´ asica usa básicamente técnicas de dominio frecuencial en el plano complejo. Se basa en metodos de transformadas y es aplicable principalmente a sistemas lineales autónomos. Los métodos de Nyquist y Bode analizan la magnitud y la fase de la respuesta en frecuencia del sistema. Tiene la gran ventaja de que la respuesta en frecuencia del sistema puede ser medida experimentalmente y se puede pues calcular la funci´ on de transferencia. No se necesita bajo este punto de vista una descripci´ on interna exacta del sistema en el dominio temporal, es decir, sólo es importante el comportamiento entrada– salida del sistema.

Sin embargo esta teor´ıa de control cl´ asica es dif´ıcil de aplicar a sistemas con varias entradas y m´ ultiples salidas. La teor´ıa de control moderna es fundamentalmente una teor´ıa de dise¯ no en el dominio temporal (a diferencial de la cl´ asica que es en el dominio frecuencial). Se requiere un modelo en el espacio de estados del sistema que se desea controlar. Esta visión moderna es la que se dará en este libro.

1.2 Algunos ejemplos f´ısicos

1.2.

3

Algunos ejemplos f´ısicos

1.2.1.

Part´ıcula en movimiento unidimensional

Consideremos una part´ıcula que se mueve en una recta y sea s(t) su distancia al origen en funci´ on del tiempo. Supongamos que sobre la part´ıcula act´ ua una on pero también fuerza por unidad de masa u1 (t) que le produce una aceleraci´ ńicas cantiact´ ua una fuerza de rozamiento u2 (t) por unidad de masa. Si las u dades cinemáticas de interés son la posición de la part´ıcula x1 (t) = s(t) y su ˙ aplicando las leyes de la mecánica clásica, el estado de velocidad x2 (t) = s(t), la part´ıcula (x1 (t), x2 (t)) verifica el sistema de ecuaciaciones diferenciales x˙ 1 = x2 , x˙ 2 = u1 (t) − u2 (t) . Por supuesto, este sistema se puede escribir en notación matricial de la forma x˙ = Ax + Bu , siendo

x=

x1 x2

, u=

u1 u2

, A=

0 1 0 0

(1.1)

, B=

0 1

0 −1

.

Se puede preguntar c´ omo puede llegar la part´ıcula a un cierto punto fijo en el menor tiempo posible o bien con un consumo m´ınimo de una cierta variable. Matem´ aticamente se trata de determinar las variables de control u1 (t) y u2 (t) con un cierto objetivo.

1.2.2.

Termostato y transferencia de calor

La temperatura interior Ti de un horno se quiere controlar variando la entrada de calor u en las paredes. Definamos las capacidades calor´ıficas de la pared y del interior del horno como cp y ci respectivamente. Denotemos por Ai y Ae las áreas interiores y exteriores de la pared del horno y por Ri y Re los coeficientes de radiaci´ on de las superficies interior y exterior del horno. Ignorando otros efectos, si las temperaturas de la pared y del exterior del horno son Tp y Te , se tiene, tomando Ti > Tp > Te , para la pared cp T˙p = −Ae Re (Tp − Te ) − Ai Ri (Tp − Ti ) + u , y para el interior del horno ci T˙i = Ai Ri (Tp − Ti ) . Suponiendo Te constante, si definimos las variables de estado x1 = Tp − Te y x2 = Ti − Te se tiene x˙ 1 x˙ 2

1 u 1 , = T˙p − T˙e = T˙p = − Ae Re x1 − Ai Ri (x1 − x2 ) + cp cp cp 1 = T˙i − T˙e = T˙i = Ai Ri (x1 − x2 ) . ci

4

Introducci´ on y Ejemplos

Estas ecuaciones se pueden escribir de la forma x˙ = Ax + bu, siendo x = (x1 , x2 )T , −(Ae Re + Ai Ri )/cp (Ai Ri )/cp 1/cp A= , b= . (Ai Ri )/ci −(Ai Ri )/ci 0 El objetivo puede ser analizar un controlador que regule, mediante una v´ alvula por ejemplo, la cantidad de calor u(t) que se debe suministrar dependiendo de la temperatura que marquen los term´ ometros.

Cap´ıtulo 2

Control Cl´ asico y Transformada de Laplace 2.1.

El problema cl´ asico de control

La teor´ıa de control cl´ asica estudia sistemas con una entrada escalar u(t) ∈ R y una salida escalar z(t) ∈ R. Si se requiere que z(t) esté lo más cerca posible (en un cierto sentido) de una cierta se¯ nal de referencia r(t), el sistema de control se llama servomecanismo. Un caso particular es cuando r es constante, entonces el sistema de control se llama regulador. El modelo de la teor´ıa de control cl´ asica es una ecuación diferencial lineal de orden n a coeficientes constantes, es decir, z(t) verifica la ecuación z (n) + k1 z (n−1) + · · · + kn−1 z + kn z = β0 u(m) + β1 u(m−1) + · · · + βm u , (2.1) donde los supra´ındices indican derivadas respecto de la variable independiente t y ki y βi son constantes.

2.2.

Transformada de Laplace y funci´ on de transferencia

Definici´ on 1. Sea z : R+ → R una funci´ on real y continua a trozos. Se define la transformada de Laplace z¯(s) de z(t) de la forma ∞ z¯(s) = L{z(t)} = z(t) exp(−st) dt . (2.2) 0

Algunas de las propiedades de esta transformaci´ on son las siguientes: n n Linealidad: L{ i=1 ci zi (t)} = i=1 ci L{zi (t)} para todo ci ∈ R. 5

6

Control Cl´ asico y Transformada de Laplace Derivaci´ on: L{z (j) (t)} = sj z¯ − sj−1 z(0) − sj−2 z (1) (0) − · · · − sz (j−2) (0) − z (j−1) (0) . (2.3)

Tomando transformadas de Laplace en ambos miembros de la ecuación (2.1) y teniendo en cuenta la propiedad (2.3) se tiene k(s)¯ z (s) = β(s)¯ u(s), siendo los polinomios k(s) = sn + k1 sn−1 + · · · + kn−1 s + kn , β(s) = β0 sm + β1 sm−1 + · · · + βm−1 s + βm .

(2.4) (2.5)

Es habitual presentarlo de la forma z¯(s) = g(s)¯ u(s) , siendo g(s) =

β(s) , k(s)

(2.6)

(2.7)

una funci´ on racional conocida como funci´ on de transferencia. Es preciso resaltar que este método sólo es aplicable a sistemas lineales con coeficientes constantes. Nosotros seguiremos ideas diferentes para atacar los problemas b´ asicos del control.

Cap´ıtulo 3

Soluci´ on de Sistemas Lineales 3.1.

Soluci´ on espectral de sistemas lineales

Consideraremos sistemas lineales sin variables de control, es decir, sistemas del tipo x˙ = Ax , (3.1) con A ∈ Mn (R) una matriz cuadrada de orden n y x ∈ Rn . Tomaremos un problema de valor inicial o de Cauchy, de este modo impondremos que x(0) = x0 . Supondremos que todos los valores propios λi con i = 1, . . . , n, de la matriz A son distintos1 . Sean wi , con i = 1, . . . , n los vectores propios de A asociados a los valores propios λi , es decir, Awi = λi wi . Puesto que λi = λj con i = j, se sabe que el conjunto de vectores propios {w1 , . . . , wn } son linealmente independientes y en particular forman base de Rn . En consecuencia se puede expresar la solución de (3.1) de la forma n ci (t)wi , x(t) = i=1

siendo ci (t) funciones escalares del tiempo. Derivando respecto de t la ecuación anterior y sustituyendo en (3.1) se obtiene n i=1

c˙i (t)wi = A

n

ci (t)wi =

n

λi ci (t)wi .

i=1

i=1

1 Esta

no es una restricci´ on fuerte en los sistemas que aparecen en las aplicaciones reales puesto que una peque¯ na perturbaci´ on en los coeficientes aij de la matriz A es suficiente para separar las ra´ıces iguales del polinomio caracter´ıstico de A. Recordemos que, habitualmente dichos coeficientes provienen de medidas experimentales y son, por lo tanto, conocidos hasta un cierto grado de precisi´ on.

7

8

Soluci´ on de Sistemas Lineales

Puesto que el conjunto de vectores propios {w1 , . . . , wn } son linealmente independientes, se obtiene c˙i = λi ci , i = 1, . . . , n , cuya solución es ci (t) = exp(λi t) ci (0) , i = 1, . . . , n . Se tiene, en definitiva, que x(t) =

n

ci (0) exp(λi t)wi .

(3.2)

i=1

Sin embargo, interesa dar la soluci´ on x(t) en funci´ on de x(0). Para ello, definamos la matriz W = [w1 , . . . , wn ] ∈ Mn (R) cuyas columnas son los vectores wi . Sea W −1 la matriz inversa de W , que siempre existe puesto que rangW = n, y definamos sus filas como los vectores fila vi con i = 1, . . . , n, es decir, ⎤ ⎡ v1 ⎥ ⎢ W −1 = ⎣ ... ⎦ ∈ Mn (R) . vn Puesto que W −1 W = In , siendo In la matriz identidad de orden n, se verifica on (3.2) por la izquierda vi wi = 1 y vi wj = 0 si i = j. Multiplicando la ecuaci´ por el vector vj y tomando finalmente t = 0 se obtiene vj x(0) = cj (0). De este modo, la soluci´ on de (3.1) viene dada por x(t) =

n

(vi x(0)) exp(λi t)wi .

(3.3)

i=1

Esta expresión se conoce como la forma espectral de la solución de (3.1).

3.1.1.

La matriz exponencial

Veamos en esta sección una forma alternativa para obtener la soluci´ on de (3.1). Esta nueva forma evita el c´ alculo de los vectores propios wi de la matriz A necesarios en la forma espectral (3.3) de la solución de (3.1). La idea se basa en generalizar el caso escalar n = 1, con lo cual A es un escalar en (3.1) y su solución viene dada por x(t) = exp(At)x0 .

(3.4)

Para conseguir dicha generalizaci´ on, definamos la matriz exponencial exp(At) =

∞ k t k=0

k!

Ak = In + tA +

t2 2 A + ··· . 2!

(3.5)

3.1 Soluci´ on espectral de sistemas lineales

9

Esta serie infinita de matrices es convergente2 para cualquier A ∈ Mn (R) y para todo t ∈ R debido al hecho de que exp(zt) converje para cualquier pareja de escalares finitos z y t. Es claro que, a partir de su definici´ on (3.5), se tiene exp(O) = In siendo O la matriz nula de orden n. Además, derivando término a término (3.5) respecto de t se verifica 1 2 3 1 2 2 d 2 exp(At) = A + tA + t A + · · · = A In + tA + t A + · · · dt 2! 2! = A exp(At) , con lo cual

d exp(At) = A exp(At) , dt de modo que podemos concluir que (3.4) representa la soluci´ on de (3.1). La solución del problema de Cauchy x˙ = Ax , x(t0 ) = x0 ,

(3.6)

se encuentra, en la literatura de control, de la forma x(t) = Φ(t, t0 )x0 ,

(3.7)

Φ(t, t0 ) = exp[A(t − t0 )] ,

(3.8)

donde es una matriz cuadrada de orden n llamada matriz de transici´ on de estados. Es f´ acil demostrar que dicha matriz verifica las siguientes propiedades: d Φ(t, t0 ) = AΦ(t, t0 ) , dt Φ(t, t) = In , Φ(t0 , t) = Φ−1 (t, t0 ), Φ(t, t0 ) = Φ(t, t1 )Φ(t1 , t0 ).

(3.9) (3.10) (3.11) (3.12)

Observar que, en particular, la propiedad (3.11) implica (exp(At))−1 = exp(−At) .

(3.13)

Resumimos a continuaci´ on algunas de las propiedades m´ as importantes de la matriz exponencial. A y exp(A) conmutan, es decir, A exp(A) = exp(A)A. 2 Una sucesi´ on de matrices {Ak } converge hacia A si l´ımk→∞ Ak − A = 0 para cualquier on de sumas parciales {Sn } norma matricial. La serie matricial ∞ k=0 Ak converge si la sucesi´ n con Sn = k=0 Ak converge cuando n → ∞. Si la funci´ on escalar f (z) se puede representar k on por una serie de potencias f (z) = ∞ k=0 ck x convergente para |z| < R entonces la funci´ k converge si todos los valores propios λ de la matriz cuadrada c A matricial f (A) = ∞ i k=0 k A verifican |λi | < R.

10

Soluci´ on de Sistemas Lineales exp(On ) = In , siendo On e In las matrices nulas e identidad de orden n. Si A y B conmutan, es decir, AB = BA entonces exp(AB) = exp(A) exp(B). (exp(A))−1 = exp(−A). (exp(A))T = exp(AT ). d exp(tA)/dt = A exp(tA). Sea P una matriz no singular tal que B = P −1 AP . Entonces exp(B) = P −1 exp(A)P . Si D = diag{λ1 , . . . , λn } es una matriz diagonal con elementos λi en la diagonal principal, entonces exp(D) = diag{exp(λ1 ), . . . , exp(λn )} es una matriz diagonal con elementos exp(λi ) en la diagonal principal.

Las dos u ´ltimas propiedades permiten calcular de forma eficaz la exponencial de una matriz A diagonalizable. En efecto, si A es diagonalizable existe una matriz de cambio de base P no singular tal que D = P −1 AP , siendo D una matriz diagonal con los valores propios de A en la diagonal principal. M´ etodos para calcular la matriz exponencial (i) Si todos los valores propios λi de A son distintos, la f´ ormula de Sylvester viene dada por n Zk exp(λk t) , (3.14) exp(At) = k=1

siendo Zk =

n A − λj In , λk − λj j=1

(3.15)

j=k

matrices constantes que sólo dependen de A y de sus valores propios λj , pero no de sus vectores propios wi . Notar la semejanza entre la fórmula de Sylvester y la fórmula de interpolaci´ on polinomial de Lagrange utilizada en métodos numéricos. (ii) Otro método alternativo para calcular la matriz exponencial de A cuando todos sus valores propios λi son diferentes es el siguiente. exp(At) = r(A) ,

(3.16)

siendo r(λ) un polinomio de grado menor o igual que n − 1 cuyos coeficientes son funciones del tiempo t obtenidos de la forma r(λi ) = exp(λi t) , i = 1, . . . , n . La demostración de este hecho es una consecuencia de la siguiente sección.

3.1 Soluci´ on espectral de sistemas lineales

11

(iii) El cálculo de la matriz exponencial es tedioso en muchas ocasiones, de modo que los manipuladores algebraicos son de gran ayuda en dicho c´ alculo. Por ejemplo, con el programa Mathematica, el comando MatrixExp[A] calcula la matriz exponencial de la matriz A.

3.1.2.

El teorema de Cayley-Hamilton y la matriz exponencial

Sea A ∈ Mn (R). En este método no es necesario suponer ni que los valores propios de A son distintos ni que A diagonaliza. En primer lugar, denotamos por k(λ) el polinomio caracter´ıstico de A, y sean λi los valores propios de A con multiplicidad algebraica mi para i = 1, . . . , k, siendo k ≤ n. Por el Teorema de Cayley-Hamilton se sabe que k(A) = 0. Se efect´ ua la descomposición en fracciones simples k

pi (λ) 1 = , k(λ) (λ − λi )mi i=1 siendo pi (λ) polinomios de grado menor o igual que mi − 1. Definiendo los on se polinomios qi (λ) = k(λ)/(λ − λi )mi , con i = 1, . . . , k, la anterior ecuaci´ reescribe de la forma k pi (λ)qi (λ) . (3.17) 1= i=1

Evaluando este polinomio en A se tiene la identidad In =

k

pi (A)qi (A) .

i=1

Por otra parte, observemos que exp(λi tIn ) exp[(A − λi In )t] = exp(λi t)In exp[(A − λi In )t] ∞ j t (A − λi In )j , = exp(λi t) j! j=0

exp(At) =

de modo que, multiplicando por la izquierda ambos miembros por qi (A) se obtiene la siguiente igualdad con suma finita qi (A) exp(At) = exp(λi t)

m i −1 j j=0

t qi (A)(A − λi In )j , j!

(3.18)

donde se ha teniendo en cuenta el Teorema de Cayley-Hamilton de modo que qi (A)(A − λi In )j = k(A)(A − λi In )j−mi = 0. Multiplicando finalmente la

12


ecuación (3.18) por la izquierda por pi (A), sumando para todo i desde 1 hasta k y teniendo en cuenta (3.17) se obtiene la siguiente expresi´ on de la matriz exponencial

exp(At) =

k

⎡ ⎣exp(λi t)pi (A)qi (A)

i=1

⎤ t (A − λi In )j ⎦ . j!

m i −1 j j=0

(3.19)

Nota 2. Observar que los grados de los polinomios pi y qi son mi − 1 y n − mi respectivamente, de modo que, en el caso particular de que todos los valores propios de A sean distintos, es decir, mi = 1, entonces exp(At) = r(A) siendo r un polinomio de grado menor o igual que n − 1 verificando adem´ as r(λi ) = exp(λi t) para i = 1, . . . , n.

3.2.

Soluci´ on de sistemas controlados

Consideremos el sistema controlado x˙ = Ax + Bu , x(0) = x0 ,

(3.20)

con A ∈ Mn (R), B ∈ Mn×m (R) matrices constantes, x(t) ∈ Rn y u(t) ∈ Rm . Multiplicando por la izquierda ambos miembros de (3.20) por la matriz exp(−At) se obtiene exp(−At)x˙ − exp(−At)Ax = exp(−At)Bu . Utilizando ahora que las matrices exp(−At) y A conmutan y teniendo en cuenta la propiedad (3.9), es decir, d [exp(−At)] = −A exp(−At) , dt se llega a que d [exp(−At)x] = exp(−At)Bu . dt Integrando dicha ecuaci´ on respecto de t y teniendo en cuenta (3.13) se tiene

t x(t) = exp(At) x0 + exp(−Aτ )Bu(τ ) dτ .

(3.21)

0

Si se utiliza la condici´ on inicial x(t0 ) = x0 , integrando entre t0 y t se obtiene

t x(t) = Φ(t, t0 ) x0 + Φ(t0 , τ )Bu(τ ) dτ . t0

(3.22)

3.3 Relaci´ on entre espacio de estados y el control cl´ asico

3.2.1.

13

Existencia y unicidad de la soluci´ on

Presentamos a continuaci´ on un resultado sobre la existencia y unicidad del problema de Cauchy asociado a una sistema lineal controlado por un control continuo a trozos de gran utilidad en aplicaciones a la ingenier´ıa. Definici´ on 3. El control u(t) : [t0 , t1 ] → Rm es continuo a trozos en el intervalo on finita de dicho intervalo t0 = s0 < s1 < · · · < sh = [t0 , t1 ] si existe una partici´ on u(t) es continua en cada subintervalo abierto (si−1 , si ) para t1 , tal que la funci´ i = 0, 1, . . . , h. Adem´ as, para cada i, l´ımt→s+ u(t) existe y es finito si i = h y i l´ımt→s− u(t) existe y es finito si i = 0. i

Definici´ on 4. Una soluci´ on del problema de Cauchy x˙ = Ax + Bu(t) con on x(t) : x(t0 ) = x0 y u(t) continua a trozos en el intervalo [t0 , t1 ] es una funci´ [t0 , t1 ] → Rn continua y con derivada continua en [t0 , t1 ] excepto en los puntos de discontinuidad de u(t) que adem´ as verifique x˙ = Ax + Bu(t) donde exista. Teorema 5. El problema de Cauchy x˙ = Ax + Bu(t) con x(t0 ) = x0 y u(t) conńica soluci´ on x(t) : [t0 , t1 ] → tinua a trozos en el intervalo [t0 , t1 ] admite una u Rn que viene dada por (3.22).

3.3.

Relaci´ on entre espacio de estados y el control cl´ asico

Consideremos la ecuación diferencial escalar lineal de orden n a coeficientes constantes que aparece en teor´ıa de control cl´ asica z (n) + k1 z (n−1) + · · · + kn z = u(t) ,

(3.23)

siendo u(t) ∈ R la variable de control. Realizando el cambio de variable w1 = z , w2 = z˙ , . . . , wn = z (n−1) ,

(3.24)

y teniendo en cuenta que w˙ i = wi+1 para i = 1, . . . , n − 1, la ecuación (3.23) se puede escribir de la forma can´ onica controlable3 w˙ = Cw + du , donde w = (w1 , . . . , wn )T , d = (0, . . . , 0, 1)T ⎛ 0 1 0 . . ⎜ 0 0 1 . . ⎜ ⎜ . . . . . C=⎜ ⎝ . . . . . −kn −kn−1 −kn−2 . .

(3.25)

and . 0 . . . . . . . −k1

⎞ ⎟ ⎟ ⎟ ∈ Mn (R) , ⎟ ⎠

(3.26)

3 En el cap´ ıtulo siguiente (ver Teorema 14) se entender´ a mejor el nombre “can´ onica controlable”.

14


es la llamada matriz “companion”. Observar que el polinomio caracter´ıstico k(λ) de C es k(λ) = det(λIn − C) = λn + k1 λn−1 + · · · + kn . Notar que este polinomio caracter´ıstico es justo el denominador de la funci´ on de transferencia que se obtiene si se aplican transformadas de Laplace a la ecuaci´ on (3.23). Hemos visto que la ecuación escalar (3.23) se puede escribir en forma matricial con las variables de estado. Es natural preguntarse si el inverso es siempre posible, es decir, nos preguntamos si cualquier sistema lineal x˙ = Ax + bu(t) con u escalar se puede escribir en la forma clásica (3.23) mediante un cambio de variables lineal. La respuesta la contiene el siguiente resultado. Teorema 6. El sistema lineal x˙ = Ax + bu con A ∈ Mn (R) matriz constante, b ∈ Rn vector columna y u(t) ∈ R puede ser transformado mediante un cambio onica controlable lineal w = T x con T ∈ Mn (R) no singular en la forma can´ (3.25) si y s´ olo si (3.27) rang[b, Ab, A2 b, . . . , An−1 b] = n .

Demostraci´ on. Demostremos la suficiencia, es decir, supongamos que se verifica (3.27) y veamos que existe la matriz T del enunciado. Realizando el cambio de variables w = T x en el sistema x˙ = Ax + bu se obtiene el sistema transformado w˙ = T AT −1 w + T bu .

(3.28)

Veamos que es posible construir la matriz T de modo que la ecuación (3.28) sea (3.25). Impongamos la siguiente forma de dicha T ⎡ ⎢ ⎢ ⎢ T =⎢ ⎢ ⎣

ξ ξA ξA2 .. .

⎤ ⎥ ⎥ ⎥ ⎥ ∈ Mn (R) , ⎥ ⎦

(3.29)

ξAn−1 siendo ξ ∈ Rn un vector fila tal que T sea no singular. De momento supondremos que tal ξ existe. Denotemos por si a la i-ésima columna de la matriz T −1 , es decir, sea T −1 = [s1 , . . . , sn ]. Consideremos la matriz ⎛ ⎜ ⎜ T AT −1 = ⎜ ⎝

ξAs1 ξA2 s1 .. .

ξAn s1

ξAs2 ξA2 s2 .. .

ξAn s2

··· ··· .. .

··· ··· .. .

···

···

ξAsn ξA2 sn .. .

ξAn sn

⎞ ⎟ ⎟ ⎟ . ⎠

3.3 Relaci´ on entre espacio de estados y el control cl´ asico

15

Comparando los coeficientes de esta matriz con los coeficientes que se obtienen de la identidad T T −1 = In , es decir, ⎞ ⎛ ξs2 ··· ··· ξsn ξs1 ⎜ ξAs1 ξAs2 ··· ··· ξAsn ⎟ ⎟ ⎜ ⎟ = In , ⎜ .. .. .. .. .. ⎠ ⎝ . . . . . ξAn−1 s1 ξAn−1 s2 · · · · · · ξAn−1 sn se deduce que la i-ésima fila de la matriz T AT −1 coincide con la i + 1-ésima fila de la matriz identidad In para i = 1, . . . , n − 1. Por lo tanto, la matriz T AT −1 coincide con la matriz companion C dada en (3.26) con ki = −ξAn sn−i+1 . Para que (3.28) coincida con (3.25) s´ olo nos falta imponer que T b = d. Sustituyendo en esta igualdad la matriz (3.29) y teniendo en cuenta que d = (0, . . . , 0, 1)T se obtiene ξb = 0 , ξAb = 0, . . . , ξAn−2 b = 0 , ξAn−1 b = 1 ,

(3.30)

o, de forma equivalente, ξ[b, Ab, . . . , An−1 b] = dT .

(3.31)

Observar que esta ecuación tiene una u ńica soluci´ on ξ debido a la condici´ on de rango máximo (3.27). ´ Unicamente falta por demostrar que la matriz T es, en efecto, no singular. Lo probaremos mostrando que sus filas son linealmente independientes. Realizamos por lo tanto la combinaci´ on lineal de las filas de T siguiente α1 ξ + α2 ξA + · · · + αn ξAn−1 = 0 , on lineal por la derecha con ciertos escalares αi . Multiplicando esta combinaci´ aloga, es decir, multiplipor b y usando (3.30) se tiene que αn = 0. De forma an´ cando por la derecha por Ab, luego A2 b y as´ı sucesivamente y usando (3.30) se llega a que αn−1 = · · · = α1 = 0. En resumen, T es no singular. Demostremos ahora la necesidad, es decir, supongamos que existe la matriz T no singular tal que el sistema (3.28) coincida con el sistema (3.25). En particular se tiene T b = d y T AT −1 = C. Definamos r = rang[b, Ab, A2 b, . . . , An−1 b]. Queremos llegar a ver que r = n. Puesto que T es no singular, es claro que r = rang[T b, T Ab, T A2 b, . . . , T An−1 b]. Reordenando se obtiene que r = rang[T b, (T AT −1 )T b, . . . , (T AT −1 )n−1 T b]. Pero es claro que entonces r = rang[d, Cd, . . . , C n−1 d]. Teniendo en cuenta que on, es fácil C tiene la forma dada en (3.26) y que d = (0, . . . , 0, 1)T por definici´ as comprobar que la matriz U = [d, Cd, . . . , C n−1 d] tiene forma triangular. M´ concretamente, si denotamos por uij a los elementos de la matriz U entonces 1 si i = j , uij = 0 si i > j . En particular el rango r de U es máximo, o sea r = n.

16


Nota 7. Observar que la demostraci´ on del Teorema 6 es constructiva en el sentido de que T se puede construir mediante (3.29) donde ξ es la soluci´ on de (3.31).

3.4.

Problemas resueltos

1. Hallar la soluci´ on general del sistema lineal x˙ 0 1 x = y˙ −2 −3 y utilizando diferentes métodos: (i) Mediante la forma espectral de la soluci´ on. (ii) Calculando la matriz exponencial mediante la f´ ormula de Sylvester. (iii) Calculando la matriz exponencial como un polinomio matricial adecuado. Soluci´ on. (i) La forma espectral de la solución general del sistema lineal z˙ = Az viene dado por z(t) =

n

(vi z(0)) exp(λi t)wi ,

i=1

siendo λi y wi los valores y vectores propios asociados a la matriz A. Además, si W = col{w1 , . . . , wn }, entonces W −1 = fil{v1 , . . . , vn }. En nuestro caso se tiene

A=

0 −2

1 −3

,

de modo que, el polinomio caracter´ıstico k(λ) asociado a A es λ −1 = λ2 − 3λ + 2 = (λ + 1)(λ + 2) . k(A) = det(λI2 − A) = 2 λ−3 Entonces, los valores propios (ra´ıces del polinomio caracter´ıstico) de A son λ1 = −1 , λ2 = −2 . Los vectores propios wi verifican Awi = λi wi o bien, de forma equivalente (A − λi I2 )wi = 0. Entonces se tiene w1 = (α1 , β1 ) verifica el sistema lineal homogéneo 0 1 1 α1 = , β1 0 −2 −2 cuya solución es α1 + β1 = 0, de modo que, w1 = (1, −1)T .

3.4 Problemas resueltos

17

Sea w2 = (α2 , β2 ). Entonces se verifica el sistema lineal homogéneo 0 2 1 α2 = , β2 0 −2 −1 cuya solución es 2α2 + β2 = 0, de modo que, w2 = (1, −2)T . Se tiene pues que

W =

por lo tanto su matriz inversa es W −1 =

1 1 −1 −2

2 −1

1 −1

, ,

de lo que v1 = (2, 1) y v2 = (−1, −1). Finalmente, aplicando la f´ ormula de la solución en forma espectral se obtiene x(t) x(0) x(0) = v1 exp(λ1 t)w1 + v2 exp(λ2 t)w2 y(t) y(0) y(0) x(0) 1 = (2, 1) exp(−t) + y(0) −1 x(0) 1 (−1, −1) exp(−2t) y(0) −2 [2x(0) + y(0)] exp(−t) − [x(0) + y(0)] exp(−2t) = . −[2x(0) + y(0)] exp(−t) + 2[x(0) + y(0)] exp(−2t) (ii) Como todos los valores propios λi de A son distintos, se puede utilizar la f´ ormula de Sylvester (3.14) dada por exp(At) =

n

Zk exp(λk t) ,

k=1

siendo

n A − λj In , Zk = λk − λj j=1 j=k

matrices constantes que sólo dependen de A y de sus valores propios λj . En nuestro problema se tiene exp(At) = Z1 exp(−t) + Z2 exp(−2t) , donde Z1

=

Z2

=

A + 2I2 2 1 = , −2 −1 −1 + 2 A + I2 −1 −1 = . 2 2 −2 + 1

18

Soluci´ on de Sistemas Lineales De este modo, la solución general del sistema lineal es x(t) x(0) = exp(At) , y(t) y(0) siendo la matriz exponencial de A 2 1 −1 −1 exp(At) = exp(−t) + exp(−2t) . −2 −1 2 2 (iii) Otro método alternativo para calcular la matriz exponencial de A cuando todos sus valores propios λi son diferentes viene dado (3.16), es decir, puesto que A ∈ M2 (R) se tiene exp(At) = r(A) con r(λ) un polinomio de grado menor o igual que 1 cuyos coeficientes son funciones del tiempo t obtenidos de la forma r(λi ) = exp(λi t) , i = 1, 2 . Sea r(λ) = r0 (t) + r1 (t)λ. Entonces, la anterior condici´ on se escribe como r0 (t) − r1 (t) = exp(−t) , r0 (t) − 2r1 (t) = exp(−2t) , cuya solución es r0 (t) = 2 exp(−t) − exp(−2t) , r1 (t) = exp(−t) − exp(−2t) . En conclusi´ on, la soluci´ on general del sistema lineal es x(t) x(0) = exp(At) , y(t) y(0) siendo la matriz exponencial de A exp(At) = =

r0 (t)I2 + r1 (t)A

1 0 0 1

[2 exp(−t) − exp(−2t)] + 0 1 [exp(−t) − exp(−2t)] . −2 −3

2. Demostrar que L{exp(At)} = (sIn − A)−1 , siendo A ∈ Mn (R) y L{.} la transformada de Laplace. Soluci´ on. Consideremos el problema de Cauchy x˙ = Ax, x(0) = x0 . Tomando transformadas de Laplace L{x} ˙ = L{Ax} y teniendo en cuenta la condici´ on inicial x(0) = x0 se obtiene s¯ x(s) − x0 = A¯ x(s) ,


19

de donde, despejando tenemos x ¯(s) = (sIn − A)−1 x0 . Antitransformando en esta ecuaci´ on se obtiene x(t) = L−1 {(sIn − A)−1 x0 } Comparando esta expresión con la soluci´ on x(t) = exp(At)x0 del problema de Cauchy x˙ = Ax, x(0) = x0 se concluye que L{exp(At)} = (sIn − A)−1 tal y como se quer´ıa probar. 3. Hallar la soluci´ on del problema de Cauchy x˙ 1 0 1 x1 0 x1 (0) 0 = + , = . x˙ 2 −1 0 x2 t 0 x2 (0) Soluci´ on. El sistema del enunciado se puede escribir como x˙ = Ax + Bu, siendo x = (x1 , x2 )T ∈ R2 , 0 1 A= , −1 0 y donde el término Bu se puede tomar de varias formas. Por ejemplo 0 B= ∈ R2 , u(t) = t ∈ R . 1 Teniendo en cuenta (3.22), la soluci´ on del problema de Cauchy x˙ = Ax + Bu con x(0) = x0 es

t x(t) = Φ(t, 0) x0 + Φ(0, τ )Bu(τ ) dτ , 0

siendo Φ(t, 0) = exp(At) la matriz de transici´ on de estados. Es fácil comprobar que dicha matriz es una matriz de rotaci´ on en el plano cos t sin t Φ(t, 0) = , − sin t cos t siendo su inversa Φ(0, t) =

cos t − sin t

sin t cos t

−1 =

cos t sin t

Se tiene pues que t cos τ cos t sin t 0 x(t) = + sin τ − sin t cos t 0 0

− sin t cos t

− sin τ cos τ

.

0 τ

dτ

de modo que, tras finalizar todos los c´ alculos involucrados se obtiene t − sin t x(t) = . 1 − cos t

,

20

Soluci´ on de Sistemas Lineales 4. Sea u(t) la fuerza total por unidad de masa que act´ ua sobre una part´ıcula. Demostrar que la posici´ on en funci´ on del tiempo x(t) viene dada por

t

x(t) = x(0) + x(0)t ˙ + 0

(t − τ )u(τ ) dτ .

Soluci´ on. La ecuación del movimiento es, seg´ un las leyes de la mecánica clásica, x ¨ = u(t). Pasando a las variables de estado o variables de fase y = x˙ se tiene x˙ = y , y˙ = u(t) . Podemos escribir este sistema en forma matricial como x˙ 0 1 x 0 = + u(t) , y˙ 0 0 y 1 es decir, definiendo z = (x, y)T , se tiene z˙ = Az + bu, siendo 0 1 0 A= , b= . 0 0 1 La solución z(t) viene expresada, seg´ un (3.21), por

t z(t) = exp(At) z0 + exp(−Aτ )bu(τ ) dτ ,

(3.32)

0

siendo z(0) = z0 . Realizamos ahora un paréntesis para calcular la matriz exp(−Aτ ). Un simple cálculo por cualquiera de los métodos explicados en teor´ıa muestra que 1 t . (3.33) exp(At) = I2 + At = 0 1 En realidad, en este caso particular es f´ acil ver que A2 = O, es decir, la matriz A es nilpotente. Por lo tanto, la expresi´ on dada de exp(At) se sigue de la propia definici´ on de matriz exponencial (3.5) donde la serie es truncada y se convierte en polinomial, es decir, exp(At) =

∞ k t k=0

k!

Ak = I2 + tA +

t2 2 A + · · · = I2 + tA . 2!

Introduciendo (3.33) en (3.32) y teniendo en cuenta las expresiones de A y b se tiene t 1 −τ 0 x(t) 1 t x(0) u(τ ) dτ . z(t) = = + 0 1 1 y(t) 0 1 y(0) 0


21

Finalmente, extrayendo de esta ecuación la primera componente x(t) y recordando que y = x˙ se obtiene el resultado deseado

t

x(t) = x(0) + x(0)t ˙ + 0

(t − τ )u(τ ) dτ .

5. Averiguar si el sistema lineal 1 −3 x1 1 x˙ 1 = + u(t) , x˙ 2 4 2 x2 1 puede ser transformado en la forma can´ onica controlable mediante un cambio lineal de coordenadas. Dar, si es posible, dicho cambio y la forma can´ onica controlable asociada al sistema. Soluci´ on. El sistema del enunciado se puede escribir en forma compacta como x˙ = Ax + bu con las definiciones usuales, es decir, 1 −3 1 A= , b= . 4 2 1 onica controlable Puesto que A ∈ M2 (R), recordemos que, la forma can´ (3.25) dada por w˙ = Cw + du con w ∈ R2 , d = (0, 1)T y C ∈ M2 (R) la matriz “companion”, no es más que la expresión de un sistema lineal que provenga mediante el cambio usual de una ecuaci´ on diferencial escalar de segundo orden a coeficientes constantes (3.23), es decir, z¨ + k1 z˙ + k2 z = u(t). Sabemos, por el Teorema 6, que existirá un cambio de variables lineal w = T x con T ∈ M2 (R) no singular que transformar´ a el sistema dado en el enunciado a la forma can´ onica controlable si y s´ olo si rang[b, Ab] = 2. En nuestro problema se tiene 1 −2 [b, Ab] = , 1 6 de modo que se verifica rang[b, Ab] = 2. Una vez se sabe que existe la transformación lineal w = T x, sólo falta hallar la matriz T ∈ M2 (R). Para ello recordemos la Nota 7, es decir, T se puede construir mediante (3.29) donde ξ es la solución de (3.31). En otras palabras,

ξ T = ∈ M2 (R) , ξA donde ξ verifica ξ[b, Ab] = dT .

22

Soluci´ on de Sistemas Lineales En nuestro problema, despejando ξ de la anterior ecuaci´ on se tiene ξ = dT [b, Ab]−1 = (0, 1)

1 −2 1 6

−1 =

1 (0, 1) 8

6 −1

2 1

=

1 (−1, 1) . 8

Finalmente, se tiene

T =

ξ ξA

=

1 8

−1 1 3 5

.

Sabemos que, ver (3.28), realizando el cambio de variables w = T x en el sistema x˙ = Ax + bu se obtiene el sistema transformado w˙ = T AT −1 w + T bu , con C = T AT −1 y T b = d = (0, 1)T . Desarrollando se tiene 0 1 w1 0 w˙ 1 = + u(t) , w˙ 2 w2 −14 3 1 que, obviamente, corresponde a una ecuaci´ on escalar de segundo orden del control clásico. Más concretamente, definiendo z = w1 = 1/8(−x1 + x2 ) se concluye que z¨ − 3z˙ + 14z = u(t) .

Cap´ıtulo 4

Sistemas de Control Lineal 4.1.

Introducci´ on

En este cap´ıtulo abordaremos algunos de los problemas propios del control. Consideremos el sistema lineal controlado x˙ = Ax + Bu ,

(4.1) n

con A ∈ Mn (R) y B ∈ Mn×m (R) matrices constantes, x(t) ∈ R y u(t) ∈ Rm .

4.2.

Controlabilidad

Definici´ on 8. El sistema (4.1) se dice que es completamente controlable si, para cualquier t0 , cualquier estado inicial x(t0 ) = x0 y cualquier estado final xf , existe un tiempo finito t1 con t1 > t0 y un control u(t) definido para t0 ≤ t ≤ t1 tal que x(t1 ) = xf . El término “completamente”significa que la anterior definici´ on se satisface “para todo”x0 y xf . Nota 9. Es habitual suponer que el control u(t) es cotinuo a trozos, es decir, continuo excepto en un n´ umero finito de tiempos del intervalo [t0 , t1 ]. Nota 10. En la Definici´ on 8 se puede tomar el estado final xf = 0 sin pérdida de generalidad. Adem´ as, puesto que A y B son matrices constantes (no dependen del tiempo) también se puede elegir el origen de tiempos de modo que t0 = 0. La Nota 10 es debida a que, la soluci´ on x(t) de (4.1) viene dada por (3.22), es decir,

t Φ(t0 , τ )Bu(τ ) dτ . x(t) = Φ(t, t0 ) x0 + t0

De este modo se tiene

xf = Φ(t1 , t0 ) x0 +

t1

t0

23

Φ(t0 , τ )Bu(τ ) dτ

,

24

Sistemas de Control Lineal

que, reagrupando términos, se puede escribir como

t1 Φ(t0 , τ )Bu(τ ) dτ , 0 = Φ(t1 , t0 ) {x0 − Φ(t0 , t1 )xf } + t0

puesto que la matriz Φ(t1 , t0 ) es no singular, siendo su inversa Φ(t0 , t1 ). Por lo tanto se observa que, si el control u(t) transfiere el estado inicial x0 en el estado final xf , entonces tambien transfiere el estado x0 − Φ(t0 , t1 )xf al origen en el mismo tiempo t1 − t0 . El siguiente resultado muestra un criterio algebraico que caracteriza a los sistemas lineales (4.1) completamente controlables. Teorema 11. El sistema x˙ = Ax + Bu dado en (4.1) es completamente controlable si y s´ olo si la matriz de controlabilidad U = [B, AB, A2 B, . . . , An−1 B] ∈ Mn×nm (R)

(4.2)

tiene rango m´ aximo, es decir, rangU = n. Demostraci´ on. Probemos en primer lugar la necesidad. Suponemos pues que el sistema (4.1) es completamente controlable (asumiendo t0 = 0 y xf = 0 debido a la Nota 10) y queremos ver entonces que rangU = n. Utilizaremos una demostración por reducci´ on al absurdo, es decir, supondremos que rangU < n y llegaremos a una contradicci´ on. Si rangU < n entonces existe una combinaci´ on lineal de filas de U que da el vector fila nulo. Dicho de otro modo, existe un vector fila q ∈ Rn con q = 0 tal que (4.3) qB = 0 , qAB = 0 , . . . , qAn−1 B = 0 . Utilizando la soluci´ on (3.21) del sistema (4.1) con las condiciones x(0) = x0 se tiene

t exp(−Aτ )Bu(τ ) dτ , x(t) = exp(At) x0 + 0

de modo que, tomando t = t1 con x(t1 ) = 0 se obtiene t1 −x0 = exp(−Aτ )Bu(τ ) dτ . 0

Teniendo en cuenta que exp(−At) = r(A) con r un polinomio de grado menor o igual que n − 1, la anterior ecuaci´ on queda de la forma t1 (r0 (τ )In + r1 (τ )A + · · · + rn−1 (τ )An−1 )Bu(τ ) dτ , −x0 = 0

on siendo ri (t) ∈ R los coeficientes del polinomio r. Multiplicando esta ecuaci´ por la izquierda por el vector fila q y utilizando (4.3) se llega a que qx0 = 0. Como el sistema (4.1) es completamente controlable, el desarrollo efectuado anon qx0 = 0 implica q = 0 teriormente sirve para cualquier x0 . Entonces la condici´

4.2 Controlabilidad

25

en contradicci´ on con la hip´ otesis inicial rangU < n. Probemos en segundo lugar la suficiencia. Para ello supondremos que rangU = n y hemos de demostrar que el sistema (4.1) es completamente controlable, es decir, que existe una funci´ on de control u(t) definida en el intervalo de tiempo 0 ≤ t ≤ t1 tal que x(t1 ) = 0. Consideremos la siguiente matriz cuadrada t1 M= exp(−Aτ )BB T exp(−AT τ ) dτ ∈ Mn (R) . (4.4) 0

Es obvio que M es una matriz constante. Además M es simétrica puesto que

t1 T T T T M = exp(−Aτ )BB exp(−A τ ) dτ

0 t1

= 0

t1

= 0

=

t1

T exp(−Aτ )BB T exp(−AT τ ) dτ

[exp(−AT τ )]T BB T [exp(−Aτ )]T dτ exp(−Aτ )BB T exp(−AT τ ) dτ = M .

0

Vamos a probar que M es no singular demostrando que es, en realidad, una matriz definida positiva. Sea α ∈ Rn un vector columna arbitrario y definamos el atica vector fila β(τ ) = αT exp(−Aτ )B. Asociamos a M la siguiente forma cuadr´ t1 t1 β(τ )β T (τ ) dτ = β(τ )2 dτ ≥ 0 , (4.5) αT M α = 0

0

de modo que M es semidefinida positiva. Vamos a ver a continuaci´ on que M es no singular por reducci´ on al absurdo. Se sabe que M es singular si y sólo si ˜ = 0. Esto implica, teniendo en existe un vector columna, α ˜ = 0 tal que α ˜T M α T ˜ cuenta (4.5), que β(τ ) = α ˜ exp(−Aτ )B ≡ 0 para todo 0 ≤ τ ≤ t1 . Utilizando la definici´ on de matriz exponencial, la anterior condici´ on se escribe como 2 τ α ˜ T In − τ A + A2 + · · · B ≡ 0 , 0 ≤ τ ≤ t1 . 2! Ello es equivalente a imponer que ˜ T AB = 0 , α ˜ T A2 B = 0, . . . , α ˜T B = 0 , α de manera que α ˜ T U = 0 siendo U la matriz de controlabilidad (4.2). Pero, como por hip´ otesis U tiene rango máximo, es decir rango n, es imposible que existe un vector α ˜ = 0 verificando α ˜ T U = 0. Hemos llegado pues a una contradicci´ on y por lo tanto M es no singular. Entonces es posible tomar el siguiente control u(τ ) = −B T exp(−AT τ )M −1 x0 , 0 ≤ τ ≤ t1 .

(4.6)

26


Sustituyendo este control en la soluci´ on (3.21) dada por

t x(t) = exp(At) x0 + exp(−Aτ )Bu(τ ) dτ , 0

del sistema (4.1) se tiene x(t1 )

t1 exp(At1 ) x0 − exp(−Aτ )BB T exp(−AT τ )M −1 x0 dτ 0 = exp(At1 ) x0 − M M −1 x0 = 0 ,

=

donde en el segundo paso se ha utilizado (4.4). Nota 12. Debido al Teorema 11, es habitual hablar de la controlabilidad del par [A, B] en lugar del sistema x˙ = Ax + Bu. Corolario 13. Si rangB = r entonces la condici´ on (4.2) del Teorema 11 se reduce a la siguiente restricci´ on rang[B, AB, A2 B, . . . , An−r B] = n .

Demostraci´ on. Definimos las siguientes matrices Uk = [B, AB, A2 B, . . . , Ak B] , k = 0, 1, 2, . . . Observar que, si rangU = rangU+1 , entonces todas las columnas de la matriz A+1 B son linealmente dependientes con respecto a las columnas de U . Además, en este caso, se verificará también que todoas las columnas de A+2 B, A+3 B, . . ., son también linealmente dependientes con respecto a las columnas de U . Se tiene pues que rangU = rangU+1 = rangU+2 = · · · . De este modo se ve que el rango de Uk aumenta al menos en 1 cuando ¿k aumenta en 1 hasta que el máximo rango de Uk se alcanza para k = . Como rangU0 = rangB = r y además rangUk ≤ n se tiene que r + ≤ n y por lo tanto ≤ n − r. El siguiente resultado relaciona el sistema de control cl´ asico (3.25) con el concepto de controlabilidad. Teorema 14. El sistema lineal x˙ = Ax + bu con A ∈ Mn (R) matriz constante, b ∈ Rn vector columna y u(t) ∈ R puede ser transformado mediante un cambio lineal en la forma can´ onica controlable (3.25) si y s´ olo si es completamente controlable. Demostraci´ on. Como b ∈ Rn es un vector columna, es decir, estamos en el caso m = 1, la matriz B reduce al vector b. Entonces, la condici´ on (4.2) del Teorema 11 es exactamente igual a la condición (3.27) del Teorema 6. El teorema se demustra pues teniendo en cuenta el Teorema 11 y el Teorema 6.

4.2 Controlabilidad

27

Nota 15. Debido al Teorema 14, el sistema cl´ asico (3.25) se suele decir que est´ a en la forma canónica controlable. Observemos que el Teorema 11 es puramente existencial, es decir, a partir de él se sabe que existe un vector control u(t) de modo que el sistema es completamente controlable. Sin embargo no determina cual es dicho control. Para responder a esta pregunta se dispone del siguiente resultado. Teorema 16. El sistema x˙ = Ax + Bu es completamente controlable si y s´ olo si la matriz de controlabilidad simétrica t1 Φ(t0 , τ )BB T ΦT (t0 , τ ) dτ ∈ Mn (R) (4.7) U (t0 , t1 ) = t0

es no singular, siendo Φ la matriz de transici´ on de estados definida en (3.8). Adem´ as, en este caso el control u(t) definido para t0 ≤ t ≤ t1 que transfiere el estado x(t0 ) = x0 al estado x(t1 ) = xf viene dado por u(t) = −B T ΦT (t0 , t)U −1 (t0 , t1 )[x0 − Φ(t0 , t1 )xf ] .

(4.8)

Demostraci´ on. Demostremos en primer lugar la suficiencia. Asumimos pues que U es no singular, de modo que el control dado por (4.8) existe. S´ olo falta por ver que, aplicando dicho control, se obtiene x(t1 ) = xf . Para ello, recordemos que la soluci´ on (3.22) del sistema x˙ = Ax+Bu con la condici´ on inicial x(t0 ) = x0 viene dada por

t x(t) = Φ(t, t0 ) x0 + Φ(t0 , τ )Bu(τ ) dτ . (4.9) t0

Se tiene pues, sustituyendo el control u(t) dado en (4.8), que t T T −1 x(t1 )

=

Φ(t1 , t0 ) x0 +

=

Φ(t1 , t0 ) x0 −

=

xf ,

t0

Φ(t0 , τ )B −B Φ (t0 , τ )U t1

t0

Φ(t0 , τ )BB T Φ(t0 , τ ) dτ

(t0 , t1 )[x0 − Φ(t0 , t1 )xf ]

U −1 (t0 , t1 )(x0 − Φ(t0 , t1 )xf )

como se quer´ıa ver. Hemos tenido en cuenta que Φ(t1 , t0 )Φ(t0 , t1 ) = In y que adem´ as, la expresión encerrada entre llaves {.} es justo la definici´ on de U (t0 , t1 ). Demostremos ahora la necesidad. Hemos de probar que, si el sistema x˙ = Ax + Bu es completamente controlable, entonces U (t0 , t1 ) es no singular. Como U (t0 , t1 ) es una matriz simétrica, tomando un vector columna arbiatica trario α ∈ Rn podemos construir la forma cuadr´ t1 t1 αT U (t0 , t1 )α = β T (τ, t0 )β(τ, t0 ) dτ = β(τ, t0 )2 dτ ≥ 0 , (4.10) t0

t0

dτ

28


donde se ha definido el vector columna β(τ, t0 ) = B T ΦT (t0 , τ )α ∈ Rn . Se tiene pues que U (t0 , t1 ) es una matriz semidefinida positiva. Veamos por reducción al absurdo que, en realidad, U (t0 , t1 ) es una matriz definida positiva y por lo tanto no singular. Supongamos pues que entonces que α = 0. Teniendo en cuenta (4.10) con existe un vector α ˆ = 0 tal que α ˆ T U (t0 , t1 )ˆ α=α ˆ se tiene que t1 ˆ t0 )2 dτ = 0 β(τ, t0

ˆ t0 ) = B T ΦT (t0 , τ )ˆ donde β(τ, α ∈ Rn . Por lo tanto, por las propiedades de la ˆ t0 ) ≡ 0 para t0 ≤ τ ≤ t1 . norma, concluimos que β(τ, Sin embargo, puesto que por hip´ otesis el sistema x˙ = Ax + Bu es completamente controlable, se sabe que existe un vector control u ˆ(t) ∈ Rn de modo que ˆ Entonces, utilizando (4.9) con la condici´ on inicial x(t1 ) = 0 cuando x(t0 ) = α. ˆ y la final x(t1 ) = 0 se tiene x(t0 ) = α α ˆ=−

t1

t0

Φ(t0 , τ )B u ˆ(τ ) dτ .

De este modo 2

T

ˆ α ˆ=− ˆ α = α

t1

T

T

T

u ˆ (τ )B Φ(t0 , τ ) α ˆ dτ = −

t0

t1

t0

ˆ t0 ) dτ = 0 , u ˆT (τ )β(τ,

contradiciendo la hip´ otesis α ˆ = 0. Es importante observar que pueden haber, en general, otros controles u(t) diferentes del control (4.8) que transfieran el estado x(t0 ) = x0 al estado x(t1 ) = on que el control (4.8) tiene una propiedad xf . Sin embargo veremos a continuaci´ especial que le confiere optimabilidad respecto del resto de posibles controles. Teorema 17. Consideremos el sistema x˙ = Ax+Bu. Sea u ˆ(t) un vector control que lleva el estado inicial x(t0 ) = x0 al estado final x(t1 ) = xf verificando u ˆ(t) ≡ u(t) para t0 ≤ t ≤ t1 , donde u(τ ) es el control (4.8). Entonces

t1

t0

2

ˆ u(τ ) dτ >

t1

t0

u(τ )2 dτ .

Demostraci´ on. Puesto que ambos controles u(t) y u ˆ(t) transfieren el estado un (3.22), inicial x(t0 ) = x0 al estado final x(t1 ) = xf , deben verificar, seg´

Φ(t1 , t0 ) x0 +

xf

=

xf

= Φ(t1 , t0 ) x0 +

t1

t0 t1 t0

Φ(t0 , τ )B u ˆ(τ ) dτ

,

Φ(t0 , τ )Bu(τ ) dτ

.

4.2 Controlabilidad

29

Restando estas expresiones se obtiene t1 0= Φ(t0 , τ )B[ˆ u(τ ) − u(τ )] dτ . t0

Multiplicando esta ecuaci´ on por la izquierda por [x0 − Φ(t0 , t1 )xf ]T [U −1 (t0 , t1 )]T , se llega a 0=

t1

t0

[x0 − Φ(t0 , t1 )xf ]T [U −1 (t0 , t1 )]T Φ(t0 , τ )B [ˆ u(τ ) − u(τ )] dτ .

Teniendo en cuenta que u(τ ) es el control (4.8), es claro que la expresión dentro de las llaves {.} es justo −uT (τ ) de modo que t1 0= uT (τ )[ˆ u(τ ) − u(τ )] dτ . (4.11) t0

Por otra parte, utilizando las propiedades del producto escalar y de la norma ˆ) = u2 + ˆ u2 − 2uT u ˆ. Entonces, se tiene que (u − u ˆ)T (u − u t1 t1 u(τ ) − u ˆ(τ )2 dτ = [(u(τ ) − u ˆ(τ ))T (u(τ ) − u ˆ(τ ))T ] dτ t0

t0 t1

=

t0 t1

=

t0

[u(τ )2 + ˆ u(τ )2 − 2uT (τ )ˆ u(τ )] dτ [ˆ u(τ )2 − u(τ )2 ] dτ ,

donde, en el u ´ltimo paso, se ha utilizado (4.11). De aqu´ı, finalmente obtenemos t1 t1 t1 2 2 2 ˆ u(τ ) dτ = [u(τ ) + u(τ ) − u ˆ(τ ) ] dτ > u(τ )2 dτ t0

t0

t0

como se quer´ıa probar. Observar que en la u ´ltima desigualdad se ha utilizado que u(t) ≡ u ˆ(t) en el intervalo t0 ≤ t ≤ t1 . Nota 18. Dado un control u(t), en muchas aplicaciones f´ısicas se interpreta la t integral t01 u(τ )2 dτ como la energ´ıa o bien el consumo necesario para llevar el sistema desde el estado inicial x(t0 ) = x0 al estado final x(t1 ) = xf . De este modo, el Teorema 17 se interpreta de manera que, de todos los posibles controles u ˆ(t) admisibles, el que minimiza la energ´ıa o el consumo es justo el control u(t) dado por (4.8). Si el sistema lineal x˙ = Ax + Bu no es completamente controlable, no es una buena terminolog´ıa llamarle incontrolable. Ello es bebido a que cuando no se verifica la Definición 8 de sistema completamente controlable sólo significa

30


que hay ciertos estados finales xf que no pueden ser alcanzados con ning´ un un control, control. Para caracterizar si un estado final xf es alcanzable con alg´ independientemente de si el sistema es completamente controlable o no y además conocer qué control u(t) es admisible se tiene el siguiente resultado. Teorema 19. Si, para un cierto estado final xf , existe un vector columna constante ξ ∈ Rn tal que

entonces el control

U (t0 , t1 )ξ = x0 − Φ(t0 , t1 )xf ,

(4.12)

u(t) = −B T Φ(t0 , t)ξ ,

(4.13)

transfiere el sistema x˙ = Ax + Bu desde el estado inicial x(t0 ) = x0 al estado final x(t1 ) = xf . Demostraci´ on. La solución del sistema x˙ = Ax + Bu con la condici´ on inicial x(t0 ) = x0 viene dada por (3.22), es decir,

t Φ(t0 , τ )Bu(τ ) dτ . x(t) = Φ(t, t0 ) x0 + t0

Sustitutendo en esta expresi´ on el control u(t) dado en (4.13) y evaluando para t = t1 se obtiene

t1 x(t1 ) = Φ(t1 , t0 ) x0 − Φ(t0 , τ )BB T Φ(t0 , τ )ξ dτ . t0

Teniendo en cuenta, seg´ un (4.7), la expresi´ on de la matriz de controlabilidad U (t0 , t1 ), la anterior ecuaci´ on se reescribe como x(t1 ) = Φ(t1 , t0 ) [x0 − U (t0 , t1 )ξ] . Imponiendo ahora la condici´ on (4.12) se obtiene x(t1 ) = Φ(t1 , t0 )Φ(t0 , t1 )xf = xf , como se quer´ıa probar. Nota 20. Si el sistema x˙ = Ax+Bu es completamente controlable, entonces, por on el Teorema 4.8, la matriz de controlabilidad U (t0 , t1 ) es invertible y la expresi´ dada en (4.13) para el control u(t) reduce a la expresi´ on (4.8). Nota 21. Se puede demostrar, ver por ejemplo la p´ agina 77 de [2], que el rec´ıproco del Teorema 19 es cierto. De este modo, s´ olo los estados finales xf para los cuales (4.12) se verifique pueden ser alcanzados con alg´ un control. Existe otra forma an´ aloga al Teorema 19 demostrada en la página 167 de [11] que es la siguiente. Teorema 22. Un cierto estado inicial x(t0 ) = x0 se puede transferir al estado final xf si ambos vectores x0 y xf pertenecen al subespacio vectorial engendrado por los vectores columnas de la matriz de controlabilidad U dada en (4.2).

4.3 Equivalencia Algebraica

4.3.

31

Equivalencia Algebraica

Definici´ on 23. Sea P ∈ Mn (R) una matriz no singular. El sistema lineal ˆx + Bu ˆ obtenido a partir del sistema lineal x˙ = Ax + Bu mediante el x ˆ˙ = Aˆ cambio de variables x ˆ(t) = P x(t) se llama sistema algebraicamente equivalente. ˆx + Bu ˆ son algebraicaProposici´ on 24. Si los sistemas x˙ = Ax + Bu y x ˆ˙ = Aˆ mente equivalentes mediante el cambio de variables x ˆ = P x entonces ˆ = PB . Aˆ = P AP −1 , B

(4.14)

Demostraci´ on. Derivando respecto de t la identidad x ˆ(t) = P (t)x(t) se obtiene x ˆ˙ (t) = P˙ x + P x˙ = P x˙ = P (Ax + Bu) = P Ax + P Bu = P AP −1 x ˆ + P Bu . ˆx + Bu ˆ se obtiene el Identificando los anteriores términos con los de x ˆ˙ = Aˆ resultado. ˆ t0 ) las matrices de transici´ Proposici´ on 25. Sean Φ(t, t0 ) y Φ(t, on de estados ˆx + Bu ˆ de los sistemas algebraicamente equivalente x˙ = Ax + Bu y x ˆ˙ = Aˆ respectivamente mediante el cambio de variables x ˆ(t) = P x(t). Entonces ˆ t0 ) = P Φ(t, t0 )P −1 . Φ(t,

(4.15)

Demostraci´ on. Recordemos que Φ(t, t0 ) es la u ´ nica matriz no singular que satisface ˙ t0 ) = AΦ(t, t0 ) , Φ(t0 , t0 ) = In . Φ(t, ˆ t0 ) es la matriz de transición de estados del sistema Para demostrar que Φ(t, ˙x ˆ ˆ ˆ = Aˆ x + Bu se ha de probar que ˆ˙ t0 ) = AˆΦ(t, ˆ t0 ) , Φ(t ˆ 0 , t0 ) = In . Φ(t, Efectivamente, ˙ t0 )P −1 = P AΦ(t, t0 )P −1 = (P AP −1 )(P Φ(t, t0 )P −1 ) = AˆΦ(t, ˆ t0 ) , ˆ˙ t0 ) = P Φ(t, Φ(t, donde en el primer paso se ha utilizado (4.15) y en el u ´ltimo (4.14). Adem´ as ˆ 0 , t0 ) = P Φ(t0 , t0 )P −1 = P P −1 = In , Φ(t con lo cual queda demostrada la proposici´ on. Teorema 26. En un sistema lineal, la propiedadad de ser completamente controlable se preserva por equivalencia algebraica.

32


Demostraci´ on. Consideremos dos sistema lineales algebraicamente equivaˆx + Bu. ˆ Hemos de demostrar que si x˙ = Ax + Bu es lente x˙ = Ax + Bu y x ˆ˙ = Aˆ ˆx + Bu ˆ también es completamente completamente controlable entonces x ˆ˙ = Aˆ controlable. ˆx + Bu ˆ viene ˆ (t0 , t1 ) para el sistema x ˆ˙ = Aˆ La matriz de controlabilidad U dada, seg´ un (4.7), por t1 ˆ 0 , τ )B ˆ T (t0 , τ ) dτ . ˆB ˆT Φ ˆ Φ(t U (t0 , t1 ) = t0

ˆ = P B y Φ(t, ˆ t0 ) = P Φ(t, t0 )P −1 , de modo Teniendo en cuenta (4.14) y (4.15), B que t1 ˆ (t0 , t1 ) = U (P Φ(τ, t0 )P −1 )(P B)(B T P T )((P −1 )T ΦT (τ, t0 )P T ) dτ t0 t1

=

t0

P Φ(τ, t0 )BB T ΦT (τ, t0 )P T dτ

=

t1

P t0

T

T

Φ(τ, t0 )BB Φ (τ, t0 ) dτ

PT .

Se verifica pues la siguiente ecuación ˆ (t0 , t1 ) = P U (t0 , t1 )P T . U

(4.16)

ˆ (t0 , t1 ) es no singular puesto que, En particular, la matriz de controlabilidad U por hip´ otesis, P y U (t0 , t1 ) son no singulares. Utilizando el Teorema 4.8 se concluye con el resultado deseado.

4.4.

Problemas resueltos

1. Averiguar si el sistema x˙ = Ax + bu, ⎛ 0 1 0 1 ⎜ 0 0 −1 0 A=⎜ ⎝ 0 0 0 1 0 0 5 0

siendo x ∈ R4 , u ∈ R, ⎞ ⎛ ⎞ 0 ⎟ ⎜ ⎟ ⎟ , b=⎜ 1 ⎟ , ⎠ ⎝ 0 ⎠ −2

es completamente controlable. Soluci´ on. Puesto que A ∈ Mn (R), la matriz de controlabilidad U del sistema lineal del enunciado es ⎛ ⎞ 0 −1 0 −8 ⎜ 1 0 2 0 ⎟ ⎟ . U = [b, Ab, A2 b, A3 b] = ⎜ ⎝ 0 −2 0 −10 ⎠ −2 0 −10 0 Es f´ acil ver que rang U = 4, de modo que, aplicando el Teorema 11, se concluye que el sistema es completamente controlable.


33

2. Consideremos un péndulo invertido cuya ecuaci´ on del movimiento es u ¨ cos θ + θ¨ = g sin θ , siendo la longitud del péndulo, g la intensidad de campo gravitatorio, θ el ´ angulo que forma el péndulo respecto de la vertical y u la posici´ on horizontal de la base del péndulo que es la variable de control. (i) Demostrar que, si θ es inicialmente peque¯ no entonces el sistema es ´ n: Utilizar la variable x de completamente controlable. Indicacio posici´ on horizontal del péndulo y demostrar que, para peque¯ nas oscilaciones, se tiene x ¨ ≈ w2 (x − u), siendo w2 := g/. (ii) Calcular la expresi´ on expl´ıcita de u(t) si se pretende llevar el péndulo de peque¯ nas oscilaciones desde el estado inicial (x(0), x(0)) ˙ = (1, 1) hasta el estado final (x(1), x(1)) ˙ = (0, 0) en una unidad de tiempo. Soluci´ on. (i) Para θ peque¯ na se tiene cos θ ≈ 1, sin θ ≈ θ y además sin θ = (x − u)/, siendo x la posición horizontal del péndulo. Se tiene pues que la ecuación del movimiento es aproximadamente u ¨ + θ¨ = g sin θ .

(4.17)

Por otra parte x = sin θ + u , x˙ = θ˙ cos θ + u˙ ≈ θ˙ + u˙ , x ¨ ≈ θ¨ + u ¨. on (4.17) obtenemos el sistema linearizado Sustituyendo en la ecuaci´ x ¨ = w2 (x − u) , siendo w2 := g/. Esta ecuación se puede reescribir en la forma y˙ = Ay + bu, siendo y = (x, x) ˙ T, 0 1 0 A= , b = . w2 0 −w2 La matriz de controlabilidad U de este sistema lineal es 0 −w2 U = [b, Ab] = , −w2 0 de modo que rang U = 2 y, aplicando el Teorema 11, se concluye que el sistema es completamente controlable. (ii) Consideramos de nuevo el sistema y˙ = Ay + bu, siendo y = (x, x) ˙ T, 0 1 0 A= , b = . w2 0 −w2

34

Sistemas de Control Lineal Aplicaremos el Teorema 16 para calcular el control u(t). En primer lugar calculamos la matriz de controlabilidad 1 Φ(0, τ )bbT ΦT (0, τ ) dτ , U (0, 1) = 0

siendo la matriz de transici´ on de estados Φ(0, τ ) = exp(−Aτ ). Es f´ acil mostrar que cosh(τ w) − sinh(τ w)/w Φ(0, τ ) = −w sinh(τ w) cosh(τ w) De este modo se obtiene 1 w2 sinh2 (τ w) −w3 sinh(2τ w)/2 U (0, 1) = dτ −w3 sinh(2τ w)/2 w4 cosh2 (τ w) 0 −w(w − sinh(w) cosh(w))/2 −w2 sinh2 (w)/2 . = −w2 sinh2 (w)/2 w3 (w + sinh(w) cosh(w))/2 Como el sistema es totalmente controlable, está claro que la matriz U (0, 1) es no singular. De cualquier forma, se tiene que det U (0, 1) =

1 4 w [−1 − 2w2 + cosh(2w)] , 8

que es no nulo para todo w > 0. Seg´ un el Teorema 16, el control que llevar´ a el sistema y˙ = Ay + bu con y = (x, x), ˙ desde el estado inicial y(0) = y0 = (1, 1)T hasta el estado final y(1) = yf = (0, 0)T viene dado por u(t) = −bT ΦT (0, t)U −1 (0, 1)[y0 − Φ(0, 1)yf ] . Realizando los cálculos pertinentes se obtiene u(t) =

2{−w cosh[(−2 + t)w] + 3w cosh[tw] + sinh[(−2 + t)w] + (−1 + 2w2 ) sinh[tw]} . w + 2w3 − w cosh(2w)

3. Consideremos un satélite de masa unidad en una ´ orbita plana. Tomando coordenadas polares (r, θ), las ecuaciones del movimiento son 2r˙ θ˙ 1 c + uθ , r¨ = rθ˙2 − 2 + ur , θ¨ = − r r r siendo c un par´ ametro de la ´ orbita y ur y uθ las componentes radial y angular de la fuerza u ejercida por el motor del satélite que son tomadas como variables de control. Notar que si se para el motor, es decir, ur = orbita circular uniforme (ypor lo tanto uθ = 0, entonces es posible una ´ con r˙ = θ¨ = 0) de radio R y velocidad angular θ˙ = w = c/R3 . Si la


35

fuerza u del motor es peque¯ na se pueden tomar las siguientes coordenadas ˙ x3 = θ − wt y x4 = θ˙ − w. Entonces de perturbaci´ on: x1 = r − R, x2 = r, se tiene x˙ ≈ Ax + Bu, siendo x = (x1 , x2 , x3 , x4 )T , u = (ur , uθ )T , ⎞ ⎛ ⎞ ⎛ 0 0 0 1 0 0 ⎜ ⎜ 3w2 0 0 2wR ⎟ 0 ⎟ ⎟ , B=⎜ 1 ⎟ . A=⎜ ⎠ ⎝ ⎝ 0 0 0 1 0 0 ⎠ 0 −2w/R 0 0 0 1/R Demostrar que el sistema es completamente controlable. Soluci´ on. Utilizaremos el Teorema 4.2 para resolver la cuestión planteada. En nuestro caso, la matriz de controlabilidad es U

= =

[B, AB, A2 B, A3 B]

⎛

0 ⎜ 1 ⎜ ⎝ 0 0

0 0 0 1/R

1 0 0 −2w/R

0 2w 1/R 0

0 −w2 −2w/R 0

2w 0 0 −4w2 /R

−w2 0 0 2w3 /R

⎞

0 −2w3 ⎟ ⎟. −4w2 /R ⎠ 0

Es f´ acil comprobar que rang U = 4, de modo que, aplicando el Teorema 4.2 se concluye que el sistema es completamente controlable.

Cap´ıtulo 5

Observabilidad Consideremos el sistema lineal controlado x˙ = Ax + Bu ,

(5.1)

con A ∈ Mn (R) y B ∈ Mn×m (R) matrices constantes, x(t) ∈ Rn y u(t) ∈ Rm . Adem´ as, supondremos que la salida y(t) ∈ Rr del sistema es y = Cx ,

(5.2)

con C ∈ Mr×n (R).

5.1.

Observabilidad

La idea de observabilidad para un sistema controlable consiste en determinar, si es posible, el estado x(t) de un sistema conociendo u ńicamente la salida y(t). De forma más precisa se tiene la siguiente definición. Definici´ on 27. El sistema (5.1) y (5.2) se dice que es completamente observable si, para cualquier t0 y cualquier estado inicial x(t0 ) = x0 , existe un tiempo finito t1 con t1 > t0 tal que el conocimiento del control u(t) y la salida y(t) para ńica. t0 ≤ t ≤ t1 es suficiente para determinar x0 de forma u Nota 28. En la Definici´ on 27 se puede tomar el control u(t) ≡ 0 para t0 ≤ t ≤ t1 sin pérdida de generalidad. Teorema 29. El sistema (5.1) y (5.2) definido por x˙ = Ax + Bu, y = Cx es completamente observable si y s´ olo si la matriz simétrica V (t0 , t1 ) =

t1

t0

ΦT (τ, t0 )C T CΦ(τ, t0 ) dτ ,

llamada matriz de observabilidad, es no singular. 37

(5.3)

38

Observabilidad

Demostraci´ on. Demostremos en primer lugar la suficiencia, es decir, supondremos que V (t0 , t1 ) es no singular y hemos de demostrar que el sistema es completamente observable. Por la Nota 28 tomaremos el control u(t) ≡ 0 para on del sistema x˙ = Ax sujeta a la condici´ on t0 ≤ t ≤ t1 . Entonces, la soluci´ inicial x(t0 ) = x0 viene dada por x(t) = Φ(t, t0 )x0 . De este modo, la salida es on por la izquierda por ΦT (t, t0 )C T y(t) = CΦ(t, t0 )x0 . Multiplicando esta ecuaci´ e integrando respecto de t se tiene t1 ΦT (τ, t0 )C T y(t) dτ = V (t0 , t1 )x0 . t0

Puesto que, por hip´ otesis, V (t0 , t1 ) es no singular, se tiene de la anterior ecuaci´ on t1 x0 = V −1 (t0 , t1 ) ΦT (τ, t0 )C T y(τ ) dτ , (5.4) t0

de modo que el sistema es completamente observable. Demostremos ahora la necesidad. Suponemos que el sistema es completamente observable y probaremos que V (t0 , t1 ) es no singular. Sea α ∈ Rn un vector columna arbitrario. Entonces, la forma cuadr´ atica t1 αT ΦT (τ, t0 )C T CΦ(τ, t0 )α dτ αT V (t0 , t1 )α = t0 t1

=

t0 t1

= t0

(CΦ(τ, t0 )α)T (CΦ(τ, t0 )α) dτ CΦ(τ, t0 )α2 dτ ≥ 0 ,

de modo que la matriz de observabilidad V (t0 , t1 ) es semidefinida positiva. on al absurdo. Veamos que en realidad V (t0 , t1 ) es definida positiva por reducci´ α= Supongamos entonces que existe un vector columna α ˆ = 0 tal que α ˆ T V (t0 , t1 )ˆ α ≡ 0 para t0 ≤ t ≤ t1 . Entonces, la salida 0. Entonces se tiene que CΦ(t, t0 )ˆ ˆ viene dada por y(t) con control nulo u(t) ≡ 0 para t0 ≤ t ≤ t1 cuando x0 = α y(t) = CΦ(t, t0 )x0 = CΦ(t, t0 )ˆ α≡0, para t0 ≤ t ≤ t1 . Esto implica que x0 no puede determinarse de forma u ńica ˆ entonces también se a partir del conocimiento de y(t) ya que si x0 = 0 = α on con tiene y(t) = 0 para t0 ≤ t ≤ t1 . Se ha llegado pues a una contradicci´ la hip´ otesis de que el sistema es completamente observable y, en consecuencia, V (t0 , t1 ) es definida positiva. En particular V (t0 , t1 ) es no singular. Nota 30. Observemos que, cuando el sistema x˙ = Ax + Bu es completamente observable, la ecuaci´ on (5.4) da una expresi´ on expl´ıcita para el estado inicial x0 independiente de la matriz B. Entonces, es habitual hablar de la observabilidad del par [A, C] en lugar del sistema x˙ = Ax + Bu, y = Cx.

5.1 Observabilidad

39

Nota 31. Como se observa de la demostraci´ on del Teorema 29, cuando u(t) ≡ 0, olo se el estado inicial x0 dado por (5.4) es independiente del tiempo final t1 . S´ requiere t1 > t0 . El siguiente resultado es muy u ´til ya que permite deducir a partir de resultados de controlabilidad el correspondiente de observabilidad y viceversa. Teorema 32 (Dualidad). El sistema x˙ = Ax + Bu, y = Cx es completamente controlable si y s´ olo si el sistema dual x˙ = −AT x + C T u , y = B T x

(5.5)

es completamente observable. Demostraci´ on. Es f´ acil demostrar que, si Φ(t, t0 ) es la matriz de transición de estados para el sistema x˙ = Ax, entonces ΦT (t, t0 ) es la matriz de transición de estados para el sistema x˙ = −AT x. De este modo, comparando las expresiones (4.7) y (5.3) de las matrices de controlabilidad U (t0 , t1 ) y de observabilidad V (t0 , t1 ), es decir, U (t0 , t1 ) = y

t0

V (t0 , t1 ) =

t1

t1 t0

Φ(t0 , τ )BB T ΦT (t0 , τ ) dτ

ΦT (τ, t0 )C T CΦ(τ, t0 ) dτ ,

se observa que la matriz de controlabilidad del par [A, B] es idéntica a la matriz de observabilidad del par [−AT , B T ]. Adem´ as, la matriz de observabilidad del par [A, C] es idéntica a la matriz de controlabilidad del par [−AT , C T ], con lo cual el teorema queda probado. Una aplicaci´ on inmediata del Teorema 32 de dualidad es el siguiente. Teorema 33. El sistema x˙ = Ax + Bu, y = Cx es completamente observable si y s´ olo si la matriz de observabilidad ⎡ ⎤ C ⎢ CA ⎥ ⎢ ⎥ ⎢ CA2 ⎥ V =⎢ (5.6) ⎥ ∈ Mnr×n (R) ⎢ .. ⎥ ⎣ . ⎦ CAn−1

tiene rango m´ aximo, es decir, rangV = n. Demostraci´ on. Por el Teorema 32 de dualidad, el sistema x˙ = Ax + Bu, y = Cx es completamente observable si y sólo si el sistema dual x˙ = −AT x + C T u, y = B T x es completamente controlable. Utilizando ahora el Teorema 11 aplicado al sistema dual se tiene el resultado deseado.

40

Observabilidad

5.2.

Problemas resueltos

1. Consideremos un circuito RLC donde la resistencia R y el condensador de capacidad C se han colocado en paralelo y la bobina de autoinducci´ on L en serie. Si el circuito se conecta a una pila de potencial V , es f´ acil demostrar que la diferencia de potencial VC en los extremos del condensador y la intensidad IL que circula por la bobina verifican el sistema 1 1 2 V˙ C − RC VC C RC = + V . 1 IL − L1 0 I˙L L Supongamos que la salida del sistema es −VC y el control V . (i) Averiguar para qué valores de R, L y C el circuito es completamente controlable y completamente observable. (ii) ¿Es posible, cuando el sistema no es completamente observable, obtener V (t) de modo que, si el circuito tiene inicialmente (VC (0), IL (0)) = (1, 5) pase a (VC (10), IL (10)) = (2, 3)? Soluci´ on. (i) El sistema se puede escribir de la forma x˙ = Ax + bV , y = Cx, donde x = (VC , IL )T y siendo 1 2 1 − RC C RC A= , C = (−1, 0) . , b= 1 − L1 0 L Puesto que A ∈ M2 (R), la matriz de controlabilidad U del sistema lineal del enunciado es 1 1 − R22C 2 RC LC U = [b, Ab] = . 1 1 − RLC L Como

1 1 − 2 , det U = 2 2 R LC L C se tiene que si R= L/C entonces det U = 0 y rang U = 1. Por el contrario, si R = L/C entonces det U = 0 y rang U = 2, de modo que, aplicando el Teorema 11, se concluye que el sistema es completamente controlable. La matriz de observabilidad V es

−1 C V = = 2 CA RC

0 − C1

,

de modo que siempre tiene rango máximo, es decir, rang V = 2. En conclusión, aplicando el Teorema 33, el sistema es siempre completamente observable.


41

(ii) Para resolver esta cuestión aplicaremos el Teorema 19. Nos situamos en el caso no completamente controlable, es decir, tomamos R = L/C. Es f´ acil calcular la matriz de transici´ on de estados Φ(0, t) = exp(−At), siendo esta ⎞ √ ⎛ L+ L/Ct t L/Ct ⎝ − L ⎠ . √C Φ(0, t) = exp L− L/Ct L t L

L

Entonces, la matriz de controlabilidad U (0, 10)

=

=

10

Φ(0, τ )bbT ΦT (0, τ ) dτ 0 √ 20 L/C −1 + exp L L/C 1 2L

1 1/ L/C

.

Por supuesto det U (0, 10) = 0 puesto que el sistema no es completamente controlable. Sea x0 = (1, 5)T y xf = (2, 3). Calculemos ahora el vector ⎞ ⎛ √ √ 20 L/C 10 L/C 30 − −2 + 1 + exp L C L ⎟ ⎜ ⎟ . √ √ x0 − Φ(0, 10)xf = ⎜ ⎠ ⎝ 10 L/C −20−3L+30 L/C 5L + exp L L Es f´ acil ver que el sistema lineal algebraico U (0, 10)ξ = x0 − Φ(0, 10)xf , es incompatible, es decir, no existe ning´ un vector ξ ∈ R2 que sea su soluci´ on. En definitiva, aplicando el Teorema 19 se concluye que no es posible llevar el circuito desde el estado inicial x0 = (VC (0), IL (0)) = (1, 5) al estado final xf = (VC (10), IL (10)) = (2, 3) en 10 unidades de tiempo. 2. Consideremos el sistema lineal x˙ = Ax, siendo x = (x1 , x2 )T ∈ R2 y −1 −1 A= , 2 −4 donde la salida viene dada por y(t) = x1 (t) + 2x2 (t). ¿Es posible hallar el estado inicial del sistema x(0) sabiendo que la salida es y(t) = −20 exp(−3t) + 21 exp(−2t)? En caso de respuesta positiva, hallar x(0). Soluci´ on. Se tiene el sistema x˙ = Ax con salida y = Cx, siendo C = (1, 2). La matriz de observabilidad V asociada al sistema es

C 1 2 V = = , CA 3 −9

42

Observabilidad que tiene rango m´ aximo, es decir, rang V = 2. Entonces, aplicando el Teorema 33, el sistema es siempre completamente observable. Ello significa que siempre es posible hallar el estado inicial del sistema x(0) conociendo la salida y(t). Como el sistema es completamente observable y además el control es u(t) ≡ 0, se sabe, ver (5.4), que el estado inicial x(0) = x0 viene dado por x0 = V −1 (0, t1 )

t1

0

ΦT (τ, t0 )C T y(t) dτ ,

donde el tiempo final t1 no juega ning´ un papel esencial (ver la Nota 31) de modo que, se toma t1 = 1 sin pérdida de generalidad. Aqu´ı, la matriz de observabilidad V (0, 1) viene dada, seg´ un (5.3), por

1

V (0, 1) =

ΦT (τ, 0)C T CΦ(τ, 0) dτ ,

0

siendo Φ(τ, 0) = exp(Aτ ) la matriz de transici´ on de estados. Es fácil calcular la matriz exponencial obteniéndose −1 + 2 exp(τ ) 1 − exp(τ ) Φ(τ, 0) = exp(−3τ ) . 2[−1 + exp(τ )] 2 − exp(τ ) La matriz de observabilidad V (0, 1) es T −1 + 2 exp(τ ) 1 − exp(τ ) 1 V (0, 1) = exp(−6τ ) 2[−1 + exp(τ )] 2 − exp(τ ) 2 0 −1 + 2 exp(τ ) 1 − exp(τ ) (1, 2) dτ 2[−1 + exp(τ )] 2 − exp(τ ) " 2 + 25+27(e−2)e 7 + −25+18(4−3e)e 1 e6 e6 ! . = −50+9(8−3e)e 1 2 + 25+27(e−2)e 6 e6 2 5+ e6

1

Finalmente, teniendo en cuenta que la salida es y(t) = −20 exp(−3t) + 21 exp(−2t), se tiene x0 = V −1 (0, 1)

1 0

ΦT (τ, 0)C T y(t) dτ =

3 −1

.

Cap´ıtulo 6

Realimentaci´ on Lineal 6.1.

Definici´ on y preliminares

Definici´ on 34. Al sistema x˙ = Ax + Bu se le aplica una realimentaci´ on lineal si cada variable de control es una combinaci´ on lineal de las variables de estado, es decir, u(t) = Kx(t) , (6.1) siendo K ∈ Mm×n (R) una matriz constante llamada matriz de ganancia. El sistema resultante x˙ = (A + BK)x , (6.2) se llama sistema de lazo cerrado. En este cap´ıtulo se utilizar´ a la siguiente notaci´ on. a un conjunto arbitrario de n Definici´ on 35. Λn = {ρ1 , . . . , ρn } ⊂ C denotar´ n´ umeros complejos tales que cualquiera de ellos que no sea real aparece emparejado con su complejo conjugado. El siguiente resultado, cuya demostraci´ on general data del a¯ no 1967 [13], proporciona una nueva v´ıa de investigaci´ on en control mediante el uso de realimentaci´ on lineal en términos de variables de estado. Teorema 36. Si el sistema x˙ = Ax + Bu es completamente controlable entonces existe una matriz real K tal que los valores propios de la matriz A + BK on 35. pertenecen al conjunto Λn de la Definici´ Nota 37. Sabemos que, la soluci´ on del sistema lineal de lazo cerrado x˙ = (A + BK)x depende de los valores propios de la matriz A+BK. Por lo tanto, si el par [A, B] es completamente controlable, del Teorema 36 se desprende que, usando realimentaci´ on lineal, es posible ejercer una gran influencia sobre la evoluci´ on temporal de las soluciones del sistema lineal de lazo cerrado. 43

44

Realimentaci´ on Lineal

Demostraci´ on del Teorema 36 con una u ńica variable de control (m = 1). Como el sistema x˙ = Ax + bu con A ∈ Mn (R) matriz constante, b ∈ Rn vector columna y u(t) ∈ R es completamente controlable, por el Teorema 14, existe cambio lineal de variables w = T x con T ∈ Mn (R) no singular que lo transforma en la forma canónica controlable (3.25), es decir, w˙ = Cw + du con C = T AT −1 y d = T b. En concreto d = (0, . . . , 0, 1)T y ⎞ ⎛ 0 1 0 . . . 0 ⎜ 0 0 1 . . . . ⎟ ⎟ ⎜ ⎜ . . . . . . . ⎟ C=⎜ ⎟ ∈ Mn (R) . ⎝ . . . . . . . ⎠ −kn −kn−1 −kn−2 . . . −k1 Aqu´ı, ki son los coeficientes del polinomio caracter´ıstico de A, es decir, det(λIn − A) = λn + k1 λn−1 + · · · + kn . Por otra parte, la realimentaci´ on lineal en el caso m = 1 que nos ocupa se escribe de la forma u = κw, con el vector fila κ = (κn , κn−1 , . . . , κ1 ) ∈ Rn . De este modo, el sistema en lazo cerrado que se obtiene es w˙ = (C + dκ)w, siendo la matriz ⎞ ⎛ 0 1 0 . . . 0 ⎜ 0 0 1 . . . . ⎟ ⎟ ⎜ ⎜ . . . . . . ⎟ C + dκ = ⎜ . ⎟+ ⎝ . . . . . . . ⎠ −kn −kn−1 −kn−2 . . . −k1 ⎛ ⎞ 0 ⎜ 0 ⎟ ⎜ ⎟ ⎜ .. ⎟ ⎜ . ⎟ (κn , κn−1 , . . . , κ1 ) ⎜ ⎟ ⎝ 0 ⎠ 1 ⎛ ⎞ 0 1 0 . . . 0 ⎜ 0 0 1 . . . . ⎟ ⎜ ⎟ . . . . . . . ⎟ = ⎜ ⎜ ⎟+ ⎝ . . . . . . . ⎠ −kn −kn−1 −kn−2 . . . −k1 ⎛ ⎞ 0 0 0 . . . 0 ⎜ 0 0 0 . . . 0 ⎟ ⎜ ⎟ ⎜ . . . . . . . ⎟ ⎜ ⎟ ⎝ . . . . . . . ⎠ κn κn−1 κn−2 . . . κ1 ⎛ ⎞ 0 1 0 . . . 0 ⎜ 0 0 1 . . . . ⎟ ⎜ ⎟ ⎜ . . . . . . ⎟ = ⎜ . ⎟ , ⎝ . . . . . . . ⎠ −γn −γn−1 −γn−2 . . . −γ1

6.1 Definici´ on y preliminares

45

siendo κi = ki − γi , i = 1, . . . , n .

(6.3)

Notemos que la matriz C + dκ tiene la misma forma que la matriz C excepto la u ´ltima fila. Puesto que los sistemas de lazo cerrado inicial x˙ = (A + bK)x y final w˙ = (C + dκ)w son algebraicamente equivalentes mediante la transformaci´ on lineal w = T x, se tiene que C + dκ = T (A + bK)T −1 , de donde κ = KT −1 y por lo tanto la matriz de ganancia K que se busca es K = κT ,

(6.4)

donde los elementos de κ están definidos por (6.3) y los elementos γi se calculan imponiendo que los valores propios de la matriz A + bK (que son los mismos que los de la matriz C + dκ) sean justo el conjunto Λn , es decir, λn + γ1 λn−1 + · · · + γn =

n

(λ − ρi ) ,

i=1

con lo cual queda demostrado. Nota 38. Observar que la demostraci´ on efectuada del Teorema 36 con una u ńica variable de control (m = 1) es constructiva en el sentido de que existe un algoritmo para calcular la matriz de ganancia K de la realimentaci´ on lineal que nos permita tener un conjunto prefijado de autovalores para la matriz del sistema en lazo cerrado. A partir del Teorema 36 y utilizando el Teorema 32 de dualidad, es f´ acil deducir el siguiente resultado de observabilidad. Corolario 39. Si el sistema x˙ = Ax+Bu, y = Cx es completamente observable, entonces existe una matriz real L tal que los valores propios de la matriz A+LC pertenecen a un conjunto prefijado Λn ⊂ C. Los métodos para controlar sistemas prefijando los autovalores de la matriz del sistema de lazo cerrado se llaman habitualmente técnicas de control modal, donde el nombre “modo”proviene de los términos exp(λi t)wi de la expresión espectral (3.3) de la solución de sistemas lineales, siendo λi y wi valores y vectores propios. Esta idea de control modal da lugar a una nueva definici´ on. Definici´ on 40. El sistema x˙ = Ax + Bu es modalmente completamente controlable si existe una matriz de ganancia K tal que los valores propios de la matriz de lazo cerrado A + BK pertenecen a un conjunto arbitrario Λn ⊂ C como el dado en la Definici´ on 35. El Teorema 36 asegura que si el sistema x˙ = Ax + Bu es completamente controlable entonces también es modalmente completamente controlable. Se puede de hecho demostrar también el rec´ıproco de manera que se tiene una total equivalencia de las dos definiciones de controlabilidad establecidas tal y como el siguiente teorema enuncia, ver por ejemplo [1]. Teorema 41. El sistema x˙ = Ax + Bu es modalmente completamente controlable si y s´ olo si es completamente controlable.

46

Realimentaci´ on Lineal

6.1.1.

Realimentaci´ on y funci´ on de transferencia

Es interesante relacionar el concepto de realimentaci´ on lineal y el de funci´ on de transferencia g(s) para el caso de una u ńica entrada (m = 1) y u ńica salida (r = 1). Recordemos que, seg´ un la teor´ıa del control cl´ asica, la funci´ on de transferencia viene dada por la funci´ on racional (2.7), es decir, g(s) =

β0 sm + β1 sm−1 + · · · + βm−1 s + βm β(s) = . k(s) sn + k1 sn−1 + · · · + kn−1 s + kn

Teorema 42. Si se tiene una u ńica entrada (m = 1) y u ńica salida (r = 1), la realimentaci´ on lineal permite elegir de forma arbitraria en el plano complejo los polos de la funci´ on de transferencia g(s) del sistema de lazo cerrado pero deja sus ceros invariantes. Demostraci´ on. Es f´ acil demostrar que ˆ n − C)−1 d , g(s) = β(sI

(6.5)

siendo d = (0, 0, . . . , 0, 1)T y βˆ = (βm , . . . , β1 , β0 , 0, . . . , 0). De este modo, g(s) es la funci´ on de transferencia del sistema ˆ . x˙ = Cx + du , y = βx Como vimos en la demostración del Teorema 36 con una u ńica variable de control (m = 1), la realimentaci´ on lineal s´ olo afecta a la u ´ltima columna de la on de matriz C, la cual contiene los coeficientes ki del denominador de la funci´ transferencia g(s). Nota 43. El Teorema 42 no es cierto si m > 1, puesto que en este caso los ceros de la funci´ on de transferencia son afectados también por la realimentaci´ on lineal.

6.2.

Realimentaci´ on frente a control precalculado

En esta sección veremos, bajo la ´ optica del clásico ejemplo del péndulo invertido, las ventajas de tener una ley de control por realimentaci´ on u = Kx frente a un control precalculado como una funci´ on del tiempo u(t). Consideremos un péndulo de masa m, formando un ańgulo θ(t) con la vertical para tiempo t. Sea u(t) un momento externo que act´ ua sobre el péndulo y podemos modificar a voluntad, dicho de otro modo, u es nuestra variable de control. Bajo la hip´ otesis de rozamiento despreciable y gravedad g constante, la ecuación del movimiento del péndulo viene dada, seg´ un las leyes de la mecánica clásica de Newton, por la siguiente ecuación diferencial no lineal de segundo

6.2 Realimentaci´ on frente a control precalculado

47

orden mθ¨ + mg sin θ = u(t). Tomemos las unidades de espacio y de tiempo de forma que m = g = 1, con lo cual el sistema a estudiar es θ¨ + sin θ = u(t) .

(6.6)

˙ = (π, 0) es un punto inestable Si u(t) ≡ 0, el punto de equilibrio (θ, θ) puesto que una peque¯ na desviación de dicho punto provocar´ıa que el péndulo no volviera al punto. El problema de control que se pretende abordar consiste en actuar sobre el péndulo mediante el momento u(t) de modo que, para valores iniciales ˙ (θ(0), θ(0)) cercanos al punto (π, 0), la evoluci´ on del péndulo sea hacia (π, 0). En otras palabras se pretende hallar u(t) de modo que el punto de equilibrio (π, 0) sea estable. Realizemos el cambio de variable ϕ = θ −π. Con este cambio, el punto cr´ıtico que se ha de analizar es justo el origen (ϕ, ϕ) ˙ = (0, 0). Además, sin θ = − sin ϕ, de modo que la ecuación a estudiar es ϕ¨ − sin ϕ = u(t) .

(6.7)

on, que se realizan peque¯ nas oscilaSupongamos, en primera aproximaci´ ciones, es decir, que ϕ(t) es un ángulo peque¯ no en valor absoluto. Realizando la aproximaci´ on habitual sin ϕ ≈ ϕ, podemos sustituir la ecuación no lineal (6.7) por la siguiente ecuación lineal ϕ¨ − ϕ = u(t) .

(6.8)

Se pretende dise¯ nar una estrategia de control u(t) de modo que, para valores iniciales (ϕ(0), ϕ(0)) ˙ no nulos pero peque¯ nos, el sistema evolucione hacia (ϕ, ϕ) ˙ = (0, 0). Supongamos que se pretende controlar el sistema linealizado (6.8) sabiendo que parte de las condiciones iniciales (ϕ(0), ϕ(0)) ˙ = (1, −2) .

(6.9)

Es f´ acil comprobar que el control u(t) = 3 exp(−2t) ,

(6.10)

es adecuado para nuestros prop´ ositos puesto que la solución de (6.8) es, en este caso, ϕ(t) = exp(−2t) de modo que está claro que (ϕ(t), ϕ(t)) ˙ → (0, 0) cuando t → ∞. Es f´ acil mostrar que un control por realimentaci´ on del tipo u(t) = −αϕ(t) − β ϕ(t) ˙ ,

(6.11)

48

Realimentaci´ on Lineal con α > 1 y β > 0 también cumple nuestro prop´ osito, es decir, las soluciones de (6.8) con condiciones iniciales (ϕ(0), ϕ(0)) ˙ suficientemente cercanas a (0, 0) verifican (ϕ(t), ϕ(t)) ˙ → (0, 0) cuando t → ∞. La demostración de este hecho proviene del an´ alisis del sistema en lazo cerrado ϕ+β ¨ ϕ+(α−1)ϕ ˙ = 0, con ecuación caracter´ıstica λ2 +βλ+α−1 = 0 cuyas ra´ıces son λ± = [−β ± β 2 − 4(α − 1)]/2. Es obvio que la parte real ˙ → (0, 0) cuando t → ∞. (λ± ) < 0 de modo que (ϕ(t), ϕ(t))

6.2.1.

Sensibilidad a las condiciones iniciales

Veamos qué ocurre si, como es habitual en la experimentaci´ on, existen incertidumbres en las medidas de las condiciones iniciales. Supongamos, por ejemplo, que las nuevas condiciones iniciales son (ϕ(0), ϕ(0)) ˙ = (1, −2 + ) ,

(6.12)

con = 0 y peque¯ no en valor absoluto. En este caso, la soluci´ on de (6.8) con el control (6.10) es

ϕ(t) = exp(−2t) + [exp(t) − exp(−t)] , 2 que verifica

l´ım ϕ(t) =

t→∞

∞ si −∞ si

>0,

<0,

de modo que ahora el control u(t) no cumple nuestro prop´ osito.

6.2.2.

Sensibilidad a perturbaciones externas

Supongamos ahora que, durante un cierto instante de tiempo, el control u(t) que efect´ ua por ejemplo un motor, se perturba por cualquier motivo mediante un valor > 0 peque¯ no como puede ser una subida de tensi´ on en la red eléctrica. Matemáticamente podemos realizar la siguiente modelización de una situaci´ on similar ϕ¨ − ϕ = u(t) + δ(t), (6.13) siendo

δ(t) =

si t ∈ [1, 2] , 0 si t ∈ [1, 2] .

(6.14)

Vamos a probar que, si se toma el control (6.10) entonces, la solución de (6.13) diverge con lo cual no se verifican nuestros prop´ ositos. En efecto, la soluci´ on de (6.13) con condiciones iniciales (6.9), es decir, (ϕ(0), ϕ(0)) ˙ = (1, −2), es ϕ(t) = exp(−2t) para t ∈ [0, 1]. Posteriormente, se resuelve (6.13) con la condición inicial (ϕ(1), ϕ(1)) ˙ =


49

(exp(−2), −2 exp(−2)) durante el intervalo de tiempo t ∈ [1, 2]. Se obtiene entonces

1 1 ϕ(t) = exp(−2t) + −1 + exp(1 − t) + exp(t − 1) . 2 2 Finalmente, se resuelve (6.13) con la condición inicial

e 1

1 −4 −4 (ϕ(2), ϕ(2)) ˙ = e + −1 + + , −2e + e− 2e 2 2 e durante el intervalo de tiempo t ≥ 2. Se obtiene entonces 1 ϕ(t) = e−2t + (e − 1) et−2 − e1−t para t ≥ 2 , 2 con lo que l´ımt→∞ ϕ(t) = ∞. De una forma similar se puede demostrar que, si se toma por control u(t) la realimentaci´ on (6.11), es decir, u(t) = −αϕ(t) − β ϕ(t) ˙ con α > 1 y β > 0, entonces la solución de (6.13) no s´ olo no diverege sino que todav´ıa mantiene la propiedad asint´ otica deseada (ϕ(t), ϕ(t)) ˙ → (0, 0) cuando t → ∞.

6.3.

Problemas resueltos

1. ¿Es invariante la observabilidad de un sistema respecto de la realimentaci´ on lineal de estados? Responder después del estudio del sistema 1 2 0 x˙ = x+ u(t) , 3 1 1 con salida y = (1, 2)x y realimentaci´ on u(t) = (3, 1)x. Soluci´ on. El sistema x˙ =

1 3

2 1

x+

0 1

u(t) ,

es completamente observable puesto que la matriz de observabilidad V es 1 2 V = , 7 4 y tiene rango m´ aximo. Tras la realimentaci´ on, el sistema en lazo cerrado es 1 2 x˙ = x , y = (1, 2)x , 0 0

50

Realimentaci´ on Lineal siendo su matriz de observabilidad 1 VK = 1

2 2

.

Es claro que el rango de VK no es máximo de modo que el sistema de lazo cerrado no es controlable. La respuesta a la pregunta ¿Es invariante la observabilidad de un sistema respecto de la realimentación lineal de estados? es NO. 2. Consideremos el sistema lineal −1 −1 x1 1 x˙ 1 = + u(t) . x˙ 2 x2 2 −4 3 Averiguar si es posible realizar una realimentaci´ on lineal de modo que el sistema en lazo cerrado tenga asociados los valores propios ρ1 = −4 y ρ2 = −5. En caso de que sea posible, hallarla. Soluci´ on. Nos basaremos en la demostración constructiva del Teorema 36 con una u ńica variable de control (que es nuestro caso). Sea x = (x1 , x2 )T ∈ R2 . El sistema a estudiar es x˙ = Ax + bu con −1 −1 1 A= ∈ M2 (R) , b = ∈ R2 . 2 −4 3 Si el sistema es completamente controlable, se sabe por el Teorema 36 que existe una realimentaci´ on lineal de modo que el sistema en lazo cerrado tenga asociados los valores propios que deseemos de un conjunto Λ2 y por acil ver lo tanto, en particular, los valores propios ρ1 = −4 y ρ2 = −5. Es f´ que el sistema del enunciado es completamente controlable. Entonces, por el Teorema 14, existe un cambio lineal de variables w = onica T x con T ∈ M2 (R) no singular que lo transforma en la forma can´ controlable (3.25), es decir, w˙ = Cw + du con d = (0, 1)T ∈ R2 y 0 1 C= ∈ M2 (R) . −k2 −k1 Aqu´ı, ki son los coeficientes del polinomio caracter´ıstico de A, es decir, det(λI2 − A) = λ2 + k1 λ + k2 = λ2 + 5λ + 6 , de modo que k1 = 5 , k2 = 6 . La matriz de ganancia es K = κT , siendo el vector fila κ = (κ2 , κ1 ) ∈ R2 donde (6.15) κi = ki − γi , i = 1, 2 ,


51

y γi se calcula imponiendo que λ2 + γ1 λ + γ2 =

2

(λ − ρi ) ,

i=1

siendo ρi los valores propios que se desean en lazo cerrado. En concreto se tiene λ2 + γ1 λ + γ2 = (λ + 4)(λ + 5) , de donde obtenemos γ1 = 9 , γ2 = 20 . Volviendo a (6.15) se tiene κ1 = 5 − 9 = −4 , κ2 = 6 − 20 = −14 , de modo que κ = (−14, −4). S´ olo falta por calcular la matriz T . Para conseguirlo, recordemos la Nota 7, es decir, T se puede construir mediante (3.29) donde ξ es la solución de (3.31). En otras palabras,

ξ T = ∈ M2 (R) , ξA donde ξ verifica

ξ[b, Ab] = dT .

En nuestro problema, despejando ξ de la anterior ecuaci´ on se tiene −1 1 −4 T −1 ξ = d [b, Ab] = (0, 1) 3 −10 1 −5 2 = (0, 1) = (−3, 1) . −3/2 1/2 2 Finalmente, se tiene

T =

ξ ξA

1 = 2

−3 5

1 −1

.

Resumiendo, la realimentaci´ on lineal ser´ a del tipo u = Kx, siendo la matriz de ganancia 1 −3 1 = (11, −5) , K = κT = (−14, −4) 5 −1 2 es decir, el control por realimentaci´ on lineal utilizado ser´ a x1 u = Kx = (11, −5) = 11x1 − 5x2 . x2

52

Realimentaci´ on Lineal Para comprobar que todo es correcto se puede ver que los valores propios de la matriz 10 −6 A + bK = 35 −19 asociada al sistema en lazo cerrado son ρ1 = −4 y ρ2 = −5.

Cap´ıtulo 7

Estabilidad 7.1.

Introducci´ on y definiciones

A diferencia de los cap´ıtulos anteriores, los conceptos de este cap´ıtulo pueden ser aplicados a sistemas más generales que los sistemas lineales. Consideraremos pues un sistema (en general no lineal) de ecuaciones diferenciales de primer orden x˙ = f (x, t) , (7.1) siendo x(t) = (x1 (t), . . . , xn (t)) ∈ Rn el vector de estados y f (x, t) = (f1 (x, t), . . . , fn (x, t)) ∈ Rn satisfaciendo las condiciones del teorema de existencia y unicidad del problema de Cauchy asociado para una condici´ on inicial x(t0 ) = x0 dada. En particular, si fi y ∂fi /∂xj son funciones continuas en un entorno del punto (t0 , x0 ) ∈ Rn+1 se verifican dichas condiciones. El sistema (7.1) se llama aut´ onomo si f no depende expl´ıcitamente del tiempo, es decir, son sistemas del tipo x˙ = f (x). En caso contrario se llaman no aut´ onomos. Si la funci´ on f es lineal y aut´ onoma, entonces se recuperan los sistemas lineales x˙ = Ax estudiados anteriormente. Definici´ on 44. El estado x∗ ∈ Rn se llama punto cr´ıtico o punto singular o estado estacionario del sistema (7.1) si verifica f (x∗ , t) ≡ 0 para todo t. Siempre supondremos que los puntos cr´ıticos son aislados, lo que significa que no hay otros puntos cr´ıticos de (7.1) en un entorno de x∗ . Nota 45. Supongamos que x∗ es un punto cr´ıtico del sistema (7.1). Entonces, on (7.1) se deduce que x(t) = x∗ para todo si x(t0 ) = x∗ , utilizando la ecuaci´ t ≥ t0 . El hecho de que las soluciones de (7.1) que comiencen en x∗ permanezcan en x∗ da el nombre de estado estacionario a x∗ . Nota 46. Sin pérdida de generalidad siempre se puede suponer que el punto cr´ıtico x∗ es el origen, es decir x∗ = 0 ∈ Rn . En caso de que que no lo sea, se realiza el cambio de variables (traslaci´ on) z = x − x∗ que lleva el punto x∗ al origen. 53

54

Estabilidad

Intuitivamente, la estabilidad de un punto cr´ıtico x∗ de (7.1) medir´ a si las soluciones x(t) de (7.1) que comiencen “cerca”de x∗ permanecerán o no “cerńica definici´ on de estabilidad. Nosotros utilizareca”para t > t0 . No existe una u mos la siguiente definici´ on dada por Liapunov. Definici´ on 47. Sea x∗ = 0 ∈ Rn punto cr´ıtico del sistema (7.1). Entonces: (i) 0 es estable si para todo > 0 existe un δ > 0 tal que x(t0 ) < δ implica x(t) < para todo t ≥ t0 . (ii) 0 es asint´ oticamente estable si es estable y adem´ as l´ımt→∞ x(t) = 0. (i) 0 es inestable si no es estable, es decir, existe un > 0 tal que para todo un δ > 0 existe un x(t0 ) verificando x(t0 ) < δ y x(t1 ) ≥ para alg´ as esto sucede para todo x(t0 ) con x(t0 ) < δ entonces 0 t1 > t0 . Si adem´ se llama completamente inestable. Si se conoce la expresión expl´ıcita de la solución general x(t) del sistema (7.1) para cualquier condici´ on inicial x(t0 ) = x0 entonces también se conoce la estabilidad de sus puntos cr´ıticos. Por ejemplo, el sistema lineal x˙ = Ax, siendo 2 −3 A= , 2 −1 con la condici´ on inicial x(0) = (x1 (0), x2 (0)) tiene por soluci´ on general (utilizando los métodos anteriormente estudiados) √ " √ " 3 t 15 15 t + √ sin t − x1 (t) = x1 (0) exp cos 2 2 2 15 √ " t 6 15 √ x2 (0) exp t , sin 2 2 15 con una expresi´ on similar para x2 (t). Obviamente, debido al término exponencial, l´ımt→∞ x1 (t) = ∞ para cualquier valor inicial x(0) de modo que el origen del sistema lineal es inestable. Sin embargo, excepto para sistemas lineales y alg´ un otro tipo simple, no se conocerá a priori la soluci´ on general x(t) del sistema de modo que este método no es de gran aplicaci´ on.

7.2.

Estabilidad en sistemas lineales

Consideremos otra vez el sistema lineal x˙ = Ax ,

(7.2)

7.2 Estabilidad en sistemas lineales

55

con x ∈ Rn y A ∈ Mn (R). Suponiendo que A sea no singular, es claro que el u ńico punto cr´ıtico del sistema (7.2) es el origen, de modo que tiene sentido hablar de la estabilidad del sistema lineal (7.2) aunque nos refiramos a la estabilidad del origen de (7.2). Definici´ on 48. Sean λk ∈ C, con k = 1, . . . , n, los valores propios de la matriz A ∈ Mn (R). Denotemos por (λk ) a la parte real de λk . Entonces A se llama matriz de estabilidad o bien de Hurwitz si (λk ) < 0 para todo k = 1, . . . , n. Teorema 49. Sean λk , con k = 1, . . . , n, los valores propios de la matriz A. El sistema (7.2) es: (i) asint´ oticamente estable si y s´ olo si A es una matriz de estabilidad. (ii) inestable si existe alg´ un k tal que (λk ) > 0. (iii) completamente inestable si (λk ) > 0 para todo k = 1, . . . , n. olo probaremos la condici´ on suficiente de (i) en el caso Demostraci´ on de (i). S´ de que todos los valores propios de A sean diferentes. La solución del sistema (7.2) con la condici´ on inicial x(0) = x0 viene dada por x(t) = exp(At)x0 . Utilizando las propiedades de las normas matriciales se tiene x(t) ≤ exp(At) x0 .

(7.3)

Realizaremos la demostración sólo para el caso en que todos los valores propios λk de A son distintos. Entonces se puede aplicar la f´ ormula de Sylvester (3.14), es decir, n Zk exp(λk t) , exp(At) = k=1

siendo Zk matrices constantes que sólo dependen de A y de sus valores propios. Tomando normas se tiene n Zk | exp(λk t)| . exp(At) ≤ k=1

Como λk = (λk ) + i (λk ), siendo i =

√

−1, se tiene que

exp(λk t) = exp( (λk )t) exp(i (λk )t) , de modo que el módulo | exp(λk t)| = exp( (λk )t). Se tiene pues que exp(At) ≤

n

Zk exp( (λk )t) .

(7.4)

k=1

Supongamos que (λk ) < 0 para todo k. Entonces, de (7.4) se deduce que exp(At) → 0 cuando t → ∞. Finalmente, teniendo en cuenta (7.3), se concluye que x(t) → 0 cuando t → ∞ y por lo tanto el sistema (7.2) es asint´ oticamente estable.

56

Estabilidad

Nota 50. Observar que tras el Teorema 49, los casos dif´ıciles en el estudio de la estabilidad del sistema (7.2) son aquellos en los que alg´ un valor propio λj de A sea o bien nulo o bien imaginario puro, es decir, (λj ) = 0. Es f´ acil demostrar el siguiente resultado de gran utilidad en aplicaciones. Proposici´ on 51. Sea el polinomio k(λ) = λ3 + aλ2 + bλ + c ∈ R3 [λ]. Todas las ra´ıces de k(λ) tienen parte real negativa si y s´ olo si a > 0, b > 0, c > 0 y ab > c.

7.3.

Teor´ıa de Liapunov de la estabilidad

En esta sección estudiaremos la estabilidad de los puntos cr´ıticos de un sistema (7.1) no necesariamente lineal pero si autónomo. Sin pérdida de generalidad se estudiar´ a en muchas ocasiones la estabilidad del origen, ver la Nota 46. Por lo tanto, el ojetivo es estudiar la estabilidad del punto cr´ıtico x∗ en el sistema (7.5) x˙ = f (x) , f (x∗ ) = 0 . La teor´ıa de Liapunov de la estabilidad para este tipo de sistemas consiste en generalizar el concepto de energ´ıa V que aparece en los sistemas conservativos en mecánica, donde es bien sabido que un punto cr´ıtico es estable si la energ´ıa tiene un m´ınimo en dicho punto. Definici´ on 52. Dada una funci´ on V : Rn → R de clase C 1 , se define la ˙ derivada orbital V de V respecto del sistema (7.5) como la siguiente derivada direccional n ∂V dV = x˙ i = ∇V.f . V˙ = dt ∂xi i=1 Nota 53. N´ otese que, dada una soluci´ on x(t) del sistema (7.5), se tiene dV (x(t)) . V˙ (x(t)) = dt on Definici´ on 54. La funci´ on V : D ⊂ Rn → R con x∗ = 0 ∈ D es una funci´ de Liapunov para el sistema (7.5) si verifica lo siguiente: (i) V ∈ C 1 (D). (ii) V es definida positiva, es decir, V (0) = 0 y V (x) > 0 para todo x ∈ D con x = 0. (iii) V˙ es semidefinida negativa, es decir, V˙ (0) = 0 y V˙ (x) ≤ 0 para todo x ∈ D. El siguiente resultado es b´ asico en esta teor´ıa y posee un fuerte caracter geométrico.

7.3 Teor´ıa de Liapunov de la estabilidad

57

Teorema 55 (Liapunov). (i) El punto cr´ıtico x∗ = 0 del sistema (7.5) es estable si existe una funci´ on de Liapunov definida en un entorno D ⊂ Rn con 0 ∈ D. oticamenete estable si (ii) El punto cr´ıtico x∗ = 0 del sistema (7.5) es asint´ existe una funci´ on de Liapunov definida en un entorno D ⊂ Rn con 0 ∈ D tal que V˙ es definida negativa, es decir, V˙ (0) = 0 y V˙ (x) < 0 para todo x ∈ D con x = 0. Demostraci´ on. Daremos una prueba s´ olo para el caso bidimensional, es decir x = (x1 , x2 ) ∈ R2 . En primer lugar exploremos algunos aspectos geométricos de las funciones definidas positivas. Est´ a claro que, si V : R2 → R es definida positiva en un 2 entorno D ⊂ R del origen, entonces V tiene un m´ınimo aislado en el origen. Esto implica que las curvas de nivel V −1 (c) = {x ∈ R2 | V (x) = c} en el plano de fases x1 − x2 son curvas cerradas rodeando al origen para c ∈ R positivo y peque¯ no. En el siguiente paso, queremos determinar c´ omo cruzan las soluciones x(t) del sistema x˙ = f (x) a las anteriores curvas de nivel de V . Para verlo, calculamos la derivada orbital de V seg´ un el sistema x˙ = f (x), es decir, ∂V ∂V (x(t))x˙ 1 + (x(t))x˙ 2 . V˙ (x(t)) = ∂x1 ∂x2 Por supuesto, esta expresión coincide con la del producto escalar del vector gradiente ∇V (x) y del vector f (x) de modo que V˙ (x) = f (x).∇V (x) = f (x) ∇V (x) cos θ , siendo θ el ángulo formado por los vectores ∇V (x) y f (x). Recordemos que el vector ∇V (x) es perpendicular a la curva de nivel de V que pasa por x y con el sentido de máxima variaci´ on de V . De este modo, si V˙ (x) < 0, entonces el angulo θ es obtuso con lo cual la órbita del sistema que pasa por x corta a la ´ curva de nivel desde la parte del plano no acotada (exterior de la curva de nivel) hacia la parte acotada por la curva de nivel (interior de la curva de nivel). De forma similar, si V˙ (x) > 0, la órbita del sistema que pasa por x corta a la curva de nivel de V que pasa por x desde su interior hacia su exterior; si V˙ (x) = 0 la órbita del sistema que pasa por x toca a la curva de nivel tangencialmente en x. Con todos estos preliminares vamos finalmente a demostrar el teorema. (i) Sea > 0 suficientemente peque¯ no de manera que el entorno del origen consistente en los puntos x con x ≤ esté contenido en D. Sea m = m´ın {V (x)} > 0 . x=

58

Estabilidad

Recordemos que m existe puesto que el subconjunto de R2 de puntos x con x = es cerrado y acotado y por lo tanto compacto. Ahora tomemos un δ tal que 0 < δ ≤ tal que V (x) < m para todo x con x ≤ δ. Dicho δ existe ya que V es una funci´ on continua con V (0, 0) = 0. Sea x(t) la solución del sistema que pasa por el punto x0 ∈ R2 para tiempo t = 0. Si se toma x0 ≤ δ entonces x(t) ≤ para todo t ≥ 0 ya que V˙ (x(t)) ≤ 0 implica V (x(t) ≤ V (x0 ) para todo t ≥ 0. Hemos probado pues la estabilidad del origen. (ii) Se demuestra de manera similar. Por supuesto, una de las principales dificultades que uno se encuentra cuando quiere utilizar la teor´ıa de estabilidad de Liapunov consiste en buscar la función de Liapunov. Una situaci´ on simple consiste en buscar formas cuadráticas definidas positivas como funciones de Liapunov. Lema 56. Una forma cuadr´ atica V (x1 , x2 ) = ax21 + 2bx1 x2 + cx22 con a, b, c ∈ R es definida positiva si y s´ olo si a > 0 y ac − b2 > 0. Demostraci´ on. Probaremos sólo la necesidad (la suficiencia sigue razonamientos análogos). Supongamos pues que V es definida positiva. Entonces, como ˆ2 = 0 V (x1 , 0) > 0 para todo x1 = 0, se debe tener a > 0. Tomemos ahora un x ˆ2 ) > 0 para todo x1 . As´ı, el polinomio p(x1 ) = V (x1 , x ˆ2 ) fijo de modo que V (x1 , x no puede tener ra´ıces reales. En consecuencia, su discriminante 4(b2 − ac) debe ser negativo.

7.3.1.

Aplicaci´ on a sistemas lineales

Consideremos de nuevo el sistema lineal x˙ = Ax ,

(7.6)

con x ∈ Rn y A ∈ Mn (R). En esta sección mostraremos cómo se puede aplicar a los sistemas lineales la teor´ıa de Liapunov. Para ello ensayaremos la siguiente “posible”funci´ on de Liapunov definida mediante una forma cuadr´ atica V = xT P x ,

(7.7)

siendo P ∈ Mn (R) simétrica. La derivada orbital V˙ de V respecto del sistema lineal (7.6) viene dada por V˙ = x˙ T P x + xT P x˙ = xT AT P x + xT P Ax = xT (AT P + P A)x , de modo que, si definimos la matriz Q mediante la llamada ecuaci´ on matricial de Liapunov (7.8) AT P + P A = −Q , se tiene que

V˙ = −xT Qx .

(7.9)

7.3 Teor´ıa de Liapunov de la estabilidad

59

Es f´ acil ver que Q es una matriz simétrica puesto que P lo es. Esto es debido a que QT = −(AT P + P A)T = −(P T A + AT P T ) = −(P A + AT P ) = Q . Por lo tanto la expresi´ on de V˙ dada en (7.9) es una forma cuadr´ atica. Es habitual fijar V˙ , es decir, fijar Q y calcular V o equivalentemente calcular una soluci´ on P de la ecuación matricial de Liapunov (7.8). Los siguientes teoremas son una herramienta teórica con grandes aplicaciones en la teor´ıa de control ´ optimo como veremos en cap´ıtulos posteriores. Teorema 57 (Liapunov). La ecuaci´ on matricial de Liapunov (7.8) verifica lo siguiente: (i) La ecuaci´ on (7.8) tiene una u ńica soluci´ on P simétrica para cada Q simétrica fijada si y s´ olo si para cada par de valores propios λi y λj de A se verifica λi + λj = 0. (ii) Si A es una matriz de estabilidad, para cada Q definida positiva existe una u ńica soluci´ on de (7.8) que es definida positiva y viene dada de forma expl´ıcita por ∞

P =

exp(AT t)Q exp(At) dt .

(7.10)

0

Nota 58. La f´ ormula (7.10) tiene interés te´ orico puesto que en muchas ocasiones es m´ as simple resolver la ecuaci´ on matricial de Liapunov (7.8) directamenete que calcular mediante (7.10). Teorema 59. La matriz A es una matriz de estabilidad si y s´ olo si para cualquier matriz Q ∈ Mn (R) simétrica y definida positiva, la matriz P ∈ on de la ecuaci´ on matricial de Liapunov (7.8) es también definida Mn (R) soluci´ positiva. Demostraci´ on. Consta de dos partes: Si P y Q son ambas definidas positivas, entonces la V dada en (7.7) es una funci´ on de Liapunov y adem´ as V˙ es definida negativa de modo que, por el Teorema 55, el origen del sistema (7.6) es asintóticamente estable. Si Q es definida positiva y P es definida negativa o indefinida entonces en ambos casos la V definida en (7.7) puede tomar valores negativos en un entorno del origen. Entonces, el origen del sistema (7.6) es inestable.

7.3.2.

Estabilidad mediante linearizaci´ on

La teor´ıa de estabilidad desarrollada para sistemas lineales puede ser aplicada en algunos casos a sistemas no lineales mediante el concepto de linearizaci´ on.

60

Estabilidad

Consideremos de nuevo la ecuación no lineal (7.5) con el punto cr´ıtico x∗ ∈ Rn , es decir, (7.11) x˙ = f (x) , f (x∗ ) = 0 . Supongamos que f es suficientemente regular de modo que podemos realizar el desarrollo de Taylor de f (x) en un entorno de x∗ , es decir, f (x) = Df (x∗ )(x− on g(x) la que contiene a los términos de orden superior x∗ )+g(x), siendo la funci´ a los lineales en el desarrollo efectuado y ⎛ ∂f ⎞ ∂f1 ∗ ∗ 1 ∂x1 (x ) · · · ∂xn (x ) ⎜ ∂f2 ∗ ⎟ ∂f ⎜ ∂x1 (x ) · · · ∂xn2 (x∗ ) ⎟ ∗ ⎜ ⎟ ∈ Mn (R) , Df (x ) = ⎜ .. .. .. ⎟ ⎝ ⎠ . . . ∂fn ∗ ∂x1 (x )

···

∂fn ∗ ∂xn (x )

la matriz jacobiana de f calculada en el punto x∗ . Definici´ on 60. El sistema lineal x˙ = Df (x∗ )x ,

(7.12)

se llama la linearización del sistema no lineal x˙ = f (x) en el punto singular x∗ . Teorema 61 (Estabilidad por linearizaci´ on). Si el sistema lineal x˙ = Df (x∗ )x es asint´ oticamente estable o inestable, entonces el punto cr´ıtico x∗ del sistema no lineal x˙ = f (x) con f (x∗ ) = 0 tiene la misma propiedad de estabilidad. Demostraci´ on. Demostraremos sólo la primera parte del teorema relacionada con la estabilidad asint´ otica. Supondremos adem´ as que x∗ = 0. T Consideremos la funci´ on V = x P x donde la matriz P satisface la ecuación matricial de Liapunov Df T (0)P + P Df (0) = −Q , con Q una matriz arbitraria real simétrica y definida positiva. Supongamos que el sistema linearizado x˙ = Df (0)x es asint´ oticamente estable. Entonces P es definida positiva por el Teorema 59. Adem´ as, la derivada orbital de V respecto del sistema x˙ = f (x) = Df (0)x + g(x) es V˙

= x˙ T P x + xT P x˙ = [Df (0)x + g]T P x + xT P [Df (0)x + g] = [xT Df T (0) + g T ]P x + xT P [Df (0)x + g] = xT [Df T (0)P + P Df (0)]x + g T P x + xT P g = −xT Qx + g T P x + xT P g .

Por otra parte

xT P g = (P g)T x = g T P T x = g T P x ,

7.4 Estabilidad y control

61

donde en el primer paso se ha utilizado la propiedad conmutativa del producto escalar y en el u ´ltimo el hecho de que P es simétrica. De este modo se tiene que V˙ = −xT Qx + 2g T P x . Puesto que g comienza en términos de al menos grado 2, es claro que el segundo término 2g T P x de V˙ comienza con términos de al menos grado 3. Entonces, para x suficientemente cercano al origen domina el primer término −xT Qx en la expresi´ on de V˙ . Esto implica que, puesto que Q es una matriz definida positiva, para todo x suficientemente cercano al origen se tiene V˙ < 0. Aplicando ahora el Teorema 55, se concluye que el origen del sistema no lineal x˙ = f (x) es asint´ oticamente estable. Nota 62. Observar que el Teorema 61 no aporta ninguna informaci´ on sobre la estabilidad del punto cr´ıtico x∗ del sistema x˙ = f (x) cuando el sistema linoticamente estable. earizado x˙ = Df (x∗ )x es estable pero no es asint´

7.4.

Estabilidad y control

Veamos alguna de las relaciones existentes entre los conceptos de estabilidad y los de control. Consideremos de nuevo el sistema lineal controlado x˙ = Ax + Bu ,

(7.13)

con A ∈ Mn (R) y B ∈ Mn×m (R) matrices constantes, x(t) ∈ Rn y u(t) ∈ Rm . Adem´ as, supondremos que la salida y(t) del sistema es y = Cx ,

(7.14)

con C ∈ Mr×n (R).

7.4.1.

Estabilidad entrada–salida

Cuando un sistema es controlado mediante una entrada, es decir un vector de control u(t), es de gran utilidad pr´ actica definir un nuevo tipo de estabilidad. Definici´ on 63. El sistema x˙ = Ax + Bu, y = Cx, se dice que es entrada acotada–salida acotada estable (o bien b.i.b.o) si cualquiere control de entrada u(t) acotado produce una salida y(t) acotada. En otras palabras, independientemente del estado inicial x(t0 ), si u(t) < c1 para todo t ≥ t0 con c1 una constante entonces existe otra constante c2 tal que y(t) < c2 para todo t ≥ t0 . Nota 64. Es habitual en la literatura encontrar a los sistems con entrada acotada–salida acotada estables mediante el nombre b.i.b.o provenientes de la abreviatura del inglés bounded input–bounded output. Teorema 65. Si el sistema lineal x˙ = Ax es asint´ oticamente estable, entonces el sistema x˙ = Ax + Bu, y = Cx es b.i.b.o. estable.

62

Estabilidad

Demostraci´ on. En primer lugar, recordemos que la soluci´ on del sistema x˙ = un (3.21) por la f´ ormula Ax + Bu con condición inicial x(0) = x0 viene dada seg´ de variaci´ on de constantes

t exp(−Aτ )Bu(τ ) dτ . x(t) = exp(At) x0 + 0

Tomando normas en la ecuación y = Cx y utilizando sus propiedades se obtiene y(t) ≤ C x(t) t exp[A(t − τ )] Bu(τ ) dτ (7.15) . ≤ C exp(At)x0 + 0

Por otra parte, si A es una matriz de estabilidad, entonces existen constantes reales positivas K y α de modo que exp(At) ≤ K exp(−αt) ≤ K , ∀ t ≥ 0 .

(7.16)

Esta propiedad, para el caso particular en que todos los valores propios λk de A son distintos, se deduce de tomar normas en la f´ ormula de Sylvester (3.14), es decir, # # n n # # # # Zk exp(λk t)# ≤ Zk exp(λk t) exp(At) = # # # ≤ ≤

k=1 n

ˆ K

ˆ K

k=1 n

k=1 n

ˆ exp(λk t) = K

exp[ (λk )t]

k=1

ˆ exp[−αt] exp[−αt] = nK

k=1

ˆ se ˆ := máx1≤k≤n {Zk } y α := m´ın1≤k≤n { (λk )}. Tomando K = nK donde K deduce (7.16). Teniendo en cuenta que la entrada de control est´ a acotada, es decir u(t) < c1 , de las ecuaciones (7.15) y (7.16) se deduce la siguiente cota t exp[A(t − τ )] dτ y(t) ≤ C Kx0 + c1 B 0 t ≤ C Kx0 + c1 BK exp[−α(t − τ )] dτ 0 1 = KC x0 + c1 B [1 − exp(−αt)] α 1 ≤ KC x0 + c1 B para todo t ≥ 0 . α En definitiva, y(t) está acotada y por lo tanto el sistema es b.i.b.o estable.

7.4 Estabilidad y control

63

El rec´ıproco del Teorema 65 no es cierto en general. Sin embargo se dispone del siguiente resultado demostrado en la p´ agina 53 de [12]. Teorema 66. Si el sistema lineal x˙ = Ax + Bu, y = Cx es completamente controlable, completamente observable y b.i.b.o. estable, entonces el sistema x˙ = Ax es asint´ oticamente estable.

7.4.2.

Estabilizaci´ on por realimentaci´ on lineal

Supongamos que el sistema de lazo abierto x˙ = Ax es inestable. Dicho de otro modo, existe al menos un valor propio de la matriz A con parte real positiva (ver Teorema 49). En estos casos, un problema de gran aplicaci´ on consiste en averiguar si se puede aplicar alg´ un control u(t) de forma que el sistema x˙ = Ax + Bu se estabilice. Sabemos que, si el sistema x˙ = Ax + Bu es completamente controlable, entonces existe una realimentación lineal u = Kx de modo que el sistema de lazo cerrado x˙ = (A + BK)x es asint´ oticamente estable, ver Teorema 36. De hecho existen infinitas matrices de ganancia K diferentes que consiguen este objetivo. Sin embargo, cuando el sistema x˙ = Ax + Bu no es completamente controlable, entonces podemos definir una propiedad m´ as débil que esta. En concreto, se dice que el sistema x˙ = Ax + Bu es estabilizable si existe una matriz real K tal que el sistema de lazo cerrado x˙ = (A + BK)x es asint´ oticamente estable.

7.4.3.

Linealizaci´ on de sistemas de control no lineales

Consideremos un sistema de control no lineal x˙ = f (x, u) ,

(7.17)

siendo x(t) ∈ Rn las variables de estado y u(t) ∈ Rm las variables de control. Aqu´ı no estamos suponiendo que la funci´ on vectorial f sea lineal. De hecho, cuando se modeliza matemáticamente un sistema f´ısico de control, habitualmente estos son no lineales. Sin embargo, existen algunos sistemas que son o bien débilmente no lineales o bien sus caracter´ısticas no lineales se dan bajo ciertas restricciones de operación. El hecho de que podamos aproximar este tipo de sistemas no lineales mediante sistemas lineales nos brinda poder atacar el sistema con toda la artiller´ıa de métodos anal´ıticos y exactos propios de este tipo de sistemas. on un proceso de linearizaci´ on Suponiendo f ∈ C 1 , describiremos a continuaci´ del sistema (7.17) en un entorno de una soluci´ on x∗ (t) correspondiente al control u∗ (t). Usualmente, x∗ es un punto singular del sistema (7.17) correspondiente a un cierto control constante u∗ , es decir, f (x∗ , u∗ ) = 0 aunque no es necesario que esto ocurra.

64

Estabilidad Realizando una expansi´ on de Taylor de la funci´ on f (x, u) = (f1 (x, u), . . . , fn (x, u))

en un entorno del punto (x∗ , u∗ ) = (x∗1 , . . . , x∗n , u∗1 , . . . , u∗m ) y truncando a primer orden se tiene n m ∂fi ∗ ∗ ∂fi ∗ ∗ x˙ i ≈ fi (x∗ , u∗ )+ (x , u )(xj −x∗j )+ (x , u )(uj −u∗j ) , i = 1, . . . , n . ∂x ∂u j j j=1 j=1

Puesto que se pretende estudiar el sistema en un entorno de (x∗ , u∗ ), definimos las variables (7.18) ∆xi = xi − x∗i , ∆ui = ui − u∗i . Como ∆x˙ i = x˙ i − x˙ ∗i = x˙ i − fi (x∗ , u∗ ), se tiene ∆x˙ i =

n m ∂fi ∗ ∗ ∂fi ∗ ∗ (x , u )∆xj + (x , u )∆uj , i = 1, . . . , n . ∂xj ∂uj j=1 j=1

Por supuesto, este sistema lineal se puede escribir en la forma matricial ˆ ˆ ∆x˙ = A∆x + B∆u ,

(7.19)

siendo ∆x = (∆x1 , . . . , ∆xn )t ∈ Rn , ∆u = (∆u1 , . . . , ∆um )t ∈ Rm , ⎞ ⎛ ∂f ∂f1 ∗ ∗ ∗ ∗ 1 ∂x1 (x , u ) · · · ∂xn (x , u ) ⎟ ⎜ ∂f2 ∗ ∗ ∂f ⎜ ∂x1 (x , u ) · · · ∂xn2 (x∗ , u∗ ) ⎟ ⎟ ∈ Mn (R) , Aˆ = ⎜ .. .. .. ⎟ ⎜ ⎠ ⎝ . . . ∂fn ∂fn ∗ ∗ ∗ ∗ (x , u ) · · · (x , u ) ∂x ∂xn ⎞ ⎛ ∂f 1 ∂f1 ∗ ∗ ∗ ∗ 1 (x , u ) · · · ∂u1 ∂um (x , u ) ⎟ ⎜ ∂f2 ∗ ∗ ∂f2 (x∗ , u∗ ) ⎟ ⎜ ∂u1 (x , u ) · · · ∂um ˆ ⎟ ∈ Mn×m (R) . ⎜ B = ⎜ .. .. .. ⎟ ⎠ ⎝ . . . ∂fn ∗ ∗ ∂u1 (x , u )

7.5.

···

∂fn ∗ ∗ ∂um (x , u )

Problemas resueltos

1. Consideremos un péndulo de longitud inmerso en un fluido que le proporciona un rozamiento proporcional a la velocidad. Es f´ acil ver que la ecuaci´ on del movimiento del péndulo viene dada por θ¨ + 2cθ˙ + Ω2 sin θ = 0 , donde c > 0 es un coeficiente de rozamiento y Ω2 = /g con g la gravedad supuesta constante. Estudiar la estabilidad de sus puntos singulares.


65

Soluci´ on. Usando el cambio de variables habitual θ˙ = w se tiene el sistema θ˙ = w , w˙ = −Ω2 sin θ − 2cw .

(7.20)

Los puntos singulares del sistema en el plano de fases son (θ, w) = (nπ, 0), con n ∈ Z. Estudiemos a continuaci´ on la estabilidad del origen del sistema (7.20) por linearizaci´ on. La matriz Jacobiana de (7.20) es 0 1 , Df (θ, w) = −Ω2 cos θ −2c de modo que, calculada en el origen es 0 Df (0, 0) = −Ω2

1 −2c

.

√ Los valores propios asociados a Df (0, 0) son λ± = −c ± c2 − Ω2 . Debido a que c > 0 y Ω > 0 se tiene que (λ± ) < 0 (independientemente de si c ≥ Ω o no). Por lo tanto, utilizando una combinaci´ on del Teorema 49 y del Teorema 61 se concluye que el origen (θ, w) = (0, 0) es asint´ oticamente estable para el sistema (7.20). Es f´ acil demostrar que si n es par, las ecuaciones que se obtienen del estudio de la estabilidad por linearizaci´ on de los puntos singulares (θ, w) = (nπ, 0) son las mismas que para el origen. De este modo, el análisis realizado para el origen se aplica de igual forma aqu´ı. En definitiva se tiene que los puntos singulares (θ, w) = (nπ, 0) son asint´ oticamente estables para el sistema (7.20). Por el contrario, si calculamos la matriz Jacobiana de (7.20) en los puntos singulares (θ, w) = (nπ, 0) con n impar se obtiene 0 1 . Df (nπ, 0) = Ω2 −2c Los √ valores propios asociados a Df (nπ, 0) con n impar son λ± = −c ± c2 + Ω2 . Se tiene ahora que (λ+ ) > 0 y (λ− ) < 0. Aplicando de nuevo el Teorema 49 y el Teorema 61 se concluye que los puntos singulares (θ, w) = (nπ, 0) con n impar son inestable para el sistema (7.20). 2. Consideremos el sistema no lineal siguiente obtenido mediante una perturbaci´ on de un oscilador arm´ onico x˙ = y + αx(x2 + y 2 ) , y˙ = −x + αy(x2 + y 2 ) , con α ∈ R\{0} un par´ ametro. (i) ¿Es posible estudiar la estabilidad del origen del sistema por linearizaci´ on?

66

Estabilidad (ii) Sea (x(t), y(t)) una soluci´ on del sistema. Utilizar la funci´ on (x(t), y(t)2 para estudiar la estabilidad del origen del sistema. Soluci´ on. (i) La matriz Jacobiana asociada al sistema y calculada en el origen es 0 1 Df (0, 0) = . −1 0 Los valores propios asociados a Df (0, 0) son λ± = ±i y por lo tanto imaginarios puros. Puesto que (λ± ) = 0, no se puede aplicar el Teorema 61 y por lo tanto la linearizaci´ on no es suficiente para analizar la estabilidad del origen del sistema. (ii) Sea V (t) = (x(t), y(t)2 = x2 (t) + y 2 (t). Calculamos la derivada orbital de V y obtenemos V˙ = 2xx+2y ˙ y˙ = 2x[y+αx(x2 +y 2 )]+2y[−x+αy(x2 +y 2 )] = 2α(x2 +y 2 )2 . Está claro que, excepto en el origen, el signo de V˙ coincide con el de α. De este modo, teniendo en cuenta el significado geométrico de la funci´ on V , concluimos que Si α < 0 entonces (x(t), y(t)) → (0, 0) cuando t → ∞ de manera monótona. En particular, el origen del sistema es asint´ oticamente estable. Notar que, aplicando el apartado (ii) del Teorema 55, se puede deducir el resultado teniendo en cuenta que V es una funci´ on de Liapunov y V˙ es definida negativa. Si α > 0 entonces todas las soluciones (x(t), y(t)) del sistema, excepto el origen, escapan hacia el infinito cuando t → ∞. En particular, el origen del sistema es inestable. 3. Consideremos un péndulo de longitud y masa m. Se sabe que, despreciando los efectos del rozamiento y tomando la gravedad g constante, la ecuaci´ on del movimiento viene dada por g θ¨ + sin θ = 0 , siendo θ el ´ angulo que forma el péndulo con la vertical. ˙ = (0, 0) del sistema (i) ¿Es posible estudiar la estabilidad del origen (θ, θ) por linearizaci´ on? ˙ = θ2 + θ˙2 como funci´ (ii) ¿Es posible utilizar la funci´ on V (θ, θ) on de Liapunov en un entorno del origen? ˙ = 1 m2 θ˙2 + mg(1 − (iii) ¿Es posible utilizar la funci´ on energ´ıa E(θ, θ) 2 cos θ) como funci´ on de Liapunov en un entorno del origen?


67

Soluci´ on. (i) Usando el cambio de variables habitual θ˙ = w se tiene el sistema g θ˙ = w , w˙ = − sin θ . Los puntos singulares del sistema en el plano de fases son (θ, w) = (nπ, 0), con n ∈ Z. La matriz Jacobiana del sistema es Df (θ, w) =

0 − g cos θ

de modo que, calculada en el origen es 0 Df (0, 0) = − g

1 0

1 0

,

.

Los valores propios λ± asociados a Df (0, 0) son imaginarios puros, de modo que (λ± ) = 0 y no se puede concluir nada sobre la estabilidad del origen por el método de linearización. (ii) Sea V (θ, w) = θ2 + w2 . La derivada orbital de V a lo largo de las órbitas del sistema es ! g V˙ = 2w θ − sin θ . Observemos que V˙ (θ, w) no tiene signo definido en ning´ un entorno del origen. Una forma de verlo es la siguiente. El desarrollo de Taylor de sin θ entorno del origen es sin θ = θ + O(θ3 ), de modo que ! $ g% + O(θ3 ) . V˙ (θ, w) = 2w θ 1 −

68

Estabilidad De esta expresión se obtiene que el signo de sign [1 − g/] si ˙ sign [V ] = sign [−(1 − g/)] si

V˙ viene dado por w > 0, 0 < θ << 1 , w < 0, 0 < θ << 1 ,

Se concluye que no es posible tomar V como función de Liapunov en un entorno del (θ, w) = (0, 0). ˙ = 1 m2 θ˙2 +mg(1−cos θ) es definida (iii) Observar que la funci´ on E(θ, θ) 2 positiva en un entorno suficientemente peque¯ no del origen puesto que , m y g son constantes positivas. Además, la derivada orbital E˙ ≡ 0. Entonces, E es una función de Liapunov para el sistema y, aplicando la parte (i) del Teorema 55, el origen es estable. 4. Considerar un oscilador amortiguado, es decir, un sistema del tipo m¨ x + cx˙ + kx = 0 , con c > 0 y k > 0 las constantes de rozamiento y el´ astica respectivamente. Demostrar, utilizando la energ´ıa mec´ anica total como funci´ on de Liapunov, que el punto cr´ıtico (x, x) ˙ = (0, 0) es estable.

Soluci´ on. En primer lugar escribimos el sistema en el plano de fases c k x− y . m m La energ´ıa mecánica E del sistema es la suma de la energ´ıa potencial elástica y de la energ´ıa cinética, es decir, x˙ = y , y˙ = −

1 1 my 2 + kx2 . 2 2 Puesto que m y k son constantes positivas, es obvio que E(x, y) > 0 para todo (x, y) = 0 y que E(0, 0) = 0 de modo que la funci´ on E es definida positiva. Además, su derivada orbital es E(x, y) =

E˙

∂E ∂E y˙ x˙ + ∂y ∂x c k = kxy + my − x − y = −cy 2 ≤ 0 , m m =


69

por lo tanto E˙ es semidefinida positiva. En definitiva E(x, y) es una funci´ on de Liapunov para el sistema en un entorno del punto singular (x, y) = (0, 0). Aplicando el apartado (i) del Teorema 55 se concluye que el punto singular (x, y) = (0, 0) es estable1 . 5. Estudiar la estabilidad de los puntos singulares del sistema x˙ = −2xy , y˙ = x2 − y 3 , utilizando una funci´ on de Liapunov del tipo V (x, y) = ax2m + by 2n con constantes adecuadas a, b, m, n. Soluci´ on. Est´ a claro que la u ńica soluci´ on del sistema de ecuaciones 0 = −2xy , 0 = x2 − y 3 , es (x, y) = (0, 0). Por lo tanto el origen es el u ńico punto singular del sistema. La derivada orbital de V es V˙

∂V ∂V x˙ + y˙ = 2max2m−1 (−2xy) + 2nby 2n−1 (x2 − y 3 ) ∂x ∂y = [−4max2m y + 2nbx2 y 2n−1 ] − 2nby 2n+2 .

=

Elejimos los valores de los par´ ametros de manera que la expresión encerrada entre corchetes se anule. Por ejemplo, mediante simple inspección, una elección es m = n = a = 1, b = 2. Se tiene entonces V (x, y) = x2 + 2y 2 que es definida positiva y V˙ (x, y) = −4y 4 que es semidefinida negativa. Concluimos que V es una funci´ on de Liapunov para el sistema y, aplicando el apartado (i) del Teorema 55 se concluye que el punto singular (x, y) = (0, 0) es estable. 6. Considerar el sistema lineal x˙ = Ax con x ∈ R2 y −1 2 A= . 0 −3 (i) Demostrar que el sistema es asint´ oticamente estable. (ii) Hallar una funci´ on de Liapunov cuadr´ atica par el origen del sistema. Soluci´ on. (i) Aplicando la parte (i) del Teorema 49 se concluye que el sistema es asintóticamente estable puesto que los valores propios de la matriz A son λ1 = −1 y λ2 = −3 y por lo tanto todos tienen parte real negativa. (ii) Consideremos la siguiente posible funci´ on de Liapunov cuadr´ atica V (x) = xT P x, siendo P ∈ M2 (R) una matriz simétrica. Sabemos que, 1 De hecho es f´ acil demostrar que el origen es asint´ oticamente estable aunque la funci´ on de Liapunov E no nos permite detectarlo.

70

Estabilidad si la derivada orbital de V respecto del sistema x˙ = Ax es V˙ = −xT Qx entonces se verifica (7.8) que es la llamada ecuación matricial de Liapunov AT P + P A = −Q . Observar que la matriz A es una matriz de estabilidad, ver el apartado anterior donde se comprueba que todos sus valores propios tienen parte real negativa. Entonces, por el apartado (ii) del Teorema 57, se sabe que, para cada Q definida positiva existe una u ńica soluci´ on de la ecuación matricial de Liapunov que es definida positiva. Sea pues a 0 Q= , 0 b con a, b ∈ R+ \{0}. Obviamente Q es una matriz simétrica y definida positiva. Sea P = (pij ) ∈ M2 (R) simétrica. Entonces, la ecuación matricial de Liapunov queda de la forma p11 p12 −1 2 a 0 −1 0 p11 p12 + =− . p12 p22 p12 p22 0 −3 0 b 2 −3 De aqu´ı se tiene el sistema lineal para las incógnitas pij siguiente 2p11 = a , p11 − 2p12 = 0 , 2p12 − 3p22 = −b , cuya solución es p11 = a/2, p12 = a/4 y p22 = (a + b)/6. En definitiva obtenemos 1 a a/2 P = . a/2 (a + b)/3 2 Nótese que P no sólo es simétrica, también es definida positiva como puede comprobarse aplicando, por ejemplo, el criterio de Sylvester y comprobando que todos los determinantes menores principales de P son positivos. Por supuesto existe otra forma de calcular P basada en la expresi´ on (7.10), es decir ∞ exp(AT t)Q exp(At) dt .

P = 0

Un simple cálculo muestra que ∞ ae−2t ae−t [e−t − e−3t ] P = dt ae−t [e−t − e−3t ] be−6t + [e−t − e−3t ]2 0 1 a a/2 = . a/2 (a + b)/3 2 En definitiva, obtenemos la siguiente familia biparamétrica de funciones de Liapunov cuadr´ aticas para el sistema 1 a a/2 x1 V (x) = xT P x = (x1 , x2 ) a/2 (a + b)/3 x2 2 1 = [3a(x21 + x1 x2 ) + (a + b)x22 ] , 6


71

siendo a, b constantes reales positivas. 7. Considerar el sistema x˙ = Ax, con x ∈ R3 y ⎛ ⎞ −1 −1 0 A = ⎝ −1 −1 2a ⎠ . b 0 −b (i) Averiguar para qué valores de los par´ ametros reales a y b el sistema x˙ = Ax es asint´ oticamente estable. (ii) Supongamos que a = 1 y b = 2. ¿Es b.i.b.o estable el sistema controlado x˙ = Ax + Bu, y = Cx para cualquier matriz constante B, C adecuadas y vector control u(t)? (iii) Supongamos que a = 1, b = −3. Averiguar si el sistema x˙ = Ax+Bu, y = Cx es b.i.b.o estable donde ⎛ ⎞ 1 0 B = ⎝ 2 1 ⎠ , C = (1, 5, 2) . 7 3 Soluci´ on. (i) El polinomio caracter´ıstico asociado a la matriz A es k(λ) = det(A − λI3 ) = λ3 + (2 + b)λ + 2bλ + 2ab . Teniendo en cuenta la Proposici´ on 51, todas las ra´ıces del polinomio k(λ) tienen parte real negativa si y s´ olo si b + 2 > 0 , b > 0 , a > 0 , b(2 + b) > ab , o, equivalentemente b + 2 > 0 , b(b + 2 − a) > 0 , ab > 0 . (ii) Seg´ un el apartado anterior, para los valores de los par´ ametros a = 1 y b = 2 se verifica que todos los valores propios de la matriz A tienen parte real negativa, es decir el sistema x˙ = Ax es asint´ oticamente estable. Entonces, por el Teorema 65, se concluye que el sistema x˙ = Ax + Bu, y = Cx es b.i.b.o estable. (iii) Con la nueva elección de par´ ametros a1, b = −3 el sistema x˙ = Ax no es asint´ oticamente estable, ver apartado (i). En consecuencia no se puede aplicar el Teorema 65 para estudiar si el sistema x˙ = Ax + Bu, y = Cx es b.i.b.o estable. Estudiemos la controlabilidad y la observabilidad del sistema x˙ = Ax+Bu, y = Cx.

72

Estabilidad La matriz de controlabilidad es ⎛ 1 0 U = [B, AB, A2 B] = ⎝ 2 1 7 3

−3 11 28

−1 7 5 18 14 63

⎞ 3 9 ⎠ ∈ M3×6 (R) 30

y es fácil comprobar que su rango es 3 por ejemplo calculando el determinante de orden tres formado por sus primeras tres columnas 1 0 −3 2 1 11 = −2 = 0 . 7 3 28 Puesto que rang U = 3 es maximal se tiene, aplicando el Teorema 4.2, que el sistema es completamente controlable. La matriz de observabilidad V es ⎡ ⎤ ⎛ ⎞ C 1 5 2 V = ⎣ CA ⎦ = ⎝ −2 −6 16 ⎠ . CA2 −30 18 36 Como det V = −1536 = 0 se tiene que V tiene rango máximo, es decir, rang V = 3. En conclusi´ on, aplicando el Teorema 33, el sistema es completamente observable. Finalmente, puesto que el sistema x˙ = Ax + Bu, y = Cx es completamente controlable y observable pero x˙ = Ax no es asint´ oticamete estable, aplicando el Teorema 66, se concluye que el sistema x˙ = Ax + Bu, y = Cx no es b.i.b.o estable. 8. Consideremos un péndulo de longitud sometido a una gravedad g constante. Supongamos que se pretende controlar el péndulo aplic´ andole una fuerza por unidad de masa u de modo que la ecuaci´ on del movimiento es θ¨ = −g sin θ + u. Suponiendo un péndulo de peque¯ nas oscilaciones de modo que sin θ ≈ θ y que se utiliza una realimentaci´ on del tipo u = −kθ, averiguar si es posible estabilizar θ a zero. Soluci´ on. Puesto que sin θ ≈ θ se tiene que las ecuaciones del movimiento son aproximadamente θ¨ = −gθ + u. En el plano de las fases θ − w, siendo w = θ˙ la velocidad angular, las ecuaciones del movimiento se escriben de la forma x˙ = Ax + bu, siendo x = (θ, w)T y A=

0 1 −g 0

, b=

0 1

.

Puesto que la realimentaci´ on es del tipo u = −kθ, se puede escribir de la forma u = Kx, siendo K = (−k, 0).


73

El sistema en lazo cerrado es x˙ = (A + bK)x, siendo la matriz del sistema 0 1 0 0 1 A + bK = + (−k, 0) = , −g 0 1 −g − k 0 √ que tiene valores propios λ± = ± −g − k. Para estabilizar θ a zero, se tiene que elegir k de manera que el sistema lineal en lazo cerrado sea asint´ oticamente estable. Dicho de otro modo, se deber´ıa tomar cualquier k tal que (λ± ) < 0. Por supuesto esto es imposible ya que si −g − k > 0 entonces (λ+ ) > 0; si −g − k ≤ 0 entonces (λ± ) = 0. Hemos demostrado pues que con una realimentaci´ on del tipo u = −kθ no es posible estabilizar θ a zero en un péndulo de peque¯ nas oscilaciones. 9. Consideremos un péndulo invertido montado en un carrito que se puede desplazar sobre unos railes horizontalmente. El carro tiene masa M y est´ a sujeto a una pared por un muelle de constante el´ astica k2 . El péndulo se ha montado sobre el carrito mediante un muelle de constante el´ astica k1 . La longitud del péndulo es 2 y su masa total m se supone distribuida de forma homogénea a lo largo de su varilla. Sea u la fuerza de control horizontal que se puede ejercer sobre el carro. Denotemos por z(t) la coordenada del centro de gravedad del carro respecto de su posici´ on de equilibrio y sea θ(t) el ´ angulo formado por el péndulo y la vertical que es la variable que puede ser medida. Aplicando las leyes de la mec´ anica cl´ asica es f´ acil ver que las ecuaciones del movimiento del sistema son: (M + m)¨ z + k2 z + m cos θθ¨ = mθ˙2 sin θ + u , 4 m cos θ¨ z + m2 θ¨ = mg sin θ − k1 θ . 3 Supondremos por simplicidad que M = m = = 1. (i) Reescribir las ecuaciones del movimiento en la forma x˙ = f (x, u), ˙ T. siendo las variables de estado x = (z, z, ˙ θ, θ) (ii) Si no act´ ua ninguna fuerza exterior sobre el carro y éste est´ a en reposo en su posici´ on natural con el péndulo también en reposo y en posici´ on vertical est´ a el sistema en equilibrio? (iii) Linealizar el sistema entorno del anterior punto cr´ıtico. Soluci´ on. (i) Es f´ acil ver que, las ecuaciones del movimiento se reescriben de la forma x˙ = f (x, u), siendo las variables de estado ˙ T ∈ R4 , x = (x1 , x2 , x3 , x4 ) = (z, z, ˙ θ, θ) y la funci´ on f (x, u) = (f1 (x, u), f2 (x, u), f3 (x, u), f4 (x, u))T definida por ⎞ ⎛ x2

⎜ Λ[4k2 x1 − 3k1 x3 cos x3 − 4x24 sin x3 + 3g cos x3 sin x3 − 4u] f (x, u) = ⎜ ⎝ x4

Λ[−3k2 x1 cos x3 + 6k1 x3 − 6g sin x3 + 3x24 cos x3 sin x3 + 3u cos x3 ]

⎟ ⎟ , ⎠

74

Estabilidad siendo Λ = 1/[3 cos3 x3 − 8]. (ii) Bajo las condiciones del enunciado el sistema se encuentra en el punto x∗ = (0, 0, 0, 0) con u∗ = 0. Notemos que f (x∗ , u∗ ) = (0, 0, 0, 0)T , de modo que el sistema se encuentra en un estado de equilibrio, es decir un punto cr´ıtico. (iii) El sistema linealizado entorno del punto (x∗ , u∗ ) = (0, 0) ∈ R4 × R viene dado por ˆ + ˆbu , x˙ = Ax siendo ⎛ ⎜ ⎜ Aˆ = ⎜ ⎜ ⎝ ⎛ ⎜ ˆb = ⎜ ⎜ ⎝

∂f1 ∗ ∗ ∂x1 (x , u ) ∂f2 ∗ ∗ ∂x1 (x , u )

∂f1 · · · ∂x (x∗ , u∗ ) 4 ∂f2 · · · ∂x4 (x∗ , u∗ ) .. .. .. . . . ∂f4 ∂f4 ∗ ∗ ∗ (x , u ) · · · (x , u∗ ) ∂x1 ∂x4 ⎞ ∂f1 ∗ ∗ ∂u (x , u ) ∂f2 ∗ ∗ ⎟ ∂u (x , u ) ⎟ ⎟ ∈ R4 . .. ⎠ .

⎞ ⎟ ⎟ ⎟ ∈ M4 (R) , ⎟ ⎠

∂f4 ∗ ∗ ∂u (x , u )

En concreto se obtiene ⎛ Aˆ =

⎜ ⎜ ⎝ ⎛

ˆb =

⎜ ⎜ ⎝

⎞ 0 1 0 0 −4k2 /5 0 3(k1 − g)/5 0 ⎟ ⎟ ∈ M4 (R) , 0 0 0 1 ⎠ 3k2 /5 0 6(g − k1 )/5 0 ⎞ 0 4/5 ⎟ ⎟ ∈ R4 . 0 ⎠ −3/5

Cap´ıtulo 8

C´ alculo de Variaciones 8.1.

Un ejemplo: braquistocrona

El nacimiento del c´ alculo de variaciones se atribuye al problema de la braquistocrona1 o curva de descenso más rápida, propuesto y resuelto por el eminente matemático suizo Johann Bernoulli en 1696, aunque también dieron soluciones otros contemporáneos suyos como Jacob Bernoulli, Leibniz y Newton. Dicho problema consist´ıa en determinar cu´ al, de entre todas las posibles trayectorias, era la que llevaba a una part´ıcula sin rozamiento en el menor tiempo posible, desde un punto A hasta otro punto B en el mismo plano vertical, sujeta sólo a la acción de la gravedad. Para resolver este problema, se deben de considerar todas las infinitas curvas que pasan por A y B, de modo que, a cada una de ellas se le asigna un tiempo (el invertido por el punto material en descender desde A hasta B) y nos hemos de quedar con la curva con menor tiempo asociado. Veamos cómo formularlo matemáticamente. Consideremos el plano vertical de coordenadas cartesianas x e y que contiene a los puntos A y B. Obviamente, la curva soluci´ on debe estar contenida en dicho plano. Tomemos el punto A como el origen de coordenadas, el eje x de abscisas con sentido positivo hacia donde apunta la gravedad y, sea B = (x1 , y1 ) con x1 > 0 e y1 ≥ 0. Consideremos una curva en dicho plano arbitraria y descrita por la gr´ afica y = y(x) con 0 ≤ x ≤ x1 . Supondremos que la curva es suave, es decir, es de clase C 1 en su dominio. Adem´ as impondremos que los puntos inicial y final de dicha curva sean A y B respectivamente, es decir, y(0) = 0 e y(x1 ) = y1 . La masa m, en su descenso, conserva su energ´ıa mecánica puesto que no hay rozamiento. Supondremos que parte con velocidad inicial nula del punto A. Entonces, la energ´ıa mecánica en A vale EA = mghA , siendo EA > 0 y hA la altura a la que se encuentra el punto A. Conservando la energ´ıa mecánica en cualquier otro punto de la trayectoria de la part´ıcula se tiene 1 mv 2 + mgh = mghA , 2 1 El

nombre braquistocrona proviene del griego: braquistos m´ınimo y chronos tiempo.

75

76

C´ alculo de Variaciones

de modo que v 2 = 2g(hA − h). Con el sistema de referencia tomado está claro que x = hA − h, de modo que, la ecuación diferencial que rige el movimiento de la masa m es ds = 2gx , v := dt siendo v el módulo de la velocidad que tiene la part´ıcula cuando se encuentra en un punto con coordenada (x, y(x)) y ds el elemento de longitud de la trayectoria. Sabemos que, ds = 1 + y (x)2 dx, de modo que se obtiene 1 + y (x)2 dx √ . dt = 2gx Integrado esta expresión se obtiene el tiempo total T que tarda la masa m en recorrer la curva y = y(x) desde el punto A hasta el punto B. En concreto se tiene que x1 & 1 1 + y (x)2 dx . T = J(y) = √ x 2g 0 Hallar la braquistocrona consiste en obtener, de entre todas las funciones y(x) definidas para 0 ≤ x ≤ x1 , con y(0) = 0 e y(x1 ) = y1 , aquella que corresponda al menor valor del funcional J(y). Cabe destacar aqu´ı que, la soluci´ on de este problema viene dado por la famosa curva llamada cicloide, cuya ecuación paramétrica es x(t) =

1 1 (1 − cos t) , y(t) = 2 (t − sin t) , k2 k

siendo k una constante adecuada, ver el Problema 1 de este cap´ıtulo. -x 1

2

3

4

5

6

y

-0.5

-1

-1.5

-2

Figura 8.1: Gr´ afica de la cicloide invertida con k = 1 para t ∈ [0, 2π].

8.2 Introducci´ on

8.2.

77

Introducci´ on

Sea E un espacio de funciones x(t) ∈ Rn de clase C 1 [ta , tb ]. Consideremos el funcional J : E → R, definido de la forma

tb

F (x(t), x(t); ˙ t) dt ,

J(x) =

(8.1)

ta

con F ∈ C 2 respecto de sus argumentos. El problema fundamental del c´ alculo de variaciones consiste en hallar una funci´ on vectorial x(t) ∈ E que sea un extremal (máximo o m´ınimo) del funcional J(x) definido en (8.1). Por ejemplo, la funci´ on x∗ (t) minimiza al funcional J si se ∗ verifica J(x ) ≤ J(x) para todo x ∈ E. Por supuesto, el principal problema con el que nos encontramos y que lo diferencia radicalmente de las técnicas utilizadas en los problemas clásicos de hallar máximos y m´ınimos de funciones de varias variables es el hecho de que la dimensión del espacio vectorial E es infinita. Una forma intuitiva de ver esto consisteen suponer que se puede representar x(t) ∞ mediante la serie de Taylor x(t) = i=0 ai ti , de modo que el funcional J(x) se puede interpretar como una funci´ on de infinitas variables J(a0 , a1 , . . . , aj , . . .).

8.3.

Ecuaciones de Euler–Lagrange

En esta sección hallaremos las llamadas ecuaciones de Euler–Lagrange que son condiciones necesarias que ha de verificar una funci´ on x∗ (t) de clase C 2 para que dé un extremo del funcional J. La existencia, unicidad y condiciones suficientes de extremo se escapan a los objetivos de estas notas. Ver por ejemplo [4]. Además, en muchas ocasiones se sabe diferenciar entre máximo y m´ınimo por simples argumentos f´ısicos.

8.3.1.

Lema fundamental del c´ alculo de variaciones

En primer lugar damos un lema, conocido como el Lema Fundamental del C´ alculo de Variaciones, que necesitaremos para obtener posteriormente las ecuaciones de Euler–Lagrange. on Lema 67. Sean ta < tb constantes reales fijas y f (t) ∈ C[ta , tb ] una funci´ dada. Si tb η(t)f (t) dt = 0 , (8.2) ta

para toda η ∈ C 1 [ta , tb ] verificando η(ta ) = η(tb ) = 0, entonces f (t) = 0 para todo t ∈ [ta , tb ] .

(8.3)

78


Demostraci´ on. La demostración se basa en mostrar la existencia de al menos una funci´ on η para la cual (8.2) no se satisface cuando (8.3) no se cumple. Supongamos pues que (8.3) no se verifica, es decir, existe un t∗ ∈ (ta , tb ) para el cual f (t∗ ) = 0. Podemos suponer sin pérdida de generalidad que f (t∗ ) > 0. Como f es una funci´ on continua, debe existir un entorno de t∗ , digamos el ∗ ∗ intervalo [t1 , t2 ] en el cual f (t) > 0 para toda t ∈ [t∗1 , t∗2 ]. Pero entonces (8.2) no se puede verificar para todo η(t) admisible. Por ejemplo, consideremos la funci´ on admisible ⎧ si t ∈ [ta , t∗1 ] , ⎨ 0 ∗ 2 ∗ 2 (t − t1 ) (t − t2 ) si t ∈ [t∗1 , t∗2 ] , η(t) = ⎩ 0 si t ∈ [t∗2 , tb ] . Obviamente η ∈ C 1 [ta , tb ] verificando η(ta ) = η(tb ) = 0. Sin embargo, la integral (8.2) es tb t∗2 η(t)f (t) dt = (t − t∗1 )2 (t − t∗2 )2 f (t) > 0 , t∗ 1

ta

puesto que f (t) > 0 para toda t ∈ [t∗1 , t∗2 ].

8.3.2.

Ecuaciones de Euler–Lagrange

Sea x∗ (t) ∈ E una funci´ on que extremiza a J(x). Consideremos la siguiente familia uniparamétrica de funciones x(t) = x∗ (t) + ξ(t) , con ∈ R un par´ ametro y ξ ∈ C 1 [ta , tb ]. De este modo x(t) ∈ E. La variaci´ on de J sobre los miembros de esta familia viene dada por dJ = d

tb

ta

d ˙ F (x∗ (t) + ξ(t), x˙ ∗ (t) + ξ(t); t) dt . d

Obviamente, una condici´ on necesaria para que x∗ (t) sea un extremo del funcional J es que dJ =0. d =0 Como dx(t) ˙ dx(t) ˙ , = ξ(t) , = ξ(t) d d se tiene que, aplicando la regla de la cadena, dF d

˙ = ∇x F (x∗ (t) + ξ(t), x˙ ∗ (t) + ξ(t); t) ξ(t) ˙ ˙ , +∇x˙ F (x∗ (t) + ξ(t), x˙ ∗ (t) + ξ(t); t) ξ(t)

8.3 Ecuaciones de Euler–Lagrange

79

donde el punto indica la operaci´ on producto escalar ordinario en Rn . Sustituyendo esta expresión en el integrando de la expresi´ on dJ/d y realizando la integraci´ on por partes

tb

ta

!tb ˙ dt = ∇x˙ F ξ(t) − ∇x˙ F ξ(t) ta

tb

ta

d(∇x˙ F ) ξ(t) dt , dt

se obtiene que

dJ d

=0

= ∇x˙ F ξ(t)

!tb ta

tb

+ ta

∇x F −

d(∇x˙ F ) dt

ξ(t) dt = 0 , (8.4)

donde los gradientes están evaluados sobre la curva x∗ (t). Se tiene a partir de ahora dos casos diferenciados: (i) Los puntos inicial y final son fijos: Supongamos que el conjunto E donde buscamos soluciones sea el espacio de funciones x(t) de clase C 1 [ta , tb ] que adem´ as toman valores fijos en los extremos, es decir, x(ta ) = xa , x(tb ) = xb con xa , xb ∈ Rn fijos. Por supuesto en este caso se debe tener ξ(ta ) = ξ(tb ) = 0 de modo que (8.4) adopta la forma

tb

ta

∇x F −

d(∇x˙ F ) dt

ξ(t) dt = 0 ,

para toda ξ(t) ∈ C 1 [ta , tb ] verificando ξ(ta ) = ξ(tb ) = 0. Entonces, por el Lema Fundamental del Cálculo de Variaciones, x∗ (t) debe verificar la llamada ecuaci´ on de Euler–Lagrange d(∇x˙ F ) − ∇x F = 0 , dt o, de forma escalar equivalente d ∂F ∂F = 0 , j = 1, . . . , n. − dt ∂ x˙ j ∂xj

(8.5)

(ii) S´ olo el punto inicial es fijo: Supongamos que el conjunto E donde buscamos soluciones sea el espacio de funciones x(t) de clase C 1 [ta , tb ] tal que x(ta ) = xa con xa ∈ Rn fijo. Se puede demostrar que, en este caso, se deben de satisfacer también las ecuaciones de Euler–Lagrange (8.5) pero se han de a¯ nadir unas condiciones de transversalidad dadas por ˙ b ); tb ) = 0, es decir, ∇x˙ F (x(tb ), x(t ∂F = 0 , j = 1, . . . , n. (8.6) ∂ x˙ j t=tb

80

C´ alculo de Variaciones Para demostrar este hecho, notar que ahora s´ olo se tiene ξ(ta ) = 0, de modo que (8.4) adopta la forma tb d(∇x˙ F ) ∇x˙ F t=t ξ(tb ) + ∇x F − ξ(t) dt = 0 . (8.7) b dt ta Esta ecuación debe ser válida para ξ verificando ξ(ta ) = 0. En particular debe ser cierta también para aquellas ξ que además verifiquen ξ(tb ) = 0. Para estas segundas ξ, aplicando el Lema Fundamental del C´ alculo de Variaciones, la ecuación (8.7) reduce a las ecuaciones de Euler–Lagrange (8.5). De este modo, a partir de ahora, (8.7) reduce s´ olo a su primer término, es decir ∇x˙ F t=t ξ(tb ) = 0 . b

Finalmente, las condiciones de transversalidad (8.6) provienen del hecho de que la anterior ecuaci´ on debe verificarse para todo ξ(tb ) ∈ Rn , de modo que, ∇x˙ F t=t = 0. b

8.3.3.

Integrales primeras en casos simples

Las ecuaciones de Euler–Lagrange (8.5) son un conjunto de n ecuaciones diferenciales de segundo orden para las funciones xi (t), con i = 1, . . . , n. Por lo tanto, en general no son resolubles de manera exacta y se requieren métodos numéricos. Sin embargo, existen casos especiales que ocurren a menudo en las aplicaciones y son sencillos de tratar. (a) Si F es independiente de xj para cierta j, entonces de (8.5) se tiene ∂F =c, ∂ x˙ j

(8.8)

con c constante. (b) Si F no depende expl´ıcitamente de t entonces se verifica ˙ =c, x(t) ˙ ∇x˙ F − F (x, x)

(8.9)

con c constante. Para demostrarlo, multipliquemos por x(t) ˙ las ecuaciones de Euler–Lagrange (8.5) en forma vectorial x(t) ˙

d(∇x˙ F ) − x(t) ˙ ∇x F = 0 , dt

que, reagrupando da d(x(t) ˙ ∇x˙ F ) −x ¨(t) ∇x˙ F − x(t) ˙ ∇x F = 0 . dt Como

dF ∂F =x ¨(t) ∇x˙ F + x(t) , ˙ ∇x F + dt ∂t

8.4 Extremos con restricciones

81

se tiene la siguiente forma alternativa de las ecuaciones de Euler–Lagrange ∂F d [x(t) ˙ ∇x˙ F − F ] + =0. dt ∂t Por supuesto de aqu´ı se deduce que, si F no depende expl´ıcitamente de t, entonces se verifica (8.9).

8.4.

Extremos con restricciones

Supongamos que se desea hallar un extremo del funcional

tb

F (x(t), x(t); ˙ t) dt .

J(x) =

(8.10)

ta

pero ahora, la funci´ on x(t) está sujeta a restricciones2 . Restricciones integrales: x(t) ha de verificar

tb

G(x(t), x(t); ˙ t) dt = k ,

K(x) =

(8.11)

ta

siendo k una cierta constante. Este problema se resuelve mediante el método de los multiplicadores de Lagrange. En otras palabras, los extremos de (8.10) sujetos a (8.11) se buscan como los extremos de

tb

F (x(t), x(t); ˙ t) + λG(x(t), x(t); ˙ t) dt ,

J(x) + λK(x) = ta

para cierto multiplicador constante λ ∈ R. Dicho de otro modo, se tienen las ecuaciones de Euler–Lagrange d [∇x˙ (F + λG)] − ∇x (F + λG) = 0 . dt

(8.12)

Restricciones no integrales: x(t) ha de verificar G(x(t), x(t); ˙ t) = 0 .

(8.13)

Este problema se resuelve mediante el método de los multiplicadores de Lagrange λ(t) ahora no constantes. En otras palabras, los extremos de (8.10) sujetos a (8.13) se buscan como los extremos de Jλ (x) =

tb

F (x(t), x(t); ˙ t) + λ(t)G(x(t), x(t); ˙ t) dt . ta

2 La justificaci´ on de los métodos empleados en esta secci´ on se puede ver en el Apéndice al final de este cap´ıtulo.

82


Dicho de otro modo, se vuelven a tener las ecuaciones de Euler–Lagrange (8.12) pero ahora se ha de tener en cuenta que λ = λ(t), es decir, d [∇x˙ (F + λ(t)G)] − ∇x (F + λ(t)G) = 0 . dt

(8.14)

Este tipo de problemas de optimizaci´ on con restricciones no integrales es fundamental en el cap´ıtulo siguiente.

8.5.

Ap´ endice: Multiplicadores de Lagrange

Es bien conocido el siguiente resultado de an´ alisis. Teorema 68 (Lagrange). Sean f, g : Rn → R funciones de clase C 1 tales on que f (x) tiene un extremo en el punto x0 sobre la hipersuperficie de ecuaci´ g(x) = c con c ∈ R. Entonces, si ∇g(x0 ) = 0, existe un λ ∈ R (llamado multiplicador de Lagrange) tal que ∇f (x0 ) = λ∇g(x0 ). Demostraci´ on. Parametrizamos la hipersuperficie de Rn de ecuación g(x) = c de la forma r(t) = (x1 (t), . . . , xn (t)) con xi (t) ∈ C 1 para i = 1, . . . , n. Restringimos ahora la funci´ on f a los puntos de la anterior hipersuperficie, es decir, definimos la funci´ on h(t) = f (r(t)). Por ser x0 ∈ Rn un extremo de f sabemos que existe un t0 ∈ R que es extremo de h(t), siendo h(t0 ) = f (r(t0 )). Esto ˙ 0 ) = 0, de modo que, aplicando la regla de la cadena se tiene implica que h(t ˙ 0 ) = ∇f (x0 ) r(t h(t ˙ 0) = 0 . Esto implica que los vectores ∇f (x0 ) y r(t ˙ 0 ) son ortogonales. Como además, ˙ 0 ) también por las propiedades del vector gradiente, los vectores ∇g(x0 ) y r(t son ortogonales, se deduce que los vectores ∇f (x0 ) y ∇g(x0 ) son paralelos, es decir, existe un λ ∈ R tal que ∇f (x0 ) = λ∇g(x0 ). Nota 69. El Teorema 68 implica que los extremos de f (x) sujetos a la restricci´ on g(x) = c son los extremos de la funci´ on F (x, λ) = f (x) + λg(x) sin restricciones. Nota 70. El Teorema 68 se puede generalizar de modo que los extremos de f (x) sujetos a las restricciones gi (x) = ci con ci consantes reales mpara i = 1, . . . , m son los extremos de la funci´ on F (x, λ1 , . . . , λm ) = f (x) + i=1 λi gi (x) sin restricciones.

8.5.1.

Problemas isoperim´ etricos

Consideremos el problema de hallar un extremo del funcional

tb

F (x(t), x(t); ˙ t) dt .

J(x) = ta

(8.15)

8.5 Ap´ endice: Multiplicadores de Lagrange

83

pero ahora, la funci´ on x(t) ∈ C 1 [ta , tb ] además de tener los extremos fijos, es decir, x(ta ) = xa y x(tb ) = xb con xa y xb fijos, está sujeta a la restricción tb G(x(t), x(t); ˙ t) dt = k , (8.16) K(x) = ta

siendo k una cierta constante. Supondremos que F y G son de clase C 2 . on con x∗ (ta ) = xa y x∗ (tb ) = xb que exSea x∗ (t) ∈ C 1 [ta , tb ] una funci´ tremiza a J con la restricción K(x∗ ) = k. Consideremos la siguiente familia biparamétrica de funciones x(t) = x∗ (t) + 1 ξ1 (t) + 2 ξ2 (t) , ametros y ξi ∈ C 1 [ta , tb ]. Como los puntos inicial y final son fijos, con i ∈ R par´ se debe tener ξi (ta ) = ξi (tb ) = 0 para i = 1, 2. Para resolver el problema, utilizamos el método de los multiplicadores de Lagrange, es decir, definimos la funci´ on tb F (x(t), x(t); ˙ t) + λG(x(t), x(t); ˙ t) dt , F( 1 , 2 ; λ) = J(x) + λK(x) = ta

donde x(t) es la familia biparamétrica introducida. Una condici´ on necesaria para que x∗ (t) extremice a J(x) con la restricción K(x) = k es ∂F ∂F = 0 , i = 1, 2, =0. ∂ i 1 =2 =0 ∂λ 1 =2 =0 Desarrollando la primera anterior condici´ on se tiene tb ∂G ∂F ∂F = +λ dt . ∂ i ∂ ∂ i i ta Como

∂x(t) ∂ x(t) ˙ = ξi (t) , = ξ˙i (t) , ∂ i ∂ i se tiene que, aplicando la regla de la cadena, ∂F ∂ i ∂F ∂ i

=

∇x F ξi (t) + ∇x˙ F ξ˙i (t) ,

=

∇x F ξi (t) + ∇x˙ F ξ˙i (t) .

Sustituyendo estas expresiones en el integrando de ∂F/∂ i y realizando la integraci´ on por partes tb !tb tb d(∇ F ) x˙ ξi (t) dt , ∇x˙ F ξ˙i (t) dt = ∇x˙ F ξi (t) − dt ta ta ta tb !tb tb d(∇ G) x˙ ξi (t) dt , ∇x˙ G ξ˙i (t) dt = ∇x˙ G ξi (t) − dt t a ta ta

84


se tiene que ∂F ∂ i

=

∇x˙ F ξi (t) λ

!tb ta

tb

+

∇x˙ G ξi (t)

!tb ta

ta

∇x F − tb

+ ta

d(∇x˙ F ) dt

ξi (t) dt +

d(∇x˙ G) ∇x G − dt

*

ξi (t) dt

.

Como ξi (ta ) = ξi (tb ) = 0, se llega a tb d(∇x˙ (F + λG) ∂F = ∇x (F + λG) − ξi (t) dt = 0 , ∂ i 1 =2 =0 dt ta an evaluados sobre la curva para toda ξi (t) admisible y donde los gradientes est´ alculo de Variaciones, se debe x∗ (t). Entonces, por el Lema Fundamental del C´ verificar la ecuación d(∇x˙ (F + λG)) − ∇x (F + λG) = 0 , dt como se quer´ıa demostrar. Recordemos el siguiente resultado. Teorema 71 (Green). Sea D ⊂ R2 un dominio simplemente conexo (es decir, sin huecos) y sean P, Q ∈ C 1 (D). Sea γ ⊂ D una curva cerrada simple y orientada en sentido antihorario que encierra un recinto R ⊂ D. Entonces + ∂Q ∂P − P (x, y)dx + Q(x, y)dy = dxdy . ∂x ∂y γ R Sea A el área del recinto R. Entonces 1 1 ∂Q ∂P − dxdy = 2dxdy = A = dxdy 2 2 ∂x ∂y R R R + 1 xdy − ydx , = 2 γ donde en el u ´ltimo paso se ha utilizado el Teorema de Green. Nota 72. El nombre “Problemas isoperimétricos”proviene del siguiente problema: hallar la curva cerrada plana de per´ımetro p fijado que encierre un ´ area m´ axima. Observar que, si suponemos que la curva es C 1 , es decir, sus puntos admiten una parametrizaci´ on del tipo (x(t), y(t)) ∈ R2 para t ∈ [t0 , t1 ] con x, y ∈ 1 on C [t0 , t1 ], entonces el problema planteado consiste en hallar la parametrizaci´ que maximiza el funcional a ´rea 1 t1 A(x, y) = (xy˙ − xy) ˙ dt , 2 t0


85

sujeto a la restricci´ on

t1

P (x, y) =

x˙ 2 + y˙ 2 dt = p .

t0

Se puede demostrar con las técnicas utilizadas en este cap´ıtulo que la soluci´ on es una circunferencia de per´ımetro p.

8.6.

Problemas resueltos

1. Considerar una colecci´ on de n part´ıculas de masa m moviéndose en el espacio. La localizaci´ on de las part´ıculas para tiempo t viene dada por el vector x(t) = (x1 (t), . . . , x3n (t)) ∈ R3n . Supongamos que las fuerzas que act´ uan sobre las part´ıculas F (x; t) = (F1 (x; t), . . . , F3n (x; t)) ∈ R3n provienen de un potencial V (x; t), es decir, F = −∇x V . Demostrar que, cuando el sistema evoluciona desde un estado inicial x(t0 ) = x0 hasta un estado final x(tf ) = xf fijos, el camino que realiza x(t) es un extremal del funcional tf

J(x) =

L(x(t), x(t); ˙ t) dt , t0

conocido como acción donde la funci´ on L llamada lagrangiana es la resta de energ´ıa cinética menos potencial, es decir, L(x(t), x(t); ˙ t) =

1 mx ˙ 2 − V (x; t) . 2

´ n: Demostrar que minimizar J equivale a la segunda ley de Indicacio Newton. Soluci´ on. Sea x(t) una funci´ on de clase C 2 verificando x(t0 ) = x0 y x(tf ) = xf y que haga J extremal. Entonces, x(t) debe verificar las ecuaciones de Euler–Lagrange (8.5), es decir, ∂L d ∂L = 0 , j = 1, . . . , 3n. − dt ∂ x˙ j ∂xj Teniendo en cuenta que L(x(t), x(t); ˙ t) = se tiene

3n 1 2 m x˙ − V (x; t) , 2 i=1 i

∂L ∂L ∂V = mx˙ j , =− . ∂ x˙ j ∂xj ∂xj

Entonces, las ecuaciones de Euler–Lagrange adoptan la forma m¨ xj = −

∂V , ∂xj

86

C´ alculo de Variaciones que no son más que las ecuaciones de Newton de la mecánica clásica teniendo en cuenta que Fj = −∂V /∂xj . 2. Teniendo en cuenta el ejemplo introductorio de este cap´ıtulo, hallar la curva de descenso m´ as r´ apida (braquistocrona) que pasa por los puntos de coordenadas (0, 0), (x1 , y1 ). Soluci´ on. Se ha visto en la introducci´ on de este cap´ıtulo que se trata de hallar la funci´ on y(x) definida para 0 ≤ x ≤ x1 , con y(0) = 0 e y(x1 ) = y1 que minimice el funcional x1 & 1 + y (x)2 J(y) = dx . x 0 El funcional J es de la forma J(y) =

x1

F (y(x), y (x); x) dx ,

0

&

siendo

1 + y (x)2 . x Notar que F es independiente de y(x) de modo que se verificar´ a (8.9), es decir, la ecuación de Euler-Lagrange reduce a

F (y(x), y (x); x) =

∂F =c, ∂y siendo c ∈ R una constante. En nuestro caso se tiene y (x)

x(1 + y (x)2 )

=c,

para todo x ∈ (0, x1 ). Despejando se llega a y (x)2 =

c2 x . 1 − c2 x

Una forma simple de resolver esta ecuación diferencial consiste en usar la parametrización x(t) =

1 (1 − cos t) , η(t) := y(x(t)) . 2c2

Entonces, η(t) debe satisfacer η˙ 2 (t) = [y (x)x] ˙2=

1 − cos t sin2 t (1 − cos t)2 c2 x 2 x ˙ = = , 1 − c2 x 1 + cos t 4c4 4c4

de donde η(t) ˙ =±

1 − cos t . 2c4


87

La condición de contorno η(0) = y(x(0)) = 0 y el hecho de que y ≥ 0 implica t − sin t . η(t) = 2c2 Nota: Para que la cicloide (x(t), η(t)) sea la braquistocrona se ha de verificar también que exista un t0 tal que x0 = x(t0 ) =

t0 − sin t0 1 . (1 − cos t0 ) , y0 = η(t0 ) = 2c2 2c2

3. Determinar la curva plana y(x) de clase C 1 con sus puntos frontera fijos (x0 , y(x0 )) y (x1 , y(x1 )) que, al girar alrededor del eje de abscisas forme una superficie de ´ area m´ınima. Soluci´ on. Sea (x, y(x)) la gr´ afica de la curva soluci´ on.

Se sabe que el área de la superficie de revoluci´ on que genera dicha curva viene dada por el funcional x1 y 1 + y 2 dx . J(y) = 2π x0

La ecuación de Euler–Lagrange para este funcional es ∂F d ∂F =0, − ∂y dx ∂y

88

C´ alculo de Variaciones siendo F (y, y ; x) = y 1 + y 2 . Como F no depende de x expl´ıcitamente, entonces se verifica (8.9), es decir, F − y

∂F =c, ∂y

con c constante. En nuestro caso se tiene y

yy 2 1 + y 2 − =c. 1 + y 2

Simplificando se llega a

y

1 + y 2

=c.

(8.17)

La forma más sencilla de integrar esta ecuación consiste en realizar el cambio y = sinh ξ, de modo que, teniendo en cuenta la relaci´ on fundamental de la trigonometr´ıa hiperb´ olica cosh2 ξ − sinh2 ξ = 1, (8.17) se reescribe como y = c cosh ξ. Se tiene pues dx =

c sinh ξ dξ dy = c dξ , = y sinh ξ

que integrando da lugar a x(ξ) = cξ + c1 , con c1 constante. Hemos llegado a que la ecuación paramétrica de la curva solución es x(ξ) = cξ + c1 , y(ξ) = c cosh ξ . Por supuesto, el par´ ametro ξ se puede eliminar de las ecuaciones anteriores y se llega a la ecuación expl´ıcita de la curva y(x) = c cosh

x − c1 , c

que es una familia biparamétrica de curvas conocidas como catenarias. Nota: Por supuesto, las constantes c y c1 se calculan imponiendo que la curva pase por los puntos prefijados (x0 , y(x0 )) y (x1 , y(x1 )). Se puede ver que, seg´ un la disposici´ on en el plano de dichos puntos, existe una, dos o ninguna soluci´ on. 4. Considerar todas las curvas planas que pasan por los puntos (0, 0) y (2, 0) y que tengan longitud π. Hallar la ecuaci´ on de la curva que encierra mayor ´ area con el eje de abscisas. Soluci´ on. Sea (x, y(x)) la gr´ afica de la curva buscada.


89

y 1 0.8 0.6 0.4 0.2

0.5

1

1.5

2

x

El ´ area encerrada viene dada por el funcional 2 A(y) = F (y, y ; x) dx , 0

de donde F (y, y ; x) = y. La longitud de la curva es 2 G(y, y ; x) dx , L(y) = 0

de donde G(y, y ; x) = 1 + y 2 . El problema consiste en hallar la funci´ on y(x) que minimiza A sujeto a la restricción L(y) = π. Utilizando las ecuaciones de Euler–Lagrange (8.12) se tiene

∂ d ∂ (F + λG) = 0 , (F + λG) − dx ∂y ∂y es decir,

d λy −1=0 . dx 1 + y (x)2

Integrando respecto de x se tiene λy

1 + y 2

=c+x ,

(8.18)

con c constante. Para simplificar, tomemos y (x) = cot ξ ,

(8.19)

donde ξ ∈ [ξ0 , ξ1 ] cuando x ∈ [0, 2]. Sustituyendo en (8.18) se obtiene λ cot ξ = c + x =⇒ λ cos ξ = c + x , 1 + cot2 ξ

90

C´ alculo de Variaciones de donde, teniendo en cuenta la condici´ on x(ξ0 ) = 0, se llega a x(ξ) = λ[cos ξ − cos ξ0 ] . Ahora, de (8.19) y aplicando la regla de la cadena se llega a dy dx dy = = cot ξ(−λ sin ξ) = −λ cos ξ , dξ dx dξ que integrando da lugar a y(ξ) = λ[sin ξ0 − sin ξ] ,

(8.20)

donde se ha tenido en cuenta la condici´ on y(ξ0 ) = 0. La restricción L(y) = π implica π

2

=

1+

(y )2

dx =

0

2

0 ξ1

= λ ξ0

1 dx = sin ξ

ξ1

ξ0

1 dx dξ sin ξ dξ

cos ξ dξ = −λ[ξ1 − ξ0 ] . sin ξ

Imponiendo las condiciones x(ξ1 ) = 2 e y(ξ1 ) = 0 en las expresiones de x(ξ) e y(ξ) obtenidas anteriormente se tiene 2 = −λ[cos ξ0 − cos ξ1 ] , sin ξ0 = sin ξ1 . En definitiva se tiene el sistema de 3 ecuaciones con 3 incógnitas siguiente π

−λ[ξ1 − ξ0 ] , −λ[cos ξ0 − cos ξ1 ] ,

=

2 = sin ξ0

= sin ξ1 ,

Estas ecuaciones se satisfacen para ξ0 = 0, ξ1 = π y λ = −1. La curva solución es x(ξ) = 1 − cos ξ , y(ξ) = sin ξ , ξ ∈ [0, π] , que no es más que la semicircunferencia superior centrada en (1, 0) y de radio la unidad puesto que los puntos (x, y) de la curva verifican la ecuaci´ on (x − 1)2 + y 2 = 1. Observemos que la anterior curva resulta en un m´ aximo del área A y no en un m´ınimo por consideraciones f´ısicas. 5. Hallar el control u(t) ∈ R de modo que se minimice el funcional

tf

J(u, x) = t0

[u2 (t) + βx2 (t)] dt ,

(8.21)


91

con la restricci´ on x˙ = αx(t) + u(t) ,

(8.22)

a resiendo x(t) ∈ R verificando x(t0 ) = x0 con x0 dado pero x(tf ) no est´ stringido a ning´ un valor. Aqu´ı, α, β ∈ R son par´ ametros constantes. Soluci´ on. Utilizaremos el método de los multiplicadores de Lagrange. Utilizamos pues el funcional aumentado tf F (u, x, u, ˙ x; ˙ t) + λ(t)G(u, x, u, ˙ x; ˙ t) dt , Jλ (u, x) = t0

siendo F (u, x, u, ˙ x; ˙ t) = u2 (t) + βx2 (t) y G(u, x, u, ˙ x; ˙ t) = x˙ − αx(t) − u(t). Las ecuaciones de Euler–Lagrange (8.14) son

∂ d ∂ (F + λ(t)G) − (F + λ(t)G) = 0 , dt ∂ x˙ ∂x

∂ d ∂ (F + λ(t)G) − (F + λ(t)G) = 0 , dt ∂ u˙ ∂u que, en nuestro caso dan lugar a λ˙ − 2βx + λα 2u − λ

= 0, = 0.

(8.23) (8.24)

Por otra parte, la condici´ on de transversalidad (8.6) es ∂(F + λ(t)G) =0, ∂ x˙ t=tf es decir, λ(tf ) = 0. Por supuesto, debido a (8.24), también se tiene u(tf ) = 0. Eliminando λ del sistema (8.23) y (8.24) se tiene u˙ = βx − αu que, junto con (8.22) da lugar al sistema diferencial lineal x˙ α 1 x = , u˙ β −α u con las condiciones iniciales x(t0 ) = x0 , u(t0 ), que es fácilmente resoluble. Realizando el proceso pertinente se obtiene u(t) = u(t0 ) cosh(wt) +

βx0 − αu(t0 ) sinh(wt) , w

siendo w2 := α2 + β 2 . Usando la condici´ on de transversalidad u(tf ) = 0 eliminamos u(t0 ) y finalmente se tiene u(t) =

βx0 sinh[w(tf − t)] . α sinh(wtf ) − w cosh(wtf )

92

C´ alculo de Variaciones Nota: Existe otra forma de resolver el problema y consiste en sustituir el valor de u en (8.21) del despeje de (8.21). De este modo el problema se convierte en un problema de extremos del funcional J sin restricciones donde el punto inicial es fijo pero el final no. Se aplican pues las ecuaciones de Euler–Lagrange (8.5) con las condiciones de transversalidad (8.6). 6. Hallar el control u(t) ∈ R y la trayectoria x(t) ∈ R que minimice el funcional ∞ [x2 (t) + 4u2 (t)] dt , (8.25) J(u, x) = 0

sujeto a las restricciones x ¨ = −x(t) ˙ + u(t) ,

(8.26)

˙ = (0, 0). con x(0) = 0, x(0) ˙ = 1 y l´ımt→∞ (x(t), x(t)) Soluci´ on. Nos encontramos ante un problema de extremos de un funcional con una restricción (8.26) que es una ecuaci´ on diferencial de segundo orden. Por lo tanto, antes de utilizar el método de los multiplicadores de Lagrange, convertiremos la restricción (8.26) en dos restricciones (ecuaciones diferenciales) de primer orden. Como es habitual en este procedimiento, ˙ de modo que (8.26) se reescribe de definimos x1 (t) = x(t) y x2 (t) = x(t), la forma x˙ 1 = x2 (t) , x˙ 2 = −x2 (t) + u(t) . Utilizamos ahora el funcional aumentado ∞ F + λ1 (t)G1 + λ2 (t)G2 dt , Jλ (u, x1 , x2 ) = 0

siendo F (u, x1 , x2 , u, ˙ x˙ 1 , x˙ 2 ; t) = x21 (t) + 4u2 (t), G1 (u, x1 , x2 , u, ˙ x˙ 1 , x˙ 2 ; t) = ˙ x˙ 1 , x˙ 2 ; t) = x˙ 2 + x2 (t) − u(t). x˙ 1 − x2 (t) y G2 (u, x1 , x2 , u, Las ecuaciones de Euler–Lagrange (8.14) son, para j = 1, 2,

d ∂ ∂ (F + λ1 (t)G1 + λ2 (t)G2 ) − (F + λ1 (t)G1 + λ2 (t)G2 ) = dt ∂ x˙ j ∂xj

∂ d ∂ (F + λ1 (t)G1 + λ2 (t)G2 ) − (F + λ1 (t)G1 + λ2 (t)G2 ) = dt ∂ u˙ ∂u

0, 0,

que, en nuestro caso dan lugar a λ˙ 1 (t) − 2x1 (t) = 0 , λ˙ 2 (t) + λ1 (t) − λ2 (t) = 0 , 8u(t) − λ2 (t) = 0 .

(8.27) (8.28) (8.29)

Eliminando de este sistema los multiplicadores de Lagrange λ1 (t) y λ2 (t) y recordando que x1 (t) = x(t) se llega a 4¨ u = 4u˙ − x .


93

Esta ecuación junto con (8.26) resulta un sistema lineal resoluble. Realizando el procedimiento habitual de resoluci´ on se obtiene √ √ √ √ x(t) = A exp(−t/ 2) + Bt exp(−t/ 2) + C exp(t/ 2) + Dt exp(t/ 2) . ˙ = Teniendo en cuenta las condiciones x(0) = 0, x(0) ˙ = 1 y l´ımt→∞ (x(t), x(t)) (0, 0) se llega a A = C = D = 0, B = 1. En definitiva se obtiene √ √ √ 1 x(t) = t exp(−t/ 2) , u(t) = (1 − 2) 1 + t exp(−t/ 2) . 2

Cap´ıtulo 9

´ Control Optimo 9.1.

Introducci´ on

En este cap´ıtulo se describirá cómo controlar un sistema de forma que se comporte de la “mejor forma posible”. Por supuesto, la estrategia del control utilizado depender´ a de cómo se entienda “mejor”. El sistema será modelizado mediante una ecuación diferencial de primer orden, en general no lineal, de la forma x˙ = f (x, u, t) , x(t0 ) = x0 ,

(9.1)

siendo x ∈ Rn las variables de estado y u ∈ Rm las de control. Supondremos siempre suficiente regularidad a la funci´ on f de modo que se verifique la existencia y la unicidad de la soluci´ on x(t) de (9.1). Definici´ on 73. Un ´ındice de actuación J ser´ a un funcional que a cada pareja (x(t), u(t)) compatible con (9.1) le asigne un escalar real que proporciona una medida a partir de la cual es posible juzgar si el sistema (9.1) se ha comportado de la “mejor forma posible”. Veamos algunos ejemplos. Problemas de tiempo m´ınimo: Se ha de elegir el control u(t) de modo que el sistema evolucione desde el estado inicial x(t0 ) = x0 hasta un estado final x(t1 ) = xf en el menor tiempo posible, es decir, minimizando el ´ındice de actuaci´ on J = t1 − t0 que escribiremos de la forma t1 J= dt . (9.2) t0

Problemas de control terminal: Se ha de elegir el control u(t) de modo que el estado final x(t1 ) = xf esté lo más cerca posible de un cierto estado de 95

´ Control Optimo

96

referencia deseado r(t1 ). En este caso, un planteamiento adecuado puede ser minimizar el ´ındice de actuación dado por la forma cuadr´ atica J = ∆T (t1 )M ∆(t1 ) ,

(9.3)

siendo ∆(t) = x(t) − r(t) y M ∈ Mn (R) simétrica y definida positiva. Observar que un caso particular de especial relevancia aparece cuando se toma M = In . En este caso se tiene J = x(t1 ) − r(t1 )2 . Problemas de esfuerzo m´ınimo: Elegir el control u(t) de modo que el estado final x(t1 ) = xf se alcance con el menor “esfuerzo”posible. Existen diversos ´ındices de actuación dependiendo de como se defina el “esfuerzo”. Por ejemplo se pueden minimizar los ´ındices de actuación siguientes: " t1 m βi |ui (t)| dt , (9.4) J(u) = t0

o bien

i=1

t1

J(u) =

uT (t)M u(t) dt ,

(9.5)

t0

siendo M ∈ Mn (R) simétrica y definida positiva. Servomecanismos y reguladores: Si el control u(t) se elige de modo que el objetivo es que el sistema siga lo más cerca posible a un cierto estado de referencia r(t) durante un cierto intervalo de tiempo t0 ≤ t ≤ t1 entonces se habla de un servomecanismo. Un posible ´ındice de actuación puede ser t1 ∆T (t)M ∆(t) dt , (9.6) J(x) = t0

siendo ∆(t) = x(t)−r(t) y M ∈ Mn (R) simétrica y definida positiva. Para el caso particular de que r(t) sea constante, los servomecanismos reciben el nombre especial de reguladores. Si el control u(t) no está acotado, el problema de minimizar (9.6) puede dar lugar a soluciones con alguna componente ui de u infinita, lo cual no es aceptable en problemas reales. Por este motivo y con el objetivo de minimizar también el esfuerzo total realizado en el control efectuado, se puede tomar el ´ındice de actuación siguiente t1 [∆T (t)M ∆(t) + uT (t)P u(t)] dt , (9.7) J(u, x) = t0

siendo ∆(t) = x(t)−r(t) y M, P ∈ Mn (R) simétricas y definidas positivas. Los ´ındices de actuación (9.5), (9.6) y (9.7) se conocen como ´ındices de actuaci´ on cuadr´ aticos. De forma más general, escribiremos un ´ındices de actuación de la forma t1

J(u, x) = φ(x(t1 ), t1 ) +

F (x(t), u(t), t) dt , t0

siendo las funciones escalares φ y F de clase C 1 .

(9.8)

9.2 Control ´ optimo: m´ etodo hamiltoniano

9.2.

97

Control o ´ptimo: m´ etodo hamiltoniano

Consideraremos el problema de minimizar el funcional ´ındice de actuación J(u) definido por (9.8) sujeto a la ecuaci´ on diferencial (9.1), es decir, x˙ = f (x, u, t) , x(t0 ) = x0 ,

(9.9)

con u(t)∈ Rm y x(t) ∈ Rn . No supondremos en esta sección que el control u(t) está sujeto a ninguna otra restricci´ on. Definici´ on 74. El control u∗ (t) es un extremal y J(u) tiene un m´ınimo relativo ∗ en u si existe un > 0 tal que J(u) − J(u∗ ) ≥ 0 para todo u verificando u − u∗ < . Una aplicaci´ on del cálculo de variaciones a la teor´ıa de control ´ optimo es la siguiente. Teorema 75. Condiciones necesarias para que u∗ (t) ∈ Rm sea un extremal del funcional t1 J(u, x) = F (x(t), u(t), t) dt , t0

sujeto a la ecuaci´ on diferencial x˙ = f (x, u, t) , x(t0 ) = x0 , siendo x = (x1 , . . . , xn )T y u = (u1 , . . . , um ) ∈ C 1 ([t0 , t1 ]), son las siguientes: λ˙ = −∇x H ,

(9.10)

siendo la funci´ on Hamiltoniana H asociada al funcional J, H(x, u, λ, t) = λ(t) f (x, u, t) − F (x, u, t) ,

(9.11)

donde λ(t) = (λ1 (t), . . . , λn (t))T son multiplicadores de Lagrange y H verifica (∇u H)u=u∗ = 0 ,

(9.12)

para todo t0 ≤ t ≤ t1 . Demostraci´ on. Consideramos el funcional t1 J(u, x) = F (x(t), u(t), t) dt , t0

con las n restricciones x˙ = f (x, u, t) . Introducimos n multiplicadores de Lagrange λi (t), con i = 1, . . . , n. Las variables xi y λi se llaman conjugadas. Definimos el vector λ(t) = (λ1 (t), . . . , λn (t))T y el funcional aumentado t1 F(x(t), u(t), λ(t), x(t); ˙ t) dt , Jλ (u, x) = t0

´ Control Optimo

98

siendo F(x(t), u(t), λ(t), x(t); ˙ t) = F (x(t), u(t), t)+λ(t) [x−f ˙ (x, u, t)], donde el punto indica el producto escalar ordinario en Rn . Hemos de hallar un extremo del funcional Jλ (u, x), pero antes de hallar las ecuaciones de Euler–Lagrange, definiremos el Hamiltoniano H(x, u, λ, t) = λ(t) f (x, u, t) − F (x, u, t) , de modo que el funcional Jλ (u, x) se puede reescribir de la forma t1 λ(t) x˙ − H(x, u, λ, t) dt . Jλ (u, x) = t0

Las ecuaciones de Euler–Lagrange (8.5) para este funcional son d(∇x˙ F) − ∇x F dt d(∇u˙ F) − ∇u F dt

= 0, = 0.

Como ∇x˙ [λ x] ˙ = λ, las anteriores ecuaciones reducen a λ˙ ∇u H

= −∇x H , = 0,

y son conocidas como las ecuaciones de Hamilton–Pontryagin. Nota 76. Es interesante resaltar que la ecuaci´ on de estado x˙ = f (x, u, t) se puede reescribir de la forma x˙ = ∇λ H. Se tiene en definitiva el sistema hamiltoniano x˙ = ∇λ H , λ˙ = −∇x H , de donde proviene el hecho de que λ sea la variable conjugada de x. Nota 77. Si el problema de extremizar el funcional J(u, x) del Teorema 75 no nadirle restringe el valor del punto final x(t1 ), entonces al Teorema 75 hay que a¯ la condici´ on de transversalidad (∇x˙ F)t=t1 = 0, es decir λ(t1 ) = 0 .

(9.13)

Si el funcional a extremizar J(u, x) viene dado por (9.8), entonces la anterior condici´ on se modifica de la forma λ(t1 ) = (∇x φ)t=t1

(9.14)

Nota 78. La ecuaci´ on de estado x˙ = f (x, u, t) junto con la ecuaci´ on adjunta (9.10) escrita en forma compacta como λ˙ = −∇x H(x, u, λ, t), definen un sistema de 2n ecuaciones diferenciales (no lineales en general) con condiciones de on expresable medicontorno x(t0 ) y λ(t1 ) que no admite habitualmente soluci´ ante funciones elementales. En estos casos, los métodos numéricos son de gran utilidad.

9.3 El regulador lineal

99

Nota 79. En muchos libros de texto, el Teorema 75 se utiliza con el hamiltoniano (totalmente equivalente) H(x, u, λ, t) = F (x, u, t) + λ(t) f (x, u, t) .

9.3.

El regulador lineal

En la Nota 78 queda claro que, en general, es imposible expresar el control ´ptimo u∗ (t) mediante funciones elementales. Sin embargo, esto si es posible o cuando el sistema es lineal. Sea x(t) ∈ Rn , u(t) ∈ Rm , A ∈ Mn (R) y B ∈ Mm×n (R) y consideremos el sistema lineal x˙ = Ax + Bu , con un ´ındice de actuación dado por el funcional cuadr´ atico 1 1 t1 T J(u, x) = xT (t1 )M x(t1 ) + [x Qx + uT Ru] dt , 2 2 t0

(9.15)

(9.16)

siendo M, Q ∈ Mn (R), R ∈ Mm (R) matrices reales constantes simétricas y definidas positivas. Los factores 1/2 han sido introducidos por conveniencia en la notaci´ on futura. Proposici´ on 80. El ´ındice de actuaci´ on cuadr´ atico (9.16) sujeto al sistema lineal (9.15) se minimiza mediante un control o ´ptimo por realimentaci´ on u∗ (t) = −1 T ∗ on −R B P (t)x (t), donde P (t) es una matriz simétrica que satisface la ecuaci´ diferencial matricial de Riccati P˙ = P BR−1 B T P − AT P − P A − Q , con la condici´ on inicial P (t1 ) = M . Demostraci´ on. El Hamiltoniano asociado al problema de optimizaci´ on planteado es, seg´ un la Nota 79, H=

1 T 1 x Qx + uT Ru + λT (Ax + Bu) . 2 2

La condición (9.12) de optimabilidad (∂H/∂ui )u=u∗ = 0, con i = 1, . . . , m se escribe como (∇u H)u=u∗ = 0 en forma vectorial mediante el operador gradiente. En nuestro caso, teniendo en cuenta la forma del Hamiltoniano H y el hecho de que R es simétrica, se obtiene Ru∗ + B T λ = 0 , de modo que, puesto que R es invertible por ser definida positiva, el control optimo viene dado por ´ (9.17) u∗ (t) = −R−1 B T λ(t) .

´ Control Optimo

100

Las ecuaciones adjuntas que han de verificar los multiplicadores de Lagrange un (9.10), de modo que, λ(t) son λ˙ = −∇x H seg´ λ˙ = −Qx∗ − AT λ .

(9.18)

Sustituyendo (9.17) en (9.15) se tiene x˙ ∗ = Ax∗ − BR−1 B T λ . Combinando esta expresi´ on con (9.18) producimos el siguiente sistema de 2n ecuaciones lineales ∗ ∗ x˙ x A −BR−1 B T . (9.19) = λ −Q −AT λ˙ Se tiene además la condición de contorno λ(t1 ) = (∇x φ)t=t1 , ver la Nota 77. on de contorno da lugar a Como en nuestro caso φ = 12 xT M x, dicha condici´ λ(t1 ) = M x∗ (t1 ) .

(9.20)

Utilizamos ahora la f´ ormula (3.7) que nos da la soluci´ on del sistema (9.19) con la condici´ on inicial colocada en tiempo t1 , es decir, ∗ ∗ x (t) x (t1 ) . = Φ(t, t1 ) λ(t) λ(t1 ) siendo Φ(t, t1 ) = (φij ) ∈ M2n (R) la matriz de transici´ on de estados del sistema (9.19). Aqu´ı se utiliza la notaci´ on por bloques φij (t) ∈ Mn (R) con i, j = 1, 2. Se tiene pues que x∗ (t) = φ11 x∗ (t1 ) + φ12 λ(t1 ) = (φ11 + φ12 M )x∗ (t1 ) . Por otra parte λ(t) = =

φ21 x∗ (t1 ) + φ22 λ(t1 ) = (φ21 + φ22 M )x∗ (t1 ) (φ21 + φ22 M )(φ11 + φ12 M )−1 x∗ (t) := P (t)x∗ (t) ,

(9.21)

donde se ha definido la matriz P (t) = (φ21 +φ22 M )(φ11 +φ12 M )−1 . Nótese que, se puede demostrar que la matriz φ11 + φ12 M es invertible para todo t ≥ t0 . En definitiva se ha demostrado que el control o´ptimo u∗ (t) viene dado por una realimentaci´ on lineal del tipo u∗ (t) = −R−1 B T P (t)x∗ (t) ,

(9.22)

donde la matriz P (t) verifica una ecuaci´ on que a continuaci´ on deducimos. Derivando (9.21) respecto del tiempo t se tiene P˙ x∗ + P x˙ ∗ − λ˙ ∗ = 0 .

9.4 Teor´ıa de Pontryagin

101

Sustituyendo aqu´ı las expresiones de x˙ ∗ y de λ˙ ∗ obtenidas de (9.19) y (9.21) respectivamente se tiene (P˙ + P A − P BR−1 B T P + Q + AT P )x∗ (t) = 0 . Como esta ecuación ha de ser cierta para todo t0 ≤ t ≤ t1 se tiene que P (t) debe verificar la llamada ecuaci´ on diferencial matricial de Riccati P˙ = P BR−1 B T P − AT P − P A − Q ,

(9.23)

con la condici´ on inicial P (t1 ) = M ,

(9.24)

que se obtienen de (9.20) y (9.21). Nota 81. Como M es simétrica, de (9.24) se deduce que P (t) es simétrica para todo t. De este modo, la ecuaci´ on diferencial matricial de Riccati (9.23) representa n(n + 1)/2 ecuaciones diferenciales cuadr´ aticas de primer orden que se han de resolver, en general, utilizando métodos numéricos. Es importante en algunas aplicaciones el caso especial en el cual t1 tiende hacia infinito y M = 0. Se demuestra que en este caso l´ımite la matriz P es independiente del tiempo como muestra el siguiente resultado, ver por ejemplo [6] para una demostraci´ on. ˆ Si ˆ la matriz de igual rango que Q tal que Q = Q ˆ T Q. Proposici´ on 82. Sea Q ˆ el sistema lineal (9.15) es completamente controlable y el par [A, Q] es completamente observable, entonces el ´ındice de actuaci´ on cuadr´ atico ∞ [xT Qx + uT Ru] dt , J(x, u) = 0

sujeto al sistema lineal (9.15) se minimiza mediante un control o ´ptimo por realimentaci´ on u∗ (t) = −R−1 B T P x∗ (t), donde P es la u ńica matriz real simétrica y definida positiva que satisface la ecuaci´ on matricial algebraica de Riccati P BR−1 B T P − AT P − P A − Q = 0 .

9.4.

Teor´ıa de Pontryagin

Una gran variedad de procesos f´ısicos pueden ser controlados de m´ ultiples formas, pero sólo una forma es la mejor bajo un cierto criterio. Se trata pues de averiguar cual es el control ´ optimo del proceso en cuestión. Quiz´ as interese realizar el objetivo del proceso en el menor tiempo posible o bien con el m´ınimo consumo de energ´ıa, etc... Supondremos que el sistema está gobernado por x˙ = f (x, u) , x(t0 ) = x0

(9.25)

´ Control Optimo

102

con x = (x1 , . . . , xn )T las variables de estado o de fase, f = (f1 , . . . , fn )T ∈ as suponemos que f Rn y u = (u1 , . . . , um )T las variables de control. Adem´ suficientemente regular de modo que, dado un control u(t) definido para t0 ≤ ńica en el intervalo t0 ≤ t ≤ t1 . En t ≤ t1 , la solución x(t) de (9.25) existe y es u particular, supondremos que fi y ∂fi /∂xj son funciones continuas. Consideremos el funcional t1 F (x, u) dt , (9.26) J(u, x) = t0

de modo que, para cada control dado u(t) definido para t0 ≤ t ≤ t1 , se tiene J(u, x) bien definido y tomando un valor real fijado. Fijados el estado inicial x0 y final xf del sistema, supongamos que existen varios controles u(t) que transfieren el sistema desde el estado inicial x(t0 ) = x0 hasta el estado final x(t1 ) = xf en tiempo t1 − t0 . Nos preguntamos si, de entre esos controles, existe uno u∗ tal que se minimize el funcional (9.26). Dicho optimo y la solución x∗ (t) de (9.25) asociado al control u∗ (t) es llamado control ´ ∗ optima. control u una trayectoria ´ Notemos que, en el problema planteado, los valores de t0 y de t1 no están en general fijados. Lo u ńico que se pide es que para tiempo t0 el sistema se encuentre en el estado x0 y para tiempo t1 en el estado xf . Es importante notar que, en los problemas aplicados, muy a menudo se tiene an que las variables de control ui no pueden tomar valores arbitrarios si no que est´ acotadas y además pueden ser discontinuas. Entonces, s´ olo tiene sentido tomar controles u pertenecientes a un cierto dominio U ⊂ Rm . As´ı, diremos que el control u(t) es admisible si u(t) ∈ U para todo t0 ≤ t ≤ t1 y además es continuo a trozos en dicho intervalo. Dicho de otro modo, pueden existir tiempos si con t0 < s1 < s2 < · · · < sn < t1 que corresponden a discontinuidades de primera especie de u(t). Un caso muy habitual de conjunto de controles admisibles U es el hipercubo de Rm siguiente , i = 1, . . . , m} ⊂ Rm . U = {u = (u1 , . . . , um ) : |ui | ≤ umax i

(9.27)

Observemos que, al a¯ nadir la condici´ on (9.27) al problema de optimizaci´ on planteado se tiene un car´ acter radicalmente diferente del que se ten´ıa con la condición (9.12), es decir, (∇u H)u=u∗ = 0 puesto que ahora el hamiltoniano H puede ser no derivable respecto de u en todos los puntos ya que u(t) puede tener puntos donde no exista la derivada. El siguiente resultado fundamental, ver su demostraci´ on en [8] por ejemplo, aporta condiciones necesarias que debe verificar un control o´ptimo u∗ donde se tiene en cuenta la condici´ on (9.27) sobre la acotación de los controles admisibles . Para conseguirlo, definimos la funci´ on de Hamilton H asociada al sistema (9.25) y al funcional (9.26) de la forma H(λ, x, u) = λT f (x, u) − F (x, u) ,

(9.28)


103

siendo λ(t) = (λ1 (t), . . . , λn (t))T ∈ Rn un vector de multiplicadores de Lagrange verificando ∂H , i = 1, . . . , n . (9.29) λ˙ i = − ∂xi Teorema 83 (Principio del m´ aximo de Pontryagin). Sea u∗ (t) definido para t0 < t < t1 un control admisible tal que la trayectoria x∗ (t) del sistema optimo para (9.25) asociada verifica x∗ (t0 ) = x0 y x∗ (t1 ) = xf . Si u∗ (t) es ´ el funcional (9.26), entonces existe un vector de multiplicadores de Lagrange as, λ(t) = (λ1 (t), . . . , λn (t))T ∈ Rn continuo y no nulo, verificando (9.29). Adem´ aximo para cualquier t0 ≤ t ≤ t1 , el hamiltoniano H(λ(t), x(t), u) alcanza su m´ on de los controles admisibles u ∈ U, es en u = u∗ (t) mirado como una funci´ decir, (9.30) H(λ(t), x∗ (t), u∗ (t)) = sup {H(λ(t), x(t), u)} . u∈U

Nota 84. Si se trata de un problema de control en tiempo ´ optimo, es decir, F (x, u) ≡ 1 en el funcional (9.26), entonces es habitual tomar la funci´ on de Hamilton (9.31) H(λ, x, u) = λT f (x, u) , en la utilizaci´ on del Teorema 83.

9.5.

Problemas resueltos

1. Minimizar el ´ındice de actuaci´ on T J(u, x) = (x(t)2 + u2 (t)) dt , 0

sabiendo que x˙ = −ax + u(t) con x(0) = x0 . Tomar a y T constantes positivas. Soluci´ on. Utilizaremos el Teorema 75 para la resolución del problema. La funci´ on Hamiltoniana H asociada al funcional J(u, x) y al sistema x˙ = −ax + u(t) viene dado, seg´ un (9.11), por H(x, u, t) = λ(−ax + u) − (x2 + u2 ) , donde el multiplicador de Lagrange λ(t) verifica, teniendo en cuenta (9.10), la ecuación λ˙ = −∂H/∂x, es decir, λ˙ = 2x∗ + aλ ,

(9.32)

donde se ha denotado por x∗ (t) a la trayectoria ´ optima. Además, si u∗ es el control óptimo, la condici´ on (9.12), es decir, (∂H/∂u)u=u∗ = 0 impone en nuestro caso la restricción −2u∗ + λ = 0 .

(9.33)

´ Control Optimo

104

Introduciendo esta condici´ on en el sistema x˙ ∗ = −ax∗ + u∗ se obtiene 1 x˙ ∗ = −ax∗ + λ . 2

(9.34)

Teniendo en cuenta que el punto final x(T ) no está fijado se tiene la condición de transversalidad (9.14), es decir, λ(T ) = 0 . Observar que, el sistema formado por las ecuaciones (9.32) y (9.34) es un sistema lineal y puede, por lo tanto, ser resuelto con los métodos ya estudiados anteriormente. En concreto, es f´ acil demostrar que, la soluci´ on general del sistema (9.32) y (9.34) es x∗ (t)

=

λ(t)

=

1 1 + a2 t) − √ (2aC1 + C2 ) sinh( 1 + a2 t) , 2 1 + a2 1 (−2C1 + aC2 ) sinh( 1 + a2 t) , C2 cosh( 1 + a2 t) + √ 1 + a2 C1 cosh(

siendo C1 y C2 constantes arbitrarias que se determinan con las condiciones de contorno x∗ (0) = x0 y λ(T ) = 0. En concreto se obtiene C1 = x0 , C2 =

a+

√

1+

a2

2x0 √ . coth( 1 + a2 T )

un (9.33), u∗ (t) = −λ(t)/2, es En definitiva, el control o´ptimo u∗ es, seg´ decir, √ x0 sinh[ 1 + a2 (t − T )] ∗ √ √ u (t) = √ . 1 + a2 cosh[ 1 + a2 T ] + a sinh[ 1 + a2 T ] 2. Una part´ıcula realiza un movimiento rectil´ıneo sujeta a la acci´ on de una fuerza por unidad de masa u(t). Supongamos que el esfuerzo realizado cuando la part´ıcula evoluciona desde tiempo cero hasta tiempo τ viene medido por τ

J(u) =

u2 (t) dt .

0

Averiguar la expresi´ on de u(t) de modo que la part´ıcula sea transferida desde el reposo en el origen de coordenadas (x = 0) hacia el reposo en x = 1 en 1 unidad de tiempo y minimizando el esfuerzo realizado. Dar también la posici´ on ´ optima en funci´ on del tiempo. Soluci´ on. Utilizaremos el Teorema 75 para la resolución del problema. Sea x(t) la posición de la part´ıcula en funci´ on del tiempo. Seg´ un las leyes de la mecánica clásica se tiene x ¨ = u(t), o de forma equivalente, x˙ = y , y˙ = u(t) ,

(9.35)


105

donde se ha definido y(t) como la velocidad de la part´ıcula. La funci´ on Hamiltoniana H asociada al sistema (9.35) y al funcional J(u) dado por 1 J(u) = u2 (t) dt , 0

es, seg´ un (9.11),

H(x, u, t) = λ1 y + λ2 u − u2 ,

donde los multiplicador de Lagrange λi (t) verifican, teniendo en cuenta (9.10), el sistema λ˙ 1 = −∂H/∂x, λ˙ 2 = −∂H/∂y, es decir, λ˙ 1 = 0 , λ˙ 2 = −λ1 .

(9.36)

La solución del sistema (9.36) es λ1 (t) = −c1 , λ2 (t) = c1 t + c2 , siendo ci constantes reales. Si u∗ es el control óptimo, la condici´ on (9.12), es decir, (∂H/∂u)u=u∗ = 0 da lugar a (9.37) −2u∗ + λ2 = 0 . De este modo, el control óptimo es u∗ (t) =

c1 t + c2 λ2 (t) = . 2 2

Introduciendo este control en la segunda ecuaci´ on (9.35) e integrando se tiene % 1 $ c1 2 t + c2 t + c3 . y ∗ (t) = u∗ (t) dt = 2 2 con c3 constante. Puesto que la part´ıcula se encuentra para tiempo 0 y 1 en reposo se tienen las condiciones y ∗ (0) = y ∗ (1) = 0, de donde obtenemos c3 = 0 y c1 = −2c2 . En definitiva c2 , y ∗ (t) = t − t2 . 2 Integrando ahora la primera ecuaci´ on de (9.35) se tiene t3 c2 t2 ∗ ∗ x (t) = y (t) dt = − + c4 , 2 2 3 con c4 constante. Puesto que la part´ıcula se en cuentra para tiempo 0 y 1 en posiciones 0 y 1 respectivamente, es decir, x∗ (0) = 0 y x∗ (1) = 1 optima es obtenemos c4 = 0 y c2 = 12. Entonces, la trayectoria ´ t ∗ ∗ 2 1 − (x (t), y (t)) = 6t , 6t[1 − t] , 2 3 y el control ´ optimo

u∗ (t) = 6(1 − 2t) .

´ Control Optimo

106

3. Hallar un control u(t) por realimentaci´ on lineal que haga m´ınimo el funcional ∞ 1 2 2 y + u dt , 10 0 sujeto a x˙ = −x + u(t), y˙ = x. Soluci´ on. Utilizaremos la Proposici´ on 82 para solucionar el problema planteado. El sistema lineal viene definido por x˙ x =A + bu , y˙ y

siendo A=

−1 0 1 0

1 0

, B=

.

El funcional a minimizar se puede escribir de la forma ∞ x T (x, y)Q + u Ru dt , y 0 con las matrices simétricas

0 0 0 1

Q=

, R=

1 . 10

Utilizando la Proposici´ on 82, el control o ´ptimo por realimentaci´ on lineal es x , u∗ (t) = −R−1 B T P y siendo la matriz simétrica

P =

p11 p12

p12 p22

,

solución de la ecuación matricial algebraica de Riccati P BR−1 B T P − AT P − P A − Q = 0 . Es f´ acil comprobar que dicha ecuaci´ on es 2[p11 + 5p211 − p12 ] p12 + 10p11 p12 − p22 0 = p12 + 10p11 p12 − p22 −1 + 10p212 0 dando lugar a la soluci´ on ⎛ √ ! 1 −1 + 1 + 2 10 ⎜ 10 & P =⎝ √1 10

√1 10 1 10

+

⎞ .

2 5

⎟ ⎠ .

En definitiva, el control o´ptimo por realimentaci´ on lineal es

. √ √ u∗ (t) = 1 − 1 + 2 10 x − 10y .

0 0

,


107

4. Considerar un sistema regido por la ecuaci´ on x ¨ = u(t), con |u(t)| ≤ 1. Calcular la expresi´ on de u(t) que minimiza el tiempo T que tarda el sistema en llegar al estado (x(T ), x(T ˙ )) = (0, 0) y su interpretaci´ on geométrica dependiendo de las condiciones iniciales x(0) y x(0). ˙ Soluci´ on. Pasamos la ecuación de segundo orden x ¨ = u(t) a un sistema de primer orden con el cambio habitual, es decir, x˙ = y , y˙ = u(t) .

(9.38)

Como se tiene un problema de tiempo óptimo, para aplicar el Principio del máximo de Pontryagin tomaremos el Hamiltoniano (9.31) dado por H = λ 1 y + λ2 u ,

(9.39)

donde λi verifican las ecuaciones (9.29), es decir, λ˙ 1 = −∂H/∂x, λ˙ 2 = −∂H/∂y. En nuestro caso se tiene λ˙ 1 = 0 , λ˙ 2 = −λ1 ,

(9.40)

de modo que, integrando estas ecuaciones respecto del tiempo se obtiene λ1 (t) = c1 , λ2 (t) = c2 − c1 t , siendo ci constantes reales. Teniendo en cuenta que |u(t)| ≤ 1, el control ´ optimo u∗ (t) que verifica la condici´ on (9.30) de maximizar el Hamiltoniano (9.39) viene dado por la funci´ on continua a trozos −1 si λ2 < 0 , u∗ (t) = sign λ2 (t) = (9.41) 1 si λ2 > 0 . Este tipo de controles u∗ (t) se conocen en la literatura como bang–bang por motivos obvios. Como u∗ (t) = sign(c2 − c1 t) es el signo de una funci´ on lineal del tiempo, est´ a claro que como mucho el control óptimo u∗ tendr´ a una discontinuidad en el intervalo 0 ≤ t ≤ T . Estudiemos las trayectorias en el plano x − y de las fases o de estados para el sistema (9.38) dependiendo del control utilizado. Si u(t) = 1, el sistema (9.38) es x˙ = y, y˙ = 1. Eliminando el tiempo se tiene dx/dy = y que es una ecuación de variables separables. Se tiene en definitiva que las trayectorias del sistema vienen dadas por las par´ abolas (9.42) x = y 2 /2 + α , con α una constante arbitraria. Adem´ as, como y˙ = 1 > 0 está claro que y(t) es una funci´ on creciente y por lo tanto se sabe cómo se recorre la par´ abola (9.42), ver Figura 9.1.

´ Control Optimo

108

Figura 9.1: Trayectorias parab´ olicas con control u(t) = 1. Si u(t) = −1, el sistema (9.38) es x˙ = y, y˙ = −1 o, equivalentemente dx/dy = −y. Integrando esta ecuación se halla que las trayectorias del sistema vienen dadas por las par´ abolas x = −y 2 /2 + β ,

(9.43)

con β una constante arbitraria. Adem´ as, como y˙ = −1 < 0 se tiene que y(t) es una funci´ on decreciente y, en consecuencia queda claro el sentido en que se recorre la parábola (9.43), ver Figura 9.2.

Figura 9.2: Trayectorias parab´ olicas con control u(t) = −1.


109

Se tienen pues las siguientes posibilidades en la secuencia de control a aplicar dependiendo de las condiciones iniciales (x(0), y(0)) = (x0 , y0 ) ∈ R2 del sistema. Si el control óptimo u∗ (t) tiene la secuencia temporal 1 si 0 ≤ t ≤ t1 < T , ∗ u (t) = −1 si t1 < t ≤ T , la trayectoria óptima (x∗ (t), y ∗ (t)) ∈ R2 con 0 ≤ t ≤ T que sigue el sistema (9.38) viene dada por una porci´ on de la par´ abola (9.42) on de la que pase por el punto inicial (x0 , y0 ) seguido de otra porci´ parábola (9.43) que pase por el origen de coordenadas, ver Figura 9.3.

Figura 9.3: Trayectorias parab´ olicas con control bang–bang inicial u(t) = 1 y posterior u(t) = −1. La trayectoria ´ optima con condici´ on inicial (x0 , y0 ) está doblemente marcada. Si el control óptimo u∗ (t) tiene la secuencia temporal −1 si 0 ≤ t ≤ t1 < T , u∗ (t) = 1 si t1 < t ≤ T , la trayectoria óptima (x∗ (t), y ∗ (t)) ∈ R2 con 0 ≤ t ≤ T que sigue el sistema (9.38) viene dada por una porci´ on de la par´ abola (9.43) on de la que pase por el punto inicial (x0 , y0 ) seguido de otra porci´ parábola (9.42) que pase por el origen de coordenadas, ver Figura 9.4.

110

´ Control Optimo

Figura 9.4: Trayectorias parab´ olicas con control bang–bang inicial u(t) = −1 y posterior u(t) = 1. La trayectoria ´ optima con condici´ on inicial (x0 , y0 ) está doblemente marcada. Si el control óptimo no tiene discontinuidades y es u∗ (t) = 1 significa que se ha dado la casualidad de que el punto inicial (x0 , y0 ) está incluido en la semi–par´ abola (9.42) que pasa por el origen de coordenadas cuando el tiempo evoluciona, ver Figura 9.5.

Figura 9.5: Trayectoria parab´ olica con control u(t) = 1. Si el control óptimo es de la forma u∗ (t) = −1 significa que el punto abola (9.43) que pasa por inicial (x0 , y0 ) está contenido en la semi–par´ el origen de coordenadas cuando avanza el tiempo, ver Figura 9.6.

Figura 9.6: Trayectoria parab´ olica con control u(t) = −1. 5. Consideremos un oscilador arm´ onico de frecuencia natural w0 = 1 sobre el que act´ ua una fuerza u(t) con |u(t)| ≤ 1. Calcular la fuerza u(t) que hay que aplicar de modo que el oscilador sea transferido al estado final de


111

posici´ on y velocidad nula en el m´ınimo tiempo posible. Interpretar geométricamente el resultado en funci´ on de la posici´ on y velocidad inicial del oscilador. Soluci´ on. Sea x(t) la posición en funci´ on del tiempo del oscilador. Se sabe que el sistema está regido por la ecuaci´ on de segundo orden x ¨ +w02 x = u(t). Tomando w0 = 1 y pasando con el cambio habitual al plano de las fases se tiene x˙ = y , y˙ = −x + u(t) . (9.44) Como se tiene un problema de tiempo óptimo, el Principio del m´ aximo de Pontryagin ser´ a aplicado con el Hamiltoniano (9.31) dado por H = λ1 y + λ2 (−x + u) ,

(9.45)

donde λi verifican las ecuaciones λ˙ 1 = −∂H/∂x, λ˙ 2 = −∂H/∂y. En nuestro caso se tiene (9.46) λ˙ 1 = λ2 , λ˙ 2 = −λ1 , de modo que las variables auxiliares λi verifican la ecuación de un oscilador arm´ onico. Se tiene pues que λ2 (t) = A sin(t − φ) , siendo A > 0 y φ constantes reales (la amplitud y la fase respectivamente). Teniendo en cuenta que |u(t)| ≤ 1, el control ´ optimo u∗ (t) que maximiza el Hamiltoniano (9.45) viene dado por la funci´ on continua a trozos u∗ (t) = sign [λ2 (t)] = sign [sin(t − φ)] , o de forma equivalente −1 si (2k + 1)π ≤ t − φ ≤ (2k + 2)π , k ∈ N ∪ {0} , ∗ u (t) = 1 si 2kπ ≤ t − φ ≤ (2k + 1)π , k ∈ N ∪ {0} . (9.47) Se tiene pues un control u∗ (t) de tipo bang–bang donde las discontinuidades tienen lugar cada intervalo de tiempo π. Estudiemos las trayectorias en el plano x − y de las fases o de estados para el sistema (9.44) dependiendo del control u(t) = u∗ (t) utilizado. Si u(t) = 1, el sistema (9.44) es x˙ = y, y˙ = −x + 1. Eliminando el tiempo se tiene dx/dy = y/(−x + 1) que es una ecuación de variables separables. Se tiene en definitiva que las trayectorias del sistema vienen dadas por (9.48) (x − 1)2 + y 2 = R2 , con R una constante arbitraria. Dicho de otro modo, las trayectorias óptimas son circunferencias de radio R centradas en el punto

´ Control Optimo

112 Y

1

0.5

0.5

1

1.5

2

X

-0.5

-1

Figura 9.7: Trayectorias circulares (x − 1)2 + y 2 = R2 con control u(t) = 1. (x, y) = (1, 0). Además, es fácil ver que las circunferencias (9.48) son recorridas (seg´ un avanza el tiempo) en el sentido de las agujas del reloj con un periodo de 2π. En particular, una semicircunferencia se recorre en el tiempo π, ver Figura 9.7. De forma totalmente análoga al caso anterior, si u(t) = −1, el sistema (9.44) es x˙ = y, y˙ = −x + −1, es decir, dx/dy = y/(−x − 1) de modo que las trayectorias del sistema vienen dadas por (x + 1)2 + y 2 = R2 ,

(9.49)

con R una constante arbitraria. As´ı, las trayectorias óptimas son circunferencias de radio R centradas en el punto (x, y) = (−1, 0). Además, al igual que en el caso anterior, las circunferencias (9.49) son recorridas en el sentido de las agujas del reloj, donde una semicircunferencia es recorrida en tiempo π, ver Figura 9.8. Y

1

0.5

-2

-1.5

-1

X

-0.5

-0.5

-1

Figura 9.8: Trayectorias circulares (x + 1)2 + y 2 = R2 con control u(t) = −1.


113

Veamos cómo son las trayectorias óptimas de este problema dependiendo de las condiciones iniciales (x(0), y(0)) = (x0 , y0 ) ∈ R2 del sistema. Tomaremos el tiempo t1 de modo que el sistema se encuentre en el estado final, en nuestro caso (x(t1 ), y(t1 )) = (0, 0). Supongamos que el u ´ltimo tramo del control o´ptimo u∗ corresponde ´ltimo intervalo de al valor u∗ = 1. Geométricamente, durante este u tiempo, el sistema evoluciona hacia el origen por la circunferencia (9.48) de radio R = 1 partiendo de un punto que denotaremos por A y está localizado en el IV cuadrante. Por supuesto, como es el tramo final y por lo tanto u∗ ya no vuelve a cambiar, el sistema llega sobre esa circunferencia al origen habiendo recorrido un arco inferior a la semicircunferencia. El sistema habr´ a llegado al punto A desplazándose, durante un intervalo de tiempo π y bajo la acci´ on del control u∗ = −1, por la semicircunferencia (9.49) de extremo final A. Llamemos B al extremo inicial de dicha semicircunferencia que, obviamente, estar´ a en el II cuadrante. Notar que, el punto B es el simétrico de A respecto del punto (−1, 0). De forma an´ aloga, la porci´ on de trayectoria ´ optima anterior a la BA, llamémosla CB, consiste en la semicircunferencia (9.48) que pasa por el punto B y el sistema evoluciona en ella durante un tiempo π. Por supuesto, el punto C estará situado en el IV cuadrante y es el simétrico de B respecto del punto (1, 0). Esta construcción geométrica continua hasta que la trayectoria o´ptima pase por el punto inicial (x0 , y0 ), ver Figura 9.9. Y 2

1

-3

-2

-1

1

2

X

-1

-2

-3

-4

Figura 9.9: Una trayectoria o´ptima con el control final u = 1, compuesta de tres arcos de circunferencias

´ Control Optimo

114

Supongamos ahora que el u ´ltimo tramo del control o´ptimo u∗ cor∗ ´ ltimo responde al valor u = −1. Geométricamente, durante este u intervalo de tiempo, el sistema evoluciona hacia el origen por la circunferencia (9.49) de radio R = 1 partiendo de un punto que denotaremos por A y está localizado en el II cuadrante. Por supuesto, como es el tramo final y por lo tanto u∗ ya no vuelve a cambiar, el sistema llega sobre esa circunferencia al origen habiendo recorrido un arco inferior a la semicircunferencia. El sistema habr´ a llegado al punto A desplazándose, durante un intervalo de tiempo π y bajo la acci´ on del control u∗ = 1, por la semicircunferencia (9.48) de extremo final A . Llamemos B al extremo inicial de dicha semicircunferencia que, obviamente, estar´ a en el IV cuadrante. Notar que, el punto B es el simétrico de A respecto del punto (1, 0). De forma an´ aloga, la porci´ on de trayectoria ´ optima anterior a la B A , llamémosla C B , consiste en la semicircunferencia (9.49) que pasa por el punto B y el sistema evoluciona en ella durante un tiempo π. a situado en el II cuadrante y es el Por supuesto, el punto C estar´ simétrico de B respecto del punto (−1, 0). Esta construcción geométrica continua hasta que la trayectoria o´ptima pase por el punto inicial (x0 , y0 ), ver Figura 9.10. Y

3

2

1

-1

1

2

3

X

-1

-2

Figura 9.10: Una trayectoria o´ptima con el control final u = −1, compuesta de tres arcos de circunferencias

Cap´ıtulo 10

Pr´ acticas Este cap´ıtulo está compuesto de varias pr´ acticas propuestas donde se aplica la teor´ıa introducida a lo largo del libro a problemas de diferente naturaleza.

10.1.

Control de la temperatura de una c´ amara

Consideremos el problema de controlar la temperatura Tc en una cámara. El objetivo consiste en analizar un controlador que regule, mediante una v´ alvula, la cantidad de calor que se debe introducir en la c´ amara conociendo la temperatura de esta mediante un term´ ometro de mercurio (Hg) introducido en la c´ amara y contenido en un tubo de cristal de vidrio. Por supuesto, la temperatura de la cámara está influenciada por la temperatura ambiente Ta del entorno. Supongamos adem´ as que la variable que puede ser medida es la altura h de la columna de Hg. Definamos los siguientes flujos de calor: qa es el flujo de calor de la cámara al ambiente. qcr es el flujo de calor de la cámara al crital del term´ ometro. qHg es el flujo de calor del cristal hacia el mercurio. Definamos también las siguientes temperaturas: Tc es la temperatura de la cámara. Ta es la temperatura ambiente. Tcr es la temperatura del cristal del termómetro. THg es la temperatura del mercurio del termómetro. Consideremos las siguientes magnitudes geométricas: 115

116

Pr´ acticas Acr es el área interna del cristal del term´ ometro. Vcr es el volumen del cristal del termómetro. VHg es el volumen del mercurio.

Supongamos que la v´ alvula puede ser regulada mediante un control u(t) de modo que, la cantidad de calor q(t) que se introduce en la cámara viene dada por q(t) = a0 u(t) . Supongamos v´ alidas las siguientes ecuaciones de transferencia de calor qa dTc dt

= a1 (Ta − Tc ) , qcr = a2 (Tc − Tcr ) , qHg = a3 (Tcr − THg ) , dTcr dTHg = b1 (qa + q) , = b2 (qcr − qHg ) , = b3 qHg , dt dt

y las siguientes relaciones de proporcionalidad debidas a las dilataciones por las temperaturas 2 3 3 Acr = c1 Tcr , Vcr = c2 Tcr , VHg = c3 THg . La altura h de la columna de mercurio viene dada por h=

VHg − Vcr . Acr

Se supone que los par´ ametros ai , bj y cj con i = 0, 1, 2, 3 y j = 1, 2, 3 son constantes positivas que dependen de la geometr´ıa y las propiedades de los materiales. Además, en un buen term´ ometro se verificar´ a c2 > c3 .

10.1.1.

Realizaci´ on de la pr´ actica

1. Hallar las ecuaciones del sistema tomando como variables de estado Tc −Ta , ametro Tcr − Ta y THg − Ta , siendo la variable de control u y Ta un par´ constante. 2. Demostrar que si no se ejerce ning´ un control sobre el sistema (es decir u = 0), entonces el u ńico punto de equilibrio x∗ del sistema corresponde al equilibrio termodin´ amico Tc = Tcr = THg = Ta y además es asintóticamente estable en concordancia con el principio cero de la termodin´ amica. 3. ¿Es el sistema completamente controlable para todo valor de ai , bi ? 4. Hallar el control u(t) de modo que el sistema evolucione desde el estado inicial (Tc (0), Tcr (0), THg (0)) = (10, 25, 15) hasta un estado final dado por (Tc (10), Tcr (10), THg (10)) = (20, 25, 25) minimizando la cantidad de calor 10 introducida en la c´ amara, es decir, minimizando 0 u2 (t) dt. Datos: Ta = 20, a0 = 1, a1 = 2, a2 = 3, a3 = 4, b1 = 1, b2 = 2, b3 = 3.

10.2 Din´ amica de sat´ elites de comunicaci´ on

10.2.

117

Din´ amica de sat´ elites de comunicaci´ on

El principal objetivo en la din´ amica de un satélite de comunicación consiste en que dicho satélite permanezca en una posición que sea fija respecto de un cierto punto estratégico de la Tierra. Este tipo de satélites se conocedn como geoestacionarios. Se utilizan, por ejemplo, en comunicaci´ on telefónica, de televisión, navegaci´ on de barcos y aviones, predicci´ on del tiempo, etc... En un principio, una o´rbita geoestacionaria puede ser obtenida mediante la fuerza que ejercen los motores instalados en el satélite. Sin embargo, de este modo, el consumo energético de estos motores ser´ıa continuo y no deseable. Surgen de este modo las siguientes preguntas: ¿Existe una orbita geoestacionaria para el satélite aunque los motores estén parados? Supondremos en este caso, para simplificar, que la u ńica fuerza externa que act´ ua sobre el satélite es la gravitatoria que ejerce la Tierra. ¿Es esta órbita estable? Si no es estable dicha órbita, ¿qué controles mantienen el satélite en dicha órbita?

Figura 10.1: La din´ amica de un satélite terrestre.

10.2.1.

Un modelo matem´ atico para la din´ amica de sat´ elites

Supondremos que sobre el satélite act´ uan las siguientes fuerzas: La fuerza inercial Fin debida al sistema de referencia no inercial. La fuerza gravitacional Fg debida a la atracci´ on de la Tierra. La fuerza Fm ejercida por los motores. Esta será la variable de control.

118

Pr´ acticas La fuerzas perturbativas Fp ejercidas por la Luna, el Sol, el viento solar, etc... Estas fuerzas son perturbativas y el objetivo de los controles será precisamente compensar sus efectos.

La posición del satélite es descrita mediante las coordenadas esféricas (r, φ, θ), siendo r la distancia entre el centro de la Tierra y el centro de masas del satélite. De este modo, si tomamos un sistema de coordenadas cartesianas (x, y, z) con origen en el centro de la Tierra y siendo el plano x − y el plano ecuatorial se tiene x = r cos φ cos θ , y = r cos φ sin θ , z = r sin φ . Estudiaremos la din´ amica del satélite proyectada sobre el plano ecuatorial x − y. En dicho plano, denotamos por er y eθ a los vectores unitarios en las direcciones radial y tangencial. De este modo r = rer . Se sabe que der /dt = dθ/dteθ y deθ /dt = −dθ/dter . De este modo dr dt d2r dt2

= =

dr dθ er + r eθ , dt dt 2 " 2 2 d r dθ d θ dr dθ − r + r + 2 e eθ . r dt2 dt dt2 dt dt

Aplicando las leyes de la din´ amica clásica de Newton, la fuerza inercial Fin viene dada por 2 " 2 d2 r dθ dr dθ d θ d2r −r er + m r 2 + 2 m 2 =m eθ , dt dt2 dt dt dt dt siendo m la masa del satélite. Por otra parte, la fuerza gravitacional Fg es m Fg = −k 2 er , r siendo k = 4 × 104 m3 /s2 el producto de la masa de la Tierra y la constante de gravitaci´ on universal. Adem´ as, definimos Fp = pr er + pθ eθ , Fm = ur er + uθ eθ , siendo ur y uθ las variables de control. En el sistema de referencia no inercial solidario al satélite la fuerza total que act´ ua sobre el satélite es nula de manera que 2 " 2 dθ d θ dr dθ m d2 r − r + m r + 2 e m eθ + k 2 er = Fm + Fp . r dt2 dt dt2 dt dt r Igualando las componentes tangencial y radial en la anterior igualdad se obtiene 2 pr dθ k ur d2 r + , (10.1) = r − 2+ 2 dt dt r m m uθ pθ 2 dr dθ d2 θ + + . (10.2) = − 2 dt r dt dt rm rm

10.3 El p´ endulo doble invertido

10.2.2.

119


1. Supongamos que los motores están parados y que se desprecia la fuerza perturbativa Fp . Demostrar que existe una trayectoria circular geoestacionaria para el satélite. Hallar el radio R de dicha trayectoria. 2. Introducir el cambio de variable θ(t) → ξ(t) = θ(t) − Ωt, siendo Ω la velocidad angular de la Tierra. Reescribir las ecuaciones del movimiento (10.1) y (10.2) del satélite con las variables (r, ξ). Verificar que (r∗ , ξ ∗ ) = (R, 0) es un punto de equilibrio del sistema cuando los motores est´ an parados y que se desprecia la fuerza perturbativa Fp . 3. Hallar el sistema de control lineal asociado al satélite en un entorno del anterior punto equilibrio si se desprecia la fuerza perturbativa Fp . 4. ¿Es el anterior sistema linearizado completamente controlable? 5. Supongamos que uθ ≡ 0, es decir, el motor sólo puede actuar sobre el satélite en la dirección radial. ¿Es posible hallar, en el sistema linearizado del apartado 3, un control por realimentaci´ on lineal de modo que el conjunto de valores propios en lazo cerrado sea el que deseemos?

10.3.

El p´ endulo doble invertido

Consideremos un péndulo doble invertido montado sobre el centro de masas de un carrito que puede desplazarse sobre una recta horizontal. Definimos las siguientes magnitudes. m es la masa del carrito. m1 y m2 son las masas de los dos péndulos. 1 y 2 son las longitudes de los dos péndulos. Definimos las variables siguientes. u es la fuerza externa horizontal que act´ ua sobre el carrito. q es la posición del centro de masas del carrito. θ1 y θ2 son los ángulos de los dos péndulos respecto de la vertical. Es f´ acil demostrar que la energ´ıa cinética K y la energ´ıa potencial U del sistema viene dada por

$ %2 $ %2 1 1 K = mx˙ 2 + m1 q˙ + 1 θ˙1 cos θ1 + 1 θ˙1 sin θ1 2 2

$ %2 1 + m2 q˙ + 1 θ˙1 cos θ1 + 2 θ˙2 cos(θ1 + θ2 ) 2 $ %2 + 1 θ˙1 sin θ1 + 2 θ˙2 sin(θ1 + θ2 ) , U

= m1 g1 cos θ1 + m2 g[1 cos θ1 + 2 cos(θ1 + θ2 )] .

120

Pr´ acticas

Figura 10.2: El péndulo doble invertido. Utilizando el Lagrangiano del sistema L = K−U , las ecuaciones del movimiento se obtienen mediante las ecuaciones de Euler-Lagrange ∂L d ∂L ∂L d ∂L d ∂L ∂L =0. − =0, − =u, − ˙ ˙ ∂θ2 dt ∂ θ2 ∂θ1 dt ∂ θ1 ∂q dt ∂ q˙ Desarrollando estas ecuaciones se obtienen las siguientes ecuaciones del movimento del doble péndulo invertido = −(1 m1 + 1 m2 )(θ˙1 )2 sin θ1 − 2 m2 θ˙1 θ˙2 sin(θ1 + θ2 ) q (10.3) −2 m2 (θ˙2 )2 sin(θ1 + θ2 ) + (m + m1 + m2 )¨ +(1 m1 + 1 m2 )θ¨1 cos θ1 + 2 m2 θ¨2 cos(θ1 + θ2 ) , 0 = −g1 (m1 + m2 ) sin θ1 − g2 m2 sin(θ1 + θ2 ) + 2 m2 q˙θ˙2 sin(θ1 + θ2 ) q cos θ1 + 21 (m1 + m2 )θ¨1 (10.4) −1 2 m2 (θ˙2 )2 sin θ2 + (1 m1 + 1 m2 )¨ +1 2 m2 θ¨2 cos θ2 , 0 = −g2 m2 sin(θ1 + θ2 ) − 2 m2 q˙θ˙1 sin(θ1 + θ2 ) + 2 m2 q¨ cos(θ1 + θ2 ) +1 2 m2 θ¨1 cos θ2 + 22 m2 θ¨2 . (10.5)

u

10.3.1.


1. Utilizar las variables de estado ˙ θ˙1 , θ˙2 )T ∈ R6 x = (x1 , x2 , x3 , x4 , x5 , x6 )T = (q, θ1 , θ2 , q, y reescribir las ecuaciones del movimiento como un sistema de primer orden de la forma x˙ = f (x; u). 2. Demostrar que x∗ = (q ∗ , θ1∗ , θ2∗ , q˙∗ , θ˙1∗ , θ˙2∗ ) = (0, 0, 0, 0, 0, 0) es un punto cr´ıtico del sistema para u = u∗ = 0. Darle un significado f´ısico. 3. Obtener el sistema de control lineal asociado al péndulo doble invertido en un entorno del punto cr´ıtico (x∗ ; u∗ ).

10.4 Estabilidad en el p´ endulo c´ onico

121

4. Averiguar para qué valores de los parámetros m, m1 , m2 , 1 y 2 es el anterior sistema linearizado completamente controlable. 5. Se pretende estabilizar el sistema linearizado usando realimentaci´ on lineal. Hallar dicha realimentaci´ on de modo que los valores propios de la matriz del sistema en lazo cerrado sean ρi = −i, con i = 1, 2, 3, 4, 5, 6. Notar que todos los valores propios tienen parte real negativa y por lo tanto se estabiliza el sistema hacia x∗ .

10.4.

Estabilidad en el p´ endulo c´ onico

Consideremos un péndulo de masa m y longitud sometido a un campo gravitatorio de intensidad constante g y suspendido de un punto fijo O. La posición de la masa m queda determinada a través de los ángulos θ y φ. θ es el ángulo que forma el péndulo respecto de la vertical que pasa por el punto O. φ es un ángulo formado por dos rectas contenidas en el plano horizontal donde se encuentra la masa m. Una recta de referencia con dirección fija y la recta que pasa por m y es perpendicular al eje vertical que pasa por O.

Figura 10.3: (a) El péndulo cónico con su masa realizando una trayectoria circular uniforme. (b) Otra trayectoria cercana a la anterior. Es fácil demostrar que la energ´ıa cinética K y la energ´ıa potencial U del sistema viene dada por K U

1 2 ˙2 ˙ 2 sin2 θ] , m [θ + phi 2 = mg[1 − cos θ] . =

122

Pr´ acticas

Utilizando el Lagrangiano del sistema L = K −U , las ecuaciones del movimiento se obtienen mediante las ecuaciones de Euler-Lagrange d ∂L ∂L − = 0 , q = θ, φ . dt ∂ q˙ ∂q Desarrollando estas ecuaciones se obtienen las siguientes ecuaciones del movimento del péndulo c´ onico g (10.6) θ¨ = − sin θ + φ˙ 2 sin θ cos θ , φ¨ = −2θ˙φ˙ cot θ . (10.7)

10.4.1.


1. Utilizar las variables de estado ˙ φ) ˙ T ∈ R3 x = (x1 , x2 , x3 )T = (θ, θ, y reescribir las ecuaciones del movimiento como un sistema de primer orden de la forma x˙ = f (x). 2. Supongamos que la masa m describe un trayectoria circular en un plano horizontal a velocidad angular w constante cuando el péndulo tiene un ángulo θ = α constante. Averiguar qué relación deben verificarse entre w y α para que dicho movimiento sea posible. ˙ φ) ˙ = (α, 0, w) es un punto cr´ıtico del péndulo 3. Demostrar que x∗ = (θ, θ, cónico si se verifica la condición obtenida en el apartado anterior. 4. Realizar una traslaci´ on de coordenadas de modo que el punto cr´ıtico x∗ se traslade al punto cr´ıtico y ∗ = 0 ∈ R3 en las nuevas coordenadas. ¿Es posible estudiar la estabilidad del origen por linearizaci´ on bajo la condici´ on del apartado 2? 5. Calculando derivadas orbitales, demostrar que el péndulo c´ onico tiene las siguientes constantes del movimiento H1 (y1 , y2 , y3 )

= [y22 + (y3 + w)2 sin2 (y1 + α)] −

H2 (y1 , y2 , y3 )

= (y3 + w) sin2 (y1 + α) .

2g cos(y1 + α) ,

Nota: Dichas constantes del movimiento provienen f´ısicamente de la conservación de la energ´ıa y del momento angular. 6. Hallar un valor de la constante λ ∈ R de modo que la funci´ on V (y1 , y2 , y3 ) = H1 (y1 , y2 , y3 ) − H1 (0, 0, 0) + λ[H2 (y1 , y2 , y3 ) − H2 (0, 0, 0)] , sea una funci´ on de Liapunov en un entorno del origen cuando se verifica la relación del apartado 2. ¿Qué se puede concluir sobre la estabilidad del movimiento del péndulo c´ onico del apartado 2?

10.5 Un modelo de gr´ ua con servomecanismo y regulador

10.5.

123

Un modelo de gr´ ua con servomecanismo y regulador

Un modelo simple de gr´ ua consiste en tomar una masa m1 realizando un movimiento rectil´ıneo que representa el punto de suspensi´ on de un péndulo simon donde ple de longitud y masa m2 . Se requiere que el punto de suspensi´ on final de modo que está m1 se traslade desde su posición inicial hasta su posici´ las oscilaciones del péndulo sean m´ınimas. Tomemos las siguientes coordenadas del sistema gr´ ua con dos grados de libertad: x(t) ∈ R posición en funci´ on del tiempo de la masa m1 a lo largo de un eje horizontal. θ(t) será el ángulo respecto de la vertical en funci´ on del tiempo que forma el péndulo. Utilizando el Lagrangiano del sistema gr´ ua L = K − U , siendo K y U la energ´ıa cinética y potencial del sistema respectivamente, es fácil mostrar que ˙ = L(x, x, ˙ θ, θ)

1 1 (m1 + m2 )x˙ 2 + m2 2 θ˙2 cos2 θ + m2 x˙ θ˙ cos θ + m2 g cos θ . 2 2

Mediante las ecuaciones de Euler-Lagrange d ∂L ∂L d ∂L ∂L − = u(t) , =0, − dt ∂ x˙ ∂x dt ∂ θ˙ ∂θ siendo u(t) la fuerza horizontal aplicada a la masa m1 , se obtienen las siguientes ecuaciones del movimento de la gr´ ua x + m2 (θ¨ cos θ − θ˙2 sin θ) , u = (m1 + m2 )¨ 0 = θ¨ + x ¨ cos θ + g sin θ .

(10.8) (10.9)

Si suponemos que el péndulo realiza peque¯ nas oscilaciones, es decir, sin θ ≈ θ y cos θ ≈ 1, y además su velocidad angular es peque¯ na con lo cual θ˙ ≈ 0, entonces las ecuaciones del movimiento se linearizan de la forma m2 x ¨, (10.10) u = m2 θ¨ + µ 0 = θ¨ + x ¨ + gθ , (10.11) siendo µ := m2 /(m1 + m2 ).

10.5.1.


1. Escribir las ecuaciones (10.10) y (10.11) de la gr´ ua en la forma de sistema de control lineal y˙ = Ay + bu tomando como variables de estado y = ˙ x, x) (θ, θ, ˙ T ∈ R4 .

124

Pr´ acticas

2. Se desea que el puente de la gr´ ua donde est´ a m1 se desplace desde el origen de coordenadas en el eje x hasta una distancia D en un tiempo t1 . Además, en el inicio y el fin del movimiento se supone velocidad nula. Suponiendo una ley polinomial para el movimiento de referencia z(t) de m1 , hallar un polinomio de grado m´ınimo adecuado. 3. El sistema de gr´ ua se compone de un servomecanismo (se desea que m1 se mueva seg´ un x(t) siguiendo lo m´ as cerca posible a z(t)) y un regulador (se quiere que el péndulo tenga m´ınima oscilación de modo que se quiere mantener θ y θ˙ tan pr´ oximo a cero como se pueda). De este modo, se puede tomar la se¯ nal de referencia r(t) para el sistema lineal del apartado T ∈ R4 y el siguiente ´ındice de actuación 1 como r(t) = (0, 0, z(t), z(t)) ˙ para minimizar 1 1 t1 T J(u, x) = ∆T (t1 )M ∆(t1 ) + [∆ (t)Q∆(t) + αu2 (t)] dt , 2 2 0 siendo ∆(t) = y(t) − r(t) ∈ R4 , α ∈ R un par´ ametro, M, Q ∈ M4 (R) simétricas y definidas positivas. Calcular el control o´ptimo u∗ (t) por realimentaci´ on lineal y la trayectoria o´ptima y ∗ (t). Utilizar para las integraciones numéricas de ecuaciones diferenciales o bien el método RK4 con longitud de paso h = 0,1 o directamente el comando N DSolve[] del programa Mathematica. Datos: t1 = 20, D = 6, = 10, m1 = m2 = 1. Tomar, para simplificar cálculos, las matrices identidad M = Q = I4 y α = 1. 4. A la vista de los resultados del apartado anterior, dar las gr´ aficas del on x(t) y la control ´ optimo u∗ (t), la de las superposiciones de la posici´ nal de referencia z(t) y z(t) ˙ y finalmemte la velocidad x(t) ˙ de m1 con la se¯ del ángulo θ(t). 5. Seg´ un las gr´ aficas del apartado anterior, decir si son verdaderas o falsas las siguientes afirmaciones: Esencialmente, se empuja a m1 en la dirección positiva del eje x durante la mitad del recorrido y se frena despues. Durante la fase en que m1 es acelerada, m2 queda detr´ as y en la segunda fase en que m2 se frena entonces m2 adelanta a m1 . 6. Sea P (t) = (pij ) ∈ M4 (R) la matriz que satisface la ecuación diferencial matricial de Riccati asociada al problema de optimizaci´ on. Dada la norma . 2 aficamente que de P (t) definida como P (t) = i,j pij (t), mostrar gr´ se mantiene prácticamente constante durante una larga parte del recorrido óptimo de la g´ ua. Esto justifica que, en muchas ocasiones se sustituya la ecuación diferencial matricial de Riccati por la ecuaci´ on matricial algebraica de Riccati aunque el proceso sea controlado durante un tiempo finito.

10.6 El regulador centr´ıfugo de Watt

10.6.

125

El regulador centr´ıfugo de Watt

El regulador centr´ıfugo, ideado por Watt a finales del siglo XVIII para ser utilizado en las máquinas de vapor, fue uno de los primeros mecanismos utilizados para el control autom´ atico.

Figura 10.4: El regulador centr´ıfugo de Watt. El regulador centr´ıfugo consiste en una barra vertical S que gira sobre su propio eje. En su extremos superior hay dos varillas articuladas L1 y L2 de la misma longitud y con una masa m de carga en cada uno de sus extremos libres. L1 y L2 se encuentran en el mismo plano vertical que S y sólo pueden separarse de su posición de equilibrio de forma simultánea y formando el mismo ángulo ϕ respecto de la vertical S. En este deslazamiento las varillas ponen en movimiento, mediante articulaciones, un manguito M colocado sobre la barra S. Consideremos una de las masas m. Seg´ un la geometr´ıa descrita anteriormente, dicha masa se desplazará sobre una esfera de radio la longitud de una varilla. Tomaremos dicha longitud como la unidad de longitud de modo que, tomando el origen de coordenadas en el punto de uni´ on entre la barra S y las varillas L1 y L2 , siendo el eje z en la dirección de la barra S y con sentido positivo hacia abajo, la posici´ on en funci´ on del tiempo (x(t), y(t), z(t)) de la carga m viene dada, en coordenadas esféricas, por (1, θ(t), ϕ(t)), siendo x = sin ϕ cos θ , y = sin ϕ sin θ , z = cos ϕ . Las fuerzas que act´ uan sobre m son básicamente las siguientes: El peso FP = mg siendo g la gravedad supuesta constante. Una fuerza de tensi´ on FT que mantiene m sujeta a la varilla.

126

Pr´ acticas Una fuerza de rozamiento Fr debida a los engranajes de la varilla que se opone a la variaci´ on de ángulo ϕ de modo que su módulo se puede tomar Fr = k ϕ˙ siendo k una constante positiva. Una fuerza, llamemosla FB , responsable del giro del regulador sobre su eje. Aplicando la segunda ley de Newton, es f´ acil ver que se verifica la siguiente ecuación mϕ¨ = mθ˙2 sin ϕ cos ϕ − mg sin ϕ − k ϕ˙ . Por otra parte, la m´ aquina de vapor es, simplificando, un volante con un momento de inercia J que se mueve debido a la fuerza del vapor y es capaz de realizar un trabajo, por ejemplo mover un peso. Si denotamos por P1 el momento que realizan las fuerzas del vapor sobre el volante y por P el momento de dicho peso, utilizando las leyes de la mecánica clásica y en particular las del s´ olido r´ıgido se tiene que J w˙ = P1 − P , siendo w la velocidad angular del volante. El volante de la m´ aquina est´ a conectado a la barra vertical S del regulador mendiante engranajes de modo que existe una relaci´ on constante θ˙ = nw entre las velocidades angulares θ˙ y w de la barra y el volante respectivamente. Además, el mangito M del regulador est´ a conectado a la válvula que administra el vapor de modo que P1 = F1 +c cos ϕ, siendo F1 y c constantes positivas. El momento P1 depende de la obertura que tenga la v´ alvula que regula la entrada de vapor. El regulador centr´ıfugo funciona de la manera siguiente: si la velocidad angular w del volante es muy alta, la barra S gira muy r´ apido y, debido a la fuerza centr´ıfuga, las varillas L1 y L2 se abren un ´ angulo ϕ superior lo cual hace que la obertura de la v´ alvula disminuya y la cantidad de vapor que entra en la m´ aquina decrezca. Si, por el contrario, la velocidad angular del volante w disminuye el efecto contrario es proporcionado. Seg´ un lo descrito anteriormente, parece natural esperar que la velocidad angular w del volante tienda a estabilizarse, y de hecho as´ı ocurri´ o hasta mediados del siglo XIX. No fue hasta que el ingeniero ruso Vichnegradski (1876) solvent´ o el problema de dicha inestabilidad en un trabajo te´ orico que reproducimos a continuaci´ on.

10.6.1.


a) Tomar las variables de estado (ϕ, Ω, w) con Ω = ϕ˙ y reescribir las ecuaciones del movimiento del regulador como un sistema de primer orden.

10.6 El regulador centr´ıfugo de Watt

127

b) Si se desea que el regulador cumpla su misi´ on debe funcionar lo m´ as cercano posible de un punto cr´ıtico del sistema con ϕ = ϕ∗ constante. Averiguar si existe dicho punto. c) Con el paso del tiempo y gracias a las mejoras técnicas, se fueron modificando diversas magnitudes del regulador. Por un lado, el aumento de la potencia de las máquinas requer´ıa válvulas más pesadas con lo que m aument´ o. La constante de rozamiento k disminuy´ o debido a que las piezas estaban cada vez más pulidas. También se redujo el tama¯ no de los volantes y, en consecuencia, J disminuy´ o. Sin embargo, parad´ ojicamente, los reguladores funcionaron peor con dichas modificaciones. Averiguar la causa mediante un an´ alisis de la estabilidad del punto cr´ıtico del apartado anterior.

Bibliograf´ıa [1] S. Barnet & R.G. Cameron. Introduction to mathematical control theory. Claredon Press. Oxford, 1985. [2] R.W. Brockett. Finite dimensional linear systems. Wiley, New York, Dover, 1970. [3] E.J. Davison & S.H. Wang. I.E.E.E. Trans. Autom. Control, 20, 516. [4] L. Elsgoltz. Ecuaciones diferenciales y c´ alculo variacional. Ed. Mir, Moscou, 1992. [5] B.C. Kuo. Sistemas de control autom´ atico. Prentice Hall, 1996. [6] H. Kwakernaak & R. Sivan. Linear optimal control theory. Wiley– Interscience, New York, 1972. [7] L.S. Pontriaguin. Ecuaciones diferenciales ordinarias. Aguilar, Madrid, 1973. [8] L.S. Pontriaguin, V. Boltianski, R. Gamkrelidze & E. Michtchenko. The mathematical theory of optimal processes. Nueva York: Wiley, 1962. [9] J.W. Ponlderman & J.C. Willens. Introduction to mathematical systems theory. Texts in Applied Mathematics 26, Springer. [10] B. Porter & R. Crossley. Modal control. Taylor and Francis, London, 1972. [11] H.H. Rosenbrock. State–space and multivariable theory. Nelson, London, 1970. [12] J.L. Willems. Stability theory of dynamical systems. Nelson. London, 1970. [13] H.H. Wonham. Linear multivariable control: a geometric approach, 2nd edn. Springer-Verlag. Berlin, 1979.

129

Teoría de Estabilidad y Control Isaac A. García

Recommend Documents