X SBAI – Simpósio Brasileiro de Automação Inteligente
18 a 21 de setembro de 2011
São João del-Rei - MG - Brasil
FILTRAGEM ESTOCÁSTICA APLICADA A SISTEMAS MAX-PLUS LINEARES
D IEGO F IGUEIRÊDO E S ILVA∗, R AFAEL S ANTOS M ENDES∗, L AURENT H ARDOUIN†, C ARLOS A NDREY M AIA‡,
B ERTRAND C OTTENCEAU†
∗
†
Faculdade de Engenharia Elétrica e de Computação - UNICAMP
DCA/FEEC/UNICAMP-C.P. 6101
13083-970 Campinas - SP - Brasil
Laboratoire d’Ingénierie des Systèmes Automatisés - ISTIA/Université d’Angers
Av. Notre Dame du Lac 62
49000 Angers - França
‡
Depto. de Engenharia Elétrica - UFMG
Av. Antônio Carlos 6627 - Pampulha
31270-010 Belo Horizonte - MG - Brasil
Emails: [email protected], [email protected], [email protected],
[email protected], [email protected]
Abstract— This paper proposes a particle filter for linear max-plus systems. A brief introduction to discrete event systems is
presented in a max-plus context. The paper describes some fundamental results in particle filtering and applies these results to the
max-plus linear systems leading to a stochastic filtering algorithm. An example illustrates the method and some conclusions are
drawn.
Keywords— Stochastic Filtering, Discrete Event Systems, Max-plus Algebra, Particle Filters
Resumo— Este trabalho propõe um filtro estocástico do tipo filtro de partículas para sistemas max-plus lineares. Uma pequena
introdução aos sistemas a eventos discretos no contexto da abordagem max-plus é apresentada. Em seguida os fundamentos de
filtros de partículas, aplicados aos sistemas max-plus, são descritos, chegando assim a um algoritmo de filtragem estocástica. Um
exemplo ilustrativo é apresentado e são esboçadas conclusões.
Palavras-chave—
1
Filtragem Estocástica, Sistemas a Eventos Discretos, Álgebra Max-plus, Filtros de Partículas
Introdução
O desenvolvimento das técnicas de filtragem estocástica (Brown and Hwang, 1997) vem permitindo os
avanços em várias áreas de pesquisa, em especial na
geração de tecnologia espacial e militar (Bar-Shalom
et al., 2001). Pode-se citar particularmente o filtro
de Kalman (Brown and Hwang, 1997), cuja presença
se tornou indispensável em equipamentos eletrônicos
de comunicação via frequência modulada (FM) e em
sistemas de rastreamento aéreo e de navegação (BarShalom et al., 2001).
Os Sistemas a Eventos Discretos (SED)
(Cassandras and Lafortune, 1999) constituem
uma classe de sistemas caracterizados pelo fato de
sua dinâmica ser dirigida pela ocorrência de eventos,
em oposição aos sistemas contínuos cuja dinâmica é
determinada pela passagem do tempo. Num SED as
mudanças de estado ocorrem única e exclusivamente
devido à ocorrência de eventos. O estudo destes sistemas é de fundamental importância em diversas áreas
como redes de computadores, sistemas de tráfego e
no planejamento de atividades de carácter logístico,
podendo-se destacar a sua aplicação nos sistemas de
manufatura, onde é utilizado para otimizar o tempo
de produção e o uso do estoque (Maia et al., 2005).
Dentre os modelos existentes para o tratamento de
tais sistemas, destacam-se as Redes de Petri (RP)
(Murata, 1989) e a abordagem baseada na álgebra de
ISSN: 2175-8905 - Vol. X
dióides, ou mais simplesmente, na álgebra max-plus
(Baccelli et al., 1992). Esta última é a abordagem
utilizada neste trabalho.
Um enfoque possível para os SED considera
como variáveis de modelagem os “datadores”, isto é,
as informações relativas aos instantes de ocorrência
de eventos. Neste contexto, podem ocorrer situações
práticas onde nem todos os datadores são diretamente
observados, de forma análoga ao que ocorre em sistemas contínuos nos quais o estado não é totalmente
observado. Além disso, é comum que os datadores
não observados tenham seus valores influenciados por
variáveis aleatórias, que combinadas com os aspectos
estruturais do sistema provocam variações aleatórias
nos instantes de ocorrência dos eventos.
O presente trabalho pretende abordar exatamente
esta situação, tratando do problema de obter estimativas ótimas, num contexto ruidoso, de instantes de
ocorrência de eventos não-observáveis em SED. Estas estimativas deverão ser obtidas de modo análogo
ao utilizado nos filtros de Kalman para sistemas contínuos, isto é, a partir de observações de algumas variáveis do sistema, obter-se-ão estimativas de outras
variáveis não-observadas sujeitas a perturbações estocásticas previamente modeladas. A técnica de filtragem conhecida como “Filtro de Partículas” será utilizada neste trabalho. Os filtros de partículas utilizam
o conceito de simulação de Monte-Carlo, isto é, produzem uma aproximação de uma distribuição proba-
1388
X SBAI – Simpósio Brasileiro de Automação Inteligente
18 a 21 de setembro de 2011
São João del-Rei - MG - Brasil
bilística através de um grande número de amostras às
quais são atribuídos pesos. Estas amostras evoluem
de acordo com a dinâmica do sistema e têm seus pesos atualizados a partir da existência de medidas das
saídas do sistema.
O artigo está organizado como segue. Na seção 2 os sistemas max-plus são brevemente descritos, definindo-se também o cenário de perturbações
que afetam o sistema. A seção 3 descreve o princípio de funcionamento do filtro de partículas e a maneira pela qual ele pode ser aplicado aos sistemas maxplus lineares com perturbação. Finalmente, na seção 4
apresentam-se alguns resultados de simulação numérica e na seção 5 esboçam-se algumas conclusões.
2
Álgebra Max-Plus
É possível modelar um sistema a eventos discretos
como uma rede de Petri (Baccelli et al., 1992). Uma
rede de Petri é um objeto abstrato constituído de um
grafo bi-partido e de uma marcação (função). No
grafo os nós podem ser de dois tipos: lugares e transições. A marcação associa a cada lugar um número
inteiro (diz-se que o lugar contém fichas). As transições estão relacionadas aos eventos do sistema e os lugares definem condições para que estes eventos ocorram. Uma transição pode disparar se os lugares acima
dela contêm fichas. Neste trabalho consideraremos RP
p-temporizadas, nas quais os lugares têm um tempo
de permanência mínimo para as fichas. O disparo de
uma transição altera a marcação da rede, redefinindo
as condições de disparo e produzindo a dinâmica do
sistema. Neste contexto, os recursos do sistema são
representados por fichas nas diversas posições do processo de produção, as durações das atividades são representadas pelo tempo de permanência das fichas nos
lugares e as ocorrências dos eventos são representadas
pelos disparos das transições, que consomem fichas
nos lugares anteriores à transição e produzem fichas
nos lugares seguintes à transição. Uma classe importante de redes de Petri são os Grafos de Eventos Temporizados (GET), caracterizados por apresentarem um
único arco de entrada e um único arco de saída em
cada lugar. Na figura 1 apresenta-se um grafo de eventos temporizado.
Figura 1: Exemplo de Grafo de Eventos Temporizado
ISSN: 2175-8905 - Vol. X
Nesta figura é importante observar a existência de
transições não condicionadas por nenhuma outra transição do sistema (transição u na figura) e transições
que não condicionam nenhuma outra transição do sistema (transição z na figura). Estas transições são chamadas respectivamente de entrada e de saída do sistema; as outras transições são chamadas de internas ou
transições de estado. Assume-se que apenas as transições de saída do sistema são observadas no contexto
de filtragem desenvolvido na próxima seção.
Pode-se associar a cada transição de um GET uma
seqüência crescente de números inteiros {x(k)}, para
k = 0, 1, 2, . . ., onde cada elemento da série representa o instante do k-ésimo disparo dessa transição.
Supondo agora conhecida a seqüência associada às
transições de entrada de um GET, é possível determinar as seqüências de disparo de todas as transições do
GET. De fato, considerando novamente o GET da figura 1, é possível escrever as seguintes relações entre
os instantes de disparo das transições:
x1 (k) = 1 + x2 (k − 1)
(1a)
x2 (k) = max{2 + x2 (k − 1); u(k)}
(1b)
z(k) = max{2 + x1 (k); 1 + x2 (k)}
(1c)
Renomeando o operador max como sendo ⊕ e o operador + como sendo ⊗, pode-se reescrever:
x1 (k) = 1 ⊗ x2 (k − 1)
(2a)
x2 (k) = 2 ⊗ x2 (k − 1) ⊕ u(k)
(2b)
z(k) = 2 ⊗ x1 (k) ⊕ 1 ⊗ x2 (k)
(2c)
ou, em forma matricial:
xk
=
A ⊗ xk−1 ⊕ B ⊗ uk
zk
=
C ⊗ xk
(3)
sendo:
−∞ 1
−∞
; B=
; C= 2 1
−∞ 2
0
T
xk = x1 (k) x2 (k) ; uk = u(k) e zk = z(k)
A=
Tem-se portanto um sistema de equações recursivas lineares numa nova álgebra. De modo geral, um
semi-anel idempotente (ou dióide) é caracterizado por
um conjunto e duas operações (soma e produto), notado (D, ⊕, ⊗), tal que a soma seja associativa, comutativa e idempotente (a ⊕ a = a), e o produto seja
associativo (mas não necessariamente comutativo) e
distributivo à esquerda e à direita em relação à soma.
Além disso, devem existir elementos neutros para ambas as operações, notados por ε (elemento nulo) e por
e (elemento unitário), e o elemento nulo deve ser absorvente em relação ao produto. Isto é, ∀a ∈ D,
a ⊕ ε = a, a ⊗ e = a, a ⊗ ε = ε. É imediato
perceber que o conjunto Z ∪ {−∞} munido das duas
operações ⊕ ≡ max e ⊗ ≡ + é um dióide, no qual
ε = −∞ e e = 0. Um dióide é completo se ele for
fechado em relação a somas infinitas e se o produto for
1389
X SBAI – Simpósio Brasileiro de Automação Inteligente
18 a 21 de setembro de 2011
São João del-Rei - MG - Brasil
distributivo em relação a somas infinitas. A estrutura
(Z ∪ {−∞} ∪ {∞}, max, +) é um dióide completo
usualmente denominado Max-plus e notado por Zmax .
Este exemplo utiliza o que se convenciona chamar
de datadores, isto é, sequências crescentes {x(k)} que
representam as datas ou instantes de ocorrência dos
disparos da transição x. Observa-se que, formalmente,
as equações 3 são muito similares às equações de estado em sistemas com dinâmica contínua.
Os sistemas abordados neste trabalho são, por definição os sistemas a eventos discretos que podem ser
descritos pelas equações 3 admitindo-se em geral que
xk ∈ Rn×1 , uk ∈ Rp×1 e zk ∈ Rq×1 e que as matrizes A, B e C tem dimensões apropriadas. Além disso,
duas hipóteses adicionais são consideradas.
Por um lado admite-se que a matriz C contém
exatamente um elemento e em cada linha, sendo todos
os outros iguais a havendo no máximo um elemento
e em cada coluna. Esta hipótese não é restritiva pois
qualquer sistema do tipo da equação 3 pode ser colocado nesta forma pela inclusão de linhas e colunas nas
matrizes A e B (Hardouin et al., 2010). Em outras
palavras assume-se que a matriz C tem a matriz identidade como sub-matriz e que algumas das variáveis
internas são diretamente observadas. Pode-se portanto
reescrever o modelo dado pela equação 3 da seguinte
forma:
xk
= A ⊗ xk−1 ⊕ B ⊗ uk
(4)
x00k
zk =
0
0
xk
A
sendo xk =
e
A
=
. A partir destas
x00k
A00
equações é possível escrever:
zk = A00 ⊗ xk−1 ⊕ B ⊗ uk
(5)
Por outro lado, assume-se que as perturbações que
afetam o sistema atuam exclusivamente nos elementos
da matriz A e da matriz B, admitindo-se que os elementos aij e bij destas matrizes são variáveis aleatórias independentes entre si com distribuições uniformes e conhecidas.
No contexto deste trabalho o problema de filtragem consiste portanto em obter estimativas para os
estados não diretamente observados x0k conhecendose uma sequência de medidas zj e de entradas uj ,
j = 1, . . . , k e a distribuições de probabilidade dos
elementos das matrizes A e B.
3
Filtros de Partículas aplicados a Sistemas
Max-Plus
Os Filtros de Partículas utilizam uma representação
por partículas da densidade de probabilidade do estado do sistema para realizar uma estimação sequencial de Monte-Carlo1 deste estado. Uma representação por partículas é um conjunto de amostras da variável que se deseja estimar, amostradas de acordo com
1 Uma estimação sequencial de Monte-Carlo é uma técnica para
a implementação de um filtro Bayesiano recursivo através de simulações de Monte-Carlo (Ristic et al., 2004)
ISSN: 2175-8905 - Vol. X
uma densidade chamada "densidade de importância".
O método de
R Monte-Carlo é baseado no seguinte fato.
Seja I = g(x) dx uma integral que se deseja avaliar e π(x) uma densidade de probabilidade tal que
g(x) = f (x) · π(x). Se for possível sortear N amostras da variável x de acordo com a densidade π(x)
PN
{xi , i = 1, . . . , N }, então IN = N1 i=1 f (xi ) é
uma estimativa não-polarizada da integral I. Além
disso, se as amostras da variável x forem sorteadas
segundo uma outra densidade q(x), similar2 a π(x),
PN
i
i
i
então IN =
i=1 f (x ) · w(x ), sendo w(x ) =
1
N
w̃(xi )
PN
j
j=1 w̃(x )
e w̃(xi ) =
π(xi )
q(xi ) .
A densidade q(x) é
a densidade de importância.
Considerando agora um sistema max-plus linear
conforme descrito na seção anterior, seja Xk = {xj },
j = 0, . . . , k o conjunto de todos os valores das variáveis de estado x até o instante k e seja Zk definido
de modo análogo. Uma abordagem comum para o
problema de filtragem consiste em realizar o cálculo
da probabilidade p(Xk |Zk ) a partir da probabilidade
p(Xk−1 |Zk−1 ) e do valor da medida zk . Contudo, levando em conta a equação 4 e que os elementos da matriz A são variáveis aleatórias independentes, concluise que, dado xk−1 , os vetores x0k e x00k = zk são independentes. Em outras palavras: p(X0 k |Xk−1 , zk ) =
p(X0 k |Xk−1 ). A melhor estimativa que se pode ter
do estado x0k depende exclusivamente da estimativa
do estado xk−1 . Por este motivo busca-se neste trabalho um procedimento recursivo para a obtenção da
probabilidade p(Xk−1 |Zk ), em função da probabilidade p(Xk−2 |Zk−1 ) e da medida zk . A representação para as densidades de probabilidade acima é feita
através do conceito de "partículas", isto é, supõe-se
a existência de um conjunto de amostras (ou partículas) {Xik−1 }, i = 1, . . . , N amostrado segundo a
densidade de importância q(Xk−1 |Zk ) com os respectivos pesos {wki } que representam a probabilidade
p(Xk−1 |Zk ), isto é:
p(Xk−1 |Zk ) ∼
=
N
X
wki · δ(Xk−1 − Xik−1 )
(6)
i=1
Conforme a discussão ao início desta seção temse que:
p(Xik−1 |Zk )
wki ∝
q(Xik−1 |Zk )
Deseja-se portanto, a partir de uma representação
de p(Xk−2 |Zk−1 ) (do tipo definido pela equação 6) e
da medida zk , obter uma representação por partículas
para p(Xk−1 |Zk ). Para obter fórmulas recursivas adequadas, deve-se escolher uma densidade de importância que satisfaça a seguinte condição de fatorização:
q(Xk−1 |Zk ) , q(xk−1 |Xk−2 , Zk ) · q(Xk−2 |Zk−1 ).
(7)
2 Uma densidade de probabilidade q(x) é similar a π(x) se ∀x :
π(x) > 0 ⇒ q(x) > 0.
1390
X SBAI – Simpósio Brasileiro de Automação Inteligente
18 a 21 de setembro de 2011
São João del-Rei - MG - Brasil
A densidade q(xk−1 |Xk−2 , Zk ) deve ser utilizada para expandir a partícula Xik−2 para Xik−1 , incluindo um novo estado xik−1 . A equação de atualização do respectivo peso wki é obtida a seguir. Tem-se:
p(Xk−1 |Zk ) =
= p(Xk−1 |zk , Zk−1 )
p(zk |Xk−1 , Zk−1 ) · p(Xk−1 |Zk−1 )
=
, (8)
p(zk |Zk−1 )
wki
mas dado que:
p(Xk−1 |Zk−1 ) =
= p(xk−1 |Xk−2 , Zk−1 ) · p(Xk−2 |Zk−1 )
∼
(9)
= p(xk−1 |xk−2 ) · p(Xk−2 |Zk−1 ),
pode-se substituir (9) em (8) obtendo-se a seguinte expressão:
p(Xk−1 |Zk ) =
p(zk |xk−1 ) · p(xk−1 |xk−2 ) · p(Xk−2 |Zk−1 )
p(zk |Zk−1 )
p(zk |xk−1 ) · p(xk−1 |xk−2 ) · p(Xk−2 |Zk−1 ).
∼
=
∝
Mas:
wki
Esta escolha permite avançar o contador k no processo de filtragem simplesmente multiplicando-se (⊗)
uma nova realização das matrizes A e B por cada partícula xik−2 , i = 1, . . . , N e pelos valores conhecidos de uk−1 , expandindo assim a partícula Xik−2 para
Xik−1 . Esta etapa é similar, em filtragem de Kalman,
ao procedimento de previsão do estado.
A substituição da equação 11 na equação 10 resulta em:
∝
p(Xik−1 |Zk )
q(Xik−1 |Zk )
∝
p(zk |xik−1 ) · p(xik−1 |xik−2 ) · p(Xik−2 |Zk−1 )
q(xik−1 |Xik−2 , Zk ) · q(Xik−2 |Zk−1 )
∝
p(zk |xik−1 ) · p(xik−1 |xik−2 )
i
· wk−1
,
q(xik−1 |Xik−2 , Zk )
i
= wk−1
· p(zk |xik−1 ).
Esta equação é utilizada para atualizar os pesos
das partículas no momento em que a medida zk se
torna disponível. Em filtragem de Kalman isto equivale ao procedimento de atualização das previsões,
isto é, o procedimento de correção da estimativa feita
no procedimento de previsão, à luz de uma nova medida. Levando-se em conta que tanto a medida zk
quanto a partícula xik−1 e as variáveis de entrada uk
são conhecidas no momento deste cálculo, conclui-se
que a atualização deve ser feita por:
wki
i
= wk−1
· V (ηk , zk ),
=
i
wk−1
p(zk |xik−1 ) · p(xik−1 |xik−2 )
.
·
q(xik−1 |Xik−2 , Zk )
Admite-se ainda:
q(xik−1 |Xik−2 , Zk ) = q(xik |xik−2 , zk )
e finalmente:
wki
=
i
wk−1
·
p(zk |xik−1 ) · p(xik−1 |xik−2 )
.(10)
q(xik−1 |xik−2 , zk )
Desde que respeitada a condição de similaridade
citada no início desta seção, em princípio a função de
densidade de importância pode ser escolhida com liberdade. Entretanto, segundo (Ristic et al., 2004) esta
escolha é um aspecto crítico em filtros de partículas.
Mostra-se em (Doucet and Andrieu, 2000) que a escolha ótima é dada por:
qopt (xik−1 |xik−2 , zk ) = p(xik−1 |xik−2 , zk ).
Visando a simplicidade do algoritmo, neste trabalho adotar-se-á uma escolha sub-ótima para a densidade de importância dada por:
q(xik−1 |xik−2 , zk ) = p(xik−1 |xik−2 ). (11)
ISSN: 2175-8905 - Vol. X
(13)
sendo a função de verossimilhança V (η, z) deduzida
T
no apêndice (equação 15) e ηkT = [xiT
k−1 | uk ].
Observa-se que o procedimento de previsão baseado na equação 11 não depende da ocorrência da késima medida zk . Portanto nada impede que seja realizado logo após a atualização da estimativa de Xk−1 ,
resultando, conforme discutido anteriormente, na melhor estimativa possível para xk dadas as medidas em
k. A estimativa de xk é portanto dada por:
Portanto:
wki
(12)
x̂k =
N
X
xik · wki .
(14)
i=1
É usual na evolução de um algoritmo de filtragem
por partículas a observação de um fenômeno de colapso de probabilidades no qual ao longo das iterações, muitas partículas adquirem pesos muito baixos
e poucas assumem peso não desprezível. A ocorrência
deste fato leva a um mal-condicionamento da representação por partículas da densidade de probabilidade
de Xk . Para evitar este problema, quando se detecta
que o número efetivo de partículas que representa a
função de densidade de probabilidade é muito baixo,
utiliza-se um procedimento de reamostragem de partículas, que pode ser resumido como segue. Cada partícula é clonada um número de vezes proporcional ao
seu peso. Como o número total de partículas deve permanecer igual a N , resulta que as patículas com peso
muito baixo não terão clones, sendo assim abandonadas. Após o procedimento de clonagem, os pesos das
partículas sobreviventes são reajustados, assumindo o
valor 1/N . Este procedimento é formalmente descrito
em (Ristic et al., 2004).
Finalmente, é possível sintetizar o algoritmo proposto para filtragem por partículas de sistemas maxplus da seguinte forma:
1391
X SBAI – Simpósio Brasileiro de Automação Inteligente
18 a 21 de setembro de 2011
São João del-Rei - MG - Brasil
1. k=0;
2. Inicializar N partículas, X0i , i = 1, . . . , N ;
3. Para cada k:
Ler a medida zk ;
i
Atualizar os pesos de Xk−1
(eq. 13);
Se necessário, reamostrar;
Gerar as partículas para k (eq. 4);
Estimar xk (eq. 14);
4. Fim
Na próxima seção um exemplo permite avaliar o
desempenho do algoritmo de filtragem proposto.
4
Resultados
Consideremos o sistema descrito pela matriz A
abaixo, proposto em (Loreto et al., 2010). Supõe-se
que as entradas sejam nulas (sistema autônomo) e que
os estados de número 3, 6 e 8 sejam diretamente observados. Os elementos (2, 1), (5, 2), e (5, 4) da matriz
A são variáveis aleatórias uniformemente distribuídas
segundo os valores entre colchetes. Os outros valores
são determinísticos.


ε
ε
4
ε
ε ε 2 ε ε
[1, 7]
ε
ε
ε
ε ε ε 3 ε


 ε
5
ε
ε
ε ε ε ε 1


 4
ε
ε
ε
ε 3 ε ε ε


[3, 5] ε [1, 3] ε ε ε ε ε
A=

 ε
 ε
ε
5
ε
4 ε ε ε ε


 ε
ε
ε
4
ε ε ε ε 3


 ε
ε
ε
ε
3 ε 5 ε ε
ε
ε
ε
ε
ε 2 ε ε ε
A este sistema foi aplicado o filtro de partículas descrito anteriormente e comparado com o observador de
estado proposto em (Hardouin et al., 2010), caracterizado por determinar um limite inferior para cada variável de estado. As figuras 2 e 3 e as tabelas comparativas dos erros médio e máximo para cada estado são
apresentadas a seguir.
Figura 2: Trajetória da transição 2
ISSN: 2175-8905 - Vol. X
Figura 3: Trajetória da transição 5
Erro Médio
Est. Filt. Obs.
x1
0,00 0,00
x2
1,19 2,61
x3
0,00 0,00
x4
0,00 0,00
x5
0,58 3,17
x6
0,00 0,00
x7
0,00 0,00
x8
0,00 0,00
x9
0,00 0,00
Erro Máximo
Est. Filt. Obs.
x1
0,00 0,00
x2
2,44 5,85
x3
0,00 0,00
x4
0,00 0,00
x5
1,25 5,66
x6
0,00 0,00
x7
0,00 0,00
x8
0,00 0,00
x9
0,00 0,00
Observa-se que, exceto os estados 2 e 5, todos os
outros são observados e filtrados com precisão. As
figuras mostram que nos trechos selecionados as estimativas produzidos pelo filtro de partículas são em
geral mais próximas do estado real que o seu limiar
inferior. Este fato é confirmado pelas estimativas de
erro dos dois algoritmos mostradas nas tabelas.
5
Conclusão
Neste trabalho desenvolveu-se um algoritmo de filtragem estocástica para sistemas max-plus lineares baseado na técnica conhecida como filtro de partículas.
Mostrou-se uma maneira simples de realizar a amostragem por importância, baseada na geração de previsões a partir de partículas do estado no valor precedente de contagem. Mostrou-se também como atualizar os pesos das partículas a partir da obtenção de uma
nova medida, utilizando unicamente a função de verossimilhança obtida a partir da probabilidade da medida condicionada ao estado. Uma expressão para esta
função de verossimilhança foi deduzida para o caso
max-plus linear em que as perturbações estocásticas
estão associadas aos elementos das matrizes A e B do
modelo dinâmico do sistema, sendo variáveis aleatórias uniformemente distribuídas. Os resultados apresentados mostram que as técnicas de filtragem de partículas aplicadas a sistemas max-plus lineares, embora
envolvendo uma escolha sub-ótima da densidade de
1392
X SBAI – Simpósio Brasileiro de Automação Inteligente
18 a 21 de setembro de 2011
São João del-Rei - MG - Brasil
importância, levam a estimativas próximas aos valores
reais do estado do sistema. Estudos estão em curso
visando a utilização da amostragem ótima, que tem
como contrapartida maior complexidade tanto na fase
de previsão como na fase de atualização.
Agradecimentos
Este trabalho é parte das atividades do projeto CAPES/COFECUB 642/09 e contou com o financiamento do CNPq através do projeto PIBIC (Programa
Institucional de Bolsas de Iniciação Científica - Unicamp)
Referências
Baccelli, F., Cohen, G., Olsder, G. and Quadrat, J.
(1992). Synchronization and Linearity: An Algebra for Discrete Event Systems, John Wiley and
Sons, New York.
Bar-Shalom, Y., Li, X. and Kirubarajan, T. (2001). Estimation with Applications to Tracking and Navigation, Artech House, Norwood.
Brown, R. and Hwang, P. (1997). Introduction to Random Signals and Applied Kalman Filtering, John
Wiley and Sons, New York.
Cassandras, C. G. and Lafortune, S. (1999). Introduction to Discrete Event Systems, Kluwer Academic Publishers.
Doucet, A., G. S. and Andrieu, C. (2000). On sequential monte-carlo sampling methods for bayesian
filtering, Statistics and Computing 10: 197–208.
Hardouin, L., C.A.Maia, Cottenceau, B. and Lhommeau, M. (2010). Observer design for (max,+)
linear systems, IEEE Trans. on Automatic Control 55-2: 538–543.
Loreto, M. D., Gaubert, S., Katz, R. D. and Loiseau,
J. (2010). Duality between invariant spaces for
max-plus linear discrete event systems, SIAM J.
on Control and Optimaztion 48-8: 5606–5628.
Maia, C., Santos-Mendes, R., Luders, R. and Hardouin, L. (2005). Estratégias de controle por modelo de referência de sistemas a eventos discretos
max-plus lineares, SBA 16(3): 263–278.
Murata, T. (1989). Petri nets : properties, analysis and applications., Proceedings of the IEEE
77(4): 541–580.
Ristic, B., Arulampalam, S. and Gordon, N. (2004).
Beyond the Kalman Filter - Particle Filters for
Tracking Applications, Artech House.
ISSN: 2175-8905 - Vol. X
Apêndice: Cálculo da função de verossimilhança
Consideremos a equação max-plus z = A00 ⊗ x ⊕
B ⊗ u, sendo x ∈ Rn×1 , z ∈ Rq×1 , u ∈ Rp×1 ,
A00 ∈ Rq×n e B ∈ Rq×p . Definindo η T = [xT | uT ],
e H = [A00 | B], tem-se z = H ⊗ η, sendo η ∈ Rr×1 ,
H ∈ Rq×r e r = n + p. Supõe-se que os elementos
da matriz H, notados por hij , sejam variáveis aleatórias independentes e uniformemente distribuídas entre
hij e hij , sendo suas funções de probabilidade acumulada e de densidade de probabilidade respectivamente
dadas por:


0
se τ ≤ hij


 τ −hij
se hij < τ ≤ hij
Fij (τ ) =
hij −hij




1
se τ > hij

1

se hij < τ ≤ hij
hij −hij
pij (τ ) =
 0
se τ ≤ hij ou τ > hij
Determina-se a seguir a densidade de probabilidade da
variável aleatória z, condicionada pelo vetor η, isto é,
pz (t|η), sendo t = [t1 . . . tq ]T ∈ Rq×1 . Em razão da
independência dos elementos de H tem-se que:
P [z ≤ t] = P [z1 ≤ t1 e . . . e zq ≤ tq ] =
q
Y
P [zi ≤ ti ]
i=1
Porém:
P [zi ≤ ti ]
= P [max(hij + ηj ) ≤ ti ]
j
= P [hi1 ≤ ti − η1 e . . . e hir ≤ ti − ηr ]
r
Y
Fij (ti − ηj )
=
j=1
Portanto a função de probabilidade acumulada conjunta do vetor z, condicionada por η, é dada por:
Fz (t|η)
=
q Y
r
Y
Fij (ti − ηj )
i=1 j=1
Derivando-se sucessivamente em relação a t1 . . . tq
obtém-se a densidade de probabilidade procurada:
pz (t|η) =
=
q
r
Y
∂ Y
(
Fij (ti − ηj ))
∂ti j=1
i=1
q X
r
r
Y
Y
∂
(
(Fij (ti − ηj ))
Fik (ti − ηk ))
∂ti
i=1 j=1
k6=j
q X
r
r
Y
Y
=
(
pij (ti − ηj )
Fik (ti − ηk ))
i=1 j=1
k6=j
Se o vetor z for conhecido, a função acima é chamada de função de verossimilhança de η. Neste caso,
substituindo-se t pelo valor conhecido de z, obtém-se:
q X
r
r
Y
Y
V (η, z) =
(
pij (zi − ηj )
Fik (zi − ηk ))
i=1 j=1
k6=j
(15)
1393