Abstract

Risolvere con Scilab una qualsiasi rete lineare in corrente continua e alternata, in regime stazionario, a partire dalla descrizione dei rami.

Introduzione

Sfruttando le potenzialita' di programmazione di Scilab, specie nel calcolo matriciale, cerchiamo di eliminare anche il lavoro di scrittura del sistema risolutivo di una rete elettrica, in corrente alternata, ma anche in corrente continua, con la semplice descrizione dei rami.

La scelta del metodo risolutivo ricade su quello "dei potenziali ai nodi" molto piu' facile da implementare rispetto a quello "delle correnti di maglia"; si tratta di trasformare ogni ramo nel circuito equivalente secondo Norton-Mayer; parallelo fra generatore di corrente Jeq e impedenza equivalente Zeq.

La necessita' di ammettere anche lati a impedenza nulla (sola f.e.m. o corto circuito) impedisce il calcolo della corrente di corto circuito di ramo (pari ella Jeq) e impone la considerazione di nuove correnti incognite oltre ai potenziali.

Cio' fatto, il metodo consiste nella semplice applicazione del principio di Kirchhoff ai nodi.

Input

L'input della matrice complessa m() che descrive gli rm rami della rete e' realizzato editando lo script di Scilab.

I rami della rete possono contenere: generatori di tensione, impedenze, generatori di corrente e loro combinazioni. In regime continuo gli elementi saranno ovviamente rappresentati da un numero reale mentre in regime alternato sinusoidale saranno rappresentati da numeri complessi a+j*b, anche in forma polare pc(modulo,fase) con fase in gradi.

Fissati arbitrariamente i versi delle correnti di ramo da I₁ a I_rm, assegnato un numero ad ogni nodo da 1 a nn e copiato lo script sotto riportato nell'editor di Scilab, si inseriscono a partire dal ramo con corrente I₁, in ordine, fino al ramo I_rm :

1) nodo iniziale,

2) nodo finale (in relazione al verso della corrente, che andra' dal primo verso il secondo),

3) forza elettromotrice (positiva se diretta verso il nodo finale),

4) impedenza (o resistenza nel caso di rete in corrente continua),

5) generatore di corrente (positivo se concorde con il verso della corrente).

N.B. Ricordo che:

a) agli elementi mancanti nel ramo si deve assegnare il valore 0,

b) il nodo numero 1 sara' il nodo di riferimento a potenziale zero V(1)=0,

c) si dovranno separare i valori inseriti con 1 o piu' spazi o con una virgola, un "punto e virgola" a fine ramo [con eccezione per l'ultimo!]

Il programma elabora ogni ramo in successione, da 1 a rm, costruendo progressivamente: la matrice dei coefficienti A() e il vettore dei termini noti b(), che per convenienza di scrittura codice hanno dimensioni

n.equazioni = n.nodi + n.rami.anomali.

Solo successivamente si elimina la prima riga e la prima colonna per risolvere il sistema

con n.equazioni = (n.nodi − 1) + n.rami.anomali

attaverso una semplice "divisione a sinistra".

.

Esempio 1 - regime periodico alternato sinusoidale (PAS)

Per un esempio applicativo usiamo la stessa rete proposta nel tutorial del Prof. Richardson, gia' considerato nell'articolo Scilab e numeri complessi di Electroportal. Riporto comunque l'indirizzo del file:

Working with Complex Numbers and Matrices in Scilab.

fig.1 "rete usata nello script"

Sono stati inseriti per primi, solo per convenienza di confronto con i risultati di Richardson, i lati esterni, dove le correnti di maglia uguagliano le correnti di ramo.

Lo Script

(versione alfa)

   //            microSpice  
  //----------------------------------------------------------- 
 //     Copyright   Renzo DF 2008 
//------------------------------------------------------------- 
   format('v',6); j=%i; k0=%pi/180;
//-------- cartesiana --> polare ----------
 function[mf]=cp(c)
   mf=[abs(c),phasemag(c)];  
 endfunction
//------ polare --> cartesiana ------------
 function [c] = pc(m,f)
  if argn(2) == 1 then
    c = m(:,1:2:$).*cos(m(:,2:2:$)*k0) + j*m(:,1:2:$).*sin(m(:,2:2:$)*k0);
  elseif argn(2) == 2 then  c = m.*cos(f*k0) + j*m.*sin(f*k0); 
  end
endfunction
     //--------------------------------------------------------
    //     INSERIRE gli elementi dei rami della rete  
   //    secondo l'ordine indicato nel commento 
m=[...  
   //-----nodo-------generatore-----impedenza-----generatore---
  // inizio| fine |   tensione    (o resistenza)   corrente
 //   q    |   w  |      E       |      Z       |     J       |
//-------------------------------------------------------------
      1        2    pc(100,-90)         0             0       ;   
      4        1       -500             0             0       ;
      2        4        0             -j*20           0       ; 
      3        4        0               20            0       ;
      3        2        0              j*12           0       ;
      1        3        0               80            0
 ];
  //-----------------------------------------------------------
 // estraggo sottomatrice (1^ e 2^ colonna di m)
  y= real(m(:,1:2)); 
// determino "numero di nodi" 
  nn=max(y);
// determino "numero di rami" (righe matrice m)  
  nr=size(m,1);
//-----------------------------------------------   
s=0;
for i=1 : nr 
// se J==0 e Z==0 incremento n.rami anomali
    if   m(i,5) == 0 then s=  s + 1 ; 
     end
end; neq=nn+s ;  // faccio spazio in A() e b()   
 //-----------------------------------------------
// azzero coefficienti,termini noti,correnti di ramo 
    A=zeros(neq,neq); b=zeros(neq,1); I=zeros(nr,1);
 //
// costruisco matrice A() e vettore b()
r = 1;
for k=1 : nr
 q=y(k,1); w=y(k,2); E=m(k,3); Z=m(k,4); J=m(k,5); 
   if J<>0 then  
        b(q)=b(q)+J;
        b(w)=b(w)-J;
    
  elseif Z<>0 then // se J==0 e Z<>0
        A(q,q)=A(q,q)-1/Z;
        A(w,w)=A(w,w)-1/Z;
        A(q,w)=A(q,w)+1/Z;
        A(w,q)=A(w,q)+1/Z;
        b(q)=b(q)+ E/Z;
        b(w)=b(w)- E/Z;
          
      else // se J==0 e Z==0
           A(nn+r,q)=-1;      //-Vi
           A(nn+r,w)=+1;     //+Vf 
           b(nn+r)  = E;    // = E
           A(q,nn+r)=-1;   //  -Ik
           A(w,nn+r)=+1;  //   +Ik
         r=r+1;
      end 
  end 
  // per comodita' pongo V(1) = 0
 // nodo n°1=nodo di riferimento
// ed elimino 1^riga e 1^colonna
   A2=A([2:$],[2:$]); b2=b([2:$],1); 
   V2 = A2\b2; // calcolo V(2)...V(nn)
   V   =[0;V2] // aggiungo V(1)=0        
  //--------------------------------------> out V()
 // Calcolo le correnti nei rami con verso
// da "nodo inizio" a -> "nodo fine"
r=1;  
for k=1:nr
 q=y(k,1); w=y(k,2); E=m(k,3); Z=m(k,4); J=m(k,5);
 if J<>0 then I(k)= J;  
 elseif Z<>0 then I(k)=((V(q)-V(w))+ E)/Z;
    // se Z<>0 uso legge Ohm generalizzata
     else I(k)=V(nn+r);//se Z=0 pesco I dal fondo di V()
     r=r+1; 
  end
end
//-------OUT correnti di ramo--------------
                Ic = I  
//-------------------------------------------- eof

(chiaramente le righe di codice relative alle funzioni cp() e pc() possono essere eliminate se e' gia' stato lanciato uno script che implementi dette trasformazioni!).

Output

L'output di Scilab e' il seguente. I primi 4 elementi di V() rappresentano ordinatamente i potenziali dei nodi 1,2,3,4 mentre gli ultimi 2 elementi di V() rappresentano le 2 correnti incognite addizionali introdotte nel sistema per la presenza dei rami 1-2 e 4-1 "anomali".

Ic() contiene, le correnti in forma cartesiana. (in ordine dal primo all'ultimo ramo)

V  =
   0
   0.000 - 100.i
   96. + 128.i
   500.
 - 14. - 17.i
 - 15.2 - 18.6i
Ic  =
 - 14. - 17.i
 - 15.2 - 18.6i
   5. - 25.i
 - 20.2 + 6.4i
   19. - 8.i
 - 1.2 - 1.6i

Come ultime considerazioni di calcolo ricordo che la differenza di potenziale puo' essere ricavata dalla differenza fra i potenziali del nodo iniziale e finale

-->V42=V(4)-V(2)
V42  =
   500. + 100.i

e se si desidera convertire le correnti Ic in forma polare Ip

-->Ip = cp(Ic)
Ip  =
   22.02  - 129.5
   24.02  - 129.3
   25.5   - 78.69
   21.19    162.4
   20.62  - 22.83
   2.     - 126.9

E quindi, per esempio

o alternativamente nel dominio del tempo

con (se f=50Hz)

(non ω = 314,16 rad/s).

Esempio 2 - regime continuo

Vista la recente introduzione del generatore di corrente, approfitto per descrivere passo-passo l'evoluzione dell'elaborazione, dall'input della rete all'output finale.

Analizzeremo la semplice rete in c.c. di fig.2, che contiene pero' tutte le diverse tipologie di ramo. Cominciamo fissando i versi delle correnti (arbitrari) e numeriamo i nodi (arbitrariamente ma ricordando che il nodo 1 viene assunto come nodo di riferimento a zero volt). Il risultato finale e' il seguente

fig.2

Lanciamo Scilab -> apriamo l'Editor -> copiamo lo Script (CTRL-C) -> incolliamolo (CTRL-V) nella finestra dell'Editor -> inseriamo i valori, riga per riga, cancellando quelli gia' presenti, a partire dal ramo relativo a I1 percorrendolo dal nodo 1 al nodo 2 -> proseguendo con il ramo della I2 e cosi' via.

Il risultato sara' il seguente (evidenziati: il ramo e la riga associati al generatore di corrente ... dal nodo 3 al nodo 2)

fig.3

che sostituendo i valori numerici diventa

.

In questa matrice m() c'e' la descrizione completa del circuito.

A questo punto vogliamo passare il compito della soluzione della rete a Scilab attraverso la scrittura di uno script (serie di istruzioni) che, partendo da questa serie di valori numerici (25 nell'esempio), ci restituisca un vettore Ic()contenente le correnti di ramo.

Per "istruire" Scilab e' ovviamente indispensabile conoscere le procedure e le leggi usate in una soluzione "manuale" del problema.

Soluzione manuale

In una soluzione manuale, si dovrebbe in primo luogo scegliere il "migliore " metodo risolutivo; riferendosi alla scelta che porta piu' facilmente al risultato finale.

Nella rete di fig.5 le alternative potrebbero essere:

a) Kirchhoff, (sistema con 2 equazioni ai nodi + 2 equazioni alle maglie),

b) correnti di maglia, (sistema a 3 equazioni),

c) potenziali ai nodi, (sistema a 3 equazioni),

d) mista Thevenin-Millman-Ohm-Kirchhoff.

La soluzione della rete di fig.5 risulterebbe davvero semplice usando la metodologia "d"; basterebbe infatti:

1) trasformare con Thevenin il blocco (E1,R2 e J3) in un semplice (Jeq=3),
2) applicare Millman e ricavare V31=(10-5)/(1/2+1/2)=5 V,
3) calcolare con Ohm I5=(V31+E5)/R5=15/2=7,5 A;  I4=V31/R4=5/2=2,5 A;
   I2=E1/R2= 5/1=5 A; I3=-J3=-10 A  e con Kirchhoff  I1=I2-I3=5+10=15 A

L'implementazione del metodo d, gia' difficile sugli "umani", risulta quasi impossibile sui "computer" e la lasciamo a chi si occupa di intelligenza artificiale.

Intendendo trovare una soluzione generale al problema la scelta si orienta sulle due alternative b e c. Nella pratica comune direi che si ricade spesso sulle "correnti di maglia" ma, considerando le due possibile alternative nell'ottica della programmazione, a mio parere, non si puo' che decidere per il metodo dei "potenziali ai nodi".

Scelta questa strada, dobbiamo ridisegnare la rete trasformando ogni ramo secondo Norton e successivamente applicare il primo principio di Kirchhoff ad ogni nodo meno uno.

Nel processo di trasformazione troviamo che:

a) il ramo relativo a I1 NON e' trasformabile -> dobbiamo mantenere I1 come incognita!

b) il ramo relativo a I2 non necessita di trasformazione,

c) il ramo relativo a I3 si trasforma semplicemente (cancellando R3),

d) il ramo relativo a I4 non necessita di trasformazione,

e) il ramo relativo a I5 e' trasformabile calcolando J5=E5/R5

fig.4

Dopo la "cura" la rete si trasforma in quella di fig.5 dove avremo come incognite i potenziali V1,V2,V3 e la corrente I1:

fig.5

A questo punto si possono scrivere:

a) le 3 equazioni ai nodi usando il primo principio di Kirchhoff, come sappiamo

    solo 2 saranno indipendenti (n.nodi-1)

b) una relazione addizionale che deriva dalla conoscenza della differenza fra V2 e V1.

Le scriviamo comunque tutte per ragioni di chiarezza

e sostituendo i valori numerici

Delle prime 3 normalmente si elimina la piu' complessa; anzi, non si scrive nemmeno!

Ci si ritrova con 3 equazioni e 4 incognite; i potenziali non saranno calcolabili se non assegnando un valore di riferimento arbitrario ad uno di essi (nel nostro caso poniamo per es. V1=0).

Note in modo univoco saranno invece le differenze di potenziale.

Non ci resta che risolvere il seguente sistema lineare

A * V = b  :

(dove notiamo che l'ultima riga del vettore V() e' in realta' una corrente!)

In passato la soluzione era tutta manuale,

oggigiorno manuale solo la scrittura di A() e b(), via computer V=A\b.

Soluzione automatica

La scopo dello script proposto in questo articolo e' proprio quello di assegnare a Scilab anche il compito di ricavare A() e b() (che scansafatiche !!!).

Riprendendo il sistema iniziale, a 4 equazioni (non indipendenti ) e 4 incognite, scritta la matrice dei coefficienti e dei termini noti in forma letterale, incominciamo a considerare i rami con resistenza diversa da zero. Notiamo come il ramo con [q=1, w=3, R=R5] compaia negli elementi A(1,1) e A(3,3) con segno negativo mentre negli elementi A(1,3) e A(3,1) con segno positivo,

fig.6

e cosi' anche per gli altri due rami resistivi.

Di conseguenza, le prime istruzioni da fornire a Scilab per formare A(), detta x la generica resistenza, q il nodo di partenza, w il nodo di arrivo e k il generico ramo, saranno:

 //se nel ramo J==0 e Z<>0 
        A(q,q)=A(q,q)-1/Z;
        A(w,w)=A(w,w)-1/Z;
        A(q,w)=A(q,w)+1/Z;
        A(w,q)=A(w,q)+1/Z;
      

Le correnti di corto circuito di ramo (nell'esempio solo J5), con la convenzione presa (positive le correnti entranti nei nodi), daranno contributo di segno opposto nel vettore b() dei termini noti (i segni sono invertiti a causa dello spostamento a secondo membro dell'equazione) e saranno implementabili con:

       b(q)=b(q)+ E/Z;
       b(w)=b(w)- E/Z;

similmente si terra' conto dei lati che contengono un gen. di corrente quando J<>0 (nel nostro caso J3)

Infine l'r-esimo ramo anomalo, se nn=numero nodi, imporra':

  //se nel ramo J==0 e x==0

a) una riga (equazione) addizionale in A() ,

        A(nn+r,q)=-1;    //- Vi
        A(nn+r,w)=+1;   // +Vf 
        b(nn+r)= E;    // = E

b) una colonna (incognita) addizionale in b(),

        A(q,nn+r)=-1;   //  -Ik
        A(w,nn+r)=+1;  //  +Ik

Solo dopo avere completato la scrittura di A() e b() viene eliminata la prima riga e la prima colonna con:

A2=A([2:$],[2:$]); b2=b([2:$],1); 
                       //2:$  significa  "dal secondo all'ultimo elemento"

Abbiamo tutto il necessario per far calcolare i potenziali V(2)...V(nn) con un semplice:

V2 = A2\b2;

---

Flow chart

Lo script puo' essere sintetizzato dal seguente diagramma di flusso :

autom 4.jpg

... to be continued