Model ARX di Gekko

Sep 15 2020

Apakah ada cara lain untuk memperkenalkan model ARX di GEKKO selain fungsi arx ()?

Inilah alasannya: Saya mencoba mengidentifikasi model sistem sebagai model ARX. Pertama saya mencoba menggunakan sysid () dan axr () (fungsi di GEKKO) untuk mengidentifikasi sistem saya dan kemudian mensimulasikan hasilnya dan melihat apakah hasilnya sesuai keinginan. Saat menggunakan sampel data kecil (10 menit dan 1 jam) identifikasi dengan sysid () sudah baik, tetapi dengan sampel besar (5 jam), hasil identifikasi kurang baik. Jadi saya mencoba mengidentifikasi sistem saya dengan kode yang saya tulis, menggunakan regresi linier dan variabel dependen tertunda untuk mengidentifikasi model ARX (saya menggunakan kumpulan data yang sama untuk sysid () dan kode saya). Masalahnya adalah jika saya menggunakan kode saya untuk mendapatkan parameter a, b dan c untuk kamus p dan kemudian menggunakan kamus ini untuk fungsi arx (p) untuk membuat simulasi, kurva suhu logis tetapi nilai suhu tidak meski hasil prediksinya bagus.

Hasil identifikasi dengan regresi linier lebih baik dibandingkan dengan identifikasi menggunakan sysid ().

Apa yang saya lakukan salah di sini?

Berikut adalah kode yang saya gunakan untuk Regresi Linier:

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_

Hasil simulasi Menggunakan Gekko dan koefisien dari regresi linier: [1]: https://i.stack.imgur.com/B2vnL.png

Kode untuk Simulasi:

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

Jawaban

1 JohnHedengren Sep 15 2020 at 18:38

Anda bisa mendapatkan sysidhasil yang setara jika Anda menggunakan opsi pred='meas'alih-alih default pred='model'dan menggunakan shift='calc'alih-alih default shift='init'. Regresi linier yang Anda lakukan dapat memberikan hasil yang bias sedangkan opsi default sysid()memberikan hasil yang tidak bias karena menggunakan bentuk kesalahan keluaran. Perbedaannya adalah bahwa nilai berikutnya y[k]diprediksi dari nilai model sebelumnya dan bukan dari pengukuran sebelumnya untuk y[k-1]. Saya memverifikasi bahwa prediksi Gekko sudah benar dengan perhitungan Excel cepat dan satu langkah.

Ini adalah respon model yang setara di Gekko tetapi dengan lebih banyak langkah.

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

Berikut adalah cara membangun model tanpa fungsi 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()

Anda dapat membuat model ARX dengan IMODE=1Python jika Anda mengelola nilai deret waktu dengan nama variabel unik di setiap titik waktu. Perhatikan bahwa label MVdan Anda CVbertukar dalam contoh yang Anda posting. Ini CVadalah variabel yang dikontrol dan merupakan nilai prediksi keluaran. Nilai MVtersebut dapat diatur secara manual oleh operator atau disesuaikan oleh pemecah.

Jika Anda melihat ke dalam fungsi sysid, Anda juga akan melihat contoh bagaimana membangun model ARX tanpa bantuan fungsi ARX tetapi untuk kasus multivariasi. Ini lebih rumit jadi saya tidak merekomendasikan menggunakan pendekatan ini.