R Двоичная целочисленная оптимизация с группами

Я пытаюсь заставить Rsolnp ограничить мои параметры двоичными целыми числами или десятичными числами, которые почти одинаковы (например, .999 достаточно близко к 1).

У меня есть три вектора равной длины (52), каждый из которых будет умножен на мой двоичный вектор параметров в моей целевой функции.

library(Rsolnp)
a <- c(251, 179, 215, 251, 63, 45, 54, 63, 47, 34, 40, 47, 141, 101, 121, 141, 47, 34, 40, 47, 94, 67, 81, 94, 47, 34, 40, 47, 157, 108, 133, 157, 126, 85, 106, 126, 126, 85, 106, 126, 110, 74, 92, 110, 110, 74, 92, 110, 63, 40, 52, 63)
b <- c(179, 251, 215, 0, 45, 63, 54, 0, 34, 47, 40, 0, 101, 141, 121, 0, 34, 47, 40, 0, 67, 94, 81, 0, 34, 47, 40, 0, 108, 157, 133, 0, 85, 126, 106, 0, 85, 126, 106, 0, 74, 110, 92, 0, 74, 110, 92, 0, 40, 63, 52, 0)
c <- c(179, 179, 118, 179, 45, 45, 30, 45, 34, 34, 22, 34, 101, 101, 67, 101, 34, 34, 22, 34, 67, 67, 44, 67, 34, 34, 22, 34, 108, 108, 71, 108, 85, 85, 56, 85, 85, 85, 56, 85, 74, 74, 49, 74, 74, 74, 49, 74, 40, 40, 27, 40)

x - мой вектор параметров и ниже, если моя целевая функция.

objective_function = function(x){
  -(1166 *  sum(x[1:52] * a) / 2000) *
    (((sum(x[1:52] * b)) / 2100) + .05) *
    (((sum(x[1:52] * c))/1500) + 1.5)
}

По сути, я хочу, чтобы 1 параметр в каждой группе из 4 равнялся 1, а остальные 0, и я не уверен, как создать для этого правильные ограничения, но я считаю, что мне нужно использовать эти ограничения суммы в сочетании с другим типом ограничения. . Вот мои ограничения:

eqn1=function(x){
  z1=sum(x[1:4])
  z2=sum(x[5:8])
  z3=sum(x[9:12])
  z4=sum(x[13:16])
  z5=sum(x[17:20])
  z6=sum(x[21:24])
  z7=sum(x[25:28])
  z8=sum(x[29:32])
  z9=sum(x[33:36])
  z10=sum(x[37:40])
  z11=sum(x[41:44])
  z12=sum(x[45:48])
  z13=sum(x[49:52])
  return(c(z1,z2,z3,z4,z5,z6,z7,z8,z9,z10,z11,z12,z13))
}

И, наконец, вот вызов моей функции:

opti<-solnp(pars=rep(1,52), fun = objective_function, eqfun = eqn1, eqB = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), LB=rep(0,52))

Вызов opti $ pars возвращает мой вектор решения:

 [1] 7.199319e-01 2.800680e-01 6.015388e-08 4.886578e-10 5.540961e-01 4.459036e-01 2.906853e-07 4.635970e-08 5.389325e-01
[10] 4.610672e-01 2.979195e-07 3.651954e-08 6.228346e-01 3.771652e-01 1.980380e-07 3.348488e-09 5.389318e-01 4.610679e-01
[19] 2.979195e-07 3.651954e-08 5.820231e-01 4.179766e-01 2.099869e-07 2.624076e-08 5.389317e-01 4.610680e-01 2.979195e-07
[28] 3.651954e-08 6.499878e-01 3.500120e-01 1.959133e-07 1.059012e-08 6.249098e-01 3.750900e-01 2.588037e-07 1.752927e-08
[37] 6.249106e-01 3.750892e-01 2.588037e-07 1.752927e-08 6.095743e-01 3.904254e-01 2.741968e-07 2.233806e-08 6.095743e-01
[46] 3.904254e-01 2.741968e-07 2.233806e-08 5.679608e-01 4.320385e-01 6.821224e-07 3.997882e-08

Как можно видеть, вес разделяется между несколькими переменными в каждой группе по 4, вместо того, чтобы быть принудительно равным 1, а остальные равны 0.

Если это невозможно с этим пакетом, может ли кто-нибудь показать мне, как преобразовать мою целевую функцию для работы с другими пакетами оптимизации? Из того, что я видел, они требуют, чтобы целевая функция была преобразована в вектор коэффициентов. Любая помощь приветствуется. Спасибо!


person throwawaydisplayname    schedule 25.08.2020    source источник
comment
Похоже на проблему с нелинейной целевой функцией, линейными ограничениями и двоичными переменными. Это называется проблемой MINLP (смешанное целочисленное нелинейное программирование).   -  person Erwin Kalvelagen    schedule 25.08.2020
comment
Возможно, вы захотите взглянуть на гуроби.   -  person Christoph    schedule 25.08.2020
comment
Пожалуйста, укажите минимум! воспроизводимый пример.   -  person Christoph    schedule 25.08.2020
comment
@Christoph Вы должны иметь возможность запустить предоставленный код. В идеале мой вектор решения будет выглядеть примерно как 0,1,0,0,1,0,0,0 ... с каждым четвертым параметром, равным 1.   -  person throwawaydisplayname    schedule 25.08.2020


Ответы (3)


в CPLEX вы можете попробовать математическое программирование, как писал Пол, но вы также можете использовать Программирование ограничений.

В OPL (язык моделирования CPLEX)

using CP;

execute
{
  cp.param.timelimit=5; // time limit 5 seconds
  
}

int n=52;
range r=1..n;

int a[r]=[251, 179, 215, 251, 63, 45, 54, 63, 47, 34, 40, 47, 141, 101, 121, 141, 47, 34, 40, 47, 94, 67, 
81, 94, 47, 34, 40, 47, 157, 108, 133, 157, 126, 85, 106, 126, 126, 
85, 106, 126, 110, 74, 92, 110, 110, 74, 92, 110, 63, 40, 52, 63];
int b[r]=[179, 251, 215, 0, 45, 63, 54, 0, 34, 47, 40, 0, 101, 141, 121, 0, 
34, 47, 40, 0, 67, 94, 81, 0, 34, 47, 40, 0, 108, 157, 133, 0, 85, 126, 106, 0, 85, 
126, 106, 0, 74, 110, 92, 0, 74, 110, 92, 0, 40, 63, 52, 0];
int c[r]=[179, 179, 118, 179, 45, 45, 30, 45, 34, 34, 22, 34, 101, 101, 67, 101, 
34, 34, 22, 34, 67, 67, 44, 67, 34, 34, 22, 34, 108, 108, 71, 108, 85, 85, 56, 85,
 85, 85, 56, 85, 74, 74, 49, 74, 74, 74, 49, 74, 40, 40, 27, 40];

// decision variable 
 dvar boolean x[r];
 
 // objective
 dexpr float obj=
 
  -(1166 *  sum(i in r) (x[i]*a[i]) / 2000) *
    (((sum(i in r) (x[i]* b[i])) / 2100) + .05) *
    (((sum(i in r) (x[i]*c[i]))/1500) + 1.5);
    
minimize obj;

subject to
{
  // one and only one out of 4 is true
  forall(i in 1..n div 4) count(all(j in 1+(i-1)*4..4+(i-1)*4)x[j],1)==1;
  
} 

дает

// solution with objective -889.3463
x = [1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1
         0 0 0 0 1 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 1 0
         0 0];

в течение 5 секунд

NB: вы можете вызвать OPL CPLEX из R или полагаться на любой другой CPLEX API.

И на питоне можно написать то же самое

from docplex.cp.model import CpoModel

n=52

r=range(0,n)

a =[251, 179, 215, 251, 63, 45, 54, 63, 47, 34, 40, 47, 141, 101, 121, 141, 47, 34, 40, 47, 94, 67, 81, 94, 47, 34, 40, 47, 157, 108, 133, 157, 126, 85, 106, 126, 126, 85, 106, 126, 110, 74, 92, 110, 110, 74, 92, 110, 63, 40, 52, 63]
b =[179, 251, 215, 0, 45, 63, 54, 0, 34, 47, 40, 0, 101, 141, 121, 0, 34, 47, 40, 0, 67, 94, 81, 0, 34, 47, 40, 0, 108, 157, 133, 0, 85, 126, 106, 0, 85, 126, 106, 0, 74, 110, 92, 0, 74, 110, 92, 0, 40, 63, 52, 0]
c =[179, 179, 118, 179, 45, 45, 30, 45, 34, 34, 22, 34, 101, 101, 67, 101, 34, 34, 22, 34, 67, 67, 44, 67, 34, 34, 22, 34, 108, 108, 71, 108, 85, 85, 56, 85, 85, 85, 56, 85, 74, 74, 49, 74, 74, 74, 49, 74, 40, 40, 27, 40]

mdl = CpoModel(name='x')

#decision variables
mdl.x = {i: mdl.integer_var(0,n,name="x"+str(i+1)) for i in r}

mdl.minimize(-1166 *  sum(mdl.x[i]*a[i] / 2000 for i in r) \
             *((sum(mdl.x[i]* b[i] / 2100  for i in r) +0.05)) \
             *((sum(mdl.x[i]*c[i]/1500  for i in r) +1.5)) )

for i in range(0,n // 4):
    mdl.add(1==sum( mdl.x[j] for j in range(i*4+0,i*4+4)))

msol=mdl.solve(TimeLimit=5)

# Dislay solution
for i in r:
    if (msol[mdl.x[i]]==1):
        print(i+1," ")

и это дает

! Best objective         : -889.3464 
 
1  
5  
9  
13  
17  
22  
25  
30  
34  
38  
41  
45  
49  
person Alex Fleischer    schedule 29.08.2020
comment
Понятно! Спасибо Алекс - person throwawaydisplayname; 31.08.2020

Я пробовал с несколькими решателями. С помощью решателей MINLP Couenne и Baron мы можем решить эту проблему напрямую. С помощью Гуроби нам нужно разделить цель на две квадратичные части. Все эти решатели дают:

----    119 VARIABLE x.L  

i1  1.000,    i5  1.000,    i9  1.000,    i14 1.000,    i17 1.000,    i21 1.000,    i25 1.000,    i29 1.000
i34 1.000,    i38 1.000,    i41 1.000,    i46 1.000,    i49 1.000


----    119 VARIABLE z.L                   =     -889.346  obj

Здесь нули не печатаются.

Я использовал GAMS (коммерческий), но если вы хотите использовать бесплатные инструменты, вы можете использовать Pyomo (Python) + Couenne. Я не уверен насчет решателей MINLP для R, но Gurobi можно использовать из R.

Обратите внимание, что групповое ограничение просто:

groups(g)..  sum(group(g,i),x(i)) =e= 1;      

где g - это группы, а group(g,i) - это 2-й набор с отображением между группами и элементами.

Для Гуроби нужно сделать что-то вроде (в псевдокоде):

z1 = 1166 *  sum(i,x(i)*a(i)) / 2000        (linear)
z2 = ((sum(i, x(i)*b(i))) / 2100) + .05     (linear) 
z3 = ((sum(i, x(i)*c(i)))/1500) + 1.5       (linear) 
z23 = z2*z3                                 (non-convex quadratic)
obj = -z1*z23                               (non-convex quadratic)

и скажите Gurobi использовать невыпуклый решатель MIQCP.

Извините, для этого нет кода R. Но это может дать вам повод задуматься.

person Erwin Kalvelagen    schedule 25.08.2020

Я настроил ноутбук R, чтобы решить (или попытаться решить) проблему как смешанную целочисленную линейную программу, используя CPLEX в качестве решателя MIP и пакет Rcplex в качестве интерфейса к нему. Результаты были не впечатляющими. После пяти минут измельчения CPLEX получил решение, несколько худшее по сравнению с тем, что получил Эрвин (-886,8748 против его -889,346) с разрывом более 146% (что, учитывая результат Эрвина, в основном, это просто верхняя граница, сходящаяся очень медленно). Я счастлив поделиться ноутбуком, в котором показана линеаризация, но для его использования вам потребуется установить CPLEX.

Обновление: у меня есть второй ноутбук, использующий пакет генетических алгоритмов GA, который постоянно приближается к решению Эрвина (а иногда и достигает его) менее чем за пять секунд. Результаты случайны, поэтому повторный запуск может быть лучше (или хуже), и нет никаких доказательств оптимальности.

person prubin    schedule 28.08.2020
comment
В настоящее время я работаю над этим на python, но если бы вы могли поделиться со мной обоими этими блокнотами, было бы здорово @prubin - person throwawaydisplayname; 30.08.2020
comment
Я обновил свой ответ с помощью кода cplex python, который я использовал - person Alex Fleischer; 30.08.2020
comment
Первый блокнот (с использованием CPLEX): rubin.msu.domains/blog/groupselect.nb.html. - person prubin; 31.08.2020
comment
Вторая записная книжка (с использованием GA): rubin.msu.domains/blog/groupselect2.nb.html. - person prubin; 31.08.2020