1. Contexto do Problema: Desafios da Incerteza de Preço na Autoprogramação de Energia

A autoprogramação refere-se ao processo pelo qual empresas de geração de energia ou operadores de microrredes formulam plenos de geração e estratégias de lances no mercado de eletricidade para maximizar lucros ou minimizar custos. Os principais desafios incluem:

• Incerteza de Preço: Os preços do mercado de energia no dia anterior são altmaente voláteis devido a flutuações na demanda, geração de energia renovável e congestionamento da rede.
• Risco de Cauda: Cenários de preços extremamente baixos podem levar a perdas financeiras substanciais para os geradores. A otimização estocástica tradicional frequentemente negligencia esses riscos de cauda.
• Ambiguidade de Distribuição: A distribuição de probabilidade real dos preços de eletricidade é desconhecida, e apenas informações de momentos (como média e covariância) podem ser estimadas a partir de dados históricos, sujeitas a erros de estimação.

Limitações das Abordagens Tradicionais:

• Otimização Estocástica: Depende de distribuições de probabilidade precisas e é sensível a erros na especificação da distribuição.
• Otimização Robusta: Tende a ser excessivamente conservadora, considerando apenas os cenários piores.
• Aplicação Isolada de CVaR: Requer uma forma de distribuição pré-definida e não consegue lidar com a ambiguidade da distribuição.

2. Modelo Central: Estrutura Híbrida DRO Baseada em Momentos e CVaR

2.1 Princípios Fundamentais da Otimização Robusta de Distribuição (DRO)

• Construção do Conjunto de Incerteza: Com base nas estimativas da média (μ) e da matriz de covariância (Σ) de preços históricos de eletricidade, define-se um conjunto de incerteza de momentos: $$ \mathcal{P} = \left\{ P \in \mathcal{D} \mid E[\boldsymbol{\xi}] = \boldsymbol{\mu}, \text{Cov}[\boldsymbol{\xi}] = \boldsymbol{\Sigma} \right\} $$

  onde $$ \boldsymbol{\xi} $$ é o vetor de preços de eletricidade aleatórios e $$ \mathcal{D} $$ é o conjunto de todas as distribuições de probabilidade admissíveis.

• Função Objetivo: Otimiza o custo esperado sob o pior cenário dentro do conjunto de incerteza: $$ \min_{\mathbf{x} \in \mathcal{X}} \max_{P \in \mathcal{P}} E_{P}[\text{custo}(\mathbf{x}, \boldsymbol{\xi})] $$

  onde $$ \mathbf{x} $$ representa as decisões de programação (por exemplo, status de desligamento/ligamento da unidade, geração de energia).

2.2 Incorporação do Valor Condicional em Risco (CVaR)

• Definição: Em um nível de confiança $$ \beta $$, o CVaR representa a perda esperada excedendo o Valor em Risco (VaR): $$ \text{CVaR}_{\beta}(L) = \min_{\eta \in \mathbb{R}} \left\{ \eta + \frac{1}{1-\beta} E[\max(0, L - \eta)] \right\} $$

  onde $$ L $$ é a função de perda (por exemplo, lucro negativo).

• Combinação com DRO: A incorporação do CVaR no objetivo ou nas restrições permite otimizar o pior CVaR sob o conjunto de incerteza $$ \mathcal{P} $$: $$ \min_{\mathbf{x} \in \mathcal{X}} \max_{P \in \mathcal{P}} \text{CVaR}_{\beta}(\text{perda}(\mathbf{x}, \boldsymbol{\xi})) $$

  Esta formulação aborda simultaneamente a ambiguidade da distribuição e o risco de cauda.

3. Resolução do Modelo: Redução da Dimensão de Otimização (ODR) e Algoritmo ADMM

3.1 Gargalo Computacional e a Abordagem ODR

• Problema: Modelos DRO baseados em momentos frequentemente levam a Problemas de Programação de Semicontinuidade (SDP), cujas dimensões crescem exponencialmente com o número de variáveis aleatórias.
• Ideia Central do ODR: Capitaliza a característica de baixo posto dos parâmetros aleatórios para decompor o SDP de alta dimensão em subproblemas de baixa dimensão.

  1. Demonstra-se que o posto superior da matriz do problema original é $$ r \ll \text{dim}(\boldsymbol{\xi}) $$.

  2. A dimensão é reduzida para $$ r $$, construindo um SDP aproximado de duas camadas (interno e externo):

     • A camada externa fornece um limite inferior (solução conservadora).
     • A camada interna fornece um limite superior (solução agressiva).

     Essas camadas convergem iterativamente para a solução global ótima.

3.2 Fluxo de Solução ADMM

Um algoritmo ADMM aprimorado é aplicado ao SDP não convexo de dimensão reduzida:

1. Separação de Variáveis: Desacopla as variáveis de decisão $$ \mathbf{x} $$ e as variáveis auxiliares $$ \mathbf{z} $$.

2. Atualização Alternada: $$ \begin{aligned} \mathbf{x}^{k+1} &= \arg \min_{\mathbf{x}} \mathcal{L}_{\rho}(\mathbf{x}, \mathbf{z}^k, u^k) \\ \mathbf{z}^{k+1} &= \arg \min_{\mathbf{z}} \mathcal{L}_{\rho}(\mathbf{x}^{k+1}, \mathbf{z}, u^k) \\ u^{k+1} &= u^k + \rho(\mathbf{x}^{k+1} - \mathbf{z}^{k+1}) \end{aligned} $$

   onde $$ \mathcal{L}_{\rho} $$ é a função Lagrangiana aumentada.

3. Convergência: Sob as restrições de momentos, o algoritmo converge para uma solução $$ \epsilon $$-ótima.

Vantagens de Desempenho:

• Redução do tempo de computação em ordens de magnitude (por exemplo, de horas para segundos).
• Diferença entre a solução obtida e a ótima global < 0.1%.

4. Modelagem da Incerteza de Preço de Eletricidade e Geração de Cenários

4.1 Modelagem da Distribuição de Probabilidade

• Paramétrica: Assume que os preços de eletricidade seguem uma distribuição t-location-scale (parâmetro de localização $$ \mu_t $$, parâmetro de escala $$ \sigma_t $$), adequada para características de cauda pesada.
• Não Paramétrica: Utiliza a estimativa de densidade de kernel (KDE) para ajustar os dados históricos, evitando desvios devido a suposições de distribuição.

4.2 Geração e Redução de Cenários

1. Amostragem de Monte Carlo: Gera $$ N $$ cenários de preço de eletricidade.

2. Redução de Kantorovich: Seleciona $$ S $$ cenários representativos (por exemplo, $$ S=30 $$), minimizando a distância para o conjunto de cenários original: $$ \min_{ \{ \hat{P}_s \}_{s=1}^S } \sum_{s=1}^S \min_{P \in \mathcal{P}} d(\hat{P}_s, P) $$

   Garantindo a eficiência computacional.

5. Estudo de Caso: Sistemas de Barramentos IEEE

5.1 Sistema IEEE de 118 Barramentos

• Estrutura: 54 geradores, 186 linhas, 99 barramentos de carga, operando em múltiplas áreas.
• Problema de Autoprogramação:

  • Variáveis de Decisão: Status de ligar/desligar $$ u_{ti} $$, geração $$ p_{ti} $$, capacidade de reserva $$ r_{ti} $$ da unidade $$ i $$ no tempo $$ t $$.

  • Função Objetivo (DRO-CVaR Combinado): $$ \min_{\mathbf{x}} \max_{P \in \mathcal{P}} \left( \sum_{t=1}^T \sum_{i \in G} c_i(p_{ti}) - \lambda_t p_{ti} \right) $$

    onde $$ c_i(p_{ti}) $$ é o custo de combustível, $$ \lambda_t $$ é o preço de eletricidade.

  • Restrições: Balanço de potência, taxas de rampa, requisitos de reserva.

5.2 Resultados Chave

• Controle de Risco: Para $$ \beta=0.95 $$, as perdas de cauda são reduzidas em 23% (em comparação com a otimização estocástica).
• Eficiência Computacional: O tempo de solução ODR-ADMM é < 10 segundos, enquanto o SDP tradicional excede 1 hora.
• Melhoria na Conservação: O custo esperado do conjunto de incerteza de momentos é 12% menor do que o da otimização robusta.

6. Comparação e Extensões

6.1 Comparação com Outras Abordagens

6.2 Extensões para Aplicações Industriais

• Licitação de Microrredes: Incorporando a incerteza de energia solar/eólica.
• Restrições de Emissões de Carbono: Introduzindo funções de custo de carbono.
• Programação Multiagente: Coordenação de múltiplos geradores usando ADMM distribuído.

Conclusão

O modelo DRO-CVaR baseado em momentos proposto neste trabalho utiliza conjuntos de incerteza para caracterizar a ambiguidade da distribuição de preços de eletricidade e CVaR para quantificar o risco de cauda. A combinação com o algoritmo ODR-ADMM garante uma solução eficiente. As validações em sistemas IEEE (incluindo o de 118 barramentos) demonstram que:

1. Em níveis de risco equivalentes, o custo esperado é reduzido em 8%-15%.
2. As perdas em cenários extremos são reduzidas em mais de 20%.
3. A eficiência computacional atende aos requisitos de tempo do mercado do dia anterior (solução em nível de segundos).

Direções futuras incluem a integração de conjuntos de incerteza de Wasserstein orientados por dados, a fusão com aprendizado profundo para previsão de parâmetros de momentos em tempo real e a extensão para sistemas de acoplamento multie-energia.

Código MATLAB de Exemplo

T = SCUC_data.totalLoad.T;  % Número de períodos
G = SCUC_data.units.N;      % Número de geradores
N = SCUC_data.baseparameters.busN;  % Número total de barramentos

all_branch.I = [ SCUC_data.branch.I; SCUC_data.branchTransformer.I ]; % Origens de todas as linhas
all_branch.J = [ SCUC_data.branch.J; SCUC_data.branchTransformer.J ]; % Destinos de todas as linhas
all_branch.P = [ SCUC_data.branch.P; SCUC_data.branchTransformer.P ]; % Limites de potência das linhas

beta_CVaR = 0.95;

% Formar a matriz de coeficientes de fluxo de potência DC
type_of_pf = 'DC';
Y = SCUC_nodeY(SCUC_data,type_of_pf);
B = -Y.B; % Para equações DC, B ignora resistência, considerando apenas reatância

% Definir constantes úteis
diag_E_T = sparse(1:T,1:T,1); % Matriz diagonal T*T com 1s na diagonal
low_triangle = diag_E_T(2:T,:) - diag_E_T(1:T-1,:);

% Definir variáveis y e z
for i=1:N    % Processar cada barramento individualmente
   if ismember(i, SCUC_data.units.bus_G)    % Se for um barramento de gerador, as variáveis de decisão incluem potência, ângulo, restrição de potência z
       y{i}.PG = sdpvar(T,1); % Criar variáveis de decisão de número real sdpvar
       y{i}.theta = sdpvar(T,1);
       y{i}.z = sdpvar(T,1);
   else
       y{i}.theta = sdpvar(T,1);
   end
end

for i=1:N    % Processar cada barramento individualmente
   if ismember(i, SCUC_data.units.bus_G)    % Se for um barramento de gerador, as variáveis de decisão incluem potência, ângulo, restrição de potência z
       z{i}.PG = sdpvar(T,1); % Criar variáveis de decisão de número real sdpvar
       z{i}.theta = sdpvar(T,1);
       z{i}.z = sdpvar(T,1);
   else
       z{i}.theta = sdpvar(T,1);
   end
end
% Fim da definição de variáveis

% Calcular estimativas de média miu_hat e covariância sigma_hat
miu_hat_G_T = zeros(G,T); % Preços de eletricidade para cada gerador em cada período são diferentes, G*T
for q = 1:q_line % q_line marca o comprimento da amostra
   miu_hat_G_T = miu_hat_G_T + reshape(lamda_q_NT(q,:,:),G,T); % Acumular preços de eletricidade do tipo matriz G*T por amostra q
end
miu_hat_G_T = 1/q_line * miu_hat_G_T;
miu_hat = reshape(miu_hat_G_T',G*T,1); % Empilhar em um vetor de coluna para todos os períodos de cada gerador

Sigma_hat = zeros(G*T,G*T);
for q = 1:q_line
   tmp1 = reshape(lamda_q_NT(q,:,:),G,T) - miu_hat_G_T;
   tmp2 = reshape(tmp1', G*T,1);
   Sigma_hat = Sigma_hat + tmp2  * tmp2';
end
Sigma_hat = 1/q_line * Sigma_hat;
Sigma_hat_neg_half = Sigma_hat^(-1/2);
% Fim do cálculo de média e covariância

% Construir o conjunto de incerteza S (vetor coluna G*T, empilhado por todos os períodos de cada G)
tmp = max(lamda_q_NT,[],1); % Pegar o valor máximo para cada período entre os geradores G de todas as amostras
tmp = reshape(tmp,G,T)';
lamda_positive = reshape(tmp, G*T,1);  % Calcular lamda positivo
tmp = min(lamda_q_NT,[],1); % Pegar o valor mínimo para cada período entre os geradores G de todas as amostras
tmp = reshape(tmp,G,T)';
lamda_negative = reshape(tmp, G*T,1);  % Calcular lamda negativo
% Fim da construção do conjunto de incerteza

% Chamar a função de partição para dividir as regiões
[PI,PINumber,PIG] = portion(N);
D = length(PINumber); % Determinar o número de partições

% Distinguir entre nós internos e externos
ext = [];
for d = 1:D
   for i = 1:size(all_branch.I,1) % Iterar por todas as linhas a partir da primeira
       left = all_branch.I(i); % Origem e destino da linha determinam a susceptância B
       right = all_branch.J(i);
       if ismember(left, PI{d}') && ~ismember(right, PI{d}') % Verificar se a linha pertence à partição atual
           ext = [ext;left;right];
       end
   end
end
ext = unique(ext);
% Definir outras variáveis
r = sdpvar(D,1);
t = sdpvar(D,1);
alpha_CVaR = sdpvar(D,1);
z_vector = [];
P_vector = [];
for g = 1:G
   z_vector = [z_vector;  y{SCUC_data.units.bus_G(g)}.z  ];
   P_vector = [P_vector;  y{SCUC_data.units.bus_G(g)}.PG ]; % Separar vetores coluna com ;
end

% Configurar parâmetros ADMM
u_A_b = [];
rho = 5;
gap_pri = 1e-2;
gap_dual = 1e-2;

M = 100; % Número de iterações
% core = D;
% p = parpool(core);
% p.IdleTimeout = 100;
for k = 1:M
   p_A = [];
   z_A = [];
   z_A_pri = [];
   u_A = [];
   p_A_k1 = [];
   u_r = 0; % Marcador para o comprimento das variáveis de partição
   % Atualização p
   for d = 1:D
       PIi = PI{d};
       PIN = PINumber{d};
       d_g = PIG{d};
       Constraints = []; % Limpar o conjunto de restrições a cada iteração
       miu_hat_d = [];
       lamda_positive_d = [];
       lamda_negative_d = [];
       y_d = [];
       zy_d = [];
       theta_d = [];
       ztheta_d = [];
       z_d = [];
       for i = PIi' % Iterar pelos nós na partição atual
           theta_d = [theta_d;y{i}.theta]; % Obter os ângulos de todos os barramentos na partição atual
           ztheta_d = [ztheta_d;z{i}.theta]; % Obter os ângulos de todos os barramentos na partição atual

           if ismember(i, SCUC_data.units.bus_G)
               y_d = [y_d;y{i}.PG]; % Obter a geração de todos os geradores na partição atual
               zy_d = [zy_d;z{i}.PG]; % Obter a geração de todos os geradores na partição atual
               z_d = [z_d;y{i}.z]; % Obter o custo do gerador na partição atual
               i_g = find(SCUC_data.units.bus_G == i); % Encontrar o índice do gerador atual
               % Limites superior e inferior de geração da unidade
               Constraints = [Constraints, ...
                   SCUC_data.units.PG_low(i_g) * ones(T,1) <= y{i}.PG <= SCUC_data.units.PG_up(i_g) * ones(T,1) ];

               % Taxa de rampa Pup e Pdown são iguais e P0 não é considerado
               Constraints = [Constraints, ...
                   -SCUC_data.units.ramp(i_g) * ones(T-1,1) <= low_triangle * y{i}.PG <= SCUC_data.units.ramp(i_g) * ones(T-1,1) ];

               % Linearização do custo de geração
               for l = 0:(L-1)
                   p_i_l =SCUC_data.units.PG_low(i_g) +  ( SCUC_data.units.PG_up(i_g) - SCUC_data.units.PG_low(i_g) ) / L * l;
                   Constraints = [Constraints, ...
                       (2* p_i_l * SCUC_data.units.gamma(i_g) +SCUC_data.units.beta(i_g) ) * y{i}.PG - y{i}.z <= (p_i_l^2 * SCUC_data.units.gamma(i_g) - SCUC_data.units.alpha(i_g)) * ones(T,1) ];
               end

               % Calcular a média de preço de eletricidade da partição
               miu_hat_m = zeros(1,T);
               for q = 1:q_line % q_line marca o comprimento da amostra
                   miu_hat_m = miu_hat_m + reshape(lamda_q_NT(q,i_g,:),1,T); % Acumular preços de eletricidade do gerador i por amostra q
               end
               miu_hat_m = 1/q_line * miu_hat_m;
               miu_hat_m = reshape(miu_hat_m',T,1); % Empilhar em vetor de coluna para todos os períodos do gerador
               miu_hat_d = [miu_hat_d ; miu_hat_m];

               % Obter valores máximos e mínimos correspondentes para o gerador
               lamda_positive_m = lamda_positive(1+T*(i_g-1):T*i_g,1);
               lamda_positive_d = [lamda_positive_d;lamda_positive_m];
               lamda_negative_m = lamda_negative(1+T*(i_g-1):T*i_g,1);
               lamda_negative_d = [lamda_negative_d;lamda_negative_m];
           end

           if ~ismember(i, ext) % Verificar se o barramento é um nó interno
               % Construir a desigualdade de fluxo de potência DC interna da partição
               % Construir termo intermediário
               cons_PF = sparse(T,1); % Uma restrição por T
               if ismember(i, SCUC_data.units.bus_G) % Se i for um barramento de gerador, o termo intermediário deve incluir PG
                   cons_PF = cons_PF +  y{i}.PG;
               end
               for j = 1:N
                   cons_PF = cons_PF -B(i,j) .* y{j}.theta;
               end

               % Construir o lado direito
               % find retorna o número do índice
               index_loadNode = find(SCUC_data.busLoad.bus_PDQR==i); % O barramento i é um nó de carga?
               if index_loadNode>0
                   b_tmp = SCUC_data.busLoad.node_P(:,index_loadNode); % Puxar a carga do nó de carga pelo índice
               else
                   b_tmp = sparse(T,1);
               end
               Constraints = [Constraints, ...
                   0 <= cons_PF <= b_tmp ];
               % Fim da construção da desigualdade de fluxo de potência da partição interna
           end
   

Referências

Algumas partes deste artigo foram extraídas da internet e as fontes são indicadas ou referenciadas. Pode haver omissões, e quaisquer imprecisões podem ser corrigidas mediante contato. (O conteúdo do artigo é apenas para referência; os resultados de execução são o critério final.)

[1] R. A. Jabr, "Robust self-scheduling under price uncertainty using conditional value-at-risk,"
IEEE Trans. Power Syst., vol. 20, no. 4, pp. 1852-1858, Nov. 2005.