We developed via ML an analytical potential that can accurately capture many-body interactions between polymer-grafted NPs in a polymer matrix and used the potential to explore NP assembly over large length and time scales. Our approach relies on computing relevant PMFs from CG MD simulations to extract effective two-and three-body interactions between NPs, which are then fitted to permutationally invariant polynomials. As the polymer is explicitly treated in the CG model, these fitted interactions implicitly capture all energetic and entropic effects of the grafted and matrix chains. We validated the many-body potential by comparing the assembly behavior of NPs treated using the fitted potential against those treated using the CG model and demonstrating that the potential can not only accurately reproduce the assembled structures but also capture the assembly pathways. Given that this potential eliminates the need to simulate the numerous degrees of freedom associated with the grafted polymer and the polymer matrix, the potential can speed up NP simulations by at least three orders of magnitude. By simulating the assembly behavior of large systems of NPs over long timescales, we could successfully capture the formation of the experimentally observed 1D string and 2D hexagonal sheet phases and illuminate the critical role of three-body interactions in the formation of these phases. We also elucidated the role of NP interactions on the formation of 3D globular aggregate and dispersed phases and showed that three-body interactions do not play a significant role in their formation. The manybody potentials developed here for specific polymer-grafted NPs also allowed us to rapidly explore assembly behavior under other conditions, e.g., NP core-core interaction strengths and NP volume fractions, without the need for generating new training data. In this manner, we were able to discover several other interesting NP phases, including string networks, small clusters, and two different varieties of gels, where one is a dense network of short 1D strings and the other is a network of partial-2D and 3D motifs. In addition, we showed the utility of Steinhardt's bond orientational order parameters in effective classification of NP assemblies.