Monitoração da qualidade de energia utilizando estimação ...coral.ufsm.br/sepoc/sepoc2018/arquivos/presentations/SEPOC_S2c.pdfFigura 3: Alocação de PMUs em um sistema de 33 barras

Introdução e Objetivos Fundamentos teóricos Metodologia proposta Testes e Resultados Considerações Finais

Monitoração da qualidade de energia utilizandoestimação de estados e medição fasorial sincronizada

José Luiz Rezende Pereira

23 de outubro de 2018

1 / 31


Sumário

1 Introdução e Objetivos

2 Fundamentos teóricos

3 Metodologia proposta

4 Testes e Resultados

5 Considerações Finais

2 / 31


Introdução

Motivações:

Aumento de distorções harmônicas;

Falta de sistemas de monitoramento para redes dedistribuição;

Averiguação e manutenção de padrões de qualidade deenergia;

Desafio de estimar o estado operativo com número reduzidomedidores;

Utilização de PMUs (Phasor Measurement Units).

3 / 31


Objetivo:

Desenvolvimento de uma metodologia para estimar o estadooperativo da rede tanto para 60Hz como para frequênciasHarmônicas

1 Utiliza-se um número reduzido de PMUs;2 considera-se curvas diárias de carga.

Formulação de problemas de otimização com restrições nãolineares;

Figura 1: Imagem ilustrativa da aplicação do estimador.

4 / 31


Estimação Estática de estados

Conceitos geraisTécnica usada para estimar grandezas elétricas em sistemas depotência;

Baseada na alocação mínima de medidores e conhecimento datopologia;

Processa conjuntamente informações de vários medidores (análisesistêmica).

Metodologia tradicional: Mínimos Quadrados Ponderados (1970)Objetiva minimizar a diferença quadrática entre valores medidos evalores estimados em função das variáveis de estados;

Variáveis de estados: tensões (magnitude e ângulos) em todos osbarramentos;

Consolidada nas redes de transmissão através do sistema SCADA(Supervisory and Control Data Acquisition).

5 / 31


Estimação Estática de estados

Aplicação nas redes de distribuiçãoEscassez de medidores alocados (observabilidade não completa);

Uso de pseudomedidas (dados históricos baseados de demanda)∗;

Uso de Medição Fasorial Sincronizada via GPS (Global PositioningSystem);

Uso de medidores inteligentes (smart meters).

Figura 2: Exemplo de curva diária de carga ∗

6 / 31


Estimação de estados harmônicos

Aplicação nas redes de distribuiçãoObjetivo: Identificação de fontes harmônicas dominantes;

Escassez de medidores instalados;

A literatura destaca uso de PMUs neste contexto;

Desenvolvimento de micro-PMUs com suporte para fasoresharmônicos;

Adoção de métodos meta-heurísticos;Tempo computacional desvantajoso;Técnicas não determinísticas;Confiabilidade necessária inerente ao processo de estimação de estados.

O tema é atual e desafiador, sendo um campo fértil de pesquisa.

7 / 31


Metodologia proposta

Considerações iniciais

Formulação de problemas de otimização;

Cada ordem harmônica h é estimada individualmente;

Execução em intervalos de tempo regulares t.

Procedimento de alocação de PMUs

Figura 3: Alocação de PMUs em um sistema de 33 barras.

Preconiza 2 PMUs por ramais do alimentador (principal e laterais);

Outros procedimentos são testados durante o trabalho.

8 / 31



Vetor de medidas

zht = [zh

(1,re), zh(1,im), z

h(2,re), z

h(2,im), ..., z

h(Nm,re), z

h(Nm,im)]

Tt (1)

Em que um fasor (zh(re)+jzh

(im)) é dividido em partes real e imaginária;

2 Nm é o número total de elementos dentro do vetor.

Vetor de variáveis de estados

xht = [IA,h

(1,re), IA,h(1,im), I

A,h(2,re), I

A,h(2,im), ..., I

A,h(Nl,re), I

A,h(Nl,im),

IB,h(1,re), I

B,h(1,im), I

B,h(2,re), I

B,h(2,im), ..., I

B,h(Nl,re), I

B,h(Nl,im),

IC,h(1,re), I

C,h(1,im), I

C,h(2,re), I

C,h(2,im), ..., I

C,h(Nl,re), I

C,h(Nl,im)]

Tt

(2)

(Is,h(b,re) + jIs,h

(b,im)): fasor de corrente ramal;

Considerando as três fases, o número total de elementos do vetor será (2 × 3 × Nl).

9 / 31



Grandezas elétricas calculadas em função das variáveis

Is,hk,re =

∑b∈Ωk

Is,hb,re (3)

Is,hk,im =

∑b∈Ωk

Is,hb,im (4)

[VABC,hk ] = [YABC,h

bar ]−1[IABC,hk ] (5)

Ps,hk = Vs,h

(k,re)Is,h(k,re) + Vs,h

(k,im)Is,h(k,im)

Qs,hk = −Vs,h

(k,re)Is,h(k,im) + Vs,h

(k,im)Is,h(k,re)

(6)

s ∈ A,B,C.10 / 31



Formulação:

min Jht (xh

t ) =12

2Nm∑j=1

zhj,t − φj,t(xh

t )

σj

2

(7)

sujeito a:

para h = 1 :(1 − ps,1

k,t )Ps,1k,(t−1) ≤ Ps,1

k,t ≤ (1 + ps,1k,t )Ps,1

k,(t−1)(1 − qs,1

k,t )Qs,1k,(t−1) ≤ Qs,1

k,t ≤ (1 + qs,1k,t )Qs,1

k,(t−1)

(8)

para h , 1 :−(1 + ps,h

k,t )Ps,hk,(t−1) ≤ Ps,h

k,t ≤ (1 + ps,hk,t )Ps,h

k,(t−1)−(1 + qs,h

k,t )Qs,hk,(t−1) ≤ Qs,h

k,t ≤ (1 + qs,hk,t )Qs,h

k,(t−1)

(9)

Restrições formadas apenas para as barras k não monitoradas;

Usadas para suprir a escassez de medições (informaçõessuplementares acerca do estado operativo da rede).

11 / 31



A função objetivo

min Jht (xh

t ) =12

2Nm∑j=1

zhj,t − φj,t(xh

t )

σj

2

(10)

zht : valor da medição;

φht (xh

t ): valor estimado correspondente;

σ: desvio padrão associado à classe do instrumento de medição(ex.: σ = 0.5%);

Jht : valor da função objetivo.

12 / 31


Restrições para h = 1

(1 − ps,1k,t )P

s,1k,(t−1) ≤ Ps,1

k,t ≤ (1 + ps,1k,t )P

s,1k,(t−1)

(1 − qs,1k,t )Q

s,1k,(t−1) ≤ Qs,1

k,t ≤ (1 + qs,1k,t )Q

s,1k,(t−1)

(11)

Ps,1k,t , Qs,1

k,t : potência ativa/reativa a ser determinada em um instantede tempo t;

Ps,1k,(t−1), Qs,1

k,(t−1): valores estimados em um intervalo de tempoanterior (t − 1);

ps,1k,t e qs,1

k,t : fatores percentuais de incerteza relativos à variaçãoentre instantes de tempo consecutivos.

Estado em t=1 Estado em t=2 Estado em t=3Informações

históricas ...

t-1 t

13 / 31


Restrições para h , 1

para h , 1 :−(1 + ps,h

k,t )Ps,hk,(t−1) ≤ Ps,h

k,t ≤ (1 + ps,hk,t )Ps,h

k,(t−1)−(1 + qs,h

k,t )Qs,hk,(t−1) ≤ Qs,h

k,t ≤ (1 + qs,hk,t )Qs,h

k,(t−1)

(12)

Ps,hk,t , Qs,h

k,t podem assumir valores positivos ou negativos (entreguesda fonte para a rede ou da rede para o consumidor).

Figura 4: Caso típico de carga não linear.

14 / 31


Resolução por Método de Pontos Interiores com barreira de segurança

min Jht (xh

t )sujeito a (s.a.):l ≤ g(xh

t ) ≤ u(13)

min Jht (xh

t ) − µn∑

j=1ln(sl,j) − µ

n∑j=1

ln(su,j)

s.a.:g(xh

t ) − sl = lg(xh

t ) + su = usl > 0, su > 0

(14)

min Jht (xh

t ) − µn∑

j=1ln(sl,j + ρ) − µ

n∑j=1

ln(su,j + ρ)

s.a.:g(xh

t ) − sl = lg(xh

t ) + su = usl ≥ 0, su ≥ 0, ρ > 0

(15)

l, u: valores limítrofes mínimos e máximos das restrições g;

sl , su : variáveis de folga associados aos limites mínimos e máximos;

ρ: parâmetro barreira de segurança. Inicialmente pequeno (0,01) e decai exponencialmente até um valor mínimoρmin = 10−8 .

15 / 31


Resolução por MPIB modificado

Uma vez formulado o problema, aplica-se as condições de optimalidade de primeiraordem e resolve-se o problema via Newton-Raphson;

Como diferencial, pode-se detectar as restrições ativas que ‘encostam’ nos seusrespectivos limites (variáveis de folga serão nulas);

Para garantir que a solução encontrada esteja dentro de limites adequados, essespoderão ser relaxados a cada iteração:

uj,it+1 = uj,it(1 + ζu,j)lj,it+1 = lj,it(1 + ζl,j)

(16)

ζu,j = ρ e

(−

su,jρmin

)ζl,j = ρ e

(−

sl,jρmin

) (17)

ρmin = 10−8 ;

Caso a variável de folga seja zero (restrição encosta no limite), o parâmetro de relaxamento será igual a ρ;

Caso contrário, será praticamente nulo;

A solução é obtida dentro de limites propostos e atualizáveis de acordo com a necessidade de determinar oestado operativo mais provável do sistema;

Critério de convergência: |∆xht | < 10−6 e niter = 300.

16 / 31


Selecionar ordem

harmônica

t =1

Fluxo de potência

harmônico

Obter medições fasoriais

para ordem harmônica h

para o intervalo de tempo t

Variação de

carga

Obter solução

ótima

t =t+1

Bloco-1

Bloco-2

Bloco-3

Bloco-4

Bloco-5

Bloco-6

Bloco-7

17 / 31


Processamento Sequencial e Paralelo

Estima ordem

harmônica h=1

Estima ordem

harmônica h=2

Estima ordem

harmônica h=3

Estima ordem

harmônica h=4Início

Início

Fim

Estima ordem

harmônica

h=1

Estima ordem

harmônica

h=2

Estima ordem

harmônica

h=3

Estima ordem

harmônica

h=4

Fim

SEQUENCIAL

PARALELO

Estima ordem

harmônica h=1

Estima ordem

harmônica h=2

Estima ordem

harmônica h=3

Estima ordem

harmônica h=4Início

Início

Fim

Estima ordem

harmônica

h=1

Estima ordem

harmônica

h=2

Estima ordem

harmônica

h=3

Estima ordem

harmônica

h=4

Fim

SEQUENCIAL

PARALELO

18 / 31


Figura 5: Sistema de 33 barras.

Fontes harmônicas (FH) inseridas para simular a resposta do sistema;

Modelagem seguindo tutorial da PES (Power Energy Society);

Os parâmetros ps,hk,t e qs,h

k,t serão adotados como 10%;

Dados históricos de demanda são obtidos por valores tabelados no apêndice datese (Dados de carregamento do sistema).

19 / 31


Figura 6: Resultado de estimação para a barra 11, fase A, frequência fundamental.

0 4 8 12 16 20 24350

400

450

500

Tempo (horas)

Pot

ênci

a A

tiva

(KW

)

Estimado

Verdadeiro

(a) Potência ativa

0 4 8 12 16 20 240

0,1

0,2

0,3

0,4

0,5

Tempo (horas)

Err

o (%

)

Erro de estimação

(b) Erros

0 4 8 12 16 20 24

0,98

0,985

0,99

0,995

1

Tempo(horas)

Mag

nitu

de (

p.u.

)

EstimadoVerdadeiro

(c) Magnitude de tensão

0 4 8 12 16 20 240

0,05

0,1

0,15

0,2

Tempo (horas)

Err

o (%

)

Erro de estimação

(d) Erros20 / 31


Figura 7: Resultado de estimação para a barra 11, fase A, h = 5.

0 4 8 12 16 20 24

−0,56

−0,54

−0,52

−0,5

−0,48

Tempo (horas)

Pot

ênci

a A

tiva

(KW

)

Estimado

Verdadeiro

(a) Potência ativa

0 4 8 12 16 20 240

0,05

0,1

0,15

0,2

Tempo (horas)

Err

o (%

)

Erro de estimação

(b) Erros

0 4 8 12 16 20 240,0068

0,007

0,0072

0,0074

0,0076

0,0078

Tempo(horas)

Mag

nitu

de (

p.u.

)

EstimadoVerdadeiro

(c) Magnitude de tensão

0 4 8 12 16 20 240

0,02

0,04

0,06

0,08

0,1

Tempo (horas)

Err

o (%

)

Erro de estimação

(d) Erros21 / 31


Figura 8: Magnitude de tensões para todas as barras em ∆t = 24:00,frequência fundamental, fase A.

5 10 15 20 25 300,980

0,985

0,990

0,995

1,000

1,005

Barras

Tens

ões

(p.u

.)

VerdadeiroEstimado

(a) Real/ Estimado

5 10 15 20 25 300,000

1,000

2,000

3,000 ·10−3

BarrasE

rros

(%)

Erro percentual

(b) Erro Percentual

22 / 31


Figura 9: Resultado de estimação de tensões harmônicas em 24:00 (fase A).

3

5

7

9

11

13

15

14

711

141719

2225

2831

0

0,01

0,02

HarmônicaBarras

Mag

nitu

de (

p.u.

)

0

0,002

0,004

0,006

0,008

0,01

0,012

0,014

0,016

(a) Verdadeiro

3

5

7

9

11

13

15

14

711

141719

2225

2831

0

0,01

0,02

HarmônicaBarras

Mag

nitu

de (

p.u.

)

0

0,002

0,004

0,006

0,008

0,01

0,012

0,014

0,016

(b) Estimado

35

79

1113

15

1 4 7 11 14 1719 22 25 28 31

0

0,0005

0,001

0,0015

Harmônica

Barra

Err

os (

%)

0

0,0002

0,0004

0,0006

0,0008

0,001

(c) Erros 23 / 31


Tabela 1: EQMN (%) máximos para a fase A.

h 1 3 5 7 9 11 13 15

Vmag 0,002 0,031 0,035 0,030 0,031 0,032 0,030 0,033

Vang 0,022 0,026 0,022 0,021 0,029 0,021 0,021 0,025

Imag 0,001 0,011 0,008 0,008 0,011 0,009 0,008 0,010

Iang 0,041 0,008 0,007 0,008 0,010 0,009 0,010 0,011

Erros bem reduzidos;

Erros maiores para as tensões, já que o vetor de variáveis é formado por correntes;

Análise não trivial: depende de diversos fatores:Número de PMUs;Dimensão da rede;Classe de exatidão dos medidores.

24 / 31


Procedimentos de alocação

Figura 10: Sistema de 33 barras com 16 PMUs (usando Programação Linear Inteira).

Tabela 2: EQMN (%) máximos para a fase A (16 PMUs).

h 1 3 5 7 9 11 13 15

Vmag 1,21x10−4 0,010 0,008 0,007 0,009 0,010 0,009 0,010

Vang 1,11x10−4 0,008 0.007 0,007 0,008 0,005 0,007 0,006

Imag 0,91x10−4 0,001 0,002 0,001 0,003 0,001 0,001 0,001

Iang 1,78x10−4 1,77x10−4 1,52x10−4 1,79x10−4 0,97x10−4 2,01x10−4 1,45x10−4 1,91x10−4

25 / 31


Figura 11: Sistema de 69 barras.

0 4 8 12 16 20 240

0,2

0,4

0,6

0,8

1

Tempo (horas)

fato

r m

ultip

licat

ivo

Geração distribuída com capacidade de gerar 80kW na barra 16 (com demanda de60kW) e um filtro passivo sintonizado em 180Hz na barra 5.

Apesar de haver FH dominantes, há fontes dispersas menores na rede;

Intuito: testar robustez do método em um cenário mais complexo.

26 / 31


Figura 12: Resultado de estimação para a barra 16, h = 1.

0 4 8 12 16 20 24−20

0

20

40

60

80

Tempo (horas)

Pot

ênci

a A

tiva

(KW

)Estimado

Verdadeiro

(a) Potência ativa, h = 1

0 4 8 12 16 20 240

0,5

1

1,5

Tempo (horas)

Err

o (%

)

Erro de estimação

(b) Erros

0 4 8 12 16 20 24

0,982

0,983

0,984

0,985

Tempo(horas)

Mag

nitu

de (

p.u.

)

0,981

Estimado

Verdadeiro

(a) Magnitude de tensão, h = 1

0 4 8 12 16 20 240

0,02

0,04

0,06

0,08

0,1

Tempo (horas)

Err

o (%

)

Erro de estimação

(b) Erros

27 / 31


Figura 13: Resultado de estimação para a barra 20.

0 4 8 12 16 20 240,016

0,0165

0,017

0,0175

0,018

0,0185

Tempo(horas)

Mag

nitu

de (

p.u.

)

Estimado

Verdadeiro

(a) Magnitude da tensão, h = 3

0 4 8 12 16 20 240

0,05

0,1

0,15

0,2

Tempo (horas)

Err

o (%

)

Erro de estimação

(b) Erros

0 4 8 12 16 20 24

6,4

6,8

7,2

7,6

Tempo (horas)

THD(%

)

Estimado

Verdadeiro

(a)THD, barra 20

0 4 8 12 16 20 240

0,1

0,2

0,3

0,4

0,5

Tempo (horas)

Err

o (%

)

Erro de estimação

(b) Erros

28 / 31


Tabela 3: EQMN (%) máximos para a fase A (ps,hk,t e qs,h

k,t iguais a 10%).

h 1 3 5 7 9 11 13 15

Vmag 0,010 0,032 0,033 0,033 0,032 0,038 0,035 0,034

Vang 0,035 0,008 0,002 0,009 0,009 0,007 0,004 0,007

Imag 0,001 0,012 0,009 0,005 0,011 0,007 0,006 0,011

Iang 0,041 0,008 0,007 0,008 0,011 0,007 0,011 0,013

Figura 14: Maiores erros para diferentes valores dos fatores ps,hk,t e qs,h

k,t .

0 10 20 30 40 500

0,05

0,1

0,15

0,2

0,25

0,3

0,35

0,4

EQ

MN

max

(%)

ps,hk,t

e qs,hk,t

(%)29 / 31


Tabela 4: Comparação entre métodos publicados.

Pesquisa Técnica Ano Barras Medidores Razão (%) Erro(%) PMU

Arefi HBM 2009 70 (mon.) 8 11,42 2,133 Sim

Arruda AG 2010 50 (tri.) 8 16 ≤5 Não

Arefi PSO 2011 70 (mon.) 6 8,57 2,103 Sim

Moradifar Fuzzy+MQP 2017 69 (mon.) 24 34,78 ≤ 2 Sim

Bevcirovic AG 2018 7 (mon.) 4 57,14 1,129 Sim

PROPOSTO MPIBS 2018 69 (trif.) 16 23,18 ≤ 1 Sim

HBM: Honey Bee Mating;

AG: Algoritmo Genético;

PSO: Particle Swarm Optimization;

MQP: Minimos quadrados ponderados;

MPIBS: Método de pontos interiores com barreira de segurança.

30 / 31


Considerações finais

ConclusõesEstados harmônicos e curvas diárias são estimados de maneirasatisfatória;

Erros de estimação são menores que 1%;

Valida-se o uso do método de pontos interiores com barreira desegurança;

Análises de sensibilidade mostram robustez do algoritmo;

Tempo computacional vantajoso usando processamento paralelo.

31 / 31

Documents

Monitoração da qualidade de energia utilizando estimação ...coral.ufsm.br/sepoc/sepoc2018/arquivos/presentations/SEPOC_S2c.pdfFigura 3: Alocação de PMUs em um sistema de 33 barras