Modelagem de Incertezas em Sistemas Elétricos

A otimização de sistemas elétricos modernos enfrenta desafios significativos devido à variabilidade inerente das fontes renováveis e a demanda flutuante dos consumidores. Este artigo apresenta técnicas avançadas de modelagem probabilística e otimização robusta implementadas em MATLAB, permitindo engineers desenvolverem estratégias de despacho mais resilientes.

1. Metodologias de Modelagem de Incertezas

As principais abordagens para representar incertezas na geração renovável e na carga de energia são:

Este estudo concentra-se nas duas abordagens predominantes: porgramação estocástica baseada em cenários e otimização robusta.

2. Técnicas de Geração e Redução de Cenários

2.1 Modelagem de Fontes de Incerteza

O código a seguir demonstra a modelagem estocástica para diferentes fontes de incerteza:

%% Modelagem de Incertezas para Vento, Solar e Carga
clc; clear all; close all;

%% Configuração de Parâmetros
horizonte_temporal = 24;          % Horizonte de 24 horas
quantidade_cenarios_inicial = 100;  % Cenários iniciais
quantidade_cenarios_finais = 10;    % Cenários após redução
indices_tempo = 1:horizonte_temporal;

%% 2.1.1 Modelagem de Geração Eólica
% Curva média típica (perfil diário)
potencia_eolica_media = 50 + 20 * sin(2*pi*(indices_tempo-6)/24);  % MW

% Desvio padrão (maior variabilidade durante o dia)
potencia_eolica_desvio = 0.2 * potencia_eolica_media;  % 20% de variação

% Geração de cenários (distribuição normal)
cenarios_eolicos = zeros(quantidade_cenarios_inicial, horizonte_temporal);
for h = 1:horizonte_temporal
    cenarios_eolicos(:, h) = potencia_eolica_media(h) + ...
        potencia_eolica_desvio(h) * randn(quantidade_cenarios_inicial, 1);
    % Limitação entre [0, capacidade máxima]
    cenarios_eolicos(:, h) = max(0, min(cenarios_eolicos(:, h), 100));
end

%% 2.1.2 Modelagem de Geração Solar Fotovoltaica
potencia_solar_media = 60 * sin(pi*(indices_tempo-6)/12);
potencia_solar_media(potencia_solar_media < 0) = 0;

potencia_solar_desvio = 0.15 * potencia_solar_media;
potencia_solar_desvio(isnan(potencia_solar_desvio)) = 0.1;

cenarios_solares = zeros(quantidade_cenarios_inicial, horizonte_temporal);
for h = 1:horizonte_temporal
    if potencia_solar_media(h) > 0
        cenarios_solares(:, h) = potencia_solar_media(h) + ...
            potencia_solar_desvio(h) * randn(quantidade_cenarios_inicial, 1);
        cenarios_solares(:, h) = max(0, min(cenarios_solares(:, h), 80));
    end
end

%% 2.1.3 Modelagem de Incerteza na Carga
carga_media = 150 + 50 * sin(2*pi*(indices_tempo-14)/24);  % MW
carga_desvio = 0.1 * carga_media;

cenarios_carga = zeros(quantidade_cenarios_inicial, horizonte_temporal);
for h = 1:horizonte_temporal
    cenarios_carga(:, h) = carga_media(h) + carga_desvio(h) * randn(quantidade_cenarios_inicial, 1);
    cenarios_carga(:, h) = max(100, cenarios_carga(:, h));
end

%% 2.1.4 Probabilidades Iniciais (iguais)
prob_iniciais = ones(quantidade_cenarios_inicial, 1) / quantidade_cenarios_inicial;

%% 2.2 Redução de Cenários via Clustering K-Means
function [cenarios_reduzidos, prob_reduzidas, indices_grupo] = ...
    reduzir_cenarios(cenarios_originais, prob_originais, n_clusters)
    
    [n_scenarios, n_horas] = size(cenarios_originais);
    
    % Aplicação do algoritmo K-means
    [indices_grupo, centroides] = kmeans(cenarios_originais, n_clusters, ...
        'MaxIter', 1000, 'Replicates', 10);
    
    cenarios_reduzidos = centroides;
    
    % Cálculo de probabilidade por grupo
    prob_reduzidas = zeros(n_clusters, 1);
    for k = 1:n_clusters
        membros_grupo = (indices_grupo == k);
        prob_reduzidas(k) = sum(prob_originais(membros_grupo));
    end
    
    % Visualização dos grupos
    figure;
    for k = 1:min(4, n_clusters)
        subplot(2,2,k);
        dados_grupo = cenarios_originais(indices_grupo == k, :);
        plot(indices_tempo, dados_grupo', 'Color', [0.8, 0.8, 0.8], 'LineWidth', 0.5);
        hold on;
        plot(indices_tempo, cenarios_reduzidos(k, :), 'r-', 'LineWidth', 2);
        title(sprintf('Grupo %d (Prob=%.3f)', k, prob_reduzidas(k)));
        xlabel('Hora');
        ylabel('Potência (MW)');
        grid on;
    end
end

% Aplicação da redução de cenários
[cenario_eolico_final, prob_eolica_final, idx_eolico] = ...
    reduzir_cenarios(cenarios_eolicos, prob_iniciais, quantidade_cenarios_finais);
[cenario_solar_final, prob_solar_final, idx_solar] = ...
    reduzir_cenarios(cenarios_solares, prob_iniciais, quantidade_cenarios_finais);
[cenario_carga_final, prob_carga_final, idx_carga] = ...
    reduzir_cenarios(cenarios_carga, prob_iniciais, quantidade_cenarios_finais);

%% 2.3 Geração de Cenários Conjuntos (considerando correlação)
function [cenarios_conjuntos, prob_conjuntas] = ...
    gerar_cenarios_conjuntos(cenarios_eolicos, prob_eolicos, ...
    cenarios_solares, prob_solares, cenarios_carga, prob_cargas)
    
    n_eolico = size(cenarios_eolicos, 1);
    n_solar = size(cenarios_solares, 1);
    n_carga = size(cenarios_carga, 1);
    
    % Assumindo independência, gerar todas as combinações
    n_conjunto = n_eolico * n_solar * n_carga;
    cenarios_conjuntos = zeros(n_conjunto, 3, 24);
    prob_conjuntas = zeros(n_conjunto, 1);
    
    idx = 1;
    for i = 1:n_eolico
        for j = 1:n_solar
            for k = 1:n_carga
                cenarios_conjuntos(idx, 1, :) = cenarios_eolicos(i, :);
                cenarios_conjuntos(idx, 2, :) = cenarios_solares(j, :);
                cenarios_conjuntos(idx, 3, :) = cenarios_carga(k, :);
                prob_conjuntas(idx) = prob_eolicos(i) * prob_solares(j) * prob_cargas(k);
                idx = idx + 1;
            end
        end
    end
    
    % Redução adicional dos cenários conjuntos
    cenarios_conjuntos_flat = reshape(cenarios_conjuntos, n_conjunto, []);
    [cenarios_reduzidos, prob_reduzidas] = reduzir_cenarios(...
        cenarios_conjuntos_flat, prob_conjuntas, 20);
    
    % Restauração da dimensão original
    cenarios_conjuntos = reshape(cenarios_reduzidos, [], 3, 24);
    prob_conjuntas = prob_reduzidas;
end

% Geração dos cenários conjuntos
[cenarios_juntos, prob_juntas] = gerar_cenarios_conjuntos(...
    cenario_eolico_final, prob_eolica_final, ...
    cenario_solar_final, prob_solar_final, ...
    cenario_carga_final, prob_carga_final);

n_juntos = size(cenarios_juntos, 1);
fprintf('Cenários conjuntos após redução: %d\n', n_juntos);

3. Modelo de Programação Estocástica

3.1 Estrutura de Dois Estágios

%% Programação Estocástica de Dois Estágios: Despacho com Incertezas
clear; clc;

%% 1. Parâmetros do Sistema
horizonte_tempo = 24;
num_geradores = 3;       % Número de unidades térmicas
num_parques_eolicos = 1; % Número de parques eólicos
num_parques_solares = 1; % Número de parques solares
num_cenarios = size(cenarios_juntos, 1);

%% 2. Parâmetros das Unidades Geradoras
% Unidade 1: Base
gerador(1).pot_min = 50;     % MW
gerador(1).pot_max = 200;    % MW
gerador(1).coef_a = 0.01;    % $/MW^2
gerador(1).coef_b = 20;      % $/MW
gerador(1).coef_c = 100;     % $
gerador(1).rampa_subida = 100;  % MW/h
gerador(1).rampa_descida = 100; % MW/h
gerador(1).tempo_min_ligada = 4;     % h
gerador(1).tempo_min_desligada = 4;  % h
gerador(1).custo_partida = 500;   % $
gerador(1).custo_parada = 100;    % $

% Unidade 2: Média
gerador(2).pot_min = 20;
gerador(2).pot_max = 100;
gerador(2).coef_a = 0.02;
gerador(2).coef_b = 30;
gerador(2).coef_c = 50;
gerador(2).rampa_subida = 80;
gerador(2).rampa_descida = 80;
gerador(2).tempo_min_ligada = 2;
gerador(2).tempo_min_desligada = 2;
gerador(2).custo_partida = 300;
gerador(2).custo_parada = 50;

% Unidade 3: Pico
gerador(3).pot_min = 0;
gerador(3).pot_max = 50;
gerador(3).coef_a = 0.05;
gerador(3).coef_b = 50;
gerador(3).coef_c = 0;
gerador(3).rampa_subida = 50;
gerador(3).rampa_descida = 50;
gerador(3).tempo_min_ligada = 1;
gerador(3).tempo_min_desligada = 1;
gerador(3).custo_partida = 200;
gerador(3).custo_parada = 0;

%% 3. Parâmetros das Renováveis
% Penalidades por corte
penalidade_corte_eolico = 50;   % $/MWh
penalidade_corte_solar = 40;    % $/MWh

% Reserva girante requerida
reserva_girante = 0.1;  % 10% da carga

%% 4. Construção do Modelo de Dois Estágios
% Estágio 1: Estado de ligado/desligado (binário)
% Estágio 2: Geração, geração efetiva renovável, cortes (contínuo)

addpath(genpath('yalmip'));
addpath(genpath('gurobi'));

% Variáveis de decisão
estado_operacao = binvar(num_geradores, horizonte_tempo, 'full');
indicador_partida = binvar(num_geradores, horizonte_tempo, 'full');
indicador_parada = binvar(num_geradores, horizonte_tempo, 'full');

% Variáveis do segundo estágio (por cenário)
potencia_gerada = sdpvar(num_geradores, horizonte_tempo, num_cenarios, 'full');
potencia_eolica_real = sdpvar(horizonte_tempo, num_cenarios, 'full');
potencia_solar_real = sdpvar(horizonte_tempo, num_cenarios, 'full');
corte_eolico = sdpvar(horizonte_tempo, num_cenarios, 'full');
corte_solar = sdpvar(horizonte_tempo, num_cenarios, 'full');

%% 5. Função Objetivo: Minimizar Custo Total Esperado
custo_total = 0;

% Custo de partida e parada (Estágio 1)
for i = 1:num_geradores
    for t = 1:horizonte_tempo
        custo_total = custo_total + gerador(i).custo_partida * indicador_partida(i,t) + ...
                     gerador(i).custo_parada * indicador_parada(i,t);
    end
end

% Valor esperado do segundo estágio
for s = 1:num_cenarios
    prob_cenario = prob_juntas(s);
    
    for t = 1:horizonte_tempo
        % Obtenção dos dados do cenário
        eolico_disponivel = cenarios_juntos(s, 1, t);
        solar_disponivel = cenarios_juntos(s, 2, t);
        demanda_carga = cenarios_juntos(s, 3, t);
        
        % Custo de geração
        for i = 1:num_geradores
            custo_total = custo_total + prob_cenario * (...
                gerador(i).coef_a * potencia_gerada(i,t,s)^2 + ...
                gerador(i).coef_b * potencia_gerada(i,t,s) + ...
                gerador(i).coef_c * estado_operacao(i,t));
        end
        
        % Penalidades por corte
        custo_total = custo_total + prob_cenario * (...
            penalidade_corte_eolico * corte_eolico(t,s) + ...
            penalidade_corte_solar * corte_solar(t,s));
    end
end

%% 6. Restrições
restricoes = [];

% 6.1 Lógica de partida/parada
for i = 1:num_geradores
    for t = 1:horizonte_tempo
        if t > 1
            restricoes = [restricoes, indicador_partida(i,t) >= estado_operacao(i,t) - estado_operacao(i,t-1)];
            restricoes = [restricoes, indicador_parada(i,t) >= estado_operacao(i,t-1) - estado_operacao(i,t)];
        else
            restricoes = [restricoes, indicador_partida(i,1) >= estado_operacao(i,1)];
            restricoes = [restricoes, indicador_parada(i,1) >= 0 - estado_operacao(i,1)];
        end
        
        restricoes = [restricoes, indicador_partida(i,t) <= estado_operacao(i,t)];
        restricoes = [restricoes, indicador_parada(i,t) <= 1 - estado_operacao(i,t)];
        restricoes = [restricoes, indicador_partida(i,t) + indicador_parada(i,t) <= 1];
    end
end

% 6.2 Tempos mínimos de operação
for i = 1:num_geradores
    for t = 1:horizonte_tempo
        % Tempo mínimo ligado
        if t + gerador(i).tempo_min_ligada - 1 <= horizonte_tempo
            restricoes = [restricoes, ...
                sum(estado_operacao(i, t:min(t+gerador(i).tempo_min_ligada-1, horizonte_tempo))) >= ...
                gerador(i).tempo_min_ligada * indicador_partida(i,t)];
        end
        
        % Tempo mínimo desligado
        if t + gerador(i).tempo_min_desligada - 1 <= horizonte_tempo
            restricoes = [restricoes, ...
                sum(1 - estado_operacao(i, t:min(t+gerador(i).tempo_min_desligada-1, horizonte_tempo))) >= ...
                gerador(i).tempo_min_desligada * indicador_parada(i,t)];
        end
    end
end

% 6.3 Restrições por cenário
for s = 1:num_cenarios
    for t = 1:horizonte_tempo
        eolico_disp = cenarios_juntos(s, 1, t);
        solar_disp = cenarios_juntos(s, 2, t);
        demanda = cenarios_juntos(s, 3, t);
        
        % Balanço de potência
        geracao_total = sum(potencia_gerada(:,t,s)) + potencia_eolica_real(t,s) + potencia_solar_real(t,s);
        restricoes = [restricoes, geracao_total == demanda];
        
        % Limites de geração
        for i = 1:num_geradores
            restricoes = [restricoes, ...
                gerador(i).pot_min * estado_operacao(i,t) <= potencia_gerada(i,t,s) <= ...
                gerador(i).pot_max * estado_operacao(i,t)];
        end
        
        % Restrições de rampa
        if t > 1
            for i = 1:num_geradores
                restricoes = [restricoes, ...
                    potencia_gerada(i,t,s) - potencia_gerada(i,t-1,s) <= gerador(i).rampa_subida];
                restricoes = [restricoes, ...
                    potencia_gerada(i,t-1,s) - potencia_gerada(i,t,s) <= gerador(i).rampa_descida];
            end
        end
        
        % Restrições de renováveis
        restricoes = [restricoes, 0 <= potencia_eolica_real(t,s) <= eolico_disp];
        restricoes = [restricoes, 0 <= potencia_solar_real(t,s) <= solar_disp];
        
        % Definição de corte
        restricoes = [restricoes, ...
            corte_eolico(t,s) == eolico_disp - potencia_eolica_real(t,s)];
        restricoes = [restricoes, ...
            corte_solar(t,s) == solar_disp - potencia_solar_real(t,s)];
        
        % Reserva girante
        reserva_rapida = 0;
        for i = 1:num_geradores
            reserva_rapida = reserva_rapida + min(gerador(i).rampa_subida, ...
                gerador(i).pot_max * estado_operacao(i,t) - potencia_gerada(i,t,s));
        end
        restricoes = [restricoes, reserva_rapida >= reserva_girante * demanda];
    end
end

%% 7. Resolução
ops = sdpsettings('solver', 'gurobi', 'verbose', 1, 'showprogress', 1);
ops.gurobi.TimeLimit = 3600;
ops.gurobi.MIPGap = 0.01;

tic;
diagnosticos = otimizar(restricoes, custo_total, ops);
tempo_solucao = toc;

if diagnosticos.problem == 0
    disp('=== Solução Encontrada ===');
    fprintf('Tempo de solução: %.2f segundos\n', tempo_solucao);
    fprintf('Custo objetivo: %.2f $\n', valor(custo_total));
    
    % Extração de resultados
    estado_otimo = valor(estado_operacao);
    potencia_otima = valor(potencia_gerada);
    eolica_otima = valor(potencia_eolica_real);
    solar_otima = valor(potencia_solar_real);
    corte_eolico_otimo = valor(corte_eolico);
    corte_solar_otimo = valor(corte_solar);
    
    % Cálculo de valores esperados
    potencia_esperada = zeros(num_geradores, horizonte_tempo);
    for i = 1:num_geradores
        for t = 1:horizonte_tempo
            for s = 1:num_cenarios
                potencia_esperada(i,t) = potencia_esperada(i,t) + prob_juntas(s) * potencia_otima(i,t,s);
            end
        end
    end
    
else
    disp('Falha na solução');
    disp(yalmiperror(diagnosticos.problem));
end

%% 8. Visualização de Resultados
visualizar_resultados(estado_otimo, potencia_esperada, eolica_otima, solar_otima, ...
    corte_eolico_otimo, corte_solar_otimo, cenarios_juntos, prob_juntas);

3.2 Função de Visualização

function visualizar_resultados(estado_opt, potencia_esp, eolica_opt, solar_opt, ...
    corte_eolico_opt, corte_solar_opt, cenarios, probabilidades)
    
    [n_geradores, n_horas] = size(potencia_esp);
    n_cenarios = size(cenarios, 1);
    tempo = 1:n_horas;
    
    figure('Position', [100, 100, 1400, 900]);
    
    % Subplot 1: Estado de operação
    subplot(3,3,1);
    imagesc(tempo, 1:n_geradores, estado_opt);
    colormap(flipud(gray));
    colorbar;
    xlabel('Tempo (h)');
    ylabel('Gerador');
    title('Estado de Operação');
    set(gca, 'YTick', 1:n_geradores);
    
    % Subplot 2: Geração esperada
    subplot(3,3,2);
    area(tempo, potencia_esp');
    legend('Gerador 1', 'Gerador 2', 'Gerador 3', 'Location', 'best');
    xlabel('Tempo (h)');
    ylabel('Potência (MW)');
    title('Geração Esperada');
    grid on;
    
    % Subplot 3: Renováveis
    subplot(3,3,3);
    hold on;
    
    eolica_esperada = zeros(1, n_horas);
    solar_esperado = zeros(1, n_horas);
    for t = 1:n_horas
        for s = 1:n_cenarios
            eolica_esperada(t) = eolica_esperada(t) + probabilidades(s) * eolica_opt(t,s);
            solar_esperado(t) = solar_esperado(t) + probabilidades(s) * solar_opt(t,s);
        end
    end
    
    eolica_disponivel_esp = zeros(1, n_horas);
    solar_disponivel_esp = zeros(1, n_horas);
    for t = 1:n_horas
        for s = 1:n_cenarios
            eolica_disponivel_esp(t) = eolica_disponivel_esp(t) + probabilidades(s) * cenarios(s,1,t);
            solar_disponivel_esp(t) = solar_disponivel_esp(t) + probabilidades(s) * cenarios(s,2,t);
        end
    end
    
    plot(tempo, eolica_esperada, 'b-', 'LineWidth', 2, 'DisplayName', 'Eólica Efetiva');
    plot(tempo, solar_esperado, 'g-', 'LineWidth', 2, 'DisplayName', 'Solar Efetivo');
    plot(tempo, eolica_disponivel_esp, 'b--', 'LineWidth', 1, 'DisplayName', 'Eólica Disponível');
    plot(tempo, solar_disponivel_esp, 'g--', 'LineWidth', 1, 'DisplayName', 'Solar Disponível');
    
    xlabel('Tempo (h)');
    ylabel('Potência (MW)');
    title('Geração Renovável');
    legend('Location', 'best');
    grid on;
    
    % Subplot 4: Cortes
    subplot(3,3,4);
    corte_eolico_esp = zeros(1, n_horas);
    corte_solar_esp = zeros(1, n_horas);
    for t = 1:n_horas
        for s = 1:n_cenarios
            corte_eolico_esp(t) = corte_eolico_esp(t) + probabilidades(s) * corte_eolico_opt(t,s);
            corte_solar_esp(t) = corte_solar_esp(t) + probabilidades(s) * corte_solar_opt(t,s);
        end
    end
    
    barras = [corte_eolico_esp; corte_solar_esp]';
    bar(tempo, barras, 'stacked');
    xlabel('Tempo (h)');
    ylabel('Energia Cortada (MWh)');
    title('Corte de Renováveis');
    legend('Corte Eólico', 'Corte Solar', 'Location', 'best');
    grid on;
    
    % Subplot 5: Balanço em cenário típico
    subplot(3,3,5);
    [~, idx_maior_prob] = max(probabilidades);
    
    carga_cenario = squeeze(cenarios(idx_maior_prob, 3, :))';
    geracao_total = sum(potencia_optima(:,:,idx_maior_prob), 1) + ...
                eolica_opt(:,idx_maior_prob)' + solar_opt(:,idx_maior_prob)';
    
    area(tempo, [potencia_optima(:,:,idx_maior_prob)', ...
                   eolica_opt(:,idx_maior_prob)', ...
                   solar_opt(:,idx_maior_prob)']);
    hold on;
    plot(tempo, carga_cenario, 'k-', 'LineWidth', 2, 'DisplayName', 'Carga');
    xlabel('Tempo (h)');
    ylabel('Potência (MW)');
    title('Balanço de Potência - Cenário Típico');
    legend('G1', 'G2', 'G3', 'Eólica', 'Solar', 'Carga', 'Location', 'best');
    grid on;
    
    % Subplot 6: Reserva
    subplot(3,3,6);
    reserva_disponivel = zeros(1, n_horas);
    reserva_necessaria = zeros(1, n_horas);
    
    for t = 1:n_horas
        carga_exp = 0;
        for s = 1:n_cenarios
            carga_exp = carga_exp + probabilidades(s) * cenarios(s,3,t);
        end
        reserva_necessaria(t) = 0.1 * carga_exp;
        
        for s = 1:n_cenarios
            for i = 1:n_geradores
                reserva_disponivel(t) = reserva_disponivel(t) + probabilidades(s) * ...
                    min(gerador(i).rampa_subida, gerador(i).pot_max*estado_opt(i,t)-potencia_optima(i,t,s));
            end
        end
    end
    
    bar(tempo, [reserva_disponivel; reserva_necessaria]', 'grouped');
    xlabel('Tempo (h)');
    ylabel('Reserva (MW)');
    title('Reserva Girante');
    legend('Disponível', 'Requerida', 'Location', 'best');
    grid on;
    
    % Subplot 7: Distribuição de probabilidade
    subplot(3,3,7);
    [prob_ord, idx_ord] = sort(probabilidades, 'descend');
    bar(1:length(prob_ord), prob_ord);
    xlabel('Cenário (ordenado)');
    ylabel('Probabilidade');
    title('Distribuição de Probabilidade');
    grid on;
    
    % Subplot 8: Incerteza da carga
    subplot(3,3,8);
    carga_min = min(squeeze(cenarios(:,3,:)), [], 1);
    carga_max = max(squeeze(cenarios(:,3,:)), [], 1);
    carga_media = mean(squeeze(cenarios(:,3,:)), 1);
    
    fill([tempo, fliplr(tempo)], [carga_max, fliplr(carga_min)], ...
        [0.8, 0.8, 0.8], 'FaceAlpha', 0.3, 'EdgeColor', 'none');
    hold on;
    plot(tempo, carga_media, 'b-', 'LineWidth', 2);
    xlabel('Tempo (h)');
    ylabel('Carga (MW)');
    title('Faixa de Incerteza da Carga');
    legend('Faixa', 'Média', 'Location', 'best');
    grid on;
    
    % Subplot 9: Custos
    subplot(3,3,9);
    rotulos_custo = {'Combustível', 'Partida/Parada', 'Penalidades'};
    valores_custo = [70, 15, 5];
    bar(valores_custo);
    set(gca, 'XTickLabel', rotulos_custo);
    ylabel('Custo (k$)');
    title('Decomposição de Custos');
    grid on;
end

4. Modelo de Otimização Robusta

4.1 Estrutura Adaptativa de Otimização Robusta

%% Otimização Robusta Adaptativa: Pior Caso de Incertezas
clear; clc;

%% 1. Parâmetros do Sistema
horizonte_tempo = 24;
num_geradores = 3;
num_parques_eolicos = 1;
num_parques_solares = 1;

%% 2. Definição do Conjunto de Incertezas
% Modelo de incerteza por intervalo
carga_nominal = [120, 115, 110, 105, 100, 110, 130, 150, 170, 180, ...
                175, 170, 165, 160, 155, 160, 170, 185, 190, 180, ...
                170, 160, 140, 125];  % MW

eolica_nominal = [45, 48, 50, 52, 55, 60, 65, 70, 75, 80, ...
                78, 75, 72, 70, 68, 65, 62, 60, 58, 55, ...
                52, 50, 48, 46];  % MW

solar_nominal = [0, 0, 0, 0, 0, 5, 15, 30, 45, 60, ...
              70, 75, 80, 75, 70, 60, 45, 30, 15, 5, ...
              0, 0, 0, 0];  % MW

% Intervalos de incerteza (±percentual)
incerteza_carga = 0.15;    % ±15%
incerteza_eolica = 0.30;   % ±30%
incerteza_solar = 0.25;    % ±25%

% Parâmetros do orçamento de incertezas
orcamento_carga = 8;
orcamento_eolica = 6;
orcamento_solar = 4;

%% 3. Modelagem de Otimização Robusta
yalmip('clear');

% Variáveis de primeiro estágio: Estado de operação
estado_oper = binvar(num_geradores, horizonte_tempo, 'full');
partida = binvar(num_geradores, horizonte_tempo, 'full');
parada = binvar(num_geradores, horizonte_tempo, 'full');

% Variáveis de segundo estágio: Geração
potencia_ger = sdpvar(num_geradores, horizonte_tempo, 'full');

% Variáveis de incerteza
delta_carga = sdpvar(1, horizonte_tempo, 'full');
delta_eolica = sdpvar(1, horizonte_tempo, 'full');
delta_solar = sdpvar(1, horizonte_tempo, 'full');

% Definição do conjunto de incertezas
conjunto_incertezas = [];
for t = 1:horizonte_tempo
    conjunto_incertezas = [conjunto_incertezas, -1 <= delta_carga(t) <= 1];
    conjunto_incertezas = [conjunto_incertezas, -1 <= delta_eolica(t) <= 1];
    conjunto_incertezas = [conjunto_incertezas, -1 <= delta_solar(t) <= 1];
end

% Restrições de orçamento
conjunto_incertezas = [conjunto_incertezas, ...
    norm(delta_carga, 1) <= orcamento_carga, ...
    norm(delta_eolica, 1) <= orcamento_eolica, ...
    norm(delta_solar, 1) <= orcamento_solar];

% Valores reais = nominal + incerteza
carga_real = carga_nominal + incerteza_carga * carga_nominal .* delta_carga;
eolica_real = eolica_nominal + incerteza_eolica * eolica_nominal .* delta_eolica;
solar_real = solar_nominal + incerteza_solar * solar_nominal .* delta_solar;

%% 4. Restrições
restricoes = [];

% 4.1 Lógica de partida/parada
for i = 1:num_geradores
    for t = 1:horizonte_tempo
        if t > 1
            restricoes = [restricoes, partida(i,t) >= estado_oper(i,t) - estado_oper(i,t-1)];
            restricoes = [restricoes, parada(i,t) >= estado_oper(i,t-1) - estado_oper(i,t)];
        else
            restricoes = [restricoes, partida(i,1) >= estado_oper(i,1)];
            restricoes = [restricoes, parada(i,1) >= 0 - estado_oper(i,1)];
        end
        restricoes = [restricoes, 0 <= partida(i,t) <= 1, 0 <= parada(i,t) <= 1];
    end
end

% 4.2 Limites de geração
for i = 1:num_geradores
    for t = 1:horizonte_tempo
        restricoes = [restricoes, ...
            gerador(i).pot_min * estado_oper(i,t) <= potencia_ger(i,t) <= gerador(i).pot_max * estado_oper(i,t)];
    end
end

% 4.3 Restrições de rampa
for i = 1:num_geradores
    for t = 2:horizonte_tempo
        restricoes = [restricoes, potencia_ger(i,t) - potencia_ger(i,t-1) <= gerador(i).rampa_subida];
        restricoes = [restricoes, potencia_ger(i,t-1) - potencia_ger(i,t) <= gerador(i).rampa_descida];
    end
end

% 4.4 Balanço de potência (robusto)
balanceamento = [];
for t = 1:horizonte_tempo
    geracao_total = sum(potencia_ger(:,t)) + eolica_real(t) + solar_real(t);
    balanceamento = [balanceamento, geracao_total == carga_real(t)];
end
restricoes = [restricoes, implica(conjunto_incertezas, balanceamento)];

% 4.5 Restrições de reserva
reservas = [];
for t = 1:horizonte_tempo
    reserva_disponivel = 0;
    for i = 1:num_geradores
        reserva_disponivel = reserva_disponivel + ...
            min(gerador(i).rampa_subida, gerador(i).pot_max*estado_oper(i,t) - potencia_ger(i,t));
    end
    
    % Pior caso de necessidade de reserva
    pior_caso_carga = carga_nominal(t) * (1 + incerteza_carga);
    necessidade_reserva = 0.1 * pior_caso_carga;
    
    reservas = [reservas, reserva_disponivel >= necessidade_reserva];
end
restricoes = [restricoes, reservas];

%% 5. Função Objetivo
% Minimizar custo no pior caso
custo_total = 0;

% Custos de partida/parada
for i = 1:num_geradores
    for t = 1:horizonte_tempo
        custo_total = custo_total + gerador(i).custo_partida * partida(i,t) + ...
                     gerador(i).custo_parada * parada(i,t);
    end
end

% Custos de geração
custo_geracao = 0;
for i = 1:num_geradores
    for t = 1:horizonte_tempo
        custo_geracao = custo_geracao + gerador(i).coef_a * potencia_ger(i,t)^2 + ...
                       gerador(i).coef_b * potencia_ger(i,t) + gerador(i).coef_c * estado_oper(i,t);
    end
end

% Penalidades por corte (pior caso)
corte_eolico_penal = 0;
corte_solar_penal = 0;

for t = 1:horizonte_tempo
    corte_eolico_penal = corte_eolico_penal + penalidade_corte_eolico * ...
        max(0, eolica_real(t) - (carga_real(t) - sum(potencia_ger(:,t)) - solar_real(t)));
    corte_solar_penal = custo_total + corte_solar_penal + penalidade_corte_solar * ...
        max(0, solar_real(t) - (carga_real(t) - sum(potencia_ger(:,t)) - eolica_real(t)));
end

custo_total = custo_total + custo_geracao + corte_eolico_penal + corte_solar_penal;

%% 6. Resolução do Problema Robusto
ops = sdpsettings('solver', 'gurobi', 'robust.lplp', 'duality');
ops.verbose = 1;
ops.gurobi.TimeLimit = 1800;

tic;
diagnosticos = otimizar([restricoes, conjunto_incertezas], custo_total, ops);
tempo_sol = toc;

if diagnosticos.problem == 0
    disp('=== Otimização Robusta Concluída ===');
    fprintf('Tempo: %.2f segundos\n', tempo_sol);
    fprintf('Custo no pior caso: %.2f $\n', valor(custo_total));
    
    % Extração
    estado_robusto = valor(estado_oper);
    potencia_robusta = valor(potencia_ger);
    
    % Análise de robustez
    verificar_robustez(estado_robusto, potencia_robusta, carga_nominal, eolica_nominal, solar_nominal, ...
        incerteza_carga, incerteza_eolica, incerteza_solar, gerador);
    
else
    disp('Falha na otimização robusta');
    disp(yalmiperror(diagnosticos.problem));
end

4.2 Função de Análise de Robustez

function verificar_robustez(estado_otimo, potencia_otima, carga_nom, eolica_nom, solar_nom, ...
    incerteza_carga, incerteza_eolica, incerteza_solar, geradores)
    
    n_horas = length(carga_nom);
    n_geradores = size(potencia_otima, 1);
    tempo = 1:n_horas;
    
    % Geração de cenários para validação
    n_monte_carlo = 1000;
    contagem_viavel = 0;
    distribuicao_custos = zeros(n_monte_carlo, 1);
    
    for mc = 1:n_monte_carlo
        % Geração aleatória de realizações
        xi_carga = 2*rand(1, n_horas) - 1;
        xi_eolica = 2*rand(1, n_horas) - 1;
        xi_solar = 2*rand(1, n_horas) - 1;
        
        % Aplicação de orçamento
        if sum(abs(xi_carga)) > 8
            xi_carga = xi_carga * 8 / sum(abs(xi_carga));
        end
        if sum(abs(xi_eolica)) > 6
            xi_eolica = xi_eolica * 6 / sum(abs(xi_eolica));
        end
        if sum(abs(xi_solar)) > 4
            xi_solar = xi_solar * 4 / sum(abs(xi_solar));
        end
        
        % Cálculo de valores reais
        carga_atual = carga_nom + incerteza_carga * carga_nom .* xi_carga;
        eolica_atual = eolica_nom + incerteza_eolica * eolica_nom .* xi_eolica;
        solar_atual = solar_nom + incerteza_solar * solar_nom .* xi_solar;
        
        % Verificação de viabilidade
        viavel = true;
        
        for t = 1:n_horas
            % Balanço de potência
            geracao_total = sum(potencia_otima(:,t)) + eolica_atual(t) + solar_atual(t);
            if abs(geracao_total - carga_atual(t)) > 1e-3
                viavel = false;
                break;
            end
            
            % Restrições de geradores
            for i = 1:n_geradores
                if potencia_otima(i,t) < geradores(i).pot_min * estado_otimo(i,t) - 1e-3 || ...
                   potencia_otima(i,t) > geradores(i).pot_max * estado_otimo(i,t) + 1e-3
                    viavel = false;
                    break;
                end
            end
            
            if ~viavel, break; end
        end
        
        if viavel
            contagem_viavel = contagem_viavel + 1;
            
            % Custo neste cenário
            custo_cenario = 0;
            for i = 1:n_geradores
                for t = 1:n_horas
                    custo_cenario = custo_cenario + ...
                        geradores(i).coef_a * potencia_otima(i,t)^2 + ...
                        geradores(i).coef_b * potencia_otima(i,t) + ...
                        geradores(i).coef_c * estado_otimo(i,t);
                end
            end
            distribuicao_custos(mc) = custo_cenario;
        end
    end
    
    taxa_viabilidade = contagem_viavel / n_monte_carlo;
    fprintf('Teste de Viabilidade Monte Carlo: %.2f%%\n', taxa_viabilidade * 100);
    
    % Visualização dos resultados
    figure('Position', [100, 100, 1200, 800]);
    
    % Subplot 1: Distribuição de custos
    subplot(2,3,1);
    histograma(distribuicao_custos(distribuicao_custos > 0), 20);
    xlabel('Custo Total ($)');
    ylabel('Frequência');
    title('Distribuição de Custos por Cenário');
    grid on;
    
    % Subplot 2: Comparação robusto vs nominal
    subplot(2,3,2);
    
    % Cálculo para caso nominal
    potencia_nominal = zeros(n_geradores, n_horas);
    for t = 1:n_horas
        carga_liquida = carga_nom(t) - eolica_nom(t) - solar_nom(t);
        restante = carga_liquida;
        for i = 1:n_geradores
            if restante > 0
                potencia_nominal(i,t) = min(restante, geradores(i).pot_max * estado_otimo(i,t));
                restante = restante - potencia_nominal(i,t);
            end
        end
    end
    
    plot(tempo, sum(potencia_otima,1), 'r-', 'LineWidth', 2, 'DisplayName', 'Robusto');
    hold on;
    plot(tempo, sum(potencia_nominal,1), 'b--', 'LineWidth', 2, 'DisplayName', 'Nominal');
    xlabel('Tempo (h)');
    ylabel('Geração Total (MW)');
    title('Despacho Robusto vs Nominal');
    legend('Location', 'best');
    grid on;
    
    % Subplot 3: Análise de reservas
    subplot(2,3,3);
    
    reservas_disp = zeros(1, n_horas);
    reserva_pior_caso = zeros(1, n_horas);
    
    for t = 1:n_horas
        for i = 1:n_geradores
            reservas_disp(t) = reservas_disp(t) + ...
                min(geradores(i).rampa_subida, geradores(i).pot_max*estado_otimo(i,t) - potencia_otima(i,t));
        end
        
        pior_carga = carga_nom(t) * (1 + incerteza_carga);
        reserva_pior_caso(t) = 0.1 * pior_carga;
    end
    
    bar(tempo, [reservas_disp; reserva_pior_caso]');
    xlabel('Tempo (h)');
    ylabel('Reserva (MW)');
    title('Análise de Reservas');
    legend('Disponível', 'Necessária Pior Caso', 'Location', 'best');
    grid on;
    
    % Subplot 4: Impacto em cenários extremos
    subplot(2,3,4);
    
    n_extremos = 5;
    desbalanceamento = zeros(n_extremos, n_horas);
    
    for s = 1:n_extremos
        xi_carga_ext = sign(randn(1, n_horas));
        xi_eolica_ext = sign(randn(1, n_horas));
        xi_solar_ext = sign(randn(1, n_horas));
        
        carga_ext = carga_nom + incerteza_carga * carga_nom .* xi_carga_ext;
        eolica_ext = eolica_nom + incerteza_eolica * eolica_nom .* xi_eolica_ext;
        solar_ext = solar_nom + incerteza_solar * solar_nom .* xi_solar_ext;
        
        for t = 1:n_horas
            geracao_ext = sum(potencia_otima(:,t)) + eolica_ext(t) + solar_ext(t);
            desbalanceamento(s,t) = geracao_ext - carga_ext(t);
        end
    end
    
    plot(tempo, desbalanceamento', 'LineWidth', 1);
    hold on;
    plot(tempo, zeros(1,n_horas), 'k-', 'LineWidth', 2);
    xlabel('Tempo (h)');
    ylabel('Desequilíbrio (MW)');
    title('Balanço em Cenários Extremos');
    grid on;
    
    % Subplot 5: Análise de conservadorismo
    subplot(2,3,5);
    
    custo_nominal = 0;
    for i = 1:n_geradores
        for t = 1:n_horas
            custo_nominal = custo_nominal + ...
                geradores(i).coef_a * potencia_nominal(i,t)^2 + ...
                geradores(i).coef_b * potencia_nominal(i,t) + ...
                geradores(i).coef_c * estado_otimo(i,t);
        end
    end
    
    custo_robusto_medio = mean(distribuicao_custos(distribuicao_custos > 0));
    custo_robusto_max = max(distribuicao_custos);
    
    barras = [custo_nominal, custo_robusto_medio, custo_robusto_max];
    bar(1:3, barras);
    set(gca, 'XTickLabel', {'Nominal', 'Robusto Médio', 'Pior Caso'});
    ylabel('Custo ($)');
    title('Análise de Conservadorismo');
    grid on;
    
    for i = 1:3
        text(i, barras(i), sprintf('%.0f', barras(i)), ...
            'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom');
    end
    
    % Subplot 6: Viabilidade vs orçamento
    subplot(2,3,6);
    
    orcamentos_teste = 0:2:10;
    taxas_viav = zeros(length(orcamentos_teste), 1);
    
    for g = 1:length(orcamentos_teste)
        orc = orcamentos_teste(g);
        cont_viavel = 0;
        
        for mc = 1:100
            xi_teste = 2*rand(1, n_horas) - 1;
            if sum(abs(xi_teste)) > orc
                xi_teste = xi_teste * orc / sum(abs(xi_teste));
            end
            
            carga_teste = carga_nom + incerteza_carga * carga_nom .* xi_teste;
            
            viavel = true;
            for t = 1:n_horas
                if carga_teste(t) > sum(potencia_otima(:,t)) + eolica_nom(t) + solar_nom(t) + 10
                    viavel = false;
                    break;
                end
            end
            
            if viavel
                cont_viavel = cont_viavel + 1;
            end
        end
        
        taxas_viav(g) = cont_viavel / 100;
    end
    
    plot(orcamentos_teste, taxas_viav*100, 'b-o', 'LineWidth', 2);
    xlabel('Orçamento de Incerteza \Gamma');
    ylabel('Taxa de Viabilidade (%)');
    title('Viabilidade vs Nível de Incerteza');
    grid on;
    ylim([0, 105]);
end

5. Recomendações Práticas

5.1 Guia de Seleção de Métodos

5.2 Diretrizes de Parametrização

• Número de cenários: 50-200 cenários geralmente são suficientes; acima disso, os ganhos diminuem
• Orçamento de incertezas (Γ): Configure entre 20-40% do número de variáveis incertas
• Penalidades de corte: Defina como 1,2-1,5 vezes o custo máximo de geração
• Requisitos de reserva: Padrão de 10-15%, aumentando para 15-20% com alta penetração de renováveis

5.3 Estratégias de Eficiência Computacional

1. Computação paralela: Utilize parfor para processar cenários simultaneamente

2. Algoritmos de decomposição: Aplique Benders ou decomposição por cenários

3. Redução de cenários: Mantenha cenários com alta probabilidade e grande diversidade

4. Simplificação de modelos: Utilize linearização ou aproximações lineares por partes

5. Considerações Finais

Este artigo apresentou uma implementação abrangente em MATLAB cobrindo:

1. Modelagem de incertezas: Geração e redução de cenários estocásticos
2. Programação estocástica: Modelo de otimização de dois estágios
3. Otimização robusta: Estrutura adaptativa para o pior caso
4. Aálise visual: Visualização completa de resultados

Aspectos fundamentais a considerar:

• A programação estocástica é mais adequada para cenários onde a econômica é prioritária
• A otimização robusta é preferível quando a segurança é crítica
• A aplicação prática frequentemente combina ambas as abordagens
• A qualidade e quantidade de cenários influencia significativamente os resultados

Direções para trabalhos futuros:

1. Integração de aprendizado de máquina para previsão de incertezas
2. Extensão com restrições de rede elétrica
3. Otimização coordenada em múltiplas escalas temporais
4. Desenvolvimento de algoritmos de otimização distribuída

O código foi desenvolvido com arquitetura modular, facilitando extensões e adaptação para sistemas elétricos reais.