Initial commit.
--- /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
--- /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()
--- /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()
--- /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
--- /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
--- /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???
--- /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()
--- /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]
--- /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)
--- /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)
--- /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)
--- /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)
--- /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)
--- /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.
+# """
+
+
--- /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
--- /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)
--- /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)
--- /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()
--- /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()
--- /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()
--- /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()
--- /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)
--- /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)
--- /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
--- /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
--- /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
--- /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
--- /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
--- /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
--- /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
--- /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
--- /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
--- /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
--- /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
--- /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)
--- /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)
--- /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
--- /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
--- /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)
--- /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
--- /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)
--- /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"
--- /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()
+
--- /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
--- /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
+
--- /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)
+
--- /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
+
--- /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
+
--- /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)
--- /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 "
+
--- /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
--- /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 "
+
--- /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
--- /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
--- /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
--- /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
--- /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
--- /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
+
--- /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
+
--- /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
+
--- /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
+
+
--- /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
--- /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
+
--- /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
--- /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);
--- /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
--- /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
--- /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
--- /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
--- /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
--- /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
--- /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
--- /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
+
+
--- /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
--- /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
+
--- /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)
+
--- /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
--- /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
--- /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
--- /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)
--- /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)
--- /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
--- /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
--- /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
--- /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;
--- /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
--- /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
+
+
+
--- /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
--- /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
+
--- /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))<eps&real(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
--- /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)
--- /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
--- /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)<a-1.0e-5)
+ 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 polsplit3(fac, a=1):
+ 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) + (pl.real(roots)<-0.05))
+ badpoly = pl.poly(roots[badindex])
+ goodindex = pl.find((pl.absolute(roots)<a-1.0e-5) * (pl.real(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
--- /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)
--- /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)
--- /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
--- /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)
--- /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)
--- /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
--- /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<pr
+ if not pl.any(bb * aa):
+ if T1.size is 0:
+ T1 = b.copy()
+ else:
+ T1 = pl.row_stack((T1, b))
+ T1rows += 1
+ if pr.size is 0:
+ pr = pl.array([j])
+ else:
+ pr = pl.vstack((pr, j))
+ while j < Srows:
+ sel[j] = 0
+ j += Frows
+ return T1, T1rows, sel, pr
+
+if __name__== "__main__":
+ pass
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/python/tfvar.py Fri May 27 14:24:59 2011 +0530
@@ -0,0 +1,21 @@
+#!/usr/bin/python
+# 11.7
+
+# N and D polynomials in z^{-1} form; discrete case
+import scipy as sp
+from l2r import l2r
+from scipy import signal
+from covar_m import covar_m
+
+def tfvar(N,dN,D,dD):
+ N,dN,D,dD = l2r(N,dN,D,dD)
+ N = N/D[0]
+ D = D/D[0]
+ LN, LD = len(N), len(D)
+ D1 = D
+ if LD<LN:
+ D1 = sp.hstack((D, sp.zeros(LN-LD)))
+ dD1 = dD+LN-LD
+ H = signal.lti(N, D1)
+ yvar = covar_m(H,1)
+ return N,dN,D,dD,yvar
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/python/tran.py Fri May 27 14:24:59 2011 +0530
@@ -0,0 +1,16 @@
+#!/usr/bin/python
+
+from polyfuncs import polsize
+import scipy as sp
+
+def transp(Q, degQ):
+ """ Function to transpose a polynomial matrix. """
+ rQ, cQ = polsize(Q, degQ)
+ rP = cQ
+ cP = rQ
+ degP = degQ
+ P = sp.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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/python/transp.py Fri May 27 14:24:59 2011 +0530
@@ -0,0 +1,16 @@
+#!/usr/bin/python
+
+from polyfuncs import polsize
+import scipy as sp
+
+def transp(Q, degQ):
+ """ Function to transpose a polynomial matrix. """
+ rQ, cQ = polsize(Q, degQ)
+ rP = cQ
+ cP = rQ
+ degP = degQ
+ P = sp.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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/python/xdync.py Fri May 27 14:24:59 2011 +0530
@@ -0,0 +1,53 @@
+#!/usr/bin/python
+
+import pylab as pl
+from rowjoin import rowjoin
+from polyfuncs import polsize
+from cindep import cindep
+from seshft import seshft
+from colsplit import colsplit
+from left_prm import left_prm
+from move import move
+from clcoef import clcoef
+
+def xdync(N, degN, D, degD, C, degC, gap=1.0e8):
+ """ """
+ C = pl.atleast_2d(C)
+ F, degF = rowjoin(D, degD, N, degN)
+ Frows, Fbcols = polsize(F, degF)
+ B, degB, A, degA, S, sel, degT1, Fbcols = left_prm(N, degN, D, degD, 3, gap)
+ Crows, Ccols = C.shape
+ Srows, Scols = pl.atleast_2d(S).shape
+ S = S[sel!=0,:]
+ T2 = pl.array([])
+
+ for i in range(Crows):
+ Saug = seshft(S, C[i,:], 0)
+ b = cindep(Saug)
+ b = move(b, pl.find(sel), Srows)
+ if T2.size != 0:
+ T2 = pl.vstack((T2, b))
+ else:
+ T2 = b.copy()
+
+ X, degX, Y, degY = colsplit(T2, degT1, Fbcols, Frows-Fbcols)
+ X, degX = clcoef(X, degX)
+ Y, degY = clcoef(Y, degY)
+
+ return Y, degY, X, degX, B, degB, A, degA
+
+
+if __name__== "__main__":
+ N = pl.array([[0, 4, 0, 1],
+ [-1, 8, 0, 3]])
+ dN = 1
+ D = pl.array([[0, 0, 1, 4, 0, 1],
+ [0, 0, -1, 0, 0, 0]])
+ dD = 2
+ C = pl.array([[1, 0, 1, 1],
+ [0, 2, 0, 1]])
+ dC = 1
+
+ print xdync(N,dN,D,dD,C,dC)
+
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/python/zpowk.py Fri May 27 14:24:59 2011 +0530
@@ -0,0 +1,8 @@
+#!/usr/bin/python
+
+import scipy as sp
+def zpowk(k):
+ zk = sp.zeros((1, k+1))
+ zk[0,k] = 1
+ dzk = k
+ return zk, dzk
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/specs/abex.py Fri May 27 14:24:59 2011 +0530
@@ -0,0 +1,27 @@
+#!/usr/bin/python
+# 7.8
+
+import os, sys
+sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"]
+
+import pylab as pl
+from xdync import xdync
+
+N = pl.convolve([0, 1], [1, 1])
+D = pl.convolve([1, -4], [1, -1])
+dN = 2
+dD = 2
+C = pl.array([1, -1, 0.5])
+dC = 2
+
+Y, dY, X, dX, B, dB, A, dA = xdync(N, dN, D, dD, C, dC)
+
+print "Y =", Y
+print "dY =", dY
+print "X =", X
+print "dX =", dX
+print "B =", B
+print "dB =", dB
+print "A =", A
+print "dA =", dA
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/specs/ant_lead.py Fri May 27 14:24:59 2011 +0530
@@ -0,0 +1,38 @@
+#!/usr/bin/python
+# 7.6
+import os, sys
+sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"]
+
+# Continuous time antenna model
+
+import pylab as pl
+from scipy import signal
+from dscr import dscr
+
+a = 0.1
+F = pl.array([[0, 1],[0, -a]])
+g = pl.array([[0],[a]])
+c = pl.array([1, 0])
+d = 0
+Ga = signal.lti(F,g,c,d)
+num,den = Ga.num, Ga.den
+Ts = 0.2
+G = dscr(Ga,Ts)
+
+#lead controller
+beta1 = 0.8
+N = pl.array([1, -0.9802])*(1-beta1)/(1-0.9802)
+Rc = pl.array([1, -beta1])
+
+# // simulation parameters using g_s_cl2.cos
+# gamm = 1; Sc = 1; Tc = 1; C = 0; D = 1;
+# st = 1; st1 = 0;
+# t_init = 0; t_final = 20;
+
+# // u1: -4 to 11
+# // y1: 0 to 1.4
+
+# [Tcp1,Tcp2] = cosfil_ip(Tc,1); // Tc/1
+# [Np,Rcp] = cosfil_ip(N,Rc); // N/Rc
+# [Scp1,Scp2] = cosfil_ip(Sc,1); // Sc/1
+# [Cp,Dp] = cosfil_ip(C,D); // C/D
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/specs/data01.py Fri May 27 14:24:59 2011 +0530
@@ -0,0 +1,26 @@
+#!/usr/bin/python
+# 7.9
+import os, sys
+sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"]
+
+import pylab as pl
+from left_prm import left_prm
+
+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
+
+B, dB, A, dA = left_prm(N, dN, D, dD)
+
+print "A =", A
+print "dA =", dA
+
+print "B =", B
+print "dB =", dB
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/specs/data05.py Fri May 27 14:24:59 2011 +0530
@@ -0,0 +1,32 @@
+#!/usr/bin/python
+import os, sys
+sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"]
+
+
+import pylab as pl
+from xdync import xdync
+
+N = pl.array([[0, 4, 0, 1],
+ [-1, 8, 0, 3]])
+dN = 1
+
+D = pl.array([[0, 0, 1, 4, 0, 1],
+ [0, 0, -1, 0, 0, 0]])
+dD = 2
+C = pl.array([[1, 0, 1, 1],
+ [0, 2, 0, 1]])
+dC = 1
+
+Y, dY, X, dX, B, dB, A, dA = xdync(N, dN, D, dD, C, dC)
+
+print "Y =", Y
+print "dY =", dY
+print "X =", X
+print "dX =", dX
+print "B =", B
+print "dB =", dB
+print "A =", A
+print "dA =", dA
+
+
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/specs/nyquist_ex1.py Fri May 27 14:24:59 2011 +0530
@@ -0,0 +1,42 @@
+#!/usr/bin/python
+# 7.2
+import os, sys
+sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"]
+
+import pylab as pl
+from scipy import signal
+from bode import freqresp, phasemag
+
+
+# def freqresp(H, omega):
+# num = pl.poly1d(H.num)
+# den = pl.poly1d(H.den)
+
+# fresp = map(lambda w: num(w*1j) / den(w*1j), omega)
+# fresp = pl.array(fresp)
+# mag = pl.absolute(fresp)
+# phase = pl.angle(fresp)
+
+# return mag, phase, omega
+
+def nyquist(sys, omega=None):
+ if (omega == None):
+ omega = pl.logspace(-2, 2, 5000)
+
+ omega, fresp = freqresp(sys, omega, 'd')
+
+ x = pl.real(fresp)
+ y = pl.imag(fresp)
+
+ pl.plot(x, y, '-');
+ pl.plot(x, -y, '--');
+
+ pl.plot([-1], [0], '+k')
+
+ pl.show()
+
+num = [1]
+den = [1, -1, 0]
+H = signal.lti(num, den)
+
+nyquist(H)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/specs/rlocus_ex1.py Fri May 27 14:24:59 2011 +0530
@@ -0,0 +1,62 @@
+#!/usr/bin/python
+# 7.1
+import os, sys
+sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"]
+
+import pylab as pl
+from scipy import signal
+
+def rlocus(H, k_vec=None):
+ """Plots the Root Locus given
+ an lti system and
+ a vector of gain values"""
+ if k_vec is None:
+ k_vec = linspace(0, 100, 10000)
+
+ num = pl.poly1d(H.num)
+ den = pl.poly1d(H.den)
+
+ roots = []
+ for k in k_vec:
+ p = den + k*num
+ roots.append(pl.roots(p))
+ roots = pl.array(roots)
+
+ # Tried to do it with array slicing. Looks BAD
+ # n = H.num
+ # d = H.den
+ # l = d.shape[0]-n.shape[0]
+ # n_vec = n[pl.newaxis,:]*k_vec[:,pl.newaxis]
+ # A = pl.zeros((k_vec.shape[0], d.shape[0]))
+ # A = A+d
+ # A[:,l:] = n_vec
+ # r = map(pl.roots, A)
+ # r = pl.array(r)
+
+ x = pl.linspace(0, 2*pl.pi, 100)
+ pl.plot(pl.cos(x), pl.sin(x), 'b.', markersize=2)
+
+ for i, k_root in enumerate(roots[0]):
+ pl.plot(pl.real(roots[:,i]), pl.imag(roots[:,i]))
+
+ xmin, xmax = pl.xlim()
+ x = pl.linspace(int(xmin), int(xmax))
+ ymin, ymax = pl.ylim()
+ y = pl.linspace(int(ymin), int(ymax))
+
+ pl.plot(x, pl.zeros_like(x), 'b.', pl.zeros_like(y), y, 'b.', markersize=2)
+
+ pl.xlim(xmin, xmax)
+ pl.ylim(ymin, ymax)
+
+ pl.xlabel('Real part')
+ pl.ylabel('Imaginary part')
+
+ pl.show()
+
+n = [1]
+d = [1, -1, 0]
+k_vec = pl.linspace(0, 20, 10000)
+H = signal.lti(n, d)
+
+rlocus(H, k_vec)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/ss/ex_comp.py Fri May 27 14:24:59 2011 +0530
@@ -0,0 +1,30 @@
+#!/usr/bin/python
+# 14.2
+import os, sys
+sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"]
+
+import scipy as sp
+from scipy import linalg
+from pend import cont_mat, delta_a, acker
+from polyfuncs import polmul
+
+A = sp.array([[1, 2], [0, 3]])
+c = sp.array([1, 0])
+p = sp.roots([1, -0.5, 0.5])
+b = sp.array([[0], [1]])
+K = acker(A, b, p)
+
+p1 = 0.1 + 0.1j
+p2 = 0.1 - 0.1j
+phi = sp.real(sp.convolve([1, -p1], [1, -p2]))
+Obs = sp.vstack((c, sp.dot(c, A)))
+alphae = sp.dot(A, A)-0.2*A+0.02*sp.eye(2,2)
+Lp = sp.dot(alphae, sp.dot(linalg.inv(Obs), sp.array([[0],[1]])))
+
+C = sp.array([[1, 0, 0.5, 2], [0, 1, -4.71, 2.8]])
+dC = 1
+
+[HD,dHD] = polmul(K,0,C,dC)
+[HD,dHD] = polmul(HD,dHD,Lp,0)
+
+print HD, dHD
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/ss/kalrun.py Fri May 27 14:24:59 2011 +0530
@@ -0,0 +1,37 @@
+#!/usr/bin/python
+# 14.4 & 14.3
+
+import pylab as pl
+
+def kal_ex(x, xline, M):
+ y = x + pl.rand()
+ Q = 0.0
+ R = 1.0
+ xhat_ = xline
+ P_ = M + Q
+ K = P_/(P_+R)
+ P = (1-K)*P_
+ xhat = xhat_ + K*(y-xhat_)
+ return xhat, P, y
+
+x = 5
+xhat = 2
+P = 1
+xvec = x
+xhat_vec = xhat
+Pvec = P
+yvec = x
+for i in range(200):
+ xline = xhat
+ M = P
+ xhat, P, y = kal_ex(x, xline, M)
+ xvec = pl.vstack((xvec, x))
+ xhat_vec = pl.vstack((xhat_vec, xhat))
+ Pvec = pl.vstack((Pvec, P))
+ yvec = pl.vstack((yvec, y))
+
+n = range(201)
+pl.plot(Pvec)
+pl.plot(n, xhat_vec, n, yvec, n, xvec)
+pl.xlabel('n')
+pl.show()
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/ss/pend.py Fri May 27 14:24:59 2011 +0530
@@ -0,0 +1,76 @@
+#!/usr/bin/python
+# 14.1
+# Uses A, B from 2.1 (pend_model.py)
+
+import os, sys
+sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"]
+
+import scipy as sp
+from scipy import signal
+from scipy import linalg
+from dscr import dscr
+from pend_model import pend_model
+
+def pol2cart(theta, r):
+ x = r*sp.cos(theta)
+ y = r*sp.sin(theta)
+ return x, y
+
+def cont_mat(a, b):
+ n = len(a)
+ con = b
+ for i in range(1, n):
+ prod = sp.identity(n)
+ for j in range(i):
+ prod = sp.dot(prod, a)
+ prod = sp.dot(prod, b)
+ con = sp.hstack((con, prod))
+ return con
+
+def delta_a(a, P):
+ n = len(a)
+ delta_a = sp.identity(n)
+ for pole in P:
+ delta_a = sp.dot(delta_a, (a - pole*sp.identity(n)))
+ return delta_a
+
+def acker(a, b, P):
+ """ Crude implementation of ackermann's formula """
+ con = cont_mat(a, b)
+ con_inv = linalg.inv(con)
+ del_a = delta_a(a, P)
+ e_n = sp.zeros(len(a))
+ e_n[-1] = 1
+
+ K = sp.dot(sp.dot(e_n, con_inv), del_a)
+ return K
+
+if __name__ == "__main__":
+ A, B = pend_model()
+ # Changed form of C, D from the scilab one to make it work with signal.lti
+ C = sp.eye(1, 4)
+ D = 0
+ Ts = 0.01
+ G = signal.lti(A, B, C, D)
+ H = dscr(G, Ts)
+ a, b, c, d = H.A, H.B, H.C, H.D
+ rise = 5
+ epsilon = 0.1
+ N = rise/Ts
+ omega = sp.pi/2/N
+ r = epsilon**(omega/sp.pi)
+ r1 = r
+ r2 = 0.9*r
+
+ x1, y1 = pol2cart(omega, r1)
+ x2, y2 = pol2cart(omega, r2)
+
+ p1 = x1+y1*1j
+ p2 = x2+y2*1j
+ p3 = x1-y1*1j
+ p4 = x2-y2*1j
+ P = [[p1], [p2], [p3], [p4]]
+
+ K = acker(a,b,P)
+
+ print K
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/ss/pend_model.py Fri May 27 14:24:59 2011 +0530
@@ -0,0 +1,43 @@
+#!/usr/bin/python
+## 2.1
+
+from numpy import zeros
+
+def pend_model():
+ Km = 0.00767
+ Kg = 3.7
+ Rm = 2.6
+ r = 0.00635
+ M = 0.522
+ m = 0.231
+ g = 9.81
+ L = 0.305
+ J = 0
+
+ D1 = (J+m*L**2)*(M+m)-m**2*L**2
+ alpha = m*g*L*(M+m)/D1
+ beta1 = m*L/D1
+ gamma1 = m**2*g*L**2/D1
+ delta = (J+m*L**2)/D1
+
+ alpha1 = Km*Kg/Rm/r;
+ alpha2 = Km**2*Kg**2/Rm/r**2
+
+ A = zeros([4,4])
+ A[0,2] = 1
+ A[1,3] = 1
+ A[2,1] = -gamma1
+ A[2,2] = -alpha2*delta
+ A[3,1] = alpha
+ A[3,2] = alpha2*beta1
+
+ B = zeros([4,1])
+ B[2] = alpha1*delta
+ B[3] = -alpha1*beta1
+
+ return A, B
+
+if __name__ == "__main__":
+ A, B = pend_model()
+ print "A", A
+ print "B", B
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/system/conv2.py Fri May 27 14:24:59 2011 +0530
@@ -0,0 +1,11 @@
+#!/usr/bin/python
+# 3.2
+
+from scipy import signal
+import scipy as sp
+
+h = sp.array([1,2,3])
+u = sp.array([4,5,6])
+y = signal.convolve(u,h)
+print "y =",
+print y
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/system/energy.py Fri May 27 14:24:59 2011 +0530
@@ -0,0 +1,10 @@
+#!/usr/bin/python
+
+from scipy import signal
+from scipy import linalg
+import scipy as sp
+
+u = sp.array([4,5,6])
+Eu = linalg.norm(u)**2
+ruu = signal.correlate(u,u)
+Eu = ruu[len(ruu)/2]