You can use couplings from regularized optimal transport: For $n=1$ (i.e. a coupling of just two marginals) there are two popular ways.

Entropic regularization: Use $\pi_k$ defined by $$\newcommand{\RR}{\mathbb{R}} \min_\pi \int_{\RR^2} cd\pi + \gamma_k\int_{\RR^2}\pi(\log(\pi)-1)d\lambda $$ where $\lambda$ denotes the Lebesgue measure, $c$ is some cost function (e.g. continuous and non-negative), and $\gamma_k$ is a sequence of regularization parameters converging to zero and the minimum is taken over all measure $\pi$ which have a density with respect to the Lebesgue-measure (and the density has finite entropy) and which have the marginals $\mu_1$ and $\mu_2$. If $\gamma_k$ tends to zero it should hold that $\pi_k$ converges weakly (in the sense of measures) to $\pi$.

Quadratic regularization: Use $\pi_k$ defined by $$\newcommand{\RR}{\mathbb{R}} \min_\pi \int_{\RR^2} c d\pi + \gamma_k\int_{\RR^2}\pi^2 d\lambda $$ where $\lambda$ denotes the Lebesgue measure, $c$ is again some cost function (e.g. square integrable and non-negative) and $\gamma_k$ is a sequence of regularization parameters converging to zero and the minimum is taken over all measures $\pi$ which have a square integrable density with respect to the Lebesgue-measure and which have the marginals $\mu_1$ and $\mu_2$. I guess that weak convergence is also true here.

For entropic regularization I recommend Computational Optimal Transport (free version here) and for quadratic regularization I only have my own paper Quadratically Regularized Optimal Transport (also on the arXiv).

If you work on compact domains instead of $\RR$ (or if you assume that all marginals have compact support), you can use couplings from regularized optimal transport: For $n=1$ (i.e. a coupling of just two marginals) there are two popular ways.

Entropic regularization: Use $\pi_k$ defined by $$\newcommand{\RR}{\mathbb{R}} \min_\pi \int_{\RR^2} cd\pi + \gamma_k\int_{\RR^2}\pi(\log(\pi)-1)d\lambda $$ where $\lambda$ denotes the Lebesgue measure, $c$ is some cost function (e.g. continuous and non-negative), and $\gamma_k$ is a sequence of regularization parameters converging to zero and the minimum is taken over all measures $\pi$ which have a density with respect to the Lebesgue-measure (and the density has finite entropy), which have their support in the product of the supports of the marginals, and which have the marginals $\mu_1$ and $\mu_2$. If $\gamma_k$ tends to zero it should hold that $\pi_k$ converges weakly (in the sense of measures) to $\pi$.

Quadratic regularization: Use $\pi_k$ defined by $$\newcommand{\RR}{\mathbb{R}} \min_\pi \int_{\RR^2} c d\pi + \gamma_k\int_{\RR^2}\pi^2 d\lambda $$ where $\lambda$ denotes the Lebesgue measure, $c$ is again some cost function (e.g. square integrable and non-negative) and $\gamma_k$ is a sequence of regularization parameters converging to zero and the minimum is taken over all measures $\pi$ which have a square integrable density with respect to the Lebesgue-measure, having their support in the product of the supports of the marginals and which have the marginals $\mu_1$ and $\mu_2$. I guess that weak convergence is also true here.

For entropic regularization I recommend Computational Optimal Transport (free version here) and for quadratic regularization I only have my own paper Quadratically Regularized Optimal Transport (also on the arXiv).