Macchine elettriche e Linee sono “descritte” da modelli elettrici costituiti da tre bipoli “ideali”: Induttanza, resistenza e capacità. La resistenza rappresenta tutto ciò che è “attivo”, ovvero dove avviene la conversione energetica. Induttanza e condensatore, invece, rappresentano tutto ciò che “accumula” energia. Le grandezze “energetiche”, ovvero quelle associate allo stato energetico del sistema, sono Variabili di Stato (V.D.S). Dal momento che le V.D.S. non ammettono salti (principio di conservazione dell’energia), qualora si cercasse di far variare qualche parametro, il sistema risponderà con dei transitori (ovvero con una certa dinamica), che garantiscano la continuità delle stesse. Va tenuto quindi presente che non vi è alcuno modo di ignorare la dinamica di un sistema elettrico. La si può rendere veloce quanto si vuole (riesce….) ma il circuito elettrico non sarà mai svincolato dalle leggi della Fisica.

Circuito RL

Dato un circuito RL, ad esempio, la variabile di stato è la corrente.

L’equazione costitutiva dell’induttanza è la seguente:

Le lettere in minuscolo sono funzioni del tempo.

Analiticamente si osserva facilmente che la corrente induttiva non può variare a gradino (essendo V.D.S.) : questo comporterebbe una derivata “virtualmente” tendente all’ infinito, portanto sarà la tensione stessa a tendere “virtualmente” all’ infinito (nella realtà ciò non avviene a causa della presenza di parametri quali capacità parassite ecc) . Una tensione “tendente” ad infinito si traduce in archi elettrici che consentono comunque il passaggio di corrente. L’arco elettrico, di fatto, è il fenomeno che garantisce la naturale continuità della corrente.

L’enegia immagazzinata dall’induttore, infatti, vale:

Se la corrente variasse “a gradino” varierebbe a gradino anche E: questo è incompatibile con la Fisica (l’apertura di un circuito induttivo è “incompatibile”).

Scrivendo l’equazione della suddetta maglia RL, si può esplicitare la variazione di corrente:

Ora, senza passare per la risoluzione analitica risolvendo l’equazione differenziale ed il classico problema di Cauchy, si può ricorrere al calcolo numerico.

Prima di fare ciò, è interessante notare come dalla suddetta equazione sia possibile ritrovare alcuni “concetti” fisici. Qualora l’inerzia del circuito, ovvero L, sia nulla, l’equazione presenterebbe una singolarità: La derivata di corrente “tenderebbe” ad infinito, ovvero si avrebbe un vero e proprio gradino di corrente. Si è già detto che questo, in un circuito induttivo, non è possibile. Il punto è che, se L fosse nulla, il circuito non sarebbe induttivo ma si tradurrebbe in un circuito resistivo (sola R). La corrente non sarebbe più variabile di stato e quindi sarebbero consentite variazioni in tempo zero.

Se si analizzasse, ad esempio, la soluzione analitica con ingresso di tensione a gradino (ovvero chiudendo l’interruttore S) :

Ove:

Si osserverebbe che l’assenza di inerzia porterebbe ad una costante di tempo nulla e quindi :

Qualora non vi fosse la resistenza R, l’equazione presenterebbe singolarità: La corrente tenderebbe, nel tempo, virtualmente ad infinito (in sostanza si tenderebbe ad un corto circuito).

Scrivendo il modello della maglia in Python per poi risolverlo con il solver “odeint()”, si può facilmente verificare che si tratta di un transitorio del primo ordine:

import numpy as np
 from scipy.integrate import odeint
 import matplotlib.pyplot as plt
 import math 
 R = 0.3 #[OHM]
 L = 0.07 #[H]

 def RL(i,t,e):
    di_dt = (1/L)*(e-R*i)
    return di_dt

 n = 10000
 t_max =2
 t = np.linspace(0,t_max,n)
 i0 = 0
 e =3

 i= odeint(RL,i0 ,t,args=(e,))

 plt.plot(t,i, color ="red") 
 plt.xlabel('Tempo [s]')
 plt.ylabel('i(t) [A]')
 plt.axis([0,7*L/R, 0, 10.5])
 plt.grid()
 plt.show() 

L’andamento della corrente (variabile di stato) in funzione del tempo è il seguente.

La tensione ai capi dell’induttore può essere ricavata risolvendo la maglia:

Notare che la tensione ai capi di L si annulla asintoticamente. Questo è coerente dal momento che, qualora non fosse così, la corrente crescerebbe ininterrottamente (la motivazione è già stata espletata).

La potenza fornita dal generatore viene in parte dissipata sulla resistenza, in parte viene trasformata in energia magnetica ed immmagazzinata sull’indottore L:

In nero la potenza dissipata sulla resistenza e in verde la potenza di magnetizzazione. L’area sottesa alle curve (per definizione) è l’energia.

Chiudendo un induttore carico (quindi sta scorrendo una corrente non nulla) su una resistenza si innescerebbe un’evoluzione libera, ovvero senza forzante (generatore).

Si osserva come al variare della resistenza si ha una variazione della costante di tempo (ovviamente): Più è alta la resistenza più velocemente viene dissipata l’energia (quindi il transitorio sarà più veloce). In sostanza, qualora la resistenza fosse nulla si avrebbe una costante di tempo tendente “virtualmente” ad infinito. Questo si traduce nel fatto che il transitorio non si estinguerebbe mai (l’energia non verrebbe dissipata):

Tracciando quindi l’andamento della corrente al crescere della resistenza si ottiene il seguente risultato:

Plottando, invece, l’andamento della potenza dissipata sulla resistenza, si osservano che l’energia dissipata (l’area sottesa alla potena) è la medesima ma la velocità con cui viene dissipata (che dipende dalla costante di tempo del sistema) varia:

Per alti valori di R, l’energia viene dissipata più lentamente.

Per concludere questo ragionamento, il valore di resistenza R è associata alla dissipazione energetica, mentre il valore dell’induttanza L rappresenta porpio la capacità di immagazzinamento di energia.

Circuito RC

Per dualismo, si tenga in considerazione un circuito capacitivo RC:

Lo stato energetico è associato alla tensione sul condensatore vc (che è v.d.s.). L’energia immagazzinata in un condensatore, infatti, vale:

Azzerando la resistenza, quando si chiude l’interruttore e la tensione vc non è pari alla tensione di alimentazione V, si avrebbe un corto circuito. Questo perchè la tensione su un condensatore è variabile di stato (e quindi non ammette variazioni a gradino). Questo è osservabile risolvendo l’equazione differenziale:

La soluzione analitica è :

Al tempo t = 0 si avrebbe:

La caduta sulla sulla resistenza è pari a (istante per istante):

All’istante iniziale tutta la tensione cade su R: qualora fosse assente si avrebbe un corto circuito per antonomasia (sostanzialmente la corrente tenderebbe, all’istante iniziale, virtualmente ad inifinito).

Il modello in Python:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math 
R = 0.3 #[OHM]
C = 0.1 #[F]

def RC(vc,t,e):
   dvc_dt = (1/(R*C))*(e-vc)
   return dvc_dt

n = 10000
t_max =2
t = np.linspace(0,t_max,n)
e =3
vc0 = 0
vc= odeint(RC,vc0 ,t,args=(e,))
 
plt.plot(t,vc, color ="blue") 
plt.xlabel('Tempo [s]')
plt.ylabel('vc(t) [V]')
plt.axis([0,7*R*C, 0,3.5])
plt.grid()
plt.show()  

Da cui:

Osservazioni

Si può quindi affermare che:

• Se all’instante della chiusura di un interruttore su un condensatore la tensione su di esso non è pari alla tensione di alimentazione, si ottiene un corto circuito;
• Se all’istante di apertura di un interruttore in un ramo induttivo non si consente il ricircolo della corrente ma si “pretenderebbe” una discontinuità, si ottengono delle sovratensioni tali da danneggiare gli isolanti tramite Archi elettrici.

Circuito RLC

In un circuito RLC si instaura un transitorio del seondo ordine, ovvero si hanno due variabili di stato: la corrente su L e la tensione ai capi di C.

In Python, una funzione che descrive il suddetto circuito è il seguente:

def RLC(v_d_s,t,e):
     i = v_d_s[0]
     vc = v_d_s[1]
     di_dt = (1/L)*(e-R*i-vc)
     dvc_dt = (1/C)*i
     v_d_s = [di_dt, dvc_dt ]
     return v_d_s 

Applicando un gradino di tensione, le due variabili di stato (partendo da una condizione di energia nulla, ovvero sia tensione che corrente iniziale sono nulle) hanno il seguente andamento:

Quindi, all’aumentare del valore della resistenza, si ha un aumento dello smorzamento (c.v.d.). Il caso limite lo si avrebbe con resistenza nulla, ovvero il transitorio non si estinguerebbe mai:

Passivando il generatore di tensione e indicando quali sono le condizioni iniziali (corrente sull’induttore e tensione sul condensatore all’istante zero), si ottiene l’evoluzione libera:

Qualora il circuito fosse alimentato in alternata, si oserverebbe un rimbalzo di energia tra le due reattanza, oltre che ad una dissipazione attiva sulla resistenza. Si ha una perfetta anaolgia con il Pendolo (meccanica classica):

Nella posizione 1 la velocità è nulla è l’energia posseduta dalla massa m è solo quella potenziale (che rappresenta la tensione o “potenziale elettrico”). Lasciato libero di cadere aumenta la propria velocità trasformando l’energia potenziale in energia cinetica . Nel punto 2, energia potenziale nulla (rispetto al riferimento scelto) e la velocità sarà massima (massima energia cinetica). Ora, se non vi è smorzamento dovuto agli attriti (componente “resistiva” che dissipa energia) dal punto 2 tornerà alla medesima altezza del punto 1: pendolo ideale. Se invece lo smorzamento è presente, oscillera fino a fermarsi in modo “smorzato”: tutta l’energia posseduta nel punto 1 è stata dissipata in calore (come nella resistenza).

Quindi, se è presenta attrito il pendolo prima o poi si fermerà (a seguito di oscillazioni smorzate). L’analogia con l’ elettrotecnica è completa: senza resistenza l’energia viene scambiata dalle reattanze ma non si trasforma in calore. La presenza della resistenza, invece, fa si che l’energia “esca” dal circuito sottoforma di calore (effetto Joule).

Matteo Gentileschi

Calcolo dei flussi di potenza

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

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

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

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

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

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

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

Facendone l’inversa:

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

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

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

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

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

Il nodo i-esimo presenta la seguente struttura:

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

Notare che, in notazione matriciale:

Quindi, si può scrivere:

Il BFS è costituito da due fasi:

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

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

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

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

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

In python il ciclo è il seguente:

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

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

Il codice completo è il seguente:

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

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

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

V , Is = BFS(Zci)

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

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

Il risultato, ad esempio, è il seguente:

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

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

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

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

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

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

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

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

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

La funzione BFS diventa:

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

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

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

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

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

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

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

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

BFS probabilistico

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

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

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

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

L’algoritmo viene così sintetizzato:

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

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

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

La CDF sarà la seguente:

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

Matteo Gentileschi

Quando si ha a che fare con un sistema fisico si deve tener presente che esso può “instabilizzarsi”. Si può affermare che uno dei parametri che quantizzi al meglio e “a colpo d’occhio” quanto tale sistema sia distante dell’instabilità è il Fattore di Sicurezza. Esso può essere definito, ad esempio, come rapporto tra “forze resistenti” e “forze instabilizzanti”:

Se vuol dire che il denominatore è maggiore del numeratore, ovvero che . In altri termini, le forze instabilizzanti sono superiori a quelle che si oppongono a tale instabilizzazione. Pertanto, a seconda dell’incertezza dei parametri, è bene mantenere tale rapporto superiore all’unità.

Notare che sarebbe possibile utilizzare anche il reciproco di tale rapporto ma, a quel punto, si deve puntare ad avere un valore inferiore all’unità (complementarità).

Tale concetto è applicabile a più ambiti. Si consideri, ad esempio, un piano inclinato.

La forza instabilizzante Fi è la componente della forza peso parallela alla superficie di scivolamento. Tale forza sarà contrastata dalla forza di attrito che assumerà le vesti di forza resistente Fr.

Tralasciando la notazione vettoriale (essendo la direzione esplicitata dal disegno), note la massa m del corpo e la geometria del problema si ha:

La forza resistente Fr sarà data dal prodotto tra la componente P (normale alla superficie) e il coefficiente di attrito statico.

La componente P viene bilanciata dalla reazione vincolare offerta dal piano (se la struttura è sufficientemente robusta, altrimente cede). Il moto lungo la direzione tangenziale del piano avviene se sussiste una forza dinamica, ovvero se la differenza tra la forza ressitente e la forza instabilizzante non è nulla. Il corpo “resterà fermo” se:

Ovvero l’angolo critico (oltre il quale avviene il moto) sarà quello che soddisfa la seguente condizione:

Rammentando la definizione di Fattore di sicurezza Fs, applicandola al caso in esame, si ha che:

Si osservi che con Fs = 1 (come accennato sopra, equilibrio perfetto tra le forze) si ottiene nuovamente la relazione che definisce l’angolo critico:

A questo punto è d’uopo un chiarimento. Qualora il fattore di sicurezza fosse maggiore di uno, matematicamente ne risulterebbe che la forza di attrito è maggiore della componente tangenziale. Questo NON corrisponde alla realtà fisica. Se così fosse, il blocco inizierebbe a muoversi verso l’alto. Prima di contattare le autorità ecclesiastiche, si rammenti che la forza di attrito, definita come forza normale per un coefficiente di attrito statico, è quella “massima”. Supponento infatti che l’angolo sia nullo.

Per muovere il blocco è necessario comunque applicare una forza dall’esterno che superi F_r, ma se ciò non avviene non si ha moto.

Tale risultato lo si ritrova plottando le componenti della forza peso (supposta, per questo calcolo, unitaria) al variare di theta, scrivendo il seguente codice:

import numpy as np
import matplotlib.pyplot as plt
import math
n= 100 #numero di campioni
us = 0.61 #coefficiente attrito statico
theta = np.linspace(0, math.pi ,n) #angolo theta_ da 0 a 90°
Fp = 1 #Forza peso
#vettori delle componenti al variare di theta
Fi = np.zeros(n)
Fr = np.zeros(n)
#calcolo componenti
for i in range(0,n):
    Fi[i] = Fp * math.sin(theta[i])
    Fr[i] = us*Fp * math.cos(theta[i])
plt.plot(theta,Fi, label="Fi")
plt.plot(theta,Fr, label="Fr")
plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc='lower left', ncol=2, mode="expand", borderaxespad=0.)
plt.xlabel('Angolo')
plt.ylabel('Forze')
plt.grid()
plt.show()

Le curve si incrociano proprio nel valore di theta di 0,54rad: oltre il quale, la forza resistente Fr non riuscirà più a contrastare la forza motrice Fi.

Quindi, idealmente, fino a che l’inclinazione non supera tale valore, non si corre il rischio che il blocco scivoli via. Il problama “reale” è che i suddetti valori non sono noti senza incertezza.

Metodo Monte Carlo per la stima del fattore di sicurezza

Fino a che si è in un contesto deterministico, ovvero il caso ideale, quanto detto è sufficiente. Nella realtà si hanno delle “incertezze” dei parametri misurati. Essi infatti non saranno definiti da uno scalare ma da una distribuzione probabilistica.

A seguire un esempio in cui i valori numerici sono scelti di autorità.

Si supponga di disporre un blocco di Alluminio su un piano di acciaio, quindi, idealmente :

Da un semplice calcolo si ha che l’angolo limite vale:

Quindi, per sicurezza,si predispone un blocco che non consenta di andare oltre i 0,5 rad. Bisogna ora verificere se tale “azione per la sicurezza” è sufficiente.

Il problema è che il dispositivo di regolazione dell’inclinazione presenta una tolleranza del 5%. Il piano, inoltre, potrebbe essere sporco, cosa che potrebbe variare il coefficiente di attrito statico del 20%.

Partendo dalle distribuzioni probabilistiche cumulative (C.D.F.) , sia del coefficiente di attrito che dell’angolo:

Dalle quali, partendo dall’estrazione casuale di un valore da una distribuzione probabilistica uniforme (0-1), si estrarranno i due valori necessari. Per ogni estrazione si calcolerà il fattore di sicurezza (relativo all’estrazione i-esmima):

Procedendo quindi con un numero di estrazioni pari a n=1000, si “contano” le volte in cui è minore ad un certo valore limite che per questioni cautelative poniamo pari a 1.05, 5% in più del necessario, per poi calcolare la frequenza con cui ciò accade.

Il codice in Python:

import numpy as np
import matplotlib.pyplot as plt
import math 
import random 

#Valori centrali
us = 0.61 #Coefficiente di attrito statico 
theta = 0.5 #angolo di inclinazione [rad]

tol_theta= 0.05 #tolleranza sull'angolo
tol_u = 0.2 #tolleranza sul coefficiente di attrito

theta_min = theta(1-tol_theta) 
theta_max = theta(1+tol_theta)
us_min = us(1-tol_u) 
us_max = us(1+tol_u)

K = 300 #numero di passi
#Parametri distribuzione normale 
sigma_us = 0.025 #Deviazione standard (sqrt)
sigma_theta = 0.007 #Deviazione standard 
(sqrt)


#Funzione per la costruzione della PDF
def PDF(x,sig,xc,k):#(x, sigma, valore centrale, passi)
     a = 1/(pow(2*np.pi*pow(sig,2),0.5))
     pdf = np.zeros(k)
     for i in range(0,k):
         pdf[i] = a*math.exp(-0.5*pow(x[i]-xc,2)/(pow(sig,2)))
     return pdf

#Funzione per la costruzione della CDF
def CDF(pdf,k):
     cdf = np.zeros(k)
     cdf[0] = pdf[0]
     for i in range(1,k):
         cdf[i] = cdf[i-1]+pdf[i]# Normalizzazione C.D.F.
     cdf = cdf /max(cdf)
     return cdf

#Costruzione PDF angolo
vec_theta = np.linspace(theta_min,theta_max,K)#vettore angoli
PDF_theta = PDF(vec_theta,sigma_theta,theta,K)

#Costruzione PDF "u"
vec_us = np.linspace(us_min,us_max,K)#vettore coefficienti
PDF_us = PDF(vec_us,sigma_us,us,K)

#Costruzione CDF 
CDF_theta = CDF(PDF_theta,K)
CDF_us = CDF(PDF_us,K)
#Estrazione valore da DISTRIBUZIONE UNIFORME

#Funzione per l'estrazione (distribuzione uniforme 0-1)
def dado(): return np.random.uniform(0,1)#Random

#Funzione di estrazione
def estrazione(cdf, k):
    r = dado()
    if(r<= cdf[0]):
        r_v = 0 
    for i in range(1,k):
        if(r<= cdf[i] and r > cdf[i-1]):
            r_v = i
     return r_v

 #ESTRAZIONI
 n = 1000 #numero di estrazioni
 Fs_min = 1.05
 f = np.zeros(K)
 for i in range(0,K):
     f[i] = Fs_min*math.tan(vec_theta[i])
 cont = 0
 for i in range(0,n):       
     us_es = estrazione(CDF_us,K)#Estrai u
     theta_es = estrazione(CDF_theta,K)#Estrai theta
     #Plotta il punto (coppia di valori) estratto
     plt.scatter(vec_theta[theta_es],vec_us[us_es],color ="red",s = 1) 
     #Verifica se il Fsi è inferiore al valore limite
     if(vec_us[us_es]/math.tan(vec_theta[theta_es])< Fs_min):
         cont += 1 #Incrementa contatore (Fs <Fs_min)


 plt.plot(vec_theta,f,color ="black") 
 plt.ylabel(" coefficiente di attrito")
 plt.xlabel("Angolo")
 plt.grid()
 plt.show()

 Prob = 100*cont/n #Calcola probabilità
 risultato = "La probabilità che Fs sia inferiore a " + str(Fs_min) + " è del: " + str(Prob) + "%"
 print(risultato)  

Il risultato ottenuto:

“La probabilità che Fs sia inferiore a 1.05 è del: 9.0%”

Quindi, seppur l’angolo è stato impostato al di sotto della soglia limite, ovvero quella che “matematicamente” garantiva la staticità del blocco, dal Metodo Monte Carlo si evince che si ha comunque un rischio elevato.

Ciò che si può fare, ora, è ridurre il valore di angolo massimo impostabile. Provando con 0,45 rad (quindi cambiando solo il valore centrale della P.D.F. di partenza):

“La probabilità che Fs sia inferiore a 1.05 è del: 0.1%”

Risultato che consente di lavorare in sicurezza. Si rammenta che, qualora fosse possibile, sia per questioni etiche che per Legge, qualora fosse possibile prendere ulteriori provvedimenti che consentano di ridurre tale probabilità, questi dovranno essere messi in campo.

Si rammenta che il rischio è definito come prodotto tra la magnitudo del “danno” che l’evento nefasto causerebbe e la probabilità che ciò accada. Seppur 0,1% può risultare un valore irrisorio, gli effetti di un blocco di qualche quintale che scivola via sono piuttosto ingenti, pertanto il rischio rimane elevato.

Matteo Gentileschi

Quando il numero di Permutazioni è eccessivo, oppure si hanno variabili di difficile stima, studiare un fenomeno in maniera analitica è praticamente impossibile. Per tale motivo si usa il “Metodo Monte Carlo“.

Si prenda ad esempio il tiro con l’arco. Conoscendo il peso e la forma della freccia, la quantità di moto iniziale ecc è possibile (grazie alla Fisica Classica) calcolarne il punto di impatto. Tra i parametri da considerare, però, vi sono l’emotività del tiratore, la stanchezza, le tolleranze dell’arco e delle frecce le condizioni atmosferiche (vento) ecc: lo studio analitico è piuttosto complesso.

Un fenomeno che dipende fortemente dalle condizioni iniziali ( si veda la “Teoria del caos“) può essere studiato con approccio statistico.

Si considerano quindi le variabili come “Aleatorie“. Un fenomeno aleatorio, a differenza di quello deterministico, non porta sempre allo stesso risultato. Il lancio di un dado ne è l’esempio per antonomasia: il risultato varia di volta in volta.

Il dado, infatti, “tira fuori” un numero secondo una distribuzione probabilistica (P.D.F.) uniforme , ovvero quell’elemento ha la stessa probabilità degli altri di essere selezionato e non dipende da chi tira, da come tira ECC.

Problemi legati alla generazione di “valori casuali“

Si tenga presente che la generazione di un “valore” casuale non è così scontata. Un processore infatti può compiere miliardi di operazioni al secondo ma solo “di tipo algebrico. Il programmatore gestisce tale potenza di calcolo impartendo “istruzioni”. Non è possibile quindi, concettualmente, che il processore fornisca un numero casuale!

Daltronde non esiste un sistema più deterministico di un processore: con gli stessi dati in ingresso restituisce (return) sempre gli stessi risultati.

Molti software, infatti, sfruttano degli algoritmi ma tirano fuori comunque valori riproducibili: si tratta quindi di un processo “Pseudorandom“. Per avere numeri realmente casuali è necessario ricorrere a “Generatori hardware di numeri casuali“, basati su “rumori” fisici (rumore termico , decadimento raioattivo ecc) che, proprio per loro natura, sono casuali.

Quelli utilizzati in questa sede sono dunque dei ” PRNG” ( pseudorandom number generator ).

Distribuzione probabilistica P.D.F. e relativa cumulativ C.D.F.

Esistono diverse distribuzioni probabilistiche. Come esempio si può richiamare la distribuzione probabilistica di Gauss, detta anche “normale”, che viene descritta come:

Ove è il e è la “varianza“.

Una semplice funzione in Python implementa tale funzione:

#Costruzione PDF
 def PDF(x,sig,xc,k):#(x, sigma, valore centrale, passi)
     a = 1/(pow(2*np.pipow(sig,2),0.5))
     pdf = np.zeros(k)
     for i in range(0,k):
         pdf[i] = a*math.exp(-pow(x[i]-xc,2)/(2*pow(sig,2)))
     return pdf

Per questioni di calcolo numerico, si è reso necessario “quantizzare” la P.D.F., ovvero la variabile x è stata discretizzata (in K passi). L’integrale della P.D.F. porta alla distribuzione probabilistica cumulativa C.D.F.. Nel discreto, un integrale si riduce ad una semplice sommatoria:

#Costruzione CDF
 def CDF(pdf,k):
         cdf = np.zeros(k)
         cdf[0] = pdf[0]
         for i in range(1,k):
             cdf[i] = cdf[i-1]+pdf[i]
         cdf = cdf/max(cdf)#Normalizzazione CDF
         return cdf

Lancio del dado

Il lancio del dado, o meglio i lanci (servono migliagli di lanci), verrà simulato grazie al modulo “random” di Python (si tratta quindi di PRNG).

import random
#Estrazione da distribuzione uniforme 0-1
def dado(): return np.random.uniform(0,1)
[…]
r = dado() #lancio del dado 
[…]

“Alea iacta est!”

Quindi si vuole avere una distribuzione di probabilità uniforme, ovvero ogni “faccia” ha la medesima probbilità di essere estratto (1/6).

Nota: si tratta di distribuzioni quantizzate. I numeri estratti, ovviemente, sono solo interi. Non è possibile lanciare un dado e ottenere 1,3….

Il modello di distribuzione è riportato in seguito.

facce = [1,2,3,4,5,6]#Facce del dado
PDF =(1/6)* np.ones(6)

A cui è associata una distribuzione probabilistica cumulativa C.D.F.

Notare che la C.D.F. deve avere come valore massimo 1: questo perchè, data la distribuzione probabilistica (P.D.F.), la somma delle probabilità che un evento accada deve coincidere con l’unità. In sostanza, si ha il 100% di probabilità che un valore venga estratto (lanciando il dato non è concepibile che rimanga fermo su uno spigolo).

Il dato di ingresso è il risultado dell’estrazione di una variabile aleatoria a distribuzione uniforme r. Entrando nella C.D.F. si ottiene il relativo valore estratto.

def estrazione(cdf):
  r = dado() #estrazione
  if(r<= CDF[0]):
     p = 0
  for i in range(1,6):
     if(r<= CDF[i] and r > CDF[i-1]):
         p = i
  return p

In p viene quindi memorizzato l’indice della faccia estratta.

Ad esempio, se venisse estratto il valore 0,6 il risultato sarebbe “4”, OPPURE SE VENISSE ESTRATTO IL VALORE 0,25 IL RISULTATO SAREBBE “2”:

Se si procedesse al’estrazione di un numero virtualmente infinito di valori e si calcolasse, grazie ai risultati, la probabilità che un evento avvenga, si otterrebbe proprio la P.D.F. Ad esempio, si conta il numero delle volte in cui è stao estratto il numero 4 e si divide per il numero N di estrazioni.

METODO MONTE CARLO - n iterazioni
contatore = np.zeros(6)
for i in range(0,n):
   estratto = estrazione(CDF)#estrazione da CDF
   contatore[estratto] +=1
pid_r = contatore/n#valore probabilità "ricostruita"

Il codice completo:

import numpy as np
import matplotlib.pyplot as plt
import math
import random

n = 10000 #numero estrazioni
#Ritorna il valore estratto
def dado(): return np.random.uniform(0,1)
#estrazione di un valore contenuto in v in funzione della
#relativa CDF
def estrazione(cdf):
  r = dado() #estrazione
  if(r<= CDF[0]): 
      p = 0 
  for i in range(1,6): 
      if(r<= CDF[i] and r > CDF[i-1]):
         p = i
  return p

facce = [1,2,3,4,5,6]#Facce del dado
PDF =(1/6)* np.ones(6)

#Costruzione CDF
CDF = np.zeros(6)
CDF[0] = PDF[0]
for i in range(1,6):
   CDF[i] = CDF[i-1]+PDF[i]

METODO MONTE CARLO - n iterazioni
def monte_carlo(n):
   contatore = np.zeros(6)#azzeramento contatore
   for i in range(0,n):
      estratto = estrazione(CDF)#estrazione da CDF
      contatore[estratto] +=1 #se estratto, incrementa!
   return contatore/n #valore probabilità "ricostruita"

pid_r = monte_carlo(n)#chiama l'algoritmo "Monte Carlo" 

Il risultato:

Questo risultato è stato ottenuto con un numero di estrazioni n = 100000. Riducendo “n”, la P.D.F. RICOSTRUITA approssimerà sempre meno quella di partenza. Di seguito alcuni esempi.

n=1000

n = 100

n = 10

Daltronde non ci si potrebbe aspettare diversamente. Portando all’estremo il ragionamento tentando un solo lancio, il valore ottenuto “mostrerebbe” una probabilità del 100%. Questo sarebbe assurdo.

n = 1

Inoltre, con “n” eccessivamente basso, il risultato varia continuamente ad ogni simulazione. Si prenda, ad esempio, la faccia numero 4.

[...] C.D.F calcolata a monte [...]
#METODO MONTE CARLO - n iterazioni
def monte_carlo(n):
   contatore = np.zeros(6)
   for i in range(0,n):
      estratto = estrazione(CDF)#estrazione da CDF
      contatore[estratto] +=1
   return contatore[4]/n #probabilità del solo numero "4"!

k = 100 #numero simulazioni per ogni valore di n
P_n_10000 = np.zeros(k)
P_n_1000 = np.zeros(k)
P_n_100 = np.zeros(k)
P_n_10 = np.zeros(k)
P_n_1 = np.zeros(k)
#nota: monte_carlo(), in questo caso, non ritorna "contatore/n"
#ma solamente "contatore[4]/n"
for i in range(0,k):
   P_n_10000[i] = monte_carlo(10000)
   P_n_1000[i] = monte_carlo(1000)
   P_n_100[i] = monte_carlo(100)
   P_n_10[i] = monte_carlo(10)
   P_n_1[i] = monte_carlo(1)

Notare come la probabilità calcolata subisca maggiori variazioni, durante i vari tentativi, per valori di n bassi. Aumentando n, infatti, si osserva che la differenza tra il risultato ottenuto e quello reale diminuisce (varianza) e le curve tendono ad appiattirsi.

Con questo si è voluto dimostrare che il “singolo dato” è decisamente inutile ai fini statistici.

Calcolo della probabilità con il Metodo Monte Carlo

Ipotizzando di avere due valori conosciuti però non in maniera deterministica ma con tolleranza. Per esempio A = 3 e B = 4. Di questi due valori si conoscono le distribuzioni probabilistiche C.D.F. (cumulative):

Volendo calcolare quale sia la probabilità che, estraendo i valori A e B dalle rispettive cumulative, il loro prodotto C rispetti la condizione:

Si simulano n lanci (es. n = 100), come in esempio:

#Funzione costruzione PDF
def PDF(x,sig,xc,k):#(x, sigma, valore centrale, passi)
     a = 1/(pow(2*np.p*ipow(sig,2),0.5))
     pdf = np.zeros(k)
     for i in range(0,k):
         pdf[i] = a*math.exp(-pow(x[i]-xc,2)/(2*pow(sig,2)))
     return pdf
#Funzuone costruzione CDF
def CDF(pdf,k):
         cdf = np.zeros(k)
         cdf[0] = pdf[0]
         for i in range(1,k):
             cdf[i] = cdf[i-1]+pdf[i]
             cdf = cdf/max(cdf)
         return cdf

def dado(): return np.random.uniform(0,1)#Ritorna il valore estratto
 #estrazione di un valore contenuto in v in funzione della
 relativa CDF
 def estrazione(cdf, k):
     r = dado()
     if(r<= cdf[0]): 
        r_v = 0  
     for i in range(1,k):  
        if(r<= cdf[i] and r > cdf[i-1]):
             r_v = i
     return r_v

 PDF_A = PDF(x,sigma_A, A, K)
 PDF_B = PDF(x,sigma_B, B, K)

 CDF_A = CDF(PDF_A,K)
 CDF_B = CDF(PDF_B,K)

 plt.plot(x,CDF_A ,color = "red",label =" C.D.F. A") 
 plt.plot(x,CDF_B ,color = "blue",label =" C.D.F. B") 
 plt.ylabel("C.D.F.")
 plt.xlabel('x')
 plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc='lower left', ncol=2, mode="expand", borderaxespad=0.)
 plt.grid()
 plt.show()

 n = 1000
 cont = 0
 for i in range(0,n):       
     a = estrazione(CDF_A,K)
     b = estrazione(CDF_B,K)
     plt.scatter(x[b],x[a],color ="red",s = 1.2) 
     c = x[b]*x[a] 
    if(c < 1.5*A*B and c > 0.9*A*B):
         cont += 1 #Volte in cui si verifica la condizione richiesta

plt.plot(x,0.9AB/x ,color = "black",label ="ab = 0.9AB")  plt.plot(x,1.5AB/x ,color = "blue",label ="ab = 1.5AB") 
 plt.ylabel("A")
 plt.xlabel("B")
 plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc='lower left', ncol=2, mode="expand", borderaxespad=0.)
 plt.grid()
 plt.show()
 Prob = 100cont/n risultato = "La probabilità che 0.9*AB < C < 1.5*A*B  " + str(Prob) + "%"
 print(risultato)  

La nuvola di punti è costituita dalle coppie (a,b) di valori estratti. Il rapporto tra i punti compresi tra le due curve (condizioni limite) e il numero di punti totali ci da la probabilità che la suddetta condizione si manifesti.

Il risultato (output) è il seguente:

“La probabilità che 0.9AB < C < 1.1AB è del 9.79%”

Applicazioni

Una delle applicazioni del Metodo Monte Carlo è la progettazione e lo studio delle reti elettriche, in particolare in presenza di generazione locale.

Considerando infatti un generatore eolico, la potenza prodotta è proporzionale alla velocità del vento, che a sua volta è un parametro aleatorio. Considerando un certo periodo dell’anno ad un certo orario, si ha una C.D.F. determinata dai dati ottenuti da prolungate campagne di misura:

Il costruttore fornisce i dati necessari a stimare la potenza generata in funzione della velocità del vento.

Quindi, supponendo di avere una rete che coinvolga differenti generatori , per determinare le correnti che passano nei vari tratti di linea, si procede con il Metodo Monte Carlo:

Si imposta un numero di iterazioni sufficientemente alto (esempio n =1000) e per ogni iterazione si estrae un valore di r (da una distribuzione probabilistica uniforme) a cui è univocamente associata (grazie alla C.D.F.) una velocità del vento.

Un calcolo più accurato, deve tener conto della disponibilità del generatore. Si ha probabilità che sia disponibile e che non lo sia: grazie alla relativa C.D.F., estraendo r si ottiene la disponibilità del generatore.

Quindi, ricapitolando, si “estrae” la disponibilità e poi si “estrae” la velocità del vento (ergo, la potenza generata). Questo in ogni iterazione e per ogni generatore, per poi calcolare le correnti sulle linee (tenendo presenti i carichi ecc. con il metodo del load flow).

La “media” di tutte le correnti risultanti è la corrente con cui si dimensionano o si verificano le linee.

Il metodo trova utilità in ambito di “sicurezza“: per calcolare il fattore di sicurezza Fs per quanto riguarda la caduta massi, ad esempio, i dati in ingresso sono “aleatori”: peso del masso, inglinazione delle superfici di scorrimenti, angolo di attrito, attrito delle superfici ecc. Tutti parametri che deve fornire un geologo. La “fisica” è nota ma i parametri e le condizioni iniziali sono dunque casuali. Si procede a numerore estrazioni di tali valori per poi calcolare ogni volta il risultato. Doi numerosi risultati si ricava la “freauenza” con cui il risultato viola i parametri di sicurezza, ovvero (nel complesso) si va sotto il noto “fattore di sicurezza”.

Matteo Gentileschi

Modello analitico in python

Questo articolo è ridondante dal momento che si è parlato di come risolvere circuiti elettrici tramite il calcolo analitico dell variabili di stato, utilizzando però Matlab. In questa sede verrà riproposto il medesimo approccio ma utilizzando uno strumento open source, ovvero Python.

Una volta importate le apposite librerie, si crea il modello della funzione e si chiama la funzione odeint.

Si consideri un circuito con mutui accoppiamenti (quindi una prima approssimazione di un trasformatore) come quello in figura:

Le cui equazioni sono le seguenti:

Si osserva subito che:

Il cambio di segno nella seconda equazione è dovuta alla convenzione di segno (fisico) adottata. Mettendo a sistema le suddette equazioni, si riesce a scrivere:

In sostanza si è arrangiato il sistema per ottenere la seguente forma:

Il modello viene così scritto in Python:

def trafo(Z,t,e):
   r1 = R1
   rl = RL
   l2 = L2
   l1 = L1
   m = M
   a = l1-(pow(m,2)/l2)
   b = m*rl/l2 
   #Condizioni iniziali (al tempo [i]) 
   i1 = Z[0] 
   i2 = Z[1]
   #Sistema di equazioni differenziali 
   di1_dt = (b/a)i2-(r1/a)*i1+e/a
   di2_dt = (- (rl/l2)-m*b/(l2*a))*i2+(m*r1/(l2*a))i1-(m/(l2*a))*e
   #Vettore soluzioni del sistema
   Z =[di1_dt ,di2_dt] 
   return Z

Le variabili a e b sono variabili di comodo (per semplificare la scrittura delle equazioni). I parametri del circuito sono stati dichiarati come variabili globali (in maiuscolo) e richiamati all’interno del modello come variabili locali (in minuscolo).

Il solver ODE viene fatto lavorare all’interno di un ciclo for: all’iterazione i-esima viene fornito il relativo valore della forzante e(t) (che quindi diventa e(i))) e le nuove condizioni iniziali che saranno i valori di i1 ed i2 all’terazione “i-1” , ovvero il risultato dell’iterazione precedente. Notare che il contatore “i” non parte da 0 ma da 1: questo perchè, all’istante 0 i valori di i1 ed i2 sono forniti dall’utente (condizioni iniziali al tempo 0).

#condizioni iniziali al tempo 0
z0 = [0,0]

la funzione forzante e il vettore “tempo” vengono generati a monte della suddetta iterazione, in funzione nel numero “n” di iterazioni:

t = np.linspace(0,t_max,n)
e = np.zeros(n)
for i in range(1,n):
    e[i] = 5 + V*math.sin(Wt[i])

Quindi si procede con l’iterazione:

for i in range(1,n):
   #passo temporale
   tspan = [t[i-1],t[i]]
   #risolvi step
   z = odeint(trafo,z0,tspan,args=(e[i],))
   i_1[i] = z[0][0]
   i_2[i] = z[0][1]
   #Aggiorna le "condizioni iniziali"
   z0 = z[1]

Notare che in tspan viene definito l’intervallo di tempo nel quale si vuole risolvere il sistema: in altre parole, questo approccio, consente una risoluzione a tratti. Aumentado il numero passi “n”, mantenendo invariato il tempo totale di cacolo (t_max), si otterranno risultati migliori.

Si è dato per assodato, dato che si simula per basse frequenze, che il solver stesso, di default, presenti passo di calcolo adeguato.

Di seguito il codice completo:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math
#Parametri
R1 = 1
RL = 3
L1 = 0.01
L2 = 0.1
M = 0.9* pow((L1L2),0.5) 
n = 1000 #Numero di campioni
V = 10 #Valore di picco tensione in ingresso
f = 10 #Frequenza
W = 2 * 3.14 * f
t_max = 0.6
#Inizzializza vettrori correnti
i_1 = np.zeros(n)
i_2 = np.zeros(n)
#Modello del circuito con accoppiamenti mutui
def trafo(Z,t,e):
   r1 = R1
   rl = RL
   l2 = L2
   l1 = L1
   m = M
   a = l1-(pow(m,2)/l2)
   b = mrl/l2 #Condizioni iniziali (al tempo [i])
   i1 = Z[0] 
   i2 = Z[1] 
   #Sistema di equazioni differenziali
   di1_dt = (b/a)i2-(r1/a)i1+e/a 
   di2_dt = (- (rl/l2)-mb/(l2a))i2+(mr1/(l2a))i1-(m/(l2a))*e
   #Vettore soluzioni del sistema
   Z =[di1_dt ,di2_dt]
   return Z
#condizioni iniziali al tempo 0
z0 = [0,0]
#time points
t = np.linspace(0,t_max,n)
#input
e = np.zeros(n)
for i in range(1,n):
   e[i] = 5 + Vmath.sin(Wt[i])
for i in range(1,n):
   #passo temporale
   tspan = [t[i-1],t[i]]
   #risolvi step
   z = odeint(trafo,z0,tspan,args=(e[i],))
   i_1[i] = z[0][0]
   i_2[i] = z[0][1]
   #Aggiorna le "condizioni iniziali"
    z0 = z[1]

Una nota sulla definizione dei parametri è necessaria. Dato un circuito la cui morfologia è quella qui presentata, è necessario che:

Questa condizione fisica deve essere rispettata per amor della Fisica: se così non fosse, l’accumulo energetico sarebbe negativo. Notare che non è una disuguaglianza in senso stretto, è ammessa anche l’uguaglianza, ma per una maggiore coerenza del modello matematico si è scritto:

M = 0.9* pow((L1L2),0.5)

Il coefficiente 0.9 è stato scelto di autorità, ma non deve mai superare l’unità (per le suddette ragioni).

Verifica bontà del modello

I parametri sono arbitrari, ma il comportamento del sistema “reale” è noto, pertanto si procede con simulazioni i cui risultati sono ben noti. Se il risultato della simulazione è coerente con la fisica del problema, allora il modello sarà valido.

Transitori impulsivi

Applicando un gradino di tensione ci si aspetta che:

a) La corrente iniziale sia nulla: le induttanze non ammettono gradini di corrente;

b) Il valore finale della corrente i1 (a regime) sia dato dal rapporto tra l’ampiezza del gradino e la resistenza R1;

c) Il valore di corrente al secondario sia nullo a transitorio estinto: non essendovi, a regime, variazioni di tensione sul primario, non può sussistere f.e.m. al secondario.

Quindi, dato il seguente gradino di tensione:

e = np.zeros(n)
for i in range(1,n):
    if i > n/10:
    e[i] =V

Il risultato è il seguente:

la tensione v1 (in viola) assume il valore di “e” all’istante iniziale, dal momento che la corrente i1 (in verde) è nulla. Negli istanti successivi la corrente i1 inizia a crescere fino a raggiungere asintoticamente il suo valore finale (che è pari a V/R1). Notare che al secondario i parametri tensione e corrente sussistono solamente durante le variazioni dei parametri al primario: la corrente i2 (in rosso) e la tensione v2 (in arancione) si evolvono secondo le dinamiche di un transitorio del primo ordine ma si estinguono una volta che sul primario i parametri diventano costani.

NOTA: la corrente i2 assume valori negativi. Questo è coerente per la convenzione di segno adottata. In caso contrario la resistenza produrrebbe energia rendendo necessario l’intervento delle autorità ecclesiastiche.

Disaccoppiando i circuiti, ponendo M=0, ci si aspetta un transitorio RL al primario (al secondario non accade nulla).

Infatti:

Chiusura di un generatore (ideale) di f.e.m. sinusoidale

Alimentando il tutto da una rete prevalente ci si aspetta la presenza di componente unidirezionale della corrente.

La tensione sarà così formata:

e = np.zeros(n)
for i in range(1,n):
   e[i] = 5 + V*math.sin(W*t[i])

Nota: la presenza di offset serve a testare ulteriormente il modello, dal momento che non la si deve ritrovare al secondario (per ovvie ragioni fisiche).

In questo caso ci si aspetta che:

a) la componente continua della tensione e(t) non abbia effetti sul secondario;

b) La componente unidirezionele di corrente deve estinguersi;

c) Le correnti siano in ritardo (carico induttivo):

Lanciando la simulazione, il risultato è il seguente:

A colpo d’occhio si può affermare che il risultato è coerente con la Fisica. Nello specifico:

La componente unidirezionale si annulla. Al secondario, a regime, non vi è traccia di offset.

Al primario è presente un offset di tensione e quindi di corrente: si osservi che, a regime, il valor medio della tensione sull’indutture (mutuamente accoppiato) è nullo. Se così in fosse, la corrente divergerebbe senza limite.

Si può affermare che il modello matematico risponde bene.

Matteo Gentileschi

Mi chiamo Matteo Gentileschi ed ebbi i natali nel lontano 1991, a Rieti. Appassionato di elettronica e programmazione, perito industriale diplomato all’ I.T.I.S. “C.Rosatelli” di Rieti, Ingegnere Elettrotecnico presso il Politecnico di Torino, coordinatore per la sicurezza in fase di progettazione (C.S.P), dichiaratamente “smanettone” e con un debole per l’ “open source“.

Ho deciso di aprire questo piccolo blog per raccogliere alcuni miei appunti e progetti, perchè ritengo che il confronto sia d’uopo, che sia il solo modo per progredire. Pertanto, si accettano critiche: solo se ben argomentate e costruttive.

Questo blog è e resterà sempre totalmente demonetizzato.

Personalmente faccio una fatica non idifferente a scrivere un articolo, ne ho sempre una decina “in bozza”. Ogni codice viene opportunatamente provato, ogni schema viene implementato e testato, ogni concetto viene “revisionato”. Sostanzialmente “garantisco” per quello che c’è scritto: non per questo sarà scevro da errori (è qui che serve il confronto).

Cerco di utilizzare prevalentemente piattaforme open source, perchè credo siano il futuro.

Matteo Gentileschi

Quanto segue è un ragionamento che ho trovato su un vecchio numero della rivista Tuttonormel (rivista del settore Impianti Elettrici). Ho sempre sostenuto che noi Elettrotecnici siamo dotati di un certo grado di goliardia, questo piccolo articolo ne fornisce prova schiacciante. Si vuole infatti dimostrare la relazione tra le natiche dei cavalli e le dimensioni dei vettori spaziali.

Nel progettare qualcosa, il contesto storico e culturale è sicuramente uno dei dati di ingresso più importanti. Disponibilità di materiali, grado di tecnologia ecc rimangono dunque vincoli importanti. Sicuramente, una delle leggi scientifiche più importanti (più importante delle leggi della Dinamica e della conservazione dell’Energia) è una delle leggi di Murphy:

“Se si necessita di N componenti per la realizzazione di un progetto , se ne avranno sempre a disposizione N-X, con X che va da 1 a N”

I razzi propulsori degli Shuttle (progettati dalla ATK Thikol nell’Utah) presentanto ingombri importanti. In realtà i progettisti avrebbero voluto farli più grandi ma hanno dovuto fare i conti con la fase di trasporto: sarebbero stati trasportati in treno.

Ovviamente, una ferrovia degna di tale nome passa attraverso qualche tunnel. Questo aspetto è il punto di partenza del ragionamento proposto.

• I tunnel ferroviari sono poco più larghi di una carrozza ferroviaria;
• La dimensione della carrozza ferroviaria è determinata dalla distanza tra i binari (scartamento);
• Le prime ferrovie americane furono costruite dagli inglesi che quindi le costruirono “a modo loro”(scartamento di 4 piedi e 8,5 pollici);
• Le prime ferrovie inglesi furono costruite dalle stesse imprese di costruzione che costruivano le tramvie;
• Le prime carrozze del tram avevano la stessa larghezza delle carrozze a cavalli;
• Le carrozze a cavalli avevano la stessa larghezza dei solchi presenti nelle vecchie strade inglesi (in caso contrario si sarebbero danneggiate più frequentemente);
• Le prime strade inglesi (ovviamente) sono state realizzate dai romani, che le percorrevano con i loro carri da guerra;
• I carri da guerra romani avevano la larghezza necessaria per “tenere” il sedere di due cavalli;

Da questi elementi si evince che:

“La larghezza dei razzi dello Shuttle è stata influenzata (e fortemente vincolata) da scelte/necessità di duemila anni prima”

Questo ci consente anche di affermare empiricamente che alcune norme tecniche sono fatte letteralmente “a culo di cavallo“.

Matteo Gentileschi

Braccio Meccanico

Un braccio meccanico è costituito da “links” (membri) e da giunti (servocomandi). Esiste una relazione analitica tra posizione angolare dei giunti e la “postura” del bracio.

Tramite semplici considerazioni geometriche è piuttosto semplice trovare la posizione (nello spazio) di un qualsiasi punto del braccio note le posizioni angolari dei singoli attuatori e le dimensioni dei singoli membri.

Meno semplice è la determinazione della posizione angolare dei giunti assegnato un obiettivo (es. posizione). Si tratta, in questo caso, della “cinematica inversa”.

LEGGE DEL COSENO

Uno dei metodi usati è quello basato sulla “Teorema del Coseno“, che mette in correlazione analitica la lunghezza dei membri e l’angolo tra due di essi:

Dall grafico precedente, partendo dalle coordinate x e y (obiettivo) si risale alle posizioni angolari come segue (codice in processing):

void calcolo_angoli(){
  //TRIGONOMETRIA 
  I = sqrt(pow(X,2)+pow(Y,2)); //lunghezza ipotenusa
  theta_t = atan(Y/X); //angolo tra ipotenusa e asse X
  //LEGGE DEL COSENO
  angolo[1] = theta_t +acos((pow(L1,2)+pow(I,2)-pow(L2,2))/(2*L1*I));
  // NOTA: la posizione angolare viene riferita all'asse X
  angolo[2] = acos((pow(L1,2)-pow(I,2)+pow(L2,2))/(2*L1*L2));
}

Appare evidente che non sempre esistono soluzioni: se ad esempio le coordinate fossero tali da avere “ipotenusa” più lunga della massima estenzione del braccio, non si avrebbero soluzioni (fisicamente non può arrivarci!).

Lo spazio W dei punti p(x,y) ammissibili (raggiungibili) è il seguente:

Notare, inoltre, che la solusione può non essere univoca.

A seguire un’animazione in Processing (vedasi come fare animazioni in Processing):

float theta_t;
float[] angolo = new float[2];

float L1 = 120,L2 = 120;

float X,Y;
float X0,Y0=0;
float I;

float t = 0,passo = 0.1;
float R = 40,r = 20;

void setup(){
  size(250,270,P3D);
  X0 = 1.2*(L1+L2)/2;
}

void draw(){
  background(155);
  translate(width/10,height/2);
  //equazione parametrica circonferenza
  X =X0+ R*cos(t);
  Y =Y0+ R*sin(t);
  t = t + passo; //incremento tempo 
  noFill();
  strokeWeight(2);
  ellipse(X0,Y0,2*R,2*R); 
  //ellipse(X0,Y0,2*R-r,2*R-r);
  //ellipse(X0,Y0,2*R+r,2*R+r);
  calcolo_angoli();//funzione calcolo angoli (legge del coseno)
  fill(0,255,0);
  strokeWeight(4);
  //costruzione "postura"
  rotateZ(-angolo[0]);
  line(0,0,L1,0);
  ellipse(0,0,r,r);
  translate(L1,0);
  rotateZ(PI-angolo[1]);
  line(0,0,L2,0);
  ellipse(0,0,r,r);
  fill(255,0,0);

  ellipse(L2,0,r,r);
}

void calcolo_angoli(){
  //TRIGONOMETRIA
  I = sqrt(pow(X,2)+pow(Y,2));
  theta_t = atan(Y/X);
  //LEGGE DEL COSENO
  angolo[0] = theta_t +acos((pow(L1,2)+pow(I,2)-pow(L2,2))/(2*L1*I));
  angolo[1] = acos((pow(L1,2)-pow(I,2)+pow(L2,2))/(2*L1*L2));
}

Le cose si complicano se il numero di membri è maggiore di due.

FABRIK

FABRIK è l’acronimo (Forward And Backward Reaching Inverse Kinematics”).

In figura viene mostrato il principio che sta alla base del FABRIK, per il cui utilizzo è necessario conoscere lo stato iniziale:

Si vuole spostare il punto V1 nel punto t (target): il tratto v1-v2 viene spostato in modo tale da cadere sulla direzione della congiungente tra t e v2. Ovviemente per garantire la continuità meccanica è necessario che anche il secondo tratto venga traslato e ruotato in modo da farlo giacere sulla direzione della congiungente v2-v3. Il procedimento si ripete per tutti i membri.

A seguire un esempio in Processing:

float L = 10; //lunghezze
int n = 50; //numero elementi del braccio 
PVector a,b,dir;//vettori di servizio 

float Xt,Yt;//coordinate "target" t

//MATRICE DELLE COORDINATE DELLE ARTICOLAZIONI 
float[][] coordinate = new float[2][n+1]; 

float t=0;
void setup(){
  size(500,500);
  
  //NECESSARIO IMPOSTARE UNO STATO INIZIALE:
  for (int i = 0; i <n; i = i+1) {
      coordinate[1][i] = 0; 
      coordinate[0][i] = (i+1)*L;
  }
  
}

void draw(){
  background(155);
  Xt = width/2 +width/3*cos(t);
  Yt = height/2 +width/3*sin(t);
  if(t>2*PI){
    t = 0;//periodico 2*PI
  }
  t = t+0.1;
  //Xt = mouseX;
  //Yt = mouseY;
  delay(100);
  coordinate[0][0] = Xt;
  coordinate[1][0] = Yt;
  
  for (int i = 0; i <n; i = i+1) {
  //per ogni segmento se ne prendono le coordinate degli estremi 
  a = new PVector(coordinate[0][i],coordinate[1][i]);
  b = new PVector(coordinate[0][i+1],coordinate[1][i+1]);
  dir = PVector.sub(a,b);//vettore nella direzione a-b
  dir.setMag(L);//settaggio lunghezza del suddetto vettore
  dir.mult(-1);//rotazione
  b = PVector.add(a,dir);//aggiunta ad "a"
  
  //aggiornamento matrice delle coordinate:
  coordinate[0][i+1] = b.x;
  coordinate[1][i+1] = b.y;
  }
  
  //tracciamento 
  for (int i = 0; i <n; i = i+1) {
  strokeWeight(5);
  line(coordinate[0][i] ,coordinate[1][i],coordinate[0][i+1] ,coordinate[1][i+1]);
  }
  
}

è possibile cambiare sia la lunghezza L dei singoli elementi che il numero totale n.

Matteo Gentileschi

Di base, il linguaggio è molto semplice, consente di concentrare l’attenzione esclusivamente sulla grafica e sugli algoritmi da implementare.

Si parte da una schermata in cui l’origine degli assi è il vertice in alto a sinistra (l’asse delle ordinate e rivolto verso il basso). Tale sistema di riferimento può essere traslato e ruotato:

translate(x,y,y); //sposta l'origine degli assi in (x,y,z);

rotateX(a);//ruota il pino perpendicolare all'asse X di "a" radianti
rotateY(b);//ruota il pino perpendicolare all'asse Y di "b" radianti
rotateZ(c);//ruota il pino perpendicolare all'asse Z di "c" radianti
//n.b. rotazione antioraria

Una volta modificato il riferimento è possibile disegnare gli oggetti tramite parametri e coordinate.

Se si è interessati alla creazione di gif animate è sufficiente “prelevare” dei fotoframe in sequenza. La seguente funzione crea uno “screenshot” della schermata e la salva nella cartella ove è salvato lo sketch in questione:

saveFrame("Nome_immagine"+contatore+".jpg");

A seguire una carrellata di alcune animazioni. Carrellata “casuale”, ovvero tale articolo verrà aggiornato costantemente ogni volta che mi capita di creare qualche piccola animazione.

Elisse

Equazione ellisse con centro Xc e Yc:

Parametrizzando le coordinate x e y si ha:

Si nota subito che se si volesse tracciare una circonferenza basta porre uguali l’asse maggiore e l’asse minore ( “a” e “b” in questo caso sono metà asse ), entrambi uguali al raggio.

Codice Processing:

float xc,yc,x,y;//centro ellisse
int a = 60,b=100;//parametri ellisse
float theta = 0; //parametro
float r = 20;//raggio pallina

void setup() {
  size(300,300);
  xc = width/2;//coordinata centro
  yc= height/2;//coordinata centro
}

void draw() {
  background(155);
  noFill();
  ellipse(xc, yc, 2*a,2*b);//disegna ellisse
  x = xc+a*cos(theta);
  y = yc+b*sin(theta);
  theta = theta+0.1;//incremento parametro 
  fill(0,0,100);
  circle(x,y,r);//disegna pallina
}

MOTO DELLA TERRA

Quello che segue è un esempio di ellisse in tre dimensioni (3D), che potrebbe rappresentare, ad esempio, un orbita. Infatti, ci si propone di “simulare” il moto relativo tra Terra e Luna.

nota: le proporzioni tra Terra e Luna non sono rispettate, come i rapporti tra le velocità angolari, per questioni di licenza artistica….

float x,y; //coordinate Luna

float theta = 0; //parametro
float rt = 60,rl =10;//raggio Terra e Luna
float dist = 100;//distanza Terra-Luna

PImage terra,luna,universo;
PShape mondo,moon;

void setup() {
  size(600,309,P3D);
  //creazione universo!!!!!!!!!
  universo = loadImage("universo.jpg");
  //creazione del Mondo!!!!!!!!!!!!!!!!!!!!!! XD
  noStroke();
  terra = loadImage("mondo.jpg");
  mondo = createShape(SPHERE, rt);//creazione shape "Mondo"
  mondo.setTexture(terra);
  //creazione Luna
  luna = loadImage("full_moon.jpg");
  moon = createShape(SPHERE, rl);//creazione shape "Luna"
  moon.setTexture(luna);
  
}

void draw() {
  background(universo);
  translate(width/2,height/2);//sposta l'origine degli assi al centro 

  rotateZ(PI/10);
  rotateY(theta);
  shape(mondo);//disegna la Terra
  //line(0,height,0,0,-height,0);//asse terrestre
  rotateY(-theta);
  theta = theta+0.03;//tempo 
  x = dist*cos(-4*theta);
  y = dist*sin(-4*theta);
  rotateX(PI/2);
  translate(x,y);
  rotateZ(theta);
  shape(moon);//disegna la Luna
  delay(70);
}

SPIRALE

Essendo il raggio della spirale parametrizzabile come segue:

Pertanto le coordinate saranno:

float a=1,b=1,theta = 10;
float x,y;
float rs = 40;

float step = 0.3;
PShape pallina;

void setup() {
  size(400,400,P3D);
  pallina = createShape(SPHERE,rs);//crea oggetto pallina
}

void draw() {
  background(155);
  //sposta al centro il sistema di riferimento:
  translate(width/2,height/2); 
  //Equazioni parametriche della spirale:
  x = (a + theta*b)*cos(theta);
  y = (a + theta*b)*sin(theta);
  theta = theta +step;//incremento parametro 
  //sposta il sistema di riferimento in (x,y):
  translate(x,y);
  shape(pallina);//disegna la pallina
  //sposta il sistema di riferimento "nella parte opposta":
  translate(-2*x,-2*y);
  shape(pallina);//disegna la pallina  
}

SEGA CIRCOLARE

Si possono creare degli oggetti e dei gruppi di oggetti per poi semplificare la scrittura del codice, richiamandoli come oggetti unici:

PShape oggetto_1, oggetto_2, complessivo; 

Void setup(){
  size(450,450,P3D);
  oggetto_1 = createShape(); //creazione
  oggetto_1.beginShape(TRIANGLES); //inizio creazione (TRIANGOLO)
  oggetto_1.vertex(a,b,c); //vertici (per coordinate)
  oggetto_1.vertex(e,f,g);
  oggetto_1.vertex(h,i,l);
  oggetto_1.endShape();// fine creazione 

  oggetto_2 = createShape(); //creazione
  oggetto_2.beginShape(TRIANGLES); //inizio creazione (TRIANGOLO)
  oggetto_2.vertex(a2,b2,c2); //vertici (per coordinate)
  oggetto_2.vertex(e2,f2,g2);
  oggetto_2.vertex(h2,i2,l2);
  oggetto_2.endShape();// fine creazione 

 complessivo ==createShape(GROUP); //creazione gruppo di oggetti 
 complessivo.addChild(oggetto_1); //aggiungi ogetto_1 al "complessivo"
 complessivo.addChild(oggetto_2); //aggiungi ogetto_2 al "complessivo"

}


void draw(){
shape(complessivo); //disegna l'oggetto "complessivo"
}

Ecco un esempio:

PShape dente_sega,a,b,c,d; 

float L = 70,K = 40;;
float H = 40;
float M = 60;

int count =0;
int num_denti = 20;
float R;//raggio del disco 
float step=0;
void setup(){
  size(450,450,P3D);
  
  //CREAZIONE "DENTE DI SEGA" 
  a = createShape(); 
  a.beginShape(TRIANGLES);
  a.vertex(-L/2,0,-H/2);
  a.vertex(-L/2,0,H/2);
  a.vertex(L/2,0,0);
  a.endShape();
  
  b = createShape(); 
  b.beginShape(TRIANGLES);
  b.vertex(-L/2,0,-H/2);
  b.vertex(-L/2,0,H/2);
  b.vertex(-L/2,-M,0);
  b.endShape();
  
  c = createShape(); 
  c.beginShape(TRIANGLES);
  c.vertex(-L/2,0,-H/2);
  c.vertex(-L/2,-M,0);
  c.vertex(L/2,0,0);
  c.endShape();
  
  d = createShape(); 
  d.beginShape(TRIANGLES);
  d.vertex(-L/2,-M,0);
  d.vertex(L/2,0,0);
  d.vertex(-L/2,0,H/2);
  d.endShape();
   
  dente_sega =createShape(GROUP);
  dente_sega.addChild(a);
  dente_sega.addChild(b);
  dente_sega.addChild(c);
  dente_sega.addChild(d);
    
  R = K/((2*PI)/num_denti);//calcolo del raggio in funzione del
  //numero di denti
}


void draw(){
  background(155);
  translate(width/2,height/2);

  rotateX(PI/5);
  rotateY(step);
  rotateZ(-2*step);
  
  if (step >2*PI){
    step=0;
   //periodicità 2*PI
  }
  else{
    step = step + 0.1;
  }
  
  //CICLO FOR PER DISEGNARE LA CORONA DENTATA
  for (int i =1 ; i < num_denti +1; i = i+1) {
    rotateZ((2*PI)/num_denti);
    shape(dente_sega,0,-R);
   }
  
  delay(100);
}

SONAR

Quando viene tracciata una linea è possibile settarne alcuni parametri grafici, come spessore, colore ecc.

Il colore viene imposto tramite il modello RGB:

stroke(Reed,Green,Blue,opacità);

I parametri Red, Green e Blue (che vanno da 0 e 255), rappresentano la “quantità” di quel colore nel colore finale che si vuole ottenere. Esistono apposite tabelle che mostrano i valori di tali parametri sulla base del colore che si vuole ottenere.

float R,x,y,t = 0;
float alpha=255;
//int cont =0;
void setup() {
  size (400, 400); 
  R = 0.8*height;
}

void draw() {
  background(0);
  translate(width/2,height/2);
  stroke(98,245,31,alpha);
  noFill();
  strokeWeight(3);
  ellipse(0,0,R,R);
  ellipse(0,0,0.8*R,0.8*R);
  ellipse(0,0,0.6*R,0.6*R);
  ellipse(0,0,0.4*R,0.4*R);
  ellipse(0,0,0.2*R,0.2*R);
  line(0,R, 0,-R);
  line(R,R,-R,-R);
  line(-R,0,R,0);
  line(R,-R,-R,R);
  t = t + 0.09;
  strokeWeight(1);
  for (int i = 0; i <alpha; i = i+1) {
  x = 0.55*R*cos(t-i*0.01);
  y = 0.55*R*sin(t-i*0.01);
  stroke(98,245,31,alpha-i);
  line(0,0, x,y);
}
  //saveFrame("sonar"+cont+".jpg");
  //cont = cont+1;
  delay(20);
 
}

Matteo Gentileschi

Un servocomando è un dispositiovo che amplifica la potenza di una funzione di comando mettendo in movimento un dispositivo meccanico. In sostanza sono dei trasduttori di posizione: si portano alla posizione richiesta.

Il servocomandi AC ricevono in ingresso un treno di impulsi (generalmente alla frequenza di 50Hz) con duty-cycle variabile all’interno di una determinata finestra.

Sono dispositivi molto versatili e piuttosto robusti. Nelle applicazioni più comuni non si ha un controllo ad anello chiuso, ovvero si da l’input ed il servocomando cercherà di raggiungere l’obiettivo. Per tale motivo, se rimanesse incastrato, continuerebbe a fornire coppia ma un motore che fornisce coppia a rotore fermo (coppia di spunto) equivale ad un cortocircuito (a transitorio estinto).

Per tale motivo è indispensabile consultare i dati forniti dalla casa produttrice. Nel caso in esame, se si volesse far muovere un peso di 1.6kg e si alimentasse con una tensione di 4.8v, il servo si danneggerà: in quelle condizioni infatti, la coppia di stallo verrebbe vinta dal carico che trascinerà l’albero.

La velocità nominale (operating speed) viene definita infatti a vuoto (senza carico): a 4.8v si riuscirà ad avere rotazioni di 60° in 0.12 secondi.

La velocità di uscita del servo non è la medesima dell’albero del motore interno ma modificata tramite “motoriduttori”: Riducendo il numero di giri a parità di potenza si ha coppia maggiore. Quello che si vuole ottenere da tali dispositivi è infatti una coppia elevata.

Calettato all’albero d’uscita vi sarà un potenziometro rotativo: l’uscita del partitore resistivo sarà pertanto proporzionale alla posizione fisica.

La differenza tra la posizione richiesta e quella fisica, chiamata “errore” di posizione, viene amplificata tramite un regolatore che ne aumenta la stabilità.

Pilotaggio di un servocomando con Arduino

Un semplice sketch in Processing può fungere da pannello di controllo per un servocomando:

import processing.serial.*; //libreria per la comunicazione seriale
Serial porta_arduino; //variabile inizializzata 

float x_sx,x_dx,y = 100,r = 50; //pulsanti

float angle = 90,L = 105,Y0 = 200;

float lim_inf = 20,lim_sup = 160;//limiti 

float arr = 10;

int delay = 10 ; //ritardo in ms
void setup(){
     size(400,250);
     x_sx = width/4;
     x_dx = width/4+width/2;
     porta_arduino = new Serial(this, "COM5", 9600);//settaggio parametri 
}
  
  
  void draw(){
      background(100);
      strokeWeight(3);  // Default
      text("UP",x_sx,150);
      text("DOWN",x_dx,150);
      text("POSIZIONE SERVOCOMANDO",width/2,30);
      circle(x_sx,y,r);
      circle(x_dx,y,r);
      text(angle,width/2,Y0+3*arr);
      textAlign(CENTER);

      translate(width/2,Y0);
      rotate(-radians(angle));
      line(0,0,L,0);
      triangle(L, 0, 0, arr, 0, -arr);
      circle(0,0,2*arr);
      rotate(radians(angle));

      rotate(-radians(lim_sup));
      line(0,0,L,0);
      rotate(radians(lim_sup-lim_inf));
      line(0,0,L,0);

      if(mousePressed){
      //se il pulsante del mouse viene premuto:
      //controlla se è all'interno del pulsante
      //up o down
          if(sqrt(pow(mouseX-x_sx,2)+pow(mouseY-y,2))<r){
            if(angle < lim_sup){angle = angle +1;}//controllo limiti 
          }
        
          if(sqrt(pow(mouseX-x_dx,2)+pow(mouseY-y,2))<r){
             if(angle > lim_inf){angle = angle -1;}//controllo limiti 
          }
          porta_arduino.write(int(angle));//invia dato ad Arduino
          delay(delay); //ritardo 
      }

}

è opportuno impostare dei limiti (in termini di posizione angolare) oltre i quali l’attuatore arriverebbe a battuta.

Per quanto riguarda la scheda Arduino invece, è possibile utilizzare la libreria Servo.h.

Una volta definito l’oggetto della classe “servo” si determina, nel setup, il pin a cui è collegato:

servo_1.attach(pin_servo); //pin a cui è collegato il servo 

Per impostare la posizione, è sufficiente chiamare la seguente funzione:

    servo_1.write(angle);//scrivi posizione angolare

Ove l’argomento (angle) rappresenta l’angolo desiderato, da 0 a 180°.

Notare che una scheda Arduino non può gestire un numero qualsiasi di servocomandi tramite tale libreria, pertanto è necessario consultare la documentazione relativa alla scheda in esame.

A seguire uno sketch che consente ad Arduino di portare il servocomando nella posizione inviatagli tramite porta seriale.

#include <Servo.h>
#define pin_servo 7
int angle;//angolo in gradi 

Servo servo_1;//oggetto della classe "Servo"

void setup() {
  servo_1.attach(pin_servo); //pin a cui è collegato il servo 
  Serial.begin(9600); 
  servo_1.write(90);  // set servo to mid-point

}

void loop() {
  if (Serial.available() >0) {
    angle = Serial.read();
    servo_1.write(angle);//scrivi posizione angolare
    //N.B. L'argomento (angolo) è in gradi
  }
  
}

Matteo Gentileschi