In this work we derive an exact formula to calculate the Expected Shortfall (ES) value for the one-factor delta-gamma approach which, to the best of our knowledge, was still missing in the literature. We then use the one-factor delta-gamma as a control variate to estimate the ES of the multi-factor delta-gamma approach. A one-factor delta-gamma approximation is used for each risk factor appearing in the problem. Since the expected values of control variates are computed by means of an exact formula, the additional effort of computation with respect to the naive estimator of the multi-factor delta-gamma can be neglected. With this method, we achieve a considerable reduction of the variance. We have established a theorem to prove that the variance is further reduced when we use all the risk factors instead of just some of them. We show that one of the main potential applications takes place in the insurance industry regulation within the Swiss solvency test framework. We perform a model risk analysis and illustrate these results with numerical experiments.