Minimizing a function with some restrictions

Question

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.

Answer 1

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