如何用Gekko构建并求解最优功率损耗优化问题？

问题描述 #

需要最小化以下目标函数：

sum[sum[Ri*(Pi² + (Qi - Qcj*Xij)²) for j in range(Nc)] for i in range(N)]

其中：

• P、Q 为常量
• Qc 是候选解列表
• X 为二元决策变量（取值0或1）

当前使用Gekko求解时，无论是否设置X的初始值，总功率损耗始终为350684.76，目标函数值为9554569.42，无法得到更优解，需排查问题并优化。

现有Gekko求解代码 #

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)

# 初始值设置（可选）
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,1],
               [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, 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, 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]])
for i in range(N):
    for j in range(Nc):
        X[i,j].value = bv[i,j]
        #m.Equation(X[i,j]==bv[i,j])

# 单位转换：P、Q转为KW，R转为千欧
for i in range(N):
    Q[i]=Q[i]*1000
    P[i]=P[i]*1000
    R[i]=R[i]*0.001

# 约束：每个i最多选一个j（X[i][j]之和≤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)

# 求解器参数设置
m.solver_options = ['minlp_gap_tol 1e-5',\
                    'minlp_maximum_iterations 10000',\
                    'minlp_max_iter_with_int_sol 2000']
m.options.SOLVER = 1  # 使用APOPT求解器
m.solve(debug=0, disp=True)

总损耗计算代码 #

Ploss=np.zeros(N)
for i in range(N):
    for j in range(Nc):
        #print(i,j,X[i][j].value[0])
        Ploss[i]= R[i]*((P[i]**2)+(Q[i]-Qc[j]*X[i][j].value[0])**2)
print('the total Ploss en KW',sum(Ploss))

问题排查与修正方案 #

1. 总损耗计算逻辑错误 #

**这是严重错误！**原计算代码中，Ploss[i]会被每个j的循环覆盖，只有最后一个j的计算结果生效，即使X[i][j]为0，这导致计算结果完全无法反映真实的损耗值。正确的计算应该只累加X[i][j]==1的对应项：

Ploss=np.zeros(N)
for i in range(N):
    for j in range(Nc):
        if X[i][j].value[0] > 0.5:  # 整数变量取值为0或1，用阈值判断
            Ploss[i] = R[i] * ((P[i]**2) + (Q[i] - Qc[j])**2)
            break  # 每个i最多对应一个j，找到后跳出循环
print('总损耗（KW）：', sum(Ploss))

2. 求解器参数优化 #

当前使用APOPT求解器（SOLVER=1），可以调整参数帮助探索更多解空间，避免陷入局部最优：

# 调整求解器参数
m.solver_options = ['minlp_gap_tol 1e-7',
                    'minlp_maximum_iterations 20000',
                    'minlp_max_iter_with_int_sol 5000',
                    'minlp_branch_method 2']  # 广度优先搜索，更易找到全局最优

3. 约束条件明确化 #

如果每个i必须选择一个Qc[j]（而非可选），可以将约束从sum(X[i][j])<=1改为sum(X[i][j])==1，缩小解空间，帮助求解器更快收敛：

for i in range(N):
    m.Equation(m.sum(X[i,:]) == 1)

4. 数值量级优化 #

单位转换后P、Q的量级达到1e3，平方后为1e6，与Qc的差值平方量级接近，但R被缩小到1e-3，导致目标函数数值跨度较大。可以考虑统一量级（比如不转换单位，直接用原始值计算），减少求解器的数值计算压力。

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