ARX model not identified

[tikull]

Dear CPCLAB-UNIPI team,
I have started using SIPPY as the python version of system identification for my Thesis. I am implementing it on raspberry pi. In the beginning, I implemented a small discrete-time linear system in MATLAB, and also identified (with 99.94% match) with ARX model of the order [1 2 0].
My system (g_sample11) looks as below,
(0.001 z - 0.001)/ ( z - 0.006738)
Now, I am again trying to identify the same system with SIPPY with below function,
Id_ARX = system_identification(Yout1, Usim, 'ARX', IC ='AIC', tsample = ts, na_ord = [1, 1], nb_ord = [2, 2], delays = [0, 0], ARX_mod = 'LLS')
and
Id_ARX = system_identification(Yout1, Usim, 'ARX', tsample = ts, ARX_orders = [na_ord, nb_ord, theta], ARX_mod = 'RLLS')
where, identification parameters are,
na = 1
nb = 2
th = 0

na_ord = [na]
nb_ord = [[nb]]
nc_ord = [nc]
theta = [[th]]

ts = 5e-6
tfin = 0.02
npts = int(old_div(tfin, ts)) + 1
Time = np.linspace(0, tfin, npts)

I am providing the input signal as, [Usim,,] = fset.GBN_seq(npts, 0.08, Range = [-5, 5]), and output signal as,
Yout1, Time, Xsim = cnt.lsim(g_sample11, Usim, Time). Also, I have checked the output from the actual system and the identified system for 4000 points.

Here are the resulting input and output plots, please provide the insight to solve this problem. What am I doing wrong?

[RBdC]

Dear tikull,
thank you very much for your interest in our program SIPPY.

Your use of the program seems appropriate, the model definition seems ok,
but the results are indeed quite strange and the performance appear pretty poor ...

Here below some observations and some questions to investigate your problem:

• 4000 point is a very long time to investigate properly these figure as .png files ..
• the p_switch parameter (the probability of switch), the second input to GBN function, does not seem consistent:
  you say 0.08 in the text, but there is a 0.03 in the title of the figure... maybe, you did not also choose an appropriate value for your system dynamics.. maybe, we have to check ..
• for which of the two callings to function system_identification the results are produced?: 1) LSS or RLSS?
• the model matching seems indeed low.. but how much low? what is the identified Transfer Function?

Please send us your example file; we will run and test and try to give you a better help.

The SIPPY Team

[tikull]

Dear Team SIPPY,
Thank you for the quick response. Please find the code written below. I have not included the noise transfer function for simplicity. Also, there are comparisons made between ARMAX and ARX. Note, changing the LSS and RLSS is not affecting much the identified system in this case. I hope the code attached will clarify the mistake.

`# -- coding: utf-8 --
"""
Created

@author: Giuseppe Armenise, revised by RBdC
example armax mimo
case 3 outputs x 4 inputs

"""
from future import division

from past.utils import old_div

Checking path to access other files

try:
from sippy import *
except ImportError:
import sys, os

sys.path.append(os.pardir)
from sippy import *

import numpy as np
import control.matlab as cnt
from sippy import functionset as fset
from gekko import GEKKO

SISO system

generating transfer functions in z-operator

#var_list = [0.001, 0.002, 0.02]
ts = 5e-6

#Test Resistor and Capacitor network, subscript s stands for serial and p
#stands for parallel component in the network
Rs = 1000 # minimum value for Rs = 1k for the transconductance identification
Rp = 1000000000
Cp = 0.000000001 #minimum value for Cp = 1nF for transconductance identification

#RC network as a system under test, and its Laplace domain transfer
#function, SUT in the time domain and SUT in discrete time domain transfer
#function

num = [1/Rs, 1/(CpRpRs)]
den = [1, (Rp+Rs)/(CpRsRp)]

SISO transfer functions (G and H)

SUT = cnt.tf(num, den)
g_sample11 = cnt.c2d(SUT, ts)
#g_sample11 = cnt.tf(num, den, ts)

#H_sample1 = cnt.tf(H1, den, ts)

time

tfin = 0.02
npts = int(old_div(tfin, ts)) + 1
Time = np.linspace(0, tfin, npts)

#INPUT#
#Usim_noise = np.zeros((1, npts))
[Usim,,] = fset.GBN_seq(npts, 0.08, Range = [-5, 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 = cnt.lsim(H_sample1, err_inputH[0, :], Time)

Noise-free output

Yout1, Time, Xsim = cnt.lsim(g_sample11, Usim, Time)

Total output

#Ytot1 = Yout1 + err_outputH1

#Ytot = Ytot1.squeeze()

identification parameters

na = 1
nb = 2
nc = 2
th = 0

na_ord = [na]
nb_ord = [[nb]]
nc_ord = [nc]
theta = [[th]]

IDENTIFICATION STAGE

#identification with gekko

#m = GEKKO()
#yp,p,K = m.sysid(Time,Usim,Yout1,na,nb,pred='meas')

mode = 'FIXED'
#mode = 'IC'

if mode == 'IC':
# use Information criterion
Id_ARX = system_identification(Yout1, Usim, 'ARX', IC ='AIC', tsample = ts, na_ord = [0, 0], nb_ord = [2, 2], delays = [0, 0], ARX_mod = 'LLS') #
Id_ARMAXi = system_identification(Yout1, Usim, 'ARMAX', IC='AIC', tsample = ts, na_ord=[1, 1], nb_ord=[2, 2],
nc_ord=[2, 2], delays=[0, 0], max_iterations=300, ARMAX_mod = 'ILLS')

elif mode == 'FIXED':
Id_ARX = system_identification(Yout1, Usim, 'ARX', tsample = ts, ARX_orders = [na_ord, nb_ord, theta], ARX_mod = 'RLLS')
Id_ARMAXi = system_identification(Yout1, Usim, 'ARMAX', tsample = ts, ARMAX_orders = [na_ord, nb_ord, nc_ord, theta],
max_iterations = 300, ARMAX_mod = 'ILLS')

output of the identified model

Yout_ARX = Id_ARX.Yid.T
Yout_ARMAXi = Id_ARMAXi.Yid.T
#Yout_ARMAXr = Id_ARMAXr.Yid.T

#Generate Fast fourier transform of the identified data object

#Datf = cnt.fft(ID_ARX)
######plot

import matplotlib.pyplot as plt

plt.close('all')
plt.figure(0)
#plt.subplot(1, 1, 1)
plt.plot(Time, Usim)
plt.grid()
plt.ylabel("Input 1 - GBN")
plt.xlabel("Time")
plt.title("Input (Switch probability=0.03)")
#plt.show()

#plt.close('all')
plt.figure(1)
#plt.subplot(2, 1, 1)
plt.plot(Time, Yout1)
plt.plot(Time, Yout_ARX)
#plt.plot(Time, Yout_ARMAXi)
#plt.plot(Time, yp)
#plt.plot(Time, Yout_ARMAXr)
plt.ylabel("y$_1$,out")
plt.grid()
plt.xlabel("Time")
plt.title("identification data")
plt.legend(['System', 'ARX', 'ARMAXi'])
#plt.legend(['System', 'GekkoARX'])

plt.show()
`

[RBdC]

Dear tikull,
thank you for your file.

We managed to run your example and reproduce your results.
ARX models do not actually seem well identified for this application.

We are now investigating the problem in depth, searching for eventual issues in our core codes.
We will let you know updates.
Best
SIPPY staff

[tikull]

Dear Team SIPPY,
Thanks for the acknowledgment. I shall wait for the updates. I hope this will improve the code.

Best
Tikull

[RBdC]

Dear tikull,
sorry for the big delay in this reply (mainly due to summer vacations)..
I hope this may be still useful for your purpose.

We examined your problem in details.
There are actually different issues involved. Here below some comments:

1. your system has the following Continous-Time TF:
   G =
   0.001 s + 0.001

---

s + 1e06

Note that, when it is turned into discrete time domain, you should pay attention to the sampling time.
ts = 5e-6 (i.e, 0.5e-5) is definitely too large for your system dynamics, since you lose the most part of the transient response.
The system has indeed impulsive response, due to its lead-lag (1pole, 1zero) behavior and its gain (around 0), and also approx. 0.1e-4 (1e-5) as settling time.. that is, it is very fast ! (being a RC network)
so that, a smaller sampling time is a better choice, e.g., ts = 1e-7, being ts usually linked to the settling time (1/60 - 1/120);
this makes also the identification much more reliable and easy to perform
(see our simple matlab file "test_on_matlab.m" intended to explain this issue)
2) In addition, GBN sequence does not seem the best input test signal for your particular system dynamics;
since it has a very fast impulsive response, the transient dynamics could be not sufficiently visible (excited).
As an alternative solution, you could consider a sum of sinusoids (at least 2 sins), by setting suitable amplitudes and frequencies.
3) So that, with ts = 1e-7 and a sum of 2 sinusoids as input signal (see the example_file.py in attach),
being the DT TF:
0.001 z - 0.001
Gd = ---------------
z - 0.9048

after identification via SIPPY, we obtain:
0.001104 z - 0.001104
Gd = ---------------
z^2 - 0.894z

which can be considered a good estimate for the parameters,
but there is indeed an issue with orders.
You may neglect this last problem, and simply consider the following identified TF:
0.001104 z - 0.001104
Gd = ---------------
z - 0.894

we have a small bug in the general codes (actually revealed only by your specific system dynamics).
We will keep investigating to fix this and let you know.
All the best
SIPPY team

test_on_matlab.m.zip

example_file.py.zip