Variational inequalities and their relation to weak formulations of PDEs

Question

this topic is continuously boggling me as I go through material of a class on nonlinear PDEs. I will illustrate the broader issue with the following example. Consider the reaction diffusion equation involving p-Laplacian: $-\Delta_p u + b(u) = 0$ in $\Omega$ and $\nabla_p\cdot \nu + d(u) = 0$ on $\partial\Omega$, where we have some usual conditions on function $b$ and $d$ (Caratheodory and suitable growth bounds). After a rather vague discussion it is being concluded that for some closed convex $M \subset V = W^{1,p}(\Omega)$ there exists a unique solution $u \in M$ to the variational inequality (V.I.) $$ 0 \leq \int_{\Omega} [\nabla_pu\cdot\nabla(v-u) + b(u)(v-u)] + \int_{\partial\Omega} [d(u)(v-u)] \text{ for all } v \in M. $$ This looks very much like the weak formulation (W.F.), which is $$ \int_{\Omega} [\nabla_pu\cdot\nabla(v) + b(u)v] + \int_{\partial\Omega} [d(u)v] = 0 \text{ for all } v \in V. $$ Here is where my lack of conceptual understanding comes in: how are solutions of the V.I. and W.F. related? There appears to be some sort of equivalence lurking around but I cannot really get my head around it. What does knowing that there exists a unique solution of V.I. tells us about the W.F. (or the PDE in the classical sense)? Does the fact that we restrict ourselves to a convex subset $M$ not make it a separate problem? I can see the equivalence more clearly when the PDE we are dealing with is linear and we define the bilinear form $B(u,v)$ and then a functional $J(v) = \frac{1}{2} B(u,v)$ (if RHS is zero) and take $M = V = H^1_0(\Omega)$ (say), but other than that I am unfortunately in the dark. In particular, I am not sure about the relation to the original PDE (or lack of thereof) if we work with a strict subset of $V$ in the linear cases. Apologies if the question is obvious/a duplicate, but I have tried googling, looking in Evans etc, but found no satisfying answer. I will be very grateful for any insight. Cheers.