Resumos
Neste trabalho apresentaremos uma técnica iterativa baseada em simulações de Monte Carlo para calcular o controle ótimo de um problema de regulador linear quadrático de horizonte infinito para um sistema linear com saltos Markovianos a tempo discreto, quando a matriz de transição de probabilidade não é conhecida. Sabemos que o controle ótimo deste problema é dado em termos da solução maximal de um conjunto de equações algébricas de Riccati acopladas entre si (EARA) a tempo discreto, que foram extensivamente estudadas nos últimos anos. Traçaremos um paralelo com a teoria do algoritmo TD(lambda) para Processos Markovianos de Decisão (PMD) para desenvolver o algoritmo TD(lambda) para o controle ótimo associado à solução maximal de uma EARA.
Simulações de monte carlo; equações algébricas de Riccati acopladas entre si; sistemas com saltos; controle ótimo
In this paper we present an iterative technique based on Monte Carlo simulations for deriving the optimal control of the infinite horizon linear regulator problem of discrete-time Markovian jump linear systems for the case in which the transition probability matrix of the Markov chain is not known. It is well known that the optimal control of this problem is given in terms of the maximal solution of a set of coupled algebraic Riccati equations (CARE), which have been extensively studied over the last few years. We trace a parallel with the theory of TD(lambda) algorithms for Markovian decision processes to develop a TD(lambda) like algorithm for the optimal control associated to the maximal solution of the CARE. Some numerical examples are also presented.
Monte carlo simulations; coupled algebraic riccati equations; jump systems; optimal control
Método de diferenças temporais aplicado às equações de Riccati acopladas entre si
Oswaldo L. V. Costa; Julio C.C. Aya
Departamento de Engenharia de Telecomunicações e Controle, Escola Politécnica da Universidade de São Paulo, CEP: 05508 900 São Paulo SP Brazil, oswaldo@lac.usp.br, julio@lac.usp.br
RESUMO
Neste trabalho apresentaremos uma técnica iterativa baseada em simulações de Monte Carlo para calcular o controle ótimo de um problema de regulador linear quadrático de horizonte infinito para um sistema linear com saltos Markovianos a tempo discreto, quando a matriz de transição de probabilidade não é conhecida. Sabemos que o controle ótimo deste problema é dado em termos da solução maximal de um conjunto de equações algébricas de Riccati acopladas entre si (EARA) a tempo discreto, que foram extensivamente estudadas nos últimos anos. Traçaremos um paralelo com a teoria do algoritmo TD(l) para Processos Markovianos de Decisão (PMD) para desenvolver o algoritmo TD(l) para o controle ótimo associado à solução maximal de uma EARA.
Palavras-chave: Simulações de monte carlo, equações algébricas de Riccati acopladas entre si, sistemas com saltos, controle ótimo.
ABSTRACT
In this paper we present an iterative technique based on Monte Carlo simulations for deriving the optimal control of the infinite horizon linear regulator problem of discrete-time Markovian jump linear systems for the case in which the transition probability matrix of the Markov chain is not known. It is well known that the optimal control of this problem is given in terms of the maximal solution of a set of coupled algebraic Riccati equations (CARE), which have been extensively studied over the last few years. We trace a parallel with the theory of TD(l) algorithms for Markovian decision processes to develop a TD(l) like algorithm for the optimal control associated to the maximal solution of the CARE. Some numerical examples are also presented.
Keywords: Monte carlo simulations, coupled algebraic riccati equations, jump systems, optimal control.
1 INTRODUÇÃO
Neste trabalho consideraremos o seguinte modelo em um espaço probabilístico (W,P,{ k}, ) apropriado, conhecido na literatura internacional como sistemas lineares com saltos Markovianos a tempo discreto (vide (Mariton, 1990)):
onde q(k) é uma cadeia de Markov tomando valores em {1,...,N} com a matriz de transição de probabilidade = (pij). Seja S = {u = (u(0),...);u(k) é k-mensurável para cada k}. Para u Î S, considere o seguinte funcional quadrático para o sistema (1):
Deseja-se minimizar (2) sob u Î S . Nos artigos ((Costa e Fragoso,1995), (Ji e Chizeck,1988), (Ji et al., 1991)) mostra-se que a solução deste problema esta associada à existência de uma solução P = (P1,...,PN), Pi > 0 i = 1,...,N, do seguinte conjunto de equações algébricas de Riccati acopladas entre si (EARA), para i = 1,...,N
onde (X) = (1(X),..., N(X)) é definido, para X = (X1,...,XN) como
Se tal solução P existe, pode-se mostrar (vide (Costa e Fragoso, 1995)) que a lei de controle ótima para o problema dado pelas equações (1), (2) e (3) é obtida pela lei de controle de realimentação
onde F = (F1,...,FN) é dado por
Associado aos problemas de horizonte infinito (1), (2), temos o problema de horizonte finito (7) no qual, para algum S = (S1,...,SN) Î n+, deseja-se encontrar um u(0),...,u(k) tal que minimize o seguinte funcional
Como é provado em (Costa e Fragoso, 1995), a solução deste problema pode ser obtida do seguinte conjunto de equações a diferença de Riccati para k = k,...,0,
onde P(k + 1) = S e
A lei de controle ótima de (7) é dada por
u(k) = Fq (k)(k)x(k)
e sob algumas condições (vide (Costa e Fragoso, 1995)), P(k) converge para uma solução semi-definida positiva P de (3).
Traçando um paralelo entre a teoria acima e a teoria de processos Markovianos de decisão (PMD), poderíamos relacionar as iterações (8),(9) com a técnica chamada de iteração de valores em PMD. Por outro lado, como será apresentado no Lema 5 abaixo (vide Observação 1), o método chamado de quasi-linearização para obter a solução de (3) pode ser visto como a técnica de iteração de estratégias para PMD, envolvendo a parte de avaliação de estratégia e melhoramento de estratégia.
O método do TD( ) tem sido aplicado para resolver problemas relacionados a PMD (veja por exemplo (Bertesekas e Tsitsilklis, 1996),(Sutton e Barto, 1998)), para o caso onde a matriz de transição de probabilidade da cadeia de Markov é desconhecida. Neste caso assumimos que se conhece o histórico do processo, isto é, que existe um conjunto de trajetórias simuladas do sistema, e que a função de custo de uma dada estratégia é calculada iterativamente através do histórico do processo.
A meta deste trabalho é aplicar o método TD () para obter o controle ótimo (5), (6) associado à solução P do EARA dado pelas equações (3) e (4) para o caso em que a matriz de transição de probabilidade é desconhecida, sendo possível simular as trajetórias para a cadeia de Markov q(k). Para o caso em que a matriz de probabilidade é conhecida, existe uma serie de algoritmos numéricos para obter a solução maximal P do EARA (3) e (4) (veja ((Abou-Kandil et al., 1995), (Ait-Rami e Ghaoui, 1996), (Blair e Sworder, 1975), (Costa e Fragoso, 1995), (Costa e Marques, 1999), (Costa et al., 1997), (Gajic e Borno, 1995), (Ji e Chizeck, 1988), (Ji et al., 1991), (Mariton, 1990), (Val et al., 1998), (Costa e Aya, 1999)). Para o caso no qual a matriz de transição de probabilidade é não conhecida, propõe-se um algoritmo semelhante ao da avaliação de estratégias de Monte Carlo que atualiza as estimativas incrementalmente para obter o controle ótimo F dado pelas equações (5), (6). Obtém-se a prova deste convergência, seguindo os mesmos argumentos que em (Bertesekas e Tsitsilklis, 1996).
Este trabalho é apresentado do seguinte modo: na Seção 2 apresentamos algumas notações e resultados preliminares. Na seção 3 descrevemos o método TD() para resolver a EARA (3). O resultado principal é o Teorema 7, que prova a convergência do método apresentado nesta seção. Na Seção 4 mostraremos um exemplo numérico, e na Seção 5 os comentários finais.
2 NOTAÇÃO E RESULTADOS PRELIMINARES
Para e espaços de Banach, definimos (, ) o espaço de Banach de todos os operadores lineares de em , com a norma uniforme induzida representada por ||·||. Por simplicidade usaremos () := (, ). O raio espectral de um operador Î () será denotado por rs (). Se é um espaço de Hilbert o produto interno será denotado por á ·;·ñ, e para Î (), * denotará o operador adjunto de . Como é usual, > 0 ( > 0 respectivamente) denotará que o operador Î () é positivo semi-definido (positivo definido). Em particular denotaremos o espaço complexo n-dimensional por n e por (n, m) o espaço linear de todas as matrizes m × n com norma limitada, e (n) := ( n, n) e o produto interno em (n) é dado por áH;Vñ = tr{H*V} para H,V Î (n). O superescrito ' denotará a transposta de uma matriz. Usamos neste trabalho ||·||2 para a norma Euclidiana em n e a norma induzida em (n).
Neste trabalho, serão ainda empregados os seguintes resultados:
Lema 1 Sejaum espaço de Hilbert eÎ (). As seguintes afirmações são equivalentes:
a) rs () < 1
b) existe umÎ (n) invertível e Î () tal que || || < 1 e
=
-1.Prova: Corolário 1.14 de (Kubrusly, 1997), páginas 31-32.
Lema 2 Sejaum espaço de Hilbert eÎ (). Se rs () < 1 então para cada Q Î existe uma única solução S(Q) Î para o sistema em X
Além disso
Prova: Veja Teorema 5.17, página 102 em (Weidman, 1980).
O conjunto
n,m é formado pelo espaço linear de todas as N-seqüência de matrizes complexas V = (V1,...,VN) com Vi Î (n, m), i = 1,...,N, e por simplicidade n = n,n. Neste trabalho consideraremos o conjunto n equipado com o seguinte produto interno: para H = (H1,...,HN), V = (V1,...,VN) Î n definimos o produto interno á·;·ñ em n como segue:
e ||H||2: = áH;Hñ (de forma que ||H||2 = ||Hi||2). Então com o produto interno acima definido (vide equação (12)), n é um espaço de Hilbert. Nos também temos a seguinte norma ||· ||max em n. Para H = (H1,...,HN) Î n, ||H||max é definido por
||H ||max := max{||Hi||2 ;i = 1,...,N}.
Definimos
n+ := {V = (V1,...,VN) Î n; Vi > 0,"i}e escrevemos que V > S onde V = (V1,...,VN) Î n e S = (S1,...,SN) Î n, se
V - S = (V1 - S1,...,VN - SN) Î n+,
e que V > S se Vi - Si > 0 para i = 1,...,N. Para, G = (G1,...,GN) Î n definimos os seguintes operadores, (·) = (1(·),..., N(·)) Î (n) e (·) = (1(·),..., N(·)) Î (n) para V = (V1,...,VN) Î n e i,j = 1,...,N,
onde o operador (·) = (1(·),..., N(·)) Î (n) é definido como em (4). É simples verificar que os operadores , , e mapeiam n+ em n+, e que rs () = rs (). No artigo (Costa e Fragoso, 1993) se prova que rs () = rs () (de fato = * de acordo com a equação (12)). Também definimos o operador Î (n) como
onde, de novo, V = (V1,...,VN) Î n e (V) = (1(V),..., N(V)). No artigo de (Costa e Fragoso, 1993) provou-se que rs () = rs () = rs ().
Assumimos no modelo (1) e na função de custo (2) que A = (A1,...,AN) Î n, B = (B1,...,BN) Î m,n, C = (C1,...,CN) Î n,p e D = (D1,...,DN) Î m,p. Como for mostrado em (Costa e Fragoso, 1993), (Mariton, 1988), o modelo (1) com u(k) = Fq(k)x(k), e Vi(k) = E(x(k)x(k)*1 {q(k) = i}), V(k) = (V1(k),...,VN(k)) Î n+ leva a
V(k + 1) = (V(k)), k = 0,1,...
onde Gi = Ai + BiFi em (14), e E(|| x(k)||2) = tr{Vi(k)}.
O próximo passo é definir o operador para o 4essimo momento de x(t). Para isto usamos o produto de Kronecker L
K Î () para L,K Î (n) e o operador vec{·}: (n) definido da forma comum (veja (Brewer, 1978)). Seja o operador Î () da seguinte forma: para S = (S1,...,SN) Î , (S) = (1(S),..., N(S)) é definido como
onde Si(k) = E(vec{x(k) x(k)*}vec{x(k)x(k)*}*1 {q(k) = i}), S(k) = (S1(k),...,SN(k)) Î +. Temos o seguinte resultado:
Lema 3 S(k + 1) = (S(k)) e E(||x(k)||4) = tr{Si(k)}.
Prova: Como x(k + 1) = Gq(k)x(k), é fácil de se verificar que
xj(k + 1)xj(k + 1)* = 1{q(k+1) = j}Gixi(k)xi (k)*
onde xi(k) := x(k)1 {q(k) = i}. Seja zi(k) = vec{xi(k)xi(k)*}, i = 1,...,N. Depois de algumas manipulações segue que
zj(k + 1)zj(k + 1)* = 1{q(k+1)=j}
e escrevemos
Si(k) = E(zi(k)zi(k)* 1{q(k) = i}).
Conseqüentemente temos que
Sj(k + 1) = pij (iGi)Si(k)(
).Finalmente, note que
tr{zi(k)zi (k)*} = ||xi(k)||2tr{xi(k)xi(k)}
=||xi(k) ||4.
Vamos definir o conceito de estabilidade que será usado nas seções seguintes.
Definição 1 Diremos que F = (F1,...,FN) Î n,m estabiliza (A,B) no sentido da média quadrática se, quando fazemos u(k) = Fq(k)x(k) no sistema (1), temos que E(||x(k)|| 2) ® 0 quando k ® ¥ para qualquer condição inicial de x(0) e q(0). Diremos que (A,B) é estabilizável na média quadrática se para algum F = (F1,...,FN) Î n,m, temos que F estabiliza (A,B) no sentido da média quadrática.
O seguinte resultado, provado em (Costa e Fragoso, 1993), mostra que F = (F1,...,FN) estabiliza o sistema (1) no sentido da média quadrática se e somente se o raio espectral do operador (14) (ou (13)) em malha fechada é menor que 1.
Lema 4 F = (F1,...,FN) Î n,mestabiliza (A,B) no sentido da média quadrática se e somente se rs () < 1, onde é como em (14) com Gi = Ai + BiFi.
Definição 2 Definimos (·) = (1(·),..., N(·)): n+®n,m, (·) = (1(·),..., N(·)): n,m ®n e = (1(·),..., RN(·)):n+ ® n+ como
onde X = (X1,...,XN) Î n+, F = (F1,...,FN) Î n,m.
A seguinte identidade será útil nas próximas seções: para qualquer F = (F1,...,FN) Î n,m, temos que
O seguinte Lema provado em (Costa e Marques, 1999), provê a existência de uma solução maximal para (3) quando (A,B) é estável na média quadrática. É baseado na técnica de quasi-linearização para a EARA (3), e paralelamente à técnica de iteração de estratégias para PMD (vide Observação 1 abaixo).
Lema 5 Suponha que (A,B) é estabilizável na média quadrática e considere F0 = (,...,) Î n,m tal que estabiliza (A,B) no sentido da média quadrática. Então para l = 0,1,2,..., existe Xl = (,...,) que satisfaz às seguintes propriedades:
a) X0>X1>...>Xl>X, para um X Î n+ arbitrário tal que X >
(X).b) rs (l) < 1, onde l (·) = ((·),...,(·)) e para i = 1,...,N,
c) Xl satisfaz Xl = l (Xl) + (Fl) e é dado por Xl = (l)k ( (l)).
Além disso existe um X+ = (,...,) Î n+tal que X+ = (X+), X+> X para qualquer X Î n+tal que X > (X), e Xl ® X+ quando l ® ¥. Mais ainda rs (+) < 1, onde + (·) = ((·),...,(·)) é definido como (·) =
i(·), para i = 1,...,N, e= i(X+)
= Ai + Bi.
Observação 1 O passo c) do Lema 5, que corresponde ao cálculo da solução do sistema linear
Xl= l(Xl) + (Fl),
pode ser visto como o passo de avaliação de estratégias na técnica de iteração de estratégias para PMD, enquanto que da identidade (18),
e o lado direito é minimizado em Fi escolhendo Fi = = i(Xl), que pode ser visto como o passo de melhoramento de estratégias.
Encerramos esta seção com o seguinte resultado que será bastante útil ao longo desta seção. Num espaço probabilístico apropriado (F,, {St},S), considere duas seqüências de processos estocásticos {W(t);t = 0,1,...}, {g(t);t = 0,1,...} tal que para cada t = 0,1,...,W(t) é uma matriz n × n e g(t) é uma variável escalar positiva St-adaptada. Assuma que -quase sempre temos que
e
Assuma também, que o termo do ruído satisfaz às seguintes condições
e
onde {A(t);t = 0,1,...} é um processo estocástico tal que para cada t = 0,1,..., A(t) é uma variável escalar positiva St-adaptada. Considere um processo estocástico {R(t);t = 0,1,...}, R(t) uma matriz n × n, dada pela seqüência
Lema 6 Suponha que as equações (19)-(22) são satisfeitas e que as seqüências {R(t);t = 0,1,...} de matrizes n × n são dadas por (23). Se a seqüência A(t) é limitada -quase sempre então R(t) converge a zero -quase sempre.
Prova: Veja Corolário 4.1, página 161 em (Bertesekas e Tsitsilklis, 1996).
3 MÉTODO DE DIFERENÇA TEMPORAL QUANDO NÃO SE CONHECE A MATRIZ
Se a matriz de transição de probabilidade é conhecida, podemos utilizar o algoritmo descrito no Lema 5, para encontrar a solução maximal P do EARA dada pelas equações (3),(4) e então obter a lei de controle ótima F dada pela equação (6). Se a matriz de transição de probabilidade não é conhecida poderíamos tentar desenvolver um algoritmo de Monte Carlo para obter a solução Xl como no Lema 5 parte c). Entretanto, isto não é suficiente para obter Fl+1 = (Xl), como se pode observar da equação (17), pois há uma dependência em através do operador (·). Uma alternativa é calcular diretamente o termo Sl: = (Xl), que nos leva à seguinte equação:
ou
Sl = l(Sl) + ((Fl)).
Escrevemos Ql = ((Fl)) para simplificar a notação. Como foi visto na seção 2, rs (l) = rs (l) < 1, e o Lema 2 pode ser usado de forma que a equação em Y
tem uma única solução Sl dada por
Um vez que se tenha calculado Sl, Fl+1 pode ser obtida da equação como
O restante desta seção será dedicado a calcular Sl, via simulações de Monte Carlo, e traçando um paralelo com o método do TD() (veja (Bertesekas e Tsitsilklis, 1996)). Para simplificar a notação, eliminaremos o superescrito l. Para qualquer Î [0,1) definimos o operador Î (n) e Z Î n da seguinte maneira
e
Obtemos então o seguinte resultado:
Proposição 1 rs () < 1.
Prova: Como foi mencionado acima, rs () = rs () < 1. Suponha que tenhamos rs () > 1. Então para algum b Î , |b| > 1, e V Î n
(V) = bV.
Neste caso,
e
de forma que
e é um autovalor de com o autovetor associado V. Porém
que é uma contradição com o fato que rs () < 1. Deste modo rs () < 1.
Fazendo iterações na equação (24) com Y = S, temos que
Temos da equação (27) que
A equação (28) evidencia o seguinte método de diferenças temporais. Seja X := {1,...,N} ¥ e para cada i = 1,...,N e t = 1,2,... considere as seguintes variáveis aleatórias Qi(t) = (qi(t,0), qi(t,1),...) Î X tal que qi(t,0) = i e {qi(t,k)} tem a mesma distribuição que {q(k)}. Para k = 0,1,... o operador afim limitado (t,k,·) é definido, para V = (V1,...,VN) Î n, em termos de i(t,k) Î n e i(t,k,·) Î (n) como
Note que
E (
i(t,k,V)| q(0) = i) = (Q) + (V) - (V)e definimos para i = 1,...,N
Considere {g(t);t = 1,2,...} satisfazendo as equações (19), (20) com probabilidade 1 e para um Y(0) = (Y1(0),...,YN(0)) Î n arbitrário defina para t = 1,2,..., a seqüência Y(t) = (Y1(t),...,YN(t)) Î n da seguinte maneira:
Yi(t + 1) = Yi(t) + g(t)
ki(t,k) i = 1,...,N,ou de uma forma mais compacta,
Esta equação será fundamental no desenvolvimento do método pois, como será demonstrado no Teorema 7, Y(t) converge para S. Conforme fica evidenciado pelo somatório na equação (29), para a implementação do algoritmo necessitaríamos de um grande número de realizações da cadeia de Markov q(t), geradas computacionalmente ou através de observações passadas, e representadas na equação (29) por qi(t,k),k = 0,1,.... Note que da equação (28) podemos escrever a equação (29) da seguinte maneira:
Denotamos por St a história do processo até o tempo t, a qual é definida como
St = s{Y(0),...,Y(t),Qi (s),
s = 1,...,t - 1,i = 1,...,N,g(s),s = 1,...,t}.
Então para cada i = 1,...,N
Na seguinte proposição considere como na equação (16) com Gi = Ai + BiFi.
Proposição 2 Se 2rs() < 1 então existem constantes a > 0, b > 0 tais que
(||Wi(t)||2| St) < a||Y(t)||2 + b.
Prova: Visto que Wi(t) é hermitiana, podemos encontrar um x Î n, ||x||2 = 1, tal que ||Wi(t)||2 = | x* Wi(t)x |. Observe que
onde
Como visto na seção 2, E(||x(k) ||2) = tr{Vi(k)}, onde Vi(k) = E(x(k)x(k)*), V(k) = (V1(k),...,VN(k)) Î n+
V(k + 1) = (V(k)), k = 0,1,...
com Gi = Ai + BiFi em (14). Como rs() < 1, podemos encontrar c0 > 0, c1Î (0,1) tais que
||
k || < c0.Dessa forma, definindo c3 = ||C||max + ||D||max ||F||max e observando que
segue que
e similarmente
Portanto para alguma constante c4 > 0,
Finalmente observe que do Lema 3 temos que
Se
2rs () < 1 então podemos encontrar c5 > 0, c6Î (0,1) tais que 2k||k|| <c5e
Deste modo para algum c8 > 0,
e o resultado segue.
Assim as equações (21) e (22) são satisfeitas para cada Wi(t), i = 1,...,N. O seguinte resultado mostra que, através de uma transformação no algoritmo dado por (30) podemos assumir que || || < 1.
Observação 2 Sem perder a generalidade, assumimos que || || < 1.
Suponha que o resultado não seja válido. Visto que rs () < 1 (Proposição 1), sabemos que do Lema 1 existe um Î () invertível e real e Î () tal que ||| < 1 e
=
-1.Considere a seguinte transformação, no algoritmo descrito pela equação (29)
Então da equação (29) segue que
a qual pode ser reescrita como:
onde
= Z, (t) = W(t).
É fácil verificar-se que (i(t) | t) = 0 e (||i(t)||2 | t) < ã ||(t) ||2 + para algum ã > 0 , > 0 onde t é definido de forma apropriada com (s) no lugar de Y(s). Assim, se |||| > 1, podemos aplicar a transformação linear acima e trabalhar com a representação do algoritmo dado pela equações (31), (32), em lugar de (29), (30).
Assim, de agora em diante, assumimos, sem perda de generalidade que, || || < 1 na equação (30). O seguinte resultado segue os mesmos argumentos apresentados em (Bertesekas e Tsitsilklis, 1996), páginas 162-167.
Proposição 3 Se 2r s() < 1 então a seqüência {Y(t);t = 1,2,...} é limitada com probabilidade 1.
Prova: Esta prova segue os mesmos passos que a prova da Proposição 4.7 em (Bertesekas e Tsitsilklis, 1996), página 162-166. Como || || < 1 podemos encontrar um G tal que G > 1 e G > , e definimos h e como || || < h := + |||| Î (0,1), := - 1 > 0. Como em (Bertesekas e Tsitsilklis, 1996), página 163, a seqüência St-adaptada não-decrescente {G(t);t = 1,2,...} é definida como segue: G(0) = max{||Y(0)||,G} e
onde
= min{s > 0;||Y(t + 1)|| < (1 + ) sG(0)}.
Conseqüentemente para todo t = 1,2,..., ||Y(t)|| < (1 + )G(t) e ||Y(t)|| < G(t) se G(t - 1) < G(t). Definimos (t) = W(t), de forma que
e para c = a(1 + )2 + b
Definindo para i = 1,...,N, t > t0> 1, Ri(t0, t0) = 0 e
Ri(t + 1,t0) = (1 - g(t))Ri(t, t0) + g(t)i(t)
temos do Lema 6, (33) e (34), que Ri(t, t0) tende para zero, quando t tende para infinito com probabilidade 1. Considere o conjunto G Ì S, tal que as equações (19), (20) são válidas e Ri(t,t0) tende a zero, quando t tende a infinito para cada t0 = 1,2,..., i = 1,...,N. É fácil checar que G tem probabilidade 1 (é a intersecção de conjuntos contáveis com probabilidade 1). Suponha por contradição que Y(t)(w), w Î G, é ilimitado, de forma que G(t) (w)® ¥ a medida que t ® ¥ e ||Y(t)(w)|| < G(t)(w) para infinitamente freqüentes t. Seja u(w) > 1 tal que g(s)(w) < 1 para s > u(w). Considere t0(w) > u(w) tal que ||Y(t0)(w) || < G(t0)(w) e ||Ri(t,0)(w) || < para i = 1,...,N, t > t0 (w) > u(w). Visto que
Ri(t,0)(w) =(1 - g(s))Ri(t0,0) (w) + Ri(t,t0)(w)
temos para t > t0(w) > u(w) e R(t,t0) (w) = (R1(t,t0)(w),...,RN (t,t0)(w)) Î n que
Além disso temos que
Falta mostrar por indução que para t > t0(w)
e
Para t = t0(w) o resultado é trivial pois R(t0,t0) (w) = 0 e por hipótese, ||Y(t0) (w)|| < G(t0)(w). Suponha que as equações (37) e (38) são válidas para t. Então (por simplicidade suprimimos w) a partir da equação (38), temos que
Y(t + 1) - R(t + 1,t0)G(t0)=
1 - g(t))(Y(t) - R(t,t0)G(t0))
+ g(t)(Z + (Y(t)) + W(t)) - g(t)(t)G(t)
= (1 - g(t))(Y(t) - R(t,t0)G(t0))
+ g(t)(Z + (Y(t))) + g(t)W(t) - g(t) W(t)
= (1 - g(t))(Y(t) - R(t,t0)G(t0)) + g(t)(Z + (Y(t)))
e das equações (36), (37), (38) que
||Y(t + 1) - R(t + 1,t0)G(t0)||
< (1 - g(t)) ||Y(t) - R(t,t0)G(t0)||
+ g(t)||Z + (Y(t)) ||
< (1 - g(t))G (t0) + g(t)G(t0)=G(t0)
provando (37) para t + 1. Da equação (35) temos que
||Y(t + 1)|| < ||Y(t + 1) - R(t + 1,t0)G(t0)||
+ ||R(t + 1,t0)||G(t0)
<G(t0) + G(t0) = G(t0)(1 + ) = G(t)(1 + )
e pela definição G(t + 1) = G(t) = G(t0) provando para t + 1. Entretanto isto torna-se uma contradição, com o fato de que G(t) tende para infinito, provando o resultado.
Obtemos agora o resultado principal desta seção, provando a convergência de Y(t) para S.
Teorema 7 Se2rs () < 1 então a seqüência {Y(t);t = 1,2,...} converge para S com probabilidade 1.
Prova: Esta prova segue os mesmos passos que a prova da Proposição 4.5 no livro (Bertesekas e Tsitsilklis, 1996), páginas 166-167. Primeiro note que, a partir da equação (28),
S = Z + (S) = (1 - g(t))S + g(t)(Z + (S))
e então definindo X(t) = Y(t) - S, temos
X(t + 1) = (1 - g(t))X(t) + g(t) ((X(t)) + W(t)).
Defina para i = 1,...,N, t > t0> 1, Ri(t0,t0) = 0 e
Ri(t + 1,t0) = (1 - g(t))Ri(t, t0) + g(t)Wi(t).
A partir do Lema 6, Proposição 2 e 3 temos que Ri(t,t0) tende a zero, à medida que t tende a infinito com probabilidade 1. Seja L Ì S um conjunto tal que a seqüência {Y(t);t = 1,2,...} é limitada e Ri(t,t0) tende a zero à medida que t tende a infinito para todo i = 1,...,N, t0> 1. Visto que este conjunto é a intersecção de conjuntos contáveis com probabilidade 1, temos que (L) = 1. Escrevemos R(t,t0) = (R1(t,t0) ,...,RN(t,t0)) Î n. Para cada w Î L existe um d(w) tal que para todo t > 1, temos que
||X(t) (w)|| < ||Y(t)(w) || + ||S|| < d(w).
Usamos n > 0 tal que |||| + n < 1 e dk+1 = (|||| + n)dk, d0 = d. Pegamos t0(w) = u(w) onde u(w) é tal que g(s) (w) < 1 para s > u(w)) e provamos que existe sempre um tk + 1(w) > tk(w) tal que
Para simplificar a notação, omitimos o w. Para t = 0 o resultado é direto. Suponha que a equação (39) é válida para k. Considere a seguinte seqüência
Mostramos por indução que
Para t = tk temos que y(tk) = dk, R(tk,tk) = 0 e o resultado segue da equação (39). Suponha que a equação (41) é válida para t. Então
X(t + 1) - R(t + 1,tk) =(1 - g(t))(X(t) - R(t,tk))
+ g(t) (X(t))
e visto que || (X(t)) || < || || ||X(t)||, temos das equações (39), (40) e (41) que
||X(t + 1) - R(t + 1,tk)|| (1 - g(t))y(t) + g(t) ||||dk
= y(t + 1)
mostrando (41) para t + 1. Desta forma a equação (41) é válida para todo t > tk. Como y(t) converge para || || dk e R(t,tk) tende a 0 à medida que t tende a infinito, podemos encontrar um tk+1> tk tal que y(t) < (|||| + )dk e ||R(t,tk) || < dk para todo t > tk+1. Deste modo da equação temos que para todo t > tk+1
||X(t)|| < ||X(t) - R(t,tk) || + ||R(t,tk)||
<y(t)+ dk (|| || + n)dk = dk+1
provando a equação (39) para k + 1. Como dk tende a zero à medida que k tende a infinito, o resultado segue.
Observação 3 Na prática devemos testar para verificar se a convergência de Y(t) para S ocorre, visto que é desconhecido, e deste modo não é possível checar de antemão se 2rs () < 1. Deve-se destacar também que o algoritmo necessita de um F0 que estabiliza (A,B) na média quadrática para iniciar as iterações.
4 EXEMPLO NUMÉRICO
Para ilustrar o uso do resultados desenvolvidos nas seções anteriores, escolhemos um sistema econômico simples baseado no modelo multiplicador-acelerador de Samuelson's (Blair e Sworder, 1975), o qual aparece na forma de equações de estado:
x(k + 1) = Aq(k)x(k) + Bq(k)u(k).
Consideramos os seguintes três modos de operação do sistema linear a saltos Markovianos a tempo discreto:
A matriz de transição de probabilidades para este sistema é assumida como:
Definimos o algoritmo para = 0,1,..., da seguinte forma:
i) é calculado de acordo com a equação (26) (exceto para = 0).
ii) é o valor estacionário da equação (29).
Conforme mencionado anteriormente, fica evidenciado pelo somatório na equação (29) que para a implementação do algoritmo necessitaríamos de um grande número de realizações da cadeia de Markov q(t), geradas computacionalmente ou através de observações passadas, e representadas na equação (29) por qi(t,k),k = 0,1,.... Deve-se notar também que, conforme mencionado na Observação 3, necessita-se de um F0 que estabilize (A,B) na média quadrática para iniciar as iterações.
Foram realizados vários experimentos com diferentes valores (vide Tabela 2). Na coluna 3 da Tabela 2 mostra-se o erro, que foi calculado de acordo com a seguinte equação:
onde é o ganho do controlador ótimo associado ao modo i.
Além dos dados observados na tabela anterior, mostraremos também, nas figuras seguintes, a evolução do valor de Fi em relação ao número de iterações , para um dos valores de (para os outros valores o comportamento é bastante semelhante).
5 CONCLUSÕES
Neste trabalho traçou-se um paralelo com o método de simulação Monte Carlo-TD () para PMD (veja por exemplo (Bertesekas e Tsitsilklis, 1996), (Sutton e Barto, 1998)), para obter o controle ótimo associado ao conjunto de equações algébricas de Riccati acopladas entre si (EARA), para controle ótimo de sistemas lineares com saltos Markovianos a tempo discreto. Assumimos que a matriz de transição de probabilidade é desconhecida, mas é possível simular as trajetórias da cadeia de Markov q(k). Relacionamos as iterações (8),(9) com aquelas associadas à técnica de iteração de valores em PMD, e o método chamado de quase-linearização apresentado no Lema 5 com a técnica de iteração de estratégias para PMD, envolvendo uma parte de avaliação de estratégias, e a parte de melhoramento de estratégias (veja Observação 1).
Na técnica de iteração de estratégias, a avaliação de estratégias é feita via simulações de Monte Carlo, usando o método de diferenças temporais. Aplicando estas idéias ao método de quase-linearização apresentado no Lema 5, obtivemos um método iterativo que traça um paralelo com os métodos de simulação TD() em PMD. Para isto foi mostrado no Teorema 7 que se o escolhido é suficientemente pequeno, a convergência do algoritmo TD() na avaliação de custo acontece com probabilidade 1. São apresentados exemplos numéricos para ilustrar os resultados.
Artigo submetido em 14/12/2000
1a. Revisão em 19/12/2002
Aceito sob recomendação do Ed. Assoc. Prof. Liu Hsu
Referências bibliográficas
- Abou-Kandil, H., Freiling, G. e Jank, G. (1995). On the solution of discrete-time markovian jump linear quadratic control problems, Automatica 31(5): 765768.
- Ait-Rami, M. e Ghaoui, L. E. (1996). LMI optimization for nonstandard riccati equations arising in stochastic control, IEEE Trans. Automatic Control 41(11): 1666 1671.
- Bertesekas, D. P. e Tsitsilklis, J. (1996). Neuro-Dynamic Programming, Athena Scientific.
- Blair, W. e Sworder, D. (1975). Feedback control of a class of linear discrete system with jump parameters and quadratic cost criteria, Int. J. Control 21: 833841.
- Brewer, W. (1978). Kronecker product and matrix calculus in system theory, IEEE Trans. Circuits and Systems 25: 772781.
- Costa, O. e Aya, J. (1999). Temporal difference methods for the maximal solution of discrete-time coupled algebraic riccati equations, Proc. American Control Conference pp. 17911795.
- Costa, O. e Fragoso, M. (1993). Stability results for discretetime linear systems with markovian jumping parameters, J. Math. Analysis and Applic 179: 154178.
- Costa, O. e Fragoso, M. (1995). Discrete-time LQ-optimal control problems for infinite markov jump parameter systems, IEEE Trans. Automat. Control 40: 20762088.
- Costa, O. e Marques, R. (1999). Maximal and stabilizing hermitian solutions for discrete-time coupled algebraic riccati equations, Mathematics of Control, Signals and Systems 12(2): 167195.
- Costa, O., Val, J. D. e Geromel, J. (1997). A convex programming approach to H2-control of discrete-time markovian jump linear systems, Int. J. Control 66: 557579.
- Gajic, Z. e Borno, I. (1995). Lyapunov iterations for optimal control of jump linear systems at steady state, IEEE Trans. Automat. Control 40(11): 481498.
- Ji, Y. e Chizeck, H. (1988). Controllability, observability and discrete-time markovian jump linear quadratic control, Int. J. Control 48: 481498.
- Ji, Y., Chizeck, H., Feng, X. e Loparo, K. (1991). Stability and control of discrete-time jump linear systems, Control Th. and Adv. Tech 7: 247270.
- Kubrusly, C. (1997). An Introduction to Models and Decompositions in Operator Theory, Springer Verlag, New York.
- Mariton, M. (1988). Almost sure and moments stability of jump linear systems, Systems and Control Letters 11: 393397.
- Mariton, M. (1990). Jump Linear Systems in Automatic Control, Marcel Dekker.
- Sutton, R. S. e Barto, A. G. (1998). Reinforcement Learning - An Introduction, MIT Press.
- Val, J. D., Geromel, J. e Costa, O. (1998). Uncoupled riccati iterations for the linear quadratic control problem of discrete-time markov jump linear systems, IEEE Trans. Automat. Control 43(12): 17271733.
- Weidman, J. (1980). Linear Operators in Hilbert Spaces, Springer Verlag, New York.
Datas de Publicação
-
Publicação nesta coleção
19 Nov 2003 -
Data do Fascículo
Set 2003
Histórico
-
Recebido
14 Dez 2000