如何用Gekko构建含二进制变量的优化问题？求解遇阻

问题描述 #

需要最小化以下目标函数：
$$\sum_{i=0}^{N-1} \sum_{j=0}^{N_c-1} R_i \left{ P_i^2 + (Q_i - Q_{c,j} \cdot X_{ij})^2 \right}$$
其中：

• $P$、$Q$ 为已知常数
• $Q_c$ 是候选解列表
• $X$ 是二进制决策变量（取值为0或1）
• 约束条件：每个索引$i$对应的所有$X_{ij}$之和 ≤1（即每个$i$最多选择一个$j$）
• 核心目标：找到使目标函数取最小值的$X$向量

现有实现代码 #

from gekko import GEKKO
import numpy as np
P=[13.10511598922975,11.2611396806742,10.103920431906348,8.199519500182628,6.411296067052755,4.753519719147589,3.8977762462825973,2.6593092284662734,1.6399999999854893]
Q=[5.06643685386732,4.4344047044589585,3.8082608015186405,3.2626022579039584,1.2568869621197523,0.6152693459109657,0.46237064874523776,0.35226399840832523,0.20000000001140983]
R=[0.1233, 0.014, 0.7463, 0.6984, 1.9831, 0.9053, 2.0552, 4.7953, 5.3434]
Qc=[150, 300, 450, 600,750, 900,1050, 1200,1350,1500,1650,1800,1950,2100,2250,2400,2550,2700,2850,3000,3150,3300,3450,3600,3750,3900,4050]
N=len(Q)
Nc=len(Qc)
m = GEKKO(remote=False)
X = m.Array(m.Var,(N,Nc),integer=True,lb=0,ub=1,value=0)
# 转换P、Q为KW单位
for i in range(N):
    Q[i]=Q[i]*1000
    P[i]=P[i]*1000
# 约束条件：每个i对应的X行和≤1
for i in range(N):
    m.Equation(m.sum([X[i][j]for j in range(Nc)])<=1)    
# 构建目标函数
b=m.sum([m.sum([R[i]*((P[i]**2)+((Q[i])-Qc[j]*X[i][j])**2) for j in range(Nc)]) for i in range(N)])
m.Minimize(b)

尝试的三种求解方法 #

方法1：直接使用APOPT求解器 #

m.options.SOLVER = 1
m.solve()

方法2：设置初始值后用APOPT求解 #

bv = np.array([[0, 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
 [0, 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1],
 [0, 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0],
 [0, 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
 [0, 0, 0, 0, 0, 0, 1, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
 [0, 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
 [0, 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
 [0, 0, 0, 0, 1, 0, 0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
 [0, 0, 0, 0, 1, 0, 0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]])
for i in range(N):
    for j in range(Nc):
        X[i,j].value = bv[i,j]
m.options.SOLVER = 1
m.solve()

方法3：先IPOPT后APOPT求解 #

m.options.SOLVER = 3
m.solve(debug=0, disp=True)
m.options.SOLVER = 1
m.solve(debug=0, disp=True)

问题 #

上述三种方法均未得到最优解，寻求可行的解决办法。

解决思路与调整方案 #

1. 优化数值尺度

原代码中$P$、$Q$乘以1000后与$Q_c$的数值差距过大（如$Q_i$可达5000+，$Q_c$最大为4050），导致求解器难以平衡目标函数中的各项权重。建议对所有变量进行归一化处理，缩小数值尺度：

scale = 1000
P = [p / scale for p in P]
Q = [q / scale for q in Q]
Qc = [qc / scale for qc in Qc]

2. 调整求解器参数

针对APOPT求解器（SOLVER=1），增大迭代次数并减小收敛容差，让求解器有足够空间寻找最优解：

m.options.MAX_ITER = 10000  # 增大迭代次数
m.options.TOL = 1e-6        # 减小收敛容差

3. 生成高质量初始解

手动设置的初始解可能不是最优方向，建议用贪心策略生成初始解：对每个$i$，计算选择每个$j$时的项值，选择使该项最小的$j$设为1，其余为0。这能帮助求解器快速收敛到更优区域。

4. 明确约束条件

原约束为$\sum_j X_{ij} \leq 1$，若需求是每个$i$必须选择一个$j$，需将约束改为$\sum_j X_{ij} = 1$，避免求解器输出全0的平凡解。

调整后的完整代码 #

from gekko import GEKKO
import numpy as np

# 原始数据
P = [13.10511598922975, 11.2611396806742, 10.103920431906348, 8.199519500182628, 6.411296067052755, 4.753519719147589, 3.8977762462825973, 2.6593092284662734, 1.6399999999854893]
Q = [5.06643685386732, 4.4344047044589585, 3.8082608015186405, 3.2626022579039584, 1.2568869621197523, 0.6152693459109657, 0.46237064874523776, 0.35226399840832523, 0.20000000001140983]
R = [0.1233, 0.014, 0.7463, 0.6984, 1.9831, 0.9053, 2.0552, 4.7953, 5.3434]
Qc = [150, 300, 450, 600, 750, 900, 1050, 1200, 1350, 1500, 1650, 1800, 1950, 2100, 2250, 2400, 2550, 2700, 2850, 3000, 3150, 3300, 3450, 3600, 3750, 3900, 4050]

# 数值归一化
scale = 1000
P = [p / scale for p in P]
Q = [q / scale for q in Q]
Qc = [qc / scale for qc in Qc]

N = len(Q)
Nc = len(Qc)

m = GEKKO(remote=False)
X = m.Array(m.Var, (N, Nc), integer=True, lb=0, ub=1, value=0)

# 约束条件：每个i最多选一个j（若必须选则改为==1）
for i in range(N):
    m.Equation(m.sum(X[i, :]) <= 1)

# 构建目标函数
obj = 0
for i in range(N):
    for j in range(Nc):
        obj += R[i] * (P[i]**2 + (Q[i] - Qc[j] * X[i,j])**2)
m.Minimize(obj)

# 设置求解器参数
m.options.SOLVER = 1
m.options.MAX_ITER = 10000
m.options.TOL = 1e-6

# 生成贪心初始解
for i in range(N):
    min_val = float('inf')
    best_j = -1
    for j in range(Nc):
        current_val = R[i] * (P[i]**2 + (Q[i] - Qc[j])**2)
        if current_val < min_val:
            min_val = current_val
            best_j = j
    if best_j != -1:
        X[i, best_j].value = 1

# 求解并输出结果
m.solve(disp=True)

print("最优X向量：")
for i in range(N):
    row = [int(round(X[i,j].value[0])) for j in range(Nc)]
    print(f"第{i+1}行：{row}")

print(f"最小目标函数值：{m.options.objfcnval}")

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