You need to use algsubs instead of subs:

eqn:= N*x*y - x*y + f(x,y) = 0;
algsubs(x*y= 1, eqn);

Both subs and eval (and even subsindets and evalindets) can only substitute for subexpressions that occur as syntactic subunits in the main expression, but algsubs can find algebraic subunits. It would require a lengthy digression into the minutiae of computer algebra for me to precisely define syntactic subunit, but if you'd like, I'll try later. I think an example will suffice as a definition for now: Let's consider the present expression, but written in functional form (I'd also call it a tree form):

`=`(`+`(`*`(N, x, y), `*`(-1, x, y), f(x,y)), 0);

The above is not just a pedantic presentation form, but is actually executable Maple code that produces the same eqn as the original. You will note that our target expression `*`(x*y) doesn't appear in the above functional form, which is why subs doesn't see it.

The above functional form is very close to (but not identical to) Maple's internal representation of the expression. In reality, the numeric coefficients are handled a little bit differently, which is why subs does find the x*y in `*`(-1, x, y).

For the purpose of checking associativity, we can ignore the phi on the left side of the equation. Then do this:

`&*`:= (u__1,u__2)-> ln(exp(phi(u__1))+exp(phi(u__2))) + Q(phi(u__1) - phi(u__2));
AssocLeft:= (u__1 &* u__2) &* u__3;
AssocRight:= u__1 &* (u__2 &* u__3);
CriticalEqn:= AssocLeft - AssocRight = 0;

So, to check any particular Q, we substitute it into the CriticalEqn, along with some arbitrary phi and u's. Here, I let Q and phi be the identity function:

eval(CriticalEqn, [Q= (x->x), phi= (x->x), u__1= 1, u__2= 0, u__3= -1]);
evalf(%);

The left side isn't even close to 0, which shows immediately that it's not associative for this Q.

Regardless of whether the equations are linear or the solutions are exact, you can achieve what you asked for like this:

eval(Vector[row](nops(SOL[0]), symbol= tau), SOL[0]);

Here is code to exactly duplicate in Maple 2019 the functionality of rept that you show in your worksheet:

rept:= (f, R::(name= range(realcons)))-> 
   ((t,a,r)-> eval(f, t= t - r*floor((t-a)/r)))((lhs, (lhs, rhs-lhs)@~rhs)(R))
: 

Your example:

F:= rept(piecewise(t<1, t/2+1/2, t>=1, -t+2), t= -1..2);
F := piecewise(t - 3*floor(t/3 + 1/3) < 1, t/2 - (3*floor(t/3 + 1/3))/2 + 1/2, 1 <= t - 3*floor(t/3 + 1/3), -t + 3*floor(t/3 + 1/3) + 2)

Note that my code gives an exact symbolic representation (in algebraic form) of the periodic extension, not merely a procedure or a function only suitable for plotting.

I think that you're saying that 1 shouldn't be a solution under RealDomain. I agree with you, and I think that this is a bug.

In order to use RealDomain properly, you should use it with with(RealDomain) so that the operations `^` and sqrt use the versions from the package. However, doing this, I still get 1 in the solution set.

The following worksheet shows DirectSearch performing superbly when applied to a black-boxed Ackley function where I've hidden the minimum at a random location to make it more challenging.
 

Download Ackley.mw

I will need to see your code to diagnose the problem. But my first guess of the cause of your problem is that you're storing large intermediate results pursuant to your complicated calculation in indexed variables. Thus, the memory that they use can't be re-used. In our vernacular we say that they aren't being garbage collected. The solution may be as simple as changing indexed variables to plain ones.

The diff_table command allows you to use indexing (not subscripting) to specify the independent variables with respect to which derivatives are taken. Indices in Maple are placed in square brackets on input, and appear as subscripts in the prettyprinted output^[*1]. So, you must do like this:

restart:
with(PDETools):
U:= diff_table(u(x,t)):
declare(U[]);
PDE:= U[x,x] - U[t] = 0;  #Note the square brackets!
Infinitesimals(PDE);

Your input  U_x,x - U_t = 0 is pretty much nonsense to Maple in a PDE context. Sometimes it's unfortunate that Maple allows you to enter those things without complaining.

[*1] As opposed to the prettyprinted output, you can view the raw form, plaintext output of any command by doing 

lprint(%);

after the command. Or you can apply it to a name, such as 

lprint(PDE);

Then you'll see the square brackets. This can often provide a "reality check" when you're baffled by some output.

The answer to your immediate problem is that you need to use expand to force through the multiplication. This case and more-general cases can be done by 

numer(e)  #numerator

where e is the expression that you want to work on. 

While it's not directly related to your specific question, I thought that it's worth pointing out that your procedure can be vastly simplified to 

MinDistPoint:= proc(
   A::[algebraic, algebraic], B::[algebraic, algebraic], P::[algebraic, algebraic]
)
description "Point on line AB at minimum distance from P";
local D:= <A-B>, R:= <0,-1;1,0>, V:= <<P>.D, -R.<A>.<B>>;
   [V.D, R.V.D] /~ (D.D)
end proc:

There are probably numerous other formulations of similar simplicity. This was just the first that I derived.

Local variables inside procedures are ordinarily only evaluated one level deep. Your H and lambda are local variables. In the last line of the procedure, one level of evaluation turns H into [3*lambda, 4*lambda-3] without substituting the value for lambda, even though it's already known. This behavior is usually harmless, but it can be annoying when the partially evaluated value is returned to the top-level. The correction is to make the last line of the procedure eval(H). 

Note that if you don't use the eval, it will happen by itself anyway the next time that the result is used. For example:

H:= PP(X, Y, Z);
                 H := [3 lambda, 4 lambda - 3]
H;
                            [57  1 ]
                            [--, --]
                            [25  25]

A place to start might be the package Interpolation:-Kriging or something else in package Interpolation. If the system to be optimized can be expressed in terms of numeric procedures, it can then be handled by the Optimization or DirectSearch packages. These packages do not require that systems be expressed algebraically.

Like this:

HasProduct:= proc(e, p::`*`)
local P;
   has(algsubs(p= P, e), P)
end proc:

My solution uses Iterator:-CartesianProduct on the sets returned by NumberTheory:-Divisors. It works on polynomials with rational coefficients; they need not be integer. 

RationalRoots:= proc(F::polynom(rational))
description "A straightforward application of the rational root theorem.";
uses It= Iterator, NT= NumberTheory;
local f:= expand(F), r;  
   {seq(
      ({-1,1} *~ `/`(seq(r)))[], 
      r= It:-CartesianProduct(      
         ((S->{-1,S[]})@NT:-Divisors@~(tcoeff,lcoeff))(ilcm(denom~([coeffs(f)])[])*f)      
      )
    ), 
    0
   }
end proc
:   

What is your ultimate goal? 

Suppose that A is a numeric Matrix, generated from your function of two variables (z,x)-> z*x. Then what's wrong with simply using A[2,3] to refer to the the 3rd item in the 2nd row?

I think that there's a better way than StringTools to do any numeric computation.

Is there anything about the following code that doesn't do what you want?

M:= (z,x)-> z*x;
(Zres, Xres):= (5,9);
Z:= Array(1..Zres, i-> i);
X:= Array(1..Xres, i-> 1+0.5*(i-1)); 
A:= Array(1..numelems(Z), 1..numelems(X), (i,j)-> M(Z[i],X[j]));