基于ARX模型与Gekko的系统辨识：室内温度异常问题排查

问题现象 #

当外界温度上升时，出现不符合预期的系统响应：

• β不变时，室内温度未同步上升，反而保持恒定
• 极端场景下，外界温度上升、β维持不变，室内温度出现下降（示例图1、图2）

已执行的排查与修正 #

• 原始数据采样间隔为2秒（每日43200条），最初采用ARX模型进行系统辨识时错误设置了300秒采样时间，现已修正为2秒采样间隔的MPC代码
• 将单日外部温度数据（data_2）拼接为连续两天的温度曲线，重新运行代码后问题仍未解决

修正后的MPC代码 #

# Import library
import numpy as np
import pandas as pd
import time
from gekko import GEKKO
from numpy import array
import matplotlib.pyplot as plt

# import data
tx = pd.read_csv(r'C:\Users\mpc\tx_b.csv')
tz = pd.read_csv(r'C:\Users\mpc\tz.csv')

data = pd.concat([tx,tz],axis=1)
data_1 = data[0:8596800]
data_2 = data[8596800:]

#%%  Initialize Model
m = GEKKO(remote=False)

# system identification
ts = 2
t = np.arange(0,len(data_1)*ts, ts)
u_id = data_1[['Tx_i','beta_i']]
y_id = data_1[['Tz_i']]

#meas : the time-series next step is predicted from prior measurements as in ARX

na=10; nb=10 # ARX coefficients
print('Identify model')
start = time.time()
yp,p,K = m.sysid(t,u_id,y_id,na,nb,objf=100,scale=False,diaglevel=0,pred='meas')
print('temps de prediction :'+str(time.time()-start)+'s')

#%% parametres 
K = array([[  0.93688819, -12.22410568]])
p = {'a': array([[ 1.08945931],
        [-0.00243571],
        [-0.00247112],
        [-0.00273341],
        [-0.00296342],
        [-0.00319516],
        [-0.00343794],
        [-0.00366398],
        [-0.00394255],
        [-0.06661506]]),
 'b': array([[[-0.05134201, -0.01035174],
         [ 0.00170311, -0.01551259],
         [ 0.00172715, -0.01178932],
         [ 0.00178147, -0.01051817],
         [ 0.00184694, -0.00821511],
         [ 0.00192371, -0.00570574],
         [ 0.00201409, -0.00344425],
         [ 0.00210016, -0.0014708 ],
         [ 0.00222189,  0.00021622],
         [ 0.03789636,  0.04235503]]]),
 'c': array([0.0266222])}


#%% I used the last day's external temperature data as a disturbance.
T_externel = data_2[["Tx_i"]].values

m = GEKKO(remote=False)
m.y = m.Array(m.CV,1)
m.u = m.Array(m.MV,2)
m.arx(p,m.y,m.u)

# rename CVs
m.T = m.y[0]

# rename MVs
m.beta = m.u[1]

# distrubance
m.d = m.u[0]

# distrubance and parametres
m.d = m.Param(T_externel[0])

m.bias = m.Param(0)
m.Tb = m.CV()
m.Equation(m.Tb==m.T+m.bias)

# steady state initialization
m.options.IMODE = 1
m.solve(disp=False)

# set up MPC
m.d.value = T_externel

m.options.IMODE   = 6 # MPC
m.options.CV_TYPE = 1 # the objective is an l1-norm (region)
m.options.NODES   = 2 # Collocation nodes
m.options.SOLVER  = 3 # IPOPT
m.time = np.arange(0,len(T_externel)*2,2) # step time = 300s

# Manipulated variables
m.beta.STATUS = 1  # calculated by the optimizer
m.beta.FSTATUS = 0 # use measured value
m.beta.DCOST = 0.0 # Delta cost penalty for MV movement
m.beta.UPPER = 1.0 # Upper bound
m.beta.LOWER = 0.0 # Lower bound
m.beta.MV_STEP_HOR = 1
m.beta.value = 0

# Controlled variables
m.Tb.STATUS = 1        # drive to set point
m.Tb.FSTATUS = 0       # receive measurement
m.Tb.SPHI = 17.5         # set point high level
m.Tb.SPLO = 16.5         # set point low level
m.Tb.WSPHI = 100       # set point high priority
m.Tb.WSPLO = 100       # set point low priority

T_MEAS = 20
m.Tb.value = T_MEAS
m.bias.value = T_MEAS - m.T.value[0]

m.options.SOLVER = 3
m.solve(disp=False)

if m.options.APPSTATUS == 1:
   # Retrieve new values
   beta  = m.beta.NEWVAL
   for i in range(43200):
       print(i)
else:
   # Solution failed
   beta  = 0.0

# Plot the results
plt.figure(figsize=(8,3.5))
plt.subplot(3,1,1)
plt.plot(m.time,m.Tb.value,'r-',label=r'$T_{int}$')
plt.plot([0,m.time[-1]],[m.Tb.SPHI,m.Tb.SPHI],'k--',label='Upper Bound')
plt.plot([0,m.time[-1]],[m.Tb.SPLO,m.Tb.SPLO],'k--',label='Lower Bound')
plt.legend(loc=1); plt.grid()
plt.ylabel('Tin (°C)')
plt.subplot(3,1,2)
plt.plot(m.time,m.d.value,'g:',label=r'$T_{ext}$')
plt.ylabel('Tex (°C)')
plt.subplot(3,1,3)
plt.step(m.time,m.beta.value,'b--',label=r'$\beta$')
plt.ylabel('optimal control')
plt.xlabel('Time (sec)')
plt.legend(loc=1); plt.grid()
plt.savefig('results7.png',dpi=300)
plt.show()

（注：已将单日外部温度数据拼接为连续两天的曲线，相关结果图见示例）

内容的提问来源于stack exchange，提问作者m_nacereddine