I've the following linear programming problem:

This is the LP representation of the uncapacitated facility location problem.
This is the dual representation of this problem:

My question is how to derive the dual representation from the primal in this problem.
I tried to follow the step by step tutorial by S´ebastien Lahaie found here but it won't work because of the summation in the constraint equations.
I'm looking for a pure mechanical way (with no intuition) to do that.
My preference is to build the Lagrangian and go from there. In this case, the Lagrangian is $$L(x,y,u,v,w,s)=\sum_{i\in F} f_iy_i + \sum_{i\in F,j\in D} c_{ij} x_{ij} - \sum_{j\in D} v_j \left( \sum_{i\in F} x_{ij} - 1 \right) - \sum_{i\in F,j\in D} w_{ij} ( y_i - x_{ij} ) - \sum_{i\in F,j\in D} s_{ij} x_{ij} - \sum_{i\in F} u_i y_i$$ Every additional summation corresponds to a constraint. The Lagrange multipliers for the equality constraints are $v$, and for the inequality constraints $w$, $s$, and $u$, respectively; these three must all be (elementwise) nonnegative.
It is important to get the signs right for the inequalities! You should be substracting the product of two nonnegative terms for each inequality: the inequality expression times its corresponding multiplier. Also, this approach assumes that the primal problem is a minimization. If it's not, negate the objective to make it so---but don't forget to negate the dual objective at the end.
Now you rearrange and separate the terms involving each primal variable, and the terms that depend on neither: $$L(x,y,u,v,w_1,w_2)=\sum_{i\in F} y_i\left(f_i-\sum_{j\in D}w_{ij}-u_i\right) + \sum_{i\in F,j\in D} x_{ij} \left( c_{ij} - v_j + w_{ij} - s_{ij} \right) + \sum_{j\in D} v_j$$ Each expression that multiplies a primal variable becomes an $=0$ equality constraint in the dual; each expression that does not contain a primal variable enters the objective. This is due to the fact that that if you take the partial derivative of a Lagrangian with respect to a primal variable, you'll get an affine expression of the dual variables; if that expression is not zero, the Lagrangian would be unbounded below. \begin{array}{lll} \text{maximize} & \sum_{j\in D} v_j \\ \text{subject to} & f_i - \sum_{j\in D} w_{1,ij} - u_i = 0 & i\in F \\ & c_{ij} - v_j + w_{ij} - s_{ij} = 0 & i\in F, ~j\in D \\ & w_{ij},s_{ij} \geq 0 & i\in F,~j\in D \\ & u_i \geq 0 & i \in F \end{array} Eliminating $u_i$ gives you $$f_i - \sum_{j\in D} w_{ij} - u_i = 0 \quad\Longrightarrow\quad \sum_{j\in D} w_{ij} \leq f_i$$ Eliminating $s_{ij}$ gives $$c_{ij} - v_j + w_{ij} - s_{ij} = 0 \quad\Longrightarrow\quad v_j - w_{ij} \leq c_{ij}$$ And that should give you the dual you offered above.