Since you asked for general assistance and I have dealt with these type of calculations, I will give you some general points as the question evolved these last days:

• The expression for the activity of a species is strange. I imagine that a single species participates in various reaction, not one. Thus, there are many extents of reactions in the numerator. Let $i$ denote the chemical species for which we have $i = 1,2,...,M$ and $j$ the chemical reaction for which $j = 1,2,...,N$ if we have a total of $N$ reactions. Then, the molar fraction of species $i$ in the gas phase is \begin{equation} y_i = \frac{n_i}{n} = \frac{n_{j0} + \sum_{j = 1}^N \nu_{ij}\xi_{j}} {\sum_{i = 0}^M \left[n_{j0} + \sum_{j = 1}^N \nu_{ij}\xi_{j}\right]} \tag1 \end{equation} where the brackets are not really needed but are added for clarity.
• The gas composition $y_i$ is a variable that belongs to the gas phase. Thus, the total amount in de denominator in Eq. (1) must not include the amounts of any solid species.
• The activity for solid phases are one, but a little care must be taken once you solve the system and hopefully reach a solution. Once you have all the extents of reactions $\xi_1,\xi_2,...,\xi_N$, you calculate the final amounts of for the solid species. If any of those are zero or negative, then something is wrong in the calculation.
• The approach that you are taking, which is to minimise the total Gibbs energy is not correct. The function you are using is wrong. Even though it is possible to formulate this problem as that, I suggest you don't. Try to solve the $N$ equations which appear with the $N$ law of mass action \begin{equation} K_j = \prod a_{ij}^{\nu_{ij}} \quad j = 1,2,...,N \tag2 \end{equation} If you want to have the value of the total Gibbs energy, then this value is \begin{equation} G = \sum_{i = 1}^M n_i G_i = \sum_{i = 1}^M n_i [\Delta_f G_i + RT\ln a_i] \end{equation} this is the quantity which acquires a minimum value. It can be positive, or negative.

Since you asked for general assistance and I have dealt with these type of calculations, I will give you some general points as the question evolved these last days:

• The expression for the activity of a gaseous species is strange. I imagine that a single species participates in various reactions, not one. Thus, there should be several extents of reactions in the numerator. Let $i$ denote the chemical species for which we have $i = 1,2,...,M$, and $j$ the chemical reaction for which we have $j = 1,2,...,N$. Then, the molar fraction of species $i$ in the gas phase is \begin{equation} y_i = \frac{n_i}{n} = \frac{n_{i0} + \sum_{j = 1}^N \nu_{ij}\xi_{j}} {\sum_{i = 1}^M \left[n_{i0} + \sum_{j = 1}^N \nu_{ij}\xi_{j}\right]} \tag1 \end{equation} where the brackets are not really needed, but are added for clarity.
• The gas composition $y_i$ is a variable that belongs to the gas phase. Thus, the total amount in the denominator of Eq. (1) only includes the amounts of the species that 'live' in the gaseous phase.
• The activity for the solid species are one, but care must be taken once you solve the system and hopefully reach a solution. Once you have all the extents of reactions $\xi_1,\xi_2,...,\xi_N$, you calculate the final amounts for all the solid species. If any of those are zero or negative, then something is wrong in the calculation or the chosen solution is wrong. Remember that a set of non-linear equations yields many solutions, but only one in this case will be physically meaningful.
• The approach that you are taking which is to minimise the total Gibbs energy is not correct, and the equation you are minimising is wrong. Try to solve the $N$ equations which appear with the $N$ laws of mass action \begin{equation} K_j = \prod a_{ij}^{\nu_{ij}} \quad j = 1,2,...,N \tag2 \end{equation} If you want to have the value of the total Gibbs energy, then this value is \begin{equation} G = \sum_{i = 1}^M n_i G_i = \sum_{i = 1}^M n_i [\Delta_f G_i + RT\ln a_i] \tag3 \end{equation}
• The value of Eq. (3) is the minimum value and has units of $\ce{J}$. It can be positive or negative.
• If you insist on the approach of minimising $G$, i.e. Eq. (3), you can extend this method using the Lagrange multipliers. It can be proved that you will have to solve even more equations, a total of $M$ (number of species) and $w$ (number of material balances). The number $w$ coincides with the number of Lagrange multipliers $\lambda$'s which are not known. Search online and you will get many papers that tried this idea. However, keep it simple and just do it like in class like I said.