You can also use command SolveEquations from DirectSearch package to solve the problem.

with(DirectSearch);
interface(rtablesize=100):
sol:=SolveEquations(a^2+b^2+c^2=d^2,{a>=1,b>=1,c>=1,d>=1,d<=100},assume=integer,AllSolutions,solutions=20);

You can disable all warnings by command: interface(warnlevel=0);

@Markiyan Hirnyk Because objective function has many local minimums we can use global optimizers GlobalOptima and (or) GlobalSearch from DirectSearch package. These optimizers take very long time to optimize, but they solve the problems:

sol:=DirectSearch:-DataFit(F,X,Y,x,optimizer=globaloptima,pointrange=[a=-3..3,b=-3..3,c=-3..3,d=-3..3,e=-3..3]);

sol:=DirectSearch:-DataFit(F,X,Y,x,optimizer=globalsearch,pointrange=[a=-3..3,b=-3..3,c=-3..3,d=-3..3,e=-3..3],solutions=5,evaluationlimit=2000);

The sum of the squared residuals being about 0.00006 for both global optimizers.

The objective function (that optimized) has many local minimums. So we need a good initial point. For example:

startpoint:=initialpoint=[a =3,b=2,c=3,d=5,e=0.7];
sol1:=DirectSearch:-DataFit(F,X,Y,x,startpoint);

Try OrthogonalExpansions package to approximate univariate and multivariate functions.

First variant.
    1. Find the directory where your Maple installed.
    2. Find the subdirectory \lib\
    3. Copy the following two files into this subdirectory:
        DirectSearch.mla – Maple Archive File,
        DirectSearch.hdb – Maple Help file.

Second variant.
    1. Create the directory, for example, C:\DirectSearchLib\
    2. Copy the following two files into this directory:
        DirectSearch.mla – Maple Archive File,
        DirectSearch.hdb – Maple Help file.
    3. In current Maple session execute the command:
        > libname:=”C:\\DirectSearchLib”, libname;

Since the objective function ObjFn depends on 6 variables and the constraints ModCons depend on 14 variables we must set all variables explicitly. So the local solution:

var:=[op({op(indets(ObjFn)),op(indets(ModCons))})];

sol:=Search(ObjFn,[op(ModCons),op(Bounds)],startpoint,variables=var,tolerances=10^(-8),checkexit=10,maximize);

eval(ModCons, [seq(var[i]=sol[2][i],i=1..nops(var))]);

Try local solution several times because I obtain each time different solutions - due to deterministic chaos, I think.

The multimodal solutions:

sol2:=GlobalSearch(ObjFn,[op(ModCons),op(Bounds)],pointrange=Bounds,solutions=5,number=10,variables=var,tolerances=10^(-8),checkexit=10,maximize);

eval(ModCons, [seq(var[i]=sol2[1][2][i],i=1..nops(var))]);

The global solution (very very slow):

sol3:=GlobalOptima(ObjFn,[op(ModCons),op(Bounds)],pointrange=Bounds,variables=var,tolerances=10^(-8),checkexit=10,maximize);

eval(ModCons, [seq(var[i]=sol3[2][i],i=1..nops(var))]);

It is hard root to find by numerical methods. For example, how to find real root of the following equation:
      piecewise(x<> Pi, 1+x^2, 0)=0
by pure numerical methods ? The pure numerical methods must fail even if time of computation is infinity, because the set of real numbers is continuum. Remedy is symbolic or mixed symbolic-numeric methods:

solve(piecewise(x<> Pi, 1+x^2, 0)=0);

with(OrthogonalExpansions):

velocity:=piecewise(t<=6, 3*sin(t*Pi/6), t>6, 0);
n:=10; # number of fourier series terms
Fourier:=FourierSeries(velocity,t=0..20,n);
plot([velocity,Fourier],t=0..20);

fsolve command gives infinity number of solutions for various b:

fsolve([2.32=a+b,pa=.4310344828*a,pb=.4310344828*b,pa+pb=1]);

{a = 2.32 - 1. b, b = b, pa = 1.000000000 - 0.4310344828 b, pb = 0.4310344828 b }

I try solve problem by DirectSearch package v.2. I cannot find any problem to maximize function f and to solve equation df=0.

with(DirectSearch);

seq(Search(eval(f,ga=m),maximize),m=0.1..0.8,0.1);

SolveEquations(df=0);

SolveEquations command from DirectSearch package finds all four solutions in one step:

DirectSearch:-SolveEquations(x^3-3*2^x+3=0, AllSolutions);

The DirectSearch package finds two solutions:

with(DirectSearch):

T1:=abs(subs(s=x+I*y,T));

sol:=GlobalSearch(T1, tolerances=10^(-14), distance=10^(-6));

or

sol:=SolveEquations(T1, AllSolutions, distance=10^(-6));

It seems there are only two solutions, see plot:

plot3d(T1,x=0.185..0.195,y=-0.005..0.005);

Firstly, when your model function is procedure the first parameters must be independent variables followed by function parameters (see Help for DataFit). So procedure FuncH:=proc(p,pi,d,ttt) ... you should replace by FuncH:=proc(ttt,p,pi,d)...

Secondly, in the first variant, when function F is used, the procedure H := PDE(p, pi, d) is computed ONE TIMES for all 80 points pkt. In the second variant, when function FuncH is used, the procedure H := PDE(p, pi, d) is computed 80 times for all 80 points pkt.

Therefore, the second variant must be approximately 80 times slowler than the first variant.

Increase Digits to round(evalhf(Digits)).