2
$\begingroup$

I need help in finding a way to minimize a function with some complex constraints. I'm not an expert in this field, but I'm still trying to understand if I can do something with my problem.

So, here's the functions:

\begin{equation} \min z = \{ds_1*[(cp_{11}*qp_1*pp_{11}+...+(cp_{n1}*qp_n*pp_{n1})]+...+ds_m*[(cp_{1m}*qp_1*pp_{1m}+...+(cp_{nm}*qp_n*pp_{nm})]\} \end{equation}

with those given constants

\begin{equation} ds_x > 0, qp_x > 0, pp_{xy} > 0 \end{equation}

and

$$ cp_{yx} = \begin{cases} 1 \\ 0 \end{cases}, \begin{equation} \sum\limits_{i=1}^m cp_{iy} = 1 \text{ (only one element can be 1 in the whole column)} \end{equation} $$

what I have to do is to find the values for those cp terms to minimize the function.

Now, I know a little bit of optimization (simplex, branch & bound and stuff like this) but I've never went that deep and this one looks really complex to me.

Being a computer scientist, I thought of making some kind of backtracking stuff, but even with a small instance (let's say n=3 and m=7) I'd have to work on

\begin{equation} 2^{3*7}=2097152 \end{equation}

different combinations, which is really a lot.

Do you have any suggestion on this?

EDIT: Based on @DumpsterDoofus's answer, this is what I tried (with some real values)

\begin{equation} d=[7, 5, 11], \text{vector with the } ds_j \text{ constants} \end{equation}

\begin{equation} p= \begin{bmatrix} 5 & 2 & 5 & 3 & 7 & 1 & 1\\2 & 8 & 3 & 2 & 7 & 2 & 3\\3 & 6 & 4 & 5 & 9 & 1 & 1 \end{bmatrix} , \text{matrix with the } pp_{ij} \text{ constants} \end{equation}

assuming that

\begin{equation} qp_i = 1, \forall i \end{equation}

we have

\begin{equation} k= \begin{bmatrix} 7*5 & 7*2 & 7*5 & 7*3 & 7*7 & 7*1 & 7*1\\5*2 & 5*8 & 5*3 & 5*2 & 5*7 & 5*2 & 5*3\\11*3 & 11*6 & 11*4 & 11*5 & 11*9 & 11*1 & 11*1 \end{bmatrix} = \begin{bmatrix} 35 & 14 & 35 & 21 & 49 & 7 & 7\\10 & 40 & 15 & 10 & 35 & 10 & 15\\33 & 66 & 44 & 55 & 99 & 11 & 11 \end{bmatrix} \end{equation}

which, translated to a vecotr, becomes

\begin{equation} k=[35,10,33,14,40,66,35,15,44,21,10,55,49,35,99,7,10,11,7,15,11] \end{equation}

(I did it manually and now I corrected it after proper explanation of the vec() function)

Now that I've got the value of k, I'm ready to run this in Mathematica:

n=3;
m=7;
j[k_]:=ConstantArray[1,k];
k={35,10,33,14,40,66,35,15,44,21,10,55,49,35,99,7,10,11,7,15,11};
A=ArrayFlatten[{IdentityMatrix[m]\[TensorProduct]j[n]}];
b=j[m]\[TensorProduct]{1,0};
Round@LinearProgramming[k,A,b]

Output: {0,1,0,1,0,0,0,1,0,0,1,0,0,1,0,1,0,0,1,0,0} [Edited after changing k to its correct value]

Changing the output back to a matrix leads to this:

\begin{bmatrix} 0 & 1 & 0 & 0 & 0 & 1 & 1\\1 & 0 & 1 & 1 & 1 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix}

which is clearly wrong as we have one column full of zeroes and another one full of ones. seems correct now.

$\endgroup$
  • $\begingroup$ I'm in the process of writing a step-by-step answer, but one question: you say $\sum\limits_{i=1}^m cp_{yi} = 1 \text{ (only one element can be 1 in the whole column)}$, but a sum over the second index is actually doing the sum horizontally (ie, summing over each row), so shouldn't it actually be that "only one element can be 1 in the whole row", instead of column? $\endgroup$ – DumpsterDoofus Aug 3 '14 at 0:39
  • $\begingroup$ In my answer, I assume that your formula has a typo, and that you really meant to write $\sum\limits_{i=1}^m cp_{iy} = 1$, ie, that only one element can be 1 in the whole column. If this is incorrect, I can modify my answer to work for rows instead. $\endgroup$ – DumpsterDoofus Aug 3 '14 at 2:01
  • $\begingroup$ The reason you're getting an incorrect answer is because you're vectorizing is incorrect; vectorization is done by stacking the columns of a matrix, and vice versa for the reverse operation, see en.wikipedia.org/wiki/Vectorization_(mathematics). I'll edit my answer to include a vec function, and run it on the values you provided. $\endgroup$ – DumpsterDoofus Aug 3 '14 at 13:48
  • $\begingroup$ Ok, now I'm really confused. You say that $m=7$, but then you have $d=[7,5,11]$ with length 3, but that contradicts the formula for $z$ you provided, which states that the $ds_i$ values range in $i$ from $1$ to $m$. So shouldn't $d$ be of length $7$, and $qp$ be of length 3, based on the formula you provided? Or is that a typo in your formula, and the $m$'s and $n$'s should be swapped? $\endgroup$ – DumpsterDoofus Aug 3 '14 at 14:06
  • 1
    $\begingroup$ I updated my answer, and I got the same result as you did. Might want to plug in a couple other permutations and check if it's actually optimal, but it's probably correct. $\endgroup$ – DumpsterDoofus Aug 3 '14 at 14:24
2
$\begingroup$

First, note that since $ds_x,qp_x,pp_{xy}$ are constants, the expression $z$ can be simplified greatly: $$z=\{ds_1[(cp_{11}qp_1pp_{11}+...+(cp_{m1}qp_mpp_{m1})]+...+ds_n[(cp_{1n}qp_1pp_{1n}+...+(cp_{mn}qp_mpp_{mn})]\} \\ =\sum_{i,j}k_{ij}cp_{ij} \\ =\mathbf{k}\cdot\mathbf{c}$$ where $$\mathbf{c}=\text{vec}(cp) \\ \mathbf{k}=\text{vec}(k)$$ where $\text{vec}()$ is the matrix-to-vector operator, and $$k_{ij}=ds_iqp_jpp_{ij}.$$

Let $cp$ be an $n\times m$ matrix (as you describe in your question). The constraint equation (sum over columns of $cp$ is always 1) can be written as $$A\mathbf{c}=\mathbf{j}_m$$ where $\mathbf{j}_m$ is a list of $m$ ones, and $$A=\text{ArrayFlatten}[\{I_m\otimes\mathbf{j}_n\}]$$ where $I_m$ is the $m\times m$ identity matrix and $\otimes$ is the tensor product.$^1$

At that point, you're basically done, since you can then invoke LinearProgramming$^2$. Here's a very concise example (for $n=5,m=2$) which randomly generates coefficients $k_{ij}$ and then optimizes the resulting system:

n = 5;
m = 2;
j[k_] := ConstantArray[1, k];
k = RandomReal[{0, 1}, n m];
A = ArrayFlatten[{IdentityMatrix[m]\[TensorProduct]j[n]}];
b = j[m]\[TensorProduct]{1, 0};
Round@LinearProgramming[k, A, b]

Output: {0, 0, 0, 0, 1, 0, 0, 0, 0, 1}

The only confusing point is the definition $\mathbf{b}=\mathbf{j}_m\otimes\{1,0\}$, which is done because we want Mathematica to use the constraint $A\mathbf{c}=\mathbf{j}_m$, rather than the default syntax LinearProgramming[k, A, j[m]], which would use $A\mathbf{c}\geq\mathbf{j}_m$.$^3$

To get the answer for your particular problem, simply replace the line k = RandomReal[{0, 1}, n m] with your actual values for the problem you're trying to solve!

Minor notes:

$^1$: As an example, for the case $cp$ is $5\times 2$, you have $$A\mathbf{c}=\left( \begin{array}{cccccccccc} 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1 & 1 \\ \end{array} \right) \left( \begin{array}{c} cp_{11} \\ cp_{12} \\ cp_{13} \\ cp_{14} \\ cp_{15} \\ cp_{21} \\ cp_{22} \\ cp_{23} \\ cp_{24} \\ cp_{25} \\ \end{array} \right)=\left( \begin{array}{c} 1 \\ 1 \\ \end{array} \right)=\mathbf{j}_2.$$

As you can see, this ensures that the sum over each column of $cp$ is 1.

$^2$: The advantage of LinearProgramming (rather than Minimize) is that it's pretty well-suited to tackling very large problems, with a large number of variables and constraints, whereas Minimize doesn't scale well with problem size.

$^3$: From the documentation page from LinearProgramming:

enter image description here

By setting the $s_i=0$ you enforce constraint equality; a sneaky way to do this syntax is to use $\mathbf{b}=\mathbf{j}_m\otimes\{1,0\}$ instead of $\mathbf{b}=\mathbf{j}_m$.

Edit:

Based on the corrections you provided, here's the result for your specific numbers:

n = 3;
m = 7;
vec[m_] := Flatten[m\[Transpose]];
mat[v_] := Partition[v, n]\[Transpose];
j[k_] := ConstantArray[1, k];
d = {7, 5, 11};
p = ( {
    {5, 2, 5, 3, 7, 1, 1},
    {2, 8, 3, 2, 7, 2, 3},
    {3, 6, 4, 5, 9, 1, 1}
   } );
q = j[m];
k = vec[p (d\[TensorProduct]q)];
A = ArrayFlatten[{IdentityMatrix[m]\[TensorProduct]j[n]}];
b = j[m]\[TensorProduct]{1, 0};
mat[LinearProgramming[k, A, b]] // MatrixForm

Output: $$\left( \begin{array}{ccccccc} 0 & 1 & 0 & 0 & 0 & 1 & 1 \\ 1 & 0 & 1 & 1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \end{array} \right)$$

$\endgroup$
  • $\begingroup$ I've upvoted your answer until I have the chance to actually try it, but it seems what I was looking for $\endgroup$ – StepTNT Aug 3 '14 at 9:07
  • $\begingroup$ I just edited the question with my work based on your solution because the result is not what I expected and I think I did something wrong $\endgroup$ – StepTNT Aug 3 '14 at 11:22
  • $\begingroup$ The edit was perfect, thanks for the reply and sorry for messing up things :) $\endgroup$ – StepTNT Aug 3 '14 at 14:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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