Gekko中的ARX模型

编程问答 2022-06-13

问题描述

除了arx（）函数以外，还有什么其他方法可以在GEKKO中引入ARX模型吗？

这是原因：我试图将系统模型识别为ARX模型。首先，我尝试使用sysid（）和axr（）（GEKKO中的函数）来识别我的系统，然后模拟结果并查看输出是否符合要求。当使用较小的数据样本（10分钟和1小时）时，使用sysid（）进行识别的效果很好，但是对于较大的样本（5小时），识别结果却不太好。因此，我尝试使用线性回归和延迟因变量来识别我编写的代码，以识别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,0)
u4 = np.append(u3,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()

解决方法

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

Step response

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

More steps 

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,0)
u4 = np.append(u3,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,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),0)
u = u5

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

time = np.linspace(0,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,i-1,i,)
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,)
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')
gekkopandaspython

相关问答

matplotlib报错：AttributeError: module 'backend_interagg' has no attribute 'FigureCanvas'. Did you mean: 'FigureCanvasAgg'?
使用本地python环境可以成功执行 import pandas as pd impor...
gitlab登录失败，报错：This challenge page was accidentally cached by an intermediary and is no longer available.
设置时间 控制面板
后端开发常见错误
错误1：Request method ‘DELETE‘ not supported 错误还原：...
docker常见错误
错误1：启动docker镜像时报错：Error response from daemon:...
idea常见错误
错误1：private field ‘xxx‘ is never assigned 按Alt...
pip安装依赖失败
报错如下，通过源不能下载，最后警告pip需升级版本 Requirem...