基于Gekko的带自变量贡献约束回归模型构建求助

需求说明 #

开发回归模型 y = a0 + a1*x1 + a2*x2，基于200个数据点，要求所有数据点的a1*x1之和满足 lb1 < sum(a1*x1) < ub1。现有Gekko代码仅能逐行施加约束，无法实现全局求和约束，需调整代码满足需求。

现有代码 #

import gekko as gk
import numpy as np
import pandas as pd

# 假设df是已加载的数据集，X是自变量列名列表，ubound_dict是自变量约束字典
m = gk.GEKKO(remote=False)

m.options.IMODE=2 #Regression mode
y = np.array(df['y']) #因变量
x = np.array(df[X]) #自变量数组

n = x.shape[1] #自变量数量
c  = m.Array(m.FV, n+1) #参数数组（含截距项）
for ci in c:
    ci.STATUS = 1 #启用参数求解
xp = [None]*n

#加载数据
xd = m.Array(m.Param,n)
yd = m.Param(value=y)
for i in range(n):
    xd[i].value = x[:,i]
    xp[i] = m.Var()
    if ubound_dict[i] >= 0:
        xp[i] = m.Var(lb=0, ub=ubound_dict[i])
    elif ubound_dict[i] < 0:
        xp[i] = m.Var(lb=ubound_dict[i], ub=0)
    
    m.Equation(xp[i]==c[i]*xd[i])
    
yp = m.Var()
m.Equation(yp==m.sum([xp[i] for i in range(n)] + [c[n]]))
#最小化预测值与实际值的平方差
m.Minimize((yd-yp)**2)

#使用APOPT求解器
m.options.SOLVER = 1
#求解
m.solve(disp=True)
#提取参数值
a = [i.value[0] for i in c]
print(a)

解决方案 #

核心思路是计算所有数据点的a1*x1全局总和，再对该总和施加范围约束。无需逐行处理，直接利用Gekko的求和函数实现全局约束：

1. 计算a1*x1的全局总和：由于sum(a1*x1) = a1 * sum(x1)，可以直接用参数c[0]（对应a1）乘以所有x1数据点的和，避免循环遍历每个数据点，提升效率。
2. 对全局总和添加上下界约束。

修改后的完整代码：

import gekko as gk
import numpy as np
import pandas as pd

# 假设df是已加载的数据集，X是自变量列名列表，ubound_dict是自变量约束字典
# 定义全局求和的上下界
lb1 = 100  # 根据实际需求设置
ub1 = 500  # 根据实际需求设置

m = gk.GEKKO(remote=False)

m.options.IMODE=2 #Regression mode
y = np.array(df['y']) #因变量
x = np.array(df[X]) #自变量数组

n = x.shape[1] #自变量数量
c  = m.Array(m.FV, n+1) #参数数组（含截距项）
for ci in c:
    ci.STATUS = 1 #启用参数求解
xp = [None]*n

#加载数据
xd = m.Array(m.Param,n)
yd = m.Param(value=y)
for i in range(n):
    xd[i].value = x[:,i]
    xp[i] = m.Var()
    if ubound_dict[i] >= 0:
        xp[i] = m.Var(lb=0, ub=ubound_dict[i])
    elif ubound_dict[i] < 0:
        xp[i] = m.Var(lb=ubound_dict[i], ub=0)
    
    m.Equation(xp[i]==c[i]*xd[i])
    
yp = m.Var()
m.Equation(yp==m.sum([xp[i] for i in range(n)] + [c[n]]))
#最小化预测值与实际值的平方差
m.Minimize((yd-yp)**2)

# --- 新增全局求和约束 ---
# 计算所有数据点a1*x1的总和
sum_a1x1 = m.Var()
# 等价于 sum(a1*x1) = a1 * sum(x1)
m.Equation(sum_a1x1 == c[0] * m.sum(xd[0]))
# 添加范围约束
m.Equation(sum_a1x1 > lb1)
m.Equation(sum_a1x1 < ub1)

#使用APOPT求解器
m.options.SOLVER = 1
#求解
m.solve(disp=True)
#提取参数值
a = [i.value[0] for i in c]
print(a)
# 查看全局求和结果
print(f"sum(a1*x1) = {sum_a1x1.value[0]}")

关键说明 #

• 若需显式遍历每个数据点计算总和（而非利用数学简化），可将求和方程替换为：

  m.Equation(sum_a1x1 == m.sum([c[0]*xd[0][j] for j in range(len(y))]))
  

• 原代码中xp[i]是每个数据点的c[i]*xd[i]值，仅用于逐行计算预测值，无法直接实现全局求和，因此需单独定义全局求和变量并添加约束。

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