A Neumann boundary condition simply gets into the system of equations through the Green's first identity when applied to the Laplacian term, i.e.:
$$ \int_L u_{xx} \, \omega \, dx = u_x(1) \omega(1) - u_x(0) \omega(0) - \int_L u_x \omega_x \, dx, $$ where $\omega$ is the test function used in your FEM approach. Since you have $u(1) = b$, you may enforce $\omega(1) = 0$, but on the other hand, for the VN b.c. just substitute it back in the previous equation, yielding:
$$ \int_L u_{xx} \, \omega \, dx = -a \, \omega(1) - \int_L u_x \omega_x \, dx, $$ which will be further simplified if you conveniently choose $\omega(x)$ such that $\omega(0) = 1$, i.e., for example choosing $\omega(x)$ as the family of piecewise linear polynomials in $[0,1]$, known as hat functions. Therefore, we would have:
$$u(x) = \sum_{i=0}^N u_n \omega_n (x) = u_N \omega_N(x) + \sum_{i=0}^{N-1} u_n \omega_n(x), $$ where $u_n$ are the unknowns. The Dirichlet boundary condition tells us:
$$u(1) = u_N \cdot 1 + 0 = b \Rightarrow u_N = b, $$ and we have reduced the order of the system by 1, as you pointed out before. Notice now that substituting the VN b.c. would tell us nothing and therefore $u_N$ remains unknown.
I hope this helps.
Cheers!
Reference
I recommend Classical Potential Theory by Armitage and Gardiner, section 5.2. Corollary 5.2.3 implies that if a bounded function is harmonic in $\Omega\setminus F$ where $F$ is a finite set, then it extends as a harmonic function in $\Omega$. In your setup, the extended function would contradict the maximum principle.
Proof
I will use Green's function, as you expected. If $n>2$, form the sum
$$
f(x) - \epsilon\sum_{i=1}^k \|x-x_k\|^{2-n}
$$
Note that $f$ tends to $-\infty$ at the points $x_i$. Letting $f(x_i)=-\infty$, we obtain a subharmonic function in all of $\Omega$. Indeed, one characterization of subharmonic functions is the local sub-mean value property: for every point $a$ there is $r_a>0$ such that the average of $f$ on every sphere $S(a,r)$ with $r<r_a$ is at least $f(a)$. This property holds at $x_i$ trivially, and at other points by the harmonicity of $f$.
Subharmonic functions satisfy the maximum principle. For a fixed $p\in \Omega\setminus \{x_i\}$ it says
$$
f(p) \le O(\epsilon)
$$
and since $\epsilon$ was arbitrarily small, $f\le 0$ in $\Omega$. The same consideration applies to $-f$; thus, $f\equiv 0$.
For $n=2$ the proof is the same, but has the logarithm in it:
$$
f(x) + \epsilon\sum_{i=1}^k \log\|x-x_k\|
$$
Other PDE
The keyword is capacity: does a one-point set have zero capacity for the given differential operator? The Laplace equation is the Euler-Lagrange equation for the functional $\int|\nabla f|^2$, so the relevant capacity of point $p$ is
$$
\inf \left\{\int|\nabla f|^2 : f=0\text{ on boundary, } f(p)=1 \right\}
$$
The infimum can be taken over all smooth functions, or over all Lipschitz functions: does not matter. The capacity of a point is zero, as you can see by letting $f$ be a scaled and truncated version of Green's function. Since the solution of boundary value problem is meant to minimize the functional for the given boundary data, the fact that the infimum is zero says there isn't a solution.
For the biharmonic equation, the relevant functional is $\int|\Delta f|^2$. Is it true that
$$
\inf \left\{\int|\Delta f|^2 : f=0\text{ on boundary, } f(p)=1 \right\} =0 \quad ?
$$
Yes in dimensions $n\ge 4$, because biharmonic Green's function is unbounded there (source). No in dimensions $n\le 3$. The reason is that the $L^2$ norm of the Laplacian controls $W^{2,2}$ norm (roughly speaking), and in dimensions $n\le 3$ the Sobolev embedding implies that $W^{2,2}$ norm controls pointwise values of a function.
So, the biharmonic equation can be used for interpolation in low dimensions. In two dimensions, this is known as a thin plate spline. More generally, there are polyharmonic splines: the higher the dimension, the more one iterates the Laplacian to keep Green's function bounded.
Best Answer
The answer is negative. This is not something I can prove, but if such decoupling was available, it would be all over the numerical PDE literature. Instead, I see:
in A highly accurate numerical solution of a biharmonic equation by M. Arad, A. Yakhot, G. Ben-Dor.
Guo Chen, Zhilin Li, and Ping Lin in A fast finite difference method for biharmonic equations on irregular domains... also comment on the difficulty of dealing with boundary conditions in your post. They develop the following method: