I need to find four orthogonal linear combinations of complicated functions, that vanish at four different points.

I use

NMinimize[f[{a, b, c, d}, 0, 0], {a, b, c, d}]

and get a result, as replacement rule (using MachinePrecision)

{3.3493*10^-17, {a -> 2.34594, b -> -1.80385, c -> 2.51873, d -> 0.2406}}

Now if I create a vector for these coefficients

P1 = {a, b, c, d}/.%[[2]]

and call the function with this vector I get a different result than by replacing directly in the function call:

f[{a,b,c,d}, 0, 0]/.%%[[2]] returns 3.3493*10^-17, while f[P1, 0, 0] returns 0.373542.

What's worse is that f[({a, b, c, d}/.%%%%[[2]]), 0, 0] also returns the inaccurate result 0.373542.

What I wanted to do is find a linear combination that minimizes at the first spot, then define a coefficient vector from that and minimize at the second spot with the constraint that the coefficient vector has to be orthogonal to the first. f[{a, b, c, d}, x, y] normalizes the coefficient vector btw.

How can I get Mathematica to assign appropriately precise values upon using the replacement rule?

I tried setting the precision to 100 instead of MachinePrecision, but with the very same results.