$\newcommand{\Si}{\Sigma}$Your first question: "Is my projection optimal (not necessarily unique)?" The answer to this is -- Of course, not. Moreover, your projection is usually not even in $S$.

Indeed, as was finally clarified in your comments, we have
\begin{equation*}
S=S_1\cap S_2,
\end{equation*}
where
\begin{equation*}
S_1=\{x\in\{0,1\}^{n\times n}\colon x_{i,j}+x_{j,i}\le1\ \forall(i,j)\in[n]\times[n]\},
\end{equation*}
\begin{equation*}
S_2=\{ x\in \{0,1\}^{n\times n}\colon\sum_{j=1}^{n}x_{i,j}=1\ \forall i\in[n]\},
\end{equation*}
and $[n]:= \{1,...,n\}$.

Let $Px$, $P_1x$, and $P_2x$ denote the projections of $x$ onto $S$, $S_1$, and $S_2$, respectively, so that your projection of $x$ is $P_{12}x:=P_2(P_1x)$.

So, your first question is this: Is it true for all $x\in[0,1]^{n\times n}$ that
\begin{equation*}
Px=P_{12}x\text{?} \tag{1}\label{1}
\end{equation*}

The (nonlinear) projection $P_1x$ of $x\in[0,1]^{n\times n}$ is described quite simply:
\begin{equation*}
(P_1x)_{i,j}=x_{i,j}\,1(x_{i,j}>x_{j,i})
\end{equation*}
for all $(i,j)\in[n]\times[n]$; here we exclude the "zero-probability" case when $x_{i,j}=x_{j,i}$ for some distinct $i$ and $j$ in $[n]$.

Note that the matrices $x\in S_2$ are in the bijective correspondence with the $n$-tuples $s$ in the set $[n]^n$, where the correspondence is given by the formula
\begin{equation*}
x_{i,j}=X(s)_{i,j}:=1(j=s_i)
\end{equation*}
for all $(i,j)\in[n]\times[n]$, so that $s_i$ is the position of the (only) $1$ in the $i$th row of the matrix $x\in S_2$. The square of the Euclidean distance $d(x,X(s))$ from a matrix $x\in[0,1]^{n\times n}$ to the matrix $X(s)$ corresponding to an $n$-tuple $s\in[n]^n$ is
\begin{equation*}
\begin{aligned}
d(x,X(s))^2&=\sum_{(i,j)\in[n]\times[n]}(x_{i,j}-X(s)_{i,j})^2 \\
&=\sum_{(i,j)\in[n]\times[n]}(x_{i,j}-1(j=s_i))^2 \\
&=\sum_{i\in[n]}\Big((x_{i,s_i}-1)^2-x_{i,s_i}^2+\sum_{j\in[n]}x_{i,j}^2\Big) \\
&=n+\sum_{(i,j)\in[n]\times[n]}x_{i,j}^2-2\sum_{i\in[n]}x_{i,s_i}.
\end{aligned}
\end{equation*}

So, projecting a matrix $x\in[0,1]^{n\times n}$ onto $S_2$ or onto $S$ is equivalent to maximizing
\begin{equation*}
\Si_x(s):=\sum_{i\in[n]}x_{i,s_i} \tag{2}\label{2}
\end{equation*}
over all $n$-tuples $s\in[n]^n$ corresponding to the matrices $x$ in $S_2$ or $S$, respectively.

The just described projections $P,P_1,P_2$ and hence $P_{12}=P_2\circ P_1$ are implemented in a Mathematica notebook whose image is shown below. In the notebook, we take a (pseudo)random matrix $x\in[0,1]^{n\times n}$ (with $n=3$) and see that $P_{12}x\notin S$; so, $P_{12}x$ cannot be the same the projection $Px$ of $x$ onto $S$. (In the notebook, we refer to the $n$-tuples $s\in[n]^n$ corresponding to matrices $X(s)\in S_1$ as *good*; clearly, $s\in[n]^n$ is good iff $s_{s_i}\ne i$ for all $i\in[n]$.)

As for your second question, "Is there any algorithm to solve this problem?", I can say the following:

Asking multiple questions is not encouraged on MathOverflow. So, please consider moving this question to another post.

It was shown above that the subset $S$ of the set $\{0,1\}^{n\times n}$ of cardinality $2^{n^2}$ can be conveniently indexed by a subset of the set $[n]^n$ of the much smaller cardinality $n^n$.

Therefore, the procedure described in the Mathematica notebook should be feasible for $n\le10$.

For $n>10$, some "greedy"-type algorithms can probably give reasonable approximations to the desired projection.

Here is an image of the mentioned Mathematica notebook:

1more comment