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.