The vacuum energy is calculated for Yang-Mills (YM) system defined in $D$ dimensional space-time of $S^1\times R^d$ ($D=d+1$), where the possibility of the YM fields to acquire the vacuum expectation values on $S^1$ is taken into account. The vacuum energy has already been obtained to the order of one-loop in many people. Here we calculate the vacuum energy in $D$ dimensions to two-loop order. With an intention to reach higher loops, an approximation method is proposed, which is especially effective in higher dimensions. By this method, we can treat the higher-loop contributions of YM interactions as easily as we treat one-loop effect. As a check, we show reproduction of the two-loop contribution ($D$-dependence of the coefficient as well as the functional form) when the coupling constant is small. This approximation method is useful not only for the Kaluza-Klein theories but also for the finite temperature-density system (as a quark-gluon plasma).