Um tutorial sobre métodos pseudo-espectrais para controle ótimo computacional

A tutorial on pseudospectral methods for computational optimal control

OTIMIZAÇÃO EM CONTROLE E AUTOMAÇÃO

Um tutorial sobre métodos pseudo-espectrais para controle ótimo computacional

A tutorial on pseudospectral methods for computational optimal control

Victor M. Becerra^I; Roberto Kawakami Harrop Galvão^II

^IUniversity of Reading School of Systems Engineering Reading RG6 6AY, UK. [e-mail]

^IIInstituto Tecnológico de Aeronáutica Divisão de Engenharia Eletrônica 12228-900, São José dos Campos, SP, Brasil. [e-mail]

RESUMO

Este artigo é um tutorial introdutório sobre controle ótimo pseudo-espectral. Em métodos pseudo-espectrais, uma função é aproximada como uma combinação linear de funções de base suaves, tipicamente escolhidas como polinômios de Legendre ou Chebyshev. A colocação de equações algébrico-diferenciais é realizada em pontos de colocação ortogonal, que são selecionados de modo a minimizar o erro de interpolação. Métodos pseudo-espectrais discretizam o problema de controle ótimo original de modo a convertê-lo em um problema de programação não-linear. Um otimizador numérico é então empregado para obter soluções localmente ótimas.

Este artigo também descreve sucintamente a funcionalidade e a implementação de um pacote computacional de código aberto escrito em C++ chamado PSOPT. Tal pacote emprega métodos de discretização pseudo-spectrais para resolver problemas de controle ótimo com múltiplas fase. O PSOPT permite a utilização de métodos de Legendre ou Chebyshev, e possui características úteis tais como diferenciação automática, detecção de esparsidade e escalonamento automático. O uso de métodos pseudo-espectrais é ilustrado em dois problemas retirados da literatura de controle ótimo computacional.

Palavras-chave: Controle Ótimo, Métodos Pseudo-Espectrais, Programação Não-Linear, Planejamento de Trajetórias.

ABSTRACT

This paper is a tutorial introduction to pseudospectral optimal control. With pseudospectral methods, a function is approximated as a linear combination of smooth basis functions, which are often chosen to be Legendre or Chebyshev polynomials. Collocation of the differential-algebraic equations is performed at orthogonal collocation points, which are selected to yield interpolation of high accuracy. Pseudospectral methods directly discretize the original optimal control problem to recast it into a nonlinear programming format. A numerical optimizer is then employed to find approximate local optimal solutions.

The paper also briefly describes the functionality and implementation of PSOPT, an open source software package written in C++ that employs pseudospectral discretization methods to solve multi-phase optimal control problems. The software implements the Legendre and Chebyshev pseudospectral methods, and it has useful features such as automatic differentiation, sparsity detection, and automatic scaling. The use of pseudo-spectral methods is illustrated in two problems taken from the literature on computational optimal control.

Keywords: Optimal Control, Pseudospectral Methods, Nonlinear Programming, Trajectory Planning.

1 INTRODUÇÃO

Métodos pseudo-espectrais foram originalmente desenvolvidos para a solução de equações diferenciais parciais, tendo se tornado uma ferramenta amplamente utilizada em fluidodinâmica computacional (Canuto et al., 1988), (Canuto et al., 2007). Nos últimos 15 anos, tais métodos também passaram a ser explorados para solução numérica de problemas de controle ótimo (Elnagar et al., 1995), (Fahroo and Ross, 2001), (Fahroo and Ross, 2002), (Ross and Fahroo, 2004), (Kang and Bedrossian, 2007), (Gong et al., 2008). Ao contrário de técnicas de diferenças finitas, que aproximam as derivadas de uma função empregando informações locais, métodos pseudo-espectrais aproximam as funções de estado e controle como uma combinação linear de funções de base globais, tipicamente escolhidas como polinômios de Legendre ou Chebyshev. A colocação de equações algébrico-diferenciais é realizada em pontos de colocação ortogonal, que são selecionados de modo a minimizar o erro de interpolação. Métodos pseudo-espectrais possuem taxa de convergência exponencial (também dita "espectral") e praticamente não apresentam erros dissipativos ou de dispersão. Além disso, permitem a obtenção de boas aproximações, mesmo empregando grades com um número relativamente pequeno de pontos (Trefethen, 2000). Nos casos em que a colocação global não se mostra apropriada (por exemplo, quando a solução pode exibir discontinuidades), técnicas pseudo-espectrais de múltiplos domínios podem ser utilizadas. Para isso, cada problema é dividido em um conjunto de sub-intervalos, dentro dos quais é feita uma colocação global (Canuto et al., 2007).

Métodos pseudo-espectrais discretizam diretamente o problema de controle ótimo original, de modo a convertê-lo em um problema de programação não-linear (PPNL). Tal problema pode então ser resolvido com o uso de otimizadores numéricos de modo a obter soluções aproximadas localmente ótimas. Pode-se mostrar que métodos pseudo-espectrais são particularmente úteis para aproximar funções suaves, bem como operações de integração e diferenciação (Canuto et al., 2006), (Trefethen, 2000), que são elementos importantes em problemas de controle ótimo. Para fins de diferenciação, as derivadas dos estados nos pontos de discretização são calculadas multiplicando-se uma matriz constante de diferenciação pelos valores dos estados nos pontos em questão. Com isso, as equações diferenciais do problema de controle ótimo são substituídas por um conjunto de equações algébricas. Empregando as conhecidas regras de quadratura de Gauss, a integração do funcional de custo pode ser aproximada por uma soma ponderada dos valores do integrando nos pontos de discretização. Vale salientar que restrições nos estados e controles podem ser facilmente incluídas na formulação.

O método pseudo-espectral de Chebyshev para problemas de controle ótimo foi originalmente proposto em (Vlassenbroeck, 1988) e (Vlassenbroeck and Dooren, 1988). Uma variante para otimização de trajetórias foi posteriormente apresentada por Fahroo and Ross (2002). Já o método pseudo-espectral de Legendre foi proposto por Elnagar et al. (1995), tendo depois recebido contribuições de vários outros pesquisadores. Por exemplo, em (Kang et al., 2007) demonstra-se que uma seqüência de soluções obtidas com o método pseudo-espectral de Legendre converge para a solução ótima do problema de controle original. O tratamento de problemas envolvendo múltiplas fases é discutido em (Ross and Fahroo, 2004). Uma extensão para o caso de horizonte infinito é proposta em (Fahroo and Ross, 2008). Uma aplicação que ganhou notoriedade consistiu na geração de trajetórias em tempo real para reorientação da Estação Espacial Internacional (Kang and Bedrossian, 2007), (Bedrossian et al., 2009).

Os métodos de Chebyshev e Legendre fazem uso de polinômios ortonormais para a discretização de problemas de controle ótimo. Vale salientar que bases de funções ortonormais têm sido amplamente utilizadas para identificação e controle de sistemas dinâmicos (Campello et al., 2007). Na identificação de sistemas lineares, a resposta a impulso é aproximada através de uma combinação linear de funções da base escolhida. No caso mais geral, mapeamentos não-lineares precisam ser empregados (Campello and Oliveira, 2007). Os modelos assim obtidos podem ser utilizados para fins de controle utilizando, por exemplo, técnicas de controle preditivo (Oliveira et al., 2000), (Oliveira and Amaral, 2000), (Oliveira et al., 2007). Nesse contexto, costuma-se empregar funções de Laguerre e Kautz (Wahlberg and Mäkilä, 1996), (Campello et al., 2002), que permitem incorporar conhecimento a priori a respeito da dinâmica dominante do sistema (Campello et al., 2007). Tais funções são casos particulares das chamadas funções ortonormais generalizadas (Heuberger et al., 1995). Nos métodos pseudo-espectrais a serem aqui discutidos, polinômios ortonormais de Legendre ou Chebyshev são empregados em um contexto diferente. Neste caso, o modelo do sistema a ser controlado é suposto conhecido e o problema consiste na determinação de aproximações polinomiais para as trajetórias ótimas de estado e controle.

Alguns pacotes comerciais para solução de problemas de controle ótimo de larga escala oferecem a possibilidade de empregar métodos pseudo-espectrais. Exemplos incluem PROPT (Rutquist and Edvall, 2009) (www.tomdyn.com) e DIDO (www.elissar.biz), que empregam a plataforma MATLAB. Contudo, tais pacotes não permitem que o usuário visualize detalhes da implementação dos algoritmos, nem possibilitam a introdução de alterações no código, o que limita sua utilidade para fins de ensino e pesquisa. Como alternativa de código aberto, pode-se mencionar o GPOPS (Rao et al., 2008) (www.gpops.org), que usa o método pseudo-espectral de Gauss (Benson, 2004). Contudo, essa ferramenta também requer o uso de MATLAB, bem como de um pacote comercial de otimização conhecido como SNOPT (www.sbsi-sol-optimize.com). Os exemplos a serem apresentados neste artigo foram gerados empregando uma ferramenta denominada PSOPT, que foi recentemente desenvolvida por um dos autores (Becerra, 2009) (www.psopt.org). O PSOPT possui código aberto e requer apenas a disponibilidade de um compilador para C++, juntamente com bibliotecas publicamente disponíveis.

Este tutorial discutirá basicamente os métodos pseudo-espectrais de Legendre e Chebyshev. Contudo, vale salientar que os fundamentos de outros métodos são similares aos aqui apresentados. O restante do texto está organizado da seguinte forma. A Seção 2 apresenta o uso de polinômios de Legendre para aproximação de funções e para o cálculo de derivadas e integrais. A seção 3 traz uma exposição similar para o uso de polinômios de Chebyshev. A Seção 4 descreve como um problema de controle ótimo com uma única fase pode ser discretizado empregando o método pseudo-espectral de Legendre. A Seção 5 discute a discretização de problemas com múltiplas fases. A Seção 6 descreve o pacote computacional PSOPT. A Seção 7 apresenta dois exemplos ilustrativos retirados da literatura de controle ótimo computacional. Comentários finais são dados na Seção 8.

2 APROXIMAÇÃO PSEUDO-ESPECTRAL DE LEGENDRE

2.1 Aproximação de funções

Seja L_N(τ) o polinômio de Legendre de ordem N, que pode ser gerado através da seguinte expressão:

Exemplos de polinômios de Legendre são

A Figura 1 apresenta os polinômios de Legendre L_N(τ) para N = 0,1,2,5,10.

Os nós de Lagrange-Gauss-Lobatto (LGL) associados ao polinômio L_N(τ) são definidos como τ₀ = -1, τ_N = 1 e τ_k correspondendo às raízes de dL_N(τ)/dτ no intervalo [-1, 1] para k = 1,2,..., N - 1. Não há fórmulas explícitas para calcular as raízes de dL_N(τ)/dτ, mas as mesmas podem ser obtidas através de algoritmos numéricos. Para ilustração, a Figura 2 apresenta os nós LGL para N = 20.

Empregando-se polinômios de Legendre L₀(τ), L₁(τ), ..., L_N(τ), uma função f(τ): [-1,1] → pode ser aproximada por uma combinação linear da forma

com coeficientes denotados por a_k, k = 0,1,...,N.

Contudo, é possível mostrar (Hesthaven et al., 2007) que tal aproximação pode ser re-escrita de forma mais conveniente como

sendo Φ_k(τ) o polinômio interpolador de Lagrange definido por

Vale notar que Φ_k(τ_j) = 1 se k = j e Φ_k(τ_j)=0 se k ≠ j, de modo que

Em comparação com (1), a expressão (2) é mais conveniente porque os coeficientes da combinação linear são dados diretamente pelos valores da função f(τ) nos nós LGL (τ_k, k = 0,1,...,N). Além disso, a derivada de f(τ) nos nós LGL pode ser obtida de forma aproximada por uma simples multiplicação de matrizes, como se verá mais adiante.

Vale salientar que, embora as funções de base empregadas em (2) sejam polinômios de Lagrange, fala-se em uma aproximação pseudo-espectral de Legendre devido ao uso dos nós LGL. A vantagem de se empregar tais nós é ilustrada na Figura 3, considerando-se como exemplo a aproximação polinomial de ordem N para a função f(τ) = 1/(1 +τ+ 15τ² ). As colunas da esquerda e da direita apresentam os resultados obtidos com N+1 pontos equidistantes e N+1 nós LGL, respectivamente. Com o aumento de N, os erros crescem exponencialmente no primeiro caso (efeito conhecido como fenômeno de Runge), e diminuem exponencialmente no segundo.

2.2 Diferenciação numérica

A derivada de f^N(τ) nos nós LGL (τ_k) pode ser obtida diferenciando a Eq. (2). O resultado pode ser expresso na forma

em que

A expressão (5) pode ser re-escrita na seguinte forma:

A matriz D com elementos {D_ki; k, i = 0,1,...,N} é conhecida como matriz de diferenciação. Por exemplo, a matriz de diferenciação de Legendre para N = 5 é dada por

A Figura 4 apresenta o resultado da diferenciação de Legendre de f(τ) = sen(5τ²) empregando N = 5, 10 e 20. Vale notar que as escalas verticais são diferentes nas três curvas de erro. A Figura 5 apresenta o erro absoluto máximo obtido na diferenciação de Legendre de f(τ) = sen(5τ²) em função de N. Como se pode observar, o erro diminui rapidamente até ser limitado pela precisão de representação numérica do computador. Esse fenômeno é conhecido como acurácia espectral.

2.3 Integração numérica

Se h(τ) é um polinômio de grau menor ou igual a 2N-1, sua integral no intervalo τ ∈ [-1,1] pode ser calculada de forma exata através da seguinte expressão:

em que τ_k, k = 0,...,N, correspondem aos nós LGL e w_k são pesos dados por

Se L(τ) é uma função suave qualquer, então para um N apropriado, sua integral no intervalo τ ∈ [-1,1] pode ser aproximada por

Considere-se, por exemplo, o cálculo da seguinte integral:

O resultado exato com sete casas decimais é 1.9334214. Empregando N = 3, tem-se

e a integral pode ser aproximada por

obtendo-se um erro na ordem de 10^-5. Se N = 5, o resultado da aproximação é 1.9334215, correspondendo a um erro na ordem de 10^-7.

3 APROXIMAÇÃO PSEUDO-ESPECTRAL DE CHEBYSHEV

3.1 Aproximação de funções

Seja T_N(τ) o polinômio de Chebyshev de ordem N, que pode ser gerado através da seguinte expressão:

em que cos^-1 denota a função arco-cosseno. Exemplos de polinômios de Chebyshev são:

Os nós de Chebyshev-Gauss-Lobatto (CGL) no intervalo [-1,1] são definidos como τ_k = -cos(πk/N) para k = 0,1,..., N.

À semelhança do método pseudo-espectral de Legendre, o método de Chebyshev aproxima funções f(τ): [-1,1] → através de uma combinação linear da forma

sendo que agora o polinômio interpolador de Lagrange φ_k(τ) é definido por

em que

Vale notar que φ_k(τ_j) = 1 se k = j e φ_k(τ_j)=0 se k ≠ j, de modo que

assim como na aproximação de Legendre.

3.2 Diferenciação numérica

A derivada de f^N(τ) nos nós CGL pode ser obtida diferenciando a Eq. (12). O resultado pode ser expresso na forma

em que

se k ≠ i e

A expressão (16) pode ser re-escrita em forma matricial assim como em (7).

3.3 Integração numérica

Se h(τ) é um polinômio de grau menor ou igual a 2N-1, sua integral ponderada no intervalo τ ∈ [-1,1] pode ser calculada de forma exata através da seguinte expressão:

em que τ_k, k = 0,1,...,N , correspondem aos nós CGL, w_k são pesos dados por

e g(τ) é uma função de ponderação da forma

Se L(τ) é uma função suave qualquer, então para um N apropriado, sua integral ponderada no intervalo τ ∈ [-1,1] pode ser aproximada por

Desse modo, tem-se

ou seja

4 DISCRETIZAÇÃO DE PROBLEMAS DE CONTROLE ÓTIMO

Seja o seguinte problema de controle ótimo envolvendo uma única fase:

Problema

₁: Encontrar as trajetórias de controle e estado, u(t), x(t), t ∈ [t₀, t_f], bem como os parâmetros estáticos p e os instantes de tempo t₀, t_fque minimizam o seguinte índice de desempenho:

sujeito às seguintes restrições diferenciais:

restrições de trajetória:

condições de contorno:

limites de excursão para os controles e estados:

e ainda:

sendo

Empregando a transformação:

é possível reformular o problema

₁ usando uma nova variável independente τ no intervalo [-1,1], como detalhado abaixo. Por simplificidade, a seguinte notação será adotada: u(τ) ≡ u(t(τ)) , x(τ) ≡ x(t(τ)) .

Problema

₂: Encontrar as trajetórias de controle e estado, u(τ),x(τ), τ ∈ [-1, 1], bem como os parâmetros estáticos p e os instantes de tempo t₀, t_f que minimizam o seguinte índice de desempenho:

sujeito às restrições diferenciais:

restrições de trajetória:

condições de contorno:

limites de excursão para os estados e controles:

e ainda:

Para se aplicar o método de Legendre ao problema

₂, o estado x(τ), τ ∈ [-1,1] é aproximado empregando polinômios de Lagrange mediante interpolação de ordem N nos pontos LGL:

Da mesma forma, o controle u(τ), τ ∈ [-1,1] é aproximado como:

Vale notar que x^N(τ_k) = x(τ_k) e u^N(τ_k) = u(τ_k), como indicado na Equação (4).

A derivada do estado é aproximada da seguinte maneira:

sendo D_ki os elementos da matriz de diferenciação D [(N+1) × (N+1)] dados por (6).

Considere que as trajetórias de controle, estado e derivada do estado nos pontos LGL sejam armazenadas nas matrizes U^N [n_u × (N+1)], X^N [n_x × (N+1)] e ^N[n_x × (N+1)] apresentadas abaixo:

De (36), conclui-se que X^N e _N são relacionadas pela seguinte expressão:

Seja agora a seguinte matriz F^N[n_x × (N+1)] formada pelo lado direito das restrições diferenciais avaliadas nos pontos LGL:

Define-se então a matriz ζ^N [n_x×(N+1)] de defeitos diferenciais nos pontos de colocação como

Defina-se ainda a seguinte matriz H^N [n_h ×(N+1)] de restrições de trajetória avaliadas nos pontos LGL:

Por fim, o índice de desempenho para o problema

₂ pode ser aproximado da seguinte forma:

sendo os pesos w_k definidos como em (9).

É possível agora expressar

₂ como um problema de programação não-linear:

Problema

₃

sujeito a:

Neste problema, o vetor de decisão y, com dimensão n_y = n_u(N+1)+n_x(N+1)+n_p+2, é composto pelos seguintes elementos:

sendo vec(U^N), vec(X^N) vetores-coluna formados através do empilhamento vertical das colunas de U^N e X^N, respectivamente. A função de custo F : → é da forma

e a função de restrições G : → , com n_g = n_x(N+1)+n_h(N+1)+n_e+1, é dada por

Os limites para G(y) são dados por

em que 0_n denota um vetor-coluna de n zeros e [x]_n representa um vetor-coluna (nm ×1) formado pelo empilhamento de n cópias do vetor-coluna x (m ×1). Os limites para o vetor de decisão y são dados por

Exemplo 1: Sistema linear com função de custo quadrática

Considere-se o seguinte problema: Encontrar u(t), t ∈ [0, t_f] de modo a minimizar o funcional de custo

sujeito a

em que A ∈ e B ∈ , sendo Q ∈ uma matriz positivo-semidefinida e R ∈ uma matriz positivo-definida. O valor inicial dos estados é fixado em ξ₀∈ .

Este problema pode ser formulado como

₁, com a seguinte particularização para a função de custo (25), restrições diferenciais (26) e condições de contorno (28):

não havendo restrições de trajetória (27), nem limites de excursão sobre os controles (29) ou estados (30).

Para o problema

₃ resultante da discretização pseudo-espectral com N+1 nós, o vetor de decisão y corresponde a

já que t₀ e t_f são fixados e não há parâmetros p adicionais a serem otimizados. Como não foram impostos limites de excursão para os estados e controles, não há restrições da forma (46). Por outro lado, a função de restrições G(y) é dada por

com limites inferior (G_L) e superior (G_U) expressos por

Uma vez que o problema não inclui restrições de trajetória, as restrições sobre G(y) expressam somente o vínculo dinâmico entre os controles e os estados (vec(ζ^N) = ), bem como a imposição do valor inicial para os estados (x^N(-1) = ξ₀). Vale lembrar que ζ^N é a matriz de defeitos diferenciais definida em (42). Cabe ainda salientar que, devido à normalização do tempo em (32), o valor de x(t) para t = t₀ corresponde a x^N(τ) para τ = -1.

Neste problema em particular, o resultado da aplicação de métodos pseudo-espectrais pode ser comparado com uma solução analítica. Com efeito, pode-se mostrar que a solução deste problema é da forma (Costa, 2007):

sendo P(t) obtida a partir da seguinte equação diferencial de Riccati:

com P(t_f) = 0.

Suponha que as matrizes A, B, Q, R, o estado inicial x(0) e o tempo final t_f sejam dados por

A Figura 6 apresenta as soluções obtida mediante integração numérica da equação diferencial de Riccati (linhas tracejadas) e mediante aplicação do método pseudo-espectral de Legendre com 20 nós LGL (marcadores sobre as linhas tracejadas). Como se pode observar, as soluções concordam entre si.

O detalhe ampliado na Figura 6 mostra que os estados não atingem a origem ao final do horizonte de tempo especificado (t_f = 10). Uma possível modificação no problema consiste em acrescentar uma restrição terminal da forma

de modo a impor que os estados sejam levados a zero. Para isso, basta alterar a função de restrições G(y) para

com limites

lembrando que n_x = 2 e N = 19 neste exemplo.

A Figura 7 apresenta a solução gerada pelo método de Legendre após essa modificação. Para facilitar a visualização, os resultados obtidos nos nós LGL estão conectados por uma linha contínua. Como se pode ver no detalhe ampliado, a restrição terminal foi de fato satisfeita.

O problema também pode ser modificado de modo a incluir restrições no controle e/ou nos estados. A Figura 8 apresenta o resultado de se impor um limitante inferior de -0.5 para o controle u(t). Na Figura 9, tem-se o resultado de se impor ainda um limitante inferior de -0.3 para o estado x₂(t). Tais limites são incluídos introduzindo-se restrições da forma y > y_L com

sendo u_L = -0.3 e x_L = [-∞ -0.3]^T. Aqui, ∞ deve ser entendido como um valor suficientemente grande para que não haja uma restrição ativa sobre a excursão de x₁(t).

5 PROBLEMAS DE CONTROLE ÓTIMO COM MÚLTIPLAS FASES

Alguns problemas de controle ótimo podem ser formulados de forma conveniente através da definição de N_p fases. Tais fases podem ser inerentes ao problema, como no caso de uma espaçonave que ejeta estágios durante uma manobra de lançamento. Fases também podem ser introduzidas pelo analista de modo a contemplar peculiaridades da solução do problema, tais como descontinuidades nas variáveis de controle. Vale notar que as fases não precisam ser seqüenciais, como ilustrado na Figura 10. Tal disposição de fases poderia corresponder, por exemplo, ao lançamento de um veículo espacial que liberasse um estágio ao final da fase 1 e um satélite ao final da fase 2. As fases 3 e 4 envolveriam as trajetórias da espaçonave e do satélite após a separação, supondo que ambas devam ser consideradas no problema de otimização.

Um problema de controle ótimo com N_p fases pode ser formulado da seguinte forma:

Problema

₄: Encontrar trajetórias de controle e estado u⁽ⁱ⁾(t),x⁽ⁱ⁾(t), t ∈ [], bem como parâmetros estáticos p⁽ⁱ⁾, e instantes de tempo , i = 1,2,..., N_p, de modo a minimizar o seguinte índice de desempenho:

sujeito às restrições diferenciais:

às restrições de trajetória:

às condições de contorno:

às restrições de conexão entre fases:

aos limites de excursão para os estados e controles:

e ainda:

para i = 1,2,...,N_p, sendo

em que U_ψ denota o domínio da função Ψ.

Uma discussão detalhada de problemas de múltiplas fases como

₄ pode ser encontrada em (Betts, 2001).

Assim como feito para

₃, o problema ₄ pode ser discretizado de modo a formar um problema de programação não-linear. Neste caso, as variáveis de decisão para o PPNL resultante são dadas por:

em que o super-índice entre parênteses indica a fase à qual pertencem as variáveis. A função de restrição G(y) é dada por:

em que Ψ corresponde às restrições de conexão associadas ao problema, avaliadas em y.

Com base nas informações do problema, pode-se definir limites (y_L, y_U) e (G_L, G_U) para as variáveis de decisão e para a função de restrição, respectivamente, de modo a completar a definição do PPNL associado a ₄.

Exemplo 2: Sistema linear com função de custo quadrática - Problema com duas fases

A Figura 11a apresenta a trajetória percorrida pelo estado do sistema no Exemplo 1 com restrição terminal x₁(t_f) = x₂(t_f) = 0. Como se pode notar, em t = 3.8, o estado já se encontra no interior de um círculo de raio unitário centrado na origem. A fim de exemplificar o tratamento de um problema com duas fases, suponha que o estado x (t) deva permanecer no exterior do círculo entre t = 0 e t = 5.0. Neste caso, a primeira fase teria = 0 e = 5.0. A função de custo, as restrições diferenciais e o valor inicial para o estado seriam definidas como no Exemplo 1, sendo incluída ainda a seguinte restrição de trajetória:

Tal restrição pode ser colocada na forma de (71) com

A segunda fase teria = 5.0 e = 10. A função de custo, as restrições diferenciais e o valor final para o estado seriam definidas como no Exemplo 1. A restrição de conexão entre as fases seria expressa como

ou seja, o estado final para a Fase 1 seria o estado inicial para a Fase 2. Tal restrição pode ser colocada na forma de (73) com

lembrando que n_x = 2 neste exemplo.

Para o problema

₄ resultante da discretização pseudo-espectral com N+1 nós, o vetor de decisão y corresponderia a

já que , , , são fixados e não há parâmetros p⁽¹⁾ e p⁽²⁾ adicionais a serem otimizados. Neste caso, a função de restrições G(y) seria dada por

com limites inferior (G_L) e superior (G_U) expressos por

em que 1_N+1 denota um vetor-coluna com N+1 elementos iguais a 1.

A Figura 11b apresenta a solução do problema assim formulado empregando 10 nós LGL em cada fase (isto é, N = 9). Como se pode observar, a restrição de trajetória é obedecida, uma vez que o estado cruza a circunferência pontilhada em t = 5.0.

6 O PACOTE COMPUTACIONAL DE CÓDIGO ABERTO PSOPT

PSOPT é um pacote computacional de código aberto escrito em C++ que emprega métodos de discretização pseudo-espectral para resolver problemas de controle ótimo. Este pacote permite a utilização de aproximações de Legendre ou Chebyshev, podendo tratar problemas com uma ou múltiplas fases, plantas com dinâmica não-linear e atrasos de transporte, condições de contorno genéricas, restrições de trajetória não-lineares (igualdades ou desigualdades) sobre os estados e/ou variáveis de controle, restrições integrais, restrições de pontos interiores e limites de excursão sobre os controles e estados. Adicionalmente, o pacote possibilita o uso de funções de custo com termos de Lagrange e Mayer, condições iniciais e finais fixas ou livres, conexões lineares ou não-lineares entre as fases, tempos inicial e final fixos ou livres e otimização de parâmetros estáticos. A implementação faz uso de métodos de escalonamento automático, diferenciação automática utilizando a biblioteca ADOL-C (ou alternativamente, diferenciação numérica por meio de diferenças finitas esparsas) e detecção automática de esparsidade da matriz Jacobiana.

Para solução dos problemas de programação não-linear resultantes da discretização pseudo-espectral, emprega-se o pacote IPOPT. IPOPT é um pacote de código aberto escrito em C++ para otimização não-linear de larga escala, utilizando métodos de pontos interiores (Wächter and Biegler, 2006). Opcionalmente, PSOPT permite o uso de SNOPT (Gill et al., 2005), que é um pacote de otimização comercial baseado em programação seqüencial quadrática esparsa.

PSOPT não requer o uso de nenhuma plataforma comercial ou biblioteca que não seja publicamente disponível. Além disso, sua velocidade de execução tende a ser maior em comparação com rotinas implementadas em um ambiente de programação interpretado (tal como o MATLAB). Por fim, PSOPT também pode ser chamado como uma sub-rotina de aplicações escritas em C++, o que oferece maior flexibilidade de uso para o usuário.

A arquitetura do pacote PSOPT encontra-se ilustrada na Figura 12.

7 EXEMPLOS

A fim de ilustrar o potencial dos métodos descritos neste tutorial, dois exemplos serão discutidos. Em particular, será utilizado o método pseudo-espectral de Legendre. Contudo, vale salientar que a aplicação do método de Chebyshev gera soluções similares. Os resultados de cada exemplo podem ser confrontados com soluções apresentadas nas referências.

7.1 Manobra terminal de um míssil

Este exemplo ilustra a otimização da trajetória de um míssil com o intuito de atingir um dado alvo em tempo mínimo a partir de condições iniciais fixadas. Todos os parâmetros do problema foram retirados de (Subchan and Zbikowski, 2009). A Figura 13 apresenta as variáveis associadas ao modelo dinâmico do míssil considerado. Nesta figura, γ é o ângulo de trajetória, α é o ângulo de ataque, V é a velocidade, χ é a posição longitudinal, h é a altitude, D é a força aerodinâmica axial, L é a força aerodinâmica normal, T é a força de propulsão e mg é o peso do míssil.

As equações de movimento do míssil são dadas por

em que

estando os parâmetros do modelo apresentados na Tabela 1. Para fins de otimização da trajetória, os controles considerados serão a propulsão T e o ângulo de ataque α. Vale salientar que a trajetória resultante poderia servir como referência para uma malha de controle interna que efetivamente comandaria os atuadores do míssil.

As condições iniciais e finais para os estados são

e as restrições do problema são dadas por

Em vista de (89), a restrição (98) pode também ser expressa como

Vale salientar que as restrições de trajetória (99) e (100) sobre a altitude h não são suaves com respeito à posição longitudinal χ. A fim de evitar dificuldades associadas à solução do PPNL com esta não-suavidade, as restrições (99) e (100) foram aproximadas por uma só restrição da forma:

sendo uma versão suavizada da função de Heaviside, definida como:

em que > 0 é um número pequeno. Uma alternativa ao uso desta aproximação consistiria em dividir o problema em duas fases, de forma similar à empregada para o Exemplo 2 na Seção 5. Neste caso, a comutação entre as fases não ocorreria em um instante fixado a priori, mas sim quando a posição longitudinal c atingisse o valor de 7500 m.

Este problema pode ser formulado como

₁, com as seguintes particularizações:

• Estado: x = [ γ V χ h ]^T

• Controle: u = [ T α ]^T

• Função de custo:

• Restrições diferenciais: (84) - (87)

• Condições de contorno: (91) - (94)

• Restrição de excursão dos estados: (95)

• Restrição de excursão dos controles: (96) - (97)

• Restrições de trajetória: (101) - (102)

não havendo parâmetros p adicionais a serem otimizados. Para o problema ₃ resultante da discretização pseudo-espectral com N+1 nós, o vetor de decisão y corresponde a

lembrando que o tempo final t_fé livre.

O problema foi inicialmente resolvido empregando 40 nós. Em seguida, a grade foi refinada para 80 nós, utilizando-se uma interpolação da primeira solução como ponto inicial de busca para a nova solução. Como resultado, obteve-se o valor de 40.912 para a função de custo, isto é, um tempo final = 40.912 s para a manobra. A Figura 14 apresenta a trajetória ótima obtida para o míssil, com marcadores para os pontos correspondentes aos nós LGL. As Figuras 15 e 16 mostram, respectivamente, a velocidade e o ângulo de ataque em função do tempo. A Tabela 2 detalha a configuração empregada para o otimizador neste exemplo.

7.2 Lançamento de veículo espacial

Este exemplo envolve o lançamento de um veículo espacial, como descrito em (Rao et al., 2008), (Benson, 2004). O lançamento envolve quatro fases, com estágios sendo ejetados ao final das fases 1, 2 e 3. O tempo final das fases 1, 2 e 3 é fixado em 75.2 s, 150.4 s e 261.0 s, respectivamente, ao passo que o tempo final da fase 4 é livre. O problema consiste em encontrar o controle, u(t), e o tempo final da fase 4, , que minimizem a função de custo

Em outras palavras, deseja-se maximizar a massa do veículo (isto é, a quantidade de combustível remanescente) ao final da fase 4. A dinâmica do problema é descrita pelas seguintes equações:

em que r(t) = [ x(t) y(t) z(t) ]^T é a posição do veículo com relação ao centro da Terra, v = [ v_x(t) v_y(t) v_z(t) ]^T é a velocidade, µ > 0 é uma constante associada à força de atração gravitacional, T é a força de propulsão, m é a massa, g₀ é a aceleração da gravidade ao nível do mar, I_sp é o impulso específico do motor, u = [ u_x u_y u_z]^T é a direção da propulsão, e D = [ D_x D_y D_z]^T é o arrasto, que é dado por

sendo C_D o coeficiente de arrasto, A_ref a área de referência, ρ a densidade do ar, e v_rel a velocidade com relação à atmosfera, dada por

em que ω é a velocidade angular da Terra com respeito a um referencial inercial. A densidade do ar é dada por

sendo ρ₀ a densidade ao nível do mar, h = ║r║-R_e a altitude, R_e o raio equatorial da Terra e h₀ um parâmetro de escala. Os valores para essas constantes podem ser encontrados em (Rao et al., 2008) e (Benson, 2004).

O lançamento é iniciado no instante t₀ ao nível do mar com o veículo em repouso com relação à Terra. As condições iniciais são

As restrições terminais estão relacionadas com a órbita desejada, que é definida em parâmetros orbitais (Bate et al., 1971) como

Tendo em vista que u é um vetor unitário que indica a direção da força de propulsão, deve-se ainda impor a seguinte restrição de trajetória para o problema:

Fisicamente, sabe-se que a posição e a velocidade não podem apresentar descontinuidades na passagem de uma fase para a seguinte. Desse modo, deve-se impor as seguintes restrições de conexão entre fases:

em que (i) representa o índice da fase. Além disso, é necessário acrescentar uma restrição de conexão que modele a discontinuidade de massa causada pela ejeção de estágios ao final das fases 1, 2 e 3:

Este problema pode ser formulado como

₄, com as seguintes particularizações:

• Estado: x⁽ⁱ⁾ = [ r⁽ⁱ⁾v⁽ⁱ⁾m⁽ⁱ⁾ ]^T, i = 1,...,4

• Controle: u⁽ⁱ⁾ = , i = 1,...,4

• Função de custo:

• Restrições diferenciais: (105) - (107) em cada fase

• Condições de contorno:

r⁽¹⁾(), v⁽¹⁾() m⁽¹⁾() dados por (111) - (113) e r⁽⁴⁾(), v⁽⁴⁾() fixados de acordo com os parâmetros orbitais (114) - (118).

• Restrição de trajetória: (119) em cada fase

• Restrições de conexão entre fases: (120) - (122)

não havendo parâmetros p adicionais a serem otimizados, nem restrições sobre a excursão dos estados. Vale salientar que a restrição de trajetória (119) já estabelece indiretamente uma restrição sobre a excursão de cada componente do controle: |u_j(t)| < 1, j = 1,2,3, t ∈ [t₀, ].

Para o problema resultante da discretização pseudo-espectral com N+1 nós, o vetor de decisão y corresponde a

lembrando que o tempo final é livre.

O problema foi inicialmente resolvido empregando 5 nós por fase (20 nós no total). Em seguida, a grade foi refinada para 15 nós por fase (60 nós no total), utilizando-se uma interpolação da primeira solução como ponto inicial de busca para a nova solução.

Como resultado, obteve-se o valor de -7529.7 para a função de custo, isto é uma massa final de 7529.7 kg ao final do lançamento. As Figuras 17, 18 e 19 apresentam a altitude, velocidade e elementos do vetor de controle obtidos como solução do problema. A Tabela 3 detalha a configuração empregada para o otimizador neste exemplo.

8 CONCLUSÕES

Este artigo apresentou um tutorial sobre o uso de métodos de aproximação pseudo-espectral para solução computacional de problemas de controle ótimo. Fundamentos da teoria de aproximação empregando os métodos pseudo-espectrais de Legendre e Chebyshev foram discutidos, mostrando-se como problemas de uma ou mais fases podem ser discretizados na forma de problemas de programação não-linear. Nesse contexto, apresentou-se ainda um pacote de código aberto denominado PSOPT. Para fins de ilustração, tal pacote foi aplicado a dois exemplos retirados da literatura de controle ótimo computacional.

AGRADECIMENTOS

Os autores agradecem o apoio da FAPESP (Processo 2006/58850-6) e CNPq.

REFERÊNCIAS

Bate, R. R., Mueller, D. D. White, J. E. 1971 . Fundamentals of Astrodynamics, Dover, New York.

Becerra, V. 2009 . PSOPT: Optimal Control Solver User Manual.

Bedrossian, N. S., Bhatt, S., Kang, W. Ross, I. M. 2009 . Zero-propellant maneuver guidance, IEEE Control Systems Magazine 29(5): 53-73.

Benson, D. A. 2004 . A Gauss Pseudospectral Transcription for Optimal Control, PhD thesis, MIT, Department of Aeronautics and Astronautics, Cambridge, Mass.

Betts, J. T. 2001 . Practical Methods for Optimal Control Using Nonlinear Programming, SIAM, Philadelphia.

Campello, R. J. G. B., Amaral, W. C. Favier, G. 2002 . Desenvolvimento ótimo de modelos de Volterra em funções ortonormais de Laguerre, Anais do XIV Congresso Brasileiro de Automática, Natal, RN, pp. 128-133.

Campello, R. J. G. B. Oliveira, G. H. C. 2007 . Modelos não-lineares, in L. A. Aguirre, A. P. A. Silva, M. F. M. Campos W. C. Amaral (eds), Enciclopédia de Automática: Controle e Automação, Vol. 3, Editora Blucher, São Paulo, pp. 104-122.

Campello, R. J. G. B., Oliveira, G. H. C. Amaral, W. C. 2007 . Identificação e controle de processos via desenvolvimentos em séries ortonormais. Parte A: Identificação, Controle & Automação 18(3): 301-321.

Canuto, C. G., Hussaini, M. Y., Quarteroni, A. Zang, T. A. 1988 . Spectral Methods in Fluid Dynamics, Springer-Verlag, Berlin.

Canuto, C., Hussaini, M., Quarteroni, A. Zang, T. 2006 . Spectral Methods: Fundamentals in Single Domains, Springer-Verlag, Berlin.

Canuto, C., Hussaini, M., Quarteroni, A. Zang, T. 2007 . Spectral Methods: Evolution to Complex Geometries and Applications to Fluid Dynamics, Springer, Heidelberg.

Costa, O. L. V. 2007 . Projeto LQG, in L. A. Aguirre, A. H. Bruciapaglia, P. E. Miyagi R. H. C. Takahashi (eds), Enciclopédia de Automática: Controle e Automação, Vol. 1, Editora Blucher, São Paulo, pp. 86-110.

Elnagar, G., Kazemi, M. A. Razzaghi, M. 1995 . The pseudospectral Legendre method for discretizing optimal control problems, IEEE Transactions on Automatic Control 40(10): 1793-1796.

Fahroo, F. Ross, I. 2001 . Costate estimation by a Legendre pseudospectral method, Journal of Guidance, Control, and Dynamics 24(2): 270-277.

Fahroo, F. Ross, I. 2002 . Direct trajectory optimization by a Chebyshev pseudospectral method, Journal of Guidance, Control, and Dynamics 25(1): 160-166.

Fahroo, F. Ross, I. M. 2008 . Pseudospectral methods for infinite-horizon optimal control problems, Journal of Guidance, Control, and Dynamics 31(4): 927-936.

Gill, P. E., Murray, W. Saunders, M. A. 2005 . SNOPT: An SQP algorithm for Large-Scale Constrained Optimization, SIAM Review 47(1): 99-131.

Gong, Q., Fahroo, F. Ross, I. M. 2008 . Spectral algorithms for pseudospectral methods in optimal control, Journal of Guidance, Control, and Dynamics 31(3): 460 - 471.

Hesthaven, J., Gottlieb, S. Gottlieb, D. 2007 . Spectral Methods for Time-Dependent Problems, Cambridge University Press, Cambridge.

Heuberger, P. S. C., Van den Hof, P. M. J. Bosgra, O. H. 1995 . A generalized orthonormal basis for linear dynamical systems, IEEE Transactions on Automatic Control 40(3): 451-465.

Kang, W. Bedrossian, N. 2007 . Pseudospectral optimal control theory makes debut flight, saves nasa $1m in under three hours, SIAM News 40(7).

Kang, W., Gong, Q., Ross, I. M. Fahroo, F. 2007 . On the convergence of nonlinear optimal control using pseudospectral methods for feedback linearizable systems, International Journal of Robust and Nonlinear Control 17(14): 1251-1277.

Oliveira, G. H. C. Amaral, W. C. 2000 . Identificação e controle preditivo de processos não-lineares utilizando séries de Volterra e bases de funções ortonormais, Anais do XIII Congresso Brasileiro de Automática, Florianópolis, SC, pp. 2174-2179.

Oliveira, G. H. C., Amaral, W. C., Favier, G. Dumont, G. A. 2000 . Constrained robust predictive controller for uncertain processes modeled by orthonormal series functions, Automatica 36(4): 563-571.

Oliveira, G. H. C., Campello, R. J. G. B. Amaral, W. C. 2007 . Identificação e controle de processos via desenvolvimentos em séries ortonormais. Parte B: Controle preditivo, Controle & Automação 18(3): 322-336.

Rao, A., Benson, D., Huntington, G. Francolin, C. 2008 . User's manual for GPOPS version 1.3: A Matlab package for dynamic optimization using the Gauss pseudospectral method.

Ross, I. Fahroo, F. 2004 . Pseudospectral knotting methods for solving optimal control problems, Journal of Guidance, Control, and Dynamics 27(3): 397-405.

Rutquist, P. E. Edvall, M. M. 2009 . PROPT Matlab Optimal Control Software, TOMLAB Optimization.

Subchan, S. Zbikowski, R. 2009 . Computational Optimal Control: Tools and Practice, John Wiley, New York.

Trefethen, L. N. 2000 . Spectral Methods in MATLAB, SIAM, Philadelphia.

Vlassenbroeck, J. 1988 . A Chebyshev polynomial method for optimal control with state constraints, Automatica 24(4): 499-506.

Vlassenbroeck, J. Dooren, R. V. 1988 . A Chebyshev technique for solving nonlinear optimal control problems, IEEE Transactions on Automatic Control 33(4): 333-340.

Wächter, A. Biegler, L. T. 2006 . On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming, Mathematical Programming 106(1): 25-57.

Wahlberg, B. Mäkilä, P. M. 1996 . Approximation of stable linear dynamical systems using Laguerre and Kautz functions, Automatica 32(5): 693-708.

Artigo submetido em 22/10/2009 (Id.: 01069)

Revisado em 30/11/2009, 02/02/2010

Aceito sob recomendação do Editor Associado Prof. Luis Fernando Alves Pereira

Datas de Publicação

• Publicação nesta coleção
  29 Jun 2010
• Data do Fascículo
  Jun 2010

Histórico

• Revisado
  30 Nov 2009
• Recebido
  22 Out 2009
• Aceito
  02 Fev 2010