User Guide
User Manual:
Open the PDF directly: View PDF 
.
Page Count: 12
UNIVERSITY OF PISA, DEPARTMENT OF CIVIL AND INDUSTRIAL ENGINEERING
Systems Identification Package for Python
(SIPPY):
User Guide
Giuseppe Armenise, Riccardo Bacci di Capaci,
Marco Vaccari, Gabriele Pannocchia
August 6, 2018
1 INTRODUCTION
This document is intended to give a guideline to the user who has to perform his own identi-
fication with collected data.
The program is written both for Python 2.7 and 3.6 (supported on Windows \ Linux \ Mac).
The following packages must be installed: NumPy, SciPy, control, math, Slycot. In order to
make the code compatible to both Python 2.7 and 3.6, the package future has to be installed
(for more details see: http://python-future.org/index.html). This package can be installed via
pip install , or searching for future among non installed packages in Anaconda Navigator
(more information can be found at https://anaconda.org/anaconda/anaconda-navigator).
1.1 INSTALLATION E PACKAGE CONTENT
In order to make the installation easier, the user can simply use the quick command
python setup.py install
in order to gather all the required packages all together.
SIPPY is distributed as packed file SIPPY.zip that contains the following items:
• setup.py: run the installation for the whole package.
1
• SIPPY/__init__.py: main file containing the function that has to be recalled to perform
the identification
• SIPPY/functionset.py: file containing most of the functions used by the identification
functions and other useful functions (See Section 7).
• SIPPY/functionsetSIM.py: additional functions used by the Subspace identification
functions and other useful functions for state space models (See Section 7).
• Examples/*.py: folder containing various examples, stressing different aspects of the
package.
Every other file not mentioned above and inside the folder "SIPPY/", i.e. SIPPY/*.py, is called
by the main file, hence the user has not to use them.
The example file created by the user can be in the same folder as SIPPY/_init_.py of can be
in a level down folder e.g. the folder "Example“. In order to run the example, few import lines
must be inserted, e.g. take as reference the example files already in the zip file. For instance,
the following line is needed to import the main function "system_identification":
from SIPPY import *
The line used to perform the identification is:
Identified_system=system_identification(INPUT_ARGUMENTS)
The input arguments depend on the used identification method, while the output is an object
containing the attributes that define the model (see Section 6).
2 NECESSARY PARAMETERS TO BE DEFINED
Using the following notation:
L: number of sampling times of the collected data,
ny: number of the controlled variables,
nu: number of the manipulated variables.
Being L>> max(ny,nu).
The three necessary input arguments to define are:
1. Output data: y (∈Rny×Lor ∈RL×ny);
2. Input data: u (∈Rnu×Lor ∈RL×nu);
3. Identification method: a string.
The identification method can be chosen between: ’ARX’, ’ARMAX’, ’N4SID’, ’MOESP’, ’CVA’,
’PARSIM-P’, ’PARSIM-S’ or ’PARSIM-K’. For instance:
Identified_system=system_identification(y,u,’ARX’)
The optional parameters are all keyword arguments to be added after the necessary param-
eters and their name cannot be changed. These arguments are written in the following Sec-
tions in typographical style.
2
3 OPTIONAL PARAMETERS COMMON TO ALL METHODS
The user can center to zero the collected input-output data respect to:
• the mean value, in this case add: centering=’MeanVal’;
• the initial condition, in this case: centering=’InitVal’;
By default centering=’None’.
The user can use one of the implemented Information Criteria to select the model order:
• Akaike information criterion, in this case add IC=’AIC’;
• Corrected Akaike information criterion, in this case add: IC=’AICc’;
• Bayesian information criterion, in this case: IC=’BIC’.
By default IC=’None’.
The optional argument tsample is the sampling time used to build (see in Python Library
the functions control.tf and control.ss) the transfer functions of the model (see Section 6). By
default tsample=1.
4 OPTIONAL ARGUMENTS IF IC=’None’
If IC=’None’, the user has to define the model orders. Anyway, there is a default option.
The optional arguments for each model set are discussed below.
4.1 ’ARX’ AND ’ARMAX’ METHODS
The user has to define:
•na: the nyoutput orders in a list containing nyinteger numbers;
•nb: the ny×nuinput orders in a list containing nylists, in each of these lists there are
nuinteger numbers (if an input has no influence on an output, set that order equal to
0);
•nc(only for ’ARMAX’): the nyerror model orders in a list containing nyinteger num-
bers;
•θ: the ny×nuinput delays in a list containing nylists, in each of these lists there are nu
integer numbers;
For the ARX identification, add: ARX_orders=[na,nb,θ].
For the ARMAX identification, add: ARMAX_orders=[na,nb,nc,θ].
The alternative is a list containing integer numbers:
• for the ARX case, by default: ARX_orders=[1,1,0];
3
Table 4.1: Example of MIMO ARMAX system
Output 1 Output 2 Output 3
Input 1 g11 =4z3+3.3z2
z5−0.3z4−0.25z3−0.021z2g21 =−85z2−57.5z−27.7
z4−0.4z3g31 =0.2z3
z4−0.1z3−0.3z2
Input 2 g12 =10z2
z5−0.3z4−0.25z3−0.021z2g22 =71z+12.3
z4−0.4z3g32 =0.821z2+0.432z
z4−0.1z3−0.3z2
Input 3 g13 =7z2+5.5z+2.2
z5−0.3z4−0.25z3−0.021z2g23 =−0.1z3
z4−0.4z3g33 =0.1z3
z4−0.1z3−0.3z2
Input 4 g14 =−0.9z3−0.11z2
z5−0.3z4−0.25z3−0.021z2g24 =0.994z3
z4−0.4z3g34 =0.891z+0.223
z4−0.1z3−0.3z2
Error model h1=z5+0.85z4+0.32z3
z5−0.3z4−0.25z3−0.021z2h2=z4+0.4z3+0.05z2
z4−0.4z3h3=z4+0.7z3+0.485z2+0.22z
z4−0.1z3−0.3z2
• for the ARMAX case, by dafault: ARMAX_orders=[1,1,1,0]
The respective lists (na,nb,ncand θ) will contain numbers equal to the one specified.
The ARMAX identification has an iterative procedure and the number of maximum iterations
can be specified, by default: ARMAX_max_iterations=100.
Denoting z−1as the backward shift operator and zas the forward shift operator, for each
output the identified (ARMAX) structure is:
(1 +a1z−1+... +anaz−na)y(t)=(b1,1z−(1+θ1)+... +b1,nb1z−(nb1+θ1))u1(t)+
+(b2,1z−(1+θ2)+... +b2,nb2z−(nb2+θ2))u2(t)+
+... +(1 +c1z−1+... +cncz−nc)e(t)
(4.1)
Where nadenotes the output order, nbi denotes the order of the i-th input, θidenotes the
delay of the i-th input, ncdenotes the error model order, ajdenotes the j-th coefficient of the
output, bi,jdenotes the j-th coefficient of the i-th input and cjdenotes the j-th coefficient of
the error model. In the ARX structure nc=0.
The transfer functions of the system are returned as fractions of polynomials in z:
y=G(z)u+H(z)e
For instance, Table 4.1 shows the identified matrices obtained for a MIMO system using the
ARMAX method, in which
na=[3, 1, 2]
nb=[[2, 1, 3, 2],[3, 2, 1, 1], [1, 2, 1, 2]]
nc=[2, 2, 3]
θ=[[1, 2, 2, 1],[1, 2, 0, 0], [0, 1, 0, 2]]
4.2 SUBSPACE IDENTIFICATION METHODS:
4
’N4SID’, ’MOESP’, ’CVA’, ’PARSIM-P’, ’PARSIM-S’ AND ’PARSIM-K’
Subspace identification methods (SIM) identify a state space model structure, that can be
represented in three forms, recalled below.
Process form:
(xk+1=Axk+Buk+wk
yk=C xk+Duk+vk
(4.2)
where: yk∈Rny,xk∈Rn,uk∈Rnu,wk∈Rnand vk∈Rnyare the system output, state, input,
state noise and output measurement noise respectively (the subscript "k" denotes the k−th
sampling time), A∈Rn×n,B∈Rn×nu,C∈Rny×n,D∈Rny×nuare the system matrices.
Innovation form:
(xk+1=Axk+Buk+K ek
yk=C xk+Duk+ek
(4.3)
Predictor form:
(xk+1=AKxk+BKuk+K yk
yk=C xk+Duk+ek
(4.4)
where the following relations hold:
AK=A−KC
BK=B−K D
The user has to define the future and past horizons (SS_f and SS_p respectively). For ’N4SID’,
’MOESP’ and ’CVA’ methods, the future and past horizons are equal, set by default SS_f=20
(integer number).
For ’PARSIM-P’, ’PARSIM-S’ and ’PARSIM-K’ methods, the future and past horizons can be
set, by default: SS_f=20,SS_p=20 (integer numbers).
After performing the singular value decomposition (SVD), scheduled for the identification,
the model order nis chosen by the settings:
• A threshold value (SS_threshold, between 0.0 and 1.0) and the maximum order
(SS_max_order, integer number);
• A fixed order (SS_fixed_order, integer number).
The threshold value imply that all singular values such that:
σi
σmax
<SS_threshold are discarded.
σma x : largest singular value, σi: i-th singular value.
The maximum order value limits the order of the model to that value.
Using a fixed order value, the maximum order and the threshold value are not taken into
account.
By default: SS_threshold=0.1, while the other ones are not specified.
The other optional arguments are by default:
•SS_D_required=False;
5
•SS_A_stability=False (only for ’N4SID’, ’MOESP’ and ’CVA’).
If SS_D_required=False, the state space model will have the matrix Dfilled with zeros; add
SS_D_required=True to have the matrix Dcalculated.
For ’N4SID’, ’MOESP’ and ’CVA’ methods, setting SS_A_stability=True will guarantee the
stability of the matrix A. After the normal evaluation of the matrix A, if the calculated Ais not
stable, the stability is forced and a warning message is shown: "Forcing A stability". This can
lead to significant errors in the identified models.
Only in PARSIM-K method, a revaluation step for Band x0can be required by the user, adding
SS_PK_B_reval=True. This revaluation step may lead to better performances of the model
in the process form.
5 OPTIONAL ARGUMENTS IF IC=’AIC’, ’AICC’OR ’BIC’
Using an information criterion, the user has to define a range for the model orders.
5.1 ’ARX’ AND ’ARMAX’ METHODS
The ARMAX_max_iterations argument is the same described in Section 4.1.
The user has to define the ranges of the model orders:
• The output order range in a list with two integer numbers, by default: na_ord=[0,5];
• The input order range in a list with two integer numbers, by default: nb_ord=[1,5];
• The error order range (only for ARMAX) in a list with two integer numbers, by default:
nc_ord=[0,5];
• The input delays range in a list with two integer numbers, by default: delays=[0,5];
For ARX and ARMAX identification, the information criteria are implemented only in the SISO
case.
5.2 SUBSPACE IDENTIFICATION METHODS:
’N4SID’, ’MOESP’, ’CVA’, ’PARSIM-P’, ’PARSIM-S’ AND ’PARSIM-K’
The following arguments are the same described in Section 4.2: SS_f, SS_p,
SS_D_required, SS_A_stability. The order is selected in the defined range SS_orders, a
list with two integer numbers, by default: SS_orders=[1,10].
Adding SS_PK_B_reval=True, the revaluation step for Band x0is performed externally to
the evaluation of the best model according to the information criterion and therefore only as
a final step.
6
6 ATTRIBUTES OF THE MODEL
Depending on the selected identification method, the object returned by the function
"system_identification" has different attributes. E.g. to recall the attribute Gof the object
Identified_system, the command line is:
Identified_system.G
6.1 ’ARX’ AND ’ARMAX’ METHODS
The attributes for these methods are:
•na,nb,nc(∗),theta: respectively output orders, input orders, error model orders, input
delays;
•G: transfer functions that relate inputs with outputs;
•H: transfer functions of the error models;
•ts: sampling time of the transfer functions;
•NUMERATOR, DENOMINATOR: lists containing coefficients respectively of the numerators
and of the denominators of G. The Gis built with the line:
G=control.tf(NUMERATOR,DENOMINATOR,ts)
•NUMERATOR_H(∗), DENOMINATOR_H(∗): lists containing coefficients respectively of the
numerators and of the denominators of H. The His built with the line:
H=control.tf(NUMERATOR_H,DENOMINATOR_H,ts)
•Vn: estimated error norm, defined as:
Vn=1
2Ltr (²²T)
where tr is the trace operator, ²is the error sequence, calculated as:
²=y−ˆ
y
being ˆ
ythe output obtained from the identified model. This norm is calculated on the
outputs normalized with respect to their standard deviation.
In the subsection 6.1, the superscript (∗) means that the attribute can be found only in the
ARMAX structure.
7
6.2 ’N4SID’, ’MOESP’, ’CVA’, ’PARSIM-P’, ’PARSIM-S’ AND ’PARSIM-K’
METHODS
The attributes for these methods are:
•A, B, C, D: the system matrices of (4.2).
•A_K, B_K: the system matrices of (4.4). These matrices for ’N4SID’, ’MOESP’ and ’CVA’
methods can be obtained only if the Kalman filter gain Kis calculated.
•K: the Kalman filter gain. This matrix for ’N4SID’, ’MOESP’ and ’CVA’ methods can be
obtained only if the package "Slycot" is full-installed.
•ts: sampling time of the transfer function G;
•G: transfer function that relates inputs with outputs, obtained as follows:
G=control.ss(A,B,C,D,ts)
•x0: initial state estimate. For ’N4SID’, ’MOESP’ and ’CVA’ methods, x0 is filled with
zeros (no estimate of the initial state is made).
•Q(∗),R(∗),S(∗): covariance matrices defined as follows:
µQ S
STR¶=E[µρw
ρv¶¡ρT
wρT
v¢]
where ρwand ρvare the estimated sequences of the state noise and output measure-
ment noise (considering the outputs normalized with respect to their standard devia-
tions). Edenotes the expected value operator.
•Vn: the same Vn of subsection 6.1. This norm for PARSIM methods is calculated using
the predictor form (4.4), while for ’N4SID’, ’MOESP’ and ’CVA’ methods using the pro-
cess form (4.2), i.e. two values of this norm can be compared only if obtained from the
same class of methods. In PARSIM-K method, if the user selected SS_PK_B_reval=True,
the norm is calculated in the process form.
In the subsection 6.2, the superscript (∗) means that the attribute exists only if the selected
method is ’N4SID’, ’MOESP’ or ’CVA’.
7 USEFUL FUNCTIONS
To use the functions contained in the file "functionset.py", the following line is needed:
from SIPPY import functionset as fset
The functions that can be useful for the user are:
8
•fset.PRBS_seq: this function returns a pseudo-random binary sequence. The input
arguments are:
1. length_arg: size of the sequence;
2. prob_switch: probability of switch (0:no switch, 1:always switch);
3. Range (optional): Upper and lower values of the sequence, example: Range=[3.,15.]
(by default Range=[-1.0,1.0]).
E.g., using the command line:
prbs=fset.PRBS_seq(100,0.1,Range=[-10.0,5.0])
prbs is an array with size=100, 10 % of switch probability, elements that can be -10 or
5.
•fset.white_noise: this function adds a white noise to a signal. The input arguments
are:
1. y: original signal
2. A_rel: relative amplitude. The standard deviation of the returned white noise is
the standard deviation of y multiplied by A_rel.
The outputs are:
1. errors: the white noise;
2. y_err: equal to errors+y
The command line is:
errors,y_err=white_noise(y,A_rel)
•fset.white_noise_var: this function generates a white noise matrix (rows with zero
mean). The input arguments are:
1. L: size (number of columns of the returned matrix);
2. Var: variance list;
e.g. the following command line:
noise=white_noise_var(100,[1,2])
returns a matrix that has 100 columns and two row vectors with variances 1 and 2.
To use the functions contained in the file "functionsetSIM.py", the following line is needed:
from SIPPY import functionsetSIM as fsetSIM
The functions that can be useful for the user are:
9
•fsetSIM.SS_lsim_process_form: this function performs a simulation in the process
form, given the system matrices, the input sequence (an array with nurows and L
columns) and the initial state estimate (array with n rows and one column). The state
sequence and the output sequence are returned. The command line is:
x,y=fsetSIM.SS_lsim_process_form(A,B,C,D,u,x0)
•fsetSIM.SS_lsim_predictor_form: this function performs a simulation in the pre-
dictor form, given the system matrices, the Kalman filter gain, the real output sequence
(array with nyrows and L columns, the input sequence (an array with nurows and L
columns) and the initial state estimate (array with n rows and one column). The state
sequence and the estimated output sequence are returned. The command line is:
x,y_hat=SS_lsim_predictor_form(A_K,B_K,C,D,K,y,u,x0)
•fsetSIM.SS_lsim_innovation_form: this function performs a simulation in the in-
novation form. This function is analogous to the previous one, using the system matri-
ces Aand Binstead of A_K and B_K:
x,y_hat=SS_lsim_innovation_form(A,B,C,D,K,y,u,x0)
In these functions, the initial state is an optional argument and if it is not given, the initial
states are taken as zeros. The user can perform a simulation with these functions recalling
the attributes of the Identified_system, e.g.:
x,y=fsetSIM.SS_lsim_process_form(Identified_system.A,
Identified_system.B, Identified_system.C,Identified_system.D,
u,Identified_system.x0)
8 METHODS SUMMARY
For each method, the considered parameters (excluding centering and tsample, being com-
mon to all methods) and the attributes of the system are reported in the Table 8.1.
9 PROBLEMS AND SOLUTIONS
Here some problems that the user can meet are presented.
•LinAlgError: SVD did not converge. The user can try to solve the problem changing
SS_f and SS_p or, if this does not work, changing method.
• For ’N4SID’, ’MOESP’ and ’CVA’ methods, if the message "Kalman filter cannot be cal-
culated" is shown, the problem can be that the Slycot package is not well installed,
otherwise "control.dare" is not able to solve the Riccati equation. To check if Slycot
works, use the state space example (attached in Identification_code.zip), where the K
should be calculated. The Slycot package is available at:
10
Table 8.1: Methods and their parameters and attributes
Method IC= Method parameters Attributes
’None’ ARX_orders G,H,na,nb,theta,ts,
’ARX’ NUMERATOR,DENOMINATOR,Vn
’AIC’, ’AICc’ or ’BIC’ na_ord,nb_ord,delays Same as above
’None’ ARMAX_orders,G,H,na,nb,nc,theta,ts,
’ARMAX’ ARMAX_max_iterations NUMERATOR,DENOMINATOR,Vn,
NUMERATOR_H,DENOMINATOR_H
’AIC’, ’AICc’ or ’BIC’ na_ord,nb_ord,nc_ord,delays, Same as above
ARMAX_max_iterations
’None’ SS_f,SS_threshold,SS_max_order,A,B,C,D,A_K,B_K,K,ts,
’N4SID’, ’MOESP’ SS_fixed_order,G,x0,Q,R,S,Vn
or ’CVA’ SS_D_required,SS_A_stability
’AIC’, ’AICc’ or ’BIC’ SS_f,SS_orders, Same as above
SS_D_required,SS_A_stability
’None’ SS_f,SS_p,SS_threshold,A,B,C,D,A_K,B_K,K,ts,
’PARSIM-P’ or SS_max_order,SS_fixed_order,G,x0,Vn
’PARSIM-S’ SS_D_required
’AIC’, ’AICc’ or ’BIC’ SS_f,SS_p,SS_orders,SS_D_required Same as above
’None’ SS_f,SS_p,SS_threshold,A,B,C,D,A_K,B_K,K,ts,
’PARSIM-K’ SS_max_order,SS_fixed_order,G,x0,Vn
SS_D_required,SS_PK_B_reval
’AIC’, ’AICc’ or ’BIC’ SS_f,SS_p,SS_orders, Same as above
SS_D_required,SS_PK_B_reval
11
https://pypi.python.org/pypi/slycot/0.2.0
Alternatively the binaries can be found at:
https://www.lfd.uci.edu/~gohlke/pythonlibs/
• The identification is not well performed (e.g. returning an unstable system in State
space models). The user has to try different options, changing the centering argu-
ment, resetting the future and past horizons, lowering the SS_threshold value, and so
on. For PARSIM methods, an advice can be to set SS_p ≥SS_f.
12
