The Variational Germano Identity [1, 2] is used to optimize the coefficients of residual-based subgrid-scale models that arise from the application of a Variational Multiscale Method [3, 4]. It is demonstrated that numerical iterative methods can be used to solve the Germano relations to obtain values for the parameters of subgrid-scale models that are nonlinear in their coefficients. Specifically, the Newton-Raphson method is employed. A least-squares minimization formulation of the Germano Identity is developed to resolve issues that occur when the residual is positive and negative over different regions of the domain. In this case a Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm is used to solve the minimization problem. The developed method is applied to the one-dimensional unsteady forced Burgers’ equation and the two-dimensional steady Stokes’ equations. It is shown that the Newton-Raphson method and BFGS algorithm generally solve, or minimize the residual of, the Germano relations in a relatively small number of iterations. The optimized subgridscale models are shown to outperform standard SGS models with respect to a L2 error. Additionally, the nonlinear SGS models tend to achieve lower L2 errors than the linear models.