update: control, GBN and validation

updated file functionset (GBN and validation)
updated models (solved issues about version of control)
updated examples accordingly

Author: RBdC

.travis.yml

@@ -19,3 +19,4 @@ script:
   - python ARMAX.py
   - python ARX_MIMO.py
   - python SS.py
+  - python example_CST.py

Examples/ARMAX.py

@@ -24,7 +24,7 @@
 # ## Define pseudo random binary sequence as input signal and white noise as noise signal
 
 switch_probability = 0.08  # [0..1]
-Usim = fset.GBN_seq(npts, switch_probability)
+[Usim,_,_] = fset.GBN_seq(npts, switch_probability)
 white_noise_variance = [0.005]
 e_t = fset.white_noise_var(Usim.size, white_noise_variance)[0]
 
@@ -139,7 +139,7 @@
 
 switch_probability = 0.07  # [0..1]
 input_range = [0.5, 1.5]
-U_valid = fset.GBN_seq(npts, switch_probability, input_range)
+[U_valid,_,_] = fset.GBN_seq(npts, switch_probability, Range = input_range)
 white_noise_variance = [0.01]
 e_valid = fset.white_noise_var(U_valid.size, white_noise_variance)[0]
 

Examples/ARMAX_MIMO.py

@@ -109,10 +109,10 @@
 #                                   #INPUT#
 Usim = np.zeros((4, npts))
 Usim_noise = np.zeros((4, npts))
-Usim[0, :] = fset.GBN_seq(npts, 0.03, [-0.33, 0.1])
-Usim[1, :] = fset.GBN_seq(npts, 0.03)
-Usim[2, :] = fset.GBN_seq(npts, 0.03, [2.3, 5.7])
-Usim[3, :] = fset.GBN_seq(npts, 0.03, [8., 11.5])
+[Usim[0, :],_,_] = fset.GBN_seq(npts, 0.03, Range = [-0.33, 0.1])
+[Usim[1, :],_,_] = fset.GBN_seq(npts, 0.03)
+[Usim[2, :],_,_] = fset.GBN_seq(npts, 0.03, Range = [2.3, 5.7])
+[Usim[3, :],_,_] = fset.GBN_seq(npts, 0.03, Range = [8., 11.5])
 
 # Adding noise
 err_inputH = np.zeros((4, npts))
@@ -207,7 +207,7 @@
 plt.grid()
 plt.ylabel("Input 1 GBN")
 plt.xlabel("Time")
-plt.title("Input (Switch probability=0.03)")
+plt.title("Input (Switch probability=0.03) (identification data)")
 
 plt.subplot(4, 1, 2)
 plt.plot(Time, Usim[1, :])
@@ -325,10 +325,10 @@
 #                                   #INPUT#
 Usim = np.zeros((4, npts))
 Usim_noise = np.zeros((4, npts))
-Usim[0, :] = fset.GBN_seq(npts, 0.03, [0.33, 0.7])
-Usim[1, :] = fset.GBN_seq(npts, 0.03, [-2., -1.])
-Usim[2, :] = fset.GBN_seq(npts, 0.03, [1.3, 2.7])
-Usim[3, :] = fset.GBN_seq(npts, 0.03, [1., 5.2])
+[Usim[0, :],_,_] = fset.GBN_seq(npts, 0.03, Range = [0.33, 0.7])
+[Usim[1, :],_,_] = fset.GBN_seq(npts, 0.03, Range = [-2., -1.])
+[Usim[2, :],_,_] = fset.GBN_seq(npts, 0.03, Range = [1.3, 2.7])
+[Usim[3, :],_,_] = fset.GBN_seq(npts, 0.03, Range = [1., 5.2])
 err_inputH = np.zeros((4, npts))
 
 err_inputH = fset.white_noise_var(npts, var_list)
@@ -470,7 +470,7 @@
 plt.grid()
 plt.ylabel("Input 1 GBN")
 plt.xlabel("Time")
-plt.title("Input (Switch probability=0.03)")
+plt.title("Input (Switch probability=0.03) (validation data)")
 
 plt.subplot(4, 1, 2)
 plt.plot(Time, Usim[1, :])

Examples/ARX_MIMO.py

@@ -21,20 +21,9 @@
     from sippy import *
 
 import numpy as np
-import control
 import control.matlab as cnt
 from sippy import functionset as fset
-from distutils.version import StrictVersion
-
-if StrictVersion(control.__version__) >= StrictVersion('0.8.2'):
-	lsim = cnt.lsim
-else:
-	def lsim(sys, U = 0.0, T = None, X0 = 0.0):
-		U_ = U
-		if isinstance(U_, (np.ndarray, list)):
-			U_ = U_.T
-		return cnt.lsim(sys, U_, T, X0)
-		
+
 
 # generating transfer functions in z-transf.
 var_list = [50., 100., 1.]
@@ -119,33 +108,33 @@ def lsim(sys, U = 0.0, T = None, X0 = 0.0):
 #                                   #INPUT#
 Usim = np.zeros((4, npts))
 Usim_noise = np.zeros((4, npts))
-Usim[0, :] = fset.GBN_seq(npts, 0.03, [-0.33, 0.1])
-Usim[1, :] = fset.GBN_seq(npts, 0.03)
-Usim[2, :] = fset.GBN_seq(npts, 0.03, [2.3, 5.7])
-Usim[3, :] = fset.GBN_seq(npts, 0.03, [8., 11.5])
+[Usim[0, :],_,_] = fset.GBN_seq(npts, 0.03, Range = [-0.33, 0.1])
+[Usim[1, :],_,_] = fset.GBN_seq(npts, 0.03)
+[Usim[2, :],_,_] = fset.GBN_seq(npts, 0.03, Range = [2.3, 5.7])
+[Usim[3, :],_,_] = fset.GBN_seq(npts, 0.03, Range = [8., 11.5])
 
 # Adding noise
 err_inputH = np.zeros((4, npts))
 
 err_inputH = fset.white_noise_var(npts, var_list)
 
 err_outputH = np.ones((3, npts))
-err_outputH1, Time, Xsim = lsim(H_sample1, err_inputH[0, :], Time)
-err_outputH2, Time, Xsim = lsim(H_sample2, err_inputH[1, :], Time)
-err_outputH3, Time, Xsim = lsim(H_sample3, err_inputH[2, :], Time)
-
-Yout11, Time, Xsim = lsim(g_sample11, Usim[0, :], Time)
-Yout12, Time, Xsim = lsim(g_sample12, Usim[1, :], Time)
-Yout13, Time, Xsim = lsim(g_sample13, Usim[2, :], Time)
-Yout14, Time, Xsim = lsim(g_sample14, Usim[3, :], Time)
-Yout21, Time, Xsim = lsim(g_sample21, Usim[0, :], Time)
-Yout22, Time, Xsim = lsim(g_sample22, Usim[1, :], Time)
-Yout23, Time, Xsim = lsim(g_sample23, Usim[2, :], Time)
-Yout24, Time, Xsim = lsim(g_sample24, Usim[3, :], Time)
-Yout31, Time, Xsim = lsim(g_sample31, Usim[0, :], Time)
-Yout32, Time, Xsim = lsim(g_sample32, Usim[1, :], Time)
-Yout33, Time, Xsim = lsim(g_sample33, Usim[2, :], Time)
-Yout34, Time, Xsim = lsim(g_sample34, Usim[3, :], Time)
+err_outputH1, Time, Xsim = cnt.lsim(H_sample1, err_inputH[0, :], Time)
+err_outputH2, Time, Xsim = cnt.lsim(H_sample2, err_inputH[1, :], Time)
+err_outputH3, Time, Xsim = cnt.lsim(H_sample3, err_inputH[2, :], Time)
+
+Yout11, Time, Xsim = cnt.lsim(g_sample11, Usim[0, :], Time)
+Yout12, Time, Xsim = cnt.lsim(g_sample12, Usim[1, :], Time)
+Yout13, Time, Xsim = cnt.lsim(g_sample13, Usim[2, :], Time)
+Yout14, Time, Xsim = cnt.lsim(g_sample14, Usim[3, :], Time)
+Yout21, Time, Xsim = cnt.lsim(g_sample21, Usim[0, :], Time)
+Yout22, Time, Xsim = cnt.lsim(g_sample22, Usim[1, :], Time)
+Yout23, Time, Xsim = cnt.lsim(g_sample23, Usim[2, :], Time)
+Yout24, Time, Xsim = cnt.lsim(g_sample24, Usim[3, :], Time)
+Yout31, Time, Xsim = cnt.lsim(g_sample31, Usim[0, :], Time)
+Yout32, Time, Xsim = cnt.lsim(g_sample32, Usim[1, :], Time)
+Yout33, Time, Xsim = cnt.lsim(g_sample33, Usim[2, :], Time)
+Yout34, Time, Xsim = cnt.lsim(g_sample34, Usim[3, :], Time)
 
 Ytot1 = Yout11 + Yout12 + Yout13 + Yout14 + err_outputH1
 Ytot2 = Yout21 + Yout22 + Yout23 + Yout24 + err_outputH2
@@ -157,7 +146,6 @@ def lsim(sys, U = 0.0, T = None, X0 = 0.0):
 Ytot[1, :] = Ytot2.squeeze()
 Ytot[2, :] = Ytot3.squeeze()
 
-Ytot = np.column_stack([Ytot, np.ones((3, 1))])
 ##identification parameters
 ordersna = [na1, na2, na3]
 ordersnb = [[nb11, nb12, nb13, nb14], [nb21, nb22, nb23, nb24], [nb31, nb32, nb33, nb34]]
@@ -169,7 +157,7 @@ def lsim(sys, U = 0.0, T = None, X0 = 0.0):
 # output of the identified model
 # you can build g11, g12, etc. separately using the NUMERATOR and DENOMINATOR attributes
 # see how in the armax_MIMO example
-Yout_id, Time, Xsim = lsim(Id_sys.G, Usim, Time)
+Yout_id, Time, Xsim = cnt.lsim(Id_sys.G, Usim, Time)
 
 ######plot
 #  

Examples/SS.py

@@ -40,7 +40,7 @@
 
 # Input sequence
 U = np.zeros((1, npts))
-U[0] = fset.GBN_seq(npts, 0.05)
+[U[0],_,_] = fset.GBN_seq(npts, 0.05)
 
 ##Output
 x, yout = fsetSIM.SS_lsim_process_form(A, B, C, D, U)

Examples/example_CST.py

@@ -17,26 +17,11 @@
 from sippy import functionsetSIM as fsetSIM
 from sippy import *
 #
-# import functionset as fset
-# import functionsetSIM as fsetSIM
-# from __init__ import *
 #
 import numpy as np
-import control
 import control.matlab as cnt
 import matplotlib.pyplot as plt
 
-from distutils.version import StrictVersion
-
-if StrictVersion(control.__version__) >= StrictVersion('0.8.2'):
-	lsim = cnt.lsim
-else:
-	def lsim(sys, U = 0.0, T = None, X0 = 0.0):
-		U_ = U
-		if isinstance(U_, (np.ndarray, list)):
-			U_ = U_.T
-		return cnt.lsim(sys, U_, T, X0)
-
 # sampling time
 ts = 1.         # [min]
 
@@ -104,13 +89,13 @@ def Fdyn(X,U):
 F_min = 0.4
 F_max = 0.6 
 Range_GBN_1 = [F_min,F_max]
-U[0,:] = fset.GBN_seq(npts, prob_switch_1, Range = Range_GBN_1)    
+[U[0,:],_,_] = fset.GBN_seq(npts, prob_switch_1, Range = Range_GBN_1)    
 # Steam Flow rate W = U[1]          [kg/min]
 prob_switch_2 = 0.05
 W_min = 20
 W_max = 40
 Range_GBN_2 = [W_min,W_max]
-U[1,:] = fset.GBN_seq(npts, prob_switch_2, Range = Range_GBN_2)
+[U[1,:],_,_] = fset.GBN_seq(npts, prob_switch_2, Range = Range_GBN_2)
 
 # disturbance inputs as RW (random-walk)
 
@@ -226,4 +211,113 @@ def Fdyn(X,U):
     plt.xlabel("Time")
     if i == 0:
         plt.title('identification')
+        
+ 
+#### VALIDATION STAGE
+
+# Build new input sequences 
+U_val = np.zeros((m,npts))
+  
+# manipulated inputs as GBN
+# Input Flow rate Fin = F = U[0]    [m^3/min]
+prob_switch_1 = 0.05
+F_min = 0.4
+F_max = 0.6 
+Range_GBN_1 = [F_min,F_max]
+[U_val[0,:],_,_] = fset.GBN_seq(npts, prob_switch_1, Range = Range_GBN_1)  
+# Steam Flow rate W = U[1]          [kg/min]
+prob_switch_2 = 0.05
+W_min = 20
+W_max = 40
+Range_GBN_2 = [W_min,W_max]
+[U_val[1,:],_,_] = fset.GBN_seq(npts, prob_switch_2, Range = Range_GBN_2)
+
+# disturbance inputs as RW (random-walk)
+# Input Concentration Ca_in = U[2]  [kg salt/m^3 solution]
+Ca_0 = 10.0                         # initial condition
+sigma_Ca = 0.02                      # variation
+U_val[2,:] = fset.RW_seq(npts, Ca_0, sigma = sigma_Ca)
+# Input Temperature T_in            [°C]
+Tin_0 = 25.0                        # initial condition
+sigma_T = 0.1                       # variation
+U_val[3,:] = fset.RW_seq(npts, Tin_0, sigma = sigma_T)
+
+#### COLLECT DATA
 
+# Output Initial conditions
+Caout_0 = Ca_0
+Tout_0 = (ro*cp*U[0,0]*Tin_0 + U[1,0]*Lam)/(ro*cp*U[0,0])
+Xo1 = Caout_0*np.ones((1,npts))
+Xo2 = Tout_0*np.ones((1,npts))
+X_val = np.vstack((Xo1,Xo2))
+
+# Run Simulation
+for j in range(npts-1):  
+    # Explicit Runge-Kutta 4 (TC dynamics is integrateed by hand)        
+    Mx = 5                  # Number of elements in each time step
+    dt = ts/Mx              # integration step
+    # Output & Input
+    X0k = X_val[:,j]
+    Uk = U_val[:,j]
+    # Integrate the model
+    for i in range(Mx):         
+        k1 = Fdyn(X0k, Uk)
+        k2 = Fdyn(X0k + dt/2.0*k1, Uk)
+        k3 = Fdyn(X0k + dt/2.0*k2, Uk)
+        k4 = Fdyn(X0k + dt*k3, Uk)
+        Xk_1 = X0k + (dt/6.0)*(k1 + 2.0*k2 + 2.0*k3 + k4)
+    X_val[:,j+1] = Xk_1
+
+# Add noise (with assigned variances)
+var = [0.01, 0.05]    
+noise_val = fset.white_noise_var(npts,var)    
+
+# Build Output
+Y_val = X_val + noise_val
+
+
+# MODEL VALIDATION   
+        
+# ARX 
+Yv_arx = fset.validation(Id_ARX,U_val,Y_val,Time)
+
+# ARMAX
+Yv_armax = fset.validation(Id_ARMAX,U_val,Y_val,Time)
+
+
+# SS
+x_ss, Yv_ss = fsetSIM.SS_lsim_process_form(Id_SS.A,Id_SS.B,Id_SS.C,Id_SS.D,U_val,Id_SS.x0)
+
+
+##### PLOTS
+
+# Input
+plt.figure(3)
+str_input = ['F [m$^3$/min]', 'W [kg/min]', 'Ca$_{in}$ [kg/m$^3$]', 'T$_{in}$ [$^o$C]']
+for i in range(m):  
+    plt.subplot(m,1,i+1)
+    plt.plot(Time,U_val[i,:])
+    # plt.ylabel("Input " + str(i+1))
+    plt.ylabel(str_input[i])
+    plt.grid()
+    plt.xlabel("Time")
+    plt.axis([0, tfin, 0.95*np.amin(U_val[i,:]), 1.05*np.amax(U_val[i,:])])
+    if i == 0:
+        plt.title('validation')
+
+# Output
+plt.figure(4)
+str_output = ['Ca [kg/m$^3$]', 'T [$^o$C]']
+for i in range(p): 
+    plt.subplot(p,1,i+1)
+    plt.plot(Time,Y_val[i,:],'b')
+    plt.plot(Time,Yv_arx[i,:],'g')
+    plt.plot(Time,Yv_armax[i,:],'r')
+    plt.plot(Time,Yv_ss[i,:],'m')
+    # plt.ylabel("Output " + str(i+1))
+    plt.ylabel(str_output[i])
+    plt.legend(['Data','ARX','ARMAX','SS'])
+    plt.grid()
+    plt.xlabel("Time")
+    if i == 0:
+        plt.title('validation')

sippy/armax.py

@@ -5,7 +5,7 @@
 @author: Giuseppe Armenise
 """
 from __future__ import absolute_import, division, print_function
-import control as cnt
+import control.matlab as cnt
 from builtins import object
 import warnings
 

sippy/armaxMIMO.py

@@ -5,7 +5,7 @@
 @author: Giuseppe Armenise
 """
 from __future__ import absolute_import, division, print_function
-import control as cnt
+import control.matlab as cnt
 import sys
 from builtins import object
 

sippy/arx.py

@@ -8,7 +8,7 @@
 
 import sys
 from builtins import object
-import control as cnt
+import control.matlab as cnt
 from .functionset import *
 # from functionset import *
 

sippy/arxMIMO.py

@@ -5,7 +5,7 @@
 @author: Giuseppe Armenise
 """
 from __future__ import absolute_import, division, print_function
-import control as cnt
+import control.matlab as cnt
 import sys
 from builtins import object