Note tecniche.

    In questa pagina ho cercato, oltre che a discolparmi, di condensare una descrizione dei dettagli implementativi, anche commentando alcune parti di codice. Consiglio di dare prima un' occhiata alla struttura del programma. (torna all'indice).


     
     
     
     
     

    Contenuti

    • Considerazioni generali
    • Sulla unit Gauss
    • Sulla Unit Newton
    • Sulla Unit Stott

     
     
     
     
     
     
     
     
     




    Considerazioni generali
    (vedi la struttura)       Torna  Su
    allocazione statica
    Eccetto un caso strettamente necessario, non ho usato l'allocazione dinamica della memoria per dimensionare tutti gli arrays che popolano il programma (nonostante il C++ sia fornito di comodissime istruzioni in proposito).
    Oltre al semplificarmi la vita risparmiandomi codice poco familiare e tempo, la scelta dell'allocazione statica è coerente con l'ingenuità dalla quale è nato questo programma, pensato, come vuole tradizione, mediante il paradigma imperativo. Un uso intelligente ed elegante dell'allocazione dinamica prevede una struttura ben diversa, che si può ottenere ragionando con il paradigma "a oggetti".
    L' uso dell'allocazione statica in questo caso è comunque ampliamente giustificabile facendo un bilancio dei costi e dei benefici: i pochi Kbyte allocati più del necessario sono paragonabili a una goccia nel mare di RAM delle macchine odierne; se questo consente di evitare le complesse operazioni di allocazione e deallocazione della memoria, ecco che la scelta si giustifica.
    Le variabili da cui dipenderebbero le dimensioni dei vari arrays sono il numero di nodi della rete (n), l'ordine del sistema non lineare (ord) e il numero di iterazioni da effettuare (numit); gli estremi superiori fissati sono rispettivamente maxnodi, 2*maxnodi emaxnumit, i cui valori sono consultabili nella pagina delle caratteristiche.
    Ecco l'inventario degli arrays sovradimensionati:

    Variabili globali sovradimensionate

    Unit Global

    • 4 arrays double di dimensione [maxnodi][maxnodi] (matrici dati ammettenze nodali)
    • 1 array  double di dimensione [4][maxnodi]       (dati grandezze nodali)
    • 2 arrays int di dimensione  [2*maxnodi]   (arrays di supporto per il calcolo dello Jacobiano)
    • 1 array  int di dimensione  [maxnodi]     (array di supporto tipo di nodo)


    Variabili locali sovradimensionate

    Unit DataIn

    • 1 array  bool  di dimensione  [4][maxnodi]   (matrice conoscenza dati nodali)
    Unit Gauss
    • 2 arrays double di dimensione [maxnodi][maxnumit] (risultati dell'elaborazione)
    Unit Newton
    • 2 arrays double di dimensione [maxnodi][maxnumit] (risultati dell'elaborazione)
    • 2 arrays double di dimensione [2*maxnodi][2*maxnodi] (matrice Jacobiana e sua inversa)
    • 2 arrays double di dimensione [2*maxnodi]         (vettore residui e incrementi)
    Unit Stott
    • 2 arrays double di dimensione [maxnodi][maxnumit] (risultati dell'elaborazione)
    • 4 arrays double di dimensione [maxnodi][maxnodi] (matrici Ba e Br e loro inverse)
    • 2 arrays double di dimensione [2*maxnodi]         (vettore residui e incrementi)


    allocazione dinamica
    In realtà sono stato costretto ad allocare dinamicamente (mediante l'istruzione di C++  new)  le matrici di cui fanno uso gli algoritmi di Newton e Stott per definire il puntatore di puntatore double da passare alla funzione globale di inversione. Per snellezza ho comunque sovradimensionato anche queste matrici, disallocando ( con delete) la memoria solo alla chiusura del programma.

    caratteristica peculiare
    Una caratteristica del programma di cui vado fiero è che gli algoritmi accettano una numerazione arbitraria dei nodi.
    Per farli funzionare senza una numerazione decisa a priori a seconda del tipo di nodo,  ho dovuto utilizzare tre arrays di supporto dei calcoli:
    Uno è l'array tiponodo[], i cui indici (nodi) sono associati a interi che rappresentano il tipo di nodo (0=Ed,1=PQ,2=PE,...); questo array è usato in maniera estensiva, anche per costruire i successivi array di supporto che a loro volta mi sono serviti per la costruzione dello jacobiano,del vettore residui e per sommare gli incrementi delle incognite  (tutto questo fa parte del metodo di Newton) ; parlerò di questi array nel capitoletto di questa stessa pagina riguardante l' algoritmo di Newton.




    Sulla Unit Gauss
    (riguardo l'algoritmo vedi le generalità e come è applicato) Torna Su

    Un pò di storia
    Il problema dell' algoritmo di Gauss, con la numerazione arbitraria dei nodi, è che non è detto che si conosca a priori la sequenza di operazioni da svolgere: ogni tipo di nodo necessita un diverso trattamento e se la numerazione dei nodi non è legata al tipo non è possibile preordinare i passi da svolgere.
    Una prima soluzione per ovviare il problema è quella di non consentire all' utente una numerazione arbitraria, ma secondo un criterio ben preciso: in questo modo si creerebbe un legame tra la sequenza e i tipi di nodo, quindi un modo per prevedere la sequenza di operazioni una volta conteggiato il numero dei vari tipi di nodo. Questa soluzione è banale e ne ho cercata un' altra.
    Scrivendo una funzione dedicata per ciascun tipo di nodo, potrei, ad ogni iterazione, creare un ciclo esteso a tutti i nodi nel quale chiamare per ciascun nodo la funzione corretta mediante un controllo del tipo di nodo, tramite un costrutto switch.
    In questo caso riuscirei ad ottenere una sequenza di operazioni che si adatta a una qualsiasi numerazione dei nodi, con l' unico onere consistente nella iniziale costruzione di un array che lega i tipi di nodo con la sequenza scelta dall' utente..

    La soluzione
    La soluzione che ho adottato è molto vicina a questa appena descritta, ma ho voluto evitare il controllo che si effettuerebbe ogni iterazione per ogni nodo adottando un array di puntatori a funzione *pfunct[] , che associa i tipi di nodo alle sei funzioni associate ai sei tipi di nodo fEd , fPQ , fPE , fPd , fQE , fQd.
    Eccettuata  fEd (che si limita ad aggiornare la tensione col valore già noto), le funzioni puntate da *pfunct[]  hanno il compito di stimare l'eventuale potenza mancante (a seconda del nodo in questione) e calcolare la nuova tensione del nodo usando la formula di GlimmStagg, comune a tutti i nodi e quindi implementata per comodità in una settima funzione, chiamata appunto GlimmStagg .
    Ad ogni iterazione, per ogni nodo viene opportunamente gestito attraverso il meccanismo del puntatore a funzione, senza effettuere alcun controllo del tipo di nodo (a patto di aver precedentemente costruito un array che associa i nodi ai corrispondenti tipi).
    Descrivendo con un disegno (che non è coerente con alcuna convenzione formale) :

    Riassumendo, grazie all' uso intelligente dell' array  tiponodo e di un array di puntatori a funzione, si ottiene la chiamata delle funzioni opportune per ogni nodo p ad ogni iterazione k.
    Andando a vedere il codice sorgente della funzione principale di Gauss-Seidel viene quasi da sorridere tanto è semplice:
    //---------------------------------------------------------------------------
    //algoritmo Gauss
    for(k=0 ; k<numit ; k++)      //ad ogni iterazione
        for(int p=0 ; p<n ; p++)   //per ogni nodo
        pfunct[tiponodo[p]](p);    //chiamata della opportuna funzione
    //Calcolo terminato;conversione delle fasi in gradi...
    //---------------------------------------------------------------------------




    Sulla Unit Newton
    (riguardo l'algoritmo vedi le generalità e come è applicato) Torna Su

    Le difficoltà: l' ordinamento...
    Anche qui, come in Gauss, la numerazione arbitraria dei nodi crea qualche problema, forse qualcuno in più.
    L' algoritmo di Newton-Raphson rappresenta infatti una sfida più interessante rispetto Gauss: qui bisogna individuare le equazioni della rete che insieme costituiscono il sistema di equazioni non lineari da risolvere; Bisogna poi individuare le incognite (il cui numero deve eguagliare l'ordine del sistema); fatto questo bisogna ad ogni iterazione calcolare il residuo del sistema, la matrice dello Jacobiano, invertirla, moltiplicarla con il vettore residuo per ottenere il vettore incrementi il quale deve essere utilizzato per aggiornare le incognite.
    Si presentano due problemi principali: il primo è far corrispondere tra loro gli elementi dello Jacobiano e dei vari vettori che popolano le formule dell' algoritmo; è necessario quindi ordinare opportunamente gli arrays.

    ...e l' inversione
    Il secondo problema, che mi ha fatto perdere un tempo considerevole, è la "banale" inversione di matrice. Ho cercato su internet un' implementazione che mi soddisfacesse, ma non l' ho trovata.Ho avuto anche la tentazione, sentendo parlare di C++, classi, metodi, overload degli operatori, di introdurre una classe "Matrici" con relativi operatori di inversione, somma, prodotto ( magari scaricandola da internet).Ho resistito a questo, adottando la soluzione molto più ingenua di una funzione void operante su un array globale.

    La soluzione per l' ordinamento...
    Ho deciso di fissare l' ordinamento dei vari arrays così come l' ho visto a lezione, ossia voglio pensare allo Jacobiano diviso ordinatamente nei quattro sottoblocchi (vedi qui).
    Coerentemente a tale ordinamento devo costruire il vettore residui in modo che abbia come elementi prima i residui attivi, seguiti da quelli reattivi. Operando in tal modo il vettore degli incrementi che otterrò dai calcoli sarà costituito prima dagli incrementi dei moduli delle tensioni, poi da quelli delle fasi.
    Per effettuare tutte le operazioni in cui è coinvolto l'ordinamento ho fatto uso di due arrays di supporto: uno è chiamato Jeq e associa ai suoi indici (interpretati come la numerazione delle equazioni, ordinate nel modo voluto) al numero del nodo che genera l'equazione; l' altro è Jinc e associa i suoi indici (interpretati come la numerazione delle incognite, ordinate nel modo voluto) al numero del nodo che fornisce l' incognita. Quindi se e è l'indice dell' equazione desiderata, Jeq [e] dà l'indice del nodo che la genera; analogamente se i è l'indice dell'incognita desiderata, Jinc [i] dà l'indice del nodo che la fornisce.
    Jinc [ ] e  Jeq [ ] si costruiscono con l'ausilio del già menzionato array tiponodo [ ].

    ...e altri problemi dell'inversione
    La funzione che implementa l' inversione si chiama CalcoloInversa( ).
    In essa ho realizzato la decomposizione LU e la risoluzione per colonne degli n sistemi lineari che si devono risolvere per determinare le colonne dell'inversa.
    Dopo il considerevole tempo necessario alla codifica, ero molto orgoglioso di questa funzione. Amara è stata la sorpresa nel constatare che la mia bella funzione non produceva affatto risultati corretti pur essendo a prima vista priva di errori semantici. Solo dopo una notte insonne sono giunto a capo del problema: non avevo considerato gli errori dovuti al Pivot, ossia all'avere elementi diagonali della matrice molto più piccoli degli elementi fuori diagonale. Ho dovuto quindi aggiungere una parte di codice per effettuare un pivoting della matrice prima dell 'inversione, ossia una permutazione delle linee atta a scambiare gli elementi diagonali più piccoli con elementi di valore maggiore.
    Credo di aver implementato ciò in maniera piuttosto intelligente, utilizzando (guarda caso) ancora un array di supporto (mi sono convinto che la mia tecnica di programmazione sia fortemente caratterizzata dall'utilizzo di array di interi di supporto).
    Questo array, chiamato ind[ ], lega i suoi indici (che rappresentano gli indici delle righe della matrice non permutata) con gli indici delle righe della matrice dopo la permutazione; nei calcoli successivi basta sostutuire la riga i con la riga ind [i] per considerare la matrice permutata anzichè la matrice non permutata.
    Attuo quindi il pivoting con la costruzione di ind [i], ottenuta confrontando ogni elemento sulla diagonale con gli elementi della stessa colonna successivi:
    //---------------------------------------------------------------------------
    //Pivoting: permutazione delle righe per non incorrere in errori numerici.
    //Costruzione array delle permutazioni
    for (i=0 ; i<ord ; i++)
       for (j=i+1 ; j<ord ; j++)
       if ( fabs(J[i][i]) < fabs(J[ind[j]][i]) )   {
                                                               temp=ind[i];
                                                               ind[i]=ind[j];
                                                               ind[j]=temp;
                                                              }
    //---------------------------------------------------------------------------
    Preciso che nella attuale versione del programma la funzione CalcoloInversa( ) non è dotata di controllo della singolarità della matrice da invertire.
    Si consulti a tal proposito la nota sull'inversione dello Jacobiano nella pagina che descrive gli algoritmi implementati.

    l 'algoritmo
    L 'algoritmo si riduce a una serie di chiamate alle funzioni void che operano sulle variabili globali della Unit:
    //----------------------------------------------------------------------------
    //algoritmo Newton
    for(k=0 ; k<numit ; k++)
       {
       CalcoloResiduo( );        /*calcolo di res*/
       CalcoloJacobiano( );     /*calcolo di J*/
       CalcoloInversa( );          /*calcolo inversa dello jacobiano Jinv*/
       CalcoloIncrementi( );   /*prodotto tra Jinv e res*/

       /*aggiornamento incognite*/
       for(i=0 ; i<nE ; i++)     E [Jinc[i]] [k+1]=E [Jinc[i]] [k]+incr[i];
       for(i=nE ; i<ord ; i++)   d [Jinc[i]] [k+1]=d [Jinc[i]] [k]+incr[i];
       }

    //Calcolo terminato;conversione delle fasi in gradi...
    //----------------------------------------------------------------------------
    Nella parte dell' aggiornamento delle incognite si può notare l' uso dell' array di supporto Jinc [ ], grazie al quale è possibile accedere agli indici dei nodi a partire dagli indici che ordinano le incognite.
    Di seguito uno schema (che non è coerente con alcuna convenzione formale) della Unit Newton, che vorrebbe dare l'idea delle chiamate che effettua il main a ogni iterazione evidenziando i quattro passaggi principali dell' algoritmo:

    visualizzazione
    Per un confronto utile con i calcoli a mano avrei dovuto prevedere la visualizzazione dello Jacobiano, dei vettori dei residui e degli incrementi ad ogni iterazione, ma non l'ho fatto perchè avrei complicato e reso meno comprensibile la finestra di visualizzazione dei risultati diversificandola da quella di Gauss; inoltre, data la complessità dei calcoli, ci sono pochi pazzi che intraprenderebbero il calcolo a mano con questo algoritmo che qui è stato implementato principalmente per uno scopo di verifica e di confronto dei risultati.




    Sulla Unit Stott
    (riguardo l'algoritmo vedi le generalità e come è applicato)  Torna Su

    Difficoltà già superate...
    Questo algoritmo, da buon figlio, eredita da Newton le stesse difficoltà implementative: l'ordinamento e l'inversione di matrice, problemi già risolti.
    La costruzione delle matrice Ba e Br e le varie operazioni matriciali, in cui è fondamentale tenere conto dell'ordinamento, sono codificabili immediatamente riutilizzando gli indispensabili "array indice" Jeq e Jinc, di cui si è già diffusamente parlato.
    L'operazione di inversione è già stata implementata; se prima era dedicata alla Unit Newton, si tratta ora di condividerla con questa Unit.

    ...e nuovi motivi di riflessione
    Tuttavia l'implementazione di Stott ha fatto emergere dei problemi legati alla scelta di accettare i tipi di nodi che io chiamo atipici, quali Pd e QE, scelta suggerita unicamente dal considerare tutte le combinazioni a due a due dei dati PQEd e priva di motivazioni fisiche (che io sappia non so se esistano sistemi reali modellizzati con questi tipi di nodi).
    Questi problemi si traducono nell'inapplicabilità di Stott e nella produzione di risultati senza senso fisico da parte di Newton quando il numero di fasi incognite non coincide con il numero di iniezioni di potenza attiva note e, corrispondentemente, il numero di moduli incogniti non coincide con il numero di iniezioni di potenza reattiva note.
    Stott (come del resto Carpenter) è inapplicabile perchè Ba e Br non risultano più quadrate e quindi non possono essere più invertite; si noti che la matrice quadrata che le contiene risulta singolare, e perciò non invertibile.
    Ciò accade quando sono presenti in numero diverso nodi Pd e QE (infatti un nodo Pd produce un'equazione in P e un'incognita E; esso viene bilanciato solo da un nodo QE, che produce un'equazione in Q e un'incognita d).
    In queste condizioni l'algoritmo di Newton, benchè matematicamente applicabile, produce risultati senza senso.
    Nell'attesa di chiarire questo fatto è stato necessario introdurre dei controlli sugli ingressi prima dell'avvio dei due algoritmi bloccando, nel caso, l'algoritmo di Stott  e applicando l'algoritmo di Newton previo avvertimento.
    L'algoritmo di Gauss, dal canto suo, sembra non essere minimamente toccato da questi ingressi particolari, producendo risultati sensati (ma non verificabili): questo potrebbe suggerire un'origine prettamente matematica del problema, che non sarebbe così generato da una inconsistenza fisica degli ingressi.

    visualizzazione
    Anche se nella attuale versione non è previsto, prossimamente sarà opportuno prevedere la visualizzazione delle matrici Ba e Br.




    Torna all' indice      Torna  Su