Gekko中的ARX模型

Question

除了arx()函数之外，有没有其他方法可以在GEKKO中引入ARX模型？

原因是:我试图将系统的一个模型标识为ARX模型。首先，我尝试使用sysid()和axr() (GEKKO中的函数)来识别我的系统，然后模拟结果，看看输出是否符合预期。当使用小样本(10min和1h)时，sysid()的识别效果较好，但在大样本(5h)时，识别结果不是很好。因此，我尝试使用我编写的代码来识别我的系统，使用线性回归和延迟因变量来识别ARX模型(我对sysid()和我的代码使用了相同的数据集)。问题是，如果我使用我的代码为字典p获取a、b和c参数，然后使用这个字典为arx(p)函数创建模拟，那么温度曲线是符合逻辑的，但是尽管预测结果很好，但温度值却不是。

线性回归的识别结果优于sysid()的识别结果。

我在这里做错了什么？

下面是我用于线性回归的代码：

import sklearn.metrics as metrics
import pandas as pd
import numpy as np
from pandas.plotting import autocorrelation_plot
from sklearn.linear_model import LinearRegression
import seaborn as sns
import matplotlib.pyplot as plt

b_dataframe = pd.read_csv("Temp.txt")
b_dataframe.columns = ["Temp"]
a_dataframe = pd.read_csv("State.txt")
a_dataframe.columns = ["State"]
df = b_dataframe.join(a_dataframe)

# autocorrelation_plot(df["T[C]"])
X = df.drop("Temp", axis=1) # Drop column T[U]
X.loc[:, "lagged_T_1"] = df["Temp"].shift(1).fillna(0)
#X.loc[:, "lagged_T_2"] = df["T[C]"].shift(2).fillna(0)
y = df["Temp"]
[![enter image description here][1]][1]
#defined a function for linear regression
lin_reg = LinearRegression()
# Train data points --> the rest is for prediction.
n_train = 2500
# just a split 
x_train, x_test = X.iloc[:n_train,:], X.iloc[n_train:,:]
y_train, y_test = y.iloc[:n_train], y.iloc[n_train:]

#model fitting/ train.
#Fit x, y values used for train to the given data.
lin_reg.fit(x_train.values,y_train.values)

# test: With the rest of data points, test the results of the prediction.
y_pred = pd.Series(lin_reg.predict(x_test.values), name="T_pred")
print(lin_reg.coef_)

plt.plot(y_pred.values)
plt.plot(y_test.values)
#plt.text(1, 1, metrics.mean_absolute_error(y_test, y_pred))
plt.legend(["Prediction", "Actual"])
plt.ylim([11.6, 15])
lin_reg.coef_, lin_reg.intercept_

使用Gekko和线性回归系数的模拟结果: 1：https://i.stack.imgur.com/B2vnL.png

模拟代码：

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

na = 1# Number of A coefficients
nb = 1 # Number of B coefficients
ny = 1 # Number of outputs
nu = 1 # Number of inputs

# A (na x ny)
# actual A,B,C values are from 5 h data
A = np.array([[0.960187147]])
# B (ny x (nb x nu))
B= np.array([[[-0.000361506092]]])
C = np.array([ 0.565842747871903])

# create parameter dictionary
# parameter dictionary p['a'], p['b'], p['c']
# a (coefficients for a polynomial, na x ny)
# b (coefficients for b polynomial, ny x (nb x nu))
# c (coefficients for output bias, ny)
p = {'a':A,'b':B,'c':C}


m = GEKKO(remote=True)
y,u = m.arx(p)

# load inputs
#tf = 719 # final time
u1 = np.append(np.ones(500),np.zeros(500),0)
u2 = np.append(u1, np.ones(500),0)
u3 = np.append(u2, np.zeros(500),0)
u4 = np.append(u3, np.ones(500),0)
u5 = np.append(u4, np.zeros(936),0)
u[0].value = u5


mv = y[0]
cv= u[0]
mv.value = 14.2

m.time = np.linspace(0,3436,3436)
m.options.imode = 4
m.options.nodes= 2
#m.options.SOLVER = 1
# simulate
m.solve()

Answer 1

如果使用选项pred='meas'而不是默认的pred='model'，并使用shift='calc'而不是默认的shift='init'，则可以获得等效的shift='calc'结果。您正在执行的线性回归可能会给出有偏差的结果，而sysid()中的默认选项则会给出无偏差的结果，因为它使用的是输出错误形式。不同之处在于，根据先前的模型值而不是y[k-1]的先前测量来预测下一个y[k]。我通过一个快速的Excel计算和一个步骤验证了Gekko的预测是正确的。

​

​

以下是Gekko中的等效模型响应，但步骤更多。

​

​

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

na = 1# Number of A coefficients
nb = 1 # Number of B coefficients
ny = 1 # Number of outputs
nu = 1 # Number of inputs

# A (na x ny)
# actual A,B,C values are from 5 h data
A = np.array([[0.960187147]])
# B (ny x (nb x nu))
B= np.array([[[-0.000361506092]]])
C = np.array([ 0.565842747871903])

# create parameter dictionary
# parameter dictionary p['a'], p['b'], p['c']
# a (coefficients for a polynomial, na x ny)
# b (coefficients for b polynomial, ny x (nb x nu))
# c (coefficients for output bias, ny)
p = {'a':A,'b':B,'c':C}


m = GEKKO(remote=True)
y,u = m.arx(p)

# load inputs
#tf = 719 # final time
u1 = np.append(np.ones(500),np.zeros(500),0)
u2 = np.append(u1, np.ones(500),0)
u3 = np.append(u2, np.zeros(500),0)
u4 = np.append(u3, np.ones(500),0)
u5 = np.append(u4, np.zeros(936),0)
u[0].value = u5

cv = y[0]
mv= u[0]
cv.value = 14.2

# for time steps of 1 use final time of 3435
m.time = np.linspace(0,3435,3436)
m.options.imode = 4
m.options.nodes= 2
#m.options.SOLVER = 1
# simulate
m.solve()

plt.subplot(2,1,1)
plt.plot(m.time,cv.value,'b-',label='CV')
plt.legend()
plt.subplot(2,1,2)
plt.plot(m.time,mv.value,'r--',label='MV')
plt.legend()

plt.show()

下面是一种不使用ARX函数构建模型的方法：

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

A = 0.960187147
B = -0.000361506092
C = 0.565842747871903

m = GEKKO(remote=True)

u1 = np.append(np.ones(500),np.zeros(500),0)
u2 = np.append(u1, np.ones(500),0)
u3 = np.append(u2, np.zeros(500),0)
u4 = np.append(u3, np.ones(500),0)
u5 = np.append(u4, np.zeros(936),0)
u = u5

cv = m.Array(m.Var,3436)

time = np.linspace(0,3435,3436)
m.options.imode = 1

m.Equation(cv[0]==14.2)
for i in range(3435):
    m.Equation(cv[i+1] == A * cv[i] + B * u[i] + C)

# simulate
m.solve()

如果您在每个时间点使用惟一的变量名管理时间序列值，则可以在Python中使用IMODE=1构建ARX模型。注意，在您发布的示例中，您的MV和CV标签被交换了。CV是控制变量，是输出预测值。MV是可以由操作员手动调整或由求解器调整值。

如果您深入研究sysid函数，您还将看到一个示例，说明如何在不使用ARX函数的情况下构建ARX模型，但适用于多变量情况。这比较复杂，所以我不推荐使用这种方法。

syid.Raw('Objects')
syid.Raw('  sum_a[1:ny] = sum(%i)'%na)
syid.Raw('  sum_b[1:ny][1::nu] = sum(%i)'%nbk)
syid.Raw('End Objects')
syid.Raw('  ')
syid.Raw('Connections')
syid.Raw('  a[1:na][1::ny] = sum_a[1::ny].x[1:na]')
syid.Raw('  b[1:nb][1::nu][1:::ny] = sum_b[1:::ny][1::nu].x[1:nb]')
syid.Raw('  sum_a[1:ny] = sum_a[1:ny].y')
syid.Raw('  sum_b[1:ny][1::nu] = sum_b[1:ny][1::nu].y')
syid.Raw('End Connections')
syid.Raw('  ')
syid.Raw('Constants')
syid.Raw('  n = %i' %n)
syid.Raw('  nu = %i'%nu)
syid.Raw('  ny = %i'%ny)
syid.Raw('  na = %i'%na)
syid.Raw('  nb = %i'%nbk)
syid.Raw('  m = %i'%m)
syid.Raw('  ')
syid.Raw('Parameters')
syid.Raw('  a[1:na][1::ny] = 0.9 !>= 0.00001 <= 0.9999999')
syid.Raw('  b[1:nb][1::nu][1:::ny] = 0')
syid.Raw('  c[1:ny] = 0')
syid.Raw('  u[1:n][1::nu]')
syid.Raw('  y[1:m][1::ny]')
syid.Raw('  z[1:n][1::ny]')
syid.Raw('  Ks[1:ny][1::nu] = 1')
syid.Raw('  ')
syid.Raw('Variables')
syid.Raw('  y[m+1:n][1::ny] = 0')
syid.Raw('  sum_a[1:ny] = 0 !<= 1')
syid.Raw('  sum_b[1:ny][1::nu] = 0')
syid.Raw('  K[1:ny][1::nu] = 0 >=-1e8 <=1e8')
syid.Raw('  ')
syid.Raw('Equations')
if pred=='model':
    # use model to predict next y (Output error)
    eqn = '  y[m+1:n][1::ny] = a[1][1::ny]*y[m:n-1][1::ny]'
else:
    # use measurement to predict next y (ARX)
    eqn = '  y[m+1:n][1::ny] = a[1][1::ny]*z[m:n-1][1::ny]'
for j in range(1,nu+1):
    eqn += '+b[1][%i][1::ny]*u[m:n-1][%i]'%(j,j,)
    for i in range(2,nbk+1): 
        eqn += '+b[%i][%i][1::ny]*u[m-%i:n-%i][%i]'%(i,j,i-1,i,j,)
if pred=='model':
    # use model to predict next y (Output error)
    seqn = '+a[%i][1::ny]*y[m-%i:n-%i][1::ny]'
else:
    # use measurement to predict next y (ARX)
    seqn = '+a[%i][1::ny]*z[m-%i:n-%i][1::ny]'        
for i in range(2,na+1): 
    eqn += seqn%(i,i-1,i,)
eqn += '+c[1::ny]'
syid.Raw(eqn)
syid.Raw('')
syid.Raw('  K[1:ny][1::nu] * (1 - sum_a[1:ny]) = Ks[1:ny][1::nu] * sum_b[1:ny][1::nu]')        
syid.Raw('  minimize %e * (y[m+1:n][1::ny] - z[m+1:n][1::ny])^2'%objf)
syid.Raw('  minimize 1e-3 * a[1:na][1::ny]^2')
syid.Raw('  minimize 1e-3 * b[1:nb][1::nu][1:::ny]^2')
syid.Raw('  minimize 1e-3 * c[1:ny]^2')