Gekko - 优化调度的不可行解决方案，与 gurobi 进行比较 [英] Gekko - infeasible solution to optimal scheduling, comparison w/ gurobi

问题描述

我对 Gurobi 有点熟悉，但过渡到 Gekko，因为后者似乎有一些优势.不过，我遇到了一个问题，我将使用我想象中的苹果园来说明.为期 5 周的收获期 (#horizo​​n: T=5) 即将来临，我的——非常微薄的——产品将是:[3.0, 7.0, 9.0, 5.0, 4.0]我自己保留的一些苹果[2.0, 4.0, 2.0, 4.0, 2.0]，剩下的产品我将在农贸市场以下列价格出售:[0.8, 0.9, 0.5,1.2, 1.5].我有可容纳 6 个苹果的存储空间，因此我可以提前计划并在最佳时机销售苹果，从而最大限度地提高我的收入.我尝试使用以下模型确定最佳时间表:

I am somewhat familiar with Gurobi, but transitioning to Gekko since the latter appears to have some advantages. I am running into one issue though, which I will illustrate using my imaginary apple orchard. The 5-weeks harvest period (#horizon: T=5) is upon us, and my - very meagre - produce will be: [3.0, 7.0, 9.0, 5.0, 4.0] Some apples I keep for myself [2.0, 4.0, 2.0, 4.0, 2.0], the remaining produce I will sell in the farmer's market at the following prices: [0.8, 0.9, 0.5, 1.2, 1.5]. I have storage space with room for 6 apples, so I can plan ahead and sell apples at the most optimal moments, hence maximizing my revenue. I try to determine the optimal schedule with the following model:

m       = GEKKO()
m.time  = np.linspace(0,4,5)
orchard   = m.Param([3.0, 7.0, 9.0, 5.0, 4.0])
demand    = m.Param([2.0, 4.0, 2.0, 4.0, 2.0]) 
price     = m.Param([0.8, 0.9, 0.5, 1.2, 1.5])

### manipulated variables
# selling on the market
sell                = m.MV(lb=0)
sell.DCOST          = 0
sell.STATUS         = 1
# saving apples
storage_out         = m.MV(value=0, lb=0)
storage_out.DCOST   = 0      
storage_out.STATUS  = 1 
storage_in          = m.MV(lb=0)
storage_in.DCOST    = 0
storage_in.STATUS   = 1

### storage space 
storage         = m.Var(lb=0, ub=6)
### constraints
# storage change
m.Equation(storage.dt() == storage_in - storage_out) 

# balance equation
m.Equation(sell + storage_in + demand == storage_out + orchard)

# Objective: argmax sum(sell[t]*price[t]) for t in [0,4]
m.Maximize(sell*price)
m.options.IMODE=6
m.options.NODES=3
m.options.SOLVER=3
m.options.MAX_ITER=1000
m.solve()

出于某种原因，这是不可行的(错误代码 = 2).有趣的是，如果将 demand[0] 设置为 3.0，而不是 2.0(即等于 orchard[0]，该模型确实会产生一个成功的解决方案.

For some reason this is unfeasible (error code = 2). Interestingly, if set demand[0] to 3.0, instead of 2.0 (i.e. equal to orchard[0], the model does produce a succesful solution.

1. 为什么会这样?
2. 即使是成功的"输出值有点奇怪:存储空间没有使用一次，并且 storage_out 在最后一个时间步中没有得到适当的约束.显然，我没有正确地制定约束条件.我应该怎么做才能获得与 gurobi 输出相当的真实结果(请参阅下面的代码)?

1. Why is this the case?
2. Even the "succesful" output values are bit weird: the storage space is not used a single time, and storage_out is not properly constrained in the last timestep. Clearly, I am not formulating the constraints correctly. What should I do to get realistic results, which are comparable to the gurobi output (see code below)?

output = {'sell'    : list(sell.VALUE),
        's_out'     : list(storage_out.VALUE),
        's_in'      : list(storage_in.VALUE), 
        'storage'   : list(storage.VALUE)}
df_gekko = pd.DataFrame(output)
df_gekko.head()

>   sell  s_out     s_in        storage
0   0.0   0.000000  0.000000    0.0
1   3.0   0.719311  0.719311    0.0
2   7.0   0.859239  0.859239    0.0
3   1.0   1.095572  1.095572    0.0
4   26.0  24.124924 0.124923    0.0 

Gurobi 模型使用 demand = [3.0, 4.0, 2.0, 4.0, 2.0] 求解.请注意，gurobi 还生成了一个具有 demand = [2.0, 4.0, 2.0, 4.0, 2.0] 的解决方案.这只会对结果产生微不足道的影响:在 t=0 出售的 n 个苹果变成 1.

Gurobi model solved for with demand = [3.0, 4.0, 2.0, 4.0, 2.0]. Note that gurobi also produces a solution with demand = [2.0, 4.0, 2.0, 4.0, 2.0]. This only has a trivial impact on the outcome: n apples sold at t=0 becomes 1.

T = 5
m = gp.Model()
### horizon (five weeks)

### supply, demand and price data  
orchard   = [3.0, 7.0, 9.0, 5.0, 4.0]
demand    = [3.0, 4.0, 2.0, 4.0, 2.0] 
price     = [0.8, 0.9, 0.5, 1.2, 1.5]

### manipulated variables
# selling on the market
sell = m.addVars(T)

# saving apples
storage_out = m.addVars(T)
m.addConstr(storage_out[0] == 0)
storage_in  = m.addVars(T)

# storage space
storage = m.addVars(T)
m.addConstrs((storage[t]<=6) for t in range(T))
m.addConstrs((storage[t]>=0) for t in range(T))
m.addConstr(storage[0] == 0)

# storage change
#m.addConstr(storage[0] == (0 - storage_out[0]*delta_t + storage_in[0]*delta_t))
m.addConstrs(storage[t] == (storage[t-1] - storage_out[t] + storage_in[t]) for t in range(1, T))

# balance equation
m.addConstrs(sell[t] + demand[t] + storage_in[t] == (storage_out[t] + orchard[t]) for t in range(T))

# Objective: argmax sum(a_sell[t]*a_price[t] - b_buy[t]*b_price[t])
obj = gp.quicksum((price[t]*sell[t]) for t in range(T))
m.setObjective(obj, gp.GRB.MAXIMIZE)
m.optimize()

输出:

    sell    storage_out storage_in  storage
0   0.0     0.0         0.0         0.0
1   3.0     0.0         0.0         0.0
2   1.0     0.0         6.0         6.0
3   1.0     0.0         0.0         6.0
4   8.0     6.0         0.0         0.0

推荐答案

您可以通过以下方式获得成功的解决方案:

You can get a successful solution with:

m.options.NODES=2

问题在于它正在使用 NODES=3 求解主节点之间的平衡方程.您的微分方程具有线性解，因此 NODES=2 应该足够准确.

The issue is that it is solving the balance equation in between the primary node points with NODES=3. Your differential equation has a linear solution so NODES=2 should be sufficiently accurate.

以下是改进解决方案的其他几种方法:

Here are a couple other ways to improve the solution:

• 对将库存移入或移出仓库设置一个小惩罚.否则，求解器可以使用 storage_in = storage_out 找到较大的任意值.
• 我使用了 m.Minimize(1e-6*storage_in) 和 m.Minimize(1e-6*storage_out).
• 因为初始条件通常是固定的，所以我在开始时使用零值只是为了确保计算第一个点.
• 如果它们以整数单位出售和存储，我也切换到整数变量.如果您需要 SOLVER=1 的整数解，则需要切换到 APOPT 求解器.

• Set a small penalty on moving inventory into or out of storage. Otherwise the solver can find large arbitrary values with storage_in = storage_out.
• I used m.Minimize(1e-6*storage_in) and m.Minimize(1e-6*storage_out).
• Because the initial condition is typically fixed, I used zero values at the beginning just to make sure that the first point is calculated.
• I also switched to integer variables if they are sold and stored in integer units. You need to switch to the APOPT solver if you want an integer solution with SOLVER=1.

 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :  0.058899999999999994 sec
 Objective      :  -17.299986
 Successful solution
 ---------------------------------------------------
 

Sell
[0.0, 0.0, 4.0, 1.0, 1.0, 8.0]
Storage Out
[0.0, 0.0, 1.0, 0.0, 0.0, 6.0]
Storage In
[0.0, 1.0, 0.0, 6.0, 0.0, 0.0]
Storage
[0.0, 1.0, 0.0, 6.0, 6.0, 0.0]

这是修改后的脚本.

from gekko import GEKKO
import numpy as np

m       = GEKKO(remote=False)
m.time  = np.linspace(0,5,6)
orchard   = m.Param([0.0, 3.0, 7.0, 9.0, 5.0, 4.0])
demand    = m.Param([0.0, 2.0, 4.0, 2.0, 4.0, 2.0]) 
price     = m.Param([0.0, 0.8, 0.9, 0.5, 1.2, 1.5])

### manipulated variables
# selling on the market
sell                = m.MV(lb=0, integer=True)
sell.DCOST          = 0
sell.STATUS         = 1
# saving apples
storage_out         = m.MV(value=0, lb=0, integer=True)
storage_out.DCOST   = 0      
storage_out.STATUS  = 1 
storage_in          = m.MV(lb=0, integer=True)
storage_in.DCOST    = 0
storage_in.STATUS   = 1

### storage space 
storage         = m.Var(lb=0, ub=6, integer=True)
### constraints
# storage change
m.Equation(storage.dt() == storage_in - storage_out) 

# balance equation
m.Equation(sell + storage_in + demand == storage_out + orchard)

# Objective: argmax sum(sell[t]*price[t]) for t in [0,4]
m.Maximize(sell*price)
m.Minimize(1e-6 * storage_in)
m.Minimize(1e-6 * storage_out)
m.options.IMODE=6
m.options.NODES=2
m.options.SOLVER=1
m.options.MAX_ITER=1000
m.solve()

print('Sell')
print(sell.value)
print('Storage Out')
print(storage_out.value)
print('Storage In')
print(storage_in.value)
print('Storage')
print(storage.value)

这篇关于Gekko - 优化调度的不可行解决方案，与 gurobi 进行比较的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！