User Guide

User Manual:

Open the PDF directly: View PDF .
Page Count: 12

Download
Open PDF In Browser	View PDF

U NIVERSITY OF P ISA , D EPARTMENT OF C IVIL AND I NDUSTRIAL E NGINEERING

Systems Identification Package for Python
(SIPPY):
User Guide
Giuseppe Armenise, Riccardo Bacci di Capaci,
Marco Vaccari, Gabriele Pannocchia
August 6, 2018

1 I NTRODUCTION
This document is intended to give a guideline to the user who has to perform his own identification 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 I NSTALLATION 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 N ECESSARY PARAMETERS TO BE DEFINED
Using the following notation:
L: number of sampling times of the collected data,
n y : number of the controlled variables,
n u : number of the manipulated variables.
Being L >> max(n y , n u ).
The three necessary input arguments to define are:
1. Output data: y (∈ Rn y ×L or ∈ RL×n y );
2. Input data: u (∈ Rnu ×L or ∈ 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 parameters and their name cannot be changed. These arguments are written in the following Sections in typographical style.

2

3 O PTIONAL 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 O PTIONAL 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:
• n a : the n y output orders in a list containing n y integer numbers;
• n b : the n y × n u input orders in a list containing n y lists, in each of these lists there are
n u integer numbers (if an input has no influence on an output, set that order equal to
0);
• n c (only for ’ARMAX’): the n y error model orders in a list containing n y integer numbers;
• θ: the n y ×n u input delays in a list containing n y lists, in each of these lists there are n u
integer numbers;
For the ARX identification, add: ARX_orders=[n a ,n b ,θ ].
For the ARMAX identification, add: ARMAX_orders=[n a ,n b ,n c ,θ ].
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

Input 1

g 11 =

4z 3 +3.3z 2
z 5 −0.3z 4 −0.25z 3 −0.021z 2

2
g 21 = −85z 4−57.5z−27.7
3

Input 2

g 12 =

10z 2
z 5 −0.3z 4 −0.25z 3 −0.021z 2

g 22 = 71z+12.3
4
3

Input 3

g 13 =

Input 4

g 14 =

Error model

h1 =

7z 2 +5.5z+2.2

z 5 −0.3z 4 −0.25z 3 −0.021z 2
−0.9z 3 −0.11z 2
z 5 −0.3z 4 −0.25z 3 −0.021z 2
z 5 +0.85z 4 +0.32z 3

z 5 −0.3z 4 −0.25z 3 −0.021z 2

Output 3
g 31 =

z −0.4z

2 +0.432z
g 32 = 0.821z
4
3
2

z −0.4z

g 23 =

z −0.1z −0.3z

−0.1z 3

g 33 =

z 4 −0.4z 3

3
g 24 = 0.994z
4
3

h2 =

z 4 −0.4z 3

0.1z 3
z 4 −0.1z 3 −0.3z 2

g 34 = 40.891z+0.223
3
2

z −0.4z

z 4 +0.4z 3 +0.05z 2

0.2z 3
z 4 −0.1z 3 −0.3z 2

z −0.1z −0.3z

h3 =

z 4 +0.7z 3 +0.485z 2 +0.22z
z 4 −0.1z 3 −0.3z 2

• for the ARMAX case, by dafault: ARMAX_orders=[1,1,1,0]
The respective lists (n a , n b , n c and θ) 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 −1 as the backward shift operator and z as the forward shift operator, for each
output the identified (ARMAX) structure is:
(1 + a 1 z −1 + ... + a n a z −n a )y(t ) = (b 1,1 z −(1+θ1 ) + ... + b 1,nb1 z −(nb1 +θ1 ) )u 1 (t )+
+ (b 2,1 z −(1+θ2 ) + ... + b 2,nb2 z −(nb2 +θ2 ) )u 2 (t )+
+ ... + (1 + c 1 z

−1

+ ... + c nc z

−n c

(4.1)

)e(t )

Where n a denotes the output order, n bi denotes the order of the i-th input, θi denotes the
delay of the i-th input, n c denotes the error model order, a j denotes the j-th coefficient of the
output, b i , j denotes the j-th coefficient of the i-th input and c j denotes the j-th coefficient of
the error model. In the ARX structure n c = 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
n a = [3, 1, 2]
n b = [[2, 1, 3, 2], [3, 2, 1, 1], [1, 2, 1, 2]]
n c = [2, 2, 3]
θ = [[1, 2, 2, 1], [1, 2, 0, 0], [0, 1, 0, 2]]

4.2 S UBSPACE 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:
(
x k+1 = Ax k + Bu k + w k
(4.2)
y k = C x k + Du k + v k
where: y k ∈ Rn y , x k ∈ Rn , u k ∈ Rnu , w k ∈ Rn and v k ∈ Rn y are the system output, state, input,
state noise and output measurement noise respectively (the subscript "k" denotes the k − t h
sampling time), A ∈ Rn×n , B ∈ Rn×nu , C ∈ Rn y ×n , D ∈ Rn y ×nu are the system matrices.
Innovation form:
(
x k+1 = Ax k + Bu k + K e k
(4.3)
y k = C x k + Du k + e k
Predictor form:
(

x k+1 = A K x k + B K u k + K y k
y k = C x k + Du k + e k

(4.4)

where the following relations hold:
AK = A − K C
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 n is 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.
σmax : 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 D filled with zeros; add
SS_D_required=True to have the matrix D calculated.
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 A is 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 B and x 0 can 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 O PTIONAL ARGUMENTS IF IC=’AIC’, ’AIC C ’ 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 S UBSPACE 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 B and x 0 is 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 G of 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 G is 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 H is built with the line:

H=control.tf(NUMERATOR_H,DENOMINATOR_H,ts)
• Vn: estimated error norm, defined as:
Vn=

1
T
2L t r (²² )

where t r is the trace operator, ² is the error sequence, calculated as:
² = y − ŷ
being ŷ the 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 K is 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:
¶
µ ¶
QS
ρw ¡ T T ¢
ρw ρv ]
= E[
ST R
ρv

µ

where ρ w and ρ v are the estimated sequences of the state noise and output measurement noise (considering the outputs normalized with respect to their standard deviations). E denotes 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 process 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 U SEFUL 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 n u rows 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 predictor form, given the system matrices, the Kalman filter gain, the real output sequence
(array with n y rows and L columns, the input sequence (an array with n u rows 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 innovation form. This function is analogous to the previous one, using the system matrices A and B instead 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 M ETHODS SUMMARY
For each method, the considered parameters (excluding centering and tsample, being common to all methods) and the attributes of the system are reported in the Table 8.1.

9 P ROBLEMS 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 calculated" 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,
NUMERATOR, DENOMINATOR, Vn

’ARX’
’AIC’, ’AICc’ or ’BIC’

na_ord, nb_ord, delays

Same as above

’None’

ARMAX_orders,

G, H, na, nb, nc, theta, ts,

ARMAX_max_iterations

NUMERATOR, DENOMINATOR, Vn,

’ARMAX’

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,

SS_max_order, SS_fixed_order,

G, x0, Vn

’PARSIM-K’

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 argument, 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



---

Source Exif Data:

File Type                       : PDF
File Type Extension             : pdf
MIME Type                       : application/pdf
PDF Version                     : 1.5
Linearized                      : No
Page Count                      : 12
Producer                        : pdfTeX-1.40.18
Creator                         : TeX
Create Date                     : 2018:08:06 18:23:34+02:00
Modify Date                     : 2018:08:06 18:23:34+02:00
Trapped                         : False
PTEX Fullbanner                 : This is pdfTeX, Version 3.14159265-2.6-1.40.18 (TeX Live 2017) kpathsea version 6.2.3

EXIF Metadata provided by EXIF.tools