Métodos pseudo-espectral (parte 1)

Recentemente um colega trouxe com um paper em que abordava o método-pseudo espectral para cálculo de trajetória de um míssil com restrições terminais. O paper é o "UM TUTORIAL SOBRE MÉTODOS PSEUDO-ESPECTRAIS PARA CONTROLE ÓTIMO COMPUTACIONAL", do Becerra e do Kawakami. Ele está disponível na internet e é relativamente fácil de fazer o download.

Eu havia estudado controle ótimo como matéria no mestrado mas não havia ouvido falar de métodos pseudo-espectral. No fim das contas, se tornou um desafio interessante, pois eu voltei a programar alguns exercícios, utilizando o método pseudo-espectral bem como a forma analítica comumente empregada nos problemas de controle ótimo. Neste espaço, eu abordarei a conceituação e trarei códigos em Matlab para a gente ir testando as coisas aos poucos.

O paper aborda dois métodos pseudo-espectral: o de Legendre e o de Chebyshev. Eu irei me atentar, inicialmente, para o de Legendre e a questão da derivada do sistema. Mas a filosofia em se utilizando do método de Chebyshev similar e foi intencionalmente omitida.

Aproximação pseudo-espectral de Legendre

O paper aborda a criação dos polinômios de Legendre de ordem N. Ele é dado por meio da seguinte equação algébrica:

Nesse universo, o que nos interessa mesmo são os nós de Lagrande-Gauss-Lobato (LGL) ao polinômio , definida no intervalor . O nó inicial é definido por , o final por e os demais , com , são obtidos por meio das raízes de .

Tá, mas por que esses polinômios são importantes? Esses polinômios podem aproximar uma função , por meio de?

em que , para Essa mesma função pode ser reescrita por meio da seguinte equação algébrica:

em que é o polinômio interpolador de Lagrange, definido por:

Vale notar que , se , e se .

Uma outra coisa que chama a atenção é diferenciação numérica da função nos nós LGL (), que resulta em:

em que:

Agrupando todas as derivadas da função nos nós LGL (), pode-se expressar todos os resultados da forma:

Eu vou rodar um exemplo, aos moldes do que está presente no paper, com a codificação realizada no Matlab para que se veja o resultado obtido. Para tal, eu criei a função fatorial abaixo:

function x_out = fatorial(x)

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

% Calcula o fatorial de um número inteiro

% não negativo

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

if ( (x == 0) || (x == 1) )

x_out = 1;

else

x_out = x * fatorial(x - 1);

end

 

end

a função polinomio_Legendre abaixo:

function out = polinomio_Legendre(N, tau_k)

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

% Esta função calcula a polinômio de Legendre no nó

% tau_k.

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

syms x

 

% Cálculo propriamente dito

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

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

 

% Polinômio de Legendre.

Ln = num/den;

 

% Substitui e envia para saída.

out = double(subs(Ln, tau_k));

 

end

a função matriz_D:

function out = matriz_D(N, nos)

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

% Obtenção da matriz Dki, doravante de-

% nominada apenas por 'D'.

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

 

% Alocação de memória.

out = zeros(N+1, N+1);

 

% Determina individualmente os elementos da matriz D.

for k = 1:(N+1)

for i = 1:(N+1)

if ( k ~= i )

num = polinomio_Legendre(N, nos(k));

den = (nos(k) - nos(i)) * polinomio_Legendre(N, nos(i));

out(k, i) = num/den;

elseif ( (k == 1) && (i == 1) )

out(k, i) = -N*(N+1)/4;

elseif ( (k == N+1) && (i == N+1) )

out(k, i) = N*(N+1)/4;

else

out(k, i) = 0;

end

end

end

end

e a função tau, para a determinação dos nós do sistema:

function out = tau(N)

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

% Determinação dos nós dos sistema de ordem

% N

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

syms x

 

% Cálculo propriamente dito

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

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

 

% Polinômio de Legendre.

Ln = num/den;

 

% Diferenciar o polinômio de Legendre.

p = diff(Ln);

 

% Encontra as raízes.

intermed = double(solve( p==0 ));

 

% Ordenação.

out = sort([-1 intermed' 1]);

 

end

para colocar tudo isso para rodar e, consequentemente fornecer o resultado apresentado para o exemplo com , tem-se o seguinte trecho de código:

clear all; close all; clc;

 

% Ordem da aproximação.

N = 5;

 

% nós no intervalo [-1, 1] .

tau_k = tau(N);

 

% matriz D.

D = matriz_D(N, tau_k);

disp(D);

-7.5000 10.1414 -4.0362 2.2447 -1.3499 0.5000 -1.7864 0 2.5234 -1.1528 0.6535 -0.2378 0.4850 -1.7213 0 1.7530 -0.7864 0.2697 -0.2697 0.7864 -1.7530 0 1.7213 -0.4850 0.2378 -0.6535 1.1528 -2.5234 0 1.7864 -0.5000 1.3499 -2.2447 4.0362 -10.1414 7.5000

No caso da aproximação funcional do exemplo do paper , a gente irá fazer o comparativo com a sua derivada. Dada a ordem da aproximação, ou seja, o valor do N, a gente encontra os nós do sistema, a matriz D e em seguida a função derivada .

% Determinação das entradas do sistema.

N = [5, 10, 20];

% N = [5];

 

% Para o plot da função de seno (para o plot bonitinho)

t = linspace(-1, 1, 100);

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

fseno(i) = sin(5*t(i)^2);

dfseno(i) = 10*t(i)*cos(5*t(i)^2);

end

 

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

% Determina os nós e a matriz D.

tau_k = tau(N(i));

D = matriz_D(N(i), tau_k);

% Função f e df nos nós.

for k = 1:max(size(tau_k))

f(k) = sin(5*tau_k(k)^2);

end

 

% Obtenção dos pontos de derivada do sistema.

dfN = D * f';

 

if (i > 1)

figure;

end

 

% Plota o resultado de f

subplot(1, 2, 1);

plot(t, fseno);

hold on;

plot(tau_k, f, '*');

grid;

legend("sin(5\tau^2)", "sin(5\tau_k^2)");

ylabel("f(\tau)");

xlabel("\tau");

% Plota o resultado de df

subplot(1, 2, 2);

plot(t, dfseno);

hold on;

plot(tau_k, dfN, '*');

grid;

legend("10\taucos(5\tau^2)", "1\tau_k.cos(5\tau_k^2)");

ylabel("df/d\tau");

xlabel("\tau");

 

 

end