Modelos ARX en Gekko

Sep 15 2020

¿Hay alguna otra forma de introducir un modelo ARX en GEKKO que la función arx ()?

Esta es la razón: estoy tratando de identificar un modelo de un sistema como un modelo ARX. Primero intenté usar sysid () y axr () (funciones en GEKKO) para identificar mi sistema y luego simular el resultado y ver si la salida es la deseada. Cuando se utilizaron muestras de datos pequeñas (10 min y 1 h), la identificación con sysid () fue buena, pero con una muestra grande (5 h), los resultados de identificación no fueron tan buenos. Entonces traté de identificar mi sistema con el código que escribí, usando regresión lineal y una variable dependiente retrasada para identificar un modelo ARX (usé el mismo conjunto de datos para sysid () y mi código). El problema es que si usé mi código para obtener los parámetros a, byc para el diccionario py luego usé este diccionario para la función arx (p) para crear una simulación, la curva de temperatura es lógica pero los valores de temperatura no lo son a pesar de los buenos resultados de la predicción.

Los resultados de la identificación con la regresión lineal son mejores que la identificación con sysid ().

¿Qué estoy haciendo mal aquí?

Aquí está el código que utilicé para la regresión lineal:

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_

Resultados de la simulación Usando Gekko y coeficiente de regresión lineal: [1]: https://i.stack.imgur.com/B2vnL.png

Código para simulació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


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

Respuestas

1 JohnHedengren Sep 15 2020 at 18:38

Puede obtener sysidresultados equivalentes si usa la opción en pred='meas'lugar de la predeterminada pred='model'y la usa en shift='calc'lugar de la predeterminada shift='init'. La regresión lineal que está realizando puede dar resultados sesgados, mientras que las opciones predeterminadas sysid()dan resultados no sesgados porque utiliza un formulario de error de salida. La diferencia es que el siguiente y[k]se predice a partir del valor del modelo anterior en lugar de la medición anterior para y[k-1]. Verifiqué que las predicciones de Gekko son correctas con un cálculo rápido de Excel y un solo paso.

Aquí hay una respuesta de modelo equivalente en Gekko pero con más pasos.

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

Aquí hay una forma de construir el modelo sin la función 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()

Puede construir el modelo ARX IMODE=1en Python si administra los valores de la serie temporal con nombres de variable únicos en cada punto de tiempo. Tenga en cuenta que sus etiquetas MVy CVse intercambian en el ejemplo que publicó. El CVes la variable controlada y es el valor predicho de salida. El MVes ese valor que puede ser ajustado manualmente por el operador o bien ajustado por el solucionador.

Si observa dentro de la función sysid, también verá un ejemplo de cómo construir el modelo ARX sin la ayuda de la función ARX pero para el caso multivariante. Esto es más complicado, por lo que no recomiendo usar este enfoque.