Closed form solutions for the flexural-torsional buckling of elastic beam-columns may only be obtained for simple end boundary conditions, and the case of uniform bending and compression. Moment gradient cases need approximate analytical or numerical methods to be used. Investigations presented in this paper deal with the analytical energy method applied for any asymmetric transverse loading case that produces a moment gradient. Part I of this paper is devoted entirely to the theoretical investigations into the energy based out-of-plane stability formulation and its general solution. For the convenience of calculations, the load and the resulting moment diagram are presented as a superposition of two components, namely the symmetric and antisymmetric ones. The basic form of a non-classical energy equation is developed. It appears to be a function dependent upon the products of the prebuckling displacements (knowfrom the prebuckling analysis) and the postbuckling deformation state components (unknowns enabling the formulation of the stability eigenproblem according to the linear buckling analysis). Firstly, the buckling state solution is sought by presenting the basic form of the non-classical energy equation in several variants being dependent upon the approximation of the major axis stress resultant M_�� and the buckling minor axis stress resultant M_z. The following are considered: the classical energy equation leading to the linear eigenproblem analysis (LEA), its variant leading to the quadratic eigenproblem analysis (QEA) and the other non-classical energy equation forms leading to nonlinear eigenproblem analyses (NEA). The novel forms are those for which the stability equation becomes dependent only upon the twist rotation and its derivatives. Such a refinement is allowed for by using the second order out-of-plane bending differential equation through which the minor axis curvature shape is directly related to the twist rotation shape. Secondly, the effect of coupling of the in-plane and out-of-plane buckling forms is taken into consideration by introducing approximate second order bending relationships. The accuracy of the classical energy method of solving FTB problems is expected to be improved for both H- and I-section beam-columns. The outcomes of research presented in this part are utilized in Part II.