# Using NMaximize to get integer results for a fractional programming problem

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]
constraints = Map[{Apply[Plus, constraintSets[[#]]] == b[[#]]} &, Range]
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[]}, {{domain[], 1, 2}}, Take[domain, {3, 10}],
{{domain[], -2, -1}}, Take[domain, {12, 21}], {domain[], 1, 2}},
Take[domain, {23, 30}], {{domain[], -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

• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. Sep 12, 2015 at 22:05
• Worth noting is that 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. Sep 12, 2015 at 23:32
• Might be better to recast as an ILP (go from fractional LP to ordinary LP, then throw in integrality constraints). Sep 13, 2015 at 20:43
• Dropping the integer constraints using the example found here (bit.ly/1gmRCX8) as a guide, I am able to transform the problem to an ordinary LP fairly easily. However, I am unsure how to go about adding integrality constraints, since transformation of the model is such that x_i = y_i / t, where x are variables in the LFP model, and y and t are variables in the LP model. See notebook here: bit.ly/1LuueQf. If I constrain y and to the integer domain, a) this does not guarantee that y / t will be an integer and b) when I test this, the solver is not able to satisfy all constraints Sep 14, 2015 at 19:18