Merge pull request #26 from rssalessio/iv

Iv

Author: rssalessio

.gitignore

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

README.md

@@ -1,4 +1,4 @@
-# PythonVRFT Library - Version 0.0.4
+# PythonVRFT Library - Version 0.0.5
 VRFT Adaptive Control Library written in Python. Aim of this library is to provide an implementation of the VRFT (Virtual Reference Feedback Tuning) algorithm.
 
 You can find the package also at the following [link](https://pypi.org/project/pythonvrft/)
@@ -25,7 +25,7 @@ pip install .
 ``` 
 
 ## Examples
-Examples are located in the examples/ folder. At the moment only 1 example is available.
+Examples are located in the examples/ folder. At the moment 2 examples are available. One without instrumental variables and one using instrumental variables.
 
 ## Tests
 To execute tests run the following command
@@ -34,12 +34,13 @@ python -m unittest
 ``` 
 
 ## Changelog
-- [**DONE - V0.0.2**][26.03.2017] Implement the basic VRFT algorithm (1 DOF. offline, linear controller, controller expressed as scalar product theta*f(z))
-- [**DONE - V0.0.3**][05.01.2020] Code refactoring and conversion to Python 3; Removed support for Python Control library.
-- [**In Progress - V0.0.4**][07.01.2020-] Add Documentation and Latex formulas
+- [**V. 0.0.2**][26.03.2017] Implement the basic VRFT algorithm (1 DOF. offline, linear controller, controller expressed as scalar product theta*f(z))
+- [**V. 0.0.3**][05.01.2020] Code refactoring and conversion to Python 3; Removed support for Python Control library.
+- [**V. 0.0.5**][07.01.2020] Add Instrumental Variables (IVs) Support
+- [**In Progress**][07.01.2020-] Add Documentation and Latex formulas
 - [**TODO**] Add MIMO Support
-- [**TODO**] Add IV 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:
@@ -50,3 +51,5 @@ If you find this code useful in your research, please, consider citing it:
 >  doi          = {},
 >  url          = { https://github.com/rssalessio/PythonVRFT }
 >}
+
+[![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0)

examples/example1.png

@@ -19,15 +19,12 @@
 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
+# Example 1
+# ------------
+# In this example we see how to apply VRFT to a simple SISO model
+# without any measurement noise.
+# Input data is generated using a square signal
+#
 
 #Generate time and u(t) signals
 t_start = 0
@@ -43,8 +40,7 @@ def generateNoise(t):
 den = [1, -0.9]
 sys = ExtendedTF(num, den, dt=t_step)
 t, y = scipysig.dlsim(sys, u, t)
-y = y[:, 0]
-#y += generateNoise(t)
+y = y.flatten()
 data = iddata(y, u, t_step, [0])
 
 
@@ -58,7 +54,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))
@@ -78,7 +74,7 @@ def generateNoise(t):
 yr = np.array(yr).flatten()
 ys = np.array(ys).flatten()
 yc = np.array(yc).flatten()
-fig, ax = plt.subplots(4, sharex=True)
+fig, ax = plt.subplots(4, sharex=True, figsize=(12,8), dpi= 100, facecolor='w', edgecolor='k')
 ax[0].plot(t, yr,label='Ref System')
 ax[0].plot(t, yc, label='CL System')
 ax[0].set_title('Systems response')

examples/1_example.py renamed to examples/example1.py

@@ -0,0 +1,114 @@
+# 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: 08th January 2020, by alessior@kth.se
+#
+
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.signal as scipysig
+from vrft import *
+
+# Example 2
+# ------------
+# In this example we see how to apply VRFT to a simple SISO model
+# with colored measurement noise (no instrumental variables)
+# Input data is generated using random normal noise
+#
+
+def generateNoise(t):
+    # Generate colored noise
+    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.random.normal(size=t.size)
+
+#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() + 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)
+
+#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, figsize=(12,8), dpi= 100, facecolor='w', edgecolor='k')
+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/1_example.png renamed to examples/example2.png

@@ -0,0 +1,114 @@
+# 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: 08th January 2020, by alessior@kth.se
+#
+
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.signal as scipysig
+from vrft import *
+
+# Example 3
+# ------------
+# In this example we see how to apply VRFT to a simple SISO model
+# with  measurement noise using instrumental variables
+# Input data is generated using random normal noise
+#
+
+#Generate time and u(t) signals
+t_start = 0
+t_end = 10
+dt = 1e-2
+t = np.array([i * dt for i in range(int(t_end/dt    ))])
+
+#Experiment
+num = [0.5]
+den = [1, -0.9]
+sys = ExtendedTF(num, den, dt=dt)
+
+def generate_data(sys, u, t):
+    t, y = scipysig.dlsim(sys, u, t)
+    y = y.flatten() + 0.5 * np.random.normal(size = t.size)
+    return iddata(y, u, dt, [0])
+
+u = np.random.normal(size=t.size)
+data1 = generate_data(sys, u, t)
+data2 = generate_data(sys, u, t)
+data = [data1, data2]
+
+#Reference Model
+refModel = ExtendedTF([0.6], [1, -0.4], dt=dt)
+
+#PI Controller
+control = [ExtendedTF([1], [1, -1], dt=dt),
+        ExtendedTF([1, 0], [1, -1], dt=dt)]
+
+#Experiment filter
+prefilter = refModel * (1 -  refModel)
+
+# VRFT method with Instrumental variables
+theta_iv, r_iv, loss_iv, C_iv = compute_vrft(data, refModel, control, prefilter, iv=True)
+
+# VRFT method without Instrumental variables
+theta_noiv, r_noiv, loss_noiv, C_noiv = compute_vrft(data1, refModel, control, prefilter, iv=False)
+
+#Obtained controller
+print('------IV------')
+print("Loss: {}\nTheta: {}\nController: {}".format(loss_iv, theta_iv, C_iv))
+print('------No IV------')
+print("Loss: {}\nTheta: {}\nController: {}".format(loss_noiv, theta_noiv, C_noiv))
+
+
+# Closed loop system
+closed_loop_iv = (C_iv * sys).feedback()
+closed_loop_noiv = (C_noiv * sys).feedback()
+
+t = t[:len(r_iv)]
+u = np.ones(len(t))
+
+_, yr = scipysig.dlsim(refModel, u, t)
+_, yc_iv = scipysig.dlsim(closed_loop_iv, u, t)
+_, yc_noiv = scipysig.dlsim(closed_loop_noiv, u, t)
+_, ys = scipysig.dlsim(sys, u, t)
+
+yr = yr.flatten()
+ys = ys.flatten()
+yc_noiv = yc_noiv.flatten()
+yc_iv = yc_iv.flatten()
+
+fig, ax = plt.subplots(4, sharex=True, figsize=(12,8), dpi= 100, facecolor='w', edgecolor='k')
+ax[0].plot(t, yr,label='Ref System')
+ax[0].plot(t, yc_iv, label='CL System - IV')
+ax[0].plot(t, yc_noiv, label='CL System - No IV')
+ax[0].set_title('CL 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, data1.y[:len(r_iv)])
+ax[2].grid(True)
+ax[2].set_title('Experiment data')
+ax[3].plot(t, r_iv)
+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')
+
+
+plt.show()