2. TCLab in the frequency domain
[1]:
import numpy
[2]:
%matplotlib inline
[3]:
tau_p = 150
K_p = 0.33
theta = 15
[4]:
omega = numpy.logspace(-3, -2)
s = omega*1j
[5]:
G = K_p/(tau_p*s + 1)*numpy.exp(-theta*s)
[6]:
import matplotlib.pyplot as plt
Let’s choose 5 logarithmically spaced points around the corner frequency
[7]:
freqs = numpy.logspace(-2.8, -2, 5)
[8]:
def plotfreqs(ax):
for freq in freqs:
ax.axvline(freq)
[9]:
def bode(omega, G, gainax=None, phaseax=None, phasecorr=0):
if gainax is None:
fig, (gainax, phaseax) = plt.subplots(2, 1, sharex=True)
gainax.loglog(omega, numpy.abs(G))
angle = numpy.angle(G)
phaseax.semilogx(omega, numpy.unwrap(angle) + phasecorr)
return gainax, phaseax
[10]:
gainax, phaseax = bode(omega, G)
plotfreqs(gainax)
2.1. Direct frequency domain tests
[11]:
Qbar = 50
deltaQ = A = 10
How long does one sine wave take to repeat?
\[P = 2\pi/\omega\]
[12]:
P = 2*numpy.pi/freqs
We’ll do a couple of repeats.
[13]:
nperiods = 3
[14]:
switch_times = numpy.concatenate([[0], numpy.cumsum(P*nperiods)])
[15]:
switch_times
[15]:
array([ 0. , 11893.26574888, 19397.409123 , 24132.20349893,
27119.65678502, 29004.61237718])
[16]:
from tclab import runexperiment, labtime
[17]:
def send_sine_wave(t, lab):
print(f'\rTime: {t} Last sleep: {labtime.lastsleep} ', end='')
step = numpy.max(numpy.nonzero(switch_times <= t))
lab.Q1(Qbar + A*numpy.sin(freqs[step]*(t - switch_times[step])))
[18]:
totaltime = sum(nperiods*P)
totaltime
[18]:
29004.612377178186
Experiment will take this many hours
[19]:
totaltime/60/60
[19]:
8.056836771438386
Looks like we’ll be running it overnight.
[57]:
%%time
experiment = runexperiment(send_sine_wave, connected=True,
plot=False, twindow=1000,
time=int(totaltime),
speedup=1,
dbfile='sinetest.db')
TCLab version 0.4.6dev
NHduino connected on port /dev/cu.wchusbserial1410 at 115200 baud.
TCLab Firmware 1.3.0 Arduino Uno.
Time: 4872.0 Last sleep: 0.8694100379943848 TCLab disconnected successfully.
CPU times: user 15.5 s, sys: 28 s, total: 43.6 s
Wall time: 1h 21min 16s
[37]:
import tclab
[38]:
h = tclab.Historian((('Q1', lambda: [0, 0, 0, 0]),
('Q2', None),
('T1', None),
('T2', None)), dbfile='sinetest.db')
[39]:
# h = experiment.historian
[40]:
h.get_sessions()
[40]:
[(2, '2018-03-06 18:45:27', 13710),
(12, '2018-03-07 14:55:43', 2001),
(15, '2018-03-07 18:43:39', 7526),
(25, '2018-03-08 05:34:09', 5523),
(27, '2018-03-08 07:10:23', 4873),
(28, '2018-03-08 12:59:29', 55),
(29, '2018-03-08 13:00:31', 116),
(30, '2018-03-08 13:02:35', 1001),
(31, '2018-03-08 13:25:17', 2001),
(32, '2018-03-08 14:30:19', 0),
(33, '2018-03-08 14:30:32', 891),
(34, '2018-03-08 14:46:08', 536),
(35, '2018-03-08 14:55:18', 132),
(36, '2018-03-08 15:02:28', 2001),
(37, '2018-03-09 04:37:56', 0),
(38, '2018-03-09 04:39:17', 0),
(39, '2018-03-09 04:45:49', 0),
(40, '2018-03-09 04:46:44', 0),
(41, '2018-03-09 08:22:48', 0),
(42, '2018-03-09 08:23:20', 0)]
[41]:
import pandas
[42]:
%matplotlib inline
[43]:
# h = experiment.historian
[44]:
sine_sessions = [2, 15, 25, 27]
[45]:
def sinefit(session, freq, Tbar=40, inphase=0, gain=0.4, phase=0):
h.load_session(session)
t = numpy.array(h.t)
Q1 = numpy.array(h.logdict['Q1'])
Q1_sine = A*numpy.sin(t*freq + inphase)
T1 = numpy.array(h.logdict['T1'])
T1_sine = A*gain*numpy.sin(t*freq + inphase + phase)
plt.plot(t, Q1 - Qbar, color='blue', alpha=0.4)
plt.plot(t, Q1_sine, color='blue')
plt.plot(t, T1 - Tbar, color='red', alpha=0.4)
plt.plot(t, T1_sine, color='red')
plt.ylim(-A, A)
print(f'Gain={gain}, Phase={phase-inphase}')
[46]:
sinefit(2, freqs[0], 40, 0, 0.3, 0, )
Gain=0.3, Phase=0
[47]:
from ipywidgets import interact
[48]:
interact(sinefit,
session=sine_sessions,
freq=freqs,
Tbar=(30., 50.),
inphase=(0, 2*numpy.pi),
gain=(0., 1., 0.01),
phase=(-2*numpy.pi, 0),)
[48]:
<function __main__.sinefit>
[124]:
gains = [0.3, 0.32, 0.25, 0.24, 0.19]
phases = [-0.28, -0.38, -0.58, -0.88, -1.28]
[158]:
gainax, phaseax = bode(omega, G)
gainax.scatter(freqs, gains, color='red')
phaseax.scatter(freqs, phases, color='red')
[158]:
<matplotlib.collections.PathCollection at 0x11a430ba8>
2.2. FFT based bode diagram
[20]:
next_values = {'time': 0, 'value': +A}
def random_noise(t, lab):
if t > next_values['time']:
next_values['time'] += numpy.random.uniform(50, 300)
next_values['value'] = -A if next_values['value'] == +A else +A
lab.Q1(Qbar + next_values['value'])
[21]:
next_values
[21]:
{'time': 0, 'value': 10}
[22]:
%matplotlib notebook
[23]:
%%time
experiment = runexperiment(random_noise, connected=True,
plot=True, twindow=2000,
time=8*60*60,
speedup=1,
dbfile='sinetest.db',
)
TCLab version 0.4.6dev
NHduino connected on port /dev/cu.wchusbserial1410 at 115200 baud.
TCLab Firmware 1.3.0 Arduino Uno.
TCLab disconnected successfully.
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
<timed exec> in <module>()
~/Documents/Development/TCLab/tclab/experiment.py in runexperiment(function, connected, plot, twindow, time, dbfile, speedup, synced)
99 """
100 with Experiment(connected, plot, twindow, time, dbfile, speedup, synced) as experiment:
--> 101 for t in experiment.clock():
102 function(t, experiment.lab)
103 return experiment
~/Documents/Development/TCLab/tclab/experiment.py in clock(self)
78 else:
79 times = range(self.time)
---> 80 for t in times:
81 yield t
82 if self.plot:
~/Documents/Development/TCLab/tclab/labtime.py in clock(period, step, tol, adaptive)
98 'Step size was {} s, but {:.2f} s elapsed '
99 '({:.2f} too long). Consider increasing step.')
--> 100 raise RuntimeError(message.format(step, elapsed, elapsed-step))
101 labtime.sleep(step - (labtime.time() - start) % step)
102 now = labtime.time() - start
RuntimeError: Labtime clock lost synchronization with real time. Step size was 1 s, but 3.22 s elapsed (2.22 too long). Consider increasing step.
[181]:
h = experiment.historian
[30]:
h.load_session(36)
[31]:
%matplotlib inline
[32]:
plt.plot(h.t, h.logdict['T1'])
[32]:
[<matplotlib.lines.Line2D at 0x118ee7a20>]
[33]:
import numpy
[72]:
startcut = 900
[73]:
Q1 = numpy.array(h.logdict['Q1'][startcut:])
T1 = numpy.array(h.logdict['T1'][startcut:])
[74]:
plt.plot(Q1)
plt.plot(T1)
[74]:
[<matplotlib.lines.Line2D at 0x1176f1278>]
[75]:
u = numpy.fft.rfft(Q1 - Qbar)
y = numpy.fft.rfft(T1 - numpy.mean(T1))
[76]:
omegafft = numpy.fft.rfftfreq(len(Q1), 1)
Gfft = y/u
[160]:
ga, pa = bode(omega, G)
bode(omegafft, Gfft, ga, pa, numpy.pi/4)
plt.xlim(1e-3, 1e-2)
ga.set_ylim(0.01, 1)
pa.set_ylim(-numpy.pi/2, 0)
[160]:
(-1.5707963267948966, 0)
[ ]:
[ ]: