Блог

Есть ли какой-либо другой способ ввести модель ARX в GEKKO, кроме функции arx()?

Вот причина: я пытаюсь идентифицировать модель системы как модель ARX. Сначала я попытался использовать sysid() и axr() (функции в GEKKO) для идентификации моей системы, а затем смоделировать результат и посмотреть, соответствует ли результат желаемому. При использовании небольших выборок данных (10 минут и 1 час) идентификация с помощью sysid() была хорошей, но при большой выборке (5 часов) результаты идентификации были не такими хорошими. Итак, я попытался идентифицировать свою систему с помощью кода, который я написал, используя линейную регрессию и отложенную зависимую переменную для идентификации модели ARX (я использовал тот же набор данных для sysid() и моего кода). Проблема в том, что если я использовал свой код для получения параметров a, b и c для словаря p, а затем использовал этот словарь для функции arx (p) для создания симуляции, температурная кривая логична, но значения температуры не соответствуют хорошим результатам прогнозирования.

Результаты идентификации с линейной регрессией лучше, чем идентификация с использованием sysid().

Что я здесь делаю не так?

Вот код, который я использовал для линейной регрессии:

 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_
  

Результаты моделирования с использованием Gekko и коэффициента линейной регрессии:
[1]:https://i.stack.imgur.com/B2vnL.png

Код для моделирования:

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

Вы можете получить эквивалентные sysid результаты, если будете использовать опцию pred='meas' вместо значения по умолчанию pred='model' и использовать shift='calc' вместо значения по умолчанию shift='init' . Выполняемая вами линейная регрессия может давать предвзятые результаты, в то время как параметры по умолчанию в sysid() дают непредвзятые результаты, поскольку она использует форму ошибки вывода. Разница в том, что следующее значение y[k] прогнозируется на основе предыдущего значения модели вместо предыдущего измерения для y[k-1] . Я проверил правильность прогнозов Gekko с помощью быстрого вычисления в Excel и одного шага.

Вот эквивалентный ответ модели в Gekko, но с большим количеством шагов.

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

Вот способ построения модели без функции 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()
  

Вы можете построить модель ARX с IMODE=1 на Python, если вы управляете значениями временных рядов с уникальными именами переменных в каждый момент времени. Обратите внимание, что ваши MV и CV метки поменялись местами в примере, который вы опубликовали. CV является управляемой переменной и является выходным прогнозируемым значением. MV — это то значение, которое может быть скорректировано вручную оператором или другим решателем.

Если вы заглянете внутрь функции sysid, вы также увидите пример того, как построить модель ARX без помощи функции ARX, но для многомерного случая. Это сложнее, поэтому я не рекомендую использовать этот подход.

 syid.Raw('Objects')
syid.Raw('  sum_a[1:ny] = sum(%i)'%na)
syid.Raw('  sum_b[1:ny][1::nu] = sum(%i)'%nbk)
syid.Raw('End Objects')
syid.Raw('  ')
syid.Raw('Connections')
syid.Raw('  a[1:na][1::ny] = sum_a[1::ny].x[1:na]')
syid.Raw('  b[1:nb][1::nu][1:::ny] = sum_b[1:::ny][1::nu].x[1:nb]')
syid.Raw('  sum_a[1:ny] = sum_a[1:ny].y')
syid.Raw('  sum_b[1:ny][1::nu] = sum_b[1:ny][1::nu].y')
syid.Raw('End Connections')
syid.Raw('  ')
syid.Raw('Constants')
syid.Raw('  n = %i' %n)
syid.Raw('  nu = %i'%nu)
syid.Raw('  ny = %i'%ny)
syid.Raw('  na = %i'%na)
syid.Raw('  nb = %i'%nbk)
syid.Raw('  m = %i'%m)
syid.Raw('  ')
syid.Raw('Parameters')
syid.Raw('  a[1:na][1::ny] = 0.9 !>= 0.00001 <= 0.9999999')
syid.Raw('  b[1:nb][1::nu][1:::ny] = 0')
syid.Raw('  c[1:ny] = 0')
syid.Raw('  u[1:n][1::nu]')
syid.Raw('  y[1:m][1::ny]')
syid.Raw('  z[1:n][1::ny]')
syid.Raw('  Ks[1:ny][1::nu] = 1')
syid.Raw('  ')
syid.Raw('Variables')
syid.Raw('  y[m 1:n][1::ny] = 0')
syid.Raw('  sum_a[1:ny] = 0 !<= 1')
syid.Raw('  sum_b[1:ny][1::nu] = 0')
syid.Raw('  K[1:ny][1::nu] = 0 >=-1e8 <=1e8')
syid.Raw('  ')
syid.Raw('Equations')
if pred=='model':
    # use model to predict next y (Output error)
    eqn = '  y[m 1:n][1::ny] = a[1][1::ny]*y[m:n-1][1::ny]'
else:
    # use measurement to predict next y (ARX)
    eqn = '  y[m 1:n][1::ny] = a[1][1::ny]*z[m:n-1][1::ny]'
for j in range(1,nu 1):
    eqn  = ' b[1][%i][1::ny]*u[m:n-1][%i]'%(j,j,)
    for i in range(2,nbk 1): 
        eqn  = ' b[%i][%i][1::ny]*u[m-%i:n-%i][%i]'%(i,j,i-1,i,j,)
if pred=='model':
    # use model to predict next y (Output error)
    seqn = ' a[%i][1::ny]*y[m-%i:n-%i][1::ny]'
else:
    # use measurement to predict next y (ARX)
    seqn = ' a[%i][1::ny]*z[m-%i:n-%i][1::ny]'        
for i in range(2,na 1): 
    eqn  = seqn%(i,i-1,i,)
eqn  = ' c[1::ny]'
syid.Raw(eqn)
syid.Raw('')
syid.Raw('  K[1:ny][1::nu] * (1 - sum_a[1:ny]) = Ks[1:ny][1::nu] * sum_b[1:ny][1::nu]')        
syid.Raw('  minimize %e * (y[m 1:n][1::ny] - z[m 1:n][1::ny])^2'%objf)
syid.Raw('  minimize 1e-3 * a[1:na][1::ny]^2')
syid.Raw('  minimize 1e-3 * b[1:nb][1::nu][1:::ny]^2')
syid.Raw('  minimize 1e-3 * c[1:ny]^2')