将matlab系统识别转换为python

Question

我在 Matlab 中有一段代码，我想将其转换为 Python。 Matlab 代码使用此处提供的系统识别工具箱：

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);

为了提供有关代码的更多解释，我将数据封装为 iddata，并将其拆分为训练和验证数据。这些拆分将用于估计识别线性 ARX 模型的最佳顺序。一旦确定，我想用这些顺序检测系统中的非线性。然后，我想构建 ARX 模型，找到初始状态，并将其转换为稳态模型。最后，我想检测任何异常行为以识别系统是否不稳定。

我开始转换到 Python 并找到了一个名为 SIPPY 的包，用于线性 ARX mdeoling。这是我用 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)

此代码将读取数据（无需 iddata 封装）并输出给定​​订单范围的最佳 ARX 模型。这意味着它将作为arxstruc(z1,z2,struc(0:2, 1:50,1:50))、nn = selstruc(V,'aic'); 和sys=arx(z1,nn); 执行。但是，在对相同数据进行测试以比较输出时，我发现 Matlab 给出的最佳命令是 [1 25 1] 而 python 返回 [2 35 1]。当我调查原因时，我发现损失值与 Matlab 和 Python 不同，并且由于输出将是实现最小损失的顺序，因此具有不同的顺序是合乎逻辑的。那么有人可以帮我解决这个问题吗？ Matlab中使用的损失函数是什么？是否有在 Matlab 中模拟系统识别并在 Python 中提供相同结果的包？

Answer 1

这里有一些代码在另一个 StackOverflow 帖子上 Model Predictive Control 和系统标识。

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()