For the vectors $X=(x_1,\cdots, x_n),~ Y=(y_1,\cdots, y_n)$ and $\alpha=(\alpha_1,\cdots,\alpha_n),~ \beta=(\beta_1,\cdots, \beta_n)\in\mathbb R^n_+$ s.t. $\sum_{k=1}^n\alpha_k~~=~~\sum_{k=1}^n\beta_k~~=~~1$, let $\mathcal P$ be the collection of matrices $p=(p_{i,j})_{1\le i\le n, 1\le j\le n}\in\mathbb R^{n\times n}_+$ satisfying

\begin{eqnarray} &&\sum_{j=1}^np_{i,j}~=~\alpha_i,~ \forall 1\le i\le n;~~~~~\sum_{i=1}^np_{i,j}~=~\beta_j,~ \forall 1\le j\le n; \\ &&\sum_{j=1}^np_{i,j}y_j~=~\alpha_ix_i,~ \forall 1\le i\le n. \end{eqnarray}

Given a matrix $c=(c_{i,j})_{1\le i\le n, 1\le j\le n}$, my question is how to **solve numerically** the following optimization problem:

\begin{eqnarray} \sup_{p=(p_{i,j})_{1\le i\le n, 1\le j\le n}\in\mathcal P}~\sum_{i=1}^n\sum_{j=1}^n p_{i,j}c_{i,j}. \end{eqnarray}

Here assume $\mathcal P\neq \emptyset$, does someone know some related algorithm? Thank you very much!