@mmcdara Thanks for your advice. However, I find some strange things.
In the first approach, because

it is enough to search for a real instance satisfying the boolean expression. Strangely, the result is not correct.
As for the second approach, I attempt to resolve the conditions via brute-force methods: RealDomani[solve](And(RealDomain[solve](eig_1 >= 0), RealDomain[solve](eig_2 >= 0), RealDomain[solve](eig_3 >= 0), …)). Strangely, the result is still incorrect. 
 

Download PSD_BUG.mw

 Is there something wrong?

@Carl Love Thanks for knowledge. The documentation of Alias gives me a feeling of being somewhat dangerous. Besides, some new commands in the ArrayTools package seem to compromise on efficiency (maybe for the sake of simplicity?).

@C_R Maple code is as follows: 

interface(printbytes = false): 
with(SignalProcessing, FFT): 
with(LinearAlgebra, HankelMatrix):
for n from 0 to 15 do
    m:=HankelMatrix(<$1+1..2**n+2**n>,datatype=complex[8],shape=[]);
    gc();print(n,time[real](FFT(m,normalization=none,inplace=true)))
od: # shape=[] is not optional

Now you can run Python directly (cf. Connectivity in Maple 2023); nevertheless, FFT is not built into Python.

from os import environ
environ['JAX_ENABLE_X64']="1"
import timeit,jax.numpy as jnp
for n in range(0x10):
	m=sum(jnp.indices((2**n,2**n),sparse=True,dtype='cdouble')+(1,1))
	print(n,timeit.timeit(lambda:jnp.fft.fft2(m),number=1),sep="，")

Here the better choice is using a standalone release.

@Carl Love "evalf[4](…)" is not enough to obtain the desired result. For example, the entry corresponding to 0.0 and 0.09 should be "0.5359" instead of "0.5358". So one has to use "evalf[4](evalf[4+1](…))" in the anonymous function.

@C_R Strangely, although the help page writes that

this tutorial (updated on 10/26/2015) says that "there are 1D and 2D versions (to compute for real or complex data) in the library".
I believe that the documentation has some typos.

@C_R Thanks for your suggestion. But the Signal Processing package only works (and hence, only optimizes) on coerced Vectors with compatible datatype, because the documentation writes that 

Strangely, although it is said that "the code underlying these commands does not support Arrays of dimension two or higher directly", the Intel^® IPP help page says that "there are 1D and 2D versions in the library".

@TechnicalSupport Many thanks! Has this problem (concerning the QuantifierElimination package) been fixed as well?

@vv I am in agreement with you.
The QuantifierElimination command (since Maple 2020) in the SemiAlgebraicSetTools subpackage of RegularChains appears to be not experimental now. Unfortunately, it seems that neither QuantifierElimination:-QuantifierEliminate nor RegularChains:-SemiAlgebraicSetTools:-QuantifierElimination can effectively eliminate quantifiers (because of the high complexity?). 

Note. In modern Mathematica, the elimination is often much faster (Why??!): 

Mathematica 13.3.0 Kernel for Microsoft Windows (64-bit)
Copyright 1988-2023 Wolfram Research, Inc.

In[1]:= EchoTiming[Resolve[ForAll[{x,y,z},(x|y|z)\[Element]NonNegativeReals\[And]x+y+z\[GreaterEqual]x*y*z-2,12*K+(x+y+z-(-1))^2\[GreaterEqual]K*(y*z+z*x+x*y)+(6-(-1))^2],WorkingPrecision\[RuleDelayed]21+1],Method\[Rule]RepeatedTiming]
>> 10.458

Out[1]= K == 4

In[2]:= EchoTiming[Resolve[ForAll[{x,y,z},(x|y|z)\[Element]NonNegativeReals\[And]x+y+z\[GreaterEqual]x*y*z-2,12*K+(x+y+z-(+2))^2\[GreaterEqual]K*(y*z+z*x+x*y)+(6-(+2))^2],WorkingPrecision\[RuleDelayed]17+1],Method\[Rule]RepeatedTiming]
>> 5.56633

Out[2]= 2 <= K <= 3

I surmise that this is because certain core components call optimized and compiled 𝙲 code rather than execute shallow native 𝙼𝚖𝚊 code ... I suppose so.

Alas, every coin has two sides; this may also result in uncommon inefficiency (as one has seen here).

@tomleslie Well, I do not believe that this is just a “feature”, and yet Maple 2023 still cannot solve this more carefully.

@vv It is fortunate to see that solve solved those two systems at last (though one has to await them for enough long (!)). Nevertheless, the output of the first system is too complicated: 

Download half_an_hour.mw

As you have expected, the certified result is simply . I cannot know why Maple did not simplify the returned value in a long time, which only leads to an unreadable (and more or less worthless) solution!

@C_R Sorry for misunderstanding. I don't mean that "a = 2/(sqrt(5) + 1) should be treated as a separate solution"; I just want to stress that the output should be “a >= 2/(sqrt(5)+1)” instead of "a > 2/(sqrt(5)+1)". (Isn't this a bug?) 
Besides, I find that this bug can even be reproduced in an old version! 

# in Maple 1x
> solve({a > 0, ln(a)+ln(1+a) >= 0}, a);
 = 
           /      /  2                    \        \ 
          { RootOf\_Z  - _Z - 1, index = 1/ - 1 < a }
           \                                       / 

It seems that Maple's developers still don't know this ….

@mmcdara Thanks for your trick. However, it appears the output is still not very accurate when using larger Digits. (And in many cases, it is even nearly unexpected.) The documentation writes that this package "takes advantage of built-in library routines provided by NAG", but I don't know why it cannot return a result with high precision.

@Axel Vogt As for the "source", it comes from a difficult algebraic inequality. But this question is another thing (as I do not need an exact supremum here).

@acer Thanks for your kind reply. You mention the Optimization package. In my view, this toolbox is much less smart than solve/float, evalf/Int, dsolve/numeric, and evalf/Sum (for some unaccountable reason). One can often find (too) many warning messages and errors during a computation. Well, there is a reasonable prospect that it is no longer actively developed (though is still being maintained (just like the algsubs command)). What a pity.

@Carl Love Many thanks! 
You said that "Maple uses dynamic scoping", but recently, I noticed the following sentence in What Has Changed in Maple V Release 5: 

Does this mean that Maple uses a mixed version instead?

@nm You say: "It goes without saying that listing all the proc names in trace([foo,....],statements = false) is not a real solution.". But the following code does work: 

foo:=proc(n,m,k)   
   local r;  
   #if DO_TRACE then
   #   option trace=0:
   #fi;
   
   boo(1);
   r:=n:
   NULL;
 end proc:

boo:=proc(n)   
   local r;  
   r:=n:
   NULL;
 end proc:

trace~([foo, boo], statements = false):
foo(1, 2, 3);
{--> enter foo, args = 1, 2, 3
{--> enter boo, args = 1
<-- exit boo (now in foo) = }
<-- exit foo (now at top level) = }

Is there any misunderstanding?