Modelos ARX em Gekko

Sep 15 2020

Existe alguma outra maneira de introduzir um modelo ARX em GEKKO além da função arx ()?

Aqui está o motivo: estou tentando identificar um modelo de sistema como um modelo ARX. Primeiro tentei usar sysid () e axr () (funções em GEKKO) para identificar meu sistema e depois simular o resultado e ver se a saída é a desejada. Ao usar pequenas amostras de dados (10 min e 1h) a identificação com sysid () foi boa, mas com uma grande amostra (5h), os resultados da identificação não foram tão bons. Então, tentei identificar meu sistema com o código que escrevi, usando regressão linear e uma variável dependente atrasada para identificar um modelo ARX (usei o mesmo conjunto de dados para sysid () e meu código). O problema é que se eu usar meu código para obter os parâmetros a, bec para o dicionário p e, em seguida, usar este dicionário para a função arx (p) para criar uma simulação, a curva de temperatura é lógica, mas os valores de temperatura não são apesar dos bons resultados de previsão.

Os resultados da identificação com a regressão Linear são melhores do que a identificação usando sysid ().

O que eu estou fazendo errado aqui?

Aqui está o código que usei para a regressão linear:

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 da simulação usando Gekko e coeficiente de regressão linear: [1]: https://i.stack.imgur.com/B2vnL.png

Código para simulação:

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

Respostas

1 JohnHedengren Sep 15 2020 at 18:38

Você pode obter sysidresultados equivalentes se usar a opção em pred='meas'vez do padrão pred='model'e usar em shift='calc'vez do padrão shift='init'. A regressão linear que você está executando pode fornecer resultados tendenciosos, enquanto as opções padrão sysid()fornecem resultados imparciais porque usa um formulário de erro de saída. A diferença é que o próximo y[k]é previsto a partir do valor do modelo anterior em vez da medição anterior para y[k-1]. Eu verifiquei que as previsões Gekko estão corretas com um cálculo rápido do Excel e uma única etapa.

Aqui está uma resposta de modelo equivalente em Gekko, mas com mais etapas.

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

Aqui está uma maneira de construir o modelo sem a função 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()

Você pode construir o modelo ARX IMODE=1em Python se gerenciar os valores da série temporal com nomes de variáveis exclusivos em cada ponto no tempo. Observe que seus rótulos MVe CVsão trocados no exemplo que você postou. O CVé a variável controlada e é o valor previsto de saída. O MVé aquele valor que pode ser ajustado manualmente pelo operador ou então ajustado pelo solucionador.

Se você olhar dentro da função sysid, você também verá um exemplo de como construir o modelo ARX sem a ajuda da função ARX, mas para o caso multivariado. Isso é mais complicado, então não recomendo usar essa abordagem.