Is there a techinque that allows us to approximate, with arbitrary precision, any positive real number $x$ by a sequence of ratios of $p$-smooth numbers, for given prime number $p>2$? If yes, is there a way to determine the "best" of such sequences?
In particular, if $(p_i)$ is the sequence of prime numbers and $p=p_N$ for some $N>1$, then I am looking for a sequence $(z_i)$ in $\mathbb{Z}^N$, with $z_i=(z_{i,1},\ldots,z_{i,N})$, such that given $x\in\mathbb{R}_{>0}$ and $\epsilon>0$, we can always find some $I\in\mathbb{N}$ such that $i>I$ implies $\lvert{x-q_i}\rvert<\epsilon$, where $$q_i=\prod_{j=1}^{N}p_j^{z_{i,j}}\;.$$ For $N=\infty$, each term $z_i$ of the "best" sequence is determined by the prime-factorisation of $q_i\in\mathbb{Q}_{>0}$, with $(q_i)$ the sequence of the convergents of the continued-fraction representation of $x$ (see https://en.wikipedia.org/wiki/Continued_fraction#Best_rational_approximations).
Here are some remarks and proposals, too lengthy for a comment.
By taking the $\log$, approximating $x_0 \in \mathbb{R}$ by a $p$-smooth ratio is equivalent to finding a $(z_i)_{i=1..N} \in \mathbb{Z}^N$ such that $|\log(x_0) - \sum_{i=1}^N z_i \log(p_i) | < \varepsilon$, for a given $\varepsilon > 0$.
A brute force technique is to list all $(z_i)_{i=1..N}$ by order of $\max_{i=1..N} |z_i|$. This will succeed in finding an approximation as good as we want, because the $\log(p_i)$ are linearly independent on $\mathbb{Q}$. It can be improved a little, e.g. by using polygonal boundaries on exponents, but that's still very inefficient. So we are looking for a technique with a lower complexity.
We can represent real numbers by hyperplanes in $\mathbb{R}^N$, where $x \in \mathbb{R}$ is represented by the hyperplane $\sum_{i=1}^N x_i \log(p_i) = \log(x)$, each prime number being on one of the axes. In this setting, the problem is equivalent to finding a point on the $\mathbb{Z}^N$ grid in $\mathbb{R}^N$, that is nearer than $\varepsilon$ from the hyperplane representing $x_0$.
This is typically an Integer Linear Programming (ILP) problem: search $(z_i) \in \mathbb{Z}^N$ that minimizes $\sum_{i=1}^N z_i \log(p_i)$, with constraints $\forall i$ $|z_i| < K$, and $\sum_{i=1}^N z_i \log(p_i) \ge x_0$. That is, best approximation from above, for exponents lower than K. (Each $z_i$ can be replaced by the difference of 2 positive integer variables, as is usual in ILP).
The general ILP problem is NP-complete. Of course here we do not deal with the general problem, but with some particular instance where coefficients are the $\log$ of the first $N$ prime numbers, and constraints are quite simple. So there could be some algebraic reason that makes the problem less complex than NP-complete.
And there could also be a complexity difference between the stated problem, and a slightly more general one where the $(p_i)_{i=1..N}$ are not the first $N$ prime numbers, but any set of $N$ prime numbers. In this second problem, the $(p_i)$ can (I believe) be chosen so that they represent any set of coefficients $(c_i)$ strictly positive $\in \mathbb{R}^N$ with a scale factor $\lambda$.
I.e. $\forall \varepsilon > 0, \forall (c_i) \in \mathbb{R}^N, c_i > 0, \exists \lambda > 0 \in \mathbb{R}, \exists (p_i)_{i=1..N}$ prime numbers, $\forall i \text{ }|c_i - \lambda p_i| < \varepsilon$. The idea would be then to solve an ILP problem with coefficients $(c_i)$ with this technique, which would prove that, the ILP problem being NP-complete, the problem with logarithms would also be NP-complete. But there is much left to prove: that the $(c_i)$ problem is NP-complete notwithstanding the simple constraints, that we can cope with $\lambda$ (which probably means we need $\lambda$ to grow slower than $1/\epsilon$)...
So if the question about $p$-smooth approximation is not only theorical but an efficient implementation is looked for, one domain to search is Integer Linear Programming and the usual algorithms: cutting planes, branch-and-bound, etc. Or use one of the dedicated solvers.
On the theoretical side, a nice question would be to see if there are some patterns in the way exponents of the $(p_i)$ behave.
I had no time to do numerical experiments, however here is a smallish, hand-computed example, where the exponents oscillate: approximating $5$ with prime numbers $2$ and $3$ gives $2 \times 3, 3^2/2, 2^4/3, 3^4/2^4$ (whose good approximation is one of the reasons for the splitting of octava into 12 semitones, incidentally), $2^{15}/3^8$, etc. The hyperplane here is simply a line in $\mathbb{R}^2$, so the approximation oscillates between its two extremities. With $N$ prime numbers, the hyperplane can be projected onto the $S^{N-2}$ sphere; it would be interesting to see how the best approximations behave on this sphere.