python - Gekko 中的 ARX 模型

Question

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

原因如下：我试图将系统模型识别为 ARX 模型。首先，我尝试使用 sysid() 和 axr()（GEKKO 中的函数）来识别我的系统，然后模拟结果并查看输出是否符合要求。当使用小数据样本（10 分钟和 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()

标签： pythonpandasgekko

Answer 1

sysid如果您使用 optionpred='meas'而不是 defaultpred='model'和 useshift='calc'而不是default ，您可以获得等效的结果shift='init'。您正在执行的线性回归可能会给出有偏的结果，而默认选项会sysid()给出无偏的结果，因为它使用输出错误形式。不同之处在于下一个y[k]是从先前的模型值而不是先前的测量值预测的y[k-1]。我通过快速 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()

IMODE=1如果您在每个时间点使用唯一变量名称管理时间序列值，则可以使用 Python 构建 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')