Riemann explicit formula for $\pi(x)$ and its evaluation

There is the easiest way to obtain convergent explicit formulas :

Since $$\psi(x) =\sum_{p^k \le x} \log p=( x-\sum_\rho \frac{x^\rho}{\rho}-\log 2\pi)1_{x > 2}$$ (sum over trivial and non-trivial zeros)

Letting $\Pi(x) = \sum_{p^k \le x} \frac{1}{k}= \int_{2-\epsilon}^x \frac{\psi'(x)}{\log x}dx$ and $1_{x > 2}\int_2^x ( \frac{d}{dt}\frac{t^\rho}{\rho}) \frac{dx}{\log x} = 1_{x > 2}\int_2^x \frac{\rho t^{\rho-1}}{\rho\log t} dt=1_{x > 2}\int_{2^\rho}^{x^\rho} \frac{u}{\log u} du= (\text{li}(x^\rho)-\text{li}(2^\rho))1_{x > 2}$

then $$\Pi(x) = (\text{li}(x)-\text{li}(2))1_{x > 2}-\sum_\rho (\text{li}(x^\rho)-\text{li}(2^\rho))1_{x > 2}$$

Now let $\pi(x) =\sum_{n=1}^\infty \frac{\mu(n)}{n} \Pi(x^{1/n})=\sum_{n=1}^{\log_2(x)} \frac{\mu(n)}{n} \Pi(x^{1/n})1_{x^{1/n}> 2}$ and $\text{R}_\rho(x) = \sum_{n=1}^\infty \frac{\mu(n)}{n} (\text{li}(x^{\rho/n})-\text{li}(2^{\rho/n}))1_{x^{1/n} > 2}$ (those are finite sums)

therefore $$\pi(x) = \text{R}_1(x)-\sum_\rho \text{R}_\rho(x)$$

There are indeed many traps while evaluating $\pi(x)$ using the $\zeta$ zeros
(I will neglect here the smooth and quickly decreasing with $x$ contribution of the trivial zeros).
This answer detailed most of these traps but let's examine your specific points :

Now, whether I use $$\tag{3.1}R(x^\rho)=1+\sum _{n=1}^{\infty } \frac{\big(\log(x^\rho)\big)^n}{nn! \zeta (n+1)}$$ or $$\tag{3.2}R(x^\rho)=1+\sum _{n=1}^{\infty } \frac{\big(\rho \log(x)\big)^n}{nn! \zeta (n+1)}$$ in $(1)$, in both cases I do not get $\pi^*(x)$. (Series seems to diverge/oscillate.)

The known non trivial zeros $\rho$ verify the RH and may thus be written $\,\displaystyle\rho:=\frac 12+it$.
Now $(3.1)$ can't indeed provide the right answer because the numerical evaluation of $x^\rho$, before applying the logarithm, will give you for $x$ real $>1\quad\displaystyle x^\rho=\exp\left(\left(\frac 12+it\right)\log(x)\right)$.
The problem is that the evaluation of the logarithm will not return $\,\left(\frac 12+it\right)\log(x)\;$ but something like $\left(\frac 12+it\right)\log(x)+2\pi ki\;$ with $k\in\mathbb{Z}$ depending of $t$ (the evaluation of the logarithm imposing a phase in $(-\pi,\pi]\;$)

Using a different phase reference (branch) for the terms of the series isn't a fine idea but what about the Gram series as used in $(3.2)$?
This formula is indeed correct and was the one I used to produce this animation for $x\in(1,100)$. But this formula too has to be handled with care because we have to evaluate something behaving like (since $\,\zeta(n)\sim 1\,$ as $\,n\to +\infty$) :$$S(z):=\sum_{n=1}^{\infty } \frac{z^n}{nn!},\quad z:=\left(\frac 12+it\right)\log(x)$$ For the first non trivial zero and (say) $\,x=100\,$ we gave $|z|>65$, for the tenth $|z|>229$, for the hundredth $|z|>1089$.
The problem is that the terms of $S(z)$ and thus $R(x^\rho)\,$ will become very large (about $\dfrac {e^z}{z\sqrt{2 \pi z}}$ near $n=z$) before decreasing again to get under unity after about $\,\lceil e\,z\rceil\,$ terms. In practice for $\,x=100\,$ and the $100$-th non trivial zero you will need a working precision of $480$ digits!

Another idea proposed by reuns is to use the Möbius $\mu$ function formula (no high precision here) :

$$\tag{1}R(x)=\sum_{n=1}^{\lceil\log_2(x)\rceil} \frac{\mu(n)}n \operatorname{li}\bigl(x^{1/n}\bigr)$$ (since $\ \displaystyle\pi^{*}(x)=\sum_{n=1}^{\infty} \frac{\mu(n)}n \Pi^*\bigl(x^{1/n}\bigr)\;$ with $\Pi^*(z)=0\;$ for $\,z<2\;$ we need only the $\,\lceil\log_2(x)\rceil\;$ first terms for all the $R(x^\rho)$ expressions ... not just $R(x)$)

But from the logarithm in the $\operatorname{li}$ function we get exactly the same problem of "phase reduction" as previously. Fortunately a simply remedy exists from the definition of the exponential integral $\operatorname{Ei}$ by using $\,\operatorname{li}\left(e^x\right)=\operatorname{Ei}(x)$ and replacing $(1)$ with : $$\tag{2}R(x)=\sum_{n=1}^{\lceil\log_2(x)\rceil} \frac{\mu(n)}n \operatorname{Ei}\left(\frac 1n\log(x)\right)$$ and $$\tag{3}R\left(x^\rho\right)=\sum_{n=1}^{\lceil\log_2(x)\rceil} \frac{\mu(n)}n \operatorname{Ei}\left(\frac {\rho}n\log(x)\right)$$

The result for $\;\pi(x),\;x\in(2,100)$ (starting with $R(x)$ evaluated using the Gram series, subtracting the exact $\;\displaystyle\frac 1{\log(x)}-\frac 1{\pi}\arctan\left(\frac{\pi}{\log(x)}\right)\;$ contribution for the trivial zeros and two times the real part of the sum of the $\,200\,$ first non trivial zeros $\,R\left(x^\rho\right)\,$ terms using $(3)\;$) should be :

Concerning the use of the $\zeta$ zeros to evaluate the prime counting function $\pi(x)$ as in your expression $(1)$ this answer may he helpful. Shortly the combinatorial Meissel-Lehmer method allowed to compute $\pi(x)$ up to $10^{23}$ while the analytic method allowed to go up to $10^{25}$ (in $2016$, I don't know the state of the art). Note that more efficient expressions than $(1)$ were used as you may find in the history by Büthe and references.