Gekko의 ARX 모델
GEKKO에서 arx () 함수 외에 ARX 모델을 도입하는 다른 방법이 있습니까?
이유는 다음과 같습니다. 시스템 모델을 ARX 모델로 식별하려고합니다. 먼저 sysid () 및 axr () (GEKKO의 함수)를 사용하여 내 시스템을 식별 한 다음 결과를 시뮬레이션하고 출력이 원하는지 확인했습니다. 작은 데이터 샘플 (10 분 및 1 시간)을 사용하는 경우 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()
답변
당신은 해당 얻을 수 있습니다 sysid당신이 옵션을 사용하면 결과를 pred='meas'대신 기본 pred='model'및 사용 shift='calc'대신에 기본적으로 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')