Calcolo dei flussi di potenza

Il calcolo dei flussi di potenza viene utilizzato per la risoluzione di una rete, ovvero consente di calcolare lo stato della stessa. Noti i parametri, quindi carichi e impedenze di linea, si procede con il calcolo delle correnti sui rami e delle tensioni ai nodi per verificare la presenza di violazioni (superamento delle portate termiche, tensioni di alimentazione fuori dalla fascia di tolleranza ecc.).

Ovviamente il BFS, basato sul calcolo fasoriale, fornisce risultati in condizione di regime.

Una rete di distribuzione, generalmente MT, presenta elementi trasversali non trascurabili (a differenza di una rete di trasmissione, quindi in AT). Gli elementi trasversali di una linea sono prevalentemente capacitivi.

Data una rete radiale, ovvero in assenza di loop, come la seguente:

Nota: i rami vengono nomenclati come il nodo di arrivo.

Si definisce tramite ispezione visiva (quindi partendo dal grafo) la matrice delle incidenze L, i cui elementi assumono i seguenti valori:

• -1 nodo di arrivo del ramo;
• +1 nodo di partenza del ramo;

Facendone l’inversa:

Letta per righe, la matrice gamma, considerando gli elementi non nulli, permette di ricostruire il percorso tra il nodo e la radice.

Ogni ramo del grafo rappresenta una linea avente parametri specifici quali lunghezza e impedenze, parametrizzato tramite il modello a “pi-greca”. La suscettanza capacitiva trasversale totale viene “divisa” tra ingresso e uscita:

La matrice Zb conterrà sulla diagonale le impedenze longitudinali (R+JX).

A questo punto si può notare che moltiplicando per la matrice gamma, quest’ultima “filtra” le longitudinali: Il risultato, di fatto, contiene solamente elementi relativi alla porzione di rete che va dal nodo considerato alla radice:

Moltiplicando, quindi, per le correnti di linea si ottiene la caduta di tensione (dal nodo alla radice).

Il nodo i-esimo presenta la seguente struttura:

Le correnti nodali Is vengono calcolate dalla tensione nodale Vi (che è incognita, sia in modulo che in fase), noti i carichi (impedenza assegnata Zci, potenza assegnata Pci+jQci e corrente assegnata Ici). In seguito possono essere calcolate, tramite la matrice , le correnti sui rami (Ib). Ora, passando alla scrittura matriciale:

Notare che, in notazione matriciale:

Quindi, si può scrivere:

Il BFS è costituito da due fasi:

• Calcolo delle correnti (Bacward): Nella fase di backard della k-esima iterazione si calcolano le correnti nodali partendo dalle tensioni nodali ottenute nell’iterazione k-1 (per la prima iterazione, la tensione di riferimento rappresenta la “condizione iniziale” impostata a monte del ciclo iterativo).

Nota: la componente trasversale della linea, in questo esempio, è consitrerata come inglobata con il carico (dal momento che si trova in parallelo con esso). Da un punto di vista matematico la struttura non cambia: qualora si volesse esplicitare Bc, è sufficiente introdurre un vettore contenente tale valore per ogni nodo, calcolare la corrente che lo attraversa in funzione della tensione nodale e sommarla ad Is.

• Calcolo delle tensioni nodali (Forward): in sostanza, note le correnti che attraversano i rami:

La tensione ai nodi nell’iterazione k si ricava calcolando la caduta di tensione dallo radice:

Il ciclo viene interrotto quando si ha convergenza, ovvero quando la differenza tra il risultato all’iterazione K e quella all’iterazione K-1 è inferiore ad un certo valore:

In python il ciclo è il seguente:

#Funzione di calcolo BFS():
def BFS(Zci):
    V = np.ones(len(L))#Condizioni iniziali
    tol = pow(10,-7)#errore minimo 
    dv = 10#In modo da entrare nel ciclo
    i = 0
    n = 1000
    while(dv > tol):
        Is = np.divide(V,Zci)#Backward
        Ib = np.dot(np.transpose(gamma),Is)
        Vaux = one - np.dot(zbus,Is)#Forward
        dv =np.max( abs( np.divide((Vaux -V ),Vaux )) )
        dv = np.max(    np.divide(abs(Vaux -V ),abs(Vaux ))                 )
        V = Vaux
        i +=1
        if(i>n):
            print("overflow!")
            break
    return V,Ib

Il BFS converge ma è previsto un break qualora l’errore richiesto sia eccessivamente piccolo rihciedendo (per convergere) un numero eccessivo di iterazioni. Notare che il calcolo delle correnti nodali viene eseguito tramite l’operatore di Hadamard.

Il codice completo è il seguente:

import numpy as np
import matplotlib.pyplot as plt
import math 
import cmath 
#Matrice delle incidenze
L = np.array([[-1,0,0,0,0,0,0,0,0],
              [1,-1,0,0,0,0,0,0,0],
              [1,0,-1,0,0,0,0,0,0],
              [0,1,0,-1,0,0,0,0,0],
              [0,1,0,0,-1,0,0,0,0],
              [0,0,1,0,0,-1,0,0,0],
              [0,0,1,0,0,0,-1,0,0],
              [0,0,0,0,0,0,1,-1,0],
              [0,0,0,0,0,0,1,0,-1]])
#Controllo radialità della rete tramite det(L)
if (np.linalg.det(L) != 0): 
     gamma =np.linalg.inv(L)
#NOTA: TUTTI I VALORI SONO IN p.u.
#impedenze longitudinali della linea: tutte uguali!
Zl =pow(10,-3)*complex(5,3)*np.ones(9) 
zb = np.diag(Zl)
#impedenze bus
zbus=np.dot(np.dot(gamma,zb),np.transpose(gamma))
V0 =1 #Profilo piatto!
one = V0*np.ones(np.size(L,1))
  
        
#Funzione di calcolo BFS():
def BFS(Zci):
    V = np.ones(len(L))#Condizioni iniziali
    tol = pow(10,-7)#errore minimo 
    dv = 10#In modo da entrare nel ciclo
    i = 0
    n = 1000
    while(dv > tol):
        Is = np.divide(V,Zci)#Backward
        Ib = np.dot(np.transpose(gamma),Is)
        Vaux = one - np.dot(zbus,Is)#Forward
        dv =np.max( abs( np.divide((Vaux -V ),Vaux )) )
        dv = np.max(    np.divide(abs(Vaux -V ),abs(Vaux ))                 )
        V = Vaux
        i +=1
        if(i>n):
            print("overflow!")
            break
    return V,Ib

x = 10
y = 0.5
load = complex(x,y) #Singolo carico

Zci = load*np.ones(len(L)) #Carichi: tutti uguali

V , Is = BFS(Zci)

Dato che le impedenze di linea, le matrici L e gamma sono dichiarate come variabili globali, è sufficiente chiamare la funzione BFS(), con argomento il vettore dei carichi, per ottenere correnti sui rami e tensioni nodali. Ad esempio:

x = 10
y = 0.5
load = complex(x,y) #Singolo carico
Zci = load*np.ones(len(L)) #Carichi: tutti uguali

Il risultato, ad esempio, è il seguente:

Si osservi che le operazioni tra matrici, piuttosto pesanti, ove possibile sono state svolte prima del ciclo iterativo (per evitare di ripeterle).

La bontà del codice è deducibile da alcune considerazioni. Ipotizzando tutti carichi uguali :

• Confrontando le correnti di linea con carichi fortemente reattivi, quindi con bassi fattori di potenza, con le correnti di linea in caso di “rifasamento”, si osserva che nel secondo caso le correnti sono minori.

• Azzerando le impedenze di linea (matrice Zb nulla) le cadute di tensioni sono nulle e le tensioni ai nodi sono tutte pari ala tensione del nodo di alimentazione (slack);

• Imponendo carichi puramente capacitivi si osserva la sussistenza dell’effetto Ferranti, più marcato nella periferia della linea.

L’effetto Ferranti si manifesta in presenza di reattanze capacitive lungo la linea. Più è lunga la linea, maggiore sarà la reattanza trasversale (capacitiva) totale. L’esempio del solo carico capacitivo non è molto distante dalla realtà: si ipotizzi, ad esempio, il distacco di un carico associato al relativo rifasataore. La rete “vede” il rifasatore come carico e quindi è a tutti gli effetti un carico capacitivo. Personalmente, ho avuto modo di visionare un profilo di carico di una grossa azienda in cui la potenza reattiva assorbita durante la notte era negativa: questo perchè i macchinari erano spenti ma i rifasatori erano comunque alimentati.

• Imponendo i carichi a potenza assegnata (quali generatori locali) le tensioni dei nodi in prossimità del generatore possono variare.

Modificando leggermente il codice, inserendo un generatore al nodo 8:

Sg = complex(0,0)*np.ones(9) #Generatori
 P =1.3
 Q=0.2
 Sg[7] = complex(P,Q)#Generatore al nodo 8

La funzione BFS diventa:

#Funzione di calcolo BFS():
 def BFS(Zci):
     V = np.ones(len(L))#Condizioni iniziali
     tol = pow(10,-7)#errore minimo 
     dv = 10#In modo da entrare nel ciclo
     i = 0
     n = 1000
     while(dv > tol):
         Is = np.divide(V,Zci) + np.conj(np.divide(Sg,V))#Backward#Backward
         Ib = np.dot(np.transpose(gamma),Is)
         Vaux = one - np.dot(zbus,Is)#Forward
         dv =np.max( abs( np.divide((Vaux -V ),Vaux )) )
         dv = np.max(    np.divide(abs(Vaux -V ),abs(Vaux ))                 )
         V = Vaux
         i +=1
         if(i>n):
             print("overflow!")
             break
     return V,Ib

Con la generazione locale, quindi, è possibile gestire la rete in modo da riportare le tensioni di alimentazione ai nodi lontani dallo slack ai valori “richiesti” .

Anche le correnti subiscono variazioni. In alcuni tratti di linea camnbiano addirittura verso: questo vuol dire che, i nodi in prossimità del generatore (in questo caso l’8) sono alimentati dallo stesso e non dallo slack.

Nel caso esposto, si tratta di un generatore talmente potente (eccessivamente) da consentire di scaricare completamente il primo tratto di linea.

Una volta definite tensioni è correnti, è noto stato energetico della rete. Ad esempio, considerando la corrente sul ramo 1, ovvero quello che collega lo slack al nodo principale, è nota la potenza fornita dal generatore (slack):

La tensione sul nodo di slack presenta solo parte reale dal momento che è stato considerato come origine delle fasi. La potenza fornita all’intera rete è pari a:

Volendo, è possibile calcolare anche il fattore di potenza, dall’argomento della stessa:

Le cose tornano dal momento che imponendo carichi puramente resistivi e azzerando la componente reattiva in Zb il fattore di potenza è unitario (la rete e i carichi rappresenterebbero per il generatore un carico puramente resistivo).

BFS probabilistico

I carichi non sono noti a priori, è necessario conoscere i profili di carico associati al tipo di carico. Qualora fosse presente, ad esempio, un parco eolico, la potenza generata dipende dalla velocità del vento, che è un parametro aleatorio.

Nella figura, viene mostrato il processo di “estrazione” della potenza generata al nodo 4 da un generatore eolico.

• Si estrae un valore casuale r che va da 0 a 1, verificandone quindi la disponibilità;
• Se disponibile, si estrae nuovamente un numero casuale (r) e dalla CDF della velocità del vento, si ricava la “velocità del vento estratta”;
• Dalla curva del costruttore, che lega velocità del vento e potenza prodotta, si risale (tramite la velocità del vento estratta nel punto precedente) alla potenza prodotta.

Il fattore di potenza viene mantenuto costante. Questo procedimento viene ripetuto ad ogni ciclo, in ogni nodo ove sono presenti generatori (fotovoltaici, eolici, turbogas ecc).

L’algoritmo viene così sintetizzato:

1. Estrazione potenza prodotta (in ogni nodo) e costruzione della rete con i parametri estratti;
2. Calcolo BFS;
3. “Aggiornamento classe”;
4. Nuova estrazione;

Praticamente, si calcola il valore della “classe” per \textbf{costruirne una PDF}.

Notare che la PDF presenta area unitaria: quindi è necessario dividere per il numero di ripetizioni e per l’ampiezza della classe.

La CDF sarà la seguente:

Dalla quale si può ricavara la probabilità che la tensione scenda sotto ad un determinato valore.

Matteo Gentileschi