2
$\begingroup$

I have a problem in vectors $\mathbf{x} = (x_1,\ldots,x_n)$ and $\mathbf{y}=(y_1,\ldots,y_n)$, where every $x_i,y_i\in\mathbb{R}_{\geq 0}$,

$$\begin{array}{ll} \text{maximize} & \displaystyle\sum_{i=1}^n {y_i}^u {x_i}^{1-u}\\ \text{subject to} & \displaystyle\sum_{i=1}^n x_i = 1\\ & \displaystyle\sum_{i=1}^n y_i = 1\\ & \mathbf{a} \leq \mathbf{x} \leq \mathbf{b}\\ & \mathbf{c} \leq \mathbf{y} \leq \mathbf{d}\end{array}$$

where $u \in (0,1)$ and nonnegative vectors $\mathbf{a}$, $\mathbf{b}$, $\mathbf{c}$, $\mathbf{d}$ are given (namely they are known and one can choose them freely as long as they satisfy the constraints) such that the following conditions are also satisfied: $$ \sum_{i=1}^n a_i<1,\ \sum_{i=1}^n c_i<1,\quad \mbox{and}\quad \sum_{i=1}^n b_i>1,\ \sum_{i=1}^n d_i>1 $$

Can I solve such a problem with Mathematica? I was thinking about using NMaximize or Maximize but I have no experience for optimization over vectors using Mathematica, thats why I have no Mathematica code.

Addendum: I think Daniel is right. Therefore I prepared the vectors $\mathbf{a}$, $\mathbf{b}$, $\mathbf{c}$, $\mathbf{d}$:

u=0.5;
a[i_] := 0.6*PDF[BinomialDistribution[20, 0.3], i]
b[i_] := 1.5*PDF[BinomialDistribution[20, 0.3], i]
c[i_] := 0.4*PDF[BinomialDistribution[20, 0.7], i]
d[i_] := 2*PDF[BinomialDistribution[20, 0.7], i]

I am still looking for answers which could work for $n>40$.

Improve this question
$\endgroup$
9
  • $\begingroup$ Try setting n to some small value (say 2, 3, or 4) and set up the Maximize function. Maybe from trying a few values of n you might see a pattern for a general answer. $\endgroup$
    – JimB
    Commented May 4, 2017 at 4:00
  • $\begingroup$ Do you want to maximize the above function for $u$ given $x$ and $y$? Or given $u$ choose the values of $x$ and $y$ that maximize the function? $\endgroup$
    – JimB
    Commented May 4, 2017 at 4:19
  • $\begingroup$ @JimBaldwin $u$ is given. one can take it as say $0.5$. The unknowns are $\mathbf{x}$ and $\mathbf{y}$ vectors. If I could have written a code for $n=10$, then this code would have also worked for every $n$.. $\endgroup$
    – Seyhmus Güngören
    Commented May 4, 2017 at 4:32
  • $\begingroup$ If $a=c=0$ and $b=d=1$, is the answer for any $u$: $x_i=y_i=1/n$ which results in the maximum always being 1 ? $\endgroup$
    – JimB
    Commented May 4, 2017 at 5:26
  • $\begingroup$ @JimBaldwin my best guess is that you are right. Moreover, I think that for any $a=c$ and $b=d$ no matter what they are equal to one would get $x_i=y_i\forall i$. This is probably the most trivial part of the problem. The non trivial part is when a,b,c,d unequal almost everywhere. $\endgroup$
    – Seyhmus Güngören
    Commented May 4, 2017 at 5:44

3 Answers 3

Reset to default
2
$\begingroup$ 
u = 0.5;
n = 20
a = Table[0.6*PDF[BinomialDistribution[20, 0.3], i], {i, n}];
b = Table[1.5*PDF[BinomialDistribution[20, 0.3], i], {i, n}];
c = Table[0.4*PDF[BinomialDistribution[20, 0.7], i], {i, n}];
d = Table[2*PDF[BinomialDistribution[20, 0.7], i], {i, n}];
X = Array[x, n];
Y = Array[y, n];

FindMaximum[{X^(1 - u).Y^u, 
  Flatten[{ Total[X] == 1 , Total[Y] == 1, 
    MapThread[#1 <= #2 <= #3 &, {a, X, b}], 
    MapThread[#1 <= #2 <= #3 &, {c, Y, d}]}]}  , Flatten[Join[{X, Y}]]]

This yields a result (0.287865) along with a warning that it is not converged, only a best estimate.

(Note for the given a,b,c,d, n must be a pretty large number )

Improve this answer
$\endgroup$
9
  • $\begingroup$ Looks good. I guess you posted while I was still writing mine up. $\endgroup$
    – Daniel Lichtblau
    Commented May 4, 2017 at 21:29
  • $\begingroup$ would it be possible to plot the vectors a, b,c ,d and the vectors x, y together? we can then compare different answers.. $\endgroup$
    – Seyhmus Güngören
    Commented May 4, 2017 at 21:39
  • $\begingroup$ what do you mean by $n$ must be a pretty large number? it is given as $20$ already. Am I missing something? $\endgroup$
    – Seyhmus Güngören
    Commented May 4, 2017 at 21:59
  • $\begingroup$ I just tried to change every $20$ with $50$ in your code and in this case. the solution does not satisfy the constraints. $\endgroup$
    – Seyhmus Güngören
    Commented May 5, 2017 at 0:14
  • $\begingroup$ sorry to be vague. I meant n can not be like 2 or 3. (due to Sum b,d >1 requirement) Beyond about 20 does you no good either because a,b,c,d all go to zero. You should get a solution but the extra terms must be zero. $\endgroup$
    – george2079
    Commented May 5, 2017 at 1:12
2
$\begingroup$

For the example provided one might proceed as follows.

u = 1/2;
a[i_] := 0.6*PDF[BinomialDistribution[20, 0.3], i]
b[i_] := 1.5*PDF[BinomialDistribution[20, 0.3], i]
c[i_] := 0.4*PDF[BinomialDistribution[20, 0.7], i]
d[i_] := 2*PDF[BinomialDistribution[20, 0.7], i]
n = 14;

Check the constraints on the totals.

In[415]:= Map[Total, {Array[a, n], Array[b, n], Array[c, n], 
  Array[d, n]}]

(* Out[415]= {0.599495482389, 1.49873870597, 0.233451668207, \
1.16725834104} *)

Create the x and y vectors and the set of constraints.

xvec = Array[x, n];
yvec = Array[y, n];
c1 = {Total[xvec] >= 1, Total[yvec] >= 1};
c2 = Thread[Array[a, n] <= xvec <= Array[b, n]];
c3 = Thread[Array[c, n] <= yvec <= Array[d, n]];
obj = xvec^u.yvec^(1 - u);
constraints = Join[c1, c2, c3];
vars = Join[xvec, yvec];

Notice that I constrain the x and y sums to be >= unity rather than strictly equal. This is to work around a problem NMinimize seems to have with imposing the equality constraint and still finding viable initial search points. Alternatively one can use FindMinimum and avoid the issue. Notice that >= suffices because any "minimum" value found that satisfies the strong inequality can be further reduced by making it satisfy the equality instead.

Now we minimize.

{min, vals} = NMinimize[{obj, constraints}, regions]

(* Out[600]= {0.113052286045, {x[1] -> 0.0102589870128, 
  x[2] -> 0.0417687898593, x[3] -> 0.107405488744, 
  x[4] -> 0.195631440342, x[5] -> 0.268294549424, 
  x[6] -> 0.141472396852, x[7] -> 0.0985572027817, 
  x[8] -> 0.0686380474937, x[9] -> 0.0392217407402, 
  x[10] -> 0.0184902491202, x[11] -> 0.00720399313413, 
  x[12] -> 0.00231556920486, x[13] -> 0.000610699578119, 
  x[14] -> 0.000130864199554, y[1] -> 6.50982441109*10^-10, 
  y[2] -> 1.44278086865*10^-8, y[3] -> 2.01986175439*10^-7, 
  y[4] -> 2.00302480394*10^-6, y[5] -> 0.0000149559109898, 
  y[6] -> 0.0000872428055708, y[7] -> 0.000407133069209, 
  y[8] -> 0.00154371284591, y[9] -> 0.00480266214711, 
  y[10] -> 0.012326832933, y[11] -> 0.0402198296988, 
  y[12] -> 0.228793477986, y[13] -> 0.328523969226, 
  y[14] -> 0.383277964372}} *)

Check constraints:

{c1, c2, c3} /. vals

(* Out[601]= {{True, True}, {True, True, True, True, True, True, True, 
  True, True, True, True, True, True, True}, {True, True, True, True, 
  True, True, True, True, True, True, True, True, True, True}} *)

FindMinimum gives a comparable result.

{min, vals} = FindMinimum[{obj, constraints}, vars]

(* Out[602]= {0.113052279616, {x[1] -> 0.0102590056662, 
  x[2] -> 0.0417688087857, x[3] -> 0.107405508307, 
  x[4] -> 0.19563146156, x[5] -> 0.268294575854, 
  x[6] -> 0.141472291161, x[7] -> 0.0985571911308, 
  x[8] -> 0.068638043823, x[9] -> 0.0392217393274, 
  x[10] -> 0.0184902485401, x[11] -> 0.00720399293769, 
  x[12] -> 0.00231556915854, x[13] -> 0.000610699558297, 
  x[14] -> 0.000130864191064, y[1] -> 6.50866427519*10^-10, 
  y[2] -> 1.4427539019*10^-8, y[3] -> 2.01985546167*10^-7, 
  y[4] -> 2.00302333266*10^-6, y[5] -> 0.0000149559075502, 
  y[6] -> 0.0000872427940428, y[7] -> 0.000407133038865, 
  y[8] -> 0.00154371277236, y[9] -> 0.00480266195846, 
  y[10] -> 0.0123268323601, y[11] -> 0.0402198257305, 
  y[12] -> 0.22879347941, y[13] -> 0.328523970434, 
  y[14] -> 0.383277965507}} *)

--- edit ---

Here is the n=20 case. I use equality constraints and FindMaximum.

u = 1/2;
a[i_] := 0.6*PDF[BinomialDistribution[20, 0.3], i]
b[i_] := 1.5*PDF[BinomialDistribution[20, 0.3], i]
c[i_] := 0.4*PDF[BinomialDistribution[20, 0.7], i]
d[i_] := 2*PDF[BinomialDistribution[20, 0.7], i]

n = 20;
Map[Total, {Array[a, n], Array[b, n], Array[c, n], Array[d, n]}]

(* Out[678]= {0.599521246402, 1.49880311601, 0.399999999986, \
1.99999999993} *)

xvec = Array[x, n];
yvec = Array[y, n];
c1 = {Total[xvec] == 1, Total[yvec] == 1};
c2 = Thread[Array[a, n] <= xvec <= Array[b, n]];
c3 = Thread[Array[c, n] <= yvec <= Array[d, n]];
obj = xvec^u.yvec^(1 - u);
constraints = Join[c1, c2, c3];
vars = Join[xvec, yvec];

{min, vals} = FindMaximum[{obj, constraints}, vars]

(* Out[688]= {0.287179170249, {x[1] -> 0.00420039062101, 
  x[2] -> 0.0168068363863, x[3] -> 0.0430675097985, 
  x[4] -> 0.0783747301246, x[5] -> 0.107512002157, 
  x[6] -> 0.162230688536, x[7] -> 0.246265726493, 
  x[8] -> 0.17156411577, x[9] -> 0.0980431454055, 
  x[10] -> 0.0462211272314, x[11] -> 0.0180081125686, 
  x[12] -> 0.00578813622943, x[13] -> 0.00152640343632, 
  x[14] -> 0.000326847202309, x[15] -> 0.0000559562031744, 
  x[16] -> 7.46270969759*10^-6, x[17] -> 7.43504383666*10^-7, 
  x[18] -> 5.41032712888*10^-8, x[19] -> 2.4407490807*10^-9, 
  x[20] -> 5.2301766015*10^-11, y[1] -> 3.2543321076*10^-9, 
  y[2] -> 7.21376950518*10^-8, y[3] -> 9.96710077075*10^-7, 
  y[4] -> 9.96943298*10^-6, y[5] -> 0.000074658665243, 
  y[6] -> 0.000435965674377, y[7] -> 0.00203522075526, 
  y[8] -> 0.00771750424298, y[9] -> 0.0240107719009, 
  y[10] -> 0.0616279225931, y[11] -> 0.130722583872, 
  y[12] -> 0.228736998482, y[13] -> 0.30013881215, 
  y[14] -> 0.0774548114721, y[15] -> 0.0716601191227, 
  y[16] -> 0.0522522505167, y[17] -> 0.0287165994173, 
  y[18] -> 0.0112102849437, y[19] -> 0.0028066714826, 
  y[20] -> 0.000387777172855}} *)

In[691]:= Map[Total, {xvec, yvec} /. vals]

(* Out[691]= {0.999999990973, 0.999999993999} *)

Also works fine with FindMinimum.

--- end edit ---

Improve this answer
$\endgroup$
11
  • $\begingroup$ would it be possible to plot the vectors a, b,c ,d and the vectors x, y together? we can then compare different answers. $\endgroup$
    – Seyhmus Güngören
    Commented May 4, 2017 at 21:39
  • $\begingroup$ the value of $n$ must be $20$. it is actually the length of the vectors $a,b,c,d$ as well as $x$ and $y$. $\endgroup$
    – Seyhmus Güngören
    Commented May 4, 2017 at 21:41
  • $\begingroup$ Where was it specified that n=20? I missed this detail. $\endgroup$
    – Daniel Lichtblau
    Commented May 4, 2017 at 22:40
  • $\begingroup$ $n$ is the length of the vectors $\mathbf{x}$ and $\mathbf{y}$. and these vectors are bounded by $\mathbf{a},\mathbf{b},\mathbf{c},\mathbf{d}$, which also have the same length. In the question their sum runs from $1$ to $n$ (see the $4$ constraints just before the end of the question), implying that their lengths are also $n$. In the given example, $\mathbf{a},\mathbf{b},\mathbf{c},\mathbf{d}$ have length $20$, so are $\mathbf{x}$ and $\mathbf{y}$. I am sorry if It was not clear enough in the question. $\endgroup$
    – Seyhmus Güngören
    Commented May 4, 2017 at 22:47
  • 2
    $\begingroup$ There are three responses with a total of two upvotes (one of which was mine). No response has been accepted. The actual problem is at best poorly elaborated and seems to keep changing, with no indication of what, if anything, has been tried by the poster. And now it is turning into a question about different programs, which is definitely outside the scope of this forum. I am of the opinion that this belongs on a different site. $\endgroup$
    – Daniel Lichtblau
    Commented May 7, 2017 at 21:41
0
$\begingroup$

I tried something but not sure the result is right

n = 2;
varx = Array[x, n];
vary = Array[y, n];
contsX = And @@ (a < # < b & /@ varx);
contsY = And @@ (c < # < d & /@ vary);
u = 0.4; a = 0; b = 1; c = 0; d = 1;
func = Sum[y[i]^u x[i]^(1 - u), {i, n}];
NMaximize[{func,  Sum[x[i], {i, n}] == 1 && Sum[y[i], {i, n}] == 1 && contsX &&    contsY}, Join[varx, vary]]

EDITED: 

varx = Array[x, n];
vary = Array[y, n];
constX = And @@ MapThread[#1 <= #2 <= #3 &, {a, varx, b}];
constY = And @@ MapThread[#1 <= #2 <= #3 &, {c, vary, d}];
func = Total[y[#]^u x[#]^(1 - u) & /@ Range@n]
FindMaximum[{func,   Total@varx == 1 && Total@vary == 1 && contsX && contsY},  Join[varx, vary]]
Improve this answer
$\endgroup$
1
  • $\begingroup$ I provided an example in the question. $a = 0; b = 1; c = 0; d = 1$ is the trivial case. $\endgroup$
    – Seyhmus Güngören
    Commented May 4, 2017 at 19:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.