Método de diferenças temporais aplicado às equações de Riccati acopladas entre si

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, [e-mail], [e-mail]

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 = (p_ij). 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 = (P₁,...,P_N), P_i > 0 i = 1,...,N, do seguinte conjunto de equações algébricas de Riccati acopladas entre si (EARA), para i = 1,...,N

onde (X) = (₁(X),..., _N(X)) é definido, para X = (X₁,...,X_N) 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 = (F₁,...,F_N) é dado por

Associado aos problemas de horizonte infinito (1), (2), temos o problema de horizonte finito (7) no qual, para algum S = (S₁,...,S_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) = F_{q (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 r_s (). 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 ⁿ e por (ⁿ, ^m) o espaço linear de todas as matrizes m × n com norma limitada, e (ⁿ) := ( ⁿ, ⁿ) e o produto interno em (ⁿ) é dado por áH;Vñ = tr{H*V} para H,V Î (ⁿ). O superescrito ' denotará a transposta de uma matriz. Usamos neste trabalho ||·||₂ para a norma Euclidiana em ⁿ e a norma induzida em (ⁿ).

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) r_s () < 1

b) existe umÎ (ⁿ) 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 r_s () < 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 = (V₁,...,V_N) com V_i Î (ⁿ, ^m), i = 1,...,N, e por simplicidade ⁿ = ^n,n. Neste trabalho consideraremos o conjunto ⁿ equipado com o seguinte produto interno: para H = (H₁,...,H_N), V = (V₁,...,V_N) Î ⁿ definimos o produto interno á·;·ñ em ⁿ como segue:

e ||H||²: = áH;Hñ (de forma que ||H||² = ||H_i||²). Então com o produto interno acima definido (vide equação (12)), ⁿ é um espaço de Hilbert. Nos também temos a seguinte norma ||· ||_max em ⁿ. Para H = (H₁,...,H_N) Î ⁿ, ||H||_max é definido por

||H ||_max := max{||H_i||₂ ;i = 1,...,N}.

Definimos

ⁿ⁺ := {V = (V₁,...,V_N) Î ⁿ; V_i > 0,"i}

e escrevemos que V > S onde V = (V₁,...,V_N) Î ⁿ e S = (S₁,...,S_N) Î ⁿ, se

V - S = (V₁ - S₁,...,V_N - S_N) Î ⁿ⁺,

e que V > S se V_i - S_i > 0 para i = 1,...,N. Para, G = (G₁,...,G_N) Î ⁿ definimos os seguintes operadores, (·) = (₁(·),..., _N(·)) Î (ⁿ) e (·) = (₁(·),..., _N(·)) Î (ⁿ) para V = (V₁,...,V_N) Î ⁿ e i,j = 1,...,N,

onde o operador (·) = (₁(·),..., _N(·)) Î (ⁿ) é definido como em (4). É simples verificar que os operadores , , e mapeiam ⁿ⁺ em ⁿ⁺, e que r_s () = r_s (). No artigo (Costa e Fragoso, 1993) se prova que r_s () = r_s () (de fato = ^* de acordo com a equação (12)). Também definimos o operador Î (ⁿ) como

onde, de novo, V = (V₁,...,V_N) Î ⁿ e (V) = (₁(V),..., _N(V)). No artigo de (Costa e Fragoso, 1993) provou-se que r_s () = r_s () = r_s ().

Assumimos no modelo (1) e na função de custo (2) que A = (A₁,...,A_N) Î ⁿ, B = (B₁,...,B_N) Î ^m,n, C = (C₁,...,C_N) Î ^n,p e D = (D₁,...,D_N) Î ^m,p. Como for mostrado em (Costa e Fragoso, 1993), (Mariton, 1988), o modelo (1) com u(k) = F_q(k)x(k), e V_i(k) = E(x(k)x(k)^*1 _{{q(k) = i}}), V(k) = (V₁(k),...,V_N(k)) Î ⁿ⁺ leva a

V(k + 1) = (V(k)), k = 0,1,...

onde G_i = A_i + B_iF_i em (14), e E(|| x(k)||²) = tr{V_i(k)}.

O próximo passo é definir o operador para o 4^essimo momento de x(t). Para isto usamos o produto de Kronecker L

K Î () para L,K Î (ⁿ) e o operador vec{·}: (ⁿ) definido da forma comum (veja (Brewer, 1978)). Seja o operador Î () da seguinte forma: para S = (S₁,...,S_N) Î , (S) = (₁(S),..., _N(S)) é definido como

onde S_i(k) = E(vec{x(k) x(k)^*}vec{x(k)x(k)^*}^*1 _{{q(k) = i}}), S(k) = (S₁(k),...,S_N(k)) Î +. Temos o seguinte resultado:

Lema 3 S(k + 1) = (S(k)) e E(||x(k)||⁴) = tr{S_i(k)}.

Prova: Como x(k + 1) = G_q(k)x(k), é fácil de se verificar que

x_j(k + 1)x_j(k + 1)^* = 1_{{q(k+1) = j}}G_ix_i(k)x_i (k)^*

onde x_i(k) := x(k)1 _{{q(k) = i}}. Seja z_i(k) = vec{x_i(k)x_i(k)^*}, i = 1,...,N. Depois de algumas manipulações segue que

z_j(k + 1)z_j(k + 1)* = 1_{q(k+1)=j}

.(iG_i)z_i(k)z_i(k)* (

)

e escrevemos

S_i(k) = E(z_i(k)z_i(k)^* 1_{{q(k) = i}}).

Conseqüentemente temos que

S_j(k + 1) = p_ij (iG_i)S_i(k)(

).

Finalmente, note que

tr{z_i(k)z_i (k)*} = ||x_i(k)||²tr{x_i(k)x_i(k)}

=||x_i(k) ||⁴.

Vamos definir o conceito de estabilidade que será usado nas seções seguintes.

Definição 1 Diremos que F = (F₁,...,F_N) Î ^n,m estabiliza (A,B) no sentido da média quadrática se, quando fazemos u(k) = F_q(k)x(k) no sistema (1), temos que E(||x(k)|| ²) ® 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 = (F₁,...,F_N) Î ^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 = (F₁,...,F_N) 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 = (F₁,...,F_N) Î ^n,mestabiliza (A,B) no sentido da média quadrática se e somente se r_s () < 1, onde é como em (14) com G_i = A_i + B_iF_i.

Definição 2 Definimos (·) = (₁(·),..., _N(·)): ⁿ⁺®^n,m, (·) = (₁(·),..., _N(·)): ^n,m ®ⁿ e = (₁(·),..., R_N(·)):ⁿ⁺ ® ⁿ⁺ como

onde X = (X₁,...,X_N) Î ⁿ⁺, F = (F₁,...,F_N) Î ^n,m.

A seguinte identidade será útil nas próximas seções: para qualquer F = (F₁,...,F_N) Î ^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 F⁰ = (,...,) Î ^n,m tal que estabiliza (A,B) no sentido da média quadrática. Então para l = 0,1,2,..., existe X^l = (,...,) que satisfaz às seguintes propriedades:

a) X⁰>X¹>...>X^l>X, para um X Î ⁿ⁺ arbitrário tal que X >

(X).

b) r_s (^l) < 1, onde ^l (·) = ((·),...,(·)) e para i = 1,...,N,

c) X^l satisfaz X^l = ^l (X^l) + (F^l) e é dado por X^l = (^l)^k ( (^l)).

Além disso existe um X⁺ = (,...,) Î ⁿ⁺tal que X⁺ = (X⁺), X⁺> X para qualquer X Î ⁿ⁺tal que X > (X), e X^l ® X⁺ quando l ® ¥. Mais ainda r_s (⁺) < 1, onde ⁺ (·) = ((·),...,(·)) é definido como (·) =

_i(·), para i = 1,...,N, e

= _i(X⁺)

= A_i + B_i.

Observação 1 O passo c) do Lema 5, que corresponde ao cálculo da solução do sistema linear

X^l= ^l(X^l) + (F^l),

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 F_i escolhendo F_i = = _i(X^l), 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,, {S_t},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 S_t-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 S_t-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 X^l como no Lema 5 parte c). Entretanto, isto não é suficiente para obter F^l+1 = (X^l), como se pode observar da equação (17), pois há uma dependência em através do operador (·). Uma alternativa é calcular diretamente o termo S^l: = (X^l), que nos leva à seguinte equação:

ou

S^l = ^l(S^l) + ((F^l)).

Escrevemos Q^l = ((F^l)) para simplificar a notação. Como foi visto na seção 2, r_s (^l) = r_s (^l) < 1, e o Lema 2 pode ser usado de forma que a equação em Y

tem uma única solução S^l dada por

Um vez que se tenha calculado S^l, F^l+1 pode ser obtida da equação como

O restante desta seção será dedicado a calcular S^l, 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 Î (ⁿ) e Z Î ⁿ da seguinte maneira

e

Obtemos então o seguinte resultado:

Proposição 1 r_s () < 1.

Prova: Como foi mencionado acima, r_s () = r_s () < 1. Suponha que tenhamos r_s () > 1. Então para algum b Î , |b| > 1, e V Î ⁿ

(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 r_s () < 1. Deste modo r_s () < 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 Q_i(t) = (q_i(t,0), q_i(t,1),...) Î X tal que q_i(t,0) = i e {q_i(t,k)} tem a mesma distribuição que {q(k)}. Para k = 0,1,... o operador afim limitado (t,k,·) é definido, para V = (V₁,...,V_N) Î ⁿ, em termos de _i(t,k) Î ⁿ e _i(t,k,·) Î (ⁿ) 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) = (Y₁(0),...,Y_N(0)) Î ⁿ arbitrário defina para t = 1,2,..., a seqüência Y(t) = (Y₁(t),...,Y_N(t)) Î ⁿ da seguinte maneira:

Y_i(t + 1) = Y_i(t) + g(t)

^k

_i(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 q_i(t,k),k = 0,1,.... Note que da equação (28) podemos escrever a equação (29) da seguinte maneira:

Denotamos por S_t a história do processo até o tempo t, a qual é definida como

S_t = s{Y(0),...,Y(t),Q_i (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 G_i = A_i + B_iF_i.

Proposição 2 Se ²r_s() < 1 então existem constantes a > 0, b > 0 tais que

(||W_i(t)||²| S_t) < a||Y(t)||² + b.

Prova: Visto que W_i(t) é hermitiana, podemos encontrar um x Î ⁿ, ||x||₂ = 1, tal que ||W_i(t)||₂ = | x^* W_i(t)x |. Observe que

onde

Como visto na seção 2, E(||x(k) ||²) = tr{V_i(k)}, onde V_i(k) = E(x(k)x(k)^*), V(k) = (V₁(k),...,V_N(k)) Î ⁿ⁺

V(k + 1) = (V(k)), k = 0,1,...

com G_i = A_i + B_iF_i em (14). Como r_s() < 1, podemos encontrar c₀ > 0, c₁Î (0,1) tais que

||

^k || < c₀.

Dessa forma, definindo c₃ = ||C||_max + ||D||_max ||F||_max e observando que

segue que

e similarmente

Portanto para alguma constante c₄ > 0,

Finalmente observe que do Lema 3 temos que

Se

²r_s () < 1 então podemos encontrar c₅ > 0, c₆Î (0,1) tais que

^2k||^k|| <c₅

e

Deste modo para algum c₈ > 0,

e o resultado segue.

Assim as equações (21) e (22) são satisfeitas para cada W_i(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 r_s () < 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)||² | _t) < ã ||(t) ||² + 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 ²r _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 S_t-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 + )² + b

Definindo para i = 1,...,N, t > t₀> 1, R_i(t₀, t₀) = 0 e

R_i(t + 1,t₀) = (1 - g(t))R_i(t, t₀) + g(t)i(t)

temos do Lema 6, (33) e (34), que R_i(t, t₀) 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 R_i(t,t₀) tende a zero, quando t tende a infinito para cada t₀ = 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 t₀(w) > u(w) tal que ||Y(t₀)(w) || < G(t₀)(w) e ||R_i(t,0)(w) || < para i = 1,...,N, t > t₀ (w) > u(w). Visto que

R_i(t,0)(w) =(1 - g(s))R_i(t₀,0) (w) + R_i(t,t₀)(w)

temos para t > t₀(w) > u(w) e R(t,t₀) (w) = (R₁(t,t₀)(w),...,R_N (t,t₀)(w)) Î ⁿ que

Além disso temos que

Falta mostrar por indução que para t > t₀(w)

e

Para t = t₀(w) o resultado é trivial pois R(t₀,t₀) (w) = 0 e por hipótese, ||Y(t₀) (w)|| < G(t₀)(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,t₀)G(t₀)=

1 - g(t))(Y(t) - R(t,t₀)G(t₀))

+ g(t)(Z + (Y(t)) + W(t)) - g(t)(t)G(t)

= (1 - g(t))(Y(t) - R(t,t₀)G(t₀))

+ g(t)(Z + (Y(t))) + g(t)W(t) - g(t) W(t)

= (1 - g(t))(Y(t) - R(t,t₀)G(t₀)) + g(t)(Z + (Y(t)))

e das equações (36), (37), (38) que

||Y(t + 1) - R(t + 1,t₀)G(t₀)||

< (1 - g(t)) ||Y(t) - R(t,t₀)G(t₀)||

+ g(t)||Z + (Y(t)) ||

< (1 - g(t))G (t₀) + g(t)G(t₀)=G(t₀)

provando (37) para t + 1. Da equação (35) temos que

||Y(t + 1)|| < ||Y(t + 1) - R(t + 1,t₀)G(t₀)||

+ ||R(t + 1,t₀)||G(t₀)

<G(t₀) + G(t₀) = G(t₀)(1 + ) = G(t)(1 + )

e pela definição G(t + 1) = G(t) = G(t₀) 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 Se²r_s () < 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 > t₀> 1, R_i(t₀,t₀) = 0 e

R_i(t + 1,t₀) = (1 - g(t))R_i(t, t₀) + g(t)W_i(t).

A partir do Lema 6, Proposição 2 e 3 temos que R_i(t,t₀) 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 R_i(t,t₀) tende a zero à medida que t tende a infinito para todo i = 1,...,N, t₀> 1. Visto que este conjunto é a intersecção de conjuntos contáveis com probabilidade 1, temos que (L) = 1. Escrevemos R(t,t₀) = (R₁(t,t₀) ,...,R_N(t,t₀)) Î ⁿ. 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 d_k+1 = (|||| + n)d_k, d₀ = d. Pegamos t₀(w) = u(w) onde u(w) é tal que g(s) (w) < 1 para s > u(w)) e provamos que existe sempre um t_k + 1(w) > t_k(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 = t_k temos que y(t_k) = d_k, R(t_k,t_k) = 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,t_k) =(1 - g(t))(X(t) - R(t,t_k))

+ 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,t_k)|| (1 - g(t))y(t) + g(t) ||||d_k

= y(t + 1)

mostrando (41) para t + 1. Desta forma a equação (41) é válida para todo t > t_k. Como y(t) converge para || || d_k e R(t,t_k) tende a 0 à medida que t tende a infinito, podemos encontrar um t_k+1> t_k tal que y(t) < (|||| + )d_k e ||R(t,t_k) || < d_k para todo t > t_k+1. Deste modo da equação temos que para todo t > t_k+1

||X(t)|| < ||X(t) - R(t,t_k) || + ||R(t,t_k)||

<y(t)+ d_k (|| || + n)d_k = d_k+1

provando a equação (39) para k + 1. Como d_k 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 ²r_s () < 1. Deve-se destacar também que o algoritmo necessita de um F⁰ 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) = A_q(k)x(k) + B_q(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 q_i(t,k),k = 0,1,.... Deve-se notar também que, conforme mencionado na Observação 3, necessita-se de um F⁰ 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 F_i 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 H₂-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