bmarechal / py-ident

Download as
Browse Code ยป

Commit f27e91828c8845706a48e89f123fdd436906fe69

Authored by mer0m
Committed by GitHub
1 parent 9259d78054
Exists in master

Add files via upload

Showing 2 changed files with 119 additions and 0 deletions Side-by-side Diff

  1 +# -*- coding: utf-8 -*-
  2 +"""
  3 +Created on Wed Mar 16 10:22:54 2016
  4 +
  5 +@author: Baptiste Marechal
  6 +"""
  7 +
  8 +import csv, time, glob, datetime, os
  9 +import matplotlib.pyplot as plt
  10 +import matplotlib.dates as md
  11 +
  12 +#from scipy import signal
  13 +#from scipy.optimize import curve_fit
  14 +
  15 +tic = time.time()
  16 +
  17 +os.chdir('/home/user/sicav_data/Manip/2016/2016-03/')
  18 +
  19 +def getColumn(filename, column):
  20 + results = []
  21 + for dat_file in sorted(glob.glob(filename)):
  22 + file_result = csv.reader(open(dat_file), delimiter='\t')
  23 + results = results + map(float,[result[column] for result in file_result])
  24 + return results
  25 +
  26 +t = getColumn('*-lakeshore.dat', 0)
  27 +T1 = getColumn('*-lakeshore.dat', 2)
  28 +T2 = getColumn('*-lakeshore.dat', 3)
  29 +T3 = getColumn('*-lakeshore.dat', 4)
  30 +T4 = getColumn('*-lakeshore.dat', 5)
  31 +
  32 +"""
  33 +n = 400
  34 +t = [t[n*i] for i in range(len(t)/n)]
  35 +T1 = [T1[n*i] for i in range(len(T1)/n)]
  36 +T2 = [T2[n*i] for i in range(len(T2)/n)]
  37 +T3 = [T3[n*i] for i in range(len(T3)/n)]
  38 +T4 = [T4[n*i] for i in range(len(T4)/n)]
  39 +
  40 +def func(U, a0, b1, b2, y0):
  41 + sys = signal.lti([a0],[b2, b1, 1])
  42 + y = sys.output(U, t, y0)
  43 + return y[1]
  44 +
  45 +print('Fitting...')
  46 +popt, cov = curve_fit(func, T1, T4)
  47 +print(popt)
  48 +Yfit = func(T1, *popt)
  49 +"""
  50 +timetamps = [datetime.datetime.fromtimestamp(ti) for ti in t]
  51 +datenums=md.date2num(timetamps)
  52 +
  53 +plt.subplots_adjust(bottom=0.35)
  54 +plt.xticks(rotation=90)
  55 +ax=plt.gca()
  56 +xfmt = md.DateFormatter('%Y-%m-%d %H:%M:%S')
  57 +ax.xaxis.set_major_formatter(xfmt)
  58 +
  59 +plt.plot(datenums, T1, label = 'Table')
  60 +plt.plot(datenums, T2, label = 'Link st.')
  61 +plt.plot(datenums, T3, label = 'PT2')
  62 +plt.plot(datenums, T4, label = 'Reg. st.')
  63 +#plt.plot(datenums, Yfit, label = 'Fit')
  64 +
  65 +plt.legend()
  66 +plt.grid()
  67 +plt.show()
  68 +
  69 +toc = time.time() - tic
  70 +print(toc)
first_order_ident.py
Diff comments   View file @ f27e918
  1 +# -*- coding: utf-8 -*-
  2 +"""
  3 +Created on Fri Mar 11 11:04:35 2016
  4 +
  5 +@author: Baptiste Marechal
  6 +"""
  7 +
  8 +from scipy import signal, linspace, pi, randn, ones
  9 +from scipy.optimize import curve_fit
  10 +import matplotlib.pyplot as plt
  11 +import time
  12 +
  13 +tic = time.time()
  14 +
  15 +'''function to optimize'''
  16 +def func(u, K, tau, y0):
  17 + sys = signal.lti(K, [tau, 1])
  18 + y = sys.output(u, t, y0)
  19 + return y[1]
  20 +
  21 +'''input square signal'''
  22 +global t
  23 +t = linspace(0, 20, 1000)
  24 +Umes = -0.5*(signal.square(2*pi*0.1*t)+ones(len(t)))+ones(len(t))+randn(len(t))/50
  25 +
  26 +'''noisy output signal'''
  27 +p = [1, 1, 0]
  28 +Ymes = func(Umes, *p)+randn(len(Umes))/50
  29 +
  30 +'''input and output signals plot'''
  31 +plt.plot(t, Umes, label = 'Umes')
  32 +plt.plot(t, Ymes, label = 'Ymes')
  33 +
  34 +
  35 +'''optimization with non-linear least squares method'''
  36 +popt, cov = curve_fit(func, Umes, Ymes)
  37 +print(p)
  38 +print(popt)
  39 +
  40 +'''estimated response plot'''
  41 +Yfit = func(Umes, *popt)
  42 +plt.plot(t, Yfit, label ='Yfit')
  43 +
  44 +plt.legend()
  45 +plt.grid()
  46 +plt.show()
  47 +
  48 +toc = time.time() - tic
  49 +print(toc)