iv part 2

Author: rssalessio

.gitignore

@@ -1,4 +1,6 @@
 *.pyc
 /dist/
 /*.egg-info
-/build/
+/build/
+/.ipynb_checkpoints
+/examples/.ipynb_checkpoints

README.md

@@ -40,6 +40,7 @@ python -m unittest
 - [**In Progress**][07.01.2020-] Add Documentation and Latex formulas
 - [**TODO**] Add MIMO Support
 - [**TODO**] Generalize to other kind of controllers (e.g., neural nets)
+- [**TODO**] Add Cython support
 
 ## Citations
 If you find this code useful in your research, please, consider citing it:

examples/1_example.py

@@ -58,7 +58,7 @@ def generateNoise(t):
 #Experiment filter
 L = refModel * (1 -  refModel)
 #VRFT
-theta, r, _, _, C = compute_vrft(data, refModel, base, L)
+theta, r, loss, C = compute_vrft(data, refModel, base, L)
 
 #Obtained controller
 print("Controller: {}".format(C))

examples/2_example.py

@@ -0,0 +1,110 @@
+# Copyright [2020] [Alessio Russo - alessior@kth.se]  
+# This file is part of PythonVRFT.
+# PythonVRFT is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, version 3 of the License.
+# PythonVRFT is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+# You should have received a copy of the GNU General Public License
+# along with PythonVRFT.  If not, see <http://www.gnu.org/licenses/>.
+#
+# Code author: [Alessio Russo - alessior@kth.se]
+# Last update: 07th January 2020, by alessior@kth.se
+#
+
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.signal as scipysig
+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
+
+#Generate time and u(t) signals
+t_start = 0
+t_end = 10
+t_step = 1e-2
+t = np.arange(t_start, t_end, t_step)
+u = np.ones(len(t))
+u[200:400] = np.zeros(200)
+u[600:800] = np.zeros(200)
+
+#Experiment
+num = [0.5]
+den = [1, -0.9]
+sys = ExtendedTF(num, den, dt=t_step)
+t, y = scipysig.dlsim(sys, u, t)
+y = y.flatten() + 0.5 * np.random.normal(size = t.size)
+#y += generateNoise(t)
+data = iddata(y, u, t_step, [0])
+
+
+#Reference Model
+refModel = ExtendedTF([0.6], [1, -0.4], dt=t_step)
+
+#PI Controller
+base = [ExtendedTF([1], [1, -1], dt=t_step),
+        ExtendedTF([1, 0], [1, -1], dt=t_step)]
+
+#Experiment filter
+L = refModel * (1 -  refModel)
+
+#VRFT
+theta, r, loss, C = compute_vrft(data, refModel, base, L, iv=True)
+
+#Obtained controller
+print("Controller: {}".format(C))
+
+L = (C * sys).feedback()
+
+print("Theta: {}".format(theta))
+print(scipysig.ZerosPolesGain(L))
+
+#Analysis
+t = t[:len(r)]
+u = np.ones(len(t))
+_, yr = scipysig.dlsim(refModel, u, t)
+_, yc = scipysig.dlsim(L, u, t)
+_, ys = scipysig.dlsim(sys, u, t)
+
+yr = np.array(yr).flatten()
+ys = np.array(ys).flatten()
+yc = np.array(yc).flatten()
+fig, ax = plt.subplots(4, sharex=True)
+ax[0].plot(t, yr,label='Ref System')
+ax[0].plot(t, yc, label='CL System')
+ax[0].set_title('Systems response')
+ax[0].grid(True)
+ax[1].plot(t, ys, label='OL System')
+ax[1].set_title('OL Systems response')
+ax[1].grid(True)
+ax[2].plot(t, y[:len(r)])
+ax[2].grid(True)
+ax[2].set_title('Experiment data')
+ax[3].plot(t, r)
+ax[3].grid(True)
+ax[3].set_title('Virtual Reference')
+
+# Now add the legend with some customizations.
+legend = ax[0].legend(loc='lower right', shadow=True)
+
+# The frame is matplotlib.patches.Rectangle instance surrounding the legend.
+frame = legend.get_frame()
+frame.set_facecolor('0.90')
+
+# Set the fontsize
+for label in legend.get_texts():
+    label.set_fontsize('large')
+
+for label in legend.get_lines():
+    label.set_linewidth(1.5)  # the legend line width
+plt.show()

examples/notebook_example_1.ipynb

@@ -24,62 +24,62 @@
 
 class TestIDData(TestCase):
     def test_type(self):
-        a = iddata(0.0, 0.0, 0.0)
+        a = iddata(0.0, 0.0, 0.0, [0])
         with self.assertRaises(ValueError):
             a.check()
 
-        a =  iddata(0.0, [1], 0.0)
+        a =  iddata(0.0, [1], 0.0, [0])
         with self.assertRaises(ValueError):
             a.check()
 
-        a = iddata(np.zeros(10), 1, 0.0)
+        a = iddata(np.zeros(10), 1, 0.0, [0])
         with self.assertRaises(ValueError):
             a.check()
 
-        a = iddata([0 for i in range(10)], [0 for i in range(10)], 1.0)
+        a = iddata([0 for i in range(10)], [0 for i in range(10)], 1.0, [0])
         self.assertTrue(a.check())
 
-        a = iddata(np.zeros(10), np.zeros(10), 1.0)
+        a = iddata(np.zeros(10), np.zeros(10), 1.0, [0])
         self.assertTrue(a.check())
 
     def test_size(self):
-        a = iddata(np.zeros(10), np.zeros(10), 0.0)
+        a = iddata(np.zeros(10), np.zeros(10), 0.0, [0])
         self.assertEqual(len(a.y), 10)
         self.assertEqual(len(a.u), 10)
         self.assertEqual(len(a.y), len(a.u))
 
-        a = iddata([0 for i in range(10)], [1 for i in range(0,10)], 0.0)
+        a = iddata([0 for i in range(10)], [1 for i in range(0,10)], 0.0, [0])
         self.assertEqual(len(a.y), 10)
         self.assertEqual(len(a.u), 10)
         self.assertEqual(len(a.y), len(a.u))
 
-        a = iddata(np.zeros(10), np.zeros(9), 0.0)
+        a = iddata(np.zeros(10), np.zeros(9), 0.0, [0])
         with self.assertRaises(ValueError):
             a.check()
 
-        a = iddata(np.zeros(8), np.zeros(9), 0.0)
+        a = iddata(np.zeros(8), np.zeros(9), 0.0, [0])
         with self.assertRaises(ValueError):
             a.check()
 
 
     def test_sampling_time(self):
-        a = iddata(np.zeros(10), np.zeros(10), 0.0)
+        a = iddata(np.zeros(10), np.zeros(10), 0.0, [0])
         with self.assertRaises(ValueError):
             a.check()
 
-        a = iddata(np.zeros(10), np.zeros(10), 1e-9)
+        a = iddata(np.zeros(10), np.zeros(10), 1e-9, [0])
         with self.assertRaises(ValueError):
             a.check()
 
-        a = iddata(np.zeros(10), np.zeros(10), -0.1)
+        a = iddata(np.zeros(10), np.zeros(10), -0.1, [0])
         with self.assertRaises(ValueError):
             a.check()
 
-        a = iddata(np.zeros(10), np.zeros(10), 0.1)
+        a = iddata(np.zeros(10), np.zeros(10), 0.1, [0])
         self.assertTrue(a.check())
 
     def test_copy(self):
-        a = iddata(np.zeros(10), np.zeros(10), 0.1)
+        a = iddata(np.zeros(10), np.zeros(10), 0.1, [0])
         b = a.copy()
         self.assertTrue(a.check())
         self.assertTrue(b.check())
@@ -90,7 +90,7 @@ def test_copy(self):
         self.assertTrue(a.ts == b.ts)
 
     def test_filter(self):
-        a = iddata(np.zeros(10), np.zeros(10), 0.1)
+        a = iddata(np.zeros(10), np.zeros(10), 0.1, [0])
         L = scipysig.dlti([1], [1], dt=0.1)
         b = a.copy()
         a.filter(L)
@@ -100,23 +100,24 @@ def test_filter(self):
         self.assertTrue(a.ts == b.ts)
 
     def test_split(self):
-        n = 10
-        a = iddata(np.random.normal(size=n), np.random.normal(size=n), 0.1)
+        n = 9
+        a = iddata(np.random.normal(size=n), np.random.normal(size=n), 0.1, [0])
 
         b, c = a.split()
+        n0 = len(a.y0)
+        n1 = (n + n0) // 2
 
         self.assertTrue(b.y.size == c.y.size)
         self.assertTrue(b.u.size == c.u.size)
-        self.assertTrue(b.y0 == c.y0)
         self.assertTrue(b.ts == c.ts)
         self.assertTrue(b.ts == a.ts)
-        self.assertTrue(np.all(b.y == a.y[:n//2]))
-        self.assertTrue(np.all(b.u == a.u[:n//2]))
-        self.assertTrue(b.y0 == a.y0)
+        self.assertTrue(np.all(b.y == a.y[:n1 - n0]))
+        self.assertTrue(np.all(b.u == a.u[:n1 - n0]))
+        self.assertTrue(np.all(b.y0 == a.y0))
 
-        self.assertTrue(np.all(c.y == a.y[n//2:]))
-        self.assertTrue(np.all(c.u == a.u[n//2:]))
-        self.assertTrue(c.y0 == a.y0)
+        self.assertTrue(np.all(c.y == a.y[n1:n]))
+        self.assertTrue(np.all(c.u == a.u[n1:n]))
+        self.assertTrue(np.all(c.y0 == a.y[n1 - n0:n1]))
 
         y0 = [-1, 2]
         a = iddata(np.random.normal(size=n), np.random.normal(size=n), 0.1, y0)
@@ -132,8 +133,8 @@ def test_split(self):
         self.assertTrue(np.all(b.u == a.u[:n1 - n0]))
         self.assertTrue(np.all(b.y0 == a.y0))
 
-        self.assertTrue(np.all(c.y == a.y[n1:]))
-        self.assertTrue(np.all(c.u == a.u[n1:]))
+        self.assertTrue(np.all(c.y == a.y[n1:n-1]))
+        self.assertTrue(np.all(c.u == a.u[n1:n-1]))
         self.assertTrue(np.all(c.y0 == a.y[n1 - n0:n1]))
 
 

examples/notebook_example_2.ipynb

@@ -53,15 +53,10 @@ def test_virtualReference(self):
         y = y[:,0]
         data = iddata(y,u,t_step,[0,0])
 
-        r, _ = virtualReference(data, num, den)
-        for i in range(len(r)):
-            self.assertTrue(np.isclose(r[i], u[i]))
-
-        data = iddata(y,u,t_step,[0,0,0])
-        r, _ = virtualReference(data, num, den)
-        for i in range(len(r)):
-            self.assertTrue(np.isclose(r[i], u[i]))
-
+        # wrong initial conditions
+        with self.assertRaises(ValueError):
+            r, _ = virtualReference(data, num, den)
+    
         #test good data, first order
         data = iddata(y,u,t_step,[0])
 

test/test_iddata.py

@@ -27,61 +27,66 @@
 class TestVRFT(TestCase):
     def test_vrft(self):
         t_start = 0
-        t_end = 10
         t_step = 1e-2
-        t = np.arange(t_start, t_end, t_step)
-        u = np.ones(len(t)).tolist()
+        t_ends = [10, 10 + t_step]
 
-        num = [0.1]
-        den = [1, -0.9]
-        sys = scipysig.TransferFunction(num, den, dt=t_step)
-        t, y = scipysig.dlsim(sys, u, t)
-        y = y[:,0]
-        data = iddata(y,u,t_step,[0])
+        expected_theta = np.array([1.93220784, -1.05808206, 1.26623764, 0.0088772])
+        expected_loss = 0.00064687904235295
+        
+        for t_end in t_ends:
+            t = np.arange(t_start, t_end, t_step)
+            u = np.ones(len(t)).tolist()
 
-        refModel = ExtendedTF([0.2], [1, -0.8], dt=t_step)
-        prefilter = refModel * (1-refModel)
+            num = [0.1]
+            den = [1, -0.9]
+            sys = scipysig.TransferFunction(num, den, dt=t_step)
+            t, y = scipysig.dlsim(sys, u, t)
+            y = y[:,0]
+            data = iddata(y,u,t_step,[0])
 
-        control = [ExtendedTF([1], [1,0], dt=t_step),
-                ExtendedTF([1], [1,0,0], dt=t_step),
-                ExtendedTF([1], [1,0,0,0], dt=t_step),
-                ExtendedTF([1, 0], [1,1], dt=t_step)]
+            refModel = ExtendedTF([0.2], [1, -0.8], dt=t_step)
+            prefilter = refModel * (1-refModel)
 
-        theta1, _, loss1, _ = compute_vrft(data, refModel, control, prefilter)
-        theta2, _, loss2, _ = compute_vrft([data], refModel, control, prefilter)
-        theta3, _, loss3, _ = compute_vrft([data, data], refModel, control, prefilter)
-        
-        self.assertTrue(np.isclose(loss1, loss2))
-        self.assertTrue(np.isclose(loss1, loss3))
-        self.assertTrue(np.linalg.norm(theta1-theta2)<1e-15)
-        self.assertTrue(np.linalg.norm(theta1-theta3)<1e-15)
+            control = [ExtendedTF([1], [1,0], dt=t_step),
+                    ExtendedTF([1], [1,0,0], dt=t_step),
+                    ExtendedTF([1], [1,0,0,0], dt=t_step),
+                    ExtendedTF([1, 0], [1,1], dt=t_step)]
+
+            theta1, _, loss1, _ = compute_vrft(data, refModel, control, prefilter)
+            theta2, _, loss2, _ = compute_vrft([data], refModel, control, prefilter)
+            theta3, _, loss3, _ = compute_vrft([data, data], refModel, control, prefilter)
+            
+            self.assertTrue(np.isclose(loss1, loss2))
+            self.assertTrue(np.isclose(loss1, loss3))
+            self.assertTrue(np.linalg.norm(theta1-theta2)<1e-15)
+            self.assertTrue(np.linalg.norm(theta1-theta3)<1e-15)
+            self.assertTrue(np.linalg.norm(theta1-expected_theta, np.infty) < 1e-5)
+            self.assertTrue(abs(expected_loss - loss1) < 1e-5)
 
     def test_iv(self):
         t_start = 0
-        t_end = 10
         t_step = 1e-2
-        t = np.arange(t_start, t_end, t_step)
-        u = np.ones(len(t)).tolist()
+        t_ends = [10, 10 + t_step]
 
-        num = [0.1]
-        den = [1, -0.9]
-        sys = scipysig.TransferFunction(num, den, dt=t_step)
-        t, y = scipysig.dlsim(sys, u, t)
-        y = y.flatten() + np.random.normal(size=t.size)
-        data = iddata(y,u,t_step,[0])
-
-        refModel = ExtendedTF([0.2], [1, -0.8], dt=t_step)
-        prefilter = refModel * (1-refModel)
+        
+        for t_end in t_ends:
+            t = np.arange(t_start, t_end, t_step)
+            u = np.ones(len(t)).tolist()
 
-        control = [ExtendedTF([1], [1,0], dt=t_step),
-                ExtendedTF([1], [1,0,0], dt=t_step),
-                ExtendedTF([1], [1,0,0,0], dt=t_step),
-                ExtendedTF([1, 0], [1,1], dt=t_step)]
+            num = [0.1]
+            den = [1, -0.9]
+            sys = scipysig.TransferFunction(num, den, dt=t_step)
+            t, y = scipysig.dlsim(sys, u, t)
+            y = y.flatten() + 1e-2 * np.random.normal(size=t.size)
+            data = iddata(y,u,t_step,[0])
 
-        #import pdb
-        #pdb.set_trace()
-        
-        #theta, _, loss, _ = compute_vrft(data, refModel, control, prefilter, iv=True)
+            refModel = ExtendedTF([0.2], [1, -0.8], dt=t_step)
+            prefilter = refModel * (1-refModel)
 
+            control = [ExtendedTF([1], [1,0], dt=t_step),
+                    ExtendedTF([1], [1,0,0], dt=t_step),
+                    ExtendedTF([1], [1,0,0,0], dt=t_step),
+                    ExtendedTF([1, 0], [1,1], dt=t_step)]
 
+            theta, _, loss, _ = compute_vrft(data, refModel, control, prefilter, iv=True)