Mô hình ARX ở Gekko

Sep 15 2020

Có cách nào khác để giới thiệu mô hình ARX trong GEKKO ngoài hàm arx () không?

Đây là lý do: Tôi đang cố gắng xác định một mô hình của hệ thống là một mô hình ARX. Đầu tiên, tôi đã thử sử dụng sysid () và axr () (các hàm trong GEKKO) để xác định hệ thống của mình và sau đó mô phỏng kết quả và xem liệu đầu ra có như mong muốn hay không. Khi sử dụng các mẫu dữ liệu nhỏ (10 phút và 1 giờ), việc nhận dạng với sysid () là tốt, nhưng với mẫu lớn (5h), kết quả nhận dạng không tốt như vậy. Vì vậy, tôi đã cố gắng xác định hệ thống của mình bằng mã tôi đã viết, sử dụng hồi quy tuyến tính và một biến phụ thuộc trễ để xác định mô hình ARX (tôi đã sử dụng cùng một tập dữ liệu cho sysid () và mã của mình). Vấn đề là nếu tôi sử dụng mã của mình để lấy các tham số a, b và c cho từ điển p và sau đó sử dụng từ điển này cho hàm arx (p) để tạo mô phỏng, thì đường cong nhiệt độ là logic nhưng giá trị nhiệt độ thì không. mặc dù kết quả dự đoán tốt.

Kết quả nhận dạng với hồi quy tuyến tính tốt hơn so với nhận dạng bằng sysid ().

Tôi làm gì sai ở đây?

Đây là mã tôi đã sử dụng cho Hồi quy tuyến tính:

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_

Kết quả mô phỏng Sử dụng Gekko và hệ số từ hồi quy tuyến tính: [1]: https://i.stack.imgur.com/B2vnL.png

Mã cho mô phỏng:

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

Trả lời

1 JohnHedengren Sep 15 2020 at 18:38

Bạn có thể nhận được sysidkết quả tương đương nếu bạn sử dụng tùy chọn pred='meas'thay vì mặc định pred='model'và sử dụng shift='calc'thay vì mặc định shift='init'. Hồi quy tuyến tính mà bạn đang thực hiện có thể cho kết quả sai lệch trong khi các tùy chọn mặc định sysid()đưa ra kết quả không thiên vị vì nó sử dụng dạng lỗi đầu ra. Sự khác biệt là y[k]giá trị tiếp theo được dự đoán từ giá trị mô hình trước thay vì phép đo trước đó cho y[k-1]. Tôi đã xác minh rằng các dự đoán của Gekko là chính xác bằng một phép tính Excel nhanh và một bước duy nhất.

Đây là một phản ứng mô hình tương đương trong Gekko nhưng với nhiều bước hơn.

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

Đây là một cách để xây dựng mô hình mà không cần chức năng 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()

Bạn có thể xây dựng mô hình ARX bằng IMODE=1Python nếu bạn quản lý các giá trị chuỗi thời gian với các tên biến duy nhất tại mỗi thời điểm. Lưu ý rằng nhãn MVvà CVnhãn của bạn được hoán đổi trong ví dụ mà bạn đã đăng. Là CVbiến được kiểm soát và là giá trị dự đoán đầu ra. Là MVgiá trị đó có thể được điều chỉnh bằng tay bởi người vận hành hoặc điều chỉnh bởi người giải quyết khác.

Nếu bạn nhìn vào bên trong hàm sysid, bạn cũng sẽ thấy một ví dụ về cách xây dựng mô hình ARX mà không cần đến sự trợ giúp của hàm ARX nhưng đối với trường hợp đa biến. Điều này phức tạp hơn nên tôi không khuyên bạn nên sử dụng cách tiếp cận này.