📊 Aula 05 — Algoritmos Numéricos em C e Python: Fasores, Matrizes e Solução de Circuitos

Navegação Rápida

📍 Disciplina: 03 Programação Aplicada a Sistemas Elétricos | 🗺️ Área: 06 Computação e Telecomunicações | 🏠 Hub Principal ⬅️ Aula Anterior: Aula 04 - Memory-Mapped I/O e ISR | 🔝 MOC da Disciplina: Hub 03 Programação Aplicada

---

📌 Sumário da Aula

1. 🧠 1. Paradigma Computacional Híbrido: C para Tempo Real e Python para Modelagem
2. ⚡ 2. Extração Fasorial Discreta de 1 Ciclo (Filtro Cosseno/Seno)
3. 📐 3. Grandezas de Potência em Tempo Real e Aritmética Complexa em C99
4. 🌐 4. Montagem Algorítmica da Matriz Ybus e Armazenamento Esparso
5. 📈 5. Solução de Fluxo de Potência por Newton-Raphson com Fatoração LU
6. 📊 6. Diagrama de Algoritmos Numéricos, Fasores e Solução de Redes
7. 📝 7. Questões de Concursos e Problemas Avançados Resolvidos
8. 🔗 8. Conexões e Referências Cruzadas

---

🧠 1. Paradigma Computacional Híbrido: C para Tempo Real e Python para Modelagem

Na engenharia de sistemas elétricos contemporânea, nenhum projeto complexo é resolvido utilizando exclusivamente uma única linguagem de programação. Adota-se a arquitetura computacional híbrida em duas camadas:

flowchart TD
    subgraph Borda["Camada de Borda / Tempo Real Estrito (Firmware em C99/C++)"]
        C1["Amostragem ADC síncrona a 3840 Hz (DMA)"]
        C2["Filtragem Fasorial DFT de 1 ciclo em ponto fixo/FPU"]
        C3["Malhas de controle PWM em tempo real (< 50 µs)"]
        C4["Algoritmos de desarme de emergência ANSI 50/51/87T"]
    end

    subgraph Centro["Camada Central / Análise e Otimização (Python Científico)"]
        P1["Montagem de grandes redes elétricas (Ybus esparsa)"]
        P2["Fluxo de Potência Não Linear (Newton-Raphson / SciPy)"]
        P3["Despacho Econômico Hidrotérmico e Previsão de Carga"]
        P4["Interface Gráfica de Engenharia e Visualização de Dados"]
    end

    Borda <-->|Telemetria Serial / Modbus / IEC 61850| Centro

• Linguagem C (ISO C99): Oferece execução determinística em tempo real estrito (Hard Real-Time), manipulação cirúrgica de registradores e latências de nanossegundos em relés de proteção IEDs;
• Linguagem Python (NumPy / SciPy / Pandapower): Proporciona altíssima produtividade analítica, suporte nativo a álgebra linear de matrizes esparsas, bibliotecas de otimização convexa e integração com ciência de dados.

---

⚡ 2. Extração Fasorial Discreta de 1 Ciclo (Filtro Cosseno/Seno)

Em relés de proteção digitais, os sinais brutos de corrente e tensão fornecidos pelos transformadores de corrente (TC) e de potencial (TP) contêm harmônicos de chaveamento, ruídos de alta frequência e uma componente exponencial contínua transitória aperiódica decrescente (DC Offset) gerada pelo fechamento de faltas em momentos de cruzamento por zero de tensão.

Para calcular a impedância aparente de linha e atuar nas funções de distância (função ANSI 21) ou diferencial (87), o microcontrolador deve extrair a componente fasorial fundamental de $60\text{ Hz}$.

O padrão universal da indústria é o Filtro da Transformada Discreta de Fourier de Ciclo Completo (Full-Cycle DFT).

2.1 Formulação Matemática Discreta

Considere um sinal senoidal discretizado em $N$ amostras uniformemente distribuídas ao longo de exatamente um ciclo fundamental de $60\text{ Hz}$ ($T=\frac{1}{60}\text{ s}\approx 16,67\text{ ms}$), com taxa de amostragem $\Delta t=T/N$:

$x[k]=x(k\cdot \Delta t),k\in \{0,1,\dots ,N-1\}$

A decomposição ortogonal do fasor fundamental $\bar{X}=X_{\text{real}}+j{X}_{\text{imag}}$ é calculada pelas convoluções discretas com os núcleos de cosseno e seno:

$\boxed{X_{\text{real}}=\frac{2}{N}\sum _{k=0}^{N-1}x[k]\cos \left(\frac{2\pi k}{N}\right)}$ $\boxed{X_{\text{imag}}=-\frac{2}{N}\sum _{k=0}^{N-1}x[k]\sin \left(\frac{2\pi k}{N}\right)}$

2.2 Extração do Módulo RMS e Ângulo de Fase

O módulo eficaz (Valor RMS da Fundamental) e o ângulo de fase angular em relação à base de tempo são obtidos por:

$\boxed{X_{\text{rms}}=\frac{|\bar{X}|}{\sqrt{2}}=\frac{\sqrt{X_{\text{real}}^2+X_{\text{imag}}^2}}{\sqrt{2}}}$ $\boxed{\theta =\text{atan2}(X_{\text{imag}},X_{\text{real}})}$

2.3 Propriedades Matemáticas de Rejeição de Ruído

1. Rejeição Perfeita de Todos os Harmônicos Inteiros: Em virtude da ortogonalidade trigonométrica das funções senoidais discretas no intervalo $[0,2\pi ]$: ${\sum }_{k=0}^{N-1}\cos \left(\frac{2\pi mk}{N}\right)\cos \left(\frac{2\pi nk}{N}\right)=\{\begin{array}{ll}N/2,& \text{se }m=n\ne 0\\ 0,& \text{se }m\ne n\end{array}$ O filtro suprime com atenuação infinita o $2^{\circ },3^{\circ },4^{\circ },\dots ,(N-1{)}^{\circ }$ harmônicos da rede elétrica;
2. Implementação em Janela Deslizante Recursiva (Recursive DFT): Ao adquirir uma nova amostra $x_{\text{novo}}$ e descartar a amostra mais antiga $x_{\text{antigo}}$ (de $N$ passos atrás), o fasor no instante $m+1$ é atualizado a partir do fasor do instante $m$ através de uma simples rotação complexa: ${\bar{X}}_{m+1}=\left({\bar{X}}_m+\frac{2}{N}(x_{\text{novo}}-x_{\text{antigo}})\right)e^{-j\frac{2\pi }{N}}$ Isso reduz a complexidade temporal de cálculo do fasor de $\mathcal{O}(N)$ para $\mathcal{O}(1)$ operações por amostra, permitindo a atualização fasorial em microcontroladores modestos a cada poucos microssegundos.

---

📐 3. Grandezas de Potência em Tempo Real e Aritmética Complexa em C99

3.1 True RMS e Potência Ativa Discreta

Quando a forma de onda da corrente é distorcida por cargas não lineares (inversores de frequência e retificadores), o valor eficaz verdadeiro (True RMS) deve considerar o conteúdo energético de todo o espectro até a frequência de Nyquist:

$\boxed{V_{\text{rms}}=\sqrt{\frac{1}{N}\sum _{k=0}^{N-1}v[k{]}^2},I_{\text{rms}}=\sqrt{\frac{1}{N}\sum _{k=0}^{N-1}i[k{]}^2}}$

A Potência Ativa Média ($P$) em Watts é a média temporal dos produtos instantâneos de tensão e corrente:

$\boxed{P=\frac{1}{N}\sum _{k=0}^{N-1}v[k]\cdot i[k]}$

A Potência Aparente ($S$) em Volt-Ampères e o Fator de Potência ($FP$) são obtidos por:

$S=V_{\text{rms}}\cdot I_{\text{rms}},FP=\frac{P}{S}\le 1,0$

A Potência Reativa ($Q$) segundo a convenção de Budeanu/IEEE 1459 para distorção moderada é dada por:

$Q=\sqrt{S^2-P^2}$

3.2 Aritmética de Números Complexos em C99 (<complex.h>)

A partir do padrão C99, a linguagem C oferece suporte nativo de primeira classe a números complexos, dispensando a criação de estruturas manuais struct { float real, imag; }:

#include <stdio.h>
#include <complex.h>
#include <math.h>
 
void exemplo_fasor(void) {
    // Declaração de fasores fasoriais em C99
    double complex V1 = 13800.0 * cexp(I * 0.0);             // 13.8 kV /_0°
    double complex Z_linha = 1.2 + 4.5 * I;                  // Impedância Z = 1.2 + j4.5 ohms
    
    // Lei de Ohm Fasorial direta com operadores aritméticos normais:
    double complex I_linha = V1 / Z_linha;
    
    // Potência Complexa S = V * conj(I)
    double complex S_carga = V1 * conj(I_linha);
    
    double P = creal(S_carga); // Potência Ativa em Watts
    double Q = cimag(S_carga); // Potência Reativa em VAr
    double S_mag = cabs(S_carga); // Módulo Aparente em VA
    
    printf("P = %.2f kW, Q = %.2f kVAr, S = %.2f kVA\n", P/1000.0, Q/1000.0, S_mag/1000.0);
}

---

🌐 4. Montagem Algorítmica da Matriz Ybus e Armazenamento Esparso

Em estudos de planejamento de transmissão e análise de estabilidade, as redes elétricas reais contam com centenas ou milhares de barras ($N_b\gg 1000$). No entanto, cada barra conecta-se tipicamente a apenas 2 a 4 linhas vizinhas.

A Matriz Nodal de Admitâncias (${\mathbf{Y}}_{\text{bus}}$) é altamente esparsa: mais de $99\%$ de seus elementos numéricos são rigorosamente iguais a zero.

Armazenar uma matriz densa de $5000\times 5000$ números complexos de dupla precisão exigiria: $5000\times 5000\times 16\text{ bytes}=400000000\text{ bytes}\approx 400\text{ MB de RAM}$ Inversões de matrizes cheias consumiriam tempo com complexidade cúbica $\mathcal{O}(N_b^3)$.

Com o armazenamento esparso no formato CSR (Compressed Sparse Row) ou LIL (List of Lists), armazenam-se exclusivamente os termos não nulos ($\approx 20000$ termos), reduzindo o consumo de memória para menos de $500\text{ kB}$ ($>99,8\%$ de economia).

4.1 Implementação Vetorizada da Matriz Ybus em Python

import numpy as np
import scipy.sparse as sp
 
def montar_ybus_esparsa(num_barras, linhas_dados):
    """
    Monta a matriz Ybus esparsa diretamente no formato LIL.
    linhas_dados: lista de tuplas (de_barra, para_barra, r, x, b_shunt)
    """
    Ybus = sp.lil_matrix((num_barras, num_barras), dtype=complex)
    
    for de, para, r, x, b_sh in linhas_dados:
        z = complex(r, x)
        y = 1.0 / z
        b_meio = 1j * (b_sh / 2.0)
        
        # Termos Fora da Diagonal (Negativo da admitância em série)
        Ybus[de, para] -= y
        Ybus[para, de] -= y
        
        # Termos da Diagonal Principal (Soma das admitâncias incidentes)
        Ybus[de, de] += y + b_meio
        Ybus[para, para] += y + b_meio
        
    return Ybus.tocsr() # Converte para CSR para operações rápidas de produto

---

📈 5. Solução de Fluxo de Potência por Newton-Raphson com Fatoração LU

O problema do Fluxo de Carga em regime permanente consiste em determinar o vetor de magnitudes de tensão $\mathbf{V}$ e ângulos de fase $\mathbf{\theta }$ em todas as barras da rede, dadas as gerações e demandas de potência ativa e reativa especificadas.

As equações de injeção de potência na barra $k$ em coordenadas polares são não lineares:

$P_k(\mathbf{\theta },\mathbf{V})=V_k{\sum }_{m=1}^{N_b}{V}_m(G_{km}\cos ({\theta }_k-{\theta }_m)+B_{km}\sin ({\theta }_k-{\theta }_m))$ $Q_k(\mathbf{\theta },\mathbf{V})=V_k{\sum }_{m=1}^{N_b}{V}_m(G_{km}\sin ({\theta }_k-{\theta }_m)-B_{km}\cos ({\theta }_k-{\theta }_m))$

Onde $Y_{km}=G_{km}+j{B}_{km}$.

5.1 O Algoritmo de Newton-Raphson

A cada iteração $v$, o sistema não linear é linearizado em torno do ponto de operação atual através da expansão em Série de Taylor truncada de primeira ordem:

$\left[\begin{array}{c}{\mathbf{\Delta P}}^{(\nu )}\\ {\mathbf{\Delta Q}}^{(\nu )}\end{array}\right]=\left[\begin{array}{cc}{\mathbf{H}}^{(\nu )}& {\mathbf{N}}^{(\nu )}\\ {\mathbf{M}}^{(\nu )}& {\mathbf{L}}^{(\nu )}\end{array}\right]\left[\begin{array}{c}{\mathbf{\Delta }\mathbf{\theta }}^{(\nu )}\\ \frac{{\mathbf{\Delta V}}^{(\nu )}}{\mathbf{V}}\end{array}\right]$

Onde os resíduos de potência são calculados por: $\Delta P_k^{(\nu )}=P_k^{\text{especificado}}-P_k({\mathbf{\theta }}^{(\nu )},{\mathbf{V}}^{(\nu )})$ $\Delta Q_k^{(\nu )}=Q_k^{\text{especificado}}-Q_k({\mathbf{\theta }}^{(\nu )},{\mathbf{V}}^{(\nu )})$

E a Matriz Jacobiana $\mathbf{J}=\left[\begin{array}{cc}\mathbf{H}& \mathbf{N}\\ \mathbf{M}& \mathbf{L}\end{array}\right]$ possui as derivadas parciais: $H_{km}=\frac{\partial P_k}{\partial {\theta }_m},N_{km}=V_m\frac{\partial P_k}{\partial V_m},M_{km}=\frac{\partial Q_k}{\partial {\theta }_m},L_{km}=V_m\frac{\partial Q_k}{\partial V_m}$

flowchart TD
    Init["1. Inicialização do Flat-Start:<br/>θ_k = 0 rad, V_k = 1.0 pu para k ≠ Slack"] --> CalcRes["2. Calcular Injeções de Potência e Resíduos:<br/>ΔP_k = P_k^esp - P_k^calc<br/>ΔQ_k = Q_k^esp - Q_k^calc"]
    CalcRes --> TesteConv{"3. Teste de Convergência:<br/>max(|ΔP|, |ΔQ|) < 10^-4 pu?"}
    TesteConv -- Sim --> Sucesso["Convergência Atingida!<br/>Calcular fluxos de potência nos ramos e perdas Joule"]
    TesteConv -- Não --> MontaJac["4. Avaliar Matriz Jacobiana Esparsa J"]
    MontaJac --> FatoracaoLU["5. Resolver Sistema Linear por Fatoração LU:<br/>J * [Δθ; ΔV/V] = [ΔP; ΔQ]"]
    FatoracaoLU --> Atualiza["6. Atualizar Variáveis de Estado:<br/>θ^(v+1) = θ^(v) + Δθ<br/>V^(v+1) = V^(v) + ΔV"]
    Atualiza --> CalcRes

O algoritmo de Newton-Raphson exibe convergência quadrática local: uma vez que o vetor de estados entra na região de atração da solução física, o número de dígitos corretos dobra a cada iteração, garantindo convergência em tipicamente 3 a 5 iterações para qualquer tamanho de sistema interligado.

---

📊 6. Diagrama de Algoritmos Numéricos, Fasores e Solução de Redes

O diagrama técnico abaixo exibe o ciclo completo do processamento numérico de sinais elétricos em C99 e a resolução iterativa de fluxo de carga em Python:

---

📝 7. Questões de Concursos e Problemas Avançados Resolvidos

Questão 1 (Banca CESGRANRIO - Engenheiro de Eletricidade / Eletrobras)

Em um relé de proteção digital de sobrecorrente temporizada, o sinal de corrente alternada de $60\text{ Hz}$ de uma fase é digitalizado com uma taxa de amostragem de exatamente $N=16\text{ amostras por ciclo}$. As 16 amostras numéricas obtidas em um ciclo fundamental de carga resistivo-indutiva equilibrada foram processadas pelo algoritmo da DFT de 1 ciclo, resultando nos seguintes termos ortogonais de Fourier: $X_{\text{real}}=141,42\text{ A},X_{\text{imag}}=-141,42\text{ A}$ Determine:

1. O valor fasorial retangular da corrente fundamental ($\bar{I}$).
2. O valor eficaz verdadeiro (RMS) da componente fundamental da corrente em ampères ($I_{\text{rms}}$).
3. O ângulo de fase da corrente em graus ($\theta$).

Resolução Passo a Passo:

1. Fasor Retangular Fundamental ($\bar{I}$): $\bar{I}=X_{\text{real}}+j{X}_{\text{imag}}=141,42-j141,42\text{ A}$

2. Cálculo do Módulo RMS da Fundamental ($I_{\text{rms}}$): O módulo do fasor de pico é: $|\bar{I}|=\sqrt{X_{\text{real}}^2+X_{\text{imag}}^2}=\sqrt{(141,42{)}^2+(-141,42{)}^2}$ $|\bar{I}|=\sqrt{20000+20000}=\sqrt{40000}=200,0\text{ A (Valor de Pico)}$

   O valor eficaz RMS é obtido dividindo por $\sqrt{2}$: $I_{\text{rms}}=\frac{|\bar{I}|}{\sqrt{2}}=\frac{200,0}{1,41421}\approx 141,42\text{ A}$

3. Cálculo do Ângulo de Fase ($\theta$): $\theta =\arctan \left(\frac{X_{\text{imag}}}{X_{\text{real}}}\right)=\arctan \left(\frac{-141,42}{141,42}\right)=\arctan (-1,0)$ Como o componente real é positivo e o imaginário é negativo, o fasor reside no 4º quadrante: $\theta =-45^{\circ }\text{ou}315^{\circ }$

$\boxed{\bar{I}=141,42-j141,42\text{ A};I_{\text{rms}}=141,42\text{ A};\theta =-45^{\circ }}$

Interpretação: A corrente de $141,42{\text{ A}}_{\text{rms}}$ está atrasada de $45^{\circ }$ em relação à referência de tempo, correspondendo perfeitamente a um circuito com fator de potência indutivo $\cos (-45^{\circ })=\frac{\sqrt{2}}{2}\approx 0,707$.

---

Questão 2 (Banca CESPE/CEBRASPE - Engenheiro Eletricista / Perito Criminal)

Considere uma linha de transmissão entre as barras 1 e 2 de um sistema de potência com impedância em série $Z_{12}=0,02+j0,08\text{ pu}$ e capacitância total de carregamento desprezível ($b_{\text{shunt}}=0$). Calcule a admitância de ramo $y_{12}$ em sua forma retangular ($g_{12}+j{b}_{12}$) e determine a contribuição dessa linha nos elementos $Y_{11}$ e $Y_{12}$ da matriz nodal de admitâncias ${\mathbf{Y}}_{\text{bus}}$.

Resolução Passo a Passo:

1. Cálculo da Admitância em Série ($y_{12}=\frac{1}{Z_{12}}$): $y_{12}=\frac{1}{0,02+j0,08}$ Multiplica-se numerador e denominador pelo complexo conjugado $(0,02-j0,08)$: $y_{12}=\frac{0,02-j0,08}{(0,02{)}^2+(0,08{)}^2}=\frac{0,02-j0,08}{0,0004+0,0064}=\frac{0,02-j0,08}{0,0068}$

   • Condutância: $g_{12}=\frac{0,02}{0,0068}\approx 2,9412\text{ pu}$
   • Susceptância: $b_{12}=-\frac{0,08}{0,0068}\approx -11,7647\text{ pu}$ $y_{12}=2,9412-j11,7647\text{ pu}$

2. Contribuição na Matriz Nodal de Admitâncias:

   • Termo fora da diagonal ($Y_{12}=-y_{12}$): $Y_{12}=-(2,9412-j11,7647)=-2,9412+j11,7647\text{ pu}$
   • Contribuição no termo da diagonal principal ($Y_{11}+=y_{12}$): $Y_{11}\text{ recebe }+y_{12}=+2,9412-j11,7647\text{ pu}$

$\boxed{y_{12}=2,9412-j11,7647\text{ pu};Y_{12}=-2,9412+j11,7647\text{ pu}}$

---

Questão 3 (Problema Analítico de Engenharia - Convergência de Newton-Raphson)

Em um fluxo de potência de uma rede de duas barras contendo uma barra de referência (Barra 1: $V_1=1,0\angle 0^{\circ }\text{ pu}$) e uma barra de carga PQ (Barra 2: $P_2^{\text{esp}}=-1,0\text{ pu},Q_2^{\text{esp}}=-0,5\text{ pu}$), interligadas por uma linha puramente reativa $Z_{12}=j0,2\text{ pu}$ ($Y_{bus}=\left[\begin{array}{cc}-j5& +j5\\ +j5& -j5\end{array}\right]$): Iniciando a partir do Flat-Start ($V_2^{(0)}=1,0\text{ pu},{\theta }_2^{(0)}=0\text{ rad}$), determine o resíduo de potência ativa $\Delta P_2^{(0)}$ e o elemento Jacobiano $H_{22}=\frac{\partial P_2}{\partial {\theta }_2}$ na primeira iteração.

Resolução Passo a Passo:

1. Equação da Potência Ativa Calculada na Barra 2 ($P_2^{\text{calc}}$): Como $G_{12}=G_{22}=0$ e $B_{21}=+5$, $B_{22}=-5$: $P_2^{\text{calc}}=V_2(V_1{B}_{21}\sin ({\theta }_2-{\theta }_1)+V_2{B}_{22}\sin (0))$ $P_2^{\text{calc}}=V_2{V}_1(5)\sin ({\theta }_2-{\theta }_1)$

2. Avaliação no Ponto Inicial Flat-Start ($V_1=1,0,V_2=1,0,{\theta }_1=0,{\theta }_2=0$): $P_2^{\text{calc}(0)}=1,0\times 1,0\times 5\times \sin (0-0)=0,0\text{ pu}$

3. Cálculo do Resíduo de Potência Ativa ($\Delta P_2^{(0)}$): $\Delta P_2^{(0)}=P_2^{\text{especificado}}-P_2^{\text{calc}(0)}=-1,0-0,0=-1,0\text{ pu}$

4. Derivada Parcial do Jacobiano ($H_{22}=\frac{\partial P_2}{\partial {\theta }_2}$): Diferenciando a expressão de $P_2$ em relação a ${\theta }_2$: $H_{22}=\frac{\partial }{\partial {\theta }_2}[5V_2{V}_1\sin ({\theta }_2-{\theta }_1)]=5V_2{V}_1\cos ({\theta }_2-{\theta }_1)$ Substituindo os valores do Flat-Start: $H_{22}^{(0)}=5\times 1,0\times 1,0\times \cos (0)=5,0\text{ pu}$

$\boxed{\Delta P_2^{(0)}=-1,0\text{ pu};H_{22}^{(0)}=5,0\text{ pu}}$

Avanço do Algoritmo: A atualização angular de primeira iteração calculada por Newton-Raphson será: $\Delta {\theta }_2^{(0)}=\frac{\Delta P_2^{(0)}}{H_{22}^{(0)}}=\frac{-1,0}{5,0}=-0,20\text{ rad}\approx -11,46^{\circ }$ Demonstrando como o fluxo de potência ativo é transmitido do ângulo mais adiantado (Barra 1 a $0^{\circ }$) para o ângulo mais atrasado (Barra 2 a $-11,46^{\circ }$).

---

🔗 8. Conexões e Referências Cruzadas

Disciplinas Vinculadas:

• Aula 02: Representação IEEE 754 e ponto flutuante.
• Aula 03: Buffers de amostragem DMA e decaimento de vetores.
• 07 Planejamento e Análise de SEP: Métodos de Newton-Raphson, Desacoplado Rápido e FPO.
• 07 Modelagem Matemática: Resolução de sistemas de equações algébricas lineares e matrizes esparsas.

Materiais Didáticos Oficiais da Disciplina:

• Arquivo de Curso: 02 - Linguagem C Recursos Básicos.pdf
• Arquivo de Curso: 03 - Linguagem C Funções.pdf
• Arquivo de Curso: 04 - Linguagem C Nos Sitemas Embarcados.pdf