Embed Size (px)
Citation preview
KARLA CRISTIANE ARSIE
METODO DE EULER E METODO PSEUDOESPECTRAL
USANDO PONTOS LEGENDRE GAUSS RADAU PARA
UMA CLASSE DE PROBLEMAS DE CONTROLE OTIMO
Curitiba
2013
KARLA CRISTIANE ARSIE
METODO DE EULER E METODO PSEUDOESPECTRAL
USANDO PONTOS LEGENDRE GAUSS RADAU PARA
UMA CLASSE DE PROBLEMAS DE CONTROLE OTIMO
Dissertacao apresentada ao Programa de Pos-
Graduacao em Matematica Aplicada da Uni-
versidade Federal do Parana, como requisito
parcial a obtencao do grau de Mestre em Ma-
tematica Aplicada.
Orientadora: Dr.a Elizabeth Wegner Karas.
Coorientador: Dr. Miguel A. Dumett Canales.
Curitiba
2013
iii
iv
v
Agradecimentos
Agradeco primeiramente a minha famılia, pelo apoio,
compreensao, ajuda e por todo carinho ao longo deste percurso.
Amo voces.
Ao meu noivo Rodrigo, que esteve sempre ao meu lado, que me
apoiou, incentivou, ajudou e participou comigo dos momentos de
angustias e alegrias. Neoqeav.
Aos meus amigos, especialmente Izabela, Janaina, Keilla,
Leonardo e Priscila, pela ajuda e pelo carinho.
Agradeco aos meus orientadores Elizabeth e Miguel, pelas
orientacoes, conversas, por estarem sempre dispostos a me ajudar.
Obrigada pelo esforco, tempo dedicado e por sempre acreditarem
no meu trabalho.
Aos membros da banca, professores Dr. Antonio Carlos Gardel
Leitao, Dr. Saulo Pomponet Oliveira e Dr.a Ailin Ruiz Zarate
Fabregas obrigada por aceitarem fazer parte deste momento e dar
suas contribuicoes valiosas.
Ao Programa de Pos-Graduacao em Matematica Aplicada da
UFPR pela oportunidade e formacao de qualidade propiciada.
A CAPES pelo apoio financeiro.
E finalmente, mas nao menos importante, a Deus pela capacitacao
concedida, pela forca espiritual e por seguir ao meu lado.
vi
“ A mente que se abre a uma nova ideia jamais
voltara ao seu tamanho original.”
Albert Einstein
vii
Resumo
O objetivo deste trabalho e discutir o metodo pseudoespectral com pontos
de colocacao de Legendre-Gauss-Radau (LGR) apresentado em [22], para
determinar solucoes numericas de algumas classes de problemas de con-
trole otimo. Nesta dissertacao revisa-se [22], e se deriva a discretizacao do
metodo pseudoespectral LGR, de problemas de controle otimo (sem res-
tricoes nas variaveis de estado e de controle) utilizando notacao tensorial.
Adicionalmente se derivam as condicoes de otimalidade de Karush-Kuhn-
Tucker (KKT) associadas ao problema.
Para avaliar a precisao do metodo em problemas de controle otimo es-
pecıficos, e necessario conhecer a solucao exata dos problemas escolhidos.
Procurando replicar os resultados em [22], trabalhou-se num primeiro exem-
plo com um Problema de Bolza (tipo LQR) sem restricoes nas variaveis de
estado e de controle. Se apresenta uma derivacao detalhada da solucao exata
deste problema quadratico, utilizando o Princıpio do Maximo de Pontrya-
gin. O problema de minimizacao resultante foi resolvido atraves da rotina
quadprog do MATLAB. A precisao do metodo pseudoespectral LGR e com-
parada, com bons resultados, com o metodo de Euler (aplicado ao problema
de otimizacao quadratico produto da discretizacao por Euler do problema
de Bolza tipo LQR original).
Para evidenciar que o metodo pseudoespectral LGR de discretizacao pode
ser aplicado a problemas de controle otimo com restricoes nas variaveis de
controle e de estado (o que nao e abordado em [22]), dois exemplos adici-
onais, apresentados em [31], sao discutidos nesta dissertacao. No segundo
exemplo a funcao custo e nao quadratica e a rotina fmincon do MATLAB e
encarregada de fazer o trabalho de otimizacao a partir das equacoes discre-
tizadas pelo metodo pseudoespectral LGR. No terceiro exemplo, o problema
de otimizcao foi resolvido pela rotina quadprog. Existem poucos problemas
nao-lineares de controle otimo (com restricoes) cujas solucoes exatas sao co-
nhecidas. Usualmente argumentos de convexidade e outros, sao necessarios
para encontrar as solucoes exatas.
Palavras-chave: Controle otimo; problema tipo Bolza; metodo pseudoes-
pectral; pontos de colocacao LGR; metodo de Euler.
viii
Abstract
The goal of this work is to present the details of the Legendre-Gauss-Radau
(LGR) pseudoespectral numerical method. This method was published in
[22] and it is utilized to find numerical solutions of certain classes of optimal
control problem. In this dissertation, we review [22], and derive the discre-
tization of the LGR pseudoespectral method of optimal control problems
(without restrictions in the state and control variables) using tensorial no-
tation. In adition, the associated Karush-Kuhn-Tucker (KKT) optimality
conditions are derived. To assess the accuracy of the LGR method in specific
optimal control problems, it is necessary to know the exact solution of the
selected problems. Aiming to replicate the results in [22], we worked initially
a Bolza problem (LQR type) without restrictions in the state and control
variables. The process of finding the exact solution of this quadratic pro-
blem is derived in detail, utilizing the Pontryagin Maximum Principle. The
LGR pseudoespectral discretization technique was successful when applied
to the Bolza problem without restrictions. The corresponding MATLAB
code is included in the Appendix, and utilizes mainly the quadprog routine.
The accuracy of the LGR pseudoespectral method is compared successfully
against Euler method (applied to the quadratic optimization problem which
is obtained from the forward Euler discretization of the LQR type original
Bolza problem). To point out that the LGR pseudoespectral discretization
method could be applied to optimal control problems with restrictions in the
state and control variables (something not attempted in [22]), and outside of
the Bolza problem class, two additional examples from [31] are presented in
this dissertation. There are few non-linear optimal control problems (with
restrictions) whose exact solutions are known. Usually, convexity arguments
and others, are necessary to find the exact solutions.
Palavras-chave: Optimal control; Bolza problem; pseudoespectral method;
LGR collocation points; Euler method.
ix
Sumario
Introducao 1
1 O problema de Controle Otimo 4
1.1 Princıpio do Maximo de Pontryagin . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Metodo Pseudoespectral usando Pontos Legendre Gauss Radau . . . . . . 6
1.2.1 Fundamentacao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2.2 Discretizacao do problema . . . . . . . . . . . . . . . . . . . . . . . 10
1.2.3 Condicoes de Karush-Kuhn-Tucker . . . . . . . . . . . . . . . . . . 12
1.2.4 Sistema Adjunto Transformado . . . . . . . . . . . . . . . . . . . . 15
1.2.5 Formulacao Integral . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.3 Metodo de Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
1.4 Outros metodos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
1.4.1 Comparacao com o metodo de Kameswaran e Biegler . . . . . . . . 21
1.4.2 Comparacao com o metodo de Fahroo e Ross . . . . . . . . . . . . 22
1.5 Problema de Bolza . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2 Exemplos 24
2.1 Exemplo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.1.1 Solucao Exata . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.1.2 Solucao aproximada pelo metodo pseudoespectral . . . . . . . . . . 28
2.1.3 Solucao aproximada pelo metodo de Euler . . . . . . . . . . . . . . 30
2.2 Exemplo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.2.1 Solucao aproximada pelo metodo pseudoespectral . . . . . . . . . . 35
2.3 Exemplo 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.3.1 Solucao aproximada pelo metodo pseudoespectral . . . . . . . . . . 38
Conclusao 41
Apendices 42
A Revisao de Conceitos 42
A.1 Matriz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
A.1.1 Produto direto de Matrizes . . . . . . . . . . . . . . . . . . . . . . . 42
x
A.1.2 Exponencial de matriz . . . . . . . . . . . . . . . . . . . . . . . . . 43
A.2 Solucao de um sistema de EDOs . . . . . . . . . . . . . . . . . . . . . . . . 45
A.3 Topicos de otimizacao contınua . . . . . . . . . . . . . . . . . . . . . . . . 46
A.3.1 Metodo de gradientes conjugados . . . . . . . . . . . . . . . . . . . 47
A.3.2 O Teorema de Karush-Kuhn-Tucker . . . . . . . . . . . . . . . . . . 49
B Codigos em Matlab 52
Referencias Bibliograficas 67
xi
Introducao
O objetivo da teoria de controle otimo e determinar o controle, que fara com que
um sistema satisfaca um conjunto de restricoes fısicas e ao mesmo tempo minimize algum
criterio de desempenho. Normalmente, os problemas de controle otimo tem restricoes
nas variaveis de estado e de controle. As solucoes para problemas de controle otimo
(sem restricoes) podem ser encontradas, por exemplo, atraves da aplicacao do calculo
variacional ([11, 20, 31]) e da aplicacao do Princıpio do Maximo de Prontryagin ([11, 18]).
Em [11] se faz uma analogia entre problemas de controle otimo e os problemas do calculo
de variacoes, mostra-se tambem uma estrategia de como o Princıpio do Maximo pode, em
alguns casos, ser utilizado para efetivamente determinar uma solucao de um determinado
problema de controle otimo. Porem, em [31, Cap. 10], se menciona que em geral nao e
uma boa ideia trabalhar com a estrategia antes mencionada, pois ela se baseia em definir
U = Y ′ e trabalhar com a variavel estado estendida W = (Y, U), nos problemas de
calculo de variacoes para tentar encontrar as solucoes de problemas de controle otimo. A
dificuldade principal esta no fato que Y ′ e conhecida explicitamente, mas U ′ nao e. Varias
tecnicas foram desenvolvidas para resolver os problemas de controle otimo. Estes metodos
numericos sao divididos em duas categorias: metodos indiretos e metodos diretos [3, 17].
Os metodos indiretos consistem em converter o problema de otimizacao em um
sistema de equacoes algebrico-diferenciais de valor no contorno a partir da aplicacao do
Princıpio do Maximo de Pontryagin. A principal vantagem desse metodo e a alta precisao e
garantias de solucoes que satisfazem as condicoes necessarias de otimalidade. No entanto,
ha uma desvantagem significativa. As condicoes necessarias de otimalidade devem ser
encontradas analiticamente, e para a maioria dos problemas esta busca nao e trivial.
Os metodos diretos utilizam a parametrizacao na variavel de controle ou a para-
metrizacao nas variaveis de estado seguida da solucao do problema de programacao nao-
linear resultante. Nestes metodos polinomios podem ser usados para aproximar a equacoes
diferenciais em pontos de colocacao ([4, 24, 28]). Estas tecnicas baseiam-se em metodos
pseudoespectrais que consistem em encontrar solucoes numericas de equacoes diferenciais.
Metodos pseudoespectrais sao uma classe de metodos de colocacao direta, onde o proble-
ma de controle otimo e transcrito para um problema de programacao nao-linear atraves
de uma parametrizacao do estado e do controle utilizando polinomios de Legendre ou
Chebyshev em conjunto com pontos de colocacao. Tal problema de programacao nao-
linear pode entao ser resolvido com uso de algoritmos de Otimizacao para obter solucoes
Introducao 2
aproximadas localmente otimas.
Neste trabalho aborda-se um metodo pseudoespectral que consiste na parame-
trizacao do estado e do controle utilizando polinomios de Legendre em conjunto com
a discretizacao utilizando pontos de colocacao LGR. O metodo fornece, tambem, uma
maneira precisa de encontrar as condicoes de otimalidade de KKT associadas ao pro-
blema. Mostra-se que o metodo apresentado esta relacionado a metodo de integracao
global implıcita. Isto e feito provando que a inversa da matriz de diferenciacao obtida
no metodo pseudoespectral coincide com a matriz associada a um metodo de integracao
global implıcito [16].
Este trabalho esta organizado da seguinte forma. No Capıtulo 1, se apresenta bre-
vemente o que e um problema de controle otimo (sem restricoes) e se enuncia o resultado
fundamental da teoria de controle otimo, o teorema do Princıpio do Maximo de Pon-
tryagin. Em seguida, se apresenta o metodo pseudoespectral LGR, seguindo a sequencia
estabelecida em [22], mas com demonstracoes diferenciadas das encontradas nessa re-
ferencia. Adicionalmente se menciona o metodo de Euler e se apresenta o problema de
Bolza.
O Capıtulo 2 inicia especificando o tipo LQR do problema de Bolza e encontrando
a solucao exata desse problema utilizando o Princıpio do Maximo de Pontryagin. Logo, a
discretizacao pseudoespectral LGR e utilizada para encontrar uma solucao numerica do
problema de Bolza LQR. A precisao do metodo e comparada com a precisao do metodo
de Euler aplicado ao mesmo problema de Bolza sem restricoes. Os codigos MATLAB
necessarios para construir os pontos de colocacao LGR e para calcular os pesos de qua-
dratura da discretizacao pseudoespectral LGR, assim como tambem a chamada da rotina
QUADPROG, para encontrar explicitamente, as solucoes numericas do estado e do con-
trole do problema de Bolza, sao apresentados no Apendice da dissertacao.
Adicionalmente, sao incluıdos dois problemas de controle otimo com restricoes
nao-lineares, para testar a discretizacao pseudoespectral LGR (isto nao e feito em [22])
Ambos os problemas, com suas respectivas solucoes exatas, sao apresentados em [31, Cap.
10]. O metodo pseudoespectral LGR e bem sucedido em encontrar uma solucao numerica.
Os codigos em MATLAB correspondentes sao dados no Apendice da dissertacao e estarao
disponıveis no site da pos-graduacao em Matematica Aplicada.
Finalmente, nos Apendices apresentam-se alguns resultados basicos da teoria
abordada no trabalho, bem como os codigos implementados em Matlab para resolver
numericamente os exemplos discutidos.
Introducao 3
Notacao.
xT transposto do vetor x
‖.‖∞ norma infinito
‖.‖ norma 2
x, x′ derivada de x
xj(τ) derivada da j-esima componente de x em relacao a τ
〈x, y〉 produto interno de x por y
ej vetor canonico com 1 na componente j e as demais componentes nulas
AT transposta da matriz A
Ak k-esima linha da matriz A
diag(α1, . . . , αn) matriz diagonal quadrada cujos elementos da diagonal sao α1, . . . , αn
AB multiplicacao da matriz A pela matriz B
∇f gradiente da funcao f : IRn → IR. A mesma notacao sera adotada
para indicar a matriz m× n cuja i-esima linha e ∇fi, quando
f : IRn → IRm
∇φ matriz m× n cujo elemento (i, j) e (∇φ(X))ij = ∂φ(X)∂Xij
.
Aqui φ : IRm×n → IR e X ∈ IRm×n
C(I) Espaco das funcoes contınuas em I.
C(I) Espaco das funcoes contınuas por partes do intervalo I.
C1(I) Espaco das funcoes continuamente diferenciaveis em I.
C1(I) Espaco das funcoes continuamente diferenciaveis por partes no intervalo I.
Capıtulo 1
O problema de Controle Otimo
Um problema de controle otimo consiste em minimizar um certo funcional que
depende de dois tipos de variaveis: variaveis de estado ou proprias (x ∈ IRn) e variaveis
de controle (u ∈ IRm). As duas variaveis estao parametrizadas por uma variavel real
independente t, usualmente associada ao tempo, isto e, x = x(t), u = u(t). Formalmente
x ∈ C1 e uma funcao x : IR → IRn diferenciavel por pedacos e u ∈ C e uma funcao
u : IR→ IRm contınua por pedacos em cada intervalo de tempo. A variavel de estado x(t)
obedece a uma dinamica dada pela seguinte equacao diferencial, denominada equacao de
estado
x(t) = f(x(t), u(t)), (1.1)
onde f : IRn×m → IRn e uma funcao contınua. A variavel de estado pode estar ainda
sujeita a condicoes x(t0) = x0 e x(tf ) = xf . A variavel de controle u(t) e uma funcao que
pertence a um certo espaco de funcoes U , chamado o conjunto de controles admissıveis.
Por exemplo, em [18] considera-se U = u : [0,∞) −→ IRm | u(.) e mensuravel.No caso que (1.1) nao depende de u, entao (1.1) vira uma equacao diferencial. E
conhecido que sob a condicao de continuidade de f , o teorema de Peano [30] garante a
existencia de solucoes de (1.1). Adicionalmente, se a funcao f e localmente Lipschitziana
em x, o teorema de Picard[30] garante a existencia de uma solucao unica de (1.1). Quando
a funcao u(t) e constante, em [30, Cap.2] encontram-se condicoes similares as de Peano e
Picard, para garantir existencia e unicidade de solucoes de (1.1).
Usualmente, o funcional a minimizar em um problema de controle otimo e com-
posto de dois termos: o primeiro deles depende do estado final x(tf ) e o segundo e da
forma ∫ tf
t0
r(x(t), u(t))dt.
Adicionalmente, sao impostas restricoes por desigualdade nas variaveis de estado
e de controle. Por exemplo, x(t) ≥ 0 e |u(t)| ≤ 1, para t0 ≤ t ≤ tf .
Condicoes gerais para existencia e unicidade de um problema de controle otimo
podem ser encontradas em [31]. Em geral, estas condicoes sao do tipo f ∈ C1, r ∈ C1, x ∈
4
5
C1, u ∈ C0 por partes. Neste trabalho assumimos a existencia de uma unica solucao. Nos
exemplos a serem considerados no proximo capıtulo, a solucao exata e conhecida.
Em geral, encontrar analiticamente uma solucao exata de um problema de contro-
le otimo com restricoes e muito difıcil e metodos numericos sao necessarios. No entanto,
quando o problema de controle otimo nao tem restricoes, o Princıpio de Maximo de
Pontryagin e uma ferramenta muito util para procurar solucoes exatas.
Como em [22], considera-se o caso de problema de controle otimo sem restricoes
para apresentar a correspondente discretizacao pseudoespectral LGR.
Pode-se formular o problema de controle otimo descrito acima, segundo [11], por
minimizar∫ tft0r(x(t), u(t))dt
sujeito a u ∈ Ux(t) = f(x(t), u(t))
x(t0) = x0
x(tf ) = xf .
(1.2)
No problema acima, os instantes de tempo inicial e final sao dados (problema de
tempo fixo), porem, ha problemas em que o tempo final e uma das incognitas no problema
de controle otimo (problemas de tempo livre) ou as vezes o tempo tf e dado mas xf (o
estado final) e desconhecido (problema de tempo fixo com estado final variavel) [18]. Este
trabalho restringe-se nos problemas de tempo fixo. Para mais detalhes consultar [31] e
[11, Cap.10].
Neste capıtulo sera apresentado um dos principais resultados na teoria de con-
trole otimo, o Princıpio do Maximo de Pontryagin, que neste trabalho sera utilizado
para obter a solucao exata de um problema de controle otimo. Em seguida, para se
obter a solucao aproximada deste problema, apresenta-se um metodo pseudoespectral,
que consiste na parametrizacao do estado e do controle utilizando polinomios de Legendre
em conjunto com a discretizacao com pontos de colocacao LGR - Legendre Gauss Radau.
Sao discutidas as condicoes de otimalidade de Karush-Kuhn-Tucker (KKT) associadas ao
problema discretizado e mostra-se que este metodo pseudoespectral pode ser visto como
um metodo de integracao global implıcito. Por fim, apresenta-se o metodo de Euler, um
metodo alternativo, para obter uma solucao aproximada.
1.1 Princıpio do Maximo de Pontryagin
O princıpio de Pontryagin estabelece um conjunto de condicoes para que uma
curva seja solucao de um problema de controle otimo. Esse princıpio se expressa em
termos da funcao H : IRn × IRn × IRm → IR definida por
H(x, p, u) = 〈f(x, u), p〉+ r(x, u),
6
onde f : IRn × IRm → IRn, r : IRn × IRm → IR (as mesmas funcoes dadas no problema do
controle otimo) e p ∈ IRn. Essa funcao e conhecida como Hamiltoniano de Pontryagin.
O teorema abaixo, que e baseado em [11, Teorema 178] e [18, Teorema 4.3], e
apresentado de uma forma conveniente para que, na obtencao da solucao exata de um
problema de controle otimo, se possa utiliza-lo diretamente.
Teorema 1.1 - Princıpio do Maximo de Pontryagin Considere um sistema de con-
trole otimo, sendo u = ω(t) o controle e x = γ(t) a trajetoria correspondente. Se
(γ(t), ω(t)) e uma solucao otima entao existe uma funcao µ : [0, tf ] → IRn tal que, para
todo t ∈ [0, tf ] vale que:
• A curva p = µ(t) satisfaz
µi = −∂H∂xi
(γ(t), µ(t), ω(t)),
• A trajetoria x = γ(t) satisfaz
γi = −∂H∂pi
(γ(t), µ(t), ω(t)),
• O controle u = ω(t) maximiza o Hamiltoniano, ou seja,
maxu∈U
H(γ(t), µ(t), u) = H(γ(t), µ(t), ω(t)).
Demonstracao. [11, Apendice B].
Utiliza-se o princıpio do maximo para determinar a solucao de um problema.
1.2 Metodo Pseudoespectral usando Pontos Legen-
dre Gauss Radau
Nesta secao estuda-se a formulacao do metodo pseudoespectral usando pontos
LGR - Legendre Gauss Radau, proposto em [22] para resolver o problema de controle
otimo (1.2). Inicialmente o problema e discretizado, tendo como referencia [8]. Em
seguida, as condicoes de KKT sao estabelecidas para o problema discretizado, seguindo
uma abordagem diferente, mas equivalente, a apresentada em [22]. Enquanto em [22] faz-
se o uso do Lagrangeano, aqui opta-se pela forma discutida no Apendice A.3.2, ou seja,
encontram-se as condicoes de otimalidade relacionando o oposto do gradiente da funcao
objetivo com o gradiente das restricoes do problema e da condicao inicial com relacao a
cada uma de suas variaveis. Cabe ressaltar que ambas as abordagens sao equivalentes.
Por fim, mostra-se que a restricao do problema tem uma formulacao equivalente atraves
7
de integrais, mostrando que este metodo pseudoespectral pode ser visto como um metodo
de integral global implıcito.
1.2.1 Fundamentacao
Pontos de colocacao LGR
Um dos objetivos deste trabalho e discretizar um problema para que se possa
encontrar sua solucao aproximada atraves de um metodo de colocacao, que consiste em
encontrar a solucao numerica de equacoes diferenciais ordinarias, equacoes diferenciais
parciais e equacoes integrais. A ideia e escolher um espaco de dimensao finita de solucoes
canditadas, que sao obtidas, geralmente, atraves de polinomios aplicados em uma quanti-
dade de pontos no domınio, chamados de pontos de colocacao. A discretizacao proposta
por [22] utiliza pontos de colocacao Legendre-Gauss-Radau (LGR).
Alem do conjunto de pontos de colocacao LGR, outros dois conjuntos sao apresen-
tados em [22] e muito utilizados: Lengendre-Gauss(LG) e Legendre-Gauss-Lobatto(LGL).
Esses tres conjuntos de pontos sao obtidos a partir de raızes de um polinomio de Legendre
e\ou combinacoes lineares desses polinomios e suas derivadas.
O polinomio de Legendre de grau n, Pn, e dado por:
Pn(x) =
[n/2]∑k=0
(−1)k(2n− 2k)!
2nk!(n− k)!(n− 2k)!xn−2k, (1.3)
onde
[n/2] =
n
2se n e par
n− 1
2se n e ımpar.
Ha algumas relacoes entre os polinomios de Legendre Pn com valores de n con-
secutivos e tambem entre os Pn e suas derivadas. Estas relacoes de recorrencia permitem
obter um novo Pn ou relacionar derivadas, a partir de outros. Apresentam-se a seguir
algumas destas relacoes de recorrencia. Os detalhes podem ser consultados em [6]. A
primeira relacao, bastante utilizada, permite obter um polinomio a partir de outros dois
de menor ordem,
nPn(u) = (2n− 1)uPn−1(u)− (n− 1)Pn−2(u). (1.4)
A proxima, relaciona o polinomio e as derivadas de polinomios consecutivos
uP ′n(u)− P ′n−1 = nPn(u). (1.5)
Por fim, tem-se a relacao de polinomio com a derivada de outros dois
(2n+ 1)Pn(u) = P ′n+1(u)− P ′n−1(u). (1.6)
8
Por (1.3), P0(u) = 1 e P1(u) = u. A relacao (1.4) permite a obtencao rapida de
P2(u) = 32u2 − 1
2, a seguir P3(u) = 5
2u3 − 3
2u, e assim por diante. Geralmente, esse e o
processo usado em computacao, para se evitar os fatoriais de ordem alta, sendo que estes
tem grande custo computacional.
A figura abaixo mostra os polinomios de Legendre de graus 0, 1, 2 e 3.
−1.5 −1 −0.5 0 0.5 1 1.5−4
−3
−2
−1
0
1
2
3
4
t
Solu
çã
o
P0(u)
P1(u)
P2(u)
P3(u)
Figura 1.1: Polinomios de Legendre Pn(u).
O conjunto de N pontos de colocacao Legendre-Gauss-Radau (LGR), Lengendre-
Gauss(LG) e Legendre-Gauss-Lobatto(LGL) sao definidos no domınio [−1, 1] como:
LG : raızes de PN ;
LGR : raızes de PN−1 + PN ;
LGL : raızes de PN juntamente com os pontos− 1 e 1,
onde PN denota a derivada do polinomio de Legendre PN . Note que os pontos de LG
nao incluem nenhum dos extremos do intervalo [−1, 1], os pontos LGR incluem um dos
seus extremos, e os pontos LGL incluem ambas as extremidades do intervalo. Outra dife-
renca e que os pontos de LGR, ao contrario dos outros dois conjuntos, sao relativamente
assimetricos em relacao a origem. Neste trabalho concentra-se a atencao nos pontos LGR.
Forma de Lagrange do Polinomio de Interpolacao
Interpolar uma funcao f consiste em aproximar essa funcao por uma outra funcao
g, escolhida entre uma classe de funcoes definida a priori e que satisfaca algumas propri-
edades. A necessidade de se fazer esta aproximacao surge em varias situacoes, como por
exemplo, quando sao conhecidos somente os valores numericos da funcao f para um con-
junto de pontos e e necessario calcular o valor da funcao em um ponto que nao se saiba
seu valor ou ainda quando a funcao f tem uma expressao tal que operacoes como a di-
ferenciacao e a integracao sao difıceis, ou impossıveis, de serem realizadas. Considera-se
9
aqui o caso de interpolacao polinomial [27] ou seja em que a funcao e aproximada por um
polinomio.
Dados n pontos (x1, f(x1)), ..., (xn, f(xn)), o objetivo e aproximar f por um po-
linomio pn−1, de grau menor ou igual a n− 1, tal que, para k = 1, 2, ..., n:
f(xk) = pn−1(xk).
Teorema 1.2 Existe um unico polinomio pn−1, de grau menor ou igual a n− 1, tal que,
para k = 1, 2, ..., n
pn−1(xk) = f(xk),
desde que xk 6= xj, j 6= k.
Demonstracao. [27, Teo.8.1].
Considere que os n pontos x1, . . . , xn pertencem a um intervalo (a, b). Considere f
contınua e sendo n vezes diferenciavel nos n pontos no intervalo (a, b). Denote yi = f(xi)
para i = 1, . . . , n. Seja pn−1 o polinomio de grau menor ou igual a n − 1 que interpola
f em x1, . . . , xn. Pode-se representar pn−1 na forma pn−1(x) = y1`1(x) + ... + yn`n(x),
onde os polinomios `k sao de grau n − 1. Para cada i, a condicao pn−1(xi) = yi deve ser
satisfeita, ou seja
pn−1(xi) = y1`1(xi) + ...+ yn`n(xi) = yi.
A forma mais simples de se satisfazer esta condicao e impor
`k(xi) =
0 se k 6= i
1 se k = i,
e assim define-se `k(x) por
`k(x) =n∏
j=1j 6=k
x− xjxk − xj
que satisfaz a condicao acima. Como o numerador de `k(x) e um produto de n−1 fatores,
entao `k e um polinomio de grau n− 1 e alem disso, para x = xi com i = 1, . . . , n tem-se
pn−1(xi) =n∑k=1
yk`k(xi) = yi`i(xi) = yi.
Entao, a forma de Lagrange para o polinomio interpolador e
pn−1(x) =n∑k=1
yk `k(x). (1.7)
10
onde yk = f(xk) e `k(x) =n∏
j=1j 6=k
x− xjxk − xj
.
1.2.2 Discretizacao do problema
O objetivo aqui e discretizar o problema de controle otimo (1.2), e para isto
utilizam-se os pontos de colocacao LGR, que sao pontos definidos no domınio [−1, 1],
como visto na Secao 1.2.1. Note que o problema (1.2) esta definido em um intervalo [t0, tf ],
o qual pode ser facilmente transformado no intervalo [−1, 1], atraves da transformacao
afim
t =tf − t0
2τ +
t0 + tf2
, (1.8)
para t ∈ [t0, tf ] e τ ∈ [−1, 1]. Derivando a igualdade (1.8) tem-se que
dt
dτ=tf − t0
2,
e pela regra da cadeia,dx
dt=dx
dτ
dτ
dt=
2
tf − t0dx
dτ. (1.9)
De (1.8) e (1.9), o problema (1.2) pode ser reescrito em funcao da nova variavel da seguinte
formaminimizar
∫ 1
−1 r(x(τ), u(τ))dτ
sujeito adx
dτ=tf − t0
2f(x(τ), u(τ))
x(−1) = x0
x(1) = xf .
(1.10)
Discretizacao da condicao diferencial
O intervalo [−1, 1] e discretizado considerando N pontos de colocacao LGR,
τ1, τ2, · · · , τN ∈ [−1, 1], com τ1 = −1 e τN < 1, vistos na Secao 1.2.1.
Utilizando (1.7), a j-esima componente de estado e aproximada por uma soma
da forma:
xj(τ) ≈N∑i=1
xij `i(τ) (1.11)
onde `i(τ) =N∏
m=1j 6=i
τ − τmτi − τm
, para i = 1, · · · , N, sao polinomios de Lagrange. Derivando a
aproximacao em relacao a τ tem-se que
dxjdτ
(τ) ≈N∑i=1
xijd`idτ
(τ). (1.12)
11
Avaliando (1.12) em τk e denotando ˙i(τk) = Dki, tem-se
xj(τk) ≈N∑i=1
Dki xij. (1.13)
A matriz DN×N = [Dki] cujo elemento (k, i) e dado por
Dki = ˙i(τk) (1.14)
e chamada de Matriz de Diferenciacao.
Agora seja XN×n a matriz formada pelos coeficientes xij de (1.11). Multiplicando
a matriz D por X, tem-se uma matriz [DX]N×n. Por (1.13) o elemento (i, j) dessa matriz
pode ser visto como
[DX]ij ≈ xj(τi). (1.15)
Seja UN×m = [uij] uma matriz cujo elemento (i, j) denota a aproximacao discreta
para a j-esima componente de controle avaliada no i-esimo ponto de colocacao, ou seja,
uij ≈ uj(τi).
A seguir um lema que permitira relacionar a linha da matriz X com x(τ).
Lema 1.3 Para todo i = 1, . . . , N , a i-esima linha Xi da matriz X contem as componen-
tes da aproximacao discreta x(τi)T .
Demonstracao. Lembrando que `i(τ) =N∏
m=1j 6=i
τ − τmτi − τm
tem-se que,
`i(τk) =
0 se i 6= k
1 se i = k
e por (1.13) conclui-se que a j-esima componente de x(τi)T e xij, que e exatamente a
j-esima componente de Xi.
Considere agora a matriz
F (X,U)N×n = [Fij(X,U)] = [fj(Xi, Ui)], (1.16)
onde fj e a j-esima componente da restricao do problema (1.10). Pelo Lema 1.3, (1.10) e
1.15 conclui-se que
DX =tf − t0
2F (X,U). (1.17)
Discretizacao da funcao custo
A funcao custo do problema (1.2) envolve uma integral. Recorre-se entao as
tecnicas de integracao numerica. Nesta secao apresenta-se a regra da quadratura, que
pode ser vista com detalhes em [9, Cap.4].
12
A ideia basica da quadratura consiste em escrever uma aproximacao da integral de
uma funcao, geralmente estabelecida como um somatorio com pesos dos valores assumidos
pela funcao em pontos especıficos dentro do domınio de integracao.
Uma regra de quadratura de N pontos e construıda para produzir um resultado
exato para polinomios de grau 2N − 1 ou menor para uma escolha adequada dos pesos
wi para i = 1, ..., N associado a cada ponto. O domınio de integracao de tal regra e
por convencao tomado como [−1, 1]. Sejam wi, com 1 ≤ i ≤ N , os pesos da quadratura
associados com os pontos LGR. Esses pesos tem a seguinte propriedade
∫ 1
−1P(τ)dτ =
N∑i=1
wiP(τi) (1.18)
para qualquer polinomio P de grau no maximo 2N − 2. Por [1], os pesos associados a
cada valor de P (τi) sao calculados por
wi =
∫ 1
−1
N∏k=1k 6=i
(τ − τkτi − τk
)dτ. (1.19)
Assim, usando (1.17), o problema (1.10) e escrito na forma discretizada como
minimizarN∑i=1
wi r(Xi, Ui))
sujeito a DX =tf − t0
2F (X,U)
X1 = x0
XN+1 = xf
(1.20)
onde XN+1 refere-se a variavel de estado no ponto τN+1 = 1.
1.2.3 Condicoes de Karush-Kuhn-Tucker
Como em [22], considere o problema na forma
minimizar Φ(x(1))
sujeito a dxdτ
= f(x(τ), u(τ)
x(−1) = x0,
(1.21)
e em sua forma discretizada
minimizar Φ(XN+1)
sujeito a DX = F (X,U)
X1 = x0.
(1.22)
13
Nesta secao sao desenvolvidas as condicoes de KKT (ver Apendice A.3.2) para o problema
(1.22) de uma maneira distinta de [22]. Considere Λ ∈ IRN×n a matriz composta pelos
multiplicadores de Lagrange associados as restricoes do problema (1.22) e seja o vetor
linha µ ∈ IRn composto pelos multiplicadores de Lagrange associados a condicao inicial
X1 = x0. Assim, tem-se
• Gradiente da funcao objetivo com relacao a X.
∇XΦ(XN+1) =(eN+1 ⊗∇Φ(XN+1)
T)(N+1).n×1 .
• Gradiente da funcao objetivo com relacao a U
∇UΦ(XN+1) = 0.
• Produto do multiplicador de Lagrange Λ ∈ IRN×n pelo gradiente da restricao
F (X,U)−DX = 0.
– Gradiente de DX em relacao a X.
(∇D(X))TΛ =N∑j=1
N∑k=1
DjkEkΛTj ,
onde Ek = ek ⊗ In.
– Gradiente de DX em relacao a U : nulo.
– Gradiente de F (X,U) em relacao a X
(∇XF (X,U))TΛ =N∑i=1
ei ⊗(∇Xf(Xi, Ui)Λ
Ti
).
– Gradiente de F (X,U) em relacao a U .
(∇UF (X,U))TΛ =N∑i=1
ei ⊗(∇Uf(Xi, Ui)Λ
Ti
).
• Produto do multiplicador µ ∈ IRn pelo gradiente de x0 −X1 = 0.
(∇XX1)Tµ = e1µ.
14
Com todas as parcelas apresentadas, pode-se expressar as condicoes de KKT do problema
(1.22):
N∑j=1
Dj1Λj = Λ1∇Xf(X1, U1)T − µ; (1.23)
N∑j=1
DjkΛj = Λk∇Xf(Xk, Uk)T , 2 ≤ k ≤ N ; (1.24)
∇φ(XN+1) = DTNΛ; (1.25)
Λk∇Uf(Xk, Uk)T = 0, 1 ≤ k ≤ N. (1.26)
Mas as igualdades (1.23) e (1.24) podem ser reescritas de forma conjunta como
DT1:NΛ = ∇X 〈Λ, F (X,U)〉 − e1µ. (1.27)
Seja Dj:k a submatriz de D formada pelas colunas j ate k e seja Xj:k a submatriz
de X formada pelas linhas de j ate k. Utilizando essa forma, as restricoes do problema
(1.22) podem ser escritas como
D2:NX2:N = F (X,U)−D1x0. (1.28)
O proximo resultado garante que a matriz D2:N , obtida pela omissao da primeira
coluna matriz D tem inversa.
Proposicao 1.4 A matriz D2:N obtida deletando a primeira coluna de D e invertıvel.
Demonstracao. Tome p′ ∈ IRN−1 nao nulo e considere p = (0, p′) ∈ IRN . Suponha agora
que Dp = 0. Mostra-se que D2:Np′ = 0 admite apenas solucao trivial e assim que D2:N e
nao-singular.
Seja P um polinomio de grau N tal que P(τk) = pk para 1 ≤ k ≤ N , onde pk e o
k-esimo elemento de p, ou seja, o polinomio escolhido deve ter a forma P(τ) =N∑k=1
pk`(τ).
Note que para i = 1, . . . , N
(Dp)i =N∑k=1
Dikpk =N∑k=1
pk ˙k(τi) = P(τi),
ou seja, cada componente de Dp e a derivada de P avaliada nos pontos de colocacao. Por
hipotese Dp = 0, e assim (Dp)i = 0 para todo i = 1, . . . , N , logo
P(τi) = 0.
Por P ser um polinomio de grauN−1 e se anular emN pontos segue que P e identicamente
nulo e portanto P e um polinomio constante. Como P(τN) = pN = 0 segue que P e um
15
polinomio identicamente nulo. Ou seja, se Dp = 0, com a primeira coordenada de p
nula, que e equivalente a ter D2:N p′ = 0 isto implica que p = 0, ou seja, p′ = 0 e assim
D2:N p′ = 0 admite apenas a solucao trivial e desta forma D2:N e nao-singular.
1.2.4 Sistema Adjunto Transformado
Agora reformula-se as condicoes de KKT do problema (1.22), utilizando os pesos
da discretizacao, de modo que se tornem condicoes transformadas para o problema (1.21).
Desta forma, tem-se as condicoes de KKT para o problema contınuo e para o problema
discreto.
Seja W ∈ RN×N matriz diagonal cujos elementos sao wi. Seja λ ∈ IRN×n definido
por
λ = W−1Λ. (1.29)
A fim de relacionar as equacoes discretas as equacoes contınuas, utiliza-se uma matriz
D+ ∈ IRN×N , uma versao modificada de D definida como:
D+11 = −D11 −
1
w1
e D+ij = −wj
wiDji. (1.30)
Usando (1.24), (1.29), (1.30), tem-se para i = 2, . . . , N ,
Λi∇Xf(Xi, Ui) =N∑j=1
DjiΛj =N∑j=1
−wiwjD+ijΛj = −
N∑j=1
wiD+ijλj.
Pela igualdade acima e (1.29), tem-se
N∑j=1
D+ijλj = −Λi
wi∇Xf(Xi, Ui) = −λi∇Xf(Xi, Ui). (1.31)
Por outro lado, usando (1.23), (1.29) e (1.30), obtem-se:
Λ1∇Xf(X1, U1)− µ =N∑j=1
Dj1Λj.
= D11Λ1 +N∑j=2
Dj1Λj
=
(−D+
11 −1
w1
)Λ1 +
N∑j=2
−w1
wjD+
1jΛj
= −D+11Λ1 − λ1 −
N∑j=2
D+1jw1λj.
16
Dividindo ambos os membros da igualdade por w1 e usando (1.29) tem-se
−D+11λ1 −
λ1w1
−N∑j=2
D+1jλj = λ1∇Xf(X1, U1)−
µ
w1
,
donde segue queN∑j=1
D+1jλj = −λ1∇Xf(X1, U1) +
1
w1
(µ− λ1), (1.32)
que e analoga a (1.23).
Para cada i = 1, . . . , N , dividindo a igualdade (1.26) por wi e usando (1.29),
resulta que
λi∇Uf(Xi, Ui) = 0. (1.33)
Agora discute-se a igualdade (1.25). Tome um polinomio P tal que P(τi) = 1,
para 1 ≤ i ≤ N . Pela definicao dos polinomios de Lagrange discutidos na Secao 1.2.1,
tem-se que P(τ) = 1 para todo τ , e assim P(τ) = 0. Considere v ∈ IRN um vetor de
componentes unitarias e (Dv)k a k-esima componente do produto Dv. Assim
(Dv)k =N∑i=1
Dkivi =N∑i=1
Dki =N∑i=1
˙i(τk) = P(τk) = 0.
Consequentemente 0 = Dv =N∑j=1
Dj =N−1∑j=1
Dj +DN . Logo
DN = −N−1∑j=1
Dj. (1.34)
Utilizando (1.34) tem-se
DTNΛ =
N∑i=1
ΛiDi,N = −N∑i=1
N∑j=1
ΛiDij.
Com (1.30) e fazendo mudanca de ındices vem que
−N∑i=1
N∑j=1
ΛiDij =Λ1
w1
+N∑i=1
N∑j=1
ΛiD+ji
wjwi
=Λ1
w1
+N∑i=1
N∑j=1
ΛjD+ij
wiwj.
17
Utiliza-se (1.29) e se isola λ1 de (1.32) daı
Λ1
w1
+N∑i=1
N∑j=1
ΛjD+ij
wiwj
= λ1 +N∑i=1
N∑j=1
wiλjD+ij
=
(−w1
N∑j=1
D+1jλj − λ1w1∇Xf(X1, U1) + µ
)+
N∑i=1
N∑j=1
wiλjD+ij .
Finalmente, com (1.31) e reorganizando as contas segue que(−w1
N∑j=1
D+1jλj − λ1w1∇Xf(X1, U1) + µ
)+
N∑i=1
N∑j=1
wiλjD+ij = −λ1w1∇Xf(X1, U1) + µ
+N∑i=2
(−λiwi∇XF (Xi, Ui))
= µ−N∑i=1
λiwi∇Xf(X1, Ui).
Tem-se, finalmente, as condicoes transformadas de KKT:
µ = ∇Φ(XN+1) +N∑i=1
wiλi∇Xf(Xi, Ui)
D+λ = −∇X 〈λ, F (X,U)〉+1
w1
e1(µ− λ1)
∇U 〈λ, F (X,U)〉 = 0.
Considere agora o seguinte resultado.
Teorema 1.5 [22, Teo.1] Considere P um polinomio de grau no maximo N − 1, com
N ≥ 1, e p ∈ IRN um vetor cuja i-esima componente e dada por pi = P(τi). Se D+
satisfaz, para todo i = 1, . . . , N ,
(D+p)i = P(τi),
entao D+ e a matriz de diferenciacao para o espaco de polinomios de grau N − 1 definida
em (1.30).
Demonstracao. Sejam P um polinomio de grau no maximo N com P(1) = 0 e Q um
polinomio de grau no maximo N − 1, com N ≥ 1. Usando integracao por partes, vale a
seguinte igualdade∫ 1
−1P(τ)Q(τ)dτ = −P(−1)Q(−1)−
∫ 1
−1P(τ)Q(τ)dτ. (1.35)
Note que, PQ e PQ sao polinomios de grau no maximo 2N − 2. Assim, pela Secao
18
1.2.2, a quadratura de Gauss e exata e consequentemente as integrais em (1.35) podem
ser substituıdas por suas quadraturas equivalentes, ou seja,
N∑j=1
wj pjqj = −p1q1 −N∑j=1
wjpj qj,
onde pj = P(τj) e pj = P(τj). De forma compacta, isto pode ser reescrito como
(Wp)T q = −p1q1 − (Wp)T q.
Substituindo p = D1:Np e q = D+q,
pTDT1:NWq + p1q1 + pTWD+q = 0
e assim,
pT(DT
1:NW +WD+ + e1eT1
)q = 0.
Como p e q foram tomados arbitrariamente, ou seja, essa igualdade deve ser satisfeita
para quaisquer p e q segue que
DT1:NW +WD+ + e1e
T1 = 0,
que implica (1.30).
1.2.5 Formulacao Integral
Nesta secao mostra-se que a discretizacao pseudoespectral da equacao de estado,
ou seja, da restricao do problema (1.22), tem uma formulacao equivalente atraves de
integrais.
Por (1.34) tem-se que D1 +D2 + · · ·+DN = 0 entao, se v e um vetor de compo-
nentes unitarias vale que
D1 = −N∑j=2
Dj = −D2:Nv. (1.36)
Pela Proposicao 1.4 a matriz D2:N e invertıvel entao
D−12:ND1 = −v. (1.37)
Seja P um polinomio qualquer de grau no maximo N . Por D ser a matriz de
diferenciacao dos polinomios de grau N tem-se que Dp = p onde pi = P(τi) e pi = P(τi)
19
para 1 ≤ i ≤ N . Entao,
p = Dp =N∑i=1
Dipi = D1p1 +D2:Np2:N ,
multiplicando por D−12:N tem-se
D−12:N p = D−12:ND1p1 + p2:N .
Utilizando (1.37), para i = 2, . . . , N vem que
pi = p1 + (D−12:N p)i. (1.38)
Agora obtem-se uma expressao diferente para p1 − pi baseado na integracao da
interpolacao da derivada. Seja `+j , para j = 1, . . . , N o polinomio de Lagrange interpolador
associado aos pontos de colocacao:
`+j =N∏
m=1m 6=j
τ − τmτj − τm
.
Dado um polinomio P de grau no maximo N , sua derivada P e um polinomio de
grau no maximo N − 1. Assim P pode ser interpolado exatamente pelos polinomios de
Lagrange `+j :
P(τ) =N∑j=1
pj`+j .
Integrando essa igualdade de −1 a τi tem-se
∫ τi
−1P(τ)dτ =
N∑j=1
pj
∫ τi
−1`+j (τ)dτ ,
pelo teorema fundamental do calculo e denotando∫ τi−1 `
+j (τ)dτ = Aij, para i = 2, . . . , N
vem que
P(τi)− P(−1) =N∑j=1
pjAij (1.39)
ou seja,
pi − p1 = (Ap)i. (1.40)
As relacoes (1.38) e (1.40) sao satisfeitas para qualquer polinomio de grau no
maximo N . Escolha um polinomio P tal que P(1) = 0, ou seja, p1 = 0. Daı por (1.38)
20
tem-se que pi = (D−12:N p)i e por (1.40) tem-se pi = (Ap)i e assim segue que
D−12:N p = Ap.
Tomando p como a i-esima coluna da matriz identidade na igualdade acima, conclui-se
que, para todo i = 1, . . . , N − 1, a i-esima coluna de D−12:N e igual a i-esima coluna da
matriz A, ou seja,
A = D−12:N .
Sabendo isso, multiplicando a equacao (1.28) por A, e usando (1.37), para i = 2, . . . , N ,
se obtem
Xi = x0 + AiF (X,U) (1.41)
onde Ai e a i-esima linha de A.
Assim, a forma diferencial da equacao de estado DX = F (X,U) e equivalente a
forma integrada (1.41), onde os elementos de A sao integrais do polinomio interpolador
de Lagrange `+j , enquanto os elementos de D sao as derivadas do polinomio de Lagrange
`i.
Para resumir, o que se tem em (1.41) e uma equacao que esta na forma de um
metodo de integracao global implıcito, enquanto a aproximacao DX = F (X,U) esta na
forma de um metodo pseudoespectral.
O fato de que a integral ou a forma diferencial poder ser usada, mostra que o
metodo de colocacao Radau pode ser visto como um metodo de integracao global implıcito
ou um metodo pseudoespectral. Em particular, usando a forma pseudoespectral dos
pontos de colocacao LGR, resulta num sistema de equacoes que nao tem qualquer perda
de informacao quando se passa para a forma integral.
1.3 Metodo de Euler
Esta secao, que tem [2, Sec. 3.1] como principal referencia, apresenta o metodo
de Euler para obtencao de uma solucao numerica do problema (1.2). O metodo de Euler
e baseado no Teorema de Taylor [26].
O intervalo de tempo [t0, tf ] e discretizado em (N − 1) subintervalos, de modo
que:
t0 ≡ t1 < t2 < · · · < tN ≡ tf .
Definindo as amplitudes de cada subintervalo i = 1, · · · , N−1 por hi = ti+1−ti, a variavel
de estado x(t) ∈ IRn e expandida, segundo Taylor, como
x(ti+1) = x(ti + hi) = x(ti) + hidx
dt(ti) +
h2i2
d2x
dt2(ςi),
para algum ςi ∈ [ti, ti+1]. Assim, para valores suficientemente pequenos de hi, tem-se a
21
seguinte aproximacaodx
dt(ti) ≈
x(ti+1)− x(ti)
hi.
Como o estado x satisfaz a equacao diferencial do problema (1.2) tem-se que
x(ti+1)− x(ti)
hi= f(x(ti), u(ti)),
para hi suficientemente pequeno.
Tomando-se amplitude h > 0 constante suficientemente pequena para todos os
subintervalos e denotando x(tk) por xk e u(tk) por uk, a condicao diferencial acima e
escrita da seguinte forma,
f(xk, uk) =xk+1 − xk
h, (1.42)
para todo k = 1, · · · , N − 1. Analogamente, segundo [2], a discretizacao da funcao custo
do problema (1.2) e dada porN−1∑k=0
r(xk, uk)h.
Assim o problema (1.2) e discretizado, segundo o metodo de Euler, como
minimizarN−1∑k=1
r(xk, uk)h
sujeito a f(xk, uk) =xk+1 − xk
hparak = 1, · · · , N − 1
x(t1) = x1
x(tN) = xN .
(1.43)
A resolucao deste problema de otimizacao leva a uma solucao numerica do problema
original. No proximo capıtulo, veremos um exemplo particular.
1.4 Outros metodos
Dois metodos de colocacao LGR sao apresentados em [14] e [15]. O metodo de
Kameswaran e Biegler em [15] concentra-se na colocacao local usando pontos LGR. O
metodo de Fahroo e Ross em [14] descreve um metodo global para resolver problemas
de horizonte infinito. Nesta secao, um breve comentario sobre a forma como o metodo
apresentado neste trabalho refere-se a estes trabalhos.
1.4.1 Comparacao com o metodo de Kameswaran e Biegler
O metodo pseudoespectral abordado neste trabalho tem semelhancas com o metodo
de Kameswaran e Biegler ([15]). A aproximacao para a variavel de estado usa tambem
os polinomios de Lagrange. E observado, no entanto, que o metodo de Kameswaran e
22
Biegler usa colocacao local, utilizando varios subintervalos. O grau dos polinomios em
cada subintervalo e fixo e a convergencia e conseguida atraves do aumento do numero de
subintervalos. O metodo aqui tratado e um metodo de colocacao global, em que ha um
unico intervalo e a convergencia e alcancada atraves do aumento do grau de polinomios,
ou seja, no aumento do N , dos pontos de colocacao, como visto pela Secao 1.2.1. O
metodo de Kameswaran e Biegler e implementado de uma forma semelhante ao metodo
de Euler tratado na Secao 1.3 (devido ao fato de que o intervalo de tempo e dividido em
subintervalos), enquanto que o metodo do presente trabalho e implementado na forma
de um metodo pseudoespectral. Nota-se que ambas as abordagens sao validas, mas a
abordagem atual, segundo [22], e usada com mais frequencia na literatura de controle.
1.4.2 Comparacao com o metodo de Fahroo e Ross
Na abordagem pseudoespectral Lobatto como descrito em [21], a variavel de es-
tado e aproximada por polinomios de grau N − 1 e a dinamica do sistema e discretizado
com N pontos de quadratura Lobatto. Para o problema de controle do horizonte infinito
estudado em [14], Fahroo e Ross propoem utilizar uma mudanca de variaveis para trans-
formar o intervalo de tempo infinito para o intervalo [−1,+1). Esta mudanca de variaveis
leva a uma singularidade na dinamica no ponto τ = +1. Desta forma, nao e possıvel
utilizar τ = +1 como um ponto de colocacao. Para lidar com essa singularidade, Fahroo
e Ross propoem discretizar nos pontos de quadratura Radau para τ = N < 1.
A diferenca fundamental entre o metodo pseudoespectral em [14] e o metodo pre-
visto neste trabalho e que, em [14], a variavel de estado e aproximada por polinomios
de grau N − 1, enquanto que neste trabalho, a variavel de estado e aproximada usando
polinomios de grau N . Esta mudanca no grau dos polinomios leva a diferencas fundamen-
tais entre os dois esquemas. Por exemplo, uma vez que os polinomios de Lagrange sao
de diferentes graus, as matrizes de diferenciacao sao completamente diferentes. A matriz
utilizada na diferenciacao em [14] e singular, enquanto que a matrize D2:N , segundo a Pro-
posicao 1.4, e invertıvel. Se o controle e o estado inicial x0 sao dados, entao a dinamica
discretizada do problema [14] e um sistema de N equacoes e de N − 1 incognitas, X2:N .
Em contraste, 1.28 e um sistema de N − 1 equacoes com N − 1 incognitas, X2:N , onde a
matriz de coeficientes D2:N e invertıvel.
Na abordagem de [14], XN+1, a estimativa da variavel de estado em τ = 1, e
removido do problema usando polinomios de grau N−1 em vez de polinomios de grau N .
Na abordagem aqui apresentada, a variavel de estado e aproximada em τi, 1 ≤ i ≤ N + 1.
Assim, XN+1, a estimativa do estado, e uma variavel incluıda no esquema pseudoespectral.
Avaliar o estado em τ = +1 e util quando a funcao objetivo depende do estado no
momento final, ou quando ha uma restricao de ponto final, como era o caso do problema
abordado neste trabalho.
23
1.5 Problema de Bolza
No proximo capıtulo resolve-se um problema que e chamado de problema de
Bolza. Para mais detalhes sobre esse tipo de problema consultar [5]. O problema de
Bolza pode ser formulado de diversas maneiras, cada uma das quais tem as suas vantagens
peculiares e desvantagens. Uma das formulacoes mais uteis pode ser descrita brevemente
como segue:
Determinar o estado x(τ) ∈ IRn, o controle u(τ) ∈ IRm, o tempo incial t0 e tempo
final tf que encontre um mınimo local do seguinte problema
minimizar J = Φ(x(−1), t0, x(1), tf ) +tf − t0
2
∫ 1
−1 g(x(τ), u(τ), τ)dτ
sujeito adx
dτ=
tf−t02f(x(τ), u(τ), τ) ∈ IRn
φ(x(−1), x(1)) = 0 ∈ IRq
C(x(τ), u(τ)) ≤ 0 ∈ IRc.
(1.44)
Quando Φ = 0 o problema e o chamado problema de Lagrange e quando g = 0
tem-se o problema geral de Mayer como formulado por Bliss em [5]. Estas tres versoes,
ou seja, o problema completo, o problema com Φ = 0 e o problema com g = 0, sao
equivalentes no sentido de que cada um pode ser transformado em qualquer um dos
outros dois tipos.
Destas tres versoes, o problema de Bolza formulado acima parece ser o mais
conveniente, pois possui todos os termos de forma explıcita. Para o presente trabalho
formula-se o problema utilizando uma aproximacao pseudoespectral LGR. E assim, o
problema de Bolza toma a seguinte forma:
minimizar Φ(X1, t0, XN+1, tf ) +tf − t0
2
N∑k=1
wkg(Xk, Uk, τk)
sujeito a DX =tf − t0
2F (X,U)
φ(X1, XN+1) = 0
C(Xk, Uk) ≤ 0 1 ≤ k ≤ N,
onde os w′ks sao dados por (1.19).
Capıtulo 2
Exemplos
Neste capıtulo consideram-se alguns problemas de controle otimo particulares que
sao resolvidos pelo metodo pseudoespectral utilizando pontos LGR discutido no capıtulo
anterior. Inicialmente, discute-se um problema de Bolza sem restricoes [2, 22] dito pro-
blema do regulador linear quadratico, conhecido pela abreviacao LQR, do ingles linear
quadratic regulator. Discute-se sua solucao exata e solucoes numericas utilizando o metodo
de Euler visto na Secao 1.3 e o metodo pseudoespectral visto na Secao 1.2. Em seguida,
sao discutidas as solucoes numericas de dois problemas de controle otimo com restricoes
nao lineares.
2.1 Exemplo 1
O primeiro exemplo e um problema do regulador linear quadratico, conhecido
pela abreviacao LQR, do ingles linear quadratic regulator que consiste em minimizar uma
funcao custo quadratica com tempos iniciais e finais fixos [2, pag. 126]. O funcional e
dado por
J =1
2x(tf )
TSx(tf ) +1
2
∫ tf
t0
(x(t)TQx(t) + u(t)TRu(t)
)dt,
sujeito a um sistema linear dinamico
dx
dt= Ax+Bu,
e as condicoes
x(t0) = x0, x(tf ) = xf .
O numero de estados e n, x(t) ∈ IRn, o numero de controles e m, u(t) ∈ IRm, tais que
S ∈ IRn×n, Q ∈ IRn×n, R ∈ IRm×m, A ∈ IRn×n e B ∈ IRn×m.
A solucao exata do problema [2] pode ser encontrada utilizando o Princıpio do
Maximo de Pontryagin, que faz recair no seguinte sistema
24
Exemplos 25
(x(t)
p(t)
)=
(A −BR−1BT
−Q −AT
)(x(t)
p(t)
),
onde p e a variavel de co-estado.
O controle otimo e definido por
u(t) = −R−1BTp(t).
O caso especıfico considerado aqui tem um unico estado e um unico controle de
modo que n = m = 1 e S,Q,R,A e B sao escalares. Em particular, tome S = 0 e
Q = R = A = B = 1. Alem disso, considere as condicoes de contorno: inicial x(t0) =√
2
e a condicao final x(tf ) = 1, com t0 = 0 e tf = 5. Assim, o problema considerado e dado
por
maximizar −∫ 5
0(x2(t) + u2(t)) dt
sujeito a x(t) = x(t) + u(t)
x(0) =√
2
x(5) = 1.
(2.1)
Aplicam-se as tecnicas discutidas no capıtulo anterior para obtencao da solucao
exata e de solucoes aproximadas para este problema.
2.1.1 Solucao Exata
A solucao exata do problema e obtida atraves da aplicacao do Princıpio do
Maximo de Pontryagin, como visto no Teorema 1.1.
Identificam-se os elementos do problema (2.1) com os elementos do problema de
controle otimo dado em (1.2) e com os elementos do Teorema 1.1. Para nao carregar a
notacao o argumento t e omitido, ou seja, escreve-se, por exemplo, x ao inves de x(t), e
assim para todas as variaveis que dependem de t. Assim, desta identificacao, tem-se
f(x, u) = x+ u,
r(x, u) = −x2 − u2,H(x, p, u) = f(x, u)p+ r(x, u).
Desta forma,
x = ∇pH(x, p, u) = f(x, u), (2.2)
p = −∇xH(x, p, u) = −∂f∂yp− ∂r
∂y= −p+ 2x, (2.3)
Hu(x, p, u) = 0 =∂f
∂up+
∂r
∂u= p− 2u. (2.4)
Exemplos 26
Das igualdades acima, tem-se o seguinte sistemax = x+ p
2
p = 2x− p
que pode ser escrito na forma matricial(x
p
)=
(1 1
2
2 −1
)︸ ︷︷ ︸
M
(x
p
).
A solucao deste sistema de equacoes diferenciais, como discutido na Secao A.2, e(x
p
)= eMt
(x0
p0
). (2.5)
Entao precisa-se calcular eMt e encontrar p0, e assim x e p estarao bem definidos. A
matriz M e uma matriz diagonalizavel, ou seja, como visto na Secao A.1.2, existe uma
matriz S tal que M = SDS−1 com D diagonal. Neste caso,
S =
(1 −1
2√
2− 2 2√
2 + 2
);D =
( √2 0
0 −√
2
);S−1 =
( √2+24
14√2√
2−24
14√2
).
Pela Proposicao A.6, eMt = SeDtS−1, onde
eDt =
(e√2t 0
0 e−√2t
).
Desta forma,
eMt =
(1 −1
2√
2− 2 2√
2 + 2
)(e√2t 0
0 e−√2t
)√
2 + 2
4
1
4√
2√2− 2
4
1
4√
2
.
Fazendo a multiplicacao das matrizes tem-se
eMt =
√
2 + 2
4e√2t +
2−√
2
4e−√2t e
√2t − e−
√2t
4√
2√2
2e√2t −√
2
2e−√2t 2−
√2
4e√2t +
2 +√
2
4e−√2t
.
Usando as condicoes iniciais x0 = x(0) =√
2 e x(5) = 1 em (2.5) tem-se que(1
p(5)
)= e5M
( √2
p0
).
Exemplos 27
Denotando
e5M =
(B11 B12
B21 B22
),
tem-se que
1 =√
2B11 + p0B12,
e assim, como B12 6= 0,
p0 =1−B11
√2
B12
.
Substituindo os valores de B11 e B12 e fazendo algumas manipulacoes algebricas conclui-se
que
p0 = p(0) =4√
2− (2√
2 + 4)e5√2 + (2
√2− 4)e−5
√2
e5√2 − e−5
√2
. (2.6)
Voltando em (2.5) e usando a condicao inicial x0 =√
2, tem-se que(x(t)
p(t)
)= eMt
( √2
p(0)
)
ou seja
x(t) =√
2
(√2 + 2
4e√2t +
2−√
2
4e−√2t
)+
(e√2t − e−
√2t
4√
2
)p(0) (2.7)
e
p(t) =√
2
(√2
2e√2t −√
2
2e−√2t
)+
(2−√
2
4e√2t +
2 +√
2
4e−√2t
)p(0),
com p(0) dado em (2.6). Da igualdade (2.4) segue que
u(t) =p(t)
2.
Assim a solucao do problema (2.1) ex(t) =
√2(√
2+24e√2t + 2−
√2
4e−√2t)
+(e√2t−e−
√2t
4√2
)p(0)
u(t) =
((2√
2− 4)e−5√2 + 2
√2− 2
2e5√2 − 2e−5
√2
)e√2t +
(−(2√
2 + 4)e5√2 + 2
√2 + 2
2e5√2 − 2e−5
√2
)e−√2t
(2.8)
A Figura 2.1 exibe graficamente o controle u e a variavel de estado y = x2.
Exemplos 28
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5−4
−3
−2
−1
0
1
2
t
Solu
ção
y(t)
u(t)
Figura 2.1: Solucao exata do problema.
2.1.2 Solucao aproximada pelo metodo pseudoespectral
Nesta secao discute-se a solucao numerica pelo metodo pseudoespectral discutido
na Secao 1.2 do problema (2.1) que escrito na forma de minimizacao e dado por,
minimizar 12
∫ 5
0(x2(t) + u(t)2) dt
sujeito a x(t) = x(t) + u(t)
x(0) =√
2
x(5) = 1.
(2.9)
Inicialmente, considera-se a mudanca de variavel (1.8) de t ∈ [−1, 1] para τ ∈ [0, 5],
t =5
2(τ + 1).
Assim,dτ
dt=
2
5e x =
dx
dt=dx
dτ
dτ
dt=
2
5
dx
dτ
e a restricao do problema (2.14) se escreve como
2
5
dx
dτ= x(τ) + u(τ).
Sejam N pontos de colocacao LGR, −1 = τ1 < τ2 · · · < τN < 1, obtidos como raızes
da soma PN−1(τ) + PN(τ) dos polinomios de Legendre, como discutido na Secao ??.
Considere D ∈ IRN×N a matriz de diferenciacao definida em (1.14). Por (1.15), tem-se
que DX ≈ dx
dτ. Assim a restricao se reescreve como
2
5DX = X + U,
Exemplos 29
com X,U ∈ IRN . Tendo calculado o vetor de pesos w ∈ IRN associados aos pontos de
colocacao por (1.19), o problema (2.14) discretizado e dado por
minimizar1
2
N∑k=1
(x2k + u2k)wk
sujeito a 25DX = X + U
x1 =√
2
xN+1 = 1.
(2.10)
Os pontos de colocacao e o vetor w ∈ IRn foram calculados utilizando a rotina
pontos lgr w em Matlab descrita no Apendice B. Este programa consiste numa adaptacao
para calculo dos pontos de colocacao LGR [23] do codigo lglnodes.m escrito por Winckel
[32] que tem como referencia [7]. O resultado obtido pelo programa para diferentes valores
de N foi comparado com os valores fornecidos em [13]. A matriz de diferenciacao D ∈IRN×N foi calculada pela rotina matriz dif.
Tendo o vetor w ∈ IRN de pesos da quadratura associados com os N pontos de
colocacao LGR e a matriz D de diferenciacao, o problema (2.10) pode ser visto como um
problema de minimizacao de uma funcao quadratica com restricoes de igualdade lineares
e condicao de caixa, na variavel (X,U) ∈ IR2N , com uma variavel adicional que deve
satisfazer a condicao xN+1 = 1. Os problemas, para diferentes valores de N e solucoes
iniciais randomicas, foram resolvidos utilizando-se a rotina quadprog do Matlab. A Figura
2.2 apresenta o grafico da solucao exata e uma solucao numerica obtida adotando-se
N = 60 pontos de colocacao. As linhas cheias representam a solucao exata no intervalo
[−1, 1]. Os sımbolos triangulares representam a solucao aproximada para a variavel de
estado y = x2 e o sımbolo asterisco e usado para representar a variavel de controle u.
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1−4
−3
−2
−1
0
1
2
τ
yexato
uexato
yaprox
uaprox
Figura 2.2: Solucao exata e solucao obtida pelo metodo pseudoespectral.
Exemplos 30
Analise de erros.
Foram consideradas 30 instancias do problema (2.10) para cada valor de N entre
10 e 100. Cada instancia difere na solucao inicial que e tomada arbitraria. Para cada
instancia foram calculadas a norma da diferenca entre a solucao exata e a solucao numerica
tanto para a variavel de estado y = x2 e a variavel de controle u, ou seja, foram calculados
os erros:
erro y = ‖yexata − yaprox‖∞erro u = ‖uexata − uaprox‖∞ ,
onde (yexata, uexata) representa a solucao exata e (yaprox, uaprox) representa a solucao apro-
ximada. Para cada valor de N foi calculada a mediana dos erros calculados para as 30
instancias. A Figura 2.3 mostra os graficos da variacao deste erro quando varia-se o valor
de N . O grafico da esquerda mostra o erro na variavel de estado, enquanto o grafico da
direita mostra o grafico do erro para a variavel de controle. O eixo vertical de ambos
os graficos esta na escala logarıtmica. As curvas em linha contınua referem-se as curvas
de regressao linear e quadratica na escala logarıtmica, cujos parametros foram obtidos
atraves do comando polyfit do Matlab. Note que os erros descrescem mais lentamente na
variavel de controle u.
10 20 30 40 50 60 70 80 90 100
10−2
10−1
N
log(e
rro
y)
erroy
4.25N−1.44
0.0979N0.691 + (−0.292 ln(N))
10 20 30 40 50 60 70 80 90 100
10−1.4
10−1.3
10−1.2
N
log(e
rro
u)
errou
0.0832N−0.165
0.465N−1.14 + (0.133 ln(N))
Figura 2.3: Variacao em funcao de N do erro na variavel de estado (esquerda) e na variavelde controle (direita).
2.1.3 Solucao aproximada pelo metodo de Euler
Nesta secao sera aplicado o metodo de Euler estudado na Secao 1.3 para resolver
o problema (2.1). O intervalo de tempo [0, 5] e inicialmente transformado no intervalo
[−1, 1], utilizando a mudanca dada em (1.8) e este e discretizado em N subintervalos de
amplitude2
N. Adotando-se t0 = −1 e tN = 1, tem-se que o (k + 1)-esimo subintervalo
e [tk, tk+1] para k = 0, . . . , N − 1. Denota-se x(tk) por xk e u(tk) por uk e estima-se a
Exemplos 31
derivada x porxk+1 − xk
hem cada intervalo. Assim, a restricao principal do problema
(2.1) se reescreve como k restricoes da forma:
xk+1 − xk =5
N(xk + uk), (2.11)
com k = 0, . . . , N − 1.
A funcao objetivo dada por uma integral e discretizada atraves do somatorio
1
2
N−1∑k=0
(x2k + u2k
)h,
e o problema (2.1) se reescreve na forma discretizada como
minimizar 12
N−1∑k=0
(x2k + u2k)h
sujeito a xk+1 − xk =5
N(xk + uk), k = 0, . . . , N − 1
x0 =√
2
xN = 1.
A condicao de contorno no tempo final foi incorporada ao problema como uma
penalizacao a funcao objetivo, ou seja, fixando-se um valor M suficientemente grande, o
problema anterior e transformado em
minimizar 12
N−1∑k=0
(x2k + u2k)h+M(xN − 1)2
sujeito a xk+1 − xk =5
N(xk + uk), k = 0, . . . , N − 1
x0 =√
2.
(2.12)
Para as implementacoes, foi tomado M = 105.
O lema a seguir, relaciona cada variavel de estado xk com a posicao inicial x0 e
com as variaveis de controle uj com j = 0, · · · , k − 1.
Lema 2.1 Para todo k = 1, . . . , N , vale a relacao
xk = (1 + h)kx0 + hk−1∑j=0
(1 + h)k−j−1uj.
Demonstracao. A prova e feita por inducao. Primeiramente verifica-se que vale para k = 1.
Note que por (2.11) tem-se
x1 − x0 = h(x0 + u0),
Exemplos 32
ou seja,
x1 = (1 + h)x0 + hu0,
mostrando assim que vale para k = 1. Assume-se que a proposicao vale para algum k,
com 1 ≤ k ≤ N − 1, ou seja,
xk = (1 + h)kx0 + h
k−1∑j=0
(1 + h)k−j−1uj.
Prova-se entao que vale para k + 1. Utilizando novamente (2.11) tem-se que
xk+1 = xk + h(xk + uk) = (1 + h)xk + huk.
Usando a hipotese de inducao, a igualdade acima se torna
xk+1 = (1 + h)
[(1 + h)kx0 + h
k−1∑j=0
(1 + h)k−j−1uj
]+ huk.
Fazendo as multiplicacoes e agrupando os termos tem-se
xk+1 = (1 + h)k+1x0 + hk−1∑j=0
(1 + h)k−juj + h(1 + h)0uk
= (1 + h)k+1x0 + hk−1∑j=0
(1 + h)(k+1)−j−1uj + h(1 + h)(k+1)−k−1uk
= (1 + h)k+1x0 + hk∑j=0
(1 + h)(k+1)−j−1uj.
Provando assim que, para qualquer k = 1, . . . , N , vale a relacao
xk = (1 + h)kx0 + h
k−1∑j=0
(1 + h)k−j−1uj. (2.13)
Cabe ressaltar que as restricoes do problema (2.12) valida para k = 0, . . . , N − 1
foram usadas fortemente na prova do Lema 2.1. O resultado provado sera usado para
escrever a funcao objetivo do problema (2.12) apenas em funcao das variaveis de controle
uk e de x0 que pela condicao inicial e fixado em√
2.
Usando (2.13) a funcao objetivo do problema (2.12) se escreve como
Exemplos 33
J =1
2
N−1∑k=0
(x2k + u2k)h+M(xN − 1)2
=1
2(x20 + u20)h+
1
2
N−1∑k=1
(x2k + u2k)h+M(xN − 1)2
=1
2(x20 + u20)h
+1
2
N−1∑k=1
((1 + h)2kx20 + 2hx0(1 + h)k
k−1∑j=0
(1 + h)k−j−1uj+
+h2k−1∑j=0
k−1∑m=0
(1 + h)2k−j−m−2ujum + u2k
)h
+M
(((1 + h)Nx0 − 1)2 + 2h((1 + h)Nx0 − 1)
N−1∑j=0
(1 + h)N−j−1uj+
+h2N−1∑j=0
N−1∑m=0
(1 + h)2N−j−m−2ujum
),
que pode ser reescrita como
J =1
2hx20
N−1∑k=0
(1 + h)2k +M(x0(1 + h)N−1 − 1)2 + h2x0
N−1∑k=1
k∑j=0
(1 + h)2k−j−2uj
+ 2Mh(x0(1 + h)N−1 − 1)N−2∑j=0
(1 + h)N−j−2uj +h3
2
N−1∑k=1
k∑j=0
k∑m=0
(1 + h)2k−j−m−2ujum
+ h2M
N−2∑j=0
N−2∑m=0
(1 + h)2N−j−m−2ujum +h
2
N−1∑k=0
u2k,
que e uma funcao quadratica na variavel de controle u, visto que x0 =√
2.
Assim (2.12) recai no problema de minimizacao da funcao quadratica J , sem
restricoes, que pode ser resolvido pelo Algoritmo A.1 do metodo de gradientes conjugados.
Tendo obtido o vetor u das variaveis de controle, obtem-se o vetor x atraves da expressao
(2.13), o que permite determinar o vetor y das variaveis de estado atraves da igualdade
y = x2.
Exemplos 34
A listagem referente a implementacao necessaria para resolucao do problema (2.1)
pelo metodo de Euler descrito nesta secao e apresentada no Apendice B.
A Figura 2.4 apresenta a solucao obtida adotando-se N = 60
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1−4
−3
−2
−1
0
1
2
τ
yexato
uexato
yaprox
uaprox
yEuler
uEuler
Figura 2.4: Solucao exata e solucoes numericas.
Analise de erros. Como para o caso da solucao numerica obtida pelo metodo pseu-
doespectral, foram consideradas 30 instancias do problema (2.10) para cada valor de N
entre 10 e 100. Cada instancia difere na solucao inicial que e tomada arbitraria. Para
cada instancia foram calculadas a norma da diferenca entre a solucao exata e a solucao
numerica tanto para a variavel de estado y = x2 e a variavel de controle u, ou seja, foram
calculados os erros:
erro y = ‖yexata − yEuler‖∞erro u = ‖uexata − uEuler‖∞ ,
onde (yexata, uexata) representa a solucao exata e (yEuler, uEuler) representa a solucao apro-
ximada obtida pelo metodo de Euler. Para cada valor de N foi calculada a mediana dos
erros obtidos para as 30 instancias. A Figura 2.5 mostra os graficos da variacao deste erro
quando varia-se o valor de N . O grafico da esquerda mostra o erro na variavel de estado,
enquanto o grafico da direita mostra o grafico do erro para a variavel de controle. O eixo
vertical de ambos os graficos esta na escala logarıtmica. A linha cheia representa a curva
da regressao linear na escala logarıtmica. Os coeficientes da regressao foram fornecidos
pelo comando polyfit do Matlab.
O erro na variavel de estado e da ordem de 10−1, enquanto para o metodo pseu-
doespectral o erro foi da ordem de 10−2. Os erros para o metodo de Euler foram maiores
que para o metodo pseudoespectral.
Exemplos 35
10 20 30 40 50 60 70 80 90 100
10−1
N
log(e
rro
y)
erroy
0.713N−0.485
10 20 30 40 50 60 70 80 90 100
10−0.6
10−0.5
10−0.4
10−0.3
10−0.2
10−0.1
N
log(e
rro
u)
errou
2.91N−0.539
Figura 2.5: Variacao em funcao de N do erro na variavel de estado (esquerda) e na variavelde controle (direita).
2.2 Exemplo 2
Considere o problema de controle otimo, discutido em [31, pag. 347], dado por:
minimizar∫ 1
0(u2(t)− y(t)3) dt
sujeito a y(t) = y(t)u(t)
y(0) = 1
y(1) = 1
−1 ≤ u(t) ≤ 0,
(2.14)
onde u ∈ C[0, 1] e y ∈ C1[0, 1], ou seja, u : [0, 1]→ IR e uma funcao contınua por partes e
y : [0, 1]→ IR e uma funcao continuamente diferenciavel. A solucao deste problema ey(t) ≡ 1,
u(t) ≡ 0.(2.15)
A seguir discute-se a resolucao numerica pelo metodo estudado na Secao 1.2.
2.2.1 Solucao aproximada pelo metodo pseudoespectral
De maneira analoga ao exemplo anterior, se faz a mudanca de variavel de t ∈ [0, 1]
para τ ∈ [−1, 1]. Por (1.8) tem-se que:
t =1
2τ +
1
2.
Assim,dτ
dt= 2 e y =
dy
dt=dy
dτ
dτ
dt= 2
dy
dτ
Exemplos 36
e a restricao do problema (2.14) se escreve como
dy
dτ=
1
2y(τ)u(τ).
Sejam N pontos de colocacao LGR, −1 = τ1 < τ2 < ... < τN < 1, obtidos como raızes
da soma dos polinomios de Legendre PN−1(τ) + PN(τ). Considere D ∈ IRN×N a matriz
de diferenciacao definida em (1.14). Por (1.15) tem-se que DY ≈ dy
dτe a restricao se
reescreve como
DY =1
2Y U,
onde Y, U ∈ IRN .
Tendo os pontos de colocacao, e sabendo que os pesos w′is, com i = 1, . . . , N ,
para a discretizacao do problema, e dado, segundo [1], por
w(i) =
∫ 1
−1
N∏j=1j 6=i
τ − τjτi − τj
dτ,
pode-se discretizar o problema (2.14), que toma a forma
minimizarN∑k=1
(u2k − y3k)wk
sujeito a 2DY = Y U
y1 = 1
yN+1 = 1.
(2.16)
Os procedimentos utilizados para resolver este problema sao analogos aos que
foram utilizados para resolver o exemplo anterior. O problema foi resolvido numerica-
mente em Matlab. O vetor w ∈ IRN de pesos da quadratura associados com os N pontos
de colocacao LGR foram obtidos pela rotina pontos lgr w. A matriz de diferenciacao
D ∈ IRN×N foi calculada pela rotina matriz dif. O problema 2.16 que consiste na mi-
nimizacao de uma funcao nao linear, nao quadratica, sujeita a restricoes de igualdade
nao lineares e condicoes de contorno na variavel (y, u) ∈ IR2N , foi resolvido pela funcao
fmincon disponıvel em Matlab. As rotinas utilizadas estao disponıveis no Apendice B.
A Figura 2.6 apresenta uma solucao do problema tomando-se N = 60. O sımbolo
triangular representa a solucao aproximada para a variavel de estado y e o sımbolo as-
terisco e usado para representar a variavel de controle u associadas com cada um dos 60
pontos de colocacao.
Exemplos 37
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1−0.2
0
0.2
0.4
0.6
0.8
1
1.2
τ
y
u
Figura 2.6: Solucao obtida pelo metodo pseudoespectral.
Analise de erros. Adotando-se solucoes iniciais arbitrarias, foram resolvidas 30 instancias
do problema (2.16) para cada valor de N entre 2 e 100. Para cada uma destas instancias,
foi calculada a norma infinito do vetor diferenca entre a solucao exata e a solucao numerica
para a variavel de estado y e para a variavel de controle u. Ou seja, foram calculados os
erros dados por
erroy = ‖yexata − yaprox‖∞ e errou = ‖uexata − uaprox‖∞.
onde (yexata, uexata) e dado em (2.15) e (yaprox, uaprox) e a solucao numerica fornecida
pela rotina fmincon. Para cada um dos valores de N considerados, foi tomada a mediana
dos erros obtidos nas 30 instancias. Os graficos da Figura 2.7 mostram a variacao do erro
na escala logarıtmica para a solucao aproximada do estado e de controle, respectivamente,
em relacao ao numero de pontos de colocacao N entre 2 e 100. Note que os erros sao
da ordem de 10−8 e que neste caso, diferentemente do caso anterior, os erros aumentam
com o valor de N . As linhas cheias correspondem a curva de regressao linear na escala
logarıtmica cujos coeficientes foram fornecidos pelo comando polyfit do Matlab.
2.3 Exemplo 3
O segundo problema consiste em determinar o estado y ∈ C1 e o controle u ∈ Ccomo solucao do seguinte problema de Bolza:
minimizar y2(2) +∫ 2
0y(t)2dt
sujeito a y(t) = u(t)
y(0) = 1
−1 ≤ u(t) ≤ 1,
(2.17)
Exemplos 38
0 20 40 60 80 10010
−13
10−12
10−11
10−10
10−9
10−8
N
log(e
rro
y)
erroy
3.27e−14N2.24
0 20 40 60 80 10010
−13
10−12
10−11
10−10
10−9
10−8
10−7
10−6
10−5
N
log(e
rro
u)
errou
1.41e−14N2.99
Figura 2.7: Variacao do erro das variaveis y e u na escala logarıtmica em relacao a N .
cuja solucao exata apresentada em [31, pag 348] e
y(t) =
1− t se t ≤ 1,
0 se t > 1,e u(t) =
−1 se t ≤ 1,
0 se t > 1.(2.18)
2.3.1 Solucao aproximada pelo metodo pseudoespectral
Inicialmente considera-se a mudanca de variavel de t ∈ [0, 2] para τ ∈ [−1, 1], ou
seja,
t = τ + 1.
Assim,dτ
dt= 1 e y =
dy
dt=dy
dτ
dτ
dt=dy
dτ
e a restricao do problema (2.17) se escreve como
dy
dτ= u(τ).
Calculam-se os N pontos de colocacao LGR, −1 = τ1 < τ2 < ... < τN < 1, o vetor de
pesos w ∈ IRN e a matriz de diferenciacao D ∈ IRN×N como discutidos nos exemplos
anteriores. Como a funcao objetivo do problema (2.17) envolve o valor y(2), considera-se,
alem das variaveis Y ∈ IRN e U ∈ IRN , uma variavel adicional denotada por YN+1. Assim
o problema (2.17) se escreve na forma discretizada como
minimizar Y 2N+1 +
N∑k=1
y2kwk
sujeito a DY = U
y1 = 1
−1 ≤ U ≤ 1.
(2.19)
Exemplos 39
Note que o problema acima consiste na minimizacao de uma funcao quadratica com
restricoes de igualdade lineares e condicoes de contorno, na variavel (Y , U) ∈ IR2N+1,
onde Y = (Y, YN+1) ∈ IRN+1 e U ∈ IRN . Deste modo este problema pode ser resolvido
pela rotina quadprog disponıvel no Matlab.
A Figura 2.8 apresenta a solucao exata e a solucao numerica, considerando-se 60
pontos de colocacao e solucao inicial arbitraria. A solucao exata (2.18) que se escreve na
variavel τ como
y(τ) =
−τ se τ ≤ 0,
0 se τ > 0,e u(τ) =
−1 se τ ≤ 0,
0 se τ > 0.
esta representada por linhas cheias. O sımbolo triangular representa a solucao aproximada
para a variavel de estado y e o sımbolo asterisco e usado para representar a variavel de
controle u.
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
τ
yexato
uexato
yaprox
uaprox
Figura 2.8: Solucao do problema 2.17.
Note que a variavel de controle u e uma funcao descontınua. A solucao aproxi-
mada assume erros maiores no ponto de descontinuidade e no extremo direito do intervalo.
No Apendice B sao apresentadas as listagens das rotinas em Matlab utilizadas.
Analise de erros. Adotando-se solucoes iniciais arbitrarias, foram resolvidas 30 instancias
do problema (2.19) para cada valor de N entre 10 e 100. Para cada uma destas instancias,
foi calculada a norma infinito do vetor diferenca entre a solucao exata e a solucao numerica
para a variavel de estado y e para a variavel de controle u. Ou seja, foram calculados os
erros dados por
erroy = ‖yexata − yaprox‖∞ e errou = ‖uexata − uaprox‖∞.
onde (yexata, uexata) e dado em (2.18) e (yaprox, uaprox) e a solucao numerica fornecida
Exemplos 40
pela rotina quadprog. Para cada um dos valores de N considerados, foi tomada a mediana
dos erros obtidos nas 30 instancias. Como ja ressaltado, tendo em vista a descontinuidade
da variavel de controle u, o erro da solucao aproximada na regiao de descontinuidade
e grande e foi desconsiderado na analise de erro. Os graficos da Figura 2.9 mostram a
variacao do erro, na escala logarıtmica para a solucao aproximada do estado e de controle,
respectivamente, em relacao ao numero de pontos de colocacao N entre 10 e 100. As linhas
cheias representam as curvas de regressao linear na escala logarıtmica cujos coeficientes
foram fornecidos pelo comando polyfit do Matlab. Enquanto o erro na variavel de estado
e da ordem de 10−3, o erro na variavel de controle e da ordem de 10−1.
10 20 30 40 50 60 70 80 90 100
10−3
N
log
(err
oy)
erroy
0.048N−0.891
10 20 30 40 50 60 70 80 90 10010
−2
10−1
N
log
(err
ou)
errou
0.472N−0.648
Figura 2.9: Variacao do erro das variaveis y e u na escala logarıtmica em relacao a N .
Conclusao
Discutiu-se neste trabalho o metodo de discretizacao pseudoespectral com pontos
de colocacao LGR. Esta tecnica e utilizada para encontrar solucoes numericas de proble-
mas de controle otimo sem restricoes (problema de Bolza tipo LQR), como foi publicado
em [22]. Apresentacao da tecnica segue a sequencia em que o metodo pseudoespectral
LGR e mostrado em [22], mas as demonstracoes sao proprias e estao dadas em notacao
tensorial. Para medir diretamente a precisao do metodo pseudoespectral LGR, e ne-
cessario ter a solucao exata do problema. A solucao exata do problema de Bolza LQR
de controle otimo foi encontrada utilizando o Princıpio de Maximo de Pontryagin. Para
comparar (positivamente) o desempenho do metodo pseudoespectral LGR, o problema
LQR de Bolza e resolvido tambem aplicando uma discretizacao de Euler da equacao de
estado e por uma quadratura tipo Euler do funcional. A discretizacao da equacao de
estado fornece uma forma de exprimir as variaveis de estado em termos das variaveis
de controle. A substituicao desta relacao na quadratura permite utilizar um metodo de
otimizacao quadratica da variavel de controle. O metodo pseudoespectral LGR e tambem
aplicado bem sucedidamente a problemas de controle otimo com restricoes nas variaveis
de estado e de controle (algo que nao e tentado em [22]).
Todos os codigos MATLAB estao incluıdos no Apendice desta dissertacao e serao
postados no site do Programa de Pos-graduacao em Matematica Aplicada da Universidade
Federal do Parana. Tambem no Apendice encontra-se uma revisao de conceitos referentes
a este trabalho.
Apendice A
Revisao de Conceitos
Neste capıtulo apresentam-se algumas definicoes, propriedades e ressaltam-se al-
guns resultados que podem ser uteis para o entendimento deste trabalho. As principais
referencias deste capıtulo sao [9, 10, 11, 12, 18, 25, 26, 27, 29].
A.1 Matriz
Nesta secao, que e baseada principalmente em [27], apresentam-se algumas de-
finicoes e resultados essenciais de matrizes que foram utilizados no trabalho. Inicialmente
apresenta-se a definicao de um tipo de matriz que e muito comum em problemas de
otimizacao.
Definicao A.1 [25, pag. 260] Seja A ∈ IRn×n uma matriz simetrica. Diz-se que A e
definida positiva quando xTAx > 0, para todo x ∈ IRn\0. Tal propriedade e denotada
por A > 0. Se xTAx ≥ 0, para todo x ∈ IRn, A e dita semidefinida positiva, fato este
denotado por A ≥ 0.
A.1.1 Produto direto de Matrizes
Em varias demonstracoes do Capıtulo 1 usa-se produto direto de matrizes. Nesta
secao apresentam-se a definicao de produto direto de matrizes, algumas propriedades deste
produto e um exemplo de como se faz este calculo.
Definicao A.2 Dadas as matrizes A ∈ IRm×n e B ∈ IRp×q define-se o produto direto de
A por B como a matriz C ∈ IRmp×nq de tal forma que
Cmp×nq = A⊗B =
a11B a12B · · · a1nB
a21B a22B · · · a2nB...
.... . .
...
am1B am2B · · · amnB
.
42
Apendice 43
Algumas propriedades interessantes do produto direto de matrizes:
(i) A⊗B 6= B ⊗ A;
(ii) Se u e v sao vetores, entao uT ⊗ v = v ⊗ uT = vuT ;
(iii) Se as dimensoes sao compatıveis
(A⊗B)(C ⊗D) = AC ⊗BD.
Exemplo A.3 Calcule o produto direto de A por B, sendo
A =
(1 2
1/2 4
)e B =
(1 2 4
6 8 10
).
Pela definicao e sendo A ∈ IR2×2 e B ∈ IR2×3 o resultado sera uma matriz
C ∈ IR4×6. Para construı-la, os seus elementos serao vistos como blocos, ou seja, Cij e o
bloco (i, j) de mesma dimensao de B. Assim,
C11 = a11B = 1 ∗
(1 2 4
6 8 10
), C12 = a12B = 2 ∗
(1 2 4
6 8 10
),
C21 = a21B = 1/2 ∗
(1 2 4
6 8 10
), C22 = a22B = 4 ∗
(1 2 4
6 8 10
).
E entao,
C =
1 2 4 2 4 8
6 8 10 12 16 20
1/2 1 2 4 8 16
3 4 5 24 32 40
.
A.1.2 Exponencial de matriz
Considere a ∈ IR. Ao resolver a equacao diferencial ordinaria (EDO) escalar
x(t) = ax,
obtem-se a seguinte formula,[10], para a solucao do problema com condicao inicial x(0) =
x0 ∈ IR:
x(t) = x0eat, t ∈ IR.
Na proxima secao, sera vista a equacao analoga quando a e uma matriz n × n com
coeficentes reais e x : IR → IRn. Para isto, precisa-se saber como calcular a exponencial
de uma matriz. As referencias sobre o assunto sao [10, 11, 29].
Apendice 44
Considera-se a funcao f : IR → IR+ tal que f(x) = ex. Essa funcao e chamada
de exponencial de x e tem a seguinte representacao, [29], por series de Taylor
ex =∞∑k=0
xk
k!.
Desta forma parece ser natural a definicao a seguir
Definicao A.4 [11, pag.330] Seja A ∈ IRn×n. A matriz definida pela serie
∞∑k=0
Ak
k!, sendo A0 = I
e denominada exponencial da matriz A. Para representa-la usa-se a notacao eA ou ainda
exp(A).
Para os problemas estudados neste trabalho, e suficiente saber como calcular uma expo-
nencial de matriz diagonalizavel.
Definicao A.5 [29, Def.3.1] Uma matriz A ∈ IRn×n e diagonalizavel quando existem
uma matriz invertıvel S, cujas colunas sao os autovetores de A e uma matriz diagonal D
formada pelos autovalores de A tal que
D = S−1AS.
Proposicao A.6 [29, Prop. 3.2] Se A e uma matriz diagonalizavel, D = S−1AS, entao
eA = eSDS−1
= SeDS−1.
Demonstracao. Multiplicando a igualdade D = S−1AS pela esquerda por S e depois pela
direita por S−1 tem-se que
A = SDS−1.
Alem disso,
(SDS−1)k = (SDS−1)(SDS−1) · · · (SDS−1)︸ ︷︷ ︸k vezes
= SDkS−1.
Logo, para n ∈ IN as somas parciais sao dadas por
Sn =n∑k=0
(SDS−1)k
k!=
n∑k=0
SDkS−1
k!= S
(n∑k=0
Dk
k!
)S−1.
Portanto,
eA = eSDS−1
= limn→∞
Sn = limn→∞
[S
(n∑k=0
Dk
k!
)S−1
]= S
(limn→∞
n∑k=0
Dk
k!
)S−1 = SeDS−1.
Apendice 45
Teorema A.7 [29, Teo. 3.4] Seja D = diag(α1, α2, . . . , αn), uma matriz diagonal n×n.
Entao eD = diag(eα1 , eα2 , . . . , eαn).
Demonstracao. Seja D uma matriz diagonal com entradas αi (i = 1, 2, . . . , n). Assim
para todo k ∈ IN,
Dk = diag(αk1, . . . , αkn).
Logo para cada n ∈ IN as somas parciais sao dadas por
Sn =n∑k=0
Dk
k!= diag
(n∑k=0
αk1k!, . . . ,
n∑k=0
αknk!
).
Assim
eD = limn→∞
Sn = limn→∞
(n∑k=0
Dk
k!
)= diag (eα1 , . . . , eαn) .
A.2 Solucao de um sistema de EDOs
O problema tratado no decorrer do trabalho esta sujeito a uma equacao diferen-
cial. O objetivo e resolver esse problema, para isso deve-se lembrar como e a solucao de
um sistema de EDOs de 1a e de 2a ordem.
Considere Ω ⊂ IRn e f : Ω → IRn uma aplicacao contınua. Seja I um intervalo
nao degenerado da reta, ou seja, um subconjunto conexo de IR nao reduzido a um ponto.
Definicao A.8 [30, pag 04.] Uma funcao diferenciavel φ : I → IRn chama-se solucao da
equacaodx
dt= f(x) (A.1)
no intervalo I se:
i) o grafico de φ em I esta contido em Ω e
ii)dφ
dt(t) = f(φ(t)) para todo t ∈ I. Se t e um ponto extremo do intervalo, a derivada
e a derivada lateral respectiva.
A equacao (A.1) chama-se equacao diferencial ordinaria de primeira ordem e e
denotada abreviamente por
x = f(x).
Sob hipoteses bem gerais sobre f , por exemplo se f e ∂f∂x
sao contınuas em Ω, existe uma
e so uma solucao φ de (A.1) num intervalo que contem t0 e tal que φ(t0) = x0. Uma tal φ
e chamada de solucao do problema com dados iniciais (t0, x0) para a equacao (A.1). Este
Apendice 46
problema pode ser denotado abreviadamente porx = f(x)
x(t0) = x0.
Para existencia e unicidade de solucoes de uma equacao diferencial consultar [30, pag. 12]
e [10].
Lembra-se agora solucoes de equacoes especıficas, que podem ser verificadas nos
Capıtulos 3 e 10 de [10].
A solucao do sistema de EDOs de primeira ordemx(t) = Mx(t)
x(0) = x0,
onde M ∈ IRn×n e x0 e o ponto inicial, e
x(t) = etMx0.
A solucao do sistema de EDOs de segunda ordemx(t) + bx(t) + cx(t) = 0
x(tf ) = k1, x(tf ) = k2
onde b, c ∈ IR e tf e o tempo final, e dada por
x(t) = c1eλ1t + c2e
λ2t
quando a equacao caracterıstica, λ2 + bλ+ c = 0, possui raızes reais distintas λ1 e λ2. As
constantes c1 e c2 sao encontradas de forma que as condicoes iniciais sejam satisfeitas.
A.3 Topicos de otimizacao contınua
Nesta secao sao introduzidos alguns topicos de otimizacao contınua que serao
uteis para o desenvolvimento deste trabalho.
Inicia-se com algumas ideias basicas sobre minimizacao de uma funcao. Em
seguida e apresentado o metodo de gradiente conjugado para minimizacao de uma funcao
quadratica. E encerra-se esta secao com um teorema de otimalidade de primeira ordem
para problemas de otimizacao com restricoes, o teorema de Karush-Kuhn-Tucker(KKT).
Para as demonstracoes dos teoremas e maiores detalhes, consultar [12].
Pode-se dizer que otimizacao consiste em encontrar pontos de mınimos ou de
maximos de uma funcao real sobre um conjunto Ω ⊂ IRn. Isto pode ser representado pelo
Apendice 47
seguinte problema
minimizar f(x)
sujeito a x ∈ Ω,(A.2)
onde f : IRn → IR e uma funcao contınua dita funcao objetivo e Ω ⊂ IRn e um conjunto
chamado de conjunto viavel. Quando Ω = IRn, o problema e dito irrestrito.
Definicao A.9 [12, pag. 25] Considere f : IRn → IR e x∗ ∈ Ω ⊂ IRn. Diz-se que x∗
e um minimizador local de f em Ω quando existe δ > 0 tal que f(x∗) ≤ f(x) para todo
x ∈ B(x∗, δ) ∩ Ω. Caso f(x∗) ≤ f(x) para todo x ∈ Ω, x∗ e dito minimizador global de f
em Ω.
Normalmente contenta-se com a obtencao de minimizadores locais do problema, onde esta
discussao sera focada.
O teorema a seguir nos da uma condicao necessaria para caracterizar um minimi-
zador local de um problema irrestrito. Os pontos que satisfazem tal condicao necessaria
sao ditos estacionarios.
Teorema A.10 Seja f : IRn → IR diferenciavel no ponto x∗ ∈ IRn. Se x∗ e um minimi-
zador local de f em IRn entao
∇f(x∗) = 0.
Demonstracao. [12, Teo. 2.9].
Em um problema de otimizacao e difıcil resolver, de forma direta, o sistema de
n equacoes e n incognitas dado por ∇f(x) = 0. Normalmente uma solucao e obtida por
meio de um processo iterativo. Considera-se entao um algoritmo que a partir de um
ponto inicial x0, gera uma sequencia de pontos (xk) desejando-se que os seus pontos de
acumulacao sejam estacionarios.
A.3.1 Metodo de gradientes conjugados
Apresenta-se nesta secao o metodo de gradientes conjugados [12, Sec.5.3] para
minimizacao de funcoes quadraticas.
Considere a funcao quadratica f : IRn → IR dada por
f(x) =1
2xTAx+ bTx+ c, (A.3)
com A ∈ IRn×n definida positiva, b ∈ IRn e c ∈ IR.
Segundo [12], uma funcao f quadratica com Hessiana definida positiva tem um
unico minimizador x∗, que e global e satisfaz, pelo Teorema A.10, a seguinte condicao
∇f(x∗) = Ax∗ + b = 0.
Apendice 48
Para entender o algoritmo sao necessarios alguns conceitos.
Definicao A.11 [12, pag. 64] Seja A ∈ IRn×n uma matriz definida positiva. Os vetores
d0, d1, ..., dk ∈ IRn\0 sao ditos A-conjugados se
(di)TAdj = 0,
para todos i, j = 0, 1, ..., k, com i 6= j.
Note que, no caso particular onde A e a matriz identidade, vetores A-conjugados
sao ortogonais no sentido usual.
Sabe-se que um metodo de gradientes conjugados minimiza um funcao quadratica
em IRn em no maximo n iteracoes.
Dado um conjunto qualquer de direcoes A-conjugadas d0, d1, ..., dn−1, define-se
uma sequencia finita da seguinte maneira: dado x0 ∈ IRn arbitrario, define-se para
k = 0, 1, ..., n− 1,
xk+1 = xk + tkdk, (A.4)
onde tk = argmint∈IR
f(xk + tdk). O escalar tk e o comprimento do passo na direcao dk a
partir do ponto corrente xk.
Como f e quadratica pode-se obter uma formula explıcita para tk. Para isso,
define-se φ : IR → IR por φ(t) = f(xk + tdk). Como, pela definicao, tk e o minimizador
de f(xk + tdk), tem-se que φ′(tk) = 0. Note que φ′(tk) = ∇f(xk + tkdk)Tdk. Por f ser
quadratica
∇f(xk + tkdk) = A(xk + tkd
k) + b = (Axk + b) + tkAdk = ∇f(xk) + tkAd
k.
Assim,
φ′(tk) =(∇f(xk)T + tkAd
k)dk = 0,
ou seja,
tk = −∇f(xk)Tdk
(dk)TAdk.
Agora mostra-se como gerar as direcoes conjugadas.
Dado x0 ∈ IRn, defina d0 = −∇f(x0) e, para k = 0, 1, ..., n− 2,
dk+1 = −∇f(xk+1) + βkdk,
onde xk+1 e dado por (A.4) e βk e tal que dk e dk+1 sejam A-conjugados, ou seja,
(dk)TA(−∇f(xk+1) + βkdk) = (dk)TAdk+1 = 0.
Apendice 49
Isolando βk da igualdade acima
βk =(dk)TA∇f(xk+1)
(dk)TAdk+1.
Agora com todas as ferramentas necessarias, apresenta-se o algoritmo de gradi-
entes conjugados.
Algoritmo A.1 Gradientes conjugados para funcoes quadraticas
Dado: x0 ∈ IRn, faca d0 = −∇f(x0)
k = 0
repita enquanto ∇f(xk) 6= 0
tk = −∇f(xk)Tdk
(dk)TAdk
xk+1 = xk + tkdk
βk =(dk)TA∇f(xk+1)
(dk)TAdk
dk+1 = −∇f(xk+1) + βkdk
k = k + 1
O algoritmo de gradientes conjugados neste formato e aplicado apenas para
funcoes quadraticas.
A.3.2 O Teorema de Karush-Kuhn-Tucker
O objetivo desta secao e apresentar as condicoes de otimalidade ou sejam, as
condicoes de KKT (Karush-Kuhn-Tucker), para o problema geral de otimizacao, que
consiste emminimizar f(x)
sujeito a g(x) ≤ 0
h(x) = 0,
(A.5)
onde f : IRn → IR, g : IRn → IRp e h : IRn → IRm sao funcoes continuamente diferenciaveis.
O problema (A.5) e o problema (A.2) em que o conjunto viavel e dado por
Ω = x ∈ IRn | g(x) ≤ 0, h(x) = 0.
Definicao A.12 [12, pag. 97] Seja x ∈ Ω. Uma restricao de desigualdade gi e dita ativa
em x, se gi(x) = 0. Caso gi(x) < 0, diz-se que gi e inativa em x.
Denota-se por A(x) o conjunto de ındices das restricoes de desigualdade ativas em um
ponto viavel x, isto e,
A(x) = i | gi(x) = 0.
Apendice 50
Teorema A.13 (KKT) Seja x∗ ∈ Ω um minimizador local do problema (A.5) e suponha
que o conjunto dos gradientes das restricoes de desigualdade ativas e das restricoes de
igualdade sao linearmente independentes, isto e, ∇gi(x∗)i∈A(x∗) ∪ ∇hj(x∗)j∈1,...,m e
LI. Entao existem µ∗ ∈ IRp e λ∗ ∈ IRm tais que
−∇f(x∗) =
p∑i=1
µ∗i∇gi(x∗) +m∑j=1
λ∗j∇hi(x∗),
µ∗i ≥ 0, i = 1, ..., p,
µ∗i gi(x∗) = 0, i = 1, ..., p.
Demonstracao. [19, Teo. 12.1]
Os vetores λ∗ e µ∗ sao chamados de multiplicadores de Lagrange.
Note que, o resultado do teorema diz que o oposto do gradiente da funcao objetivo
pode ser expresso pela combinacao linear dos gradientes das restricoes do problema.
Interpretacao geometrica
Para representar geometricamente o que o teorema de KKT afirma, considere o
seguinte problema
minimizar f(x) = (x1 − 2)2 + (x2 − 1)2
sujeito a g1(x) = x1 + x2 − 2 ≤ 0
g2(x) = x21 − x2 ≤ 0
x1 ≥ 0
x2 ≥ 0,
(A.6)
onde f : IR2 → IR, g1 : IR2 → IR e g2 : IR2 → IR.
Primeiramente note que o ponto x∗ =
(1
1
)e solucao global do problema.
Na figura abaixo, estao representados algumas curvas de nıvel da funcao, o con-
junto viavel e a solucao x∗.
Apendice 51
Figura A.1: Ilustracao do problema (A.6).
Na figura abaixo, estao representados os gradientes das restricoes ativas no ponto(1
1
)e o oposto do gradiente da funcao objetivo.
Figura A.2: Ilustracao da relacao dos gradientes do problema (A.6).
Algebricamente tem-se que ∇f(x∗) =
(−2
0
), ∇g1(x∗) =
(1
1
), ∇g2(x∗) =(
2
−1
). Assim
−∇f(x∗) =2
3∇g1(x∗) +
2
3∇g2(x∗),
ou seja, o oposto do gradiente da funcao objetivo no ponto x∗ e combinacao linear, com
coeficientes µ = λ = 23, dos gradientes das restricoes no ponto x∗.
Apendice B
Codigos em Matlab
52
Cálculo dos Pontos de Colocação e pesos de Quadratura
function [x,w]=pontos_lgr_w(n)
%======================================================================== % Karla Arsie - 2012 % Calcula os n pontos LGR como raizes de P_n-1(x)+P_n(x). % Polinomio de Legendre P_n = P(:,n+1) % % Foi baseado no programa disponivel em: % http://www.mathworks.com/matlabcentral/fileexchange/4775-legende-gauss-
lobatto-nodes-and-weights/content/lglnodes.m. % conforme descricao abaixo e na referencia % F. B. Hildebrand , "Introduction to Numerical Analysis," Section 8.11 % Dover 1987 %======================================================================== % lglnodes.m % % Computes the Legendre-Gauss-Lobatto nodes, weights and the LGL Vandermonde % matrix. The LGL nodes are the zeros of (1-x^2)*P'_N(x). Useful for numerical % integration and spectral methods. % % Reference on LGL nodes and weights: % C. Canuto, M. Y. Hussaini, A. Quarteroni, T. A. Tang, "Spectral Methods % in Fluid Dynamics," Section 2.3. Springer-Verlag 1987 % % Written by Greg von Winckel - 04/17/2004 % Contact: [email protected] % % Truncation + 1 % N1=n; % N=n-1;
N1=n; N=n-1;
%Ponto inicial para pontos LGR x=-cos(2*pi*(0:N)/(2*N+1))';
% Matriz de Vandermonde P=zeros(N1,N1+1);
% Encontra P_(N), utilizando a relação de recorrencia
% Calcula suas derivadas de primeira e segunda ordem
% X: atualização usando o método de Newton-Raphson xold=2;
free=2:N1;
while max(abs(x-xold))> eps
xold=x;
P(1,:)=(-1).^(0:N1);
P(free,1)=1; P(free,2)=x(free);
53
Apendice 54
for k=2:N1 P(free,k+1)=( (2*k-1)*x(free).*P(free,k)-(k-1)*P(free,k-1) )/k; end
x(free)=xold(free)-((1-xold(free))/N1).*(P(free,N1)+P(free,N1+1))... ./(P(free,N1)-P(free,N1+1)); end
P=P(1:N1,1:N1);
% Calculo dos pesos w=zeros(N1,1); w(1)=2/N1^2; w(free)=(1-x(free))./(N1*P(free,N1)).^2;
Apendice 55
Matriz de Diferenciação
function D=matriz_dif(n,x)
%============================================================ % Karla, Miguel, Elizabeth % 2012 % Fornece a matriz de diferenciacao % Adaptacao do programa legsrddiff escrito por Wang Li-Lian disponivel em: % http://www1.spms.ntu.edu.sg/~lilian/bookcodes/legen/legsrddiff.m %============================================================= % D=legsrddiff(n,x) returns the first-order differentiation matrix of size % n by n, associated with the Legendre-Gauss-Radau points x, which may be
%computed by % x=legsrd(n). Note: x(1)=-1. % See Page 110 of the book: J. Shen, T. Tang and L. Wang, Spectral Methods: % Algorithms, Analysis and Applications, Springer Series in Compuational % Mathematics, 41, Springer, 2011. % Use the function: lepoly() % Last modified on August 31, 2011 %=============================================================
if n==0, D=[]; return; end; xx=x; nx=size(x);
% Encontra L_n-1+L_n e sua derivada de primeira ordem [dy1,y1]=lepoly(n-1,xx); [dy2,y2]=lepoly(n,xx); dy=dy1+dy2; y=y1+y2; if nx(2)>nx(1), y=y'; dy=dy'; xx=x'; end; %% y é o vetor coluna de L_n-1(x_k) D=(xx./dy)*dy'-(1./dy)*(xx.*dy)'; %% encontra dy(x_j) (x_k-x_j)/dy(x_k); % 1/d_kj for k not= j (see (3.204)) D=D+eye(n); % add a matriz identidade para que 1./D pode ser
% operada D=1./D; D=D-eye(n); xx=xx(2:end); D=D+diag([-(n+1)*(n-1)/4; xx./(1-xx.^2)+ n*y1(2:end)./((1-xx.^2).*dy(2:end))]);
%atualiza as diagonais return;
Apendice 56
function [varargout]=lepoly(n,x)
% lepoly Legendre polynomial of degree n % y=lepoly(n,x) is the Legendre polynomial % The degree should be a nonnegative integer % The argument x should be on the closed interval [-1,1]; % [dy,y]=lepoly(n,x) also returns the values of 1st-order % derivative of the Legendre polynomial stored in dy % Last modified on August 30, 2011 % Verified with the chart in http://keisan.casio.com/has10/SpecExec.cgi
if nargout==1, if n==0, varargout1=ones(size(x)); return; end; if n==1, varargout1=x; return; end; polylst=ones(size(x)); poly=x; % L_0(x)=1, L_1(x)=x for k=2:n, % Three-term recurrence relation: polyn=((2*k-1)*x.*poly-(k-1)*polylst)/k;
% kL_k(x)=(2k-1)xL_k-1(x)-(k-1)L_k-2(x) polylst=poly; poly=polyn; end; varargout1=polyn; end;
if nargout==2, if n==0, varargout2=ones(size(x)); varargout1=zeros(size(x));
return;
end; if n==1, varargout2=x; varargout1=ones(size(x));
return;
end;
polylst=ones(size(x)); pderlst=zeros(size(x));poly=x; pder=ones(size(x)); % L_0=1, L_0'=0, L_1=x, L_1'=1 for k=2:n, % Three-term recurrence relation: polyn=((2*k-1)*x.*poly-(k-1)*polylst)/k;
% kL_k(x)=(2k-1)xL_k-1(x)-(k-1)L_k-2(x) pdern=pderlst+(2*k-1)*poly; % L_k'(x)=L_k-2'(x)+(2k-1)L_k-1(x) polylst=poly; poly=polyn; pderlst=pder; pder=pdern; end; varargout2=polyn; varargout1=pdern; end;
return
Apendice 57
Solução Exata e pelo método Pseudoespectral do Problema do Exemplo 1
go_ex1
%=====================================================
% Dissertacao de Mestrado de Karla Arsie % 17/08/2012 % Problema: maximizar -int_0^5(x^2(t)+u^2(t)) dt % sujeito dx/dt=x+u % x(0)=sqrt2,x(5)=1 % Primeiro Exemplo do Capitulo 2 da dissertacao %=====================================================
clear all
N=60; % Numero de pontos de colocacao
X1=sqrt(2); Xend=1; % Condicoes de contorno
% Solucao aproximada [tau,VETORW]=pontos_lgr_w(N); % retorna os pontos de colocacao e os pesos da
discretizacao.
MATRIZD=matriz_dif(N,tau); %retorna matriz de diferenciacao
%ponto inicial x0=X1*rand(N,1); x0(1)=X1; u0=rand(N,1); xu0 = [x0;u0];
% Dados do problema para ser resolvido pelo quadprog nas variaveis
%(x1,..xN,u1,..uN) xu0 = [x0;u0]; H=diag([VETORW;VETORW]); f=zeros(2*N,1); id=eye(N,N); Aeq = [(2*MATRIZD/5-id) -id]; Beq =zeros(N,1); % Limites apos varios testes ----------------- LB = [X1;zeros(N-2,1);1;-inf*ones(N,1)]; UB = [X1;inf*ones(N-2,1);1;inf*ones(N,1)];
[xu,fval,EXITFLAG] = quadprog(H,f,[],[],Aeq,Beq,LB,UB,xu0);
x = [X1;xu(2:N);Xend]; u = xu(N+1:2*N);
Apendice 58
% --------------------------------------------------------- % Solucao exata s=X1; x0=s; p0=(4*s-exp(5*s)*(2*s+4)-exp(-5*s)*(4-2*s))/(exp(5*s)-exp(-5*s)); x_exata=(exp(s*((5/2)*tau+5/2))*((1-exp(-5*s))/(exp(5*s)-exp(-5*s)))+ exp(-
s*((5/2)*tau+5/2))*(s-(1-exp(-5*s))/(exp(5*s)-exp(-5*s))))'; %no intervalo [-1,1] y_exata=x_exata.^2; p=(((s/2)*exp(s*((5/2)*tau+5/2))-(s/2)*exp(-s*((5/2)*tau+5/2)))*x0+(((2-
s)/(4))*exp(s*((5/2)*tau+5/2))+((2+s)/(4))*exp(-s*((5/2)*tau+5/2)))*p0)';
%no intervalo [-1,1] u_exata=p/2;
% Graficos %----------------------------------------------------------- figure(1) clf axis([-1 5 -4 2]) %solucao exata plot(tau,y_exata,'b-',tau,u_exata,'r-') hold on xlabel('\tau','fontsize',20) % Solucao aproximada --------------------------------------- plot([tau;1],x.*x,'^b') plot(tau(1:N-1),u(1:N-1),'*r') legend('y_exato','u_exato','y_aprox','u_aprox','Location','SouthEast')
%----------------------------------------------------------- figure(2) clf axis([-1 5 -4 2]) %solucao exata plot(tau,y_exata,'b-',tau,u_exata,'r-') hold on xlabel('\tau','fontsize',20) % Solucao aproximada --------------------------------------- plot([tau;1],x.*x,'^b') plot(tau(1:N-1),u(1:N-1),'*r') euler_ex1 legend('y_exato','u_exato','y_aprox','u_aprox','y_Euler',
'u_Euler','Location','SouthEast') %-----------------------------------------------------------
% Erro na restricao %ndifEx=norm(MATRIZD*y_exata-u_exata); ndifAp=norm(2*MATRIZD/5*x(1:N)-x(1:N)-u); fprintf(' %4d | %7.3f | %2d \n',N,ndifAp,EXITFLAG)
Apendice 59
Solução pelo método de Euler do Problema do Exemplo 1 - euler_ex1
%======================================================= % Dissertacao de Mestrado de Karla Arsie % 25/05/2012 % Problema: maximizar -int_0^5(x^2(t)+u^2(t)) dt % sujeito dx/dt=x+u % x(0)=sqrt2,x(5)=1 % Primeiro Exemplo do Capitulo 2 da dissertacao %resolvido pelo método de Euler %=====================================================
n=N; M=10^3; h=5/(n-1);
[Q b c]=gra(n,X1); % gera a quadratica a ser minimizada
% somando o termo M*(Xend-1)^2 beta=X1*(1+h)^(n-1)-Xend; c=c+M*beta^2; for i=1:n-1 b(i)=b(i)+2*M*beta*h*(1+h)^(n-i-1); for j=1:n-1 Q(i,j)=Q(i,j)+2*M*h^2*(1+h)^(2*n-j-i-2); end end
u0=rand(n-1,1);
% Solucao aproximada Euler ---------------- u_euler=GC(c,b,Q,u0); % usa Gradiente Conjugado para minimizacao
x_euler=zeros(n,1); x_euler(1)=X1; for k=2:n x_euler(k)=(1+h)^(k-1)*X1; for j=1:k-1 x_euler(k)=x_euler(k)+h*(1+h)^(k-j-1)*u_euler(j); end end y_euler=x_euler.^2; % plotado y=x^2
t=-1:2/(N-1):1;t=t'; plot(t,y_euler,'og','linewidth',2) plot(t(1:N-1),u_euler,'dk','linewidth',2)
Apendice 60
Funções Auxiliares
function [Q b c]=gra(n,x1)
%Gera c,b e Q da função J=1/2sum(xk^2+uk^2)=c+b^Tu+u^TQu/2;
v=zeros(n,n-1); for i=1:n-1 for j=i+1:n v(j,i)=1; end end
%encontrar c c=zeros(n-1,1); h=1/(n-1); h=5*h; s=0; x1=sqrt(2); for k=1:n-1 s=s+(1+h)^(2*k-2)*x1^2; end c=h/2*s;
%encontrar b z=zeros(n-1,n); for k=2:n-1 for i=1:k-1 z(i,k)=h^2*(1+h)^(2*k-i-2)*x1; end end
for i=1:n-1 s=0; for k=1:n s=s+z(i,k)*v(k,i); end b(i)=s; end b=b';
%encontrar Q w=zeros(n-1,n-1,n); for k=2:n-1 for i=1:k-1 for j=1:k-1 w(i,j,k)=h^3*(1+h)^(2*k-i-j-2); end end end
Apendice 61
for i=1:n-1 for j=1:n-1 s=0; if i==j
for k=1:n s=s+w(i,i,k)*v(k,i); end s=s+h; elseif i>j for k=1:n s=s+w(i,j,k)*v(k,i); end else for k=1:n s=s+w(i,j,k)*v(k,j); end end Q(i,j)=s; end end
end
Apendice 62
% Gradiente Conjugado % função quadratica f(x)=c+bx+x'Qx/2
function x=GC(c,b,Q,x0); x=x0; n=length(x0)+1; % dimensao do problema. gradfx=b+Q*x; fx=c+b'*x+x'*Q*x/2;
d=-gradfx; % direcao d0 ngradf=norm(gradfx); % norma do grad
k=1; kmax=3*n; %n; epsilon=1e-6;
while (ngradf>=epsilon) & (k<=kmax) Qd=Q*d; dQd=d'*Qd; alfa=-(d'*gradfx)/dQd; x=x+alfa*d;
gradfx=b+Q*x; ngradf=norm(gradfx); fx=c+b'*x+x'*Q*x/2;
beta=(d'*Q*gradfx)/dQd;
d=-gradfx+beta*d; if k>n d=-gradfx; end k=k+1; end
if ngradf>epsilon | k>=kmax disp('GC - Nao resolveu') fprintf('norma do gradiente: %12f \n',ngradf) pause(0.1) end
Apendice 63
Solução pelo método Pseudoespectral do Problema do Exemplo 2 - go_ex2
%=========================================================== % Dissertacao de Mestrado de Karla Arsie % 10/03/2013 % Problema: minimizar int_0^1(u^2(t)-y^3(t)) dt % sujeito dy/dt=yu % y(0)=1,y(1)=1 % -1=<u(t)=<0 % Primeiro Exemplo Nao linear do Capitulo 2 da dissertacao %========================================================== clear all
global VETORW MATRIZD Y1 Yend
N=60; % numero de pontos de colocacao
Y1=1; Yend=1; % Condicoes de contorno
% Solucao aproximada [tau,VETORW]=pontos_lgr_w(N); % pontos de colocacao e pesos da discretizacao
MATRIZD=matriz_dif(N,tau)/2; %retorna matriz de diferenciacao
%ponto inicial y0=Y1*rand(N,1); y0(1)=Y1; u0=rand(N,1); yu0 = [y0;u0];
% Dados do problema para ser resolvido pelo fmincon A = []; B = []; Aeq = []; Beq = []; LB = [ones(N,1);-ones(N,1)]; UB = [Inf*ones(N,1);zeros(N,1)];
% Resolucao pelo fmincon OPTIONS = OPTIMSET('Display','iter','GradObj','on'); [yu,fval,EXITFLAG] = FMINCON(@(yu0) funcao_ex2(yu0),yu0,A,B,Aeq,Beq,LB,UB,@(yu0)
nonlin_ex2(yu0),OPTIONS);
y = [Y1;yu(2:N);Yend]; u = yu(N+1:2*N);
%plota graficos da solucao aproximada figure(3) hold on plot([tau;1],y,'^b') plot(tau,u,'*r') legend('y','u',4) xlabel('\tau','fontsize',16)
Apendice 64
function [J,grad_J]=my_functional_J(YU0)
%gera a funcao a ser minimziada
global VETORW MATRIZD Y1 Yend
n0 = (length(YU0))/2; % assumed to be n Y0 = YU0(1:n0); % size n
U0 = YU0(n0+1:2*n0); % size n
J=VETORW'*(Y0.^3+U0.^2);
grad_J =[3*(Y0.^2).*VETORW;2*U0.*VETORW];
end
function [c,ceq]=nonlin(YU0) %gera as restricoes do problema
global VETORW MATRIZD Y1 Yend
% YU0 is [y(2:n-1),u]
c=[]; n0 = length(YU0)/2; % assumed to be n
Y0 = max(0,YU0(1:n0)); %evitar Y <0 % size n
%SY0 = diag(sqrt(Y0)); % size n x n U0 = YU0(n0+1:2*n0); % size n
ceq = MATRIZD*Y0 - Y0.* U0; % size n
end
Apendice 65
Solução pelo método Pseudoespectral do Problema do Exemplo 3 - go_ex3
%========================================================================== % go_ex3 % Segundo Exemplo - Problema nao linear % Referencia: J. L. Troutman. Variational Calculus and Optimal Control: % Optimization with Elementary Convexity. Springer-Verlag, New York, 1996. % Solucao exata com t in [0,2]: y(t) =1-t se t<=1, 0 c.c. % u(t)=1 se t<=1, 0 c.c. % Solucao exata com tau in [-1,1]: y(tau) = -tau se tau <= 0, 0 c.c. % u(tau) = 1 se tau <= 0, 0 c.c. % Dissertacao de Mestrado - Karla Arsie % Marco/2013 %========================================================================== clear all
N=60; % Numero de pontos de colocacao
Y1=1; Yend=0;% Condicoes de contorno
% Solucao aproximada [tau,VETORW]=pontos_lgr_w(N); % retorna os pontos de colocacao
%e os pesos da discretizacao MATRIZD=matriz_dif(N,tau); % retorna a matriz de diferenciacao
%ponto inicial y0=[Y1*rand(N,1);Yend]; y0(1)=Y1; u0=rand(N,1);
% Dados do problema para ser resolvido pelo quadprog nas variaveis
(y1,..yN,yN+1,u1,..uN) yu0 = [y0;u0]; H=2*diag([VETORW;1;zeros(N,1)]); f=zeros(2*N+1,1); Aeq = [MATRIZD zeros(N,1) -eye(N,N)]; Beq =zeros(N,1); LB = [1;-inf*ones(N,1);-ones(N,1)]; UB = [1;inf*ones(N,1);ones(N,1)];
% Resolucao pelo quadprog OPTIONS = OPTIMSET('Algorithm','interior-point-convex',
'Display','iter','GradObj','on'); yu = quadprog(H,f,[],[],Aeq,Beq,LB,UB,yu0,OPTIONS);
y = [Y1;yu(2:N);Yend]; u = yu(N+2:end);
Apendice 66
%plota graficos para comparar solucao exata e solucao aproximada figure(1); clf % Solucao Exata plot([-1 0 1],[1 0 0],'b-','linewidth',2) hold on plot([-1 0],[-1 -1],'r-','linewidth',2) plot([tau;1],y,'^b') plot(tau,u,'*r') axis([-1 1 -1 1]) xlabel('\tau','fontsize',24) grid on legend('y_exato','u_exato','y_aprox','u_aprox',4) plot([0 1],[0 0],'r-') plot([tau;1],y,'^b')
% Solucao exata nos pontos de colocacao para analise do erro y_exata=-tau.*(tau<=0); u_exata=-(tau<=0);
Referencias Bibliograficas
[1] D. N. Arnold. A concise introduction to numerical analysis. School of Mathematics,
University of Minnesota, 1991.
[2] D. Benson. A gauss pseudospectral transcription for optimal control. Technical
report, Massachusetts Institute of Technology, 2005.
[3] J. T. Betts. Survey of numerical methods for trajectory optimization. Technical
report, Journal of Guidance, Control and Dynamics, Washington, Marco-Abril 1998.
[4] J. T. Betts. Practical Methods for Optimal Control and Estimation Using Nonlinear
Programming. Society for Industrial & Applied Mathematics; 2nd edition, Washing-
ton, 2009.
[5] G.A. Bliss. The problem of Bolza in the calculus of variations. Annals of Mathema-
tics(2), 33:261–274, 1932.
[6] M. L. Boas. Mathematical Methods in the Physical Sciences. John Wiley and Sons,
New York, 1983.
[7] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Tang. Spectral Methods in
Fluid Dynamics. Springer-Verlag, 1987.
[8] D. Garg e A.V. Rao C.L. Darby. Costate estimation using multiple-interval pseudos-
pectral methods. Journal of spacecraft and rockets, 48(5), 2011.
[9] P. J. Davis and P. Rabinowitz. Methods of Numerical Integration. Dover Publications,
INC, New York, 2a edition, 1984.
[10] C. I. Doering and A. O. Lopes. Equacoes Diferenciais Ordinarias. IMPA, Rio de
Janeiro, RJ, 2005.
[11] J. Baumeister e A. Leitao. Introducao a teoria de controle e programacao dinamica.
Instituto de Matematica Pura e Aplicada, Rio de Janeiro, Brasil, 1nd edition, 2008.
[12] A. A. Ribeiro e E. W. Karas. Um curso de Otimizacao. Cengage Learning, 2013. A
aparecer.
67
Referencias Bibliograficas 68
[13] S. Islam e G. Saha. Applications of gauss-radau and gauss-lobatto numerical inte-
grations over a four node quadrilateral finite element. Bangladesh J. Sci. Ind. Res.,
43(3):377–386, 2008.
[14] F. Fahroo e I. M. Ross. Pseudospectral methods for infinite-horizon nonlinear optimal
control problems. Journal of Guidance, Control, and Dynamics,, 31:927–936, 2008.
[15] S. Kameswaran e I. T. Biegler. Convergence rates for direct transcription of optimal
control problems using collocation at radau points. Computational Optimization and
Applications, 41:81–126, 2008.
[16] B. A. Finlayson e L. E. Scriven. The method of weighted residuals - a review. Applied
Mechanics Reviews, 19(9):735–748, 1966.
[17] O. Stryk e R. Bulirsch. Direct and indirect methods for trajectory optimization.
Annals of Operations Research, 37:357–373, 1992.
[18] N. Barron e R. Jensen. The pontryagin maximum principle from dynamic program-
ming and viscosity solutions to first-order partial differential equations. Trans. Amer.
Math. Soc., 298:635–641, 1986.
[19] J. Nocedal e S. J. Wright. Numerical Optimization. Springer-Verlag, Springer Series
in Operations Research, USA, 2nd edition, 2006.
[20] F. Lewis e V. Syrmos. Optimal Control. John Wiley & Sons, INC, 3nd, New Jersey,
USA, 1995.
[21] M. Kazemi e M. Razzaghi G. Elnagar. The pseudospectral legendre method for
discretizing optimal control problems. IEEE Transactions on Automatic Control,
40:1793–1796, 1995.
[22] D. Garg, M. A. Patterson, C. Francolin, C. L. Darby, G. T. Huntington, W. W. Hager,
and A. V. Rao. Direct trajectory optimization and costate estimation of finite-horizon
and infinite-horizon optimal control problems using a Radau pseudospectral method.
Comput. Optim. Appl., 49:335–358, 2011.
[23] F. B. Hildebrand. Introduction to Numerical Analysis. Dover Publications, Dover,
2a edition, 1987.
[24] D. Kraft. On converting optimal control problems into nonlinear programming codes.
NATO ASI Series, Computational Mathematical Programming, ed. K. Schittkowski,
Springer, F15:261–280, 1985.
[25] S. J. Leon. Algebra linear com Aplicacoes. Editora S.A., Rio de Janeiro, Brasil, 4nd
edition, 1998.
Referencias Bibliograficas 69
[26] E. L. Lima. Curso de Analise volume 2. Instituto de Matematica Pura e Aplicada,
Rio de Janeiro, Brasil, 1981.
[27] A. Quarteroni, R. Sacco, and F. Saleri. Numerical Mathematics. Springer, New York,
2000.
[28] J. Stoer e K. H. Well (eds.) R. Bulirsch, A. Miele. Optimal control - calculus of
variations, optimal control theory and numerical methods. International Series of
Numerical Mathematics, 111:129–143, 1993.
[29] F. A. Silva. Sistemas de equacoes diferenciais e exponencial de matrizes. Cadernos
PET - Matematica, 1:169–195, 2007.
[30] J. Sotomayor. Licoes de Equacoes Diferenciais. Projeto Euclides, IMPA, Rio de
Janeiro, Brasil, 1979.
[31] J. L. Troutman. Variational Calculus and Optimal Control: Optimization with Ele-
mentary Convexity. Springer-Verlag, New York, 1996.
[32] G. V. Winckel. Codigo em Matlab. Disponıvel em
http://www.mathworks.com/matlabcentral/fileexchange/4775-legende-gauss-
lobatto-nodes-and-weights/content/lglnodes.m, consultado em junho de 2012.
