Controle linear (estado-espaço)
Controle linear (estado-espaço)
No artigo sobre controladores PID, tentamos mostrar como funciona o controle "analógico", ou "clássico", por feedback. Um simples aquecedor de água foi nossa cobaia ao longo de todo o texto. Uma leitura cuidadosa daquele artigo permite tirar algumas conclusões.
O controle por feedback é heurístico: funciona mesmo quando sabemos muito pouco a respeito do sistema a ser controlado. Porém é preciso testar o controle de alguma maneira, principalmente no caso de controladores PI e PID que dependem de calibragem para funcionar bem.
Uma forma de testar o controle é desenvolvendo um simulador, e neste caso precisamos conhecer bem o sistema para que o simulador fique bom. Outra opção é testar empiricamente, direto no sistema ou máquina "real". Isto nega parte da vantagem inicial, que era justamente não precisar conhecer o sistema.
Outra limitação importante do controle PID é: apenas uma entrada e apenas uma saída. No caso do aquecedor de exemplo, a entrada era a temperatura da água (na verdade, a diferença entre a temperatura desejada e a real), e a saída era o volume de gás liberado para queima.
Para controlar duas saídas (por exemplo, se o aquecedor possuir um ventilador além do queimador) precisaríamos de dois controles PID e duas entradas independentes. Cada um dos controles comporta-se como se o outro controle, e a outra entrada, não existissem. (O ventilador do aquecedor poderia ser acionado de acordo com a temperatura do ar.)
Mas isto pode não ser suficiente; talvez as duas saídas dependam das duas entradas simultaneamente. A velocidade do ventilador do aquecedor deveria variar também conforme a quantidade de gás queimado, porque sua função não é só esfriar o aquecedor, mas também remover o CO2.
Ou talvez haja muito mais entradas que saídas. O aquecedor poderia ser equipado com sensor de CO2, sensor de fluxo, sensor de temperatura de entrada da água, etc. Mas não podemos aproveitar estas novas entradas para melhorar o desempenho de nosso controle PID.
É claro, sempre podemos fazer "gambiarras" no controle, para escapar das limitações do PID. Isto é fácil se ele for implementado em software. Porém estas gambiarras ficam totalmente de fora do modelo teórico. Isto dificulta muito a prova rigorosa de que a gambiarra funcionará bem, e cada desenvolvedor implementará a gambiarra de forma diferente.
Finalmente, a análise do controle PID é difícil, devido ao uso de transformadas de Laplace. Manipular fórmulas é complicado, e os métodos numéricos (ou seja, por força bruta) são limitados. Toda essa fundamentação teórica baseia-se em funções contínuas (analógicas) porém a maioria dos controladores é digital, só entende funções discretas.
A representação por estado-espaço aborda o problema do controle de forma completamente diferente.
Uma curiosidade história: a representação estado-espaço foi estudada em primeiro lugar na antiga União Soviética, e foi um dos fatores que a permitiu largar na frente na Corrida Espacial. A despeito da Guerra Fria, os cientistas russos e ocidentais compartilharam informação e desenvolveram coletivamente a teoria, que era uma necessidade em comum para a indústria aeroespacial.
Em primeiro lugar, é necessário modelar o sistema de forma muito fiel à coisa real. Para fazer isso é preciso conhecer o sistema a fundo, inclusive as leis da Física que o afetam. Compare isso com o controle por feedback, cuja maior vantagem era justamente poder ignorar o funcionamento interno do sistema.
Por outro lado, é importante lembrar que no artigo sobre controle PID nós tivemos de desenvolver um simulador em separado. Isto não será mais necessário aqui. Uma vez que as fórmulas estejam desenvolvidas, a simulação está disponível de graça.
Em segundo lugar, o estado do sistema é diretamente expresso por uma equação diferencial, no domínio do tempo. (Nos sistemas com controle PID, o sistema é expresso por funções de transferência, no domínio da freqüência.)
A forma canônica da representação estado-espaço de qualquer sistema é a seguinte:
x' = Ax + Bu y = Cx + Du x: estado do sistema (matriz S x 1) S: número de estados u: controle ou entrada do sistema (matriz E x 1) E: número de entradas ou controles y: saída do sistema (matriz T x 1) T: número de saídas A: matriz S x S B: matriz S x E C: matriz T x S D: matriz T x E
Pois é, as equações acima envolvem matrizes, o que é um incômodo inicial para todos nós que não prestamos muita atenção às aulas de álgebra linear no segundo grau...
As matrizes são utilizadas porque representam um "pacote" de números de forma muito compacta. No modelo acima, pode haver qualquer número de entradas, estados e saídas no sistema, e todos podem ser interdependentes — um notório avanço em relação ao controle PID.
Antes de falar qualquer coisa a mais, vamos modelar um sistema simples, porém completo: o aquecedor de água que foi nossa cobaia para o controle PID.
O estado do sistema (matriz x) é o conjunto de variáveis que podemos medir e queremos controlar. No caso do aquecedor, o único estado é a temperatura de saída da água (qout) em Celsius. Por simplicidade, consideramos esta ser a temperatura de toda a água acumulada dentro do aquecedor.
A entrada ou controle do sistema (matriz u) é o conjunto de variáveis que podem influir no estado. No nosso aquecedor, são três: temperatura de entrada da água (qin), temperatura desejada pelo usuário (qtgt), e a potência do aquecedor (w) em kcal/min.
Pode ser difícil decidir se algo é estado (x), controle (u), entrada (u), saída (y) ou característica do sistema (A ou B). Temperatura de entrada não pode ser controlada, mas certamente influencia a temperatura de saída. A potência do aquecedor atua sobre o estado (pois aquece a água), mas também é uma variável que desejamos controlar. Afinal, potência é estado ou controle? Na minha opinião, encaixou melhor como controle.
A saída do sistema é uma transformação simples das variáveis para o mundo físico. Por exemplo, o controle de potência (w) em kcal/min precisa ser traduzido para kBTU/h ou Watts. Para facilitar nosso exemplo, vamos presumir que se trata de um aquecedor elétrico e usaremos Watts. A equação de saída não é diferencial, e seu resultado não influencia o estado.
Mais algumas grandezas que precisamos encaixar no modelo:
Os itens acima têm de ser encaixados nas matrizes A e B, pois são características do sistema; não são controles nem estados. "Ah, mas a vazão vai mudar ao longo do tempo!" Verdade, mas todos os elementos da equação de estado-espaço também podem mudar com o tempo.
Os elementos x, u, A, B, C e D são na verdade funções no tempo: x(t), u(t), A(t), B(t), C(t) e D(t). Então não há problema em colocar características variáveis dentro das matrizes A a D. Existem sistemas onde as matrizes A a D são de fato constantes, denominados LTI (linear time-invariant), mas não é o caso do nosso aquecedor.
Acreditando que descrevemos o modelo suficientemente bem, seguem as equações, uma por uma.
x' = A x + B u
q' = [ -(l+f/z) ] q + [ l+f/z, 0, 1/z ] [ qin ]
[ qtgt ]
[ w ]
q' Variação (derivada) da temperatura, Celsius/min
q Temperatura da água, Celsius
qin Temperatura de entrada da água, Celsius
qtgt Temperatura desejada da água, Celsius
l Perda térmica, % por grau C de diferença
f Vazão de saída de água quente, litros/min
z Volume de água quente no aquecedor, litros
w Potência do aquecedor, kcal/min
A equação acima simula o aquecedor de forma muito próxima à realidade. A temperatura sobe em função da potência do aquecedor, e diminui (tendendo a qin) conforme água quente é consumida, além da perda pelo esfriamento natural. A temperatura desejada pelo usuário (qtgt) ainda não influi no estado, mas resolveremos isto depois.
A equação ainda é contínua. É preciso "discretizá-la" para permitir seu uso numa planilha ou programa. Há duas formas de fazer isto:
No nosso aquecedor, a unidade de tempo é o minuto, então se fizermos dt=0.1, estamos simulando o estado a cada seis segundos, o que é precisão suficiente. É importante lembrar de escalar (multiplicar) a derivada por dt — se a derivada do estado resultou 5, devemos adicionar 5x0.1=0.5 ao estado, pois 5 é a taxa de variação para um minuto inteiro.
Conhecer o volume de água dentro do aquecedor permite fazer uma previsão muito boa do "atraso" entre o aquecedor ligar e a água atingir a temperatura desejada. A simulação fica fiel porque a relação entre volume e vazão é perfeitamente conhecida. No controle PID, nós tínhamos utilizado uma simples média móvel.
Agora, a equação de saída:
y = C x + D u
y = 0 q + [ 0 0 4186 ] [ qin ]
[ qtgt ]
[ w ]
Conforme dito antes, a equação de saída é uma mera transdução de calorias para Watts. O resultado y comanda diretamente o aquecedor elétrico.
O único defeito da modelagem acima é que não estamos controlando a potência do aquecedor (w). O modelo está ok para uma potência fixa, ou se o usuário controlasse manualmente a potência, mas obviamente desejamos que isto seja automático.
Vamos apelar ao tradicional controle proporcional por feedback (P) para uma solução simples e ilustrativa. O feedback pode ser implementado recalculando-se "u" a cada passo, conforme a seguinte equação:
u = r - Kx u Matriz com valores de controle r Matriz com valores desejados K Matriz do controlador de feedback
A matriz "r" pode mudar com o tempo, mas ela não reage em função do estado do sistema; ela se altera conforme o usuário muda a temperatura desejada e conforme a temperatura ambiente oscila.
Para modelar um aquecedor funcional, falta então especificarmos as matrizes "r" e "K". É o que faremos a seguir.
u = r - K x
u = [ qin ] - [ 0 ] q
[ qtgt ] [ 0 ]
[ qtgt . P ] [ P ]
P: ganho do controlador P
Ou seja,
u = [ qin ]
[ qtgt ]
[ P.(qtgt - q) ]
A cada passo, nós recalculamos a potência do aquecedor, que torna-se proporcional à diferença entre a temperatura atual (q) e a desejada (qtgt). Faltaria ainda limitar a potência a valores positivos, já que nosso aquecedor não pode refrigerar, mas são detalhes.
É interessante notar que precisamos estabelecer um valor "artificial" para a terceira linha da matriz r, de modo que a conta feche. Esta é uma característica das equações de estado-espaço — às vezes precisamos manipular valores para encaixar nosso sistema dentro da forma canônica.
Este controlador P tem as mesmas deficiências que já estudamos no artigo sobre controle PID. A única melhoria real é que nosso modelo funciona com fluxo zero; o aquecedor esquenta a água contida dentro dele, e desliga quando atinge a temperatura desejada, tal qual uma caldeira.
A seguinte manipulação de fórmulas é muito interessante:
x' = Ax + Bu u = r - Kx x' = Ax + B(r - Kx) x' = Ax + Br - BKx x' = (A - BK)x + Br
A forma final omite completamente a matriz de controle "u". Determinamos o comportamento do sistema unicamente em função da configuração "r". Esta forma é muito útil para determinar se o sistema controlado por feedback é estável. Mas, na implementação real, ainda precisamos calcular "u" porque é nele que está o controle da potência.
Agora podemos usar as equações diretamente num simulador! Pressione qualquer botão para iniciar a simulação.
|
Vazão |
0.0 L/min |
Temperatura de saída --- °C Potência --- W |
|
Temp. água fria
|
20.0°C | |
| Temp. alvo | 40.0°C | |
| Peso P | ||
| Peso I | ||
Os pesos P e I podem ser manipulados. Ainda não abordamos a implementação do controle PI em estado-espaço; isto será explicado mais à frente. Para experimentar com o controle PI, digite um peso I diferente de zero.
O código-fonte pode ser consultado na página. Utilizamos a biblioteca Sylvester para fazer as operações com matrizes.
Sistemas descritos pelas equações de espaço-estado são denominados sistemas lineares, porque o uso de matrizes impõe uma relação linear entre controles, estados e saídas.
A totalidade das máquinas e dispositivos da vida real é não-linear. Nem o aquecimento da água é linear; a definição de caloria vale apenas à temperatura de 19,5C. Mas o desvio é suficientemente pequeno para ser ignorado. Em outros sistemas, a coisa pode não ser tão simples.
Em geral, é mais desejável modelar e controlar o sistema como se ele fosse linear, então existem diversas técnicas de linearização para compatibilizar as variáveis não-lineares.
Podemos dizer muito sobre um sistema analisando diretamente as matrizes A, B, C e D, sem precisar apelar à transformada de Laplace.
Esta é a vantagem de adotar a forma "canônica" da equação de estado, em vez de outra qualquer. Compêndios já foram escritos sobre o assunto, e a reutilização das técnicas de análise fica infinitamente mais fácil se todo mundo usa a mesma forma.
Um sistema é considerado estável quando todos os autovalores da matriz A são negativos. (Autovalores de uma matriz são aparentados com o determinante, e também só existem para matrizes quadradas. Uma matriz 3x3 tem 3 autovalores, enquanto o determinante é sempre único. Autovalor e determinante de uma matriz 1x1 são iguais a seu elemento.)
No caso do aquecedor, a matriz A é 1x1 e igual a -(l+f/z). Como l, f e z são sempre positivos, o autovalor é sempre negativo e o sistema é estável.
Não é difícil enxergar porque isto é verdade. Supondo a equação diferencial
x' = Ax
a derivada de x só tende a zero — e x só tende a um valor finito — se A for negativo.
Esquecemos um detalhe: nosso aquecedor é controlado por feedback! Então não faz sentido analisar a matriz "A" isoladamente, porque, na presença de feedback,
x' = (A - BK)x + Br
Precisamos verificar se A-BK é igualmente negativa.
A - B K
-(l+f/z) - [ l+f/z, 0, 1/z ] [ 0 ]
[ 0 ]
[ P ]
-(l+f/z) - P/z
Felizmente, não há problemas de estabilidade com o feedback, porque o ganho P do controlador também é sempre positivo.
O feedback é capaz inclusive de estabilizar um sistema originalmente instável, porque a composição A-BK toma o lugar da matriz A.
Um sistema é considerado controlável quando é possível conduzir o(s) estado(s) para qualquer valor, num tempo finito, dado um conjunto de valores de controle (matriz u). Nosso aquecedor é controlável porque basta manipular a temperatura da água fria e a potência para atingir qualquer temperatura final. A controlabilidade é facilmente testada com base nas matrizes A e B.
Um sistema com mais estados que controles pode ser controlável; depende do sistema. Por exemplo, uma locomotiva sobre trilhos tem no mínimo dois estados importantes (posição e velocidade) e apenas um controle (aceleração). É possível atingir qualquer combinação de posição e velocidade variando apenas a aceleração, como qualquer maquinista sabe. Já um avião não poderia ser controlado apenas com o acelerador (embora às vezes seja preciso tentar!).
Um sistema é observável se podemos determinar o estado atual do sistema com base apenas nas saídas e nas entradas/controles. Nosso aquecedor é observável desde que conheçamos a vazão, usando alguma espécie de medidor. Por outro lado, se medirmos a temperatura de saída, podemos calcular a vazão sem medidor (e um medidor de temperatura é mais barato).
Nessas alturas do texto, você talvez esteja pensando: controlar o aquecedor usando PID era tão mais fácil! Bastava medir a temperatura de saída. Já aqui, também precisamos saber vazão e temperatura de entrada, o que parece conduzir a um sistema bem mais caro.
(Ao menos, o volume do aquecedor e o índice de dissipação são valores fixos, fáceis de determinar, então a presença destes não encarece o controlador.)
A boa notícia é que não precisamos realmente implementar as equações completas no do aquecedor. O estado (x) vai ser conhecido o tempo todo, pois usaremos um sensor de temperatura. Só precisamos das equações para conferir a qualidade do controle.
Basta implementar as equações de feedback e de saída, que também podem ser "enxugadas" ainda mais. (A equação de feedback ainda é necessária porque é através dela que determinamos a potência do aquecedor.)
u = [ qin ] - [ 0 ] q
[ qtgt ] [ 0 ]
[ qtgt . P ] [ P ]
u = P (qtgt - q)
y = 0 q + [ 0 0 4186 ] [ qin ]
[ qtgt ]
[ w ]
y = 4186 . u
Como pode ser visto acima, o sistema ficou reduzido a um controlador P, que baseia-se unicamente na temperatura escolhida pelo usuário qtgt e a temperatura atual q. A equação de saída é uma mera conversão de unidade.
É claro, esta é a forma mais barata possível de implementar o sistema, não necessariamente a melhor. Poderíamos manter todas as equações, e comparar a temperatura estimada com a temperatura no sensor — e fazer alguma coisa útil com a diferença. Voltaremos a esta possibilidade mais ao final do texto.
Para implementar um controlador no estilo PI (proporcional-integral), que não apresenta o "droop" ou queda de rendimento do aquecedor, precisamos encaixar a integral do erro nas equações.
A forma mais apropriada é considerar esta integral como parte do estado do sistema.
x' = A x + B u
[ q' ] = [ -(l+f/z), 0 ] [ q ] + [ l+f/z, 0, 1/z ] [ qin ]
[ e' ] [ -1 , 0 ] [ e ] [ 0 , 1, 0 ] [ qtgt ]
[ w ]
e: integral do erro
= somatório de (qtgt - q)
Na equação acima, podemos ver que a variação do erro é a diferença instantânea entre q e qtgt, portanto o erro em si é cumulativo.
A equação de saída do nosso aquecedor não muda pois a potência de saída é função do controle (u), que não mudou. As equações de feedback têm de mudar já que precisamos dar utilidade ao erro integral:
u = r - K x
u = [ qin ] - [ 0, 0 ] [ q ]
[ qtgt ] [ 0, 0 ] [ e ]
[ qtgt . P ] [ P, -I ]
I: ganho integral do controlador PI
u = [ qin ]
[ qtgt ]
[ P.(qtgt - q) + I.e ]
O sistema PI ficou quase tão simples quanto o P, e podemos igualmente "enxugá-lo" nos mesmos moldes que fizemos com o P, se quisermos barateá-lo.
As matrizes permitem expressar diversas equações paralelas e independentes, então sistemas mais complexos, descritos por múltiplas equações diferenciais, continuam cabendo na forma. Basta fazer "x" uma matriz de diversas linhas, e ajustar as demais matrizes de acordo.
A forma usa uma equação diferencial de primeiro grau (derivada simples). Se o sistema tem derivada dupla, tripla, etc. é preciso fazer alguma manipulação para encaixá-las na forma. Ou usamos variáveis substitutas, ou adicionamos as derivadas intermediárias no estado, formando uma cadeia.
Por exemplo, num sistema onde queremos controlar a posição de um objeto fazendo uso de força (aceleração), a aceleração é a segunda derivada da posição. Precisamos incluir a velocidade, que é a derivada simples da posição:
x' = A x + B y
| velocidade | = | 0 1 | | posição | + | 0 | força
| aceleração | | 0 0 | | velocidade | | 1/massa |
velocidade = velocidade
aceleração = força / massa
Ao atualizar x com base em x':
posição = posição + velocidade
velocidade = velocidade + aceleração
Conforme dito antes, não precisamos abolir a equação de estado-espaço em nossa implementação de aquecedor. Podemos usá-la em conjunto com um sensor de temperatura para melhorar o sistema.
O estado calculado pelas equações não possui ruído, porém é apenas uma estimativa, e vai desviar do estado real. É como um barco navegando às cegas, estimando sua posição com base na velocidade e rumo. O cálculo é exato mas vai descolando cada vez mais da realidade, conforme passa o tempo.
Já o estado obtido por um sensor é "a coisa real", mas então haverá o ruído de medição, e outras limitações. Por exemplo, o GPS embute um erro aleatório, e pode demorar até 10 minutos para fornecer a posição.
Num sistema ideal, podemos combinar os dois estados (estimado e medido) de modo a obter um estado de qualidade mais alta: mais exato, com menos ruído e com o mínimo possível de atraso.
Para combinar os dois valores, podemos usar uma média, ou média móvel, mas médias simples não tiram proveito dos pontos fortes de cada medida (a ausência de ruído da estimativa, e a exatidão do sensor). A ferramenta ideal é o filtro de Kalman.
Em sua essência, o filtro de Kalman é uma média móvel, mas o ganho desta média móvel é regulado dinamicamente, de modo a aumentar a "qualidade" do estado. A fórmula básica é:
x[t+dt] = K[t].z[t] + (1 - K[t]).x[t] x[t+dt] nova estimativa do estado K[t] ganho de Kalman (entre 0 e 1) no instante t z[t] estado medido (e.g. com sensor) no instante t x[t] estimativa anterior do estado
O filtro de Kalman é naturalmente discreto, e seu mecanismo recursivo usa pouca memória num computador.
Para calcular o ganho de Kalman a cada passo, com alguma simplificação em relação à fórmula original:
Pm = Ad.P.Adt + Q K = Pm / (Pm + R) (calcular o novo estado neste ponto!) P = (I - K).Pm P: covariância estimada do sistema, no momento Q: covariância da estimativa R: covariância da medida (e.g. sensor) K: ganho de Kalman Ad: matriz A da equação de estado em forma integral Adt: Ad transposta I: Matriz identidade
As fórmulas intimidam um pouco mas vamos ilustrar com um exemplo, usando de novo o aquecedor com controle P como cobaia. Como seu estado é uma matriz 1x1, todas as matrizes de Kalman viram simples números.
Vamos estimar algumas características do nosso aquecedor para testar o filtro de Kalman:
Desvio-padrão da estimativa: 0.4C Covariância = desvio-padrão ao quadrado Q = 0.4 x 0.4 = 0.16 Desvio-padrão do sensor de temperatura: 2C R = 2 x 2 = 4 Vazão de água: zero Perda térmica (l): 1% por minuto por grau A = -(l+f/z) = -0.01 Taxa de amostragem: 6 segundos (dt=0.1 minuto) Ad ~= 1 + A.dt = 1 - (0.01).(0.1) = 0.999
A covariância da estimativa e do sensor significam coisas bem diferentes. A covariância do sensor é uma medida de ruído mesmo, devido a imperfeições do sensor. Já a covariância da estimativa mede o quanto ela pode se movimentar. A movimentação torna a estimativa "menos confiável", não porque podemos errar alguma conta, mas sim porque o modelo não é perfeito.
Tivemos de fazer uma estimativa de vazão e perda térmica porque, no nosso modelo de aquecedor, estes valores influem no valor da matriz A, que é utilizada (na forma integral) no filtro de Kalman. Para calcular a versão integral de A, utilizamos a aproximação 1+A.dt em vez do exponencial, que é aceitável quando dt é relativamente pequeno (ou seja, a taxa de amostragem é relativamente alta).
Mais algumas variáveis iniciais:
Temperatura real da água: 40C Temperatura estimada: 40.5C P inicial: 0.5
Também vamos ignorar a dinâmica do controle por feedback, embora ele certamente reagiria quando o filtro de Kalman atualizasse o estado.
Primeira rodada:
Pm = Ad.P.Adt + Q Pm = 0.999 * 0.5 * 0.999 + 0.16 Pm = 0.659 K = Pm / (Pm + R) K = 0.659 / (0.659 + 4) K = 0.1414 z = 42.5047 (leitura do sensor-termômetro) x = K.z + (1 - K).x x = 0.1414 * 42.5047 + (1 - 0.1414) * 40.5 x = 40.7834 (novo estado estimado) P = (1 - K).Pm P = (1 - 0.1414) * 0.659 P = 0.5658
Segunda rodada:
Pm = Ad.P.Adt + Q Pm = 0.999 * 0.5658 * 0.999 + 0.16 Pm = 0.7246 K = Pm / (Pm + R) K = 0.7246 / (0.7246 + 4) K = 0.15336 z = 39.7925 (leitura do sensor-termômetro) x = K.z + (1 - K).x x = 0.15336 * 39.7925 + (1 - 0.15336) * 40.7834 x = 40.6314 (novo estado estimado) P = (1 - K).Pm P = (1 - 0.15336) * 0.7246 P = 0.6134
Terceira rodada:
Pm = Ad.P.Adt + Q Pm = 0.999 * 0.6134 * 0.999 + 0.16 Pm = 0.772 K = Pm / (Pm + R) K = 0.772 / (0.772 + 4) K = 0.16177 z = 38.8512 (leitura do sensor-termômetro) x = K.z + (1 - K).x x = 0.16177 * 38.8512 + (1 - 0.16177) * 40.6314 x = 40.3434 (novo estado estimado) P = (1 - K).Pm P = (1 - 0.16177) * 0.772 P = 0.647
Devido ao ruído do sensor, a primeira rodada produziu uma estimativa ainda pior que a original, porém a terceira rodada já produziu um estado de qualidade melhor que a original.
Mais rodadas, simuladas com um script Python:
t=0, Pm=0.659001, K=0.141447, z=42.504700, x=40.783558, P=0.565787 t=1, Pm=0.724656, K=0.153378, z=39.792500, x=40.631552, P=0.613510 t=2, Pm=0.772284, K=0.161827, z=38.851200, x=40.343443, P=0.647307 t=3, Pm=0.806013, K=0.167709, z=38.981818, x=40.115086, P=0.670837 t=4, Pm=0.829496, K=0.171756, z=39.364760, x=39.986213, P=0.687025 t=5, Pm=0.845652, K=0.174518, z=38.033263, x=39.645389, P=0.698071 t=6, Pm=0.856675, K=0.176391, z=39.438286, x=39.608857, P=0.705565 t=7, Pm=0.864155, K=0.177658, z=42.294169, x=40.085924, P=0.710631 t=8, Pm=0.869210, K=0.178512, z=40.801745, x=40.213706, P=0.714046 t=9, Pm=0.872619, K=0.179086, z=39.490366, x=40.084166, P=0.716345
Devido à simplicidade do sistema, o ganho de Kalman tende logo a seu valor final (em torno de 0.18) e comporta-se de forma idêntica a uma média móvel com peso 0.18. Mas pelo menos fomos poupados de estimar este peso.
A maior vantagem, que este exemplo não mostra, é que o filtro de Kalman continua funcionando mesmo num sistema altamente dinâmico, onde todas as variáveis (inclusive as taxas de ruído) mudam o tempo todo. É o tipo de recurso que coloca o controle linear num patamar superior aos controles PID "analógicos".
O filtro de Kalman é o "pão com manteiga" de um receptor GPS, cuja leitura é muito variável (ou seja, muito ruído). O filtro permite obter uma leitura estável, e mais precisa conforme passa o tempo. Unidades GPS mais sofisticadas usam velocidade e rumo para estimar a mudança de posição, melhorando ainda mais a qualidade da posição reportada, mesmo quando o aparelho passa por lugares com perda de sinal como túneis ou ruas cercadas de prédios altos.
Como dito antes, o filtro de Kalman é apenas uma média móvel que determina seu próprio peso ou ganho. Podemos aplicá-lo em situações onde não há controle, apenas observação.
Por exemplo, suponha um termômetro de temperatura ambiente, cujo erro de medição seja conhecido. Uma média móvel ajuda a obter amostras de melhor qualidade. Mas qual seria o peso ideal dessa média móvel? Em vez de ficar "chutando" valores até chegar no peso ideal, podemos usar Kalman.
Relembrando a fórmula:
Pm = Ad.P.Adt + Q P: covariância estimada do sistema, no momento Ad: matriz A da equação de estado em forma integral Q: covariância da estimativa
Trata-se de um sistema estável — a temperatura ambiente daqui um minuto tende a ser a mesma que a de agora. Mesmo que modelássemos as variações diárias e sazonais, o fato é que elas são bastante lentas. Portanto, a matriz A de uma equação estado-espaço tende a zero. Usando a aproximação Ad=I+A.dt para pequenos valores de dt, Ad=I+0=I.
Neste exemplo, temos apenas um sensor e uma métrica a estimar, então todos os elementos são matrizes 1x1, equivalentes a números comuns.
O valor inicial de P é a covariância inicial do sistema. Ela reflete a incerteza da nossa estimativa inicial da temperatura.
Somos obrigados a começar com uma estimativa qualquer, nem que seja zero. Mas é melhor começar com algo mais provável, entre 20 e 25 graus. Obviamente, o erro dessa estimativa *vai* ser grande em relação à realidade. Usando um desvio-padrão de 10 graus, temos P=10x10=100.
Q é a covariância da estimativa. Parece o mesmo que P, mas não é. Q expressa a incerteza da próxima estimativa, supondo que nossa estimativa atual fosse perfeita.
Temperatura ambiente muda muito devagar. Vamos estimar que o desvio-padrão da temperatura seja 1 grau por hora, e vamos coletar uma amostra por minuto, então
desvio por minuto = 1 / 60 = 0.016666 variância = 0.016666 * 0.016666 = 0.0002777 P(inicial) = 100 Q = 0.0002777 Ad = 1 Pm = 1 x 0.0002777 x 1 + 0.0002777 = 100.0002777
Vamos ainda supor que o desvio-padrão do sensor-termômetro seja de 2 graus. Com isto, temos tudo que precisamos para terminar o cálculo da primeira rodada.
R = 2 x 2 = 4 K = Pm / (Pm + R) P = (I - K).Pm R: covariância da medida (e.g. sensor) K: ganho de Kalman K = 100.0002777 / (100.0002777 + 4) = 0.961538 (Recalcule a média móvel da temperatura neste ponto.) P = (1 - 0.961538) * 100.0002777 = 3.84615
O ganho de Kalman K é enorme nessa primeira rodada: 96%. Isto significa que a nova média móvel será quase toda baseada no valor obtido diretamente do sensor.
O novo valor de P é a nova incerteza da estimativa, que diminuiu bastante com apenas uma leitura do sensor (de 100 para 3.84). Conforme as rodadas seguemm, os valores convergem:
n 1 Pm 100.000278 K 0.961539 P 3.846154 n 2 Pm 3.846432 K 0.490214 P 1.960857 n 3 Pm 1.961134 K 0.328987 P 1.315947 n 4 Pm 1.316225 K 0.247586 P 0.990346 n 5 Pm 0.990623 K 0.198497 P 0.793988 n 6 Pm 0.794265 K 0.165670 P 0.662680 n 7 Pm 0.662957 K 0.142175 P 0.568701 n 8 Pm 0.568979 K 0.124531 P 0.498124 ... n 426 Pm 0.033528 K 0.008312 P 0.033250 n 427 Pm 0.033528 K 0.008312 P 0.033249 n 428 Pm 0.033527 K 0.008312 P 0.033248 n 429 Pm 0.033526 K 0.008312 P 0.033247
Depois de centenas de interações, o valor de K estaciona em torno de 0.0083, ou seja, a média móvel atribui peso de 0.8% para o sensor e 99.2% para a estimativa anterior (lembrando que o recálculo ocorre a cada minuto).
A pergunta é a seguinte: uma vez que se trata de um sistema estável, podemos pular a implementação do filtro de Kalman e usar apenas uma média móvel com peso constante 0.0082?
Não, não podemos, pois ainda existe o problema da estimativa inicial ser ruim. Um peso de 0.8% fixo só funciona quando a estimativa inicial é muito boa. Neste ponto, um amador poderia emendar a média móvel com algumas gambiarras, do tipo usar a primeira leitura do sensor como estimativa inicial. Mas aí o resultado final seria um código a) mais complicado que o filtro de Kalman e b) com qualidade menor que o filtro de Kalman.
Então, vale a pena fazer um esforço para entender e implementar as fórmulas, que não são nem tão complicadas assim.
Por último, vale a pena citar a questão da fusão de sensores. Por exemplo, utilizar dois sensores de temperatura com o objetivo de obter mais precisão e acurácia.
O filtro de Kalman "baunilha" não trata esse caso diretamente. Ele espera sensores que meçam grandezas diferentes, cuja relação pode ser descrita numa equação estado-espaço. Por exemplo, podemos parear um sensor de posição e um de aceleração para estimar velocidade.
Existem diversas alternativas para lidar com o caso da fusão de sensores. A mais fácil é usar alternadamente um ou outro sensor, a cada rodada, lembrando de usar a variância R correspondente a cada sensor (supondo que sejam diferentes).
Outra alternativa é utilizar o filtro estendido de Kalman. Este artigo possui um exemplo completo, que vamos resumir aqui. O filtro de Kalman adiciona um fator adicional C nas fórmulas para lidar com sinais não-lineares, e vamos utilizá-lo para fazer a fusão dos sensores.
Pm = Ad.P.Adt + Q K = (Pm.Ct) / (C.Pm.Ct + R) (*) x = K.z + (1 - K.C).x P = (I - K.C).Pm (*) divisão A/B definida como produto de A pelo inverso de B
Vamos considerar brevemente o caso de estimativa de uma única grandeza (temperatura) usando dois sensores. Neste caso, nas fórmulas acima, o fator C é uma matriz 2x1 (2 linhas, uma coluna) que representa os pesos dos dois sensores — normalmente, ambos iguais a 1. Essa matriz e sua versão transposta Ct (1x2) também fazem a adaptação entre matrizes de diferentes tamanhos.
A covariância R tem de ser uma matriz 2x2, com as variâncias na diagonal principal, e a covariância entre sensores (normalmente zeros) na diagonal secundária. As leituras dos sensores são acomodadas em z(n) na forma de uma matriz 2x1.
Uma vez que a estimativa continua sendo um número singelo (temperatura), x(n), P e Pm continuam sendo matrizes 1x1, equivalentes a números comuns. O ganho K tem tamanho 1x2. As regrinhas de multiplicação de matrizes garantem que tudo se encaixa no final.
A característica mais marcante de um sistema oscilatório, tipo massa-mola, pêndulo, etc. é uma equação diferencial de segunda ordem, cuja solução tem de envolver trigonometria. Vide a equação elementar abaixo:
v'' = -v v = cos(t)
A equação estado-espaço para esta cossenóide, na ausência de qualquer controle, é:
x' = A x | v' | = | 0.0, 1.0 | | v | | v'' | | -1.0, 0.0 | | v' |
É quase óbvio que este sistema é marginalmente estável, porque ele oscila para sempre, sem amortecer nem aumentar. (Na prática, um oscilador digital assim implementado seria instável devido a pequenos erros de arredondamento que se acumulam. Este artigo demonstra uma técnica melhor.)
É interessante verificar os autovalores da matriz A, que já sabemos ser marginalmente estável. Vamos aplicar uma fórmula para calcular autovalores, que pode ser usada diretamente porque nossa matriz é pequena:
determinante(A - λI) = 0 λ = algum autovalor da matriz A I = matriz identidade
Vamos calcular A - λI e achar a equação do determinante:
| 0 1 | - | λ 0 | = | -λ 1 | | -1 0 | | 0 λ | | -1 -λ | determinante: λλ + 1 autovalores satisfazem a equação: λλ + 1 = 0
As soluções são números imaginários: λ=-j,+j. A "assinatura" de um sistema oscilatório é a presença de pares conjugados complexos entre os autovalores, de forma análoga aos pólos complexos da transformada de Laplace. Autovalores ou pólos puramente imaginários denotam um sistema marginalmente estável.