Pyomo调优SARIMAX超参数时约束不满足的问题求助

我正在开发一款SARIMAX（带外生回归项的季节性自回归积分滑动平均）模型超参数调优求解器，核心目标是确保模型训练集与测试集R²均≥0.8，同时最大化两者的R²之和。目前使用Pyomo构建优化模型，尝试了glpk、ipopt、cbc等求解器，求解器均返回成功，但最终得到的参数始终违反R²≥0.8的约束条件。

已尝试的排查方向：

• 重新评估超参数的取值边界
• 检查目标函数的计算公式
• 调整约束条件的实现方式

期望得到满足训练集与测试集精度阈值的稳健SARIMAX模型。

复现代码 #

import numpy as np
import pandas as pd
from pyomo.environ import *
from statsmodels.tsa.statespace.sarimax import SARIMAX
from sklearn.metrics import r2_score

# Step 1: Create a synthetic time series dataset
np.random.seed(42)
n = 120  # Length of time series data

# Generate a seasonal pattern with trend and noise
time = np.arange(n)
seasonal_pattern = 10 * np.sin(2 * np.pi * time / 12)
trend = 0.1 * time
noise = np.random.normal(0, 1, n)
data = 20 + trend + seasonal_pattern + noise

# Create a DataFrame
df = pd.DataFrame(data, columns=["value"])
df['date'] = pd.date_range(start="2010-01-01", periods=n, freq="D")
df.set_index("date", inplace=True)

# Step 2: Split into train and test sets
train_size = int(len(df) * 0.8)
train, test = df.iloc[:train_size], df.iloc[train_size:]

# Step 3: Define the Pyomo model
model = ConcreteModel()

# Define integer variables for SARIMA parameters
model.p = Var(within=NonNegativeIntegers, bounds=(0, 10), initialize=1)
model.d = Var(within=NonNegativeIntegers, bounds=(0, 10), initialize=1)
model.q = Var(within=NonNegativeIntegers, bounds=(0, 10), initialize=1)
model.P = Var(within=NonNegativeIntegers, bounds=(0, 10), initialize=1)
model.D = Var(within=NonNegativeIntegers, bounds=(0, 10), initialize=1)
model.Q = Var(within=NonNegativeIntegers, bounds=(0, 10), initialize=1)
model.seasonality = Var(within=NonNegativeIntegers, bounds=(2, 12), initialize=2)
print("Displaying P value ",model.p.display())

# Define parameters to store R² scores
model.train_r2 = Var(within=NonNegativeReals)
model.test_r2 = Var(within=NonNegativeReals)
model.test_train_sum = Var(within=NonNegativeReals)


# Step 4: Define a function to fit SARIMAX and calculate R² values
def sarima_r2_model(p, d, q, P, D, Q, seasonality):
    try:
        # Fit SARIMAX model with the given parameters
        model = SARIMAX(
            train['value'],
            order=(p, d, q),
            seasonal_order=(P, D, Q, seasonality),
            enforce_stationarity=False,
            enforce_invertibility=False
        )
        results = model.fit(disp=False)

        # Calculate R² on both train and test sets
        train_pred = results.predict(start=train.index[0], end=train.index[-1])
        test_pred = results.predict(start=test.index[0], end=test.index[-1])
        train_r2_value = r2_score(train['value'], train_pred)
        test_r2_value = r2_score(test['value'], test_pred)
        return train_r2_value, test_r2_value
    except:
        print("========== If model fitting fails, return low R² values")
        return -10, -10

# Define objective to maximize sum of R² scores
def objective_rule(m):
    p = int(m.p.value)
    d = int(m.d.value)
    q = int(m.q.value)
    P = int(m.P.value)
    D = int(m.D.value)
    Q = int(m.Q.value)
    seasonality = int(m.seasonality.value)
    m.train_r2.value, m.test_r2.value = sarima_r2_model(p, d, q, P, D, Q, seasonality)
    
    # Store R² values in model variables for constraints

    print("train_r2_value ==========",m.train_r2.value)
    print("test_r2_value ==========",m.test_r2.value)
    model.test_train_sum = m.train_r2.value + m.test_r2.value 
    return model.test_train_sum


# Set objective function
model.objective = Objective(rule=objective_rule, sense=maximize)



# Step 5: Add constraints for minimum R² threshold
r2_threshold = 0.8
model.train_r2_constraint = Constraint(expr=model.train_r2 >= r2_threshold)
model.test_r2_constraint = Constraint(expr=model.test_r2>=r2_threshold)


# Step 6: Solve the optimization problem with a suitable solver
solver = SolverFactory('cbc',executable="/usr/bin/cbc")  # Choose solver like 'cbc' or 'glpk' if available


results = solver.solve(model, tee=True)
print(results.write())
print(results.solver.status)
print("optimal",results.solver.termination_condition)

print("# will show you the evaluation of the OBJ function")
model.objective.display()  
print(model.objective.expr()) 

# Step 7: Retrieve the best parameters
best_params = {
    'p': int(model.p.value),
    'd': int(model.d.value),
    'q': int(model.q.value),
    'P': int(model.P.value),
    'D': int(model.D.value),
    'Q': int(model.Q.value),
    'seasonality': int(model.seasonality.value)
}
print("Optimized parameters:", best_params)

# Step 8: Train the final model with optimized parameters
best_model = SARIMAX(
    train['value'],
    order=(best_params['p'], best_params['d'], best_params['q']),
    seasonal_order=(best_params['P'], best_params['D'], best_params['Q'], best_params['seasonality']),
    enforce_stationarity=False,
    enforce_invertibility=False
)
best_results = best_model.fit(disp=False)

# Step 9: Final evaluation on train and test sets
train_pred = best_results.predict(start=train.index[0], end=train.index[-1])
test_pred = best_results.predict(start=test.index[0], end=test.index[-1])

print("Final Train R²:", r2_score(train['value'], train_pred))
print("Final Test R²:", r2_score(test['value'], test_pred))

核心问题 #

求解器返回成功，但最终参数对应的训练集/测试集R²无法达到预设的0.8阈值，约束条件未被实际生效。

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