  We introduce the \emph{transport energy} functional $\mathcal E$ (a variant
of the Bouchitt\'e-Buttazzo-Seppecher shape optimization functional) and we
prove that its unique minimizer is the optimal transport density $\mu^*$, i.e.,
the solution of Monge-Kantorovich equations. We study the gradient flow of
$\mathcal E$ showing that $\mu^*$ is the unique global attractor of the flow.
  We introduce a two parameter family $\{\mathcal
E_{\lambda,\delta}\}_{\lambda,\delta>0}$ of strictly convex functionals
approximating $\mathcal E$ and we prove the convergence of the minimizers
$\mu_{\lambda,\delta}^*$ of $\mathcal E_{\lambda,\delta}$ to $\mu^*$ as we let
$\delta\to 0^+$ and $\lambda\to 0^+.$
  We derive an evolution system of fully non-linear PDEs as gradient flow of
$\mathcal E_{\lambda,\delta}$ in $L^2$, showing existence and uniqueness of
solutions. All the trajectories of the flow converge in $W^{1,p}_0$ to the
unique minimizer $\mu_{\lambda,\delta}^*$ of $\mathcal E_{\lambda,\delta}.$
  Finally, we characterize $\mu_{\lambda,\delta}^*$ by a non-linear system of
PDEs which is a perturbation of Monge-Kantorovich equations by means of a
p-Laplacian.
