Skip to content
Snippets Groups Projects
Verified Commit fbb3c76d authored by Christian Rohlfing's avatar Christian Rohlfing
Browse files

Added video for RLC system

parent 6a8faef3
No related branches found
No related tags found
1 merge request!6New DFT demo and minor changes in Laplace/z-Transform demo
%% Cell type:code id: tags:
``` python
# Copyright 2020 Institut für Nachrichtentechnik, RWTH Aachen University
%matplotlib widget
from ipywidgets import interact, interactive
import ipywidgets as widgets
from IPython.display import clear_output, display, HTML
import matplotlib.pyplot as plt
import numpy as np
import cmath # for sqrt(-1)
import rwth_nb.plots.mpl_decorations as rwth_plots
from rwth_nb.misc.signals import *
plt.close('all')
```
%% Cell type:markdown id: tags:
<div class="inline-block">
<img src="figures/rwth_ient_logo@2x.png" style="float: right;height: 5em;">
</div>
# RLC-System
Zum Starten: Im Menü: Run <span class="fa-chevron-right fa"></span> Run All Cells auswählen.
%% Cell type:markdown id: tags:
## Einleitung
Es wird folgendes RLC-System betrachtet:
![Blockdiagramm](figures/rlc_system_block_diagram.png)
%% Cell type:code id: tags:
``` python
# Exemplary values
R = 16 # Ohm
R = 16 # Ohm, Ω
L = 1.5E-3 # Henry, mH
C = 1E-6 # Farad, myF
```
%% Cell type:markdown id: tags:
## Berechnung Laplace-Übertragungsfunktion
Die Übertragungsfunktion des Systems kann im Laplace-Bereich mittels Spannungsteiler berechnet werden:
$$
H(p)
=
\frac{U_2(p)}{U_1(p)}
=
\frac{1/(Cp)}{R+Lp+1/(Cp)}
=
\frac{1}{LCp^2+RCp+1}
=
\frac{1/(LC)}{p^2+(R/L)p + 1/(LC)}
$$
mit $p = \sigma + \mathrm{j}\omega = \sigma + \mathrm{j}2 \pi f$.
%% Cell type:code id: tags:
``` python
# Laplace transfer function
H_laplace = lambda p: 1/(L*C*p**2+R*C*p+1);
```
%% Cell type:markdown id: tags:
### Polstellen
Die zugehörigen Polstellen berechnen sich zu
$$p_{\mathrm{P},1,2} =
- \underbrace{\frac{R}{2L}}_{=a} \pm \underbrace{\sqrt{\frac{R^2}{4L^2} - \frac{1}{LC}}}_{=b}
$$
mit $a=\frac{R}{2L}$ und $b=\sqrt{\frac{R^2}{4L^2}-\frac{1}{LC}}$.
%% Cell type:code id: tags:
``` python
# Poles
a = R/(2*L)
b = cmath.sqrt(R**2/(4*L**2)-1/(L*C))
p_p1 = -a+b
p_p2 = -a-b
# Print out the numbers
print("a={0:.3f}".format(a))
print("b={0:.3f}".format(b))
if R**2/(4*L**2) >= 1/(L*C): print("b reell")
else: print("b imaginär")
print("\nPolstellen:")
print("p_p1={0:.3f}, p_p2={0:.3f}".format(p_p1, p_p2))
```
%% Cell type:markdown id: tags:
### Pol-Nulstellendiagramm
Nachfolgend wird das Pol-Nulstellendiagramm geplottet. Es enthält die beiden konjugiert komplexen Polstellen, den Konvergenzbereich und das zugehörige $H_0$.
Da der Konvergenzbereich die imaginäre Achse beinhaltet, ist das System stabil.
%% Cell type:code id: tags:
``` python
beta = np.imag(b) # Imaginary part of the poles
pp = np.array([p_p1, p_p2]); pz = np.array([]) # Zeros # Poles and Zeros
ord_p = np.array([1, 1]); ord_z = np.array([]) # Poles' and Zeros' orders
roc = np.array([np.max(np.real(pp)), np.inf]) # region of convergence
# Plot
fig, ax = plt.subplots()
ax.plot(np.real(pp), np.imag(pp), **rwth_plots.style_poles); ax.plot(np.real(pp), -np.imag(pp), **rwth_plots.style_poles); rwth_plots.annotate_order(ax, pp, ord_p);
ax.plot(np.real(pz), np.imag(pz), **rwth_plots.style_zeros); ax.plot(np.real(pz), -np.imag(pz), **rwth_plots.style_zeros); rwth_plots.annotate_order(ax, pz, ord_z);
rwth_plots.plot_lroc(ax, roc, 500, np.imag(p_p1)); ax.text(-1000,ax.get_ylim()[1]*0.8,'$H_0 = 1/(LC)$',bbox = rwth_plots.wbbox);
ax.set_xlim(ax.get_xlim()[0], 300); rwth_plots.annotate_xtick(ax, r'$-a$', -a,0,'k');
ax.set_xlabel(r'$\rightarrow \mathrm{Re}\{p\}$'); ax.set_ylabel(r'$\uparrow \mathrm{Im}\{p\}$'); rwth_plots.grid(ax); rwth_plots.axis(ax);
```
%% Cell type:markdown id: tags:
## Fourier-Übertragungsfunktion
Aus der Laplaceübertragungsfunktion kann die Fourierübertragungsfunktion berechnet werden, indem $p = \mathrm{j}2\pi f$ gesetzt wird:
$$H(p = \mathrm{j}2\pi f) \quad\text{mit}\quad \omega_0 = \frac{1}{\sqrt{LC}}$$
%% Cell type:code id: tags:
``` python
# Fourier transfer function
H_fourier = lambda f: H_laplace(1j*2*np.pi*f)
f = np.linspace(0, 10000, 10000) # frequency axis
# Resonance frequency
omega0 = 1/np.sqrt(L*C)
f0 = omega0/(2*np.pi)
# Print f0
print("f0 = {0:.2f} Hz".format (f0))
# Plot
fig,ax = plt.subplots()
ax.plot(f/1000, np.abs(H_fourier(f))); rwth_plots.axis(ax)
ax.set_xlabel(r'$\rightarrow f$ [kHz]'); ax.set_ylabel(r'$\uparrow |H(f)|$');
rwth_plots.annotate_xtick(ax, r'$f_0 = \omega_0/(2\pi)$', omega0/(2*np.pi)/1000,-0.25,'rwth:black');
ax.axvline(f0/1000,0,1, color='k',linestyle='--');
```
%% Cell type:markdown id: tags:
## Impulsantwort
Die Impulsantwort kann mittels Partialbruchzerlegung (siehe Vorlesung) bestimmt werden zu
$$
h(t) = \frac{\mathrm{e}^{-at}}{bLC}
\frac{\mathrm{e}^{bt}-\mathrm{e}^{-bt}}{2}
\cdot \varepsilon(t)
$$
mit $a=\frac{R}{2L}$ und $b=\sqrt{\frac{R^2}{4L^2}-\frac{1}{LC}}$.
Der nachfolgende Plot zeigt diese Impulsantwort.
%% Cell type:code id: tags:
``` python
# Impulse response (time-domain)
h = lambda t: np.real(np.exp(-a*t)/(b*L*C)*(np.exp(b*t)-np.exp(-b*t))/2*unitstep(t))
t = np.linspace(-0.001, 0.01, 10000) # time axis
# Plot
fig,ax = plt.subplots()
ax.plot(t*1000, h(t), 'rwth:blue'); rwth_plots.axis(ax); ax.set_xlim([-0.1, 5])
ax.set_xlabel(r'$\rightarrow t$ [ms]'); ax.set_ylabel(r'$\uparrow h(t)$');
#np.mean(np.abs(h(t) - 1/(beta*L*C)*np.exp(-a*t)*np.sin(beta*t)*unitstep(t))**2)
```
%% Cell type:markdown id: tags:
## Interaktive Demo
In dieser interaktiven Demo kann das Verhalten des Systems für variable Werte von $R$ betrachtet werden. Über den Schieberegler kann der Wert für $R$ geändert werden, entsprechend sieht man die Änderungen für die Fourier-Übertragungsfunktion und die Impulsantwort.
%% Cell type:code id: tags:
``` python
fig0,axs = plt.subplots(2, 1, **rwth_plots.landscape)
@widgets.interact(Rsel=widgets.FloatSlider(min=0, max=200, step=1, value=R, description='$R$ [$\Omega$]'))
def update_plot(Rsel):
H_laplace = lambda p: 1/(L*C*p**2+Rsel*C*p+1);
H_fourier = lambda f: H_laplace(1j*2*np.pi*f);
a = Rsel/(2*L)
b = cmath.sqrt(Rsel**2/(4*L**2)-1/(L*C))
h = lambda t: np.real(np.exp(-a*t)/(b*L*C)*(np.exp(b*t)-np.exp(-b*t))/2*unitstep(t))
if not axs[0].lines: # Call plot() and decorate axes. Usually, these functions take some processing time
ax = axs[0]; ax.plot(f/1000, np.abs(H_fourier(f)), 'rwth:blue'); rwth_plots.axis(ax)
ax.set_xlabel(r'$\rightarrow f$ [kHz]'); ax.set_ylabel(r'$\uparrow |H(f)|$');
ax.axvline(f0/1000,0,1, color='rwth:black',linestyle='--');
ax = axs[1]; ax.plot(t*1000, h(t), 'rwth:blue'); rwth_plots.axis(ax); ax.set_xlim([-.225, 5])
ax.set_xlabel(r'$\rightarrow t$ [ms]'); ax.set_ylabel(r'$\uparrow h(t)$');
else: # If lines exist, replace only xdata and ydata since plt.plot() takes longer time
axs[0].lines[0].set_ydata(np.abs(H_fourier(f)))
axs[1].lines[0].set_ydata(h(t))
tmp = np.max(np.abs(H_fourier(f))); rwth_plots.update_ylim(axs[0], np.abs(H_fourier(f)), 0.1*tmp, tmp+10 )
tmp = np.max(np.abs(h(t))); rwth_plots.update_ylim(axs[1], (h(t)), 0.1*tmp, tmp+10 )
display(HTML('{}<br />{}'.format('a={0:.3f}'.format(a), 'b={0:.3f}'.format(b))))
```
%% Cell type:markdown id: tags:
### Experiment
Mit $R=0 \,\Omega$, $L=47\,\mathrm{mH}$ und $C=0,68\,\mu\mathrm{F}$ liegt die entsprechende Grenzfrequenz bei $f_0 \approx 900\,Hz$. Rauchentwicklung ab Sekunde 20 ;)
%% Cell type:code id: tags:
``` python
from IPython.display import Video
Video("figures/RLC_Experiment.mp4", width=480, height=270 )
```
%% Cell type:markdown id: tags:
___
This notebook is provided as [Open Educational Resource](https://en.wikipedia.org/wiki/Open_educational_resources) (OER). Feel free to use the notebook for your own purposes. The code is licensed under the [MIT license](https://opensource.org/licenses/MIT).
Please attribute the work as follows:
*Christian Rohlfing, Übungsbeispiele zur Vorlesung "Grundgebiete der Elektrotechnik 3 - Signale und Systeme"*, gehalten von Jens-Rainer Ohm, 2020, Institut für Nachrichtentechnik, RWTH Aachen University.
*Christian Rohlfing, Übungsbeispiele zur Vorlesung "Grundgebiete der Elektrotechnik 3 - Signale und Systeme"*, gehalten von Jens-Rainer Ohm, 2021, Institut für Nachrichtentechnik, RWTH Aachen University.
......
File added
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment