Is the following optimization problem feasible (modified optimal transport)? The primary decision variable appears in the objective function as $\mathbf{x}$, whereas the second decision variable $\mathbf{y}\in \mathbb{R}^K$ is introduced in the constraint vector $\mathbf{b}$: This variable will be optimized as an optimal weight vector for multiplication against a multivariate data matrix $\mathbf{M}\in \mathbb{R}^{N\times K}$ as shown in the last bullet point. Product $\mathbf{M} \cdot \mathbf{y}$ therefore is an optimally weighted version of the original matrix, or optimal source distribution, but of size $N\times 1$.
$$ \begin{array}{rrcl} \arg \min_{\mathbf{x}, \, \mathbf{y}} \ & \mathbf{c}^\top \mathbf{x} & & \ \\ \mathrm{s.t.} \ & \mathbf{A} \mathbf{x} & = & \ \mathbf{b} \\ \ & \mathbf{x} & \geq &\ \mathbf{0} \end{array}$$
- $\mathbf{c}$ is a vectorization of the distance matrix $\boldsymbol C$,
- $\mathbf{x}$ is a vectorization of the transport matrix $\boldsymbol\Gamma$,
- while the source and target distributions are $\mathbf{b} = \begin{bmatrix} \mathbf{M} \cdot \mathbf{y}\\ \beta\end{bmatrix}$
Any comments on this optimization set-up with two decision variables would be helpful.
(Next issue is how to make standard solvers like scipy.optimize.linprog admit the existence of a second decision variable, because the documentation of most solvers only addresses one decision variable only. but I ask this part at stack instead)
Your objective function can be written as $c^T x + 0^T y$. Your constraints can be written in block-matrix form with constant RHS as follows, where the subscripts indicate the sizes: $$ \begin{pmatrix} A'_{\ell \times n} & -M_{\ell \times K}\\ A''_{\ell \times n} & 0_{\ell \times K}\\ \end{pmatrix} \begin{pmatrix} x_{n \times 1}\\ y_{K \times 1}\\ \end{pmatrix} = \begin{pmatrix} 0_{\ell \times 1}\\ \beta_{\ell \times 1}\\ \end{pmatrix} $$
If your solver requires a single vector of decision variables, then rename $y$, as described by @MichalAdamaszek. Explicitly, $$(x_1,\dots,x_n,y_1,\dots,y_K)$$ becomes $$(x_1,\dots,x_n,x_{n+1},\dots,x_{n+K}).$$