Modèles ARX dans Gekko

Sep 15 2020

Existe-t-il un autre moyen d'introduire un modèle ARX dans GEKKO que la fonction arx ()?

Voici la raison: j'essaie d'identifier un modèle de système comme un modèle ARX. J'ai d'abord essayé d'utiliser sysid () et axr () (fonctions dans GEKKO) pour identifier mon système, puis simuler le résultat et voir si la sortie est comme vous le souhaitez. Lors de l'utilisation de petits échantillons de données (10 min et 1h), l'identification avec sysid () était bonne, mais avec un grand échantillon (5h), les résultats d'identification n'étaient pas aussi bons. J'ai donc essayé d'identifier mon système avec le code que j'avais écrit, en utilisant la régression linéaire et une variable dépendante retardée pour identifier un modèle ARX (j'ai utilisé le même jeu de données pour sysid () et mon code). Le problème est que si j'ai utilisé mon code pour obtenir les paramètres a, b et c du dictionnaire p, puis utilisé ce dictionnaire pour la fonction arx (p) pour créer une simulation, la courbe de température est logique mais les valeurs de température ne le sont pas malgré les bons résultats de prédiction.

Les résultats d'identification avec la régression linéaire sont meilleurs que l'identification utilisant sysid ().

Qu'est-ce que je fais de mal ici?

Voici le code que j'ai utilisé pour la régression linéaire:

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_

Résultats de la simulation à l'aide de Gekko et du coefficient de régression linéaire: [1]: https://i.stack.imgur.com/B2vnL.png

Code pour la simulation:

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

Réponses

1 JohnHedengren Sep 15 2020 at 18:38

Vous pouvez obtenir des sysidrésultats équivalents si vous utilisez l'option pred='meas'au lieu de la valeur par défaut pred='model'et utilisez à la shift='calc'place de la valeur par défaut shift='init'. La régression linéaire que vous effectuez peut donner des résultats biaisés tandis que les options par défaut de sysid()donnent des résultats non biaisés car elle utilise un formulaire d'erreur de sortie. La différence est que le suivant y[k]est prédit à partir de la valeur du modèle précédente au lieu de la mesure précédente pour y[k-1]. J'ai vérifié que les prédictions Gekko sont correctes avec un calcul Excel rapide et une seule étape.

Voici une réponse de modèle équivalente dans Gekko mais avec plus d'étapes.

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

Voici un moyen de créer le modèle sans la fonction 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()

Vous pouvez créer le modèle ARX avec IMODE=1en Python si vous gérez les valeurs de séries temporelles avec des noms de variables uniques à chaque instant. Notez que vos étiquettes MVet CVsont permutées dans l'exemple que vous avez publié. Le CVest la variable contrôlée et est la valeur prédite de sortie. C'est MVcette valeur qui peut être ajustée manuellement par l'opérateur ou bien ajustée par le solveur.

Si vous regardez à l'intérieur de la fonction sysid, vous verrez également un exemple de la façon de construire le modèle ARX sans l'aide de la fonction ARX mais pour le cas multivarié. C'est plus compliqué, je ne recommande donc pas d'utiliser cette approche.