ac-electricity

acelectricity module

The aim of this project is to understand the rudiments of linear AC electrical circuits and their functioning.

Prerequisites

Python : version 3
All operating systems
Libraries :

• numpy
• matplotlib.pyplot
• matplotlib.widgets

Installing acelectricity python module

From Pypi repository :
https://pypi.org/project/ac-electricity

pip install ac-electricity

Basic electrical quantities

• Voltage (volt)
• Current (ampere)
• Impedance Z (ohm)
• Admittance Y (siemens)
• Active power P (watt)
• Reactive power Q (var)
• Apparent power S (VA)

Electrical laws

• Kirchhoff’s current law
• Kirchhoff’s voltage law
• Ohm's law
• Impedances, admittances in series and parallel
• Voltage divider and current divider
• Millman's theorem
• Joule's first law (S=I²Z)
• Power law (S=V.I*)

Others features

• Bode plot (with cursors and sliders)
• User-defined transfer function H(jω)
• Digital filter frequency response

AC Circuit diagram

Circuit should only contain :

• independent sine voltage source
• independent sine current source
• resistors, inductors, capacitors

Example circuit

                  VL
           <---------------------
                          L1
            +------+    _  _  _
      --->--| R1   |---/ \/ \/ \-----+-------+
        IL  +------+                 |       |
                                     |       |
    ^                             IR v       v IC
    |                                |       |      ^
    |                              +---+     |  C1  |
    |                              |   |   -----    |
Vin |                              |R2 |   -----    | Vout
    |                              |   |     |      |
    |                              +---+     |      |
    |                                |       | 
                                     |       |
      -------------------------------+-------+

Datas :
Vin = 5 Vrms
Inductor : L1 = 100 mH ; R1 = 180 Ω
R2 = 2.2 kΩ ; C1 = 330 nF

What are the IL, IR, IC currents ?
What are the VL, Vout voltages ?
What is the frequency response Vout/Vin ?

>>> from acelectricity import *
>>> Vin = Voltage(5)  # Vrms
>>> Zr1 = Impedance(r=180)
>>> Zl1 = Impedance(l=0.1)
>>> Zr2 = Impedance(r=2200)
>>> Zc1 = Impedance(c=330e-9)

>>> Zeq1 = 1/(1/Zr2 + 1/Zc1)  # impedances in parallel
>>> # or Zeq1 = Zr2//Zc1
>>> Zeq1
Complex impedance (Ω) : 100.88-460.173j @ 1.0 kHz
>>> Zeq = Zr1 + Zl1 + Zeq1  # impedances in series
>>> Zeq(2000)  # @ 2000 Hz
(206.11818300356995+1018.3560443060462j)
>>> IL = Vin/Zeq  # Ohm's law
>>> IL.properties(2000)
Frequency (Hz) : 2000
Angular frequency (rad/s) : 12566.4

Complex current : 0.000954663-0.00471665j
Amplitude (Arms) : 0.00481229
Amplitude (A) : 0.00680561
Amplitude (dBA ref 1 Arms) : -46.3529621108357
Phase (degrees) : -78.5577
Phase (radians) : -1.37109
i(t) = 0.00680561×sin(12566.4×t -1.371091)

>>> Vout = IL*Zeq1  # Ohm's law
>>> VL = Vin-Vout  # Kirchhoff’s voltage law
>>> IC = Vout/Zc1  # Ohm's law
>>> IR = IL-IC  # Kirchhoff’s current law
>>> H = Vout/Vin  # transfer function

>>> # draw Bode plot and save datas
>>> H.bode(title='Vout/Vin transfer function', filename='h.csv')
>>> IL.bode(magnitude_unit='default', yscale='linear', title='IL current')
>>> Zeq.bode(yscale='log', title='Zeq frequency response')
>>> Zeq.bode_imag(title='Zeq reactance frequency response')
>>> show()

Zoom and data cursors :

Impedances and admittances

>>> Yr2 = 1/Zr2
>>> Yc1 = 1/Zc1
>>> Yc1
Complex admittance (S) : 0+0.00207345j @ 1.0 kHz
>>> 1/Yc1
Complex impedance (Ω) : 0-482.288j @ 1.0 kHz
>>> 1/(Yr2+Yc1)
Complex impedance (Ω) : 100.88-460.173j @ 1.0 kHz
>>> Zr2*(Zc1/(Zr2+Zc1))
Complex impedance (Ω) : 100.88-460.173j @ 1.0 kHz

Law class

>>> law = Law()
>>> # voltage divider
>>> Vout = law.VoltageDivider(vtotal=Vin, z=Zr2//Zc1, z2=Zr1+Zl1)
>>> Vout
Complex voltage (Vrms) : -2.28808-6.8219j @ 1.0 kHz
>>> Vin*(Zeq1/Zeq)
Complex voltage (Vrms) : -2.28808-6.8219j @ 1.0 kHz
>>> # Millman's theorem
>>> gnd = Voltage(0)
>>> Vout = law.Millman(v_z=[(Vin, Zr1+Zl1), (gnd, Zr2), (gnd, Zc1)])
>>> Vout
Complex voltage (Vrms) : -2.28808-6.8219j @ 1.0 kHz
>>> (Vin/(Zr1+Zl1))/(1/(Zr1+Zl1) +1/Zr2 +1/Zc1)
Complex voltage (Vrms) : -2.28808-6.8219j @ 1.0 kHz
>>> Vout/Vin
Ratio : -0.457615-1.36438j @ 1.0 kHz
>>> # current divider
>>> IC = law.CurrentDivider(itotal=IL, z=Zc1, z2=Zr2)
>>> IC
Complex current (Arms) : 0.0141449-0.00474421j @ 1.0 kHz
>>> IL*Zr2/(Zr2 +Zc1)
Complex current (Arms) : 0.0141449-0.00474421j @ 1.0 kHz

Electrical power, Joule's first law

>>> S = IL*Vin  # input source complex power
>>> S.properties(1000)
Frequency (Hz) : 1000
Angular frequency (rad/s) : 6283.19

Complex power : 0.0655242+0.0392254j

Active power P (W) : +0.0655242
Reactive power Q (var) : +0.0392254
Apparent power S (VA) : 0.0763678

Phase (degrees) : +30.9064
Phase (radians) : +0.539419
Power factor PF : 0.8580073402000852

Active power (dBW ref 1 W) : -11.835985322667693

>>> Sr1 = law.Joule(z=Zr1, i=IL)
>>> Sr1
Complex power (W): 0.0419907+0j @ 1.0 kHz
>>> Zr1*IL*IL
Complex power (W): 0.0419907+0j @ 1.0 kHz
>>> Sr2 = law.Joule(z=Zr2, v=Vout)
>>> Sr2
Complex power (W): 0.0235334 @ 1.0 kHz
>>> Sl1 = Zl1*IL*IL
Complex power (W): 0+0.146575j @ 1.0 kHz
>>> Sc1 = Zc1*IC*IC
Complex power (W): 0-0.10735j @ 1.0 kHz
>>> Sr1 +Sr2 +Sc1 +Sl1
Complex power (W): 0.0655242+0.0392254j @ 1.0 kHz

Parameters analysis

>>> Vin.RMS = 10
>>> Vin.phase = 45
>>> Zr1.r = 100
>>> Zl1.l = 0.22
>>> Zr2.r = 1000
>>> Zc1.c = 100e-9
>>> H.bode(title='Vout/Vin transfer function', filename='h2.csv')
>>> IL.bode(magnitude_unit='default', yscale='linear', title='IL current')
>>> show()

>>> for Zr2.r in [10, 100, 1e3, 1e4]:
        H.bode(title="Vout/Vin with R2={} Ω".format(Zr2.r))
>>> show()

Add matplotlib widgets : slider and button

>>> fig, ax1, ax2, mag, ph = H.bode(title='Vout/Vin transfer function')
>>> fig.subplots_adjust(left=0.15, bottom=0.35, right=0.85, top=0.9)

>>> ax_r1 = plt.axes([0.15, 0.20, 0.4, 0.03])
>>> ax_l1 = plt.axes([0.15, 0.15, 0.4, 0.03])
>>> ax_r2 = plt.axes([0.15, 0.10, 0.4, 0.03])
>>> ax_c1 = plt.axes([0.15, 0.05, 0.4, 0.03])

>>> slider_r1 = Slider(ax_r1, 'r1 (Ω)', 10, 300,
    valinit=Zr1.r, valstep=1)
>>> slider_l1 = Slider(ax_l1, 'l1 (H)', 0.01, 1,
    valinit=Zl1.l, valstep=0.01)
>>> slider_r2 = Slider(ax_r2, 'r2 (Ω)', 1e3, 1e4,
    valinit=Zr2.r, valstep=10)
>>> slider_c1 = Slider(ax_c1, 'c1 (nF)', 100, 1000,
    valinit=Zc1.c*1e9, valstep=1)

>>> def update(val):
        Zr1.r = slider_r1.val
        Zl1.l = slider_l1.val
        Zr2.r = slider_r2.val
        Zc1.c = slider_c1.val*1e-9
        xvalues = mag.get_xdata()
        # H.db() method, according to H.bode() parameters :
        # magnitude_unit='db' (default)
        # help(H) for more information
        magnitudes = [H.db(x) for x in xvalues]
        # H.phase_deg() method, according to H.bode() parameters :
        # phase_unit='degrees' (default)
        phases = [H.phase_deg(x) for x in xvalues]
        mag.set_ydata(magnitudes)
        ph.set_ydata(phases)

>>> slider_r1.on_changed(update)
>>> slider_l1.on_changed(update)
>>> slider_r2.on_changed(update)
>>> slider_c1.on_changed(update)

>>> ax_button = plt.axes([0.7, 0.10, 0.15, 0.06])
>>> button = Button(ax_button, 'Autoscale')

>>> def autoscale(event):
        ax1.relim()
        ax1.autoscale_view()
        ax2.relim()
        ax2.autoscale_view()

>>> button.on_clicked(autoscale)
>>> show()

User-defined transfer function

Example 1 : second order band-pass filter

>>> from acelectricity import *
>>> # static gain, damping value, normal angular frequency
>>> a, z, wn = 10, 0.1, 1000*2*math.pi
>>> H = Ratio(fw=lambda w: a*(2*z*1j*w/wn)/(1+2*z*1j*w/wn-(w/wn)**2))
>>> H.bode(filename='H.csv')
>>> a, z, wn = 100, 0.5, 10000
>>> H.bode(filename='H2.csv')
>>> show()

or :

>>> a, z, wn = 10, 0.1, 1000*2*math.pi
>>> H = Ratio.transfer_function(numerator=[0, a*(2*z*1j/wn)],
    denominator=[1, 2*z*1j/wn, -1/wn**2])
>>> H.bode(filename='H3.csv')
>>> a, z, wn = 100, 0.5, 10000
>>> # new instance
>>> H = Ratio.transfer_function(numerator=[0, a*(2*z*1j/wn)],
    denominator=[1, 2*z*1j/wn, -1/wn**2])
>>> H.bode(filename='H4.csv')
>>> show()

Example 2 : cascaded series, parallel filters

>>> from acelectricity import *
>>> # first order low-pass filter
>>> wn = 10000
>>> Hlp = Ratio(fw=lambda w: 1/(1+1j*w/wn))
>>> # first order high-pass filter
>>> Hhp = Ratio(fw=lambda w: 1/(1+1j*wn/w))
>>> Hs = Hlp*Hhp  # cascaded series filters
>>> Hs.bode()
>>> Hp = Hlp+Hhp  # parallel filters
>>> Hp.bode()
>>> show()

Example 3 : linear control system

          +      +------+
     -->---(X)---| G(w) |----+--->--
          - |    +------+    |
            |                |
            |    +------+    |
            +----| H(w) |-<--+
                 +------+

>>> from acelectricity import *
>>> # feedforward transfer function
>>> # first order low-pass filter
>>> wn = 10000
>>> G = Ratio.transfer_function([1], [1, 1j/wn])
>>> # feedback transfer function
>>> H = Ratio.transfer_function([10])  # constant
>>> # open-loop transfer function
>>> Hopenloop = G*H
>>> Hopenloop.bode()
>>> # closed-loop transfer function
>>> Hcloseloop = G/(1+Hopenloop)
>>> Hcloseloop.bode()
>>> show()

Digital filter frequency response

y(n) = 0.1x(n) +1.6y(n-1) -0.7y(n-2)

>>> from acelectricity import *
>>> fs = 100000  # sampling rate (Hz)
>>> H = Ratio.digital_filter(fs=fs, b=[0.1], a=[1, -1.6, 0.7])
>>> H.bode(xmin=0, xmax=fs/2, xscale='linear',
    title='IIR digital filter')
>>> show()

y(n) = (x(n)+x(n-1)+...+x(n-7))/8

>>> from acelectricity import *
>>> fs = 10000  # sampling rate (Hz)
>>> H = Ratio.digital_filter(fs=fs, b=[1/8]*8)
>>> H.bode(xmin=0, xmax=fs/2, xscale='linear',
    magnitude_unit='default', title='FIR digital filter')
>>> show()

Custom default frequency

>>> from acelectricity import *
>>> Yc = Admittance(c=220e-6)
>>> Yc
Complex admittance (S) : 0+1.3823j @ 1.0 kHz
>>> ElectricalQuantity.DEFAULT_FREQUENCY = 50
>>> Yc
Complex admittance (S) : 0+0.069115j @ 50 Hz
>>> 1/Yc
Complex impedance (Ω) : 0-14.4686j @ 50 Hz

Goodies

Ideal filters

>>> from acelectricity import *
>>> def lp(w):
        wn = 10000
        return 1 if w < wn else 0.001
>>> def hp(w):
        wn = 1000
        return 10 if w > wn else 0.001
>>> def bp(w):
        wn1, wn2 = 1000, 10000
        return 1 if wn2 > w > wn1 else 0.001
>>> def notch(w):
        wn1, wn2 = 4000, 5000
        return 0.001 if wn2 > w > wn1 else 1
>>> def allpass(w):
        wn = 10000
        return 1j if w > wn else -1j
>>> Hlp = Ratio(fw=lp)
>>> Hhp = Ratio(fw=hp)
>>> Hbp = Ratio(fw=bp)
>>> Hnotch = Ratio(fw=notch)
>>> Hallpass = Ratio(fw=allpass)
>>> Hallpass.bode(title='Ideal all-pass filter')
>>> (Hhp+Hlp).bode(title='Ideal filter Bode plot')
>>> show()

Note : remember that ideal filters are not realizable (not causal).

Asymptotic Bode diagram

>>> def lp1(w):
        # first order low-pass filter 1/(1+1j*w/wn)
        # asymptotic approximation
        wn = 1000
        return 1 if w < wn else 1/(1j*w/wn)
>>> def lp2(w):
        wn = 30000
        return 1 if w < wn else 1/(1j*w/wn)
>>> Hlp1 = Ratio(fw=lp1);Hlp2 = Ratio(fw=lp2)
>>> (Hlp1*Hlp2).bode(title='Asymptotic Bode plot')
>>> show()

Advantages and limitations

This module manages basic arithmetic operations + - * / as well as // which designates two impedances in parallel.

The dimensional homogeneity is checked :

>>> V3 = V1 - V2 + I3  # V+A -> Error
TypeError : Voltage expected
>>> I1 = Current(2)
>>> I = I1 + 0.5  # A+number -> Error
TypeError : Current expected
>>> I2 = Current(0.5, phase=30)
>>> I = I1 + I2
>>> I
Complex current (Arms) : 2.43301+0.25j @ 1.0 kHz
>>> I = 5*I2 - V1/Z1 + I3

The result of any operation must give a quantity whose unit is one of : V, A, Ω, S, W (or ratio).
Otherwise, you will get an error :

>>> Z1/V1  # Ω/V -> 1/A -> Error
TypeError
>>> U2*(Z3/(Z2+Z3))  # V*(Ω/Ω) -> V*() -> V
>>> U2*Z3/(Z2+Z3)  # V*Ω -> Error
TypeError
>>> S = V1*(V1/Z1)  # V*(V/Ω) -> V*A -> W
>>> S = V1*V1/Z1    # V*V -> Error
TypeError

See also

https://pypi.org/project/dc-electricity