ads:
Universidade Federal de Uberlândia
Faculdade de Engenharia Química
Programa de Pós-Graduação em
Engenharia Química
Controle Ótimo de Sistemas
Algébrico-Diferenciais com
Flutuação do Índice Diferencial
Adriene Artiaga Pfeifer
Uberlândia - MG
2007
ads:
Livros Grátis
http://www.livrosgratis.com.br
Milhares de livros grátis para download.
Universidade Federal de Uberlândia
Faculdade de Engenharia Química
Programa de Pós-Graduação em Engenharia Química
Controle Ótimo de Sistemas
Algébrico-Diferenciais com
Flutuação do Índice Diferencial
Adriene Artiaga Pfeifer
Uberlândia - MG
2007
ads:
Universidade Federal de Uberlândia
Faculdade de Engenharia Química
Programa de Pós-Graduação em Engenharia Química
Controle Ótimo de Sistemas
Algébrico-Diferenciais com
Flutuação do Índice Diferencial
Adriene Artiaga Pfeifer
Dissertação de Mestrado apresentada ao Pro-
grama de Pós-Graduação em Engenharia
Química da Universidade Federal de Uber-
lândia como parte dos requisitos necessários
à obtenção do título de Mestre em Engenha-
ria Química, área de concentração Desenvol-
vimento de Processos Químicos.
Uberlândia - MG
2007
Dados Internacionais de Catalogação na Publicação (CIP)
P525c
Pfeifer, Adriene Artiaga, 1982-
Controle ótimo de sistemas algébrico-diferenciais com flutuação do
índice diferencial / Adriene Artiaga Pfeifer. - 2007.
125 f. : il.
Orientador: Valéria Viana Murata.
Dissertação (mestrado) – Universidade Federal de Uberlândia, Progra-
ma de Pós-Graduação em Engenharia Química.
Inclui bibliografia.
1. Controle de processos químicos
- Teses. I. Murata, Valéria Viana. II.
Universidade Federal de Uberlândia. Programa de Pós-Graduação em En-
genharia Química. III. Título.
CDU: 66.012
Elaborada pelo Sistema de Bibliotecas da UFU / Setor de Catalogação e Classificação
Agradecimentos
Agradeço aos meus pais, Randolfo e Eleni, que sempre me incentivaram e
aconselharam em todos os momentos da minha vida. As minhas irmãs, meus sobrinhos
e à minha avó Orávia que sempre me apoiaram.
À minha orientadora Valéria Viana Murata pela amizade, orientação e pelos
incentivos em todos os momentos.
Aos Professores Luís Cláudio Oliveira Lopes, Adilson José de Assis e Darci
Odloak pela atenção dispensada e opiniões relevantes na conclusão deste trabalho.
Ao meu amigo Fran Sérgio por toda a ajuda e orientação dada para a reali-
zação deste trabalho.
Aos meus amigos e companheiros de curso Ricardo, Davi, Sandra, Marcos,
Ricardo Pires, Andréia, Emília, Alaine, José Luiz, Fabiano, Janaína, Líbia, Gislaine,
Juliana, Patrícia, Danylo pelos momentos de descontração, alegria e incentivo.
À todos o funcionários da Faculdade de Engenharia Química pelo apoio.
Aos professores da Faculdade de Engenharia Química.
Ao CNPq pela concessão de bolsa de estudo.
“Ah... mas quem sou eu senão uma formiguinha das
menores, que anda pela terra cumprindo sua
obrigação."
Chico Xavier
SUMÁRIO
Lista de Figuras iv
Lista de Tabelas vi
Lista de Quadros vii
Lista de Abreviaturas viii
Resumo x
Abstract xii
1 Introdução 1
2 Problemas de Controle Ótimo Algébrico-Diferenciais (PCOADs) 5
2.1 Conceitos gerais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2 Condições necessárias de otimalidade para PCOADs . . . . . . . . . . . 8
2.3 Métodos de Redução do Índice Diferencial . . . . . . . . . . . . . . . . 20
2.3.1 O Código ALGO (UNGER et al., 1995) . . . . . . . . . . . . . . . 20
2.3.2 O Código PALGO (UNGER et al., 1995) . . . . . . . . . . . . . . 21
Sumário ii
2.4 Consistência de inicialização . . . . . . . . . . . . . . . . . . . . . . . . 24
2.5 Função Identificadora de Fases (FIF) . . . . . . . . . . . . . . . . . . . 25
2.6 Métodos Numéricos de solução de PCOADs . . . . . . . . . . . . . . . 26
2.6.1 Métodos Diretos . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.6.2 Métodos Indiretos . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.6.3 Métodos Mistos ou Híbridos . . . . . . . . . . . . . . . . . . . . 29
2.7 Sistemas Híbridos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.7.1 Definição . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.7.2 Sistemas Chaveados . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.7.3 Analogia entre PCOADs e Sistemas Chaveados . . . . . . . . . 38
3 Desenvolvimento de uma interface de integração de ferramentas para
o tratamento de PCOADs 39
3.1 Aspectos Gerais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.2 Pacote Maplets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.3 Descrição de ferramentas de tratamento de PCOADs . . . . . . . . . . 41
3.3.1 INDEX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.3.2 ACIG . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.3.3 OTIMA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.3.4 EVENTS: geração do sistema equivalente parametrizado pelos
eventos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.4 Aplicações da interface OpCol . . . . . . . . . . . . . . . . . . . . . . . 44
3.4.1 Caso 1 - Sistema de Índice 1 . . . . . . . . . . . . . . . . . . . . 45
3.4.2 Caso 2 - Reator CSTR com reação exotérmica de primeira ordem 48
3.4.3 Caso 3 - O Problema de Newton formulado como EDOs . . . . 51
3.4.4 Caso 4 - Sistema Chaveado . . . . . . . . . . . . . . . . . . . . 53
4 Solução numérica de PCOADs com flutuação do índice diferencial 58
4.1 Método de Solução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
Sumário iii
4.2 Estudo de Casos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.2.1 Caso 1 - Sistema Chaveado . . . . . . . . . . . . . . . . . . . . 59
4.2.2 Caso 2 - Reator semi-batelada isotérmico com reações paralelas
e restrição de seletividade . . . . . . . . . . . . . . . . . . . . . 64
4.2.3 Caso 3 - Reator Semi-Batelada Isotérmico . . . . . . . . . . . . 69
5 Conclusões e Sugestões 77
5.1 Conclusões . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
5.2 Sugestões . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
Referências Bibliográficas 80
A Manual - OpCol 84
B DIRCOL 91
B.1 Aspectos Gerais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
B.1.1 Subrotinas de entrada do DIRCOL 2.1 . . . . . . . . . . . . . . 92
B.1.2 Arquivos de saída padrões do DIRCOL 2.1 . . . . . . . . . . . . 93
B.1.3 Arquivos de saída opcionais do DIRCOL 2.1 . . . . . . . . . . . 93
B.2 Arquivo de entrada DATLIM . . . . . . . . . . . . . . . . . . . . . . . 93
B.3 Arquivo de entrada DATDIM . . . . . . . . . . . . . . . . . . . . . . . 95
B.4 Subrotina USER.f . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
LISTA DE FIGURAS
2.1 Grafos do código PALGO para o problema do pêndulo. . . . . . . . . 23
2.2 Dinâmica do estado contínuo definido em n
c+1
fases. . . . . . . . . . . . 32
3.1 Fluxograma do código EVENTS. . . . . . . . . . . . . . . . . . . . . . 44
3.2 Problema Isaac Newton. . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.1 Caso 1 - Variável de estado x
1
. . . . . . . . . . . . . . . . . . . . . . . . 61
4.2 Caso 1 - Variável de estado x
2
. . . . . . . . . . . . . . . . . . . . . . . . 61
4.3 Caso 1 - Variável de co-estado λ
1
. . . . . . . . . . . . . . . . . . . . . . 62
4.4 Caso 1 - Variável de co-estado λ
2
. . . . . . . . . . . . . . . . . . . . . . 62
4.5 Trajetória de estado. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.6 Caso 1 - Variável de controle. . . . . . . . . . . . . . . . . . . . . . . . 63
4.7 Caso 2 - Evolução das concentrações. . . . . . . . . . . . . . . . . . . . 67
4.8 Caso 2 - Função Objetivo. . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.9 Perfil ótimo da taxa de alimentação. . . . . . . . . . . . . . . . . . . . 68
4.10 Caso 2 - Perfil do Volume. . . . . . . . . . . . . . . . . . . . . . . . . . 69
4.11 Perfis das Concentrações. . . . . . . . . . . . . . . . . . . . . . . . . . . 73
4.12 Perfil ótimo da taxa de alimentação. . . . . . . . . . . . . . . . . . . . 74
Lista de Figuras v
4.13 Caso 3 - Perfil da Temperatura. . . . . . . . . . . . . . . . . . . . . . . 75
4.14 Caso 3 - Perfil do Volume. . . . . . . . . . . . . . . . . . . . . . . . . . 75
A.1 Tela de Entrada. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
A.2 Tela de Entrada do Código INDEX. . . . . . . . . . . . . . . . . . . . . 85
A.3 Tela de Entrada do Código OTIMA. . . . . . . . . . . . . . . . . . . . 85
A.4 Tela de Entrada do Código ACIG. . . . . . . . . . . . . . . . . . . . . . 86
A.5 Tela de Entrada do Código EVENTS. . . . . . . . . . . . . . . . . . . . 86
A.6 Entrada via arquivo do usuário. . . . . . . . . . . . . . . . . . . . . . . 87
A.7 Tela de seleção - INDEX. . . . . . . . . . . . . . . . . . . . . . . . . . . 88
A.8 Tela de Opções. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
A.9 Tela de seleção - OTIMA. . . . . . . . . . . . . . . . . . . . . . . . . . 89
A.10 Tela de seleção - EVENTS. . . . . . . . . . . . . . . . . . . . . . . . . . 89
A.11 Tela de geração de resultados. . . . . . . . . . . . . . . . . . . . . . . . 90
B.1 Estrutura do DIRCOL 2.1. . . . . . . . . . . . . . . . . . . . . . . . . . 92
LISTA DE TABELAS
4.1 Resultados do Caso 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.2 Variações nas estimativas iniciais. . . . . . . . . . . . . . . . . . . . . . 60
4.3 Parâmetros e condições Operacionais do Modelo (Caso 2). . . . . . . . 65
4.4 Caso 2 - Função Identificadora de Fases. . . . . . . . . . . . . . . . . . 66
4.5 Resultados do Caso 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
4.6 Variações nas estimativas iniciais. . . . . . . . . . . . . . . . . . . . . . 66
4.7 Parâmetros e condições Operacionais do Modelo (Caso 3). . . . . . . . 71
4.8 Caso 3 - Função Identificadora de Fases. . . . . . . . . . . . . . . . . . 71
4.9 Resultados do Caso 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.10 Variações nas estimativas iniciais. . . . . . . . . . . . . . . . . . . . . . 72
4.11 Variações nos parâmetros. . . . . . . . . . . . . . . . . . . . . . . . . . 72
LISTA DE QUADROS
3.1 Comandos do Maple. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.2 Saída do Código INDEX para o Caso 1. . . . . . . . . . . . . . . . . . 45
3.3 Saída do Código ACIG para o Caso 1. . . . . . . . . . . . . . . . . . . 47
3.4 Saída do Código INDEX para o Caso 2. . . . . . . . . . . . . . . . . . 48
3.5 Saída do Código ACIG para o Caso 2. . . . . . . . . . . . . . . . . . . 50
3.6 Conjunto de equações para a primeira fase. . . . . . . . . . . . . . . . . 54
3.7 Conjunto de equações para a segunda fase. . . . . . . . . . . . . . . . . 55
3.8 Conjunto de equações para a terceira fase. . . . . . . . . . . . . . . . . 56
LISTA DE ABREVIATURAS
CI - Condição Inicial
COLDAE - Collocation Differential Algebraic Equation Method
CSTR - Continuous Stirred Tank Reactor
DIRCOL - Direct Collocation Method
EADs - Equações Algébrico-Diferenciais
EDOs - Equações Diferenciais Ordinárias
EDPs - Equações Diferenciais Parciais
NLP - Programação Não Linear
PCOs - Problema de Controle Ótimo
PCOADs - Problemas de Controle Ótimo Algébrico-Diferencial
PMP - Princípio Mínimo de Pontryagin
PVI - Problema de Valor Inicial
PVC - Problema de Valor no Contorno
QP - Programação Quadrática
ix
SQP - Programação Quadrática Seqüencial
FIF - Função Identificadora de Fases (Switching Structure)
TPBV - Two-Point Boundary Value
Resumo
Os Problemas de Controle Ótimo, também chamados Problemas de Otimização Dinâmi-
ca, são formados por uma Função Objetivo a ser maximizada ou minimizada, associada
a conjuntos de equações algébricas e diferenciais que incluem restrições de igualdade
e de desigualdade nas variáveis de estado e de controle que caracterizam um sistema
de Equações Algébrico-Diferenciais (EADs). A extensão do ponto de vista algébrico-
diferencial de solução numérica aos PCOs, já amplamente utilizado na simulação de
processos devido à garantia de atendimento às restrições algébricas originais e implí-
citas na formulação e à eliminação das manipulações necessárias para transformar o
problema original num sistema de equações puramente diferenciais, caracteriza o cha-
mado Problema de Controle Ótimo Algébrico-Diferencial (PCOAD). Uma categoria
de PCOAD de especial interesse é a dos que incluem restrições de desigualdade, devido
à necessidade de conhecimento prévio da seqüência de ativações e desativações destas
restrições ao longo da trajetória e também dos instantes em que elas ocorrem, cha-
mados Eventos. As ativações/desativações das restrições causam flutuações no índice
diferencial e no número de graus de liberdade dinâmicos do PCOAD, exigindo técnicas
especiais de redução deste índice até um e o em prego de métodos numéricos eficientes
que garantam a convergência e estabilidade da solução.
Estes PCOADs com restrições de desigualdade são equivalentes a uma classe de pro-
blemas de otimização dinâmica híbridos, que associam comportamentos contínuos e
discretos (FEEHERY, 1998). Um tipo particular de PCO híbrido é aquele cujo estado
contínuo não apresenta saltos nos Eventos, chamado PCO Chaveado, para o qual Xu
e Antsaklis (2004) propõem uma metodologia de solução baseada na parametrização
dos Eventos com a especificação prévia da seqüência de subsistemas ativos, resultando
na solução de um problema de valor no contorno algébrico-diferencial em dois pontos,
formado pelas equações de estado, co-estado e de estacionariedade, condições de con-
torno e de continuidade e suas diferenciações, chamadas equações de sensibilidade.
Neste trabalho, esta abordagem indireta empregada para PCO Chaveados foi esten-
dida para PCOAD com restrições de desigualdade, com o objetivo de estimar também
xi
os Eventos, além das variáveis de controle, de estado e adjuntas. A abordagem de-
senvolvida por Xu e Antsaklis (2004) para PCO Chaveados foi implementada num
código específico utilizando o Maple 9.5, chamado EVENTS, com o objetivo de gerar
simbolicamente as equações baseadas na parametrização dos Eventos. Este código foi
incorporado a uma interface chamada OpCol, que reúne ferramentas de caracterização
de sistemas de E AD e de geração das condições de otimalidade segundo o Princípio
de Pontryagin estendidas para PCOAD de diferentes classes. As ferramentas de ca-
racterização são o INDEX de Murata (1996) que identifica simbolicamente o índice,
a resolubilidade e a consistência das condições iniciais e o ACIG de Cunha e Murata
(1999) que implementa o algoritmo de Gear para a redução do índice e geração do
sistema equivalente de índice 1. O OTIMA (GOMES, 2000; LOBATO, 2004) gera as
equações de Euler-Lagrange para PCOAD. Estas ferramentas foram inicialmente im-
plementadas em diferentes versões do Maple e todas foram atualizadas para a versão
9.5 utilizando o pacote Maplets que permite a entrada de dados através de janelas
interativas com o usuário, exigindo dele pouco conhecimento da sintaxe Maple. A
interface OpCol foi testada para quatro casos e para cada ferramenta foi criado um
banco de exemplos com problemas típicos da literatura que auxiliam o usuário na sua
utilização. Além disto, o método direto implementado no código DIRCOL estendido
para formulações multifásicas com estimativa dos Eventos e o método indireto com
Parametrização dos Eventos e abordagem algébrico-diferencial implementado num có-
digo MATLAB foram utilizados na solução numérica de três estudos de casos: um
PCO chaveado e 2 PCOAD de reatores batelada onde a variável de controle é a taxa
de alimentação do componente B: o primeiro tem reações paralelas e restrições de
seletividade com 3 fases de índices 1, 3 e 1 e o segundo restrições de segurança com 2
fases de índices 2 e 1 e respectivamente e foram descritos por Srinivasan e t al. (2003).
A mesma metodologia utilizada por estes autores foi empregada na obtenção de ex-
pressões analíticas para a variável de controle em cada fase necessárias no m étodo
indireto, compondo as chamadas Funções Identificadoras de Fase (FIF), a partir das
condições de otimalidade baseadas no Princípio de Pontryagin - especificamente a par-
tir da condição de estacionariedade e da identificação da restrição ativa que permitirá
a determinação da variável de controle - e da análise física do problema de modo a
descartar seqüências de ativação/desativação não apropriadas.
Os resultados obtidos pelo método indireto e pelo método direto são comparados entre
si para os 3 problemas citados, mostrando a viabilidade tanto da formulação multi-
fásica empregando o DIRCOL quanto o desempenho satisfatório do método indireto
com estimativa de Eventos, além da utilidade das ferramentas de caracterização de
EADs, de obtenção das condições de otimalidade e de parametrização dos eventos
disponibilizadas na interface Opcol.
Palavras-chave: Problemas de Controle Ótimo, Sistemas Híbridos, Sistemas Chave-
ados, Equações Algébrico-Diferenciais, DIRCOL.
Abstract
Optimal Control Problems (OCP), also known as Dynamic Optimization Problems,
consist of an Objective Function to be maximized or minimized, associated with a set of
differential and algebraic equations which include equality and inequality constraints
in the state or control variables and characterize a system of Differential-Algebraic
Equations (DAE). The differential-algebraic approach of numerical solution widely
used in process simulation due the guarantee of attendance of the implicit algebraic
constraints in the original formulation and the elimination of the necessary manipulati-
ons to transform the original problem into a purely differential system,was extended to
OCP characterizing the called Differential-Algebraic Optimal Control Problem (DA-
OCP). A category of DAOCP of special interest includes inequality constraints, due
the necessity of previous knowledge of the activations and deactivations sequence of
these constraints along the trajectory and also of the instants where they occur, named
Events.
This DAOCPs with inequality constraints is equivalent to a class of hybrid dynamic
optimization problems, where continuous and discrete behaviors are associated (FE-
EHERY, 1998). A particular type of hybrid OCP is that one where continuous state
does not present jumps in the Events, called Switched OCP, for which Xu e Antsaklis
(2004) considers a solution methodology based on the parameterization of Events with
a previous specification of active subsystems sequence, resulting in the solution of a
two-point boundary value differential-algebraic problem, formed by the state, co-state
and stationarity equations, boundary and continuity conditions and its differentiati-
ons, called sensitivity equations.
In this work, this indirect approach for Switched OCP was extended for DAOCP with
inequality constraints, with the objective to estimate the Events, along the control,
state and adjoint variables. The developed approach for Switched OCP described by
Xu e Antsaklis (2004) was implemented in a specific co de using Maple 9.5, called
EVENTS, with the objective to symbolically generate the equations based on the pa-
rameterization of Events. This code was incorporated in a interface named OpCol,
xiii
that collect characterization to ols of DAE systems and generation of the optimality
conditions extended Pontryagin’s Principle for PCOAD of different types. The cha-
racterization tools are the INDEX of Murata (1996) that symbolically identifies the
index, the resolubility and the consistency of initial conditions and the ACIG of Cunha
e Murata (1999) that implements the Gear’s algorithm for the index reduction and the
index 1 equivalent system generation. The OTIMA (GOMES, 2000; LOBATO, 2004) ge-
nerates the Euler-Lagrange equations for DAOCP. These tools had been implemented
initially in different versions of Maple and all had been update to 9.5 version using the
Maplets package that allows the data entry through interactive windows with the user,
demanding a little knowledge of the Maple syntax. The OpCol interface was tested for
four cases and for each tool a example data bank with typical problems of literature
was created to assist the user in its use. Moreover, the direct method implemented in
DIRCOL code was extended for multi-phases formulation with estimates of Events and
the indirect method with Events Parameterization and differential-algebraic approach
implemented in a Matlab code had been used for the numerical solution of three cases:
a switched OCP and 2 DAOCP of batch reactors where the control variable is the feed
rate of the component B - the first one has parallel reactions and selectivity constraints
with 3 phases of index 1, 3 and 1 and the second a safety constraint with 2 phases of
index 2 and 1 respectively and had been described by Srinivasan et al. (2003). The
methodology used by this authors was applied to attained analytical expressions for
the control variable in each phase necessary in indirect method, composing the called
Switching Functions, from the optimality conditions based in the Pontryagin’s Princi-
ple - specifically from the stationarity condition and the active constraint identification
that will allow the control variable determination - and of the physical analysis of the
problem in order to discard not appropriate activations/deactivations sequences.
The results obtained by indirect and direct methods are compared for the 3 cited pro-
blems, showing the viability as much of the multiphase formulation using the DIRCOL
and also the satisfactory performance of the indirect method with estimates of Events,
beyond the utility of the tools of characterization of EADs, of attainment of optimality
conditions and parameterization of Events available in Opcol interface.
Keywords: Optimal Control Problems, Differential-Algebraic Equations, Hybrid Sys-
tems, Switched Systems, DIRCOL.
CAPÍTULO 1
Introdução
Os Problemas de Controle Ótimo (PCOs), também chamados Problemas de Otimiza-
ção Dinâmica, são formados por uma Função Objetivo a ser maximizada ou minimi-
zada, associada a conjuntos de equações algébricas e diferenciais que incluem restrições
de igualdade e de desigualdade nas variáveis de estado e de controle que caracterizam
um sistema de Equaçõ es Algébrico-Diferenciais (EADs). A extensão do ponto de vista
algébrico-diferencial de solução numérica aos PCOs, já amplamente utilizado na simu-
lação de processos devido à garantia de atendimento às restrições algébricas originais
e implícitas na formulação e à eliminação das manipulações n ecess árias para transfor-
mar o problema original num sistema de equações puramente diferenciais, caracteriza
o chamado Problema de Controle Ótimo Algébrico-Diferencial (PCOAD).
Uma categoria de PCOAD de especial interesse é a dos que incluem restrições de
desigualdade, devido à necessidade de conhecimento prévio da seqüência de ativações
e desativações destas r estrições ao longo da trajetória e também dos instantes em que
elas ocorrem, chamados Eventos. As ativações/desativações das restrições causam flu-
tuações no índice diferencial e no número de graus de liberdade dinâmicos do PCOAD,
exigindo técnicas especiais de redução deste índice até um e o emprego de métodos
numéricos eficientes que garantam a convergência e estabilidade da solução.
Estes métodos podem ser diretos, indiretos ou mistos. Os métodos diretos uti-
lizam a parametrização na variável de controle ou a paramet rização nas variáveis de
controle e de estado seguida da solução do problema de Programação Não Linear
(NLP) resultante. A discretização nas variáveis de estado e de controle usando fun-
ções polinomiais lineares para o controle e cúbicas para o estado e definidas por partes,
seguida da solução de uma seqüência de NLP de dimensões crescentes por métodos
do tipo Programação Quadrática Seqüencial (SQP) para sistemas esparsos e densos
2
é utilizada no código DIRCOL de Stryk (1999), que calcula estimativas das variáveis
adjuntas e das funções multiplicadoras das restrições de igualdade e de desigualdade
e também dos E ventos. Os métodos indiretos aplicam o Princípio do Mínimo de Pon-
tryagin para gerar sistemas de EADs de valor no contorno definidas para cada fase
identificada pelos Eventos. Os métodos mistos combinam as vantagens dos métodos
diretos e dos métodos indiretos.
O estudo da abordagem algébrico-diferencial em PCOs na Faculdade de Enge-
nharia Química da Universidade Federal de Uberlândia iniciou-se com o trabalho de
Gomes (2000) através da aplicação do método indireto na solução de P COs e a criação
de um código que gera as equações adjuntas e a condição de mínimo através da aplica-
ção do Princípio de Pontryagin (OTIMA). Lobato (2004) apresentou uma abordagem
mista para a solução de problemas de otimização dinâmica, unindo as vantagens do
método direto e indireto de solução aplicados a problemas com restrições de desi-
gualdade. Vale ressaltar o trabalho de iniciação científica de Cunha e Murata (1999)
com a criação de um código de redução do índice diferencial de sistemas de equações
algébrico-diferenciais (ACIG). Com os estudos pode ser verificado que a solução numé-
rica de PCOs sujeitos a restrições de desigualdade apresenta desafios especiais porque
possui as seguintes características (FEEHERY, 1998):
• Exigência do conhecimento da seqüência e do número de ativações e desativações
das restrições ao longo da trajetória;
• Natureza combinatória para problemas com grande dimensão;
• A presença de restrições nas variáveis de estado, resultantes de limitações físicas,
econômicas, de segurança, geram EADs de Índice Diferencial Superior (> 1).
Nos trabalhos realizados nesta Faculdade, nota-se a necessidade de determinar
os momentos da ativação e desativação das restrições (Eventos). Com esse propó-
sito, a busca para métodos que estime esses valores para um melhor entendimento
da solução de PCO sujeitos a restrições de desigualdade, se faz necessária. Segundo
Feehery (1998) existe uma equivalência entre PCOs com restrições de desigualdade
nas variáveis de estado e os PCOs com comportamento contínuo e discreto acopla-
dos, chamados Sistemas Híbridos. De acordo com Barton et al. (1998), um fenômeno
híbrido pode ser classificado em chaveamentos ( switches) ou saltos (jumps) . Um cha-
veamento refere-se a mudanças discretas na forma funcional de f ( ˙z, z, u, t) = 0 como
conseqüência de eventos que ocorrem instantaneamente em um determinado ponto no
tempo. Estas diferentes formas possíveis do funcional são conhecidas com o “modos"de
um sistema híbrido. O momento da ocorrência dos eventos pode ser definido a priori
ou definido implicitamente pelo estado do sistema satisfazendo alguma condição. Um
salto refere-se a descontinuidades no estado da dinâmica do sistema como conseqüência
dos eventos. Em geral, saltos e chaveamentos podem ocorrer simultaneamente em um
evento. As equações (1.1) e (1.2) representam chaveamentos e saltos, respectivamente.
A primeira refere-se a um modelo que representa a influência de um vertedouro na
3
dinâmica de um tanque pulmão, onde h é o nível do tanque e h
weir
é a altura do
vertedouro e indica que enquanto o nível do líquido está acima do vertedouro, ele
escoa para fora do tanque e enquanto estiver abaixo, não ocorrerá a saída deste lí-
quido. Na segunda equação, x
−
denota lim
t→t
∗
−
x(t) e x
+
denota lim
t→t
∗
+
x(t), etc., e
j juntamente com f( ˙z, z, u, t) = 0 define o salto e pode representar, por exemplo, uma
certa quantidade de reagente alimentado num reator batelada num espaço de tempo
muito curto, que pode ser visto como um salto instantâneo na variável que representa
o número de moles do reagente no reator.
F
out
= k
weir
(h − h
weir
)
1,5
∀ t ∈ [0, t
f
] : h(t) > h
weir
F
out
= 0 ∀ t ∈ [0, t
f
] : h(t) ≤ h
weir
(1.1)
j( ˙x
−
, x
−
, u
−
, ˙x
+
, x
+
, u
+
, t
∗
) = 0 ∀ t
∗
∈ [0, t
f
] : s( ˙x
−
, x
−
, u
−
, t
∗
) = 0 (1.2)
Na definição apresentada por Xu e Antsaklis (2004), um sistema chaveado é um
tipo particular de sistema híbrido, mas cujo estado contínuo não apresenta saltos
nos eventos. Os autores propõem uma metodologia de transcrição de um sistema
chaveado geral em um PCO equivalente parametrizado p elos eventos, no qual uma lei
de mudança entre um sistema e outro é especificada. O problema equivalente tem a
propriedade de que os eventos são fixos em relação à nova variável independente do
tempo.
Considerando agora o PCO com flutuações de índice decorrente das ativações e
desativações de restrições ao longo da trajetória que provocam alterações na forma
funcional das equações do problema e admitindo como seqüência de mudanças entre
as fases, definidas pelos eventos, a Função Identificadora de Fases baseada na inter-
pretação física do problema, o objetivo desta dissertação é utilizar um algoritmo de
parametrização com relação aos eventos desenvolvido por Xu e Antsaklis (2004) para
sistemas chaveados, para resolver PCOs algébrico-diferenciais (PCOADs), utilizando
métodos de solução apropriados para problemas de índice 1. Desta forma, torna-se
possível a determinação dos eventos e a verificação do atendimento das condições de
ótimo. Para viabilizar o tratamento destes problemas conforme proposto, foram aper-
feiçoadas e integradas ferramentas automáticas de identificação e de redução do índice
diferencial, de obtenção das condições de otimalidade e de geração do sistema parame-
trizado que utilizam a computação algébrica. Além disto, os resultados obtidos pelo
procedimento de parametrização são comparados aos obtidos através do método direto
implementado no código DIRCOL, utilizando a opção de implementação multifásica
do problema.
Esta dissertação possui a estrutura conforme segue. O Capítulo 2 apresenta con-
ceitos fundamentais sobre otimização dinâmica, a abordagem algébrico-diferencial e
métodos de solução numérica, divididos em métodos diretos, indiretos, m istos e abor-
dagem híbrida. O Capítulo 3 trata da caracterização de sistemas de EADs, redução
4
do índice, geração das condições necessárias para o ótimo e a metodologia de solução
adotada aplicando a interface criada que une códigos que facilitam o estudo da abor-
dagem algébrico-diferencial em problemas de controle ótimo. O Capítulo 4 apresenta
os resultados obtidos pela abordagem adotada para problemas típicos. As conclusões
e sugestões para trabalhos futuros são apresentadas no Capítulo 5.
CAPÍTULO 2
Problemas de Controle Ótimo
Algébrico-Diferenciais (PCOADs)
2.1 Conceitos gerais
Um Problema de Controle Ótimo (PCO) consiste na determinação dos perfis das va-
riáveis de controle que maximizem ou minimizem uma medida de desempenho. O
aumento significativo de sua aplicação na indústria ocorreu na década passada devido
a popularidade de ferramentas de simulação dinâmica, associado a mercados compe-
titivos que demandam melhor desempenho da operação do processo sujeito a mais
restrições. Algumas aplicações de PCO em diversas áreas são apresentadas a seguir:
• Determinação de condições de operação ótimas para plantas químicas sujeitas a
restrições de segurança, de condições operacionais e ambientais (FEEHERY, 1998;
MODAK et al., 1986);
– Obtenção do perfil de temperatura que maximiza a seletividade de um de-
terminado produto;
– Maximização da quantidade de produto formado sujeito a restrições de
pressão etc;
2.1. Conceitos gerais 6
• Determinação da trajetória de robôs mecânicos (STRYK, 2000).
Para reunir conceitos fundamentais necessários para a melhor compreensão desta
dissertação, é apresentado a seguir um pequeno glossário de termos e definições.
Sistema de Equações Algébrico-Diferenciais (EADs) - Uma EAD é um sis-
tema de equações que pode ser escrito como
f(˙z, z, u, t) = 0 (2.1)
Nas EADs existem restrições algébricas na variável de estado z, que podem
aparecer explicitamente como em
g( ˙x, x, y, t) = 0 (2.2)
h(x, y, t) = 0 (2.3)
onde z=(x,y), ou podem aparecer implicitamente devido à singularidade de
df
dz
quando existem linhas não nulas.
Índice Diferencial de EADs - É o número mínimo de vezes que o sistema de
EADs ou parte dele deve ser diferenciado em relação a t para se determinar
˙z como uma função contínua de z (PANTELIDES, 1988).
EADs de Índice Superior e Redução de Índice - EADs de índice superior são
as equações com índice maior do que 1. Do ponto de vista da solução numérica é
desejável que o índice das EADs seja o menor possível devido à dificuldade asso-
ciada na solução dessas EADs, que pode ser comparada à solução de EDOs com
rigidez numérica. Entretanto, esta redução obtida através da simples diferenci-
ação das restriçõ es pode não satisfazer as restrições originais de maneira exata,
com sérias implicações quando elas envolvem propriedades físicas importantes.
Portanto, devem ser consideradas form as de reintroduzir restrições perdidas no
sistema, chamadas invariantes.
Princípio do Mínimo de Pontryagin (PMP) - É decorrente da imposição de que
o Hamiltoniano de um sistema contínuo sujeito a restrições de desigualdade nas
variáveis de controle deve ser minimizado para qualquer conjunto possível destas
variáveis de controle. É aplicável a problemas com variações fortes e restrições
de fim e foi estabelecido por Pontryagin em 1962 (BRYSON; HO, 1975).
Sistema Aumentado Correspondente - Conjunto de equações formado pelo sis-
tema original de EADs acoplado às equações decorrentes da aplicação do PMP.
Sistema Estendido Correspondente - Conjunto de equações formado pelo sis-
tema original de EADs acoplado às equações escondidas decorrentes da redução
do índice.
2.1. Conceitos gerais 7
Restrições terminais - Condições que os estados iniciais e finais devem satisfazer.
Restrições de trajetória - Condições que devem ser satisfeitas em todos os pontos
da trajetória.
Restrições interiores - Condições em pontos particulares ao longo da trajetória.
Arcos Singulares - Arcos onde a matriz de derivadas segundas do Hamiltoniano
com relação as variáveis de controle é singular. Os problemas de otimização
lineares podem apresentar impulsos nas variáveis de controle se elas não forem
submetidas a limites. Quando estes limites são impostos, a solução apresentará
períodos com a variável de controle nos limites de restrição e arcos singulares.
As variáveis de controle freqüentemente apresentam descontinuidades (cantos)
quando passam do arco restrito para o arco singular ou vice-versa (BRYSON; HO,
1975). Devido a estes arcos singulares, o perfil de controle não tem influência
direta sobre as condições de otimalidade do Hamiltoniano e portanto a deter-
minação da variável de controle exige condições adicionais que podem ser de
difícil manuseio (CUTHREL; BIEGLER, 1989). Os problemas de reatores batelada
alimentada, a manipulação de vazões em reatores químicos, colunas de destila-
ção, extratores e trocadores de calor são casos típicos de problemas com controle
singular (MODAK et al., 1986).
Eventos - Pontos isolados no tempo onde os comportamentos contínuos e discretos
interagem entre si. Nesses pontos podem ocorrer mudanças na forma funcional
das EADs e/ou nas trajetórias das variáveis de controle em cada fase (jumps e/ou
switches) decorrentes de descontinuidades nas trajetórias de estado e das variá-
veis adjuntas. Do ponto de vista algébrico-diferencial, estas mudanças causam a
flutuação do Índice Diferencial ao longo das fases.
Com o desenvolvimento da abordagem algébrico-diferencial na simulação de
processos químicos, uma nova perspectiva de solução pode ser adotada para os
PCOs. Adotando a perspectiva algébrico-diferencial, um Problema de Controle Ótimo
Algébrico-Diferencial (PCOAD) pode ser definido como (LOGSDON; BIEGLER, 1989):
min
u(t),t
f
J = Ψ(z(t
f
), t
f
) +
t
f
t
0
Φ(z, u, t)dt (2.4)
sujeito a:
g(˙z, z, u, t) = 0 (2.5)
h(˙z, z, u, t) ≤ 0 (2.6)
z
m
´
in
≤ z(t) ≤ z
m´ax
(2.7)
u
m
´
in
≤ u(t) ≤ u
m´ax
(2.8)
2.2. Condições necessárias de otimalidade para PCOADs 8
z ∈ R
n
z
(2.9)
u ∈ R
n
u
(2.10)
t ∈ R
1
(2.11)
Ψ : R
n
z
→ R
1
(2.12)
Φ : R
n
z
+n
u
+1
→ R
1
(2.13)
g : R
n
z
+n
u
+1
→ R
n
g
(2.14)
h : R
n
z
+n
u
+1
→ R
n
h
(2.15)
onde z(t) é o vetor das variáveis de estado, u(t) é o vetor das variáveis de controle,
o primeiro termo da função objetivo (J) é a condição de fim e o segundo termo é a
integral de um funcional dos vetores de controle e estado, h é o vetor de rest rições
de desigualdade e g é o vetor de restrições de igualdade. Os superscritos mín e máx
identificam, respectivamente, os limites inferior e superior das variáveis.
2.2 Condições necessárias de otimalidade para PCO-
ADs
As condições necessárias para o ótimo para problemas nos quais o sistema dinâmico
é descrito somente através de equações diferenciais ordinárias são bem estabelecidas
na literatura (BRYSON; HO, 1975). Nesta seção, estas condições serão generalizadas
para sistemas descritos por EADs, conforme demostrado por Renfro (1986) e Feehery
(1998).
Para o PCO considerado nesta seção, o vetor inicial das variáveis de estado é
conhecido (isto é, o vetor inicial não é determinado pela otimização), as trajetórias
de controle não estão sujeitas a restrição e as trajetórias de estado são restritas por
EADs. Este problema pode ser matematicamente expresso como (LOBATO, 2004;
GOMES, 2000):
min
u(t),t
f
J = Ψ(z(t
f
), t
f
) +
t
f
t
0
L(z, u, t)dt (2.16)
sujeito ao sistema de EADs
f(˙z, z, u, t) = 0 (2.17)
com condições iniciais dadas por
ϕ(˙z(t
0
), z(t
0
), t
0
) = 0 (2.18)
2.2. Condições necessárias de otimalidade para PCOADs 9
onde J(.), L(.), Ψ(.) → R, f(.), ϕ(.) → R
m
x
, z ∈ R
m
x
e u ∈ R
m
u
.
As variáveis de estado z nesta formulação incluem as variáveis de estado diferen-
ciais e algébricas.
A função Ψ na Eq.(2.16) pode ser expressa como:
Ψ(z(t
f
), t
f
) = Ψ(z(t
0
), t
0
) +
t
f
t
0
dΨ
dt
dt (2.19)
Desde que seja admitido que o tempo inicial t
0
e as condições de estado z(t
0
) são
fixos, a função objetivo pode ser expressa como
J =
t
f
t
0
¯
L(˙z, z, u, t)dt (2.20)
onde:
¯
L(˙z, z, u, t) =
dΨ
dt
+ L =
∂Ψ
∂t
+
∂Ψ
∂t
T
˙z + L (2.21)
Uma função objetivo aumentada é formada adicionando as restrições à função
objetivo através do uso das variáveis adjuntas λ(t):
¯
J =
t
f
t
0
¯
L(˙z, z, u, t) + λ
T
(t)f(˙z, z, u, t)
dt (2.22)
É conveniente definir o Hamiltoniano como:
H(˙z, z, u, λ, t) =
¯
L(˙z, z, u, t) + λ
T
(t)f(˙z, z, u, t) (2.23)
Para obter as condições necessárias para a otimalidade, é necessário definir a
variação do funcional. Para o funcional:
¯
J =
t
f
t
0
H(˙z, z, u, λ, t)dt (2.24)
2.2. Condições necessárias de otimalidade para PCOADs 10
o incremento do funcional é:
∆
¯
J =
t
f
t
0
[H(˙z + δ ˙z, z + δz, u + δu, λ + δλ, t) − H(˙z, z, u, λ, t)] dt +
t
f
+δt
f
t
f
H(˙z, z, u, λ, t)dt
(2.25)
Expandindo o incremento em Séries de Taylor ao redor do ponto (˙z(t), z(t), u(t))
e extraindo os termos que são lineares tem-se a variação de
¯
J:
δ
¯
J =
t
f
t
0
∂H
∂ ˙z
δ ˙z +
∂H
∂z
δz +
∂H
∂u
δu +
∂H
∂λ
δλ
dt + H(˙z, z, u, λ, t)δt
f
(2.26)
a qual pode ser simplificada integrando o primeiro termo por partes para se obter:
δ
¯
J =
t
f
t
0
∂H
∂z
−
d
dt
∂H
∂ ˙z
δz +
∂H
∂u
δu +
∂H
∂λ
δλ
dt+
+
∂H
∂ ˙z
t=t
f
δz(t
f
) + H(˙z, z, u, λ, t)δt
f
(2.27)
Usando a relação definida a seguir:
δz(t
f
) = δz
f
− ˙zδt
f
(2.28)
e substituindo a Eq.(2.28) na Eq.(2.27) fornece:
δ
¯
J =
∂H
∂ ˙z
t=t
f
δz
f
+
H −
∂H
∂ ˙z
˙z
t=t
f
δt
f
+
+
t
f
t
0
∂H
∂z
−
d
dt
∂H
∂ ˙z
δz +
∂H
∂u
δu +
∂H
∂λ
δλ
dt (2.29)
As condições necessárias de primeira ordem para o ótimo podem ser encontradas
2.2. Condições necessárias de otimalidade para PCOADs 11
fazendo a variação de
¯
J igual a zero. As condições são:
∂H
∂z
−
d
dt
∂H
∂ ˙z
= 0 (2.30)
∂H
∂u
= 0 (2.31)
∂H
∂λ
= 0 (2.32)
∂H
∂ ˙z
δz
t=t
f
+
H −
∂H
∂ ˙z
˙z
t=t
f
δt
f
= 0 (2.33)
As condições (2.30-2.33) definem um sistema de EADs de valor no contorno. Estas
condições podem ser simplificadas, expandindo os termos que incluem Ψ na Eq.(2.30).
δ
¯
J =
∂
∂z
∂Ψ
∂z
˙z +
∂Ψ
∂t
−
d
dt
∂
∂ ˙z
∂Ψ
∂z
˙z
=
∂
2
Ψ
∂z∂t
+
∂
2
Ψ
∂z
2
˙z
+
∂
2
Ψ
∂t∂z
+
∂
2
Ψ
∂z
2
˙z
= 0
(2.34)
Admitindo que as derivadas parciais de segunda ordem são contínuas, a ordem de
diferenciação pode ser mudada e a expressão inteira igualada a zero. Substituindo as
Eq.(2.23) e Eq.(2.34) e as Eq.(2.30)-Eq.(2.33) podem ser simplificadas:
∂L
∂z
+ λ
T
∂f
∂z
−
˙
λ
T
∂f
∂ ˙z
− λ
T
d
dt
∂f
∂ ˙z
= 0 (2.35)
∂L
∂u
+ λ
T
∂f
∂u
= 0 (2.36)
f(˙z, z, u, t) = 0 (2.37)
λ
T
∂f
∂ ˙z
+
∂Ψ
∂z
t=t
f
δz
f
+
∂Ψ
∂t
+ L + λ
T
f − λ
T
∂f
∂ ˙z
˙z
t=t
f
δt
f
= 0 (2.38)
Estas condiçõ es são uma generalização das condições necessárias, também chama-
das de equações de Euler-Lagrange, para a otimização dinâmica de sistemas de EADs.
A seguir estas equações serão mostradas para casos particulares que envolvem proble-
mas com tempos finais fixos ou livres e problemas com variáveis de estado espec ificadas
no tempo final.
a) Problemas com o Tempo Final Fixo
Se t
f
é fixo então δt
f
é igual a zero na Eq.(2.38). Se a variável de estado não for
2.2. Condições necessárias de otimalidade para PCOADs 12
especificada no tempo final, as condições no ponto final têm que satisfazer a:
λ
T
∂f
∂ ˙z
+
∂Ψ
∂z
t=t
f
δz
f
= 0 (2.39)
como δz
f
é arbitrário, δz
f
= 0 e para satisfazer a Eq.(2.38):
λ
T
∂f
∂ ˙z
+
∂Ψ
∂z
t=t
f
= 0 (2.40)
b) Problemas com o Tempo Final Livre
Como t
f
é livre, a suposição de que δt
f
= 0 não pode ser feita. Neste caso além
das condições dadas pelas equações Eq.(2.35) - Eq.(2.37), para o caso em que se tem as
variáveis de estado fixas no tempo final ou variáveis de estado livre, o sistema deverá
atender a seguinte condição:
∂Ψ
∂t
+ L + λ
T
f − λ
T
∂f
∂ ˙z
˙z
t=t
f
= 0 (2.41)
c) Algumas variáveis de estado especificadas no tempo final fixo
Seja o problema de otimização definido pelas Eq.(2.16) - Eq.(2.18), com algumas
variáveis de estado especificadas em t = t
f
. Se z
i
, o i-ésimo componente do vetor de
estado z, é definido em t = t
f
, então como a variação δz
i
(t
f
) na Eq.(2.38) não pode ser
nula, é necessário que a Eq.(2.39) seja satisfeita. A Eq.(2.36), ∂H/∂u = 0 necessita de
uma condição adicional para o problema com restrição final. No presente caso, δu(t)
não é completamente arbitrário, o conjunto admissível de δu(t) é sujeito às restrições:
δz
i
(t
f
) = 0 i = 1, ..., q (2.42)
Um conjunto admissível de δu(t) pode ser definido como os valores de δu(t) que
satisfazem todas as restrições do problema, com o por exemplo a Eq.(2.42).
Desde que z
i
(t
f
) especificado para i = 1, . . . , q, é consistente considerar que:
ϕ = ϕ(z
q+1
, . . . , z
n
)
t=t
f
(2.43)
As Eq.(2.35)-Eq.(2.38) permanecem inalteradas para este caso, apenas a condição
2.2. Condições necessárias de otimalidade para PCOADs 13
de contorno em t = t
f
passa a ser dada por:
λ
i
(t
f
) =
0 j=1,. . . ,q
∂φ
∂z
j
t=t
f
j=q+1,. . . ,n
(2.44)
d) Sistemas com funções de variáveis de estado especificadas no tempo
final fixo
Seja o problema de otimização definido pelas Eq.(2.16) - Eq.(2.18) sujeito a restri-
ção definida pela Eq.(2.45), de dimensão q, função das variáveis de estado e com o
valor definido no tempo final.
ϕ(z(t
f
), t
f
) = 0 (2.45)
A Eq.(2.45) pode ser adicionada à função objetivo através de multiplicadores de
Lagrange ν, um vetor de dimensão q.
J = ψ(z(t
f
), t
f
) + ν
T
ϕ(z(t
f
), t
f
) +
t
f
t
0
L(z, u, t)dt (2.46)
Definindo:
Ψ = ψ(z(t
f
), t
f
) + ν
T
ϕ(z(t
f
), t
f
) (2.47)
O conjunto de parâmetros ν devem ser escolhidos para satisfazer a Eq.(2.45).
Portanto as condições necessárias são dadas pelas Eq.(2.35-2.38) e por:
λ
T
(t
f
) =
∂Ψ
∂z
+ ν
T
∂ϕ
∂z
t=t
f
(2.48)
A Eq.(2.36) determina o vetor u(t), as Eq.(2.35) - Eq.(2.40) formam um sistema de
EADs de valor no contorno com q parâmetros ν para serem determinados na Eq.(2.48)
tal que a Eq.(2.45) seja satisfeita.
e) Problemas com Restrição de Trajetória
Nesta seção são considerados os problemas de otimização com restrições de traje-
2.2. Condições necessárias de otimalidade para PCOADs 14
tória, que se aplicam a pontos intermediários ou sobre toda a trajetória.
i) Res trições de Igualdade na Variável de Controle
O problema é definido como:
min
u(t),t
f
J = Ψ(z(t
f
), t
f
) +
t
f
t
0
L(z, u, t)dt (2.49)
ϕ(z(t
f
), t
f
) = 0 (2.50)
f(
˙
z, z, u, t) = 0 (2.51)
C(u, t) = 0 (2.52)
onde C(u, t) é o vetor de restrições de igualdade na variável de controle. Neste caso
u(t) é um vetor de variáveis de controle de dimensão m ≥ 2 e C é uma função escalar.
O Hamiltoniano (Eq.(2.23)) será redefinido como:
H(
˙
z, z, u, λ, t) =
L(
˙
z, z, u, t) + λ
T
(t)f(
˙
z, z, u, t) + µC (2.53)
As condições necessárias Eq.(2.35, 2.37 e 2.38) permanecem inalteradas, sendo
que Eq.(2.36) será redefinida como:
∂L
∂u
+ λ
T
∂f
∂u
+ µ
∂C
∂u
= 0 (2.54)
O sistema formado pelas condições necessárias e pela Eq.(2.52) representa as m+1
condições para determinar o m-ésimo componente do vetor de controle u(t) e a função
escalar µ(t).
ii) Restrições de igualdade nas variáveis de controle e estado
Neste caso a Eq.(2.52) será substituída por:
C(z, u, t) = 0 (2.55)
Aqui as condições obtidas para a seção anterior po dem ser aplicadas, só que um
novo termo será acrescentado a Eq.(2.35), sendo reescrita como:
∂L
∂z
+ λ
T
∂f
∂z
+ µ
T
∂C
∂z
−
˙
λ
T
∂f
∂ ˙z
− λ
T
d
dt
∂f
∂ ˙z
= 0 (2.56)
2.2. Condições necessárias de otimalidade para PCOADs 15
iii) Res trições de Igualdade na Variável de Estado
Se a restrição não tiver dependência explícita na variável de controle, ocorre uma
complexidade adicional. Seja a restrição:
S(z, t) = 0 (2.57)
Se esta restrição é aplicada sobre todo o intervalo t
0
≤ t ≤ t
f
, a derivada temporal
da restrição é nula ao longo da trajetória:
dS
dt
=
∂S
∂t
+
∂S
∂z
˙
z =
∂S
∂z
f(z, u, t) = 0 (2.58)
A Eq.(2.58) pode ou não revelar a dependência da variável de controle u. Se a
Eq.(2.58) revelar a dependência de u, então pode ser tratada como uma restrição do
tipo da Eq.(2.55). Para isto, deve-se eliminar um componente de z como função dos
n-1 componentes remanescentes usando a Eq.(2.57) como uma condição de contorno
em t = t
0
ou t = t
f
. Se a Eq.(2.58) não revelar a variável de controle explicitamente,
deve-se repetir o processo de diferenciar a equação até que a variável de controle u seja
revelada explicitamente. Surge então o conceito de ordem da restrição de igualdade
na variável de controle, que é definida como o número de vezes que a restrição deve ser
diferenciada para que se obtenha a dep endência da variável de controle u. A q-ésima
derivada temporal da restrição da Eq.(2.57), é representada por:
S
q
(z, u, t) = 0 onde S
q
(z, u, t) =
d
q
S
dt
q
(2.59)
Neste caso, q componentes de z devem ser eliminados manuseando os (n-q) com-
ponentes remanescentes, usando as q relações:
S(z, t)
S
1
(z, t)
.
.
.
S
q−1
(z, t)
= 0 (2.60)
ou adicionando a Eq.(2.60) como um conjunto de condições de contorno em t = t
0
ou
t = t
f
.
A existência de restrições de igualdade nas variáveis de estado num PCOAD pode
aumentar o índice diferencial. É interessante observar que uma restrição de igualdade
2.2. Condições necessárias de otimalidade para PCOADs 16
pode surgir quando um PCO formado por Equações Diferenciais Implícitas do tipo:
F (
˙
x, x, u, p, t) = 0 (2.61)
onde u é a variável de controle, p é o parâmetro e x é a variável de estado, é reescrito
definindo uma nova variável de controle u
∗
= [u, v] tal que o Sistema Aumentado
definido como:
˙
x = v (2.62)
G(x, u, p, v, t) = 0 (2.63)
torna-se um sistema de EADs onde a equação Eq.(2.63) passa a ser a nova restrição
algébrica.
iv) Restrições de desigualdade na variável de controle
Neste caso o problema de ot imização considerado para simplificar a análise, é
de tempo fixo e sem restrição definida no ponto final, sujeito a uma restrição de
desigualdade do tipo:
C(u, t) ≤ 0 (2.64)
Se o Hamiltoniano for definido por: H
0
= λf + L (Eq.(2.2)), considerando que
δz
i
= 0, ∂H/∂λ = 0 e para simplificações em que os coeficientes de δz são iguais a
zero, pode ser reescrita como:
δJ =
t
f
t
0
H
0
u
δudt
t
f
t
0
H
0
u
(z, λ, u, t)dt (2.65)
As condições necessárias para este problema são as Eq.(2.35) - Eq.(2.38). Para
u(t) ser minimizado, δJ ≥ 0 para todo o conjunto admissível de u(t). Isto implica
que δH
0
≥ 0 para todo t e todo conjunto admissível u(t). Os pontos onde ocorrem os
valores ótimos de u(t) tem a seguinte propriedade:
δH
0
= H
0
u
δu ≥ 0
δC = C
u
δu ≤ 0 (2.66)
v) Se o Hamiltoniano for definido por:
2.2. Condições necessárias de otimalidade para PCOADs 17
H = λ
T
f + L + µ
T
C (2.67)
Aqui as condições necessárias são dadas pelas equações Eq.(2.35) - Eq.(2.38), com
uma condição adicional:
µ =
≥ 0 se C = 0
= 0 se C < 0
(2.68)
A exigência do multiplicador ser positivo quando C = 0, pode ser interpretada
como uma condição para que o gradiente H
u
≡ λ
T
f
u
+L
u
seja obtido somente violando
as restrições.
Quando a restrição de desigualdade torna-se ativa em algumas porções da traje-
tória, o problema de otimização apresenta arcos com restrição e sem restrição. Nos
pontos de junção, entre os arcos restritos e não restritos, a variável de controle, pode
ou não ser contínua. Se u(t) for descontínuo, o ponto é chamado de canto (corner). Os
Corners podem ocorrer em qualquer ponto da trajetória, mas não são mais prováveis
de ocorrer nos pontos de junção que no meio do arco sem restrição. A priori não
existe método para determinar a existência dos cantos. Se u(t) for contínuo nos pon-
tos de junção, seguido da continuidade de λ, ∂H/∂u, H então µ(t) também é contínuo.
vi) Restrições de Desigualdade nas Variáveis de Controle e de Estado
Neste caso a restrição de desigualdade é dada por:
C(z, u, t) ≤ 0 (2.69)
O problema é manuseado da mesma forma que o problema de funções de variáveis
de estado espe cificadas no tempo final fixo.O Hamiltoniano é definido como Eq.(2.67):
H = λ
T
f + L + µ
T
C (2.70)
Aqui
µ =
>0 se C = 0
= 0 se C < 0
(2.71)
e as equações de Euler-Lagrange são assim definidas:
2.2. Condições necessárias de otimalidade para PCOADs 18
˙
λ
T
= −
∂H
∂z
−L
z
− λ
T
f
z
− µC
z
se C = 0
−L
z
− λ
T
f
z
se C < 0
(2.72)
A condição que determ ina u(t) é dada por:
H
u
≡ L
u
+ λ
T
f
u
+ µC
u
(2.73)
Quando C < 0, µ = 0 a Eq.(2.74) determina u(t). Para C = 0, a Eq.(2.71) e
Eq.(2.73) determinam u(t) e µ(t) simultaneamente.
vii) Restrições de Desigualdade nas Variáveis de Estado
Seja a restrição de desigualdade dada por:
S(z, t) ≤ 0 (2.74)
Considere que S e u são escalares. A derivada temporal da Eq.(2.74) e a subs-
tituição de
˙
z devem ser realizadas até que seja revelada a dependência explícita de
u. Se forem necessárias q derivadas temporais, a Eq.(2.74) é denominada restrição de
desigualdade na variável de estado de q-ésima ordem. Aqui S
q
(z, u, t) representa a
q-ésima derivada temporal total de S. O Hamiltoniano é definido como:
H = L + λ
T
f + µS
q
(2.75)
onde a restrição está ativa se:
S
q
= 0 ⇒ S = 0 (2.76)
A restrição não está ativa se:
µ = 0 ⇒ S < 0 (2.77)
As condições necessárias são dadas pelas equações Eq.(2.35) - Eq.(2.38), substi-
tuindo C por S
q
. A condição necessária para µ(t) se a restrição está ativa é
µ(t) ≥ 0 ⇒ S = 0 (2.78)
2.2. Condições necessárias de otimalidade para PCOADs 19
Neste caso também ocorre o aparecimento de arcos restritos e arcos não restritos.
Os arcos restritos devem ser tangentes aos arcos não restritos nos pontos de junção, o
que faz com que apareçam descontinuidades nos pontos de entrada e saída de qualquer
arco. Daqui surgem as restrições de tangência, que são denominadas restrições de
contorno em pontos interiores, definidas por:
N(z, t)
S(z, t)
S
1
(z, t)
.
.
.
S
q−1
(z, t)
= 0 (2.79)
Pode-se escolher o ponto de entrada em vez do ponto de saída para satisfazer estas
restrições interiores. Então os λ e o H serão descontínuos no ponto de entrada t = t
1
e contínuos no ponto de saída. O vetor N(z, t) (Eq.(2.79)) representa as condições de
salto. Bryson e Ho (1975), através de manipulações semelhantes, mostraram que as
condições de salto no ponto de entrada podem ser obtidas por:
λ
T
(t
−
1
) = λ
T
(t
+
1
) + π
T
∂N
∂z(t
1
)
(2.80)
H
T
(t
−
1
) = H
T
(t
+
1
) + π
T
∂N
∂t
1
(2.81)
onde t
−
1
significa o tempo anterior a t
1
e t
+
1
o tempo logo após o t
1
, π é um vetor
de multiplicadores de Lagrange de dimensão q usados para adicionar as condições de
junção Eq.(2.81) à função objetivo, que são determinados de tal forma que atendam
a estas condiçõ e s.
A solução de PCOs com Restrições de Desigualdade apresenta desafios especiais,
porque exige o conhecimento da seqüência e do número de ativações e desativações
ao longo da trajetória. Quando a quantidade de restrições é reduzida, geralmente é
possível determinar esta seqüência examinando a solução do problema sem restrições.
Entretanto, a presença de um número maior de restrições torna o problema de natureza
combinatória (FEEHERY, 1998).
Uma fonte de flutuação do índice ocorre quando a variável de controle não pode
ser explicitada em termos das variáveis de estado e das variáveis adjuntas a partir
da condição estacionária, provocando arcos singulares durante a solução ( LOGSDON;
BIEGLER, 1989). É importante ressaltar que, se o número de restriçõe s de desigualdade
for maior que o número de variáveis de controle, resulta em EADs de índice superior
independente da restrição de desigualdade estar ativa ou não (FEEHERY, 1998).
2.3. Métodos de Redução do Índice Diferencial 20
2.3 Métodos de Redução do Índice Diferencial
A redução do índice diferencial é necessária para a solução numérica de EADs gerais
de índice superior. Segundo Gear e Petzold (1984) apenas os sistemas de índice menor
ou igual a 1 e as EADs lineares com coeficiente constantes de índice qualquer podem
ser resolvidos por métodos para EDOs. Para as EADs lineares em geral e EADs não
lineares de índice superior não existem métodos numéricos que as resolvam incondici-
onalmente. Isto ocorr e porque o conjunto de soluções das EDOs subordinadas obtidas
pela redução é maior do que o conjunto de soluções das EADs originais. Como as
equações algébricas das EADs se tornam apenas implícitas nas EDOs subordinadas,
a solução numérica freqüentemente causa flutuações nestas restrições algébricas (que
não se preservam com a discretização exceto se são lineares), e levam a instabilida-
des. Por isto, vários algoritmos de redução do índice usam as chamadas técnicas de
estabilização de restrições.
De modo geral, os algoritmos de redução p odem ter ou não etapas de diferenciação,
e fornecerem ou não sistemas de maior dimensão. Os algoritmos de Gear e Petzold
(1984), Gear (1988) e Bachmann et al. (1990) utilizam a diferenciação das restrições
algébricas com relação ao tempo e a linearidade das novas equações introduzidas pela
diferenciação com relação às derivadas de ordem superior, para eliminar os termos
onde estas derivadas de ordem superior aparecem. Assim, são obtidas novas equações
que são funções das variáveis de estado e das variáveis algébricas. As novas equações
algébricas obtidas durante a r edução são incorporadas ao sistema nos algoritmos de
Gear (1988) e Bachmann et al. (1990).
Faz-se necessária uma análise prévia do problema original que inclua a determi-
nação do índice diferencial e a consistência da inicialização. Para isso, dois algoritmos
de redução do índice que utilizam propriedades estruturais para caracterizar EADs são
brevemente apresentados na seqüência.
2.3.1 O Código ALGO (UNGER et al., 1995)
O Código ALGO é uma ferramenta estrutural, implementada na linguagem FOR-
TRAN, e é baseada na representação estrutural do algoritmo de Gear (GEAR, 1988)
proposto para a redução do índice diferencial. Este algoritmo aplica-se a EADs
gerais F(z(t), z
′
(t), t) = 0 de dimensão n tendo como critério de resolubilidade o
det(λ(F
z
′
) + F
z
) ser não identicamente nulo para todo λ. O índice diferencial do
sistema é igual ao número de iterações deste algoritmo para transformar as EADs
num sistema de EDOs.
O ALGO é composto pelas seguintes subrotinas:
INPUT Lê a dimensão das EADs. Lê os padrões das matrizes de derivadas pad(F
z
′
)
2.3. Métodos de Redução do Índice Diferencial 21
e de algébricas pad(F
z
). Estes dados podem ser fornecidos pelo usuário por
entrada manual, via arquivo e p ode ser gravado em um arquivo especificado pelo
usuário.
SOLVAB Determina o posto estrutural do pad(F
z
′
)+ pad(F
z
) utilizando RANKDET
e MC21A. A sub-rotina MC21A coloca o máximo número de elementos não
nulos na diagonal principal e retorna à RANKDET que determina a soma destes
elementos. Se esta soma for igual a dimensão do sistema, a matriz estrutural
será não singular e o sistema poderá ser resolvido.
PARTITION Executa a partição do sistema.
INVERSION Realiza a Eliminação Gaussiana Estrutural.
SUBK1X1 Executa a substituição.
IND2CAUS Nomeia as variáveis e as equações que poderiam ser a causa do í ndice
superior.
ALLONIN Cria um relatório para EADs de índice maior que um.
SUBK1CDG Executa a substituição.
DIFFSUBX Faz a diferenciação parcial conforme a especificação da linearidade for-
necida pelo usuário.
FAIC Determina se um conjunto de variáveis especificadas para ter seus valores ini-
ciais arbitrariamente escolhidos permite a inicialização consistente do sistema.
Utiliza as subrotinas FEA SRANK e MC21A.
2.3.2 O Código PALGO (UNGER et al., 1995)
A ferramenta estrutural PALGO é baseada no algoritmo estrutural de Pantelides
(1988). Segundo Pantelides, um sistema semi-explícito geral não apresenta proble-
mas de inicialização consistente se a matriz Jacobiana do sistema de equações for não
singular. As restrições adicionais ao sistema de equações seguem o princípio de que um
sub-conjunto de equações deve ser diferenciado se e somente se o grupo de linhas da
matriz correspondentes a estas equações forem linearmente dependentes e se existirem
linhas linearmente independentes dentro de sub-conjuntos deste grupo. Para o estudo
do algoritmo algumas definições são importantes:
• Nós: São as formas de representação das variáveis e equações no gráfico “bipar-
tite".
• E-nós: Conjunto de nós referentes às EADs.
2.3. Métodos de Redução do Índice Diferencial 22
• V-nós: Conjunto de nós referentes às variáveis x’, y.
• Nó exposto: Nó que não aparece em nenhuma linha combinada.
• Associação: Conjunto de linhas (i-j), chamadas linhas combinadas, cujos nós i e
j aparecem em apenas uma das linhas.
• Associação completa: Associação que não deixa nenhum E-nós exposto.
• Associação parcial: Ass ociação que deixa E-nós expostos.
• Caminho estendido: Caminho com nós expostos nas duas extremidades.
• Sub-conjunto estrutural e minimamente singular (ems): Subconjunto de equa-
ções singular que não tem nenhum de seus próprios subconjuntos estruturalmente
singular.
O algoritmo de Pantelides detecta o número mínimo de equações que devem ser
diferenciadas para que se obtenha resolução estruturalmente consistente do sistema.
Para tanto, este algoritmo utiliza algumas das seguintes listas de associação:
• A(.) - expressa a relação entre as variáveis e suas derivadas. A[1] = 5 indica que
a quinta variável é a derivada da primeira.
• B(.) - relaciona as equações criadas por diferenciação com as suas criadoras. B[5]
= 6 indica que a sexta equação foi criada pela diferenciação da quinta equação.
• ASSIGN(.) - contém a equação a qual c ada variável está associada. ASSIGN(5)
= 1 indica que a variável 5 está ligada pelas linhas do grafo à primeira equação.
As linhas do grafo proposto por Pantelides representam a existência das variáveis
em cada uma das equações. Os números apresentados em E-nós e V-nós representam
a posição dos vetores de equações e variáveis, respectivamente (Fig. 2.1).
Exemplo: Pêndulo de Pantelides.
F
1
= x
′
1
− x
3
= 0
F
2
= x
′
2
− x
4
= 0
F
3
= x
′
3
− x
1
x
5
= 0
F
4
= x
′
4
− x
2
x
5
+ G = 0
F
5
= x
2
1
+ x
2
2
− L
2
= 0
onde os termos L e G são constantes no tempo.
Equações = [F
1
, F
2
, F
3
, F
4
, F
5
]
Variáveis = [x
1
, x
2
, x
3
, x
4
, x
5
, x
′
1
, x
′
2
, x
′
3
, x
′
4
]
O PALGO é composto pelas seguintes subrotinas:
2.3. Métodos de Redução do Índice Diferencial 23
1
2
3
4
5
1
2
3
4
5
6
7
8
9
E-nós
V-nós
Figura 2.1: Grafos do código PALGO para o problema do pêndulo.
NEWMO - Examina as equações e estabelece os eixos entre os respectivos nós; os
eixos são criados somente com as variáveis de maior derivada; desta forma esta
subrotina consegue eliminar os A(.)=0 do algoritmo de Pantelides.
AUGMP - Através de uma modificação da subrotina MC21A, chamada AUGAD
que retorna o PATHFOUND como falso ou verdadeiro, determina qual das equa-
ções devem ser diferenciadas. A AUGAD é responsável por criar as associações
entre variáveis e equações (ASSIGN(.)) e indicar (colorir), através destas asso-
ciações, quais são as variáveis e equações que devem ser diferenciadas. Usando
o mesmo critério de MC21A para o cálculo do posto estrutural, a AUGAD de-
termina as permutações necessárias para que o PATHFOUND seja verdadeiro.
DICOLSUB - Se o PATHFOUND retornado por AUGMP é falso, então existe um
nó exposto indicando que existe uma equação a ser diferenciada. Esta subrotina
insere a equação resultante da diferenciação estrutural do grafo e retorna o sis-
tema alterado à AUGMP. Isto é feito até que o sistema não tenha nenhum nó
exposto, ou seja, até que PATHFOUND seja verdadeiro.
As subrotinas INPUT e SOLVAB, descritas para ALGO com base no algoritmo de
Gear, também são utilizadas pelo PALGO.
2.4. Consistência de inicialização 24
2.4 Consistência de inicialização
As condições iniciais de um sistema de EADs na sua forma geral:
F(x,
˙
x, y, t) = 0 (2.82)
são definidas pelo vetor (x
0
,
˙
x
0
, y
0
) em t
0
, e não por (x
0
, y
0
).
Uma condição necessária mas não suficiente para que estas condições iniciais sejam
consistentes é que elas satisfaçam ao sistema:
F(x
0
,
˙
x
0
, y
0
, t
0
) = 0 (2.83)
Por exemplo, é necessário que as condições iniciais (x
1
0
, x
2
0
, ˙x
1
0
, ˙x
2
0
) satisfaçam
as seguintes EADs de índice 1 (PANTELIDES, 1988):
˙x
1
+ ˙x
2
= a(t)
x
1
+ x
2
2
= b(t)
(2.84)
onde a(t) e b(t) são funções contínuas e diferenciáveis.
Entretanto, para que exista a consistência, elas devem satisfazer também à equa-
ção resultante da diferenciação da segunda equação:
˙x
1
+ 2x
2
˙x
2
=
˙
b(t) (2.85)
Portanto, mesmo os sistemas de índice 1, podem apresentar problemas de inicia-
lização, que são inerentes a todos os sistemas de índice superior. Em outras palavras,
enquanto a inicialização da maioria dos sistemas de índice 1 é similar à das EDOs,
todos os sistemas de índice superior e alguns sistemas de índice 1 apresentam dificul-
dades, como as mostradas. Além disto, é preciso localizar e evitar as diferenciações de
equações que não trazem informações adicionais. O que se busca, portanto, são for-
mas de extrair dos sistemas informações fundamentais que estão nele contidas, mas que
não são aparentes. Os atuais códigos numéricos de solução de EADs freqüentemente
falham ou se tornam extrem amente ineficientes se os valores iniciais são inconsistentes.
2.5. Função Identificadora de Fases (FIF) 25
2.5 Função Identificadora de Fases (FIF)
Funções Identificadoras de Fases (Switching Functions) são funções que indicam quando
uma restrição que está ativa, torna-se inativa e vice-versa.
Um caso particular de grande interesse é quando a variável de controle aparece
linearmente na função Hamiltoniano. Em geral, nenhum mínimo existirá para tais
problemas a menos que restrições de desigualdade nas variáveis de estado e/ou de
controle sejam especificadas. Se as restrições de desigualdade são lineares na variável
de controle, é razoável esperar que a solução mínima, se existir, sempre exigirá que a
variável de controle esteja localizada em um ponto ou outro do limite da região viável
de controle (BRYSON; HO, 1975).
Seja o seguinte sistema de equações
˙
x = F(x) + g(x)u (2.86)
x(0) = x
o
(2.87)
com variável de controle escalar dado por
u
min
≤ u ≤ u
max
(2.88)
A função Hamiltoniano é definida como
H = λ
T
(F (x) + g(x)u) (2.89)
Para essa classe de controle tem-se
u(t) =
u
m´ax
se λ
T
g < 0
Função Identificadora de Fase se λ
T
g = 0
u
m
´
in
se λ
T
g > 0
(2.90)
Seywald et al. (1994) obtiveram trajetórias ótimas de um foguete em um plano
vertical aplicando o PMP. Para obtenção destes resultados eles formularam FIFs associ-
adas a doze possíveis lógicas de controle baseadas na ativação/desativação de restrições
e nos sinais das variáveis adjuntas e que necessariamente tivessem significado do ponto
de vista da engenharia e não somente significado matemático.
Santos et al. (2005) analisaram a FIF para um PCO da fermentação alcoólica em
reator batelada alimentada para maximação do produto com restrições na produti-
vidade, na variável de controle e no volume do fermentador e operação com tempo
livre e que apresenta arcos singulares. A estratégia de controle é definida em 5 fases,
sendo 2 com controle singular e 3 com controle bang-bang, baseada em aspectos físicos
2.6. Métodos Numéricos de solução de PCOADs 26
do problema. O índice diferencial flutua entre 1 e 3 ao longo destas fases e o PCO é
resolvido por fase como um sistema de EDOs de valor inicial.
2.6 Métodos Numéricos de solução de PCOADs
2.6.1 Métodos Diretos
Com o desenvolvimento dos algoritmos para a solução dos problemas de Programação
Não Linear (NLP) (BIEGLER, 1984; RENFRO et al., 1987; CERVANTES; BIEGLER, 1998;
BIEGLER et al., 2002), a abordagem direta para a solução dos PCOs é a mais utilizada
atualmente e a que tem recebido maior atenção dos pesquisadores. Estes métodos
transformam o PCO num NLP através de parametrização das variáveis de controle
e/ou de estado. Eles são divididos em seqüencial e simultâneo. A estratégia seqüencial
consiste na parametrização da variável de controle. Neste caso as condições iniciais e o
conjunto de parâmetros de controle são conhecidos e o sistema de EADs é discretizado
utilizando-se uma aproximação polinomial e em seguida é resolvido por um NLP. Esse
procedimento determina o valor da função objetivo e das restrições encontrando valores
ótimos de coeficientes na parametrização do controle. A estratégia simultânea consiste
na parametrização tanto da variável de controle quanto da variável de estado. Como
resultado dessa parametrização as EADs são resolvidas somente uma vez, no ponto
ótimo, o que evita soluções intermediárias que não existem ou esforço computacional
excessivo.
Na literatura podem ser encontradas várias formas de discretização baseadas em
aproximações por colocação ortogonal (OH; LUUS, 1977; BIEGLER, 1984), por expan-
sões em séries (ÇELIK et al., 2003) e por colocação ortogonal em elementos finitos
(RENFRO, 1986; RENFRO et al., 1987; CUTHREL; BIEGLER, 1987b, 1987a; LOGSDON;
BIEGLER, 1989, 1992). Em particular, no uso de aproximações em elementos finitos, o
tamanho e número de elementos pode ser determinado automaticamente, permitindo
um melhor controle do erro de discretização (LOGSDON; BIEGLER, 1989; VASANTHA-
RAJAN; BIEGLER, 1990). As técnicas de discretização têm natureza local (métodos
das diferenças finitas e volumes finitos) ou global (colocação ortogonal padrão e em
elementos finitos). As técnicas de natureza local usualmente requerem um grande
número de pontos de discretização para que a s olução final obtida aproxime bem a
solução real procurada. Isto, em geral, significa que sistemas de dezenas ou centenas
de equações têm que ser resolvidos simultaneamente. O processo de geração das equa-
ções, no entanto, é simples e pode ser automatizado com certa facilidade. Por sua
vez, as técnicas de natureza global podem permitir considerável redução do tamanho
da malha de discretização, embora o processo de geração das equações discretizadas
seja significativamente mais complexo do que o caso anterior. Os Métodos Diretos
podem fornecer soluções sub-ótimas e resultados menos precisos quando comparados
2.6. Métodos Numéricos de solução de PCOADs 27
aos Méto dos Indiretos. Possui como vantagem a solução rápida e como desvantagem
a dificuldade na determinação da otimalidade.
Um código que resolve PCO através da abordagem do método direto é o código
DIRCOL, que é um conjunto de subrotinas desenvolvido por Stryk (1999) e imple-
mentado utilizando o Fortran com o objetivo de resolver PCO descr ito por equa-
ções diferenciais de primeira ordem sujeito a restrições de igualdade e/ou desigualdade
nas variáveis de controle e/ou variáveis de estado. A DIRCOL utiliza um método de
colocação direto que discretiza as variáveis de estado e de controle por aproximações
polinomiais por partes transformando o PCO num NLP, que é resolvido através de
Programação Quadrática Seqüencial densa (SQP) pelo código NPSOL ou pelo código
SNOPT apropriado para sistemas esparsos (GILL et al., 1997, 1986). O código informa
estimativas das variáveis adjuntas e dos eventos, resolve problemas com multi-fases,
através do particionamento do sistema com estimativa dos eventos em cada fase e
possui estratégia de refinamento da malha.
Cuthrel e Biegler (1987b) resolveram PCOs sujeito a restrições de desigualdade
nos perfis de estado e de controle e com dependência linear no controle, que causa
arcos singulares ou perfis liga-desliga. Nestes arcos, o perfil de controle não influen-
cia diretamente as condiçõ es de otimalidade do Hamiltoniano e a determinação do
controle requer condições adicionais, que podem ser de difícil manuseio. Os autores
incorporaram as equações algébricas com coeficientes desconhecidos, referentes aos re-
síduos oriundos da aplicação da colocação ortogonal em elementos finitos utilizando
os polinômios de Lagrange como base, como restriçõ es de igualdade ao NLP e os co-
eficientes como variáveis de decisão. Foram utilizadas estratégias de minimização do
erro de aproximação além dos chamados superelementos, que são um nível extra de
elementos utilizados na alocação adaptativa dos elementos.
Feehery e Barton (1998) utilizaram a parametrização no controle em PCOs sujei-
tos a restriçõ es nas variáveis de estado. Os autores mostram que as EADs se tornam
de índice superior durante as partes restritas pelas variáveis de estado (flutuação de
índice durante a trajetória) e estabelecem a equivalência entre o PCO com restrição de
desigualdade e um problema de otimização híbrido (discreto/contínuo) que contém o
fenômeno de mudança de fase. O algoritmo desenvolvido detect a a ativação e a desa-
tivação das restrições durante a solução do Problema de Valor Inicial (PVI) e resolve
as EADs de índice superior resultante e as equações de sensibilidade usando o método
das Derivadas Mudas.
2.6.2 Métodos Indiretos
Os métodos indiretos surgiram com o desenvolvimento do Cálculo Variacional, permi-
tindo a dedução das condições necessárias e suficientes para a solução de problemas de
otimização dinâmica (BRYSON; HO, 1975). Os trabalhos de Bryson e Ho (1975), Lynn
et al. (1970), Lynn e Zahradnik (1987), aplicaram o princípio de Pontryagin com a ge-
2.6. Métodos Numéricos de solução de PCOADs 28
ração das equações diferenciais adjuntas e a condição necessária para a minimização da
função Hamiltoniano, transformando o problema original em um problema de valor no
contorno em dois pontos (TPBV). Estes TPBV podem ser resolvidos por métodos de
discretização como colocação em elementos finitos e diferenças finitas, “chute"simples
e múltiplo.
Atualmente, os métodos indiretos podem ser utilizados de modo eficiente devido
ao desenvolvimento dos programas de álgebra computacional, que permite a obten-
ção automática das equações diferenciais adjuntas e demais condições de otimalidade.
Desta forma, podem ser obtidas soluções altamente precisas e testes feitos a posteriori
podem excluir as soluções subótimas (CHUDEJ; GÜNTHER, 1999).
As condições necessárias para o ótimo para problemas nos quais o sistema dinâ-
mico é descrito somente através de equações diferenciais ordinárias são bem estabe-
lecidas na literatura (BRYSON; HO, 1975). As condições generalizadas para sistemas
descritos por EADs, conforme demonstrado por Renfro (1986) e Feehery (1998) fo-
ram apresentadas na Seção 2.2. Os Métodos Indiretos têm uma faixa de convergência
restrita e apresentam dificuldades de convergência, apesar de garantir o atendimento
às restrições do problema. Um código que aborda o Método Indireto é o COLDAE
(ASCHER; SPITERI, 1994) que é aplicado a EADs de valor no contorno não lineares
semi-explícitas de índice no máximo dois e totalmente implícitos de índice um. A
COLDAE é uma subrotina geral para a solução de EDOs de valor no contorno. O
método implementado é a colo cação polinomial por partes em pontos gaussianos, aco-
plado a um método de projeção de Ascher-Petzold.
Lynn et al. (1970) e Lynn e Zahradnik (1987) aplicaram o PMP com a geração das
equações adjuntas e da condição necessária para a minimização da função Hamiltoni-
ano. No primeiro trabalho, o problema de valor no contorno puramente diferencial,
resultante da otimização da conversão de reaçõ es consecutivas em um reator tubular
com dispersão axial, foi resolvido por um método de resíduos ponderados. Esta técnica,
batizada como aproximação da trajetória para controle quase ótimo, foi estendida para
sistemas distribuídos r epresentados por equações diferenciais parciais hiperbólicas de
primeira ordem no segundo trabalho.
San e Stephanopoulos (1984) estudaram a cinética de crescimento da biomassa
e a estratégia de alimentação para a otimização da biomassa final produzida em um
reator batelada. Para esse caso particular a ativação dos arcos singulares foi detectada,
determinando assim as condições ao longo e na trajetória final do arco singular. Os
mesmos autores San e Stephanopoulos (1989) mostraram que a escolha da variável de
controle é de fundamental importância para garantir a otimalidade e a viabilidade da
implementação. Eles encontraram soluções sub-ótimas para a operação do reator com
taxa de alimentação e concentração de substrato fixa como uma variável de controle
devido à inabilidade do controlador manter baixa concentração de substrato no fim
da fase de produção. Por outro lado, o uso da concentração de substrato como uma
variável de controle pela manipulação da concentração mostrou ser uma alternativa
efetiva. Eles propuseram uma solução em dois passos, que determina primeiro a con-
2.6. Métodos Numéricos de solução de PCOADs 29
centração ótima de substrato no fermentador e depois calcula o perfil de concentração
de substrato na alimentação que resulta no perfil ótimo do fermentador.
Modak et al. (1986) apresentaram a formulação do problema de otimização do
processo de fermentação batelada alimentada e a obtenção e aplicação das condições
de ótimo, considerando que o controle singular é factível conforme as características
das funções taxas expecíficas de crescimento celular e de formação de produto. Fo-
ram consideradas três combinações diferentes de taxas representadas por funções não
monotônicas que exibem máximos. No tipo mais comum, o crescimento celular não
sofre inibição ou repressão e tem o comportamento expresso pela Equação de Mono d
enquanto a formação de produto sofre inibição. Em todos os casos, considerou-se que
a vazão da corrente alimentada pode alcançar o seu valor máximo, mínimo ou singular
de modo a permitir inicialmente o crescimento celular e na sequência fazer com que
as células produzam o produto desejado de forma ótima. Além da dependência das
funções das taxas, as seqüências ótimas são influenciadas pelas condições iniciais as-
sociadas. Uma vez identificadas as seqüências, o problema de otimização é reduzido a
um problema de determinações dos eventos, do tempo final e da taxa de alimentação
singular.
Gomes (2000) estudou a abordagem algébrico-diferencial na solução de PCOs
através do método indireto, utilizando o código COLDAE, e propôs um código que
gera as equações necessárias para o ótimo a partir da aplicação o Princípio do Mínimo
de Pontryagin. Est e código é chamado OTIMA que será abordado no Capítulo 3.
2.6.3 Métodos Mistos ou Híbridos
Os Métodos Mistos ou Híbridos são uma combinação de métodos diretos e métodos
indiretos. Nesse método, os Métodos Diretos são aplicados a problemas m ais simplifica-
dos e os resultados servem de estimativas para os Métodos Indiretos, com refinamento
da solução.
Stryk e Schlemmer (1994) resolveram um problema de minimização de energia
e tempo de um robô industrial utilizando uma combinação entre o m étodo direto e
o indireto. O método direto é aplicado com estimativas inici ais sub-ótimas da solu-
ção x(t), u(t). A solução sub-ótima fornece estimativas confiáveis das variáveis de
estado e variáveis adjuntas e da função identificadora de fases das restrições no es-
tado e no controle. Com essas estimativas, o método dos chutes múltiplos é aplicado
para um problema de valor no contorno em múltiplos pontos resultante das condições
necessárias para a otimalidade.
Chudej e Günther (1999) desenvolveram uma ab ordagem baseada no espaço de
estados, que transforma o PCO com restrições no estado num PCO com restrições no
controle, aplicando-a na solução de um problema de otimização da trajetória de um
avião supersônico num plano vertical sujeito a restrições, que possui quatro variáveis de
2.6. Métodos Numéricos de solução de PCOADs 30
estado e 2 variáveis de controle. Desta forma, o sistema de EADs de valor no contorno
de índice superior é transformado num sistema de EDOs de valor no contorno com
atendimento automático das restrições no estado através da eliminação de um dos
componentes do vetor estado no arco do contorno através da utilização da restrição
no estado, levando a uma manipulação analítica considerável das equações. A partir
de uma FIF pré-definida, os autores formulam um PCO substituto sujeito a condições
definidas em pontos interiores e nos contornos, para cada uma das fases identificadas
pelo sinal da restrição de estado. O PCO substituto é resolvido através de um método
híbrido onde estimativas iniciais dos estados e das variáveis adjuntas obtidas através
do código DIRCOL são utilizadas no algoritmo de chute múltiplo Mumus.
Oberle e Sothmann (1999) resolveram um PCO de tempo final livre, relativo
a um processo de fermentação batelada alimentada que representa a biossíntese de
penicilina, utilizando parametrização na variável de controle e solução do NLP por um
método SQP padrão, no qual admite-se a continuidade do controle. Aplicando a teoria
de controle ótimo e resolvendo o PVC resultante associado à função identificadora de
fase (switching function) por técnicas de chute múltiplo associadas a relaxação dos
parâmetros característicos das equações de taxa e métodos de continuação, os autores
verificaram que o problema resolvido pelo método indireto satisfez a todas as condições
necessárias de primeira e segunda ordens estabelecidas pela teoria, apresentando uma
solução similar à obtida pela abordagem direta mas com pequenos sub-arcos bang-bang
no controle.
Lobato (2004) resolveu PCOs inicialmente pelo Dircol, buscando a melhor solu-
ção possível fornecida pelo código, que fornece também a estimativa dos Eventos e do
perfil das variáveis adjuntas. Com a estimativa dos Eventos e análise dos res ultados
obtidos, são construídas as Funções Identificadoras de Fases. Os Sistemas Aumenta-
dos obtidos simbolicamente, relativos a cada fase e aos quais foi aplicado o Método
da Variável de Folga, no caso da presença de restrições de desigualdade, são então
caracterizados pelos códigos estruturais e as reduções de índices por eles indicadas são
executadas simbolicamente. Os problemas de valor no contorno de índice 1, referentes
a cada fase, são resolvidos seqüencialmente pelo código Coldae. Neste trabalho, não
é utilizada a formulação multifásica como entrada do código DIRCOL, ou seja, o pro-
blema é resolvido com apenas uma fase. O autor demonstra que este método possui
aplicabilidade a problemas com restrições nas variáveis de estado e com controle linear
nas restrições, que provo cam flutuações no índice do PCOADs.
2.7. Sistemas Híbridos 31
2.7 Sistemas Híbridos
2.7.1 Definição
Os Sistemas Híbridos consistem em sistemas dinâmicos onde os comportamentos con-
tínuos e discretos interagem entre si. A investigação de sistemas híbridos está criando
uma interface nova e fascinante entre a engenharia de controle, a matemática e a infor-
mática. Eles estão presentes na análise e no projeto de sistemas de controle inteligentes
e tem se disseminado com o uso difundido de máquinas digitais (XU, 2001).
Um exemplo comum é um sistema de aquecimento de um quarto. Este sistema
possui dois m odos de operação: aquecimento ligado (ON) ou desligado (OFF). As
dinâmicas contínuas associadas a esses dois modos são diferentes. Se é detectado que
a temperatura do quarto está abaixo de um determinado nível, o aquecedor é ligado e
a dinâmica de aquecimento fica ativa. Por outro lado, se a temperatura do quarto está
em um nível acima do desejado, o aquecedor é desligado e outra dinâmica contínua
fica ativa (XU, 2001).
Outro exemplo interessante é o de uma máquina do Lehrstuhl B für Mechanik da
Technische Universität München provida de 6 “pernas", para a qual o problema consiste
em determinar simultaneamente os controles contínuos ótimos para as “pernas"e o
modo ótimo de caminhar, definido por um número e seqüência discretos de passos, de
posições e tipos de contatos com o solo para cada perna, para minimizar o consumo
total de energia da máquina ou o tempo necessário para a caminhada. A estrutura das
equações do movimento para esta máquina depende da situação de contato da perna
com o solo (variáveis discretas), que pode ser de 3 tipos: nenhum contato, contato fixo e
contato deslizante. As mudanças na estrutura das equações, relacionadas às mudanças
na situação de contato, resultam em descontinuidades nas equações de movimento nos
tempos intermediários (eventos), que separam o problema em fases com dinâmicas
específicas (STRYK, 2000).
Segundo Stryk (2000), um Problema de Controle Ótimo Híbrido é definido como
a minimização de um índice de desempenho híbrido dado por:
J[u, v] =
n
c+1
i=1
ϕ
(i)
(x(t
c
i
− 0), x(t
c
i
+ 0), q(t
c
i
− 0), q(t
c
i
+ 0), t
c
i
)+
+
n
c+1
i=1
t
c
i
t
c
i−1
L
(i)
(x(t), u(t), q(t), v(t), t)dt (2.91)
2.7. Sistemas Híbridos 32
Figura 2.2: Dinâmica do estado contínuo definido em n
c+1
fases.
sujeito a um sistema dinâmico não-linear
˙x(t) = f
(i)
(x(t), u(t), q(t), v(t), t), t
c
i−1
< t < t
c
i
, i = 1, ..., n
c+1
, (2.92)
que descreve o estado dinâmico contínuo em n
c
+ 1 fases t(
c
i
), t(
c
i+1
), i = 0, ..., n
c
(Fig.2.2). As variáveis especificadas denotam:
x : [0, t
f
] → X ⊂ R
n
x
estado contínuo, q : [0, t
f
] → Q ⊂ Z
n
q
estado discreto,
u : [0, t
f
] → U ⊂ R
n
u
controle contínuo, v : [0, t
f
] → V ⊂ Z
n
v
controle discreto.
A transição de uma fase i para uma fase i+1 ocorre em um momento desconhecido
t(
c
i
), i = 0, ..., n
c
, isto é, nos Eventos.
0 = t
c
0
< t
c
1
< ... < t
c
n
c
< t
c
n
c+1
= t
f
(2.93)
Aqui, x é continuamente diferenciável por partes e u é uma função contínua
por partes do tempo. Freqüentemente o número total de eventos n
c
é desconhecido
inicialmente e deve ser determinado.
As trajetórias de estado e controle são comumente sujeitas a restrições de igual-
dade e desigualdade
0 ≤ g
(i),j
(x(t), u(t), q(t), v(t), t), j = 1, ..., n
g(i)
, (2.94)
0 = h
(i),j
(x(t), u(t), q(t), v(t), t), j = 1, ..., n
h(i)
. (2.95)
O número e o tipo de restrições podem diferir de uma fase para outra e depender
2.7. Sistemas Híbridos 33
dos valores atuais de q e v.
As condições implícitas devem satisfazer nos tempos iniciais, finais e nos eventos
0 = r
(0),j
(x(0), q(0), x(t
f
), q(t
f
), t
f
), j = 1, ..., n
r(0)
, (2.96)
0 = r
(i),j
(x(t
c
i
− 0), q(t
c
i
− 0), x(t
c
i
+ 0), q(t
c
i
+ 0), t
c
i
), j = 1, ..., n
r(i)
, i = 1, ..., n
c
. (2.97)
Os eventos devem ser definidos nos zeros das funções identificadoras de fases r
(i),j
.
Condições explícitas no tempo inicial e final são adicionadas
x
j
1
(0) = x
0,j
1
, q
j
2
(0) = q
0,j
2
,
x
k
1
(t
f
) = x
f,k
1
, q
k
2
(t
f
) = q
f,k
2
,
(2.98)
assim como condições de transição de fases explícitas
x
l
1
(t
c
i
+ 0) = R
q
(i),l
1
(x(t
c
i
− 0), q(t
c
i
− 0), t
c
i
),
q
l
2
(t
c
i
+ 0) = R
q
(i),l
2
(x(t
c
i
− 0), q(t
c
i
− 0), t
c
i
),
(2.99)
podem ser adotadas. Aqui, j
1
,k
1
,l
1
são elementos dos sub-conjuntos de 1, 2, ..., n
x
e
j
2
,k
2
,l
2
são elementos dos sub-conjuntos de 1, 2, ..., n
q
. As constantes reais x
0,j
1
,x
f,k
1
e
as constantes inteiras q
0,j
2
,q
f,k
2
são dadas. As funções que aparecem na definição do
problema de controle ótimo híbrido são, ao menos localmente, contínuos e continua-
mente diferenciável em relação aos argumentos contínuos x, u e t.
2.7.2 Sistemas Chaveados
Um Sistema Chaveado é um tipo particular de Sistema Híbrido, mas cujo est ado
contínuo não apresenta saltos nos eventos. Se uma lei de mudança for especificada,
então o Sistema Chaveado pode ser descrito como (XU, 2001):
˙x(t) = f
i(t)
(x(t), u(t), t) (2.100)
i(t) = j(x(t), i(t
−
), t) (2.101)
onde i : R → I é uma função constante por partes que determina o subsistema ativo
em cada instante t.
Para um sistema chaveado S, a seqüência de mudança σ em [t
0
, t
f
] é definida como
σ = ((t
0
, i
0
), (t
1
, e
1
), (t
2
, e
2
)...(t
K
, e
K
)), com 0 ≤ K < ∞, t
0
≤ t
1
≤ t
2
≤ ... ≤ t
K
≤ t
f
,
e i
0
∈ I, e
k
= (i
k−1
, i
k
) ∈ E para k = 1, 2, ..., K.
2.7. Sistemas Híbridos 34
A seqüência de mudança σ definida anteriormente indica que o subsistema i
k
está
ativo em [t
k
, t
k+1
] se t
k
< t
k+1
([t
K
, t
f
] se k = K), e i
k
passa para o instante t
k
se
t
k
= t
k+1
, ou seja, o sistema passa do subsistema i
k−1
para o subsistema i
k
e então do
subsistema i
k+1
no instante t
k
.
Bemporad et al. (2002) resolveram problemas de controle ótimo para sistemas
chaveados onde os eventos e a seqüência de mudança são variáveis de decisão. A solução
destes problemas é feita através de dois algoritmos que alternam entre si. O algoritmo
Master que encontra a seqüência ótima de mudança e o algoritmo Slave que encontra os
eventos ótimos. De acordo com os autores, para melhorar o desempenho do algoritmo
uma abordagem heurística pode ser adotada. Egerstedt et al. (2003) propõem uma
nova forma do gradiente do custo funcional descrita por Xu (2001) para facilitar o
uso de algoritmos baseados no método do gradiente descendente (gradient descent).
Esta nova estrutura do gradiente permite que o número de eventos torne-se parte do
controle em vez de ser uma constante dada. Xu e Antsaklis (2004) transcreveram um
sistema chaveado em um PCO equivalente onde os eventos desconhecidos no problema
original fazem parte do estado contínuo. O problema equival ente tem a propriedade
de que os eventos são fixos em relação a nova variável independente do tempo.
A metodologia descr ita por Xu e Antsaklis (2004) para sistemas chaveados defi-
nidos pelas Eqs.(2.102 e 2.103) e aplicada a um sistema com apenas um evento para
simplicidade de notação é descrita em detalhes na seqüência. Considere um sistema
chaveado dado por:
˙x(t) = f
1
(x, u, t), t
0
≤ t < t
1
(2.102)
˙x(t) = f
2
(x, u, t), t
1
≤ t ≤ t
f
(2.103)
O objetivo é encontrar um evento t
1
e uma entrada continua ótima u(t), t ∈ [t
0
, t
f
]
tal que o custo funcional dado pela Eq.(2.104) seja minimizado, dados t
0
, t
f
e x(t
0
) =
x
0
.
J = Ψ(x(t
f
)) +
t
f
t
0
L(x, u, t)dt (2.104)
É introduzida uma nova variável de estado t
s1
correspondente ao evento t
1
, tal
que:
∂t
s1
∂t
= 0 (2.105)
t
s1
(0) = t
1
(2.106)
Uma variável independente do tempo τ é também introduzida e uma relação linear
2.7. Sistemas Híbridos 35
por partes entre t e τ é estabelecida como se segue:
t =
t
0
+ (t
s1
− t
0
)τ, 0 ≤ τ < 1
t
s1
+ (t
f
− t
s1
)(τ − 1), 1 ≤ τ ≤ 2
(2.107)
Claramente, pode-se perceber que τ = 0 corresponde a t = t
0
, τ = 1 a t = t
1
e
τ = 2 a t = t
f
. Introduzindo t
s1
e τ, o problema original dado pelas Eqs.(2.102) e
(2.103) pode ser transcrito no seguinte problema equivalente definido por fases:
dx(τ)
dτ
= (t
s1
− t
0
)f
1
(x, u) (2.108)
dt
s1
dτ
= 0 (2.109)
no intervalo τ ∈ [0, 1) e
dx(τ)
dτ
= (t
f
− t
s1
)f
2
(x, u) (2.110)
dt
s1
dτ
= 0 (2.111)
no intervalo τ ∈ [1, 2]. O problema é achar t
s1
e u(τ) ótimos, τ ∈ [0, 2] tal que o custo
funcional:
J = Ψ(x(2)) +
1
0
(t
s1
− t
0
)L(x, u)dτ +
2
1
(t
f
− t
s1
)L(x, u)dτ (2.112)
seja minimizado. São dados t
0
, t
f
e x(0) = x
0
.
Definindo:
˜
f
1
(x, u, t
s1
)
∆
= (t
s1
− t
0
)f
1
(x, u) (2.113)
˜
f
2
(x, u, t
s1
)
∆
= (t
f
− t
s1
)f
2
(x, u) (2.114)
˜
L
1
(x, u, t
s1
)
∆
= (t
s1
− t
0
)L(x, u) (2.115)
˜
L
2
(x, u, t
s1
)
∆
= (t
f
− t
s1
)L(x, u) (2.116)
e definindo a função Hamiltoniana parametrizada como:
H(x, p, u, t
s1
)
∆
=
˜
L
1
(x, u, t
s1
) + p
T
˜
f
1
(x, u, t
s1
), τ ∈ [0, 1)
˜
L
2
(x, u, t
s1
) + p
T
˜
f
2
(x, u, t
s1
), τ ∈ [1, 2]
(2.117)
2.7. Sistemas Híbridos 36
Aplicando as condições necessárias, as equações de estado, co-estado e de estaci-
onariedade e as condições de contorno e continuidade, são fornecidas:
∂x
∂τ
=
∂H
∂p
T
=
˜
f
k
(x, u, t
s1
) (2.118)
∂p
∂τ
= −
∂H
∂x
T
= −
∂
˜
f
k
∂x
T
p −
∂
˜
L
k
∂x
T
(2.119)
0 =
∂H
∂u
T
=
∂
˜
f
k
∂u
T
p −
∂
˜
L
k
∂u
T
(2.120)
x(0, t
s1
) = x
0
(2.121)
p(2, t
s1
) =
∂Ψ
∂x
(x(2, t
s1
))
T
(2.122)
p(1−, t
s1
) = p(1+, t
s1
) (2.123)
As Eqs.(2.118) - (2.123) formam um sistema de EADs de valor no contorno para-
metrizado pelo evento. Assim, o valor ótimo de J
1
como uma função dos parâmetros
t
s1
é dado pela Eq.(2.124).
J
1
(t
s1
) = ψ(x(2, t
s1
)) +
1
0
˜
L
1
(x, u, t
s1
)dτ +
1
0
˜
L
2
(x, u, t
s1
)dτ (2.124)
Diferenciando J
1
em relação à t
s1
fornece a seguinte expressão:
∂J
1
∂t
s1
=
∂ψ(x(2, t
s1
))
∂x
∂x(2, t
s1
)
∂t
s1
+
1
0
L(x, u) + (t
s1
− t
0
)
∂L
∂x
∂x
∂t
s1
+
∂L
∂u
∂u
∂t
s1
dτ+
+
2
1
−L(x, u) + (t
f
− t
s1
)
∂L
∂x
∂x
∂t
s1
+
∂L
∂u
∂u
∂t
s1
dτ (2.125)
É necessário obter as funções
∂
∂τ
∂x
∂t
s1
e
∂
∂τ
∂u
∂t
s1
para se obter o valor
dJ
1
dt
s1
.
Diferenciando as equações de estado, co-estado e estacionariedade em re lação a t
s1
,
obtém-se:
∂
∂τ
∂x
∂t
s1
=
∂
∂t
s1
∂x
∂τ
= f
1
+ (t
s1
− t
0
)
∂f
1
∂x
∂x
∂t
s1
+
∂f
1
∂u
∂u
∂t
s1
(2.126)
2.7. Sistemas Híbridos 37
∂
∂τ
∂p
∂t
s1
= −
∂
∂t
s1
∂p
∂τ
= −
∂f
1
∂x
T
p −
∂L
∂x
T
− (t
s1
− t
0
) (2.127)
×
∂f
1
∂x
T
∂p
∂t
s1
+
p
T
∂
2
f
1
∂x
2
∂x
∂t
s1
T
+
p
T
∂
2
f
1
∂x∂u
∂u
∂t
s1
T
+
∂
2
L
∂x
2
∂x
∂t
s1
+
∂
2
L
∂x∂u
∂u
∂t
s1
0 =
∂f
1
∂u
T
p −
∂L
∂u
T
− (t
s1
− t
0
) (2.128)
×
∂f
1
∂u
T
∂p
∂t
s1
+
p
T
∂
2
f
1
∂u∂x
∂x
∂t
s1
T
+
p
T
∂
2
f
1
∂u
2
∂u
∂t
s1
T
+
∂
2
L
∂u∂x
∂x
∂t
s1
+
∂
2
L
∂u
2
∂u
∂t
s1
para τ ∈ [0, 1) e
∂
∂τ
∂x
∂t
s1
=
∂
∂t
s1
∂x
∂τ
= −f
2
+ (t
f
− t
s1
)
∂f
2
∂x
∂x
∂t
s1
+
∂f
2
∂u
∂u
∂t
s1
(2.129)
∂
∂τ
∂p
∂t
s1
= −
∂
∂t
s1
∂p
∂τ
=
∂f
2
∂x
T
p +
∂L
∂x
T
− (t
f
− t
s1
) (2.130)
×
∂f
2
∂x
T
∂p
∂t
s1
+
p
T
∂
2
f
2
∂x
2
∂x
∂t
s1
T
+
p
T
∂
2
f
2
∂x∂u
∂u
∂t
s1
T
+
∂
2
L
∂x
2
∂x
∂t
s1
+
∂
2
L
∂x∂u
∂u
∂t
s1
0 = −
∂f
2
∂u
T
p −
∂L
∂u
T
+ (t
f
− t
s1
) (2.131)
×
∂f
2
∂u
T
∂p
∂t
s1
+
p
T
∂
2
f
2
∂u∂x
∂x
∂t
s1
T
+
p
T
∂
2
f
2
∂u
2
∂u
∂t
s1
T
+
∂
2
L
∂u∂x
∂x
∂t
s1
+
∂
2
L
∂u
2
∂u
∂t
s1
para τ ∈ [1, 2].
Diferenciando as condições de contorno e de continuidade em relação à t
s1
:
dx(0, t
s1
)
dt
s1
= 0 (2.132)
dp(2, t
s1
)
dt
s1
=
∂
2
Ψ(x(2, t
s1
))
∂x
2
∂x(2, t
s1
)
∂t
s1
(2.133)
dp(1−, t
s1
)
dt
s1
=
dp(1+, t
s1
)
dt
s1
(2.134)
2.7. Sistemas Híbridos 38
O sistema de equações formado pelas Eqs.(2.118) - (2.123) e pelas Eqs.(2.126) -
(2.131) e as condições dadas pelas Eqs.(2.132) - (2.134), forma um problema algébrico-
diferencial de valor no contorno em dois pontos (BVPDAE) parametrizado pelo evento
t
s1
. Resolvendo-as e substituindo o valor na função objetivo J
1
, pode-se obter
dJ
1
dt
s1
.
Para cada evento, a EAD pode ser resolvida através de métodos numéricos.
2.7.3 Analogia entre PCOADs e Sistemas Chaveados
Segundo Feehery (1998), a solução de PCO com restrição de desigualdade apresenta
desafios especiais, porque demanda o conhecimento da seqüência e do número de ati-
vações e desativações ao longo da trajetória. Segundo o autor existe uma equivalência
entre PCOAD com restrições de desigualdade nas variáveis de estado e os PCOs com
comportamento contínuo e discreto acoplados, chamados Sistemas Híbridos.
Os com portamentos contínuos e discretos interagem em p ontos isolados no tempo,
conhecidos como Eventos. Nesses pontos, devido à descontinuidade no estado e/ou nas
variáveis adjuntas (saltos), podem ocorrer mudanças na forma funcional das EADs
e/ou nas trajetórias da variável de controle em cada fase. Como conseqüência, nes-
ses problemas, o índice diferencial do sistema pode flutuar ao longo da trajetória de
solução, exigindo estratégias numéricas específicas aplicadas em cada fase. Por isso,
é necessário o conhecimento prévio da história da ativação e desativação das restri-
ções, feito através das chamadas Funções Identificadoras de Fases (FIF), para resolver
os problemas adequadamente. Para formular a FIF, Srinivasan et al. (2003), Lobato
(2004) e Stryk (2000) propõem a solução de um PCO sem restrições inicialmente, de-
pois a obtenção da solução por um método direto que estime as mudanças de fases e
os eventos.
Nesta dissertação, as variações discretas na forma funcional das EADs como uma
conseqüência dos Eventos que ocorrem instantaneamente num determinado tempo são
pré-definidas através da formulação prévia das FIF enquanto o tempo de ocorrência
dos Eventos é definido implicitamente pelo estado do sistema, utilizando o algoritmo
proposto por Xu e Antsaklis (2004) descrito anteriormente. Desta forma, enquanto
em Lobato (2004) a identificação dos Eventos era feita através do acompanhamento
da ativação das restrições ao longo da solução e a solução pelo método direto era feita
fase a fase seqüencialmente, neste trabalho a identificação de eventos é feita numerica-
mente, utilizando a abordagem algébrico-diferencial. A formulação do problema com
inclusão das equações de sensibilidade é feita simbolicamente, através de algoritmo
especialmente desenvolvido. A solução pelo método direto usando o código DIRCOL
através da formulação multifásica foi realizada, com estimativa dos Eventos.
CAPÍTULO 3
Desenvolvimento de uma interface de
integração de ferramentas para o
tratamento de PCOADs
3.1 Aspectos Gerais
A necessidade de se caracterizar sistemas de EA Ds e reduzir o índice de sistemas de
índice superior, ou analisar a consistência de condições in iciais, motivaram a atualiza-
ção dos códigos INDEX e ACIG, além da atualização do código que gera as equações
adjuntas e condição de ótimo através do Princípio de Pontryagin (OTIMA) que facilita
a utilização de métodos indiretos de solução de PCOADs. Para ampliar a utilização
destes códigos à outros usuários do núcleo, não os tornando obsoletos e esquecidos e
melhorar o estudo de PCOADs, foi desenvolvida uma interface amigável para o usuário
que exija pouco conhecimento da sintaxe Maple. Esta união de códigos em uma inter-
face foi denominada OpCol. Uniu-se a esta interface o código EVENTS que através
do algoritmo desenvolvido por Xu e Antsaklis (2004) faz a parametrização de sistemas
chaveados em relação aos eventos.
3.2. Pacote Maplets 40
Neste capítulo são descritas as características da interface criada, para a atualiza-
ção, extensão e agrupamento de ferramentas de manipulação, tratamento e geração de
condições de otimalidade para PCOADs usando o Maple
R
9.5 (MAPLESOFT, 2003).
Neste desenvolvimento, foi necessário fazer a atualização de funções e comandos para
o Maple 9.5, já que foram envolvidas neste novo código ferramentas que utilizavam
o Maple V3, o Maple 4.0 e o Maple 8.0. Além disto, foram utilizados comandos do
pacote Maplets para a criação de janelas, diálogos e outras interfaces visuais que in-
teragem com o usuário. Algumas características do pacote Maplets serão descritas na
próxima seção.
3.2 Pacote Maplets
O pacote Maplets do Maple é uma interface gráfica que p erm ite que o usuário combine
pacotes e procedimentos com janelas e caixas de diálogo interativas. Este pacote pode
ser utilizado para criar aplicações customizadas que exija o mínimo de conhecimento
da sintaxe do Maple.
É composto por três sub-pacotes: Elementos (que podem ser janelas, b otões, etc),
Exemplos (contendo diversos exemplos de aplicações do Maplets) e Ferramentas (que
contém rotinas para manipulação e interação com os elementos do Maplets), além de
uma função para execução do Maplet, Display.
Para a atualização dos códigos INDEX, ACIG e OTIMA e a criação do código
EVENTS foram necessárias algumas mudanças de sintaxes e comandos do Maple para
se utilizar o pacote Maplets. Alguns deles estão descritos abaixo.
O comando para impressão no Worksheet do Maplet, print, lprint, printf entre
outros, não podem ser utilizados na interface Maplets. Para a impressão em uma tela
de resultados composto por um TextBox, é necessário o comando
Maplets:-Tools:-Set(’xB2’(’appendline’)="Exemplo");
onde xB2 é a r eferência do elemento TextBox na interface. O comando appendline
é utilizado para adicionar uma nova linha de saída. Também pode ser utilizado o
comando append onde a saída será escrita na mesma linha anterior. Um outro comando
é o Get que deve ser utilizado somente em procedimentos, ou seja, não pode ser
utilizado na criação da interface. Este comando recupera o valor de um elemento de
uma aplicação ativa. Por exemplo, o comando
opta:=Maplets:-Tools:-Get(’teste’);
3.3. Descrição de ferramentas de tratamento de PCOADs 41
atribui o valor de teste à nova variável opta, onde teste é a referência de um elemento
na interface. O comando SetOption, permite que valores de certas opções mudem
enquanto uma aplicação do Maplet e stá em funcionamento. Por exemplo, a linha
abaixo, muda o valor do elemento iC3 de verdadeiro para falso.
SetOption(’iC3’(’value’)=false)
Na atualização dos códigos alguns comandos utilizados na versão original tornaram-se
obsoletos ou tiveram sua entrada modificada, devendo ser alterados para o funciona-
mento. Dentre eles, podem-se destacar (Quadro 3.1):
Quadro 3.1: Comandos do Maple.
Anterior Atual
traperror/lasterror try/catch
submatrix matadd
RETURN return
codegen CodeGeneration
‘type/ratpoly‘() type(X,ratpoly)
Uma dificuldade encontrada foi a impressão de matrizes dentro de um TexBox e
a utilização do comando CodeGeneration, que deve ser melhorada.
3.3 Descrição de ferramentas de tratamento de PCO-
ADs
As ferramentas integradas foram o INDEX (MURATA, 1996) desenvolvido para a car ac-
terização de EADs, o ACIG ( CUNHA; MURATA, 1999) que executa a redução simbólica
do índice diferencial baseado no Algoritmo de Gear e o OTIMA, desenvolvido inicial-
mente por Gomes (2000) e generalizado e aperfeiçoado por Lobato (2004), que gera o
Sistema Equivalente utilizado na solução de PCOADs pelo Método Indireto. A nova
ferramenta associada a esta nova interface, chamada EVENTS, é o algoritmo de Xu
e Antsaklis (2004) que gera o Sistema Equivalente parametrizado pelos eventos. Na
seqüência, são apresentadas as principais características de cada uma destas ferramen-
tas.
3.3. Descrição de ferramentas de tratamento de PCOADs 42
3.3.1 INDEX
O código simbólico INDEX foi desenvolvido por Murata (1996) usando o Maple V3
como linguagem de programação. O usuário deve fornecer como dados de entrada
um arquivo do tipo texto com a dimensão do sistema e as equações a analisar, usando
nomes exclusivos para as variáveis e suas derivadas, para as equações e para a dimensão
do sistema. A partir destas informações, o código faz a partição do sistema, verifica
se ele é semi-explícito, linear, linearmente implícito ou não-linear, identifica sistemas
implícitos de índice 1 ou superior ou sistema semi-explícitos de Hessenberg de índices
2 e 3, verifica a resolubilidade de sistemas lineares com coeficientes constantes e alerta
quanto à existência de problemas de inicialização em sistemas de índice 1.
O algoritmo é baseado na determinação da est rutura das equações e na depen-
dência das variáveis. O algoritmo simbólico pode ser usado para comprovar a caracte-
rização obtida através dos códigos estruturais, que fornecem apenas um valor inferior
para o índice e um valor superior de graus de liberdade.
3.3.2 ACIG
Os códigos de caracterização de equaçõ es algébrico-diferenciais, ALGO e PALGO,
possuem limitações que po dem ser quanto ao índice estrutural obtido através destes
códigos, que é apenas um limite inferior do índice diferencial e o número estrutural de
graus de liberdade dinamicamente definidos que é um limite superior do número de
graus de lib erdade dinamicamente definidos. Por causa destas limitações, foi desen-
volvido por Cunha e Murata (1999) uma ferramenta denominada ACIG (Algoritmo
de Análise de Consistência e Redução do Índice Diferencial baseado no Algoritmo de
Gear) para ser acoplada ao código INDEX (MURATA, 1996), já que o INDEX apenas
identifica sistemas de EADs com problemas de inicialização, mas não fornece infor-
mações ao usuário quanto ao número de graus de liberdade das equações ou quanto
ao sistema estendido correspondente. O ACIG é composto por três blocos principais:
Entrada de dados, Corpo do Programa e Saída de dados.
• Entrada de Dados
A entrada de dados pode ser feita através de banco de dados próprio da interface
desenvolvida e através de arquivo de dados do próprio usuário.
• Corpo do Programa
A princípio é analisada a resolubilidade do s istema de EADs, calculando-se a
soma da matriz de derivadas parciais das equações em relação as variáveis di-
ferenciais e em relação às algébricas. Para que o sistema possa ser resolvido o
determinante desta matriz não deve ser nulo.
O programa é baseado principalmente no cálculo do posto simbólico da matriz
de derivadas do sistema de EADs. O posto de uma matriz é o número máximo
3.3. Descrição de ferramentas de tratamento de PCOADs 43
de linhas ou colunas de uma matriz que definem um sistema linearmente inde-
pendente. Definida a m atriz de derivadas, o posto p ode avaliar a existência de
variáveis diferenciais no sistema, fazendo assim a partição em blocos de variáveis
diferenciais e de variáveis algébricas. O mesmo acontece com as equações do
sistema. Esse procedimento é iterativo e tem como critério de parada a igual-
dade entre o posto e o número de equações avaliadas. Resolve-se então o sistema
diferencial para as variáveis diferenciais e diferenciam-se o restante das equações.
Efetua-se então a substituição das variáveis diferenciais no sistema de equações
diferenciado anteriormente. O posto é novamente analisado. Quando a igual-
dade acontece, garante-se que o sistema foi reduzido a um sistema de EDOs.
Enquanto a diferença entre o posto e o número de equações avaliadas no laço
persistir, as diferenciações e manipulações algébricas deverão ocorrer.
• Saída de dados
Ao final da execução, o programa apresenta o número de graus de liberdade di-
namicamente definidos, o índice do sistema original e os sistemas de índice zero
e um são gerados. O programa também gera uma análise de consistência dos
graus de liberdade definidos durante a execução apresentando as condições inici-
ais consistentes. O tempo de execução também é mostrado ao final da execução.
O ACIG apresenta grande versatilidade e bons resultados na caracterização de
EADs, embora não seja aplicável a alguns tipos de dados, tais como: Equações dife-
renciais parciais (EDPs) em sua forma original, EADs de grande dimensão e altamente
não linear.
3.3.3 OTIMA
A solução numérica de PCOs através de métodos indiretos normalmente exige um
tratamento analítico considerável na obtenção das equações adjuntas, que pode se
tornar trabalhoso e conduzir a expressões complexas, tornando-os desvantajosos com
relação a outros métodos.
Para eliminar esta desvantagem, foi desenvolvido por Gomes (2000) e estendido
por Lobato (2004) o código OTIMA, implementado inicialmente em Maple V R4 e
atualizado para V8, que através da aplicação do Princípio de Pontryagin, faz o cál-
culo simbólico das condições de ótimo, através da geração das equações adjuntas e a
condição de mínimo. O OTIMA é capaz de lidar com problemas de controle ótimo
com restrições de ponto final, restrições de desigualdade e igualdade e com condições
laterais.
3.4. Aplicações da interface OpCol 44
3.3.4 EVENTS: geração do sistema equivalente parametrizado
pelos eventos
A obtenção das equaçõ es de sensibilidade a partir da diferenciação das variáveis em
relação a cada evento conforme algoritmo proposto por Xu e Antsaklis (2004), des-
crito em detalhes no Capítulo 2, é feita automaticamente em ambiente Maple 9.5 e
incorporada à ferramenta OpCol, juntamente com as equações de estado, co-estado e
estacionariedade, para sistemas chaveados. O sistema resultante é resolvido usando o
software MATLAB.
Entrada
Variáveis adjuntas
Condição de estacionariedade
Parametrização dos eventos
Diferenciação das equações
em relação aos eventos
Sistema
chaveado
Custo
Funcional
Saída
Sistema
Equivalente
Figura 3.1: Fluxograma do código EVENTS.
3.4 Aplicações da interface OpCol
Na seqüência, a interface desenvolvida é aplicada a quatro estudo de casos, que permi-
tem identificar suas principais características. Os dois primeiros exemplos, retirados
de Pantelides (1988) e de Murata (1996), mostram as ações de identificação de índice
superior, redução deste índice executadas pelo INDEX e pela ACIG. Os dois últimos
exemplos, retirados de Lobato (2004) e de Xu (2001) mostram respectivamente, as
gerações do Sistema Estendido Corresp ondente através da aplicação do PMP, e do
Sistema Parametrizado pelos Eventos, executadas pelo OTIMA e pelo EVENTS.
3.4. Aplicações da interface OpCol 45
3.4.1 Caso 1 - Sistema de Índice 1
Este exemplo é descrito por Pantelides (1988), que apresenta um sist ema linear de
EADs de índice 1.
˙z
1
+ z
1
+ 2y
1
= 0 (3.1)
z
1
+ 2z
2
= u1 (3.2)
Para ser feita a caracterização do problema foi utilizado o c ódigo INDEX e sua
saída é apresentada no Quadro (3.2).
Quadro 3.2: Saída do Código INDEX para o Caso 1.
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ +
+ SYMBOLIC ANALYSIS OF DIFFERENTIAL-ALGEBRAIC SYSTEMS +
+ +
+ INDEX 1.0 gmscp-PEQ/COPPE +
+ Revision FEQUI/UFU/2006 +
+ +
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
From Pantelides
==============================================
SYSTEM OF DIFFERENTIAL-ALGEBRAIC EQUATIONS
OF SEMI-EXPLICIT TYPE WITH DIMENSION 2
dx1+x1+2*y1 = 0
x1+2*y1 = u1
dim(f)=1
dim(g)=1
dim(x)=1
dim(y)=1
===============================================
Analysing the system’s linearity - Wait...
3.4. Aplicações da interface OpCol 46
LINEAR SYSTEM.
Cpu time for linearity’s analysis:0.sec.
=====================================================
Generating the pattern of matrizes Fz’and Fz - Wait...
FdzS = [[1, 0], [0, 0]]
FzS = [[1, 1], [1, 1]]
=====================================================
Analysing the system’s index - Wait...
Cpu time for index analysis:0.sec.
INDEX ONE SYSTEM.
================================================
Generating the Subroutines with Functions to DASSL and RADAU5 - Wait...
Cpu time for functions subroutine =0.sec.
RES
[YPRIME[1]+Y[1]+2*Y[2]]
[Y[1]+2*Y[2]-u1]
################## The End ######################
Nesta saída do INDEX, pode-se obter os padrões utilizados nos códigos ALGO
e PALGO (FdzS e FzS), além disso mostra se o sistema é de índice 1 ou superior.
Utilizando o código ACIG, pode-se gerar o sistema de índice 1 ou 0. Com esse código,
os sistemas de índice superior podem ter o sistema estendido correspondente gerado
facilitando a tarefa de diferenciação do sistema. A saída do ACIG é apresentado no
Quadro (3.3).
3.4. Aplicações da interface OpCol 47
Quadro 3.3: Saída do Código ACIG para o Caso 1.
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ +
+ SYMBOLIC ANALYSIS OF DIFFERENTIAL-ALGEBRAIC SYSTEMS +
+ +
+ ACIG 1.0 +
+ +
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
From Pantelides
==============================================
Original DAE system
Differential
diff(x1(t),t)+x1(t)+2*x2(t)
Algebraic
x1(t)+2*x2(t)-u1
==============================================
Degrees of Freedom
1
==============================================
ODE System
{diff(x2(t),t) = 1/2*x1(t)+x2(t), diff(x1(t),t) = -x1(t)-2*x2(t)}
==============================================
Index One DAE system
{diff(x1(t),t) = -x1(t)-2*x2(t)}
[x1(t)+2*x2(t)-u1]
==============================================
Cpu time: .155sec
################# The End ###################
3.4. Aplicações da interface OpCol 48
3.4.2 Caso 2 - R eator CSTR com reação exotérmica de pri-
meira ordem
Considere o modelo de operação transiente de um reator do tipo CSTR onde ocorre
uma reação exotérmica de primeira ordem A → B, e o calor de reação é removido
por uma camisa de refrigeração. Este sistema não-linear apresenta índice superior e é
dado por (MURATA, 1996):
˙x
1
=
F
Q
(Co − x
1
) −
ko
ρ
exp
−E
Rx
2
x
1
(3.3)
˙x
2
=
F
Q
(T o − x
2
) +
−∆Hko
ρCp
exp
−E
Rx
2
x
1
−
U
ρCpV
(x
2
− y
1
) (3.4)
x
2
− u
1
= 0 (3.5)
A saída do código INDEX, que é apresentada no Quadro (3.4), mostra que o
sistema apresenta índice superior (2) e através do código ACIG pode ser feita a redução
para índice 1 e índice 0. A saída do ACIG é apresentada no Quadro (3.5).
Quadro 3.4: Saída do Código INDEX para o Caso 2.
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ +
+ SYMBOLIC ANALYSIS OF DIFFERENTIAL-ALGEBRAIC SYSTEMS +
+ +
+ INDEX 1.0 gmscp-PEQ/COPPE +
+ Revision FEQUI/UFU/2006 +
+ +
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
CSTR with control of the temperature
==============================================
SYSTEM OF DIFFERENTIAL-ALGEBRAIC EQUATIONS
OF SEMI-EXPLICIT TYPE WITH DIMENSION 3
dx1-FF/Q*(Co-x1)+exp(ln(K)-E1/R1/x2)*x1 = 0
3.4. Aplicações da interface OpCol 49
dx2-FF/Q*(To-x2)+HH*exp(ln(K)-E1/R1/x2)*x1/Rho/Cp+U*(x2-y1)/Q/Rho/Cp = 0
x2-u1 = 0
dim(f)=2
dim(g)=1
dim(x)=2
dim(y)=1
===============================================
Analysing the system’s linearity - Wait...
NONLINEAR SYSTEM.
Cpu time for linearity’s analysis:.62e-1sec.
===============================================
Generating the pattern of matrizes Fz’and Fz - Wait...
FdzS = [[1, 0, 0], [0, 1, 0], [0, 0, 0]]
FzS = [[1, 1, 0], [1, 1, 1], [0, 1, 0]]
===============================================
Analysing the system’s index - Wait...
HIGHER INDEX SYSTEM.
Looking for a special structural form...
Cpu time for index analysis:.31e-1sec.
INDEX TWO HESSENBERG FORM SYSTEM:
X’=F1(X,Y)0=F2(x)
dim(f)= 2
dim(g)=1
=============================================
Generating the Subroutines with Functions to DASSL - Wait...
Cpu time for functions subroutine =0.sec.
RES
[YPRIME[1]-FF/Q*(Co-Y[1])+exp(ln(K)-E1/R1/Y[2])*Y[1]]
[YPRIME[2]-FF/Q*(To-Y[2])+HH*exp(ln(K)-E1/R1/Y[2])*Y[1]/Rho/Cp+
U*(Y[2]-Y[3])/Q/Rho/Cp]
[Y[2]-u1]
######################### The End ###############################
3.4. Aplicações da interface OpCol 50
Quadro 3.5: Saída do Código ACIG para o Caso 2.
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ +
+ SYMBOLIC ANALYSIS OF DIFFERENTIAL-ALGEBRAIC SYSTEMS +
+ +
+ ACIG 1.0 +
+ +
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
CSTR with control of the temperature
==============================================
Original DAE system
Differential
diff(x1(t),t)-(FF*Co-FF*x1(t)-exp(-(-ln(K)*R1*x2(t)+E1)/R1/x2(t))*
x1(t)*Q)/Q
diff(x2(t),t)-(-U*x2(t)+FF*Rho*Cp*To-FF*Rho*Cp*x2(t)-HH*exp(-(-ln(K)*
R1*x2(t)+E1)/R1/x2(t))*x1(t)*Q+U*x3(t))/Q/Rho/Cp
Algebraic
x2(t)-u1
(-U*x2(t)+FF*Rho*Cp*To-FF*Rho*Cp*x2(t)-HH*exp(-(-ln(K)*R1*x2(t)+
E1)/R1/x2(t))*x1(t)*Q+U*x3(t))/Q/Rho/Cp
==============================================
Degrees of Freedom
1
==============================================
ODE System
{diff(x1(t),t)=(FF*Co-FF*x1(t)-exp(-(-ln(K)*R1*x2(t)+E1)/R1/x2(t))
*x1(t)*Q)/Q,
diff(x2(t),t)=(-U*x2(t)+FF*Rho*Cp*To-FF*Rho*Cp*x2(t)-HH*exp(-(-ln(K)
*R1*x2(t)+E1)/R1/x2(t))*x1(t)*Q+U*x3(t))/Q/Rho/Cp,
3.4. Aplicações da interface OpCol 51
diff(x3(t),t)=-(U^2*R1*x2(t)^3-U*R1*x2(t)^2*FF*Rho*Cp*To+2*U*R1*
x2(t)^3*FF*Rho*Cp+U*R1*x2(t)^2*HH*exp(-(-ln(K)*R1*x2(t)+E1)/R1/x2(t))
*x1(t)*Q-U^2*R1*x2(t)^2*x3(t)-FF^2*R1*x2(t)^2*Rho^2*Cp^2*To+FF^2*R1*
x2(t)^3*Rho^2*Cp^2+2*FF*R1*x2(t)^2*Rho*Cp*HH*exp(-(-ln(K)*R1*x2(t)+E1)
/R1/x2(t))*x1(t)*Q-FF*R1*x2(t)^2*Rho*Cp*U*x3(t)+HH*E1*exp(-(-ln(K)*R1
*x2(t)+E1)/R1/x2(t))*x1(t)*Q*U*x2(t)-HH*E1*exp(-(-ln(K)*R1*x2(t)+E1)/
R1/x2(t))*x1(t)*Q*FF*Rho*Cp*To+HH*E1*exp(-(-ln(K)*R1*x2(t)+E1)/R1/
x2(t))*x1(t)*Q*FF*Rho*Cp*x2(t)+HH^2*E1*exp(-2*(-ln(K)*R1*x2(t)+E1)/
R1/x2(t))*x1(t)^2*Q^2-HH*E1*exp(-(-ln(K)*R1*x2(t)+E1)/R1/x2(t))*x1(t)
*Q*U*x3(t)-HH*exp(-(-ln(K)*R1*x2(t)+E1)/R1/x2(t))*R1*x2(t)^2*Q*Rho*Cp
*FF*Co+HH*exp(-2*(-ln(K)*R1*x2(t)+E1)/R1/x2(t))*R1*x2(t)^2*Q^2*Rho*Cp
*x1(t))/U/R1/x2(t)^2/Q/Rho/Cp}
==============================================
Index One DAE system
{diff(x1(t),t)=(FF*Co-FF*x1(t)-exp(-(-ln(K)*R1*x2(t)+E1)/R1/x2(t))*
x1(t)*Q)/Q,
diff(x2(t),t)=(-U*x2(t)+FF*Rho*Cp*To-FF*Rho*Cp*x2(t)-HH*exp(-(-ln(K)*
R1*x2(t)+E1)/R1/x2(t))*x1(t)*Q+U*x3(t))/Q/Rho/Cp}
[(-U*x2(t)+FF*Rho*Cp*To-FF*Rho*Cp*x2(t)-HH*exp(-(-ln(K)*R1*x2(t)+E1)
/R1/x2(t))*x1(t)*Q+U*x3(t))/Q/Rho/Cp]
==============================================
CPU Time: .283sec
################# The End ###################
A tarefa de diferenciação do sistema original para a obtenção do sistema de índice
1 ou zero é realizado pelo código ACIG com sucesso, ajudando o usuário na caracte-
rização de um sistema de EADs.
3.4.3 Caso 3 - O Problema de Newton formulado como EDOs
Este problema foi um dos primeiros resolvidos pelo Cálculo das Variações, por volta
de 1686, por Isaac Newton. O objetivo é a minimização do arraste na extremidade de
3.4. Aplicações da interface OpCol 52
um cone de raio da base r(t) e comprimento l, axialmente simétrico, num escoamento
hipersônico (BRYSON; HO, 1975; FEEHERY, 1998; LOBATO, 2004). O problema pode
ser escrito matematicamente como:
l
x
r
fluxo
Figura 3.2: Problema Isaac Newton.
min
u(x)
1
2
(r(l))
2
+
l
0
ru
3
1 + u
2
dx (3.6)
dr
dx
+ u = 0 (3.7)
r(0) = r
0
(3.8)
onde x é a distância axial, r é o raio do cone e a variável de controle é a tangente do
ângulo entre a direção da velocidade e a tangente local à parede do cone em função de
x. Aplicando o PMP, através do código OTIMA, obtemos:
dr
dx
+ u = 0 (3.9)
dλ
dx
−
u
3
1 + u
2
= 0 (3.10)
λ −
ru
2
(3 + u
2
)
(1 + u
2
)
2
= 0 (3.11)
3.4. Aplicações da interface OpCol 53
com as seguintes condições de contorno
r(0) = r
0
(3.12)
λ
l
= −r(l) (3.13)
Este problema que foi apresentado por Lobato (2004), possui as mesmas condições de
ótimo encontradas em Bryson e Ho (1975), exceto pelo sinal de λ e por um fator igual
a 4, que nesta referência multiplica o termo integral da função objetivo.
3.4.4 Caso 4 - Sistema Chaveado
O problema a seguir é descrito em Xu (2001) e através da interface criada, gera-se
o sistema parametrizado que é apresentado nos Quadros (3.6-3.8). Este problema é
resolvido no Capítulo 4.
J =
1
2
(x
1
(3) − e
2
)
2
+
1
2
(x
2
(3) − e
2
)
2
+
1
2
3
0
u
2
(t)dt (3.14)
Sub − sistema 1
˙x
1
= −x
1
+2x
1
˙x
2
= x
2
+x
2
u
(3.15)
Sub − sistema 2
˙x
1
= x
1
−3x
1
u
˙x
2
= 2x
2
−2x
2
u
(3.16)
Sub − sistema 3
˙x
1
= 2x
1
+x
1
u
˙x
2
= −x
2
+3x
2
u
(3.17)
com condições iniciais x
1
(0) = 1 e x
2
(0) = 1.
3.4. Aplicações da interface OpCol 54
Quadro 3.6: Conjunto de equações para a primeira fase.
Primeira Fase
0 ≤ τ < 1
u
1
= −2y(1)y(3) − y(2)y(4)
dx1dtau=ts1*(-y(1)+2*u1*y(1))
dx2dtau=ts1*(y(2)+u1*y(2))
dl1dtau=-ts1*(-y(3)+2*y(3)*u1)
dl2dtau=-ts1*(y(4)+u1*y(4))
dx1dts1=-y(1)+2*y(1)*(-2*y(1)*y(3)-y(2)*y(4))+ts1*(-y(5)+2*y(5)*(-2*y(1)*
y(3)-y(2)*y(4))+2*y(1)*(-2*y(5)*y(3)-2*y(1)*y(7)-y(6)*y(4)-y(2)*y(8)))
dx2dts1=y(2)+y(2)*(-2*y(1)*y(3)-y(2)*y(4))+ts1*(y(6)+y(6)*(-2*y(1)*y(3)-
y(2)*y(4))+y(2)*(-2*y(5)*y(3)-2*y(1)*y(7)-y(6)*y(4)-y(2)*y(8)))
dl1dts1=y(3)-2*y(3)*(-2*y(1)*y(3)-y(2)*y(4))-ts1*(-y(7)+2*y(7)*(-2*y(1)*
y(3)-y(2)*y(4))+2*y(3)*(-2*y(5)*y(3)-2*y(1)*y(7)-y(6)*y(4)-y(2)*y(8)))
dl2dts1=-y(4)-y(4)*(-2*y(1)*y(3)-y(2)*y(4))-ts1*(y(8)+y(8)*(-2*y(1)*y(3)-
y(2)*y(4))+y(4)*(-2*y(5)*y(3)-2*y(1)*y(7)-y(6)*y(4)-y(2)*y(8)))
dx1dts2=ts1*(-y(9)+2*y(9)*(-2*y(1)*y(3)-y(2)*y(4))+2*y(1)*(-2*y(9)*y(3)-
2*y(1)*y(11)-y(10)*y(4)-y(2)*y(12)))
dx2dts2=ts1*(y(10)+y(10)*(-2*y(1)*y(3)-y(2)*y(4))+y(2)*(-2*y(9)*y(3)-2*y(1)
*y(11)-y(10)*y(4)-y(2)*y(12)))
dl1dts2=-ts1*(-y(11)+2*y(11)*(-2*y(1)*y(3)-y(2)*y(4))+2*y(3)*(-2*y(9)*y(3)
-2*y(1)*y(11)-y(10)*y(4)-y(2)*y(12)))
dl2dts2=-ts1*(y(12)+y(12)*(-2*y(1)*y(3)-y(2)*y(4))+y(4)*(-2*y(9)*y(3)-2*y(1)
*y(11)-y(10)*y(4)-y(2)*y(12)))
3.4. Aplicações da interface OpCol 55
Quadro 3.7: Conjunto de equações para a segunda fase.
Segunda Fase
1 ≤ τ < 2
u
2
= 3y(1)y(3) + 2y(2)y(4)
dx1dtau(ts2-ts1)*(y(1)-3*u2*y(1))
dx2dtau=(ts2-ts1)*(2*y(2)-2*u2*y(2))
dl1dtau=-(ts2-ts1)*(y(3)-3*u2*y(3))
dl2dtau=-(ts2-ts1)*(2*y(4)-2*u2*y(4))
dx1dts1=-y(1)+3*y(1)*(3*y(1)*y(3)+2*y(2)*y(4))+(ts2-ts1)*(y(5)-3*y(5)*(3*y(1
)*y(3)+2*y(2)*y(4))-3*y(1)*(3*y(5)*y(3)+3*y(1)*y(7)+2*y(6)*y(4)+2*y(2)*y(8)))
dx2dts1=-2*y(2)+2*y(2)*(3*y(1)*y(3)+2*y(2)*y(4))+(ts2-ts1)*(2*y(6)-2*y(6)*(
3*y(1)*y(3)+2*y(2)*y(4))-2*y(2)*(3*y(5)*y(3)+3*y(1)*y(7)+2*y(6)*y(4)+2*y(2)*
y(8)))
dl1dts1=y(3)-3*y(3)*(3*y(1)*y(3)+2*y(2)*y(4))-(ts2-ts1)*(y(7)-3*y(7)*(3*y(1)*
y(3)+2*y(2)*y(4))-3*y(3)*(3*y(5)*y(3)+3*y(1)*y(7)+2*y(6)*y(4)+2*y(2)*y(8)))
dl2dts1=2*y(4)-2*y(4)*(3*y(1)*y(3)+2*y(2)*y(4))-(ts2-ts1)*(2*y(8)-2*y(8)*(3*y(
1)*y(3)+2*y(2)*y(4))-2*y(4)*(3*y(5)*y(3)+3*y(1)*y(7)+2*y(6)*y(4)+2*y(2)*y(8)))
dx1dts2=y(1)-3*y(1)*(3*y(1)*y(3)+2*y(2)*y(4))+(ts2-ts1)*(y(9)-3*y(9)*(3*y(1)*
y(3)+2*y(2)*y(4))-3*y(1)*(3*y(9)*y(3)+3*y(1)*y(11)+2*y(10)*y(4)+2*y(2)*y(12)))
dx2dts2=2*y(2)-2*y(2)*(3*y(1)*y(3)+2*y(2)*y(4))+(ts2-ts1)*(2*y(10)-2*y(10)*(
3*y(1)*y(3)+2*y(2)*y(4))-2*y(2)*(3*y(9)*y(3)+3*y(1)*y(11)+2*y(10)*y(4)+2*y(2)
*y(12)))
dl1dts2=-y(3)+3*y(3)*(3*y(1)*y(3)+2*y(2)*y(4))-(ts2-ts1)*(y(11)-3*y(11)*(3*y(1
)*y(3)+2*y(2)*y(4))-3*y(3)*(3*y(9)*y(3)+3*y(1)*y(11)+2*y(10)*y(4)+2*y(2)*
y(12)))
dl2dts2=-2*y(4)+2*y(4)*(3*y(1)*y(3)+2*y(2)*y(4))-(ts2-ts1)*(2*y(12)-2*y(12)*(3
*y(1)*y(3)+2*y(2)*y(4))-2*y(4)*(3*y(9*y(3)+3*y(1)*y(11)+2*y(10)*y(4)+2*y(2)*
y(12))))
3.4. Aplicações da interface OpCol 56
Quadro 3.8: Conjunto de equações para a terceira fase.
Terceira Fase
2 ≤ τ < 3
u
3
= −y(1)y(3) − 3y(2)y(4)
dx1dtau=(tf-ts2)*(2*y(1)+u3*y(1))
dx2dtau=(tf-ts2)*(-y(2)+3*u3*y(2))
dl1dtau=-(tf-ts2)*(2*y(3)+u3*y(3))
dl2dtau=-(tf-ts2)*(-y(4)+3*u3*y(4))
dx1dts1=(tf-ts2)*(2*y(5)+y(5)*(-y(1)*y(3)-3*y(2)*y(4))+y(1)*(-y(5)*y(3)-y(1)*
y(7)-3*y(6)*y(4)-3*y(2)*y(8)))
dx2dts1=(tf-ts2)*(-y(6)+3*y(6)*(-y(1)*y(3)-3*y(2)*y(4))+3*y(2)*(-y(5)*y(3)-
y(1)*y(7)-3*y(6)*y(4)-3*y(2)*y(8)))
dl1dts1=-(tf-ts2)*(2*y(7)+y(7)*(-y(1)*y(3)-3*y(2)*y(4))+y(3)*(-y(5)*y(3)-y(1)*
y(7)-3*y(6)*y(4)-3*y(2)*y(8)))
dl2dts1=-(tf-ts2)*(-y(8)+3*y(8)*(-y(1)*y(3)-3*y(2)*y(4))+3*y(4)*(-y(5)*y(3)-
y(1)*y(7)-3*y(6)*y(4)-3*y(2)*y(8)))
dx1dts2=-2*y(1)-y(1)*(-y(1)*y(3)-3*y(2)*y(4))+(tf-ts2)*(2*y(9)+y(9)*(-y(1)*
y(3)-3*y(2)*y(4))+y(1)*(-y(9)*y(3)-y(1)*y(11)-3*y(10)*y(4)-3*y(2)*y(12)))
dx2dts2=y(2)-3*y(2)*(-y(1)*y(3)-3*y(2)*y(4))+(tf-ts2)*(-y(10)+3*y(10)*(-y(1)
*y(3)-3*y(2)*y(4))+3*y(2)*(-y(9)*y(3)-y(1)*y(11)-3*y(10)*y(4)-3*y(2)*y(12)))
dl1dts2=2*y(3)+y(3)*(-y(1)*y(3)-3*y(2)*y(4))-(tf-ts2)*(2*y(11)+y(11)*(-y(1)*
y(3)-3*y(2)*y(4))+y(3)*(-y(9)*y(3)-y(1)*y(11)-3*y(10)*y(4)-3*y(2)*y(12)))
dl2dts2=-y(4)+3*y(4)*(-y(1)*y(3)-3*y(2)*y(4))-(tf-ts2)*(-y(12)+3*y(12)*(-y(1)
*y(3)-3*y(2)*y(4))+3*y(4)*(-y(9)*y(3)-y(1)*y(11)-3*y(10)*y(4)-3*y(2)*y(12)))
Onde: dx
1
/dτ = y(1), dx
2
/dτ = y(2), dλ
1
/dτ = y(3), dλ
2
/dτ = y(4), dx
1
/dt
s1
=
y(5), dx
2
/dt
s1
= y(6),dλ
1
/dt
s1
= y(7), dλ
2
/dt
s1
= y(8), dx
1
/dt
s2
= y(9), dx
2
/dt
s2
=
3.4. Aplicações da interface OpCol 57
y(10), dλ
1
/dt
s2
= y(11) e dλ
2
/dt
s2
= y(12).
Neste Capítulo, o intuito foi mostrar a interface cr iada para várias problemas den-
tro da abordagem algébrico-diferencial e esta apresentou-se bastante favorável para a
caracterização de sistemas, pois apresenta um conjunto de códigos que auxiliam na
pré-solução numérica de PCOADs, com o interesse de exigir do usuário pouco conhe-
cimento da sintaxe utilizada no software MAPLE. Um Manual do Usuário simplificado
para a ferramenta Opcol é apresentado no Apêndice A.
A interface criada auxilia no estudo de PCOADs, desde a caracterização do sis-
tema de EADs até a geração das equações adjuntas e condição de mínimo, exigidos
na solução por métodos indiretos. Para a solução de um PCOAD recomenda-se as
seguintes etapas: a caracterização do sistema através do código INDEX para a verifi-
car se o sistema possui índice superior e problemas de inicialização. Caso o problema
seja de índice superior, através do código ACIG se faz a redução do índice diferencial,
gerando o sistema estendido correspondente. A partir do sistema de índice 1 po de ser
utilizado o có digo OTIMA para geração das equações adjuntas e condição de mínimo.
Com este prévio estudo, há a possibilidade de se evitar erros desnecessários. O código
EVENTS, nesta interface, só pode ser aplicado a sistemas chaveados.
CAPÍTULO 4
Solução numérica de PCOADs com
flutuação do índice diferencial
4.1 Método de Solução
As soluções dos problemas investigados foram obtidas através de uma rotina imple-
mentada no ambiente MATLAB
R
, utilizando a função bvp4c e comparadas com a
solução de problemas sujeitos a restrições de desigualdade através do código DIRCOL,
implementado por fases, com estimativa dos eventos.
• Solução através do MATLAB
1. Inicia-se com estimativas para os eventos;
2. Calcula-se a variável de controle;
3. Calcula-se a função objetivo modificada e suas respectivas diferenciações
em relação a cada evento;
4. Atualiza os eventos usando o Método de Newton até atingir o critério de
parada.
• Solução através do código DIRCOL
4.2. Estudo de Casos 59
O segundo método de solução utilizado é baseado no método de colocação direta
implementada pelo código DIRCOL (STRYK, 1999), que resolve um PCO sujeito
a restriçõ es gerais de igualdade ou desigualdade nas variáveis de controle e/ou
de estado. Descontinuidades nas equações diferenciais podem ser tratadas como
problema de fases múltiplas. Pela discretização da variável de estado e controle,
o PCO é transcrito em uma seqüência de problemas não-lineares de otimização
com restrições (NLP). O código DIRCOL também fornece estimativas confiáveis
das variáveis adjuntas, das funções multiplicadoras das restrições de estado e da
função identificadora de fases.
4.2 Estudo de Casos
A metodologia apresentada foi aplicada a 3 estudo de casos que serão descritos a seguir.
4.2.1 Caso 1 - Sistema Chaveado
Este problema foi apresentado por Xu (2001), onde o objetivo é minimizar um custo
funcional J sujeito as restrições dadas pelos subsistemas.
J =
1
2
(x
1
(3) − e
2
)
2
+
1
2
(x
2
(3) − e
2
)
2
+
1
2
3
0
u
2
(t)dt (4.1)
Sub − sistema 1
˙x
1
= −x
1
+2x
1
˙x
2
= x
2
+x
2
u
(4.2)
Sub − sistema 2
˙x
1
= x
1
−3x
1
u
˙x
2
= 2x
2
−2x
2
u
(4.3)
Sub − sistema 3
˙x
1
= 2x
1
+x
1
u
˙x
2
= −x
2
+3x
2
u
(4.4)
São dadas as condições iniciais do problema: x
1
(0) = 1 e x
2
(0) = 1.
Através da interface desenvolvida, foram obtidas as equações que compõem este
sistema, descritas no Capítulo 3. Este problema foi resolvido utilizando tolerâncias
relativa e absoluta iguais a 10
−6
e escolhidos valores iniciais para t
s1
= 1, 2 e t
s2
= 2, 2.
A Tabela (4.1) apresenta os resultados do Caso 1.
Neste caso, não foi utilizado o código DIRCOL para comparação, mas sim sua
solução analítica. Pode-se p erceber que os resultados obtidos neste trabalho apresenta-
4.2. Estudo de Casos 60
Tabela 4.1: Resultados do Caso 1.
t
S1
) t
S2
Função Objetivo N
◦
de Iterações
Xu (2001) 1,0050 1,9993 0,0000027599 18
Solução Analítica
1,0 2,0 0
Este Trabalho
0,999991 2,000005 0,0000000000 18
ram valores semelhantes à solução analítica e à solução númerica obtida na literatura.
Os valores iniciais para os eventos foram os mesmos utilizados por Xu (2001). Foram
feitas variações nas condições iniciais dos eventos para verificar a dificuldade de con-
vergência do método. Os resultados destes test es são apresentados na Tabela (4.2).
Tabela 4.2: Variações nas estimativas iniciais.
Tol. 10
−6
t
S1
t
S2
N
◦
de Iterações t
S1
t
S2
J
1 1,2 2,2 18 0,999991 2,000005 0,0000000000
2
0,5 1,5 33 0,999991 2,000005 0,0000000000
3
1,0 2,0 13 0,999991 2,000005 0,0000000000
Tol. 10
−5
t
S1
t
S2
N
◦
de Iterações t
S1
t
S2
J
1 1,2 2,2 12 0,999987 2,000006 0,0000000001
2
0,5 1,5 28 0,999990 2,000008 0,0000000003
3
1,0 2,0 8 0,999990 2,000008 0,0000000002
Os valores iniciais para os eventos influenciaram na convergência do método de
solução, aumentando o número de iterações necessárias para a obtenção da solução
numérica, admitindo uma tolerância de 10
−6
. Reduzindo a tolerância para 10
−5
, houve
a redução no número de iterações necessárias, mas a função objetivo e os valores dos
eventos não sofreram alterações significativas. As trajetórias de estado e controle são
mostradas nas Figuras (4.1) - (4.6).
4.2. Estudo de Casos 61
0 0.5 1 1.5 2 2.5 3
0
1
2
3
4
5
6
7
8
τ
x
1
(τ)
Este trabalho
Solução analítica
Figura 4.1: Caso 1 - Variável de estado x
1
.
0 0.5 1 1.5 2 2.5 3
0
5
10
15
20
25
τ
x
2
(τ)
Este trabalho
Solução analítica
Figura 4.2: Caso 1 - Variável de estado x
2
.
4.2. Estudo de Casos 62
0 0.5 1 1.5 2 2.5 3
0
0.5
1
1.5
2
2.5
3
3.5
4
x 10
−7
τ
λ
1
(τ)
Figura 4.3: Caso 1 - Variável de co-estado λ
1
.
0 0.5 1 1.5 2 2.5 3
−1.5
−1
−0.5
0
x 10
−7
τ
λ
2
(τ)
Figura 4.4: Caso 1 - Variável de co-estado λ
2
.
4.2. Estudo de Casos 63
0 1 2 3 4 5 6 7 8
0
5
10
15
20
25
x
1
(τ)
x
2
(τ)
Este trabalho
Solução Analítica
Figura 4.5: Trajetória de estado.
0 0.5 1 1.5 2 2.5 3
−1.5
−1
−0.5
0
0.5
1
1.5
2
2.5
3
x 10
−7
τ
u(τ)
Figura 4.6: Caso 1 - Variável de controle.
4.2. Estudo de Casos 64
4.2.2 Caso 2 - Reator semi-batelada isotérmico com reações
paralelas e restrição de seletividade
O caso proposto é descrito em Srinivasan et al. (2003). A descrição do problema é:
1. Características:
• Reação A + B → C
2B → D
• Condições: Semi-batelada, reação exotérmica;
• Objetivo: Maximizar a produção de C (tempo conhecido);
• Variável Manipulada: Taxa de alimentação de B;
• Restrições: Limites na variável de controle, nas concentrações dos compo-
nentes B e D no tempo final;
2. Formulação do Problema:
O problema pode ser descrito como:
min
u(t)
J = V (t
f
)c
C
(t
f
) (4.5)
sujeito a:
dc
A
dt
= −k
1
c
A
c
B
−
u
V
c
A
(4.6)
dc
B
dt
= −k
1
c
A
c
B
− 2k
2
c
2
B
+
u
V
(c
Bin
− c
B
) (4.7)
dV
dt
= u (4.8)
c
C
=
c
A0
V
0
− c
A
V
V
(4.9)
c
D
=
(c
A
+ c
Bin
− c
B
) V − (c
A0
+ c
Bin
− c
B0
) V
0
2V
(4.10)
c
B(t
f
)
≤ c
Bf max
(4.11)
c
D(t
f
)
≤ c
Df max
(4.12)
u
min
≤ u ≤ u
max
(4.13)
c
A
(0) = c
A0
(4.14)
c
B
(0) = c
B0
(4.15)
V (0) = V
0
(4.16)
4.2. Estudo de Casos 65
c
C
(0) = c
C0
(4.17)
c
D
(0) = c
D0
(4.18)
onde c
X
são as concentrações das espécies A, B, C e D, respectivamente, u é
taxa de alimentação do componente B, V é o volume do reator, k
1
e k
2
são
parâmetros cinéticos.
3. Parâmetros e Condições Operacionais do Modelo: Os parâmetros e condições
operacionais do modelo são mostrados na Tabela (4.3)
Tabela 4.3: Parâmetros e condições Operacionais do Modelo (Caso 2).
k
1
0,053 L/(mol min)
k
2
0,128 L/(mol min)
c
Ao
0,72 mol/L
c
Bo
0,05 mol/L
c
Co
0 mol/L
c
Do
0 mol/L
c
Bin
5 mol/L
c
Bfmax
0,025 mol/L
c
Df max
0,15 mol/L
u
min
0 L/min
u
max
0,001 L/min
V
0
1 L
De acordo com Srinivasan et al. (2003), o problema pode ser descrito em 3 fases,
devido às restrições de desigualdade nas variáveis de estado e controle. Estas
fases são obtidas através da análise física do problema e podem ser descritas
conforme abaixo:
4. Solução Ótima:
• Primeira Fase: devemos ter u = u
max
, para aumentarmos B o mais rápido
quanto possível.
• Segunda Fase: limitar a quantidade de D produzida.
u =
c
B
V (k
1
c
A
(2c
Bin
− c
B
) + 4k
2
c
B
c
Bin
)
2c
Bin
(c
Bin
− c
B
)
• Terceira Fase: devemos ter u = u
min
para alcançarmos c
B
em t
f
.
Com essa descrição das fases do problema, nota-se que o problema original possui
índice um, mas com a ativação da restrição na segunda fase o problema muda de
índice, passando para índice 3, descrita na Tabela (4.4). A função identificadora
de fases foi obtida conforme descrita em Srinivasan et al. (2003) que obtém a
4.2. Estudo de Casos 66
expressão analítica para a variável de controle, determina da na segunda fase.
Com a determinação desta expressão o sistema possui índice 1.
Tabela 4.4: Caso 2 - Função Identificadora de Fases.
Função Identificadora de Fases Índice
Primeira Fase u = u
max
1
Segunda Fase
u =
c
B
V (k
1
c
A
(2c
Bin
−c
B
)+4k
2
c
B
c
Bin
)
2c
Bin
(c
Bin
−c
B
)
3
Terceira Fase
u = u
min
1
5. Resultados:
O problema foi resolvido utilizando-se tolerância relativa e absoluta iguais à 10
−7
,
após 16 iterações partindo-se de [t
s1
, t
s2
]=[19, 202].
Tabela 4.5: Resultados do Caso 2.
t
S1
(min) t
S2
(min) n
c
(t
f
)
Srinivasan et al. (2003) 20,25 205 0,43
DIRCOL
20,249041 205,00547 0,43058
Este trabalho
20,2495 205,0013 0,43059
Pode-se perceber que os resultados obtidos neste trabalho apresentaram valores
semelhantes à solução pelo código DIRCOL e a solução obtida n a literatura.
Alterando os valores iniciais para os eventos, para verificar a convergência do
método, obtém-se como resultados, os valores da Tabela (4.6).
Tabela 4.6: Variações nas estimativas iniciais.
t
S1
(min) t
S2
(min) t
S1
(min) t
S2
(min) n
c
(t
f
)
1 19 202 20,2495 205,0013 0,43059
2
12 180 20,2495 205,0013 0,43059
3
25 220 20,2505 204,9987 0,43059
As variações nos valores iniciais dos eventos influenciaram na convergência do
método de solução, aumentando o número de iterações necessárias para a ob-
tenção da solução numérica, mas a função objetivo e os valores dos eventos não
sofreram alterações significativas. As trajetórias de estado e controle são mos-
tradas nas Figuras ( 4.7) - (4.10).
4.2. Estudo de Casos 67
Figura 4.7: Caso 2 - Evolução das concentrações.
Figura 4.8: Caso 2 - Função Objetivo.
4.2. Estudo de Casos 68
Figura 4.9: Perfil ótimo da taxa de alimentação.
4.2. Estudo de Casos 69
Figura 4.10: Caso 2 - Perfil do Volume.
4.2.3 Caso 3 - Reator Semi-Batelada Isotérmico
O caso proposto também é apresentado em Srinivasan et al. (2003). A descrição do
problema é:
1. Características:
• Reação A + B → C
• Condições: Semi-batelada, reação exotérmica;
• Objetivo: Minimizar o tempo necessário para a produção de C (valor co-
nhecido);
• Variável Manipulada: Taxa de alimentação de B;
• Restrições: Limites na variável de controle, de temperatura da camisa e de
volume;
2. Formulação do Problema:
O problema pode ser descrito como:
min
t
f
,u(t)
J = t
f
(4.19)
4.2. Estudo de Casos 70
sujeito a:
dc
A
dt
= −kc
A
c
B
−
u
V
c
A
(4.20)
dc
B
dt
= −kc
A
c
B
−
u
V
(c
Bin
− c
B
) (4.21)
dV
dt
= u ( 4.22)
c
C
=
c
A0
V
0
+ c
C0
V
0
− c
A
V
V
(4.23)
T
cf
= T (t) + min (c
A
, c
B
)
(−∆H)
ρc
p
(4.24)
T
cf
(t) ≤ T
max
(4.25)
V (t
f
) ≤ V
max
(4.26)
n
C
(t
f
) ≥ n
Cdes
(4.27)
u
min
≤ u ≤ u
max
(4.28)
c
A
(0) = c
A0
(4.29)
c
B
(0) = c
B0
(4.30)
V (0) = V
0
(4.31)
c
C
(0) = c
C0
(4.32)
onde c
X
são as concentrações das espécies A, B e C, respectivamente, u é taxa
de alimentação do componente B, V é o volume do reator, k é um parâmetro
cinético, T , a temperatura do reator, T
cf
é a temperatura do resfriador, ∆H a
entalpia da reação, ρ a densidade e cp a capacidade calorífica.
3. Parâmetros e Condições Operacionais do Modelo:
Os parâmetros e condições operacionais do modelo são mostrados na Tabela (4.7)
4. Condições Operacionais:
Inicialmente sabemos que o número de mols de B é sempre menor que o de A,
então implica que c
B
(t) ≤ c
A
(t). Com T
cf
≤ T
max
, implica que c
B
≤ c
Bmax
, com
c
Bmax
= ρcp(T max − T )/(−∆H). Portanto, inicialmente temos que ter a maior
quantidade de B possível.
5. Solução Ótima:
• Primeira Fase: devemos ter c
B
= c
Bmax
, para mantermos a restrição na
temperatura ativa, o que nos revela a expressão para a variável de controle
(T
cf
= T
max
).
u =
kc
A
c
B
V
c
Bin
− c
B
c
B
=c
B max
4.2. Estudo de Casos 71
Tabela 4.7: Parâmetros e condições Operacionais do Modelo (Caso 3).
k 0,0482 L/(mol h)
u
min
0 L/h
u
max
0,1 L/h
c
A0
2 mol/L
c
B0
0,63 mol/L
c
C0
0 mol/L
n
Cdes
0,6 mol
c
Bin
2 mol/L
∆H -60000 J/mol
ρ 900 g/L
T 70
o
C
T
max
80
o
C
cp 4,2 J/gK
V
max
1 L
V
0
0,7 L
• Segunda Fase: Quando a restrição do volume fica ativa, implica que u =
u
min
.
O t empo final é obtido quando a composição desejada para o componente
C é alcançada (n
C
= n
Cdes
).
Com essa descrição das fases do problema, nota-se que na segunda fase o sistema
passa de índice um para índice 2, descrita na Tabela (4.8).
Tabela 4.8: Caso 3 - Função Identificadora de Fases.
Função Identificadora de Fases Índice
Primeira Fase u =
kc
A
c
B
V
c
Bin
−c
B
c
B
=c
B max
2
Segunda Fase
u = u
min
1
n
c
= n
Cdes
6. Resultados:
Para se resolver este problema, foi feita uma mudança na função objetivo deste
problema, que é de tempo final livre. Em vez de minimizar o tempo necessário
para a produção do componente C, o tempo final da reação foi obtido ao alcan-
çar o número de mols do desejado para o componente C (nc
des
= 0, 6). Esta
mudança foi necessária pois o algoritmo proposto por Xu e Antsaklis (2004) ne-
cessita do tempo inicial(t
0
) e do tempo final(t
f
). Com essa mudança, o problema
foi resolvido utilizando-se tolerância relativa e absoluta iguais à 10
−6
, após 13 ite-
rações partindo-se de t
s1
igual a 10h. Os resultados são apresentados na Tabela
(4.9).
4.2. Estudo de Casos 72
Tabela 4.9: Resultados do Caso 3.
t
S1
(h) t
f
(h)
Srinivasan et al. (2003) 11,44 19,80
DIRCOL
11,44481 19,800916
Este Trabalho
11,44478 19,801017
Pode-se perceber que os resultados obtidos neste trabalho apresentaram valores
semelhantes à solução pelo código DIRCOL e a solução obtida na literatura. Fo-
ram feitos testes variando as estimativas iniciais para os eventos, e os resultados
obtidos são apresentados na Tabela (4.10). Testes variando os parâmetros do
modelo (k e c
Bin
) utilizando o código DIRCOL foram realizados e os resultados
obtidos são apresentados na Tabela (4.11).
Tabela 4.10: Variações nas estimativas iniciais.
t
S1
(h) t
f
(h) t
S1
(h) t
f
(h)
1 10 25 11,44478137 19,80101779
2
5 20 11,44478077 19,80101806
3
12 20 11,44485676 19,80095945
Tabela 4.11: Variações nos parâmetros.
k c
bin
J IFAIL
1 0,0482 2,0 19,800916 0
2
0,00482 2,0 21,000000 3
3
0,482 2,0 20,997092 3
4
0,0482 1,0 21,000000 3
As variações nos valores iniciais dos eventos influenciaram na convergência do
método de solução, aumentando o número de iterações necessárias para a ob-
tenção da solução numérica, mas a função objetivo e os valores dos eventos não
sofreram alterações significativas. As variações nos parâmetros, retornaram o
valor de IFAIL = 3 o que significa, de acordo com o código DIRCOL, que a solu-
ção não convergiu pois não foi p ossível achar uma solução ótima que atendesse
as restrições. As trajetórias de estado e controle são apresentadas nas Figuras
(4.11) - (4.12).
4.2. Estudo de Casos 73
Figura 4.11: Perfis das Concentrações.
4.2. Estudo de Casos 74
Figura 4.12: Perfil ótimo da taxa de alimentação.
4.2. Estudo de Casos 75
Figura 4.13: Caso 3 - Perfil da Temperatura.
Figura 4.14: Caso 3 - Perfil do Volume.
4.2. Estudo de Casos 76
Neste capítulo foram realizados três estudo de casos que apresentaram resultados
satisfatórios em relação aos encontrados na literatura. Foram considerados para os
casos resolvidos valores iniciais para os eventos de acordo com a literatura. Testes para
a análise da convergência do método foram realizados e para os casos considerados,
os valores iniciais para os eventos podem aumentar o número de iterações necessárias
para a solução numérica do problema. Um teste alterando os valores dos parâmetros
do último caso foi feito para avaliar se possui sensibilidade paramétrica na solução dos
problema.
A solução numérica através do código DIRCOL para os casos descritos por Srini-
vasan et al. (2003) foi realizada utilizando a formulação multi-fásica com estimativas
iniciais para os eventos, e mostrou-se bastante eficaz para a solução dos problemas
através do método direto. No apêndice B, é apresentada a formulação no código
DIRCOL para o caso 3.
A solução de PCO utilizando a parametrização dos eventos requer a sequência
de mudanças entre as fases, fazendo com que este método torne-se de uso restrito.
Neste estudo de casos considerou-se problemas com apenas uma variável manipulada,
devido a dificuldade de estender esta abordagem para sistema s com múltiplas entradas,
apenas analisando fisicamente o problema e considerando as ativações e desativações
das restrições ao longo da trajetória.
CAPÍTULO 5
Conclusões e Sugestões
5.1 Conclusões
Nesta dissertação a abordagem algébrico-diferencial da solução de Problemas de Con-
trole Ótimo sujeitos a restrições de desigualdade foi associada com a técnica de para-
metrização de Eventos desenvolvida originalmente para Problemas Chaveados, com o
objetivo de estimar os Eventos utilizando o Método Indireto de solução. Para isto, foi
utilizada a determinação de Funções Identificadoras de Fase que expressam a variável
de controle analiticamente em cada fase, baseada na condição de otimalidade e na aná-
lise física do problema para distinguir seqüências apropriadas de ativação/desativação
das restrições, que causam flutuações no índice diferencial das equações.
A geração do sistema equivalente parametrizado composto pelas equações de es-
tado, co-estado e estacionariedade, condições de contorno e continuidade além das di-
ferenciações das variáveis de estado e controle em relação a cada evento, feita através
das equações de sensibilidade foram geradas simbolicamente a partir de uma interface
específica desenvolvida em Maple
R
9.5 para sistemas chaveados. Para um sistema
chaveado de pequena dimensão, composto de 3 sub-sistemas com 2 equações em cada
um, gera-se um sistema parametrizado correspondente composto de 36 equações. Esta
tarefa trabalhosa de geração deste sistema foi eliminada com a implementação do có-
digo EVENTS, viabilizando a aplicabilidade desta metodologia para sistemas de maior
dimensão.
5.2. Sugestões 78
A atualização de ferramentas automáticas para o estudo de PCOADs foi rea-
lizada, unindo códigos que fazem a caracterização das EADs em relação ao índice,
resolubilidade e condições iniciais consistentes (INDEX), a implementação simbólica
do algoritmo de Gear que reduz o índice e gera o sistema de índice 1 correspondente,
feito pelo código ACIG e a geração das condições de otimalidade estendidas, atra-
vés das equações de Euler-Lagrange, obtidas pela aplicação do Princípio do Mínimo
de Pontryagin através do código OTIMA. Estas ferramentas, juntamente com o có-
digo EVENTS foram associadas numa interface amigável com o usuário, denominada
OpCol.
A solução numérica dos PCOADs foi feita através de um método direto com para-
metrização nas variáveis de controle e nas variáveis de estado implementados no código
DIRCOL utilizando uma formulação multifásica, com estimativa das variáveis adjun-
tas e dos Eventos. A solução obtida foi comparada com a obtida com um algoritmo
implementado em MATLAB
R
6.0, que resolve um problema de valor no contorno que
corresponde ao sistema parametrizado pelos eventos.
Os PCOADs que compuseram os estudos de casos desta dissertação foram escolhi-
dos para representar casos típicos com flutuação de Índice que pudessem ser associados
à metodologia de parametrização utilizada originalmente para os PCO chaveados. A
comparação entre os resultados obtidos pelo méto do direto e pelo método indireto as-
sociado à parametrização mostra que nos casos estudados os dois métodos apresentam
desempenho satisfatório, ressaltando que questões ligadas ao refinamento da discreti-
zação, estimativa das variáveis adjuntas, dos Eventos, do estado e do controle em cada
fase no caso da DIRCOL e sensibilidade do método indireto às estimativas iniciais dos
Eventos e a dificuldade de obtenção de Funções Identificadoras de Fase analiticamente
para sistemas de maior dimensão devem ser cuidadosamente consideradas na solução
de novos PCOADs.
Uma contribuição significativa desta Dissertação foi demonstrar através da im-
plementação prática a equivalência entre sistemas híbridos com comportamento con-
tínuo/discreto acoplados e os PCOAD com restrições de desigualdade que causam
flutuação no índice diferencial, identificada pioneiramente por Feehery (1998).
5.2 Sugestões
• Avaliar o desempenho dos métodos direto e indireto para sistemas de grande
dimensão;
• Estudar algoritmos de aceleração da convergência do método numérico de solução
do sistema parametrizado;
• Avaliar o efeito da seqüência de ativação e desativação das restrições pré-definida
sobre os resultados obtidos;
5.2. Sugestões 79
• Aprimorar a interface Opcol desenvolvida, c riando alternativas de entrada de
dados e de impressão de resultados e promovendo o acoplamento entre as ferra-
mentas que a compõem;
• Estender o método indireto com parametrização de Eventos para problemas com
múltiplas entradas e restrições.
REFERÊNCIAS BIBLIOGRÁFICAS
ASCHER, U.; SPITERI, R. Collocation software for boundary-value differential-
algebraic equations. SIAM J. Scient. Comput., v. 15, p. 938–952, 1994.
BACHMANN, R. et al. On methods for reducing the index of differential algebraic
equations. Computers and Chemical Engineering, v. 14, p. 1271–1273, 1990.
BARTON, P. I. et al. Dynamic optimization in a discontinuous world. Ind. Eng.
Chem. Res., v. 37, p. 966–981, 1998.
BEMPORAD, A.; GIUA, A.; SEATZU, C. An iterative algorithm for the optimal
control of continuous-time switched linear systems. IEEE Computer Society,
2002.
BIEGLER, L. T. Solution of dynamic optimization problems by sucessive quadratic
programming and orthogonal collocation. Computers and Chemical Engineering,
v. 8, p. 243–248, 1984.
BIEGLER, L. T.; CERVANTES, A. M.; WACHTER, A. Advances in simultaneous
strategies for dynamic process optimization. Computers and Chemical
Engineering, v. 57, p. 575–593, 2002.
BRYSON, A. E.; HO, Y. C. Applied optimal control. [S.l.]: Hemisphere Publishing
Washington, 1975.
CERVANTES, A.; BIEGLER, L. Large scale DAE optimization using a simultaneous
NLP formulation. AIChE Journal, v. 44, n. 5, p. 1038–1050, 1998.
CHUDEJ, K.; GÜNTHER, M. Global state space approach for the efficient numerical
solution of state-constrained trajectory optimization problems. Journal of
Optimization, Theory and Applications, v. 103, p. 75–93, 1999.
Referências Bibliográficas 81
CUNHA, J.; MURATA, V. V. Desenvolvimento de uma ferramenta de caracterização
automática e de redução do índice diferencial de modelos matemáticos de
processos de engenharia química do tipo algébrico diferencial de valor inicial.
1999. Relatório de Iniciação Científica, Uberlândia, UFU.
CUTHREL, J.; BIEGLER, L. T. On the optimization of differential-algebraic process
systems. AIChE Journal, v. 33, n. 8, p. 1257–1270, 1987.
CUTHREL, J.; BIEGLER, L. T. Simultaneous s olution and optimization of process
flowsheets with differential equation models. Chem. Eng. Res. Des., v. 64, p.
341–346, 1987.
CUTHREL, J.; BIEGLER, L. T. Simultaneous optimization and solution methods for
batch reactor control profiles. Computers Chem. Engng., v. 13, p. 49–62, 1989.
EGERSTEDT, M.; WARDI, Y.; DELMOTTE, F. Optimal control of switching times
in switched dynamical systems. 2003.
ÇELIK, E.; KARADUMAN, E.; BAYRAM, M. Numerical solution of chemical
differential-algebraic equations. Applied Mathematics and Computation, v. 139,
p. 259–264, 2003.
FEEHERY, W. F. Dynamic optimization with path constraints. Tese (Doutorado) —
Massachusetts Institute of Technology, 1998.
FEEHERY, W. F.; BARTON, P. I. Dynamic optimization with state variable path
constraints. Computers Chem, Eng, v. 22, n. 9, p. 1241–1256, 1998.
GEAR, C. W. Differential algebraic equation index transformations. SIAM J. Sci.
Stat. Comput., v. 9, p. 39–47, 1988.
GEAR, C. W.; PETZOLD, L. R. Ode method for the solution of differential/algebraic
systems. SIAM J. Numer. Anal., v. 21, p. 716–728, 1984.
GILL, P. E.; MURRAY, W.; SAUNDERS, M. A. User’s Guide for SNOPT 5.3: a
Fortran Package for Large-Scale Nonlinear Programming. 1997. Department of
Mathematics, University of California, San Diego.
GILL, P. E. et al. User’s Guide for NPSOL: a Fortran package for nonlinear
programming. 1986. Department of Mathematics, University of California, San
Diego.
GOMES, E. O. Abordagem algébrico diferencial na solução de problemas de controle
ótimo. Dissertação (Mestrado) — Universidade Federal de Uberlândia, 2000.
LOBATO, F. S. Abordagem mista para problemas de otimização dinâmica. Dissertação
(Mestrado) — Universidade Federal de Uberlândia, 2004.
LOGSDON, J. S.; BIEGLER, L. T. Accurate solution of diferential-algebraic
optimization problems. Ind Eng. Chem. Res., v. 28, p. 1628–1639, 1989.
Referências Bibliográficas 82
LOGSDON, J. S.; BIEGLER, L. T. Decomposition strategies for large-scale dynamic
optimization problems. Chemical Engineering Science, v. 47, n. 4, p. 851–864,
1992.
LYNN, L. L.; PARKIN, E. S.; ZAHRADNIK, R. L. Near-optimal control by
trajectory approximation. I&EC Fundamentals, v. 9, n. 1, p. 58–63, 1970.
LYNN, L. L.; ZAHRADNIK, R. L. The use of orthogonal polynomials in the
near-optimal control of distributed systems by trajectory approximation. Int. J.
Control, v. 12, n. 6, p. 1079–1087, 1987.
MAPLESOFT. Maple 9.5 - Introductory Programming Guide. 2003. M apleSoft.
MODAK, J. L.; LIM, H. C.; TAYLEB, Y. J. General characteristics of optimal
feed rate profiles for various fed-batch fermentation process. Biotechnology and
Bioengineering, XXVIII, p. 1396–1407, 1986.
MURATA, V. V. C aracterização simbólica de equações algébrico diferenciais por um
sistema de álgebra computacional com aplicações na engenharia química. Tese
(Doutorado) — COPPE/UFRJ, 1996. Rio de Janeiro.
OBERLE, H. J.; SOTHMANN, B. Numerical computation of optimal feed rates for a
fed-batch fermentation mo del. Journal of Optimization Theory and Applications,
v. 100, p. 1–13, 1999.
OH, S. H.; LUUS, R. Use of orthogonal collocation method in optimal control
problems. Int. J. Control, v. 26, p. 657–673, 1977.
PANTELIDES, C. The consistent initialization of differential algebraic systems.
SIAM J. SCI. STAT. Compt., v. 9, p. 213–231, 1988.
RENFRO, J. G. Computational studies in the optimization of systems described by
differential-alge braic equations. Tese (Doutorado) — University of Houston,
1986.
RENFRO, J. G.; MORDSHEDI, A. M.; ASBJORNSEN, O. A. Simultaneous
optimization solution of systems described by differential-algebraic equations.
Comput. and Chem. Engng., v. 11, p. 503–517, 1987.
SAN, K. Y.; STEPHANOPOULOS, G. A note on the optimality criteria for maximum
biomass production in a fed-batch fermenter. Biotechnology and Bioengineering,
v. 26, p. 1261–1264, 1984.
SAN, K. Y.; STEPHANOPOULOS, G. Optimization of fed-batch penicillin
fermentation: a case of singular optimal control with state constraints.
Biotechnology and Bioengineering, v. 34, p. 72–78, 1989.
SANTOS, K. G. et al. Controle ótimo da fermentação alcoólica em reator batelada
alimentada. 2005. Congresso Brasileiro de Engenharia Química - Iniciação
Científica.
Referências Bibliográficas 83
SEYWALD, H.; CLIFF, E. M.; WELL, K. H. Range optimal trajectories for an
aircraft flying in the vertical plane. Journal of Guidance, Control and Dynamics,
v. 17, p. 389–398, 1994.
SRINIVASAN, B.; PALANKI, S.; BONVIN, D. Dynamic optimization of batch
processes: I. characterization of the nominal solution. Computers and Chemical
Engineering, v. 27, p. 1–26, 2003.
STRYK, O. v.; SCHLEMMER, M. Optimal control of the industrial robot manutec
r3. International Series of Numerical Mathematics, v. 1, p. 367–382, 1994.
STRYK, O. von. User’s Guide for DIRCOL - a Direct Collocation method for the
numerical solution of optimal control problems. 1999. Fachgebiet Simulation und
Systemoptimierung (SIM). Technische Universitat Darmstadt.
STRYK, O. von. Numerical hyb rid optimal control and related topics. 2000. Technische
Universität München.
UNGER, J.; KRONER, A.; MARQUARDT, W. Structural analysis of differential
algebraic equations systems: theory and applications. Computers Chem. Engng,
v. 19, p. 867–882, 1995.
VASANTHARAJAN, S.; BIEGLER, L. T. Simultaneous strategies for optimization
of differential-algebraic systems with enforcement error criteria. Comput. Chem.
Engng., v. 14, n. 10, p. 1083–1100, 1990.
XU, X. Analysis and design of switched systems. Tese (Doutorado) — Department of
Electric Engineering, University of Notre Dame, 2001.
XU, X.; ANTSAKLIS, P. J. Optimal control of switched systems based on
parameterization of the switching instants. IEEE Transactions on Automatic
Control, v. 49, p. 1–16, 2004.
APÊNDICE A
Manual - OpCol
Neste apêndice é apresentado um manual para o usuário conhecer a interface, a entrada
de dados e geração de resultados do OpCol. A tela principal do programa é mostrada
na Figura A.1. Nest a tela é possível selecionar qual código o usuário deseja utilizar.
A seleção é feita clicando no botão de cada código. As Figuras A.2 - A.5 apresentam
a tela de entrada do código INDEX, OTIMA, ACIG e EVENTS, respectivamente.
Figura A.1: Tela de Entrada.
85
Figura A.2: Tela de Entrada do Código INDEX.
Figura A.3: Tela de Entrada do Código OTIMA.
Cada tela possui as seguintes opçõ es:
•Examples - Apresenta um banco de dados para cada código.
•Via own data files - O usuário pode entrar com um arquivo próprio através desta
opção, de acordo com a sintaxe utilizada em cada código (Figura A.6).
Para o código INDEX, a opção Examples possui o seguinte banco de dados, des-
crito em (MURATA, 1996):
86
Figura A.4: Tela de Entrada do Código ACIG.
Figura A.5: Tela de Entrada do Código EVENTS.
•Index 0
1.Van der Pol - Hairer e Wanner (1993)
•Index 1
1.Adsorption column - Von Meien e Biscaia (1994)
2.From Pantelides - Pantelides (1988)
3.Finite Bath - Mendonça e Biscaia (1995)
87
Figura A.6: Entrada via arquivo do usuário.
•Index 2
1.From Pantelides - Pantelides (1988)
2.Trajectory prescribed path control - Brenan et al (1989)
3.CSTR with temperature control - McLellan (1994)
4.Condenser with fixed volume - Unger et al (1995)
•Index 3
1.CSTR with temperature control - McLellan (1994)
2.The pendulum model - Brenan et al (1989)
que pode ser selecionado como most ra a Figura A.7. Ao selecionar um exemplo, ativa-
se a tela de opções (Figura A.8).
O código ACIG possui o mesmo banco de dados do código INDEX, mas não
apresenta a tela de opções. A tela de seleção de exemplos do código OTIMA é mostrada
na Figura A.9.
O código EVENTS apresenta como banco de dados os exemplos descritos em
(XU, 2001) e podem ser selecionados como mostra a Figura A.10. Ao selecionar um
exemplo, dos códigos ACIG, OTIMA e EVENTS, ativa-se a tela de resultados (Figura
A.11).
Para se fazer a geração de resultados, basta clicar no botão Run. Para selecionar
outro exemplo, basta selecionar a opção Restart e para finalizar a interface, seleciona-se
o botão Exit.
88
Figura A.7: Tela de seleção - INDEX.
Figura A.8: Tela de Opções.
89
Figura A.9: Tela de seleção - OTIMA.
Figura A.10: Tela de seleção - EVENTS.
90
Figura A.11: Tela de geração de resultados.
APÊNDICE B
DIRCOL
B.1 Aspectos Gerais
O código DIRCOL 2.1 foi desenvolvido por Stryk (1999) e utiliza um método direto
de colocação para a discretização das variáveis de controle e de estado do problema de
controle ótimo de dimensão infinita, que é transformado numa seqüência de problemas
de otimização sujeito a restrições não lineares de dimensão finita (Problema de Pro-
gramação Não Linear). Este código não exige a aplicação da teoria de Controle Ótimo
ou o desenvolvimento das equações diferenciais adjuntas para a sua aplicação. O pro-
blema de Programação Não Linear (NLP) resultante pode ser resolvido pelo mé todo
de Programação Quadrática Seqüencial implementado no có digo NPSOL, aplicável a
sistemas densos, ou pelo código SNOPT aplicado a sistemas esparsos. O DIRCOL 2.1
pode tratar sistemas multifásicos e problemas com descontinuidades no lado direito das
equações diferenciais, além de fornecer estimativas confiáveis das variáveis adjuntas e
das funções multiplicadoras das restrições de estado (STRYK, 1999).
A Figura (B.1) mostra esquematicamente a estrutura do código DIRCOL 2.1 . A
subrotina USER.F é a subrotina que deve ser fornecida pelo usuário, específica para o
problema tratado.
A seguir são descritas as subrotinas que compõem o código DIRCOL e os arquivos
de entrada e saída de dados.
B.1. Aspectos Gerais 92
Figura B.1: Estrutura do DIRCOL 2.1.
B.1.1 Subrotinas de entrada do DIRCOL 2.1
DIRCOM Esta subrotina é chamada primeiro pelo programa principal.
USROBJ Esta subrotina informa a função objetivo do problema de controle ótimo
reescrita, se necessário, na forma de Mayer.
USRDEQ Esta subrotina calcula o lado direito das equações diferenciais.
USRNBC Calcula as condições no tempo inicial, os pontos de mudança caso exis-
tentes e o tempo final.
USRNIC Calcula as funções das restrições de desigualdade não lineares no tempo.
USRNEC Calcula as funções das restrições de igualdade não lineares no tempo.
USRSTV Informa os valores iniciais para o estado, controle e eventos.
USRREV Informa a existência de um sistema de referência para as variáveis de
controle e estado.
USRREA Informa a existência de um sistema de referência para as variáveis adjuntas.
B.2. Arquivo de entrada DATLIM 93
B.1.2 Arquivos de saída padrões do DIRCOL 2.1
•DATRES: dados gerais da otimização.
•DATSKA: constantes das transformações lineares usadas para o escalonamento
interno de todas as variáveis e funções.
•DATGIT (final): posições normalizadas dos p ontos interiores da malha de dis-
cretização utilizados, definidos para cada fase no final da otimização.
•GDATX: variáveis de estado discretizadas calculadas (funções cúbicas definidas
por partes).
•GDATU: variáveis de controle discretizadas calculadas (funções lineares defini-
das por partes).
•GDATD: erros locais calculados para verificação de precisão.
B.1.3 Arquivos de saída opcionais do DIRCOL 2.1
•GDATL: valores estimados das variáveis adjuntas (funções lineares definidas por
partes).
•GDATM: estimativas das funções multiplicadoras das restrições (funções linea-
res definidas por partes).
•GDATS: estimativas dos novos pontos de mudanças ou eventos.
B.2 Arquivo de entrada DATLIM
Nesta seção é apresentado o arquivo de entrada DATLIM para o Caso 3 (Reator semi-
batelada isotérmico) descrito no Capítulo 4
B.2. Arquivo de entrada DATLIM 94
************************************************************************
* file DATLIM *
* (prescribed values at initial time, final time and switching points, *
* lower and upper bounds for all variables X, U, P, E) *
************************************************************************
*
* the NX values of X(1) through X(NX) at E(1)=T0, E(M)=TF are
1 , 0, 2.D0
1 , 0, 0.63D0
1 , 0, 0.7D0
*
* the LU values of U(1) through U(LU) at E(1)=T0, E(M)=TF are
0 , 0
0 , 1
*
* 1. switching point E(2):
* -----------------------
* the NX values of X(1) through X(NX) at the switching point are
0 , 1
0 , 1
0 , 1
* the LU values of U(1) through U(LU) at the switching point are
0 , 0
0 , 0
* the lower and upper bounds of the events E(2),...,E(M)=TF are
* -------------------------------------------------------------
* MIN , MAX
10.E+0 , 15.E+0
18E+0 , 21.E+0
*
* 1. phase:
* ---------
B.3. Arquivo de entrada DATDIM 95
* the NX lower and upper bounds of the state variables X are
* X(I)MIN , X(I)MAX
0.0E+0 , +2.0E+0
0.63E+0 , +0.63E+0
0.0E+0 , +2.0E+0
*
* the LU lower and upper bounds of the control variables U are
* U(K)MIN , U(K)MAX
+0.E+0 , +0.035E+0
0.0E+0 , +1.0E+0
*
* 2. Phase:
* ---------
* the NX lower and upper bounds of the state variables X are
* X(I)MIN , X(I)MAX
0.0E+0 , +2.0E+0
0.0E+0 , +2.0E+0
1.0E+0 , +1.0E+0
*
* the LU lower and upper bounds of the control variables U are
* U(K)MIN , U(K)MAX
+0.0E+0 , +0.0E+0
0.0E+0 , +2.0E+0
*
*23456789012345678901234567890123456789012345678901234567890123456789012
B.3 Arquivo de entrada DATDIM
Nesta seção é apresentado o arquivo de entrada DATDIM para o Caso 3 (Reator
semi-batelada isotérmico) descrito no Capítulo 4
B.3. Arquivo de entrada DATDIM 96
************************************************************************
* file DATDIM *
* (dimensions of the parameterized optimal control problem) *
************************************************************************
* NAME of the OPTIMAL CONTROL PROBLEM
* -----------------------------------
*2345678901234567890123456789* (<-- max. length of name)
Concentration
*
* iAction:
* -------
* - OPTIMIZATION using NPSOL .................................... (0)
* - a check of all dimensions of feasibility .................... (1)
* - a check of subroutines & computation of starting trajectory . (2)
* or computation of a FEASIBLE TRAJECTORY by
* - objective min-max1 / use NPOPT .............................. (3)
* - objective min-max1 / use NPSOL .............................. (4)
* - objective min-max2 / use NPSOL .............................. (5)
* or actions involving SNOPT:
* - OPTIMIZATION using NPOPT ............................. (6)
* - OPTIMIZATION using SNOPT (dense Jacobian)............. (7)
* - OPTIMIZATION using SNOPT (sparse Jacobian)............ (8)
* - FEASIBLE TRAJECTORY using SNOPT (sparse Jacobian)............ (9)
*
* iAction, MajItL = ?,?
*
0, -10
* Optional SQP-Parameters:
* -----------------------
* Optimality Tolerance EPSOPT = ?
B.3. Arquivo de entrada DATDIM 97
1.0E-6
*
* Nonlinear Feasibility Tolerance EPSNFT = ?
1.0E-6
*
* Major Print Level (0, 5 or 10) = ?
5
*
* which SCALINGS and DIFFERENCE APPROXIMATIONS are to be used:
* -----------------------------------------------------------
* iScale:
* ------
* - automatic scaling (but for X, U, E in each phase the same) (0)
* - read scalings from file ’DATSKA’ (1)
* - use no scaling (2)
* - automatic scaling (X, U, E in each phase different) (3)
* - automatic scaling (X, U in each phase the same, but E different)(4)
*
* iDiff:
* -----
* - forward difference approximations of DIRCOL (default) (0)
* - internal difference approximation of NPSOL or SNOPT (-1)
*
* iScale, iDiff = ?,?
0, 0
*
* NUMBER of STATE VARIABLES ( NX ),
* ------ of CONTROL VARIABLES ( LU ),
* of CONTROL PARAMETERS ( LP ),
B.3. Arquivo de entrada DATDIM 98
* NX, LU, LP = ?
3, 2, 0
*
* NUMBER of PHASES M1 = ?
* ------------------------
2
*
* NUMBERS of NONLINEAR IMPLICIT BOUNDARY CONSTRAINTS
* --------------------------------------------------
* NRNLN(1)
* ...
* NRNLN(M1)
0
0
* NUMBERS of NONLINEAR INEQUALITY and EQUALITY CONSTRAINTS
* --------------------------------------------------------
* in phases 1 through M1:
* NGNLN(1) , NHNLN(1)
* ...
* NGNLN(M1), NHNLN(M1)
0, 1
0, 1
*
* NUMBER of GRID POINTS in phases 1 through M1 ( NG(K) >= 3 ):
* ------------------------------------------------------------
* NG(1)
* ...
* NG(M1)
10
10
*
B.3. Arquivo de entrada DATDIM 99
* GRID POINTs parameters:
* iStartGrid | iOptGrid (during optimization):
* ---------- | --------
* (starting positions): | - fixed grid points (0)
* - equidistant (0) | - movable (collocation error) (1)
* - as in file DATGIT (1) | - movable (variation) (2)
* - as Cebysev points (2) | - movable (no add. eq. cons.) (3)
*
* iStartGrid, iOptGrid = ?,?
0, 0
*
* STARTING VALUES of X(t), U(t), P, and E:
* --------------------------------------------
* - as specified in subroutine USRSTV (0)
* - as in files GDATX, GDATU (unchanged number of phases) (1)
* - X, U, P as in files GDATX, GDATU and
* E as specified in USRSTV (changed number of phases) (2)
0
*
* ESTIMATES of the ADJOINT VARIABLES and of
* -----------------------------------------
* the SWITCHING STRUCTURES of state and control constraints
* - are NOT required (0)
* - are required (1)
1
* NAMES of the NX state variables:
* ---------------------------------
* X(1)_Name
* ...
* X(NX)_Name
*2345678901234* (<-- max. length of name)
B.4. Subrotina USER.f 100
Ca(t)
Cb(t)
V(t)
*
* NAMES of the LU control variables:
* -----------------------------------
* U(1)_Name
* ...
* U(LU)_Name
*2345678901234* (<-- max. length of name)
U(t)
Cc(t)
* the I-th STATE VARIABLE (I = 1,.., NX) is an UNCONSTRAINED ANGLE
* and varies only in [ -PI, PI [ : 1 (if yes) or 0 (if not)
*
0
0
0
* the K-th control variable (K = 1,.., LU) is an UNCONSTRAINED ANGLE
* and varies only in [ -PI, PI [ : 1 (if yes) or 0 (if not)
*
0
0
B.4 Subrotina USER.f
Nesta seção é apresentada a subrotina USER.f para o Caso 3 (Reator semi-batelada
isotérmico) descrito no Capítulo 4. Esta subrotina fornece informações do problema
que será resolvido. Algumas subrotinas não aparecem nesta descrição pois não foram
necessárias para este estudo, devendo mantê-las em sua forma original.
B.4. Subrotina USER.f 101
C - file user.f of problem: SrPB03
SUBROUTINE DIRCOM()
C***********************************************************************
C Entrada dos parâmetros do modelo
C***********************************************************************
IMPLICIT NONE
C--------BEGIN---PROBLEM------------------------------------------------
C**** REAL
DOUBLE PRECISION
+ k,Delta,rho,cp,cbin,Tmax,Vmax,ncdes,cbmax
COMMON /USRCOM/ k,Delta,rho,cp,cbin,Tmax,Vmax,ncdes,cbmax
k = 0.0482D0
Delta = -60000.D0
rho = 900.D0
cp = 4.2D0
cbin = 2.D0
Tmax = 80.D0
Vmax = 1.D0
ncdes = 0.6D0
cbmax = 0.63D0
C---------END----PROBLEM------------------------------------------------
RETURN
END
SUBROUTINE USRSTV( IPHASE, NX, LU, TAU, X, U, IFAIL)
C***********************************************************************
C Estimativas iniciais para as variáveis de controle e estado, para os
C eventos, tempo inicial e tempo final. Os valores utilizados para o
C tempo inicial, evento e tempo final foram 0, 12, 20, respectivamente.
C***********************************************************************
IMPLICIT NONE
C
INTEGER IPHASE, NX, LU, IFAIL
B.4. Subrotina USER.f 102
C**** REAL
DOUBLE PRECISION
+ T,TAU, X(NX), U(LU)
DOUBLE PRECISION
+ k,Delta,rho,cp,cbin,Tmax,Vmax,ncdes,cbmax
c
COMMON /USRCOM/ k,Delta,rho,cp,cbin,Tmax,Vmax,ncdes,cbmax
C
C--------BEGIN---PROBLEM------------------------------------------------
C
IF (IPHASE .GT. 0) THEN
C
C ------ Initial estimates of X and U at TAU
C
T =TAU
X(1) = 1.2D0
X(2) = 0.5D0
X(3) = 1.D0
U(1) = 0.015D0
U(2) = 0.5D0
C
ELSE IF (IPHASE .EQ. 0) THEN
C
C ------ Initial estimates of the events E
C
U(1) = 0.D0
U(2) = 12.D0
U(3) = 20.D0
ELSE IF (IPHASE .LT. 0) THEN
C
C---------END----PROBLEM------------------------------------------------
RETURN
C --- End of subroutine USRSTV
END
B.4. Subrotina USER.f 103
SUBROUTINE USROBJ( NR, NX, LU, LP, XL, UL, P, ENR, FOBJ,
& XR, UR, TF, IFAIL)
C***********************************************************************
C* Função Objetivo do Problema. Neste caso o objetivo é minimizar o tempo
C final (problema de tempo final livre)
C***********************************************************************
C
IMPLICIT NONE
C
INTEGER NX, LU, LP, NR, IFAIL
DOUBLE PRECISION
C**** REAL
+ ENR, XL(NX), UL(LU), P(LP), FOBJ, XR(NX), UR(LU), TF
C
C--------BEGIN---PROBLEM------------------------------------------------
C
IF (NR .EQ. 2) THEN
FOBJ = TF
END IF
IF (NR .EQ. 1) THEN
FOBJ = 0.
END IF
C---------END----PROBLEM------------------------------------------------
RETURN
C --- End of subroutine USROBJ
END
B.4. Subrotina USER.f 104
SUBROUTINE USRDEQ( IPHASE, NX, LU, LP, X, U, P, T, F, IFAIL)
C***********************************************************************
C* Lado direito das equações diferenciais
C***********************************************************************
C
IMPLICIT NONE
C
INTEGER IPHASE, NX, LU, LP, IFAIL
C**** REAL
DOUBLE PRECISION
+ X(NX), U(LU), P(LP), T, F(NX)
C
DOUBLE PRECISION
+ k,Delta,rho,cp,cbin,Tmax,Vmax,ncdes,cbmax
c
COMMON /USRCOM/ k,Delta,rho,cp,cbin,Tmax,Vmax,ncdes,cbmax
C
C--------BEGIN---PROBLEM------------------------------------------------
C
F(1) = -k*X(1)*X(2)-U(1)/X(3)*X(1)
F(2) = -k*X(1)*X(2)+U(1)/X(3)*(cbin-X(2))
F(3) = U(1)
C
C---------END----PROBLEM------------------------------------------------
C
RETURN
C --- End of subroutine USRDEQ
END
B.4. Subrotina USER.f 105
SUBROUTINE USRNBC( IKIND, NRNLN, NX, LU, LP, XL, UL, P, EL, RB,
1 XR, UR, ER, IFAIL)
C***********************************************************************
C* Condições de continuidade
C***********************************************************************
C
IMPLICIT NONE
C
INTEGER IKIND, NRNLN, NX, LU, LP, IFAIL
C**** REAL
DOUBLE PRECISION
+ XL(NX), XR(NX), UL(LU), UR(LU), P(LP), EL, ER, RB(NRNLN)
C
DOUBLE PRECISION
+ k,Delta,rho,cp,cbin,Tmax,Vmax,ncdes,cbmax
c
COMMON /USRCOM/ k,Delta,rho,cp,cbin,Tmax,Vmax,ncdes,cbmax
c
C--------BEGIN---PROBLEM------------------------------------------------
C
IF (IKIND .EQ. -2) THEN
XR(1) = XL(1)
XR(2) = XL(2)
XR(3) = XL(3)
END IF
IF (IKIND .EQ. -1) THEN
UR(2) = 0.6D0
END IF
C
C---------END----PROBLEM------------------------------------------------
RETURN
C --- End of subroutine USRNBC
END
B.4. Subrotina USER.f 106
SUBROUTINE USRNEC(IPHASE, NHNLN, NEEDH, NX, LU, LP, X, U, P, T, H,
& IFAIL)
C***********************************************************************
C* Restrições de igualdade não lineares
C***********************************************************************
IMPLICIT NONE
C
INTEGER NHNLN, NX, LU, LP, IPHASE, NEEDH(NHNLN), IFAIL
C**** REAL
DOUBLE PRECISION
+ T, X(NX), U(LU), P(LP), H(NHNLN)
C
DOUBLE PRECISION
+ k,Delta,rho,cp,cbin,Tmax,Vmax,ncdes,cbmax
c
COMMON /USRCOM/ k,Delta,rho,cp,cbin,Tmax,Vmax,ncdes,cbmax
C
C--------BEGIN---PROBLEM------------------------------------------------
H(1) = U(2)-1/X(3)*(1.4D0-X(1)*X(3))
C---------END----PROBLEM------------------------------------------------
RETURN
C --- End of subroutine USRNEC
END
Livros Grátis
( http://www.livrosgratis.com.br )
Milhares de Livros para Download:
Baixar livros de Administração
Baixar livros de Agronomia
Baixar livros de Arquitetura
Baixar livros de Artes
Baixar livros de Astronomia
Baixar livros de Biologia Geral
Baixar livros de Ciência da Computação
Baixar livros de Ciência da Informação
Baixar livros de Ciência Política
Baixar livros de Ciências da Saúde
Baixar livros de Comunicação
Baixar livros do Conselho Nacional de Educação - CNE
Baixar livros de Defesa civil
Baixar livros de Direito
Baixar livros de Direitos humanos
Baixar livros de Economia
Baixar livros de Economia Doméstica
Baixar livros de Educação
Baixar livros de Educação - Trânsito
Baixar livros de Educação Física
Baixar livros de Engenharia Aeroespacial
Baixar livros de Farmácia
Baixar livros de Filosofia
Baixar livros de Física
Baixar livros de Geociências
Baixar livros de Geografia
Baixar livros de História
Baixar livros de Línguas
Baixar livros de Literatura
Baixar livros de Literatura de Cordel
Baixar livros de Literatura Infantil
Baixar livros de Matemática
Baixar livros de Medicina
Baixar livros de Medicina Veterinária
Baixar livros de Meio Ambiente
Baixar livros de Meteorologia
Baixar Monografias e TCC
Baixar livros Multidisciplinar
Baixar livros de Música
Baixar livros de Psicologia
Baixar livros de Química
Baixar livros de Saúde Coletiva
Baixar livros de Serviço Social
Baixar livros de Sociologia
Baixar livros de Teologia
Baixar livros de Trabalho
Baixar livros de Turismo
