1. Programa Principal (programa_principal.m)

%% Estrutura principal para otimização multiobjetivo com GWO e NSGA-III
clear; clc; close all;

%% 1. Preparação dos dados
fprintf('1. Carregando e processando os dados...\n');
load('conjunto_dados_multisaida.mat'); % Carrega dados multivariáveis
% Dados de entrada e saída: entradas_treinamento, saidas_treinamento, entradas_teste, saidas_teste

% Normalização dos dados
[entradas_treino_norm, config_entrada] = mapminmax(entradas_treinamento', 0, 1);
[saidas_treino_norm, config_saida] = mapminmax(saidas_treinamento', 0, 1);
entradas_teste_norm = mapminmax('apply', entradas_teste', config_entrada);
saidas_teste_norm = mapminmax('apply', saidas_teste', config_saida);

entradas_treino_norm = entradas_treino_norm';
saidas_treino_norm = saidas_treino_norm';
entradas_teste_norm = entradas_teste_norm';

%% 2. Configuração do problema de otimização multiobjetivo
fprintf('2. Definindo o problema de otimização...\n');
num_variaveis = 6; % Número de variáveis a otimizar
% Definição das variáveis: [taxa_aprendizado, dimensao_oculta, num_cabecas, num_camadas, taxa_dropout, tamanho_lote]
limites_inferiores = [1e-4, 32, 2, 1, 0.1, 16];
limites_superiores = [1e-2, 256, 8, 4, 0.5, 128];

% Função objetivo múltipla
funcao_objetivo = @(x) funcao_multiobjetivo(x, entradas_treino_norm, saidas_treino_norm, ...
                                            entradas_teste_norm, saidas_teste_norm, config_entrada, config_saida);

%% 3. Execução da otimização com NSGA-III melhorado
fprintf('3. Iniciando a otimização multiobjetivo com NSGA-III melhorado...\n');
configuracoes = struct('TamanhoPopulacao', 50, ...
                       'MaximoGeracoes', 30, ...
                       'TaxaCrossover', 0.8, ...
                       'TaxaMutacao', 0.1, ...
                       'ExibirProgresso', 'iter');

[fronteira_pareto, conjunto_pareto, informacoes_otimizacao] = ...
    nsga3_melhorado(funcao_objetivo, num_variaveis, limites_inferiores, limites_superiores, configuracoes);

%% 4. Análise dos resultados
fprintf('4. Analisando os resultados e gerando visualizações...\n');

% Visualização da fronteira de Pareto
figura = figure('Position', [100, 100, 800, 600]);
scatter(fronteira_pareto(:,1), fronteira_pareto(:,2), 50, 'filled', ...
    'MarkerFaceColor', [0.3 0.5 0.7], 'MarkerEdgeColor', 'k');
xlabel('Objetivo 1: Erro de Previsão (RMSE)');
ylabel('Objetivo 2: Complexidade do Modelo (Número de Parâmetros)');
title('Fronteira de Pareto Obtida');
grid on;

% Seleção da solução com menor erro de previsão
[~, indice] = min(fronteira_pareto(:,1));
solucao_otima = conjunto_pareto(indice, :);

%% 5. Treinamento do modelo final com a solução ótima
fprintf('5. Treinando o modelo final com os parâmetros otimizados...\n');
parametros_finais = decodificar_parametros(solucao_otima, limites_inferiores, limites_superiores);
[modelo_final, previsoes, desempenho] = ...
    treinar_transformer_final(parametros_finais, entradas_treino_norm, saidas_treino_norm, ...
                              entradas_teste_norm, saidas_teste_norm, config_entrada, config_saida);

%% 6. Avaliação do desempenho do modelo
fprintf('6. Avaliando o desempenho do modelo...\n');
avaliar_desempenho_modelo(previsoes, saidas_teste_norm, config_saida, saidas_teste);

%% 7. Salvando os resultados
save('resultados_otimizacao.mat', 'fronteira_pareto', 'conjunto_pareto', ...
     'solucao_otima', 'modelo_final', 'desempenho', 'informacoes_otimizacao');

2. Algoritmo NSGA-III Melhorado (nsga3_melhorado.m)

function [fronteira_pareto, conjunto_pareto, informacoes] = nsga3_melhorado(funcao_objetivo, ...
                                                                            num_variaveis, limites_inf, limites_sup, configuracoes)
    % Algoritmo NSGA-III melhorado com integração do GWO
    
    % Extrair parâmetros das configurações
    tamanho_populacao = configuracoes.TamanhoPopulacao;
    maximo_geracoes = configuracoes.MaximoGeracoes;
    prob_crossover = configuracoes.TaxaCrossover;
    taxa_mutacao = configuracoes.TaxaMutacao;
    
    % Inicializar a população
    populacao = inicializar_populacao(tamanho_populacao, num_variaveis, limites_inf, limites_sup);
    
    % Avaliar a população inicial
    [objetivos, ~] = avaliar_populacao(populacao, funcao_objetivo);
    
    % Gerar pontos de referência para o NSGA-III
    pontos_referencia = criar_pontos_referencia(size(objetivos, 2), 10);
    
    fronteira_pareto = [];
    conjunto_pareto = [];
    
    % Processo evolutivo
    for geracao = 1:maximo_geracoes
        fprintf('Processando a geração %d...\n', geracao);
        
        % 1. Classificação não-dominada
        niveis_fronteira = classificacao_nao_dominada(objetivos);
        
        % 2. Cálculo da distância de aglomeração (métrica melhorada)
        distancias_aglomeracao = calcular_distancia_aglomeracao(objetivos, niveis_fronteira);
        
        % 3. Seleção ambiental com base nos pontos de referência
        [populacao, objetivos] = selecao_ambiental(...
            populacao, objetivos, niveis_fronteira, distancias_aglomeracao, ...
            pontos_referencia, tamanho_populacao);
        
        % 4. Geração de descendentes usando estratégia GWO
        descendentes = gerar_descendentes(populacao, objetivos, ...
                                         limites_inf, limites_sup, prob_crossover, taxa_mutacao, geracao);
        
        % 5. Avaliação dos descendentes
        [objetivos_descendentes, ~] = avaliar_populacao(descendentes, funcao_objetivo);
        
        % 6. Combinação das populações
        populacao_combinada = [populacao; descendentes];
        objetivos_combinados = [objetivos; objetivos_descendentes];
        
        % 7. Seleção da próxima geração
        [populacao, objetivos] = selecionar_proxima_geracao(...
            populacao_combinada, objetivos_combinados, tamanho_populacao, pontos_referencia);
        
        % Atualizar a fronteira de Pareto
        indices_pareto = find(niveis_fronteira == 1);
        fronteira_pareto = objetivos(indices_pareto, :);
        conjunto_pareto = populacao(indices_pareto, :);
        
        % Exibir informações de progresso
        if mod(geracao, 5) == 0
            fprintf('  Geração: %d, Tamanho da fronteira de Pareto: %d\n', ...
                    geracao, size(fronteira_pareto, 1));
        end
    end
    
    % Extrair a fronteira de Pareto final
    niveis_finais = classificacao_nao_dominada(objetivos);
    indices_pareto_finais = find(niveis_finais == 1);
    fronteira_pareto = objetivos(indices_pareto_finais, :);
    conjunto_pareto = populacao(indices_pareto_finais, :);
    
    informacoes.numero_geracoes = maximo_geracoes;
    informacoes.tamanho_populacao = tamanho_populacao;
    informacoes.tamanho_pareto_final = size(fronteira_pareto, 1);
end

function novos_individuos = gerar_descendentes(populacao, objetivos, limites_inf, limites_sup, ...
                                               prob_crossover, taxa_mutacao, geracao_atual)
    % Gerar descendentes usando a estratégia do GWO
    tamanho_pop = size(populacao, 1);
    dimensoes = size(populacao, 2);
    novos_individuos = zeros(tamanho_pop, dimensoes);
    
    % Identificar líderes via classificação não-dominada (Alpha, Beta, Delta)
    niveis = classificacao_nao_dominada(objetivos);
    idx_alpha = find(niveis == 1, 1, 'first');
    idx_beta = find(niveis <= 2, 2, 'first');
    idx_delta = find(niveis <= 3, 3, 'first');
    
    % Parâmetro adaptativo do GWO
    valor_a = 2 * (1 - geracao_atual / 100);
    
    for i = 1:tamanho_pop
        if rand() < prob_crossover
            % Atualização usando o GWO
            r1 = rand(); r2 = rand();
            A = 2 * valor_a * r1 - valor_a;
            C = 2 * r2;
            
            D_alpha = abs(C * populacao(idx_alpha, :) - populacao(i, :));
            X1 = populacao(idx_alpha, :) - A * D_alpha;
            
            % Operação de crossover
            outro_pai = populacao(randi(tamanho_pop), :);
            ponto_corte = randi(dimensoes - 1);
            novos_individuos(i, 1:ponto_corte) = X1(1:ponto_corte);
            novos_individuos(i, ponto_corte+1:end) = outro_pai(ponto_corte+1:end);
        else
            novos_individuos(i, :) = populacao(i, :);
        end
        
        % Operação de mutação
        if rand() < taxa_mutacao
            ponto_mutacao = randi(dimensoes);
            novos_individuos(i, ponto_mutacao) = limites_inf(ponto_mutacao) + ...
                (limites_sup(ponto_mutacao) - limites_inf(ponto_mutacao)) * rand();
        end
        
        % Garantir que os valores estejam dentro dos limites
        novos_individuos(i, :) = max(min(novos_individuos(i, :), limites_sup), limites_inf);
    end
end

3. Classe do Modelo Transformer (ModeloTransformer.m)

classdef ModeloTransformer < handle
    % Modelo Transformer para regressão multivariável
    
    properties
        num_camadas
        num_cabecas
        dimensao_oculta
        taxa_dropout
        taxa_aprendizado
        tamanho_lote
        rede
        opcoes_treinamento
    end
    
    methods
        function obj = ModeloTransformer(parametros)
            % Inicializar os parâmetros do modelo
            obj.taxa_aprendizado = parametros.aprendizado;
            obj.num_camadas = parametros.camadas;
            obj.num_cabecas = parametros.cabecas;
            obj.dimensao_oculta = parametros.dimensao;
            obj.taxa_dropout = parametros.dropout;
            obj.tamanho_lote = parametros.lote;
            
            % Construir a arquitetura do modelo
            obj.construir_modelo();
        end
        
        function construir_modelo(obj)
            % Definir a estrutura do Transformer
            camadas = [
                sequenceInputLayer([], 'Name', 'camada_entrada')
                selfAttentionLayer(obj.num_cabecas, obj.dimensao_oculta, ...
                                  'Dropout', obj.taxa_dropout)
                fullyConnectedLayer(obj.dimensao_oculta, 'Name', 'camada_oculta1')
                reluLayer('Name', 'ativacao_relu')
                dropoutLayer(obj.taxa_dropout, 'Name', 'camada_dropout')
                fullyConnectedLayer(obj.dimensao_oculta, 'Name', 'camada_oculta2')
                fullyConnectedLayer(1, 'Name', 'camada_saida')
                regressionLayer('Name', 'camada_regressao')
            ];
            
            grafico_camadas = layerGraph(camadas);
            
            % Configurar as opções de treinamento
            opcoes = trainingOptions('adam', ...
                'MaxEpochs', 100, ...
                'MiniBatchSize', obj.tamanho_lote, ...
                'InitialLearnRate', obj.taxa_aprendizado, ...
                'LearnRateSchedule', 'piecewise', ...
                'LearnRateDropFactor', 0.85, ...
                'LearnRateDropPeriod', 15, ...
                'GradientThreshold', 1, ...
                'Verbose', false, ...
                'Plots', 'none');
            
            obj.rede = grafico_camadas;
            obj.opcoes_treinamento = opcoes;
        end
        
        function [modelo_treinado, historico] = treinar(obj, entradas, saidas)
            % Treinar o modelo
            [modelo_treinado, historico] = trainNetwork(entradas, saidas, ...
                                                        obj.rede, obj.opcoes_treinamento);
        end
        
        function previsoes = prever(obj, modelo, dados_teste)
            % Realizar previsões
            previsoes = predict(modelo, dados_teste);
        end
        
        function complexidade = calcular_complexidade(obj)
            % Estimar a complexidade do modelo (número de parâmetros)
            complexidade = obj.dimensao_oculta * obj.num_camadas * obj.num_cabecas * 500;
        end
    end
end

4. Função Multiobjetivo (funcao_multiobjetivo.m)

function objetivos = funcao_multiobjetivo(x, entradas_treino, saidas_treino, ...
                                         entradas_teste, saidas_teste, config_entrada, config_saida)
    % Função objetivo múltipla: minimizar erro de previsão e complexidade do modelo
    
    % Decodificar os parâmetros da solução
    parametros = struct();
    parametros.aprendizado = x(1);
    parametros.dimensao = round(x(2));
    parametros.cabecas = round(x(3));
    parametros.camadas = round(x(4));
    parametros.dropout = x(5);
    parametros.lote = round(x(6));
    
    % Inicializar e treinar o modelo Transformer
    modelo_transformer = ModeloTransformer(parametros);
    [modelo_treinado, ~] = modelo_transformer.treinar(entradas_treino, saidas_treino);
    
    % Realizar previsões nos dados de teste
    previsoes = modelo_transformer.prever(modelo_treinado, entradas_teste);
    
    % Desnormalizar as previsões e os valores reais
    previsoes_desnorm = mapminmax('reverse', previsoes', config_saida)';
    saidas_teste_desnorm = mapminmax('reverse', saidas_teste', config_saida)';
    
    % Objetivo 1: Erro quadrático médio (RMSE)
    rmse = sqrt(mean((previsoes_desnorm - saidas_teste_desnorm).^2, 'all'));
    
    % Objetivo 2: Complexidade do modelo
    complexidade = modelo_transformer.calcular_complexidade();
    
    % Retornar os objetivos (ambos devem ser minimizados)
    objetivos = [rmse, complexidade];
    
    % Exibir os resultados atuais
    fprintf('  RMSE: %.4f, Complexidade: %.0f\n', rmse, complexidade);
end

5. Funções Auxiliares (funcoes_auxiliares.m)

%% Funções auxiliares para o algoritmo NSGA-III

function populacao = inicializar_populacao(tamanho, dimensoes, limites_inf, limites_sup)
    % Gerar população inicial aleatória dentro dos limites
    populacao = zeros(tamanho, dimensoes);
    for idx = 1:tamanho
        populacao(idx, :) = limites_inf + (limites_sup - limites_inf) .* rand(1, dimensoes);
    end
end

function [objetivos, modelos] = avaliar_populacao(populacao, funcao_obj)
    % Avaliar todos os indivíduos da população
    tamanho_pop = size(populacao, 1);
    objetivos = zeros(tamanho_pop, 2);
    modelos = cell(tamanho_pop, 1);
    
    parfor i = 1:tamanho_pop
        objetivos(i, :) = funcao_obj(populacao(i, :));
    end
end

function niveis = classificacao_nao_dominada(objetivos)
    % Classificação não-dominada rápida para múltiplos objetivos
    [num_individuos, num_obj] = size(objetivos);
    niveis = zeros(num_individuos, 1);
    
    dominados_por = cell(num_individuos, 1);
    contagem_dominacao = zeros(num_individuos, 1);
    
    for p = 1:num_individuos
        dominados_por{p} = [];
        contagem_dominacao(p) = 0;
        for q = 1:num_individuos
            if eh_dominado(objetivos(p, :), objetivos(q, :))
                dominados_por{p} = [dominados_por{p}, q];
            elseif eh_dominado(objetivos(q, :), objetivos(p, :))
                contagem_dominacao(p) = contagem_dominacao(p) + 1;
            end
        end
        if contagem_dominacao(p) == 0
            niveis(p) = 1;
        end
    end
    
    contador = 1;
    fronteira_atual = find(niveis == contador);
    while ~isempty(fronteira_atual)
        proxima_fronteira = [];
        for p = fronteira_atual
            for q = dominados_por{p}
                contagem_dominacao(q) = contagem_dominacao(q) - 1;
                if contagem_dominacao(q) == 0
                    niveis(q) = contador + 1;
                    proxima_fronteira = [proxima_fronteira, q];
                end
            end
        end
        contador = contador + 1;
        fronteira_atual = proxima_fronteira;
    end
end

function dominado = eh_dominado(a, b)
    % Verificar se o vetor a é dominado por b
    nao_pior = all(a <= b);
    melhor = any(a < b);
    dominado = nao_pior && melhor;
end

function distancias = calcular_distancia_aglomeracao(objetivos, niveis)
    % Calcular distância de aglomeração para cada indivíduo
    [num_individuos, num_obj] = size(objetivos);
    distancias = zeros(num_individuos, 1);
    
    maximo_nivel = max(niveis);
    for nivel = 1:maximo_nivel
        indices_nivel = find(niveis == nivel);
        tamanho_nivel = length(indices_nivel);
        
        if tamanho_nivel > 0
            objs_nivel = objetivos(indices_nivel, :);
            [objs_ordenados, indices_ordenados] = sortrows(objs_nivel);
            
            for obj_idx = 1:num_obj
                valores_obj = objs_ordenados(:, obj_idx);
                amplitude = valores_obj(end) - valores_obj(1);
                
                if amplitude > 0
                    distancias(indices_nivel(indices_ordenados(1))) = inf;
                    distancias(indices_nivel(indices_ordenados(end))) = inf;
                    
                    for k = 2:tamanho_nivel-1
                        idx_original = indices_nivel(indices_ordenados(k));
                        distancias(idx_original) = distancias(idx_original) + ...
                            (valores_obj(k+1) - valores_obj(k-1)) / amplitude;
                    end
                end
            end
        end
    end
end

function pontos_ref = criar_pontos_referencia(num_objetivos, divisoes)
    % Gerar pontos de referência uniformemente distribuídos
    pontos_ref = [];
    % Implementação simplificada para dois objetivos
    if num_objetivos == 2
        pontos_ref = linspace(0, 1, divisoes+1)';
        pontos_ref = [pontos_ref, 1 - pontos_ref];
    end
end

6. Avaliação de Desempenho (avaliar_desempenho_modelo.m)

function avaliar_desempenho_modelo(previsoes, saidas_teste_norm, config_saida, saidas_teste_original)
    % Avaliar o desempenho do modelo treinado
    
    % Desnormalizar as previsões
    previsoes_desnorm = mapminmax('reverse', previsoes', config_saida)';
    
    % Calcular métricas de desempenho
    erro_rmse = sqrt(mean((previsoes_desnorm - saidas_teste_original).^2, 'all'));
    erro_mae = mean(abs(previsoes_desnorm - saidas_teste_original), 'all');
    r_quadrado = 1 - sum((previsoes_desnorm - saidas_teste_original).^2) / ...
                 sum((saidas_teste_original - mean(saidas_teste_original)).^2);
    
    % Exibir os resultados
    fprintf('\n========= Relatório de Desempenho =========\n');
    fprintf('RMSE: %.4f\n', erro_rmse);
    fprintf('MAE:  %.4f\n', erro_mae);
    fprintf('R²:   %.4f\n', r_quadrado);
    fprintf('============================================\n');
    
    % Gerar gráficos de análise
    figura = figure('Position', [100, 100, 1200, 400]);
    
    subplot(1, 3, 1);
    scatter(saidas_teste_original(:,1), previsoes_desnorm(:,1), 20, 'filled');
    hold on;
    plot([min(saidas_teste_original(:,1)), max(saidas_teste_original(:,1))], ...
         [min(saidas_teste_original(:,1)), max(saidas_teste_original(:,1))], 'r--', 'LineWidth', 2);
    xlabel('Valor Real');
    ylabel('Valor Previsto');
    title('Comparação entre Valores Reais e Previstos');
    grid on;
    
    subplot(1, 3, 2);
    residuos = previsoes_desnorm - saidas_teste_original;
    histogram(residuos(:,1), 20);
    xlabel('Resíduos');
    ylabel('Frequência');
    title('Distribuição dos Resíduos');
    grid on;
    
    subplot(1, 3, 3);
    plot(saidas_teste_original(1:100,1), 'b-', 'LineWidth', 1.5);
    hold on;
    plot(previsoes_desnorm(1:100,1), 'r--', 'LineWidth', 1.5);
    xlabel('Índice da Amostra');
    ylabel('Valor');
    title('Série Temporal: Real vs. Previsto');
    legend('Valores Reais', 'Valores Previstos', 'Location', 'best');
    grid on;
end

Instruções de Uso

Para uitlizar esta estrutura:

1. Prepare os dados multivariáveis em um arquivo MATLAB com variáveis para treinamento e teste.
2. Ajuste os parâmetros de otimização no programa principal, como limites das variáveis e configurações do algoritmo.
3. Execute o programa principal para iniciar o processo de otimização multiobjetivo.
4. Analise a fronteira de Pareto gerada para selecionar o compromisso desejado entre erro e complexidade.
5. Use o modelo final treinado para realizar previsões em novos dados.