# HG changeset patch # User Puneeth Chaganti # Date 1306486499 -19800 # Node ID 0efde00f9229164ccb0df22bc9f0ad9126d83c05 Initial commit. diff -r 000000000000 -r 0efde00f9229 .hgignore --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/.hgignore Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,34 @@ +# use glob syntax. +syntax: glob + +*.aux +*.dvi +*.log +*.nav +*.snm +*.toc +*.pdf +*.vrb +*.out +*.sty +*.pyc +*.zip +*~ +.project +.pydevproject +app.yaml +build +tests/coverageResults +*,cover +tests/.coverage +*.git +*.egg-info +eggs +comments +parts +.installed.cfg +bin +develop-eggs +.gitignore +.DS_Store +index.yaml \ No newline at end of file diff -r 000000000000 -r 0efde00f9229 Z-trans/aconv1.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Z-trans/aconv1.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,13 @@ +#!/usr/bin/python +# 4.1 + +import scipy as sp +from pylab import * +a = 0.9 +n = sp.arange(-10,21) +y = a**n*(n>=0) +stem(n,y,markerfmt='.') +xlabel('Time(n)') +ylabel(r'$0.9^n1(n)$') + +show() diff -r 000000000000 -r 0efde00f9229 Z-trans/aconv2.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Z-trans/aconv2.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,12 @@ +#!/usr/bin/python +# 4.2 + +import scipy as sp +from pylab import * +a = 0.9 +n = sp.arange(-10,21) +y = -(a**n*(n<=-1)) +stem(n,y,markerfmt='.') +xlabel('Time(n)') +ylabel(r'$-(0.9)^n1(-n-1)$') +show() diff -r 000000000000 -r 0efde00f9229 Z-trans/disc1.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Z-trans/disc1.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,23 @@ +#!/usr/bin/python +## 4.4 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import scipy as sp +from scipy import signal +from dscr import dscr + +F = sp.array([[0,0],[1,-0.1]]) +G = sp.array([[0.1],[0]]) +C = sp.array([0,1]) +D = 0 +Ts = 0.2 + +H = signal.lti(F, G, C, D) + +H = dscr(H,Ts) +n, d = H.num, H.den + +print "Numerator: ", n +print "Denominator: ", d +print "Sampling Time: ", Ts diff -r 000000000000 -r 0efde00f9229 Z-trans/disc2.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Z-trans/disc2.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,20 @@ +#!/usr/bin/python +## 8.1 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import scipy as sp +from scipy import signal +from dscr import dscr + +n = 10 +d = [5, 1] +Ts = 0.5 +H = signal.lti(n, d) + +H = dscr(H,Ts) +n, d = H.num, H.den + +print "Numerator: ", n +print "Denominator: ", d +print "Sampling Time: ", Ts diff -r 000000000000 -r 0efde00f9229 Z-trans/division.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Z-trans/division.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,15 @@ +#!/usr/bin/python +## 4.11 + +import scipy as sp +from scipy import signal + +num = sp.array([11,-15,6]) +den = signal.convolve([1,-2],signal.convolve([1,-1],[1,-1])) +u = sp.zeros(5) +u[0] = 1 +y = signal.lfilter(num,den,u) + +print num, den + +#Plot of impulse response??? diff -r 000000000000 -r 0efde00f9229 Z-trans/pz.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Z-trans/pz.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,21 @@ +#!/usr/bin/python +# 4.3 + +import pylab as pl + +z = pl.array([[0],[5./12]]) +p = pl.array([[1./2],[1./3]]) + +def zplane(z, p): + pl.figure(figsize=(8,8)) + x = pl.linspace(-1.5, 1.5) + pl.plot(pl.zeros_like(x), x, 'b.', x, pl.zeros_like(x), 'b.', markersize=2) + x = pl.linspace(0, 2*pl.pi, 100) + pl.plot(pl.cos(x), pl.sin(x), 'b.', markersize=2) + pl.plot(z, pl.zeros_like(z), 'ro', p, pl.zeros_like(p), 'rx', markersize=8) + pl.xlabel('Real part') + pl.ylabel('Imaginary part') + +zplane(z, p) + +pl.show() diff -r 000000000000 -r 0efde00f9229 Z-trans/respol.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Z-trans/respol.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,9 @@ +#!/usr/bin/python +## 4.5 + +from scipy import signal + +def respol(num,den): + if num[-1] == 0: + num = num[:-1] + return signal.residuez(num,den)[:-1] diff -r 000000000000 -r 0efde00f9229 Z-trans/respol1.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Z-trans/respol1.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,12 @@ +#!/usr/bin/python +## 4.6 + +from respol import respol + +import scipy as sp +from scipy import signal + +num = sp.array([2,2,0]) +den = sp.array([1,2,-3]) + +res, pol = respol(num,den) diff -r 000000000000 -r 0efde00f9229 Z-trans/respol2.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Z-trans/respol2.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,12 @@ +#!/usr/bin/python +## 4.7 + +from respol import respol + +import scipy as sp +from scipy import signal + +num = sp.array([1,1,0]) +den = signal.convolve([1,-1],signal.convolve([1,-1],[1,-1])) + +res, pol = respol(num,den) diff -r 000000000000 -r 0efde00f9229 Z-trans/respol3.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Z-trans/respol3.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,12 @@ +#!/usr/bin/python +## 4.8 + +from respol import respol + +import scipy as sp +from scipy import signal + +num = sp.array([11,-15,6]) +den = signal.convolve([1,-2],signal.convolve([1,-1],[1,-1])) + +res, pol = respol(num,den) diff -r 000000000000 -r 0efde00f9229 Z-trans/respol5.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Z-trans/respol5.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,12 @@ +#!/usr/bin/python +## 4.9 + +from respol import respol + +import scipy as sp +from scipy import signal + +num = sp.array([1,2,0]) +den = signal.convolve(signal.convolve([1,1],[1,1]), [1,-2]) + +res, pol = respol(num,den) diff -r 000000000000 -r 0efde00f9229 Z-trans/respol6.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Z-trans/respol6.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,12 @@ +#!/usr/bin/python +## 4.10 + +from respol import respol + +import scipy as sp +from scipy import signal + +num = sp.array([3, -5/6.]) +den = signal.convolve([1,-1/4.],[1,-1/3.]) + +res, pol, coeffs = signal.residuez(num,den) diff -r 000000000000 -r 0efde00f9229 chap2/mat_exp.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/chap2/mat_exp.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,16 @@ +#!/usr/bin/python +## 2.2 + +import scipy as sp +from scipy import linalg + +F=sp.array([[-1,0],[1,0]]) +print linalg.expm(F) + +# """ +# 1. linalg.expm computes the matrix exponential using Pade approximation. +# 2. linalg.expm2 computes the matrix exponential using eigenvalue decomposition. +# 3. linalg.expm3 computes the matrix exponential using Taylor series. +# """ + + diff -r 000000000000 -r 0efde00f9229 chap2/zoh1.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/chap2/zoh1.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,25 @@ +#!/usr/bin/python +## 2.3 + +import os, sys +curdir = os.getcwdu() +n = curdir.rfind("/") +curdir[:n]+"/python" +sys.path += [curdir[:n]+"/python"] + +import scipy as sp +from scipy import signal +from scipy import linalg + +from dscr import dscr + + +F = sp.array([[-1,0],[1,0]]) +G = sp.array([[1],[0]]) +C = sp.array([0,1]) +D = 0 +Ts = 1 +H = signal.lti(F, G, C, D) + +h = dscr(H, Ts) +print h.A, h.B, h.C, h.D, Ts diff -r 000000000000 -r 0efde00f9229 freq/derv_bode.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/freq/derv_bode.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,19 @@ +#!/usr/bin/python +# 5.4 + +import scipy as sp +from pylab import * + +w = arange(0, pi, 0.01) +G = 1-exp(-1j*w) +subplot(2,1,1) +loglog(w,abs(G)) +ylabel('Magnitude') + +subplot(2,1,2) +semilogx(w, 180*angle(G)/pi) +xlabel('w (rad/s)') +ylabel('Phase') + +show() +#axis tight, label('',18,'w (rad/s)','Phase',18) diff -r 000000000000 -r 0efde00f9229 freq/filter1.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/freq/filter1.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,26 @@ +#!/usr/bin/python +# 5.2 + +import scipy as sp +from pylab import * + +omega = linspace(0, pi, 100) +g1 = 0.5 / (cos(omega)-0.5+1j*sin(omega)) +mag1 = abs(g1) +angle1 = angle(g1) * 180/pi +g2 = (0.5+0.5*cos(omega)-1.5*1j*sin(omega)) \ + * 0.25 / (1.25-cos(omega)) +mag2 = abs(g2) +angle2 = angle(g2) * 180/pi +subplot(2,1,1) +plot(omega,mag1,omega,mag2,'--') + +ylabel('Magnitude') + +subplot(2,1,2) +plot(omega,angle1,omega,angle2,'--') +xlabel('w (rad/s)') +ylabel('Phase') + +show() +#axis tight, label('',18,'w (rad/s)','Phase',18) diff -r 000000000000 -r 0efde00f9229 freq/incr_freq.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/freq/incr_freq.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,17 @@ +#!/usr/bin/python +# 5.1 + +import scipy as sp +from pylab import * + +n = sp.arange(17) +subplot(2,2,1), stem(n,cos(n*pi/8)) +grid,xlabel('n'),ylabel('cos(n*pi/8)') +subplot(2,2,2), stem(n,cos(n*pi/4)) +grid,xlabel('n'),ylabel('cos(n*pi/4)') +subplot(2,2,3), stem(n,cos(n*pi/2)) +grid,xlabel('n'),ylabel('cos(n*pi/2)') +subplot(2,2,4), stem(n,cos(n*pi)) +grid,xlabel('n'),ylabel('cos(n*pi)') + +show() diff -r 000000000000 -r 0efde00f9229 freq/lead_exp.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/freq/lead_exp.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,21 @@ +#!/usr/bin/python +# 7.3 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import pylab as pl +from scipy import signal +from bode import bode + +pol1 = [1, -0.9] +pol2 = [1, -0.8] +G1 = signal.lti(pol1,[1, 0]) +G2 = signal.lti(pol2,[1, 0]) +w = pl.linspace(0.001,0.5,1000) +bode(G1, w) +bode(G2, w) +pl.figure() +G = signal.lti(pol1,pol2) +bode(G,w) + +pl.show() diff -r 000000000000 -r 0efde00f9229 freq/lead_lag.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/freq/lead_lag.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,23 @@ +#!/usr/bin/python +# 7.4 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import pylab as pl +from scipy import signal +from bode import bode + +w = pl.linspace(0.001,pl.pi,10000) +a = pl.linspace(0.001,0.999,100) + +zero = a +pole = 0.9*a +A = pl.vstack((pl.ones_like(zero), -zero)) +A = A.T +B = pl.vstack((pl.ones_like(pole), -pole)) +B = B.T +sys = map(signal.lti, A, B) + + + +#pl.show() diff -r 000000000000 -r 0efde00f9229 freq/lead_vfy.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/freq/lead_vfy.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,13 @@ +#!/usr/bin/python +# 7.5 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import pylab as pl +from scipy import signal +from bode import bode + +w = pl.linspace(0.0001, 0.5, 1000) +G = signal.lti([1, -0.8], [1, -0.24]) +bode(G, w) +pl.show() diff -r 000000000000 -r 0efde00f9229 freq/ma_bode.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/freq/ma_bode.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,18 @@ +#!/usr/bin/python +# 5.3 + +import scipy as sp +from pylab import * + +w = arange(0, pi, 0.01) +subplot(2,1,1) +loglog(w,abs(1+2*cos(w))/3) +ylabel('Magnitude') + +subplot(2,1,2) +semilogx(w,angle(1+2*cos(w))*180/pi) +xlabel('w (rad/s)') +ylabel('Phase') + +show() +#axis tight, label('',18,'w (rad/s)','Phase',18) diff -r 000000000000 -r 0efde00f9229 freq/nmp.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/freq/nmp.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,26 @@ +#!/usr/bin/python +# 5.2 + +import scipy as sp +from pylab import * + +omega = linspace(0, pi, 100) +ejw = exp(-1j*omega) +G1 = 1.5*(1-0.4*ejw) +mag1 = abs(G1) +angle1 = angle(G1)*180/pi +G2 = -0.6*(1-2.5*ejw) +mag2 = abs(G2) +angle2 = angle(G2)*180/pi + +subplot(2,1,1) +plot(omega,mag1,omega,mag2,'--') +ylabel('Magnitude') + +subplot(2,1,2) +plot(omega,angle1,omega,angle2,'--') +xlabel('w (rad/s)') +ylabel('Phase') + +show() +#axis tight, label('',18,'w (rad/s)','Phase',18) diff -r 000000000000 -r 0efde00f9229 ident/LS_ex.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/ident/LS_ex.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,19 @@ +#!/usr/bin/python +# 6.1 + +import scipy as sp +from pylab import * +from scipy import linalg + +Mag = 10 +V = 10 +No_pts = 100 +theta = 2 +Phi = Mag * (1-2*random([No_pts,1])) +E = V * (1-2*random([No_pts,1])) +Z = Phi*theta + E +LS, res, rank, s = linalg.lstsq(Phi,Z) +Max = max(Z/Phi) +Min = min(Z/Phi) + +print LS, Max, Min diff -r 000000000000 -r 0efde00f9229 ident/autocov.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/ident/autocov.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,6 @@ +from scipy import signal + +def autocov(y): + l = len(y) + nlags = l + x = signal.correlate diff -r 000000000000 -r 0efde00f9229 imc/delay.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/imc/delay.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,40 @@ +#!/usr/bin/python +# 10.1 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import scipy as sp +from zpowk import zpowk +from pp_im import pp_im + +Ts = 1 +B = 0.63 +A = sp.array([1, -0.37]) +k = int(raw_input('Enter the delay as an integer: ')) + +if k<=0: + k = 1 + +zk, dzk = zpowk(k) + +phi = sp.array([1, -0.5]) +delta = 1 + +Rc, Sc, Tc, gamm = pp_im(B, A, k, phi, delta) + +print "Rc =", Rc +print "Sc =", Sc +print "Tc =", Tc +print "gamm =", gamm + +#// simulation parameters for stb_disc.cos +#// y1: 0 to 1; u1: 0 to 1.2 +st = 1.0 #// desired change in setpoint +t_init = 0 #// simulation start time +t_final = 20 # // simulation end time + +#// simulation parameters for stb_disc_10.1.cos +N_var = 0 +C = 0 +D = 1 +N = 1 diff -r 000000000000 -r 0efde00f9229 imc/imc_stable.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/imc/imc_stable.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,36 @@ +#!/usr/bin/python +# 10.9 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import scipy as sp +from imcsplit import imcsplit +from polyfuncs import poladd +from zpowk import zpowk + +def imc_stable(B, A, k, alpha): + """ Designs Discrete Internal Model Controller + for transfer function z^{-k}B(z^{-1})/A(z^{-1}) + Numerator and Denominator of IMC HQ are outputs + Controller is also given in R,S form. """ + Kp, d, Bg, Bnmp, Bm = imcsplit(B, 1) + Bg = Kp * Bg + Bnmpr = Bnmp[::-1] + Bms = sum(Bm) + HiN = A + HiD = Bms * sp.convolve(Bg, Bnmpr) + k = k+d + + zk, dzk = zpowk(k) + zk = sp.squeeze(zk) + Bf = 1-alpha + Af = sp.array([1, -alpha]) + S = sp.convolve(Bf, A) + R1 = sp.convolve(Af, sp.convolve(Bnmpr, Bms)) + R2 = sp.convolve(zk, sp.convolve(Bf, sp.convolve(Bnmp, Bm))) + + R, dR = poladd(R1, len(R1)-1, -R2, len(R2)-1) + R = sp.squeeze(R) + R = sp.convolve(Bg, R) + + return k, HiN, HiD, R, S diff -r 000000000000 -r 0efde00f9229 imc/imc_stable1.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/imc/imc_stable1.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,19 @@ +#!/usr/bin/python +# 10.4 + +import scipy as sp +from imcsplit import imcsplit + +def imc_stable1(B, A, k, alpha): + """ Designs Discrete Internal Model Controller + for transfer function z^{-k}B(z^{-1})/A(z^{-1}) + Numerator and Denominator of IMC HQ are outputs + Controller is also given in R,S form. """ + Kp, d, Bg, Bnmp, Bm = imcsplit(B, 1) + Bg = Kp * Bg + Bnmpr = Bnmp[::] + Bms = sum(Bm) + HiN = A + HiD = Bms * sp.convolve(Bg, Bnmpr) + k = k+d + return k, HiN, HiD diff -r 000000000000 -r 0efde00f9229 imc/imcsplit.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/imc/imcsplit.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,28 @@ +#!/usr/bin/python +# 10.3 + +import scipy as sp + +def imcsplit(B, polynomial): + """ Splits a polynomial B into good, nonminimum with + positive real & with negative real parts. + All are returned in polynomial form. + Gain is returned in Kp and delay in k.""" + k = 0 + Kp = 1 + if polynomial: + roots = sp.roots(B) + Kp = sum(B)/sum(sp.poly(roots)) + else: + roots = B + Bg, Bnmp, Bm = sp.array([1]), sp.array([1]), sp.array([1]) + for root in roots: + if root == 0: + k += 1 + elif abs(root)<1 and root.real>=0: + Bg = sp.convolve(Bg, [1, -root]) + elif abs(root)>=1 and root.real>=0: + Bnmp = sp.convolve(Bnmp, [1, -root]) + else: + Bm = sp.convolve(Bm, [1, -root]) + return Kp, k, Bg, Bnmp, Bm diff -r 000000000000 -r 0efde00f9229 imc/lewin_imc1.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/imc/lewin_imc1.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,29 @@ +#!/usr/bin/python +# 10.7 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import scipy as sp +from scipy import signal +from imc_stable1 import imc_stable1 +from zpowk import zpowk +from myc2d import myc2d + +num = sp.array([1]) +den = sp.array([250, 35, 1]) +G = signal.lti(num, den) + +Ts = 3 +B, A, k = myc2d(G, Ts) + +alpha = 0.9 +k, GiN, GiD = imc_stable1(B, A, k, alpha) +print k, GiN, GiD + +zk, dzk = zpowk(k) +Bp, Ap = B, A + +t0 = 0 +tf = 10 +st = 1 +Nvar = 0 diff -r 000000000000 -r 0efde00f9229 imc/smith.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/imc/smith.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,41 @@ +#!/usr/bin/python +# 10.2 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import scipy as sp +from zpowk import zpowk +from pp_im import pp_im +from polyfuncs import poladd + +Ts = 1 +B = 0.63 +A = sp.array([1, -0.37]) +k = 3 +Bd = sp.convolve(B, [0,1]) +kd = k - 1 + +zkd, dzkd = zpowk(kd) +mzkd, dmzkd = poladd(1, 0, -zkd, dzkd) + +phi = sp.array([1, -0.5]) +delta = 1 + +Rc, Sc, Tc, gamm = pp_im(B, A, k, phi, delta) +print "Rc =", Rc +print "Sc =", Sc +print "Tc =", Tc +print "gamm =", gamm + + +#// simulation parameters for stb_disc.cos +#// y1: 0 to 1; u1: 0 to 1.2 +st = 1.0 #// desired change in setpoint +t_init = 0 #// simulation start time +t_final = 20 # // simulation end time + +#// simulation parameters for stb_disc_10.1.cos +N_var = 0 +C = 0 +D = 1 +N = 1 diff -r 000000000000 -r 0efde00f9229 imc/vande_imc.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/imc/vande_imc.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,29 @@ +#!/usr/bin/python +# 10.10 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import scipy as sp +from scipy import signal +from imc_stable import imc_stable +from zpowk import zpowk +from myc2d import myc2d + +num = sp.array([-1.117, 3.1472]) +den = sp.array([1, 4.6429, 5.3821]) +G = signal.lti(num, den) + +Ts = 0.1 +B, A, k = myc2d(G, Ts) + +alpha = 0.5 +print B, A, k, alpha +k, HiN, HiD, R, S = imc_stable(B, A, k, alpha) +print "k =", k +print "HiN =", HiN +print "HiD =", HiD +print "R =", R +print "S =", S + +zk, dzk = zpowk(k) +Bp, Ap = B, A diff -r 000000000000 -r 0efde00f9229 imc/vande_imc1.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/imc/vande_imc1.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,31 @@ +#!/usr/bin/python +# 10.7 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import scipy as sp +from scipy import signal +from imc_stable1 import imc_stable1 +from zpowk import zpowk +from myc2d import myc2d + +num = sp.array([-1.117, 3.1472]) +den = sp.array([1, 4.6429, 5.3821]) +G = signal.lti(num, den) + +Ts = 0.1 +B, A, k = myc2d(G, Ts) + +alpha = 0.9 +k, GiN, GiD = imc_stable1(B, A, k, alpha) +print "k =", k +print "GiN =", GiN +print "GiD =", GiD + +zk, dzk = zpowk(k) +Bp, Ap = B, A + +t0 = 0 +tf = 10 +st = 1 +Nvar = 0 diff -r 000000000000 -r 0efde00f9229 imc/visc_imc1.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/imc/visc_imc1.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,23 @@ +#!/usr/bin/python +# 10.6 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import scipy as sp +from imc_stable1 import imc_stable1 +from zpowk import zpowk + +B = sp.array([0.51, 1.21]) +A = sp.array([1, -0.44]) +k = 1 +alpha = 0.5 +k, GiN, GiD = imc_stable1(B, A, k, alpha) +print "k =", k +print "GiN =", GiN +print "GiD =", GiD + +Bp, Ap = B, A +Ts = 0.1 +t0 = 0 +tf = 20 +Nvar = 0.01 diff -r 000000000000 -r 0efde00f9229 lqg/lqg.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lqg/lqg.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,61 @@ +#!/usr/bin/python +# 13.4 + +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import scipy as sp +from polyfuncs import polmul, putin, ext +from xdync import xdync +from specfac import specfac + +def lqg(A,degA,B,degB,C,degC,k,rho,V,degV,W,degW,F,degF): + """ LQG controller design by method of Ahlen and Sternad. """ + r, b, degb = specfac(A,degA,B,degB,rho,V,degV,W,degW,F,degF) + print r, b, degb, "specfac" + WFA = sp.convolve(A, sp.convolve(F,W))[::-1] + dWFA = degW + degF + degA + rhs1, drhs1 = polmul(W, degW, WFA, dWFA) + rhs1, drhs1 = polmul(rhs1, drhs1, C, degC) + rhs1 = rho * rhs1 + rhs2 = sp.convolve(C, sp.convolve(V, sp.convolve(B,V)[::-1])) + drhs2 = degC + 2*degV + degB + for i in range(degb - degB - degV): + rhs2 = sp.convolve(rhs2, [0,1]) + drhs2 = drhs2 + degb - degB - degV + C1 = sp.zeros((1,2)) + + C1, degC1 = putin(C1,0,rhs1,drhs1,0,0) + C1, degC1 = putin(C1,degC1,rhs2,drhs2,0,1) + rbf = r * b[::-1] + D1 = sp.zeros((2, 2)) + D1, degD1 = putin(D1,0,rbf,degb,0,0) + for i in range(k): + rbf = sp.convolve(rbf, [0, 1]) + D1, degD1 = putin(D1,degD1,rbf,degb+k,1,1) + N = sp.zeros((1, 2)) + N, degN = putin(N,0,-B,degB,0,0) + + AF, dAF = polmul(A,degA,F,degF) + N, degN = putin(N,degN,AF,dAF,0,1) + print N, degN, "N, degN" + Y, degY, X, degX = xdync(N, degN, D1, degD1, C1, degC1)[:4] + + R, degR = ext(X,degX,0,0) + S, degS = ext(X,degX,0,1) + X = Y[::-1] + + return R, degR, S, degS + +if __name__ == "__main__": + A, dA = sp.array([1, -0.44]), 1 + B, dB = sp.array([0.51, 1.21]), 1 + F, dF = sp.array([1, -1]), 1 + C, dC = sp.array([1, -0.44]),1 + k = 1 + V, W = 1, 1 + dV, dW = 0, 0 + rho = 1 + int1 = 1 + + print lqg(A,dA,B,dB,C,dC,k,rho,V,dV,W,dW,F,dF) diff -r 000000000000 -r 0efde00f9229 lqg/lqg_as.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lqg/lqg_as.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,61 @@ +#!/usr/bin/python + +import sys +sys.path += ['/media/all/work/digital_control/pycode/python/'] + +import scipy as sp +from polyfuncs import polmul, putin, ext +from xdync import xdync +from specfac import specfac + +def lqg_as(A,degA,B,degB,C,degC,k,rho,V,degV,W,degW,F,degF): + """ LQG controller design by method of Ahlen and Sternad. """ + r, b, degb = specfac(A,degA,B,degB,rho,V,degV,W,degW,F,degF) + WFA = sp.convolve(A, sp.convolve(F,W))[::-1] + dWFA = degW + degF + degA + rhs1, drhs1 = polmul(W, degW, WFA, dWFA) + print rhs1, "rhs1", drhs1 + print C, "C", degC + rhs1, drhs1 = polmul(rhs1, drhs1, C, degC) + rhs1 = rho * rhs1 + rhs2 = sp.convolve(C, sp.convolve(V, sp.convolve(B,V)[::-1])) + drhs2 = degC + 2*degV + degB + for i in range(degb - degB - degV): + rhs2 = sp.convolve(rhs2, [0,1]) + drhs2 = drhs2 + degb - degB - degV + C1 = sp.zeros((1,2)) + + C1, degC1 = putin(C1,0,rhs1,drhs1,0,0) + C1, degC1 = putin(C1,degC1,rhs2,drhs2,0,1) + rbf = r * b[::-1] + D1 = sp.zeros((2, 2)) + D1, degD1 = putin(D1,0,rbf,degb,0,0) + for i in range(k): + rbf = sp.convolve(rbf, [0, 1]) + D1, degD1 = putin(D1,degD1,rbf,degb+k,1,1) + N = sp.zeros((1, 2)) + N, degN = putin(N,0,-B,degB,0,0) + + AF, dAF = polmul(A,degA,F,degF) + N, degN = putin(N,degN,AF,dAF,0,1) + + Y, degY, X, degX = xdync(N, degN, D1, degD1, C1, degC1)[:4] + + R, degR = ext(X,degX,0,0) + S, degS = ext(X,degX,0,1) + X = Y[::-1] + + return R, degR, S, degS + +if __name__ == "__main__": + A, dA = sp.array([1, -0.44]), 1 + B, dB = sp.array([0.51, 1.21]), 1 + F, dF = sp.array([1, -1]), 1 + C, dC = sp.array([1, -0.44]),1 + k = 1 + V, W = 1, 1 + dV, dW = 0, 0 + rho = 1 + int1 = 1 + + print lqg_as(A,dA,B,dB,C,dC,k,rho,V,dV,W,dW,F,dF) diff -r 000000000000 -r 0efde00f9229 lqg/lqg_as1.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lqg/lqg_as1.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,23 @@ +#!/usr/bin/python +# 13.8 +# Solves Example 3.1 of Ahlen and Sternad in Hunt's book + +import scipy as sp +from lqg import lqg +from cl import cl + +A, dA = sp.array([1, -0.9]), 1 +B, dB = sp.array([0.1, 0.08]), 1 +F, dF = 1, 0 +C, dC = 1, 0 +k = 2 +V, dV = 1, 0 +W, dW = sp.array([1, -1]), 1 +rho = 0.1 + +R1,dR1,Sc,dSc = lqg_simple(A,dA,B,dB,C,dC,k,rho,V,dV,W,dW,F,dF) + +print "R1", R1 +print "dR1", dR1 +print "Sc", Sc +print "dSc", dSc diff -r 000000000000 -r 0efde00f9229 lqg/lqg_julien1.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lqg/lqg_julien1.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,22 @@ +#!/usr/bin/python + +import scipy as sp +from lqg import lqg +from cl import cl + +A, dA = sp.convolve([1, -0.7], [1, 0.9]), 2 +B, dB = sp.convolve([0.079, 0.221], [1, 0.9]), 2 +C, dC = sp.array([1, -0.7]), 1 +k = 2 +int1 = 1 +F, dF = sp.array([1, -1]), 1 +V, dV = 1, 0 +W, dW = 1, 0 +rho = 0.1 + +R1, dR1, Sc, dSc = lqg (A, dA, B, dB, C, dC, k, rho, V, dV, W, dW, F, dF) +#print R1, dR1, Sc, dSc + +Nu, dNu, Du, dDu, Ny, dNy, Dy, dDy, yvar, uvar = \ + cl(A, dA, B, dB, C, dC, k, Sc, dSc, R1, dR1, int1) +print Nu, dNu, Du, dDu, Ny, dNy, Dy, dDy, yvar, uvar diff -r 000000000000 -r 0efde00f9229 lqg/lqg_julien1_loop.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lqg/lqg_julien1_loop.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,36 @@ +#!/usr/bin/python + +import scipy as sp +from lqg import lqg +from cl import cl + +A, dA = sp.convolve([1, -0.7], [1, 0.9]), 2 +B, dB = sp.convolve([0.079, 0.221], [1, 0.9]), 2 +C, dC = sp.array([1, -0.7]), 1 +k = 2 +int1 = 1 +F, dF = sp.array([1, -1]), 1 +V, dV = 1, 0 +W, dW = 1, 0 +u_lqg = [] +y_lqg =[] +uy_lqg = [] +no_points = 101 +rhovector = sp.logspace(-1.63, 1.2, no_points) + +for i, rho in enumerate(rhovector): + print i, rho + R1, dR1, Sc, dSc = lqg (A, dA, B, dB, C, dC, k, rho, V, dV, W, dW, F, dF) + print R1, dR1, Sc, dSc + Nu, dNu, Du, dDu, Ny, dNy, Dy, dDy, yvar, uvar = \ + cl(A, dA, B, dB, C, dC, k, Sc, dSc, R1, dR1, int1) + uvar = sp.atleast_1d(uvar.squeeze()) + yvar = sp.atleast_1d(yvar.squeeze()) + u_lqg = sp.concatenate((u_lqg, uvar)) + y_lqg = sp.concatenate((y_lqg, yvar)) + if sp.size(uy_lqg) == 0: + uy_lqg = sp.array([rho, uvar[0], yvar[0]]) + else: + uy_lqg = sp.vstack((uy_lqg, sp.array([rho, uvar[0], yvar[0]]))) + +plot(u_lqg, y_lqg) diff -r 000000000000 -r 0efde00f9229 lqg/lqg_mac1.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lqg/lqg_mac1.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,24 @@ +#!/usr/bin/python +# 13.5 + +import scipy as sp +from lqg import lqg +from cl import cl + +A, dA = sp.array([1, -1.4, 0.45]), 2 +C, dC = sp.array([1, -0.5]), 1 +B, dB = 0.5 * sp.array([1, -0.9]), 1 +k = 1 +int1 = 0 +F, dF = 1, 0 +V, dV = 1, 0 +W, dW = 1, 0 +rho = 1 + +R1, dR1, Sc, dSc = lqg (A, dA, B, dB, C, dC, k, rho, V, dV, W, dW, F, dF) + +Nu, dNu, Du, dDu, Ny, dNy, Dy, dDy, yvar, uvar = \ + cl(A, dA, B, dB, C, dC, k, Sc, dSc, R1, dR1, int1) + +print R1, dR1, Sc, dSc +print Nu, dNu, Du, dDu, Ny, dNy, Dy, dDy, yvar, uvar diff -r 000000000000 -r 0efde00f9229 lqg/lqg_simple.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lqg/lqg_simple.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,34 @@ +#!/usr/bin/python +# 13.4 + +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import scipy as sp +from polyfuncs import polmul, putin, ext +from xdync import xdync +from zpowk import zpowk +from specfac import specfac + +def lqg_simple(A,dA,B,dB,C,dC,k,rho,V,dV,W,dW,F,dF): + """ LQG controller design by method of Ahlen and Sternad. """ + r, b, db = specfac(A,dA,B,dB,rho,V,dV,W,dW,F,dF) + D, dD = polmul(A, dA, F, dF) + zk, dzk = zpowk(k) + N, dN = polmul(zk, dzk, B, dB) + RHS, dRHS = polmul(C, dC, b, db) + S, dS, R1, dR1 = xdync(N, dN, D, dD, RHS, dRHS)[:4] + return R1, dR1, S, dS + +if __name__ == "__main__": + A, dA = sp.array([1, -0.44]), 1 + B, dB = sp.array([0.51, 1.21]), 1 + F, dF = sp.array([1, -1]), 1 + C, dC = sp.array([1, -0.44]),1 + k = 1 + V, W = 1, 1 + dV, dW = 0, 0 + rho = 1 + int1 = 1 + + print lqg_simple(A,dA,B,dB,C,dC,k,rho,V,dV,W,dW,F,dF) diff -r 000000000000 -r 0efde00f9229 lqg/lqg_visc.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lqg/lqg_visc.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,23 @@ +#!/usr/bin/python +# 13.6 +# Viscosity control problem of MacGregor + +import scipy as sp +from lqg import lqg +from cl import cl + +A, dA = sp.array([1, -0.44]), 1 +B, dB = sp.array([0.51, 1.21]), 1 +F, dF = sp.array([1, -1]), 1 +C, dC = sp.array([1, -0.44]),1 +k = 1 +V, W = 1, 1 +dV, dW = 0, 0 +rho = 1 +int1 = 1 + +R1,dR1,Sc,dSc = lqg(A,dA,B,dB,C,dC,k,rho,V,dV,W,dW,F,dF) + +Nu,dNu,Du,dDu,Ny,dNy,Dy,dDy,yvar,uvar = cl(A,dA,B,dB,C,dC,k,Sc,dSc,R1,dR1,int1) +print Nu,dNu,Du,dDu,Ny,dNy,Dy,dDy,yvar,uvar +print "Nu,dNu,Du,dDu,Ny,dNy,Dy,dDy,yvar,uvar" diff -r 000000000000 -r 0efde00f9229 lqg/lqg_visc_loop.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lqg/lqg_visc_loop.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,36 @@ +#!/usr/bin/python +# 13.6 +# Viscosity control problem of MacGregor + +import scipy as sp +from pylab import plot, show +from lqg import lqg +from cl import cl + +A, dA = sp.array([1, -0.44]), 1 +B, dB = sp.array([0.51, 1.21]), 1 +F, dF = sp.array([1, -1]), 1 +C, dC = sp.array([1, -0.44]),1 +k = 1 +V, W = 1, 1 +dV, dW = 0, 0 +u_lqg = [] +y_lqg =[] +uy_lqg = [] +int1 = 1 + +for rho in sp.arange(0.001, 3, 0.1): + R1,dR1,Sc,dSc = lqg(A,dA,B,dB,C,dC,k,rho,V,dV,W,dW,F,dF) + Nu,dNu,Du,dDu,Ny,dNy,Dy,dDy,yvar,uvar = cl(A,dA,B,dB,C,dC,k,Sc,dSc,R1,dR1,int1) + uvar = sp.atleast_1d(uvar.squeeze()) + yvar = sp.atleast_1d(yvar.squeeze()) + u_lqg = sp.concatenate((u_lqg, uvar)) + y_lqg = sp.concatenate((y_lqg, yvar)) + if sp.size(uy_lqg) == 0: + uy_lqg = sp.array([rho, uvar[0], yvar[0]]) + else: + uy_lqg = sp.vstack((uy_lqg, sp.array([rho, uvar[0], yvar[0]]))) + +plot(u_lqg, y_lqg) +show() + diff -r 000000000000 -r 0efde00f9229 lqg/spec.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lqg/spec.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,38 @@ +#!/usr/bin/python +# 13.1 + +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import scipy as sp +from polyfuncs import polmul, poladd + +def spec(A, dA, B, dB, rho): + AA = rho * sp.convolve(A, A[::-1]) + BB = sp.convolve(B, B[::-1]) + diff1 = dA - dB + dBB = 2 * dB + for i in range(diff1): + BB, dBB = polmul(BB, dBB, sp.array([0, 1]), 1) + rbbr, drbbr = poladd(AA, 2*dA, BB, dBB) + rbbr = rbbr.squeeze() + rts = sp.roots(rbbr) + rtsin = rts[dA:2*dA+1] + b = 1 + for i in range(dA): + b = sp.convolve(b, sp.array([1, -rtsin[i]])) + bbr = sp.convolve(b, b[::-1]) + r = rbbr[0] / bbr[0] + return r,b,rbbr + +if __name__ == "__main__": + A = sp.convolve([-0.5, 1], [-0.9, 1]) + dA = 2 + B = 0.5*sp.array([-0.9, 1]) + dB = 1 + rho = 1 + + r, beta1, sigma = spec(A, dA, B, dB, rho) + print "sigma", sigma + print "beta1", beta1 + print "r", r diff -r 000000000000 -r 0efde00f9229 lqg/spec_ex.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lqg/spec_ex.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,17 @@ +#!/usr/bin/python +# 13.2 + +import scipy as sp +from spec import spec + +A = sp.convolve([-0.5, 1], [-0.9, 1]) +dA = 2 +B = 0.5*sp.array([-0.9, 1]) +dB = 1 +rho = 1 + +r, beta1, sigma = spec(A, dA, B, dB, rho) +print "sigma", sigma +print "beta1", beta1 +print "r", r + diff -r 000000000000 -r 0efde00f9229 lqg/specfac.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lqg/specfac.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,44 @@ +#!/usr/bin/python +# 13.3 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import scipy as sp +from polyfuncs import polmul, poladd + +def specfac(A, dA, B, dB, rho, V, dV, W, dW, F, dF): + """Implements the spectral factorization for use with LQG + control by design method of Ahlen and Sternad""" + + AFW = sp.convolve(A, sp.convolve(W, F)) + dAFW = dA + dF + dW + AFWWFA = rho * sp.convolve(AFW, AFW[::-1]) + BV = sp.convolve(B, V) + dBV = dB + dV + BVVB = sp.convolve(BV, BV[::-1]) + diff = dAFW - dBV + dBVVB = 2*dBV + for i in range(diff): + BVVB, dBVVB = polmul(BVVB, dBVVB, sp.array([0, 1]), 1) + rbb, drbb = poladd(AFWWFA, 2*dAFW, BVVB, dBVVB) + rbb = rbb.squeeze() + rts = sp.roots(rbb) + rtsin = rts[dAFW:2*dAFW+1] + b = 1 + for i in range(dAFW): + b = sp.convolve(b, sp.array([1, -rtsin[i]])) + b = sp.real(b) + bbr = sp.convolve(b, b[::-1]) + r = rbb[0] / bbr[0] + return r, b, dAFW + +if __name__ == "__main__": + A, dA = sp.array([1, -0.44]), 1 + B, dB = sp.array([0.51, 1.21]), 1 + F, dF = sp.array([1, -1]), 1 + V, W = 1, 1 + dV, dW = 0, 0 + rho = 1 + int1 = 1 + print specfac(A,dA,B,dB,rho,V,dV,W,dW,F,dF) + diff -r 000000000000 -r 0efde00f9229 minv/mv.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/minv/mv.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,20 @@ +#!/usr/bin/env python +# 11.5 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import pylab as pl +from xdync import xdync +from polyfuncs import polmul + +def mv(A, dA, B, dB, C, dC, k, int1): + zk = pl.zeros(k+1) + zk[-1] = 1 + if int1 >= 1: + A, dA = polmul([1, -1], 1, A, dA) + Fj, dFj, Ej, dEj = xdync(zk, k, A, dA, C, dC)[:4] + Gk, dGk = polmul(Ek, dEk, B, dB) + S, dS = Fk, dFk + R, dR = Gk, dGk + return S, dS, R, dR + diff -r 000000000000 -r 0efde00f9229 minv/mv_nm.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/minv/mv_nm.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,23 @@ +#!/usr/bin/env python +# 11.9 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import pylab as pl +from xdync import xdync +from zpowk import zpowk +from polyfuncs import polsplit3, polmul + +def mv(A, dA, B, dB, C, dC, k, int1): + zk, dzk = zpowk(k) + Bzk, dBzk = polmul(B, dB, zk, dzk) + + Bg, Bb = polsplit3(B) + Bbr = Bb[::-1] + + RHS = pl.convolve(C, pl.convolve(Bg, Bbr)) + dRHS = len(RHS) - 1 + + Sc, dSc, Rc, dRc = xdync(Bzk, dBzk, A, dA, RHS, dRHS)[:4] + return Sc, dSc, Rc, dRc + diff -r 000000000000 -r 0efde00f9229 minv/mv_visc.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/minv/mv_visc.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,19 @@ +#!/usr/bin/env python +# 11.11 +# Viscosity control problem of MacGregor +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import pylab as pl +import mv_nm as mv_nm +from cl import cl + +A, dA = pl.array([1, -0.44]), 1 +B, dB = pl.array([0.51, 1.21]), 1 +C, dC = pl.array([1, -0.44]), 1 +k, int1 = 1, 1 + +Sc, dSc, Rc, dRc = mv_nm(A, dA, B, dB, C, dC, k, int1) + +[Nu, dNu, Du, dDu, Ny, dNy, Dy, dDy, yvar, uvar] = \ +cl(A, dA, B, dB, C, dC, k, Sc, dSc, Rc, dRc, int1) diff -r 000000000000 -r 0efde00f9229 minv/pm_10.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/minv/pm_10.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,18 @@ +#!/usr/bin/env python +# 11.3 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import pylab as pl +from xdync import xdync + +C, dC = pl.array([1, 0.5]), 1 +A, dA = pl.array([1, -0.6, -0.16]), 2 +j = 2 + +zj = pl.zeros(j+1) +zj[-1] = 1 +Fj, dFj, Ej, dEj = xdync(zj, j, A, dA, C, dC)[:4] +print Fj, dFj, Ej, dEj +print "Fj, dFj, Ej, dEj " + diff -r 000000000000 -r 0efde00f9229 minv/recursion.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/minv/recursion.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,48 @@ +#!/usr/bin/env python +# 11.1 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import pylab as pl +from polyfuncs import poladd + +def recursion(A, dA, C, dC, j): + Fo, dFo = C, dC + Eo, dEo = 1, 0 + A_z = A[1:dA+1] + dA_z = dA-1 + zi, dzi = 1, 0 + for i in range(j-1): + if dFo == 0: + Fn1 = 0 + else: + Fn1 = Fo[1:dFo+1] + dFn1 = max(dFo-1, 0) + Fn2 = -Fo[0]*A_z + dFn2 = dA-1 + Fn, dFn = poladd(Fn1, dFn1, Fn2, dFn2) + Fn = Fn.squeeze() + zi = pl.convolve(zi,[0,1]) + dzi += 1 + En2 = Fn[0]*zi + dEn2 = dzi + En, dEn = poladd(Eo, dEo, En2, dEn2) + En = En.squeeze() + Eo, dEo = En, dEn + Fo, dFo = Fn, dFn + + if dFo == 0: + Fn1 = 0 + else: + Fn1 = Fo[1:dFo+1] + + dFn1 = max(dFo-1, 0) + Fn2 = -Fo[0]*A_z + dFn2 = dA-1 + Fn, dFn = poladd(Fn1, dFn1, Fn2, dFn2) + Fn = Fn.squeeze() + Ej, dEj = Eo, dEo + Fj, dFj = Fn, dFn + + + return Fj, dFj, Ej, dEj diff -r 000000000000 -r 0efde00f9229 minv/recursion_ex1.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/minv/recursion_ex1.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,13 @@ +#!/usr/bin/env python +# 11.2 +import pylab as pl +from recursion import recursion + +C, dC = pl.array([1, 0.5]), 1 +A, dA = pl.array([1, -0.6, -0.16]), 2 +j = 2 + +Fj, dFj, Ej, dEj = recursion(A, dA, C, dC, j) +print Fj, dFj, Ej, dEj +print "Fj, dFj, Ej, dEj " + diff -r 000000000000 -r 0efde00f9229 mpc/gpc_N.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mpc/gpc_N.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,43 @@ +#!/usr/bin/env python +# 12.2 + +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import pylab as pl +from xdync import xdync +from polyfuncs import polmul + +def gpc_N(A, dA, B, dB, k, N1, N2, Nu, rho): + D, dD = pl.array([1, -1]), 1 + AD, dAD = pl.convolve(A, D), dA+1 + zj, dzj = 1, 0 + for i in range(N1+k-1): + zj = sp.convolve(zj, [0,1]) + dzj += 1 + G = pl.zeros((N2-N1+1, Nu+1)) + H1 = pl.zeros((N2-N1+1, k-1+dB)) + H2 = pl.zeros((N2-N1+1, dA+1)) + + for j in range(k+N1, k+N2+1): + zj = pl.convolve(zj, [0,1]) + dzj = dzj + 1 + Fj, dFj, Ej, dEj = xdync(zj, dzj, AD, dAD, 1, 0)[:4] + Gj, dGj = polmul(B, dB, Ej, dEj) + if j-k >= Nu: + G[j-(k+N1),:Nu+1] = Gj[range(j-k, j-k-Nu-1, -1)] + else: + G[j-(k+N1),:j-k+1] = Gj[j-k::-1] + H1[j-(k+N1),:k-1+dB] = Gj[j-k+1:j+dB] + H2[j-(k+N1),:dA+1] = Fj + K = pl.dot(pl.inv(pl.dot(G.T, G) + rho*pl.eye(Nu+1, Nu+1)), G.T) + KH1 = pl.dot(K, H1) + KH2 = pl.dot(K, H2) + R1 = pl.concatenate(([1], KH1[0,:])) + dR1 = len(R1)-1 + Sc = KH2[0,:] + dSc = len(Sc)-1; + Tc = K[0,:] + dTc = len(Tc)-1; + + return K,KH1,KH2,Tc,dTc,Sc,dSc,R1,dR1 diff -r 000000000000 -r 0efde00f9229 mpc/gpc_Nc.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mpc/gpc_Nc.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,50 @@ +#!/usr/bin/env python +# 12.2 + +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import pylab as pl +from xdync import xdync +from polyfuncs import polmul, poladd + +def gpc_Nc(A, dA, B, dB, C, dC, k, N1, N2, Nu, rho): + D, dD = pl.array([1, -1]), 1 + AD, dAD = pl.convolve(A, D), dA+1 + zj, dzj = 1, 0 + for i in range(N1+k-1): + zj = sp.convolve(zj, [0,1]) + dzj += 1 + M = 2*k+N2-2+dB + P = max(k+N2+dA-1, dC-1) + G = pl.zeros((N2-N1+1, Nu+1)) + H1 = pl.zeros((N2-N1+1, M)) + H2 = pl.zeros((N2-N1+1, P+1)) + + for j in range(k+N1, k+N2+1): + zj = pl.convolve(zj, [0,1]) + dzj = dzj + 1 + Fj, dFj, Ej, dEj = xdync(zj, dzj, AD, dAD, C, dC)[:4] + Nj, dNj, Mj, dMj = xdync(zj, dzj, C, dC, 1, 0)[:4] + Gj, dGj = polmul(Mj, dMj, Ej, dEj) + Gj, dGj = polmul(Gj, dGj, B, dB) + Pj, dPj = polmul(Mj, dMj, Fj, dFj) + Pj, dPj = poladd(Nj, dNj, Pj, dPj) + if j-k >= Nu: + G[j-(k+N1),:Nu+1] = Gj[range(j-k, j-k-Nu-1, -1)] + else: + G[j-(k+N1),:j-k+1] = Gj[j-k::-1] + H1[j-(k+N1),:j+k-2+dB] = Gj[j-k+1:2*j+dB-1] + dPj = max(j-1+dA, dC-1) + H2[j-(k+N1),:dPj+1] = Pj + K = pl.dot(pl.inv(pl.dot(G.T, G) + rho*pl.eye(Nu+1, Nu+1)), G.T) + KH1 = pl.dot(K, H1) + KH2 = pl.dot(K, H2) + R1 = pl.concatenate(([1], KH1[0,:])) + dR1 = len(R1)-1 + Sc = KH2[0,:] + dSc = len(Sc)-1; + Tc = K[0,:] + dTc = len(Tc)-1; + + return K, KH1, KH2, Tc, dTc, Sc, dSc, R1, dR1 diff -r 000000000000 -r 0efde00f9229 mpc/gpc_bas.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mpc/gpc_bas.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,41 @@ +#!/usr/bin/env python +# 12.2 + +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import pylab as pl +from xdync import xdync +from polyfuncs import polmul + +def gpc_bas(A, dA, B, dB, N, k, rho): + D, dD = pl.array([1, -1]), 1 + AD, dAD = pl.convolve(A, D), dA+1 + Nu = N+1 + zj, dzj = 1, 0 + G = pl.zeros((Nu, 1)) + H1 = pl.zeros((Nu, k-1+dB)) + H2 = pl.zeros((Nu, dA+1)) + + for j in range(Nu): + zj = pl.convolve(zj, [0,1]) + dzj = dzj + 1 + Fj, dFj, Ej, dEj = xdync(zj, dzj, AD, dAD, 1, 0)[:4] + Gj, dGj = polmul(B, dB, Ej, dEj) + m, n = G.shape + G = pl.column_stack((G, pl.zeros((m, dGj-n)))) + G[j, :dGj] = Gj[dGj-1::-1] + H1[j,:k-1+dB] = Gj[dGj:dGj+k-1+dB] + H2[j,:dA+1] = Fj + + K = pl.dot(pl.inv(pl.dot(G.T, G) + rho*pl.eye(Nu,Nu)), G.T) + KH1 = pl.dot(K, H1) + KH2 = pl.dot(K, H2) + R1 = pl.concatenate(([1], KH1[0,:])) + dR1 = len(R1)-1 + Sc = KH2[0,:] + dSc = len(Sc)-1; + Tc = K[0,:] + dTc = len(Tc)-1; + + return K, KH1, KH2, Tc, dTc, Sc, dSc, R1, dR1 diff -r 000000000000 -r 0efde00f9229 mpc/gpc_col.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mpc/gpc_col.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,46 @@ +#!/usr/bin/env python +# 12.2 + +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import pylab as pl +from xdync import xdync +from polyfuncs import polmul, poladd + +def gpc_col(A, dA, B, dB, C, dC, N, k, rho): + D, dD = pl.array([1, -1]), 1 + AD, dAD = pl.convolve(A, D), dA+1 + Nu = N+1 + zj, dzj = 1, 0 + G = pl.zeros((Nu, 1)) + H1 = pl.zeros((Nu, 2*k+N-2+dB)) + H2 = pl.zeros((Nu, k+N+dA)) + + for j in range(Nu): + zj = pl.convolve(zj, [0,1]) + dzj = dzj + 1 + Fj, dFj, Ej, dEj = xdync(zj, dzj, AD, dAD, C, dC)[:4] + Nj, dNj, Mj, dMj = xdync(zj, dzj, C, dC, 1, 0)[:4] + Gj, dGj = polmul(Mj, dMj, Ej, dEj) + Gj, dGj = polmul(Gj, dGj, B, dB) + Pj, dPj = polmul(Mj, dMj, Fj, dFj) + Pj, dPj = poladd(Nj, dNj, Pj, dPj) + + m, n = G.shape + G = pl.column_stack((G, pl.zeros((m, j+1-n)))) + G[j, :j+1] = Gj[j::-1] + H1[j,:dGj-j] = Gj[j+1:dGj+1] + H2[j,:dPj+1] = Pj + + K = pl.dot(pl.inv(pl.dot(G.T, G) + rho*pl.eye(Nu,Nu)), G.T) + KH1 = pl.dot(K, H1) + KH2 = pl.dot(K, H2) + R1 = pl.concatenate(([1], KH1[0,:])) + dR1 = len(R1)-1 + Sc = KH2[0,:] + dSc = len(Sc)-1; + Tc = K[0,:] + dTc = len(Tc)-1; + + return K, KH1, KH2, Tc, dTc, Sc, dSc, R1, dR1 diff -r 000000000000 -r 0efde00f9229 mpc/gpc_ex11.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mpc/gpc_ex11.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,39 @@ +#!/usr/bin/env python +# 12.1 +# Camacho and Bordon's GPC example; model formation +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import pylab as pl +from xdync import xdync +from polyfuncs import polmul + +A, dA = pl.array([1, -0.8]), 1 +B, dB = pl.array([0.4, 0.6]), 1 +N=3 +k=1 +D, dD = pl.array([1, -1]), 1 +AD, dAD = pl.convolve(A,D), dA+1 +Nu = N+k +zj, dzj = 1, 0 +G = pl.zeros((Nu, 1)) +H1 = pl.zeros((Nu, k-1+dB)) +H2 = pl.zeros((Nu, dA+1)) + +for j in range(Nu): + zj = pl.convolve(zj, [0,1]) + dzj = dzj + 1 + Fj, dFj, Ej, dEj = xdync(zj, dzj, AD, dAD, 1, 0)[:4] + Gj, dGj = polmul(B, dB, Ej, dEj) + m, n = G.shape + G = pl.column_stack((G, pl.zeros((m, dGj-n)))) + G[j, :dGj] = Gj[dGj-1::-1] + H1[j,:k-1+dB] = Gj[dGj:dGj+k-1+dB] + H2[j,0:dA+1] = Fj + +print "G" +print G +print "H1" +print H1 +print "H2" +print H2 diff -r 000000000000 -r 0efde00f9229 mpc/gpc_ex12.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mpc/gpc_ex12.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,27 @@ +#!/usr/bin/env python +# 12.3 +# Camacho and Bordon's GPC example; model formation +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import pylab as pl +from gpc_bas import gpc_bas + +A, dA = pl.array([1, -0.8]), 1 +B, dB = pl.array([0.4, 0.6]), 1 +N=3 +k=1 +rho = 0.8 + +K, KH1, KH2, Tc, dTc, Sc, dSc, R1, dR1 = gpc_bas(A, dA, B, dB, N, k, rho) + +print "K =", K +print "KH1 =", KH1 +print "KH2 =", KH2 +print "Tc =", Tc +print "dTc =", dTc +print "Sc =", Sc +print "dSc =", dSc +print "R1 =", R1 +print "dR1 =", dR1 + diff -r 000000000000 -r 0efde00f9229 mpc/gpc_ex2.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mpc/gpc_ex2.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,28 @@ +#!/usr/bin/env python +# 12.7 + +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import pylab as pl +from gpc_col import gpc_col + +A, dA = pl.array([1, -0.44]), 1 +B, dB = pl.array([0.51, 1.21]), 1 +N = 2 +k = 1 +C, dC = pl.array([1, -0.44]), 1 +rho = 1 + +K, KH1, KH2, Tc, dTc, Sc, dSc, R1, dR1 = gpc_col(A, dA, B, dB, C, dC, N, k, rho) + +print "K =", K +print "KH1 =", KH1 +print "KH2 =", KH2 +print "Tc =", Tc +print "dTc =", dTc +print "Sc =", Sc +print "dSc =", dSc +print "R1 =", R1 +print "dR1 =", dR1 + diff -r 000000000000 -r 0efde00f9229 mpc/gpc_wt.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mpc/gpc_wt.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,26 @@ +#!/usr/bin/env python +# 12.4 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import pylab as pl +from gpc_N import gpc_N + +A, dA = pl.array([1, -0.8]), 1 +B, dB = pl.array([0.4, 0.6]), 1 +k=1 +rho = 0.8 +N1, N2, Nu = 0, 3, 2 + +K, KH1, KH2, Tc, dTc, Sc, dSc, R1, dR1 = gpc_N(A, dA, B, dB, k, N1, N2, Nu, rho) + +print "K =", K +print "KH1 =", KH1 +print "KH2 =", KH2 +print "Tc =", Tc +print "dTc =", dTc +print "Sc =", Sc +print "dSc =", dSc +print "R1 =", R1 +print "dR1 =", dR1 + diff -r 000000000000 -r 0efde00f9229 mpc/gpc_wtc.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mpc/gpc_wtc.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,28 @@ +#!/usr/bin/env python +# 12.7 + +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import pylab as pl +from gpc_Nc import gpc_Nc + +A, dA = pl.array([1, -0.44]), 1 +B, dB = pl.array([0.51, 1.21]), 1 +C, dC = pl.array([1, -0.44]), 1 +k = 1 +rho = 1 +N1, N2, Nu = 0, 2, 0 +K, KH1, KH2, Tc, dTc, Sc, dSc, R1, dR1 = gpc_Nc(A, dA, B, dB, C, dC, k, N1, N2, Nu, rho) + +print "K =", K +print "KH1 =", KH1 +print "KH2 =", KH2 +print "Tc =", Tc +print "dTc =", dTc +print "Sc =", Sc +print "dSc =", dSc +print "R1 =", R1 +print "dR1 =", dR1 + + diff -r 000000000000 -r 0efde00f9229 pid/gmvc_pid.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pid/gmvc_pid.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,43 @@ +#!/usr/bin/env python +# 11.17 + +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import pylab as pl +from xdync import xdync +from zpowk import zpowk +from polyfuncs import polmul + +def gmvc_pid(A, B, k, T, Ts): + A = pl.atleast_1d(A) + B = pl.atleast_1d(B) + dA, dB = len(pl.atleast_1d(A)) - 1, len(pl.atleast_1d(B)) - 1 + dT = len(T) - 1 + if dA > 2: + print 'degree of A cannot be more than 2' + exit(1) + elif dB>1: + print 'degree of B cannot be more than 1' + exit(1) + delta, ddelta = pl.array([1, -1]), 1 + Adelta, dAdelta = polmul(A, dA, delta, ddelta) + Q, dQ, P, dP = xdync(Adelta, dAdelta, B, dB, T, dT)[:4] + P = pl.atleast_1d(P.squeeze()) + PAdelta = P[0]*Adelta + + zk, dzk = zpowk(k) + E, degE, F, degF = xdync(PAdelta, dAdelta, zk, dzk, P, dP)[:4] + E = pl.atleast_1d(E.squeeze()) + F = pl.atleast_1d(F.squeeze()) + nu = P[0]*E[0]*B[0] + Kc = -1/nu*(F[1]+2*F[2]) + tau_i = -(F[1]+2*F[2])/(F[0]+F[1]+F[2])*Ts + tau_d = -F[2]/(F[1]+2*F[2])*Ts + L = pl.empty(3) + L[0] = 1+Ts/tau_i+tau_d/Ts + L[1] = -(1+2*tau_d/Ts) + L[2] = tau_d/Ts + L = Kc * L.T + + return Kc, tau_i, tau_d, L diff -r 000000000000 -r 0efde00f9229 pid/gpc_pid.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pid/gpc_pid.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,78 @@ +#!/usr/bin/env python +# 12.11 + +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import pylab as pl +from xdync import xdync +from zpowk import zpowk +from polyfuncs import polmul, poladd + +def gpc_pid(A, dA, B, dB, C, dC, N1, N2, Nu, lambda1, gamm, gamma_y): + Adelta = pl.convolve(A,[1, -1]) + G = pl.array([]) + Gtilda1 = pl.array([]) + for i in range(N1, N2+1): + zi=zpowk(i)[0] + E, dE, F, dF = xdync(Adelta, dA+1, zi, i, C, dC)[:4] + Gtilda, dGtilda, Gbar, dGbar = xdync(C, dC, zi, i, E*B, dE+dB)[:4] + Gtilda = pl.atleast_1d(Gtilda.squeeze()) + Gtilda1 = pl.empty(i) + for j in range(i): + Gtilda1[j] = Gtilda[i-j-1] + if i <= Nu-1: + Gtilda2 = pl.zeros(Nu) + Gtilda2[:i] = Gtilda1 + if pl.size(G) is not 0: + G = pl.vstack((G, Gtilda2)) + else: + G = Gtilda2 + else: + G = pl.vstack((G, Gtilda1[:Nu])) + + es=sum(pl.atleast_1d(C))/sum(pl.atleast_1d(A)) + gs=sum(pl.atleast_1d(B))/sum(pl.atleast_1d(A)) + F_s=es * A + G_s=[] + + for i in range(1, Nu+1): + if (Nu-i) == 0: + row= gs * pl.ones(i) + else: + row = gs * pl.ones(i) + row = pl.column_stack((row, pl.zeros((Nu-i,Nu-i)))) + row = row.squeeze() + if pl.size(G_s) is not 0: + G_s = pl.row_stack((G_s, row)) + else: + G_s = row + + lambda_mat = lambda1 * pl.identity(Nu) + gamma_mat = gamm * pl.identity(Nu) + gamma_y_mat = gamma_y * pl.identity(N2-N1+1) + mat1 = pl.inv(pl.dot(G.T, pl.dot(gamma_y_mat, G))+lambda_mat+pl.dot(G_s.T, pl.dot(gamma_mat,G_s))) + mat2 = pl.dot(mat1, pl.dot(G.T, gamma_y_mat)) + mat2_s = pl.dot(mat1, pl.dot(G_s.T, gamma_mat)) + h_s= sum(mat2_s[0,:]) + h = mat2[0,:] + T = C + R = C * (sum(h)+h_s) + S = 0 + for i in range (N1, N2+1): + zi=zpowk(i)[0] + E, dE, F, dF = xdync(Adelta, dA+1, zi, i, C, dC)[:4] + Gtilda, dGtilda, Gbar, dGbar = xdync(C, dC, zi, i, E*B, dE+dB)[:4] + S += F*h[i-1] + S += F_s*h_s + S = S.squeeze() + if len(A) == 3: + Kp = S[0] - R - S[2] + Ki = R + Kd = S[2] + else: + Kp = S[1] - R + Ki = R + Kd = 0 + return Kp, Ki, Kd + diff -r 000000000000 -r 0efde00f9229 pid/gpc_pid_test.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pid/gpc_pid_test.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,22 @@ +#!/usr/bin/env python +# 12.10 + +import pylab as pl +from gpc_pid import gpc_pid + +A = pl.array([1, -1.95, 0.935]) +B = -0.015 +C = 1 +degA = 2 +degB = 0 +degC = 0 +N1 = 1 +N2 = 5 +Nu = 2 +gamm = 0.05 +gamma_y = 1 +lambda1 = 0.02 + +Kp, Ki, Kd = gpc_pid(A, degA, B, degB, C, degC, N1, N2, Nu, lambda1, gamm, gamma_y) + +print Kp, Ki, Kd diff -r 000000000000 -r 0efde00f9229 pid/miller.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pid/miller.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,27 @@ +#!/usr/bin/env python +# 11.15 +# GMVC PID tuning of example given by Miller et al. + +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import pylab as pl +from ch_pol import ch_pol +from gmvc_pid import gmvc_pid + +# Model +A = pl.array([1, -1.95, 0.935]) +B = -0.015 +k = 1 +Ts = 1 + +# Transient specifications +N = 15 +epsilon = 0.1 +T = ch_pol(N, epsilon)[:1] + +# Controller Design +Kc, tau_i, tau_d, L = gmvc_pid(A, B, k, T, Ts) + +# L1 = filtval(L,1); +# zk = zpowk(k); diff -r 000000000000 -r 0efde00f9229 place/ball_basic.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/place/ball_basic.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,62 @@ +#!/usr/bin/python +# 9.1 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import scipy as sp +from scipy import signal +from myc2d import myc2d +from pp_basic import pp_basic +from desired import desired + +# Magnetically suspended ball problem +# Operating conditions + +M = 0.05 +L = 0.01 +R = 1 +K = 0.0001 +g = 9.81 + +#Equilibrium conditions +hs = 0.01 +i = sqrt(M*g*hs/K) + + +# State space matrices +a21 = K*i**2/M/hs**2 +a23 = - 2*K*i/M/hs +a33 = - R/L +b3 = 1/L +a1 = sp.array([[0, 1, 0], [a21, 0, a23], [0, 0, a33]]) +b1 = sp.array([[0], [0], [b3]]) +c1 = sp.array([1, 0, 0]) +d1 = 0 + +# Transfer functions +G = signal.lti(a1, b1, c1, d1) +Ts = 0.01 +B, A, k = myc2d(G, Ts) + +#polynomials are returned +num, den = G.num, G.den + +# Transient specifications +rise = 0.15 +epsilon = 0.05 +phi = desired(Ts, rise, epsilon) + +# Controller design +Rc, Sc, Tc, gamm = pp_basic(B, A, k, phi) + +# Setting up simulation parameters for basic +st = 0.0001 # desired change in h, in m. +t_init = 0 # simulation start time +t_final = 0.5 # simulation end time + +# Setting up simulation parameters for c_ss_cl +N_var = 0 +xInitial = [0 0 0] +N = 1 +C = 0 +D = 1 diff -r 000000000000 -r 0efde00f9229 place/ball_im.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/place/ball_im.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,60 @@ +#!/usr/bin/python +# 9.9 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import scipy as sp +from scipy import signal +from myc2d import myc2d +from pp_im import pp_im +from desired import desired + +# Magnetically suspended ball problem +# Operating conditions + +M = 0.05 +L = 0.01 +R = 1 +K = 0.0001 +g = 9.81 + +#Equilibrium conditions +hs = 0.01 +i = sqrt(M*g*hs/K) + + +# State space matrices +a21 = K*i**2/M/hs**2 +a23 = - 2*K*i/M/hs +a33 = - R/L +b3 = 1/L +a = sp.array([[0, 1, 0], [a21, 0, a23], [0, 0, a33]]) +b = sp.array([[0], [0], [b3]]) +c = sp.array([1, 0, 0]) +d = 0 + +# Transfer functions +G = signal.lti(a, b, c, d) +Ts = 0.01 +B, A, k = myc2d(G,Ts) + +# Transient specifications +rise = 0.15 +epsilon = 0.05 +phi = desired(Ts, rise, epsilon) + +# Controller design +Delta = sp.array([1, -1]) +Rc, Sc, Tc, gamm = pp_im(B, A, k, phi) + +# Setting up simulation parameters for basic +st = 0.0001 # desired change in h, in m. +t_init = 0 # simulation start time +t_final = 0.5 # simulation end time + +# Setting up simulation parameters for c_ss_cl +N_var = 0 +xInitial = [0 0 0] +N = 1 +C = 0 +D = 1 diff -r 000000000000 -r 0efde00f9229 place/dof_choice.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/place/dof_choice.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,38 @@ +#!/usr/bin/env python +# 9.17 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import scipy as sp +from pp_pid import pp_pid +from pp_im import pp_im +from zpowk import zpowk + +# test problem to demonstrate benefits of 2_dof + +Ts = 1 +k = 1 +B = sp.convolve([1, 0.9], [1, -0.8]) +A = sp.convolve([1, -1], [1, -0.5]) + +# closed loop characteristic polynomial +phi = [1, -1, 0.5] + +Delta = 1 # Choice of internal model of step +control = 1 +if control == 1: #/ 1-DOF with no cancellation + Rc, Sc = pp_pid(B, A, k, phi, Delta) + Tc = Sc + gamm = 1 +else: #2-DOF + Rc, Sc, Tc, gamm = pp_im(B, A, k, phi, Delta) + +# simulation parameters for stb_disc +zk, dzk = zpowk(k) +st = 1 # desired step change +t_init = 0 # simulation start time +t_final = 20 # simulation end time +xInitial = [0, 0] +C = 0 +D = 1 +N_var = 0 diff -r 000000000000 -r 0efde00f9229 place/ibm_pp.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/place/ibm_pp.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,36 @@ +#!/usr/bin/python +# 9.10 + +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import scipy as sp +from zpowk import zpowk +from desired import desired +from pp_im import pp_im + +# Control of IBM lotus domino server +# Transfer function + +B = 0.47 +A = sp.array([1, -0.43]) +k = 1 +zk, dzk = zpowk(k) + +# Transient specifications +rise = 10 +epsilon = 0.01 +Ts = 1 +phi = desired(Ts,rise,epsilon) + +# Controller design +Delta = sp.array([1, -1]) # internal model of step used +Rc, Sc, Tc, gamm = pp_im(B, A, k, phi, Delta) + +# Simulation parameters for stb_disc.cos +st = 1 # desired change +t_init = 0 # simulation start time +t_final = 40 # simulation end time +C = 0 +D = 1 +N_var = 0 diff -r 000000000000 -r 0efde00f9229 place/motor.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/place/motor.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,40 @@ +#!/usr/bin/python +# 9.11 + +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import scipy as sp +from scipy import signal +from myc2d import myc2d + +# Motor control problem +# Transfer function + +a1 = sp.array([[-1, 0], [1, 0]]) +b1 = sp.array([[1],[0]]) +c1 = sp.array([0, 1]) +d1 = 0 +G = signal.lti(a1, b1, c1, d1) +Ts = 0.25 +B, A, k = myc2d(G,Ts) + +# Transient specifications +rise = 3 +epsilon = 0.05 +phi = desired(Ts,rise,epsilon) + +# Controller design +Delta = 1 # No internal model of step used +Rc, Sc, Tc, gamm = pp_im(B, A, k, phi, Delta) + +# simulation parameters +st = 1 # desired change in position +t_init = 0 # simulation start time +t_final = 10 # simulation end time + +xInitial = [0, 0] # initial conditions +N = 1 +C = 0 +D = 1 +N_var = 0 diff -r 000000000000 -r 0efde00f9229 place/motor_pd.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/place/motor_pd.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,58 @@ +#!/usr/bin/python +# 9.21 + +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import scipy as sp +from myc2d import myc2d +from scipy import signal +from desired import desired +from pd import pd + +# Motor control problem +# Transfer function + +a = sp.array([[-1, 0], [1, 0]]) +b = sp.array([[1],[0]]) +c = sp.array([0, 1]) +d = 0 +G = signal.lti(a,b,c,d) +Ts = 0.25 +B, A, k = myc2d(G,Ts) +num, den = G.num, G.den + +# Transient specifications +rise = 3 +epsilon = 0.05 +phi = desired(Ts,rise,epsilon) + +# Controller design +Delta = 1 # No internal model of step used +Rc, Sc = pp_pid(B, A, k, phi, Delta) + +# continuous time controller +K, taud, N = pd(Rc, Sc, Ts) +numb = K*sp.array([1, taud*(1+1/N)]) +denb = sp.array([1, taud/N]) +numf = 1 +denf = 1 + +# simulation parameters +st = 1 # desired change in position +t_init = 0 # simulation start time +t_final = 20 # simulation end time +st1 = 0 + +# continuous controller simulation: g_s_cl3.cos +num1 = 0 +den1 = 1 + +# discrete controller simulation: g_s_cl2.cos +# u1: -0.1 to 0.8 +# y1: 0 to 1.4 +C = 0 +D = 1 +N = 1 +gamm = 1 +Tc = Sc diff -r 000000000000 -r 0efde00f9229 place/pd.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/place/pd.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,9 @@ +#!/usr/bin/python +# Both Rc and Sc have to be degree one polynomials + +def pd(Rc,Sc,Ts): + K = (Sc[0]+Sc[1])/(1+Rc[1]) + N = (Sc[1]-Sc[0]*Rc[1])/Rc[1]/(Sc[0]+Sc[1]) + taudbyN = -Ts*Rc[1]/(1+Rc[1]) + taud = taudbyN * N + return K, taud, N diff -r 000000000000 -r 0efde00f9229 place/sigurd.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/place/sigurd.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,57 @@ +# Updated(2-8-07) +# 9.13 + +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import scipy as sp +from scipy import signal +from myc2d import myc2d +from dscr import dscr +from zpowk import zpowk +from desired import desired +from pp_im2 import pp_im2 + +num = 200 +den = sp.convolve([0.05, 1], [0.05, 1]) +den = sp.convolve([10, 1], den) +G = signal.lti(num,den) +Ts = 0.025 +num = G.num +den = G.den +B, A, k = myc2d(G, Ts) +zk, dzk = zpowk(k) + +# Transient specifications +a = 0.9 +rise = 0.24 +epsilon = 0.05 +phi = desired(Ts, rise, epsilon) + +# Controller design +Delta = sp.array([1, -1]) # internal model of step is present +Rc, Sc, Tc, gamm = pp_im2(B, A, k, phi, Delta, a) + +# margin calculation +Lnum = sp.convolve(Sc, sp.convolve(B,zk)) +Lden = sp.convolve(Rc, A) +L = signal.lti(Lnum, Lden) +L = dscr(L, Ts) +Gm = g_margin(L) +Pm = p_margin(L) + +num1 = 100 +den1 = [10, 1] +Gd = signal.lti(num1, den1) +C, D, k1 = myc2d(Gd, Ts) +zk, dzk = zpowk(k) +C = sp.convolve(C, zk) + +# simulation parameters g_s_cl2 +N = 1 +st = 1 # desired change in setpoint +st1 = 0 # magnitude of disturbance +t_init = 0 # simulation start time +t_final = 1.5 # simulation end time + + diff -r 000000000000 -r 0efde00f9229 place/type_2DOF.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/place/type_2DOF.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,27 @@ +#!/usr/bin/env python +# 9.16 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import scipy as sp +from zpowk import zpowk +from pp_im import pp_im + +B = 1 +A = [1 -1] +k = 1 +zk = zpowk(k) +Ts = 1 +phi = sp.array([1, -0.5]) + +Delta = 1 # Choice of internal model of step +Rc, Sc, Tc, gamm = pp_im(B, A, k, phi, Delta) + +# simulation parameters for stb_disc +st = 1 # desired step change +t_init = 0 # simulation start time +t_final = 20 # simulation end time +xInitial = [0, 0] +C = 0 +D = 1 +N_var = 0 diff -r 000000000000 -r 0efde00f9229 place/unstb.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/place/unstb.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,39 @@ +#!/usr/bin/python +# 9.7 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import scipy as sp +from scipy import signal +from zpowk import zpowk +from pp_basic import pp_basic +from desired import desired + +Ts = 1 +B = sp.array([1, -3]) +A = sp.array([1, 2, -8]) +k = 1 + +# Since k=1, tf is of the form z^-1 +zk, dzk = zpowk(k) + + +# Transient specifications +rise = 10 +epsilon = 0.1 +phi = desired(Ts, rise, epsilon) + +# Controller design +Rc, Sc, Tc, gamm = pp_basic(B, A, k, phi) + +# Setting up simulation parameters for basic +st = 1.0 # desired change in h, in m. +t_init = 0 # simulation start time +t_final = 1000 # simulation end time + +# Setting up simulation parameters for c_ss_cl +N_var = 0 +N = 1 +C = 0 +D = 1 + diff -r 000000000000 -r 0efde00f9229 python/bode.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/bode.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,41 @@ +import pylab as pl + +def freqresp(H, omega, sys_type='c'): + if (omega == None): + omega = pl.logspace(-2, 2); + + num = pl.poly1d(H.num) + den = pl.poly1d(H.den) + if sys_type == 'd': + fresp = map(lambda w: num(w) / den(w), pl.exp(2*pl.pi*omega*1j)) + else: + fresp = map(lambda w: num(w) / den(w), 2*pl.pi*omega*1j) + fresp = pl.array(fresp) + return omega, fresp + +def phasemag(fresp): + mag = abs(fresp) + mag = 20*pl.log10(mag) + + phase = pl.angle(fresp)*180/pl.pi + + return mag, phase + +def bode(sys, omega=None): + omega, fresp = freqresp(sys, omega, 'd') + mag, phase = phasemag(fresp) + + pl.subplot(211) + pl.grid(linestyle='--') + pl.semilogx(omega, mag) + xmin, xmax = pl.xlim() + pl.xlim(1e-3, xmax) + + pl.subplot(212) + pl.grid(linestyle='--') + pl.semilogx(omega, phase) + xmin, xmax = pl.xlim() + pl.xlim(1e-3, xmax) + + return (211, 212) + diff -r 000000000000 -r 0efde00f9229 python/ch_pol.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/ch_pol.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,15 @@ +#!/usr/bin/env python + +# function [phi,psi] = ch_pol(N,epsilon) +# Returns desired characteristic polynomial and numerator +# N = rise time in number of sample times +# epsilon = overshoot as a fraction of ss gain + +import pylab as pl + +def ch_pol(N,epsilon): + omega = pl.pi/2/N + r = epsilon**(omega/pl.pi) + phi = pl.array([1, -2*r*pl.cos(omega), r**2]) + psi = pl.array([1-r*pl.cos(omega), r**2-r*pl.cos(omega)]) + return phi, psi diff -r 000000000000 -r 0efde00f9229 python/cindep.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/cindep.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,37 @@ +#!/usr/bin/python + +import pylab as pl +from makezero import makezero + +def cindep(S, gap=None): + E = pl.finfo(float).eps + + if gap == None: + gap = 1e8 + + r, c = pl.atleast_2d(S).shape + + if r>c: + ind = 0 + else: + sigma = pl.svd(S)[1] + l = len(sigma) + if(sigma[-1]/sigma[0] < (E*max(1, c))): + ind = 0 + else: + if pl.any(sigma[:-1]/sigma[1:] >= gap): + ind = 0 + else: + ind = 1 + if ind: + b = pl.array([]) + + else: + A, B = S[r-1,:], S[:r-1,:] + b = pl.lstsq(B.T, A.T)[0].T + b = makezero(b, gap) + b = b.squeeze() + return pl.atleast_2d(b) + +if __name__== "__main__": + pass diff -r 000000000000 -r 0efde00f9229 python/cl.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/cl.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,32 @@ +#!/usr/bin/python +# 11.6 + +from zpowk import zpowk +from polyfuncs import polmul, poladd +from tfvar import tfvar + +def cl(A,dA,B,dB,C,dC,k,S,dS,R,dR,int1): + """ + int>=1 means integrated noise and control law: + delta u = - (S/R)y + Evaluates the closed loop transfer function and + variances of input and output + """ + + zk, dzk = zpowk(k) + BS, dBS = polmul(B,dB,S,dS) + zBS, dzBS = polmul(zk,dzk,BS,dBS) + RA, dRA = polmul(R,dR,A,dA) + + if int1>=1: + RA, dRA = polmul(RA, dRA, [1, -1], 1) + D, dD = poladd(RA, dRA, zBS, dzBS) + D = D.squeeze() + Ny, dNy = polmul(C,dC,R,dR) + Nu, dNu = polmul(C,dC,S,dS) + + + Nu, dNu, Du, dDu, uvar = tfvar(Nu, dNu, D, dD) + Ny, dNy, Dy, dDy, yvar = tfvar(Ny, dNy, D, dD) + + return Nu, dNu, Du, dDu, Ny, dNy, Dy, dDy, yvar, uvar diff -r 000000000000 -r 0efde00f9229 python/clcoef.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/clcoef.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,36 @@ +#!/usr/bin/python + +import pylab as pl +from polyfuncs import polsize + +def clcoef(Q, degQ): + """Function to clear zero leading coefficient matrices + of a Polynomial Matrix""" + rQ, cQ = polsize(Q, degQ) + + if pl.allclose(Q, pl.zeros_like(Q)): + P = pl.zeros((rQ,cQ)) + degP = 0 + else: + P, degP = Q, degQ + rP, cP = rQ, cQ + j = degP + while j>=0: + if pl.norm(P[:, -cP:], pl.inf) < 1e-8 *pl.norm(P, pl.inf): + P = P[:, :j*cP] + degP = degP-1 + else: + j=0 + j -= 1 + return P, degP + +if __name__== "__main__": + s = """P = pl.array([[1, 2, 3, 4, 0, 0, 0, 0], + [9, 10, 11, 12, 0, 0, 0, 0], + [5, 6, 7, 8, 0, 0, 0, 0]])""" + t = """degP = 3""" + print s, t + exec(s) + exec(t) + + print clcoef(P, degP) diff -r 000000000000 -r 0efde00f9229 python/colsplit.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/colsplit.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,40 @@ +#!/usr/bin/python +# Produces two polynomial matrices P1, P2 +# P1 consists of first p1 columns of P +# P2 consists of remaining p2 columns of P + +import pylab as pl +from polyfuncs import polsize + +def colsplit(P, degP, p1, p2): + P = pl.atleast_2d(P) + if pl.atleast_2d(P).size is 0: + P1 = pl.atleast_2d([]) + P2 = pl.atleast_2d([]) + return P1, None, P2, None + + rP, cP = polsize(P, degP) + + if not 0 <= p1 <= cP or not 0 <= p2 <= cP or p1+p2 != cP: + print 'colsplit: Inconsistent numbers of columns' + return None + exit(1) + + P1 = pl.atleast_2d(pl.zeros((rP, (degP+1)*p1))) + P2 = pl.atleast_2d(pl.zeros((rP, (degP+1)*p2))) + + + for i in range(degP+1): + P1[:,i*p1:(i+1)*p1] = P[:,i*cP:i*cP+p1] + P2[:,i*p2:(i+1)*p2] = P[:,i*cP+p1:(i+1)*cP] + + return P1, degP, P2, degP + +if __name__== "__main__": + s = """P = pl.array([[1, 2, 3, 4, 1, 2, 3, 4], + [9, 10, 11, 12, 9, 10, 11, 12], + [5, 6, 7, 8, 5, 6, 7, 8]])""" + print s + exec(s) + print """colsplit(P, 3, 1, 1)""" + print colsplit(P, 3, 1, 1) diff -r 000000000000 -r 0efde00f9229 python/covar_m.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/covar_m.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,65 @@ +#!/usr/bin/python + +import pylab as pl +from scipy.linalg import schur +from dscr import dscr +from scipy import signal + +def dlyap(a, b): + n = len(a) + x = pl.zeros_like(a) + s, u = schur(a) + b = pl.dot(u.T, pl.dot(b,u)) + j = n-1 + while j>=0: + k = j + ## Check for Schur block. + if j==0: + blksiz = 1 + elif s[j, j-1]!=0: + blksiz = 2 + j = j - 1 + else: + blksiz = 1 + Ajj = pl.kron(s[j:k+1,j:k+1], s) - pl.eye(blksiz*n) + rhs = pl.reshape(b[:,j:k+1].T, (blksiz*n, 1)) + if (k < n-1): + rhs2 = pl.dot(s, pl.dot(x[:,k+1:n], s[j:k+1, k+1:n].T)) + rhs = rhs + pl.reshape(rhs2, (blksiz*n, 1)) + v = -pl.solve(Ajj, rhs) + x[:,j] = v.squeeze()[:n] + if(blksiz == 2): + x[:, k] = v[n:blksiz*n].squeeze() + j = j - 1 + + ## Back-transform to original coordinates. + x = pl.dot(u, pl.dot(x, u.T)) + return x + + +def covar_m(H, W): + """ + User defined equivalent function to Matlab covar function + For discrete time domain only + Uses Lyapunov's equation for computation + W: noise intensity (scalar) + """ + a = pl.roots(H.den) + if pl.any(abs(a) > 1): +# print "Warning: System being unstable has infinite covariance" + P = pl.inf + return P + else: + A, B, C, D = H.A, H.B, H.C, H.D + # Sylvester and Lyapunov solver + Q1 = pl.dot(-B, pl.dot(W, B.T)) + Q = dlyap(A, -Q1) + # Q = linmeq(2,A,Q1,[1, 0],1) + # A*X*A' - X + B*W*B' = 0, (2b) + # Discrete time Lyapunov equation; A is general form. Hessenberg-Schur method. + # linmeq(2, A, C, [1,0], 1) + # A*X*A' - X = C, (2b) + # + P = pl.dot(C, pl.dot(Q,C.T)) + pl.dot(D, pl.dot(W,D.T)) + + return P diff -r 000000000000 -r 0efde00f9229 python/desired.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/desired.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,10 @@ +#!/usr/bin/python +# 9.4 +import scipy as sp +def desired(Ts,rise,epsilon): + Nr = rise/Ts + omega = sp.pi/2/Nr + rho = epsilon**(omega/sp.pi) + phi = sp.array([1 -2*rho*sp.cos(omega), rho**2]) + dphi = len(phi)-1 + return phi, dphi diff -r 000000000000 -r 0efde00f9229 python/dscr.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/dscr.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,15 @@ +import scipy as sp +from scipy import linalg +from scipy import signal + +def dscr(H,Ts): + """ Given a Continuous system and a Time step, returns a discretized system.""" + (n,m) = H.B.shape + S = sp.zeros((n+m,n+m)) + S[:n,:n] = H.A + S[:n,n:] = H.B + s = linalg.expm(S*Ts) + f = s[:n,:n] + g = s[:n,n:] + Hd = signal.lti(f, g, H.C, H.D) + return Hd diff -r 000000000000 -r 0efde00f9229 python/fitval.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/fitval.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,11 @@ +#!/usr/bin/python + +# finds the value of a polynomial in powers of z^{-1} +# function Y = filtval(P,z) + +def fitval(P, z): + N = len(P) - 1 + return Y + +Q = polyno(P,'x'); +Y = horner(Q,z)/z^N; diff -r 000000000000 -r 0efde00f9229 python/indep.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/indep.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,43 @@ +#!/usr/bin/python +# Determines the first row that is dependent on the previous rows of S. +# Returns the coefficients of dependence + +import pylab as pl +from makezero import makezero + +def indep(S, gap=1e8): + S = pl.atleast_2d(S) + E = pl.finfo(float).eps + r, c = S.shape + ind = 1 + i = 2 + while ind and i<=r: + shortS = S[:i,:] + sigma = pl.svd(pl.dot(shortS,shortS.T))[1] + if(sigma[-1]/sigma[0] < (E*max(i, c))): + ind = 0 + else: + shsig = sigma.copy() + shsig[:-1] = sigma[1:] + if pl.any(sigma/shsig > gap): + ind = 0 + else: + ind = 1 + i += 1 + if ind: + b = pl.atleast_2d([]) + + else: + A, B = S[i-1,:], S[:i-1,:] + # c = A/B, + # c = (B.T\A.T).T + c = pl.lstsq(B.T, A.T)[0].T +# c = pl.solve(B.T, A.T).T + c = makezero(c, gap) + c = c.squeeze() + b = pl.hstack((-c, 1)) + b = pl.atleast_2d(b) + return b + +if __name__== "__main__": + pass diff -r 000000000000 -r 0efde00f9229 python/l2r.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/l2r.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,19 @@ +#!/usr/bin/python +from transp import transp +from left_prm import left_prm + +def l2r(N, degN, D, degD): + # given Num and Den polynomial matrices in left form, + # not necessarily coprime, finds right coprime factorisation. + N, degN = transp(N, degN) + D, degD = transp(D, degD) + Rnum, Rnumdeg, Rden, Rdendeg = left_prm(N, degN, D, degD) + Rnum, Rnumdeg = transp(Rnum, Rnumdeg) + Rden, Rdendeg = transp(Rden, Rdendeg) + + Rnum = Rnum.squeeze() + Rden = Rden.squeeze() + return Rnum, Rnumdeg, Rden, Rdendeg + + + diff -r 000000000000 -r 0efde00f9229 python/left_prm.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/left_prm.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,92 @@ +#!/usr/bin/python + +import pylab as pl +from rowjoin import rowjoin +from polyfuncs import polsize +from t1calc import t1calc +from seshft import seshft +from colsplit import colsplit +from clcoef import clcoef + +def left_prm(N, degN, D, degD, job=1, gap=None): + """ """ + if gap == None: + gap = 1.0e8 + + F, degF = rowjoin(D, degD, N, degN) + Frows, Fbcols = polsize(F, degF) + Fcols = Fbcols * (degF+1) + T1 = pl.array([]) + pr = pl.array([]) + degT1, T1rows, shft = 0, 0, 0 + S=F + sel = pl.ones((Frows,1)) + T1bcols =1 + abar = pl.arange(Fbcols, Frows) + if T1rows: + T1 = pl.hstack((T1, pl.zeros((T1rows, Frows)))) + if T1.size == 0: + T1rows = 0 + while T1.size == 0 or T1rows < (Frows - Fbcols): + if T1.size == 0: + T1rows = 0 + Srows = Frows*T1bcols + T1, T1rows, sel, pr = t1calc(S, Srows, T1, T1rows, sel, pr, Frows, Fbcols, abar, gap) + T1rows, T1cols = pl.atleast_2d(T1).shape + if T1.size == 0: + T1rows = 0 + if T1rows < (Frows - Fbcols): + if T1.size != 0: + T1 = pl.hstack((T1, pl.zeros((T1rows,Frows)).squeeze())) + T1bcols += 1 + degT1 += 1 + shft += Fbcols + S = seshft(S,F,shft); + sel = pl.hstack((sel,sel[Srows-Frows:Srows+1])) + rowvec = pl.arange((T1bcols-1)*Frows+(Fbcols), (T1bcols * Frows)) + abar = pl.hstack((abar, rowvec)) + + B, degB, A, degA = colsplit(T1, degT1, Fbcols, Frows-Fbcols) + B, degB = clcoef(B, degB) + B = -B; + A, degA = clcoef(A, degA) + if job == 1: + return B, degB, A, degA + + if job == 2: + S = S[sel!=0,:] + redSrows, Scols = S.shape + C = pl.hstack(pl.eye(Fbcols), pl.zeros(Fbcols,Scols-Fbcols)) + T2 = pl.dot(C, pl.inv(S)) + T2 = makezero(T2,gap); + T2 = move(T2, pl.find(sel), Srows); + X, degX, Y, degY = colsplit(T2, degT1, Fbcols, Frows - Fbcols); + X, degX = clcoef(X, degX); + Y, degY = clcoef(Y, degY); + elif job == 3: + Y = S + degY = sel + X = degT1 + degX = Fbcols + + else: + print 'Message from left_prm:no legal job number specified' + exit() + + return B, degB, A, degA, Y, degY, X, degX + + +if __name__== "__main__": + D = pl.array([[1, 0, 0, 0, 0, 0], + [0, 1, 0, 1, 0, 0], + [0, 0, 1, 1, 1, 0]]) + + N = pl.array([[1, 0, 0], + [0, 1, 0], + [0, 0, 1]]) + + dD, dN = 1, 0 + + result = left_prm(N, dN, D, dD) + print result + pass diff -r 000000000000 -r 0efde00f9229 python/makezero.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/makezero.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,28 @@ +#!/usr/bin/python + +import pylab as pl + +def makezero(B, gap=None): + B = pl.array(B) # B is a vector! + B = B.squeeze() + if len(B.shape)!=1: + print "B must be a vector!" + exit() + if gap == None: + gap = 1e8 + temp = B[pl.find(B)] + temp = -pl.sort(-abs(temp)) + ratio = temp[:-1]/pl.float32(temp[1:]) + if pl.find(ratio>gap).size != 0: + min_ind = min(pl.find(ratio>gap)) + our_eps = temp[min_ind+1] + zeroind = pl.find(abs(B)<=our_eps) + B[zeroind] = 0 + return pl.atleast_2d(B) + +if __name__== "__main__": + B = pl.array([2,3,.000000000004,5,.00006]) + print B + B = makezero(B) + print B + diff -r 000000000000 -r 0efde00f9229 python/margin.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/margin.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,19 @@ +import pylab as pl +def p_margin(H, type='d', Ts=1): + // |h(e^(i*w*dt))|=1 <-- h(e^(i*w*dt))*h(e^(-i*w*dt)) + z=poly(0,varn(h.den)); + den = H.den + sm=simp_mode();simp_mode(%f);hh=h*horner(h,1/z)-1;simp_mode(sm) + //find the numerator roots + z=roots(hh.num); + z(abs(abs(z)-1)>eps)=[];// retain only roots with modulus equal to 1 + w=log(z)/(%i*dt); + ws=real(w(abs(imag(w))0)); //frequency points with unitary modulus + if ws==[] then phm=%inf,fr=[],return,end + f=horner(h,exp(%i*ws*dt)); + end + return phm, fr + +def g_margin(H, type='d', Ts=1): + + return gm, fr diff -r 000000000000 -r 0efde00f9229 python/move.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/move.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,27 @@ +#!/usr/bin/python + +import pylab as pl + +def move(b, nonred, maxi): + """ """ + br, bc = pl.atleast_2d(b).shape + b = pl.hstack((b, pl.zeros((br, len(nonred)-bc)))) + maxi = max(max(nonred)+1, maxi) + result = pl.zeros((br, maxi)) + nonred = pl.array(nonred) + result[:,nonred.T]=b + return result + +if __name__== "__main__": + + s = """b = pl.array([[1, 2], + [9, 10]])""" + t = """nonred = pl.array([7, 2])""" + u = """maxi = 10""" + print s, t, u + exec(s) + exec(t) + exec(u) + + print """move(b, nonred, maxi)""" + print move(b, nonred, maxi) diff -r 000000000000 -r 0efde00f9229 python/myc2d.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/myc2d.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,18 @@ +#!/usr/bin/python + +from dscr import dscr +import scipy as sp + +def myc2d(G,Ts): + """ Produces numerator and denominator of discrete transfer + function in powers of z^{-1} + G is continuous transfer function; time delays are not allowed + Ts is the sampling time, all in consistent time units. """ + H = dscr(G,Ts) + num1, den1 = H.num, H.den + A = den1[::-1] + num2 = num1[::-1] + nonzero = sp.find(num1) + B = num2(nonzero) + k = len(den1) - len(num1) + return B,A,k diff -r 000000000000 -r 0efde00f9229 python/polyfuncs.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/polyfuncs.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,164 @@ +#!/usr/bin/python +# Contains all the polynomial functions + +import pylab as pl + +def polsize(Q, degQ): + """ Determines dimensions of a polynomial matrix. """ + rQ, cQ = pl.atleast_2d(Q).shape + cQ = cQ/float(degQ+1) + if abs(round(cQ)-cQ) > 1e-6: + print "Degree of input inconsistent with number of columns" + return + else: + cQ = int(round(cQ)) + return rQ, cQ + +def polmul(A, degA, B, degB): + A = pl.atleast_2d(A) + B = pl.atleast_2d(B) + rA, cA = polsize(A, degA) + rB, cB = polsize(B, degB) + if cA != rB: + print "polmul: Inconsistent dimensions of input matrices" + return + degC = degA + degB + C = [] + for k in range(0, degC+1): + mi = 0 + if k-degB > mi: + mi = k-degB + ma = degA + if k < ma: + ma = k + Ck = pl.zeros((rA,cB)) + for i in range(mi, ma+1): + Ck = Ck + pl.dot(A[..., i*cA:(i+1)*cA], B[..., (k-i)*cB:(k-i+1)*cB]) + Ck = pl.squeeze(Ck) + C = pl.hstack((C, Ck)) + return C, degC + +def poladd(A, degA, B, degB): + A = pl.atleast_2d(A) + B = pl.atleast_2d(B) + rA, cA = polsize(A, degA) + rB, cB = polsize(B, degB) + if cA != rB: + print "polmul: Inconsistent dimensions of input matrices" + return + degC = max(degA, degB) + + if degC >= degA: + A = pl.hstack((A, pl.zeros((rA,(degC-degA)*cA)))) + + if degC >= degB: + B = pl.hstack((B, pl.zeros((rB,(degC-degB)*cB)))) + + C = A + B + return C, degC + +def polsplit2(fac, a=1-1e-5): + fac = pl.atleast_1d(fac) + if a>1: + print "good polynomial also is unstable" + return + roots = pl.roots(fac) + + # extract good and bad roots + badindex = pl.find(pl.absolute(roots)>=a-1.0e-5) + badpoly = pl.poly(roots[badindex]) + goodindex = pl.find(pl.absolute(roots)1: + print "good polynomial also is unstable" + return + roots = pl.roots(fac) + + # extract good and bad roots + badindex = pl.find((pl.absolute(roots)>=a-1.0e-5) + (pl.real(roots)<-0.05)) + badpoly = pl.poly(roots[badindex]) + goodindex = pl.find((pl.absolute(roots)=-0.05)) + goodpoly = pl.poly(roots[goodindex]) + # scale by equating the largest terms + index = pl.absolute(fac).argmax() + goodbad = pl.convolve(goodpoly, badpoly) + factor = fac[index]/goodbad[index] + goodpoly = goodpoly * factor + badpoly = pl.atleast_1d(badpoly) + goodpoly = pl.atleast_1d(goodpoly) + return goodpoly, badpoly + +def putin(A, degA, B, degB, i, j): + from clcoef import clcoef + A = pl.atleast_2d(A) + B = pl.atleast_2d(B) + rA, cA = polsize(A,degA) + if degB > degA: + A = pl.hstack((A, pl.zeros((rA,(degB-degA)*cA)))) + degA = degB + + for k in range(degB+1): + A[i,(k*cA)+j] = B[0,k] + + if degA > degB: + for k in range(degB+1,degA+1): + A[i,(k*cA)+j] = 0 + A, degA = clcoef(A,degA) + return A, degA + + +def ext(A, degA, k, l): + from clcoef import clcoef + rA, cA = polsize(A, degA) + degB = degA + B = pl.zeros((1, degB+1)) + for m in range(degB+1): + B[0, m] = A[k, (m*cA)+l] + B,degB = clcoef(B, degB) + return B, degB + + +def transp(Q, degQ): + """ Function to transpose a polynomial matrix. """ + rQ, cQ = polsize(Q, degQ) + rP = cQ + cP = rQ + degP = degQ + P = pl.zeros((rP, (degP+1)*cP)) + for i in range(degP+1): + P[:, i*cP:(i+1)*cP] = Q[:, i*cQ:(i+1)*cQ].T + + return P, degP + + +if __name__== "__main__": + + # print "Test for polsize" + # print polsize([1, 2, 1],4) + + # print "Test for polmul" + # C = pl.array([[1, 0, 0.5, 2], [0, 1, -4.71, 2.8]]) + # A = pl.array([0.5, 3.5]) + # print polmul(A, 0, C, 1) + + # print "Test for polsplit3" + # print polsplit3([1, -0.37]) + + print "Test for putin" + A = pl.array([0,0]) + B = pl.array([0.44, -1.6, 1.6, -0.44]) + print putin(A, 0, B, 3, 0, 0) + + + pass diff -r 000000000000 -r 0efde00f9229 python/pp_basic.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/pp_basic.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,47 @@ +#!/usr/bin/python +# 9.8 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import scipy as sp +from polyfuncs import polmul +from polyfuncs import polsplit2 +from zpowk import zpowk +from xdync import xdync + +def pp_basic(B, A, k, phi): + # Setting up and solving Aryabhatta identity + Ag, Ab = polsplit2(A) + dAb = len(Ab) - 1 + Bg, Bb = polsplit2(B) + dBb = len(Bb) - 1 + + zk, dzk = zpowk(k) + + N, dN = polmul(Bb, dBb, zk, dzk) + dphi = len(phi)-1 + + S1, dS1, R1, dR1 = xdync(N, dN, Ab, dAb, phi, dphi)[:4] + + # Determination of control law + R1 = sp.squeeze(R1) + S1 = sp.squeeze(S1) + Rc = sp.convolve(Bg, R1) + Sc = sp.convolve(Ag, S1) + Tc = Ag + gamm = phi.sum()/Bb.sum() + return Rc, Sc, Tc, gamm + +if __name__ == "__main__": + Ts = 1 + B = 0.63 + A = sp.array([1, -0.37]) + k = 1 + if k<=0: + k = 1 + + zk, dzk = zpowk(k) + + phi = sp.array([1, -0.5]) + + print pp_basic(B, A, k, phi) diff -r 000000000000 -r 0efde00f9229 python/pp_im.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/pp_im.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,52 @@ +#!/usr/bin/python +# 9.8 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import scipy as sp +from polyfuncs import polmul +from polyfuncs import polsplit3 +from zpowk import zpowk +from xdync import xdync + +def pp_im(B, A, k, phi, Delta): + Delta = sp.atleast_1d(Delta) + # Setting up and solving Aryabhatta identity + Ag, Ab = polsplit3(A) + dAb = len(Ab) - 1 + Bg, Bb = polsplit3(B) + dBb = len(Bb) - 1 + + zk, dzk = zpowk(k) + + N, dN = polmul(Bb, dBb, zk, dzk) + dDelta = len(Delta)-1 + D, dD = polmul(Ab, dAb, Delta, dDelta) + dphi = len(phi)-1 + + S1, dS1, R1, dR1 = xdync(N, dN, D, dD, phi, dphi)[:4] + + # Determination of control law + R1 = sp.squeeze(R1) + S1 = sp.squeeze(S1) + Rc = sp.convolve(Bg, sp.convolve(R1, Delta)) + Sc = sp.convolve(Ag, S1) + Tc = Ag + gamm = phi.sum()/Bb.sum() + return Rc, Sc, Tc, gamm + +if __name__ == "__main__": + Ts = 1 + B = 0.63 + A = sp.array([1, -0.37]) + k = int(raw_input('Enter the delay as an integer: ')) + + if k<=0: + k = 1 + + zk, dzk = zpowk(k) + + phi = sp.array([1, -0.5]) + delta = 1 + + print pp_im(B, A, k, phi, delta) diff -r 000000000000 -r 0efde00f9229 python/pp_im2.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/pp_im2.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,38 @@ +#!/usr/bin/python +# 9.8 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import scipy as sp +from polyfuncs import polmul +from polyfuncs import polsplit3 +from zpowk import zpowk +from xdync import xdync + +def pp_im2(B, A, k, phi, Delta, a=1): + Delta = sp.atleast_1d(Delta) + # Setting up and solving Aryabhatta identity + Ag, Ab = polsplit3(A, a) + dAb = len(Ab) - 1 + Bg, Bb = polsplit3(B, a) + dBb = len(Bb) - 1 + + zk, dzk = zpowk(k) + + N, dN = polmul(Bb, dBb, zk, dzk) + dDelta = len(Delta)-1 + D, dD = polmul(Ab, dAb, Delta, dDelta) + dphi = len(phi)-1 + + S1, dS1, R1, dR1 = xdync(N, dN, D, dD, phi, dphi)[:4] + + # Determination of control law + R1 = sp.squeeze(R1) + S1 = sp.squeeze(S1) + Rc = sp.convolve(Bg, sp.convolve(R1, Delta)) + Sc = sp.convolve(Ag, S1) + Tc = Ag + gamm = phi.sum()/Bb.sum() + # Total characteristic polynomial + phit = sp.convolve(phi, sp.convolve(Ag, Bg) + return Rc, Sc, Tc, gamm, phit diff -r 000000000000 -r 0efde00f9229 python/pp_pid.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/pp_pid.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,43 @@ +#!/usr/bin/python +# 9.8 +import os, sys +sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"] + +import scipy as sp +from polyfuncs import polmul +from zpowk import zpowk +from xdync import xdync + +def pp_pid(B, A, k, phi, Delta): + Delta = sp.atleast_1d(Delta) + B = sp.atleast_1d(B) + # Setting up and solving Aryabhatta identity + dA = len(A) - 1 + dB = len(B) - 1 + + zk, dzk = zpowk(k) + + N, dN = polmul(B, dB, zk, dzk) + dDelta = len(Delta)-1 + D, dD = polmul(A, dA, Delta, dDelta) + dphi = len(phi)-1 + + Sc, dSc, R, dR = xdync(N, dN, D, dD, phi, dphi)[:4] + + R = sp.squeeze(R) + Rc = sp.convolve(R, Delta) + + return Rc, Sc + +if __name__ == "__main__": + Ts = 1 + B = 0.63 + A = sp.array([1, -0.37]) + k = 1 + + zk, dzk = zpowk(k) + + phi = sp.array([1, -0.5]) + delta = 1 + + print pp_pid(B, A, k, phi, delta) diff -r 000000000000 -r 0efde00f9229 python/rowjoin.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/rowjoin.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,32 @@ +#!/usr/bin/python +# Superposes two polynomial matrices. + +import pylab as pl +from polyfuncs import polsize + +def rowjoin(P1, degP1, P2, degP2): + + rP1, cP1 = polsize(P1, degP1) + rP2, cP2 = polsize(P2, degP2) + if not cP1 == cP2: + print 'rowjoin: Inconsistent numbers of columns' + exit() + + rP = rP1 + rP2 + cP = cP1 + degP = max(degP1, degP2) + + P = pl.atleast_2d(pl.zeros((rP,(degP+1)*cP))) + P[0:rP1,0:(degP1+1)*cP1] = P1 + P[rP1:rP,0:(degP2+1)*cP2] = P2 + + return P, degP + +if __name__== "__main__": + s = """P = pl.array([[1, 2, 3, 4, 1, 2, 3, 4], + [9, 10, 11, 12, 9, 10, 11, 12], + [5, 6, 7, 8, 5, 6, 7, 8]])""" + print s + exec(s) + print """rowjoin(P, 3, P, 3)""" + print rowjoin(P, 3, P, 3) diff -r 000000000000 -r 0efde00f9229 python/seshft.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/seshft.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,38 @@ +#!/usr/bin/python +# given A and B arrays, returns C = [<-A-> 0 +# 0 <-B->] with B shifted east by N cols + +import pylab as pl +from scipy import signal + +def seshft(A, B, N): + Arows, Acols = pl.atleast_2d(A).shape + Brows, Bcols = pl.atleast_2d(B).shape + + if N>=0: + n = Arows+Brows + m = max(Acols, Bcols+N) + C = pl.zeros((n, m)) + C[:Arows, :Acols] = A + C[Arows:, N:N+Bcols] = B + + else: + N = -N + n = Arows+Brows + m = max(Acols+N, Bcols) + C = pl.zeros((n, m)) + C[:Arows, N:N+Arows] = A + C[Arows:, :Bcols] = B + + return C + +if __name__== "__main__": + a = "A = pl.eye(3)" + b = "B = pl.eye(10)" + n = "N = 0" + exec(a) + exec(b) + exec(n) + print a, b, n + C = seshft(A, B, N) + print C diff -r 000000000000 -r 0efde00f9229 python/t1calc.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/python/t1calc.py Fri May 27 14:24:59 2011 +0530 @@ -0,0 +1,40 @@ +#!/usr/bin/python + +import pylab as pl +from indep import indep +from move import move + +def t1calc(S, Srows, T1, T1rows, sel, pr, Frows, Fbcols, abar, gap=None): + b = pl.array([1]) + sel = sel.squeeze() + while (T1rows < Frows - Fbcols) and pl.any(sel==1) and b.size!=0: + b = indep(S[pl.logical_and(sel, sel),:], gap) + if b.size != 0: + b = move(b, pl.find(sel), Srows) + b = b.squeeze() + j = b.shape[-1]-1 + while not (b[j] and pl.any(abar==j)): + j -= 1 + if j == -1: + print 'Message from t1calc, called from left_prm' + print 'Denominator is noninvertible' + exit() + aa = pl.remainder(pr, Frows) == pl.remainder(j, Frows) + bb = j