Dando continuidade ao método pseudo-espectral (discutido em detalhes aqui e aqui). A formulação do problema de controle ótimo anteriormente proposto, pode ser reescrita como um problema de programação não-linear, da seguinte forma:
sujeito a:
Neste problema, o vetor de decisão y, com dimensão 
, composto pelos seguintes elementos:
, composto pelos seguintes elementos:
sendo 
e 
vetor-coluna formados através do empilhamento vertical das colunas de 
e 
, respectivamente. A função de custo F pode ser reescrita por meio de:
e vetor-coluna formados através do empilhamento vertical das colunas de e , respectivamente. A função de custo F pode ser reescrita por meio de:
e a função de restrições G, com 
, dada por:
, dada por:
Os limites inferiores e superiores (
e 
, respectivamente) para 
são dados por:
e , respectivamente) para são dados por:
em que 
denota um vetor-coluna de n zeros e 
representa um vetor-coluna 
formado pelo empilhamento de n cópidas do vetor-coluna x 
. Os limites inferiores e superiores (
e 
, respectivamente) para o vetor de decisão y são dados por:
denota um vetor-coluna de n zeros e representa um vetor-coluna formado pelo empilhamento de n cópidas do vetor-coluna x . Os limites inferiores e superiores ( e , respectivamente) para o vetor de decisão y são dados por:
Pode ter ficado confuso, mas nas secões a seguir a gente fará alguns exemplos numéricos com exemplos comumente apresentados no formato de um problema de controle ótimo e vamos solucioná-los.
Exemplo 1
Deseja-se encontrar um controle 
, 
de modo a minimizar o funcional de custo:
, de modo a minimizar o funcional de custo:
sujeito a:
em que 
e 
, sendo 
uma matriz positivo-semidefinida e 
uma matriz positivo-definida. O valor inicial dos estados é 
. Reescrevendo o problema como um problema de programação não linear resultante da discretização pseudo-espectral com 
nós, o vetor de decisão corresponde a:
e , sendo uma matriz positivo-semidefinida e uma matriz positivo-definida. O valor inicial dos estados é . Reescrevendo o problema como um problema de programação não linear resultante da discretização pseudo-espectral com nós, o vetor de decisão corresponde a:
Isso se dá porque não há parâmetro adicional, o tempo inicial 
e final 
são fixados (o valor de será especificado a seguir). Não há limites de excursão para os estados e controles (são variáveis livres). E por outro lado, a função de restrições é dada por:
e final são fixados (o valor de será especificado a seguir). Não há limites de excursão para os estados e controles (são variáveis livres). E por outro lado, a função de restrições é dada por:
com os limites inferiores e superiores dados por:
Isso decorre pelo fato de que não há restrições de trajetória (
) para os estados e controle, e há restrição inicial para o estado no instante inicial (
). Suponha que as matrizes do sistema dado (
) , o estado inicial 
e o tempo final sejam dados por:
) para os estados e controle, e há restrição inicial para o estado no instante inicial (). Suponha que as matrizes do sistema dado () , o estado inicial e o tempo final sejam dados por:
Vamos resolver por meio de código em Matlab. Para os dados propostos, tem-se que:
clear all; close all; clc;
% Matrizes do problema
A = [0 1;
0 0];
B = [0;
1];
Q = eye(2);
R = 10;
t0 = 0;
tf = 10;
Definição da quantidade de nós N para a discretização pseudo-espectral:
N = 39;
Parâmetros derivados a serem utilizados na formulação da otimização pseudo-espectral (já foram definidos anteriormente em aqui e aqui):
% Parâmetros derivados a partir de 'N'
nos = tau(N);
D = matriz_D(N, nos);
wk = pesos_wk(N, nos);
Para a execução do restante do código, é fundamental que sejam definidas as seguintes funções de custo:
function J = funcao_custo(y, Q, R, N, peso, t0, tf)
%%% Alocação inicial de memória.
persistent u x1 x2 X
%%% Criação das dimensões dos vetores (uma única vez)
if isempty(u) || length(u) ~= (N+1)
u = zeros(1, N+1);
x1 = zeros(1, N+1);
x2 = zeros(1, N+1);
X = zeros(2, N+1);
end
%%% Separa as variáveis para uso (mantendo área alocada)
u(:) = y(1, 0*(N+1) + 1 : 1*(N+1)); % Controle u
x1(:) = y(1, 1*(N+1) + 1 : 2*(N+1)); % Estado x1
x2(:) = y(1, 2*(N+1) + 1 : 3*(N+1)); % Estado x2
% Vetor X atualizado
X(1,:) = x1;
X(2,:) = x2;
%%% Função de custo
integral = 0;
for k = 1:(N+1)
Lk = X(:,k)' * Q * X(:,k) + u(k)' * R * u(k);
integral = integral + Lk * peso(k);
end
J = ((tf-t0)/2) * integral;
end
e a função de constraints (restrições):
function [c, ceq] = constraints(y, N, D, A, B)
%%% Para alocação de memória
persistent u x1 x2 X FN chiN
if isempty(u) || length(u) ~= (N+1)
u = zeros(1, N+1);
x1 = zeros(1, N+1);
x2 = zeros(1, N+1);
X = zeros(2, N+1);
FN = zeros(2, N+1);
chiN = zeros(2, N+1);
end
%%% Separa as variáveis para uso (mantendo área alocada)
u(:) = y(1, 0*(N+1) + 1 : 1*(N+1)); % Controle u
x1(:) = y(1, 1*(N+1) + 1 : 2*(N+1)); % Estado x1
x2(:) = y(1, 2*(N+1) + 1 : 3*(N+1)); % Estado x2
X = [x1; x2];
% Cálculo do FN
FN = matriz_FN(y, N, A, B);
% Determinação do vetor chi^(N)
chiN = X * D' - FN;
% Restrições de desigualdade
rest_0 = [x1(1); x2(1)] - [1; 1];
c = [ chiN(:);
-chiN(:);
rest_0(:);
-rest_0(:);];
% Restrições de igualdade
ceq = [];
end
Dentro dessas funções de constraints, é chamada a função matriz_FN, que é dada por:
function FN = matriz_FN(y, N, A, B)
%%% Alocação de memória
persistent u x1 x2 X X_dot
if isempty(u) || length(u) ~= (N+1)
u = zeros(1, N+1);
x1 = zeros(1, N+1);
x2 = zeros(1, N+1);
X = zeros(2, N+1);
X_dot = zeros(2, N+1);
end
%%% Separa as variáveis para uso (mantendo área alocada)
u(:) = y(1, 0*(N+1) + 1 : 1*(N+1)); % Controle u
x1(:) = y(1, 1*(N+1) + 1 : 2*(N+1)); % Estado x1
x2(:) = y(1, 2*(N+1) + 1 : 3*(N+1)); % Estado x2
X(1,:) = x1;
X(2,:) = x2;
% Tempo inicia e final
t0 = 0;
tf = 10;
%%% Inicialização dos parâmetros a se resolverem.
X_dot = zeros(2, N+1);
for k = 1:(N+1)
X_dot(:,k) = A * X(:,k) + B * u(k);
end
FN = ((tf-t0)/2) * X_dot;
end
De posse dessas funções, a gente consegue configurar o restante do código para emprego da função fmincon, do Matlab:
% Instância da função de custo
Fy = @(y)funcao_custo(y, Q, R, N, wk, t0, tf);
% Instância de restrições
nonlcon = @(y)constraints(y, N, D, A, B);
% Demais parâmetros para a função fmincon
Aa = []; % Não há restrição do tipo 'Aa.x <= bb';
bb = []; % Não há restrição do tipo 'Aa.x <= bb';
Aeq = []; % Não há restrição do tipo 'A.x = b';
beq = []; % Não há restrição do tipo 'A.x = b';
lb = []; % Não há, inicialmente, restrições do tipo 'lb <= x'; e
lu = []; % Não há, inicialmente, restrições do tipo 'x <= ub'.
% Opções do otimizador
options = optimoptions('fmincon', ...
'Display', 'iter', ...
'Algorithm','SQP', ...
'MaxFunEvals', 5000, ...
'TolX', 1e-6);
% Chute inicial para as variáveis do otimizador.
u = zeros(1, N+1);
x1 = zeros(1, N+1);
x2 = zeros(1, N+1);
y0 = [u, x1, x2];
% Roda o otimizador.
solution = fmincon(Fy, y0, Aa, bb, Aeq, beq, lb, lu, nonlcon, options);
Depois da resolução do problema, a gente consegue plotar o resultado por meio dos seguintes códigos:
% Extrai os valores do sistema
t_original = ((tf - t0) / 2) * nos + ((tf + t0) / 2);
u(:) = solution(1, 0*(N+1) + 1 : 1*(N+1)); % Controle u
x1(:) = solution(1, 1*(N+1) + 1 : 2*(N+1)); % Estado x1
x2(:) = solution(1, 2*(N+1) + 1 : 3*(N+1)); % Estado x2
% Plota os resultados das variáveis de estados.
subplot(2,1,1);
plot(t_original, x1); hold on;
plot(t_original, x2, 'r');
ylim([-0.5 2]);
grid;
xlabel("t");
legend('x_1(t)', 'x_2(t)');
% Variáveis de controle.
subplot(2,1,2);
plot(t_original, u);
grid;
ylim([-1.5 0.5]);
xlabel("t");
ylabel("u(t)");
