You need to solve a system of nonlinear equations, a.k.a. a multivariate root-finding problem: If $m_j$ are the sample moments and $\theta_1, ...,\theta_k$ the parameters, you're solving a set of equations of the form $f_j^*(\theta_1, ...,\theta_k)=0$ where in this case $f_j^*(\theta_1, ...,\theta_k)=f_j(\theta_1, ...\theta_k)-m_j$.
Is this the correct/standard approach in this setting?
Certainly one way to attempt to solve a system of nonlinear equations is to turn it into an optimization problem (as you have done in one particular way) and attempt to solve that, though it's not the only way to approach these problems.
It's an approach that should work reasonably well in a variety of situations; whether it works well in your specific case is not clear.
Note, however, that in general there may not necessarily be a unique solution to a moment problem, and even if there is a unique solution, it may not necessarily be "nice" enough to always yield good results with this particular approach (for example, the surface of your quadratic minimization problem may not be nice and unimodal in the parameter-space. It may instead have multiple local minima, or may have ridges or saddle points, and so on. One thing that's important to check with such an approach is that a minimum that you find is actually a solution (or very near to being a solution) to the original system of equations.
Likelihood-based methods are also infeasible, since there's no clean way to write the likelihood in terms of the parameters.
There's no need to write anything in any simple/"closed" form.
If I give you a set of parameters and an observation, can you evaluate the density/pmf?
If so, then it's potentially solvable. Some optimization algorithms don't require more than this. If you can evaluate derivatives, that may help but it's not necessarily required.