gekko - 为什么这段代码使用不同的程序会返回不同的结果？

Question

此代码在 python gekko 中的目标值为 568，但此代码的 MATLAB 版本为 70。我不明白为什么。也许关于求解器选项？

from gekko import GEKKO    
import numpy as np

m = GEKKO(remote=False)
capacity=np.array([70, 55, 51, 43, 41, 80])
demand=np.array([60, 57, 62, 38, 70])
cost=np.array([[5, 4, 5, 7, 2],
   [2, 9, 2, 6, 3],
   [6, 5, 1, 7, 9],
    [7, 3, 9, 3, 5],
    [4, 8, 7, 9, 7],
    [2, 5, 4, 2, 1]])

x = m.Array(m.Var,(6,5),lb=0,integer=True)

for i in range(6):
     for j in range(5):
           m.Minimize(cost[i,j]*x[i,j])

for i in range(6):
    m.Equation(m.sum(x[i,:])<=capacity[i])    


for j in range(5):
    m.Equation(m.sum(x[:,j])==demand[j])     


m.options.solver = 1
m.solve()
print('Objective Function: ' + str(m.options.objfcnval))
print(x)

这是给出 fval=70 的 MATLAB 代码版本。也许我写错了代码，但看起来大部分是一样的，我不明白为什么。感谢帮助。

clear
clc
capacity_points=6;
demand_points=5;
capacity=[70 55 51 43 41 80];
demand=[60 57 62 38 70];
cost=[5 4 5 7 2;
    2 9 2 6 3;
    6 5 1 7 9;
    7 3 9 3 5;
    4 8 7 9 7;
    2 5 4 2 1];
x=optimvar('x',capacity_points,demand_points,'Type','integer','LowerBound',0);
%expr=optimexpr;
expr=0;
for i=1:capacity_points
    for j= 1:demand_points
        expr=expr+cost(i,j)*x(i,j);
    end
end

%const1=optimconstr(capacity_points);
for i=1:capacity_points
    const1=sum(x(i,:)) <=capacity(i);
end
%const2=optimconstr(demand_points);
for j=1:demand_points
    const2=sum(x(:,j)) ==demand(j);
end

prob=optimproblem;
prob.Objective=expr;
prob.Constraints.const1=const1;
prob.Constraints.const2=const2;
sol=solve(prob)
[sol,fval] = solve(prob)

Answer 1

看来 Gekko 中的解决方案是正确的。您可能还想发布 MATLAB 代码，以便我们可以查看问题陈述之间是否存在任何差异。

当您想要确定解决方案是否最优时，这里有一些一般的故障排除策略:

• 确定问题是否是非凸的。此问题是一个符合 minimize c x 形式的线性规划 (LP) 问题受制于 A x = b和 A x < b . LP 是凸的，因此局部解也是全局解。如果问题是非凸的，那么您可以尝试不同的初始猜测以查看解决方案是否发生变化，或者使用全局优化器为您自动完成。

• 尝试不同的求解器。您可以设置一个循环并使用 m.options.SOLVER 更改求解器.如果您不想更改代码，请使用 m.options.SOLVER = 0尝试所有求解器。这是关于 SOLVER=0 问题的输出.

 Solver           Objective  Solution Time  Status
 -------------- ------------ ------------- ---------
 APOPT (v1.0)    5.68000E+02         0.036 Success  
 BPOPT (v1.0)    5.68001E+02         0.016 Success  
 IPOPT (v3.12)   5.68000E+02         0.017 Success  
 IPOPT (v2.3)    0.00000E+00         0.000 Skip     
 SNOPT (v6.1)    0.00000E+00         0.000 Skip     
 MINOS (v5.5)    0.00000E+00         0.000 Skip     
 -------------- ------------------------------------

免费提供的求解器都得出了相同的解决方案，因此这表明多种求解器方法都达成了相同的共识。

• 对于除小规模问题(<100 个变量)以外的任何问题，都可能难以理解优化解决方案。另一种策略是降低问题的复杂性，直到您可以得出解析解或直到解易于验证。对于您的情况，列必须等于 demand而该行必须位于 capacity 下方约束。

capacity=np.array([70, 55, 51, 43, 41, 80])
demand  =np.array([60, 57, 62, 38, 70])
cost    =np.array([[5, 4, 5, 7, 2],
                   [2, 9, 2, 6, 3],
                   [6, 5, 1, 7, 9],
                   [7, 3, 9, 3, 5],
                   [4, 8, 7, 9, 7],
                   [2, 5, 4, 2, 1]])

如果我用np.min(cost,axis=0)为每一列挑出最低成本的项目来满足需求, 它是 [2 3 1 2 1] .如果我将其乘以 demand与 np.dot(np.min(cost,axis=0),demand)它给出了 499 的目标函数值.这是没有 capacity 的最小目标函数约束。它表明 70 的目标除非 demand 否则不可能解决此问题或 cost减少了。

虽然这些策略特定于您的问题，但它们可以应用于其他优化问题，以诊断和解决非直观的解决方案。如果求解器说它找到了解决方案，那么 Karush-Kuhn-Tucker conditions满足最优性，它至少是一个局部解决方案。

对编辑的回应

MATLAB 代码的问题是只有最后一个约束被强制执行，因为 const1 和 const2 在每个循环中都被重新定义。这给出了解决方案 sol.x .

     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    70

您需要在 MATLAB 中包含所有等式和不等式约束。

clear
clc
capacity_points=6;
demand_points=5;
capacity=[70 55 51 43 41 80];
demand=[60 57 62 38 70];
cost=[5 4 5 7 2;
    2 9 2 6 3;
    6 5 1 7 9;
    7 3 9 3 5;
    4 8 7 9 7;
    2 5 4 2 1];
x=optimvar('x',capacity_points,demand_points,'Type','integer','LowerBound',0);
%expr=optimexpr;
expr=0;
for i=1:capacity_points
    for j= 1:demand_points
        expr=expr+cost(i,j)*x(i,j);
    end
end

%const1=optimconstr(capacity_points);
for i=1:capacity_points
    const1(i)=sum(x(i,:)) <=capacity(i);
end
%const2=optimconstr(demand_points);
for j=1:demand_points
    const2(j)=sum(x(:,j)) ==demand(j);
end

prob=optimproblem;
prob.Objective=expr;
prob.Constraints.const1 = const1;
prob.Constraints.const2 = const2;
sol=solve(prob)
[sol,fval] = solve(prob)

这给出了相同的目标函数值 568作为壁虎。