Monitorizarea riscului declanșării instabilității de tensiune [615542]

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
060042 București, Splaiul I ndependenței, nr. 313, sector 6

PROIECT DE DIPLOMĂ

Monitorizarea riscului declanșării instabilității de tensiune

Autor: Victor SAV

Cadru didactic îndrumător: Prof. dr. ing. Constantin BULAC

București
Iulie 2018

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

2

Proiect de diplomă

Monitorizarea riscului declan șării instabilit ății de
tensiune

prezentat la

Universitatea POLITEHNICA din București
Facultate de Energetică

pentru obținerea titlului de
inginer

Specializarea : Informatică industrială

de către

Victor SAV
(absolvent: [anonimizat])

sub îndrumarea
Prof. dr. ing. Constantin BULAC

Susținut la data de 4 iulie 2018 , în fa ța comisiei de examinare

Conf. dr. ing. Radu PORUMB Președinte
Prof. dr. ing. Constantin BULAC Membru
Prof. dr. ing. Virgil DUMBRAVĂ Membru
Prof. dr. ing. Ion TRIȘTIU Membru
Conf. dr. ing. George Cristian LĂZĂROIU Membru
Ș.l. dr. ing. Nicoleta ARGHIRA Membru
S.l. dr. ing. Alexandru MANDIȘ Secretar

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

3

Această lucrare a fost pregătită în cadrul Facultății de Energetică a UPB

(model de cuprins)

Pag.
INTRODUCERE
1.
2.

CAPITOLUL n
n.1 DENUMIRE PARAGRAF
n.1.1. Denumire subparagraf de ordin 1
……………………………………..
n.2. Denumire paragraf
n.2.2. Denumire subparagraf
……………………………………………………………….
n.2.2.1. Denumire subparagraf de ordin 2
CONCLUZII

BIBLIOGRAFIE

ANEXE

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

4

Instrucțiuni pentru elaborarea proiectului de diplomă

Informații generale
• Proiectul de diplomă se redactează pe format A4, la un rând, „Times New Roman”, mărimea
literei 12. Marginile stânga și dreapta 2,5 cm, marginile sus și jos 2 cm, paginarea lucrării la
mijlocul paginii – jos.
• In “Header” se trece titlul proiectului de dip lomă, eventual prescurtat, astfel încât să încapă pe un
singur rând, cu litere „Times New Roman”, dimensiunea 10.
• Denumirile figurilor, relațiile și tabelele se scriu cu „Times New Roman”, dimensiunea 11 și se
numerotează în raport cu fiecare capitol (Exe mplu pentru capitolul 2: Fig. 2.1, Fig. 2.2, … ; relația
(2.1), (2.2), …..; Tabelul 2.1, …..).
• Între Copertă și Cuprins se introduce prima pagină a proiectului conform modelului din prezenta
anexă . Cuprinsul lucrării începe cu pagina 3, apoi lista d e notații (dacă este cazul) și introducerea.
•
• Fiecare capitol al lucrării începe cu pagină nouă.
• DENUMIREA CAPITOLULUI (Centrat, majuscule, bold, dimensiune 14, la 5 cm de marginea
de sus). Textul începe după titlu la 2 cm.
• Titlurile paragrafelor și sub-paragrafelor se scriu după modelul din Cuprins (caractere, felul
caracterului, mărimea, alinierea).
• Ecuații
▪ Ecuațiile se numerotează aliniat la dreapta și se scriu cu editorul de ecuații;
▪ Setări:
▪ size: 12, 10, 7, 18,12;
▪ style: variabile italice,
▪ LC Greek: italice.
• Tabele
▪ Tabelele se scriu cu caractere de11 pct. Titlul tabelului se scrie centrat, bold, deasupra tabelului
cu „Times New Roman” dimensiunea 11.
• Figuri
▪ Figurile trebuie plasate centrat, de preferință la începutul paginii. Denumirea fi gurii se scrie
sub figură, cu caractere de 11 pct. (Exemplu: Fig. 2.1. Schema instalației ______)
• Bibliografia , scrisă cu „Times New Roman” dimensiunea 11, se numerotează cu cifre arabe în
ordinea citării în text sau în ordine alfabetică. Mențiunea în textul lucrării a unei referințe
bibliografice se face cu paranteze drepte (Exemplu pentru marcarea lucrării cu nu mărul de ordine
12: [12]). Numele și prenumele autorilor se scriu cu caractere italice.
Se vor indica mai jos câteva exemple de citare ale unor lucrări, pe categorii: cărți, lucrări în reviste,
lucrări la conferințe ( congrese, simpozioane etc.)
[1] Cristina Saviana , Embedding Data and Task Parallelism in Image Processing Applications, PhD
Thesis, Technische Universiteit Delft, 2003.
[2] A. Mauthe, D. Hutchison, G. Coulson and S. Namuye , “Multimedia Group Communications Towards
New Services”, in Dist ributed Systems Eng., vol. 3 , no. 3, Sept. 1996, pp. 197 -210.
[3] V. I. Arnold , Metodele matematice ale mecanicii clasice, Editura Științifică și Enciclopedică,
București, 1980
[4] V. Gioncu, M. Ivan, Teoria comportării critice și postcritice a structuri lor elastice, Editura Academiei,
București, 1984.
[5] *** COSMOS/M – Finit Element System, User Guide, 2005.

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

5

INTRODUCERE
3. Sistemul energetic
Sistemul energetic (SE) este definit ca fiind ansamblul instalațiilor/elementelor care
participă la procesele de producere, transport, distribuție și utilizare a energiei.
Sistemul electroenergetic (SEE) este un subsistem al sistemului energetic și este definit
ca fiind totalitatea instalațiilor/elementelor electrice care participă la procesele de producere,
transport, distribuție și utilizare a energiei electrice. [1]
Elementele de bază ale structurii unui SEE sunt următoarele:

Fig. 1.1 Structura un ui SEE

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

6

4. Sistemul energetic național
Structura rețelei electrice de transport a SEN este următoarea:
Fig 1.2 Sistemul energetic national
Sistemul energetic național este constituit din:
• Subsisteme de producere
• Subsisteme de transport
• Subsisteme de distribuție
Aceste subsisteme reprezintă cele mai extinse și complexe sisteme dinamice create de
către om, deoarece rețeaua electrică de transport cuprinde sute și chiar mii de noduri la care sunt
racordate centralele electrice și/sau zonele de consum, ia r generatoarele centralelor electrice
aflate la distanțe mari (zeci, sute sau chiar mii de kilometrii ) unele de altele, trebuie să
funcționeze în sincron și să asigure în orice moment echilibrul dintre puterea produsă și cea
consumată. [1]
Ținând cont de cele menționate rezultă că, aspectele tehnico -economice legate de
controlul, planificarea și exploatarea sistemelor electroenergetice sunt extrem de complexe și
vizează pe de o parte problemele fiecărui subsistem în parte, dar mai ales interacțiunile dintre
acestea care definesc domeniul de activitate al ingineriei de sistem. [1]

Fig 1. 3 Interconexiunea sistemelor

Subsistemul de producere
Subsistemul de transport
Subsistemul de distribuție
Ingineria de sistem ( control,
planificarea si exploatarea SEE )

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

7

Proiectarea și exploatarea corespunzătoare a SEE trebuie să satisfacă următoarele
cerințe:
• Sistemul trebuie să fie capabil să satisfacă în orice moment cererea consumatorilor
de putere activă și reactiva. Deoarece energia electrică nu poate fi stocată în cantități
mari, trebuie ca sistemul să dispună de rezerve corespunzătoare și controlabile în
orice moment pentru a acoperii cererea.
• Energia furnizată de sistem către consumatori trebuie să aibă un cost și un impact
ecologic minim
• Calitatea energiei furnizate trebuie să satisfacă standardele impuse de lege (
menținerea frecvenței, nivelului de tensiune și nivelului de fiabilitate )
O soluție pentru a satisface aceste cerințe este controlarea pe niveluri ierarhice a
sistemelor electro energetice prin folosirea dispozitivelor automate de comandă și control,
menținând frecvența, tensiuniile și alte variabile în limitele admisibile. De a semenea ele au un
impact major asupra performanțelor dinamice ale sistemului și a abilității acestuia de a face față
perturbațiilor. [1]
Controlul subsistemului de producere are ca scop principal asigurarea balanței dintre
puterea activă produsă și cea co nsumată ( inclusiv pierderile ), astfel încât frecvența să fie
menținută la valoarea impusă iar circulația de putere activă pe liniile de interconexiune cu
sistemele vecine să fie menținută la valorile stabilite cu aceștia.
Pentru subsi stemul de transpor t, controlul are ca principal obiectiv menținerea unui
nivel corespunzător al tensiunii în întregul sistem precum și controlul circulațiilor de puteri.
Ținând cont de cele evocate anterior, capacitatea SEE de a face față la evenimente care
vizează stabil itatea și integritatea acestuia fără a întrerupe alimentarea consumatorilor, dar și
evaluarea stării de echilibru din punct de vedere al încadrării în limitele admisibile ale
tensiunilor, frecvenței și încărcării elementelor rețelei electrice. [1]
Din punct de vedere al securității, SEE se pot afla în una din următoarele cinci stări de
funcționare:
• Normală – În care toți parametrii de stare ai SE sunt în limitele admisibile, iar
sistemul este capabil să facă față oricărei contingențe fără să apăra depășiri a le
limitelor admisibile [1]
• Alertă – Toți parametrii sunt în limitele admisibile și toate restricțiile sunt satisfăcute,
însă datorită nivelului scăzut de securitate, o contingență poate duce la
supraîncărcarea unor echipamente. În funcție de gradul de sev eritate a perturbației
sistemul poate trece în starea de urgență sau extremă urgență În această stare se iau
măsuri preventive pentru a readuce sistemul în starea normală de funcționare [1]
• Urgență – Sistemul poate fi readus în starea de alertă prin aplica rea măsurilor în caz
de urgență ( eliminarea defectului, controlul excitației, acționarea rapidă a vanelor
de admisie a agentului primar, pornirea unor grupuri, etc…) [1]
• Extremă urgență – Dacă măsurile luate în starea de urgență nu au efect atunci se
ajunge în starea de extremă urgență, care este caracterizată de deconectări în cascadă
de echipamente și prăbușirea unor zone importante a sistemului. În această situație
se impun următoarele acțiuni : deconectarea de sarcini, insularizarea sistemului (
pentru limitarea extinderii avariilor ) [1]
• Restaurare – După eliminarea defectelor se face reconectarea echipamentelor
deconectate și restabilirea alimentării consumatorilor afectați.

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

8

Capitolul 1
Calcul ul regimului permanent folosind metoda Newton -Raphson.
Studii de caz pe un sistem test.

1.1. Regimul permanent
Regimul permanent al SEE poate fi definit ca regimul simetric sinusoidal de secvență
pozitivă față de fazele a,b,c ale sistemului trifazat. Realizarea unui astfel de regim are o
probabilitate scăzută în practică. Studiul acestui regim are o importanță ridi cată pentru că
oferă informații privind starea de echilibru a sistemului la un moment dat precum și
condițiile inițiale pentru determinarea stabilității SEE. [1]
Pentru calculul regimului permanent se vor avea în vedere următoarele particularități:
• SEE su nt constituite din elemente trifazate între care nu există cuplaje electrice și
magnetice.
• Elementele trifazate fiind construite simetric față de fazele a,b,c și pământ se
utilizează pentru studiul regimurilor simetrice și echilibrate, cum este și regimul
permanent.
• Elementele active (Generatoare sau consumatori) ale rețelei de succesiune pozitivă
sunt conectate la un conductor fictiv de nul.

1.2. Modelarea elementelor pentru calculul regimului permanent
Modelarea liniilor electrice -în practică pentru calculul regimurilor de funcționare ale
SEE, liniile electrice se reprezintă printr -un cuadripol în Π [3]

Fig 1.2.1 Schema echivalenta a unei linii electrice
Calculul parametrilor liniei se face cu următoarele formule:
• Impedanța longitudinală a liniei:
𝑍𝑖𝑘=(𝑟0+𝑗𝑥0)𝐿𝑛 [Ω] (1.2.1)
• Admitanța longitudinală a liniei:
𝑦𝑖𝑘=1
𝑍𝑖𝑘𝑍𝑏𝑛 [𝑢.𝑟.] (1.2.2)
• Admitan ța transversal ă a liniei:
𝑦𝑖𝑘0=𝑦𝑘𝑖0=1
2(𝑔0+𝑗𝑏0)𝐿𝑛𝑍𝑏𝑛 [𝑢.𝑟.] (1.2.3)
Dacă linia este dublu circuit atunci pentru fiecare linie se va considera câte un cuadripol.
Modelarea transformatoarelor electrice -transformatoarele electrice se reprezintă printr –
un cuadripol în pi (Π) sau în gamma (Γ) în serie cu un transformator ideal reprezentat printr -un
operator de transformare. [3]
𝑦𝑖𝑘
𝑦𝑖𝑘0
𝑦𝑘𝑖0

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

9

Deoarece în cazul sistemului " TEST 2 " transformatoarele sunt de tip coborâtor,
reprezentarea acestora se va face printr -un cuadripol în Γ

În realizarea acestei scheme am ținut cont de următoarele reguli:
• Admitanța transversală ce modelează pierderile în fier ( 𝑦𝑖𝑘0 ) se va reprezenta
întotdeauna pe înfășurarea primară, indiferent de tipul transformatorului. [1]
• Mai întâi, parametrii transformatorului sunt raportați la tensiunea înfășurării fixe (
cea cu tensiune mai mic ă ). [1]
Relații de calcul pentru parametrii transformatoarelor:
• Rezistența transformatorului
𝑅𝑇=∆𝑃𝑠𝑐𝑛𝑜𝑚𝑈𝑛𝑓2
𝑆𝑛𝑇2 [Ω] (1.2.4)
∆𝑃𝑠𝑐𝑛𝑜𝑚 [𝑘𝑊] – pierderile de putere activa la incercarea de scurtcircuit a
transformatorului
𝑈𝑛𝑓 [𝑘𝑉] – tensiunea nominală a transformatorului pe înfășurarea fixă
𝑆𝑛𝑇 [𝑀𝑉𝐴] – puterea aparentă nominală a transformatorului
• Modulul impedanței longitu dinale a transformatorului
𝑍𝑇=𝑢𝑠𝑐
100𝑈𝑛𝑓2
𝑠𝑛𝑡 [Ω] (1.2.5)
𝑢𝑠𝑐 [%] – tensiunea de scurtcircuit a transformatorului
• Reactanța transformatorului
𝑋𝑇=√𝑍𝑇2−𝑅𝑇2 [Ω] (1.2.6)
• Impedanța longitudinală a transformatorului
𝑍𝑇=𝑅𝑇+𝑗𝑋𝑇 [Ω] (1.2.7)
• Admitanța longitudinală a transformatorului in unități relative
𝑦𝑇=1
𝑍𝑇𝑍𝑏 [𝑢.𝑟.] (1.2.8)
• Conductanța transversal ă a transformatorului
𝑦𝑖𝑘
𝑦𝑖𝑘0
𝑁𝑖𝑘
𝑖
𝑘
Fig 1.2. 2 Schemă echivalentă a unui transformator coborâtor în Γ

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

10

𝐺𝑇0=∆𝑃0
𝑈𝑛𝑓2 [𝜇𝑆] (1.2.9)
∆𝑃0 [𝑘𝑊] – pierderile de putere la încercarea de mers în gol a transformatorului
• Modulul admitanței transversale a transformatorului
𝑌𝑇0=𝑖0
100𝑆𝑛𝑇
𝑈𝑛𝑓 [𝜇𝑆] (1.2.10)
𝑖0 [%] – curentul nominal la mersul in gol
• Susceptanța inductiva a transformatorului
𝐵𝑇0=√𝑌𝑇02−𝐵𝑇02 [𝜇𝑆] (1.2.11)
• Admitanța transversal ă a transformatorului
𝑌𝑇0=𝐺𝑇0−𝑗𝐵𝑇0 [𝜇𝑆] (1.2.12)
• Admitanța transversal ă a transformatorului in unități relative
𝑦𝑇0=𝑌𝑇0𝑍𝑏 [𝑢.𝑟.] (1.2.13)
• Raportul de transformare
𝑁𝑇=𝑈𝑛𝑓
𝑈𝑛𝑟(1+𝑛𝑝∆𝑢
100)𝑈𝑏1
𝑈𝑏2 (1.2.14)

1.3. Determinarea mărimilor electrice de stare în nod – Ecuația puterilor
nodale
Luăm în considera re un nod al REȚ care este conectat cu nodurile învecinate prin laturi
de legătura ( linii și transformatoare electrice ) și în același nod mai avem un generator electric
care injectează putere în rețea și un consumator care este reprezentat în acest caz d e rețeaua de
distribuție. [1]

Puterea nodală netă complexă este diferența dintre puterea injectată în nod și cea
transmisă către consumator :
𝑆𝑖=𝑃𝑖+𝑗𝑄𝑖=(𝑃𝑔𝑖−𝑃𝑐𝑖)+[𝑄𝑔𝑖−𝑄𝑐𝑖]=𝑆𝑔𝑖−𝑆𝑐𝑖 (1.3.1)
Sau mai poate fi scrisă ca produsul dintre tensiunea complexă nodală și curentul injectat
în nod :
𝑆𝑖=𝑈𝑖𝐼𝑖∗ (1.3.2)
Puterea care circula pe laturi se calculează cu formula:
Legături cu
nodurile vecine
𝑃𝑐𝑖,𝑄𝑐𝑖
𝑃𝑔𝑖,𝑄𝑔𝑖
𝑃𝑖,𝑄𝑖
𝑖
𝑘
𝑆𝑡,𝑖
Fig 1. 3.1 . Distribuția puterilor in nod

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

11

𝑆𝑡,𝑖=𝑈𝑖𝐼𝑡,𝑖∗ (1.3.3)
𝐼𝑡,𝑖∗ – suma algebric ă a curenților care circul ă pe laturile incidente nodului i
Ținând cont de relația Kirchhoff I vom scrie relația curentului nodal sub forma
matriceal ă:
[𝐼𝑛]=[𝑌𝑛𝑛][𝑈𝑛] (1.3.4)
𝐼𝑖=∑𝑦𝑖𝑘𝑈𝑘𝑛
𝑘=1 (1.3.5)
De unde rezult ă puterea nodal ă netă:
𝑆𝑖=𝑈𝑖∑𝑦𝑖𝑘∗𝑈𝑘∗ 𝑛
𝑘=1 (1.3.6)
Puterea activ ă și reactiv ă se calculează cu ajutorul următoarelor relații :
𝑃𝑖=∑𝑈𝑖𝑈𝑘[𝐺𝑖𝑘cos(𝜃𝑖−𝜃𝑘)+𝐵𝑖𝑘sin(𝜃𝑖−𝜃𝑘)]𝑛
1 (1.3.7)
𝑄𝑖=∑𝑈𝑖𝑈𝑘[𝐺𝑖𝑘sin(𝜃𝑖−𝜃𝑘)−𝐵𝑖𝑘cos(𝜃𝑖−𝜃𝑘)]𝑛
1 (1.3.8)

1.4. Metoda Newton -Raphson in coordonate polare
Din punct de vedere matematic metoda Newton -Raphson este una din cele mai eficiente
metode pentru rezolvarea ecuaț iilor neliniare. [1]

Principiul metodei
Fie, fi(x)=0,i=1,n un sistem de ecuații neliniare cu n necunoscute. Soluția unui
astfel de sistem se realizează in mod iterativ pornind de la o estimare inițială a soluției .[1]

x(0)=[ x1(0),x2(0)……… xi(0)……… xn(0)]T (1.4.1)
∆x(0)=[ ∆x1(0),∆ x2(0)……… ∆xi(0)……… ∆xn(0)]T (1.4.2)
𝑥(1)=𝑥(0)+∆x(0) (1.4.3 )
𝑥(1) = Solutia sistemului de ecuații (1.4.1 )
𝑓𝑖(x(0)+∆x(0))=𝑓𝑖(𝑥(1))=0 (1.4.4 )
Dezvoltam in serie Taylor relația (1.4.4 ) și reținem termenii de ordinul întâi
𝑓𝑖(x(0)+∆x(0))=𝑓𝑖(𝑥(0))+𝜕𝑓𝑖
𝜕𝑥1∆𝑥1……+𝜕𝑓𝑖
𝜕𝑥𝑛∆𝑥𝑛=0 𝑖=1,𝑛 (1.4.5 )
𝑓(x(0))+𝐽∆x(0)=0 (1.4.5’ )
𝐽=[𝜕𝑓𝑖
𝜕𝑥1𝜕𝑓𝑖
𝜕𝑥𝑛
𝜕𝑓𝑛
𝜕𝑥1𝜕𝑓𝑛
𝜕𝑥𝑛] (1.4.7)
Sistemul de ecuații (1.4.5’ ) se rezolva prin metoda eliminării în raport cu necunoscuta
∆x(0) și se obține o aproximare a solutiei. [1]
𝑥(1)=𝑥(0)+∆x(0) (1.4.8 )
Procesul iterativ se desfășoară astfel:
1) Se calculează 𝐽(𝑝),𝑓(𝑥(𝑝))
2) Se rezolv ă sistemul −𝐽∆x(p)=𝑓(𝑥(𝑝))

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

12

3) Se corecteaz ă soluția 𝑥(𝑝+1)= 𝑥(𝑝)+∆x(p)
Procesul se oprește atunci când eroarea corecțiilor ∆𝑥 este mai mic ă sau egal ă cu eroarea
impus ă.[1]
max |∆x(n)|≤𝜀 (1.4.9)
Pentru aplicarea metodei Newton -Raphson la calculul SEE se vor avea in vedere
următoarele aspecte.
1) Starea de regim permanent este complet definit ă de tensiunile nodale complexe
respectiv de modulele și argumentele acestora. [1]
𝑉𝑆=[𝑈1,𝑈2,……𝑈𝑛,𝜃1,𝜃2,……𝜃𝑛]𝑇 (1.4.10)
2) Având în vedere gruparea nodurilor în noduri de tip PU,PQ și echilibru rezultă
că dintre cele "n" componente ale vectorului de stare VS , necunoscute sunt
argumentele tensiunilor .[1]
𝜃=[𝜃1,𝜃2,……𝜃𝑛−1]𝑇 (1.4.11)
În nodul de echilibru 𝜃𝑛𝑒=𝜃𝑛=0 deoarece fazorul tensiunii în acest nod este ales ca
referință .[1]
𝑈=[ 𝑈1,𝑈2,……𝑈𝑛𝑃𝑄]𝑇 (1.4.12)
Prin urmare vectorul de stare va fi c onstituit din argumentele tensiunilor în toate
nodurile mai puțin nodul de echilibru și modulele tensiunilor din nodurile de tip PQ. [1]
𝑋=[ 𝜃,𝑈]𝑇 (1.4.13)
Dimensiunea acestui vector este 𝑛−1+𝑛𝑃𝑄, pentru a determina necunoscutele este
nevoie de un sistem de ecuații în număr de 𝑛−1+𝑛𝑃𝑄, ecuațiile utilizate sunt [1]:
• Ecuația bilanțului de puteri active in toate nodurile mai puțin nodul de echilibru
∆𝑃𝑖≜𝑃𝑖𝑚𝑝𝑖−𝑃𝑐𝑎𝑙𝑐𝑖=(𝑃𝑔𝑖+𝑃𝑐𝑖)−∑𝑈𝑖𝑈𝑘[𝑛
𝑘=1𝐺𝑖𝑘cos(𝜃𝑖−𝜃𝑘)+ 𝐵𝑖𝑘sin(𝜃𝑖−𝜃𝑘)]
=𝑓𝑃,𝑖(𝜃,𝑈)=0,𝑖=1,2…𝑛−1 (1.4.14 .a)
Aceste ecuații se asociază cu necunoscutele 𝜃1 în număr de n-1
• Ecuația bilanțului de puteri reactive în toate nodurile de tip PQ
∆𝑄𝑖≜𝑄𝑖𝑚𝑝𝑖−𝑄𝑐𝑎𝑙𝑐𝑖=(𝑄𝑔𝑖+𝑄𝑐𝑖)−∑𝑈𝑖𝑈𝑘[𝑛
𝑘=1𝐺𝑖𝑘sin(𝜃𝑖−𝜃𝑘)−𝐵𝑖𝑘cos(𝜃𝑖−𝜃𝑘)]=
𝑓𝑄,𝑖(𝜃,𝑈)=0,𝑖=1,2…𝑛𝑃𝑄 (1.4.14 .b)
Ecuațiile (1.4.14 .a) si ( 1.4.14 .b) reprezintă sistemul de ecuații ce trebuie rezolvat pentru
a determina componentele necunoscute X ale vectorului de stare VS. [1]
Din cele 2n ecuații de bilanț de putere s -au utilizat:
• n-1 ecuații de bilanț de putere activ ă
• 𝑛𝑃𝑄 ecuații de bilanț de putere reactiv ă

Aplicând metoda Newton -Raphson și ținând cont de partiționarea necunoscutelor
obținem:
−𝐽[∆𝜃
∆𝑈]=[∆𝑃
∆𝑄] (1.4.15 )
J – matricea Jacobian
∆𝑃, ∆𝑄 – valoarea funcțiilor 𝑓𝑃,𝑖 si 𝑓𝑄,𝑖 calculate folosind estimarea curentă și
reprezentând abaterile de putere în nodurile rețelei și abaterile care trebuie să fie 0

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

13

Mai departe vom deriva termenii matricei J in raport cu modulele tensiunii și se
înmulțesc cu valoarea tensiunii și în același timp corecțiile ∆𝑈 se impart la aceeași valoare, în
urma acestor calcule de simplificare, sistemul ( 1.4.15 ) se scrie sub forma [1]:
[𝐻𝐾
𝑀𝐿][∆𝜃
∆𝑈
𝑈]=[∆𝑃
∆𝑄] (1.4.16 )
Este necesar să calculăm derivatele funcțiilor 𝑓𝑃,𝑖 si 𝑓𝑄,𝑖 în raport cu modulele și
argumentele tensiunii pentru a determina expresiile submatricelor H,K,M,L. [1]
Pentru a stabilii aceste expresii scriem (6.a) și (6.b) sub forma:
𝑓𝑃,𝑖(𝜃,𝑈)=𝑃𝑖𝑚𝑝𝑖−𝑃𝑐𝑎𝑙𝑐𝑖(𝜃,𝑈)=0
𝑃𝑐𝑎𝑙𝑐𝑖=𝑃𝑡𝑖=𝐺𝑖𝑖𝑈𝑖2+∑𝑈𝑖𝑈𝑘[𝐺𝑖𝑘cos(𝜃𝑖−𝜃𝑘)+𝐵𝑖𝑘sin(𝜃𝑖−𝜃𝑘)]𝑛
𝑘=1
𝑘≠𝑖 (1.4.17 .a)
𝑓𝑄,𝑖(𝜃,𝑈)=𝑄𝑖𝑚𝑝𝑖−𝑄𝑐𝑎𝑙𝑐𝑖(𝜃,𝑈)=0
𝑄𝑐𝑎𝑙𝑐𝑖=𝑄𝑡𝑖=−𝐵𝑖𝑖𝑈𝑖2+∑𝑈𝑖𝑈𝑘[𝐺𝑖𝑘sin(𝜃𝑖−𝜃𝑘)−𝐵𝑖𝑘cos(𝜃𝑖−𝜃𝑘)]𝑛
𝑘=1
𝑘≠𝑖 (1.4.17 .b)
Calculul termenilor submatricelor H,K,M,L
• Se va ține cont de faptul c ă 𝑃𝑖𝑚𝑝𝑖 si 𝑄𝑖𝑚𝑝𝑖 sunt constante
𝐻𝑖𝑖≜−𝜕𝑓𝑃𝑖
𝜕𝜃𝑖=𝜕𝑃𝑐𝑎𝑙𝑐𝑖
𝜕𝜃𝑖=∑𝑈𝑖𝑈𝑘[−𝐺𝑖𝑘sin(𝜃𝑖−𝜃𝑘)+𝐵𝑖𝑘cos(𝜃𝑖−𝜃𝑘)]𝑛
𝑘=1
𝑘≠𝑖
𝐻𝑖𝑖=−𝑄𝑐𝑎𝑙𝑐𝑖−𝐵𝑖𝑖𝑈𝑖2 (1.4.18 )
𝐻𝑖𝑘≜−𝜕𝑓𝑃𝑖
𝜕𝜃𝑘=𝜕𝑃𝑐𝑎𝑙𝑐𝑖
𝜕𝑄𝑘=𝑈𝑖𝑈𝑘[𝐺𝑖𝑘sin(𝜃𝑖−𝜃𝑘)−𝐵𝑖𝑘cos(𝜃𝑖−𝜃𝑘)] (1.4.19 )
𝐾𝑖𝑖≜−𝜕𝑓𝑃𝑖
𝜕𝑈𝑖𝑈𝑖=𝜕𝑃𝑐𝑎𝑙𝑐𝑖
𝜕𝑈𝑖𝑈𝑖
=[2𝐺𝑖𝑖𝑈𝑖+∑𝑈𝑘[𝐺𝑖𝑘cos(𝜃𝑖−𝜃𝑘)+𝐵𝑖𝑘sin(𝜃𝑖−𝜃𝑘)]𝑛
𝑘=1
𝑘≠𝑖]𝑈𝑖

𝐾𝑖𝑖=𝑃𝑐𝑎𝑙𝑐𝑖+𝐺𝑖𝑖𝑈𝑖2 (1.4.20 )
𝐾𝑖𝑘≜−𝜕𝑓𝑃𝑖
𝜕𝑈𝑘𝑈𝑘=𝑈𝑖𝑈𝑘[𝐺𝑖𝑘cos(𝜃𝑖−𝜃𝑘)+𝐵𝑖𝑘sin(𝜃𝑖−𝜃𝑘)] (1.4.21 )
𝑀𝑖𝑖≜−𝜕𝑓𝑄𝑖
𝜕𝜃𝑖=𝜕𝑄𝑐𝑎𝑙𝑐𝑖
𝜕𝜃𝑖=∑𝑈𝑖𝑈𝑘[𝐺𝑖𝑘cos(𝜃𝑖−𝜃𝑘)+𝐵𝑖𝑘sin(𝜃𝑖−𝜃𝑘)]𝑛
𝑘=1
𝑘≠𝑖

𝑀𝑖𝑖=−𝐺𝑖𝑖𝑈𝑖2+𝑃𝑐𝑎𝑙𝑐𝑖 (1.4.22 )
𝑀𝑖𝑘≜−𝜕𝑓𝑄𝑖
𝜕𝜃𝑘=𝜕𝑄𝑐𝑎𝑙𝑐𝑖
𝜕𝜃𝑘=𝑈𝑖𝑈𝑘[−𝐺𝑖𝑘cos(𝜃𝑖−𝜃𝑘)−𝐵𝑖𝑘sin(𝜃𝑖−𝜃𝑘)]
𝑀𝑖𝑘=−𝑈𝑖𝑈𝑘[𝐺𝑖𝑘cos(𝜃𝑖−𝜃𝑘)+𝐵𝑖𝑘sin(𝜃𝑖−𝜃𝑘)] (1.4.23 )
𝐿𝑖𝑖≜−𝜕𝑓𝑄𝑖
𝜕𝑈𝑖𝑈𝑖=𝜕𝑄𝑐𝑎𝑙𝑐𝑖
𝜕𝑈𝑖𝑈𝑖
=[−2𝐵𝑖𝑖𝑈𝑖+∑𝑈𝑘[𝐺𝑖𝑘sin(𝜃𝑖−𝜃𝑘)−𝐵𝑖𝑘cos(𝜃𝑖−𝜃𝑘)]𝑛
𝑘=1
𝑘≠𝑖]𝑈𝑖
𝐿𝑖𝑖=−𝐵𝑖𝑖𝑈𝑖2+𝑄𝑐𝑎𝑙𝑐𝑖 (1.4.24 )

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

14

𝐿𝑖𝑘≜−𝜕𝑓𝑄𝑖
𝜕𝑈𝑘𝑈𝑘=𝑈𝑘[𝑈𝑖[𝐺𝑖𝑘sin(𝜃𝑖−𝜃𝑘)−𝐵𝑖𝑘cos(𝜃𝑖−𝜃𝑘)]]
𝐿𝑖𝑘=𝑈𝑖𝑈𝑘[𝐺𝑖𝑘sin(𝜃𝑖−𝜃𝑘)−𝐵𝑖𝑘cos(𝜃𝑖−𝜃𝑘)] (1.4.25 )
Se constată faptul c ă 𝑀𝑖𝑘=−𝐾𝑖𝑘 iar 𝐿𝑖𝑘=𝐻𝑖𝑘 ceea ce va ajuta la reducerea efortului
de calcul în etapa de formare a matricei Jacobian [1]. Calculul se efectueaza efectuează în mod
iterativ pornind de la un set de valori inițiale pentru modulele tensiunilor și argumentele
acestora, alese astfel :
{𝜃1(0)=𝜃2(0)=⋯𝜃𝑛−1(0)=𝜃𝑛=0
𝑈1(0)=𝑈1(0)=⋯𝑈1(0)=1
𝑈𝑛𝑃𝑄+1=𝑈𝑖𝑚𝑝𝑛𝑃𝑄+1=⋯𝑈𝑛=𝑈𝑖𝑚𝑝𝑛
Valorile argumentelor tensiunilor la pasul inițial sunt egale cu 0, adică sunt egale cu
argumentele tensiunilor în nodul de echilibru.
Modulele tensiunii în nodurile de ti pul PQ se aleg egale cu tensiunea nominală a
nodurilor respective ( adică 1 în unități relative ).
Modulele tensiunii în nodurile de tipul PU și în nodul de echilibru sunt deja cunoscute,
valorile acestora fiind egale cu valoarea impusă.
Fiind cunoscute expresiile inițiale:
𝑋(0)=[𝜃(0),𝑈(0)]𝑇
Se vor calcula ∆𝑃, ∆𝑄 cu ajutorul matricei Jacobian. Pentru asta se rezolv ă următorul
sistem de ecua ții:
[𝐻𝐾
𝑀𝐿][∆𝜃
∆𝑈
𝑈]=[∆𝑃
∆𝑄]
Din acestea rezult ă vectorul curenților
∆𝑋(0)=[∆𝜃(0),∆𝑈(0)
𝑈(0)]𝑇
(1.4.26)
Noua soluție îmbunătățită va fi:
𝜃(1)=𝜃(0)+∆𝜃(0) (1.4.27)
𝑈(1)=𝑈(0)+(∆𝑈(0)
𝑈(0)+1) (1.4.28)
Procesul de calcul se repet ă iterati v pan ă când abaterile devin acceptabile:
𝑚𝑎𝑥{|∆𝑃|,|∆𝑄|}≤𝜀 (1.4.29)
Algoritmul Newton -Raphson nu poate fi utilizat dacă trebuie să se țină cont de limitele
impuse puterii reactive la nodurile de tip PU. Pentru a face posibilă folosirea acestuia trebuie
adăugată o procedură de tratare a nodurilor de tip generator [1].
Pentru aceasta, în cadrul fiecărei iterații după determinarea argumentelor și modulelor
tensiunilor nodale pentru fiecare nod de tipul PU, se calcul ează puterea reactivă [1]:
𝑄𝑐𝑎𝑙𝑐𝑖(𝑝)=−𝐵𝑖𝑖𝑈𝑖+∑𝑈𝑖𝑈𝑘[𝐺𝑖𝑘sin(𝜃𝑖−𝜃𝑘)−𝐵𝑖𝑘cos(𝜃𝑖−𝜃𝑘)]𝑛
𝑘=1
𝑘≠𝑖 (1.4.30 )
După determinarea puterii reactive se verifică dacă aceasta se încadrează în limitele
admisibile. În funcție de valoarea obț inută din relația (1 .4.30 ) putem deosebii 3 cazuri [1]:
a) 𝑄𝑖𝑚𝑖𝑛≤𝑄𝑐𝑎𝑙𝑐𝑖(𝑝)≤𝑄𝑖𝑚𝑎𝑥 rezult ă că nodul rămâne de tip PU
b) 𝑄𝑐𝑎𝑙𝑐𝑖(𝑝)<𝑄𝑖𝑚𝑖𝑛 rezult ă că pentru iterația următoare nodul este trecut in categoria
nodurilor de tip PQ iar 𝑄𝑖𝑚𝑝𝑖=𝑄𝑖𝑚𝑖𝑛
c) 𝑄𝑐𝑎𝑙𝑐𝑖(𝑝)>𝑄𝑖𝑚𝑎𝑥 rezult ă că pentru iterația următoare nodul este trecut in categoria
nodurilor de tip PQ iar 𝑄𝑖𝑚𝑝𝑖=𝑄𝑖𝑚𝑎𝑥

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

15

In iterația următoare înafara de verificarea încadrării în interval, este necesar ca pentru
nodurile declarate inițial de tip PU si care au fost trecute în categoria nodurilor de tip PQ s ă se
verifice îndeplinirea condițiilor de revenire la nod de tip PU. [1]
Pentru a decide dacă se revine sau nu, se ține cont de faptul că între puterea reactivă și
tensiune există o dependentă directă, adică :
• Dac ă 𝑄↑ atunci 𝑈↑
• Dac ă 𝑄↓ atunci 𝑈↓
Ținând cont de acestea rezult ă că:
a) Dac ă nodul a fost trecut în categoria PQ cu 𝑄𝑖𝑚𝑝𝑖=𝑄𝑖𝑚𝑖𝑛 si rezult ă 𝑈𝑖<𝑈𝑖𝑖𝑚𝑝
atunc i pentru a creste tensiunea la valoarea impus ă, este necesar ă creșterea lui Q
de unde rezult ă revenirea nodului in categoria PU
b) Dacă nodul a fost trecut în categoria PQ cu 𝑄𝑖𝑚𝑝𝑖=𝑄𝑖𝑚𝑎𝑥 și rezultă 𝑈𝑖>𝑈𝑖𝑖𝑚𝑝
atunci pentru a scădea tensiunea la valoarea impusă, este necesară scăderea lui
Q de unde rezultă revenirea nodului în categoria PU
Prin utilizarea îndelungată a metodei Newton -Raphson s -a demonstrat că pentru a
asigura o bună convergență, procedura de tratare a nodurilor PQ sau PU nu trebu ie aplicată în
primele 2 iterații. [1]
În cazul sistemelor electro energetice complexe pentru a asigura o bună convergență,
procedura de tratare a nodurilor de tip PU se folosește dacă s -a asigurat o anume abatere a
valorii puterii reactive:
𝑚𝑎𝑥{|∆𝑄𝑖|}≤𝜀1, 𝜀1>𝜀−𝑝𝑟𝑒𝑐𝑖𝑧𝑖𝑎 𝑑𝑒 𝑐𝑎𝑙𝑐𝑢𝑙
În condiții normale de funcționare, calculul regimului permanent prin metoda Newton –
Raphson necesită între 4 și 10 iterații.
Deșii este o metodă cu un număr redus de iterații, efortul de calcul este unul rid icat
deoarece este necesar să se calculeze matricea Jacobian și să se rezolve sistemul de ecuații
liniare printr -o tehnică de factorizare în cadrul fiecărei iterații. [1]

1.5. Metoda Newton -Raphson in coordonate rectangulare

Pentru a utiliza metoda Newton -Raphson la calculul regimului permanent trebuie
exprimate tensiunile [4]. În cazul variantei rectangulare tensiunea se exprimă sub forma:
𝑈𝑖=𝑈𝑖′+𝑗𝑈𝑖′′ (1.5.1)
Ca urmare puterile nodale se vor scrie sub forma:
𝑃𝑖=𝐺𝑖𝑖(𝑈𝑖′2+𝑈𝑖′′2)+∑𝑈𝑖′(𝐺𝑖𝑘𝑈𝑘′−𝐵𝑖𝑘𝑈𝑘′′)𝑛
𝑘=1
𝑘≠𝑖+∑𝑈𝑖′′(𝐺𝑖𝑘𝑈𝑘′′+𝐵𝑖𝑘𝑈𝑘′)𝑛
𝑘=1
𝑘≠𝑖 (1.5.2)
𝑄𝑖=−𝐵𝑖𝑖(𝑈𝑖′2+𝑈𝑖′′2)+∑𝑈𝑖′′(𝐺𝑖𝑘𝑈𝑘′−𝐵𝑖𝑘𝑈𝑘′′)𝑛
𝑘=1
𝑘≠𝑖−∑𝑈𝑖′(𝐺𝑖𝑘𝑈𝑘′′+𝐵𝑖𝑘𝑈𝑘′)𝑛
𝑘=1
𝑘≠𝑖 (1.5.3)
Plecând de la sistemul de ecuații de forma fi(x)=0,i=1,n bilan țul puterilor la noduri
se scrie sub forma:
𝑃𝑖−𝑃𝑖𝑖𝑚𝑝=∆𝑃𝑖′=0 𝑖=1,2,…𝑛;𝑖≠𝑒 (1.5.4)
𝑄𝑖−𝑄𝑖𝑖𝑚𝑝=∆𝑄𝑖′=0 𝑖=1,2,…𝑛;𝑖≠𝑒 (1.5.5)
Aplicând liniarizarea [−𝐹𝑖(𝑝)]=[𝐽(𝑝)][∆𝑋(𝑝)] sistemul de ecua ții devine [4]:
−∆𝑃𝑖′=∑𝜕𝑃𝑖
𝜕𝑥𝑘𝑛
𝑘=1∆𝑥𝑘 𝑖=1,2,…𝑛;𝑖≠𝑒 (1.5.6)
−∆𝑄𝑖′=∑𝜕𝑄𝑖
𝜕𝑥𝑘𝑛
𝑘=1∆𝑥𝑘 𝑖=1,2,…𝑛;𝑖≠𝑒 (1.5.7)

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

16

Din acestea rezultă abaterile de putere:
∆𝑃𝑖=−∆𝑃𝑖′=𝑃𝑖𝑖𝑚𝑝−𝑃𝑖=∑𝜕𝑃𝑖
𝜕𝑥𝑘𝑛
𝑘=1∆𝑥𝑘 (1.5.8)
∆𝑄𝑖=−∆𝑄𝑖′=𝑄𝑖𝑖𝑚𝑝−𝑄𝑖=∑𝜕𝑄𝑖
𝜕𝑥𝑘𝑛
𝑘=1∆𝑥𝑘 (1.5.9)
Pentru nodurile de tip consumator există două relații pentru abaterile puterilor impuse
în nod [4]:
∆𝑃𝑖=𝜕𝑃𝑖
𝜕𝑈𝑖′𝑈𝑖′+𝜕𝑃𝑖
𝜕𝑈𝑖′∆𝑈𝑖′′+∑𝜕𝑃𝑖
𝜕𝑈𝑘′∆𝑈𝑘′ 𝑛
𝑘=1
𝑘≠𝑖,𝑒+∑𝜕𝑃𝑖
𝜕𝑈𝑘′′∆𝑈𝑘′′ 𝑛
𝑘=1
𝑘≠𝑖,𝑒 (1.5.10)
∆𝑄𝑖=𝜕𝑄𝑖
𝜕𝑈𝑖′𝑈𝑖′+𝜕𝑄𝑖
𝜕𝑈𝑖′∆𝑈𝑖′′+∑𝜕𝑄𝑖
𝜕𝑈𝑘′∆𝑈𝑘′ 𝑛
𝑘=1
𝑘≠𝑖,𝑒+∑𝜕𝑄𝑖
𝜕𝑈𝑘′′∆𝑈𝑘′′ 𝑛
𝑘=1
𝑘≠𝑖,𝑒 (1.5.11)
Pentru nodurile de tip generator vom scrie relația pentru abaterea puterii active [4]:
∆𝑃𝑖=𝜕𝑃𝑖
𝜕𝑈𝑖′𝑈𝑖′+𝜕𝑃𝑖
𝜕𝑈𝑖′∆𝑈𝑖′′+∑𝜕𝑃𝑖
𝜕𝑈𝑘′∆𝑈𝑘′ 𝑛
𝑘=1
𝑘≠𝑖,𝑒+∑𝜕𝑃𝑖
𝜕𝑈𝑘′′∆𝑈𝑘′′ 𝑛
𝑘=1
𝑘≠𝑖,𝑒 (1.5.12)
Iar pentru variația puterii reactive vom introduce o ecuație ce exprimă constanta
modulului tensiunii [4]:
𝑈𝑖2=𝑈𝑖′2+𝑈𝑖′′2=𝑐𝑡 (1.5.13)
Această relație se mai poate scrie și sub forma:
∆(𝑈𝑖2)=𝑈𝑖2−(𝑈𝑖𝑖𝑚𝑝)2 (1.5.14)
Din acestea vom obține sistemul liniar de ecuații , de forma:
[∆𝑃
∆𝑄
∆(𝑈2)]=[𝐽1𝐽2
𝐽3𝐽4
𝐽5𝐽6][∆𝑈′
∆𝑈′′] (1.5.15)
Sau
[∆𝑃
∆𝑄
∆(𝑈2)]=[𝐽][∆𝑈′
∆𝑈′′] (1.5.16)
[𝐽] – matricea Jacobian a sistemului, elementele acesteia se calculează cu relațiile [4]:
𝐽1,𝑖𝑖=𝜕𝑃𝑖
𝜕𝑈𝑖′=2𝑈𝑖′𝐺𝑖𝑖+∑(𝐺𝑖𝑘𝑈𝑘′−𝐵𝑖𝑘𝑈𝑘′′)𝑛
𝑘=1
𝑘≠𝑖 (1.5.17)
𝐽1,𝑖𝑘=𝜕𝑃𝑖
𝜕𝑈𝑘′=𝑈𝑖′𝐺𝑖𝑘+𝑈𝑖′′𝐵𝑖𝑘 (1.5.18)
𝐽2,𝑖𝑖=𝜕𝑃𝑖
𝜕𝑈𝑖′′=2𝑈𝑖′′𝐺𝑖𝑖+∑(𝐺𝑖𝑘𝑈𝑘′′+𝐵𝑖𝑘𝑈𝑘′′)𝑛
𝑘=1
𝑘≠𝑖 (1.5.19)
𝐽2,𝑖𝑘=𝜕𝑃𝑖
𝜕𝑈𝑘′′=−𝑈𝑖′𝐵𝑖𝑘+𝑈𝑖′′𝐺𝑖𝑘 (1.5.20)

𝐽3,𝑖𝑖=𝜕𝑄𝑖
𝜕𝑈𝑖′=−2𝑈𝑖′𝐵𝑖𝑖+∑(𝐺𝑖𝑘𝑈𝑘′′+𝐵𝑖𝑘𝑈𝑘′)𝑛
𝑘=1
𝑘≠𝑖 (1.5.21)

𝐽3,𝑖𝑘=𝜕𝑄𝑖
𝜕𝑈𝑘′=𝑈𝑖′′𝐺𝑖𝑘−𝑈𝑖′𝐵𝑖𝑘 (1.5.22)
𝐽4,𝑖𝑖=𝜕𝑄𝑖
𝜕𝑈𝑖′′=−2𝑈𝑖′′𝐵𝑖𝑖+∑(𝐺𝑖𝑘𝑈𝑘′−𝐵𝑖𝑘𝑈𝑘′′)𝑛
𝑘=1
𝑘≠𝑖 (1.5.23)
𝐽4,𝑖𝑘=𝜕𝑄𝑖
𝜕𝑈𝑘′′=−𝑈𝑖′′𝐵𝑖𝑘−𝑈𝑖′𝐺𝑖𝑘 (1.5.24)
𝐽5,𝑖𝑖=2𝑈𝑖′ (1.5.25)
𝐽5,𝑖𝑘=0 (1.5.26)
𝐽6,𝑖𝑖=2𝑈𝑖′′ (1.5.27)
𝐽6,𝑖𝑘=0 (1.5.28)

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

17

Pentru nodurile de tip consumator se va utiliza vectorul ∆𝑄 cu submatricele 𝐽3 și 𝐽4 iar
pentru nodurile de tip generator se va utiliza vectorul ∆(𝑈2) cu submatricele 𝐽5 și 𝐽6.[4]
Corecțiile de tensiune se vor calcula cu relațiile :
𝑈𝑖′(𝑝+1)=𝑈𝑖′(𝑝)+∆𝑈𝑖′(𝑝) (1.5.29)
𝑈𝑖′′(𝑝+1)=𝑈𝑖′′(𝑝)+∆𝑈𝑖′′(𝑝) (1.5.30)
Iar noul set de valori pentru tensiune nodale cu relația :
𝑈𝑖(𝑝+1)=𝑈𝑖′(𝑝+1)+𝑗𝑈𝑖′′(𝑝+1) (1.5.31)
Procesul iterativ se repetă, calculând noi valori ale puterilor nodale, până când abaterile
∆𝑃 și ∆𝑄 iau valori m ai mici decât eroarea stabilită. [4]

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

18

Capitolul 2
Stabilitatea de tensiune. Aspecte fizice.
Interacțiunea dintre rețeaua electrică și sarcini.
Metoda celei mai mici valori singulare, indicatorul global VSI.
Studii de caz pe sistemul test.

2.1. Stabilitatea de tensiune in SEE
Stabilitatea de tensiune reprezintă capacitatea sistemului electroenergetic de a menține
un nivel de tensiune în limitele admisibile în toate nodurile atât în condiții normale, cât și în
urma unor perturbații.
Se consideră că un sistem electroenergetic ( SE ) se află într -o stare de instabilitate atunci
când o perturbație determină o scădere progresivă și necontrolabilă a nivelului de tensiune într –
un nod, într -o zonă sau în tot sistemul. [5]
Cauza p rincipală a declanșării fenomenului de instabilitate de tensiune este reprezentată
de căderile de tensiune, datorate modificărilor circulației de putere prin elementele inductive
ale rețelei de transport [5] în urma:
• Creșterii consumului de putere care duc e la un deficit de putere reactivă locală sau
zonală [5]
• Incidentelor care slăbesc controlul local de tensiune ( declanșarea anumitor grupuri
generatoare ), slăbesc rețeaua de transport ( declanșarea anumitor linii de transport,
transformatoare sau autotransformatoare ) sau măresc puterea tranzitată prin rețeaua
de transport ( separarea rețelelor electrice ) [5]
• Funcționării defectuoase a reglajului sub sarcină a prizelor transformatoarelor [5]
Stabilitatea de tensiune se clasific ă in:
• Stabilitate de tensiune la mici perturbații
• Stabilitate de tensiune la mari perturbații

2.2. Caracteristicile rețelei de transport

Considerăm cazul în care alimentarea unui consumator se realizează prin intermediul
unei linii electrice scurte. Aplicând teorema a II -a a lui Kirchoff vom obține relația de legătură
dintre tensiunile fază -fază de la capetele liniei [5]:
𝑈1=𝑈2+√3𝑍𝐼=𝑈2+∆𝑈+𝑗𝛿𝑈 (2.2.1 )

Fig 2.2.1. Schema echivalent ă monofazat ă a unei linii scurte ( a) și diagrama fazorial ă (b)

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

19

Dacă alegem ca origine de faz ă tensiunea la capătul receptor 𝑈2=𝑈2 și ținem cont c ă
𝑆2=𝑃2+𝑗𝑄2=𝑆2(cos𝜑+𝑗sin𝜑)=√3𝑈2𝐼∗, atunci rela ția (2.2.1 ) devine [5]:
𝑈1=𝑈2+𝑅𝑃2+𝑋𝑄2
𝑈2+𝑗𝑋𝑃2−𝑅𝑄2
𝑈2 (2.2.1')
𝑅=𝑍cos𝛽 – rezisten ța corespunzătoare impedan ței echivalente 𝑍 a liniei de transport
𝑋=𝑍sin𝛽 – reactan ța corespunzătoare impedan ței echivalente 𝑍 a liniei de transport
Prin proiectarea relației (1') pe axa real ă și axa imaginar ă se obțin ecuațiile :
𝑈1𝑈2cos𝜃=𝑈22+𝑅𝑃2+𝑋𝑄2 (2.2.1'')
𝑈1𝑈2sin𝜃=𝑋𝑃2−𝑅𝑄2 (2.2.1''')
Vom ridica la pătrat și vom însuma cele doua ecuații pentru a elimina variabila 𝜃[5]:
𝑈12𝑈22(cos2𝜃+sin2𝜃)=𝑈24+2𝑈22(𝑅𝑃2+𝑋𝑄2)+(𝑅𝑃2+𝑋𝑄2)2+(𝑋𝑃2−𝑅𝑄2)2
Ținând cont de faptul c ă (𝑅𝑃2+𝑋𝑄2)2+(𝑋𝑃2−𝑅𝑄2)2=𝑍2(𝑃22+𝑄22)=𝑍2𝑆22 iar
sin2𝜃+cos2𝜃=1 vom ob ține:
𝑓(𝑈2,𝑆2,𝜑)=𝑈24+[2(𝑅𝑃2+𝑋𝑄2)−𝑈12]𝑈22+𝑍2𝑆22=0 (2.2.2)
Având in vedere faptul c ă 𝑃2=𝑆2cos𝜑, iar 𝑄2=𝑆2sin𝜑, rela ția (2.2.2) devine:
𝑓(𝑈2,𝑆2,𝜑)=𝑈24+[2𝑆2(𝑅cos𝜑+𝑋sin𝜑)−𝑈12]𝑈22+𝑍2𝑆22=0 (2.2.2')
Relațiile (2.2.2) și (2.2.2') arată legătura dintre tensiunea 𝑈2 la capătul receptor al liniei
electrice de transport și puterile transportate ( 𝑃2 și 𝑄2 sau 𝑆2 și 𝜑 ). Vom exprima valoarea
tensiunii 𝑈2 ca o funcție explicită [5]:
𝑈2=𝑔(𝑃2,𝑄2)=𝑔(𝑆2,𝜑) (2.2.3)
Aceast ă funcție definește caracteristicile rețelei de transport reprezentat ă in figur a 2.2.2

Fig 2.2.2 Relația dintre tensiunea 𝑈2 și puterile 𝑃2 și 𝑄2, pentru o structur ă simpl ă
de tipul generator -linie-consumator
Pentru ca funcția explicit ă (2.2.3) să existe este necesar ca ecuațiile (2.2.2) și (2.2.2') să
admit ă soluții reale [5]. Aceast ă condiție este îndeplinit ă dacă:
[2𝑆2(𝑅cos𝜑+𝑋sin𝜑)−𝑈12]2−4𝑍2𝑆22≥0 (2.2.4)
Dacă impunem condiția de egalitate cu zero in relația (2.2.4), obținem :
2𝑆2(𝑅cos𝜑+𝑋sin𝜑)−𝑈12=±2𝑍𝑆2
sau
2𝑆2(𝑅cos𝜑+𝑋sin𝜑±𝑍)=𝑈12
Ținând cont de faptul c ă 𝑅=𝑍cos𝛽 și 𝑋=𝑍sin𝛽 rezult ă relația :
𝑆2=𝑈12
2(𝑅cos𝜑+𝑋sin𝜑±𝑍)=𝑈12
2𝑍[cos(𝛽−𝜑)±1] (2.2.5)

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

20

Având in vedere faptul c ă 𝑆2≥0 rezult ă din relația (2.2.5) că valoarea puterii aparente
maxime ce poate fi tranzitat ă către zona de consum este [5]:
𝑆2𝑚𝑎𝑥=𝑈12
2𝑍[cos(𝛽−𝜑)+1]=𝑈12
4𝑍cos2(𝛽−𝜑
2) (2.2.6)
Iar puterea maxim ă se calculează cu relația [5]:
𝑃2𝑚𝑎𝑥=𝑆2𝑚𝑎𝑥cos𝜑=𝑈12cos𝜑
4𝑍cos2(𝛽−𝜑
2) (2.2.7)

2.3. Aspecte statice ale stabilității de tensiune. Puncte și zone de
funcționare.

Pentru a putea face analiza punctelor din zona hașurată din figur a 2.2.2 , vom ține cont
de relația de legătura dintre 𝑃2 și 𝑄2, cu ajutorul factorului de putere cos𝜑, de unde rezultă
ecuația [5]:
𝑈24+[2𝑃2(𝑅+𝑋tan𝜑)−𝑈12]𝑈22+𝑍2(𝑃2
cos𝜑)2
=0 (2.3.1)
Dacă consider ăm 𝑈1 și cos𝜑 constante, putem determina o rela ție de forma 𝑈2=𝑓(𝑃2)
pentru 𝑃2∈[0,𝑃2𝑚𝑎𝑥], aceast ă relație este denumit ă caracteristica 𝑈2−𝑃2 a re țelei
electrice [5].

Fig 2.3.1 Caracteristica 𝑈2−𝑃2. Puncte și zone de funcționare .
Din analiza graficului, putem observa c ă pentru o valoare 𝑃2<𝑃2𝑚𝑎𝑥 exist ă două
puncte de func ționare:
• Punctul A – este caracterizat de o tensiune ridicată și o funcționare normală a
sistemului, de unde rezultă că este un punct stabil de funcționare [5]
• Punctul B – este caracterizat de o valoare a tensiunii mai mică decât valoarea critică
admisă 𝑈2𝑐𝑟𝑖𝑡 și de o funcționare anormală a sistemului, acesta fiind un punct instabil
de funcționare [5]
Pentru a putea definii mai bine zonele de funcționare, vom ține cont de efectele
reglajului de tensiune la capătul receptor al liniei de transport :
• Efectul compensării puterii reactive transportate [5]

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

21

Prin menținerea valorii tensiunii la capătul sursă al l iniei de transport 𝑈1=𝑐𝑡 și prin
variația valorii lui cos𝜑 vom ob ține urm ătoarele caracteris tici 𝑈2−𝑃2[5]:

Fig 2.3.2 Efectul compensării puterii reactive asupra caracteristicii 𝑈2−𝑃2
Din analiza graficului se observă că în cazul punctului de funcționare A valoarea
tensiunii la capătul receptor al liniei de transport ( 𝑈2 ) crește atunci când factorul de putere
trece de la o valoare inductivă (cos𝜑𝑖𝑛𝑑) la o valoare capacitivă (cos𝜑𝑐𝑎𝑝). Pentru punctul de
funcționare B efectul compensării puterii reactive este invers față de cel în cazul punctului A,
valoarea tensiunii scăzând atunci când factorul de putere trece de la o valoare inductivă la una
capacitivă .
• Efectul modificării tensiunii la capătul sursă al liniei de transport [5]
De data aceasta vom menține valoarea factorului de putere constantă (cos𝜑=𝑐𝑡) și
vom da valori diferite tensiunii de la capătul sursă (𝑈1) al liniei de tran sport ( prin modificarea
excitației generatorului ) pentru a obține o noua caracteristică 𝑈2−𝑃2[5]

Fig 2.3.3 Efectul modificării tensiunii 𝑈1 asupra caracteristicii 𝑈2−𝑃2

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

22

Se poate observa că și în acest caz, prin modificarea valorii tensiunii la capătul sursă,
valoarea tensiunii corespunzătoare punctului de funcționare A crește odată cu creșterea tensiunii
sursei în timp ce tensiunea corespunzătoare punctului B scade.
În concluzie prin reglarea tensiunii la capătul sursă sau a factorului de putere, vom
obține efectul dorit (de creștere al tensiunii) doar în punctul de funcționare A, în timp ce prin
aplicarea acelorași măsuri punctului B obținem un efect complet opus (de scădere al tensiunii).
Pentru a înțelege mai bine aceste fenomene vom ține cont de sensibilitățile tensiunii la
variația puterilor activă și reactivă. Pentru determinarea acestora vom pornii de la relația [5]:

𝑈22−𝑈1𝑈2+𝑅𝑃2+𝑋𝑄2=0 (2.3.2)
Din această relație vom scoate valorile tensiunilor punctelor de funcționare A și B [5]:
𝑈2𝐴≅𝑈1+√𝑈12−4(𝑅𝑃2+𝑋𝑄2)
2 (2.3.3)
𝑈2𝐵≅𝑈1−√𝑈12−4(𝑅𝑃2+𝑋𝑄2)
2 (2.3.3’)
Exprimăm condiția de controlabilitate a sistemului electroenergetic prin relațiile [5]:
{ 𝑑𝑈𝑖
𝑑𝐺𝑖|
𝐵𝑖=𝑐𝑡<0
𝑑𝑃𝑖
𝑑𝐺𝑖|
𝐵𝑖=𝑐𝑡>0
𝑑𝑈𝑖
𝑑𝐵𝑖|
𝐺𝑖=𝑐𝑡<0
𝑑𝑄𝑖
𝑑𝐵𝑖|
𝐺𝑖=𝑐𝑡>0
Ținând cont de aceste condiții și de relațiile precedente, putem calcula sensibilitățile
tensiunii în raport cu puterea activă și reactivă în cele două puncte [5]:
𝑑𝑈2𝐴
𝑑𝑃2|
𝑄2=𝑐𝑡=−𝑅
√𝑈12−4(𝑅𝑃2+𝑋𝑄2)<0
𝑑𝑈2𝐴
𝑑𝑄2|
𝑃2=𝑐𝑡=−𝑋
√𝑈12−4(𝑅𝑃2+𝑋𝑄2)<0
𝑑𝑈2𝐵
𝑑𝑃2|
𝑄2=𝑐𝑡=𝑅
√𝑈12−4(𝑅𝑃2+𝑋𝑄2)>0
𝑑𝑈2𝐵
𝑑𝑄2|
𝑃2=𝑐𝑡=𝑋
√𝑈12−4(𝑅𝑃2+𝑋𝑄2)>0

Având în vedere cele menționate anterior, se pot definii următoarele zone de
funcționare :
• Zona controlabil ă ( zona stabil ă de funcționare ) – Aceasta este definit ă de punctele
de funcționare pentru care 𝑍𝑐>𝑍, iar sensibilitățile date de relațiile sunt
negative [5].
• Zona de securitate – aceasta este definit ă între valorarea maxim ă și minim ă de
tensiune si este caracterizat ă de valori negative foarte sc ăzute ( apropiate de zero )
ale sensibilităților [5].
• Zona critic ă – este definit ă intre valorile 𝑈𝑚𝑖𝑛 si 𝑈𝑐𝑟 ale tensiunii și este caracterizat ă
la fel ca zona de securitate prin valori negative ale sensibilităților dar mari in valori

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

23

absolute. O func ționare îndelungată in aceast ă zonă este periculoas ă deoarace poate
sa apar ă fenomenul de instabilitate datorat unei cereri suplimentare de putere sau
acționare a transformatoarelor cu reglaj sub sarcin ă[5].
• Zona necontrolabil ă (zona instabil ă de funcționare ) – aceasta este definit ă de
punctele de funcționare pentru care 𝑍𝑐<𝑍, iar sensibilitățile au valori negative [5].

2.4. Metoda celei mai mici valori singulare. Indicatorul global VSI

Deoarece sistemul electroenergetic are o complexitate foarte ridicată, această metodă
are scopul de a determina punctele vulnerabile ale sistemului și de a identifica măsurile ce sunt
necesare a fi luate pentru a crește siguranța în exploatare a acestuia [5].
Pentru a calcula valoarea indicatorului VSI vom pornii prin a liniariza sistemul de
ecuații [5]:
∑𝑈𝑖𝑈𝑘[𝐺𝑖𝑘cos(𝜃𝑖−𝜃𝑘)+𝐵𝑖𝑘sin(𝜃𝑖−𝜃𝑘)]−𝑃𝑖𝑛
𝑘=1 =0 𝑖=2,3,…𝑛 (2.4.1)
∑𝑈𝑖𝑈𝑘[𝐺𝑖𝑘sin(𝜃𝑖−𝜃𝑘)−𝐵𝑖𝑘cos(𝜃𝑖−𝜃𝑘)]−𝑄𝑖𝑛
𝑘=1 =0 𝑖=𝑛𝑔+1,…𝑛 (2.4.2)
Vom obține modelul liniar:
[𝐽(𝑝)][∆𝜃(𝑝)
∆𝑈(𝑝)]=[∆𝑃(𝑝)
∆𝑄(𝑝)] (2.4.3 )
[𝐽(𝑝)]=[𝐽𝑃𝜃𝐽𝑃𝑈
𝐽𝑄𝜃𝐽𝑄𝑈] – Matricea Jacobian calculata cu valorile tensiunilor la pasul (p)
Descompunând după valorile singulare ale matricei Jacobian în jurul punctului de
echilibru analizat vom obține:
[𝐽]=[𝑈][Σ][𝑉]𝑇=∑𝜎𝑖[𝑈𝑖][𝑉𝑖]𝑇𝑛
𝑖=1(2.4.4)

Considerăm că matricea Jacobian nu este singulară și îi calculăm inversa:
[𝐽]−1=[𝑈][Σ]−1[𝑉]𝑇=∑𝜎𝑖−1[𝑈𝑖][𝑉𝑖]𝑇 𝑛
𝑖=1 (2.4.5)
Deci putem determina efectul modificărilor puterilor active și reactive asupra modulelor
și argumentelor tensiunilor cu relația [5]:
[∆𝜃
∆𝑈]=[𝐽]−1[∆𝑃
∆𝑄]=∑𝜎𝑖−1𝑈𝑖𝑉𝑖𝑡[∆𝑃
∆𝑄]𝑛
𝑖=1 (2.4.6)
Daca consider ăm [∆𝑃
∆𝑄]=[𝑉𝑛] și ținând cont de faptul c ă matricea [𝑉] este ortonormat ă
atunci relația devine [5]:
[∆𝜃
∆𝑈]=𝜎𝑛−1[𝑈𝑛] (2.4.7 )
Conform relației ( 2.4.7 ) și a literaturii de specialitate, cele mai mici valori singulare a
matricei Jacobian și a vectorilor asociați , au următoarele semnificații :
• 𝜎𝑚𝑖𝑛(𝐽) este un indicator global al stabilității de tensiune a sistemului
electroenergetic, acesta poate să determine dacă sistemul este sau nu în apropierea
unui punct critic [5].
• Componentele vectorului [𝑈𝑛] indică sensibilitățile modulelor și argumentelor
tensiunilor la variația puterilor active și reactive [5].
• Componentele vectorului [𝑉𝑛] indică cele mai sensibile direcții de variație ale
puterilor [5].

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

24

Ținând cont de cuplajul puternic dintre U și Q în practică se recomandă utiliza rea
matricei Jacobian redusă [5]. Pentru a putea să o calculăm, pornim prin rescrierea relației ( 2.4.3 )
sub forma :
[∆𝑃]=[𝐽𝑃𝜃][∆𝜃]+[𝐽𝑃𝑈][∆𝑈] (2.4.8)
[∆𝑄]=[𝐽𝑄𝜃][∆𝜃]+[𝐽𝑄𝑈][∆𝑈] (2.4.9)
Daca cuplajul U – P este slab atunci vom considera ∆𝑃=0 de unde rezult ă:
[∆𝜃]=−[𝐽𝑃𝜃]−1[𝐽𝑃𝑈][∆𝑈] (2.4.10)
Introducând relația (2.4.10) în (2.4.9 ) obținem :
[∆𝑄]=([𝐽𝑄𝑈]−[𝐽𝑄𝜃][𝐽𝑃𝜃]−1[𝐽𝑃𝑈])[∆𝑈] (2.4.11)
Rezultând matricea Jacobian redus ă:
[𝐽𝑅]=[𝐽𝑄𝑈]−[𝐽𝑄𝜃][𝐽𝑃𝜃]−1[𝐽𝑃𝑈] (2.4.12)
Aceasta a fost obținută în ipoteza că nu există probleme de stabilitate unghiulară în
sistemul analizat. Conform teoremei lui Schur putem trage concluz ia că matricea Jacobian este
determinată de existența unor probleme de stabilitate de tensiune [5].
De unde rezultă faptul că cea mai mică valoare singulară a matricei Jacobian redusă
devine un indicator global al stabilității de tensiune [5].
Indicatorul de stabilitate de tensiune VSI poate fi scris ca raportul dintre cea mai mică
valoare singulară 𝜎𝑚𝑖𝑛(𝐽𝑅) și valoarea minimă singulară pentru regimul de mers in gol
𝜎0𝑚𝑖𝑛(𝐽𝑅)[5]:
𝑉𝑆𝐼=𝜎𝑚𝑖𝑛(𝐽𝑅)
𝜎0𝑚𝑖𝑛(𝐽𝑅) (2.4.13)
Dacă valoarea indicatorului VSI este apropiată de 0 atunci sistemul se află în apropierea
frontierei critice, iar dacă acesta are o valoarea apropiată de 1 sistemul se află departe de
aceasta [5].

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

25

Capitolul 3
Utilizarea RNA de tip MLP pentru monitorizarea riscului declanșării
instabilității de tensiune pentru o zon ă de consum.
Studii de caz pe un sistem test.

3.1. Rețele neuronale artificiale – RNA

Rețelele neuronale artificiale sunt sisteme dinamice de prelucrare a informației, formate
dintr -un număr mare de structuri elementare de prelucrare, puternic interconectate, numite
neuroni [6].

3.1.1. Caracteristicile principale ale unei RNA
Acestea pot prelucra informații în para lel ceea ce ne permite construirea unor arhitecturi
performante care pot oferii răspunsuri rapide, în timp real, la probleme complexe [6].
Învățarea din exemple a fenomenului pentru care se antrenează, procesul de învățare
poate fi supervizat sau nesupervi zat.
Permit modelarea unor procese complexe ale căror legi de funcționare sunt fie prea
complicate, fie nu pot fi algoritmizate [6].
Capacitatea de a asocia datele de intrare cu setul sau seturile de date de antrenare cu
ajutorul cărora poate oferii un răspuns cât mai bun. Această capacitate asigură un comportament
bun chiar în condițiile unor seturi de date de intrare incomplete sau p arțial greșite [6].
În cazul "defectării" anumitor neuroni, memorarea distribuită a informației asigură o
funcționare relativ corectă .
3.1.2. Arhitectura general ă a unei RNA
Rețelele neuronale artificiale sunt formate din neuroni care sunt structuri elementare de
procesare legați prin conexiuni sinaptice organizați pe straturi succesive [6].

Fig 3.1.2.1 Arhitectura unei RNA
Stratul de intrare este alcătuit din neuroni liniari care nu realizează nici un proces de
calcul, rolul acestora fiind de a dirija componentele 𝑥1,……𝑥𝑖,……𝑥𝑛 ale vectorului de intrare
X către neuronii de calcul din primul strat [6].
Stratul de ieșire poate fi format din unul sau mai mulți neuroni, în general liniari, care
furnizează răspunsul corespunzător intrării curente X [6].

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

26

O RNA poate să conțină unul sau mai multe straturi ascunse ( intermediare ) care sunt
constituite din neuroni neliniari (sigmoidali) care formează o reprezentare internă ce poate fi
considerată o memorie asociativă [6].
Legăturile dintre straturile RNA se realizează prin conexiuni sinaptice ponderate. În
cazul rețelei prezentată în figură, fiecare neuron dintr -un strat este legat cu toți neuronii din
stratul următor iar între două straturi neco nsecutive nu există legături [6]. O astfel de rețea
neuronală se numește complet conectată și se numește perceptron multistrat sau Mulți Layer
Perceptron -MLP .
Structura prezentată în figura anterioară are următoarele proprietăți :
• Fiecar e neuron acționează independent de ceilalți neuroni din același strat ( ieșirea
unui neuron este influențată doar de informația primită prin conexiunile sinaptice de
intrare) [6].
• Activarea fiecărui neuron depinde numai de informații cu caracter local. Nu es te
necesară cunoașterea informațiilor neuronilor care nu au legături directe cu neuronul
considerat [6].
• Un număr cât mai mare de conexiuni asigură un grad cât mai ridicat de rezervare și
ușurează distribuirea informației în scopul creșterii vitezei de calc ul[6].

3.2. Multi Layer Perceptron – MLP
Rețelele care fac parte din această categorie pot avea un număr nelimitat de straturi. Toate
straturile cu excepția primului și ultimului sunt straturi ascunse. Principalele aplicații ale acestui
tip de rețea au în vedere probleme de clasificare sau de aproximare euristică a unor funcții
matematice [6].
3.2.1. Antrenarea MLP

Vom considera o rețea neuronală artificială de tip MPL cu un singur strat ascuns :

Fig 3.2.1.1 RNA cu un singur strat ascuns

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

27

Se vor face următoarele notații :
𝑛𝑖 – Num ărul neuronilor de pe stratul de intrare
𝑛𝑗 – Num ărul neuronilor de pe stratul ascuns
𝑛𝑘- Num ărul neuronilor de pe stratul de ieșire
W – matricea ponderilor neuronilor de pe stratul ascuns ale c ărei elemente 𝑤𝑗,𝑖
sunt ponderile conexiunilor sinaptice dintre neuronul j de pe stratul ascuns și neuronu l i
de pe stratul de intrare
V – matricea ponderilor de pe stratul de ieșire ale cărei elemente 𝑣𝑘,𝑗 sunt
ponderile conexiunilor sinaptice dintre neuronul k de pe stratul de ieșire și neuronul j de
pe stratul ascuns
𝑛𝑒 – numărul de exemple ( num ărul d e perechi ) ( 𝑋𝑚,𝑑𝑚) care alc ătuiesc
mulțimea de antrenare
𝑧𝑗(𝑚)=𝑔(𝑛𝑒𝑡𝑗(𝑚)) – ieșirea neuronului j de pe stratul ascuns când la intrare se
prezintă exemplul m
𝑔:𝑅→𝑅 – functia de activare a neuronilor de pe stratul ascuns
𝑦𝑘(𝑚)=𝑓(𝑛𝑒𝑡𝑘(𝑚)) – ieșirea neuronului k de pe stratul de ieșire când la intrare se
prezintă exemplul m
𝑓:𝑅→𝑅 – functia de activare a neuronilor de pe stratul de ie șire
O rețea neuronală de tip MLP poate fi antrenată pentru aproximarea unei funcții
ℎ:𝑅𝑛𝑖→𝑅𝑛𝑘 sau pentru clasificarea formelor de intrare.
Procesul de antrenare constă în modificarea iterativă a ponderilor și pragurilor
neuronilor de pe straturile de prelucrare în scopul minimizării funcției de performanță a
rețelei [6]. Această funcție poate fi:
• abaterea pătratică parțială
𝐸(𝑊,𝑉)=∑(𝑑𝑘(𝑚)−𝑦𝑘(𝑚))2 𝑛𝑘
𝑘=1 (3.1.2.1)
• abaterea pătratică medie parțială
𝐸(𝑊,𝑉)=1
𝑛𝑘∑(𝑑𝑘(𝑚)−𝑦𝑘(𝑚))2 𝑛𝑘
𝑘=1 (3.1.2.2)
• abaterea pătratică totală
𝐸(𝑊,𝑉)=∑∑(𝑑𝑘(𝑚)−𝑦𝑘(𝑚))2 𝑛𝑘
𝑘=1𝑛𝑒
𝑚=1 (3.1.2.3)
• abaterea pătratica medie total ă
𝐸(𝑊,𝑉)=1
𝑛𝑒𝑛𝑘∑∑(𝑑𝑘(𝑚)−𝑦𝑘(𝑚))2 𝑛𝑘
𝑘=1𝑛𝑒
𝑚=1 (3.1.2.4)
Pentru a implementa algoritmul de antrenare putem proceda în următoarele moduri :
• modul incremental în care ponderile sunt corectate după fiecare exemplu, iar funcția
de performanță este abaterea pătratică parțială sau abaterea pătratică medie
parțială [6].
• modul global în care ponderile sunt corectate după ce toate exemplele din mulțimea
de antrenare au fost prezentate rețelei ( la sfârșitul unei epoci ), iar funcția de
performanță este abaterea pătratică totală sau abaterea pătratică medie totală [6].
Teoretic, pentru determinarea matricelor ponderilor W și V ar trebui impuse condițiile :
𝜕𝐸
𝜕𝑤𝑗𝑖=0 𝑠𝑖 𝜕𝐸
𝜕𝑣𝑘𝑗=0,𝑖=1,𝑛𝑖,𝑗=1,𝑛𝑗,,𝑘=1,𝑛𝑘 (3.1.2.5 )
Și rezolvând sistemele de ecuații rezultate .
Din păcate în practică nu este adecvat să folosim această abordare deoarece relațiile
(3.1.2.5 ) se îndeplinesc în toate punctele de extrem ale hipersuprafeței definită de funcția de

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

28

performanță ceea ce face dificilă găsirea unui minim global în mod direct [6]. Din acest motiv,
algoritmele de antrenare ale rețelelor de tip MPL folosesc tehnici de opt imizare având ca
obiectiv minimizarea funcției de performanță .

3.2.2. Algoritme de antrenare ale rețelelor de tip MLP
3.2.2.1. Metoda gradientului descendent

Cea mai ușoară și simplă metodă de determinare a punctului de extrem local al unei
funcții constă în anularea derivatei [6]. Se va considera funcția 𝑦=𝑓(𝑥)
𝑓′(𝑥)=𝑑𝑓
𝑑𝑥=0 (3.2.2.1. 1)
În cazul în care funcția este neliniară și are o formă complexă, rezolvarea directă a
acesteia nu este posibilă și va fi nevoie să aplicăm o metodă iterativă. Pentru a determina
minimul funcției se va pornii de la o aproximație inițială 𝑥0 și prin folosirea relației de recurență
(3.2.2.1.3) se va genera un șir de aproximații succesive [6].
𝑥𝑡+1=𝑥𝑡−𝜂𝑑𝑓
𝑑𝑥|
𝑥=𝑥𝑡=𝑥𝑡−𝜂𝑓′(𝑥𝑡) (3.2.2.1.2 )
η – are o valoare pozitivă și este utilizat la amplificarea sau atenuarea deplasării în lungul
direcției 𝑓′(𝑥)=𝑑𝑓
𝑑𝑥 . În cazul în care acesta are valori mari, punctul de minim poate fi depășit,
iar dacă valorile acestuia sunt mici, găsirea punctului de minim o să fie prea lentă [6].

Fig 3.2.2.1.1 Principiul metodei gradientului descendent
Aplicând cele menționate anterior funcției de evaluare a performantei unei RNA, relația
(3.2.2.1.2) devine [6]:
𝑋𝑡+1=𝑋𝑡− 𝜂∇𝑓(𝑋𝑡) (3.2.2.1.3)
Cu ajutorul acesteia se determină relațiile pentru calcularea ponderilor conexiunilor
sinaptice ale neuronilor de pe stratul ascuns și de pe stratul de ieșire [6]:
𝑊𝑡+1=𝑊𝑡−𝜂∇𝐸(𝑊𝑡) (3.2.2.1.4.a)
𝑉𝑡+1=𝑉𝑡−𝜂∇𝐸(𝑉𝑡) (3.2.2.1.4.b)
𝜂 – viteza de învățare

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

29

3.2.2.2. Algoritmul backpropagation

Algoritmul BP sau algoritmul de propagare înapoi a erorii a fost introdus de Rumelhart
și membrii grupului "Parallel Distributed Processing" în 1986 și este primul algoritm propus
pentru antrenarea rețelelor de tip MLP, fapt care a contribuit la relansarea calculului neuronal
în inteligența artificială [6].
Deoarece acest algorit m folosește o metodă de gradient pentru minimizarea funcției de
performanță, atunci pe tot domeniul de definiție funcțiile de activare sau transfer ale neuronilor
trebuie să fie continue și derivabile, acest lucru putând fi realizat cu ajutorul unor funcți i
sigmoidale sau liniare [6].
Generarea rețelei MLP:
• etapa de învățare (parcurgere directă a rețelei) -pe baza mulțimii de antrenare se
determină ieșirile generate de rețea ( ponderile și valorile pragurilor de activare
ale neuronilor )[6].
• etapa de testare (parcurgere înapoi a rețelei) -ieșirile calculate de rețea se
compară cu ieșirile dorite pentru fiecare exemplu și pe baza acestora se
determină eroarea care se utilizează la actualizarea ponderilo r[6].
Pentru a modifica ponderile din cadrul etapei de învăț are vom folosii regula delta
generalizată :
Vom lua urm ătorul exemplu (𝑋𝑚,𝑑𝑚),𝑚=1,𝑛𝑒 din mulțimea de antrenare corec ția
unui neuron oarecare j conectat cu neuronul i din stratul precedent este proporțional ă cu un
termen de eroare 𝛿𝑗(𝑚) asociat neuronului j [6]:
∆𝑤𝑗𝑖(𝑚)=𝜂𝛿𝑗(𝑚)𝑦𝑖(𝑚) (3.2.2.2.1)
𝜂 – rata de învățare
𝑦𝑖(𝑚) – ieșirea neuronului i din stratul precedent
În cazul în care structura dată are un strat ascuns, atunci pentru a aplica regula delta se
va ține cont de poziția neuronului [6].
Corecția ponderilor de pe stratul de ieșire se face cu ajutorul următoarelor relații :
𝛿𝑘(𝑚)=(𝑑𝑘(𝑚)−𝑦𝑘(𝑚))𝑓′(𝑛𝑒𝑡𝑘(𝑚)) (3.2.2.2.2)
∆𝑣𝑘𝑗(𝑚)=𝜂𝛿𝑘(𝑚)𝑧𝑗(𝑚)=𝜂(𝑑𝑘(𝑚)−𝑦𝑘(𝑚))𝑓′(𝑛𝑒𝑡𝑘(𝑚))𝑧𝑗(𝑚) (3.2.2.2.3)
𝛿𝑘(𝑚) – termenul de eroare
∆𝑣𝑘𝑗(𝑚) – corec ția ponderii
(𝑑𝑘(𝑚)−𝑦𝑘(𝑚)) – diferen ța dintre valoarea dorit ă 𝑑𝑘(𝑚) și cea calculat ă de re țeaua
neuronal ă 𝑦𝑘(𝑚) pentru exemplu l m
𝑓′(𝑛𝑒𝑡𝑘(𝑚)) – derivata func ției de activare a neuronului k în raport cu intrarea net ă a
acestuia 𝑛𝑒𝑡𝑘(𝑚) pentru exemplul m
Corecția ponderilor de pe stratul ascuns se face cu ajutorul următoarelor relații :
𝛿𝑗(𝑚)=(∑𝑣𝑘𝑗(𝑚)𝛿𝑘(𝑚) 𝑛𝑘
𝑘=1 )𝑔′(𝑛𝑒𝑡𝑗(𝑚)) (3.2.2.2.4)
∆𝑤𝑗𝑖(𝑚)=𝜂𝛿𝑗(𝑚)𝑥𝑖(𝑚)=𝜂(∑𝑣𝑘𝑗(𝑚)𝛿𝑘(𝑚) 𝑛𝑘
𝑘=1)𝑔′(𝑛𝑒𝑡𝑗(𝑚))𝑥𝑖(𝑚) (3.2.2.2.5)
𝛿𝑗(𝑚) – termenul de eroare
∆𝑤𝑗𝑖(𝑚) – corec ția ponderii

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

30

Un dezavantaj major al algoritmului BP este convergența relativ lentă și dependența
acesteia de alți parametrii ( valorile inițiale ale pragurilor și ponderilor, viteza de învățare, tipul
funcțiilor de activare etc… )[6].
3.2.2.3. Algoritmul backpropagation cu moment

Ținând cont de metoda prezentată anteri or, rezultă că algoritmul BP de antrenare a
rețelelor MLP cuprinde două etape, procesul de calcul este iterativ iar corecțiile ponderilor unui
neuron oarecar e j conectat cu un neuron i aflat pe stratul precedent se calculează cu relația :
∆𝑤𝑗𝑖(𝑡)=𝜂𝛿𝑗(𝑡)𝑦𝑖(𝑡) (3.2.2.3.1)
Valoarea vitezei de învățare 𝜂 trebuie s ă fie :
• suficient de mică pentru a asigură convergența algoritmulu i
• suficient de mare pentru a obține un proces de învățare rapi d
Pentru a îmbunătății viteza de convergență a algoritmului BP vom introduce un termen
suplimentar în regula de corectare a ponderilor, termen denumit moment [6]:
∆𝑤𝑗𝑖(𝑡)=𝜂𝛿𝑗(𝑡)𝑦𝑖(𝑡)+∝∆𝑤𝑗𝑖(𝑡−1) (3.2.2.3.2)
Ținând cont de aceste condiții , corecția unei ponderi la pasul curent 𝑡𝑐 devine:
∆𝑤𝑗𝑖(𝑡)=𝜂∑∝𝑡𝑐−𝑡𝛿𝑗(𝑡)𝑦𝑖(𝑖)𝑡𝑐
𝑡=0 =𝜂∑∝𝑡𝑐−𝑡 𝑡𝑐
𝑡=0𝜕𝐸(𝑡)
𝜕𝑤𝑗𝑖(𝑡) (3.2.2.3.3)
• dacă în iterații consecutive derivatele parțiale ale funcției de performanță au același
semn atunci ponderea se va modifica cu o valoare din ce în ce mai mare deoarece
corecția crește în amplitudine, rezultând în accelerarea procesului de învățar e[6]
• dacă în iterații consecutive derivatele parțiale ale funcției de performanță au semne
diferit e atunci ponderea se va modifica cu o valoare din ce în ce mai mică deoarece
corecția scade în amplitudine, rezultând în încetinirea procesului de în vățare[6]
În concluzie, introducerea termenului moment în legea de corecție a ponderilor are ca
efect stabi lizarea procesului oscilatoriu de învățare .

3.2.2.4. Metoda Newton de ordinul II

Această metodă oferă o convergență rapidă a procesului iterativ de calcul către punctul
de minim căutat atunci când se cunosc detalii suplimentare privind funcția 𝑦=𝑓(𝑋), care să
permită estimarea unei valori optime pentru viteza de învățare η[6].
În cazul unei funcții de o singură variabilă 𝑦=𝑓(𝑋), dac ă se cunoaște o estimare
oarecare 𝑥𝑡 a punctului de extrem căutat se pune problema determinării corecției ∆𝑥𝑡 astfe l încât
𝑥=𝑥𝑡+1=𝑥𝑡+∆𝑥𝑡 să constituie soluția problemei [6], adică :
𝑓′(𝑥)=𝑓′(𝑥𝑡+∆𝑥𝑡)=0 (3.2.2.4. 1)
Dezvoltând f(x) in serie Taylor in jurul punctului 𝑥𝑡 și reținând termenii până la ordinul
II obținem :
𝑓(𝑥)=𝑓(𝑥𝑡)+𝑓′(𝑥𝑡)(𝑥−𝑥𝑡)+1
2𝑓′′(𝑥𝑡)(𝑥−𝑥𝑡)2 (3.2.2.4. 2)
Prin derivarea relației (3.2.2.4.2 ) obținem :
𝜕𝑓
𝜕𝑥=𝑓′(𝑥𝑡)+𝑓′′(𝑥𝑡)(𝑥−𝑥𝑡)=0 (3.2.2.4. 3)
𝑥=𝑥𝑡+1=𝑥𝑡−𝑓′(𝑥𝑡)
𝑓′′(𝑥𝑡) (3.2.2.4. 4)
Din metoda gradientului descendent rezultă că 1
𝑓′′(𝑥𝑡) reprezintă o estimare a vitezei de
învățare în punctul 𝑥𝑡.

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

31

În cazul unei funcții de mai multe variabile f:R'' →R, y=f(X) relația ( 3.2.2.4.4 ) capătă
formă vectorială :
𝑋𝑡+1=𝑋𝑡−[𝐻(𝑋𝑡)]−1∇𝑓(𝑋𝑡) (3.2.2.4.5 )
în care X = [ 𝑥1,𝑥2,……𝑥𝑛]𝑇
𝐻(𝑋𝑡) – matricea hessian
∇𝑓(𝑋𝑡) – gradientul func ției în punctul 𝑋 = 𝑋𝑡
Prin aplicarea metodei Newton la antrenarea rețelelor MLP rezultă următoarele relații
de recurență pentru modificarea ponderilor :
𝑊𝑡+1=𝑊𝑡−[𝐻(𝑊𝑡)]−1∇𝐸(𝑊𝑡) (3.2.2.4. 6.a)
𝑉𝑡+1=𝑉𝑡−[𝐻(𝑉𝑡)]−1∇𝐸(𝑉𝑡) (3.2.2.4. 6.b)
Pentru utilizarea relațiilor ( 3.2.2.4. 6.a) și ( 3.2.2.4.6 .b) va trebui să se realizeze
transformarea matricelor W și V în vectori de form a:
𝑊=[𝑤11,… ,𝑤1𝑛𝑖,… ,𝑤𝑛𝑗1,… ,𝑤𝑛𝑗𝑛𝑖]𝑇
𝑉=[𝑣11,… ,𝑣1𝑛𝑖,… ,𝑣𝑛𝑘1,… ,𝑣𝑛𝑘𝑛𝑖]𝑇
Adaptarea ponder ilor se realizează folosind viteze de învățare distincte pentru fiecare
pondere, ponderi ce se modifică de la o iterație la alta [6].
Metoda Newton are o aplicabilitate practică restrânsă datorită complexității extrem de
ridicate a calculelor ce trebuie e fectuate în cadrul fiecărei iterații, însă are ca avantaje un număr
extrem de redus de epoci pentru a converge [6].

3.2.2.5. Metoda Levenberg -Marquardt

Această metodă este o îmbunătățire adusă metodei Newton prin introducerea
următoarelor aproximații :
• Dacă matricea H este singulară, aceasta este înlocuită cu o matrice simetrică și
nesingulară [6]:
𝐻=𝐻+𝜇𝐼 (3.2.2.5.1)
• Deoarece funcția de performanță a unei RNA are forma unei sume de pătrate
matricea hessian poate fi a proximată prin [6]:
𝐻≅𝐽𝑇𝐽 (3.2.2.5.2)
Iar gradientul se determin ă cu relația [6]:
∆𝑓=𝐽𝑇𝑒 (3.2.2.5.3)

Unde 𝑒=[𝑒1(𝑋),…𝑒𝑘(𝑋),…𝑒𝑛(𝑋)]𝑇 și 𝐽 = matricea Jacobian
Ținând cont de aceste condiții, relația de bază din metoda Newton devine:
𝑋𝑡+1=𝑋𝑡−[𝐻(𝑋𝑡)+𝜇𝐼]−1∇𝑓(𝑋𝑡)=[𝐽𝑇(𝑋𝑡)𝐽(𝑋𝑡)+𝜇𝐼]−1𝐽𝑇(𝑋𝑡)𝑒(𝑋𝑡) (3.2.2.5.4)
Din relația prezentată anterior rezultă că :
• Dacă 𝜇 este zero, obținem metoda Newton în care hessianul este aproximat cu relația
𝐻≅𝐽𝑇𝐽 [6].
• Dacă 𝜇 are valori foarte ridicate, obținem metoda gradientului descendent cu un pas
foarte mic [6].
• Pentru a utiliza metoda la antrenarea RNA de tip MLP este necesară vectorizarea
simultană a matricelor ponderilor sinaptice [6].
𝑈=[𝑤11,…𝑤1𝑛𝑖,…𝑤𝑛𝑗1,…𝑤𝑛𝑗𝑛𝑖,𝑣11,…𝑣1𝑛𝑗,…𝑣𝑛𝑘1,…𝑣𝑛𝑘𝑛𝑗]𝑇

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

32

În concluzie metoda Levenberg -Marquardt este una din cele mai rapide metode de
antrenare a rețelelor de tip ML P de dimensiuni moderate și face parte din categoria metodelor
numite cvasi -Newton ( metode derivate din metoda Newton ) .

3.2.2.6. Metoda gradientului conjugat

Ținând cont de metoda gradientului descendent unde căutarea punctului de minim al
unei funcții cu mai multe variabile se face prin deplasări succesive în direcția antigradientului,
folosind următoarele relații :
𝑋𝑡+1=𝑋𝑡+𝜂𝑑𝑡 (3.2.2.6.1)
𝑑𝑡=−∇𝑓(𝑋𝑡) (3.2.2.6.2)
Deși funcția scade cel mai rapid în lungul antigradientului, nu este mereu indicată
deplasarea după această direcție pentru că nu conduce în mod sigur la cea mai rapidă
convergență, ceea ce ne determină să adoptăm un drum mai bun de căutare, format dintr -o
combinație liniară între direcția antigradientului și direcția anterioară de căutare [6].
𝑑𝑡=−∇𝑓(𝑋𝑡)+𝛽𝑡𝑑𝑡−1 (3.2.2.6.3 )
𝑑𝑡 și 𝑑𝑡−1 trebuie să fie conjugați pentru a nu fi modificată componența gradientului de –
a lungul direcției anterioare [6].
Se numesc canonic conjugați , doi vectori 𝑋,𝑌∈ 𝑅𝑛, dacă îndeplinesc relația :
𝑋𝑇𝑌=0 (3.2.2.6.4)
Generalizând relația anterioară , pentru o matrice H de ordinul n, dac ă vectorii 𝑋,𝑌∈
𝑅𝑛 satisfac rela ția 𝑋𝑇𝐻𝑌=0 atunci ace știa se numesc H conjuga ți.
Pentru a calcula minimul funcției 𝑓:𝑅𝑛→𝑅,𝑦=𝑓(𝑋) cu ajutorul metodei gradientului
conjugat, vom pornii din punctul inițial 𝑋0, alegem ca sens de deplasare sensul antigradientului
și determinăm un punct nou [6]:
𝑋1=𝑋0−𝜂0∇𝑓(𝑋0) (3.2.2.6.5)
Procesul se desfășoară conform relației de recurență î n n pași după n direcții H
conjugate:
𝑋𝑡+1=𝑋𝑡+𝜂𝑡𝑑𝑡 (3.2.2.6.6)
Pentru fiecare pas nou, direcția de căutare se stabilește conform relației ( 3.2.2.6.3 ) în
care 𝛽𝑡 se determină astfel încât vectorii să fie H conjugați [6]:
𝑑𝑡𝑇𝐻𝑑𝑡−1=0 (3.2.2.6.7)
Ținând cont de relația (3.2.2.6.3 ) relația anterioară poate fi scrisă sub forma:
[−∇𝐹(𝑋𝑡)+𝛽𝑡𝑑𝑡−1]𝑇𝐻𝑑𝑡−1=−[∇𝑓(𝑋𝑡)]𝑇𝐻𝑑𝑡−1+𝛽𝑡𝑑𝑡−1𝑇𝐻𝑑𝑡−1=0 (3.2.2.6.8)
De unde rezult ă factorul 𝛽𝑡
𝛽𝑡=[𝑑𝑡−1𝑇𝐻𝑑𝑡−1]−1[∇𝑓(𝑋𝑡)]𝑇𝐻𝑑𝑡−1 (3.2.2.6.9)
Pentru algoritmele care implementează această metodă este necesară resetarea periodică
a direcției de căutare la direcția antigradientului, resetarea se face atunci când număr ul de
direcții conjugate generate devine egal cu dimensiunea n a vectorului X [6].
Pentru a crește eficiența metodei gradientului și implicit a algoritmilor bazați pe aceasta
se va folosi ideea propusă de Powell și Beale în care resetarea la direcția antigradientului se
face dacă există o foarte mică ortogonalitate între gradientul curent și cel anterior [6].
Direcția de căutare este resetată la direcția antigradientului dacă inegalitatea se
îndeplinește [6]:
|[∇𝑓(𝑋𝑡−1)]𝑇∇𝑓(𝑋𝑡−1)|≥0.2‖∇𝑓(𝑋𝑡)‖2 (3.2.2.6.10)
Pentru a putea utiliza metoda gradientului conjugat la antrenarea RNA de tip MLP este
necesară vectorizarea simultană a matricelor ponderilor sinaptice W și V [6]:
𝑈=[𝑤11,…𝑤1𝑛𝑖,…𝑤𝑛𝑗1,…𝑤𝑛𝑗𝑛𝑖,𝑣11,…𝑣1𝑛𝑗,…𝑣𝑛𝑘1,…𝑣𝑛𝑘𝑛𝑗]𝑇

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

33

Pentru a implementa această metodă mai este necesară :
• Selectarea vitezei de învățare
• Calculul factorului 𝛽𝑡
Selectarea vitezei de învățare ( 𝜂𝑡 ) – o metodă de determinare a acesteia constă în
aproximarea variației funcției de performanță cu 𝜂 printr -un polinom de gradul II [6].
𝐸(𝜂)=𝐴𝜂2+𝐵𝜂+𝐶 (3.2.2.6.11)
Anulând derivată vom obține valoarea optimă
𝜂∗=−𝐵
2𝐴 (3.2.2.6.12 )
Pentru fiecare pas de calcul, coeficienții A,B și C se determină evaluând funcția de
performanță în trei puncte distincte [6]:
• Pentru 𝜂=0 rezult ă 𝑈1=𝑈𝑡+𝜂𝑑𝑡=𝑈𝑡 respectiv 𝐶=𝐸(𝑈1)=𝐸1
• Pentru 𝜂=0.5 rezult ă 𝑈2=𝑈𝑡+0.5𝑑𝑡 respectiv 1
4𝐴+1
2𝐵+𝐶=𝐸(𝑈2)=𝐸2
• Pentru 𝜂=1 rezult ă 𝑈3=𝑈𝑡+𝑑𝑡 respectiv 𝐴+𝐵+𝐶=𝐸(𝑈3)=𝐸3
Rezolvând sistemul de ecuații format din relațiile anterioare obținem coeficienții A,B și
C, iar din relația ( 3.2.2.6.12 ) valoarea vitezei de învățare .
Selectarea matricei H pentru calculul factorului 𝛽𝑡 – matrice a H din expresia factorului
𝛽𝑡 este hessianul funcției de performanță E[6]. Deoarece calculul acestei matrice este un proces
complex computațional s -au căutat metode pentru reducerea timpului de calcul. Aceste metode
permit calculul aproximativ al factorului 𝛽𝑡 folosind numai valorile gradientului fără a se
calcula matrice a Hessian.
Cele mai folosite metode sunt:
• Metoda Fletcher -Reeves în cadrul căreia factoru l 𝛽𝑡 se calculeaz ă ca raportul dintre
pătratul normei gradientului la pasul curent t și pătratul normei gradientului la pasul
precedent t-1
𝛽𝑡=‖∇𝐸(𝑈𝑡)‖2
‖∇𝐸(𝑈𝑡−1)‖2=[∇𝐸(𝑈𝑡)]𝑇[∇𝐸(𝑈𝑡)]
[∇𝐸(𝑈𝑡−1)]𝑇[∇𝐸(𝑈𝑡−1)] (3.2.2.6.13)
• Metoda Polak -Ribiere în cadrul căreia factorul 𝛽𝑡 se calculeaz ă ca raptorul dintre
produsul scalar al variațiilor componentelor gradientului la pasul precedent cu
gradientul la pasul curent t și pătratul normei gradientului la pasul precedent t-1
𝛽𝑡=∆[∇𝐸(𝑈𝑡−1)]𝑇∇𝐸(𝑈𝑡)
‖∇𝐸(𝑈𝑡−1)‖2 (3.2.2.6.14)

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

34

Capitolul 4
Studii de caz pe un sistem test

4.1. Studiul regimului permanent folosind metoda Newton -Raphson în
coordonate polare pentru rețeaua TEST2.

Acest studiu de caz are în vedere analiza regimului permanent în mai multe puncte de
funcționare, pentru rețeaua TEST2 din literatura de specialitate.
Schema rețelei este prezentată în figura 4.1.1:

1110
13121234
8
9765
CHE-RO2
2x170MWCTE-RO1
2x330MWCTE-ALF A
2x200MWCTE-GAMA
3x100MWCTE-BETA
2x200MW

Fig 4.1.1 Schema monofilară a rețelei TEST2

Datele nodale
Nod Tip Un
[kV] Pg
[MW] Pc
[MW] Qc
[MVAr] Uimp
[kV] Qmin
[MVAr] Qmax
[MVAr]
1 PQ 220 0 250 155 – – –
2 PU 220 255 0 0 230 0 190
3 PQ 220 0 60 35 – – –
4 PQ 220 0 190 130 – – –
5 PU 220 240 0 0 225 0 180
6 PQ 220 0 220 135 – – –
7 PU 220 240 0 0 225 0 180
8 PQ 220 0 65 35 – – –
9 PQ 220 0 130 70 – – –
10 PQ 220 0 200 140 – – –
11 PU 220 165 0 0 233 -60 160
12 EC 220 – 0 0 235 0 300
13 PQ 220 0 150 90 – – –

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

35

Datele liniilor electrice
Linia Nodul i Nodul k R
[] X
[] B0
[S] Imax
[A]
L1 1 2 9.90 60.6 416 875
L2 2 3 5.30 32.3 222 875
L3 2 4 7.40 45.2 311 875
L4 4 5 4.00 24.2 167 875
L5 5 6 3.30 20.9 544 1750
L6 6 7 6.65 37.9 243 805
L7 7 8 4.60 28.3 195 875
L8 7 9 4.00 24.2 167 875
L9 9 10 5.30 32.3 222 875
L10 1 10 5.50 31.5 202 805
L11 10 11 7.90 48.5 333 875
L12 11 13 5.30 32.3 222 875
L13 11 12 3.60 22.2 153 875
L14 12 13 5.50 33.5 230 875
L15 1 12 4.90 31.4 812 1750

Datele generatoarelor din centralele electrice
Centrala Sng
[MVA] Png
[MW] Ung
[kV] Reactanțe (raportate la Sng)
[%] Constante de timp
[s]
dX

qX
'
dX
'
qX
'
0dT
'
0qT H
RO1 388 330 24.00 215 215 36.4 26.9 6.13 0.724 6.1
RO2 190 175 15.75 114 84 38.0 28.0 6.60 0.210 8.1
ALFA 235 210 15.75 185 185 27.0 27.0 6.40 0.284 5.8
BETA 235 200 15.75 192 192 27.6 27.6 6.40 0.284 6.7
GAMA 125 100 10.50 206 206 20.5 20.5 6.00 0.25 7.9

Datele transformatoarelor bloc și ale celor din stațiile de alimentare a consumatorilor
Stația Sn
[MVA] Unf
[kV] Unr
[kV] usc
[%]
nom
scP
[kW] i0
[%]
0P
[kW] Reglaj
[%]
RO1 2×400 24 231 14,1 1150 0.45 280 ±12×1,33
RO2 2×190 15.75 242 11.7 500 0.6 150 ±1×5
ALFA 2×250 15.75 242 10.9 500 0.8 150 –
BETA 2×250 15.75 242 14.0 947 1.2 282 ±1×1,40
GAMA 3×125 10.5 231 13.5 375 1.2 107 ±13×1,25
C1 2×100 121 231 10.5 275 0.8 65 ±1×1,25
C2 2×200 121 231 10.5 485 0.8 105 ±12×1,25
C3 2×200 121 231 10.5 485 0.8 105 ±12×1,25
C4 2×100 121 231 10.5 275 0.8 65 ±12×1,25
C5 2×100 121 231 10.5 275 0.8 65 ±12×1,25
C6 2×200 121 231 15.5 708 0.7 147 ±11×1,45
C7 2×200 121 231 15.5 708 0.7 147 ±11×1,45
C8 2×200 121 231 15.5 708 0.7 147 ±11×1,45

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

36

Analizele au fost realizate în MATLAB cu ajutorul programului E PSA.
Meniul programului este prezentat în figura 4.1.2:

Fig.4.1.2 Meniu program EPSA.

Pentru selectarea rețelei pe baza cărora se vor efectua calculele se utilizează comanda 1
„Selecție rețea”.
Pasul următor este de a calcula matricea admitanțelor nodale Ynn prin utilizarea
comenzii 2 „Calcul Ynn”.
Pentru calculul regimului permanent se va utiliza comanda 3 „Calcul regim permanent”,
care va deschide fereastr a de opțiuni din figura 4.1.3 unde se poate alege metoda de calcul a
regimului permanent.

Fig 4.1.3 Meniu calcul regim permanent.

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

37

Într-o primă etapă am analizat regimul mediu de bază (RMB) prin metoda Newton –
Raphson în coordonate polare, rezultatele fiind cuprinse în tabelul de mai jos:

Rezultatele analizei RMB

U Pc Qc Pg Qg
0.966 2.500 1.550 0.000 0.000
1.045 0.000 0.000 2.550 1.245
1.019 0.600 0.350 0.000 0.000
0.979 1.900 1.300 0.000 0.000
1.023 0.000 0.000 2.400 1.460
0.978 2.200 1.350 0.000 0.000
1.023 0.000 0.000 2.400 1.718
0.998 0.650 0.350 0.000 0.000
0.964 1.300 0.700 0.000 0.000
0.944 2.000 1.400 0.000 0.000
1.059 0.000 0.000 1.650 1.160
1.028 1.500 0.900 0.000 0.000
1.068 0.000 0.000 3.910 1.729RMB
11
12
135
6
7
8
9
104Nr Nod
1
2
3

În urma efectuării calculelor regimului permanent pentru rețeaua TEST2 care
funcționează în regimul mediu de bază, am constatat că nodul cu gradul de risc de instabilitate
cel mai ridicat, în raport cu perturbațiile din rețea este nodul 10 .
Constatăm că valoarea tensiunii în nodul 10 (0.944 u.r.) este apropiată de valo area
minimă admisibilă (0.9 u.r.) a tensiunii pentru rețelele de transport .

Fig 4.1.4 Graficul valorilor tensiunilor nodale.
0.880.90.920.940.960.9811.021.041.061.08
1 2 3 4 5 6 7 8 9 10 11 12 13U [ u.r ]
NodRMB

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

38

În etapa a ÎI -a, am analizat regimul de sarcină minimă (RSMin) pentru care puterile
active și reactive consumate au fost reduse la 60% din valorile lor nominale, în scopul
evidențierii fa ptului că tensiunile nodale cresc odată cu scăderea puterii consumate .

Rezultatele analizei RSMin

U Pc Qc Pg Qg
1.005 1.500 0.930 0.000 0.000
1.045 0.000 0.000 2.550 0.607
1.032 0.360 0.210 0.000 0.000
1.002 1.140 0.780 0.000 0.000
1.023 0.000 0.000 2.400 0.481
0.999 1.320 0.810 0.000 0.000
1.023 0.000 0.000 2.400 0.759
1.009 0.390 0.210 0.000 0.000
0.991 0.780 0.420 0.000 0.000
0.99 1.200 0.840 0.000 0.000
1.059 0.000 0.000 1.650 0.323
1.045 0.900 0.540 0.000 0.000
1.068 0.000 0.000 -1.213 1.513RSMin
11
12
135
6
7
8
9
104Nr Nod
1
2
3

În urma efectuării analizei în regimul de sarcină minimă (RSMin), se poate observa
efectul scăderii puterilor active și reactive consumate asupra tensiunii nodale .
Constatăm că, nodul cu cel mai mare risc de instabilitate a devenit nodul 13, a cărui
valoare (1.068 u.r.) este apropiată de pragul maxim admisibil (1.1 u.r.) al tensiunii în rețeaua de
transport .

Fig 4.1.5 Graficul valorilor tensiunilor nodale în cazul RSMin 0.850.90.9511.051.1
1 2 3 4 5 6 7 8 9 10 11 12 13U [ u.r.]
NodRSMin (0.6 din RMB )
Tensiunea nodala in cazul RSMin Tensiunea nodala in cazul RMB

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

39

În etapa a III -a, am analizat regimul de sarcină maximă (RSMax) pentru care puterile
active și reactive consumate au fost crescute la 115% din valorile lor nominale, în scopul
evidențierii faptului că tensiunile nodale scad odată cu creșterea puterii consumate .

Rezultatele analizei RSMax

U Pc Qc Pg Qg
0.939 2.875 1.782 0.000 0.000
1.045 0.000 0.000 2.550 1.675
1.014 0.690 0.403 0.000 0.000
0.969 2.185 1.495 0.000 0.000
1.023 0.000 0.000 2.400 1.876
0.969 2.530 1.552 0.000 0.000
1.022 0.000 0.000 2.400 2.250
0.993 0.748 0.403 0.000 0.000
0.945 1.495 0.805 0.000 0.000
0.913 2.300 1.610 0.000 0.000
1.058 0.000 0.000 1.650 1.700
1.022 1.725 1.035 0.000 0.000
1.068 0.000 0.000 5.994 2.198RSMax
11
12
135
6
7
8
9
104Nr Nod
1
2
3

În urma efectuării analizei în regimul de sarcină maximă (RSMax), se poate observa
efectul creșterii puterilor active și reac tive consumate asupra tensiunii nodale.
Constatăm că, nodul cu cel mai mare risc de instabilitate a devenit nodul 10, a cărui
valoare (0.913 u.r.) este foarte apropiată de pragul minim admisibil (0.9 u.r.) al tensiunii în
rețeaua de transport.
Pentru acest regim de funcționare au fost realizate mai multe analize din care s -a
constatat faptul că pentru valori mai mari de 1.15 din puterile nominale cerute de consumatori,
este foarte probabilă apariția efectului de instabilitate de tensiune .

Fig 4.1.6 Graficul valorilor tensiunilor nodale în cazul RSMax 0.80.850.90.9511.051.1
1 2 3 4 5 6 7 8 9 10 11 12 13RSMax (1.15 din RMB)
Tensiunea nodala in cazul RSMax Tensiunea nodala in cazul RMB

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

40

În etapa a IV -a, am analizat regimul mediu de bază (RMB) în ipoteza în care linia ce
face legătura între nodurile 10 și 11 a fost deconectată, pentru a identifica efectele unei astfel
de situații .

Rezultatele analizei unei contingente

U Pc Qc Pg Qg
0.900 2.500 1.550 0.000 0.000
1.045 0.000 0.000 2.550 1.816
1.019 0.600 0.350 0.000 0.000
0.978 1.900 1.300 0.000 0.000
1.023 0.000 0.000 2.400 1.780
0.965 2.200 1.350 0.000 0.000
0.990 0.000 0.000 2.400 2.250
0.965 0.650 0.350 0.000 0.000
0.890 1.300 0.700 0.000 0.000
0.821 2.000 1.400 0.000 0.000
1.059 0.000 0.000 1.650 -0.058
1.028 1.500 0.900 0.000 0.000
1.068 0.000 0.000 4.130 3.149RMB
11
12
135
6
7
8
9
104Nr Nod
1
2
3

În urma analizei rețelei TEST2 cu o linie indisponibilă (linia care face legătura între
nodurile 10 și 11) am constatat faptul că nodul 10 are cel mai ridicat grad de vulnerabilitate la
instabilitate de tensiune.
De asemenea, valoarea tensiunii în nod ul 10 (0.821 u.r.) este mult mai scăzută decât
valoarea minimă acceptată (0.9 u.r) în rețelele de transport .

Fig 4.1.7 Graficul valorilor tensiunilor nodale în cazul RMB cu contingență 0.0000.2000.4000.6000.8001.0001.200
1 2 3 4 5 6 7 8 9 10 11 12 13U [ u.r ]
NodRMB cu contingenta
Tensiunea nodala in cazul RMB cu contingenta Tensiunea nodala in cazul RMB fara contingenta

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

41

4.2. Studiul stabilității de tensiune pentru rețeaua TEST2. Calculul
indicatorului global VSI.

Acest studiu de caz are în vedere calculul indicatorului global VSI, în vederea
determinării dacă nodul analizat este stabil sau nu .
Pentru această analiză am ales nodul 9 din figura 4.2.1:

1110
13121234
8
9765
CHE-RO2
2x170MWCTE-RO1
2x330MWCTE-ALF A
2x200MWCTE-GAMA
3x100MWCTE-BETA
2x200MW

Fig 4.2.1. Schema monofilară a rețelei TEST2

Aceast ă schem ă poate fi redus ă la:

Fig 4.2.2. Schema monofilar ă echivalent ă de calcul pentru studiul stabilității nodului 9 din
rețeaua TEST2

Pentru determinarea indicatorului global VSI, am realizat în MATLAB un program de
calcul care are ca inputuri (intrări) datele despre nod, datele despre transformatorul care face
legătura cu consumatorul, datele despre linie și tensiunea de bază .

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

42

Într-o primă etapă am întocmit schema echi valentă prin reprezentarea liniei electrice
printr -un cuadripol Π, iar transformatoarele printr -un cuadripol Γ.

Fig 4.2.3. Schema echivalent ă pentru studiul stabilității de tensiune.

Pentru această schemă echivalentă programul calculează parametrii longitudinali și
transversali, ținând cont că plotul de funcționare al transformatorului este -4.

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

43

În etapa a II -a, am realizat schema pentru calculul echivalentului Thevenin.

Fig 4.2.4 Schema de calcul redusa pentru echivalentul Thevenin.

Pentru calculul echivalentului Thevenin s -au ținut cont de notațiile făcute anterior,
acestea fiind:
𝑌1=𝑦13
𝑌2=𝑁𝑦23
𝑌10=𝑦130

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

44

𝑌20=𝑁(𝑁−1)𝑦23
𝑌30=𝑦310+𝑦320+(1−𝑁)𝑦23
𝐸𝑇ℎ=𝑌23𝑌31
𝑌22𝑌33−𝑌23𝑌32𝑈1
𝑌𝑇ℎ=𝑌20+𝑌1(𝑌1+𝑌30)
𝑌1+𝑌2+𝑌30
𝑍𝑇ℎ=1
𝑌𝑇ℎ

După calculul echivalentului lui Thevenin, folosind formulele din literatura de
specialitate referite în subcapitolul 2.4 din Capitolul 2, vom calcula indicatorul global VSI

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

45

După rularea programului de calcul in MATLAB am obținut următoarele rezultate
prezentate în tabelul de mai jos:

Rezultate calcul VSI

0 0 1.0838 59.59 0 – 1
0.91 0.49 1.0724 58.472 0.0175 -120.335 0.9813
1.3 0.7 1.0673 57.9754 0.0252 -118.0608 0.9729
1.56 0.84 1.0639 57.6381 0.0303 -116.5505 0.9673
1.82 0.98 1.0605 57.2956 0.0355 -115.0445 0.9615Punct de functionare A Punct de functionare BVSI 𝑃2𝑄2𝑈2𝐴 𝑈2𝐵 𝐽𝑅 𝐽𝑅

Din analiza acestora, se poate observa efectul creșterii puterii active și reactive asupra
indicatorului global VSI.
Cu cât acestea cresc mai mult, valoarea indicatorului tinde să scadă spre 0 indicând
faptul că sistemul se destabilizează.
Indicatorul atinge valoarea 0 în cazul în care puterile active și reactive ajung la valorile
lor maxime .

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

46

4.3. Utilizarea RNA de tip MLP pentru analiza stabilității nodului 9 din
rețeaua TEST2.

Pentru antrenarea rețelei neuronale am realizat un script în MATLAB care să genereze
aleator 3000 de valori (date de intrare) necesare calculării indicatorului global VSI, care va fi
target -ul rețelei .

Am ales ca structura rețelei să fie constituită din 3 intrări, 15 neuroni pe stratul ascuns
și o ieșire .

Fig 4.3.1 Structura MLP

Prin încercări repetate am ajuns la concluzia că din cele 3000 de date de intrare cea mai
bună distribuție a acestora este 2400 în antrenare 300 în validare și 300 în testare .

Fig 4.3.2. Distribuția datelor de intrare

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

47

În urma rulării rețelei neuronale artificiale, aceasta a generat cel mai bun răspuns după
114 epoci, eroare fiind de 10^(-5).

Fig 4.3.3 Eroare răspuns MLP

Fig 4.3.4. Histograma erorilor la antrenarea MLP

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

48

In figura 4.3.5 se poate observa un răspuns foarte bun al rețelei neuronale la antrenare,
însemnând ca aceasta a fost antrenat ă cu succes.

Fig 4.3.5 Răspunsul rețelei la antrenare

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

49

Concluzii

➢ Pentru monitorizarea riscului producerii instabilității de tensiune în SEE folosind RNA de tip
MLP am dezvoltat un program de calcul în MATLAB.
➢ Ca etapă intermediară, am rulat programul EPSA pentru a analiza diferite regimuri de
funcționare a rețelei TEST2 și a pune în evidență efectul modificării puterilor consumate asupra
tensiunii nodale și a încărcării liniilor .
➢ Am ales pentru calculul regimului permanent metoda Newton -Raphs on în coordonate polare
în detrimentul metodei Newton -Raphson în coordonate rectangulare, deoarece prima oferă
un avantaj în calul indicatorului global VSI folosit în cadrul analizei instabilității de tensiune .
➢ În urma rulării programului dezvoltat am pus în evidență efectul modificării puterilor active și
reactive asupra stabilității tensiunii în noduri .
➢ În urma analizei datelor obținute pentru nodul 9 din rețeaua TEST2 care funcționează în
regimul mediu de bază (RMB), am constatat că valoarea indicatorul ui VSI (0.9729) este foarte
apropiată de 1 fapt ce, relevă că nodul este stabil .
➢ Rularea rețelei neuronale artificiale (RNA) de tip MLP a demonstrat faptul că există
posibilitatea implementării unei rețele de acest tip pentru monitorizarea riscului produce rii
instabilității de tensiune în SEE .
➢ În final, având în vedere complexitatea procesului dezvoltat în cadrul rețelei neuronale
artificiale de tip MLP și beneficiile aduse prin identificarea timpurie a riscului producerii unei
instabilități de tensiune în SEN, mă voi preocupa pe viitor să îmi dezvolt cunoștințele în
domeniu, pentru aplicarea rețelelor artificiale la monitorizarea SEN .

Universitatea POLITEHNICA din București
FACULTATEA DE ENERGETICĂ
Departamentul Sisteme Electroenergetice
060042 București, Splaiul I ndependenței, nr. 313, sector 6 <Sigla
departamentului
care coordonează
programul de
licență>

50

Bibliografie

[1] C.Bulac – Curs Teoria și modelarea SEE.
[2] P.Kundur – Power System Stability and Control, McGraw -Hill, Inc. New York, 1994.
[3] I. Tri știu – Curs Rețele Electrice
[4] M. Eremia, H. Cri șciu, B. Ungureanu, C. Bulac – Analiza asistată de calculator a regimurilor
sistemelor electroenegetice. Editura Tehnic ă, București 1985.
[5] C. Bulac, M. Eremia – Dinamica sistemelor electroenergetice. Editura Printech, 2006.
[6] M. Eremia, Gh. Car țină, D. Petricic ă, A-I Bulac, C. Bulac, I. Tri știu, Gh. Grigora ș – Tehnici de
inteligentă artificială în conducerea sistemelor electroenergetice, Editura AGIR, Bucure ști 2006.
[7] J. Machowski, J.W. Bialek, J.R. Bumby – Power System Dynamics Stability and Control, John Wiley
& Sons, United Kingdom 2008 .
[8] V.M. da Costa, N. Martins, J.L.R. Pereira, Transactions on Power Systems, 1999.