Gekko - 最优调度的不可行解决方案，与 gurobi 的比较
Gekko - infeasible solution to optimal scheduling, comparison w/ gurobi
我对 Gurobi 比较熟悉，但是转用 Gekko，因为后者似乎有一些优势。不过，我 运行 陷入了一个问题，我将使用我想象中的苹果园来说明这一问题。为期 5 周的收获期 (#horizon: 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]。我的储物空间 space 可容纳 6 个苹果，因此我可以提前计划并在最佳时机销售苹果，从而最大限度地提高收入。我尝试使用以下模型确定最佳时间表：
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] to 3.0, instead of 2.0（即等于 orchard[0]，该模型确实产生了成功的解决方案。
- 为什么会这样？
- 即使是“成功”的输出值也有点奇怪：存储 space 一次都没有使用，并且
storage_out 在最后一个时间步中没有得到适当的约束。显然，我没有正确地制定约束条件。我应该怎么做才能获得与 gurobi 输出相当的真实结果（请参阅下面的代码）？
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
用 demand = [3.0, 4.0, 2.0, 4.0, 2.0] 求解的 Gurobi 模型。请注意，gurobi 还会生成 demand = [2.0, 4.0, 2.0, 4.0, 2.0] 的解决方案。这对结果的影响很小：在 t=0 时售出的 n 个苹果变成了 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
您可以通过以下方式获得成功的解决方案：
m.options.NODES=2
问题是它正在用NODES=3求解主节点之间的平衡方程。您的微分方程有线性解，因此 NODES=2 应该足够准确。
以下是改进解决方案的其他几种方法：
- 对将库存移入或移出存储设置一个小惩罚。否则，求解器可以使用
storage_in = storage_out. 找到较大的任意值
- 我用了
m.Minimize(1e-6*storage_in)和m.Minimize(1e-6*storage_out)。
- 因为初始条件通常是固定的，所以我在开始时使用零值只是为了确保计算出第一个点。
- 如果它们以整数单位出售和存储，我也切换到整数变量。如果你想要
SOLVER=1. 的整数解，你需要切换到 APOPT 求解器
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)
我对 Gurobi 比较熟悉，但是转用 Gekko，因为后者似乎有一些优势。不过，我 运行 陷入了一个问题，我将使用我想象中的苹果园来说明这一问题。为期 5 周的收获期 (#horizon: 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]。我的储物空间 space 可容纳 6 个苹果，因此我可以提前计划并在最佳时机销售苹果，从而最大限度地提高收入。我尝试使用以下模型确定最佳时间表：
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] to 3.0, instead of 2.0（即等于 orchard[0]，该模型确实产生了成功的解决方案。
- 为什么会这样？
- 即使是“成功”的输出值也有点奇怪：存储 space 一次都没有使用，并且
storage_out在最后一个时间步中没有得到适当的约束。显然，我没有正确地制定约束条件。我应该怎么做才能获得与 gurobi 输出相当的真实结果（请参阅下面的代码）？
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
用 demand = [3.0, 4.0, 2.0, 4.0, 2.0] 求解的 Gurobi 模型。请注意，gurobi 还会生成 demand = [2.0, 4.0, 2.0, 4.0, 2.0] 的解决方案。这对结果的影响很小：在 t=0 时售出的 n 个苹果变成了 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
您可以通过以下方式获得成功的解决方案：
m.options.NODES=2
问题是它正在用NODES=3求解主节点之间的平衡方程。您的微分方程有线性解，因此 NODES=2 应该足够准确。
以下是改进解决方案的其他几种方法：
- 对将库存移入或移出存储设置一个小惩罚。否则，求解器可以使用
storage_in = storage_out. 找到较大的任意值
- 我用了
m.Minimize(1e-6*storage_in)和m.Minimize(1e-6*storage_out)。 - 因为初始条件通常是固定的，所以我在开始时使用零值只是为了确保计算出第一个点。
- 如果它们以整数单位出售和存储，我也切换到整数变量。如果你想要
SOLVER=1. 的整数解，你需要切换到 APOPT 求解器
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)