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]));

restart:

#PDE system specification:
f:= x-> 1+x:
alpha:= t-> 1+t:
t0:= 0:
x0:= 0:
pde:= {diff(u(t,x),t)= -diff(u(t,x),x) + 2 + t + x}:
ibc:= {u(t0,x) = f(x), u(t,x0) = alpha(t)}:

#Solve the PDE, as a basis of comparison:
pdsolve(pde union ibc);
                          u(t, x) = (1 + x) (1 + t)

#Picard iteration:
U:= eval(u(t0,x), ibc);  #Initialize iteration.
to Order do  #Order is an envirnoment variable that defaults to 6; it can be changed.
   U1:= eval(u(t0,x), ibc) + 
      int(eval(eval(eval(diff(u(t,x),t), pde), u(t,x)= U), t= s), s= t0..t);
   U1:= U1 - (eval(U1, x= x0) - eval(u(t,x0), ibc));
   if U=U1 then break fi;  #Fixed point reached.
   U:= U1  #No need to save the iterates.
od;
                                 U := 1 + x

                                          1  2      
                        U1 := 1 + x + t + - t  + x t
                                          2  
       
                            U1 := t x + t + x + 1
                            U := t x + t + x + 1
                            U1 := t x + t + x + 1
                            U1 := t x + t + x + 1

RemoveEveryNth:= proc(L::list, n::posint)::list;
description "Return the list L with every nth entry removed.";
   L[remove(k-> irem(k,n)=0, [seq(1..nops(L))])]
end proc:

Integrating a constant from -1 to 1 is equivalent to multiplying it by 2. That the constant is a matrix makes no difference mathematically.

Looking at your PDE system (as output), I see some terms like this:

svt1(x, t)(x, t) - sht1(x, t)(x, t)

Note the extra (x,t). This means that somewhere you got functions and expressions mixed up.

Consider these two lines from your code:

mubt1:= (x,t)-> (-5/3)*(3*ks+4*mus)/(3*ks*(5*Jirt1(x,t)-2*(1-phi0))+4*mus*(5*Jirt1(x,t)-3*(1-phi0)));
#...
Fht1:= ((svt1(x,t)+2*sht1(x,t))*Kbt1-(svt1(x,t)-sht1(x,t))*mubt1);

I have omitted two intervening lines that aren't related to this error.  As you can see from the parts that I highlighted, the first line makes mubt1 a procedure, but the second line uses it as if it were an algebraic expression. This is the type of thing that led to your error message. You may have done this elsewhere in the code. I've only highlighted here the first instance that I noticed.

You might notice these errors yourself if you didn't insist on ending every line with a colon.