R 优化最大值 buy/sell 取决于库存水平
R optimisation max buy/sell dependent on stock level
我想找到优化问题的解决方案。目的是通过低价买入和高价卖出来实现利润最大化。存在诸如最大库存水平和最大 buy/sell 单位数量等限制。此外，买卖限制取决于库存水平。我问了一个类似的问题，尽管这里没有最后一个条件。
这是一个例子：
price = c(12, 11, 12, 13, 16, 17, 18, 17, 18, 16, 17, 13)
capacity = 25
max_units_buy_30 = 4 # when inventory level is lower then 30% it is possible to buy 0 to 4 units
max_units_buy_65 = 3 # when inventory level is between 30% and 65% it is possible to buy 0 to 3 units
max_units_buy_100 = 2 # when inventory level is between 65% and 100% it is possible to buy 0 to 2 units
max_units_sell_30 = 4 # when inventory level is lower then 30% it is possible to sell 0 to 4 units
max_units_sell_70 = 6 # when inventory level is between 30% and 70% it is possible to sell 0 to 6 units
max_units_sell_100 = 8 # when inventory level is between 70% and 100% it is possible to sell 0 to 8 units
这里发生了很多事情。
- 描述
描述中似乎有问题。 “最大卖出价/价格 取决于库存水平。”这似乎是错误的。从数据来看，价格似乎是恒定的，但买卖限制取决于库存水平。
- 时间
把握好时机很重要。通常，我们将buy和sell视为t期间发生的事情（我们称它们为流量变量） . inv 是一个 股票变量 ，在期末 t 进行测量。说 sell[t] 和 buy[t] 依赖于 inv[t] 有点奇怪（我们在时间上倒退）。当然，我们可以对其建模并求解（我们求解联立方程，所以我们可以做这些事情）。但是，这在现实世界中可能没有意义。可能我们应该查看 inv[t-1] 以更改 buy[t] 和 sell[t]。
- 细分库存水平。
我们需要将库存水平分成几个部分。我们有以下部分：
0%-30%
30%-65%
65%-70%
70%-100%
我们将一个二进制变量与每个段相关联：
inventory in [0%-30%] <=> δ[1,t] = 1, all other zero
[30%-65%] δ[2,t] = 1
[65%-70%] δ[3,t] = 1
[70%-100%] δ[4,t] = 1
因为我们需要在所有时间段都这样做，所以我们增加了一个额外的索引 t。警告：我们将 δ[k,t] 与期初的库存关联起来，即 inv[t-1]。我们可以 link δ[k,t] 到 inv[t-1] 通过根据我们所在的段更改下限和上限。
- 边界 buy/sell
类似于库存的界限，我们有以下买卖上限：
segment buy sell
0%-30% 4 4
30%-65% 3 6
65%-70% 2 6
70%-100% 2 8
第一步是建立数学模型。这里发生的事情太多了，我们可以立即编写代码。数学模型就是我们的"design"。所以我们开始吧：
有了这个，我们可以开发一些R代码。这里我们使用CVXR作为建模工具，GLPK作为MIP求解器。
> library(CVXR)
>
> # data
> price = c(12, 11, 12, 13, 16, 17, 18, 17, 18, 16, 17, 13)
> capacity = 25
> max_units_buy = 4
> max_units_sell = 8
>
> # capacity segments
> s <- c(0,0.3,0.65,0.7,1)
>
> # corresponding lower and upper bounds
> invlb <- s[1:(length(s)-1)] * capacity
> invlb
[1] 0.00 7.50 16.25 17.50
> invub <- s[2:length(s)] * capacity
> invub
[1] 7.50 16.25 17.50 25.00
>
> buyub <- c(4,3,2,2)
> sellub <- c(4,6,6,8)
>
> # number of time periods
> NT <- length(price)
> NT
[1] 12
>
> # number of capacity segments
> NS <- length(s)-1
> NS
[1] 4
>
> # Decision variables
> inv = Variable(NT,integer=T)
> buy = Variable(NT,integer=T)
> sell = Variable(NT,integer=T)
> delta = Variable(NS,NT,boolean=T)
>
> # Lag operator
> L = cbind(rbind(0,diag(NT-1)),0)
>
> # optimization model
> problem <- Problem(Maximize(sum(price*(sell-buy))),
+ list(inv == L %*% inv + buy - sell,
+ sum_entries(delta,axis=2)==1,
+ L %*% inv >= t(delta) %*% invlb,
+ L %*% inv <= t(delta) %*% invub,
+ buy <= t(delta) %*% buyub,
+ sell <= t(delta) %*% sellub,
+ inv >= 0, inv <= capacity,
+ buy >= 0, sell >= 0))
> result <- solve(problem,verbose=T)
GLPK Simplex Optimizer, v4.47
120 rows, 84 columns, 369 non-zeros
0: obj = 0.000000000e+000 infeas = 1.200e+001 (24)
* 23: obj = 0.000000000e+000 infeas = 0.000e+000 (24)
* 85: obj = -9.875986758e+001 infeas = 0.000e+000 (2)
OPTIMAL SOLUTION FOUND
GLPK Integer Optimizer, v4.47
120 rows, 84 columns, 369 non-zeros
84 integer variables, 48 of which are binary
Integer optimization begins...
+ 85: mip = not found yet >= -inf (1; 0)
+ 123: >>>>> -8.800000000e+001 >= -9.100000000e+001 3.4% (17; 0)
+ 126: >>>>> -9.000000000e+001 >= -9.100000000e+001 1.1% (9; 11)
+ 142: mip = -9.000000000e+001 >= tree is empty 0.0% (0; 35)
INTEGER OPTIMAL SOLUTION FOUND
> cat("status:",result$status)
status: optimal
> cat("objective:",result$value)
objective: 90
> print(result$getValue(buy))
[,1]
[1,] 3
[2,] 4
[3,] 4
[4,] 3
[5,] 3
[6,] 1
[7,] 0
[8,] 0
[9,] 0
[10,] 4
[11,] 0
[12,] 0
> print(result$getValue(sell))
[,1]
[1,] 0
[2,] 0
[3,] 0
[4,] 0
[5,] 0
[6,] 0
[7,] 8
[8,] 6
[9,] 4
[10,] 0
[11,] 4
[12,] 0
> print(result$getValue(inv))
[,1]
[1,] 3
[2,] 7
[3,] 11
[4,] 14
[5,] 17
[6,] 18
[7,] 10
[8,] 4
[9,] 0
[10,] 4
[11,] 0
[12,] 0
> print(result$getValue(delta))
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,] 1 1 1 0 0 0 0 0 1 1 1 1
[2,] 0 0 0 1 1 0 0 1 0 0 0 0
[3,] 0 0 0 0 0 1 0 0 0 0 0 0
[4,] 0 0 0 0 0 0 1 0 0 0 0 0
>
所以，我认为有人为此欠我一瓶好白兰地。
我想找到优化问题的解决方案。目的是通过低价买入和高价卖出来实现利润最大化。存在诸如最大库存水平和最大 buy/sell 单位数量等限制。此外，买卖限制取决于库存水平。我问了一个类似的问题，尽管这里没有最后一个条件
这是一个例子：
price = c(12, 11, 12, 13, 16, 17, 18, 17, 18, 16, 17, 13)
capacity = 25
max_units_buy_30 = 4 # when inventory level is lower then 30% it is possible to buy 0 to 4 units
max_units_buy_65 = 3 # when inventory level is between 30% and 65% it is possible to buy 0 to 3 units
max_units_buy_100 = 2 # when inventory level is between 65% and 100% it is possible to buy 0 to 2 units
max_units_sell_30 = 4 # when inventory level is lower then 30% it is possible to sell 0 to 4 units
max_units_sell_70 = 6 # when inventory level is between 30% and 70% it is possible to sell 0 to 6 units
max_units_sell_100 = 8 # when inventory level is between 70% and 100% it is possible to sell 0 to 8 units
这里发生了很多事情。
- 描述
描述中似乎有问题。 “最大卖出价/价格 取决于库存水平。”这似乎是错误的。从数据来看，价格似乎是恒定的，但买卖限制取决于库存水平。
- 时间
把握好时机很重要。通常，我们将buy和sell视为t期间发生的事情（我们称它们为流量变量） . inv 是一个 股票变量 ，在期末 t 进行测量。说 sell[t] 和 buy[t] 依赖于 inv[t] 有点奇怪（我们在时间上倒退）。当然，我们可以对其建模并求解（我们求解联立方程，所以我们可以做这些事情）。但是，这在现实世界中可能没有意义。可能我们应该查看 inv[t-1] 以更改 buy[t] 和 sell[t]。
- 细分库存水平。
我们需要将库存水平分成几个部分。我们有以下部分：
0%-30%
30%-65%
65%-70%
70%-100%
我们将一个二进制变量与每个段相关联：
inventory in [0%-30%] <=> δ[1,t] = 1, all other zero
[30%-65%] δ[2,t] = 1
[65%-70%] δ[3,t] = 1
[70%-100%] δ[4,t] = 1
因为我们需要在所有时间段都这样做，所以我们增加了一个额外的索引 t。警告：我们将 δ[k,t] 与期初的库存关联起来，即 inv[t-1]。我们可以 link δ[k,t] 到 inv[t-1] 通过根据我们所在的段更改下限和上限。
- 边界 buy/sell
类似于库存的界限，我们有以下买卖上限：
segment buy sell
0%-30% 4 4
30%-65% 3 6
65%-70% 2 6
70%-100% 2 8
第一步是建立数学模型。这里发生的事情太多了，我们可以立即编写代码。数学模型就是我们的"design"。所以我们开始吧：
有了这个，我们可以开发一些R代码。这里我们使用CVXR作为建模工具，GLPK作为MIP求解器。
> library(CVXR)
>
> # data
> price = c(12, 11, 12, 13, 16, 17, 18, 17, 18, 16, 17, 13)
> capacity = 25
> max_units_buy = 4
> max_units_sell = 8
>
> # capacity segments
> s <- c(0,0.3,0.65,0.7,1)
>
> # corresponding lower and upper bounds
> invlb <- s[1:(length(s)-1)] * capacity
> invlb
[1] 0.00 7.50 16.25 17.50
> invub <- s[2:length(s)] * capacity
> invub
[1] 7.50 16.25 17.50 25.00
>
> buyub <- c(4,3,2,2)
> sellub <- c(4,6,6,8)
>
> # number of time periods
> NT <- length(price)
> NT
[1] 12
>
> # number of capacity segments
> NS <- length(s)-1
> NS
[1] 4
>
> # Decision variables
> inv = Variable(NT,integer=T)
> buy = Variable(NT,integer=T)
> sell = Variable(NT,integer=T)
> delta = Variable(NS,NT,boolean=T)
>
> # Lag operator
> L = cbind(rbind(0,diag(NT-1)),0)
>
> # optimization model
> problem <- Problem(Maximize(sum(price*(sell-buy))),
+ list(inv == L %*% inv + buy - sell,
+ sum_entries(delta,axis=2)==1,
+ L %*% inv >= t(delta) %*% invlb,
+ L %*% inv <= t(delta) %*% invub,
+ buy <= t(delta) %*% buyub,
+ sell <= t(delta) %*% sellub,
+ inv >= 0, inv <= capacity,
+ buy >= 0, sell >= 0))
> result <- solve(problem,verbose=T)
GLPK Simplex Optimizer, v4.47
120 rows, 84 columns, 369 non-zeros
0: obj = 0.000000000e+000 infeas = 1.200e+001 (24)
* 23: obj = 0.000000000e+000 infeas = 0.000e+000 (24)
* 85: obj = -9.875986758e+001 infeas = 0.000e+000 (2)
OPTIMAL SOLUTION FOUND
GLPK Integer Optimizer, v4.47
120 rows, 84 columns, 369 non-zeros
84 integer variables, 48 of which are binary
Integer optimization begins...
+ 85: mip = not found yet >= -inf (1; 0)
+ 123: >>>>> -8.800000000e+001 >= -9.100000000e+001 3.4% (17; 0)
+ 126: >>>>> -9.000000000e+001 >= -9.100000000e+001 1.1% (9; 11)
+ 142: mip = -9.000000000e+001 >= tree is empty 0.0% (0; 35)
INTEGER OPTIMAL SOLUTION FOUND
> cat("status:",result$status)
status: optimal
> cat("objective:",result$value)
objective: 90
> print(result$getValue(buy))
[,1]
[1,] 3
[2,] 4
[3,] 4
[4,] 3
[5,] 3
[6,] 1
[7,] 0
[8,] 0
[9,] 0
[10,] 4
[11,] 0
[12,] 0
> print(result$getValue(sell))
[,1]
[1,] 0
[2,] 0
[3,] 0
[4,] 0
[5,] 0
[6,] 0
[7,] 8
[8,] 6
[9,] 4
[10,] 0
[11,] 4
[12,] 0
> print(result$getValue(inv))
[,1]
[1,] 3
[2,] 7
[3,] 11
[4,] 14
[5,] 17
[6,] 18
[7,] 10
[8,] 4
[9,] 0
[10,] 4
[11,] 0
[12,] 0
> print(result$getValue(delta))
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,] 1 1 1 0 0 0 0 0 1 1 1 1
[2,] 0 0 0 1 1 0 0 1 0 0 0 0
[3,] 0 0 0 0 0 1 0 0 0 0 0 0
[4,] 0 0 0 0 0 0 1 0 0 0 0 0
>
所以，我认为有人为此欠我一瓶好白兰地。