I have an optimization problem that I am trying to model in mathematica (wolfram cloud, to be precise).
Link to notebook for code described here: https://www.dropbox.com/s/4jet4t76fe5pysk/ReFormulate.nb?dl=0
Here is the problem: given two sets of 10 variables, which can each take a value of -1,0, or 1 (later to be extended to higher integer values), and that each variable has two associated properties given by coefficients V and T, pick 2 variables from each set such that 1 variable from each set = 1 and 1 variable from each set = -1, and all other others = 0, to maximize the following objective function (where x_i represents the entire domain of variables): obj = (-Abs(Sum(x_i * V_i))) / Sum(x_i * T_i)
. Further, assume / constrain Sum(x_i * T_i) > 0
.
V and T coefficients represented by the following lists:
VX = {158.725,158.987,149.428,145.506,141.073,136.134,130.737,124.923,112.176,105.27}
VY = {130.136,127.017,123.847,120.74,117.654,17.3212,12.2567,15.1606,18.091,18.8637}
TX = {-69.3633,-68.3942,-60.192,-56.6941,-53.149,-49.5749,-46.0265,-42.5295,-35.7162,-32.3862}
TY = {-65.2871,-64.4114,-63.4106,-62.5911,-61.6879,-14.9112,-11.1885,-13.3237,-15.4669,-16.0089}
For my model, I have created 2 variables for each set of variables, one to represent a positive value and one to represent a negative value:
domain = Join[Table[Subscript[xpos,i],{i,10}],Table[Subscript[xneg,i],{i,10}],Table[Subscript[ypos,i],{i,10}],Table[Subscript[yneg,i],{i,10}]]
To constraint the problem as described, I created the following constraints
- Sum of all xpos == 1
- Sum of all xneg == -1
- Sum of all ypos == 1
- Sum of all yneg == -1
- All xpos, ypos >= 0
- All xneg, yneg <= 0
- Sum(domainvar_i * T_i) >= 0
Here is the code I used to do this:
A = {Join[Table[1,{i,10}],Table[0,{i,30}]],Join[Table[0,{i,10}],Table[1,{i,10}], Table[0,{i,20}]], Join[Table[0,{i,20}], Table[1,{i,10}],Table[0,{i,10}]], Join[Table[0,{i,30}], Table[1,{i,10}]]}
b = {1,-1,1,-1}
c = Join[VX, VX,VY,VY];
d = Join[TX,TX,TY,TY];
constraintSets = Map[A[[#]] * domain &, Range[4]]
constraints = Map[{Apply[Plus, constraintSets[[#]]] == b[[#]]} &, Range[4]]
conPos = Map[# >=0&,Join[Take[domain,10],Take[domain,{21,30}]]]
conNeg = Map[# <=0&,Join[Take[domain,{11,20}],Take[domain,{31,40}]]]
conDenom = {Apply[Plus,d*domain] >0}
Then, I use NMaximize
to get optimized result:
NMaximize[
Join[obj, constraints, conPos, conNeg, conDenom, {domain ∈ Integers}],
domain]
This is where I run into issues. Instead of returning an answer with integer results, I get the following error:
Maximize::incst: NMaximize was unable to generate any initial points satisfying the inequality constraints ....The initial region specified may not contain any feasible points. Changing the initial region or specifying explicit initial points may provide a better solution.
I tried changing the constraints to bound the variables to a finite range. I get an integer result in this case that satisfies most but not all the constraints:
boundNeg = conNeg =
Map[-5 <= # <=0&, Join[Take[domain,{11, 20}], Take[domain, {31, 40}]]]
boundPos =
Map[5 >= # >= 0&, Join[Take[domain, 10], Take[domain, {21, 30}]]]
NMaximize[
Join[obj, constraints, boundPos, boundNeg, conDenom, {domain ∈ Integers}],
domain]
This results in the error above, followed by a result that satisfies all but 3 constraints:
Maximize::nosat: Obtained solution does not satisfy the following constraints within Tolerance -> 0.001: {-xpos_1 ≤ 0, -ypos_1 ≤ 0,yneg_1 ≤0}
I further tried running the above model with bounded constraints with different working precision parameters. While the results varied, none produced a result that satisfied all constraints.
I am puzzled by this, as it seems to me fairly trivial to pick an instance that would satisfy the constraints. For example:
- xpos_2 = 1
- xneg_1 = -1
- ypos_2 = 1
- yneg_1 = -1
- all other variables = 0
- check: -TX_1 + TX_2 - TY_1 + TY_2 = 69.3633 - 68.3942 + 65.2871 - 64.4114 = 1.8448 > 0
I have tried to specify an initial region using the 4 values above as starting values:
initialRegion =
Join[
{domain[[1]]}, {{domain[[2]], 1, 2}}, Take[domain, {3, 10}],
{{domain[[11]], -2, -1}}, Take[domain, {12, 21}], {domain[[22]], 1, 2}},
Take[domain, {23, 30}], {{domain[[31]], -2, -1}}, Take[domain, {32, 40}]]
NMaximize[
Join[obj, constraints, boundPos, boundNeg, conDenom, {domain ∈ Integers}],
initialRegion]
However, this returns the same message:
Maximize::incst: NMaximize was unable to generate any initial points satisfying the inequality constraints....The initial region specified may not contain any feasible points. Changing the initial region or specifying explicit initial points may provide a better solution.
Maximize::nosat: Obtained solution does not satisfy the following constraints within Tolerance -> 0.001: {-xpos_1 ≤0, -ypos_1 ≤0, xneg_1 ≤0}
I am looking for some insight into how to solve this problem. Is this an issue with my model, or the way I am utilizing NMaximize
? Would this be a problem better suited to expressing with LinearProgramming
?
Link to notebook for code described here: https://www.dropbox.com/s/4jet4t76fe5pysk/ReFormulate.nb?dl=0
NMaximize
is using differential evolution here, which is meant for real-valued problems. That it does or doesn't give reasonable results for integer-valued problems is really more down to chance than design. Forcing the parameters to be integer does make the problem more difficult from the perspective of differential evolution, so you may have to play with the algorithmic parameters, which are to some (large) extent problem-specific. I described them here. $\endgroup$