Given $n^2$ constants $a_{11},a_{12},\ldots,a_{1n},a_{21},\ldots,a_{nn}$ and $n^2$ non-negative variables $x_{11},x_{12},\ldots,x_{1n},x_{21},\ldots,x_{nn}$. Find the minimum value of $$\sum_{i=1}^n \left(\sum_{j=1}^n x_{ij}\right)^2 + \sum_{j=1}^n \left(\sum_{i=1}^n x_{ij}\right)^2 + \sum_{i=1}^n \sum_{j=1}^n a_{ij} x_{ij}.$$
Even for $n=2$ the problem is not trivial at all for me:
Find the minimum value of $$ f(x,y,z,t) = (x+y)^2 + (z+t)^2 + (x+z)^2 + (y+t)^2 + ax + by +cz +dt$$ where $a,b,c,d$ are constants and $x,y,z,t\ge 0$ are variables.
Thank you in advance for your discussion.
The presence of non-negativity constraints on $x$ suggests that there is no way to guarantee that there will be an analytic solution to this problem, at least not without further information about the values $a_{ij}$. Or did the person who posed this question to you believe that they can make such a guarantee?
What you have here is $$\begin{array}{ll} \text{minimize} & \|Xe\|_2^2 + \|X^Te\|_2^2 + \langle A,X \rangle \\ \text{subject to} & x_{ij} \geq 0, ~ i,j\in\{1,2,\dots,n\} \end{array}$$ where $e$ is the vector of all ones. This is actually a convex simple quadratic program, but with some Kronecker structure that makes it a bit more difficult to express. If you use a convex modeling framework like YALMIP or CVX (disclaimer: the latter is mine) you can express it rather simply. For example, in CVX, the model is
EDIT: For grins I decided see how far one could get with an analytic solution. It's not easy working with matrix variables! This is the Lagrangian: $$L(X,Z) = \|Xe\|_2^2+\|X^Te\|_2^2+\langle A,X \rangle - \langle Z,X\rangle$$ where $Z\in\mathbb{R}^{n\times n}$ contains the Lagrange multipliers for the nonnegativity constraints; it is elementwise nonnegative. The gradient of the Lagrangian is $$2ee^T X + 2 X ee^T + A - Z = 0$$ To see how I arrived at this, note that $$\|Xe\|_2^2=e^TX^TXe=\mathop{\textrm{Tr}}(X^TXee^T)=\langle X, Xee^T \rangle \quad\Longrightarrow\quad \nabla_X \|Xe\|_2^2 = 2Xee^T$$ $$\|X^Te\|_2^2=e^TXX^Te=\mathop{\textrm{Tr}}(X^Tee^TX) = \langle X, ee^TX \rangle \quad\Longrightarrow\quad \nabla_X \|X^Te\|_2^2 = 2ee^TX$$ I believe this is what we get if we look at an individual term $(i,j)$: $$2\sum_p X_{pj}+2\sum_{q} X_{iq} + A_{ij} - Z_{ij} = 0, \quad i,j\in\{1,2,\dots,n\}$$ It's interesting that the optimality conditions depend only on the $2n$ row and column sums. I am sure that can be exploited somehow. Also, note that if I treat $Z_{ij}$ as a slack variable, I get $$\sum_p X_{pj}+\sum_{q} X_{iq} \leq \tfrac{1}{2} A_{ij}, \quad i,j\in\{1,2,\dots,n\}$$ I'm going to be honest, I'm not yet sure what I'd do from here---if anything. Since existing optimization frameworks handle this problem simply, I'm not sure I'd bother to try for an analytic solution. Still, even if an analytic solution does not exist, exploiting the structure of this problem could help you arrive at much faster approaches to solving the problem numerically.