UNIVERSIDADE FEDERAL DE UBERLÂNDIA
4ª Semana do Servidor e 5ª Semana Acadêmica
2008 – UFU 30 anos
IMPLEMENTAÇÃO E AVALIAÇÃO DE TÉCNICAS DE IDENTIFICAÇÃO
DE SISTEMAS LINEARES USANDO SOFTWARE LIVRE
Rosiane Ribeiro Rocha1
Universidade Federal de Uberlândia
Faculdade de Engenharia Química
Av. João Naves de Ávila, 2121, 38408-100
[email protected]
Luís Cláudio Oliveira Lopes2
[email protected]
Valéria Viana Murata2
[email protected]
Resumo: Neste trabalho são aplicados algoritmos de identificação de um sistema que representa
uma caldeira de aquecimento, conforme modelo proposto por Vasquez et al (2008). O modelo
descreve a variação da pressão de vapor que sai da caldeira, decorrente da variação no fluxo de
combustível alimentado. A caracterização e avaliação da dinâmica do processo foi feita aplicandose vários tipos de perturbações na entrada, concluindo-se sobre a linearidade do modelo. Os
modelos do tipo ARMAX foram obtidos através de um código implementado usando o software
livre Scilab. A seleção do modelo foi feita através de índices e critérios estruturais apropriados,
como MSE, AIC, FPE e BIC. A sintonia do modelo foi executada através da evolução diferencial.
Os parâmetros do modelo estimados representaram adequadamente o processo, demonstrando a
adequação do procedimento de estimação.
Palavras-chave: Identificação de sistema, ARMAX, sintonia de modelo.
1. INTRODUÇÃO
A modelagem matemática é a área do conhecimento que estuda maneiras de desenvolver e
implementar modelos matemáticos de sistemas reais. Há várias formas de classificar as técnicas de
modelagem. Uma delas agrupa os métodos em três categorias, denominadas modelagem caixa
branca, modelagem caixa preta e modelagem caixa cinza.
Na modelagem caixa branca é necessário conhecer bem o sistema em estudo, bem como as
leis físicas que descrevem o sistema a ser modelado. Por essa razão esse tipo de modelagem é
também conhecido como modelagem conceitual.
A identificação de sistemas é uma área do conhecimento que estuda técnicas alternativas à
modelagem caixa branca. Uma das características dessas técnicas é que pouco ou nenhum
conhecimento prévio do sistema é necessário e, consequentemente, tais métodos são também
referidos como modelagem ou identificação caixa preta ou modelagem empírica.
Em muitos casos será preferível usar técnicas de identificação para se obter modelos que
descrevem o comportamento de um sistema. O que se pretende descrever com tais modelos são as
relações de causa e efeito entre as variáveis de entrada e de saída. Nesse caso, o tipo de modelos, as
técnicas usadas e os requisitos necessários são bastante distintos dos correspondentes na
modelagem pela natureza do processo.
1 - Acadêmico do curso de Engenharia Química; 2 – Orientadores.
Uma categoria de técnicas intermediárias entre a modelagem pela natureza do processo e a
identificação caixa preta é chamada identificação caixa cinza. As técnicas desse grupo
caracterizam-se por usar informações auxiliares que não estão incluídas no conjunto de dados
utilizados durante a identificação. O tipo de informações auxiliares e a forma são usadas variam
muito entre as diversas técnicas disponíveis. Sendo assim, existem métodos de identificação caixa
cinza mais “claros”, que usam mais informação auxiliar, e mais “escuros”.
Esse trabalho apresenta um estudo de caso com aplicação de procedimentos para uma
identificação de sistema. Através de informações obtidas de um artigo desenvolvido por Vasquez et
al (2008), realizou-se diversas etapas para análise e aplicação de funções do software livre Scilab,
que simulam e identificam um sistema.
2. PROCEDIMENTOS PARA A IDENTIFICAÇÃO DE UM SISTEMA
A identificação de processos pode ser dividida nas seguintes etapas:
a) Testes dinâmicos e coleta de dados: os arquivos de dados do processo têm o mesmo papel
das equações constitutivas na modelagem teórica, pois fornecem as bases específicas para o
desenvolvimento de modelos para processos específicos. Como o modelo obtido por identificação é
totalmente baseado nos dados experimentais, é importante ter em mente que a informação que não
está contida nos dados não pode aparecer num passe de mágica no modelo, da mesma forma que
não é razoável esperar que uma equação constitutiva não especificada contribua para a qualidade do
modelo teórico final. As funções de entrada típicas empregadas na identificação de processos são:
degrau, impulso, pulso retangular ou arbitrário, ondas senoidais, ruído branco e seqüências binárias
pseudo-randômicas;
b) Escolha correta da estrutura dos modelos: determinação dos termos que devem compor os
modelos, preferencialmente de forma automática, através do reconhecimento da importância destes
diversos termos, utilizando os chamados dados de identificação e evitando a sobreparametrização
que ocorre quando são utilizados mais termos do que o necessário;
c) Estimação de parâmetros utilizando métodos numéricos adequados;
d) Verificação da capacidade dos modelos em representar o processo estudado.
3. O PROCESSO
3.1 Descrição do Processo
Com a finalidade de descrever o comportamento da variação da pressão de vapor em uma
caldeira com aquecimento através de gases, Vasquez et al (2008) coletaram dados de uma caldeira
em escala laboratorial cujo esquema está apresentado na Figura 1.
Figura 1: Esquematização do Processo (Vasquez)
2
Estudos preliminares sobre o comportamento dinâmico deste processo revelam que existem
diferentes interações entre as variáveis. Inicialmente define-se a demanda de vapor que ajusta o
acendimento do queimador até obter a produção de vapor desejada. Isto requer por sua vez que o
ventilador de ar forneça uma quantidade adequada de ar para a combustão do combustível. Ao
mesmo tempo, água deve ser alimentada na caldeira para se combinar com o vapor produzido.
Todos esses fenômenos são acoplados, e interagem intensamente afetando assim a dinâmica da
troca de calor dentro da caldeira. Como a pressão de vapor dentro da caldeira é a variável
dependente e o fluxo de combustível que será inflamado no queimador uma variável independente,
já que pode ser controlada, um modelo matemático poderá ser obtido usando-se a pressão de vapor
como a saída y(t) e o fluxo de combustível como a entrada u(t). Os autores apresentam uma análise
do ganho estático da variação da pressão de vapor dentro da caldeira que permite determinar a
região operacional onde o processo se comporta linearmente. Os resultados desta análise são
reproduzidos na Figura 2, onde observa-se que a região de comportamento aproximadamente linear
ocorre para abertura da válvula de combustível na faixa de 10 a 25% (Vasquez).
Figura 2: Relação entre a pressão de vapor e a abertura da válvula numa caldeira (reproduzido de
Vasquez et al (2008))
3.2 Dinâmica do Processo
O modelo que representa o comportamento do processo é dado pela Equação 1, conforme
apresentado por Vasquez et al (2008):
T1.T2 .
d 2 ΔP ( t )
dt 2
+ (T1 + T2 ) .
d ΔP ( t )
dt
+ ΔP ( t ) = K .Δ ( t − τ )
(1)
sendo T1 e T2 constantes do tempo, P(t) a saída do processo, u(t) a entrada, τ o atraso do tempo e K
o ganho do processo. Os valores dos parâmetros são: T1 ≈ 33,92 s, T2 ≈ 1,81 s, τ ≈ 35 s e K ≈ 11,3
kPa.
Um experimento numérico com o sinal de entrada foi feito para caracterizar a dinâmica do
processo e obter uma primeira estimativa da ordem e dos parâmetros do modelo. O sinal de entrada
u(t) é aplicado ao fluxo de combustível que entra no queimador e o comportamento da pressão de
vapor P(t) dentro da caldeira é medido.
Através do software livre Scilab (disponível no site: www.scilab.org), uma simulação do
processo foi executada. As condições iniciais para a presente simulação foram: para um tempo
inicial de 150 segundos, o ∆P(t) eram [0;0]
3
O gráfico que descreve o comportamento do processo ao longo do intervalo de tempo de 150
a 450 segundos é apresentado na Figura 3. Tal gráfico foi gerado pela resolução da equação
diferencial ordinária de 2º ordem, representada pela Equação 1, através do método de Runge Kutta,
para uma abertura da válvula de 10%.
Figura 3: Variação dinâmica da pressão de vapor
O estado estacionário do processo (Pvapor = 565 kPa) foi estimado através da função fsolve
do Scilab para uma abertura de 50% da válvula.
3.2.1 Perturbação Degrau
A caracterização da dinâmica do processo foi feita através de perturbações degrau de ± 10,
± 20 e ± 30 a partir do estado estacionário. Os resultados são mostrados na Figura 4.
Figura 4: Resultados das perturbações degraus a partir do estado estacionário
O modelo pode ser considerado linear, pois perturbações simétricas geraram respostas
simétricas da saída, conforme apresentado no gráfico.
3.2.2 Perturbação Randômica
Agora com uma perturbação mais intensa, mas ainda com o objetivo de avaliar a dinâmica
do processo, a perturbação randômica foi aplicada e dados da saída coletados para posterior análise
e avaliação. Esta perturbação foi desenvolvida através da função grand do software Scilab, que gera
um sinal totalmente randômico e diferente em cada execução. A Figura 5 apresenta a perturbação e
os resultados obtidos.
4
(a)
(b)
Figura 5: (a) Perturbação randômica e (b) Resposta da pressão de vapor à perturbação
randômica
A perturbação randômica possibilita excitar a planta em vários pontos, permitindo
posteriormente obter a partir desses dados, um modelo muito mais condizente com a realidade.
3.2.3 Perturbação PRBS
Aplicando um outro tipo de perturbação, do tipo PRBS, apresentada na Figura 6(a), obtémse os resultados representados na Figura 6(b). O sinal PRBS, além de ser de fácil geração, constrói
uma perturbação de grande controle. As principais características desse sinal são:
a) Pode tomar dois valores (+V e –V);
b) Só comuta entre níveis em instantes discretos t = 0, 1h, 2h...;
c) Os instantes de comutação são pre-determinados, i.e. o sinal é determinístico;
d) É periódico com período T0 = N.h, onde N é um número inteiro ímpar;
e) Em cada período existem (N+1)/2 intervalos a um nível e (N-1)/2 intervalos no outro
nível;
f) A função de autocorrelação, num período, é semelhante a um impulso (Lemos).
(b)
(a)
Figura 6: (a) Perturbação PRBS e (b) Resultados obtidos com a perturbação PRBS
Pela análise gráfica, nota-se que a planta responde de forma rápida a uma perturbação de
sinal PRBS, ou seja, o atraso do processo é pequeno.
4. IDENTIFICAÇÃO DO SISTEMA
4.1 Representação em Tempo Discreto: Modelo ARMAX
Considere o seguinte modelo geral:
A(q)y(k) = B(q)u(k) + C(q)v(k)
F(q)
D(q)
(2)
sendo q-1 o operador de atraso, de forma que y(k)q-1 = y(k-1), v(k) ruído branco e A(q), B(q), C(q),
D(q) e F(q) os polinômios definidos a seguir:
5
A(q) = 1 – a1q-1 - ... - anyq-ny
(3)
B(q) = b1q-1 + ... + bnuq-nu
(4)
C(q) =1 + c1q-1 + ... + cnvq-nv
(5)
D(q) = 1 + d1q-1 + ... + dndq-nd
(6)
F(q) = 1 + f1q-1 + ... + fnfq-nf
(7)
O modelo auto-regressivo com média móvel e entradas exógenas (ARMAX, do inglês
Autoregressive Moving Average with Exogenous Inputs) pode ser obtido a partir do modelo geral,
tomando-se D(q) = F(q) =1 e A(q), B(q) e C(q) polinômios arbitrários, resultando em
(8)
A(q)y(k) = B(q)u(k) + C(q)v(k)
ou, alternativamente,
(9)
y(k) = B(q)u(k) + C(q)v(k)
A(q)
A(q)
(10)
y(k) = H(q)u(k) + e(k)
sendo e(k) não branco. O modelo ARMAX pertence à classe de modelos de erro na equação e é
representado na Figura 7. No presente caso o erro na equação é modelado como um processo de
média móvel (MA), e o ruído adicionado à saída, e(k), é modelado como ruído branco filtrado pelo
filtro ARMA, C(q)/A(q) (Aguirre, 2007).
Figura 7: Representação esquemática do modelo ARMAX (reproduzido de Aguirre (2007))
4.2 Identificação do caso da caldeira
O modelo ARMAX foi utilizado para representar o sistema cujos dados foram gerados pela
perturbação randômica. Os possíveis modelos ARMAX que descrevem o processo foram obtidos
pelo software Scilab e avaliados por critérios adequados. Dos 1150 dados foram usados 575 para a
identificação e 575 para a validação.
A avaliação de um modelo pode ser executado através de índices e critérios estruturais,
como MSE (Erro Quadrático Médio), AIC (critério de informação de AKaike), FPE (do inglês,
Final Prediction Error) e BIC (do inglês, B Information Criterion). O índice MSE é definido pela
Equação 11:
n
MSE =
∑( y − y )
i =1
i
n
i
2
(11)
6
sendo n o número de amostras, yi o valor de referência e ŷi o valor previsto pelo modelo para a iésima amostra (mencionado por Morgano et al, 2008).
O critério AIC é obtido através da função de Máxima Verossimilhança, a partir dos
parâmetros ajustados para os modelos conforme os métodos. Este critério foi desenvolvido a partir
da distância ou informação de Kulback e Leibler (1951). Esta distância é uma medida de
discrepância entre as linhas do modelo verdadeiro e o modelo aproximado. Akaike (1983)
relacionou a distância de Kulback e Leibler com a Máxima Verossimilhança, surgindo o AIC. O
AIC é definido pela Equação 12:
AIC = - 2log(L) + 2K
(12)
sendo L a Verossimilhança Maximizada do modelo candidato e K, o número de parâmetros deste
modelo (mencionado por Shirotai et al, 2002).
Já o critério do erro de predição final (FPE), proposto por Akaike (1969) é dado pela
Equação 13:
⎛ N + ( p + 1)⎞
⎟⎟
FPE(p) = pˆ p .⎜⎜
⎝ N − ( p + 1)⎠
(13)
na qual N é o número de amostras, p é a ordem do modelo e p̂ p é a estimativa da variância do ruído
branco de entrada para o modelo de ordem p (mencionado por Shirotai et al, 2002).
E finalmente, o índice BIC é o critério B de informação, que usa conceitos estatísticos de
Bayes (Akaike, 1977; Akaike, 1978; Schwartz, 1978). A função a ser minimizada pela ordem ótima
é dada pela Equação 14:
⎧⎪
⎛ pˆ x
⎞⎫⎪
p⎞
⎛
BIC(p) = N . ln( pˆ p ) − (N − p ). ln⎜⎜1 − ⎟⎟ + p. ln( N ) + p. ln ⎨ p −1 .⎜⎜ − 1⎟⎟⎬
⎪⎩
⎝ N⎠
⎝ pˆ p ⎠⎪⎭
(14)
na qual p̂ x é a variância da saída do modelo AR (mencionado por Shirotai et al, 2002).
Através da avaliação de índices e critérios estruturais já mencionados, realizada com ordens
de 1 a 9 do modelo ARMAX, variando o auto-regressor r de 1 a 3 e o auto-regressor s de 1 a 3,
conclui-se que o de ordem 4 (r = 2 e s = 1) é o que melhor representa o processo em questão. Para
chegar a essa conclusão, além de verificar qual ordem que apresenta os menores valores dos índices
e critérios citados, foi levado em consideração a ordem que apresenta o menor esforço
computacional. A Figura 8 apresenta as avaliações executadas.
7
Figura 8: Avaliação da qualidade do modelo com índice MSE e critérios estruturais AIC, FPE e
BIC.
Os melhores valores de MSE, AIC, FPE e BIC estão na Tabela 1.
Tabela 1 – Melhores valores de MSE, AIC, FPE e BIC
MSE
0,0265110
AIC
-2252,051
FPE
-2252,0508
BIC
-2234,6335
A identificação e a validação do modelo ARMAX é mostrada na Figura 9.
Figura 9: Identificação e validação do modelo ARMAX
O melhor modelo ARMAX que ajusta os dados, fornecido pelo Scilab é descrito pelas
Equações 15 a 17:
8
A(q) = 1 – 1.5465337q1 + 0,5588348q2,
(15)
B(q) = - 0,139168 – 0,1391585q1 e
(16)
C(q) = 1.
(17)
Essas matrizes substituem seus correspondentes no modelo geral ARMAX descrito pela
Equação 8.
5. SINTONIA DO MODELO UTILIZANDO EVOLUÇÃO DIFERENCIAL
A Evolução Diferencial foi desenvolvida por Storn e Price em meados da década de noventa
e surgiu de tentativas de resolver o problema de ajuste polinomial de Chebychev (mencionado por
Arantes et al, 2006). Para caso em questão, foi utilizado o algoritmo fornecido pelo Matlab®, mas
adaptado para uso no software livre Scilab.
O algoritmo é iniciado criando uma população inicial escolhida aleatoriamente devendo
cobrir todo o espaço de busca. Geralmente, é criada por uma distribuição de probabilidade
uniforme, quando não há nenhum conhecimento sobre o problema.
A idéia principal da evolução diferencial é gerar novos indivíduos, denotados vetores
modificados ou doadores, pela adição da diferença ponderada entre dois indivíduos aleatórios da
população a um terceiro indivíduo. Esta operação é chamada mutação.
As componentes do indivíduo doador são misturadas com as componentes de um indivíduo
escolhido aleatoriamente (denotado vetor alvo), para resultar o chamado vetor tentativa, ou vetor
experimental. O processo de misturar os parâmetros é referido freqüentemente como "cruzamento"
na comunidade dos algoritmos evolutivos.
Se o vetor experimental resultar um valor da função objetivo menor que o vetor alvo, então
o vetor experimental substitui o vetor alvo na geração seguinte. Esta última operação é chamada
seleção. O procedimento é finalizado através de algum critério de parada.
Os parâmetros T1 (constante do tempo) e K (ganho do processo) do modelo descrito na
Equação 1 foram obtidos através de 10 execuções do algoritmo de Evolução Diferencial. Os
melhores valores encontrados para os parâmetros foram: T1 = 33,780974 s e K = 11,345646 s.
A Tabela 2 apresenta uma comparação entre os parâmetros teóricos e os identificados com
respectivos erros relativos.
Tabela 2: Erros relativos entre os parâmetros identificados e os teóricos.
Parâmetros Teóricos
Parâmetros Identificados
Erro relativo
T1 = 33,92 s
T1 = 33,780974 s
0,41%
K = 11,3 s
T2 = 11,345646 s
0,40%
Pela análise da tabela acima, nota-se que a sintonia foi excelente e os parâmetros foram
identificados de forma satisfatória.
6. AGRADECIMENTOS
Agradeço à FAPEMIG pela oportunidade e incentivo financeiro. Aos professores Valéria
Viana e Luís Cláudio Lopes, pela dedicação e apoio.
9
7. REFERÊNCIAS
Aguirre, L. A. (2007). Introdução à Identificação de Sistemas: Técnicas Lineares e Não Lineares
Aplicadas a Sistemas Reais, 3 edn, Editora UFMG, Belo Horizonte,BH.
Arantes, M. B., Oliveira, G.T.S. e Saramago, S.F.P., Evolução Diferencial aplicada à solução de
alguns problemas de Engenharia de Produção, FAMAT em revista, nº 6, Maio de 2006, Uberlândia,
MG.
Ken Price, Rainer Storn, and Jouni Lampinen (2005): Differential Evolution - A Practical Approach
to Global Optimization . Springer.
Lemos, J.M. e Bernardino, A., Modelação, Identificação e Controle Digital – Métodos não
paramétrico,
ISTSecção
de
Sistemas
e
Controle,
disponível
no
site
<users.isr.ist.utl.pt/~alex/micd0506/micd3b.pdf>. Acessado em 20 de agosto de 2008.
Morgano, M.A., Faria, C.G., Ferrão, M.F., Bragagnolo, N. e Ferreira, M.M.C., 2008, Determinação
de umidade em café cru usando espectroscopia NIR e regressão multivariada, Ciência e Tecnologia
de Alimentos, vol.28 nº1, Campinas, SP.
Shirotai, C. e Itiki, C., 2002, Determinação da ordem na modelagem auto-regressiva de sinais
eletromiográficos, Produção em Iniciação Científica da EPUSP, São Paulo, SP.
Vasquez, J. R. R., Perez, R. R., Moriano, J. S., & Gonzalez, J. R. P., System identification of steam
pressure in a fire-tube boiler, Computers and Chemical Engineering (2007),
doi:10.1016/j.compchemeng.2008.01.010.
IMPLEMENTATION AND EVALUATION OF TECHNIQUES OF
IDENTIFICATION OF LINEAR SYSTEMS USING FREE SOFTWARE
Rosiane Ribeiro Rocha
Federal University of Uberlândia
School of Chemical Engineering
2121, Avenue João Naves de Ávila, 38408-100
[email protected]
Luís Cláudio Oliveira Lopes
[email protected]
Valéria Viana Murata
[email protected]
Abstract: In this article, identification algorithms are applied to a system that represents a fire-tube
boiler, as modeled by Vasquez et al. (2008). The model describes the change in the pressure of
steam that leaves the boiler, resulting from the change in the fuel flow fuel supplied. The
characterization and evaluation of the dynamics of the process were made according to various
kinds of disturbances at the entrance, concluding about the linearity of the model. The models of the
ARMAX type were obtained through a code implemented using the free software Scilab.. The
selection of the model was made through indices and structural appropriate criteria, such as MSE,
AIC, FPE and BIC. Model tuning was addressed by using the Differential Evolution. The
parameters of the estimated model properly accounted for the process behavior, demonstrating the
suitability of the estimation procedure.
Keywords: Systems Identification, ARMAX, model tuning.
10