We derive and analytically solve renormalization group (RG) equations of gauge invariant non-local Wilson line operators. This work aims to achieve a systematic resummation of subleading power infra-red effects in event shapes. The RG equations are solved using a novel approach that combines the use of gauge invariant non-local Wilson lines with the standard RG flow equations. The resulting solution is a resummed expression for the event shape distribution that includes subleading power corrections. This work has important implications for the study of event shapes in high-energy collisions, where subleading power effects can significantly impact the accuracy of theoretical predictions.