Métodos pseudo-espectral (parte 2)

Dando continuação ao registro comentado, essa seção procurará abordar a integrala numérica utilizando os métodos de Legendre. No post passado eu utilizei apresentei o método de Legendre e comentei sobre a derivação numérica. Para maiores detalhes, só ir neste link.

Seja uma função um polinômio de grau menor ou igual a . Sua integral no intervalo pode ser calculada de forma exata por meio da seguinte expressão algébrica:

em que , para , correspondem aos nós LGL e são os pesos, dados por:

Se é uma função suave qualquer, então para um N apropriado, sua integral no intervalo pode ser aproximada por:

Ou seja, de posse de uma função a integral pode ser representada por meio de uma soma "relativamente simples". De modo a poder calcular os pesos , eu criei a seguinte função em Matlab:

function wk = pesos_wk(N, nos_tau)

% ====================================================

% (*) Encontrar os parâmetros de pesos 'w(k)', da eq 9

% do paper. Para isso, utilizarei variáveis simbó-

% licas.

%

% Entradas:

%

% 1. N ...........: ordem da aproximação;

% 2. tau .........: nós para o cálculo de L(tau_k)

%

% Saída:

%

% 1. out .........: todos os w_k calculados.

% ====================================================

% Variável simbólica

syms x

 

for k = 1:(N+1)

% Cálculo propriamente dito

num = diff((x^2 - 1)^N, N);

den = (2^N) * fatorial(N);

% Polinômio de Legendre.

Ln = num/den;

% Cálculo de LN(tau_k)

p = double(subs(Ln, nos_tau(k)));

% Saída

wk(k) = 2/(N*(N+1)*p^2);

end

 

end

A título de exemplo, considere o interesse em calcular a seguinte integral:

Ela pode ser facilmente calculada por meio das ferramentas já disponívels no Matlab, conforme código abaixo:

clear all; close all; clc;

 

syms t

 

% Calcula a integral de modo simbólico

F = int(exp(t)*cos(t), t, -1, 1);

 

% Converte para numérico

val = double(F);

vpa(val, 8)

ans = 

1.9334215

 

Ou pode ser calculada utilizando o método pseudo-espectral de Legendre:

clear all; close all; clc;

 

% Ordem do sistema.

N = 8;

 

% Nós tau na ordem N.

nos_tau = tau(N);

 

% Pesos w_k.

pesos = pesos_wk(N, nos_tau);

 

soma = 0;

for i = 1:max(size(nos_tau))

L_tau = exp(nos_tau(i))*cos(nos_tau(i));

soma = soma + L_tau * pesos(i);

end

vpa(soma, 8)

ans = 

1.9334215

Tá, mas beleza. Para quê todos esses passos foram deduzidos, programados e explicados? A ideia é utilizar o método pseudo-espectral para resolução de problemas de controle ótimo. Desse modo, a gente irá abordar a característica geral do problema e como transformar ele num formato em que se empregue as ferramentas pseudo-espectral.

Seja um problema de controle ótimo em que se deseja encontrar as trajetórias de controle e estado bem como os parâmetros estáticos p e os instantes de tempo que minimizam o seguinte índice de desempenho:

sujeito às restrições diferenciais para cara um dos estados estados :

restrições de trajetória, com limite inferior () e limite superior ():

condições de contorno, com limite inferior () e limite superior ():

limites superiores (Upper) e inferiores (Lower) de excursão para os controle () e estados ():

e ainda, restrições para o parâmetros p e os parâmetros temporais e .

Esse é tipicamente a formulação de um problema de controle ótimo comumente encontrado na literatura. A gente consegue, por meio de algumas transformações, converter o problema acima apresentado em algo que possa ser utilizado utilizando o método pseudo-espectral

Inicialmente, será realizada a seguinte transformação algébrica,

onde agora e possível utilizar uma variável independente . Perceba agora que para o novo problema, os extremos temporais serão -1 e 1. Já está enxergando onde a gente quer chegar, né?

Reformulando o problema de controle ótimo tendo por base a nova variável independente τ, tenose que deseja-se agora encontrar as trajetórias de controle e estado bem como os parâmetros estáticos p e os instantes de tempo que minimizam o seguinte índice de desempenho:

sujeito às restrições diferenciais para cara um dos estados estados :

restrições de trajetória, com limite inferior () e limite superior ():

condições de contorno, com limite inferior () e limite superior ():

limites superiores (Upper) e inferiores (Lower) de excursão para os controle () e estados ():

e ainda:

Para aplicar o método de Legendre, o estado , e o controle, é aproximado empregando polinômios de Lagrange mediante interpolação de ordem N nos pontos LGL, conforme:

A derivada do estado, conforme já abordado no passo 1 deste registro, tem-se que:

sendo os elementos os elementos da matriz de diferenciação obtido no passo 1 (veja lá, meu brother).

No paper, ele propões a organização das trajetórias de controle, estado e derivadas do estado nos pontos LGL seguindo a seguinte notação notação , e :

As matrizes e estão relacionadas por meio de:

Seja, ainda, a matriz formada pelo lado direito das restrições diferenciais avaliadas nos pontos LGL:

e a matriz com as restrições de trajetória avaliadas nos pontos LGL:

Nesse ínterim, defini-se a matriz de defeitos diferenciais nos pontos de colocação por meio de:

e com isso redefin-se o índice de desempenho inicial por meio de:

Agora sim é possível expressar toda a formulação apresentada anteriormente como um problema de programação linear.