Оптимизация R максимальная покупка / продажа в зависимости от уровня запасов

Здесь много чего происходит.

1. Описание

Похоже, в описании есть проблема. «Максимальная продажа / цена зависит от уровня запасов». Кажется, это неверно. Судя по данным, цена кажется постоянной, но ограничения на продажу и покупку зависят от уровня запасов.

2. Время

Важно правильно выбрать время. Обычно мы рассматриваем buy и sell как события, происходящие в период t (мы называем их переменными потока). inv - это биржевая переменная, которая измеряется в конце периода t. Сказать, что sell[t] и buy[t] зависят от inv[t], немного странно (мы идем назад во времени). Конечно, мы можем его смоделировать и решить (мы решаем как одновременные уравнения, поэтому мы можем делать это). Но это может не иметь смысла в реальном мире. Вероятно, нам следует посмотреть на inv[t-1], чтобы изменить buy[t] и sell[t].

3. Сегментирование уровней запасов.

Нам нужно разделить уровни запасов на сегменты. У нас есть следующие сегменты:

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] с инвентарем в начале периода t, то есть inv[t-1]. Мы можем связать δ[k,t] с inv[t-1], изменив нижнюю и верхнюю границы в зависимости от того, в каком сегменте мы находимся.

4. Границы на покупку / продажу

Как и в случае с ограничениями запасов, у нас есть следующие верхние границы для покупки и продажи:

     segment     buy   sell
     0%-30%       4     4 
     30%-65%      3     6
     65%-70%      2     6
     70%-100%     2     8

Первый шаг - разработать математическую модель. Здесь происходит слишком много всего, и мы можем сразу же все запрограммировать. Математическая модель - это наш «дизайн». Итак, начнем:

Благодаря этому мы можем разработать некоторый код 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
> 

Так что, думаю, кто-то должен мне за это хорошую бутылку коньяка.