Modele ARX w Gekko

Sep 15 2020

Czy istnieje inny sposób wprowadzenia modelu ARX w GEKKO niż funkcja arx ()?

Oto powód: próbuję zidentyfikować model systemu jako model ARX. Najpierw próbowałem użyć sysid () i axr () (funkcje w GEKKO), aby zidentyfikować mój system, a następnie zasymulować wynik i sprawdzić, czy wyjście jest zgodne z oczekiwaniami. Przy użyciu małych próbek danych (10 min i 1 godz.) Identyfikacja za pomocą sysid () była dobra, ale przy dużej próbce (5 godz.) Wyniki identyfikacji nie były tak dobre. Próbowałem więc zidentyfikować mój system za pomocą napisanego przez siebie kodu, używając regresji liniowej i opóźnionej zmiennej zależnej do identyfikacji modelu ARX (użyłem tego samego zestawu danych dla sysid () i mojego kodu). Problem polega na tym, że gdybym użył mojego kodu do uzyskania parametrów a, b i c dla słownika p, a następnie użył tego słownika dla funkcji arx (p) do stworzenia symulacji, krzywa temperatury jest logiczna, ale wartości temperatury nie są pomimo dobrych wyników prognoz.

Wyniki identyfikacji za pomocą regresji liniowej są lepsze niż identyfikacja za pomocą sysid ().

Co ja tu robię źle?

Oto kod, którego użyłem do regresji liniowej:

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_

Wyniki symulacji przy użyciu Gekko i współczynnika z regresji liniowej: [1]: https://i.stack.imgur.com/B2vnL.png

Kod do symulacji:

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

Odpowiedzi

1 JohnHedengren Sep 15 2020 at 18:38

Możesz uzyskać równoważne sysidwyniki, jeśli użyjesz opcji pred='meas'zamiast domyślnej pred='model'i użyjesz shift='calc'zamiast domyślnej shift='init'. Regresja liniowa, którą wykonujesz, może dawać wyniki obciążone, podczas gdy domyślne opcje sysid()dają obiektywne wyniki, ponieważ używa formy błędu wyjściowego. Różnica polega na tym, że następny y[k]jest prognozowany na podstawie wartości poprzedniego modelu, a nie poprzedniego pomiaru dla y[k-1]. Sprawdziłem, czy prognozy Gekko są poprawne, wykonując szybkie obliczenia w programie Excel i wykonując jeden krok.

Oto równoważna odpowiedź modelu w Gekko, ale z większą liczbą kroków.

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

Oto sposób na zbudowanie modelu bez funkcji 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()

Możesz zbudować model ARX IMODE=1w Pythonie, jeśli zarządzasz wartościami szeregów czasowych z unikalnymi nazwami zmiennych w każdym punkcie czasowym. Zwróć uwagę, że w opublikowanym przykładzie Twoje MVi CVetykiety są zamienione. Jest CVto zmienna kontrolowana i jest przewidywaną wartością wyjściową. Jest MVto wartość, która może być dostosowywana ręcznie przez operatora lub przez solvera.

Jeśli zajrzysz do wnętrza funkcji sysid, zobaczysz również przykład, jak zbudować model ARX bez pomocy funkcji ARX, ale dla przypadku wielowymiarowego. Jest to bardziej skomplikowane, więc nie polecam stosowania tego podejścia.