# Minimizing a function with some restrictions

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:

$$\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})]\}$$

with those given constants

$$ds_x > 0, qp_x > 0, pp_{xy} > 0$$

and

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

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

$$2^{3*7}=2097152$$

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)

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

$$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}$$

assuming that

$$qp_i = 1, \forall i$$

we have

$$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}$$

which, translated to a vecotr, becomes

$$k=[35,10,33,14,40,66,35,15,44,21,10,55,49,35,99,7,10,11,7,15,11]$$

(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.

• 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? Aug 3, 2014 at 0:39
• 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. Aug 3, 2014 at 2:01
• 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. Aug 3, 2014 at 13:48
• 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? Aug 3, 2014 at 14:06
• 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. Aug 3, 2014 at 14:24

## 1 Answer

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:

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)$$

• I've upvoted your answer until I have the chance to actually try it, but it seems what I was looking for Aug 3, 2014 at 9:07
• 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 Aug 3, 2014 at 11:22
• The edit was perfect, thanks for the reply and sorry for messing up things :) Aug 3, 2014 at 14:27