A key problem in making precise perturbative QCD predictions is to set the proper renormalization scale of the running coupling. The extended renormalization group equations, which express the invariance of physical observables under both the renormalization scale- and scheme-parameter transformations, provide a convenient way for estimating the scale- and scheme-dependence of the physical process. In this paper, we present a solution for the scale-equation of the extended renormalization group equations at the four-loop level. Using the principle of maximum conformality (PMC) / Brodsky-Lepage-Mackenzie (BLM) scale-setting method, all non-conformal $\{\beta_i\}$ terms in the perturbative expansion series can be summed into the running coupling, and the resulting scale-fixed predictions are independent of the renormalization scheme. Different schemes lead to different effective PMC/BLM scales, but the final results are scheme independent. Conversely, from the requirement of scheme independence, one not only can obtain scheme-independent commensurate scale relations among different observables, but also determine the scale displacements among the PMC/BLM scales which are derived under different schemes. In principle, the PMC/BLM scales can be fixed order-by-order, and as a useful reference, we present a systematic and scheme-independent procedure for setting PMC/BLM scales up to NNLO. An explicit application for determining the scale setting of $R_{e^{+}e^-}(Q)$ up to four loops is presented. By using the world average $\alpha^{\bar{MS}}_s(M_Z) =0.1184 \pm 0.0007$, we obtain the asymptotic scale for the 't Hooft associated with the $\bar{MS}$ scheme, $\Lambda^{'tH}_{\bar{MS}}= 245^{+9}_{-10}$ MeV, and the asymptotic scale for the conventional $\bar{MS}$ scheme, $\Lambda_{\bar{MS}}= 213^{+19}_{-8}$ MeV. Show Partial Abstract