Pucheta Acse 2009

UNIVERSIDAD NACIONAL DE CÓRDOBA FACULTAD DE CIENCIAS EXACTAS, FÍSICAS Y NATURALES

CONTROL ÓPTIMO Y SISTEMAS ESTOCÁSTICOS

Filminas para apuntes de clases a cargo de

Prof. Dr. Ing. Julián A. Pucheta

2009

Control óptimo y procesos estocásticos

Contenidos 1. INTRODUCCIÓN .......................................................................................................................................................... 3 1.1. Modelo en el espacio de estado ......................................................................................................................... 3 1.2. Diseño de controladores de estado lineales ....................................................................................................... 6 1.3. Esquema básico del controlador lineal de estado .............................................................................................. 7 1.4. Metodologías de diseño más utilizadas ............................................................................................................. 8 1.5. Diseño del controlador mediante asignación de polos ...................................................................................... 8 1.6. Controlador de tiempo finito ............................................................................................................................. 9 2. CONTROL ÓPTIMO EN SISTEMAS LINEALES ................................................................................................................. 9 2.1. Motivación ........................................................................................................................................................ 9 3. REGULADOR ÓPTIMO LINEAL EN TIEMPO CONTINUO ................................................................................................ 10 3.1. Formulación del problema............................................................................................................................... 10 3.2. Estabilidad en el sentido de Lyapunov ............................................................................................................ 10 3.3. Problema de control óptimo cuadrático........................................................................................................... 12 4. REGULADOR ÓPTIMO LINEAL EN TIEMPO DISCRETO ................................................................................................. 14 4.1. Formulación del problema............................................................................................................................... 14 4.2. Formulación del problema de estado estacionario........................................................................................... 14 4.3. Problema de control óptimo lineal de continuo a discreto............................................................................... 19 5. REGULADOR ÓPTIMO LINEAL EN EL TRANSITORIO .................................................................................................... 28 5.1. Formulación del problema en el transitorio..................................................................................................... 28 6. CONTROL ÓPTIMO BASADO EN PROGRAMACIÓN DINÁMICA ...................................................................................... 33 6.1. Principio de optimalidad de Bellman .............................................................................................................. 33 7. PROGRAMACIÓN DINÁMICA ..................................................................................................................................... 40 7.1. Versión simbólica: Ecuación de Hamilton-Jacobi-Bellman............................................................................ 40 7.2. Versión numérica: Ecuación de Bellman ........................................................................................................ 44 7.3. Problema básico .............................................................................................................................................. 44 7.4. La política óptima de decisiones ..................................................................................................................... 45 7.5. Programación dinámica regresiva ................................................................................................................... 45 7.6. Algunos funcionales típicos ............................................................................................................................ 49 7.7. Programación Dinámica iterativa .................................................................................................................... 51 7.8. Programación dinámica aproximada ............................................................................................................... 52 7.9. Discusión y comentario final........................................................................................................................... 61 8. CONTROL DIGITAL ESTOCÁSTICO ................................................................................................................. 62 8.1. Modelo matemático estocástico de señales reales. .......................................................................................... 62 8.2. Ecuaciones diferenciales estocásticas.............................................................................................................. 63 8.3. Modelos de Estado para Sistemas Estocásticos de Tiempo continuo.............................................................. 65 8.4. Modelos de Estado para Sistemas Estocásticos de Tiempo Discreto. ............................................................. 71 8.5. Diseño de Controladores de Estado para Sistemas Estocásticos Lineales....................................................... 78 9. BIBLIOGRAFÍA .......................................................................................................................................................... 87

J. A. Pucheta (www.labimac.blogspot.com)

2

Control óptimo y procesos estocásticos

1. Introducción 1.1. Modelo en el espacio de estado La representación de sistemas en el espacio de estado constituye una herramienta de gran utilidad para el análisis y diseño de sistemas de control en el dominio temporal. En particular resulta de gran significación para el tratamiento de los sistemas multivariable. Esta forma de representación fue desarrollada para el tratamiento de modelos continuos y fue extendida posteriormente a los modelos discretos en razón de los requerimientos impuestos por el control digital. Se puede dar informalmente para un sistema la siguiente definición de estado dinámico del sistema. El estado de un sistema causal, es la información mínima que es necesario conocer en un instante t=t0 para que conjuntamente con el valor de las entradas definidas en todo tiempo a partir de t≥t0; se pueda determinar el comportamiento del sistema para cualquier t≥t0. El estado dinámico de un sistema constituye una información instantánea que se va modificando con la evolución temporal del sistema. Las variables que son necesarias para definir el estado se denominan variables de estado. Se puede dar la siguiente definición. Las variables de estado constituyen el conjunto más pequeño de variables, tales que el conocimiento de las mismas en t=t0, conjuntamente con las entradas para t≥t0, determinan el comportamiento del sistema para cualquier tiempo t≥t0. De igual modo se puede definir el vector de estado como: Un vector de estado de dimensión n es aquél cuyas componentes están constituidas por las n variables de estado. Finalmente se define al espacio de estado de la siguiente manera. Espacio de estado es el espacio geométrico n-dimensional donde se puede representar cualquier estado por un punto. Con el objeto de asociar estas definiciones a la modelación de un sistema físico, se toma como ejemplo un circuito elemental RLC; representado en la Fig. 1-1.

Fig. 1-1. (a) Circuito RLC; (b) Entrada-Salida del circuito RLC

Se toma u=ve(t) como señal de entrada al sistema y la tensión vr(t) sobre el resistor R como salida. Por relaciones físicas es conocido que la evolución de las distintas variables físicas en este circuito, tales como tensiones y corrientes, quedará definida en un futuro si se conoce para un instante de tiempo t=t0, la corriente que fluye en el inductor L, la tensión que exista sobre el capacitor C y la tensión de entrada desde t0 en adelante. En base a la definición que se ha dado de variables de estado es posible elegir a la corriente en el J. A. Pucheta (www.labimac.blogspot.com)

3

Control óptimo y procesos estocásticos

circuito y a la tensión sobre el capacitor como variables de estado, ya que éstas definen el estado dinámico del circuito. La evolución futura del estado dinámico para t≥t0 se podrá determinar si se conoce para t=t0 las variables de estado i(t), vc(t) y además la tensión de entrada ve(t) para t≥t0. Para analizar la evolución del circuito se pueden plantear las ecuaciones diferenciales del mismo. di R 1 1 = − i − vc + ve dt L L L (1-1) dvc = 1 i dt C Las Ec. (1-1) se pueden expresar en una ecuación matricial-vectorial.  R  di  −  dt   L  dv  =   c  1  dt    C

1 −  L  i  1   +  L  [ ].   v    ve 0  c   0 

(1-2)

Definiendo a i, vc como variables de estado y a x como vector de estado, la Ec. (1-2) se tiene x& = A x(t) + b u(t) con

 - R/L - 1/L  1/L  A=  , b =  . 0  1/C  0

(1-3)

La matriz A se denomina matriz del sistema y b vector de entrada. La variable de salida y=vR puede obtenerse también a partir del vector de estado mediante y = cT x(t)

(1-4)

con cT=[R 0]. El vector C se denomina vector de salida. De esta forma el circuito RLC de la Fig. 1-1 queda modelado en el espacio de estado por x& (t) = Ax(t) + b u(t)

(1-5)

y(t) = c T x(t)

con T x (t) = [i v c ] u(t) = ve (t) y(t) = vr (t). estando A, b, c definidas en las Ec. (1-3) y Ec. (1-4).

1.1.1. Representación de sistemas multivariables Cuando se consideran varias entradas y varias salidas del sistema simultáneamente, se recurre a la J. A. Pucheta (www.labimac.blogspot.com)

4

Control óptimo y procesos estocásticos

representación mostrado en la Fig. 1-2, en el cual existen interacciones múltiples de las e entradas con las s salidas. Si se desea modelar con ecuaciones diferenciales, conduce a un sistema de s×e ecuaciones diferenciales, de distinto orden que contemplan las relaciones dinámicas de todas las entradas con las distintas salidas. La de mayor orden define el orden n del sistema multivariable. Además, el orden del sistema está dado por el número mínimo de variables de estado necesarias para describir la evolución del sistema.

Fig. 1-2. Sistema Multivariable.

El sistema de las s×e ecuaciones diferenciales transformadas al dominio de la frecuencia en variable compleja s permiten modelar el sistema multivariable a través de la matriz de transferencia G(s),

y(s ) = G(s ) u (s )

(1-6)

donde y(s) es el vector de salida de dimensión s, u(s) es el vector de entrada de dimensión e, y G(s) es la matriz de transferencia de dimensión s×e. Cada elemento de la matriz G(s) representa la Función de Transferencia Gij(s) de la entrada uj(s) respecto de la salida yi(s). De la misma forma que para el caso monovariable, aunque con un mayor grado de complejidad resulta posible a través de una adecuada elección de las variables de estado, transformar todas las ecuaciones diferenciales en conjuntos de ecuaciones diferenciales de primer orden, y compactar la notación para obtener una ecuación diferencial matricial-vectorial de primer orden de la misma forma que las Ec. (1-5),  x& (t) = A(t) x(t) + B(t) u(t)   y(t) = C(t) x(t) + D(t) u(t).

(1-7)

Donde A es la matriz del sistema, B es la matriz de entrada, C es la matriz de salida y D es la matriz de transferencia directa, todas expresadas como función del tiempo para señalar la dependencia temporal en el caso de ser necesario. Para determinar la correcta dimensión de las distintas matrices componentes de la Ec. (1-7), resulta útil representar los vectores y matrices de la Ec. (1-7) por rectángulos cuyas longitudes de lados representan la dimensión considerada. Las Ec. (1-7) pueden representarse esquemáticamente para un sistema multivariable con e entradas y s salidas como en la Fig. 1-3.

J. A. Pucheta (www.labimac.blogspot.com)

5

Control óptimo y procesos estocásticos

Fig. 1-3. Representación esquemática de las ecuaciones de estado.

Se observa que para un sistema multivariable la matriz de entrada B toma la dimensión n×e, la matriz de salida C la dimensión s×n, la matriz de transferencia directa D la dimensión s×e y la matriz de entrada A, la dimensión n×n, igual que para el caso monovariable.

1.2. Diseño de controladores de estado lineales

Existen diversos esquemas de control, basados en la teoría de Entrada-Salida y en la de variables de estado. A continuación se muestran los esquemas más difundidos. 1.1.2. Entrada-salida Se realimenta el error de control, definido como e(k)=yd(k)-y(k). Los esquemas más difundidos son los del tipo Proporcional Integral Derivativo PID, con sus diversas variantes, por ejemplo, Modificado, con predictor, con anti-wind up, auto sintonía, etc.

yd ek Controlador

uk

Proceso

yk

-

Fig. 1-4. Esquema de control en la representación de sistemas Entrada-Salida.

J. A. Pucheta (www.labimac.blogspot.com)

6

Control óptimo y procesos estocásticos

1.1.3. Espacio de estados Se realimenta el estado del proceso, x(k). yk

uk Proceso

rk -

ek Controlador Fig. 1-5. Esquema de control basado en realimentación de estados.

1.3. Esquema básico del controlador lineal de estado Se modifica el funcionamiento del proceso mediante el controlador, implementado como u(k)=-k(k)x(k). Modelo de proceso lineal en el espacio de estados con realimentación lineal x(k + 1) = (A (k ) − B(k )K (k ))x(k )  y (k ) = (C(k ) − D(k )K (k ))x(k )

D(k) x(0) u(k)

x(k+1) B(k)

Iq

-1

y(k)

x(k) C(k)

A(k)

-K(k)

J. A. Pucheta (www.labimac.blogspot.com)

7

Control óptimo y procesos estocásticos

1.4. Metodologías de diseño más utilizadas En el ámbito de diseño de controladores de estado lineales, los más difundidos son:

Asignación de polos

a) Ecuación característica dada b) De tiempo finito

Espacio Funcional de costo de estados Regulador óptimo lineal

1.5. Diseño del controlador mediante asignación de polos El objetivo es trasladar los polos de det(zI-A)=0 a un lugar deseado mediante det (zI-A+BK)=0. Ecuación característica dada Para el caso monovariable, se hacen transformaciones lineales en el sistema para obtener la forma canónica controlable 0  0 1 0 K   0  0 0 1 K ~  A =  0 0 0 K 0  . . . K .   - a n . . K - a1

0    0  ~   ~T b = 0 c = [b m K bo 0 K 0]    .    1

1 0 K 0  0   0  0 1 K 0  0 .  w (k + 1) =  0 0 0 K 0 w (k ) +   u (k ) .  .   . . K .   1  - a n − a n −1 . K - a1

Con u(k)=-K.w(k) 0   . w (k + 1) =   .  − a n − k n

La ecuación característica es

J. A. Pucheta (www.labimac.blogspot.com)

1 . . − a n −1 − k n −1

   w (k )  . .  ... − a 1 − k 1  ...

0

1

0

(a n + k n ) + (a n −1 + k n −1 )z + ... + (a 1 + k 1 )z n

=0. 8

Control óptimo y procesos estocásticos

Las raíces de ésta ecuación característica será

(z − p 1 )(z − p 2 )...(z − p n ) = 0 . El diseño del controlador comienza ubicando a los polos pi, que son los polos de lazo cerrado. Tiene como ventaja de que es un método simple y directo. La desventaja es que no se tiene en cuenta el efecto conjunto de los polos en el comportamiento del sistema, ni tampoco la magnitud de las acciones de control.

1.6. Controlador de tiempo finito Es un caso particular del anterior, donde se sitúa a todos los polos de lazo cerrado en el origen del plano complejo. Se llega a ki=ai por lo tanto u será u (k ) = [a n

a n −1 ... a 1 ]x(k ) .

La ventaja de éste método es que es rápido, simple, y considera el efecto del conjunto de los polos de lazo cerrado. La desventaja es que las acciones de control son muy elevadas.

2. Control óptimo en sistemas lineales 2.1. Motivación Para el caso del diseño de controladores en tiempo discreto, siempre que el sistema sea controlable, los polos de lazo cerrado pueden ubicarse en cualquier punto del plano complejo, pero el límite de las respuestas está dado por las acciones de control. La velocidad de respuesta del proceso y la magnitud de las acciones de control están inversamente relacionadas. Una solución se encuentra proponiendo un funcional de costo que incluya estos elementos y luego realizar su minimización. El uso del criterio de minimización es ampliamente utilizado, incluso en controladores EntradaSalida para optimizar los parámetros del PID, Fletcher-Powell propone J = J (e, u ) =

∑ [e ( ) M

2

k

k =0

+ r∆u (k )

2

]

La solución analítica es posible solamente para controladores de bajo orden. Se puede minimizar a prueba y error numéricamente. Para el diseño en el espacio de estados, se empleará el funcional J = J (e, u ) =

∑ [e ( ) M

k

k =0

J. A. Pucheta (www.labimac.blogspot.com)

2

+ ru (k )

2

] 9

Control óptimo y procesos estocásticos

El funcional propuesto es convexo y continuo de sus argumentos ek y uk. El controlador será óptimo en el sentido de éste funcional.

3. Regulador óptimo lineal en tiempo continuo 3.1. Formulación del problema Dado el sistema lineal determinístico en tiempo continuo x& t = A ⋅ x t + B ⋅ u t  y t = C ⋅ x t

(3-1)

u t = −K ⋅ x t

(3-2)

se desea encontrar una ley de control ut que haga evolucionar al proceso desde x(0)≠0 a x(∞)=0 minimizando el funcional de costo J (x, u ) =

∫ (x

∞

T

)

(3-3)

Qx + u T Ru dt

0

con Q simétrica y semidefinida positiva y R simétrica y definida positiva.

Para diseñar el controlador en espacio de estados en el domino del tiempo continuo, se usará el segundo método de Lyapunov, porque no requiere resolver las ecuaciones diferenciales del sistema a controlar. En general, el sistema se define como x& = f (x , t ), x (0 ) = x 0

con la solución

φ(t 0 ; x 0 , t 0 ) = x 0 .

3.2. Estabilidad en el sentido de Lyapunov Definición de equilibrio: Es un punto del espacio de estado, xe, donde

f (x e , t ) = 0 ∀t ≥ 0.

Los sistemas lineales tienen el origen como único punto de equilibrio si A es no singular. Se definen las esferas S(δ) y S(ε) alrededor del punto de equilibrio, ∀t≥t0. mediante x 0 − x e ≤ δ,

φ(t; x 0 , t 0 ) − x e ≤ ε,

respectivamente. El sistema será estable en el sentido de Lyapunov si para cada esfera S(ε) existe una esfera S(δ) tal que las trayectorias que empiezan en S(δ) no salen de S(ε) con t→∞. Si δ no depende de t0, el equilibrio es uniformemente estable. El sistema será inestable si para algún ε>0 y cualquier δ>0 siempre existirá un x0 en S(δ) tal que las trayectorias que allí comienzan se salen de S(ε).

J. A. Pucheta (www.labimac.blogspot.com)

10

Control óptimo y procesos estocásticos

3.2.1. Teorema de Lyapunov Sea la función escalar definida positiva V(x) una función de energía (ficticia) que depende de x y de t. Si la derivada temporal de V(x) es definida negativa entonces el punto de equilibrio xe en el origen es uniformemente asintóticamente estable y la función V(x) se denominará función Lyapunov. 3.2.2. Teorema de Krasovskii Sea el sistema modelado mediante las ecuaciones

x& = f (x , t ), x (0) = x 0 ,

n

donde x∈R . Para éste sistema, se define la matriz Jacobiano  ∂f1  ∂x ∂f1...∂f n  1 F(x ) = = M ∂x 1...x n  ∂f n  ∂x1 

Además, se define la matriz

∂f1  ∂x n   O M . ∂f n  L ∂x n   L

Fˆ(x ) = F T (x ) + F(x ),

se establece una función Lyapunov para éste sistema haciendo V (x ) = f T (x ) ⋅ f (x ),

debido a que

& (x ) = f T (x ) ⋅ Fˆ ⋅ f (x ). V

3.2.3. Aplicación en sistemas lineales Sea el sistema modelado mediante las ecuaciones lineales homogéneas x& = Ax, x (0) = x 0 .

(3-4)

Se elige la función candidata de Lyapunov V(x ) = x T ⋅ P ⋅ x ,

(3-5)

siendo P una matriz simétrica y definida positiva, la derivada temporal de V(x) es

(

)

& (x ) = x T ⋅ A T P + PA ⋅ x V

(3-6)

donde se requiere que la matriz entre paréntesis sea definida negativa para que la candidata V(x) propuesta sea función de Lyapunov. Por lo tanto debe cumplirse que − Q = A T P + PA

(3-7)

donde Q debe ser definida positiva. Para verificar la existencia de P, se hace el estudio sobre

(A

J. A. Pucheta (www.labimac.blogspot.com)

T

)

P + PA = − I, 11

Control óptimo y procesos estocásticos

igualando a la matriz identidad. Nótese la relación existente entre la función Lyapunov V(x) y su derivada temporal, las expresiones (3-5) y (3-6) muestran que ∂ T x ⋅ P ⋅ x = x T ⋅ A T P + PA ⋅ x = − x T ⋅ Q ⋅ x. ∂t

(

)

(

)

(3-8)

3.3. Problema de control óptimo cuadrático Se usará el segundo método de Lyapunov para resolver el problema del control óptimo formulado. Primero se fijan las condiciones de estabilidad y luego se diseña el controlador dentro de ésas condiciones. El método supone que el sistema es controlable. Reemplazando (3-2) en la primer ecuación del sistema (3-1), se tiene x& t = (A − BK ) ⋅ x t

(3-9)

donde se asume que (A-BK) es estable, es decir, tiene todos los autovalores con parte real negativa. Sustituyendo (3-2) en (3-3), J (x, u ) =

∫ (x

∞

T

)

(

∞

)

Qx + (Kx )T RKx dt = ∫ x T Q + K T RK xdt.

0

(3-10)

0

Usando el Teorema de Lyapunov, para que la candidata V(x ) = x T ⋅ P ⋅ x,

(3-11)

sea Función Lyapunov, entonces su derivada temporal deberá ser definida negativa. Derivando en el tiempo a la V(x), & (x ) = x& T Px + x T Px& V

(3-12)

reemplazando la Ec (3-9) en la derivada temporal de x, se tiene que

(

)

(

)

& (x ) = x T A T P + PA x − x T K T BT P + PBK x. V

(3-13)

El primer término está definido en la Ec (3-7). El segundo término, debe ser definido positivo, y se hará la igualdad K T RK = K T BT P + PBK

(3-14)

teniendo en cuenta la (3-10). Así, la derivada temporal de V(x) deberá cumplir con

(

)

& (x ) = − x T ⋅ Q + K T RK ⋅ x. V

(3-15)

Igualando las expresiones, se tiene que ∂ T x ⋅ P ⋅ x = − x T ⋅ Q + K T RK ⋅ x. ∂t

(

)

(

)

(3-16)

Derivando respecto a t, se tiene T x T Px + x T Px = x T (A − BK ) Px + x T P(A − BK )x

(

)

(

)

x T (A − BK )T P + P(A − BK ) x = − x T ⋅ Q + K T RK ⋅ x ,

(3-17)

que debe resolverse en P simétrica y definida positiva. Como la condición (3-17) debe cumplirse para todo x∈Rn, se resuelve la igualdad a partir de igualar las matrices de ponderación de la forma cuadrática. Por lo tanto, J. A. Pucheta (www.labimac.blogspot.com)

12

Control óptimo y procesos estocásticos

(

)

− (A − BK )T P + P(A − BK ) = Q + K T RK,

(3-18)

corresponde al argumento del funcional de costos a minimizar, y definiendo a la función ξ como ξ = Q + K T RK + (A − BK )T P + P(A − BK ).

(3-19)

Para hallar K, se minimiza la expresión (3-19) respecto de K, teniendo en cuenta las reglas de derivación matricial e igualando a 0 el resultado verificando que la derivada segunda de (3-19) sea positiva para que el extremo sea un mínimo. Dada una matriz X y dos vectores x, y se verifica

(

)

(

)

(

)

∂ xT ⋅ X ⋅ y ∂ xT ⋅ X ⋅ y ∂ xT ⋅ X ⋅ x = X⋅y ; = X T ⋅x ; = 2X ⋅ x , ∂x ∂y ∂x

(3-20)

pero la tercera propiedad sólo es válida si X es simétrica. Derivando la (3-19) respecto a K, se tiene ∂ξ = BT P + BT P + 2RK ∂K

y la derivada segunda de (3-19) es ∂ 2ξ = 2R T ∂K 2

que es definida positiva, lo que indica que el extremo de la derivada primera es un mínimo. Por lo tanto, igualando a cero la derivada primera y despejando K se tiene que K = R −1B T P.

(3-21)

u t = −R −1BT P ⋅ x t .

(3-22)

La ley de control será, entonces,

Reemplazando el valor de K en la igualdad (3-18) se obtiene el valor mínimo de la función implícita. Así, 0 = K T RK + (A − BK )T P + P(A − BK ) + Q = K T RK − (BK )T P + A T P + PA − PBK + Q,

reemplazando, entonces, K por la Ec. (3-21), se tiene

(

)

(

T

) ) P + A P + PA − PBR T

0 = R −1B T P RR −1B T P − BR −1B T P P + A T P + PA − PBR −1B T P + Q,

( ) (R ) B P − (B P ) (BR

0 = BT P

T

−1 T

T

T

T

−1 T

T

−1 T

B P + Q,

que operando, se llega a A T P + PA − PBR −1 B T P + Q = 0,

(3-23)

que es la Ecuación de Riccati reducida. Evaluando a J de la Ec. (3-10) con el ut de la Ec. (3-22) se obtiene ∞

∫ (

)

J (x, u ) = x T Q + K T RK xdt. = − x T Px

∞ 0

(

= − x ∞ Px ∞ − x 0 Px 0 T

T

)

(3-24)

0

donde se ha usado la igualdad (3-8) para resolver la integral. Para determinar el valor en la Ec. (3-24), se considera que los autovalores de (A-BK) tienen parte real negativa, entonces x(t)→0 con t→∞. Por lo tanto la Ec. (3-24) resulta J (x, u ) = x 0 Px 0 . T

J. A. Pucheta (www.labimac.blogspot.com)

(3-25) 13

Control óptimo y procesos estocásticos

Dado el caso en que se diseñe al funcional de costos en términos de la salida y, del sistema de Ec. (3-1), ∞ J (x, u ) =

∫ (y

T

)

Qy + u T Ru dt

(3-26)

0

se reemplaza y por la segunda fila de la Ec. (3-1), quedando J (x, u ) =

∫ (x

∞

)

C QCx + u T Ru dt ,

T T

(3-27)

0

y se emplea CTQC en lugar de Q. Para el diseño del controlador óptimo cuadrático, una vez formulado el problema, se debe resolver la Ecuación de Riccati (3-23) con respecto a P verificando que (A-BK) sea estable. Para calcular el controlador usando Matlab, se dispone del comando K=LQR(A,B,Q,R) o bien [K,P,E]=LQR(A,B,Q,R) donde E contiene a los autovalores de (A-BK).

4. Regulador óptimo lineal en tiempo discreto 4.1. Formulación del problema La formulación del problema de control para el Regulador Óptimo lineal en tiempo discreto es la siguiente. Dado el sistema lineal determinístico x k +1 = A k x k + B k u k  y k = C k x k + D k u k

(4-1)

se desea encontrar una ley de control uk que haga evolucionar al proceso desde x(0)≠0 a x(N)=0 minimizando el siguiente funcional de costo J (x , u ) =

∑ [x Tk Qx k + u Tk Ru k ]+ x TN Sx N

N -1

(4-2)

k =0

con S y Q simétricas y semidefinidas positivas y R simétrica y definida positiva. Para encontrar la ley de control uk, existen diversos métodos, entre los más difundidos están los basados en el principio de optimalidad de Bellman y los que emplean los multiplicadores de Lagrange. Para el caso en que N tienda a infinito en la definición del funcional (4-2), se tiene una formulación del problema conocida como de estado estacionario donde pierde sentido el término xTNSxN ya que al ser estable el sistema controlado siempre será nulo, la cual admite un procedimiento de cómputo basado en la Teoría de Lyapunov. 4.2. Formulación del problema de estado estacionario Se propone formular el problema de control óptimo para emplear la Teoría de Lyapunov, que considera un sistema dinámico en estado estacionario. Dado el modelo dinámico de la Ec. (4-1), se desea encontrar una ley de control uk u k = −K ⋅ x k J. A. Pucheta (www.labimac.blogspot.com)

(4-3) 14

Control óptimo y procesos estocásticos

que haga evolucionar el sistema para k=0 hasta k=∞, minimizando el funcional de costos J (x, u ) =

∞

∑ x Tk Qx k + u Tk Ru k .

(4-4)

k =0

donde Q es simétrica y semidefinida positiva, y R es simétrica y definida positiva. Para resolver éste problema, se empleará el Teorema de estabilidad de Lyapunov. 4.2.1. Estabilidad en tiempo discreto El análisis de estabilidad en el sentido de Lyapunov sirve para sistemas lineales o no lineales de tiempo discreto, variantes o invariantes en el tiempo. Se basa en el segundo método de Lyapunov. Teorema Sea el sistema en tiempo discreto x (k +1)T = f (x kT ),

(4-5)

donde x∈Rn, f(x) ∈Rn, f(0)=0, y T período de muestreo. Se emplea una función que contempla la energía del sistema, y de ésta función se calcula la diferencia temporal, es decir, que dada V(x kT ) = funcion candidata

(4-6)

y la función diferencia para dos intervalos de muestreo se define como

(

)

∆V(x kT ) = V x (k +1)T − V (x kT ).

(4-7)

Si existe una función escalar continua V(x) tal que 1. V(x)>0 ∀x≠0, 2. ∆V(x)<0 ∀x≠0, ∆V(x)=V(f(xkT))-V(xkT), 3. V(0)=0, 4. V(x)→∞ con ||x||→∞, entonces el estado de equilibrio x=0 es asintótica-globalmente estable y V(x) es una función de Lyapunov. Nótese que 2 puede ser reemplazado por ∆V(x)≤0 ∀x, y ∆V(x) no se hace cero para toda secuencia {xkT} solución de (4-5). Suponiendo que en el sistema de la Ec. (4-1) se hace uk=0, se propone la siguiente función candidata de Lyapunov V(x k ) = x Tk Px k

(4-8)

donde P es simétrica y definida positiva. Entonces, se calcula ∆V (x k ) = V (Ax k ) − V (x K ).

J. A. Pucheta (www.labimac.blogspot.com)

15

Control óptimo y procesos estocásticos

Reemplazando la (4-8) se tiene,

∆V (x k ) = (Ax k )T PAx k − x Tk Px k ,

(

)

∆V (x k ) = x Tk A T PA − P x k .

(4-9)

Para asegurar estabilidad asintótica, se impone que la (4-9) sea definida negativa, y se puede escribir que ∆V(x k ) = − x Tk Qx k ,

(4-10)

donde Q es definida positiva, da la condición suficiente para estabilidad asintótica − Q = A T PA − P.

(4-11)

Es conveniente especificar Q simétrica y definida positiva, y luego verificar que P - determinada por la (4-11)- es definida positiva o no. Si P es definida positiva, entonces la V(x) propuesta por (4-8) es función de Lyapunov y se demuestra estabilidad. Por otro lado, nótese que de la Ec. (4-8), se calcula su diferencia temporal como ∆V (x k ) = V (x k +1 ) − V (x k )

de donde resulta que ∆V (x k ) = x Tk +1Px k +1 − x Tk Px k

(4-12)

y a su vez se iguala al lado derecho de la Ec. (4-10), y da − x Tk Qx k = x Tk +1 Px k +1 − x Tk Px k .

(4-13)

4.2.2. Problema de control óptimo discreto estacionario Reemplazando la ley de control (4-3) en la expresión del funcional de costos (4-4), se puede escribir que J (x, u ) =

∑ x Tk (Q + K T RK )x k . ∞

(4-14)

k =0

Ahora se busca la solución al problema de control en tiempo discreto en estado estacionario, que se basa en la Ec. (4-13), pero si la acción de control no es nula, aparece la modificación de incorporar al controlador. Así, para el caso en que uk≠0 en la entrada al sistema (4-1), se reemplaza en las ecuaciones del sistema la expresión de la ley de control (4-3), y se tiene que la Ec. (4-13) se transforma en

(

)

− x Tk Q + K T RK x k = x Tk +1Px k +1 − x Tk Px k .

(4-15)

si se hace coincidir a la matriz Q Ec. (4-13) con el argumento de ponderación de la forma cuadrática de la Ec. (4-14). De aquí que minimizando la (4-14) con respecto a K, se encuentra la ley de control óptima (4-3). Para ello, se opera en la igualdad (4-15) que puede escribirse como

(

)

(

)

− x Tk Q + K T RK x k = x Tk (A − BK )T P(A − BK ) − P x k .

(4-16)

Ésta última igualdad debe cumplirse para todo valor de xk, por lo tanto, se tiene que minimizando Q + K T RK = −(A − BK )T P(A − BK ) + P

(4-17)

con respecto a K es lo mismo que hacerlo en la (4-14). Para ello, se define una función ζ a partir de J. A. Pucheta (www.labimac.blogspot.com)

16

Control óptimo y procesos estocásticos

la Ec. (4-17) como, ς = Q + K T RK + (A − BK )T P(A − BK ) − P.

(4-18)

Se procede minimizando a la Ec. (4-18) respecto a K, derivando miembro a miembro, considerando las reglas (3-20) y que la regla de la cadena en operaciones con matrices es

(

∂ F(Tx ) ⋅ X ⋅ F(x ) ∂x

) = ∂F( ) ⋅ ∂F( ) ⋅ X ⋅ F( ) . T x

x

∂x

(4-19)

x

∂X

Para minimizar a la expresión de la Ec. (4-18), se deriva respecto a K, se iguala a cero y se despeja K en el extremo. Para demostrar que es un mínimo, se hace la derivada segunda de la Ec. (4-18) respecto a K, que deberá ser positiva. Entonces, derivando respecto a K a la Ec. (4-18), se tiene ∂ς = 2RK + (− B)T ⋅ 2P(A − BK ) . ∂K

La derivada segunda será

{

}

(

)

T ∂ 2 ς ∂ 2RK − 2 ⋅ BT PA + 2 ⋅ BT PBK = = 2R T + 2 BT PB 2 ∂K ∂K

cuyo resultado es una matriz definida positiva. Por lo tanto, el extremo en cuestión es un mínimo. Operando para despejar K de la derivada primera, RK = B T PA − B T PBK ,

es decir RK + B T PBK = B T PA

y de aquí

(

K = R + B T PB

)

−1

B T PA.

(4-20)

Reemplazando el valor de K óptimo encontrado en la Ec. (4-17), se llega a la Ecuación de Riccati de estado estacionario. Para ello, se reemplaza entonces la (4-20) en la (4-18). Entonces, operando primero a la Ec. (4-18), Q + A T PA + K T RK − A T PBK − (BK )T PA + (BK )T PBK = P,

(4-21)

donde el lado derecho y los primeros dos términos del lado izquierdo no dependen de K, pero sí los términos desde el tercero hasta el sexto. Reemplazando K de la Ec. (4-20) en el tercer término, se tiene

(

)

(

T

)

−1 −1 K T RK =  R + B T PB B T PA  R R + B T PB B T PA  

(

) (

)

T

(

)

T −1 −1 = B T PA  R + B T PB  R R + B T PB B T PA.  

Haciendo lo mismo con el cuarto término,

(

A T PBK = B T PA

) (R + B PB) T

T

−1

B T PA .

Con el quinto,

(BK )T PA =  B(R + B T PB)

−1



J. A. Pucheta (www.labimac.blogspot.com)

T

(

) (

)

T

−1 T B T PA  PA = B T PA  B R + B T PB  PA   

17

Control óptimo y procesos estocásticos

(

) (

(

) (

)

T

−1 T = B T PA  B R + B T PB  PA  

)

T

T −1 = B T PA  R + B T PB  B T PA  

y con el sexto término, se tiene

(BK )T PBK =  B(R + B T PB)

−1



(

= (PA )T  B R + B T PB 

(

)

−1

T

(

)

−1 B T PA  PB R + B T PB B T PA  T

(

)

T

(

)

(

)

B T  PB R + B T PB 

) (

)

−1

B T PA

T −1 −1 = B T PA  B R + B T PB  PB R + B T PB B T PA  

(

) (

)

T

T −1 −1 = B T PA  R + B T PB  B T PB R + B T PB B T PA.  

Nótese que todos los términos se expresaron como una forma cuadrática de BTPA, por lo tanto, se puede hacer la suma de los términos y agruparlos reemplazándolos en la Ec. (4-21), entonces operando se tiene,

(

) (

)

(

)

(

) (

)

(

)

T −1 T −1 −1 −1 T P = Q + A T PA + B T PA   R + B T PB  R R + B T PB − R + B T PB −  B R + B T PB  +      −1 T −1  +  R + B T PB  B T PB R + B T PB B T PA     ,

(

)

que operando con el primer y cuarto término dentro del paréntesis, para agrupar como forma cuadrática, se tiene

(

) (

)(

T

) (

)

(

)

T

 R + B T PB −1  B T PB + R R + B T PB −1 − R + B T PB −1 −  R + B T PB −1          ,

y ahora se puede simplificar ya que las inversas existen

(

)

T

(

)

(

)

(

T

)

 R + B T PB −1  − R + B T PB −1 −  R + B T PB −1  = − R + B T PB −1         ,

que será el término medio de la forma cuadrática en BTPA, finalmente, la Ec. (4-21) queda

(

P = Q + A T PA − B T PA

) (R + B PB) T

−1

T

B T PA

(4-22)

que es la Ecuación de Riccati de estado estacionario en tiempo discreto. A su vez, operando se llega a

(

)

A

(4-23)

P = Q + A T P −1 + B T R −1B T

)

A.

(4-24)

P = Q + A T P I + B T R −1B T P

−1

y también a,

(

−1

Para el diseño del controlador, debe resolverse la Ecuación de Riccati en P. El funcional de costo se puede evaluar, usando las igualdades de la Ec. (4-14) y la Ec. (4-15), se llega a J (x, u ) =

∞

∑

k =0 J. A. Pucheta (www.labimac.blogspot.com)

x Tk Qx k + u Tk Ru k =

∞

∑ x Tk +1Px k +1 − x Tk Px k = x T0 Px 0 . k=0

18

Control óptimo y procesos estocásticos

(4-25) 4.3. Problema de control óptimo lineal de continuo a discreto Se estudiará el problema de control óptimo que se plantea en el tiempo continuo pero se implementa en tiempo discreto. Formulación del problema Dado el sistema modelado mediante las Ec. (3-1), se desea encontrar la ley de control ut de la forma u t = u kT , t ∈ [kT , (k + 1)T ]

(4-26)

que minimice el funcional de costo definido por J (x, u ) =

(

)

1 T 1 tf x t f Sx t f + ∫ x Tt Qx t + u Tt Ru t dt 2 20

(4-27)

con Q simétrica y semidefinida positiva y R simétrica y definida positiva. Se supone que el sistema representado por las Ec. (3-1) puede ser reemplazado por x (k +1)T = A T ⋅ x kT + B T ⋅ u kT  y kT = C ⋅ x kT

(4-28)

La solución general del sistema de Ec. (3-1) ∀t∈[t0,t] es t

(4-29)

x t = e A (t − t 0 ) x t 0 + ∫ e A (t −s ) Bu s ds. t0

Si se expresa la solución en el intervalo ∀t∈[kT,(k+1)T] y considerando que us en ése intervalo es constante según se expresa en la Ec. (4-26), se tiene que  t  x t = e A (t − kT ) x kT +  ∫ e A (t −s )Bds u kT .  kT 

(4-30)

Así, la solución xt ∀t∈[kT,(k+1)T] puede expresarse como x t = F1 ⋅ x kT + F2 ⋅ u kT ,

(4-31)

F1= e A (t −kT )  t  A ( t −s ) Bds. F2= ∫ e  kT

(4-32)

donde

Reemplazando la Ec. (4-31) y la Ec. (4-26) en la expresión del funcional de costo J de la Ec. (4-27) J (x, u ) =

1 T 1 N −1(k +1)T x NT Sx NT + ∑ ∫ (F1 ⋅ x kT + F2 ⋅ u kT )T Q(F1 ⋅ x kT + F2 ⋅ u kT ) + u TkT Ru kT dt (4-33) 2 2 k =0 kT

{

}

Operando el término entre llaves, considerando que todas las cantidades son escalares, y que u y x J. A. Pucheta (www.labimac.blogspot.com)

19

Control óptimo y procesos estocásticos

son constantes en el intervalo de integración, se llega a la expresión de J de la forma J (x, u ) =

1 T 1 N −1 x NT Sx NT + ∑ x Tk Q1x k + 2 x Tk M1u k + u Tk R 1u k 2 2 k =0

(4-34)

donde (k +1)T

Q1=

∫

F1T QF1dt M1=

kT

(k +1)T

∫

kT

F1T QF2 dt R 1=

(k +1)T

∫

(F QF T 2

2

)

+ R dt .

kT

Nótese que el funcional en tiempo continuo de la Ec. (4-27) no queda expresado como J (x, u ) =

1 T 1 N −1 x N Sx N + ∑ x Tk Qx k + u Tk Ru k , 2 2 k =0

(4-35)

sino que se agrega un término cruzado que involucra a xkT y ukT y se modifican Q y R.

J. A. Pucheta (www.labimac.blogspot.com)

20

Control óptimo y procesos estocásticos

Ejemplo de aplicación Para la planta de la Fig. 6, suponiendo que el flujo entre los tanques es proporcional a la diferencia de nivel de líquido en los mismos, y aplicando el principio de conservación de la masa (o del volumen si el fluido es incompresible), las ecuaciones que gobiernan al sistema son dh1 1 = (h 2 - h1) + u t dt R1 dh 2 1 1 =(h 2 - h1) A2 h2 dt R1 R2 A1

donde ut es el caudal de líquido que entra al tanque 1.

Fig. 6. Planta hidráulica.

Se considera que A1 = A 2 = 1 ; R1 =

1 1 ; R 2 = y el sistema de medidas es el MKS. Se pretende 2 3

controlar la altura del líquido en el tanque 2 controlando el caudal de líquido entrante al tanque 1. Se pide, para los dominios en tiempo continuo y discreto: 1. Deducir las ecuaciones de estado del sistema. 2. Calcular el controlador óptimo K siguiendo el procedimiento visto, donde deben resolverse las Ecuaciones de Riccati en P. 3. Graficar diferentes respuestas del sistema para distintas matrices Q y R. Verificar los resultados con las igualdades Jmin=xT(0)Px(0). 4. Proponer una estrategia de control para referencia distinta de cero, por ejemplo 1, 10, y 100.

Solución 1. Utilizando las ecuaciones que gobiernan al sistema A1

J. A. Pucheta (www.labimac.blogspot.com)

dh1 1 (h 2 − h1 ) + u ( t ) = dt R1

21

Control óptimo y procesos estocásticos

A2

dh 2 1 1 = − (h 2 − h1 ) − h2 dt R1 R2

Valuando con A1=A2=1, R1=0,5 y R2=0,333, se tiene: dh1 = 2(h 2 − h1 ) + u ( t ) dt dh 2 = −2(h 2 − h1 ) − 3h 2 dt

O lo que es lo mismo, para tiempo continuo se tiene: h& 1 = −2h1 + 2h 2 + u ( t ) & h 2 = 2h1 − 5h 2

Expresando las ecuaciones con matrices  x& 1  − 2 2   x1  1  x&  =   ⋅   +   ⋅ u.  2   2 − 5  x 2   0  x  y = [0 1] ⋅  1  x 2 

(4-36)

2. Para hallar la matriz P a partir de la Ec. de Riccati reducida, que es el caso del dominio del tiempo continuo, se opta por la iteración siguiente

(

)

Pg = - Q - A T P + PBR -1BT P A -1   T T P = P − γ P − Pg donde γ=0,01 si se hacen 500 iteraciones.

(

)

(4-37)

Desde el programa Prog_TC_01.m pueden ensayarse diferentes casos para hallar la matriz P. Nótese que el método no sirve así planteado si A es singular. Para hallar la matriz P a partir de la Ecuación de Riccati en estado estacionario, se propone hacer la siguiente operación recursiva

(

)

-1

Pk +1 = Q + A T Pk A - (BT Pk A)T R + BT Pk B BT Pk A

(4-38)

que en 10 iteraciones llega al P invariante. Desde el programa Prog_TD_02.m pueden ensayarse diferentes casos para hallar la matriz P. 3. En el caso de tiempo continuo, la evaluación del funcional de costo se implementa mediante la evaluación de la Ec. (3-3), donde el estado se calcula mediante el método de integración de Euler. Debe coincidir con el valor de la Ec. (3-25) para referencia nula. Ensayar al programa Prog_TC_01.m.

J. A. Pucheta (www.labimac.blogspot.com)

22

Control óptimo y procesos estocásticos

En el caso del tiempo discreto, la evaluación del funcional de costos (4-4) se implementa directamente, y debe coincidir con el valor según la Ec. (4-25) para referencia nula. Ensayar el programa Prog_TD_02.m. 4. Cuando la referencia es distinta de cero, puede usarse el mismo controlador que se diseñó, pero hay que modificar a la acción de control. Una propuesta, se fundamenta en una ley de control de la forma u t = − Kx t + Grt

(4-39)

donde rt es la referencia y tiene la misma dimensión que yt, G es una ganancia de prealimentación de la referencia.

x0 rt

ut G

x&

B

∫ ⋅ dt

xt

yt

C

A

-K

Fig. 7. Esquema de control óptimo en tiempo continuo con referencia distinta de cero.

Como la referencia se alcanza en el estado estacionario del sistema, se hace el análisis para diseñar a G usando Laplace. De la Ec. (3-1) se tiene sx s = A ⋅ x s + B ⋅ u s  y s = C ⋅ x s

(4-40)

Reemplazando la transformada de Laplace de la Ec. (4-39) en la primera de la Ec. (4-40), se tiene sx s = A ⋅ x s + B(− Kx s + Grs ) = (A − BK )x s + BGrs (sI − (A − BK ))x s = BGrs

x s = (sI − (A − BK ))−1 BGrs

y la salida y se obtiene con la segunda de la Ec. (4-40)

y s = C(sI − (A − BK ))−1 BGrs

de donde se puede definir una ganancia en s, como H s = C(sI − (A − BK ))−1 BG

(4-41)

y así J. A. Pucheta (www.labimac.blogspot.com)

23

Control óptimo y procesos estocásticos

y s = H s ⋅ rs .

(4-42)

La referencia rs será un escalón unitario, ya que se trata de un problema de regulación, y se emplea el teorema del valor final, que establece que lim s ⋅ {y s } = lim y t .

s→0

t →∞

Aplicando la igualdad a la Ec. (4-42) con rs escalón unitario, se tiene que la salida será I = C(− (A − BK ))−1 BG

y despejando G de aquí, se tiene que

[

G = − C(A − BK )−1 B

]

−1

.

(4-43)

En el programa Prog_TC_01.m se tiene una implementación de ésta propuesta. Para el caso de tiempo discreto, se procede de manera análoga, sólo que cambia el dominio. La acción de control propuesta será u k = − Kx k + Grk

(4-44)

donde nuevamente la ganancia de prealimentación de la referencia G es la incógnita. Reemplazando ésta uk en la Ec. (4-1) en el dominio de Z, se tiene zx = Ax z + B(− Kx z + Grk )  y z = Cx z

despejando xz de la primera, se tiene y por lo tanto yz es

(4-45)

x z = (zI − A + BK )−1 BGrz

y z = C(zI − A + BK )−1 BGrz

(4-46)

en donde se define que H z = C(zI − A + BK )−1 BG.

(4-47)

Para hallar G, se supone que yz alcanza la referencia rz en el estado estacionario, y por lo tanto se usa el teorema del valor final (TVF), que establece que z −1 ⋅ {y z } = lim y k . z →1 z k →∞ lim

Aplicando el TVF a yz cuando rz es un escalón unitario, se tiene que I = C(I − A + BK )−1 BG

y por lo tanto, despejando G se tiene que

[

]

−1

G = C(I − A + BK )−1 B .

(4-48)

En el programa Prog_TD_02.m se tiene una implementación de ésta propuesta. Nótese que en los dos casos vistos las ganancias de prealimentación dejan al sistema a lazo abierto con respecto a la referencia. Para ello se proponen sendas estrategias que incorporan un integrador J. A. Pucheta (www.labimac.blogspot.com)

24

Control óptimo y procesos estocásticos

del error de control.

Incorporación de un integrador Incorporando un término de integración en el diseño en tiempo continuo, se tiene que la ley de control ut es u t = − Kx t + K1ξ t

(4-49)

ξ& t = rt − y t = rt − Cx t

(4-50)

donde

definiéndose el nuevo estado ξ, como la salida de un integrador cuando a la entrada está presente el error de control rt-yt. El esquema de control se muestra en la Fig. 8. x0 ξ& t

rt

∫ ⋅ dt

ξt

x&

ut K1

∫ ⋅ dt

B

yt

xt C

-

A

-K

Fig. 8. Esquema de control óptimo en tiempo continuo con un integrador en el lazo para referencia distinta de cero.

El sistema de orden incrementado puede escribirse como x& t   A 0  x t  B 0  ⋅   +   u t +   rt . ξ&  =   1   t   − C 0  ξ t   0 

(4-51)

Para el estado estacionario, se tiene con t→∞, x& ∞   A 0  x ∞  B 0  ξ&  = − C 0 ⋅  ξ  +  0  u ∞ + 1  r ∞ .     ∞    ∞ 

(4-52)

Suponiendo una referencia constante, escalón, se pueden restar las ecuaciones (4-51) y (4-52) obteniendo  A 0  B e& =  ⋅ e +  u e   − C 0 0  J. A. Pucheta (www.labimac.blogspot.com)

(4-53) 25

Control óptimo y procesos estocásticos

que determinan la dinámica del error, donde x − x ∞  e= t , ξ t − ξ∞ 

u e = − Kx e + K 1ξ

(4-54)

Haciendo el diseño del controlador con las matrices  A 0  B Aa =  , B a =  0  − C 0    

(4-55)

se llevará el error de de control de la Ec. (4-54) a cero. El algoritmo está implementado en el programa Prog_TC_Integrador.m para Matlab.

Para el caso del tiempo discreto, si se incorpora un integrador en el lazo del controlador, se tendría un esquema como el de la Fig. 9.

x0 rk

v

xk+1

uk K1

Iq-1

B

xk

yk C

vk-

q-1 A

-K

Fig. 9. Esquema de control óptimo en tiempo discreto con un integrador en el lazo para referencia distinta de cero.

La ley de control uk es u k = − Kx k + K1v k

(4-56)

donde la variable vk se define como la salida de un integrador del error de control rk-yk, v k = v k −1 + rk − y k

(4-57)

Operando sobre la señal vk, se tiene

v k +1 = v k + rk +1 − y k +1 = v k + rk +1 − C(Ax k + Bu k ) v k +1 = −CAx k + v k − CBu k + rk +1

de lo cual se tiene que 0  x k   B   x k +1   A 0   v  = − CA 1  ⋅  v  + − CB u k + 1  rk +1.   k      k +1   J. A. Pucheta (www.labimac.blogspot.com)

26

Control óptimo y procesos estocásticos

(4-58) Asignando a rk un escalón, y haciendo k→∞, se tiene

0  x ∞   B  x ∞   A 0  v  = − CA 1  ⋅  v  + − CB u ∞ + 1  r∞ .   ∞      ∞ 

(4-59)

restando a la Ec. (4-58) la Ec. (4-59), se llega a 0  A  B  e k +1 =  ⋅ ek +    u ek . − CA 1 − CB

(4-60)

x − x ∞  ek =  k , u ek = −[K K 1 ]e k . v k − v∞ 

(4-61)

donde

La expresión (4-60) determina la dinámica del error de control. Por lo tanto, diseñando el controlador considerando las matrices de orden incrementado 0  A Aa =  , y  − CA 1 

 B  Ba =   − CB

(4-62)

se obtendrá el controlador óptimo con un integrador del error de control. En el programa para Matlab Prog_TD_Integrador.m está implementado el algoritmo propuesto.

J. A. Pucheta (www.labimac.blogspot.com)

27

Control óptimo y procesos estocásticos

5. Regulador óptimo lineal en el transitorio Se pretende deducir la expresión del controlador óptimo en el dominio del tiempo para que haga evolucionar a un proceso modelado como un sistema lineal, desde cualquier estado inicial en el instante k=0 hasta el valor 0 en el instante k=N. 5.1. Formulación del problema en el transitorio El procedimiento que va a implementarse se conoce como minimización de una función sujeto a restricciones, las cuales están impuestas por el modelo dinámico del proceso. El concepto del método es útil para considerar restricciones arbitrarias en la formulación del problema de control óptimo. Las restricciones se agregan al funcional de costo mediante los multiplicadores de Lagrange (ML) λ, que son vectores con la misma dimensión que el vector de estado. Se minimiza J de la Ec (4-2), cuando está sujeta a las restricciones especificadas por la Ec (4-1), para una condición inicial x(0)≠0. Al emplear un conjunto de ML {λ(1), λ(2),.... λ(N)} que forman un vector adjunto o covector, se re define al funcional J(x,u) de la Ec (4-2) como el funcional de costo aumentado Ja(x,u), dado por J a (x , u ) =

∑ {x Tk Qx k + u Tk Ru k + λTk +1[Ax k + Bu k − x k +1 ] + [Ax k + Bu k − x k +1 ]

N -1

T

λ k +1

k =0

}

(5-1)

+ x TNSx N

donde los ML se escriben así para mantener la dimensión escalar de L. Minimizar (5-1) es equivalente que minimizar (4-2) sujeto a (4-1). Para minimizar Ja(x,u), se busca diferenciar a Ja(x,u) respecto de cada componente dadas como x, u y λ e igualar a 0 las expresiones. Suele ser conveniente trabajar con los conjugados de cada componente, pero aquí se trabajará con las componentes originales ya que sólo se hará tratamiento simbólico. Por lo tanto, se pretende hallar las siguientes igualdades ∂J a (x , u ) ∂J (x , u ) ∂J (x , u ) ∂J (x , u ) = 0, a = 0, a = 0, a = 0, ∂x k ∂x N ∂u k ∂λ k

(5-2)

con k de 1 a N para xk y λk, y de 1 a N-1 para u. Para facilitar el procedimiento, se explicitará al funcional de costo Ja(x,u) de la Ec. (5-1) J a (x , u ) =

∑ {⋅ ⋅ ⋅ + x Tk −1Qx k -1 + u Tk -1Ru k -1 + λTk [Ax k −1 + Bu k −1 − x k ] + [Ax k −1 + Bu k −1 − x k ]

N -1

k =0

T

λk

+ x Tk Qx k + u Tk Ru k + λTk +1[Ax k + Bu k − x k +1 ] + [Ax k + Bu k − x k +1 ]T λ k +1

+ x TN −1Qx N -1 + u TN -1Ru N -1 + λTN [Ax N −1 + Bu N −1 − x N ] + [Ax N −1 + Bu N −1 − x N ]T λ N + x TNSx N .

J. A. Pucheta (www.labimac.blogspot.com)

28

}

(5-3)

Control óptimo y procesos estocásticos

Ahora, se visualiza claramente dónde están las funciones de las variables correspondientes para ejecutar el procedimiento de (5-3). Recordar las igualdades de derivación matricial expresadas en la Ec. (3-20). Se obtiene,

∂J a (x, u ) = 2Qx k + A T λ k +1 + A T λ k +1 − λ k − λ k = 2 Qx k + A T λ k − λ k ∂x k

{

}

de donde se deduce ∂J a (x , u ) = 0 ⇒ Qx k + A T λ k +1 − λ k = 0, k = 1,2,...., N - 1. ∂x k

(5-4)

∂J a (x , u ) = 2Sx N − λ N − λ N = 2{Sx N − λ N } ∂x N

de donde se obtiene ∂J a (x , u ) = 0 ⇒ Sx N − λ N = 0. ∂x N

(5-5)

∂J a (x , u ) = 2Ru k + Bλ k +1 + Bλ k +1 = 2{Ru k + Bλ k +1} ∂u k

por lo tanto ∂J a (x , u ) = 0 ⇒ Ru k + BT λ k +1 = 0, k = 1,2,...., N - 1. ∂u k

(5-6)

∂J a (x , u ) = Ax k + Bu k −1 − x k + Ax k + Bu k −1 − x k = 2{Ax k + Bu k −1 − x k } ∂λ k

por último

∂J a (x , u ) = 0 ⇒ Ax k −1 + Bu k −1 − x k = 0, k = 1,2,..., N. ∂λ k

(5-7)

Nótese que la Ec. (5-7) es la ecuación de estados del sistema (4-1). Para verificar que el extremo encontrado sea un mínimo, puede hacerse la derivada segunda de (5-3) respecto a sus variables y se verá que es positiva en los casos de xk, xN, y uk, pero en λk el método no decide. Para hallar la expresión del controlador, hay que operar entre las ecuaciones (5-4), (5-5), (5-6) y (5-7) para independizarse de λ. Entonces, despejando λ de (5-4), λ k = Qx k + A T λ k +1 , k = 1,2,...., N - 1.

(5-8)

tiene como condición final a la Ec. (5-5), es decir λ N = Sx N .

(5-9)

u k = − R −1B T λ k +1 , k = 1,2,...., N - 1.

(5-10)

De la expresión (5-6) se despeja uk,

Sustituyendo la expresión de uk de (5-10) en la ecuación de estado del sistema (4-1), se tiene para J. A. Pucheta (www.labimac.blogspot.com)

29

Control óptimo y procesos estocásticos

A y B invariantes x k +1 = Ax k − BR −1BT λ k +1 , x (0 ) = x 0 .

(5-11)

Para resolver el problema de optimización, deben resolverse simultáneamente las expresiones de la Ec. (5-8) y la Ec. (5-11). Las condiciones de borde serán λN y x(0). De la Ec. (5-9), se define λ k = Pk x k , k = 0,1,...N - 1

(5-12)

donde la matriz Pk es simétrica y real, y PN=S. Sustituyendo (5-12) en las expresiones (5-8) y (5-11), se tiene Pk x k = Qx k + A T Pk +1x k +1 , k = 1,2,...., N - 1.

y análogamente para (5-11)

(5-13)

x k +1 = Ax k − BR −1B T Pk +1x k +1 , x (0 ) = x 0 .

(5-14)

Las expresiones (5-13) y (5-14) se denominan transformaciones de Riccati, por no contener a λ. De la Ec. (5-14) se tiene que

(I + BR

−1

)

B T Pk +1 x k +1 = Ax k , x (0 ) = x 0 .

(5-15)

Se demuestra que el paréntesis tiene inversa si P es al menos semidefinida positiva. Por lo tanto

(

x k +1 = I + BR −1 B T Pk +1

)

−1

Ax k , x (0 ) = x 0 .

(5-16)

Sustituyendo la Ec. (5-16) en la Ec. (5-13)

(

Pk x k = Q + A T Pk +1 I + BR −1B T Pk +1 

)

−1

A  x k , k = 1,2,...., N - 1. 

(5-17)

Multiplicando ambos miembros de la Ec. (5-17) por xkT se tiene que

(

x Tk Pk x k = x Tk Q + A T Pk +1 I + BR −1B T Pk +1 

)

−1

A  x k , k = 1,2,...., N - 1. 

(5-18)

Como la igualdad debe cumplirse para todo x, entonces

(

Pk = Q + A T Pk +1 I + BR −1B T Pk +1

)

−1

A,

k = 1,2,...., N - 1.

(5-19)

que es la Ecuación de Riccati. Sabiendo que PN=S, se resuelve PN-1 y así sucesivamente hasta P0. Despejando λk+1 de a Ec. (5-8) se tiene que

( )

λ k +1 = A T

−1

(λ k − Qx k ),

(5-20)

k = 1,2,...., N - 1.

El vector de control óptimo se escribe entonces reemplazando en (5-10) la Ec.(5-20), se obtiene

( )

u k = − R −1B T A T J. A. Pucheta (www.labimac.blogspot.com)

−1

(λ k − Qx k ),

k = 1,2,...., N - 1.

(5-21) 30

Control óptimo y procesos estocásticos

Reemplazando λk de la Ec. (5-12) en la (5-21), se tiene

( )

u k = − R −1B T A T

−1

(Pk − Q )x k

= − K k x k , k = 1,2,...., N - 1.

(5-22)

de donde se define el controlador Kk como

( )

K k = R −1B T A T

−1

(Pk − Q )

k = 1,2,...., N - 1.

(5-23)

Evaluando el funcional de costo de la Ec. (4-2) con la uk de la (5-22) se encontrará el Jmin. Para simplificar el procedimiento, multiplicando ambos miembros de (5-13) por xkT, se tiene x Tk Pk x k = x Tk Qx k + x Tk A T Pk +1x k +1 , k = 1,2,...., N - 1.

(5-24)

y reemplazando el término xkTAT por su equivalente de la Ec. (5-15)

(

)

(5-25)

)

(5-26)

T

x Tk Pk x k = x Tk Qx k + x Tk +1 I + BR −1B T Pk +1 Pk +1 x k +1 , k = 1,2,...., N - 1.

y despejando el término en xkTQxk se tiene

(

T

x Tk Qx k = x Tk Pk x k − x Tk +1 I + BR −1B T Pk +1 Pk +1 x k +1 , k = 1,2,...., N - 1.

En la Ec. (5-10) se reemplaza λk+1 partiendo de la Ec. (5-12) que está valuada en k, u k = − R −1B T Pk +1 x k +1 , k = 1,2,...., N - 1.

Así,

ukTRuk

(5-27)

resulta

(

) (

)

T

u Tk Ru k = R −1B T Pk +1 x k +1 R R −1B T Pk +1 x k +1 , k = 1,2,...., N - 1.

(5-28)

Sumando las expresiones de la Ec. (5-26) con la Ec. (5-28) se tiene que

(

)

T

x Tk Qx k + u Tk Ru k = x Tk Pk x k − x Tk +1 I + BR −1B T Pk +1 Pk +1 x k +1 +

(R

−1

) ( T

)

B T Pk +1 x k +1 R R −1B T Pk +1 x k +1 ,

(5-29)

k = 1,2,...., N - 1.

y operando se llega a x Tk Qx k + u Tk Ru k = x Tk Pk x k − x Tk +1Pk +1 x k +1 ,

k = 1,2,...., N - 1.

(5-30)

Sustituyendo la Ec. (5-30) en la Ec. (4-2)

[(

) (

)

]

J min = x T0 P0 x 0 − x1T P1 x1 + x1T P1 x 1 − x T2 P2 x 2 + L + x TN −1PN −1 x N −1 − x TN PN x N + x TN Sx N (5-31)

pero como PN=S, se tiene que el Jmin resulta J min = x T0 P0 x 0 . J. A. Pucheta (www.labimac.blogspot.com)

31

Control óptimo y procesos estocásticos

(5-32) En el Ejemplo de aplicación visto, se pretende implementar un control con 3, 4 y 5 etapas con el tiempo de muestreo de 0,104 seg. Agregar como ítem N° 5, con la salvedad de que se implementará en tiempo discreto únicamente. Como ítem N°6, obtener simulaciones con la referencia distinta de cero. Tiene sentido incorporar un término integrador?. En el programa para Matlab Prog_TD_N_03.m está codificado el algoritmo.

%Uso del control óptimo en la planta hidráulica %Control óptimo de tiempo discreto en el transitorio. %Autor JAP. %29 de noviembre de 2007. clear; Ts=.1*2*pi/6;kmax=5; Ac=[-2 2;2 -5]; Bc=[1;0];Cc=[0 1];Dc=0; SYSc=ss(Ac,Bc,Cc,Dc); [SYS,G] = c2d(SYSc,Ts,'zoh'); %Acá SYS está discretizado con el To=Ts. [A,B,C,D] = ssdata(SYS); Q=[1 0;0 10];R=.01;S=[1 0;0 100]; %Cálculo del LQR P=zeros(2); K=zeros(kmax,2); P=S; %condición inicial de P for hi=kmax:-1:1 P= Q + A'*P*inv(eye(2)+B*inv(R)*B'*P)*A; K(hi,:)=inv(R)*B'*inv(A')*(P-Q); end %Con referencia distinta de 0. r=0*pi;x=[1;1]; Jmin=x'*P*x; Jxu(1)=0;Jxx(1)=0; for k=1:kmax x1(k)=x(1); x2(k)=x(2); y(k)=C*x; Gj=inv(C*inv(eye(2)-A+B*K(k,:))*B); uq=Gj*r; u=-K(k,:)*x+uq; Jxu(k+1)=Jxu(k)+x'*Q*x+u'*R*u; Jxx(k+1)=Jxx(k)+x'*(Q+K'*R*K)*x; x=A*x+B*u; end Jxu(k+1)=Jxu(k+1)+x'*S*x; Jinf=Jxu(k+1); Jmin/Jinf t=0:Ts:Ts*(kmax-1); figure; subplot(2,2,1),plot(t,x1);title(['h_1(t)']);

J. A. Pucheta (www.labimac.blogspot.com)

32

Control óptimo y procesos estocásticos

6. Control óptimo basado en programación dinámica Se ha visto el procedimiento de resolver el problema de control óptimo para sistemas en tiempo contínuo y tiempo discreto, basados en el Teorema de Lyapunov y usando multiplicadores de Lagrange. Ahora se empleará un principio de optimización que sirve para resolver el mismo problema de control pero que además permite utilizar cualquier funcional de costo y contemplar las restricciones en las variables de estado y en las de control de manera natural. 6.1. Principio de optimalidad de Bellman

c • b• • d

a•

Fig. 10. Trayectoria óptima desde a hasta d en línea continua. En línea de trazos, otra posible trayectoria óptima.

Sea la trayectoria óptima mostrada en la Fig. 10. Suponiendo que la primer decisión, hecha en a, resulta en el tramo a-b con costo ga-b y las decisiones siguientes tienen un costo Jb, correspondiente al segmento b-d es decir, desde el punto b hasta el final. El costo mínimo desde a hasta d es J *a = g ab + J b

Afirmación: Si el tramo a-b-d es una trayectoria óptima de a hasta d, entonces b-d es la trayectoria óptima desde b hasta d. Prueba: Suponiendo que el tramo b-c-d sea la trayectoria óptima, como muestra en línea de trazos la Fig. 10, el costo desde b hasta d será entonces J 1b = g bc + J c g ab + J1b < g ab + J b = J *a

pero puede cumplirse únicamente violando la condición que a-b-d es la trayectoria óptima desde a hasta d. Así se prueba la afirmación. El principio de optimalidad puede enunciarse: Una secuencia óptima de control (política óptima) tiene la propiedad de que cualquiera sea el par (estado, acción) inicial las decisiones restantes deben constituir una secuencia óptima de control (política óptima) con respecto al estado resultante de la primera acción de control.

J. A. Pucheta (www.labimac.blogspot.com)

33

Control óptimo y procesos estocásticos

6.1.1. Formulación del problema Formulación del problema de control para el Regulador Óptimo lineal determinístico: dado el sistema lineal determinístico descrito por las Ec. (4-1) se desea encontrar una ley de control uk que haga evolucionar al proceso desde x0≠0 a xN=0 minimizando el funcional de costo definido en la Ec. (4-2). Se resuelve utilizando un método basado en Programación dinámica. 6.1.2. Ejemplo de tres etapas Se presenta un problema de control óptimo definido mediante el funcional de costo J (x , u ) =

∑ [x 2k + u 2k ]+ x 32 , 2

k =0

siendo la ecuación del sistema x (k + 1) = x (k ) + u (k ) .

Nótese que el funcional de costo define que la evolución del proceso está determinada en 3 etapas. Solución: Se tiene que J * ( x ,3) = x 3 2 .

Se comienza calculando J*(x,2)

{

}

{

}

J * ( x ,2) = min x 2 2 + u 2 2 + J * (x ,3) = min x 2 2 + u 2 2 + x 3 2 . u

u

Reemplazando mediante la expresión del modelo dinámico para poner todos los términos en función de las variables del instante k=2,

{

}

J * (x ,2 ) = min x 2 2 + u 2 2 + (x 2 + u 2 )2 . u

Diferenciando la cantidad entre llaves con respecto a u e igualando el resultado a cero, se obtiene

{

}

d x 2 2 + u 2 2 + (x 2 + u 2 )2 = 2u 2 + 2(x 2 + u 2 ) = 2x 2 + 4u 2 = 0 . du Usando este resultado y la convexidad de la función x 2 2 + u 2 2 + (x 2 + u 2 )2 se deduce que el valor 1 de u que la minimiza es u o2 = − x 2 . Sustituyendo se obtiene 2 2

2

1  3  1   J * ( x ,2) = x 22 +  − x 2  +  x 2 − x 2  = x 2 2 . 2  2  2  

Para obtener J*(x,1), la ecuación es entonces 3   J * ( x ,1) = min  x12 + u 12 + (x1 + u 1 )2  . u  2  3 Siguiendo el mismo procedimiento analítico se encuentra que u 1o = − x . Sustituyendo se obtiene 5

8 J * ( x ,1) = x 2 . 5 * Finalmente la ecuación recursiva para J (x,0) es J. A. Pucheta (www.labimac.blogspot.com)

34

Control óptimo y procesos estocásticos

8   J * ( x ,0) = min x 02 + u 02 + (x 0 + u 0 )2  , u  5  8 minimizando se obtiene u o0 = − x 0 y la función de costo resultante es 13 21 J * ( x ,0) = x 02 . 13

La especificación completa de J*(x,k) y de uo(x,k), k=0,1,2 se indica en la Tabla 6-1. k 0

J*(x,k) uo(x,k) 2 1,615x (0 ) − 8 x (0 ) 13

1

1,600 x (1)

2

1,500 x (2 )

3

x (3 )

2

2

3 − x (1) 5 1 − x (2 ) 2

2

Tabla 6-1. Solución explícita para la función de costo mínimo y la secuencia óptima de control del ejemplo

En la sección siguiente, se va a generalizar la metodología para un planteo con referencia variante en el tiempo distinta de cero. 6.1.3. Ejemplo de tres etapas con referencia rk Se presenta un problema de control óptimo definido mediante el funcional de costo J (x , u ) =

∑ [(x k − rk ) 2

2

k =0

]

+ u 2k + (x 3 − r3 )2 ,

siendo la ecuación del sistema x (k + 1) = x (k ) + u (k ) .

Nótese que el funcional de costo define que la evolución del proceso está determinada en 3 etapas. Solución: Se tiene que

J* ( x ,3) = (x 3 − r3 )2 .

Se comienza calculando J*(x,2)

{

}

{

}

J* ( x ,2) = min J (x , u ) = (x 2 − r2 )2 + u 22 + J* (x ,3) = min (x 2 − r2 )2 + u 22 + (x 3 − r3 )2 . u

u

Reemplazando mediante la expresión del modelo dinámico para poner todos los términos en función de las variables del instante k=2,

{

}

J* (x ,2) = min (x 2 − r2 )2 + u 2 2 + (x 2 + u 2 − r3 )2 . u

Diferenciando la cantidad entre llaves con respecto a u e igualando el resultado a cero, se obtiene

{

}

d (x 2 − r2 )2 + u 2 2 + (x 2 + u 2 − r3 )2 = 2u 2 + 2(x 2 + u 2 − r3 ) = 2x 2 + 4u 2 − 2r3 = 0 . du Usando este resultado y la convexidad de la función (x 2 − r2 )2 + u 22 + (x 2 + u 2 − r3 )2 se deduce J. A. Pucheta (www.labimac.blogspot.com)

35

Control óptimo y procesos estocásticos

que el valor de u que la minimiza es u o2 = −0,5x 2 + 0,5r3 .

Sustituyendo se obtiene J * ( x , 2) = + (x 2 − r2 )2 = x 2 2 − 2x 2 r2 + r22

[

+ [0,5(− x 2 + r3 )]2 = 0,25 (x 2 )2 − 2 ⋅ x 2 ⋅ r3 + r32

]

= 0,25x 2 2 − 0,5x 2 r3 + 0,25r32

[

+ (x 2 − 0,5x 2 + 0,5r3 − r3 )2 = [0,5(x 2 − r3 )]2 = 0,25 (x 2 )2 − 2 ⋅ x 2 ⋅ r3 + r32

]

.

= 0,25x 2 2 − 0,5x 2 ⋅ r3 + 0,25r32 = 1,5x 2 2 − x 2 (2r2 + r3 ) + r22 + 0,5r32

J * ( x ,2) = 1,5x 2 2 − x 2 (2r2 + r3 ) + r22 + 0,5r32

Para obtener J*(x,1), la ecuación es entonces

{

}

J * ( x ,1) = min (x1 − r1 )2 + u12 + 1,5x 2 2 − x 2 (2r2 + r3 ) + r22 + 0,5r32 . u

{

d (x1 − r1 )2 + u12 + 1,5(x1 + u1 )2 − (x1 + u1 )(2r2 + r3 ) + r22 + 0,5r32 du1

}

= 2u1 + 1,5 ⋅ 2(x1 + u 1 ) − (2r2 + r3 ) = 0 2u1 + 3(x1 + u1 ) − (2r2 + r3 ) = 0

Siguiendo el mismo procedimiento analítico se encuentra que u1o = −0,6x1 + 0,2(2r2 + r3 ) .

Sustituyendo en

J * ( x ,1) = (x 1 − r1 )2 + u 12 + 1,5x 2 2 − x 2 (2r2 + r3 ) + r22 + 0,5r32

J* ( x ,1) = (x1 − r1 )2 + u12 + 1,5(x1 + u1 )2 − (x1 + u1 )(2r2 + r3 ) + r22 + 0,25r32 = + (x1 − r1 )2 = x12 − 2 x1r1 + r12 + (− 0,6 x1 + 0,2(2r2 + r3 ))2 = (0,6 x1 )2 − 2 ⋅ 0,6 x1 0,2(2r2 + r3 ) + 0,04(2r2 + r3 )2 = 0,36x12 − 0,24x1 (2r2 + r3 ) + 0,04(2r2 + r3 )2

[

+ 1,5(0,4 x1 + 0,2(2r2 + r3 ))2 = 1,5 (0,4 x1 )2 + 2 ⋅ 0,4 x1 0,2(2r2 + r3 ) + 0,04(2r2 + r3 )2 = 0,24x12

+ 0,24x1 (2r2 + r3 ) + 0,06(2r2 + r3 )

]

2

− (0,4 x1 + 0,2(2r2 + r3 )) ⋅ (2r2 + r3 ) = −0,4 x1 (2r2 + r3 ) − 0,2(2r2 + r3 )2 + r22 + 0,25r32

J* ( x ,1) = 1,6 x12 − x1[2r1 + 0,4(2r2 + r3 )] + r12 − 0,1(2r2 + r3 )2 + r22 + 0,25r32

se obtiene J* (x ,1) = 1,6 x12 − x1[2r1 + 0,4(2r2 + r3 )] + r12 − 0,1(2r2 + r3 )2 + r22 + 0,5r32 J. A. Pucheta (www.labimac.blogspot.com)

36

Control óptimo y procesos estocásticos

. *

Finalmente la ecuación recursiva para J (x,0) es

(

J* ( x ,0) = min  x 0 − r0 u 

) +u 2

2 0

+ J* (x ,1) , 

minimizando respecto de u0,

1,6(x 0 + u 0 )2 − (x 0 + u 0 )[2r1 + 0,4(2r2 + r3 )] + r12 − 0,1(2r2 + r3 )2 + r22 + 0,5r32 2u 0 + 2(1,6 )(x 0 + u 0 ) + [− 2r1 − 0,4(2r2 + r3 )] = 0

5,2u 0 + 3,2 x 0 − [2r1 + 0,4(2r2 + r3 )] = 0 u0 =

− 3,2 x 0 + [2r1 + 0,4(2r2 + r3 )] 5,2

se obtiene

u 0 = −0,6153x 0 + 0,1923[2r1 + 0,4(2r2 + r3 )]

y la función de costo resultante se obtiene reemplazando a u0 en la expresión de J1*. En la sección siguiente, se va a generalizar la metodología para un sistema multivariable y con N etapas de evolución. En el programa para Matlab Sol_Analitica_ref.m está implementado el código que genera la operatoria expuesta y da los resultados directamente. En los programas Prog_3_rk_04.m y Prog_3_PD_05.m está implementada la evolución del sistema con éste controlador. 6.1.4. Solución al problema del regulador óptimo lineal determinístico. Siguiendo el principio de Optimalidad, el funcional de Ec. (4-2) se puede descomponer en dos partes J0 =

∑ [x iT Qx i + u iT Ru i ]+ x TN-1Qx N-1 + u TN-1Ru N-1 + x TN Sx N .

N -2

(6-1)

i =0

El funcional de (6-1) será minimizado respecto de uN-1, por cuanto la primer sumatoria no depende de uN-1. Entonces se obtiene

{

}

min J N-1 = min x TN -1Qx N -1 + u TN-1Ru N-1 + x TN Sx N . u N −1

u N −1

(6-2)

Luego se reemplaza el valor obtenido de u en la expresión del funcional para obtener el costo mínimo, expresado como J*N-1, que estará en función del modelo dinámico del proceso y de las matrices de ponderación. Haciendo nuevamente el mismo procedimiento, se obtiene la acción u para el instante uN-2 lo que permitirá obtener el J*N-2. Repitiendo éste procedimiento hasta el instante k, se obtiene la expresión de la ley de control óptima que generará el costo mínimo J*. Entonces, definiendo al funcional de costos en el instante k como función del funcional en el instante k+1 se tiene, J k = J( x k , u k ) = x(k) Q x(k) + uT (k) R u(k) + J k +1

(6-3)

T

Para k+1=N se tiene J N−1 = x(N - 1) Q x(N - 1) + u T (N - 1) R u(N - 1) + x T ( N )Sx ( N ) T

J. A. Pucheta (www.labimac.blogspot.com)

(6-4)

37

Control óptimo y procesos estocásticos

donde JN = xT(N) S x(N), y el estado x(N) es x(N) = A x(N - 1) + B u(N - 1).

(6-5)

Reemplazando la Ec. (6-5) en la Ec. (6-3) y haciendo min J = J(x (N - 1), u °(N - 1)) . ( ) N−1 u N −1

Tener en cuenta las condiciones (3-20). Se tiene

∂J N −1 = 2 BT (N - 1) S [A x(N - 1)] + ∂u(N - 1)

(6-6)

+ 2 BT S B u(N - 1) + 2 R u(N - 1) = 0.

Despejando u se obtiene el uo(N-1) -1

u °(N - 1) = - [R + BT S B] BT S A x(N - 1).

(6-7)

reemplazando el estado y la acción de control para k=N-1, que son las Ecs. (6-5) y (6-7) en la Ec. (6-4) y reordenando min J N -1 = [A x(N - 1)]T H(N - 1) ⋅ [A x(N - 1)] + x T (N - 1) S x(N - 1),

(6-8)

H(N - 1) = S - S B [ R + BT S B ]-1 BT S

(6-9)

u ( N −1)

con

agrupando los términos en x(N-1), y definiendo P(N - 1) = Q + AT H(N - 1) A

(6-10)

se obtiene min J = x T (N - 1) P(N - 1) x(N - 1). ( ) N -1

(6-11)

u N -1

Se procede a continuación a minimizar a la Ec (6-3) en k=N-2, J N-2 = xT (N - 2) Q x(N - 2) + u T (N - 2) R u(N - 2) + xT (N - 1) P(N - 1) x(N - 1)

(6-12)

con x(N - 1) = A x(N - 2) + B u(N - 2)

(6-13)

como las Ecs. (6-12) y (6-13) son análogas a las Ecs. (6-4) y (6-5), se opera repitiendo los pasos desde Ec. (6-6) hasta Ec. (6-11) obteniéndose u °(N - 2) = -[R + BT P(N - 1) B ]-1 ⋅

(6-14)

⋅ BT P(N - 1) A x(N - 2) H(N - 2) = P(N - 1) - P(N - 1) B ⋅

[

⋅ R + BT P(N - 1) B J. A. Pucheta (www.labimac.blogspot.com)

]-1 BT P(N - 1)

(6-15) 38

Control óptimo y procesos estocásticos

P(N - 2) = Q + AT H(N - 2) A

(6-16)

P(N - 2) = Q + A T P(N - 2) A - A T P(N - 1) B ⋅

(6-17)

⋅ [R + BT P(N - 1) B ]-1 BT P(N - 1) A.

donde la diferencia está en cambiar a S por P(N-1). Siguiendo con la inducción, para un k cualquiera,

[

u °(k) = - R + BT (k) P(k + 1) B

]-1 ⋅ BP(k + 1)x(k),

(6-18)

H(k) = P(k + 1) - P(k + 1) B [R + BT P(k + 1) B ]-1 BT P(k + 1)

(6-19)

P(k + 1) = Q + A T H(k + 1)A P(k) = Q + A T P(k + 1) A - AT P(k + 1)B

[R + BT P(k + 1) B]

−1

(6-20)

BT P(k + 1) A.

(6-21)

J*k = xT (k) P(k) x(k).

la Ec. (6-20) s la ecuación matricial de Riccati y se debe resolver para obtener el vector de control. La condición inicial para P es P(N) = S . (6-22) Se ha obtenido así, la expresión para

u o ( k ) = −K (k )x (k ) . Con el proceso lineal, variante o invariante, donde K(k) es

[

K(k) = R + BT P(k + 1) B

]-1 BT P(k + 1) A.

(6-23)

Si el sistema es invariante, entonces se puede obtener P en estado estacionario, esto se cumple haciendo

P = lim P(k) k →∞

y por lo tanto se obtiene

(

(6-24)

)

-1

P = Q + A T PA - A T PB R + B T PB B T PA

(6-25)

que es la ecuación matricial de Riccati en estado estacionario. Nótese que igual a la Ec. (4-22). Entonces la ecuación del controlador será

(

)

-1

u °(k) = - R + B T PB B T PAx k .

(6-26)

Para el caso del proceso invariante, se puede obtener la ecuación matricial de Riccati de estado J. A. Pucheta (www.labimac.blogspot.com)

39

Control óptimo y procesos estocásticos

estacionario tomando límite sobre P(k). De esta manera, el controlador de estado estacionario queda u ( k ) = − Kx (k ) . Donde el controlador es sub óptimo en las transiciones, pero óptimo al alcanzar la referencia.

6.1.5. Extensión a los sistemas multivariable La metodología de diseño del controlador desarrollada para los sistemas mono variable, puede ser extendida a los sistemas multivariable definidos en la Sección 1.1.1. Para ello, se formula el problema de control considerando que el sistema (1-7) y expresando a la acción de control mediante u t = −K x t

(6-27)

donde K∈Rm×n, u y x son funciones temporales con magnitudes u∈Rm y x∈Rn, respectivamente. Para ello, se define b = Bh

(6-28)

siendo B∈Rn×e y h∈Re es un vector que se elige de manera tal que el par (A, b) sea controlable. Ahora se tiene un sistema mono variable con entrada ut x& t = A x t + b u t

(6-29)

y se busca el vector kT de realimentación para obtener a ut ut = kT xt

(6-30) T

de manera que los autovalores de A-BK sean los mismos que los de A-bk . Igualando

[A − B K ] = [A − B h k T ]

(6-31)

K = h kT

(6-32)

se concluye por comparación que Se ha convertido el problema multivariable en un problema de una simple entrada, y por lo tanto se puede resolver con las ecuaciones del caso para la obtención de kT. El vector h no es único, por lo que habrá más de una matriz de realimentación que cumpla con la asignación de polos requerida, la única restricción para h es que el par (A, Bh) sea completamente controlable.

7. Programación dinámica 7.1. Versión simbólica: Ecuación de Hamilton-Jacobi-Bellman La Ecuación de HJB es la versión en el dominio del tiempo contínuo de la Ecuación de Bellman mostrada en la Ec. (7-37) y Ec. (7-38), y es una ecuación diferencial no lineal en derivadas parciales. Aquí, se formula el problema del control óptimo en el dominio del tiempo continuo proponiendo un funcional de costos de la forma

(

)

tf

J (x , u ) = h x t f , t f + ∫ g (x τ , u τ , τ )dτ

(7-1)

t0

J. A. Pucheta (www.labimac.blogspot.com)

40

Control óptimo y procesos estocásticos

El objetivo de control consiste en hallar la ley de control que haga evoluciona al proceso modelado mediante la expresión (7-2) minimizando al funcional de costos (7-1). El modelo del sistema es (7-2) x& = a [x , u , t ], t ∈ [ t , t ] t

t

t

0

f

donde x(t0) tiene un valor fijo C. La minimización del funcional (7-2) para todo t∈[t0,tf], será tf   J * (x t , t ) = min h x t f , t f + ∫ g (x τ , u τ , τ )dτ u τ∈[ t , tf ]   t

(

)

(7-3)

Subdividiendo el intervalo en ∆t tf t + ∆t J * (x t , t ) = min  ∫ g (x τ , u τ , τ )dτ + ∫ g (x τ , u τ , τ )dτ + h x t f , t f u τ∈[ t , t + ∆t ]  t t + ∆t donde se aplicará el Principio de optimalidad de Bellman.

(

(7-4)



∆t •

t

•



)

•

t0 •

tf

Fig. 11. Partición de la trayectoria óptima desde t hasta tf, en un ∆t.

Se obtiene, entonces t + ∆t  (7-5) J * (x t , t ) = min  ∫ g (x τ , u τ , τ )dτ + J * (x t + ∆t , t + ∆t ). u τ∈[ t , t + ∆t ]  t  Asumiendo que existe derivada segunda de J* en la Ec. (7-5), y que es acotado, se hace la expansión mediante serie de Taylor en el punto (xt,t). Se tiene,

t + ∆t  ∂J * (x t , t )  J * (x t , t ) = min  ∫ g (x τ , u τ , τ )dτ + J * (x t , t ) +   ∆t + u τ∈[ t , t + ∆t ] ∂t    t

(7-6)

T   ∂J * (x t , t ) [ x − x ] + tér min os de orden superior + . t + ∆ t t  ∂x   

Haciendo ∆t pequeño, se tiene que,

{

J * (x t , t ) = min g (x t , u t , t )∆t + J * (x t , t ) + J *t (x t , t )∆t + ut

+

J *x

(7-7)

(x t , t ) [a (x t , u t , t )]∆t + términos de orden superior}.

J. A. Pucheta (www.labimac.blogspot.com)

T

41

Control óptimo y procesos estocásticos

donde se definen

J *x

T

∂J *  ∂J * = = ∂x  ∂x 1

 ∂J * ∂J * . L , y J *t = ∂x 2 ∂t 

Los términos J* y Jt* pueden sacarse de las llaves de minimización porque no dependen de ut. Se tiene,

J * (x t , t ) = J * (x t , t ) + J *t (x t , t )∆t + min{g (x t , u t , t )∆t + ut

+

J *x

(7-8)

(x t , t ) [a (x t , u t , t )]∆t + términos de orden superior}. T

Se cancelan los términos en J*,

0 = J *t (x t , t )∆t + min{g (x t , u t , t )∆t

(7-9)

ut

+ J *x (x t , t ) [a (x t , u t , t )]∆t + términos de orden superior }. T

Dividiendo por ∆t y tomando límite cuando ∆t→0, considerando que términos de orden sup erior lim = 0, ∆t → 0 ∆t

Se llega a

0 = J *t (x t , t ) + min{g (x t , u t , t ) + J *x (x t , t ) [a (x t , u t , t )] T

(7-10)

ut

con la condición de frontera para t=tf,

(

J * (x t , t f ) = h x t f , t f

)

(7-11)

Definiendo al Hamiltoniano H como

(

)

H x t , u t , J *x , t = g (x t , u t , t ) + J *x (x t , t ) [a (x t , u t , t )] y también se define a

(

T

)

{(

H x t , u*, J *x , t = min H x t , u t , J *x , t

(7-12)

)}

(7-13)

que debido a la acción minimizante dependerá de x, J*x y t. Así, se llega a la Ecuación de Hamilton Jacobi Bellman, expresada en función del Hamiltoniano

(

0 = J *t (x t , t ) + H x t , u*, J *x , t

)

(7-14)

que es la versión en tiempo continuo de la Ecuación de Bellman, expresada en la Ec. (7-37).

7.1.1. Ejemplo de aplicación de la Ecuación de HJB Sea el sistema (7-15)

x& t = x t + u t

hallar la ley de control ut que minimiza al funcional de costos T

J. A. Pucheta (www.labimac.blogspot.com)

1 1 J = x 2 (T ) + ∫ u 2t dt. 4 4 0

(7-16) 42

Control óptimo y procesos estocásticos

Aquí,

1 2 ut 4

(7-17)

a = xt + ut

(7-18)

g=

y también

entonces se construye el Hamiltoniano como 1 2 u t + J *x (x t , u t ) 4

(7-19)

∂H 1 = u t + J *x (x t , t ) ∂u 2

(7-20)

H=

la derivada respecto de u es

y la derivada segunda de H es

∂ 2H ∂u 2

(7-21)

1 >0 2

=

por lo que igualando a 0 la Ec. (7-20) la u que se obtendrá corresponderá a un mínimo de J de la Ec (7-16). Se obtiene así, que (7-22) u *t = −2J *x (x t , t ). Ahora se construye la Ec HJB

resulta

0 = J *t +

[

1 − 2J *x 4

] + [J ]x 2

* x

[ ] + [J ]x

0 = J *t − J *x que para t=T se tiene

2

* x

J * (x T , T ) =

[ ]

− 2 J *x

t

2

(7-23)

(7-24)

t

1 2 xT. 4

(7-25)

1 K t x 2t 2

(7-26)

Se propone una solución, por ejemplo J * (x t , t ) =

donde Kt es una función temporal, incógnita.

J *x (x t , t ) = K t x t

(7-27)

u *t = −2K t x t

(7-28)

que por la Ec. (7-22) se tiene que

Tomando KT=1/2, J *t (x t , t ) =

1 & 2 Ktxt 2 y mediante la HJB Ec (7-24), y debido a que debe resolverse para todo xt, se llega a 0=

Una solución es

(7-29)

1 & K t − K 2t + K t . 2

(7-30)

exp(T − t ) (7-31) exp(T − t ) + exp(t − T ) Nótese que si T>>t, K es aproximadamente 1, y el sistema es estable ya que reemplazando en la Ec Kt =

J. A. Pucheta (www.labimac.blogspot.com)

43

Control óptimo y procesos estocásticos

(7-15) se tiene x& t = − x t .

(7-32)

Nótese que la solución para un problema con un sistema modelado con una variable de estado, y un funcional de costo cuadrático es muy laborioso y tedioso, incrementando la dificultad si el sistema es modelado como no lineal y el funcional de costo propuesto no es cuadrático. 7.2. Versión numérica: Ecuación de Bellman Hasta aquí no se han considerado las restricciones, ni tampoco las no linealidades en el modelo del proceso a controlar ni funcionales de costo no cuadráticos. Sin embargo, en determinadas situaciones es conveniente tener en cuenta las restricciones en el proceso para el diseño del controlador. Cuando se trata con procesos cuyos modelos no son del tipo lineal o se deben considerar saturaciones en los actuadores o en las variables de estado, sucede que la expresión analítica cerrada de la solución al problema de control óptimo no siempre puede hallarse, por lo que se hace uso de la aproximación numérica de la ley de control mediante la cuantificación de los estados de la planta. 7.3. Problema básico Para formular el problema, se presentan las expresiones del modelo del proceso, las restricciones en las variables y el funcional de costo a minimizar. Se considera el problema de minimizar la función de costos separable J = ∑ L[x (k ), u (k ), k ] N

(7-33)

k =0

donde x(0) tiene un valor fijo C y deben ser satisfechas la ecuación del sistema

x ( k + 1) = f [x ( k ), u ( k ), k ], k = 0,1,... N − 1 y las restricciones

(7-34)

x∈X ⊂ Rn ,

(7-35) (7-36)

u∈U ⊂ R . m

Por simplicidad se supone que f y L son funciones acotadas y continuas de sus argumentos, y que x y u pertenecen a subconjuntos cerrados y acotados de Rn y Rm, respectivamente. Entonces, el teorema de Weierstrass asegura que existe una política de minimización. u(k)=µ(x(k),k)

Sistema x(k+1)=f(x(k),u(k),k)

x(k)

µ(x(k),k)

Fig. 12. Implementación del controlador basado en programación dinámica numérica.

Se desea hallar una función µ(x (k ), k ) : ℜn → ℜm , que haga evolucionar al proceso modelado J. A. Pucheta (www.labimac.blogspot.com)

44

Control óptimo y procesos estocásticos

mediante la Ec (7-34) desde cualquier condición inicial hasta el estado final x(N) cumpliendo con las restricciones (7-35) - (7-36), y que minimice al funcional de costo (7-33). Aplicando el principio de optimalidad, se obtiene

{

}

min J(x, k) = min L(x, u(k), k) + J*(f (x, u(k), k ), k + 1) , u(k )

u(k )

(7-37)

denominada Ecuación de Bellman., y por lo tanto la acción de control u será

{

}

uo (x, k ) = arg min L(x, u(k), k ) + J*(f (x, u(k), k ), k + 1) . u(k )

(7-38)

7.4. La política óptima de decisiones Se desea obtener la ley de control mediante el cómputo de

{

}

u o (x , k ) = arg min L(x, u (k ), k ) + J * (f (x, u ( k ), k ), k + 1) . u(k )

Para obtener la ley de control se proponen tres métodos, conocidos como Programación dinámica regresiva (Clásica), Programación dinámica Iterativa, y Programación dinámica aproximada. 7.5. Programación dinámica regresiva Mediante éste método, se obtiene una tabla de valores con dos entradas: x y k, conociendo que J * = min{L(x , u ( N ), N )} u(N) y luego se resuelve numéricamente mediante la programación dinámica. En la metodología del ejemplo 6.1.2, cuyos resultados están en la Tabla 6-1, puede verse que hay dos pasos: uno es para calcular la secuencia óptima de decisiones y otra para calcular la trayectoria óptima del estado del proceso para cada valor numérico de x(0).

7.5.1. Procedimiento de evaluación numérica Sea el funcional de costo J=

4

∑ [2 + u (k)] ⋅ e− x (k ) + x (5) − 1. ,

k =0

el sistema expresado por

5 1   x (k + 1) = x (k ) + 2 − 2 ⋅ x (k ) + ⋅ x 2 (k ) − ⋅ x 3 (k ) ⋅ u (k ) 4 4   con las siguientes restricciones 0 ≤ x ≤ 3, − 1 ≤ u ≤ 1. utilizando una cuantificación uniforme con ∆x=1 y ∆u=1. Se desea: -Hallar la solución completa por programación dinámica. -Generar la solución para el estado inicial x(0)=2. J. A. Pucheta (www.labimac.blogspot.com)

45

Control óptimo y procesos estocásticos

-Aplicar interpolación lineal para cualquier interpolación requerida. Solución: Para los incrementos de cuantificación indicados, el conjunto de estados admisibles es x={0,1,2,3} y el conjunto de decisiones admisibles es u={-1,0,1}. Con el objeto de facilitar el análisis es útil calcular la función de transición de estados y la función de costo para cada etapa en cada estado cuantificado en función de la variable de decisión u. Dado que, en este ejemplo, ambas funciones son invariantes en el tiempo, ahorra tiempo de computarlas al principio y almacenarlas para referencia. Los resultados se muestran en las tablas Tabla 7-1 y Tabla 7-2. x(k) 0 1 2 3

x(k+1) 2 ⋅ u (k ) 1 + u (k ) 2 + u (k ) 3 + 0,5 ⋅ u (k )

Tabla 7-1. Función de transición de estados.

x(k) 0 1 2 3

L[x(k),u(k),k] 2+u(k) 0,36788 ⋅ [2 + u (k )] 0,13534 ⋅ [2 + u (k )] 0,04979 ⋅ [2 + u (k )]

Tabla 7-2. Índice de desempeño para cada etapa en cada estado cuantificado.

Debe notar que en el estado x = 3 la decisión u = -1 lleva a un próximo estado igual a 2,5. Dado que este estado no es uno de los estados cuantificados de X será necesaria una interpolación para obtener el costo mínimo de dicho próximo estado. Los cálculos comienzan especificando L[x,5] en los estados cuantificados como se muestra en la Fig. 7-13 donde para cada estado cuantificado de esa etapa se tiene, L[3,5] = |3-1| = 2 L[2,5] = |2-1| = 1 L[1,5] = |1-1| = 0 L[0,5] = |0-1| = 1. x=3

2

x=2

1

x=1

0

x=0 k=0

1

2

3

4

5

1

Fig. 7-13. Condiciones de borde.

Para ilustrar el procedimiento de interpolación, se consideran los tres controles aplicados al estado J. A. Pucheta (www.labimac.blogspot.com)

46

Control óptimo y procesos estocásticos

x=3 en la etapa k=4, como se muestra en la Tabla 7-3.

u 1 0 -1

g(x,u,k) 3,5 3 2,5

J*(g,k+1) * 2 1,5

L(x,u,k) * 0,09958 0,04979

Costo Total * 2,09958 1,54979

Tabla 7-3. Cálculos en x=3, k=4.

Para u = +1 el próximo estado es x = 3,5, el cual viola las restricciones de estados. Para u = 0 el próximo estado es x = 3. El costo mínimo en el próximo estado es I(3,5) = 2 y el costo de la etapa es L(3,0,4)=0,09958 lo que da un costo total de 2,09958. Para u = -1 el próximo estado es x = 2,5, un valor intermedio en la mitad entre los estados cuantificados x = 2 y x = 3. Una interpolación lineal entre el costo I(2,5)=1 e I(3,5)=2 da un costo de I(2,5,5) = 1,5. El costo de la etapa es L(3,-1,4) = 0,04979 lo que da un costo total de 1,54979. Este último costo es, claramente, el valor mínimo. La solución completa se muestra en la Fig. 7-14. x=3

0 1

x=2

-1

-1

-1

0

0

0

1

1

1

0

2.05099 1.78033 1.50965 1.23898 0.73576

x=0 1

1

1

1

k=1

k=2

k=3

1 0

0

3.67669 3.40601 3.13534 3.000 k=0

2

-1

1.14363 0.94735 0.67669 0.40601 0.13534 1

x=1

-1

0.83720 0.73762 0.69898 0.89236 1.54979

k=4

1

k=5

Fig. 7-14. Solución completa utilizando interpolación lineal.

k 0 1 2 3 4 5

X 2 3 2,50 2,05 2 1

u 1 -1 -0,50 -0,05 -1 -

L(x,u,k) 0,40602 0,04979 0,12313 0,24792 0,13534 0,0 0,96220

Tabla 7-4. Solución óptima para x(0) = 2.

Para el estado inicial x(0) = 2 la solución se muestra en la Tabla 7-4. Cabe destacar que esta solución requiere interpolación para obtener la decisión óptima para k ≥ 2. La decisión óptima para x(0) = 2 es u = 1 la cual lleva al próximo estado x(1) = 3 con un estado de etapa L(2,1,0) = 0,40602. La decisión óptima para x(1) = 3 es u=-1 la cual lleva al próximo estado x(2) = 2,5 y un costo de etapa de L(3,-1,1) = 0,04979. A partir de este estado se requiere realizar interpolaciones dado que no es uno de los estados cuantificados para los cuales ya se han calculado las decisiones óptimas. Las decisiones óptimas para los dos estados más cercanos son uˆ (2,2) = 0 y uˆ (3,2) = −1 . La ecuación general de interpolación lineal es J. A. Pucheta (www.labimac.blogspot.com)

47

Control óptimo y procesos estocásticos

uˆ (x , k ) = uˆ (a ⋅ ∆x, k ) +

para a ⋅ ∆x ≤ x ≤ (a + 1) ⋅ ∆x . En este caso x = 2.5 y ∆x = 1 , a = 2

uˆ[(a + 1)∆x, k ] − uˆ (a∆x, k ) ⋅ (x − ∆x ) ∆x

uˆ (a∆x , k ) = uˆ (2,2) = 0 uˆ[(a + 1)∆x , k ] = uˆ (3,2) = −1

Reemplazando en la ecuación general de interpolación lineal se obtiene uˆ (2.5,2) = 0 +

(− 1 − 0) (2.5 − 2) = −0,5 . 1

Estos valores del estado y de la decisión óptima pueden ahora introducirse en las ecuaciones originales del sistema y del costo de etapa. Dichos cálculos indican que el próximo estado es x(3)=2,05 y que el costo de la etapa es 0,12313. 5 1   x 3 = 2.5 + 2 − 2.5 ∗ 2 + ∗ 2.5 2 − ∗ 2.5 3 (− 0.5) ≅ 2,05 4 4   L(2.5,−0.5,3) = [2 + (− 0.5)] ∗ e −2.5 = 0,12313

La interpolación entre uˆ (2,2) = 0 y uˆ (3,3) = −1 da uˆ (2.05,3) = 0 +

−1 − 0 (2.05 − 2) = −0,05 . 1

La sustitución en las ecuaciones originales del sistema y del costo muestra que x(4)=2 y que el costo es 0,24792. A continuación se obtiene directamente que la decisión óptima uˆ (2,4) = −1 ,el costo en k=4 es 0,13534 y el estado final x(5)=1, con una penalidad (o costo) terminal de 0. Es dable destacar que el costo total a lo largo de esta trayectoria es de 0,96220 y no 1,14363 que es el valor calculado en el armado de la grilla. Esta discrepancia muestra que para ecuaciones del sistema y funciones de costo no lineales, como las aplicadas en este ejemplo, es a veces necesario utilizar intervalos de cuantificación más finos y/o interpolaciones de mayor orden para obtener resultados precisos. En el programa para Matlab PDNumerica_01.m se presenta una implementación del algoritmo propuesto.

J. A. Pucheta (www.labimac.blogspot.com)

48

Control óptimo y procesos estocásticos Estados

Estados

3 2 1 0

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

3

3.5

4

4.5

5

Costo 2 Costos

1.5 1 0.5 0

0

0.5

1

1.5

2

2.5

Acciones de control

Acción de control 1 0.5

Nmax =4 . Mmax =3

0

Nmax =45 . Mmax =45

-0.5 -1 0

0.5

1

1.5

2

2.5 Etapas

3

3.5

4

4.5

5

Fig. 15. Evolución usando la política tabulada sin interpolación.

7.6. Algunos funcionales típicos En la formulación del problema general de optimización en múltiples etapas o discreto, se distinguen tres elementos esenciales: La descripción de la planta o proceso en consideración por medio de una ecuación dinámica discreta como la expresada mediante la Ec (7-34). La presencia de k en los argumentos de g indica que la función puede, en general, variar con el tiempo (o etapa). Las dinámicas del sistema están fijadas por la física del problema. Las restricciones en los estados y en las acciones de control, expresadas mediante las Ec (7-35) (7-36), están fijadas por la física del proceso y por el ingeniero. El funcional de costo (a minimizar o índice de desempeño a maximizar), expresado mediante la Ec (7-33), es elegido exclusivamente por el ingeniero. El valor del índice L(x(k),u(k),k) para k=N normalmente es una expresión que no depende explícitamente de u(N), por lo que se puede expresar mediante una función que es una función del tiempo final N y del estado en el tiempo final N, es decir L(x (k ), u (k ), k ) := φ(x ( N ), N ) ∀k = N.

(7-39)

Representa el costo o penalidad que el usuario le designa a cada estado final admisible. Si se desea que el sistema llegue, en N etapas a un determinado estado final, se asigna costo cero a dicho estado y un costo o penalización elevada a los otros estados. De esta manera se puede lograr que el costo total de la trayectoria sea el mínimo, cuando el estado final es el estado al cual el ingeniero desea que llegue dicha trayectoria óptima en las N etapas. La función L[x(k),u(k),k] es una función que puede variar con la etapa k (a menudo el tiempo), J. A. Pucheta (www.labimac.blogspot.com)

49

Control óptimo y procesos estocásticos

diseñada por el ingeniero para alcanzar una determinada respuesta o comportamiento por parte del sistema. Para alcanzar diferentes objetivos de control se seleccionan diferentes tipos de índices de comportamiento. Algunos de los índices comúnmente utilizados se describen a continuación. 7.6.1. Problemas de tiempo mínimo Suponiendo que se desea encontrar la secuencia de control uk que lleve el sistema desde el estado inicial dado x(0) a un estado final deseado x(N) en el tiempo mínimo. Entonces se podría seleccionar el funcional N −1

J = φ[x ( N)] + N = φ[x ( N)] + ∑1 k =0

y especificar la condición de borde x(N)=xd. En este caso es L=1, y N puede o no ser una variable a minimizar. 7.6.2. Problemas de Mínimo Consumo de Combustible Cuando se desea encontrar la secuencia de control {uk} para llevar el sistema desde x(0) a un estado final deseado x(N), en un tiempo fijo N utilizando el mínimo combustible, el funcional a utilizar es N −1

J = ∑ u (k ) + φ[x ( N)] k =0

debido a que el combustible que se quema es proporcional a la magnitud del vector de control. En este caso L=u(k). Por ejemplo si el control es proporcional a la diferencia de temperaturas TD-TA entre la temperatura deseada y la temperatura ambiente, al variar TA la diferencia puede ser positiva o negativa, indicando la necesidad de aplicar calefacción o refrigeración. En ambos casos existe consumo de combustible. 7.6.3. Problemas de mínima energía Este funcional se utiliza si se desea encontrar u(k) para minimizar la energía del estado final y de todos los estados intermedios y también del control. Nuevamente suponiendo fijado el tiempo final N, se puede utilizar el funcional J=

1 N−1 T 1 ∑ [x (k ) Q x (k ) + u T (k ) R u (k )] + x T ( N) S x ( N) 2 k =0 2

donde Q, R y S son matrices de ponderación definidas positivas. 1 1 En este caso φ = x T ( N) S x ( N) y L = [x T (k ) Q x (k ) + u T (k ) R u (k )] , ambas son funciones 2 2 cuadráticas. Minimizar la energía corresponde, en cierto sentido, a mantener el estado y el control cerca de cero. Si se considera más importante que los estados intermedios sean pequeños entonces se podrá elegir qi grande para pesar los estados fuertemente en J que es el funcional que se trata de minimizar. Si es J. A. Pucheta (www.labimac.blogspot.com)

50

Control óptimo y procesos estocásticos

más importante que sea pequeña la energía de control, entonces se elegiría un valor grande de ri. Si interesa más que el estado final sea pequeño, entonces S debería ser grande. El problema de control óptimo se caracteriza por compromisos y ajustes con diferentes factores de peso en J que resultan en diferentes equilibrios entre objetivos de desempeño y magnitud de las acciones óptimas requeridas. En la práctica, es usualmente necesario hacer un diseño de control con un funcional J estimado, computar la secuencia de control óptimo uk y correr una simulación en computadora para ver como responde el sistema a esta secuencia de acciones de control. Si la respuesta no es aceptable, se repite la operación usando otro J con diferentes pesos en los estados y controles. Después de varias repeticiones para encontrar una secuencia uk aceptable, esta versión final de uk se aplica al sistema real. Las ventajas de la PDR pueden destacarse en los problemas de baja dimensionalidad, ya que es una metodología que encuentra un mínimo global, y la ley de control resultante reside en una tabla de valores. Sin embargo, para ampliar el ámbito de aplicación de la PDR a sistemas de alta complejidad y de grandes dimensiones existen dos alternativas: la Programación dinámica iterativa y la Programación dinámica aproximada.

7.7. Programación Dinámica iterativa Fue propuesta por Luus en 1990, para procesos químicos. Intenta ampliar el campo de aplicación hacia los problemas de ingeniería de grandes dimensiones. Propone implementar el algoritmo de la PDR sobre una región determinada del espacio de estados. Para aplicar la PDR, se aproxima el problema de control óptimo buscando una política de control constante por partes, generando una secuencia de decisiones como política de control que varía en forma continua, sobre P etapas de tiempo, cada una de longitud L, tal que L=

tf . P

Por lo tanto, en el intervalo tk ≤ t < tk+1, se tiene el control constante: ut = uk . El problema entonces es encontrar la secuencia u0, u1, ...., uP – 1 que minimiza el funcional de costos. Se define L = t k +1 − t k

con tP = tf y t0 = 0. 7.7.1. Algoritmo de la PDI Paso 1: Dividir el intervalo de tiempo tf en P etapas de tiempo, cada una de longitud L. Paso 2: Elegir el número N de puntos de la grilla x y el número M de valores admisibles para cada una de las variables de control uj. Paso 3: Elegir la región rj para cada una de las variables de control. Paso 4: Eligiendo Nu (impar) valores de control dentro de la región admisible, evaluar la Ec. (7-34) del modelo dinámico de proceso Nu veces para generar la grilla x en cada etapa de tiempo. J. A. Pucheta (www.labimac.blogspot.com)

51

Control óptimo y procesos estocásticos

Paso 5: Comenzando en la última etapa de tiempo P, correspondiente al tiempo (tf – L), para cada punto de la grilla x evaluar la Ec. (7-34) desde (tf – L) hasta tf para todos los Mum valores admisibles de control. Elegir el control u que minimiza el funcional de costo y almacenar el valor del control para usarlo en el paso 6. Paso 6: Retroceder a la etapa P – 1, correspondiente al tiempo (tf – 2L) y evaluar a la Ec. (7-34) desde (tf – 2L) hasta (tf – L) para cada punto de la grilla x con los Mm valores admisibles de control. Para continuar la integración desde (tf – L) hasta tf elegir el control del paso 5 que corresponde al punto de la grilla más cercano al x resultante en (tf – L). Comparar los Mm valores del funcional de costo y almacenar el valor de control que da el mínimo valor. Paso 7: Continuar el procedimiento hasta alcanzar la etapa 1 correspondiente al tiempo inicial t=0. Almacenar la política de control que minimiza el funcional de costo y almacenar la trayectoria x correspondiente. Paso 8: Reducir la región de valores de control admisibles por un factor ε, o sea rj+1 = (1 − ε )rj

donde j es el índice de iteración. Usar la trayectoria x óptima del paso 7 como puntos medios para la grilla x en cada etapa de tiempo, y usar la política de control óptima del paso 7 como puntos medios para los valores admisibles del control u en cada etapa de tiempo. Paso 9: Incrementar el índice de iteración j en 1 y vuelva al paso 4. Continúe iterando por un número especificado de iteraciones tal como 20 y verificar los resultados. 7.7.2. Comentario Como ventaja, se logró obtener una secuencia de acciones de control empleando menores recursos computaciones que en el caso de la PDR. Sin embargo, ésta técnica permite encontrar una solución para el problema de control óptimo dependiente del estado inicial del proceso, por lo cual en la operación en línea requiere de un equipamiento capaz de implementar la PDI en línea y en el tiempo de muestreo impuesto por el proceso real y la especificación de control. 7.8. Programación dinámica aproximada Existen varias posibilidades de disminuir la cantidad de recursos según la dimensionalidad y naturaleza del problema. A continuación se van a describir dos métodos muy difundidos, que son los métodos de Iteración de política aproximada aproximando el funcional de costos y aproximando el funcional y la ley de control. 7.8.1. Aprendizaje Q Es un método computacional alternativo que puede usarse cuando no se dispone de un modelo explícito del proceso. Se supone que es dato el valor del costo asociado a las transiciones de estados. Se define para un par estado acción (i,u) el factor Q mediante Q µ (i, u ) = L(i, u , j) + J µ ( j), J. A. Pucheta (www.labimac.blogspot.com)

52

Control óptimo y procesos estocásticos

(7-40) donde i es el estado actual, j es el estado alcanzado una vez aplicado la acción u, µ es la ley o política de control para las transiciones desde el estado j hasta el final del proceso y Jµ es el costo asociado formado por la sumatoria de los costos parciales. A partir de éstos valores, se mejora la política µ mediante µ (i ) = arg min Q µ (i, u ). u∈U (i )

(7-41)

A ésta política mejorada, se le vuelve a computar los costos asociados y se calculan nuevamente los factores Q mediante la Ec. (7-40). Repitiendo éste procedimiento, se llega a la política óptima cuando no hay más cambio en los factores Q o en la mejora de la política. Se actualizan los valores de los factores Q asociados a una política, evitando evaluar la política en el proceso. Se definen los factores Q óptimos Q*(i,u) correspondientes al par (i,u) mediante Q * (i, u ) = L(i, u , j) + J * ( j),

(7-42)

donde la Ecuación de Bellman puede escribirse como J * (i ) = min Q * (i, u ). u∈U (i )

(7-43)

Combinando las ecuaciones Ec. (7-42) y la Ec. (7-43), se tiene Q * (i, u ) = L(i, u , j) + min Q * ( j, v ). v∈U ( j)

(7-44)

Los factores Q óptimos Q*(i,u) son la solución única del sistema de la Ec. (7-44). El algoritmo se escribe como Q (i, u ) = L (i, u , j) + min Q( j, v ). v∈U ( j)

(7-45)

y en una forma más general,   Q n +1 (i, u ) = (1 − γ )Q n (i, u ) + γ  L (i, u , j) + min Q n ( j, v ). v∈U ( j)  

(7-46)

La función γ cambia de una iteración a otra para el mismo par (i,u). Se demuestra que la convergencia estará asegurada si se cumple que ∞

∞

n =0

n =0

∑ γ n (i, u ) = ∞, ∑ γ 2n (i, u ) < ∞,

∀i, u ∈ U (i ),

(7-47)

entonces Qn(i,u) converge a Q*(i,u). Se supone que existe una política óptima y que Q está acotado para todo su dominio. En el programa para Matlab PD_Q_731.m se implementó un algoritmo que calcula la política de control óptima para el Ejemplo 7.3.1, el mismo se muestra en la Fig. 16. La función γ tiene la forma J. A. Pucheta (www.labimac.blogspot.com)

53

Control óptimo y procesos estocásticos

γn =

10 10 + n

donde n indica la cantidad de veces que se actualiza al par estado acción (i,u). En la Fig. 17 se muestran gráficamente los resultados obtenidos.

J. A. Pucheta (www.labimac.blogspot.com)

54

Control óptimo y procesos estocásticos

% Programación dinámica. % Apredizaje Q. % Para el Ejemplo 7.3.1. % x(k+1)=x(k)+(2-2*x(k)+5/4*x(k)^2-1/4*x(k)^3)*u(k); % con el funcional de costo J=sum((2+u(k))*exp(-x(k))); %Autor JAP %06 12 07 clear,clc,close all; TM=200; Mmax=15; color='.-k'; tic; du=Mmax; etapas=6; xmin=0; xmax=3;umin=-1; umax=1; %%%Carga de datos rand('state',0); equis1=3*(rand(TM,1)); tiempo=ceil((etapas-1)*rand(TM,1)); M_est = [tiempo,equis1]'; Au=(umax-umin)/(Mmax-1); for i=1:Mmax uf(i)=umin+Au*(i-1); end CI=2;vfinal=1; Q = zeros(TM,du); J=zeros(1,TM); sal(1)=CI;costo(1)=0; Ya=zeros(1,TM); for k=1:etapas-1 entrada = [k; sal(k)]; consigna(k) = pol_tab_mu1(entrada,M_est,Ya); sal(k+1)=mopdm(k,sal(k),consigna(k)); costo(k+1)=costo(k)+indice(k,sal(k),consigna(k)); end costo(k+1)=costo(k+1)+abs(sal(k+1)-vfinal); evoluc(1)=costo(etapas); costo(etapas) m=zeros(TM,du); for iterac=1:10 for iq=1:TM k=M_est(1,iq); x=M_est(2,iq); %Recorre todas las acciones for acc=1:du xy=mopdm(k,x,uf(acc)); m(iq,acc)=m(iq,acc)+1; gama=10/(10+m(iq,acc)); if k<etapas-1 [aux lugar]=min(abs(((k+1)-M_est(1,:))/max(M_est(1,:)))+abs((xyM_est(2,:))/max(M_est(2,:)))); Q(iq,acc)=(1-gama)*Q(iq,acc)+gama*(indice(k,x,uf(acc))+J(lugar)); else Q(iq,acc)=(1-gama)*Q(iq,acc)+gama*(indice(k,x,uf(acc))+abs(xy-vfinal)); end end end for iq=1:TM [val lugar]=min(Q(iq,:)); J(iq)=val; Ya(iq)=uf(lugar); end sal(1)=CI; costo=0; for k=1:etapas-1 entrada = [k; sal(k)]; consigna(k) = pol_tab_mu1(entrada,M_est,Ya); sal(k+1)=mopdm(k,sal(k),consigna(k));

Fig. 16. Código para implementar el algoritmo de aprendizaje Q en Matlab.

J. A. Pucheta (www.labimac.blogspot.com)

55

Control óptimo y procesos estocásticos

Estados

Costo

3

4

2.5 3 Costos

Estados

2 1.5

2

1 1 0.5 0

0

1

2

3

4

Acción de control

0

1

2

3

4

5

8

10

Evolucion del J(0)

10

10

1 Acciones de control

0

5

0.5

5

10 0

0

10 -0.5

-1

-5

0

1

2

3

4

5

10

0

2

Etapas

4 6 Iteraciones

Fig. 17. Resultados obtenidos mediante el aprendizaje Q para el Ejemplo 7.3.1.

7.8.2. Aproximación del funcional de costo La ecuación de Bellman de la programación dinámica indica que

{

}

J * ( x , k ) = min L(x , u k , k ) + J * (f (x , u k , k ), k + 1) , uk

(7-48)

donde f es el modelo dinámico del proceso x k +1 = f (x k , u k , k ),

k = 0,1,2...N - 1.

(7-49)

Como muestra la Ec. (7-48), si se dispone de la función J* en el estado siguiente al xk, se podría encontrar la acción de control óptima uk. Sin embargo, como es bien sabido la minimización de la Ec. (7-48) es laboriosa y demanda recursos computacionales cuando la dimensionalidad del problema aumenta. Una alternativa al problema de la dimensionalidad consiste en utilizar una función los valores de la Ec. (7-48) en un dominio compacto. Así, se obtiene una representación compacta del costo asociado a cada estado del proceso, el esquema de la Fig. 18 muestra esta característica.

J. A. Pucheta (www.labimac.blogspot.com)

56

Control óptimo y procesos estocásticos

Aproximación del

i

~ J (i )

~ costo J µ (i, r ) Fig. 18. Expresión

compacta de aproximador de costo

Como se observa en la Fig. 18, encontrando el valor del vector r adecuado, se dispone del valor aproximado del costo mínimo en el que se incurre para llegar al estado de costo nulo partiendo del estado actual i, y si se dispone del modelo del sistema se podrá encontrar la política de control mediante

{

}

~ J * ( x , k ) = min L(x , u k , k ) + J * (f (x , u k , k ), k + 1, r ) . uk

(7-50)

7.8.3. Procedimiento de ajuste Para encontrar r, se divide en dos tareas el proceso de búsqueda que encuentra la función de política, como se muestra en la Fig. 19, una es la evaluación de una política estacionaria definida a partir de la que se calculan los costos para todos los estados del proceso, y la otra tarea es la mejora de la política. Ambas tareas se realizan de manera aproximada con respecto al sistema original, porque se utiliza una función que ajusta su comportamiento. La función de aproximación podrá ser una red de neuronas y por lo tanto la aproximación en la evaluación y mejora de la política se debe a su utilización. Actualización de política

Jµ

µ

µ

Evaluación aproximada de política

Fig. 19. Esquema

del proceso de búsqueda de la política óptima

Implementación ~

~

Se dispone de un conjunto de datos representativos S y para cada estado i ∈ S se calculan los valores de costo c(i), de la función Jµ(i), luego minimizando en r

(

)

2 ~ min ∑ J (i, r ) − c(i ) r

~ i∈S

(7-51)

se obtiene una función de aproximación para el costo asociado a la política evaluada. Minimizando la expresión (7-51) se encuentra el vector de parámetros r. La iteración del gradiente incremental es

(

)

~ ~ r := r + γ∇ J (i, r ) c(i ) − J (i, r )

(7-52)

~

para todo i que pertenece a S , donde γ cumple con las condiciones de (7-47). J. A. Pucheta (www.labimac.blogspot.com)

57

Control óptimo y procesos estocásticos

Luego se calculan los costos asociados a cada par estado acción, mediante ~ ~ Q(i, u ) = L(i, u ) + J ( j, r ).

(7-53)

Entonces, se obtiene la política mejorada mediante

(

)

~ µ = arg min L(i, m ) + J ( j, r ) ∀i ) µ∈U (i~ para cada estado correspondiente al conjunto S .

(7-54)

Una vez que se dispone de J µ (i, r ) , se obtiene µ a partir de la ecuación (7-54). Luego se evalúan los costos asociados a cada estado, simbolizados por c(i), mediante la Ec. (7-52) se ajustan los ~ parámetros de r obteniendo una nueva versión de la función de aproximación J µ (i, r ) . Luego, se procede a efectuar la de mejora de política, en la que se obtiene una nueva política de control expresada como µ . Se comienza con el cálculo de los costos para cada estado, y entre cada iteración se efectúa la actualización de la función γ. ~

Consideraciones Dada la política estacionaria µ, los valores de Jµ no se calcularán exactamente, porque se tiene ~ una aproximación J µ (i, r ) , donde r es el parámetro de ajuste. Las fuentes de error, en la iteración de política aproximada son principalmente dos: ~ — La estructura de J µ (i, r ) puede que no sea lo suficientemente potente, por ejemplo si se disponen de pocos parámetros de ajuste. — El adecuado ajuste de los parámetros r cuando el algoritmo de sintonía de r no está bien ajustado o es deficiente. La política inicial del método requiere que sea tan buena como sea posible, caso contrario deberá fijarse el vector de parámetros r.

En el programa para Matlab PDA_J_It_Politica_731.m se implementó el algoritmo descrito. El algoritmo usado para resolver la Ec. (7-52) es el de Levenberg-Marquardt. En la Fig. 20 se muestran los resultados obtenidos. La función γ utilizada es γn =

J. A. Pucheta (www.labimac.blogspot.com)

5 5 + 1,5n

.

58

Control óptimo y procesos estocásticos

Estados

Costo

3

4

2.5 3 Costos

Estados

2 1.5

2

1 1 0.5 0

0

1

2

3

4

5

0

0

1

2

Acción de control

4

5

Evolucion del J(0) 2.5

1 Acciones de control

3

0.5

2

0

1.5

-0.5

1

-1 0

1

2

3

4

5

0.5

0

Etapas

5

10 Iteraciones

15

20

Fig. 20. Desempeño del controlador para el Ejemplo 7.3.1 mediante iteración de política aproximada.

7.8.4. Función compacta para generar el control La implementación de cualquier esquema de control con realimentación de estados, consiste en obtener la acción de control a aplicar al sistema partiendo del valor del vector de estados, se puede escribir como x k +1 = f ( x k , u k , k )

u k = µ(x k )

(7-55)

donde f (x k , u k , k ) es el modelo del proceso y µ (i ) es la función que representa la política óptima de decisiones a partir de los estados. La función µ (i ) es la solución analítica del planteo de control óptimo, que mediante la técnica de

~ (i, v ) de la función µ , donde v iteración de política aproximada se obtiene una aproximación µ (i ) ~ es el vector de parámetros de ajuste. Para encontrar la función aproximadora µ(i, v ) , se lleva a cabo la metodología que se detalla a continuación.

Se calculan las µ (i ) dentro del conjunto Sˆ mediante ~ (i, v ) − µ (i ))2 min ∑ (µ

(7-56)

v i∈S

~ (i, v ) , siendo v el vector de parámetros de ajuste. donde la ley de decisiones se representa por µ J. A. Pucheta (www.labimac.blogspot.com)

59

Control óptimo y procesos estocásticos

Una solución la Ec (7-56) se obtiene mediante el método del gradiente incremental. Cada iteración es representada por ~ (i, v )(µ (i ) − µ ~ (i, v )) v := v + γ∇µ

(7-57)

para todo i que pertenece a Sˆ , donde γ cumple con las condiciones de (7-47). Se resuelven dos problemas de aproximación al mismo tiempo: ~ -Dado µ, se simula y evalúa para encontrar J de Jµ. ~ -Dado J , calcular la política mejorada µ para algunos de los estados y luego encontrar la nueva ~ (⋅, v ) . política µ Aproximación de

i

u

política µ~(i, r ) Fig. 21. Expresión

compacta de la (función política) o ley de control.

~ (⋅, v ) , se obtienen las acciones de control como se indica en la Una vez disponible la función µ Fig. 21. El esquema de control es el que se muestra en la Fig. 22, es conocido como neurocontrolador debido a que una función de aproximación implementada mediante una red de neuronas genera las acciones de control.

f (x , u , k )

u Aproximación de política µ~(i, v ) Fig. 22. Expresión

i

compacta de la función política como controlador

En el programa para Matlab PDA_J_It_Politica_mu_731.m está implementado el algoritmo descrito para el Ejemplo 7.3.1. La función algoritmo usado para resolver la Ec.(7-52) y la Ec. (7-57) es el de Levenberg-Marquardt. En la Fig. 23 se muestran los resultados obtenidos. La función γ utilizada es γn =

J. A. Pucheta (www.labimac.blogspot.com)

100 . 100 + 1,5n

60

Control óptimo y procesos estocásticos

Estados

Costo

3

4

2.5 3 Costos

Estados

2 1.5

2

1 1 0.5 0

0

1

2

3

4

0

5

0

1

Acción de control

2

3

4

5

80

100

Evolucion del J(0)

Acciones de control

1

0.5

0

-0.5 0

10 -1 0

1

2

3 Etapas

4

5

0

20

40 60 Iteraciones

Fig. 23. Evolución del sistema del Ejemplo 7.3.1 cuando se controla con un neurocontrolador.

7.9. Discusión y comentario final Se han visto diferentes métodos para implementar un controlador que resuelva el problema de control óptimo formulado mediante las Ec. (7-33), (7-34), (7-35), y (7-36); también conocido como el problema básico del control óptimo. Cuando la implementación del controlador pretende hacerse mediante una forma compacta, al estilo de la Ec. (7-55), debe tenerse presente que para la adecuada implementación de la técnica descripta se requiere de gran conocimiento previo del proceso a controlar. Como sugerencias generales, se pueden mencionar: - El algoritmo es fuertemente dependiente de las condiciones iniciales, es decir, de la política inicial y de los estados usados empleados para calcular el neurocontrolador representados por los ~ conjuntos S y Sˆ . - La velocidad de ajuste de los parámetros se fija mediante la función γ, y el método es muy sensible a éste parámetro. Es común hacer los primeros intentos dejando γ=1 constante, con pocas iteraciones, y luego comenzar a modificarla para que se desvanezca con las iteraciones, siempre verificando que el desempeño del controlador mejore a largo plazo. - La cantidad de parámetros de ajuste en cada función de aproximación depende de la complejidad de los datos, en general se acondicionan implementando alguna técnica como normalización o extracción de características.

J. A. Pucheta (www.labimac.blogspot.com)

61

Control óptimo y procesos estocásticos

8. CONTROL DIGITAL ESTOCÁSTICO En el tratamiento convencional de los sistemas de control se considera que las señales actuantes sobre los mismos y las perturbaciones a las que están sometidos son modeladas con expresiones matemáticas determinísticas. Cuando se aplican criterios de optimización, los resultados que se obtienen dependen de las señales utilizadas. Eventualmente un sistema de control diseñado en forma óptima para señales determinísticas, podrá ser subóptimo cuando el mismo está sometido a perturbaciones reales. En la mayoría de los casos los resultados obtenidos del diseño determinístico son aplicables a sistemas reales, afectados por señales reales. No obstante, cuando se aumentan las exigencias de control, los controladores deben diseñarse no sólo basándose en la dinámica del sistema a controlar, sino también teniendo en cuenta las características estocásticas de las señales actuantes. 8.1. Modelo matemático estocástico de señales reales. Un PROCESO ESTOCÁSTICO ES UNA ABSTRACCIÓN MATEMÁTICA DE un proceso empírico cuyo desarrollo está gobernado por leyes probabilísticas. El proceso empírico se estudia como un modelo probabilístico que evoluciona en el tiempo y genera secuencias de valores numéricos. Experimento

Ley probabilística Evento A

P(A)

P(B)

Evento B

Espacio muestral Ω Conjunto de resultados

Fig. 8-1. Componentes del modelo probabilístico de un proceso incierto.

8.1.1. Modelo probabilístico de un experimento Los elementos del modelo probabilístico se pueden observar en la Fig. 8-1. Se simboliza mediante la terna (Ω, F, P) donde Ω es el conjunto de todas las posibles salidas del experimento, F es una colección de subconjuntos de Ω denominado eventos, que incluye a Ω, con las propiedades siguientes: (1) si A es un evento, entonces Ac={ω∈Ω | ω∉A} es también un evento, y se incluye Ωc ya que es el conjunto vacío ∅. J. A. Pucheta (www.labimac.blogspot.com)

62

Control óptimo y procesos estocásticos

(2) si A1, A2, ...Ak,... son eventos, entonces U∞k=1 A k y I∞k=1 A k son también eventos. P es una función que asigna a cada evento A un número no negativo P(A) denominado probabilidad del evento A que cumple con (1) P(Ω)=1. ∞ (2) P(U∞k =1 A k ) = ∑k =1 P(A k ) para cada secuencia mutuamente disjunta A1, A2, ...Ak,... P es la medida de probabilidad. 8.1.2. Variable aleatoria Se define sobre el espacio probabilístico (Ω, F, P) la variable aleatoria x, si para cada escalar λ el conjunto {ω∈Ω | x(ω)≤λ}∈F siendo así, un evento. Además, tendrá una función de distribución de probabilidad F(z)=P{ω∈Ω | x(ω)≤z}.

8.1.3. Proceso estocástico Un proceso estocástico en tiempo contínuo es una colección de variables aleatorias, definida mediante, {x(ω,t) | t∈ℜ}. Paralelamente, un proceso estocástico en tiempo discreto, se define mediante {x(ω,t) | t∈Z}.

8.2. Ecuaciones diferenciales estocásticas La definición de Ecuación diferencial estocástica es análoga a la de ED ordinaria, salvo que la función solución es un Proceso estocástico. La motivación para estudiar a ésta clase de ED consiste en resolver un problema real con determinada precisión. Por ejemplo, sea la EDO lineal de un sistema de primer orden sin excitación externa, dX (8-1) = a (t ) ⋅ X , X (0 ) = X 0 cte. dt Si a(t) es una función determinística escalar variante en el tiempo, X también lo será. Pero, si a(t) se describe mediante a (t ) = rt + ruido

(8-2)

entonces X será un proceso estocástico, y la Ec. (8-1) será una Ecuación diferencial estocástica (EDE). Para hallar una solución de la Ec. (8-1), se debe modelar a la señal “ruido”, es decir, proponer características que permitan limitar el alcance del planteo. Suponiendo que “ruido” sea un Proceso estocástico W ruido blanco con las siguientes propiedades: (i) ∀t1≠t2 entonces Wt1 y Wt2 son independientes. (ii) {Wt} es estacionario, es decir que la distribución conjunta de {Wt1+t,···,Wtk+t} no depende de t. (iii) E[Wt]=0 ∀t. Con ésta definición de W, se puede escribir a la Ec. (8-1) como

J. A. Pucheta (www.labimac.blogspot.com)

63

Control óptimo y procesos estocásticos

(8-3) dX = rt ⋅ X + σ t ⋅ W t ⋅ X dt donde σt es un escalar función del tiempo. Con mayor generalidad, las funciones temporales rt y σt pueden ser función de t y de X. Así, la Ec. (8-33) se convierte en dX (8-4) = r (t k , X k ) ⋅ X + σ (t k , X k ) ⋅ W t ⋅ X dt Resolviendo, se define el intervalo de integración [0, t] y en el mismo una partición 0=t0
(8-6)

donde Xj=X(tj), Wk=W(tk), ∆tk=tk+1-tk. Reemplazando el producto

Wk ⋅ ∆t k por

∆Vk = V(t k +1 ) − V(t k ) donde {Vt,t≥0} sugiere que sea un movimiento Browniano. Se demuestra porque es el único PE con trayectorias continuas y tales características en sus incrementos. Entonces, poniendo al PE Vt=Bt, se tiene que k −1

(

)

k −1

(

)

X k = X 0 + ∑ r t j , X j ⋅ X j ⋅ ∆t j + ∑ σ t j , X j ⋅ X j ⋅ ∆B j . j= 0

j= 0

(8-7)

Haciendo en la Ec. (8-7) ∆tk→0, si es que existe límite en algún sentido, se puede escribir a la solución X como t

t

X t = X 0 + ∫ r (s , X s ) ⋅ X s ⋅ ds + ∫ σ (s , X s ) ⋅ X s ⋅ dB s 0

(8-8)

0

donde dBt es un mB que inicia en el origen. Ahora la Ec. (8-38) puede escribirse incluso en la forma compacta diferencial dX t = r (t , X t ) ⋅ X t ⋅ dt + σ (t , X t ) ⋅ X t ⋅ dB t , X (0 ) = X 0 cte.

(8-9)

que es la versión estocástica general de la Ec. (8-1). 8.2.1. Solución de la ecuación diferencial estocástica Para hallar la solución de la EDE (8-38), una de las herramientas más poderosas es la fórmula de Itô. En ésta sección sólo se desarrollará la aplicación, como motivación del estudio de las EDEs y su utilidad para representar procesos reales. Suponiendo que r,σ∈R, en la Ec. (8-38), se tiene que dX t = r ⋅ dt ⋅ X t + σ ⋅ X t ⋅ dB t

(8-10)

donde se pueden agrupar las variables de la forma dX t = r ⋅ dt + σ ⋅ dB t Xt integrando ambos miembros, se tiene J. A. Pucheta (www.labimac.blogspot.com)

(8-11)

64

Control óptimo y procesos estocásticos t

t

t

dX ∫ X ss ⋅ ds = ∫ r ⋅ ds + ∫ σ ⋅ dB s = r ⋅ t + σ ⋅ B t . 0 0 0

(8-12)

Para hallar la solución del lado izquierdo, se evaluará mediante la fórmula de Itô la función ln(Xt). Por lo tanto, se tiene que dX t 1 1 1  1 (8-13) d (ln X t ) = dX t +  − 2  ⋅ (dX t )2 = − ⋅ (dX t )2 2 Xt 2  Xt  Xt 2X t Reemplazando el dXt de la Ec. (8-10) d (ln X t ) =

dX t 1 − ⋅ (r ⋅ dt ⋅ X t + σ ⋅ X t ⋅ dB t )2 2 Xt 2X t

(8-14)

de donde finalmente, usando las reglas de derivación de Itô, se obtiene d (ln X t ) =

dX t σ 2 − ⋅ dt . Xt 2

(8-15)

Integrando ambos miembros respecto del tiempo, se tiene dX s σ2 ∫ d (ln X s ) ⋅ ds = ∫ X s ⋅ ds − ∫ 2 ⋅ ds. 0 0 0 t

t

t

(8-16)

donde el segundo término del lado derecho es la igualdad de la Ec. (8-12) (8-17) Xt 1 = r ⋅ t + σ ⋅ B − σ2 ⋅ t X0 2 despejando Xt se tiene la expresión del proceso estocástico solución de la EDE de la Ec. (8-10), ln

La expectativa es

  1  X t = X 0 ⋅ exp   r − σ 2  ⋅ t + σ ⋅ B . 2   

(8-18)

E [X t ] = E [X 0 ] ⋅ exp (r ⋅ t ).

(8-19)

Nótese que la Ec. (8-32) tiene igual estructura que las soluciones a EDOs lineales determinísticas. 8.3. Modelos de Estado para Sistemas Estocásticos de Tiempo continuo Sea el sistema RLC representado en la Fig. 1-1, con la carga del circuito modelada mediante la EDO d 2Q dQ 1 dQ L 2 +R + Q = Ve , Q (0 ) = Q 0 , = I0 (8-20) dt C dt dt cuyas condiciones iniciales corresponden a la corriente de la bobina y la carga del capacitor. Si se considera que la tensión de entrada ve tiene una componente de ruido, la EDO de la Ec. (8-20) se convierte en una EDE. Más aún, si se pretende medir la carga q y se considera un ruido en la medición, el problema se transforma en un problema de estimación de q, cuya metodología se conoce como filtro Kalman-Bucy. Es un procedimiento que consiste en estimar el estado de un sistema que cumple una ODE ruidosa basada en series de mediciones ruidosas. Entonces, expresando a la tensión de entrada como una función determinística más una componente de ruido, se tiene ve = G + α ⋅ W (8-21) J. A. Pucheta (www.labimac.blogspot.com)

65

Control óptimo y procesos estocásticos

y haciendo la asignación  x1   Q  X =   =  dQ .  x 2   dt 

(8-22)

Reemplazando en la Ec. (8-20) se obtiene  dx 1  dt = x 2  dx 1  L 2 = − Rx 2 − x 1 + G + α ⋅ W C  dt

(8-23)

y en notación matricial, se puede escribir dX = A ⋅ X ⋅ dt + H ⋅ dt + K ⋅ dB t

 0  dx1   1 , A = donde Bt es un mB unidimensional, dX =   dx − LC  2

(8-24)

1   0  0 R , H =  1 , K =  α . −   L G   L  L

Reescribiendo la Ec. (8-24) para agrupar a X, dX − A ⋅ X ⋅ dt = H ⋅ dt + K ⋅ dB t

(8-25)

y pre multiplicando por exp(-A·t) ambos miembros, exp (− A ⋅ t ) ⋅ dX − exp (− A ⋅ t ) ⋅ A ⋅ X ⋅ dt = exp (− A ⋅ t ) ⋅ [H ⋅ dt + K ⋅ dB t ]

(8-26)

donde se intentará relacional al lado izquierdo con d(exp(-A·t)·X). Para ello, se usa la versión multidimensional de Itô. Se definen dos funciones g1,g2 tal que g:[0,∞)×R2→R2, dada por x  g (t , x 1 , x 2 ) = exp (− A ⋅ t ) ⋅  1  x 2 

(8-27)

Aplicando Itô multidimensional, se tiene que d (exp (− A ⋅ t ) ⋅ dX ) = (− A ) exp (− A ⋅ t ) ⋅ dX ⋅ dt + exp (− A ⋅ t ) ⋅ dX .

(8-28)

Sustituyendo en la Ec. (8-26) se tiene que t

t

0

o

exp (− A ⋅ t ) ⋅ X − X 0 = ∫ exp (− A ⋅ s ) ⋅ H ⋅ ds + ∫ exp (− A ⋅ s ) ⋅ K ⋅ dB s

(8-29)

que integrando por partes siguiendo el método Itô X = exp (A ⋅ t ) ⋅ [X 0 + exp (− A ⋅ t ) ⋅ K ⋅ B t +

(8-30)

t  + ∫ exp (− A ⋅ s ) ⋅ [H + A ⋅ K ⋅ B s ] ⋅ ds   0 Nótese la complejidad de las soluciones que expresan matemáticamente el estado del sistema en el dominio del tiempo continuo, cuando se desea encontrar una solución exacta para todo t∈[0, ∞). Sin embargo, es una metodología muy poderosa que merece mencionarse.

J. A. Pucheta (www.labimac.blogspot.com)

66

Control óptimo y procesos estocásticos

8.3.1. Procesos Estocásticos de Tiempo Discreto. Siguiendo con la definición anterior, dado un experimento aleatorio cuyos elementos θ pertenecen a un espacio muestral Ω, se puede definir una función ξ que asocia a cada elemento de dicho espacio una función de tiempo discreto xi(k), donde k es la notación simplificada de KT, siendo T el período de muestreo. Al conjunto imagen χ que contiene todas las funciones vinculadas al experimento aleatorio se lo define como Proceso Aleatorio o Estocástico. Cada función es una realización del proceso. Esto se puede representar como (8-31)

ξ:Ω → χ

con χ = {x i (k, θ), i = 1, 2, ...}.

(8-32)

De esta manera para un instante determinado n, queda definida una variable aleatoria xn. Si en lugar del tiempo se fija un elemento cualquiera θj del experimento aleatorio queda definida una función. Las propiedades estadísticas del proceso estocástico se pueden describir por su función de densidad de probabilidad y por las funciones de probabilidad conjunta de todas las variables aleatorias que se pueden definir en el proceso. Cuando estas caracterizaciones estadísticas son funciones del tiempo, entonces se dice que el proceso es variante. Por el contrario si son independientes del instante de observación se dice que es invariante. Esto implica que los parámetros estadísticos son iguales para todas las variables aleatorias que se puedan definir en el proceso. Cuando en un conjunto de señales aleatorias, que conforman un proceso estocástico estacionario, cualquiera de las realizaciones del proceso se puede obtener por desplazamiento temporal de alguna de ellas, el mismo se denomina Ergódico. En este caso los parámetros estadísticos obtenidos para una variable aleatoria, en un instante definido, sobre todas las funciones o realizaciones, son iguales a los que se obtienen a lo largo del tiempo de evolución de una de las señales. 8.3.2. Parámetros estadísticos de procesos estocásticos. La caracterización estadística de un proceso estocástico se realiza teniendo en cuenta las variables aleatorias que se pueden definir en el proceso. Si el proceso es ergódico se cumple lim

N→∞

1 N

N

∑ k =1

(8-33)

1 M ∑ Mi xi (k) M →∞ M i =1

x(k) = lim

donde N es el número total de muestras de una realización, Mi es la cantidad de veces que tiene lugar el suceso xi(k) y M el número total de funciones o realizaciones. El primer miembro de la Ec. (8-33) se denomina valor medio temporal. El segundo representa el valor medio de las muestras para un tiempo dado k. El poder establecer esta igualdad es de gran importancia ya que permite aplicar todas las herramientas matemáticas de la estadística a las señales temporales. El valor medio, expectativa o esperanza matemática se expresa como m x = E { x(k) } = lim

N →∞

J. A. Pucheta (www.labimac.blogspot.com)

1 N

(8-34)

N

∑ x(k). k =1

67

Control óptimo y procesos estocásticos

Los parámetros estadísticos como valor medio, varianza, dispersión y momentos de distinto orden, son valores que dan información sobre las señales estocásticas, pero que son independientes de la variable tiempo. Esto es, no tiene información sobre cuan rápido cambia la señal o cual es su contenido armónico. Una función que provee información sobre la coherencia interna de una señal o sobre su velocidad de cambio es la Función de Autocorrelación, definida como la esperanza matemática del producto de la señal consigo misma desplazada en el tiempo, m intervalos de muestreo. Φxx (m) = E{x(k) x(k + m)}. (8-35) En el dominio discreto se tiene: 1 N (8-36) Φ xx (m) = E {x(k) x(k + m)} = lim ∑ x(k) x(k + m). N →∞ N k =1 Se observa de la definición que el valor medio de la señal influye en la correlación. Si se desea obtener una función que solamente tenga en cuenta las fluctuaciones respecto del valor medio, se pueden substraer éstos de la señal. En este caso se obtiene la Función de Autocovarianza o simplemente Covarianza, definida como: R xx (m) = cov {x, m} = E {[x(k) - m x ] [x(k+ m) - m x ]} = Φ xx (m) - m x .

(8-37)

2

Para m=0, se obtiene la varianza de la señal 1 N 2 2 σx = R xx (0) = lim ∑ [x(k) - m x ] . N →∞ N k =1

(8-38)

Si la señal tiene distribución gaussiana, las características estadísticas estáticas y dinámicas quedan totalmente definidas por la esperanza y la función de autocovarianza. Si estos valores son independientes del tiempo, la señal es estacionaria en sentido amplio. La dependencia estadística entre dos señales estocásticas estacionarias x(k) e y(k) está definida por la Función de Intercorrelación. En tiempo discreto está definida por 1 N (m) = E { x(k) y(k + m) } = lim x(k) y(k + m). Φ xy ∑ (8-39) N→∞ N k =1 De igual modo, cuando se eliminan los valores medios de las señales se obtiene la denominada Función de Intercovarianza. R xy (m) = cov {x, y, m } = E {[x(k) - m x ][y(k+ m) - m y ]} = Φ xy (m) - m x m y .

(8-40)

Dos variables aleatorias x e y se denominan incorreladas si

R xy (m) = 0 .

(8-41)

Si además es mx my = 0 , se dicen ortogonales y son independientes si p(x, y) = p(x) p(y)

(8-42)

donde p(x), p(y) son las funciones de densidad de probabilidad de x, y respectivamente, mientras que p(x, y) es la función de densidad de probabilidad conjunta. Un ruido blanco es una señal para la cual cualquier valor es independiente de todos los valores pasados. Esto implica que la función de covarianza será idénticamente nula para cualquier m, excepto para m=0 donde tomará un valor definido. Si el ruido blanco tiene una distribución gaussiana, el mismo quedará totalmente definido por su valor medio mx y la función de covarianza 2 Rxx (m) = σx (m) J. A. Pucheta (www.labimac.blogspot.com)

(8-43) 68

Control óptimo y procesos estocásticos

siendo 1 para m = 0 2 σx (m) =  0 para m ≠ 0. Cuando una señal está definida por un conjunto de n valores en cada instante, se tiene una señal estocástica vectorial o vector señal de orden n definida como x T (k) = [ x1 (k)

x 2 (k) ...

x n (k)] .

(8-44) Este vector siempre se podrá descomponer en sus n señales escalares. Si cada señal es estacionaria su valor medio estará definido por T T m x = E {x (k)} = [ m x 1 m x 2 ... m xn ] .

(8-45) La coherencia interna del vector señal se define mediante el valor esperado del producto del vector señal, sin su valor medio, con su traspuesto desplazado en el tiempo. Esto queda expresado por la Función Matricial de Covarianza. T R xx (m) = E { [ x(k) - m x ] [ x(k + m) - m x ] } o bien

 R x x (m) R x x (m)   R x x (m) R x x (m)   . . R xx (m) =  . .   . .  R (m) R (m) xx  xx 1

1

1

2

2

1

2

2

n

1

n

K R x x (m)   K R x x (m)   K .  K .  K . K R x x (m) 

2

1

n

2

n

n

(8-46)

n

Esta función matricial queda formada por todas las funciones de covarianza de las distintos componentes del vector señal de orden n. Sobre la diagonal se encuentran las n funciones de autocovarianza de las n componentes de la señal vectorial y todos los otros elementos son las funciones de intercovarianza entre las componentes de la señal. Cuando las componentes del vector señal son señales escalares incorreladas, su intercovarianza será cero y la función matricial será una matriz diagonal.

0  R x x (m)  0 R x x (m)    . 0 R xx (m) =  . .   . .   0 0  1

K

1

2

K

2

K K K K

0  0  .  .  0 R x x (m) n

(8-47)

n

8.3.3. Procesos Gaussianos. Una variable aleatoria x se dice que es gaussiana cuando está caracterizada por la función de densidad de probabilidad

 1  (x - m )  2  x  p(x) = exp −    1/2  2  σx  (2π ) σ x   1

J. A. Pucheta (www.labimac.blogspot.com)

(8-48)

69

Control óptimo y procesos estocásticos

propuesta por Gauss y de allí su nombre. Se puede hacer una extensión para n variables aleatorias, organizadas en forma vectorial, obteniéndose  1  −1 / 2 T -1 (8-49) p[x (k )] = 2 π-n/2 [det R xx (k )] exp − [x (k ) - m x (k )] R xx (k ) [x (k ) - m x (k )] 2   donde m x (k ) = E{x (k )} (8-50)

 R x1x1 (k ) R x1x 2 (k )   R x 2 x1 (k ) R x 2 x 2 (k )   . . R xx (k ) =  . .   . .    R x n x1 (k ) R x n x 2 (k )

K K K K K K

R x1x n (k )  R x 2 x n (k )  .  .  .  R x n x n (k )

(8-51)

R x i x j (k ) = E{[xi (k ) - m xi (k )] ⋅ [x j (k ) − m xj (k )]}.

(8-52)

De esta manera, un proceso aleatorio estacionario caracterizado por la Ec. (8-49), se dice que es gaussiano. 8.3.4. Proceso de Markov. Cuando la probabilidad condicional de un evento de valor x(k) en un proceso estocástico depende solamente de su último valor x(k-1) y de ninguno anterior, se denomina Proceso de Markov de primer orden. Esto se expresa mediante la función de densidad de probabilidad condicionada tal que p {x(k) / [x(k -1), x(k -2), ..., x(0) ] }= p {x(k) / x(k -1)} . (8-53) Se puede decir entonces que en un proceso de Markov, el próximo valor de la señal x(k+1) estará influenciado solamente por el valor presente x(k). En una ecuación escalar en diferencias de primer orden: x (k + 1) = a ⋅ x(k ) + b ⋅ v(k ) (8-54) el valor futuro depende de los valores presente de x(k) y v(k). En el caso que v(k) sea una señal estadísticamente independiente, como un ruido no correlacionado, la función de densidad de probabilidad de x(k+1) dependerá sólo de x(k) y será por lo tanto un proceso de Markov. Si la ecuación en diferencias es de mayor orden como x (k + 1) = a n x (k ) + a n -1 x (k − 1) + L + a1 x (k − n − 1) + b ⋅ v(k) (8-55) siempre se puede transformar la ecuación que describe el proceso, en varias ecuaciones de primer orden mediante una apropiada elección de variables auxiliares, tal que x(k ) = x1 (k ) x(k + 1) = x1 (k + 1) = x 2 (k ) (8-56) .

. .

xn -1 (k + 1) = x n (k ). J. A. Pucheta (www.labimac.blogspot.com)

70

Control óptimo y procesos estocásticos

El proceso podrá expresarse como una ecuación matricial en diferencias de primer orden 1 0 K 0  x1 (k + 1) 0      x 2 (k + 1)  0 0 1 K 0     . . K . .  .   = . . . K . .    .  0 0 0 K 1    x n (k + 1) a1 a 2 . K a n 

 x1 (k )  0       x 2 (k )  .     .  .    +  . v(k ). .     .  .    x n (k ) b   

(8-57)

En forma sintética x (k + 1) = A x (k ) + b v(k ).

(8-58)

Si la matriz A y el vector b son constantes, cada elemento de x(k+1) depende sólo del valor presente de x(k) y es por lo tanto un proceso de Markov. En la Fig. 8-2 se representa la Ec. (8-58).

Fig. 8-2. Modelo de Markov para x(k).

8.4. Modelos de Estado para Sistemas Estocásticos de Tiempo Discreto. 8.4.1. Extensión del modelo de estado determinístico para sistemas estocásticos. Frecuentemente los sistemas se ven expuestos a perturbaciones que no se pueden describir analíticamente y por lo tanto no es posible conocer su evolución futura como en el caso de los sistemas determinísticos estudiados. En esta circunstancia, las señales aleatorias de perturbación deberán ser tratadas estadísticamente. Por otra parte, el estado también se verá afectado aleatoriamente y por lo tanto tendrá una función de distribución de probabilidad. Una extensión natural del concepto de estado para sistemas determinísticos a los sistemas estocásticos, es exigir que la función de distribución de probabilidad del vector de estado, en el futuro, sólo dependa del valor actual del estado. Una forma de transformar la ecuación de estado determinística en una ecuación estocástica, es suponer que x(k+1) es una variable aleatoria que depende de x(k) y de e(k). Es decir x (k + 1) = f [x (k ), u (k ), k ] + e[x (k ), k ] (8-59) donde f es el valor medio de la función de distribución de probabilidad condicional de x y e es una secuencia de vectores aleatorios con distribución gaussiana con valor medio me [x (k ), k ] = E{e[x (k ), k ]} = 0 (8-60) y matriz de covarianza T R ee[x(k - m), k - m] = E{e[x(k ), k ] e [x(m), m]}.

(8-61)

Si las componentes de e son estadísticamente independientes, entonces la matriz de covarianza Ree, para m=0, toma la forma 2 2 R ee[x(k ), k ] =diag σe1 [x (k ), k ] ,L, σen [x(k), k ] J. A. Pucheta (www.labimac.blogspot.com)

{

}

71

Control óptimo y procesos estocásticos

(8-62) donde σ ei es la varianza de la componente i-ésima del vector aleatorio e y proporciona una medida de la potencia de la señal aleatoria o de ruido que afecta el sistema. 2

El vector e se puede modelar de manera que provenga de otro vector de ruido que tenga valor medio nulo y matriz de covarianza unitaria, esto es e[x(k), k ] = F[x (k ), k] ⋅ v[x(k ), k ]

se cumple con

(8-63)

m v [x (k ), k ] = E{v[x (k ), k ]} = 0

(8-64)

T R vv [x (k ), k ] = E{v[x (k ), k ] ⋅ v [x(k ), k ]} = I

(8-65)

reemplazando la Ec. (8-63) en la Ec. (8-60) y la Ec. (8-61)

me [x (k ), k ] = F[x (k ), k ] E{v[x(k), k ]}

(8-66)

T R ee[x(k), k ] = F[x(k ), k ] E{v[x(k ), k ] v [x(k ), k ]}FT [x(k), k ].

(8-67)

Comparando las Ecs. (8-66) y (8-67) con las Ecs. (8-64) y (8-65) se deduce que me [x(k), k ] = F[x (k ), k ] m v [x (k ), k ] T R ee [x (k ), k ] = F[x (k ), k ] R vv [x (k ), k ] F [x (k ), k ].

(8-68) (8-69)

De esta forma se obtienen los parámetros estadísticos del vector e en función de los del vector v. Reemplazando los valores definidos para éstos, las ecuaciones quedan me [x (k ), k ] = 0

(8-70)

R ee [x(k), k ] = F [x (k ), k ] FT [x (k ), k ].

(8-71)

Por lo tanto conocida la matriz de covarianza del vector aleatorio que afecta al sistema, éste se puede modelar de manera que el vector aleatorio de entrada tenga componentes estadísticas con distribución gaussiana normal (0, 1). La matriz F se obtiene descomponiendo Ree, para el caso particular dado por Ec. (8-62) tomando la forma F [x (k ), k ] =diag {σe1 [x (k ), k ], σe2 [x(k), k ], L , σen [x (k ), k ]}. (8-72) Para que la función de distribución de probabilidad condicional x(k+1) dado x(k), sea estadísticamente independiente de los valores pasados de x, la función de distribución condicional de e[x(k), k] dado x(k), tampoco debe depender del pasado de x. Sólo cumpliendo estas condiciones la Ec. (8-59) recibe el nombre de ecuación matricial-vectorial estocástica en diferencias y es un proceso aleatorio de Markov. Si la función vectorial f en la Ec. (8-59) es lineal en x y en u, y el vector de ruido e no depende de x, la Ec. (8-59) se puede escribir como x (k + 1) = A(k ) ⋅ x (k ) + B(k ) u (k ) + F(k ) v(k ) con parámetros estadísticos m v = E{v(k)} = 0

{

(8-73)

}

T R vv = E v(k ) ⋅ v (k ) = I

para el vector e será

J. A. Pucheta (www.labimac.blogspot.com)

me (k ) = F(k ) m v (k ) = 0 T T R ee(k ) = F(k )R vv(k ) F (k ) = F(k )F (k ). 72

Control óptimo y procesos estocásticos

De igual modo que para el vector de estado, la ecuación determinística para la salida, se puede transformar en una ecuación estocástica. El vector de salida estocástico estará dado por y(k ) = g[x (k ), u (k ), k ] + z[y(k ), k ] (8-74) donde z es un vector aleatorio independiente de e, que afecta a la salida y cuyos parámetros estadísticos son mz [y(k ), k ] = E{z[y(k ), k ]} = 0

(8-75)

T R zz [y (k ), k ] = E {z [y (k ), k ] z [y (k ), k ]}.

(8-76)

De igual manera, como se trató al vector aleatorio de estado e, el vector aleatorio de salida z se puede modelar dando como resultado. z[y(k ), k ] = G[y(k ), k ] w [y(k ), k ] donde w es un vector aleatorio con

(8-77)

m w [y(k ), k ] = E{w [y(k ), k ]} = 0 T R ww[y(k ), k ] = E w[y(k ), k ] w [y(k ), k ] = I de manera que para z se obtiene

{

(8-78)

}

(8-79)

mz [z(k ), k ] = G[y(k ), k ]m w [y(k ), k ]

(8-80)

T R zz[y(k ), k ] = G[y(k ), k ]R ww[y(k ), k ]G [y(k ), k ]

(8-81)

y reemplazando los valores dados en las Ecs. (8-78) y (8-79) en las Ecs. (8-80) y (8-81) mz [y(k ), k ] = 0 T R zz[y(k ), k ] = G[y(k), k ]G [y(k ), k ]. Si las componentes de z son estadísticamente independientes, la matriz Rzz tiene la forma

[

(8-82) (8-83)

]

2 2 R zz[y(k ), k ] = diag σz1[y(k), k ],L, σzn [y(k), k ]

(8-84)

y por lo tanto G

G[y(k ), k ] = diag [σz1 [y(k ), k ], L , σzn [y(k ), k ]] si G es lineal en x y en u, y z es independiente de y, la Ec. (8-74) se puede escribir

(8-85)

y(k ) = C(k ) x (k ) + D(k ) u (k ) + G (k ) w (k ). Con parámetros estadísticos m w (k ) = E{w (k )} = 0

(8-86)

T R ww (k ) = E{w (k ) ⋅ w (k )} = I

para el vector z

m z (k ) = G (k ) m w (k ) = 0 T T R zz (k ) = G (k ) R ww (k ) G (k ) = G (k ) G (k ). En la Fig. 8-3 se muestra una representación en diagrama de bloques de las Ecs. (8-73) y (8-86).

J. A. Pucheta (www.labimac.blogspot.com)

73

Control óptimo y procesos estocásticos

Fig. 8-3. Diagrama de bloques del sistema lineal estocástico de tiempo discreto.

Si el sistema es invariante y los vectores aleatorios son estacionarios, las Ecs. (8-73) y (8-86) se transforman en x (k + 1) = Ax (k ) + B u (k ) + F v(k ) (8-87) y(k ) = C x (k ) + D u (k ) + G w (k )

con

(8-88)

m v = 0; R vv = I m w = 0; R ww = I.

8.4.2. Solución de la ecuación en diferencias de estado para sistemas estocásticos. La solución de una ecuación en diferencias para sistemas determinísticos, consiste en obtener el valor de estado x en todo instante. En forma semejante resolver la ecuación en diferencias para sistemas estocásticos significa conocer la distribución de probabilidad conjunta de las variables de estado, en todo instante. Suponiendo conocida la distribución de probabilidad conjunta de x(0), x(1) (condiciones iniciales); se puede calcular la distribución condicional de x(1) dado x(0); puesto que x(1) está dado como una función de x(0) y v(0). Si la función de distribución de x(0) es conocida, la distribución conjunta se puede encontrar aplicando el teorema de Bayes. p [x (k ), x (k - 1), L , x (0 )] = p [x(k) / x(k -1) ]L p [x (1) / x(0) ] p[x (0 )]

(8-89)

Admitiendo que el vector de estado inicial x(0) tiene una distribución gaussiana, entonces x es un proceso gaussiano, puesto que se forma por combinación lineal de variables aleatorias que tienen esa distribución; lo mismo sucede con el vector de salida y formado por la combinación lineal de variables aleatorias con la misma distribución. Entonces los procesos estocásticos x, y quedan completamente caracterizados por su función de valor medio y su función de covarianza. Para la función de valor medio se tiene m x (k + 1) = E{x (k + 1)}

(8-90)

reemplazando la Ec. (8-73) en la Ec. (8-90) m x (k + 1) = E{A (k )x (k ) + B(k ) u (k ) + F(k ) v(k )} = = A(k ) E{x (k )} + B(k ) E{u (k )} + F(k )E {v(k )}

(8-91)

m x (k + 1) = A (k ) m x (k ) + B(k ) mu (k ) + F(k ) m v (k ).

J. A. Pucheta (www.labimac.blogspot.com)

74

Control óptimo y procesos estocásticos

Si el sistema es invariante y el vector v es estacionario m x (k + 1) = A m x (k ) + B mu (k ) + F mv (k ) cuya solución es k -1

m x (k ) = A m x (0 ) + ∑ A k

k -i -1

(8-92)

B m u (i ) + F m v (k )

(8-93)

i=0

con valor inicial y para el vector de salida

m x (0) = E{x (0 )} = m0

(8-94)

m y (k ) = E {y (k )}

(8-95)

Reemplazando la Ec.(8-86) en la Ec. (8-95) m y (k ) = C(k ) E{x (k )} + D(k ) E{u (k )} + G (k )E{w (k )} , y por lo tanto queda

m y (k ) = C (k ) m x (k ) + D (k ) m u (k ) + G (k ) m w (k ) Si el sistema es invariante y el proceso estacionario m y (k) = C m x (k) + D m u (k) + G m w (k)

(8-96) (8-97)

k -1

m y (k ) = C A k m x (0 ) + D m u (k ) + G m w (k ) + ∑ C A k -i -1 B m u (i ) + CF m v (k )

(8-98)

i=0

Las Ecs. (8-91), (8-92), (8-93), (8-96), (8-97), y (8-98) son análogas a las ecuaciones correspondientes a sistemas determinísticos 23, 24, 28 y 30. Cuando el valor medio de los vectores aleatorios v y w son nulos quedan mx (k + 1) = A(k ) mx (k ) + B(k ) mu (k ). Si el sistema es invariante

(8-99)

m x (k + 1) = A m x (k ) + B mu (k ) k -1

(8-100)

i=0

(8-101)

k k -i -1 B m u (i ). m x (k ) = A m x (0 ) + ∑ A

De igual modo

m y (k ) = C (k ) m x (k ) + D (k ) m u (k ).

Si el sistema es invariante

(8-102)

m y (k ) = C m x (k ) + D m u (k )

(8-103)

k -1

m y (k ) = C A k m x (0 ) + D m u (k ) + ∑ C A k -i -1 B m u (i ).

(8-104)

i=0

En la Fig. 8-4 se muestra una representación del modelo, en diagrama de bloques de las Ecs. (8-91) y (8-96), referidas a los valores medios, el cual es análogo al mostrado en la Fig. 8-3, referido al modelo de estado estocástico.

J. A. Pucheta (www.labimac.blogspot.com)

75

Control óptimo y procesos estocásticos

Fig. 8-4. Diagrama de bloques que representa a las Ecs. (8-91) y (8-96).

Para la función de covarianza del vector de estado se tiene

{

}

T R xx (k + 1) = E [x (k + 1) - mx (k + 1)] [x (k + 1) - mx (k + 1)] .

Reordenando

(8-105)

R xx (k + 1) = E{x (k + 1) x T (k + 1)}- m x (k + 1) mTx (k + 1).

(8-106)

reemplazando la Ec. (8-73) en la Ec. (8-106) y considerando además que v es estadísticamente independiente de u y x; y que mv(k)=0 se obtiene T T R xx (k + 1) = A(k ) R xx (k ) A (k ) + A(k ) R xu (k ) B (k ) + + B(k ) R ux (k ) AT (k ) + B(k )R uu (k ) BT (k ) +

(8-107)

+ F(k ) R vv(k ) FT (k ) - m x (k + 1) mTx (k + 1)

para sistemas invariantes y vector aleatorio v estacionario se tiene

R xx (k + 1) = A R xx (k ) AT + A R xu (k ) BT + B R ux (k ) AT + + B R uu (k ) BT + F R vv FT

(8-108)

con valor inicial

T T R xx (0 ) = E{x (0 ) x (0 )}- m0 m0 = R 0 . Suponiendo que el valor de las entradas es cero, lo que implica R xu (k ) = R ux (k ) =R uu (k ) = 0 , y utilizando la Ec. 5.4-13, la Ec. 5.4-49 se reduce a

con el valor inicial

R xx (k + 1) = A(k ) R xx (k ) AT (k ) +R ee(k ).

(8-109)

R xx (0 ) = E{x (0 ) x T (0 )} = R 0. . Si el sistema es invariante y el vector aleatorio e es estacionario, la Ec. (8-109) se transforma en R xx (k + 1) = A R xx (k ) AT + R ee. (8-110) Desarrollando para k=0, 1, 2, …

J. A. Pucheta (www.labimac.blogspot.com)

76

Control óptimo y procesos estocásticos T T R xx (1) = A R xx (0) A + R ee = A R 0 A + R ee k = 1, R xx (2) = A R xx (1) AT + R ee = AA R xx (0) AT AT + A R ee AT + R ee

k = 0,

k = 2,

T T T T R xx (3) = A R xx (2) A + R ee = AAA R xx (0) A A A + + AA R ee AT AT + A R ee AT + R ee . .

.

.

para k genérico k -1

∑ Ak-i-1 R ee (AT )k-i-1 .

T k

k

R xx (k) = A R xx (0) (A ) +

(8-111)

i=0

Si todos los autovalores de la matriz A tienen módulo estrictamente menor que 1 entonces la serie Ec. (8-110) converge a un valor finito R xx = lim R xx (k) k →∞

por lo tanto la Ec. (8-110) toma la forma R xx = A R xx AT + R ee .

(8-112)

La Ec. (8-110) se denomina ecuación matricial en diferencias de Liapunov y la Ec. (8-112) ecuación matricial de Liapunov en estado estacionario o ecuación algebraica de Liapunov. De esta manera se observa que conociendo la matriz A del sistema y la matriz de covarianza del vector aleatorio que afecta al sistema, se puede determinar el valor medio y la matriz de covarianza del estado, en cualquier instante. Con ello queda completamente definida la caracterización estadística del mismo.

{

Para el vector de salida se define T R yy(k) = E [y(k) - m y (k)] [y(k) - m y (k)]

{

T

}

= E y(k) y (k) - m y (k)

}

(8-113)

mTy (k)

reemplazando las Ecs. (8-86) y (8-96) en la Ec. (8-113)5.4-55 R yy(k + 1) = C(k) R xx (k) CT (k) + C(k) R xu (k) DT (k) + + D(k) R ux (k) CT (k) + D(k) R uu (k) DT (k) + + G(k) R ww (k) G T (k) - m y (k) - mTy (k) con valor inicial

(8-114)

R xx (0) = E{x(0) x T (0)}- m0 mT0 = R 0.

Para sistemas invariantes y vector aleatorio w estacionario se tiene T T T R yy(k + 1) = C R xx (k) C + C R xu (k) D + D R ux (k) C +

(8-115)

+ D R uu (k) DT + G R ww G T .

Suponiendo que el valor de las entradas es cero, lo que implica R xu (k) = R ux (k) = R uu (k) = 0 , J. A. Pucheta (www.labimac.blogspot.com)

77

Control óptimo y procesos estocásticos

y utilizando la Ec. (8-83), la Ec. (8-115) queda T R yy(k) = C R xx (k) C + R zz

(8-116)

reemplazando la Ec. 5. (8-111) en la Ec. (8-116) k R yy(k) = C Ak R xx (0) (AT ) CT +

k -1

∑ C Ai

i R ee (AT ) CT + R zz

(8-117)

i =0

si los valores de A tienen módulo estrictamente menor que 1 entonces la Ec. (8-116) converge a un valor finito T T R yy = C A R xx AT C + C R ee C + R zz . (8-118) La Ec. 5.4-60 constituye la matriz de covarianza de la salida, en estado estacionario. Las Ecs. (8-107), (8-108), (8-111), (8-114), (8-115) y (8-117) son análogas en estructura a las Ecs. 23, 24, 28, 30 obtenidas para los sistemas determinísticos. 8.5. Diseño de Controladores de Estado para Sistemas Estocásticos Lineales. 8.5.1. Formulación del problema de control para sistemas estocásticos Si se toma el funcional de costos N −1

J (x, u ) = x T (N )Sx (N ) + ∑ x T (k )Qx (k ) + u T (k )Ru (k ) k =0

para obtener una ley de control óptimo, cuando se trata de un sistema estocástico lineal de tiempo discreto, modelado por las ecuaciones

x (k + 1) = A(k ) x (k ) + B(k) u(k) + F(k ) v(k )

(8-119)

y(k ) = C(k ) x (k ) + D(k ) u (k ) + G (k ) w (k )

(8-120)

donde v(k), w(k) son vectores aleatorios con E{v(k )} = E{w (k )} = 0

E{v(k ) vT (k )} = E{w (k ) w T (k )} = I.

E{v(k ) vT (n )} = E{w (k ) w T (n )} = 0 ∀ n ≠ k

E{v(k ) w T (k )} = E{x (k ) vT (k )} = E{v(k ) w T (k )} = 0 para todo valor de k; y con parámetro estadístico para el estado inicial, E{[x (0)]} = m0

{

}

E [x (0) - m0] ⋅ [x (0) - m0]T = R 0

se presentará el inconveniente, que como el vector de estado es aleatorio, entonces el funcional también lo es, y por lo tanto no hay forma de establecer el valor mínimo del mismo. Por este motivo se define un nuevo funcional de costo como algún parámetro estadístico, por ejemplo el valor medio, J. A. Pucheta (www.labimac.blogspot.com)

78

Control óptimo y procesos estocásticos N -1  T  { ( ) } ( ) ( ) = E J x , u = E N Sx N + Jm x ∑ xT (k ) Q x (k ) + uT (k ) R u (k ) k =0  

(8-121)

con S y Q simétricas y semidefinidas positivas y R simétrica y definida positiva. De esta manera el problema del Regulador Óptimo Lineal Estocástico (LQGR) se puede formular como sigue. Dado un sistema estocástico lineal de tiempo discreto modelado en el espacio de estado con Ecs. (8-119) y (8-120); se debe encontrar una ley de control u, que modifique el estado llevándolo desde un valor inicial x(0) a uno final x(N)=0 y que simultáneamente haga mínimo el funcional de costos de Ec. (8-121). 8.5.2. Solución al problema del regulador óptimo lineal estocástico. Suponiendo que se ha medido el vector de estado y se tiene x(0), x(1), …, x(k) se debe determinar el vector de entrada u(k). Puesto que la Ec. (8-119) es una ecuación estocástica matricial-vectorial en diferencias, la función de distribución de probabilidad condicional de los valores futuros del estado, dado los valores pasados del mismo, sólo dependen del valor actual x(k). Por lo tanto es suficiente obtener u(k), como función de x(k) y no de los valores pasados x(k-1), …, x(0). Siguiendo la misma metodología que para el caso determinístico, el funcional de Ec. (8-121) se puede descomponer en dos partes de la siguiente manera N -1   J m = E{J(x , u )} = E x T (N ) S x (N ) + ∑ x T (i )Qx (i ) + u T (i ) R u (i ) + (8-122) i=k    k -1  + E ∑ x T (i ) Q x (i ) + u T (i ) R u (i )  i =0  donde el segundo término no es necesario minimizar puesto que no depende de u(k), u(k+1), …, u(N-1). Entonces se obtiene min J m = min E{J (x, u )} = min E{x T (N ) S x (N ) + u u u (8-123) N -1  T T + ∑ x (i ) Q x (i ) + u (i ) R u (i ) i= k  o Si se indica con u la función vectorial que hace mínimo el funcional de costo; y se considera que esta tiene un mínimo único se puede escribir

E{J( x, u )} ≥ E{J(x, u°)} = min[E{J (x, u )}].

(8-124)

u

Por otra parte también se cumple

(

)

J (x , u ) ≥ J x, u o = min J (x, u ).

(8-125)

u

Tomando esperanza matemática en todos los miembros

{(

)}

{ [ ( )]}

E{J (x , u )} ≥ E J x, u o = E min J x, u o . J. A. Pucheta (www.labimac.blogspot.com)

(8-126)

u

79

Control óptimo y procesos estocásticos

Minimizando esta expresión respecto de u y teniendo en cuenta que el segundo y tercer término son constantes se tiene

)}

{

u

Comparando las Ecs. (8-124) y (8-127) se tiene

{

u

{(

}

(8-127)

}

(8-128)

min [E{J (x, u )}] ≥ E J x, u o = E min [J (x, u )] . u

min[E{J (x, u )}] = E min[J (x, u )] . u

Esta última expresión indica que cuando el índice de desempeño tiene un mínimo único, las operaciones de minimización y de cálculo de esperanza matemática son conmutables. Entonces la Ec. (6-2) toma la forma

{

}

min J m = E min[J (x, u )] = u

u

k -1 (8-129)     = E min  x T (N ) S x (N ) + ∑ x T (i ) Q x (i ) + u T (i ) R u (i ) .  u  i=0   Es decir que, minimizar el funcional de costo estocástico dado por Ec. (6-1) es equivalente a minimizar el funcional determinístico y luego calcular la esperanza matemática. De esta manera se ha trasladado el problema de minimización de un funcional estocástico al caso de minimización de un funcional determinístico. Por lo tanto, aplicando el principio de optimización de Bellman y tomando como punto de partida la ecuación

N -1

J (x , u ) = x (N ) S x (N ) + ∑ x T (i ) Q x (i ) + u T (i ) R u (i )

(8-130)

T

i=k

ésta se puede descomponer en

J (x k , u k ) = x T (N ) S x (N ) + x T (N - 1) Q x (N - 1) + u T (N - 1) R u (N - 1) +

(8-131)

N -2

+ ∑ x (i ) Q x (i ) + u (i ) R u (i ). T

T

i=k

Definiendo al funcional de costos en el instante k como función del funcional en el instante k+1,

J k = J (x k , u k ) = x (k )T Q x (k ) + u T (k ) R u (k ) + J k +1

(8-132)

Para k+1=N se tiene

J N −1 = x (N - 1)T Q x (N - 1) + u T (N - 1) R u (N - 1) + x T (N ) S x (N )

(8-133)

donde JN = xT(N) S x(N), y el estado x(N) es

x (N ) = A(N - 1) x (N - 1) + B(N - 1) u (N - 1) + F(N - 1) v(N - 1)

(8-134)

Reemplazando la Ec. (6-5) en la Ec. (6-3) y haciendo min J (x , u ) = J x, u o ,

(

)

(Tener en cuenta que, dada una matriz X y dos vectores x, y se verifica

∂ ( x T X y) ∂ ( x T X y) ∂(xT X x) =X y ; = XT x ; = 2X x ∂x ∂y ∂x J. A. Pucheta (www.labimac.blogspot.com)

80

Control óptimo y procesos estocásticos

(8-135)

la tercera propiedad sólo es válida si X es simétrica) se tiene ∂J N −1 = 2 BT (N - 1) S [A(N - 1) x(N - 1) + F(N - 1) v(N - 1)] + ∂u(N - 1)

(8-136)

+ 2 B (N - 1) S B(N - 1) u (N - 1) + 2 R u (N - 1) = 0 T

Despejando u se tiene

u°(N - 1) = - [R + BT (N - 1) S B(N - 1)] BT (N - 1) S [A(N - 1) x (N - 1) + + F(N - 1) v(N - 1)]. -1

(8-137)

reemplazando el estado y la acción de control para k=N-1, que son las Ecs. (6-5) y (6-7) en la Ec. (6-4) y reordenando min J N-1 = [A(N - 1) x (N - 1) + F(N - 1) v(N - 1)] H(N - 1) ⋅ T

u ( N −1)

(8-138)

⋅ [A(N - 1) x (N - 1) + F(N - 1) v(N - 1)] + x T (N - 1) S x (N - 1),

con

[

]−1 BT (N - 1) S

H(N - 1) = S - S B(N - 1) R + BT (N - 1) S B(N - 1)

(8-139)

agrupando los términos en x(N-1) y v(N-1), definiendo P(N - 1) = Q + AT (N - 1) H(N - 1) A(N - 1)

(8-140)

y aplicando esperanza matemática en la Ec. (6-8) se obtiene E  min J N -1  = E {x T (N ) S x (N )}+ E {v T (N - 1) F T (N - 1) H (N - 1) F (N - 1) v (N - 1)}+  u ( N −1 )  + 2 E {x T (N - 1) A T (N - 1) H (N - 1) F (N - 1) v (N - 1)}

(8-141)

Reemplazando la Ec. (6-9) en la Ec. (6-10) se obtiene P(N - 1) = Q + AT (N - 1) S A(N - 1) - AT (N - 1) S B(N - 1)

[R + B (N - 1) S B(N - 1)] T

-1

(8-142)

B (N - 1) S A(N - 1). T

Por otra parte, se puede demostrar que dada las formas bilineal y cuadrática zTXy y zTXz respectivamente, con matriz X simétrica; la esperanza matemática es

J. A. Pucheta (www.labimac.blogspot.com)

E{zT X y}= mTz X m y + tr[X R zy ]

(8-143)

E{zT X z}= mTz X mz + t r[X R zz ]

(8-144)

81

Control óptimo y procesos estocásticos

Aplicando estas propiedades al segundo y tercer término de la Ec. (8-141) y teniendo en cuenta que Rvv(k)=I y mv(k)=0 se obtiene E  min J N -1  = E{x T (N - 1) P (N - 1) x (N - 1)}+ u ( N -1)  + t r FT (N - 1) H (N - 1) F(N - 1)

[

(8-145)

]

Se procede a continuación a minimizar para u(N-2), se aplica entonces la propiedad expresada en Ec. (8-128) a la Ec. (6-11)   T min E x T (N - 2 )Qx (N - 2 ) + u (N - 2 ) Ru (N - 2 ) + min J N -1  = u ( N − 1 )  

u ( N − 2 ) 

  T E  min  x T (N - 2 )Qx (N - 2 ) + u (N - 2 ) Ru (N - 2 ) + min J( x, u )  = u ( N − 2 ) u ( N − 1 )   

[

]

[

]

E  min J d N − 2  + tr F T (N - 1) H (N - 1) F(N - 1) u ( N − 2 )  donde

J d N-2 = x T (N - 1) P(N - 1) x (N - 1) + x T (N - 2) Q x (N - 2) + u T ( N - 2 ) R u (N - 2 )

(8-146)

con

x (N - 1) = A(N - 2) x (N - 2) + B(N - 2) u (N - 2) + F(N - 2) v(N - 2)

(8-147)

como las Ecs. (6-12) y (6-13) son análogas a las Ecs. (6-4) y (6-5), se opera repitiendo los pasos desde Ec. (6-6) hasta Ec. (6-11) obteniéndose

[

]-1 BT (N - 2) P(N - 1).

u °(N - 2 ) = - R + BT (N - 2 ) P(N - 1) B(N - 2 ) .[A(N - 2 ) x (N - 2 ) + F(N - 2 ) v(N - 2 )]

(8-148)

H(N - 2) = P(N - 1) - P(N - 1) B(N - 2).

[

] -1 BT (N - 2) P(N - 1)

. R + BT (N - 2) P(N - 1) B(N - 2)

(8-149)

P(N - 2) = Q + AT (N - 2) H(N - 2) A(N - 2)

(8-150)

P(N - 2 ) = Q + AT (N - 2 ) P(N - 2 ) A(N - 2 ) - AT (N - 2 ) P(N - 1) B(N - 2 ) .

[

(8-151)

]−1 BT (N - 2) P(N - 1) A(N - 2)

. R + BT (N - 2 ) P(N - 1) B(N - 2 )

[

]

E  min J d N-2 = E{x T (N - 2 ) P(N - 2 ) x (N - 2 )}+ tr F T (N - 2 ) H(N - 2 ) F(N - 2 ) + (8-152) u ( N − 2 ) 

[

]

+ tr F T (N - 1) H(N - 1) F(N - 1)

Continuando con la inducción para u(N-3, u(N-4, ···, u(k+1), u(k)

[

]-1

u °(k) = - R + BT (k ) P(k + 1) B(k )

(8-153)

B(k ) P(k + 1) [A(k ) x (k ) + F(k ) v(k )]

[

]−1 BT (k ) P(k + 1)

H(k ) = P(k + 1) - P(k + 1) B(k ) R + BT (k ) P(k + 1) B(k ) J. A. Pucheta (www.labimac.blogspot.com)

P(k + 1) = Q + AT (k ) H(k + 1)A(k )

P(k ) = Q + A (k )P(k + 1) A(k ) - A (k ) P(k + 1)B(k ) T

T

[R + T (k ) P(k + 1) B(k )]−1 T (k ) P(k + 1) A(k ).

(8-154) 82

Control óptimo y procesos estocásticos

(8-155)

N -1  T  ( ) ( ) ( ) = E k P k x k + Jm x ∑ tr FT (i ) H(i ) F(i )  i=k  

[

]

(8-156)

min J k = m Tx (k ) P(k ) m x (k ) + tr P(k ) R xx (k ) +

(8-157)

uk

N -1

+ ∑ tr F T (i ) H(i ) F(i ). la Ec. (6-20) s la ecuación matricial de Riccati y se debe resolver para obtener el vector de control. i=k Comparando las Ecs. (8-142), (6-17) y (6-20) se deduce que la condición inicial para P es

P (N ) = S .

(8-158)

De la Ec. (6-18) se concluye que la ley para obtener el vector de entrada óptimo, está formado por dos partes; la primera es una realimentación lineal del vector de estado tal como se obtiene para el regulador óptimo lineal determinístico y además una ley de prealimentación del vector aleatorio de perturbación. En la Fig. 8-5 se observa un diagrama de bloques del sistema con controlador. Es importante destacar que se necesita conocer el vector de perturbación aleatorio; en caso que esto no sea posible, se debe hacer una predicción del mismo. Cuando la predicción no es factible se obtiene una ley incompleta con un incremento de Jm. La Ec. (6-18) se puede escribir u °(k ) = -K x (k ) x (k ) - K v (k ) v(k )

[

(8-159)

]-1

T T K x (k ) = R + B (k ) P(k + 1) B(k ) B (k ) P(k + 1) A(k ).

[

]

K v (k ) = R + B T (k ) P(k + 1) B(k )

-1

(8-160)

B T (k ) P(k + 1) F(k )

(8-161)

la Ec. (8-160) es la misma que para el problema determinístico (6-23). Si el sistema es invariante, entonces se puede obtener P en estado estacionario, esto se cumple haciendo P = lim P(k) (8-162) k →∞ y por lo tanto se obtiene

[

]

-1

P = Q + A T PA - A T PB R + B T P B B T P A

(8-163)

que es la ecuación matricial de Riccati en estado estacionario. Entonces la ecuación del controlador será

[

]

u °(k ) = - R + B T P B B T P[A x (k ) + F v(k )].

J. A. Pucheta (www.labimac.blogspot.com)

-1

(8-164)

83

Control óptimo y procesos estocásticos

Si F no se conoce y se elimina de la Ec. (6-26), se obtiene el controlador determinístico. Esto significa que dicho controlador determinístico es también la solución óptima obtenible para el caso estocástico cuando no se puede modelar el ruido.

Fig. 8-5. Sistema con ley óptima de realimentación y prealimentación del vector de perturbación.

Para la Fig. 8-5 se ha tomado

donde

u°(k ) = - K (k )[A(k ) x (k ) + F(k ) v(k )]

[

(8-165)

]

K (k ) = R + B T (k ) P(k + 1) B(k ) B T (k ) P(k + 1). -1

(8-166)

Por último si se calcula la derivada segunda en la Ec. (6-6) se tiene

[

]

∂ 2 J( x, u ) = 2 R + B T (N - 1) S B(N - 1) (8-167) 2 ∂ u (N - 1) de manera que para que el funcional de costo sea mínimo, se debe asegurar que S o R sean definidas positivas. Si se asegura a R como definida positiva, S puede ser semidefinida positiva, pero tanto S como Q nunca pueden ser definidas negativas porque entonces el índice Jm puede resultar negativo, lo que contradice la definición de cualquier funcional de costo. 8.5.3. Controladores de Estado de Varianza Mínima. En la sección anterior se estudió el problema del regulador óptimo estocástico, como una extensión del regulador óptimo lineal determinístico. Se presenta a continuación otro criterio de diseño basado en el concepto de varianza mínima. Dado un sistema dinámico estocástico lineal de tiempo discreto con ecuaciones x(k + 1) = A(k) x(k) + B(k) u(k) + F(k) v(k)

(8-168)

y(k) = C(k) x(k) + D(k) u(k) + G(k) w(k)

(8-169)

donde v(k), w(k), x(k) son vectores aleatorios con distribución gaussiana con parámetros estadísticos J. A. Pucheta (www.labimac.blogspot.com)

84

Control óptimo y procesos estocásticos

E{v(k)} = E{w(k)} = 0

(8-170)

E{v(k) w T (k)}= E{v(k) x T (k)} = 0

(8-171)

E{v(k) vT (k)}= E{w(k) w T (k)} = I

(8-172)

E{v(k) vT (n)}= E{w(k) w T (n)} = 0, ∀n ≠ k.

(8-173)

Se define el funcional de costo como la varianza del estado J m (k + 1) = E{J(x, u )} = E{x 2 (k + 1)}.

(8-174)

Para que el vector de control no resulte con amplitudes demasiado grandes conviene incorporarlo en el índice de manera que J m (k + 1) = E { J( x, u ) } = E { x 2 (k + 1) + r u 2 (k)} (8-175) puesto que x, u son vectores J m (k + 1) = E{J( x, u )} = E{x T (k + 1) x(k + 1) + r u T (k) u(k)}.

(8-176)

Por último es importante disponer de la posibilidad de modificar cuanto influye cada componente de x y de u sobre el índice, por ello se incluyen dos matrices de peso R, Q, donde Q es simétrica y semidefinida positiva y R es simétrica y definida positiva. J m (k + 1) = E{J(x, u )} = E{x T (k + 1)Qx(k + 1) + u T (k)Ru(k)}

(8-177)

donde el factor r se ha incluido en la matriz R. Se puede hacer entonces la formulación del problema. Dado un sistema dinámico estocástico lineal de tiempo discreto modelado en el espacio de estado con Ecs. (8-168) y (8-169) se debe encontrar una ley de control u que minimice el funcional de costo dado por la Ec. (8-177) o lo que es lo mismo se minimiza la varianza del vector de estado. Para resolver el problema planteado se reemplaza la Ec. (8-168) en la Ec. (8-177); obteniéndose T J m (k + 1) = E { [A(k) x(k) + B(k) u(k) + F(k) v(k) ] Q [ A(k) x(k) + B(k) u(k) +

+ F(k) v(k) ] + u T (k) R u(k) } T J m ( x, u) = [ A(k) x(k) + B(k) u(k) + F(k) v(k) ] Q [ A(k) x(k) + B(k) u(k) + + F(k) v(k) ] + u T (k) R u(k)

(8-178)

operando en la Ec. (8-177) y minimizando con respecto a u se tiene

{

}

min J m (k + 1) = E min J( x, u ) = E{J m (x, u°)}. u

(8-179)

donde uo se obtiene de J. A. Pucheta (www.labimac.blogspot.com)

85

Control óptimo y procesos estocásticos

∂ J m ( x, u ) =0 ∂ u(k) teniendo en cuenta que, dada una matriz A y dos vectores x, y es ∂ ( xT A y) ∂ ( xT A y) ∂ ( xT A x) =A y ; = A T x; = 2Ax ∂ x ∂ y ∂x se obtiene

∂ J m (u, x ) = 2 BT (k)Q[A(k) x(k) + F(k) v(k)] + ∂ u(k)

(8-180)

+ 2 BT (k) Q B(k) u(k) + 2 R u(k) = 0

despejando u se tiene

[

u °(k) = - R + 2 BT (k) Q B(k)

]-1

T B (k) Q[A(k) x(k) + F(k) v(k)]

(8-181)

llamando u °(k) = - K x (k) x(k) - K v (k) v(k)

(8-182)

[

]-1

T B (k) Q A(k)

(8-183)

[

]-1

T B (k) Q F(k).

(8-184)

K x (k) = R + 2 BT (k) Q B(k)

K v (k) = R + 2 BT (k) Q B(k)

Las Ecs. (8-182), (8-183) y (8-184) son equivalentes en estructura que las Ecs. (8-159), (6-23) y (8-161) obtenidas para el controlador óptimo lineal estocástico. En este caso la matriz P que satisface la ecuación matricial de Riccati es sustituida por la matriz de peso Q. Para el caso que el sistema sea lineal e invariante y el vector aleatorio v estacionario, la ecuación del controlador será -1 u °(k) = - R + 2 BT Q B BT Q[A x(k) + F v(k)]. (8-185)

[

]

Para la Fig. 8-6 se ha tomado

u°(k) = K(k)[A(k) x(k) + F(k) v(k)]

(8-186)

donde

[

K(k) = R + 2 BT (k) Q B(k)

]-1

T B (k) Q.

(8-187)

De esta manera la ley de control óptima, en el sentido de la obtención de una varianza mínima, posee dos partes, una ley de realimentación de estado similar al caso determinístico y una ley por prealimentación del vector aleatorio de ruido que afecta al sistema.

J. A. Pucheta (www.labimac.blogspot.com)

86

Control óptimo y procesos estocásticos

Es importante destacar que se necesita conocer el vector de perturbación aleatorio; en caso que esto no sea posible, el mismo se debe predecir. Cuando la predicción no es factible se obtiene una ley incompleta con un incremento del índice J. La Ec. (8-185) se puede obtener de la Ec. (6-18) si se reemplaza P(k+1) por Q. Por este motivo la estructura de los controladores son iguales como se puede observar comparando la Fig. 8-5 y la Fig. 8-6.

Fig. 8-6. Sistema de control con ley óptima de realimentación de estado y ley de prealimentación del vector de perturbación para índice de mínima varianza.

9. Bibliografía Anderson B., Moore J., 1971. Linear optimal control. Prentice-Hall International Inc., London. Bellman R. and S. Dreyfus, 1962. Applied dynamic programming. Princenton University Press. Bertsekas D. and J. Tsitsiklis, 1996. Neuro-dynamic programming. Athena scientific. MIT. Bertsekas D. and J. Tsitsiklis. Notas del curso Introduction to Probability. MIT. Kirk, Donald E.. Optimal Control Theory: An Introduction, Dover Publications, 2004. Luus, R. Iterative Dynamic Programming.CRC Press Ogata, Katsuhiko. Ingeniería de Control Moderna. Prentice Hall. Ogata, Katsuhiko. Sistemas de Control en Tiempo Discreto 2da Ed. Prentice Hall. Oksendal B., 2003. Stochastic differential equations 6ed., Springer.

J. A. Pucheta (www.labimac.blogspot.com)

87