MODELAGEM DE COLUNAS DE DESTILAÇÃO ATRAVÉS DE
MODELOS AUTO-REGRESSIVOS
ADELSON SIQUEIRA CARVALHO
UNIVERSIDADE ESTADUAL DO NORTE FLUMINENSE - UENF
CAMPOS DOS GOYTACAZES – RJ
MARÇO – 2008
ii
MODELAGEM DE COLUNAS DE DESTILAÇÃO ATRAVÉS DE
MODELOS AUTO-REGRESSIVOS
ADELSON SIQUEIRA CARVALHO
Tese apresentada ao Centro de Ciências e
Tecnologia, da Universidade Estadual do Norte
Fluminense, como parte das exigências para
obtenção do titulo de Mestre em Engenharia de
Produção.
Orientador: Luís Humberto Guillermo Felipe, D.Sc. UENF/LEPROD
CAMPOS DOS GOYTACAZES – RJ
MARÇO – 2008
iii
MODELAGEM DE COLUNAS DE DESTILAÇÃO ATRAVÉS DE
MODELOS AUTO-REGRESSIVOS
Autor:
ADELSON SIQUEIRA CARVALHO
Tese apresentada ao Centro de Ciências e Tecnologia, da Universidade Estadual do
Norte Fluminense, como parte das exigências para obtenção do titulo de Mestre em
Engenharia de Produção.
Aprovado em ____ de ____________ de 2008
Comissão Examinadora:
_________________________________________
Prof. Geraldo Galdino de Paula Júnior (Doutor, Engenharia de Sistemas e
Computação) – UENF/LEPROD
_________________________________________
Prof. José Ramón Arica Chávez (Doutor, Engenharia de Sistemas e Computação) –
UENF/LEPROD
_________________________________________
Prof. Carlos Antonio Abanto Valle (Doutor, Ciências Estatísticas) – UFRJ/IM
_________________________________________
Prof. Luis Humberto Guillermo Felipe (Doutor, Ciências da Engenharia) –
UENF/LEPROD (orientador)
iv
“Dedico este trabalho a meus
pais e entes queridos”
v
AGRADECIMENTOS
Agradeço em primeiro lugar a Deus por encaminhar a minha vida e por me
permitir alcançar meus objetivos.
A minha família e amigos por compreenderem minhas falhas e apoiarem as
minhas escolhas.
A Aline por ser uma grande companheira e inspiração para minha trajetória.
Ao meu orientador, professor Luis Humberto Guillermo Felipe, por toda
atenção e apoio.
vi
RESUMO
Sistemas dinâmicos são em sua grande maioria não lineares e de dinâmica não
totalmente conhecida, desta forma existe certa distância entre o sistema real e o
gerado a partir da modelagem analítica. Modelos mais próximos da realidade podem
ser obtidos a partir de métodos que utilizem um conjunto de dados oriundo do
sistema real, dentre estes os modelos auto-regressivos merecem destaque. O
presente trabalho apresenta a modelagem de uma coluna de destilação piloto com
esta metodologia, destacando suas particularidades. Os modelos em questão são
capazes de fornecer representação linear discreta do sistema dinâmico. Os modelos
obtidos são validados com índices de desempenho e análise dos resíduos da
modelagem. Os resultados indicam o potencial da metodologia aplicada para
sistemas de porte industrial devido à fidelidade do sistema de aquisição de dados do
objeto de estudo. Como principal contribuição, a validação do modelo com dados de
outros quatro testes dinâmicos que não o utilizado para estimação dos parâmetros.
vii
ABSTRACT
The most of dynamic systems are non-linear and it has dynamic that is not fully
known, this way there is a distance between the real system and the generated from
analytical modeling. Models closer to reality can be obtained from methods that use a
set of data derived from real system, among these, auto-regressive can be
emphasized. This work presents a modeling of a pilot distillation column using this
methodology, highlighting its details. These models are able to provide linear discrete
representation of dynamic system. Models gotten are validated with performance
indexes and modeling residue analysis. Results indicate potential of the applied
methodology to industrial systems due to accuracy of data acquisition system of
object of research. Like main contribution, to validate the model with other four
dynamic tests different from the dynamic test used to parameters estimation.
viii
LISTA DE SIGLAS E ABREVIATURAS
AIC
Akaike Information Criteria
CEFET- Campos
Centro Federal de Educação Tecnológica de Campos
o
Grau Gay Lussac
GL
DFI
Device Fieldbus Interface
MIMO
Multiple-Input and Multiple-Output
ARX
Auto-Regressive with eXogenous input
ARIMA
Auto-Regressive Integrated Moving Average
AR
Auto-Regressive
OE
Output Error
BJ
Box-Jenkins
ARMAX
Auto-Regressive Moving Average with eXogenous input
CC
Continuous Current
MATLAB
MATrix LABoratory
LABVIEW
Laboratory Virtual Instrumentation Engineering Workbench
PRBS
Pseudo Random Binary Signal
FF
Foundation Fieldbus
NARX
Non linear Auto-Regressive with eXogenous input
NARMAX
Non linear Auto-Regressive Moving Average with eXogenous
input
FPE
Final Prediction Error
CCF
Craquamento Catalítico em leito Fluidizado
MRSE
Mean Regularizated Square Error
MSE
Mean Square Error
MAE
Mean Absolute Error
MAPE
Mean Absolute Percentage Error
ix
LISTA DE FIGURAS
Figura 1 - Esquemas com colunas A, B e C para os processos de destilação,
retificação e desidratação (Rasovsky, 1973)........................................................7
Figura 2 - Fracionador de pratos – coluna ou torre fracionadora (Rasovsky, 1973)...8
Figura 3 - Coluna de destilação piloto (Crespo, 2000). ..............................................9
Figura 4 - Coluna de destilação piloto instalada (Naegele, 2000). ...........................10
Figura 5 - Instrumentos interligados em rede fieldbus.............................................11
Figura 6 - Coluna de destilação piloto com a instrumentação atual. .......................12
Figura 7 - Trocadores de calor da coluna de destilação...........................................13
Figura 8 - Base da coluna de destilação. .................................................................13
Figura 9 - Ponto de alimentação da coluna de destilação........................................14
Figura 10 - Sistema de armazenamento da mistura binária. ....................................15
Figura 11 - Método dos mínimos quadrados. ...........................................................26
Figura 12 - Ciclo de identificação com auxílio computacional. Os retângulos são de
responsabilidade do computador e as elipses responsabilidade do usuário
(Adaptado de Ljung (1994)). ..............................................................................30
Figura 13 - Covariância cruzada das temperaturas nos estágios e a composição
(Zanata, 2005)....................................................................................................32
Figura 14 - Funções de correlação cruzada (a, b) entre a variável de saída e duas
variáveis candidatas a serem entradas, (c) entre (outras) duas variáveis
candidatas a serem entradas (Aguirre, 2000). ...................................................33
Figura 15 - Dados utilizados para a seleção das variáveis do modelo. ....................35
Figura 16 - Correlação cruzada entre a vazão e a temperatura no topo. .................36
Figura 17 - Correlação cruzada entre a temperatura na base e a temperatura no
topo. ...................................................................................................................36
Figura 18 - Correlação cruzada entre a pressão no topo e a temperatura no topo. .37
Figura 19 - Correlação cruzada entre o nível da base e a temperatura no topo. .....37
Figura 20 - Correlação cruzada entre a vazão e a temperatura da base. ................38
Figura 21 - Correlação cruzada entre a pressão e a vazão......................................39
Figura 22 - Correlação cruzada entre a temperatura da base e a pressão. .............39
Figura 23 - Sinal aplicado na entrada do sistema para realizar a identificação linear,
variável vazão de gás (Dallagnol, 2005). ...........................................................41
x
Figura 24 - Estímulo na forma de degrau positivo na vazão de produto da coluna de
destilação, para o primeiro teste dinâmico. ........................................................42
Figura 25 - Resposta da variável de saída (composição) a uma perturbação na
variável de entrada do sistema. Resposta com menor constante de tempo
(Zanata, 2005)....................................................................................................44
Figura 26 - Dados adequados para a nova taxa de amostragem.............................50
Figura 27 - Conjunto de estimação. .........................................................................56
Figura 28 - Primeiro conjunto de validação. .............................................................56
Figura 29 - Dados do segundo teste dinâmico utilizado para validação...................57
Figura 30 - Dados do terceiro teste dinâmico utilizado para validação.....................58
Figura 31 - Dados do quarto teste dinâmico utilizado para validação. .....................58
Figura 32 - Dados do quinto teste dinâmico utilizado para validação.......................59
Figura 33 - Valores do AIC calculado através das iterações da rotina. (Carvalho e
Guillermo, 2007).................................................................................................60
Figura 34 - Curva do AIC calculado através das iterações da rotina atual. ..............61
Figura 35 - Saída do modelo comparada com a saída do sistema, para o primeiro
conjunto de validação.........................................................................................64
Figura 36 - Saída do modelo comparada com a saída do sistema, para o segundo
conjunto de validação.........................................................................................65
Figura 37 - Saída do modelo comparada com a saída do sistema, para o terceiro
conjunto de validação.........................................................................................65
Figura 38 - Saída do modelo comparada com a saída do sistema, para o quarto
conjunto de validação.........................................................................................66
Figura 39 - Saída do modelo comparada com a saída do sistema, para o quinto
conjunto de validação.........................................................................................66
Figura 40 - Saída do modelo comparada com a saída do sistema, para o primeiro
conjunto de validação.........................................................................................68
Figura 41 - Saída do modelo comparada com a saída do sistema, para o segundo
conjunto de validação.........................................................................................69
Figura 42 - Saída do modelo comparada com a saída do sistema, para o terceiro
conjunto de validação.........................................................................................69
Figura 43 - Saída do modelo comparada com a saída do sistema, para o quarto
conjunto de validação.........................................................................................70
xi
Figura 44 - Saída do modelo comparada com a saída do sistema, para o quinto
conjunto de validação.........................................................................................70
Figura 45 - Saída do modelo comparada com a saída do sistema, para o primeiro
conjunto de validação.........................................................................................72
Figura 46 - Saída do modelo comparada com a saída do sistema, para o segundo
conjunto de validação.........................................................................................72
Figura 47 - Saída do modelo comparada com a saída do sistema, para o terceiro
conjunto de validação.........................................................................................73
Figura 48 - Saída do modelo comparada com a saída do sistema, para o quarto
conjunto de validação.........................................................................................73
Figura 49 - Saída do modelo comparada com a saída do sistema, para o quinto
conjunto de validação.........................................................................................74
Figura 50 - MSE calculado para o erro de predição de k passos à frente para o
primeiro conjunto de validação...........................................................................79
Figura 51 - Comparação entre saída do modelo e saída do sistema para a predição
de 20 passos à frente para o primeiro conjunto de validação. ...........................79
Figura 52 - MSE calculado para o erro de predição de k passos à frente para o
segundo conjunto de validação..........................................................................80
Figura 53 - Comparação entre saída do modelo e saída do sistema para a predição
de 20 passos à frente para o segundo conjunto de validação. ..........................80
Figura 54 - MSE calculado para o erro de predição de k passos à frente para o
terceiro conjunto de validação............................................................................81
Figura 55 - Comparação entre saída do modelo e saída do sistema para a predição
de 20 passos à frente para o terceiro conjunto de validação. ............................81
Figura 56 - MSE calculado para o erro de predição de k passos à frente para o
quarto conjunto de validação. ............................................................................82
Figura 57 - Comparação entre saída do modelo e saída do sistema para a predição
de 20 passos à frente para o quarto conjunto de validação. ..............................82
Figura 58 - MSE calculado para o erro de predição de k passos à frente para o
quinto conjunto de validação..............................................................................83
Figura 59 - Comparação entre saída do modelo e saída do sistema para a predição
de 20 passos à frente para o quinto conjunto de validação. ..............................83
Figura 60 - Autocorrelação dos resíduos para o conjunto de estimação..................84
xii
Figura 61 - Correlação cruzada entre os resíduos e a entrada para o conjunto de
estimação...........................................................................................................85
Figura 62 - Autocorrelação dos resíduos para o primeiro conjunto de validação. ....85
Figura 63 - Correlação cruzada entre os resíduos e a entrada para o primeiro
conjunto de validação.........................................................................................86
Figura 64 - Autocorrelação dos resíduos para o segundo conjunto de validação. ...86
Figura 65 - Correlação cruzada entre os resíduos e a entrada para o segundo
conjunto de validação.........................................................................................87
Figura 66 - Autocorrelação dos resíduos para o terceiro conjunto de validação. .....87
Figura 67 - Correlação cruzada entre os resíduos e a entrada para o terceiro
conjunto de validação.........................................................................................88
Figura 68 - Autocorrelação dos resíduos para o quarto conjunto de validação........88
Figura 69 - Correlação cruzada entre os resíduos e a entrada para o quarto conjunto
de validação. ......................................................................................................89
Figura 70 - Autocorrelação dos resíduos para o quinto conjunto de validação. .......89
Figura 71 - Correlação cruzada entre os resíduos e a entrada para o quinto conjunto
de validação. ......................................................................................................90
Figura 72 - Autocorrelação dos resíduos para o conjunto de estimação..................91
Figura 73 - Correlação cruzada entre os resíduos e a entrada para o conjunto de
estimação...........................................................................................................91
Figura 74 - Autocorrelação dos resíduos para o primeiro conjunto de validação. ....92
Figura 75 - Correlação cruzada entre os resíduos e a entrada para o primeiro
conjunto de validação.........................................................................................92
Figura 76 - Autocorrelação dos resíduos para o segundo conjunto de validação. ...93
Figura 77 - Correlação cruzada entre os resíduos e a entrada para o segundo
conjunto de validação.........................................................................................93
Figura 78 - Autocorrelação dos resíduos para o terceiro conjunto de validação. .....94
Figura 79 - Correlação cruzada entre os resíduos e a entrada para o terceiro
conjunto de validação.........................................................................................94
Figura 80 - Autocorrelação dos resíduos para o quarto conjunto de validação........95
Figura 81 - Correlação cruzada entre os resíduos e a entrada para o quarto conjunto
de validação. ......................................................................................................95
Figura 82 - Autocorrelação dos resíduos para o quinto conjunto de validação. .......96
xiii
Figura 83 - Correlação cruzada entre os resíduos e a entrada para o quinto conjunto
de validação. ......................................................................................................96
xiv
LISTA DE TABELAS
Tabela 1 - Tabela de polinômios retirados do modelo geral......................................23
Tabela 2 - Tabela que relaciona os modelos ao tipo de situação a qual devem ser
utilizados. ...........................................................................................................46
Tabela 3 – Valores calculados do AIC para modelos ARX de diferentes ordens. .....62
Tabela 4 - Índices calculados para os cinco conjuntos de validação referente ao
modelo obtido com o AIC. ..................................................................................67
Tabela 5 - Índices calculados para os cinco conjuntos de validação referente ao
primeiro modelo obtido com a validação cruzada. .............................................71
Tabela 6 - Índices calculados para os cinco conjuntos de validação referente ao
segundo modelo obtido com a validação cruzada..............................................74
Tabela 7 - Coeficientes estimados para o modelo ARX [20 20 0]. ............................75
Tabela 8 - Coeficientes estimados para o modelo ARX [16 15 15]. ..........................76
Tabela 9- Coeficientes estimados para o modelo ARX [6 12 12]. .............................77
Tabela 10 - MSE calculado para o erro de predição de k passos à frente do modelo
ARX [6 12 12] para os cinco conjuntos de validação. ........................................78
xv
SUMÁRIO
AGRADECIMENTOS.................................................................................................vi
RESUMO...................................................................................................................vii
ABSTRACT..............................................................................................................viii
LISTA DE SIGLAS E ABREVIATURAS ....................................................................ix
LISTA DE FIGURAS...................................................................................................x
LISTA DE TABELAS ................................................................................................xv
CAPÍTULO I ................................................................................................................1
INTRODUÇÃO ............................................................................................................1
1.1- APRESENTAÇÃO ............................................................................................1
1.2- MOTIVAÇÃO ....................................................................................................2
1.3- OBJETIVO ........................................................................................................3
1.4- ORGANIZAÇÃO ...............................................................................................4
CAPÍTULO II ...............................................................................................................6
DESTILAÇÃO: UM BREVE HISTÓRICO ...................................................................6
2.1- COLUNAS DE DESTILAÇÃO...........................................................................6
2.2 – OBJETO DE ESTUDO....................................................................................9
2.3 – CARACTERÍSTICAS FÍSICAS .....................................................................11
2.4 – O PROCESSO DE OBTENÇÃO DE ÁLCOOL HIDRATADO NA COLUNA
PILOTO..................................................................................................................15
2.5 – ALGUNS PROBLEMAS DO SISTEMA.........................................................16
2.6 – MODELAGEM DE COLUNAS DE DESTILAÇÃO.........................................17
CAPÍTULO III ............................................................................................................21
FERRAMENTAL TEÓRICO......................................................................................21
3.1- MODELOS AUTO-REGRESSIVOS................................................................21
3.1.1- ESTADO DA ARTE DA APLICAÇÃO DE MODELOS AUTOREGRESSIVOS .................................................................................................24
3.2- ESTIMAÇÃO DE PARÂMETROS...................................................................25
3.2.1- MÉTODO DOS MÍNIMOS QUADRADOS................................................26
CAPÍTULO IV............................................................................................................29
METODOLOGIA .......................................................................................................29
4.1- INTRODUÇÃO................................................................................................29
xvi
4.2 – ESCOLHA DAS VARIÁVEIS DO MODELO..................................................30
4.3 – ESCOLHA DOS ESTÍMULOS DE ENTRADA ..............................................40
4.4 – ESCOLHA DO PERÍODO DE AMOSTRAGEM ............................................42
4.5 – SELEÇÃO DA ESTRUTURA DO MODELO .................................................45
4.6 – SELEÇÃO DA ORDEM DO MODELO ..........................................................46
4.7 – ESTIMAÇÃO E VALIDAÇÃO ........................................................................47
4.8 – ÍNDICES DE DESEMPENHO .......................................................................51
4.9 – ANÁLISE DE RESÍDUOS .............................................................................53
CAPÍTULO V.............................................................................................................55
RESULTADOS..........................................................................................................55
5.1- DADOS UTILIZADOS.....................................................................................55
5.2- MÉTODOS DE SELEÇÃO DA ORDEM DO MODELO...................................59
5.3- RESULTADOS GRÁFICOS E NUMÉRICOS DOS MODELOS OBTIDOS.....63
5.4- COEFICIENTES ESTIMADOS DOS MODELOS OBTIDOS...........................74
5.5- RESULTADOS DOS MODELOS VARIANDO O HORIZONTE DE PREDIÇÃO
...............................................................................................................................77
5.6- RESULTADOS DA ANÁLISE DOS RESÍDUOS .............................................84
CAPÍTULO VI............................................................................................................98
CONCLUSÃO ...........................................................................................................98
REFERÊNCIAS.......................................................................................................100
xvii
1
CAPÍTULO I
INTRODUÇÃO
1.1- APRESENTAÇÃO
O problema abordado neste trabalho decorre da dificuldade ao se modelar
sistemas dinâmicos1 e não lineares para otimizá-los, haja vista que estes não
possuem explicitadas as relações entre suas variáveis e as não linearidades destas
relações.
Os sistemas dinâmicos encontrados na prática são, em última análise, não
lineares. É bem verdade que em alguns casos aproximações lineares são suficientes
para aplicações práticas (Aguirre, 2000).
Sistemas dinâmicos estão presentes na maioria dos processos de
transformação, onde as variáveis que influenciam na qualidade do produto final são
contínuas. Dentre as quais se destacam: temperatura, vazão, nível e pressão. As
indústrias petroquímicas, álcool, celulose e de cimento se enquadram nesta classe.
Nestas indústrias, cuja produção é contínua, e que possuem como essência a
separação de componentes da matéria-prima, a destilação é um processo de
predominância. Sobretudo nas indústrias petroquímica e do álcool é o principal
processo de transformação de matéria prima em produto final.
Colunas de destilação ou torres de destilação estão presentes em grande
número nas usinas de produção de álcool e também em refinarias de petróleo, e
unidades de processamento de gás-natural. Estes processos industriais são
perfeitos exemplos de sistemas dinâmicos, multivariáveis e de comportamento não
linear.
Na indústria é comum que a matéria prima seja destilada inicialmente em uma
coluna A (fracionadora) obtendo assim uma gama de derivados. Alguns destes por
sua vez seguem para uma coluna B (retificadora) para a extração de derivados de
1
Sistemas dinâmicos são sistemas cujas variáveis possuem determinado comportamento que sofre
influência do fator tempo.
2
maior valor agregado. Este procedimento está presente tanto na indústria do
petróleo quanto na do álcool.
As colunas de destilação industriais são sistemas dinâmicos, multivariáveis, e
não lineares; a coluna de destilação existente no CEFET-Campos é um exemplo.
Esses sistemas para serem otimizados primeiro devem ser modelados. O problema
está, então, em encontrar uma classe de modelos para representá-los.
Carvalho, et al. (2004) utilizaram redes neurais artificiais do tipo MLP (Multi
Layer Perceptron) para modelar a coluna de destilação do CEFET-Campos. Os
dados utilizados foram registrados por uma placa de aquisição de dados do tipo CIODAS08/jr-AO, com resolução de 12 bits. As variáveis com dados disponíveis são a
vazão de alimentação e a temperatura de topo, únicas possíveis de mensuração. A
rede neural MLP contém três entradas: vazão no instante atual, vazão no instante
anterior, temperatura no instante anterior. A variável de saída selecionada foi a
temperatura de topo no instante atual. Os resultados do trabalho foram registrados
apenas de forma gráfica, comparando a saída do modelo com a saída observada do
sistema.
A proposta deste trabalho é utilizar modelos auto-regressivos para modelar a
relação entre as principais variáveis da coluna de destilação: vazão de alimentação,
temperatura da base, nível da base, temperatura de topo e pressão no topo, levando
em consideração as restrições construtivas do equipamento e do sistema de
medição.
Os modelos auto-regressivos possuem vasto histórico de aplicações na
modelagem de sistemas industriais e de séries temporais.
1.2- MOTIVAÇÃO
Pesquisa operacional é o nome dado a um conjunto de modelos e algoritmos
destinados a determinar o melhor curso das ações que visam garantir o
funcionamento ótimo de sistemas, sob restrições de recursos escassos (Paula
Junior, 1998).
Neste processo de busca pelo “ótimo”, a pesquisa operacional passa a ser um
método com duas etapas bem distintas, a modelagem do sistema como um
3
problema de otimização ou modelo de programação matemática, e a resolução
deste modelo através de algum algoritmo.
A modelagem dos sistemas ou processos que serão otimizados é um misto
de arte e ciência; sua fundamentação é matemática, mas a visão do sistema varia de
modelador para modelador, em suma, pode-se obter modelagens diferentes para o
mesmo problema, e ambas estarem corretas.
Modelagem caixa branca é também conhecida como modelagem pela física
ou modelagem conceitual. Contudo, devido ao conhecimento e tempo necessário
para modelar um sistema partindo do equacionamento dos fenômenos envolvidos,
nem sempre é viável seguir este procedimento (Aguirre, 2000).
Sistemas não lineares possuem descontinuidades, ou seja, diferentes
relações entre as mesmas variáveis tratadas, à medida que se caminha pelas
diferentes partes do domínio.
Os modelos auto-regressivos podem ser ajustados aos dados amostrados de
um sistema industrial, com o intuito de fornecer a representação matemática
necessária para o objetivo do pesquisador.
Tanto a otimização de colunas de destilação quanto o controle das variáveis
de processo são de grande interesse do ponto de vista da produtividade.
Considerando que as duas áreas trabalham em conjunto, os valores operacionais
ideais para as variáveis de processo podem ser obtidos por ferramentas de pesquisa
operacional e a garantia da manutenção destas variáveis nos valores ideais é
fornecido pela teoria de controle.
Nos dois casos a necessidade de modelagem do sistema em estudo é de
primeira necessidade. A técnica proposta constitui-se como ferramenta para
profissionais que necessitam modelar sistemas de natureza dinâmica, considerando
o fato de que para tanto, não necessitam de modelo fenomenológico, mas sim um
conjunto de dados representativos do comportamento deste.
1.3- OBJETIVO
O objetivo principal deste trabalho é a modelagem da coluna de destilação do
CEFET-Campos, suprindo a ausência de modelo matemático que reproduza seu
comportamento dinâmico.
Como objetivos intermediários podem-se destacar:
4

Domínio de ferramental teórico necessário, acerca dos modelos auto-
regressivos, estes selecionados para representar a coluna de destilação;

Documentação dos resultados dos modelos aplicados;
O objetivo do trabalho não é encontrar a solução ótima para o problema de
produção alcoólica da coluna de destilação, mas sim, obtenção de um modelo
matemático que possa ser utilizado adiante para este e outros fins.
1.4- ORGANIZAÇÃO
A dissertação está organizada da seguinte forma:
Capitulo I – Introdução
Esse capítulo fornece uma introdução a respeito dos assuntos que serão
abordados ao longo da monografia, descrevendo de forma sucinta os objetivos, a
motivação e como estão organizados os capítulos.
Capítulo II – Destilação
O segundo capítulo trata de questões relativas ao problema abordado, ou
seja, o processo de destilação e a modelagem de colunas de destilação. Um breve
histórico, técnicas utilizadas para a modelagem e descrição do objeto de estudo.
Capítulo III – Ferramental Teórico
No terceiro capítulo, são apresentados conceitos relativos aos modelos autoregressivos, classe de modelos escolhida para representar matematicamente o
sistema. Este capítulo contém explanação sobre os pontos mais importantes de tal
paradigma e também o estado-da-arte.
Neste capítulo, é apresentada também a fundamentação do método dos
estimadores mínimos quadrados. É também exposto o método quando aplicado ao
problema de estimação dos coeficientes do modelo auto-regressivo.
5
Capítulo IV – Metodologia
A partir do ferramental teórico apresentado nos capítulos anteriores,
desenvolve-se neste a metodologia utilizada para aplicação do ferramental teórico
aos dados do problema real. Os artifícios e detalhes técnicos são apresentados e os
critérios utilizados nas etapas da metodologia para obtenção do modelo.
Capítulo V – Resultados Obtidos
Neste capítulo são apresentados os resultados obtidos por meio da aplicação
da metodologia. O desempenho do modelo ajustado será verificado por índices de
desempenho e seus resultados serão sumariados em tabelas e a aderência da
predição do modelo aos dados observados do sistema será contemplada de forma
gráfica. Os resíduos do modelo serão analisados.
Capítulo VI – Conclusões e Recomendações
O último capítulo apresenta as conclusões em relação ao trabalho
desenvolvido. Sugestões de trabalhos futuros também são apresentadas.
6
CAPÍTULO II
DESTILAÇÃO: UM BREVE HISTÓRICO
2.1- COLUNAS DE DESTILAÇÃO
Os processos de separação de diversos líquidos tem origem no século II da
era cristã, quando, no Egito, Zozime e Hermes foram considerados os mestres na
arte da destilação. Dos séculos VII ao XII, os árabes dedicam-se ao processo de
destilação nos chamados alambiques. A partir do século XVI o avanço foi maior e no
século passado (XIX) conseguiram desenvolver o alambique e chegar ao ponto de
partida dos equipamentos de destilação (colunas) que produzem álcool de forma
contínua (Rasovsky, 1973).
Para prosseguir no entendimento do processo de destilação e suas
particularidades faz-se necessário definir a destilação.
O princípio da destilação consta do fenômeno de fracionamento dos líquidos,
onde os mais voláteis, com pontos de ebulição mais baixos separam-se em primeiro
lugar, seguidos pelos outros componentes em seqüência correspondente às suas
respectivas volatilidades (Rasovsky, 1973).
O processo de destilação tem por objetivo a separação de componentes de
uma fase líquida através de sua vaporização parcial. Estes vapores oriundos da
vaporização são ricos de componentes mais voláteis do que o líquido, o que permite
a separação (Gomide, 1988).
As colunas de destilação são os mais importantes equipamentos para a
separação de uma mistura de líquidos miscíveis em seus componentes na indústria
química e petroquímica. Esta separação é realizada aproveitando-se o fato de os
elementos constituintes da mistura terem diferentes temperaturas de ebulição. Assim
através do fornecimento de calor à mistura consegue-se preferencialmente vaporizar
as substâncias mais voláteis, que são condensadas no topo da coluna, enquanto as
menos voláteis tendem a permanecer na fase líquida do fundo da coluna (Campos e
Teixeira, 2006).
As operações de destilação em geral baseiam-se no mesmo princípio, porém,
difere entre si, principalmente pelo tipo de equipamento usado.
7
De todos os processos, é na destilação fracionada que se obtém um maior
enriquecimento do vapor produzido, isto acontece em um equipamento chamado de
coluna ou torre de fracionamento2. Esta coluna é composta por pratos ou bandejas
horizontais perfuradas.
Na Figura 1 podem ser visualizados os três tipos de colunas utilizadas na
destilação do álcool.
Figura
1 - Esquemas com colunas A, B e C para os processos de destilação, retificação e
desidratação (Rasovsky, 1973).
Ao ser aquecido no refervedor3 o vapor sobe pela coluna em contra-fluxo à
mistura que desce pela mesma, através da vazão de refluxo. O refluxo estabelece,
portanto uma corrente líquida que desce de prato em prato na coluna (Gomide,
1988).
Na Figura 2 é apresentada o fracionador de pratos, uma representação
gráfica para uma coluna de destilação contendo refluxo e refervedor.
2
3
Equipamento industrial capaz de separar por vaporização componentes de uma mistura.
Trocador de calor situado na base da torre com a finalidade aquecer a mistura.
8
Figura 2 - Fracionador de pratos – coluna ou torre fracionadora (Rasovsky, 1973).
O álcool produzido do vinho (oriundo do processo de fermentação do caldo da
cana de açúcar) considerado um sistema multicomponente, é obtido usando-se as
colunas de fracionamento. De acordo com sua graduação alcoólica dividem-se em
álcool bruto ou fraco ( 52 a 92
o
GL), álcool retificado ou hidratado (94 a 96
o
GL) e
álcool anidro ou absoluto (99,95 o GL).
O álcool fraco é produzido em uma primeira coluna de fracionamento (coluna
A); este é retificado em uma segunda coluna (coluna B) e finalmente para o
processo de obtenção do álcool anidro (coluna C), que usa um terceiro agente
chamado de arrastador4 ou outro elemento apropriado, eliminando-se do álcool
retificado a parcela d’água nele contida que não se consegue separar mais pelo
fracionamento (Rasovsky, 1973).
4
Tipo de equipamento utilizado para realizar destilação por arraste.
9
2.2 – OBJETO DE ESTUDO
A partir de um sistema real instalado no laboratório de pesquisa e automação
do CEFET-Campos, podem-se coletar dados das variáveis: temperatura, pressão
vazão e nível de uma coluna de destilação piloto, e a partir deles ajustar uma classe
de modelos paramétricos que reproduza de maneira satisfatória o comportamento do
sistema real.
Visando o desenvolvimento de pesquisa e aprimoramento no processo de
destilação, o professor Luiz Paulo Miranda Vaillant, do curso de Química da então
Escola Técnica Federal de Campos, hoje CEFET-Campos, projetou uma coluna de
destilação para fins didáticos. Em 1989, como resultado de um intercâmbio entre a
Refinaria Nacional de Sal – Sal Cisne – localizada em Cabo Frio – RJ e a instituição,
o curso técnico de Instrumentação recebeu como doação a coluna de destilação
(Crespo, 2000).
De um modo geral as destilações na prática industrial envolvem misturas
multicomponentes, porém a compreensão dos princípios da destilação de misturas
binárias é particularmente importante, pois ela constitui a base de operações mais
complexas (Gomide, 1988).
Para efeito didático usa-se neste trabalho uma mistura binária composta por
água e etanol comercial (álcool hidratado) a uma graduação em torno de 6 a 10o GL.
A coluna de destilação, em questão neste trabalho, é composta das colunas A
e B sobrepostas, diferenciando-a das destilarias industriais convencionais (Figura 3).
Figura 3 - Coluna de destilação piloto (Crespo, 2000).
10
Diferente dos sistemas industriais a coluna de destilação em questão não
apresenta subsistemas como referverdor e refluxo de topo5. A fonte térmica é uma
resistência elétrica de 1000W, alojada no interior da base da coluna. Todo o
destilado obtido se condensa e é coletado em uma proveta, portanto a temperatura
de topo é controlada totalmente através da vazão na entrada, uma vez que a
resistência elétrica não varia.
Na Figura 4 é apresentada a coluna de destilação piloto instalada e com a
sua instrumentação.
Figura 4 - Coluna de destilação piloto instalada (Naegele, 2000).
Como forma de reproduzir em laboratório o ambiente encontrado na indústria,
a coluna de destilação possui um sistema de aquisição de dados de porte industrial.
Os instrumentos de medição das variáveis do processo são interligados em uma
rede Foundation Fieldbus6 - vide Figura 5. Através de uma DFI (Device Fieldbus
Interface) a rede de instrumentos se integra a uma rede ponto-a-ponto com o
computador, onde são disponibilizadas as informações oriundas do processo e estas
podem ser monitoradas e modificadas através do Syscon®7.
5
Retorno de parte do produto destilado, para favorecer a produção de álcool e controlar a
temperatura.
6
Sistema de comunicação de instrumentos em rede, aplicação industrial.
7
Software de configuração da rede fieldbus.
11
Figura 5 - Instrumentos interligados em rede fieldbus.
Através de uma integração entre os seguintes softwares: Syscon®, OPC link®
e InTouch 6.0®8, um programa foi criado para registrar os valores numéricos das
variáveis ao longo do tempo em arquivo texto. Este arquivo é manipulado
posteriormente para que seja montada a matriz de dados para formar os conjuntos
de estimação e validação utilizados na modelagem.
Existem duas diferenças importantes entre, a coluna didática e uma coluna
industria: a primeira está no sistema de aquecimento da mistura. Em colunas
industriais o aquecimento faz-se em um equipamento chamado de refervedor, cuja
função é vaporizar a mistura oriunda do fundo da coluna através de transferência de
calor, realizada por vapor de água que passa através dos tubos. O aquecimento da
mistura é realizado no interior da base, através de uma resistência elétrica de aço
inoxidável e potência máxima de 1000 W sob uma tensão de 127 V;
A segunda está no sistema de retorno do destilado à coluna, chamado de
refluxo, cuja finalidade é enriquecer a composição do destilado de forma a obter-se
um álcool entre 92 a 96º GL. A coluna didática também não dispõe deste sistema.
2.3 – CARACTERÍSTICAS FÍSICAS
A coluna de destilação é constituída de 5 gomos de 150 mm de altura e 150
mm de diâmetro
cada e 10 pratos perfurados sem calotas, funcionando como
coluna de destilação e uma coluna de retificação, sobreposta à primeira, com 5
8
Software de supervisão de plantas de processo na indústria, propriedade da Invensys Systems, Inc.
12
gomos de 80 mm de altura e 150 mm de diâmetro cada e 10 pratos perfurados sem
calotas, funcionando como coluna de retificação e uma base, também chamada de
panela, com 200 mm de altura, perfazendo 1370 mm de altura total e 150 mm de
diâmetro, trocadores de calor, instrumentos industriais para a medição e atuação do
processo, tubulações, tanques de preparação e armazenamento da mistura binária e
sistema digital de aquisição de dados. O sistema de instrumentação atual, instalado
na coluna de destilação pode ser visto na Figura 6.
Figura 6 - Coluna de destilação piloto com a instrumentação atual.
Os trocadores de calor tipo “casco e tubo”, projetados e construídos em aço
inoxidável 316, possuem, cada um, 300 mm de comprimento e 60 mm de diâmetro
do casco e 6 mm de diâmetro do tubo, foram projetados para serem utilizados como
aquecedores da mistura que será destilada e como condensadores dos vapores
alcoólicos que saem do topo da torre e podem ser vistos na Figura 7.
13
Figura 7 - Trocadores de calor da coluna de destilação.
As conexões, tanto da torre quanto dos trocadores de calor, são de ¼ de
polegada. Por esse motivo, o diâmetro de ¼ de polegada foi adotado como medida
padrão de todas as tubulações metálicas. A base da coluna de destilação é onde ser
deposita a mistura a ser destilada – vide Figura 8.
Figura 8 - Base da coluna de destilação.
14
Na Figura 9 é mostrado o ponto de alimentação da coluna de destilação, e
sua válvula de controle de vazão.
Figura 9 - Ponto de alimentação da coluna de destilação.
O sistema de aquecimento da mistura é feito na base da coluna através de
uma resistência de 1000W de potência (2X500W) em 127V de tensão alternada. A
resistência foi fabricada em aço inoxidável para evitar a corrosão devido à presença
de álcool na mistura e da temperatura de uso.
Para a preparação e armazenamento de mistura foram utilizadas bombonas
plásticas de 50 litros cada – vide Figura 10. São interligadas duas a duas através de
tubos e conexões plásticas, e a graduação alcoólica desejada obtida através do
acréscimo de água ou álcool na mistura. Ao lado da parte inferior de uma delas foi
instalada uma bomba elétrica de corpo plástico, similar às utilizadas em máquinas de
lavar roupas, com a função de transferir toda a mistura para o tanque de
armazenamento da mistura, construído da mesma forma que o outro. Este tanque foi
instalado num suporte em madeira e com uma altura de 800 mm acima do ponto de
alimentação, de forma a permitir o fluxo da mistura por gravidade do tanque para
coluna.
15
Figura 10 - Sistema de armazenamento da mistura binária.
O processo de destilação possui dois circuitos líquidos distintos: um
constituído de mistura hidroalcoólica, tanto com graduação baixa (no caso da
mistura e do refugo) como graduação alta (no caso do destilado) e outro de água
para refrigeração. Desta forma, utilizou-se tubos e conexões em aço inoxidável para
o primeiro caso, principalmente devido ao alto poder de corrosão e da temperatura
do álcool, e tubos e conexões de cobre, para o segundo caso.
2.4 – O PROCESSO DE OBTENÇÃO DE ÁLCOOL HIDRATADO NA COLUNA
PILOTO
O processo de obtenção de álcool hidratado na coluna de destilação pode ser
descrito da seguinte forma:
a) A mistura binária a ser destilada é produzida no tanque de preparo de
mistura, com percentual de álcool que se deseja obter. Tal percentual é medido
através de análise química da mistura;
b) A mistura é bombeada para os tanques de armazenamento de mistura
binária;
c) A mistura é então enviada para coluna através de uma tubulação. Os
tanques de armazenamento estão instalados numa plataforma 800 mm acima do
ponto de entrada na coluna. Portanto, a entrada da mistura é feita por gravidade.
16
Antes de ser admitida na coluna, porém, a mistura é aquecida nos tubos internos de
dois pré-aquecedores, que recebem o refugo que será descartado. Esse refugo, na
mesma temperatura da base, funciona como meio quente no processo. A mistura
passa pelo medidor magnético de vazão;
d) A mistura pré-aquecida entra na coluna através da entrada do prato
número 10, equivalente ao topo da coluna de destilação. Neste instante inicia-se o
processo da destilação, pois na base da torre está instalada a resistência elétrica
para gerar calor e permitir que a temperatura na base fique em torno de 99 ºC. Dois
fluxos internos são estabelecidos: um, partindo da base da coluna em forma de
vapores que vão sendo enriquecidos na medida em que ascendem até o topo e
outro em forma de líquido que, iniciado no topo, devido à temperatura do topo não
ser suficiente para a vaporização do componente água, desce tornando-se cada vez
mais pobre do componente álcool e aumentando o ponto de ebulição, até que chega
à base, praticamente isenta de álcool.
O álcool inicialmente destilado na coluna de destilação é desidratado na
coluna de retificação, saindo sob forma de vapor rico em álcool no topo, este é
enviado através da tubulação ao condensador. Esse condensador recebe no casco,
em contra-fluxo, água na temperatura ambiente. Dessa forma, o álcool é
condensado e enviado ao recipiente de destilado.
O nível desejado para a base da coluna é de 75% da faixa de medição, sendo
de grande importância este controle que objetiva evitar que a resistência elétrica
fique fora do meio líquido, o que provocaria seu rompimento. O valor desejado
descoberto empiricamente para a temperatura de topo é de 78 ºC.
2.5 – ALGUNS PROBLEMAS DO SISTEMA
Segundo Crespo (2000):
a) A pressão no topo da coluna foi extremamente crítica quanto à
regularidade.
Observou-se uma variação muito rápida, de uma pressão relativa positiva de
200 mm H2O a uma pressão relativa negativa em torno de 80 mm H2O. Este fato foi
resultante das interrupções da fonte térmica, que causaram a condensação do álcool
em forma de vapor no último prato (topo da torre, prato número 20), fazendo-o
descer para os pratos anteriores e conseqüentemente causando vácuo;
17
b) A vazão de entrada, por conseguinte, variou de acordo com a pressão
interna da coluna, visto que a alimentação foi feita por gravidade;
c) A temperatura na base da coluna manteve-se num valor abaixo de 99,5 ºC
durante a fase em que a temperatura da vazão de entrada equilibrou-se num valor
em torno de 55 ºC. Isso é explicável pela mistura ter entrado no processo já
aquecida.
Assim, o refervedor e o refluxo são de grande importância na destilação em
processos contínuos. O primeiro, responsável pela maior capacidade volumétrica da
mistura na base, permitindo um tempo maior de vaporização do álcool, resultando
em um rendimento maior e em um menor teor de álcool no refugo.
O segundo, responsável pelo aumento da concentração do álcool e pela
regularidade da temperatura no topo da coluna.
Os problemas acima mencionados foram solucionados durante a pesquisa de
Crespo (2000). Quanto ao item C, a injeção da mistura passou a ser feita em
temperatura ambiente.
2.6 – MODELAGEM DE COLUNAS DE DESTILAÇÃO
O modelo matemático de um sistema é definido como um conjunto de
equações usado para representar o sistema físico (Phillips e Harbor, 19969, apud
Crespo, 2000).
Do ponto de vista da Engenharia de Controle, define-se modelagem como:
A atividade de representação dos principais fenômenos que ocorrem no
processo, por equações e correlações entre suas variáveis mais significativas e que
tenham um papel importante nos projetos de suas malhas de controle (Valdman,
199910, apud Crespo, 2000).
A modelagem fenomenológica usa leis físicas e correlações para descrever o
sistema. Um processo pode ser caracterizado por suas variáveis de estado que
descrevem a quantidade de massa, energia e momento linear do sistema. As
variáveis típicas que são escolhidas são: posições e velocidades (sistemas
mecânicos), tensões e correntes (sistemas elétricos), níveis e vazões (sistemas
9
Phillips, C. L., Harbor, R. D. (1996) Sistemas de controle e Realimentação. São Paulo: Makron
Books, 558p.
10
Valdman, B. (1999) Dinâmica e Controle de Processos. 1. ed. Santiago: Belkis valdman, 216p.
18
hidráulicos), e temperaturas, pressões e concentrações (sistemas químicos, térmicos
e de reações). A relação entre os estados é determinada usando-se balanços
(princípios de conservação) de momento linear, massa, energia e também outras
equações constitutivas (correlações) (Seborg e Mellichamp, 198911, apud Campos e
Teixeira, 2006).
Como alternativa para modelagem de sistemas sem a exigência da
modelagem analítica, surge a identificação de sistemas.
A identificação de sistemas se propõe a obter um modelo matemático que
explique, pelo menos em parte e de forma aproximada, a relação de causa e efeito
presente nos dados (Aguirre, 2000).
Em linhas gerais existem duas metodologias para a identificação de modelos
em sistemas físicos e químicos.
1-Metodologia Empírica  Conhecida como identificação experimental de
processos, consiste em submeter à entrada do processo a uma perturbação de um
dos tipos:
a)
degrau12
b)
pulso13
c)
senoidal14
analisar as respostas ao teste escolhido, determinando as equações e correlações
do modelo.
2-Metodologia Analítica  É baseada no levantamento de equações
diferenciais e algébricas que compõem o modelo, é fundamentada nas leis da física,
química e físico-química. Esta metodologia analítica também é denominada de
modelo matemático ou fenomenológico do sistema.
Segundo Campos e Teixeira (2006), a identificação do processo inclui os
seguintes passos:
 Planejamento e execução experimental;
 Seleção da estrutura do modelo (linear ou não);
 Estimação dos parâmetros do modelo e
 Validação do modelo.
11
Seborg, Melichamp (1989) Process Dynamics and Control, Wiley.
Variação de amplitude de um sinal entre patamares e de forma abrupta.
13
Variação de amplitude de um sinal partindo de um patamar, alcançando outro e retornando de
forma imediata.
14
Variação de amplitude de um sinal descrevendo trajetória semelhante a do gráfico da função seno.
12
19
Aguirre (2000) é mais completo e apresenta a identificação de sistemas com
mais uma etapa:
 Testes dinâmicos e coleta de dados;
 Escolha da representação matemática;
 Determinação da estrutura do modelo;
 Estimação dos parâmetros e
 Validação do modelo.
Sistemas e colunas de destilação são foco de pesquisas na área de
modelagem e identificação de sistemas. A motivação vem do fato de serem sistemas
multivariáveis e não lineares, o que requisita do modelador capacidade técnica e
metodologia aprimorada.
Hovd, et al. (1997) projetam e desenvolvem uma estratégia de controle
preditivo baseada em modelo para uma coluna de destilação a vácuo. Para a
modelagem escolhem uma representação no espaço de estados, considerado semirigoroso, pois as restrições do sistema real inviabilizam uma modelagem mais
rigorosa, sobretudo pela gama variável de produtos produzidos pela unidade de
refino.
Schneider, et al. (2001) realizam a modelagem dinâmica e simulação de uma
coluna de destilação reativa de lotes. As relações dinâmicas das equações de
transporte foram modeladas através de equações diferenciais, e incluem a
transferência de massa, calor e da reação química. Foram utilizados para validação
do modelo dados de uma coluna de destilação piloto.
Doma, et al. (2001) utilizam uma metodologia de modelagem em malha
fechada, ou seja, a dinâmica do sistema é modelada enquanto o mesmo é
controlado. As vantagens são significativas uma vez que a perturbação na forma de
degrau em malha aberta leva o sistema a condições operacionais improdutivas. O
modelo identificado é do tipo MIMO (Multiple-Input and Multiple-Output), e foi
utilizado como base para um controlador preditivo para controle de uma coluna de
destilação de uma refinaria de petróleo.
Belincanta (2004) realizou testes dinâmicos em uma coluna de destilação do
tipo etanol-água com o objetivo de modelar as relações existentes entre velocidade
superficial do vapor, a fração de área livre do escoamento e a concentração de
etanol nas condições hidrodinâmicas e de eficiência do prato. Seus resultados
demonstram que a velocidade superficial é muito influenciada pela fração de área
20
livre do escoamento e apresenta uma correlação empírica para a predição da altura
da dispersão em função das propriedades físicas e operacionais da planta.
Moraes (2004) apresenta uma abordagem com modelos auto-regressivos
lineares com variáveis exógenas (ARX) para modelar os diversos subsistemas de
uma unidade de fracionamento de Nafta, que grosso modo, nada mais é que uma
coluna de destilação.
Marangoni (2005) estudou o comportamento de uma coluna de destilação e
realizou testes dinâmicos para levantamento do modelo em regime estacionário da
planta. A partir de uma perturbação na composição da mistura de alimentação da
planta, foram calculados os parâmetros de funções de transferência contínuas para
sistema de primeira ordem com tempo morto. A partir destes modelos uma
estratégia de controle distribuída foi implementada e seus resultados expostos.
Zanata (2005) realizou testes dinâmicos em uma coluna de destilação de
petróleo para obtenção de um modelo dinâmico que relacionasse temperatura nos
pratos com a composição do produto de topo da coluna. O modelo utilizado para
ajuste aos dados foi uma rede neural artificial.
21
CAPÍTULO III
FERRAMENTAL TEÓRICO
3.1- MODELOS AUTO-REGRESSIVOS
A metodologia proposta neste trabalho faz uso de representações lineares,
para a modelagem de sistemas dinâmicos, através de equações de recorrência
conhecidas como modelos auto-regressivos.
Estes modelos foram utilizados inicialmente por pesquisadores da área de
modelagem de séries temporais tendo como principais divulgadores Box e Jenkins
(1994) que apresentaram o modelo ARIMA. Eles utilizaram os modelos ARIMA para
solução de problemas de previsão de séries temporais e controle.
Os pesquisadores da área de identificação de sistemas passaram a utilizar
variações dos modelos auto-regressivos para modelagem de sistemas dinâmicos da
indústria. Dando origem a novas configurações de modelos auto-regressivos que
incluíam no modelo séries temporais de diversas variáveis, componentes para
modelagem de não linearidades do sistema e outras variações. Estes modelos são
também conhecidos como modelos de predição de erro, devido ao critério de
seleção baseado no erro de predição do modelo obtido.
Estes modelos são aplicados sobre seqüências de observações de uma ou
mais variáveis do sistema em estudo. Este conjunto de observações da variável
pode ser representado por um sinal yt   y t , yt 1 , , y t  k . Onde t é o instante
referente a um dado valor de y e k é o valor do atraso considerado para o sinal.
Os modelos auto-regressivos possuem diversas variações e formatos em que
podem ser apresentados. Abaixo são apresentadas algumas destas variações em
sua forma polinomial compacta e estendida.
O modelo auto-regressivo AR é o mais comum destes modelos paramétricos
e sua formulação é a seguinte:
A(q 1 ) y (t )  e(t )
O operador q  na toma o valor da função y t  em um instante anterior
yt  n a  . Ao selecionar a ordem n a do operador é determinado o número de valores
22
atrasados de y t  serão utilizados para determinar o valor no instante atual. O erro
cometido ao se tentar modelar yt  em função dos seus valores atrasados é dado
por et .
Desta forma:
q 1 y t   y t  1
q 2 y t   y t  2 
q  na yt   yt  n a 
O termo
A(q 1 ) é um polinômio na variável q e possui a seguinte
representação:
A( q 1 )  1  a1 q 1  ...  a na q  na
Sendo assim a versão expandida do modelo AR fica sendo:
y(t )  a1 y (t  1)  ...  a na y (t  na )  e(t )
Se, ao modelo AR for acrescentado uma entrada externa u (t ) , que nada mais
é que uma série temporal de uma variável tida como entrada que também deve ser
admitida como parte que explica o comportamento de y(t ) , ou seja, o valor de y(t )
não pode ser explicado somente pelos seus valores atrasados, mas também pelos
valores atrasados de outra variável. Chega-se ao seguinte modelo ARX:
A( q 1 ) y(t )  B ( q 1 )u (t  nk )  e(t ) , onde:
B(q 1 )  b1  b2 q 1  ...  bnb q  nb 1
O modelo ARX expandido seria:
y (t )  a1 y(t  1)  ...  a na y (t  n a )  b1u (t  n k )  b2 u (t  n k  1)  ...  bnb u (t  nk  nb  1)  e(t )
onde na e nb são as ordens dos polinômios A( q 1 ) e B ( q 1 ) , enquanto que nk é o
número de atrasos da saída para a entrada. Se acrescentada uma componente de
média móvel (MA) ao modelo passamos a um modelo ARMAX, descrito por:
A( q 1 ) y (t )  B ( q 1 )u (t  n k )  C (q 1 )e(t ) , onde:
C (q 1 )  1  c1q 1  ...  cnc q  nc , é possível notar que o polinômio C (q 1 ) é similar
ao polinômio A( q 1 ) .
Aplicando os polinômios expandidos no modelo ARMAX, tem-se:
y(t )  a1 y (t  1)  ...  a na y (t  na )  b1u (t  n k )  b2 u (t  n k  1)  ...  bnb u (t  n k  nb  1) 
 e(t )  c1e(t  1)  ...  c nc e(t  n c )
23
Existem ainda outras variações dos modelos auto-regressivos que possuem
denominações próprias tais como os modelos OE (Output Error) e BJ (Box &
Jenkins).
A estrutura do modelo OE é a seguinte:
y (t ) 
B(q 1 )
u (t  nk )  e(t ) , onde:
F (q 1 )
F (q 1 )  1  f 1q 1  ...  f n f q
n f
, similar ao polinômio A( q 1 ) ,
a estrutura do modelo BJ é a seguinte:
y (t ) 
B(q 1 )
C (q 1 )
u
(
t

nk
)

e(t ) , onde:
F (q 1 )
D(q 1 )
D( q 1 )  1  d 1q 1  ...  d nd q  nd , similar ao polinômio A( q 1 )
De uma forma geral, existe um modelo que representa qualquer um dos
modelos auto-regressivos anteriormente mencionados, ou seja, qualquer modelo
linear e discreto no tempo.
A(q 1 ) y (t ) 
B(q 1 )
C (q 1 )
u
(
t

nk
)

e(t ) , podendo ser utilizado para todos os
F (q 1 )
D(q 1 )
casos anteriores, sendo que estes seriam casos especiais onde, as ordens dos
polinômios desnecessários seriam zeradas da seguinte forma:
Tabela 1 - Tabela de polinômios retirados do modelo geral.
Modelo
Grau do polinômio a ser zero
AR
B ( q 1 ), F (q 1 ), C ( q 1 ), D (q 1 ) ou seja, nb , n f , nc e nd  0
ARX
F ( q 1 ), C ( q 1 ), D( q 1 ) ou seja, n f , nc e nd  0
ARMAX
F ( q 1 ), D( q 1 ) ou seja, n f e n d  0
OE
A(q 1 ), C (q 1 ), D(q 1 ) ou seja, n a , nc e n d  0
BJ
A(q 1 ) ou seja, na  0
Dependendo do sistema a ser modelado ele pode possuir mais de uma
entrada externa gerando a seguinte estrutura:
Bnu (q 1 )
B1 (q 1 )
C (q 1 )
A(q ) y (t ) 
u
(
t

n
)

...

u
(
t

nk
)

e(t )
1
k1
n
nu
F1 (q 1 )
Fnu (q 1 ) u
D (q 1 )
1
Demais estruturas podem ser encontradas na literatura, mas para este
trabalho explanar sobre estas, é o suficiente.
24
3.1.1- ESTADO DA ARTE DA APLICAÇÃO DE MODELOS AUTO-REGRESSIVOS
A modelagem de colunas de destilação através de modelos auto-regressivos
é consolidada sobre pesquisas relacionadas a diversas áreas de conhecimento, tais
como: sistemas dinâmicos, engenharia de controle e processamento de sinais,
sempre com intuito de obtenção do modelo dinâmico do processo para fundamentar
análises e projetos que venham a melhorar seu desempenho em operação.
Sob a ótica da pesquisa operacional, uma coluna de destilação pode ser
enxergada como um problema de otimização, dotado de função objetivo e restrições.
Encontrar as relações implícitas aos dados oriundos deste sistema é dar suporte
para a modelagem da função-objetivo do problema. A técnica se sustenta no
potencial dos modelos auto-regressivos propostos para a modelagem de sistemas
dinâmicos multivariáveis.
Aguirre et al. (1998) aplicam modelos NARMAX (Non linear Auto-Regressive
Moving Average with eXogenous inputs), para modelar um forno elétrico com
elemento de aquecimento interno de 200 watts.
Barroso et al. (2002) apresentam a identificação de um conversor CC-CC
15
BUCK , que possuía faixa restrita para os valores de entrada.
Furtado et al. (2002) experimentam modelos NARMAX para identificação de
um forno a arco elétrico, modelados e gerados os dados através de simulação,
apenas para verificação da eficácia do método.
Jurado e Carpio (2003) aplicam os modelos NARX para identificação de uma
micro-turbina composta pelos demais subsistemas: combustível, controle de
velocidade e de temperatura.
Moraes (2004) apresenta uma abordagem com modelos auto-regressivos
lineares com variáveis exógenas (ARX) para modelar os diversos subsistemas de
uma unidade de fracionamento de Nafta, que grosso modo, nada mais é que uma
coluna de destilação.
Bravo (2006) utiliza modelos ARX para identificação da resposta dinâmica de
um motor DC, utilizando dois softwares distintos: MATLAB®16 e LABVIEW®17.
Júnior e Barbabela (2006) utilizam modelos auto-regressivos para predição de
valores futuros dos preços de derivados de petróleo, onde modelavam séries
15
Dispositivo eletrônico utilizado para reduzir a tensão aplicada na entrada.
MATrix LABoratory – software proprietário da THE MATWORKS INC.
17
Software de simulação computacional, propriedade da National Instruments.
16
25
temporais com histórico superior a dez anos e através de correlação cruzada,
definiram a ordem dos modelos auto-regressivos.
Carvalho e Guillermo (2007) ajustam um modelo auto-regressivo do tipo ARX
a partir de dados de uma coluna de destilação piloto do CEFET-Campos e validam o
modelo com dados de cinco testes dinâmicos realizados no sistema real.
3.2- ESTIMAÇÃO DE PARÂMETROS
Os modelos auto-regressivos devem ter seus coeficientes estimados a partir
de métodos específicos de estimação de parâmetros, pois somente depois de
estimados os coeficientes do modelo a partir dos dados é que se tem o modelo
capaz de representar o comportamento do sistema identificado.
Estimação pode ser definida como a forma de obtenção do provável valor de
determinado parâmetro, coeficiente ou variável. Todo valor estimado possui um erro
implícito.
A estimação de parâmetros é utilizada em diversas áreas de pesquisa, mas
em todas há necessidade de obtenção de modelos mais próximos do sistema real.
Para distribuição de probabilidades de variáveis aleatórias utiliza-se a estimação
para obtenção de seus coeficientes. Exemplos destas distribuições são a normal,
weibull, exponencial, etc.
Modelos baseados nas relações físicas que governam os sistemas podem ter
coeficientes estimados. No presente trabalho os coeficientes estimados pertencem a
uma classe de modelos pré-formatados e os valores destes coeficientes determinam
o quanto este modelo pode explicar o sistema modelado.
26
3.2.1- MÉTODO DOS MÍNIMOS QUADRADOS
Figura 11 - Método dos mínimos quadrados.
O método dos mínimos quadrados baseia-se na estimação de parâmetros de
um modelo a partir da minimização da função de erro existente entre o valor da
saída do modelo f xi  para um dado valor de entrada xi e o valor da saída
observada yi obtida por amostragem do sistema real. É aplicado tipicamente
quando se deseja aproximar um modelo padrão devido ao desconhecimento da
função real.
Trata-se de obter valores para a e b que façam com que f xi  y i seja o
menor possível para todo xi . Isto é conseguido através da minimização da função
    f xi  y i 2 sendo que esta é uma função que depende de a e b ,
 a, b    axi  b  y i 2 .
Para minimizar a função  , faz-se:


0 e
0
a
b
É provável que a reta não consiga passar por todos os pontos, portanto o
método não garante que  seja zero, mas que seja o menor possível.
Dos modelos auto-regressivos apresentados na seção anterior, o escolhido
para ajustar-se aos dados do sistema foi o ARX (auto-regressivo com entradas
externas).
27
O modelo ARX da seção anterior pode ser visto como um sistema de
equações e ter seus coeficientes estimados através do método dos mínimos
quadrados. O método pode ser resumido da seguinte forma.
A( q 1 ) y(t )  B (q 1 )u (t  n k )  e(t ) , onde:
A(q 1 )  1  a1q 1  ...  ana q  na
B ( q 1 )  b1  b2 q 1  ...  bnb q  nb 1
O modelo ARX expandido seria:
y (t )  a1 y (t  1)  ...  a na y (t  n a )  b1u (t  n k )  b2 u (t  n k  1)  ...  bnb u (t  n k  nb  1)  e(t )
Modificando sua forma para aplicação do método dos mínimos quadrados:
y(t )   T (t  1)  e(t )
 T (t  1)   y(t  1), , y (t  n a ), u (t  1),, u (t  1  nb )

 T  a1 , , a na , b0 , , bnb

Onde  T (t  1) é a matriz linha de valores coletados do processo para as
variáveis envolvidas no modelo, neste caso considera-se um modelo relacionando
uma variável de saída y(t ) com uma variável de entrada u (t ) . A matriz coluna de
coeficientes a serem estimados através do método dos mínimos quadrados  . Um
erro obtido na modelagem é representado por e(t ) .
O método de mínimos quadrados é largamente utilizado para estimação de
parâmetros do modelo de sistemas lineares discretos no tempo.
De início será assumido que se conhece o valor estimado do vetor de
parâmetros,  , e que é cometido um erro e ao se tentar explicar o valor observado
de y a partir do vetor de regressores  T e de  , ou seja,
y   T  e
A partir daí o que se tem a fazer é determinar o índice de desempenho do
método, sendo que J MQ é a função custo dependente de  a ser minimizada:
N
J MQ    (i ) 2   T   
2
i 1
J MQ quantifica a qualidade de ajuste de
 T  ao vetor de dados y . Portanto é
tido como vantajoso que o valor estimado de  minimize J MQ .
28
Determinando-se e
N
em
y   T   e e substituindo-se o resultado em
2
J MQ   e(i) 2  e T e  e , tem-se:
i 1
J MQ  ( y   T  ) T ( y   )  y T y  y T    T  T y   T  T 
A fim de minimizar a função custo J MQ com respeito à  , é necessário
resolver J MQ /   0 . Fazendo-se isso, tem-se:
J MQ

 ( y T  ) T   T y  ( T    T  )   T y   T y  2 T 
Igualando-se a última equação a zero tem-se:
  [ T  ]1 T y
Para que  seja o único ponto de mínimo, é necessário verificar que:
 2 J MQ

2
 2 T   0
A inequação é verdadeira, pois 2 T  é positiva definida por construção
(Aguirre, 2000). Portanto, a equação   [ T  ] 1  T y é o estimador que fornece o
valor de  que minimiza o somatório dos quadrados dos erros. Resumindo:
 MQ  arg  min J MQ  [ T  ] 1  T y
Caso o estimador de mínimos quadrados seja adequado, o vetor de erros e
será do tipo ruído branco, com média nula e variância  e2 .
29
CAPÍTULO IV
METODOLOGIA
4.1- INTRODUÇÃO
Segundo Box e Jenkins (1994) a metodologia para a identificação de sistemas
através de modelos auto-regressivos passa pelas seguintes etapas:
Identificação: é a etapa de definição dos filtros mais adequados para a série
analisada. Utilizar filtros em excesso na descrição de uma série não promove
previsões melhores, além de complicar desnecessariamente seu processamento
matemático.
Estimação: escolhidos os filtros, estimam-se as ordens e os parâmetros que
compõem os modelos. Novamente, é importante manter o modelo em sua forma
mais simples, com a menor ordem em cada filtro. Elevar a ordem, ou seja, usar
dados referentes ao passado distante que não exercem influencia sobre o
comportamento presente não melhora o desempenho do modelo.
Diagnóstico: determinar a acurácia do modelo estimado através da
comparação dos dados previstos com os dados obtidos a partir da série real. Nessa
etapa é realizada avaliação para verificar se o modelo é capaz de prever de forma
segura o comportamento futuro do mercado. Caso não seja obtido sucesso no seu
uso como “preditor”, retornam-se as etapas anteriores para a busca de um novo
modelo.
Aplicação: após as etapas anteriores terem obtido êxito, a ferramenta é
colocada em uso.
Ljung (1994) sugere o ciclo necessário para realizar a identificação de
sistemas dinâmicos através de ready-made models, que são a classe à qual
pertencem os modelos auto-regressivos – vide Figura 12.
30
Figura 12 - Ciclo de identificação com auxílio computacional. Os retângulos são de responsabilidade
do computador e as elipses responsabilidade do usuário (Adaptado de Ljung (1994)).
Para aplicação desta metodologia a um sistema real faz-se necessário um
maior detalhamento, e explicitar as etapas seguidas por pesquisadores que
utilizaram-na para identificação de sistemas similares.
4.2 – ESCOLHA DAS VARIÁVEIS DO MODELO
Como primeiro passo para a modelagem da coluna de destilação, pode-se
considerar a escolha das variáveis do modelo. Existem diversos critérios utilizados
para esta etapa, mas em todos os casos por se tratar de um sistema real, deve-se
verificar a viabilidade da medição destas variáveis ao longo do tempo, ou seja, estes
valores devem ser amostrados, através de algum sistema de aquisição de dados.
Neste ponto surge a primeira forma de seleção das variáveis envolvidas no
modelo, pois se a modelagem será realizada a partir de dados obtidos do sistema,
só faz sentido selecionar as variáveis que podem ter seus valores amostrados pelo
sistema de aquisição de dados da planta. São elas: vazão de entrada, nível da base,
temperatura da base, pressão do topo e temperatura do topo. Portanto, o modelo a
31
ser ajustado deve conter no máximo estas cinco variáveis cujas medições estão
disponíveis.
Esta restrição do processo a ser modelado, é uma primeira seleção das
variáveis do modelo, porém é considerada mais como restrição do sistema do que
como critério de seleção.
Diversos métodos de seleção das variáveis a serem consideradas na
modelagem tem sido propostos em trabalhos científicos.
A fim de se reduzir o número de variáveis selecionadas para o treinamento da
rede, fez-se um estudo comparativo baseado na covariância das temperaturas de
cada prato com relação à variável de saída que se quer estimar (composição na
saída da torre) (Zanata, 2005).
Esta é uma outra forma de selecionar as variáveis para o modelo, chama-se
covariância cruzada e pode ser obtida através da equação:
Cov x, y k  
1 N
  xt  x  y t  k  y 
N t 1
onde: N é o número de elementos dos vetores;
x e y são respectivamente os vetores de entrada e saída;
x e y são os valores médios dos respectivos vetores;
k é o número de atrasos analisados para o sinal.
A figura 13 apresenta o gráfico da covariância cruzada entre as temperaturas
dos estágios (pratos) e a composição na saída da torre.
32
Figura 13 - Covariância cruzada das temperaturas nos estágios e a composição (Zanata, 2005).
Uma das formas mais difundidas é a utilização da função correlação cruzada.
A função de correlação cruzada apresenta-se como uma ferramenta capaz de
determinar se há correlação significativa entre duas variáveis candidatas a compor
um modelo (Aguirre, 2000).
Segundo Lara (2005) a função de correlação cruzada discreta entre dois
sinais x t  e yt  com a mesma propriedade é definida por:
1
N  N
Rxy k   lim
N
 xt  yt  k 
t 1
onde: N é o número de elementos dos vetores;
x e y são respectivamente os vetores de entrada e saída;
k é o número de atrasos analisados para o sinal.
Sendo assim as variáveis possíveis de mensuração pelo sistema de aquisição
de dados do processo podem ser selecionadas pela significância da correlação que
exista entre elas. Um método bem prático da utilização deste critério de seleção é
visto em Aguirre (2000).
Neste exemplo foram utilizados dados reais coletados na unidade de águas
ácidas de uma refinaria de petróleo. Para a seleção das variáveis envolvidas na
modelagem, um primeiro critério foi utilizado para selecionar a variável de saída do
33
modelo. Este critério foi o objetivo da modelagem. Desejava-se obter um modelo
matemático que explicasse a temperatura no topo da torre de águas ácidas, portanto
ela foi selecionada para a modelagem como sendo a única variável de saída do
modelo.
Restavam então outras sete variáveis do processo as quais deveriam ser
selecionadas ou descartadas para integrar o modelo como variáveis de entrada.
Após a seleção da variável de saída as demais foram selecionadas como
variáveis de entrada de acordo com a significância da correlação que possuíam com
a variável de saída. Para medir esta significância foi utilizado o correlograma gerado
a partir da utilização da função correlação cruzada.
Além da investigação da correlação entre cada uma das candidatas a
variáveis de entrada e a variável de saída, também foi investigada a correlação entre
elas com o intuito de evitar a inclusão de variáveis redundantes no modelo, o que
aumentaria de maneira desnecessária a complexidade.
Na Figura 14 (a) (b) e (c) são apresentados os gráficos das funções de
correlação cruzada para os casos relatados.
Figura
14 - Funções de correlação cruzada (a, b) entre a variável de saída e duas variáveis
candidatas a serem entradas, (c) entre (outras) duas variáveis candidatas a serem entradas (Aguirre,
2000).
A utilização de procedimento puramente matemático para a seleção de
variáveis de um modelo sem sentido e de difícil interpretação (Silva, 2006).
34
No objeto de estudo deste trabalho existem restrições quanto ao número de
variáveis monitoradas e possíveis de amostragem, o que influencia na seleção das
variáveis do modelo. Além disto é comum nos sistemas industriais considerar como
variáveis de entrada aquelas que podem ser manipuladas, e que provocam
modificação
no
estado
das
demais
variáveis
e
influem
diretamente
no
comportamento da variável de saída.
Também é comum serem consideradas como variáveis de saída do sistema,
as variáveis que se deseja controlar, com o objetivo da continuidade operacional e
produtiva do sistema. No caso específico deste trabalho estas considerações levam
a
inclusão
das
variáveis:
vazão
de entrada
e
temperatura
de
topo,
respectivamente como entrada e saída do modelo.
Desta forma cabe salientar que restrições operacionais, o objetivo da
modelagem e heurística devem nortear a formulação de modelos matemáticos que
serão utilizados em sistemas de qualquer natureza inclusive os processos
industriais.
O método de utilização da correlação cruzada entre as variáveis mesuradas
foi utilizado para determinação das variáveis que devem fazer parte do modelo. Os
dados utilizados para o método são referentes ao primeiro teste dinâmico realizado
na planta e que será usado para a estimação de parâmetros do modelo, e podem
ser vistos na Figura 15.
35
Comportamento das variáveis ao longo do tempo
100
90
80
Valores das variáveis
70
60
50
40
temp topo
vazao
temp base
pressao
nivel base
30
20
10
0
0
500
1000
1500
2000
2500
Tempo em segundos
3000
3500
4000
Figura 15 - Dados utilizados para a seleção das variáveis do modelo.
Baseado no objetivo da modelagem, a temperatura de topo foi definida como
variável de saída do modelo, ou seja, ela será determinada em função das demais
variáveis consideradas como entradas do sistema.
Partindo desta prerrogativa, o método de seleção das variáveis para compor o
modelo segue o mesmo que em Aguirre (2000), verificação da correlação cruzada
entre as variáveis de entrada e a de saída e também entre as de entrada para
verificação de redundâncias.
Nas Figuras 16, 17, 18 e 19 são apresentados os gráficos para a correlação
cruzadas das variáveis de entrada (vazão, temperatura da base, pressão e nível
respectivamente) com a variável de saída temperatura de topo.
36
Correlação cruzada entre vazão e temperatura de topo
0.6
Valor de correlação cruzada
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-400
-300
-200
-100
0
Atrasos
100
200
300
400
Figura 16 - Correlação cruzada entre a vazão e a temperatura no topo.
Correlação cruzada entre a temperatura da base e a temperatura de topo
Valor de correlação cruzada
1
0.5
0
-0.5
-400
-300
-200
-100
0
Atrasos
100
200
300
Figura 17 - Correlação cruzada entre a temperatura na base e a temperatura no topo.
400
37
Correlação cruzada entre pressão e temperatura de topo
0.4
0.3
Valor de correlação cruzada
0.2
0.1
0
-0.1
-0.2
-0.3
-0.4
-0.5
-400
-300
-200
-100
0
Atrasos
100
200
300
400
300
400
Figura 18 - Correlação cruzada entre a pressão no topo e a temperatura no topo.
Correlação cruzada entre nivel e temperatura de topo
0.15
Valor de correlação cruzada
0.1
0.05
0
-0.05
-0.1
-0.15
-0.2
-400
-300
-200
-100
0
Atrasos
100
200
Figura 19 - Correlação cruzada entre o nível da base e a temperatura no topo.
38
Através da contemplação das correlações existentes entre as variáveis
candidatas a entrada e a variável de saída pode-se perceber que a única que não
possui correlação significativa – valor de correlação acima de 0,3 – com a
temperatura de topo é o nível da base, portanto esta variável não compõe o modelo.
Nas Figuras 20, 21 e 22 têm-se os gráficos da função correlação cruzada
entre as variáveis de entrada restantes.
Correlação cruzada entre vazão e temperatura da base
0.6
0.5
Valor de correlação cruzada
0.4
0.3
0.2
0.1
0
-0.1
-0.2
-0.3
-0.4
-400
-300
-200
-100
0
Atrasos
100
Figura 20 - Correlação cruzada entre a vazão e a temperatura da base.
200
300
400
39
Correlação cruzada entre pressão e vazão
0.4
Valor de correlação cruzada
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
-400
-300
-200
-100
0
Atrasos
100
200
300
400
Figura 21 - Correlação cruzada entre a pressão e a vazão.
Correlação cruzada entre temperatura da base e pressão
0.5
0.4
Valor de correlação cruzada
0.3
0.2
0.1
0
-0.1
-0.2
-0.3
-0.4
-0.5
-400
-300
-200
-100
0
Atrasos
100
200
Figura 22 - Correlação cruzada entre a temperatura da base e a pressão.
300
400
40
De acordo com as correlações cruzadas entre as candidatas a variáveis de
entrada do sistema, pode-se verificar que todas possuem correlação significativa
entre si, indicando redundância. Desta forma, foi escolhida a vazão como variável de
entrada do sistema, e as variáveis, temperatura da base, e pressão foram
descartadas do modelo.
Esta decisão tomou como base o fato de a vazão ser uma variável
diretamente manipulada pela abertura da válvula, e as demais serem conseqüência
de variações dela e de outras variáveis. Em diversos tipos de aplicações a que
possa ser submetido o modelo obtido neste trabalho, a possibilidade de intervenção
direta na variável de entrada do sistema é de grande valia.
4.3 – ESCOLHA DOS ESTÍMULOS DE ENTRADA
Uma parte importante da elaboração do experimento para coleta de dados é a
escolha dos sinais que serão utilizados para estimular as dinâmicas do sistema.
Trata-se de selecionar o tipo de variação realizada nas variáveis de entrada
do sistema de forma a produzir algum efeito na variável de saída. Através da
observação do comportamento da variável de saída às variações na entrada, é
possível subsidiar algumas decisões posteriores para a modelagem do sistema.
Na literatura sobre identificação de sistemas dinâmicos na indústria três testes
são muito utilizados como estímulos de entrada para o sistema: impulso, degrau,
rampa e PRBS (Pseudo Random Binary Signal).
Segundo Aguirre (2000) os sinais de entrada devem possuir as seguintes
características:
 A persistência de excitação tem interpretações numéricas e dinâmicas;
 Sinais aleatórios e pseudo-aleatórios são comuns;
 No caso de sistemas com múltiplas entradas, além de tais entradas
precisarem ser persistentemente excitantes, não devem estar correlacionadas
entre si;
 Na identificação de dinâmica linear o importante é a largura do
espectro do sinal de entrada;
 Na identificação de sistemas não lineares o perfil de amplitude torna-se
importante também.
41
O sucesso do método de identificação depende do estímulo de entrada ser
persistentemente excitante (Ljung, 1994).
Um sinal PRBS é considerado persistentemente excitante, pelo fato de
levantar as dinâmicas ocultas do sistema.
Zanata (2005) utiliza estímulos do tipo degrau e rampa, de 10 e 20 % na
quantidade de calor fornecido à coluna de destilação. Para a seleção das variáveis
do modelo foram selecionados estímulos de entrada aleatórios com amplitude de
 5% em cinco variáveis de entrada do sistema e manipuláveis.
Dallagnol (2005) utiliza um sinal PRBS constante como estímulo de entrada
para realizar a identificação linear. Ele comenta que o sinal deve ser suficientemente
pequeno para não excitar as dinâmicas não lineares presentes no sistema, conforme
o gráfico da Figura 23.
Figura 23 - Sinal aplicado na entrada do sistema para realizar a identificação linear, variável vazão de
gás (Dallagnol, 2005).
Neste trabalho a utilização de testes como rampa e PRBS foram inviáveis do
ponto de vista operacional. O sistema de manipulação da variável de entrada, vazão
de produto para a coluna, é realizado através de uma válvula de controle que recebe
o sinal de controle de um conversor de sinal FF (Foundation Fieldbus) para sinal de
corrente elétrica padronizado 4 a 20 mA. Devido a restrições técnicas do
equipamento com operação em manual os testes realizados foram do tipo degrau
(Figura 24), que corresponde a uma mudança entre dois patamares de vazão de
produto em função da abertura da válvula.
42
Vazão ao longo do tempo
13
Valor de vazão em litros por hora
12
11
10
9
8
7
6
5
4
0
500
1000
1500
2000
2500
Tempo em segundos
3000
3500
4000
Figura 24 - Estímulo na forma de degrau positivo na vazão de produto da coluna de destilação, para
o primeiro teste dinâmico.
Apesar de não ser um sinal tão persistentemente excitante quanto o PRBS, o
degrau é utilizado em identificação de sistemas por diversos trabalhos científicos,
vide Zanata (2005).
4.4 – ESCOLHA DO PERÍODO DE AMOSTRAGEM
A taxa de amostragem utilizada para o registro dos valores das variáveis do
sistema, é um importante parâmetro a ser definido para a execução dos testes
dinâmicos realizados.
Existem diversos métodos para definição da taxa de amostragem a partir do
período ou intervalo de amostragem e alguns deles estão relacionados com a
constante de tempo do sistema.
A constante de tempo é determinada a partir da reação da variável de saída a
uma perturbação na forma de degrau na variável de entrada. Considera-se o
instante em que a variável de saída começa a reagir ao estímulo até o instante em
que a mesma atinge 63% da amplitude de sua variação total.
43
Primeiro método
Aguirre (2000) propõe o seguinte método:
Assume-se que o sinal original y * t  foi registrado utilizando-se um tempo de
amostragem muito pequeno, ou seja, muito menor do que necessário. A questão
passa a ser a definição de uma taxa pela qual o sinal observado y * t  será
decimado de forma a gerar o sinal de trabalho y t , devidamente amostrado.
Deseja-se determinar    de forma que yt   y * t .
Determinam-se as seguintes funções:

k   E y


ry* k   E y * t  y * t  y * t  k  y * t 
r y * 2
*2
t  y *2 t y *2 t  k  y *2 t 
e seus primeiros mínimos, k y* k y*2 , respectivamente. O menor desses mínimos


passará a ser o valor de trabalho, ou seja, k m*  min k y* k y*2 . Deseja-se escolher 
de forma que as funções de autocovariância do sinal decimado yt   y * t 
satisfaçam
10  k m  20 (Aguirre, 2000).
Sendo que os limites, inferior e superior podem ser relaxados para 5 e 25,
respectivamente.
Zanata (2005) aplica o método e obtém um dos três períodos de amostragem
obtidos em seu trabalho, os outros dois períodos foram encontrados através de
outros dois métodos distintos.
Segundo método
Zanata (2005) apresenta um método para a escolha da taxa de amostragem
baseado em estimativa, realizando o período de amostragem como sendo um
décimo da menor constante de tempo do sistema.
Para tal, realiza estímulos na forma de degrau de +5% e -5% nas variáveis
selecionadas como entrada do sistema e observa o comportamento da variável de
saída a cada uma destas perturbações.
44
A menor constante de tempo do sistema é obtida na reação mais rápida da
variável de saída, ou seja, é aquela na qual a variável de saída se estabiliza mais
rapidamente. A resposta da variável de saída pode ser vista na Figura 25.
Figura 25 - Resposta da variável de saída (composição) a uma perturbação na variável de entrada do
sistema. Resposta com menor constante de tempo (Zanata, 2005).
Antes da aplicação do degrau: x 2 0,500  0,98245
Estabilização após o degrau: x 2 0,515  0,98365
  x 2 0,515 x 2 0,500  0,0012
Temos que 63% do   0,000756
No gráfico corresponde ao valor de 0,98321, o instante de tempo de 0,5035h.
Logo a menor constante de tempo é:  m  0,0035 h.
Finalmente, o período de amostragem escolhido é: T  0,0004 h ou 1,5s
(Zanata, 2005).
Neste trabalho o período de amostragem foi escolhido de acordo com o
segundo método exposto – (Zanata, 2005). A constante de tempo utilizada para o
cálculo do período de amostragem foi escolhida a partir dos trabalhos realizados por
Naegele (2000), eles encontraram duas constantes de tempo diferentes, uma para a
reação da temperatura de todo a uma perturbação na forma de degrau negativo na
vazão. Outra constante foi encontrada a partir de uma perturbação na forma de
degrau positivo. Seus valores são respectivamente 436 e 475 segundos. Para evitar
que o período escolhido deixasse escapar alguma dinâmica importante do sistema,
45
foi escolhida uma constante de tempo de 400 segundos e aplicando o segundo
método, o período de amostragem escolhido foi 40 segundos.
Este passo foi utilizado em uma rotina elaborada no software MATLAB®, cujo
objetivo é montar o novo conjunto de dados a partir do original que possuía registros
das variáveis com período de amostragem de 1 segundo.
4.5 – SELEÇÃO DA ESTRUTURA DO MODELO
Identificação é a etapa de definição dos filtros mais adequados para a série
analisada. Utilizar filtros em excesso na descrição de uma série não promove
previsões melhores, além de complicar desnecessariamente seu processamento
matemático (Box e Jenkins, 1994).
Os filtros citados se referem às componentes do modelo auto-regressivo que
será utilizado para a modelagem das séries temporais referentes às variáveis do
sistema. A estrutura do modelo auto-regressivo a ser utilizado vai depender da
necessidade de inclusão destas componentes.
Por exemplo, se o objetivo da modelagem é identificar a relação existente
entre a variável de saída do modelo, e as demais variáveis tidas como entradas do
sistema, este modelo auto-regressivo selecionado deve conter a componente X que
se refere a entradas exógenas ou entradas externas. Série(s) temporal (ais) da(s)
variável (eis) de entrada que explica o comportamento da variável de saída, neste
caso temperatura de topo.
Outra componente a ser acrescentada no modelo é a MA que se refere à
média móvel. É utilizada quando existe dinâmica relevante no ruído presente na
série temporal, quando o ruído possuir as características de ruído branco gaussiano,
tais dinâmicas devem ser desconsideradas e a componentes MA é descartada, o
ruído será acrescido como erro na equação. A Tabela 2 resume para os principais
modelos as circunstâncias para a inclusão de componentes.
A principal característica do modelo ARX é utilizar a entrada do sistema para
modelar a saída do mesmo. Já, o modelo ARMAX utiliza a média móvel do erro para
predizer o estado futuro (Cruz, 2006).
46
Tabela 2 - Tabela que relaciona os modelos ao tipo de situação a qual devem ser utilizados.
Modelo
Denominação
AR
ARX
Auto-Regressive
Uma única variável
Auto-Regressive with
Outras variáveis devem
eXogenous inputs
ser consideradas
Auto-Regressive Moving
ARMAX
Circunstância
Average with eXogenous
inputs
Dinâmicas relevantes no
ruído
NARX
Non linear ARX
Modelo não linear
NARMAX
Non linear ARMAX
Modelo não linear
4.6 – SELEÇÃO DA ORDEM DO MODELO
Além de selecionar a estrutura do modelo auto-regressivo faz-se necessário
determinar a ordem de tal modelo, ou seja, quantos e quais valores atrasados das
variáveis devem ser incluídos na modelagem.
A regra mais utilizada em séries temporais é o chamado critério de
informação de Akaike (Akaike, 1974), denotado por AIC. A definição mais
comumente utilizada é:


2
n   2  n ,
AIC n   N  ln  erro
2
n  é a variância do erro de
sendo que N é o número de dados,  erro
modelagem (erro de predição de um passo a frente ou resíduos) e n é o número de
parâmetros do modelo (Aguirre, 2000).
O uso do AIC pressupõe que existe uma ordem predefinida para inclusão dos
termos candidatos seqüencialmente no modelo. Primeiramente, assume-se que o
número de termos do ruído é fixo, portanto, com os termos do ruído no modelo, uma
ordem de inclusão dos demais termos seria: y(t  1), u (t  1), y (t  2), u (t  2), e assim
por diante. Quando AIC passar por um mínimo o procedimento pode ser terminado e
a ordem terá sido determinada (Aguirre, 2000).
O índice AIC n  normalmente atinge um mínimo para um determinado
número de parâmetros do modelo n* . Do ponto de vista do critério usado, esse
número é ótimo. Deve ser lembrado que tal critério é fundamentalmente estatístico e
47
não garante necessariamente que o modelo com o número “ótimo” de termos seja
um modelo válido (Aguirre, 2000).
Logo se pode inferir que apesar de utilizado o critério de informação de
Akaike, a etapa de validação do mesmo não deve ser descartada, sob pena de
utilização de um modelo que não possua uma boa resposta para novos dados.
Ljung (1994) apresenta não só o AIC como método estatístico seleção da
ordem
min
d ,
do
1 d / N 1
1 d / N N
modelo,
mas
também
o
FPE
(Final
Prediction
Error).
N
  t,  
2
t 1
Onde N é o número de dados utilizados para a estimação, d é a dimensão
do vetor de parâmetros  e  2 é a variância do erro de predição.
Os modelos resultantes de métodos estatísticos para seleção da ordem são
encontrados em diversos trabalhos científicos.
4.7 – ESTIMAÇÃO E VALIDAÇÃO
Segundo Souza e Camargo (1996), a análise de séries temporais – e
conseqüentemente os sistemas geradores – são realizadas normalmente com quatro
objetivos possíveis:
a) Descrição: inicia-se a análise de séries temporais com a construção de um
gráfico mostrando como o fenômeno evolui no tempo. Obtém-se também medidas
descritivas simples de suas principais características.
b) Explicação: tendo-se um conjunto de observações de duas ou mais
variáveis, pode ser possível explicar o comportamento de uma em função das
demais.
c) Previsão: a previsão do comportamento futuro da variável constitui um dos
principais objetivos da análise de séries temporais.
d) Controle: quando uma série temporal mede qualidade, o objetivo da análise
pode ser de controlar o processo gerador. Uma estratégia de controle sofisticada é
ajustar um modelo à série de dados, previsões e, então, tomar, medidas corretivas
nas séries de entrada para evitar que a qualidade se afaste de um nível
estabelecido.
48
A modelagem ou identificação de sistemas utilizando modelos autoregressivos se inspira na análise de séries temporais. Todavia, diferem nas
características das variáveis envolvidas.
Na maioria das vezes a modelagem das séries temporais se utiliza de apenas
uma série de valores de apenas uma variável, o modelo que utiliza observações
anteriores desta variável são aplicados à previsão do valor futuro desta para um
dado horizonte de predição k .
Já a identificação de sistemas não se restringe a modelagem de uma série
temporal, mas sim das relações entre as séries temporais das variáveis entrada e
variáveis de saída de um dado sistema. E nem sempre o objetivo da modelagem é a
previsão.
Sobretudo os modelos auto-regressivos são aplicados com sucesso nos dois
casos o que sugere que a mesma metodologia pode ser aplicada para ambos os
casos, considerando as devidas adequações.
É necessário saber como o modelo será utilizado de forma a poder julgar se
ele incorpora, ou não, as características requeridas. Essa necessidade advém do
fato que nenhum modelo, por definição, representará o sistema real em todos os
aspectos (Aguirre, 2000).
A validação do modelo está então associada ao objetivo do mesmo e as
características do sistema modelado que se deseja reproduzir. Em outras palavras, o
modelo é tão bom quanto for sua capacidade de reproduzir a característica mais
importante do sistema para o objetivo do mesmo.
Comparar
a
simulação
do
modelo
obtido
com dados
medidos
é
provavelmente a forma mais usual de se validar um modelo. Nesse caso, deseja-se
saber se o modelo reproduz ao longo do tempo os dados observados. Em outras
palavras, deseja-se saber se o modelo serve para explicar um outro conjunto de
dados observado do mesmo sistema. É comum, portanto, referir-se à capacidade de
generalização do modelo (Aguirre, 2000).
A característica em questão é a reconstrução dinâmica da saída do sistema
pelo modelo, com o intuito de utilizar o modelo para simulação do sistema real. A
simulação de um sistema real através de modelos é conveniente para diversas áreas
de pesquisa, tais como: controle, detecção de falhas, otimização, etc.
A vantagem destes modelos em sistemas complexos é que este pode ser o
método mais rápido e prático de se obter um modelo da dinâmica do processo, A
49
desvantagem é que este modelo tem uma validade apenas local, isto é em torno de
um ponto de operação, não permitindo grandes extrapolações. (Campos e Teixeira,
2006).
Como praticamente todos os sistemas reais são não-lineares, é importante
que os dois testes sejam efetuados com o sistema operando em condições
semelhantes, principalmente se os modelos identificados forem lineares. (Aguirre,
2000).
Portanto, modelos obtidos a partir de dados coincidem com o sistema na faixa
de valores utilizadas para se realizar a estimação de seus parâmetros. E caso seja
possível realizar-se dois testes dinâmicos no sistema, um deve ser utilizado para
estimação e outro para validação do modelo. Caso isso não seja possível o mesmo
lote de dados do teste realizado será dividido para a formação dos conjuntos de
estimação e validação.
Neste trabalho ambos as formas de validação foram utilizadas. O dados
coletados do primeiro teste dinâmico realizado foram separados em conjunto de
estimação e validação.
Em um segundo momento novos testes dinâmicos foram realizados, com
diferentes estímulos de entrada, mas na mesma faixa de operação e foram utilizados
também para validação do modelo.
Cruz (2006) realiza uma identificação dinâmica de uma microturbina a gás,
utilizando modelos auto-regressivos lineares e não lineares. Para validação dos
modelos utiliza dados de vários testes dinâmicos realizados no sistema. Os índices
de desempenho são calculados para predições de um passo a frente e também de
mais passos.
Dallagnol (2005) utiliza também modelos lineares e não lineares para a
modelagem de um sistema de injeção de gás em um poço de petróleo, realizam
testes dinâmicos em um simulador em diferentes pontos de operação do processo.
Validaram o modelo com um conjunto de dados de validação também gerado pelo
simulador.
Lara (2005) valida os modelos auto-regressivos através de dados de um
conjunto de validação em malha aberta. Também realiza a validação do modelo em
malha fechada comparando sua resposta com a do sistema ambos submetidos a
perturbações na entrada. Um controlador preditivo foi utilizado no sistema de
controle em malha fechada.
50
Carvalho e Guillermo (2007) utilizaram para validação do modelo, quatro
testes dinâmicos realizados no sistema (coluna de destilação), e também 33% do
número de amostras do primeiro teste dinâmico, os outros 66% foram utilizados para
estimação dos parâmetros do modelo.
Neste trabalho os dados utilizados para compor os conjuntos de estimação e
validação foram registrados com período de amostragem de 40 segundos entre
observações resultando em uma matriz de dados de dimensão 89x5, e podem ser
vistos na Figura 26.
Valor das variáveis ao longo do tempo
100
90
80
Valores das variáveis
70
60
50
temp topo
vazao
temp base
pressao
nivel base
40
30
20
10
0
0
10
20
30
40
50
60
Tempo em segundos
70
80
90
Figura 26 - Dados adequados para a nova taxa de amostragem.
O conjunto de estimação foi montado com as primeiras 63 amostras da
matriz, e o conjunto de validação com as últimas 26 amostras. Resultando em uma
divisão de 70% dos dados para estimação e 30% para validação, uma divisão
sugerida em trabalhos científicos sobre modelos e validação.
51
4.8 – ÍNDICES DE DESEMPENHO
A validação dos modelos obtidos através da estimação de parâmetros de
modelos auto-regressivos possui um grande número de índices de desempenho
citados em trabalhos da área.
São índices que não possuem sua utilidade associada somente a modelos
auto-regressivos, mas sim a todos os tipos de modelos que visam realizar a
reconstrução dinâmica do sistema modelado.
Fazem parte desta classe não só os modelos auto-regressivos, mas também
as redes neurais artificiais, modelos de wiener, splines, e outros.
Os modelos mencionados possuem a característica de reproduzir a saída do
sistema a partir de dados de entrada sua validação através dos índices é obtida a
partir da comparação entre a saída calculada pelo modelo e a saída real observada
do sistema.
Durante a etapa de estimação dos parâmetros destes tipos de modelo, muitas
vezes, o desvio entre estas “duas saídas” é minimizado e serve de critério para
obtenção de parâmetros ótimos dos modelos.
É ideal, portanto que para a validação do modelo os índices sejam calculados
a partir da comparação entre a saída real do sistema e a obtida pelo modelo para os
dados de validação e não para os de estimação.
Os índices na sua maioria visam medir a capacidade preditiva dos modelos
para um horizonte de k passos adiante, onde k  1, ,  .
Uma forma de fazer predição é, no instante t  1 , montar um vetor de
regressores  com observações obtidas do conjunto de dados e simplesmente
calcular yˆ (t )   T (t  1)ˆ . O valor predito para o instante seguinte, ou seja, k , é yˆ (t ) .
Por esta razão, tal predição é denominada de predição de um passo à frente, sendo
que  T (t  1)ˆ é o preditor de um passo à frente. Uma outra forma de simular um
modelo é reutilizar predições passadas para compor o vetor de regressores a fim de
continuar fazendo predição. O sinal simulado desta maneira é chamado de
simulação livre ou predição de infinitos passos à frente (Aguirre, 2000).
A validação do modelo pode ser então determinada a partir da seqüência de
valores calculados com determinado índice de desempenho à medida que se faz
predições de mais passos à frente.
52
Neste trabalho a validade do modelo sob este critério será verificada, depois
de definida a ordem do modelo a partir de dois métodos de seleção: AIC e validação
cruzada. Os modelos resultantes serão submetidos à validação com os demais
testes dinâmicos e seu desempenho calculado variando o horizonte de predição.
Borjas e Garcia (2004) apresentam a modelagem de uma unidade de CCF
(Craquamento Catalítico em leito Fluidizado), sistema multivariável e não linear,
através de identificação por sub-espaços e o método de predição de erro. Os dados
utilizados para a identificação foram obtidos através da excitação com sinais pseudoaleatórios multinível, o desempenho dos modelos encontrados são verificados
através de dois índices: média relativa do erro quadrático (MRSE) e média da
variância relativa (MVAF).
N
n
MRSE(%) 
1

n k 1
 y  yˆ 
2
t 1
N
 100
 y 
2
t 1
MVAF (%) 
1 n  var iânciay  yˆ 
  100
  1 
n k 1 
var iânciay  
Os índices calculados são o R2, MAE (mean absolute error), MSE (mean
square error). Todos são índices de desempenho de modelos citados em literatura.
O R2 é um índice que mede o quanto o modelo encontrado se ajusta aos dados
observados do sistema varia entre 0 e 1. O MAE retorna a média dos valores
absolutos dos erros encontrados quando se compara saída do modelo e saída do
sistema. O MSE é semelhante ao MAE, porém é mais preciso, pois não mascara a
variância, ao utilizar o quadrado dos erros, problema ocorrido no MAE por utilizar a
função modular.
N
R2  1
 y(t)  yˆ (t )
2
t 1
N
 y(t)  y(t )
2
t 1
N
 y (t )  yˆ (t )
MAE 
t 1
N
53
N
2
y(t)  yˆ(t)
MSE
t 1
N
Sendo y(t ) os valores observados da saída do sistema, yˆ (t ) os valores
calculados na saída do modelo e N o número de amostras utilizadas (Carvalho e
Guillermo, 2007).
Trautwein (2004) utiliza redes neurais artificiais do tipo multicamadas para
previsão do consumo de água em curtíssimo prazo. Para verificação do
desempenho dos modelos obtidos são usados os seguintes índices: MSE, MAE,
MAPE, MRSE, VAR e R2.
Zanata (2005) desenvolve um sensor virtual para medição por inferência da
composição do produto de uma coluna de destilação, utiliza redes neurais artificiais,
e para medir o ajuste utiliza os índices: MSE, MAE e erro máximo.
Uma forma de validar a capacidade preditiva dos modelos preditores é utilizar
a relação que compara o desempenho do preditor de k-passos à frente com o
desempenho do preditor de 1 passo à frente, i.e., calculando os erros de predição
dos diferentes modelos com os dados medidos até o instante t :
N
2
 yˆ t / t  k  yt 
k
Rk 
t 1
N
 yˆ t / t  1 yt 
2
1
t 1
Note que quando Rk  1 o erro de predição k passos à frente será menor que
o erro de predição um passo à frente, implicando que o modelo preditor possui
melhor precisão e que este pode ser considerado como um preditor candidato (Lara,
2005).
4.9 – ANÁLISE DE RESÍDUOS
Uma outra forma de validação de modelos bastante utilizada é a análise dos
resíduos do modelo sob critérios estatísticos.
Os testes de simulação não indicam se o modelo possui falhas corrigíveis. A
análise de resíduos designa um conjunto de testes que são efetuados para verificar
se os resíduos são aleatórios ou não (Aguirre, 2000).
54
Logo se faz necessário realizar a análise de resíduos nos modelos obtidos por
este trabalho, com o intuito de encontrar um modelo que seja válido sob dois
métodos: índices de desempenho e análise dos resíduos, e que ainda assim possua
baixa complexidade, ou seja, número reduzido de parâmetros. Esta última
característica é encontrada nos textos como parcimônia.
Ao se testar um vetor de resíduos for verificado que se trata de uma variável
aleatória branca, isso significa que não há informação útil nos resíduos, ou seja, o
modelo explicou tudo que era possível explicar. Infelizmente isso não quer dizer que
a simulação livre do modelo será boa, mas simplesmente que a simulação de um
passo a frente do modelo será boa. Um vetor de valores e(t ) é branco se a sua
função de autocorrelação for nula para todos os valores de atraso maiores ou iguais
a um (Aguirre, 2000).
Portanto, a metodologia utilizada para a validação dos modelos encontrados
caminhará no sentido de realizar testes que venham verificar se os resíduos da
modelagem se comportam como ruído branco.
Mesmo que a autocorrelação dos resíduos seja nula para todo t  0 , este
resultado é obtido a partir de um conjunto de dados específico, e o que se deseja do
modelo é que seu desempenho seja o mesmo para outras situações. De forma a
verificar o quão geral é o modelo procede-se como a utilização da correlação
cruzada entre a entrada e o vetor de resíduos do modelo para diferentes conjuntos
de dados.
A implicação deste resultado é que as predições de um passo à frente do
modelo terão características semelhantes se calculadas para um outro conjunto de
dados (Aguirre, 2000).
Dallagnol (2005) utiliza a análise de resíduos através da autocorrelação, de
maneira intercalada com a estimação de parâmetros dos modelos. Desta forma ele
define a ordem da componente MA dos modelos ARMAX e NARMAX utilizados para
a modelagem do sistema.
Os resultados da validação pela análise dos resíduos são apresentados no
capítulo V.
55
CAPÍTULO V
RESULTADOS
Os resultados deste trabalho serão apresentados de forma numérica e gráfica
e seguirão uma seqüência de apresentação semelhante à seqüência das etapas da
metodologia aplicada. Em cada etapa da aplicação da metodologia, resultados
intermediários são obtidos, até que se alcance os resultados finais que estão
associados à etapa de validação dos modelos obtidos.
5.1- DADOS UTILIZADOS
Já foi mencionado no capítulo de metodologia que os dados utilizados para
formar o conjunto de estimação e o primeiro conjunto de validação foram originados
do primeiro teste dinâmico realizado na coluna de destilação. Os dados podem ser
vistos na sua forma original (período de amostragem de 1 segundo) na Figura 15 e
adequados ao novo período de amostragem na Figura 26.
Para a divisão dos dados em conjuntos de estimação e validação, foram
selecionadas as primeiras 63 amostras da matriz de dados adequada para compor o
conjunto de estimação e as 26 últimas para o conjunto de validação. Os valores
percentuais dos dois conjuntos são respectivamente 70% e 30%. Nas Figuras 27 e
28 são apresentados os dados referentes aos conjuntos de estimação e validação
respectivamente.
56
Valores das variáveis ao longo do tempo
100
90
80
Valores das variáveis
70
60
50
temp topo
vazao
temp base
pressao
nivel base
40
30
20
10
0
0
10
20
30
40
Tempo em segundos
50
60
70
Figura 27 - Conjunto de estimação.
Valores das variáveis ao longo do tempo
90
80
Valores das variáveis
70
60
50
40
temp topo
vazao
temp base
pressao
nivel base
30
20
10
0
0
5
Figura 28 - Primeiro conjunto de validação.
10
15
20
Tempos em segundos
25
30
57
Apesar de apresentadas nas figuras 27 e 28 todas as variáveis disponíveis
para a modelagem, através da correlação cruzada apenas a temperatura de topo e a
vazão foram selecionadas. Isto é válido para os demais conjuntos de dados
utilizados para validação do modelo.
Como praticamente todos os sistemas reais são não-lineares, é importante
que os dois testes sejam efetuados com o sistema operando em condições
semelhantes, principalmente se os modelos identificados forem lineares. (AGUIRRE,
2000).
Os demais conjuntos serão utilizados com a amostragem original para a
validação dos modelos obtidos, e originam-se de outros quatro testes dinâmicos
realizados na coluna de destilação, sob condições operacionais semelhantes às do
primeiro teste. Nas Figuras 29, 30, 31 e 32 são apresentados os dados do segundo,
terceiro, quarto e quinto testes dinâmicos.
Valores das variáveis ao longo do tempo
100
90
80
Valores das variáveis
70
60
50
temp topo
vazao
temp base
pressao
nivel topo
40
30
20
10
0
0
500
1000
1500
Tempo em segundos
Figura 29 - Dados do segundo teste dinâmico utilizado para validação.
2000
2500
58
Vaores das variáveis ao longo do tempo
100
90
80
Valores das variáveis
70
60
50
temp topo
vazao
temp base
pressao
nivel base
40
30
20
10
0
0
200
400
600
800
Tempo em segundos
1000
1200
1400
Figura 30 - Dados do terceiro teste dinâmico utilizado para validação.
Valores das variáveis ao longo do tempo
120
Valores das variáveis
100
80
60
temp topo
vazao
temp base
pressao
nivel base
40
20
0
0
200
400
600
800 1000 1200
Tempo em segundos
Figura 31 - Dados do quarto teste dinâmico utilizado para validação.
1400
1600
1800
2000
59
Valores das variáveis ao longo do tempo
100
90
80
Valores das variáveis
70
60
temp topo
vazao
temp base
pressao
nivel base
50
40
30
20
10
0
0
200
400
600
800
Tempo em segundos
1000
1200
1400
Figura 32 - Dados do quinto teste dinâmico utilizado para validação.
5.2- MÉTODOS DE SELEÇÃO DA ORDEM DO MODELO
Os métodos utilizados para seleção da ordem do modelo, foram o AIC, e
validação cruzada (duas formas).
Em um primeiro instante foi elaborada uma rotina no MATLAB® que testasse
todas as possibilidades de combinações dos parâmetros n a , nb e n k para estimação
dos coeficientes an e bn de um modelo ARX – vide capítulo 3 seção 3.1. Os
parâmetros mencionados representam respectivamente o número de atrasos na
saída, número de atrasos na entrada e número de atrasos da saída para a entrada,
que serão utilizados como sinais de entrada do modelo.
A rotina acima citada realiza um processo iterativo de estimação dos
coeficientes do modelo ARX tendo como argumento de entrada os valores de n a , nb
e nk e também os dados do conjunto de estimação. Os valores dos parâmetros são
elevados gradualmente, na seguinte seqüência: primeiro na , depois nb e por último
nk , com isso todas as possibilidades de combinações destes parâmetros foram
60
percorridas e para comparar a eficiência dos modelos gerados foi utilizado o MSE
sobre o erro de predição do modelo quando comparado ao conjunto de estimação e
também o primeiro conjunto de validação ambos com amostragem de 40 segundos.
Este método de verificação do MSE tanto para o conjunto de estimação quanto para
o de validação é conhecido como validação cruzada e é utilizado para verificar em
que ponto o modelo deixa de ser geral para se tornar específico aos dados de
estimação. Os valores selecionados para as ordens estavam compreendidos nos
seguintes intervalos: na  1, , 16 , nb  0,  , 15 e nk  0,  , 15 .
O segundo método, AIC, foi utilizado testando apenas modelos com ordem
n k  0 . E as demais ordens foram elevadas da mesma forma que no método anterior
da validação cruzada. A comparação entre os modelos gerados foi feita através do
AIC calculado para cada um deles.
Segundo Aguirre (2000) o critério de informação de Akaike possui um valor
mínimo que quando alcançado, indica que o aumento da ordem do modelo não mais
implica em melhoria significativa do ajuste.
Carvalho e Guillermo (2007) utilizaram o AIC seguindo este mesmo
procedimento e obtiveram para os mesmos conjuntos de dados aqui apresentados,
o gráfico apresentado na Figura 33:
Valores de AIC variando a ordem dos modelos
0
-0.5
Valores de AIC calculados
-1
-1.5
-2
-2.5
-3
-3.5
-4
1
1.5
2
2.5
3
3.5
Ordem dos modelos
4
4.5
5
Figura 33 - Valores do AIC calculado através das iterações da rotina. (Carvalho e Guillermo, 2007).
61
Percebe-se que o ponto no qual o AIC passa por um mínimo é onde as
ordens na e nb eram iguais a quatro. A rotina utilizada elevava ambas as ordens ao
mesmo tempo.
Já a curva do AIC calculada neste trabalho, é apresentada na Figura 34.
Valores de AIC ao longo das iterações
1
0
Valores de AIC calculados
-1
-2
-3
-4
-5
-6
-7
-8
-9
0
5
10
15
20
25
Iterações da rotina
30
35
40
Figura 34 - Curva do AIC calculado através das iterações da rotina atual.
É possível notar que diferente do gráfico da figura 33, o AIC atinge o mínimo
na última iteração de sua rotina, e sugere que seu decréscimo é possível com o
aumento da ordem do modelo que neste ponto é na  20 e nb  20 . A rotina atual
aumenta as ordens do modelo de forma separada, conforme mencionado.
Desta forma é percebido que o AIC sugere um modelo sobre-parametrizado
se comparado ao modelo sugerido no trabalho de Carvalho e Guillermo (2007).
62
Tabela 3 – Valores calculados do AIC para modelos ARX de diferentes ordens.
Ordem na
Ordem nb
1
1
2
2
3
3
4
4
5
5
6
6
7
7
8
8
9
9
10
10
11
11
12
12
13
13
14
14
15
15
16
16
17
17
18
18
19
19
20
20
0
1
1
2
2
3
3
4
4
5
5
6
6
7
7
8
8
9
9
10
10
11
11
12
12
13
13
14
14
15
15
16
16
17
17
18
18
19
19
20
AIC
0.418
0.4375
0.4684
-2.2217
-3.7464
-3.8239
-3.9802
-4.3563
-4.3692
-4.3582
-4.5619
-4.5455
-4.5337
-4.5432
-4.5659
-4.5497
-4.524
-4.4957
-4.6346
-4.6616
-4.646
-4.6265
-4.6706
-4.6476
-4.9032
-4.8753
-4.9984
-4.9687
-4.9811
-4.9761
-5.2085
-5.2723
-5.3663
-5.4002
-5.4755
-5.9752
-6.7593
-6.7823
-7.5572
-8.0313
No primeiro método que utiliza a validação cruzada o índice de desempenho
para comparação dos modelos foi o MSE, calculado para o erro de predição de um
passo a frente, sobre o conjunto de estimação e sobre o primeiro conjunto de
validação.
63
As variações testadas da ordem do modelo incluíam a variação de nk ,
portanto o número de combinações possíveis foi bem maior que no método que
utilizou o AIC, mesmo que o limite superior para as ordens tenha sido menor. Os
resultados conseguidos por este método sugerem um modelo também sobreparametrizado, em especial o modelo ARX com as seguintes ordens: n a  16 ,
nb  15 e nk  15 , o que possui ordem mais elevada entre os testados, foi o modelo
que obteve menor MSE tanto para o conjunto de estimação quanto validação.
Algumas diferenças existem em relação aos resultados obtidos no trabalho de
Carvalho e Guillermo (2007) e o presente trabalho, logo, serão relacionadas:
 O número de amostras presentes nos conjuntos de estimação e validação
neste trabalho são inferiores devido à escolha de um novo período de amostragem.
 A rotina utilizada neste trabalho eleva as ordens do modelo de forma
separada, conforme sugerido em Aguirre (2000).
 O número de variáveis de entrada utilizadas na modelagem atual é menor.
Apenas a variável vazão foi selecionada como entrada do modelo, devido às
funções de correlação cruzada, calculadas durante a aplicação da metodologia,
enquanto que no trabalho anterior foram utilizadas todas as variáveis medidas como
entradas do sistema, com exceção da temperatura de topo que é saída do modelo
nos dois casos.
As diferenças pesam, sobretudo no modelo sugerido pelo AIC obtido que
neste trabalho encontra-se sobre-parametrizado.
Os modelos obtidos pelos dois métodos apresentaram um bom ajuste para o
conjunto de estimação e o primeiro conjunto de validação considerando o MSE e os
demais índices de desempenho, porém para os outros quatro conjuntos de validação
o desempenho não foi o mesmo.
5.3- RESULTADOS GRÁFICOS E NUMÉRICOS DOS MODELOS OBTIDOS
A seguir serão apresentados os resultados gráficos obtidos com o modelo
resultante do método que utiliza o AIC. Trata-se de um modelo ARX com as ordens
n a  20 , nb  20 e n k  0 .
64
Os resultados do modelo ARX [20 20 0] para o primeiro, segundo, terceiro,
quarto e quinto conjuntos de validação são apresentados nas Figuras 35, 36, 37, 38
e 39 respectivamente.
Temperatura de topo ao longo do tempo
53
sistema
modelo
Valores da temperatura de topo
52
51
50
49
48
47
46
0
5
10
15
20
Tempo em segundos
25
30
Figura 35 - Saída do modelo comparada com a saída do sistema, para o primeiro conjunto de
validação.
65
Temperatura de topo ao longo do tempo
80
sistema
modelo
Valores da temperatura de topo
75
70
65
60
55
50
45
0
500
1000
1500
Tempo em segundos
2000
2500
Figura 36 - Saída do modelo comparada com a saída do sistema, para o segundo conjunto de
validação.
Temperatura de topo ao longo do tempo
80
sistema
modelo
Valores da temperatura de topo
78
76
74
72
70
68
66
64
62
0
200
400
600
800
Tempo em segundos
1000
1200
1400
Figura 37 - Saída do modelo comparada com a saída do sistema, para o terceiro conjunto de
validação.
66
Temperatura de topo ao longo do tempo
81
Valores da temperatura de topo
80
79
sistema
modelo
78
77
76
75
74
73
72
Figura
0
200
400
600
800 1000 1200
Tempo em segundos
1400
1600
1800
2000
38 - Saída do modelo comparada com a saída do sistema, para o quarto conjunto de
validação.
Temperatura de topo ao longo do tempo
80
sistema
modelo
Valores da temperatura de topo
79
78
77
76
75
74
73
0
Figura
200
400
600
800
Tempo em segundos
1000
1200
1400
39 - Saída do modelo comparada com a saída do sistema, para o quinto conjunto de
validação.
67
Tabela 4 - Índices calculados para os cinco conjuntos de validação referente ao modelo obtido com o
AIC.
2
Validação
MSE
MSE(%)
MAE
MAE(%)
MSE(R)
ERRO(%)
AMPVAR
DISCREP
R
1
0.0456
4.5635
0.0987
9.8657
0.0411
0.1997
0.1316
0.0022
0.9841
2
8.8615
886.1522
2.8710
287.0999
7.9754
-4.5629
0.1549
0.0240
0.8058
3
18.6094
1860.9
4.2364
423.6439
16.7485
-5.7198
0.2051
0.0299
-0.241
4
33.6543
3365.4
5.7669
576.6855
30.2888
-7.2510
1.3196
0.0378
-96.94
5
29.5411
2954.1
5.3887
538.8742
26.5870
-6.7709
46.6130
0.0353
-62.18
A Tabela 4 resume os resultados numéricos do modelo obtido através do
método que utiliza o AIC. Os valores dos índices foram calculados sobre o erro de
predição de um passo a frente do modelo para os cinco conjuntos de validação.
Os índices de desempenho anteriormente citados constam na tabela 4 bem
como outros índices: MSE percentual, MAE percentual, MSE regularizado ou
normalizado, erro percentual, amplitude da variância, coeficiente de discrepância.
erro % 
yˆ t  y t 
y t 
 100
Onde ŷ t  é a média dos dados da saída do modelo e yt  é a média dos
dados da saída observados no sistema.
Amp var 
ˆ  

Onde ˆ é a variância dos dados da saída do modelo e  é a variância dos
dados da saída observados no sistema.
N
2
 yˆ (t )  y(t )
Discrep 
t 1
N
 yˆ (t) 
t 1
2
N

 y(t )
2
t 1
Desta forma, os demais índices de desempenho são apenas variações dos
anteriormente explicados como os índices MSE e MAE percentuais e o MSE(R) ou
MSE normalizado apresentado na seção 4.8 como MRSE.
A partir de agora os resultados gráficos apresentados são referentes ao
modelo obtido com a primeira rotina de validação cruzada, que utilizava apenas o
conjunto de estimação e o primeiro conjunto de validação. O modelo obtido foi o
ARX com as ordens na  16 , nb  15 e nk  15 .
68
Os resultados do modelo ARX [16 15 15] para o primeiro, segundo, terceiro,
quarto e quinto conjuntos de validação são apresentados nas Figuras 40, 41, 42, 43
e 44 respectivamente.
Temperatura de topo ao longo do tempo
53
sistema
modelo
Valores da temperatura de topo
52
51
50
49
48
47
46
0
5
10
15
20
Tempo em segundos
25
30
Figura 40 - Saída do modelo comparada com a saída do sistema, para o primeiro conjunto de
validação.
69
Temperatura de topo ao longo do tempo
85
sistema
modelo
Valores da temperatura de topo
80
75
70
65
60
55
50
0
500
1000
1500
Tempo em segundos
2000
2500
Figura 41 - Saída do modelo comparada com a saída do sistema, para o segundo conjunto de
validação.
Temperatura de topo ao longo do tempo
80
sistema
modelo
Valores da temperatura de topo
78
76
74
72
70
68
66
0
200
400
600
800
Tempo em segundos
1000
1200
1400
Figura 42 - Saída do modelo comparada com a saída do sistema, para o terceiro conjunto de
validação.
70
Temperatura de topo ao longo do tempo
80.5
Valores da temperatura de topo
80
79.5
sistema
modelo
79
78.5
78
77.5
77
76.5
76
Figura
0
200
400
600
800 1000 1200
Tempo em segundos
1400
1600
1800
2000
43 - Saída do modelo comparada com a saída do sistema, para o quarto conjunto de
validação.
Temperatura de topo ao longo do tempo
80
Valores da temperatura de topo
79.5
sistema
modelo
79
78.5
78
77.5
77
76.5
0
Figura
200
400
600
800
Tempo em segundos
1000
1200
1400
44 - Saída do modelo comparada com a saída do sistema, para o quinto conjunto de
validação.
71
Na Tabela 5, são apresentados os resultados do modelo ARX [16 15 15] para
diversos índices de desempenho, considerando os cinco conjuntos de validação.
Tabela 5 - Índices calculados para os cinco conjuntos de validação referente ao primeiro modelo
obtido com a validação cruzada.
2
Validação
MSE
MSE(%)
MAE
MAE(%)
MSE(R)
ERRO(%)
AMPVAR
DISCREP
R
1
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
1.0000
2
0.8240
82.3957
0.8733
87.0
0.7416
-1.3720
0.0245
0.0072
0.9844
3
1.2488
124.8753
1.0936
109.0
1.1239
-1.4755
0.0301
0.0076
0.9317
4
1.7455
174.5510
1.3096
131.0
1.5710
-1.6467
0.0736
0.0084
-11.72
5
1.6535
165.3549
1.2701
13096.0
1.4882
-1.5958
3.4752
0.0081
-36.63
Visando a obtenção de um modelo mais parcimonioso – relação número de
parâmetros/ajuste razoável – e com melhores resultados para todos os conjuntos de
validação, foi elaborada uma nova rotina para a validação cruzada. Desta vez, o
critério para escolha do modelo foi a manutenção do MSE dentro de uma dada
tolerância para os cinco conjuntos de validação.
E por último são apresentados os resultados gráficos do segundo modelo
obtido pelo método que utiliza a validação cruzada, selecionado pelo MSE calculado
para os cinco conjuntos de validação. O modelo sugerido foi o ARX com as ordens
n a  6 , nb  12 e n k  12 .
Os resultados do modelo ARX [6 12 12] para o primeiro, segundo, terceiro,
quarto e quinto conjuntos de validação são apresentados nas Figuras 45, 46, 47, 48
e 49 respectivamente.
72
Temperatura de topo ao longo do tempo
53
sistema
modelo
Valores da temperatura de topo
52
51
50
49
48
47
46
0
5
10
15
20
Tempo em segundos
25
30
Figura 45 - Saída do modelo comparada com a saída do sistema, para o primeiro conjunto de
validação.
Temperatura de topo ao longo do tempo
80
sistema
modelo
Valores da temperatura de topo
75
70
65
60
55
50
0
500
1000
1500
Tempo em segundos
2000
2500
Figura 46 - Saída do modelo comparada com a saída do sistema, para o segundo conjunto de
validação.
73
Temperatura de topo ao longo do tempo
80
sistema
modelo
Valores da temperatura de topo
78
76
74
72
70
68
66
0
200
400
600
800
Tempo em segundos
1000
1200
1400
Figura 47 - Saída do modelo comparada com a saída do sistema, para o terceiro conjunto de
validação.
Temperatura de topo ao longo do tempo
80.4
Valores da temperatura de topo
80.2
80
79.8
79.6
sistema
modelo
79.4
79.2
79
78.8
78.6
Figura
0
200
400
600
800 1000 1200
Tempo em segundos
1400
1600
1800
2000
48 - Saída do modelo comparada com a saída do sistema, para o quarto conjunto de
validação.
74
Temperatura de topo ao longo do tempo
80.1
sistema
modelo
80
Valores da temperatura de topo
79.9
79.8
79.7
79.6
79.5
79.4
79.3
79.2
79.1
Figura
0
200
400
600
800
Tempo em segundos
1000
1200
1400
49 - Saída do modelo comparada com a saída do sistema, para o quinto conjunto de
validação.
Na Tabela 6, são apresentados os resultados do modelo ARX [6 12 12] para
diversos índices de desempenho, considerando os cinco conjuntos de validação.
Tabela 6 - Índices calculados para os cinco conjuntos de validação referente ao segundo modelo
obtido com a validação cruzada.
2
Validação
MSE
MSE(%)
MAE
MAE(%)
MSE(R)
ERRO(%)
AMPVAR
DISCREP
R
1
0.0002
0.0242
0.0042
0.4
0.0002
0.0085
0.0065
0.0002
0.9999
2
0.0252
2.5227
0.0418
4.2
0.0227
0.0143
0.0073
0.0013
0.9995
3
0.0130
1.2956
0.0799
8.0
0.0117
0.0944
0.0120
0.0008
0.9993
4
0.0243
2.4275
0.1540
15.4
0.0218
0.1933
0.0701
0.0010
0.8469
5
0.0171
1.7096
0.1285
540.4
0.0154
0.1610
0.3458
0.0008
-0.293
5.4- COEFICIENTES ESTIMADOS DOS MODELOS OBTIDOS
Os três modelos ARX comparados através de seus resultados gráficos e
numéricos são apresentados na sua forma compacta e os coeficientes dos
polinômios estimados através do método dos mínimos quadrados. Nas Tabelas 7, 8
e 9 são sumariados os valores dos coeficientes dos três modelos ARX obtidos.
75
A(q 1 ) y (t )  B(q 1 )u (t )  e(t )
20
A(q 1 )  1   a n q  n
n 1
19
B1(q 1 )  b1   bn q n 1
n 2
Tabela 7 - Coeficientes estimados para o modelo ARX [20 20 0].
 
 
A q 1
B q 1
a1
- 0.6562
b1
-0.1057
a2
0.06746
b2
- 1.197
a3
- 0.5522
b3
2.066
a4
0.07402
b4
- 1.302
a5
0.3253
b5
1.233
a6
- 0.2144
b6
- 0.9118
a7
0.3787
b7
1.233
a8
- 0.4303
b8
- 2.201
a9
0.1074
b9
0.8053
a10
0.1076
b10
0.8458
a11
0.1171
b11
- 0.1641
a12
- 0.3068
b12
- 0.1804
a13
0.6003
b13
0.4141
a14
- 0.08385
b14
- 0.7116
a15
- 0.6896
b15
- 0.06598
a16
0.303
b16
1.031
a17
0.1898
b17
- 0.2554
a18
- 0.3674
b18
- 0.4943
a19
0.1941
b19
0.4668
a 20
- 0.082
b20
- 0.175
76
A(q 1 ) y (t )  B (q 1 )u (t  15)  e(t )
16
A(q 1 )  1   a n q  n
n 1
15
B1(q 1 )   bn q  n1
n 1
Tabela 8 - Coeficientes estimados para o modelo ARX [16 15 15].
 
 
A q 1
B q 1
a1
- 1.173
b1
1.656
a2
0.9362
b2
- 2.148
a3
- 0.6055
b3
1.436
a4
- 0.5531
b4
- 0.4922
a5
0.7945
b5
- 0.4431
a6
- 0.5493
b6
0.005149
a7
0.4292
b7
- 0.04209
a8
- 0.4666
b8
0.03544
a9
0.5344
b9
- 0.01602
a10
- 0.3969
b10
- 0.01787
a11
0.09228
b11
0.06666
a12
0.09135
b12
- 0.0748
a13
0.0936
b13
0.06738
a14
- 0.3927
b14
- 0.03392
a15
0.4539
b15
0.03034
a16
- 0.2713
77
A(q 1 ) y (t )  B (q 1 )u (t  12)  e(t )
6
A(q 1 )  1   a n q  n
n 1
12
B1(q 1 )   bn q  n1
n 1
Tabela 9- Coeficientes estimados para o modelo ARX [6 12 12].
 
 
A q 1
B q 1
a1
- 1.021
b1
0.4096
a2
0.003319
b2
- 0.4365
a3
- 0.2333
b3
- 0.01006
a4
0.03998
b4
0.03309
a5
0.133
b5
- 0.004944
a6
0.07535
b6
- 0.02207
b7
0.003886
b8
- 0.03118
b9
0.009069
b10
- 0.02352
b11
0.03605
b12
0.01308
5.5- RESULTADOS DOS MODELOS VARIANDO O HORIZONTE DE PREDIÇÃO
Os modelos obtidos pelos métodos expostos, serão validados quanto ao sua
capacidade de predição de k passos à frente, para os cinco conjuntos de validação.
O MSE será calculado para k  1,  , 20 . O resumo dos valores de MSE calculados
são apresentados na Tabela 10.
78
Tabela 10 - MSE calculado para o erro de predição de k passos à frente do modelo ARX [6 12 12]
para os cinco conjuntos de validação.
k passos
1º conjunto
2º conjunto
3º conjunto
4º conjunto
5º conjunto
1
0.0002
0.0252
0.0130
0.0243
0.0171
2
0.0006
0.0539
0.0391
0.0988
0.0694
3
0.0010
0.0758
0.0773
0.2262
0.1589
4
0.0010
0.0973
0.1386
0.4568
0.3207
5
0.0010
0.1361
0.2318
0.8289
0.5817
6
0.0010
0.1876
0.3606
1.3452
0.9437
7
0.0010
0.2438
0.5279
2.0194
1.4163
8
0.0010
0.2920
0.7447
2.8905
2.0268
9
0.0010
0.3385
1.0154
3.9768
2.7878
10
0.0010
0.4157
1.3358
5.2811
3.7012
11
0.0010
0.5201
1.7097
6.8179
4.7770
12
0.0010
0.6367
2.1429
8.6056
6.0280
13
0.0010
0.7431
2.6410
10.6533
7.4606
14
0.0010
0.8700
3.2045
12.9665
9.0783
15
0.0010
1.0236
3.8367
15.5543
10.8878
16
0.0010
1.1874
4.5375
18.4258
12.8951
17
0.0010
1.3540
5.3043
21.5873
15.1042
18
0.0010
1.5233
6.1385
25.0440
17.5188
19
0.0010
1.7046
7.0448
28.8021
20.1431
20
0.0010
1.9154
8.0267
32.8676
22.9811
Nas Figuras 50, 52, 54, 56 e 58 são apresentados os gráficos do MSE para
um horizonte de predição variando de 0 a 20 passos à frente, para os cinco
conjuntos de validação.
Nas Figuras 51, 53, 55, 57 e 59 são apresentados os gráficos de comparação
da saída do modelo com a saída do sistema para um horizonte de predição de 20
passos à frente, para os cinco conjuntos de validação.
79
-4
12
MSE calculado variando o horizonte de predição
x 10
MSE
11
10
MSE calculado
9
8
7
6
5
4
3
2
0
2
4
6
8
10
12
Horizonte de predição
14
16
18
20
Figura 50 - MSE calculado para o erro de predição de k passos à frente para o primeiro conjunto de
validação.
Temperatura de topo ao longo do tempo
53
sistema
modelo
Valores da temperatura de topo
52
51
50
49
48
47
46
0
5
10
15
20
Tempo em segundos
25
30
Figura 51 - Comparação entre saída do modelo e saída do sistema para a predição de 20 passos à
frente para o primeiro conjunto de validação.
80
MSE calculado variando o horizonte de predição
2
MSE
1.8
1.6
MSE calculado
1.4
1.2
1
0.8
0.6
0.4
0.2
0
0
2
4
6
8
10
12
Horizonte de predição
14
16
18
20
Figura 52 - MSE calculado para o erro de predição de k passos à frente para o segundo conjunto de
validação.
Temperatura de topo ao longo do tempo
85
sistema
modelo
Valores da temperatura de topo
80
75
70
65
60
55
50
0
500
1000
1500
Tempo em segundos
2000
2500
Figura 53 - Comparação entre saída do modelo e saída do sistema para a predição de 20 passos à
frente para o segundo conjunto de validação.
81
MSE calculado variando o horizonte de predição
9
MSE
8
7
MSE calculado
6
5
4
3
2
1
0
0
2
4
6
8
10
12
Horizonte de predição
14
16
18
20
Figura 54 - MSE calculado para o erro de predição de k passos à frente para o terceiro conjunto de
validação.
Temperatura de topo ao longo do tempo
84
sistema
modelo
Valores da temperatura de topo
82
80
78
76
74
72
70
68
66
0
200
400
600
800
Tempo em segundos
1000
1200
1400
Figura 55 - Comparação entre saída do modelo e saída do sistema para a predição de 20 passos à
frente para o terceiro conjunto de validação.
82
MSE calculado variando o horizonte de predição
35
MSE
30
MSE calculado
25
20
15
10
5
0
0
2
4
6
8
10
12
Horizonte de predição
14
16
18
20
Figura 56 - MSE calculado para o erro de predição de k passos à frente para o quarto conjunto de
validação.
Temperatura de topo ao longo do tempo
87
sistema
modelo
Valores de temperatura de topo
86
85
84
83
82
81
80
79
78
0
200
400
600
800 1000 1200
Tempo em segundos
1400
1600
1800
2000
Figura 57 - Comparação entre saída do modelo e saída do sistema para a predição de 20 passos à
frente para o quarto conjunto de validação.
83
MSE calculado variando o horizonte de predição
25
MSE
MSE calculado
20
15
10
5
0
0
2
4
6
8
10
12
Horizonte de predição
14
16
18
20
Figura 58 - MSE calculado para o erro de predição de k passos à frente para o quinto conjunto de
validação.
Temperatura de topo ao longo do tempo
86
Valores da temperatura de topo
85
sistema
modelo
84
83
82
81
80
79
0
200
400
600
800
Tempo em segundos
1000
1200
1400
Figura 59 - Comparação entre saída do modelo e saída do sistema para a predição de 20 passos à
frente para o quinto conjunto de validação.
84
5.6- RESULTADOS DA ANÁLISE DOS RESÍDUOS
Os modelos serão agora submetidos aos testes de análise de resíduos para
verificar se existe correlação significativa – duas vezes o desvio padrão – nos
resíduos e também entre os resíduos e a entrada. Serão comparados dois dos
modelos obtidos utilizando o conjunto de estimação e os cinco conjuntos de
validação.
Nas Figuras 60, 62, 64, 66, 68 e 70 são apresentados os gráficos da função
de autocorrelação dos resíduos do modelo ARX [20 20 0], para o conjunto de
estimação e para os cinco conjuntos de validação.
Já nas Figuras 61, 63, 65, 67, 69 e 71 são apresentados os gráficos da
função correlação cruzada entre a entrada e os resíduos do modelo ARX [20 20 0],
para o conjunto de estimação e para os cinco conjuntos de validação.
Autocorrelação dos resíduos do modelo
Valor de autocorrelação calculado
1
0.5
0
-0.5
0
2
4
6
8
10
Atrasos
12
14
Figura 60 - Autocorrelação dos resíduos para o conjunto de estimação.
16
18
20
85
Correlação cruzada entre entrada e resíduos
0.3
Valor de correlação cruzada calculado
0.2
0.1
0
-0.1
-0.2
-0.3
-0.4
-20
-15
-10
-5
0
Atrasos
5
10
15
20
Figura 61 - Correlação cruzada entre os resíduos e a entrada para o conjunto de estimação.
Autocorrelação dos resíduos do modelo
1
Valor de autocorrelação calculado
0.8
0.6
0.4
0.2
0
-0.2
-0.4
0
2
4
6
8
10
Atrasos
12
14
16
Figura 62 - Autocorrelação dos resíduos para o primeiro conjunto de validação.
18
20
86
Correlação cruzada entre entrada e resíduos
1
Valor de correlação cruzada calculado
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-20
-15
-10
-5
0
Atrasos
5
10
15
20
Figura 63 - Correlação cruzada entre os resíduos e a entrada para o primeiro conjunto de validação.
Autocorrelação dos resíduos do modelo
Valor de autocorrelação calculado
0.8
0.6
0.4
0.2
0
-0.2
0
2
4
6
8
10
Atrasos
12
14
16
Figura 64 - Autocorrelação dos resíduos para o segundo conjunto de validação.
18
20
87
Correlação cruzada entre entrada e resíduos
0.1
Valor de correlação cruzada calculado
0
-0.1
-0.2
-0.3
-0.4
-0.5
-0.6
-20
-15
-10
-5
0
Atrasos
5
10
15
20
Figura 65 - Correlação cruzada entre os resíduos e a entrada para o segundo conjunto de validação.
Autocorrelação dos resíduos do modelo
Valor de autocorrelação calculado
0.8
0.6
0.4
0.2
0
-0.2
0
2
4
6
8
10
Atrasos
12
14
16
Figura 66 - Autocorrelação dos resíduos para o terceiro conjunto de validação.
18
20
88
Correlação cruzada entre entrada e resíduos
0.1
Valor de correlação cruzada calculado
0
-0.1
-0.2
-0.3
-0.4
-0.5
-0.6
-0.7
-0.8
-20
-15
-10
-5
0
Atrasos
5
10
15
20
Figura 67 - Correlação cruzada entre os resíduos e a entrada para o terceiro conjunto de validação.
Autocorrelação dos resíduos do modelo
Valor de autocorrelação calculado
0.8
0.6
0.4
0.2
0
-0.2
0
2
4
6
8
10
Atrasos
12
14
16
Figura 68 - Autocorrelação dos resíduos para o quarto conjunto de validação.
18
20
89
Correlação cruzada entre entrada e resíduos
0.1
Valor de correlação cruzada calculado
0
-0.1
-0.2
-0.3
-0.4
-0.5
-0.6
-20
-15
-10
-5
0
Atrasos
5
10
15
20
Figura 69 - Correlação cruzada entre os resíduos e a entrada para o quarto conjunto de validação.
Autocorrelação dos resíduos do modelo
Valor de autocorrelação calculado
0.8
0.6
0.4
0.2
0
-0.2
0
2
4
6
8
10
Atrasos
12
14
Figura 70 - Autocorrelação dos resíduos para o quinto conjunto de validação.
16
18
20
90
Correlação cruzada entre entrada e resíduos
0.2
Valor de correlação cruzada calculado
0.1
0
-0.1
-0.2
-0.3
-0.4
-0.5
-20
-15
-10
-5
0
Atrasos
5
10
15
20
Figura 71 - Correlação cruzada entre os resíduos e a entrada para o quinto conjunto de validação.
Nas Figuras 72, 74, 76, 78, 80 e 82 são apresentados os gráficos da função
de autocorrelação dos resíduos do modelo ARX [6 12 12], para o conjunto de
estimação e para os cinco conjuntos de validação.
Nas Figuras 73, 75, 77, 79, 81 e 83 são apresentados os gráficos da função
correlação cruzada entre a entrada e os resíduos do modelo ARX [6 12 12], para o
conjunto de estimação e para os cinco conjuntos de validação.
91
Autocorrelação dos resíduos do modelo
1
Valor de autocorrelação calculado
0.8
0.6
0.4
0.2
0
-0.2
-0.4
0
2
4
6
8
10
Atrasos
12
14
16
18
20
Figura 72 - Autocorrelação dos resíduos para o conjunto de estimação.
Correlação cruzada entre entrada e resíduos
0.3
Valor de correlação cruzada calculado
0.2
0.1
0
-0.1
-0.2
-0.3
-0.4
-20
-15
-10
-5
0
Atrasos
5
10
15
20
Figura 73 - Correlação cruzada entre os resíduos e a entrada para o conjunto de estimação.
92
Autocorrelação dos resíduos do modelo
1
Valor de autocorrelação calculado
0.8
0.6
0.4
0.2
0
-0.2
-0.4
0
2
4
6
8
10
Atrasos
12
14
16
18
20
Figura 74 - Autocorrelação dos resíduos para o primeiro conjunto de validação.
Correlação cruzada entre entrada e resíduos
0.4
Valor de correlação cruzada calculado
0.3
0.2
0.1
0
-0.1
-0.2
-0.3
-0.4
-0.5
-20
-15
-10
-5
0
Atrasos
5
10
15
20
Figura 75 - Correlação cruzada entre os resíduos e a entrada para o primeiro conjunto de validação.
93
Autocorrelação dos resíduos do modelo
1
Valor de autocorrelação calculado
0.8
0.6
0.4
0.2
0
-0.2
-0.4
0
2
4
6
8
10
Atrasos
12
14
16
18
20
Figura 76 - Autocorrelação dos resíduos para o segundo conjunto de validação.
Correlação cruzada entre entrada e resíduos
0.16
Valor de correlação cruzada calculado
0.14
0.12
0.1
0.08
0.06
0.04
0.02
0
-0.02
-0.04
-20
-15
-10
-5
0
Atrasos
5
10
15
20
Figura 77 - Correlação cruzada entre os resíduos e a entrada para o segundo conjunto de validação.
94
Autocorrelação dos resíduos do modelo
Valor de autocorrelação calculado
0.8
0.6
0.4
0.2
0
-0.2
0
2
4
6
8
10
Atrasos
12
14
16
18
20
Figura 78 - Autocorrelação dos resíduos para o terceiro conjunto de validação.
Correlação cruzada entre entrada e resíduos
0.3
Valor de correlação cruzada calculado
0.25
0.2
0.15
0.1
0.05
0
-0.05
-0.1
-20
-15
-10
-5
0
Atrasos
5
10
15
20
Figura 79 - Correlação cruzada entre os resíduos e a entrada para o terceiro conjunto de validação.
95
Autocorrelação dos resíduos do modelo
Valor de autocorrelação calculado
0.8
0.6
0.4
0.2
0
-0.2
0
2
4
6
8
10
Atrasos
12
14
16
18
20
Figura 80 - Autocorrelação dos resíduos para o quarto conjunto de validação.
Correlação cruzada entre entrada e resíduos
0.7
Valor de correlação cruzada calculado
0.6
0.5
0.4
0.3
0.2
0.1
0
-0.1
-20
-15
-10
-5
0
Atrasos
5
10
15
20
Figura 81 - Correlação cruzada entre os resíduos e a entrada para o quarto conjunto de validação.
96
Autocorrelação dos resíduos do modelo
Valor de autocorrelação calculado
0.8
0.6
0.4
0.2
0
-0.2
0
2
4
6
8
10
Atrasos
12
14
16
18
20
Figura 82 - Autocorrelação dos resíduos para o quinto conjunto de validação.
Correlação cruzada entre entrada e resíduos
0.7
Valor de correlação cruzada calculado
0.6
0.5
0.4
0.3
0.2
0.1
0
-0.1
-20
-15
-10
-5
0
Atrasos
5
10
15
20
Figura 83 - Correlação cruzada entre os resíduos e a entrada para o quinto conjunto de validação.
97
A partir dos gráficos obtidos, é possível perceber que na maioria dos casos as
funções de autocorrelação dos resíduos não apresentam o comportamento
desejado. Em contrapartida os resíduos gerados pelo modelo se distribuem
normalmente com média nula e variância  e2 .
98
CAPÍTULO VI
CONCLUSÃO
A partir dos resultados expostos neste trabalho foi possível concluir que a
modelagem de colunas de destilação, bem como a de outros sistemas industriais,
pode ser efetuada a contento através de modelos auto-regressivos.
Os modelos aqui obtidos não utilizam os recursos e técnicas mais avançadas
dentro da teoria de identificação de sistemas, mas mesmo com modelos lineares
reproduz de forma bastante aproximada o comportamento dinâmico do sistema
pesquisado.
Estes modelos são representações discretas lineares do sistema dinâmico em
estudo. As amostras utilizadas para o ajuste destes modelos representam o
funcionamento do sistema em condições normais de operação. Portanto os modelos
aqui expostos visam representar de maneira fiel o sistema próximo a este ponto
operacional.
A metodologia utilizada segue passos comuns aos mencionados em trabalhos
científicos da área. A partir de sua aplicação foram obtidos três modelos através de
dois métodos distintos – o critério de informação de Akaike (AIC) e a validação
cruzada – sendo eles: ARX [20 20 0], ARX [16 15 15] e ARX [6 12 12], tendo este
último atingido melhor desempenho em termos de minimização do MSE para os
cinco conjuntos de validação.
Tanto o modelo ARX [20 20 0] como o ARX [16 15 15] coincidiram em ser os
modelos de ordem mais elevada nas possíveis combinações. Apesar de ambos
obterem bons resultados, em termos de ajuste, para os dados de estimação e um
primeiro conjunto de dados de validação, não obtiveram o mesmo desempenho para
os demais conjuntos de validação, além de irem contra o princípio da parcimônia,
que orienta a escolha do modelo não só na sua capacidade de ajuste, como também
na quantidade reduzida de parâmetros.
O modelo que se mostrou mais eficiente em minimizar o índice MSE para os
cinco conjuntos de validação, foi o ARX [6 12 12]. O desempenho deste modelo foi
superior aos outros dois e além de ser o que possuía um menor número de
parâmetros, atendendo o princípio da parcimônia.
99
Os modelos obtidos foram ainda analisados através do cálculo de outros
índices de desempenhos como o MAE, MAPE, R2, etc. O desempenho dos modelos
foi semelhante à comparação com o MSE.
O desempenho preditivo dos modelos foi avaliado com o índice MSE para um
horizonte de predição variável. A análise destes resultados mostra que, à medida
que o horizonte de predição aumenta o MSE também aumenta, e mais uma vez os
melhores resultados são conseguidos pelo segundo modelo obtido pela validação
cruzada.
Os resíduos do modelo foram analisados com o artifício das funções de
autocorrelação para os resíduos e correlação cruzada entre os resíduos e a variável
de entrada. É sabido que uma boa modelagem gera resíduos com comportamento
de ruído branco, ou seja, a autocorrelação é nula para os atrasos superiores ou
iguais a um. Apesar de na maioria dos casos as funções de autocorrelação dos
resíduos não apresentam o comportamento desejado, os resíduos gerados pelo
modelo se distribuem normalmente com média nula e variância  e2 .
A contribuição deste trabalho foi a aplicação de uma metodologia já
consolidada para obtenção de um modelo matemático que represente de forma pelo
menos parcial o comportamento dinâmico do sistema estudado. E também a
utilização de cinco diferentes conjuntos de dados para a validação do modelo.
Colunas de destilação são equipamentos industriais que influenciam a
produção de diversos derivados necessários para a sociedade. Usinas de refino de
álcool combustível e álcool hidratado, usinas de refino de petróleo, ou seja, as
principais indústrias do setor energético e de combustíveis, indústrias que
fundamentam o processamento da matéria-prima em separação de seus
componentes e comercialização dos componentes de maior valor agregado.
Dentre as sugestões para trabalhos futuros fica a modelagem do sistema
utilizando modelos não lineares (NARX e NARMAX) e validação da capacidade
preditiva destes modelos para o mesmo horizonte utilizado neste trabalho.
Cabe destacar também a utilização dos modelos obtidos para subsidiar o
desenvolvimento de técnicas e métodos para a otimização e controle do sistema,
como por exemplo, desenvolvimento de controladores preditivos.
Modelos matemáticos dinâmicos, como os obtidos neste trabalho, podem ser
utilizados como base para outros estudos, que objetivem controle, estabilidade,
otimização e simulação de sistemas industriais.
100
REFERÊNCIAS
Aguirre, L. A. (2000) Introdução à identificação de sistemas: técnicas lineares e não
lineares aplicadas a sistemas reais. 1. ed. Belo Horizonte-MG: Editora UFMG, 554p.
Akaike, H. (1974) A new look at the statistical model identification. IEEE Transactions
on Automatic Control. 19:6 716-723.
Barroso, M. F. S., Saldanha R. R., Aguirre, L. A. (2002) Comparação de métodos
mono-objetivo em identificação caixa-cinza, TEMA – Tendências em Matemática
Aplicada e Computacional, publicação da Sociedade Brasileira de Matemática
Aplicada e Computacional. 3(2):43-52.
Belicanta, J. (2004) Coluna de para-destilação: análise das características
hidrodinâmicas e de eficiência de Murprhree. Dissertação (Mestrado em Engenharia
Química) – Campinas-SP - Universidade Estadual de Campinas – UNICAMP, 124p.
Borjas, S. D. M., Garcia, C. (2004) Modelagem de FCC usando métodos de
identificação por predição de erro e por sub-espaços. IEEE Latin America
Transactions. 2.
Box, G. E. P., Jenkins, G. M., Reinsel, G. C. (1994) Time Series Analysis:
Forecasting and Control. 3. ed. New Jersey: Prentice Hall, 592p.
Bravo, M. A. S. (2006) Métodos Matemáticos de Sistemas mediante Identificación
Paramétrica. III Seminario Internacional de Aplicaciones de la Matemática, Lima,
Peru.
Campos, M. C. M. M., Teixeira, H. C. G (2006) Controles típicos de equipamentos e
processos industriais. 1. ed. São Paulo-SP: Blücher, 416 p.
Carvalho, A. S., Guillermo, L. H. F. (2007) Modelagem de colunas de destilação
através de modelos auto-regressivos. IV Simpósio de Excelência em Gestão e
Tecnologia – SEGeT. Resende-RJ: AEDB.
101
Carvalho, A. S., Braga, A. G., Figueiredo, H. A. de (2004) Identificação por RNA da
relação vazão de alimentação por temperatura de topo de uma coluna de destilação
didática. Monografia (Graduação de Tecnologia em Automação Industrial) –
Campos-RJ – Centro Federal de Educação Tecnológica de Campos – CEFETCampos, 82p.
Crespo, L. S. (2000) Montagem, verificação de comportamento dinâmico e teste de
operacionalidade de uma torre de destilação piloto. Dissertação (Mestrado em
Engenharia Elétrica) – Vitória-ES – Universidade Federal do Espírito Santo – UFES,
120p.
Cruz, T. V. G. (2006) Identificação experimental de um modelo dinâmico de uma
microturbina a gás com câmara de combustão com baixa emissão de NOx.
Dissertação (Mestrado em Engenharia Mecânica) - Brasília-DF - Faculdade de
Tecnologia da Universidade de Brasília, 152p.
Dallagnol, V. A. F. (2005) Identificação de modelos ARMAX e NARMAX para um
poço de petróleo operando com por injeção contínua de gás. Dissertação (Mestrado
em Engenharia Elétrica) - Florianópolis-SC - Universidade Federal de Santa Catarina
- UFSC, 138p.
Doma, M. J., Taylor, P. A., Vermeer, P. J. (1996) Closed loop identification of MPC
models for MIMO processes using genetic algorithms and dithering one variable at a
time: application to an industrial distillation tower. Computers Chemical. Engineering.
20: SI035-SI040.
Furtado, E. C., Mendes, E. M. A. M., Nepomuceno, E. G., Silva, V. V. R. E. (2002)
Identificação de sistemas dinâmicos não-lineares contínuos utilizando modelos
NARMAX: estudo de caso de um forno a arco elétrico. Anais do XIV Congresso
Brasileiro de Automática, Natal-RN.
Gomide, R. (1988) Operações Unitárias Volume IV. 1. ed. São Paulo-SP: Edição do
Autor,187p.
102
Hovd, M., Michaelsen, R., Montin, T. (1997) Model Predictive Control of a Crude Oil
Distillation Column. Computers & Chemical Engineering. 21:S893-S897.
Júnior, C. A. V., Barbabela, B. A., (2006) Análise exploratória da aplicação de
modelos auto-regressivos na previsão do comportamento dos preços de derivados
de petróleo. Anais do IX SEMEAD – Seminários em Administração FEA-USP. São
Paulo:USP.
Jurado, F., Carpio, J. (2003) Modelado de la micro-turbina em el sistema de
distribución mediante estructuras NARX. Revista IEEE América Latina, 1(1).
Lara, J. M. V. (2005) Identificação de modelos para controle preditivo: aplicação a
uma planta de lodos ativados. Tese (Doutorado em Engenharia Elétrica) –
Campinas-SP - Universidade Estadual de Campinas - UNICAMP, 215p.
Ljung, L., Glad, T. (1994) Modeling of dynamic systems. 1. ed. New Jersey: Prentice
Hall, 368p.
Marangoni, C. (2005) Implementação de uma estratégia de controle com ação
distribuída em uma coluna de destilação. Tese (Doutorado em Engenharia Química)
– Florianópolis-SC – Universidade Federal de Santa Catarina - UFSC, 151p.
Moraes, C. A. S. (2004) Modelagem, controle e minimização do consumo de energia
de uma unidade de fracionamento de nafta. Dissertação (Mestrado em Engenharia
Elétrica) – Campinas-SP – Universidade Estadual de Campinas – UNICAMP, 193p.
Naegele, E. F. da S. (2000) Proposta para uma coluna de destilação didática:
mistura etanol-água. Dissertação (Mestrado em Engenharia Elétrica) – Vitória-ES –
Universidade Federal do Espírito Santo – UFES, 135p.
Rasovsky, E. M. (1973) Álcool: Destilarias. 1. ed. Rio de Janeiro-RJ: MIC, 384p.
Schneider, R., Noeres, C., Kreul, L.U., Górak, A. (2001) Dynamic modeling and
simulation of reactive batch distillation. Computers & Chemical Engineering. 25:169176.
103
Silva, C. A. M. (2006) Exploração de métodos de seleção de variáveis pela técnica
de regressão logística para análises de dados epidemiológicos. Dissertação
(Mestrado em Ciência Médicas) – Campinas-SP – Universidade Estadual de
Campinas – UNICAMP, 66p.
Souza, R. C., Camargo, M. E. (1996) Análise e previsão de séries temporais: os
modelos ARIMA. 1. ed. Santa Maria: Edição do autor, 242p.
Trautwein, B. J. (2004) Avaliação de métodos para previsão de consumo de água
para curtíssimo prazo: um estudo de caso em empresa de saneamento. Dissertação
(Mestrado em Engenharia de Produção e Sistemas) - Curitiba-PR - Universidade
Católica do Paraná, 123p.
Zanata, D. R. P. (2005) Desenvolvimento de um sensor virtual empregando redes
neurais para a medição da composição em uma coluna de destilação. Dissertação
(Mestrado em Engenharia) – São Paulo-SP - Escola Politécnica da Universidade de
São Paulo, 245p.