In this paper we investigate a Yukawa gravity modification of the Newtonian gravitational potential in a weak field approximation. For that purpose we derived the corresponding equations of motion and used them to perform two-body simulations of the stellar orbits. In 2020 the GRAVITY Collaboration detected the orbital precession of the S2 star around the supermassive black hole (SMBH) at the Galactic Center (GC) and showed that it is close to the general relativity (GR) prediction. Using this observational fact, we evaluated parameters of the Yukawa gravity (the range of Yukawa interaction Λ and universal constant δ) with the Schwarzschild precession of the S-stars assuming that the observed values as indicated by the GRAVITY Collaboration will have a small deviation from GR prediction [1]. GR provides the most natural way to fit observational data for S-star orbits, however, their precessions can be fitted by Yukawa gravity. Our main goal was to study the possible influence of the strength of Yukawa interaction, i.e. the universal constant δ, on the precessions of S-star orbits. We analyze S-star orbits assuming different strength of Yukawa interaction δ and find that this parameter has strong influence on range of Yukawa interaction Λ. For that purpose we use parameterized post-Newtonian (PPN) equations of motion in order to calculate the simulated orbits of S-stars in GR and Yukawa gravity. Using MCMC simulations we obtain the best-fit values and uncertainties of Yukawa gravity parameters for S-stars. Also, we introduce a new criterion which can be used for classification of gravitational systems in this type of gravity, according to their scales. We demonstrated that performed analysis of the observed S-stars orbits around the GC in the frame of the Yukawa gravity represent a tool for constraining the Yukawa gravity parameters and probing the predictions of gravity theories.