Let
$$F(x):=\mu((-\infty,x])=\sum_i a_i\,1(x_i\le x)
=\sum_{j=1}^n s_j\,1(x_{n:j}\le x<x_{n:j+1}),$$
where $x_{n:1}<\cdots<x_{n:n}$ are the values $x_1,\dots,x_n$ put in the increasing order (with $x_{n:n+1}:=\infty$),
$$s_j:=\sum_{i=1}^j a_{n:i},$$
and $a_{n:1},\dots,a_{n:n}$ are the values $a_1,\dots,a_n$ put in the increasing order of the $x_k$'s, so that, if $x_{n:i}=x_k$ for some $k$, then $a_{n:i}=a_k$.
So, $F$ is the cdf of the probability measure $\mu$. Similarly considered is the function $G$ defined as the cdf of the probability measure $\nu$.
Consider then the generalized inverse/quantile function $F^{-1}\colon(0,1)\to\mathbb R$ defined by
$$F^{-1}(u):=\inf\{x\in\mathbb R\colon F(x)\ge u) \\
=\max\{x\in\mathbb R\colon F(x)\ge u) \\
=\sum_{j=1}^n x_{n:j}\,1(s_{j-1}<u\le s_j)$$
for $u\in(0,1)$, with the similarly defined and considered $G^{-1}$. Let $U$ be a random variable (r.v.) uniformly distributed on $(0,1)$. Then the distributions of the r.v.'s $X:=F^{-1}(U)$ and $Y:=G^{-1}(U)$ will be $\mu$ and $\nu$, respectively. Finally, let
$$M(x_1,\dots,x_n,y_1,\dots,y_n,a,b):=M(\mu,\nu) \\
:=(E|X-Y|^p)^{1/p}
=\Big(\int_0^1|F^{-1}(u)-G^{-1}(u)|^p\,du\Big)^{1/p}.$$
Then
$$W_p(\mu,\nu)\le M(x_1,\dots,x_n,y_1,\dots,y_n,a,b)$$
and
$$M(x_1,\dots,x_n,x_1,\dots,x_n,a,a)=0,$$
as desired.
Remark: The upper bound $M(\mu,\nu)$ given above is actually the exact value of the Wasserstein distance for $p\ge1$, according to the last sentence of Theorem 2.1 -- thank you alesia for this reference.
