基于Python GEKKO的ARX模型参数复用与输入外推问题

我已经用GEKKO的sysid函数创建了ARX模型，完成了训练数据集的参数识别，现在想复用这些参数对其他输入做外推计算。我目前的核心代码片段是：

m = GEKKO(remote=False)
na = n_a
nb = n_b
yp_train, p, K = m.sysid(t_train, u_train, y_train, na, nb)

我尝试的完整代码如下：

def arxouriço2(df, inputs, outputs, n_a, n_b, train_slice=None, test_slice=None):
    datafreime_arx = df.copy()
    if isinstance(inputs, pd.Series):
        inputs = pd.DataFrame(inputs)
    if isinstance(outputs, pd.Series):
        outputs = pd.DataFrame(outputs)
    colunas_arx = list(inputs.columns) + list(outputs.columns)
    mean_col = []
    std_dev_col = []

    for i in colunas_arx:
        mean = np.nanmean(datafreime_arx[i])
        mean_col.append(mean)
        datafreime_arx[i] -= mean

        std_dev = np.nanstd(datafreime_arx[i])
        std_dev_col.append(std_dev)
        datafreime_arx[i] /= std_dev

    if train_slice is not None:
        train_data = datafreime_arx.iloc[train_slice]
        t_train = train_data['Tempo']
        u_train = inputs.iloc[train_slice]
        y_train = outputs.iloc[train_slice]
    else:
        t_train = datafreime_arx['Tempo']
        u_train = inputs
        y_train = outputs

    if test_slice is not None:
        test_data = datafreime_arx.iloc[test_slice]
        t_test = test_data['Tempo']
        u_test = inputs.iloc[test_slice]
        y_test = outputs.iloc[test_slice]
    else:
        t_test = datafreime_arx['Tempo']
        u_test = inputs
        y_test = outputs

    m = GEKKO(remote=False)
    na = n_a
    nb = n_b
    yp_train, p, K = m.sysid(t_train, u_train, y_train, na, nb)

    for i, col in enumerate(colunas_arx):
        datafreime_arx[col] /= std_dev_col[i]
        datafreime_arx[col] -= mean_col[i]

    resis = y_train - yp_train
    resis2 = resis ** 2
    print('Sum of squared residuals: ' + str(np.sum(resis2)))

    if test_slice is not None:
        yp_test = m.sysid(t_test, u_test, y_test, na, nb)[0]
        plt.figure(figsize=(8, 5))
        plt.subplot(2, 1, 1)
        plt.plot(t_train, u_train)
        plt.legend(inputs.columns)
        plt.ylabel('Temperature (ºC)')
        plt.subplot(2, 1, 2)
        plt.plot(t_test, y_test, label=outputs.columns)
        plt.plot(t_test, yp_test, label=outputs.columns + ' predicted')
        plt.legend()
        plt.ylabel('Rotation (º)')
        plt.xlabel('Time')
        plt.tight_layout()
        plt.show()

    return p, K, m

arxouriço2(df_c3_avg5_c2,df_c3_avg5_c2[['STH3_avg','STH-04']],df_c3_avg5_c2['C14-Y'],10,10,range(465, 829),range(829, 903))

请问有没有办法复用已经得到的参数，应用到不同输入上？

问题分析 #

你当前代码在测试阶段重新调用m.sysid(t_test, u_test, y_test, na, nb)，这会重新识别参数而非复用训练好的模型做预测，这是核心问题。要复用ARX参数，需手动基于sysid返回的p和K构建预测逻辑。

解决方案 #

1. 理解ARX模型结构 #

ARX模型的数学形式为：
$$y(t) = -a_1 y(t-1) - ... -a_{n_a} y(t-n_a) + b_1 u(t-1) + ... + b_{n_b} u(t-n_b) + K$$
其中p包含识别得到的$a_1...a_{n_a}$和$b_1...b_{n_b}$系数，K是模型稳态增益偏移。

2. 修改代码实现参数复用 #

替换测试阶段的sysid调用，用训练得到的参数直接计算预测值：

def arxouriço2(df, inputs, outputs, n_a, n_b, train_slice=None, test_slice=None):
    datafreime_arx = df.copy()
    if isinstance(inputs, pd.Series):
        inputs = pd.DataFrame(inputs)
    if isinstance(outputs, pd.Series):
        outputs = pd.DataFrame(outputs)
    
    # 拆分输入输出的标准化参数，避免混淆
    input_cols = inputs.columns
    output_cols = outputs.columns
    input_mean = []
    input_std = []
    output_mean = []
    output_std = []

    # 输入数据标准化
    for col in input_cols:
        mean = np.nanmean(datafreime_arx[col])
        input_mean.append(mean)
        std = np.nanstd(datafreime_arx[col])
        input_std.append(std)
        datafreime_arx[col] = (datafreime_arx[col] - mean) / std
    
    # 输出数据标准化
    for col in output_cols:
        mean = np.nanmean(datafreime_arx[col])
        output_mean.append(mean)
        std = np.nanstd(datafreime_arx[col])
        output_std.append(std)
        datafreime_arx[col] = (datafreime_arx[col] - mean) / std

    # 拆分训练/测试数据
    if train_slice is not None:
        train_data = datafreime_arx.iloc[train_slice]
        t_train = train_data['Tempo'].values
        u_train = train_data[input_cols].values
        y_train = train_data[output_cols].values
    else:
        t_train = datafreime_arx['Tempo'].values
        u_train = datafreime_arx[input_cols].values
        y_train = datafreime_arx[output_cols].values

    if test_slice is not None:
        test_data = datafreime_arx.iloc[test_slice]
        t_test = test_data['Tempo'].values
        u_test = test_data[input_cols].values
        y_test = test_data[output_cols].values
    else:
        t_test = datafreime_arx['Tempo'].values
        u_test = datafreime_arx[input_cols].values
        y_test = datafreime_arx[output_cols].values

    # 训练ARX模型，获取参数
    m = GEKKO(remote=False)
    na = n_a
    nb = n_b
    yp_train, p, K = m.sysid(t_train, u_train, y_train, na, nb)

    # 计算训练残差（反标准化后）
    yp_train_denorm = yp_train * output_std[0] + output_mean[0]
    y_train_denorm = y_train * output_std[0] + output_mean[0]
    resis = y_train_denorm - yp_train_denorm
    print('Sum of squared residuals: ' + str(np.sum(resis**2)))

    if test_slice is not None:
        # 复用参数进行外推预测
        n_test = len(t_test)
        yp_test = np.zeros(n_test)
        # 用测试集前na个实际值作为预热初始化
        yp_test[:na] = y_test[:na]
        
        for t in range(na, n_test):
            # ARX公式计算：y(t) = -Σ(a_i*y(t-i)) + Σ(b_i*u(t-i)) + K
            y_term = -np.dot(p[:na], yp_test[t-na:t])
            # 单输入场景，多输入需循环处理每个输入的b系数
            u_term = np.dot(p[na:na+nb], u_test[t-nb:t, 0])
            yp_test[t] = y_term + u_term + K
        
        # 反标准化恢复原始尺度
        yp_test_denorm = yp_test * output_std[0] + output_mean[0]
        y_test_denorm = y_test * output_std[0] + output_mean[0]
        u_train_denorm = u_train * input_std[0] + input_mean[0]
        u_test_denorm = u_test * input_std[0] + input_mean[0]

        # 绘图展示
        plt.figure(figsize=(8, 5))
        plt.subplot(2, 1, 1)
        plt.plot(t_train, u_train_denorm)
        plt.legend(input_cols)
        plt.ylabel('Temperature (ºC)')
        plt.subplot(2, 1, 2)
        plt.plot(t_test, y_test_denorm, label=output_cols[0])
        plt.plot(t_test, yp_test_denorm, label=output_cols[0] + ' predicted')
        plt.legend()
        plt.ylabel('Rotation (º)')
        plt.xlabel('Time')
        plt.tight_layout()
        plt.show()

    return p, K, m

关键说明 #

• 单独拆分输入输出的标准化参数，避免反标准化时出现错误。
• 预测阶段用ARX数学公式直接计算，完全复用训练得到的p和K参数。
• 前na个时间步用测试集实际值作为预热，保证预测序列的连续性。
• 多输入场景下，需调整u_term的计算逻辑，对应每个输入的b系数组。

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