@Markiyan Hirnyk 

What about the OP's: " there is the obvious solution X=[e_3,e_1,e_2] where (e_i) is the canonical basis of the space".
Corollary ..., emptyness ...
 

@Markiyan Hirnyk 

The existence was known. It's about the method; it does not work.
If you like your solution (including warnings and/or errors) , it's ok.
 

@Markiyan Hirnyk 

Solid people should check their own "product". Anyway, for your convenience:

restart; with(LinearAlgebra):
X := Matrix(3, 3, symbol = x):
Y := Transpose(X) . X-IdentityMatrix(3):
F := NULL: for i to 3 do for j from i to 3 do F := F, Y[i, j] = 0 end do end do:
sol := [solve([F], explicit)]:
nops(sol);
                               64
solve({x[1, 1] < 1/2, 3/10 < x[2, 3], sol[62][]});
Warning, solutions may have been lost
 ...

@Markiyan Hirnyk 
It produces several "Warning, solutions may have been lost".  And lost they are.

@petit loup 

In my opinion solve or  SolveTools:-SemiAlgebraic  (which is called by solve) are not good tools for your problem. 
Just for fun I have reduced tour toy example to a system

[0 <= 3*a^2+3*b^2-c^2-1, 0 <= 3*a^2-20*a*b+3*b^2+3*c^2-20*c+3]

and  both commands failed for it (I had to interrupt them).

If your specific problem is to find an orthogonal matrix in a neighborhood of a given matrix in GL(n,R)  then I'd suggest to use the Cayley transform (see wiki). This way, the problem reduces to find an antisymmetric matrix in the nbd of a given matrix, which is a simple task (for any n).

@Markiyan Hirnyk 

What explanation are you expecting? We have a good approximate solution and you have checked it. Of course ANY approximate solution could be invalid (in very rare situations).
BTW, you probably know that solve itself uses approximations to isolate roots etc.

@Markiyan Hirnyk 

You should know the difference between 1d and 2d input. You have copied my 1d code in 2d mode.

@Markiyan Hirnyk 

I know the difference between numeric and symbolic solutions. But try X:=Matrix(n,symbol=x)   for n=4 or 5. The numeric solution still works. 
Integer solutions? They were not requested for this problem.

@Markiyan Hirnyk 

Of course it's infinite. The problem was about the existence of real solutions. Good luck too.

@petit loup 

For now we must be happy with a numerical solution. It is easy and fast:

restart;
Digits:=30:
X:=Matrix(3,symbol=x):
eqs:={entries((X.X^+ - X^0) =~ 0, nolist)}[]:
Optimization:-NLPSolve(0, {eqs, x[1,1] <= 1/2, 3/10 <= x[2,3]});

[0., [x[1, 1] = 0.794146063666383512119841132959e-1, x[1, 2] = .620068305692267153579926862637, x[1, 3] = .780518171839421696958016866613, x[2, 1] = .822632510153236016792344665016, x[2, 2] = -.482986286802209366562016187596, x[2, 3] = .300000000000000000000000000000, x[3, 1] = -.563000065306051208096993503723, x[3, 2] = -.618255241010486777840153327601, x[3, 3] = .548444512624975566069664290363]]

@Muhammad Usman 

Replace exp(k2*eta)   by  Z^k2. The problem reduces to find the coefficients of a polynomial in the variables x, y, t, Z and you can use coeff as above.

@Carl Love 

I don't understand. You said (correctly) "The problem is not caused by using a variable whose name begins with underscore" and then you provide a workaround exactly for such variables (only).

restart;
assume(Z1::integer);
sys:={x <> (  (1/2)*Pi+Z1*Pi ), x < infinity, -infinity < x};
expand(thaw(solve(subsindets(sys, suffixed(_), freeze@``) ,{x}))); #error

Also, the freezed variable is without assumptions; so, a simpler woraround would be to simply remove the assumptions which are anyway ignored now.

@Carl Love 

Yes, you are right; and for other types too, e.g. positive.

@Preben Alsholm

@digerdiga 

Using such strange names for variables (for an outsider) your code looks as being obfuscated and very hard to follow.
Why don't you isolate the problem and use normal (short) variables?