It's clear that
$$
\sup _{\Phi_{c} \cap C_{b}} J(\varphi, \psi) \leq \sup _{\Phi_{c} \cap L^{1}} J(\varphi, \psi) \leq \inf _{\Pi(\mu, \nu)} I[\pi].
$$
So it remains to prove
$$
\inf _{\Pi(\mu, \nu)} I[\pi] \le \sup _{\Phi_{c} \cap C_{b}} J(\varphi, \psi) .
$$
To simplify notations, let $\varphi \oplus \psi: (x,y) \mapsto \varphi(x)+\psi(y)$. We have
$$
\inf _{\pi \in \Pi(\mu, \nu)} I[\pi]=\inf _{\pi \in M_{+}(X \times Y)}\left(I[\pi]+\left\{\begin{array}{l}
0 \text { if } \pi \in \Pi(\mu, \nu) \\
+\infty \text { else }
\end{array}\right)\right.
$$
with $M_+(X \times Y)$ the space of non-negative Borel measures on $X\times Y$. Also,
$$
\left\{\begin{array}{l}
0 \text { if } \pi \in \Pi(\mu, \nu) \\
+\infty \text { else }
\end{array}\right\}=\sup \left \{ \int \varphi d \mu+\int \psi d \nu-\int \varphi \oplus \psi d \pi \,\middle\vert\, (\varphi, \psi) \in C_b(X) \times C_b(Y) \right \},
$$
It follows that
$$
\begin{aligned}
&\inf _{\pi \in \Pi(\mu, \nu)} I[\pi]\\
=&\inf _{\pi \in M_{+}(X \times Y)} \sup \left \{ \int_{X \times Y} c d \pi +\int_{X} \varphi d \mu+\int_{Y} \psi d \nu - \int_{X \times Y} \varphi \oplus \psi d \pi \,\middle\vert\, (\varphi, \psi) \in C_b(X) \times C_b(Y) \right\}.
\end{aligned}
$$
- Let us first assume that $X, Y$ are compact.
Let $E:=C_{b}(X \times Y)$ be the set of all bounded continuous functions on $X \times Y$, equipped with its usual supremum norm $\|\cdot\|_{\infty}$.
Because $X,Y$ are compact, so is $X\times Y$. So $E$ coincides with the space $C_{0}(X \times Y)$ of all continuous functions vanishing at infinity on $X \times Y$.
By Riesz' theorem, the topological dual $E^*$ of $E$ can be identified with the space of regular Radon measures, $M(X \times Y)$, normed by total variation.
Moreover, a nonnegative linear form on $E$ corresponds with a regular nonnegative Borel measure.
Then we introduce
$$
\begin{gathered}
\Theta: u \in E \longmapsto\left\{\begin{array}{l}
0 &\text {if } u \geq-c, \\
+\infty &\text {else}.
\end{array}\right.\\
\text{}\\
\Xi: u \in E \longmapsto\left\{\begin{array}{l}
\int_{X} \varphi d \mu+\int_{Y} \psi d \nu &\begin{align*}
&\text {if } u = \varphi \oplus \psi \text{ for}\\
& \text{some } (\varphi, \psi) \in C_b(X) \times C_b(Y)
\end{align*},\\
+\infty &\text {else. }
\end{array}\right.
\end{gathered}
$$
It's easy to verify that $\Theta, \Xi$ satisfy the condition of Fenchel-Rockafellar duality, so
$$
\inf _{u\in E}[\Theta(u)+\Xi(u)] = \max _{\pi \in E^{*}}\left[-\Theta^{*}\left(-\pi\right)-\Xi^{*}\left(\pi\right)\right].
$$
It's clear that
$$
\begin{align*}
\inf _{u\in E}[\Theta(u)+\Xi(u)] &= \inf \left\{\int_{X} \varphi d \mu+\int_{Y} \psi d \nu \,\middle\vert\, (\varphi, \psi) \in C_b(X) \times C_b(Y) \text{ s.t. } \varphi \oplus \psi \geq -c \right\} \\
&=-\sup \left\{J(\varphi, \psi) \mid (\varphi, \psi) \in \Phi_{c}\cap C_{b}\right\}.
\end{align*}
$$
Next, we compute the Legendre-Fenchel transforms of $\Theta, \Xi$. First, for any $\pi \in E^*$,
$$
\begin{aligned}
\Theta^{*}(-\pi) =\sup _{u \in E}\left\{-\int u d \pi \,\middle\vert\, u \geq-c\right\}
= \sup _{u \in E}\left\{\int u d \pi \,\middle\vert\, u \leq c\right\} .
\end{aligned}
$$
- If $\pi$ is not nonnegative, then there exists a positive function $v \in E$ such that $\int v d \pi<0$. Then, the choice $u=\lambda v$, with $\lambda \rightarrow-\infty$, shows that the supremum is $+\infty$.
- On the other hand, if $\pi$ is nonnegative, then the supremum is clearly $\int c d \pi$. This is because lower semi-continuous and bounded from below function is a limit of an increasing sequence of Lipschitz continuous functions.
Thus
$$
\Theta^{*}(-\pi) = \begin{cases}
\int c d \pi &\text {if } \pi \in M_{+}(X \times Y) \\
+\infty &\text {else}.
\end{cases}
$$
We also have
$$
\begin{align*}
\Xi^{*}(\pi) &= \sup_{u\in E} \left \{ \int ud\pi - \int \varphi d \mu- \int \psi d\nu \,\middle\vert\, u = \varphi \oplus \psi \text{ for some } (\varphi, \psi) \in C_b(X) \times C_b(Y) \right \} \\
&= \begin{cases}
0 &\text {if } \quad \forall(\varphi, \psi) \in C_b(X) \times C_b(Y) : \int \varphi \oplus \psi d \pi= \int \varphi d \mu+\int \psi d \nu \\
+\infty & \text {else}
\end{cases} \\
&= \begin{cases}
0 &\text {if } \quad \pi \in \Pi(\mu, \nu) \\
+\infty & \text {else}.
\end{cases}
\end{align*}
$$
It follows that
$$
\max _{\pi \in E^{*}}\left[-\Theta^{*}\left(-\pi\right)-\Xi^{*}\left(\pi\right)\right] = \max _{\pi \in \Pi(\mu, \nu) \cap M_{+}(X \times Y)} -\int c d \pi = - \min _{\pi \in \Pi(\mu, \nu)} \int c d \pi.
$$
Hence
$$
-\sup \left\{J(\varphi, \psi) \mid (\varphi, \psi) \in \Phi_{c} \cap C_{b}\right\} = - \min _{\pi \in \Pi(\mu, \nu)} \int c d \pi.
$$
- $c$ is bounded.
Fix
$$
\pi_* \in \operatorname{argmin}_{\pi \in \Pi(\mu, \nu)}I [\pi].
$$
Then $\pi_*$ is tight, so there are compact sets $X_n \subset X$ and $Y_n \subset Y$ such that $\mu(X_n^c), \nu(Y_n^c) \le \delta$. Let $Z_n := X_n \times Y_n$. Then $\pi (Z_n^c) \le 2\delta$. We define $\pi_n \in \mathcal P(Z_n)$ by
$$
\pi_n := \frac{1_{Z_n}\pi_*}{\pi_* (Z_n)}.
$$
Let $\mu_n := P_\sharp^{X_n} \pi_{n} \in \mathcal P(X_n)$ and $\nu_n := P_\sharp^{Y_n} \pi_{n} \in \mathcal P(Y_n)$ be marginals of $\pi_n$. Let
$$
I_n[\pi] := \int_{Z_n} c d \pi \quad \forall \pi \in \Pi(\mu_n, \nu_n).
$$
Fix
$$
\tilde \pi_{n} \in \operatorname{argmin}_{\pi \in \Pi(\mu_n, \nu_n)} I_n [\pi].
$$
Define $\pi_{*n}$ by
$$
\pi_{*n} := \pi_* (Z_n) \tilde \pi_{n} + 1_{Z_n^c} \pi_*.
$$
It is easy to verify that $\pi_{*n} \in \Pi(\mu, \nu)$. We have
$$
\begin{align}
I[\pi_*] &\le I[\pi_{*n}] \\
&= \int_{X\times Y} c d \pi_{*n} \\
&= \pi_*(Z_n)\int_{Z_n} c d \tilde \pi_n + \int_{Z_n^c} c d \pi_* \\
&\le \pi_*(Z_n) I_n [\tilde \pi_n] + 2 \delta \|c\|_\infty \\
& \le I_n [\tilde \pi_n] + 2\delta \|c\|_\infty.
\end{align}
$$
Now we introduce
$$
J_n: L^{1}(\mu_{n}) \times L^{1}(\nu_{n}) \to \mathbb R, (\varphi, \psi) \mapsto \int_{X_{n}} \varphi d \mu_{n}+\int_{Y_{n}} \psi d \nu_{n}.
$$
Let $\Phi_{c,n}$ be the set of all $(\varphi, \psi) \in L^{1}(\mu_n) \times L^{1}(\nu_n)$ satisfying $\varphi(x)+\psi(y) \leq c(x, y)$ for $\mu_n$-almost all $x \in X_n$ and $\nu_n$-almost all $y \in Y_n$. By step 1, we know that
$$
\inf_{\pi \in \Pi(\mu_n, \nu_n)} I_n [\pi] = \sup_{(\varphi, \psi) \in \Phi_{c,n}} J_n (\varphi, \psi).
$$
In particular, there is $(\varphi_n, \psi_n) \in \Phi_{c,n}$ such that
$$
J_n (\varphi_n, \psi_n) \ge \sup_{(\varphi, \psi) \in \Phi_{c,n}} J_n (\varphi, \psi) - \delta.
$$
To ensure $\varphi_n (x)+\psi_n (y) \leq c(x, y)$ for all $x \in X_n$ and all $y \in Y_n$. We redefine $\varphi_n, \psi_n$ such that $\varphi_n, \psi_n$ take value $-\infty$ on their null sets respectively. WLOG, we assume $\delta <1$. Then $J_n (\varphi_n, \psi_n) \ge -1$. Then there is $(x_n, y_n) \in Z_n$ such that
$$
\varphi_n(x_n)+\psi_n(y_n) \ge -1.
$$
If we replace $(\varphi_n, \psi_n)$ by $(\varphi_n +s, \psi_n-s)$ for some real number $s$, we do not change the value of $J_n (\varphi_n, \psi_n)$, so the resulting couple is still admissible. By a proper choice of $s$, we can ensure
$$
\varphi_n (x_n) \ge -\frac{1}{2} \quad \text{and} \quad\psi_n (y_n) \ge -\frac{1}{2}.
$$
This implies for all $(x,y) \in Z_n$,
$$
\varphi_n(x) \le c(x, y_n) - \psi_n(y_n) \le c(x, y_n) + 1/2,\\
\psi_n(y) \le c(x_n, y) - \varphi_n(x_n) \le c(x_n, y) + 1/2.
$$
Define $\psi_n^c$ by
$$
\psi_n^c (x) := \inf_{y\in Y_n} [c(x, y) - \psi_n(y)] \quad \forall x \in X_n.
$$
Moreover,
$$
\inf_{y\in Y_n} [c(x,y)-c(x_n, y)] - 1/2 \le \psi_n^c (x) \le c (x, y_n) - \psi_n (y_n) \le c(x, y_n) + 1/2.
$$
It follows from $c$ is bounded that $\psi_n^c$ is bounded. Hence $\psi_n^c \in L^1(\mu)$. It follows from $\varphi_n(x)+\psi_n(y) \leq c(x, y)$ that $\varphi_n \le \psi_n^c$ on $X_n$. This implies $J_n (\psi_n^c , \psi_n) \ge J_n (\varphi_n, \psi_n)$. Similarly, we define $\psi_n^{cc}$ by
$$
\psi_n^{cc} (y) := \inf_{x\in X_n} [c(x, y) - \psi_n^c (x)] \quad \forall y \in Y_n.
$$
For all $y\in Y_n$, we have
$$
\inf_{x\in X_n} [c(x,y)-c(x, y_n)] - 1/2 \le \psi_n^{cc} (y) \le c(x_n, y) - \psi_n^c (x_n) \le c(x_n, y) - \varphi_n (x_n) \le c(x_n,y)+1/2.
$$
So $\psi_n^{cc}$ is also bounded and thus $\psi_n^{cc} \in L^1(\nu)$. Moreover,
$$
\psi_n^c (x) + \psi_n^{cc} (y) = \psi_n^c (x) + \inf_{z \in X} [c(z, y) - \psi_n^c(z)] \le \psi_n^c(x) + [c(x, y) - \psi_n^c(x)] = c(x,y).
$$
Then $(\psi_n^c, \psi_n^{cc}) \in \Phi_c$. Also
$$
\psi_n^{cc} (y) = \inf_{x\in X_n} [c(x, y) - \psi_n^c(x)] \ge \inf_{x\in X_n} [c(x, y) - [c(c,y) - \psi_n (y)]] = \psi_n (y) \quad \forall y \in Y_n.
$$
This implies
$$
J_n (\psi_n^c , \psi_n^{cc}) \ge J_n(\psi_n^c, \psi_n) \ge J_n (\varphi_n, \psi_n).
$$
In particular,
$$
\psi_n^c (x) \ge - \|c\|_\infty -1/2 \quad \text{and} \quad \psi_n^{cc} (y) \ge - \|c\|_\infty -1/2.
$$
Indeed,
$$
\begin{align*}
J(\psi_n^c, \psi_n^{cc})
&= \int_{X} \psi_n^c d \mu+\int_{Y} \psi_n^{cc} d \nu \\
&= \int_{X \times Y}[\psi_n^c (x) + \psi_n^{cc} (y)] d \pi_{*n}(x, y) \\
&= \pi_{*}[Z_n] \int_{Z_n}[ \psi_n^c (x) + \psi_n^{cc} (y)] d \tilde \pi_{n}(x, y) +\int_{Z_n^{c}}[\psi_n^c(x)+\psi_n^{cc}(y)] d \pi_{*}(x, y) \\
&\ge (1-2 \delta) \left (\int_{X_n} \psi_n^c d \mu_n+\int_{Y_n} \psi_n^{cc} d \nu_n \right )-(2\|c\|_{\infty}+1) \pi_{*}[Z_n^c] \\
&\geq(1-2 \delta) J_n( \psi_n^c , \psi_n^{cc} )-2(2\|c\|_{\infty}+1) \delta \\
&\geq(1-2 \delta) J_n({\varphi}_n, {\psi}_n)-2(2\|c\|_{\infty}+1) \delta \\
&\geq(1-2 \delta) \left ( \inf_{\pi \in \Pi(\mu_n, \nu_n)} I_n [\pi] - \delta \right)-2(2\|c\|_{\infty}+1) \delta \\
&= (1-2 \delta)( I_n [\tilde \pi_n] -\delta)-2(2\|c\|_{\infty}+1) \delta \\
&\ge (1-2 \delta)(I[\pi_{*}] -(2\|c\|_{\infty}+1) \delta)-2(2\|c\|_{\infty}+1) \delta.
\end{align*}
$$
The result then follows by taking the limit $\delta \to 0^+$.
- Let us assume that $X, Y$ are not compact, and that $c$ is lower semi-continuous on $X \times Y$.
Fix $(\varphi, \psi) \in \Phi_c$.
- There are a $\mu$-null set $N_x$ and $\nu$-null net $N_y$ such that $\varphi(x)+\psi(y) \le c(x,y)$ for all $(x,y) \in N_x^c \times N_y^c$.
- We re-define $\varphi, \psi$ by letting them take value $-\infty$ on $N_x, N_y$ respectively. In this way, $\varphi(x)+\psi(y) \le c(x,y)$ for all $(x,y) \in X \times Y$.
We define $\varphi^c$ by
$$
\varphi^c (y) := \inf_{x\in X} [c(x,y) - \varphi(x)].
$$
- The infimum of a collection of extended real-valued measurable functions is again a measurable function, so $\varphi^c$ is measurable.
- There is $x_0 \in N^c_x$, so $\varphi^c (y) \le c(x_0,y)-\varphi(x_0) <+\infty$ for all $y \in Y$.
- Given $y\in Y$, $\varphi(x)+\psi(y) \le c(x,y)$ for all $x \in X$, so $\varphi^c (y) \ge \psi(y)$ for all $y\in Y$.
We define $\varphi^{cc}$ by
$$
\varphi^{cc} (x) := \inf_{y\in Y} [c(x,y) - \varphi^c (y)].
$$
For all $x\in X$,
$$
\begin{align}
\varphi^{cc} (x) &= \inf_{y\in Y} \left [c(x,y) - \inf_{z\in X} [c(z,y) - \varphi(z)] \right] \\
&= \inf_{y\in Y} \left [c(x,y) + \sup_{z\in X} [-c(z,y) + \varphi(z)] \right] \\
&= \inf_{y\in Y} \sup_{z\in X} [c(x,y) -c(z,y) + \varphi(z)] \\
&\ge \inf_{y\in Y} \varphi(x) =\varphi(x) \quad \text{by picking} \quad z=x.
\end{align}
$$
There is $y_0 \in N_y^c$. Then
$$
\varphi^{cc} (x) \le c(x, y_0) - \varphi^c (y_0) \le c(x, y_0) - \psi(y_0) <+\infty \quad \forall x\in X.
$$
For all $(x,y) \in X \times Y$,
$$
\varphi^{cc} (x) + \varphi^c (y) = \inf_{z\in Y} [c(x,z) - \varphi^c(z)] + \varphi^c (y) \le [c(x,y) - \varphi^c (y)] + \varphi^c (y)= c(x,y).
$$
- There are a $\mu$-null set $M_x$ and $\nu$-null set $M_y$ such that $c (x, y) \le c_X(x)+c_Y(y)$ for all $(x, y) \in M_x^c \times M_y^c$.
- We re-define $c_X, c_Y$ such that $c_X (x) = c_Y(y) := +\infty$ for all $(x, y) \in M_x \times M_y$. In this way, $c(x,y) \le c_X(x)+c_Y(y)$ for all $(x,y) \in X \times Y$.
Let
$$
a := \inf_{y\in Y} [c_Y(y) - \varphi^c (y)].
$$
First,
$$
a \le c_Y(y_0) - \varphi^c (y_0) \le c_Y(y_0) - \psi (y_0) < +\infty.
$$
There is $x_1 \in (N_x \cup M_x)^c$, so
$$
\begin{align}
c_Y(y) - \varphi^c (y) &= c_Y(y) - \inf_{x\in X} [c(x,y) - \varphi(x)] \\
&= \sup_{x\in X} [c_Y(y)-c(x,y) + \varphi(x)] \\
&\ge \sup_{x\in X} [\varphi(x)-c_X(x)] \\
&\ge \varphi(x_1)-c_X(x_1) \\
&> -\infty.
\end{align}
$$
It follows that $a \in \mathbb R$. Clearly, $a \le c_Y(y) - \varphi^c (y)$ for all $y\in Y$, so $\varphi^c+a \le c_Y$. We have
$$
\begin{align}
(\varphi^{cc} (x)-a) - c_X(x) &= \inf_{y\in Y} [c(x,y) - \varphi^c (y)]-a - c_X(x) \\
&= \inf_{y\in Y} [(c(x,y)-c_X(x)) - \varphi^c (y)] -a \\
&\le \inf_{y\in Y} [c_Y(y) - \varphi^c (y)] -a \\&=0.
\end{align}
$$
Clearly, $(\varphi', \psi') := (\varphi^{cc}-a, \varphi^c+a)$ satisfies the requirement.
Best Answer
Let $c_X = c_Y := \frac{1}{2} |\cdot|^2$. We need the following Lemma
Let $\pi^\dagger$ minimize $\mathbb K$ over $\Pi(\mu, \nu)$. By Lemma, there is $(\varphi, \varphi^c)$ with $\varphi^{cc} = \varphi$ maximizing $\mathbb J$ over $\Phi_c$. Because $c(x, y) = c_X(x) + c_Y(y) - \langle x, y\rangle$, we have $(\tilde \varphi, \tilde \psi) := (c_X - \varphi, c_Y - \varphi^c)$ minimizes $\mathbb J$ over $\tilde \Phi$. It's can be shown that $\tilde \psi$ is the convex conjugate of $\tilde \varphi$, i.e., $\tilde \psi = \tilde \varphi^*$, and that $\tilde \varphi$ is the convex conjugate of $\tilde \varphi^*$, i.e., $\tilde \varphi = \tilde \varphi^{**}$. By this result, $\tilde \varphi$ is convex l.s.c. By Lemma, $$ \int_X \tilde \varphi (x) d \mu (x) + \int_Y \tilde \psi (y) d \nu (y)= \int_{X \times Y} \langle x, y\rangle d\pi^\dagger(x, y), $$ or equivalently $$ \int_{X \times Y} [\tilde \varphi (x) + \tilde \varphi^* (y) - \langle x, y\rangle] d\pi^\dagger(x, y) = 0. $$
We have $\tilde \varphi (x) + \tilde \varphi^* (y) - \langle x, y\rangle \ge 0$ for all $(x,y) \in X \times Y$, so $\tilde \varphi (x) + \tilde \varphi^* (y) - \langle x, y\rangle = 0$ for $\pi^\dagger$-a.e. $(x, y) \in X \times Y$. By this result, $y \in \partial \tilde \varphi (x)$ for $\pi^\dagger$-a.e. $(x, y) \in X \times Y$. It follows from $\tilde \varphi = \tilde \varphi^{**}$ that $\tilde \varphi$ is convex l.s.c.
Let $\pi^\dagger \in \Pi(\mu, \nu)$ and $\tilde \varphi \in L_1 (\mu)$ be such that $y \in \partial \tilde \varphi (x)$ for $\pi^\dagger$-a.e. $(x, y) \in X \times Y$. Then $\tilde \varphi (x) + \tilde \varphi^* (y) - \langle x, y\rangle = 0$ for $\pi^\dagger$-a.e. $(x, y) \in X \times Y$, so $$ \int_{X \times Y} [\tilde \varphi (x) + \tilde \varphi^* (y) - \langle x, y\rangle] d\pi^\dagger(x, y) = 0. $$
By this result, $\langle \cdot, \cdot \rangle$ is $\pi^\dagger$-integrable, so the map $(x, y) \mapsto \tilde \varphi (x) + \tilde \varphi^* (y)$ is also $\pi^\dagger$-integrable. Because $\tilde \varphi$ is $\mu$-integrable, so $\tilde \varphi^* \in L_1 (\nu)$. It follows that $$ \int_X\tilde \varphi d \mu + \int_Y \tilde \varphi^* d \nu = \int_{X \times Y}\langle x, y\rangle d\pi^\dagger(x, y). $$
On the other hand, $$ \int_X \varphi' d \mu + \int_Y \psi' d \nu \ge \int_{X \times Y} \langle x, y\rangle d\pi^\dagger(x, y) \quad \forall (\varphi', \psi') \in \tilde \Phi. $$
It follows that $(\tilde \varphi, \tilde \varphi^*)$ minimizes $\mathbb J$ over $\tilde \Phi$. Because $c(x, y) = c_X(x) + c_Y(y) - \langle x, y\rangle$, we have $(\varphi, \psi) := (c_X - \tilde\varphi, c_Y - \tilde \varphi^*)$ minimizes $\mathbb J$ over $\Phi_c$. This completes the proof.