add deconvolve test; refactor code part 1

Author: rssalessio

examples/1_example.py

@@ -4,14 +4,14 @@
 from vrft import *
 
 def generateNoise(t):
-	omega = 2*np.pi*100
-	xi = 0.9
-	dt = t[1] - t[0]
-	noise =  np.random.normal(0,0.1,t.size)
-	tf = scipysig.TransferFunction([10*omega**2], [1, 2*xi*omega, omega**2])
-	# Second order system
-	_, yn, _ = scipysig.lsim(tf, noise, t)
-	return yn
+    omega = 2*np.pi*100
+    xi = 0.9
+    dt = t[1] - t[0]
+    noise =  np.random.normal(0,0.1,t.size)
+    tf = scipysig.TransferFunction([10*omega**2], [1, 2*xi*omega, omega**2])
+    # Second order system
+    _, yn, _ = scipysig.lsim(tf, noise, t)
+    return yn
 
 #Generate time and u(t) signals
 t_start = 0
@@ -28,7 +28,7 @@ def generateNoise(t):
 sys = ExtendedTF(num, den, dt=t_step)
 t, y = scipysig.dlsim(sys, u, t)
 y = y[:, 0]
-y += generateNoise(t)
+#y += generateNoise(t)
 data = iddata(y, u, t_step, [0])
 
 
@@ -37,7 +37,7 @@ def generateNoise(t):
 
 #PI Controller
 base = [ExtendedTF([1], [1, -1], dt=t_step),
-	    ExtendedTF([1, 0], [1, -1], dt=t_step)]
+        ExtendedTF([1, 0], [1, -1], dt=t_step)]
 
 #Experiment filter
 L = refModel * (1 -  refModel)
@@ -47,10 +47,10 @@ def generateNoise(t):
 #Obtained controller
 print("Controller: {}".format(C))
 
-L = C * sys
-L = feedback(L)
+L = (C * sys).feedback()
 
 print("Theta: {}".format(theta))
+print(scipysig.ZerosPolesGain(L))
 
 #Analysis
 t = t[:len(r)]

test/test_tf.py renamed to test/test_utils.py

@@ -1,8 +1,25 @@
 from unittest import TestCase
 from vrft.utilities.utils import *
+from vrft.vrft.vrft_algo import virtualReference
 import numpy as np
+import scipy.signal as scipysig
+
+class TestUtils(TestCase):
+    def test_deconvolve(self):
+        t_start = 0
+        t_end = 10
+        t_step = 1e-2
+        t = np.arange(t_start, t_end, t_step)
+        sys = ExtendedTF([0.5], [1, -0.9], dt=t_step)
+        u = np.random.normal(size=t.size)
+        _, y = scipysig.dlsim(sys, u, t)
+        y = y[:, 0]
+        data = iddata(y, u, t_step, [0])
+        r1, _ = virtualReference(data, sys.num, sys.den)
+        r2 = deconvolve_signal(sys, data.y, data.ts)
+        self.assertTrue(np.linalg.norm(r2-r1[:r2.size], np.infty) <  1e-3)
+
 
-class TestTF(TestCase):
     def test_checkSystem(self):
         a = [1, 0, 1]
         b = [1, 0, 2]

vrft/utilities/utils.py renamed to vrft/extended_tf.py

@@ -2,83 +2,11 @@
 
 import numpy as np
 import scipy.signal as scipysig
-from .iddata import iddata
 
 from scipy.signal.ltisys import TransferFunction as TransFun
 from numpy import polymul, polyadd
 
-def Doperator(p: int, q: int, x: float) -> np.ndarray:
-    D = np.zeros((p * q, q))
-    for i in range(q):
-        D[i * p:(i + 1) * p, i] = x
-    return D
 
-def checkSystem(num: np.ndarray, den: np.ndarray) -> bool:
-    try:
-        M, N = systemOrder(num, den)
-    except ValueError:
-        raise
-
-    if (N < M):
-        raise ValueError("The system is not causal.")
-
-    return True
-
-def systemOrder(num: np.ndarray, den: np.ndarray) -> tuple:
-    den = den if isinstance(den, np.ndarray) else np.array([den]).flatten()
-    num = num if isinstance(num, np.ndarray) else np.array([num]).flatten()
-
-    if num.ndim == 0:
-        num = np.expand_dims(num, axis=0)
-
-    if den.ndim == 0:
-        den = np.expand_dims(den, axis=0)
-
-    N = den.size
-    M = num.size
-
-    denOrder = -1
-    numOrder = -1
-
-    for i, d in enumerate(den):
-        if not np.isclose(d, 0):
-            denOrder = N - i - 1
-            break
-
-    for i, n in enumerate(num):
-        if not np.isclose(n, 0):
-            numOrder = M - i - 1
-            break
-
-    if (denOrder == -1):
-        raise ValueError("Denominator can not be zero.")
-
-    if (numOrder == -1):
-        raise ValueError("Numerator can not be zero.")
-
-    return (numOrder, denOrder)
-
-
-def filter_iddata(data: iddata, L: scipysig.dlti) -> iddata:
-    t_start = 0
-    t_step = data.ts
-    t_end = len(data.y) * t_step
-
-    t = np.arange(t_start, t_end, t_step)
-
-    if (L != None):
-        t, y = scipysig.dlsim(L, data.y, t)
-        t, u = scipysig.dlsim(L, data.u, t)
-        return iddata(y[:, 0], u[:, 0], data.ts, data.y0)
-    return data
-
-def deconvolve_signal(T: scipysig.dlti, x: np.ndarray, dt: float) -> np.ndarray:
-    impulse = scipysig.dimpulse(T)[1][0].flatten()
-    idx1 = np.argwhere(impulse != 0)[0].item()
-    idx2 = np.argwhere(np.isclose(impulse[idx1:], 0.) == True)
-    idx2 = -1 if not idx2 else idx2[0].item()
-    recovered, _ = scipysig.deconvolve(x, impulse[idx1:idx2])
-    return recovered[np.argwhere(impulse != 0)[0].item():]
 
 class ExtendedTF(scipysig.ltisys.TransferFunctionDiscrete):
     def __init__(self, num, den, dt):
@@ -152,12 +80,15 @@ def __rsub__(self,other):
                 denom = polymul(self.den,other.den)
             return ExtendedTF(numer,denom, dt=self._dt)
 
+    def feedback(self):
+        num = self.num
+        den = self.den
+        den = polyadd(num, den)
+        self = ExtendedTF(num, den, dt=self.dt)
+        return self
+
     # sheer laziness: symmetric behaviour for commutative operators
     __rmul__ = __mul__
     __radd__ = __add__
 
-def feedback(L):
-    num = L.num
-    den = L.den
-    den = polyadd(num, den)
-    return ExtendedTF(num, den, dt=L.dt)
+

vrft/utilities/iddata.py renamed to vrft/iddata.py

@@ -0,0 +1,81 @@
+from __future__ import division
+
+import numpy as np
+import scipy.signal as scipysig
+from .iddata import iddata
+
+from scipy.signal.ltisys import TransferFunction as TransFun
+from numpy import polymul, polyadd
+
+def Doperator(p: int, q: int, x: float) -> np.ndarray:
+    D = np.zeros((p * q, q))
+    for i in range(q):
+        D[i * p:(i + 1) * p, i] = x
+    return D
+
+def checkSystem(num: np.ndarray, den: np.ndarray) -> bool:
+    try:
+        M, N = systemOrder(num, den)
+    except ValueError:
+        raise
+
+    if (N < M):
+        raise ValueError("The system is not causal.")
+
+    return True
+
+def systemOrder(num: np.ndarray, den: np.ndarray) -> tuple:
+    den = den if isinstance(den, np.ndarray) else np.array([den]).flatten()
+    num = num if isinstance(num, np.ndarray) else np.array([num]).flatten()
+
+    if num.ndim == 0:
+        num = np.expand_dims(num, axis=0)
+
+    if den.ndim == 0:
+        den = np.expand_dims(den, axis=0)
+
+    N = den.size
+    M = num.size
+
+    denOrder = -1
+    numOrder = -1
+
+    for i, d in enumerate(den):
+        if not np.isclose(d, 0):
+            denOrder = N - i - 1
+            break
+
+    for i, n in enumerate(num):
+        if not np.isclose(n, 0):
+            numOrder = M - i - 1
+            break
+
+    if (denOrder == -1):
+        raise ValueError("Denominator can not be zero.")
+
+    if (numOrder == -1):
+        raise ValueError("Numerator can not be zero.")
+
+    return (numOrder, denOrder)
+
+
+def filter_iddata(data: iddata, L: scipysig.dlti) -> iddata:
+    t_start = 0
+    t_step = data.ts
+    t_end = len(data.y) * t_step
+
+    t = np.arange(t_start, t_end, t_step)
+
+    if (L != None):
+        t, y = scipysig.dlsim(L, data.y, t)
+        t, u = scipysig.dlsim(L, data.u, t)
+        return iddata(y[:, 0], u[:, 0], data.ts, data.y0)
+    return data
+
+def deconvolve_signal(T: scipysig.dlti, x: np.ndarray, dt: float) -> np.ndarray:
+    impulse = scipysig.dimpulse(T)[1][0].flatten()
+    idx1 = np.argwhere(impulse != 0)[0].item()
+    idx2 = np.argwhere(np.isclose(impulse[idx1:], 0.) == True)
+    idx2 = -1 if idx2.size == 0 else idx2[0].item()
+    recovered, _ = scipysig.deconvolve(x, impulse[idx1:idx2])
+    return recovered[np.argwhere(impulse != 0)[0].item():]

vrft/utils.py

@@ -1,5 +1,5 @@
 from vrft.utilities.iddata import iddata
-from vrft.utilities.utils import systemOrder, checkSystem, filter_iddata
+from vrft.utilities.utils import systemOrder, checkSystem, filter_iddata, deconvolve_signal
 import numpy as np
 import scipy.signal as scipysig
 
@@ -94,14 +94,22 @@ def control_response(data: iddata, error: np.ndarray, control: list):
     return phi
 
 def compute_vrft(data: iddata, refModel: scipysig.dlti, control: list, L: scipysig.dlti):
+    # import pdb
+    # import matplotlib.pyplot as plt
+    # pdb.set_trace()
     data = filter_iddata(data, L)
     r, n = virtualReference(data,
                          refModel.num,
                          refModel.den)
+    # r2=deconvolve_signal(refModel, data.y, data.ts)
+    # plt.plot(r)
+    # plt.plot(r2)
+    # plt.show()
 
-    phi = control_response(data, np.subtract(r, data.y[:-n]), control)
+    phi = control_response(data, np.subtract(r, data.y[:n]), control)
     theta, phi = calc_minimum(data, phi)
     loss = compute_vrft_loss(data, phi, theta)
-
+    print(phi)
+    #theta[0] = -theta[0]
     final_control = np.dot(theta, control)
     return theta, r, phi, loss, final_control