โมเดล ARX ใน Gekko

Sep 15 2020

มีวิธีอื่นในการแนะนำโมเดล ARX ใน GEKKO นอกเหนือจากฟังก์ชัน arx () หรือไม่

นี่คือเหตุผล: ฉันกำลังพยายามระบุรูปแบบของระบบเป็นแบบจำลอง ARX ก่อนอื่นฉันลองใช้ sysid () และ axr () (ฟังก์ชันใน GEKKO) เพื่อระบุระบบของฉันจากนั้นจำลองผลลัพธ์และดูว่าผลลัพธ์เป็นไปตามที่ต้องการหรือไม่ เมื่อใช้ตัวอย่างข้อมูลขนาดเล็ก (10 นาทีและ 1 ชม.) การระบุด้วย sysid () เป็นสิ่งที่ดี แต่สำหรับตัวอย่างขนาดใหญ่ (5 ชม.) ผลการระบุตัวตนไม่ดีนัก ดังนั้นฉันจึงพยายามระบุระบบของฉันด้วยโค้ดที่ฉันเขียนโดยใช้การถดถอยเชิงเส้นและตัวแปรตามแบบล่าช้าเพื่อระบุโมเดล ARX (ฉันใช้ชุดข้อมูลเดียวกันสำหรับ sysid () และรหัสของฉัน) ปัญหาคือถ้าฉันใช้รหัสของฉันเพื่อรับพารามิเตอร์ a, b และ c สำหรับพจนานุกรม p แล้วใช้พจนานุกรมนี้สำหรับฟังก์ชัน arx (p) เพื่อสร้างการจำลองเส้นโค้งอุณหภูมิจะเป็นตรรกะ แต่ค่าอุณหภูมิไม่ใช่ แม้จะมีผลการทำนายที่ดี

ผลการระบุด้วยการถดถอยเชิงเส้นดีกว่าการระบุโดยใช้ sysid ()

ฉันทำอะไรผิดที่นี่?

นี่คือรหัสที่ฉันใช้สำหรับ Linear Regression:

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

คำตอบ

1 JohnHedengren Sep 15 2020 at 18:38

คุณจะได้รับเทียบเท่าsysidผลถ้าคุณใช้ตัวเลือกpred='meas'แทนการเริ่มต้นpred='model'และใช้แทนการเริ่มต้นshift='calc' shift='init'การถดถอยเชิงเส้นที่คุณกำลังดำเนินการอาจให้ผลลัพธ์ที่เอนเอียงในขณะที่ตัวเลือกเริ่มต้นsysid()ให้ผลลัพธ์ที่เป็นกลางเนื่องจากใช้รูปแบบข้อผิดพลาดเอาต์พุต ความแตกต่างก็คือว่าต่อไปเป็นที่คาดการณ์จากค่ารุ่นก่อนแทนการวัดก่อนสำหรับy[k] y[k-1]ฉันตรวจสอบแล้วว่าการคาดคะเน Gekko ถูกต้องด้วยการคำนวณ Excel อย่างรวดเร็วและขั้นตอนเดียว

นี่คือการตอบสนองแบบจำลองที่เทียบเท่าใน 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()

คุณสามารถสร้างโมเดล ARX ด้วยIMODE=1ใน Python หากคุณจัดการค่าอนุกรมเวลาด้วยชื่อตัวแปรที่ไม่ซ้ำกันในแต่ละช่วงเวลา โปรดทราบว่าป้ายกำกับMVและCVป้ายกำกับของคุณจะสลับกันในตัวอย่างที่คุณโพสต์ CVเป็นตัวแปรที่ควบคุมและมีมูลค่าการส่งออกที่คาดการณ์ไว้ MVคุ้มค่าที่สามารถปรับเปลี่ยนได้ด้วยตนเองโดยผู้ประกอบการหรือปรับอื่นโดยการแก้

หากคุณดูภายในฟังก์ชัน sysid คุณจะเห็นตัวอย่างวิธีการสร้างโมเดล ARX โดยไม่ต้องใช้ฟังก์ชัน ARX แต่สำหรับกรณีหลายตัวแปร สิ่งนี้ซับซ้อนกว่าดังนั้นฉันไม่แนะนำให้ใช้แนวทางนี้