  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 $\Lambda$ and
universal constant $\delta$) 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. 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 $\delta$, on
the precessions of S-star orbits. We analyze S-star orbits assuming different
strength of Yukawa interaction $\delta$ and find that this parameter has strong
influence on range of Yukawa interaction $\Lambda$. For that purpose we use 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.
