There is not, generally, a universal result for which you simply plug numbers into a formula. However, there is a good strategy that results in an upper bound for the achievable fidelity, with conditions on when that maximum fidelity can be achieved. I use it quite extensively in these two papers: paper1, paper2.
First, be aware that there are several different figures of merit that you could use. Your post seems to suggest that you're interested in the global fidelity, i.e. maximising the overlap of the output state $|\phi\rangle$ with the overall target state $|\psi\rangle|\psi\rangle$, if $|\psi\rangle$ was the state to be cloned. In that case, there is an optimal result, as stated by Craig Gidney in his answer (I think that result is only for "universal" cloning where the set of possible input states is any pure state of the correct dimension, with equal likelihood).
There is also the single-copy fidelity, where you want the density matrix of each individual qubit to have as large a fidelity as possible with a single copy of $|\psi\rangle$. From my awareness, this is the more commonly studied case, although that's probably just because the global fidelity case is entirely solved.
So, the strategy that I suggest is to calculate the matrix
$$
R=\sum_ip_i|\psi_i\rangle\langle\psi_i|^T\otimes |\psi_i\rangle\langle\psi_i|\otimes |\psi_i\rangle\langle\psi_i|.
$$
Let $\lambda$ be the maximum eigenvalue of $R$, and let $d$ be the dimension of the Hilbert space spanned by $S$. Then the achievable global fidelity is bounded by $F\leq d\lambda$.
If you want the single-copy fidelity, instead you'd use the matrix
$$
R=\frac{1}{2}\sum_ip_i|\psi_i\rangle\langle\psi_i|^T\otimes\left( |\psi_i\rangle\langle\psi_i|\otimes I+I\otimes |\psi_i\rangle\langle\psi_i|\right).
$$
In either case, you can decide if the bound can be achieved by looking at the corresponding eigenvector. If the maximum eigenvalue is unique and the eigenvector is maximally entangled across the partition of the first system and the other two, you can achieve the maximum. In the case where there is degeneracy, it is sufficient to find a mixture of those eigenstates such that the reduced density matrix on the first system is maximally mixed.
What I don't know about is what you can do if it turns out that the upper bound is not saturated - in all the standard cases, there's sufficient symmetry in the set of input states that the upper bound can be achieved.