2

I have a piece of code in Matlab that I want to convert to Python. The Matlab code is using the system identification toolbox which is provided here:

Ts = 1;
Znl=iddata(Xdati(:,2),Xdati(:,1),Ts);

z=iddata(Xdati(:,1),Xdati(:,2),Ts); 
z1=z(1:floor(length(Xdati(:,1))/2));  
z2=z(floor(length(Xdati(:,1))/2)+1:1:floor(2*length(Xdati(:,1))/2));  

V = arxstruc(z1,z2,struc(0:2, 1:50,1:50)); % Find the best structure of ARX model that can be 
with degrees between 1 and 50.
nn = selstruc(V,'aic'); 
[NLHyp,NLValue,NLRegs,NoiseSigma,DetectRatio] = isnlarx(Znl,nn); 

if 2*max(nn)<length(z1.y)

sys=arx(z1,nn);
x0=findstates(sys,z);
ssmodel=idss(sys);
Unstable_System=[];
Unstable_System=find(abs(eig(ssmodel.A))>1);

To provide more explanation about the code, I have data where I encapsulate it as iddata, and split it to train and validation data. these splits will be used in order to estimate the best order to identify a linear ARX model. Once identified, I want to detect nonlinearity in the system with these orders. Then, I want to construct the ARX model, find the initial states, and convert it to a steady-state model. Finally, I want to detect any abnormal behavior to identify if the system was unstable.

I started the conversion to Python and found a package called SIPPY for linear ARX mdeoling. Here is the code I wrote in Python:

T = pd.read_excel('./test_data.xlsx')
input_0 = np.array(T.iloc[:, 0])
output_0 = np.array(T.iloc[:, 1])

loss = []
na = list(range(0, 3))
nb = list(range(1, 51))
nk = list(range(1, 51))

final_model = system_identification(output_0, input_0, 'ARX', IC='AIC', na_ord=na, nb_ord=nb, delays=nk)

print(final_model.G)
print(final_model.Vn)
print(final_model.Yid)

This code will read the data (no need for iddata encapsulation) and output the best ARX model for a given range of orders. This means it will perform as arxstruc(z1,z2,struc(0:2, 1:50,1:50)), nn = selstruc(V,'aic');, and sys=arx(z1,nn);. However, when testing both on the same data to compare the output, I found that the best orders given by Matlab were [1 25 1] While python returns [2 35 1]. When I investigated the reason I found that the loss value is different from Matlab than Python, and since the output would be the order that achieves the minimum loss, it is logical to have different orders. So can anyone help me with that problem? What is the loss function used in Matlab? and is there a package that simulates the system identification in Matlab and provide the same results in Python?

abdullatif alnajim
  • 33
  • 3
  • 1
    Here is a link to the `sysid` (System Identification) function in Python Gekko: https://apmonitor.com/do/index.php/Main/ModelIdentification Once a model is identified, Python Gekko is can use it to create a Model Predictive Controller https://apmonitor.com/pds/index.php/Main/ARXTimeSeries You may also be interested in the Seeq add-on that uses Gekko: https://github.com/BYU-PRISM/Seeq – John Hedengren Dec 20 '21 at 02:59
  • 1
    @JohnHedengren I will read about it, thank you – abdullatif alnajim Dec 20 '21 at 12:12

1 Answers1

0

Here is some code that is on another StackOverflow post on Model Predictive Control and system identification.

system identification

from gekko import GEKKO
import pandas as pd
import matplotlib.pyplot as plt

# load data and parse into columns
url = 'http://apmonitor.com/do/uploads/Main/tclab_dyn_data2.txt'
data = pd.read_csv(url)
t = data['Time']
u = data[['H1','H2']]
y = data[['T1','T2']]

# generate time-series model
m = GEKKO(remote=False) # remote=True for MacOS

# system identification
na = 2 # output coefficients
nb = 2 # input coefficients
yp,p,K = m.sysid(t,u,y,na,nb,diaglevel=1)

plt.figure()
plt.subplot(2,1,1)
plt.plot(t,u)
plt.legend([r'$u_0$',r'$u_1$'])
plt.ylabel('MVs')
plt.subplot(2,1,2)
plt.plot(t,y)
plt.plot(t,yp)
plt.legend([r'$y_0$',r'$y_1$',r'$z_0$',r'$z_1$'])
plt.ylabel('CVs')
plt.xlabel('Time')
plt.savefig('sysid.png')
plt.show()
TexasEngineer
  • 684
  • 6
  • 15