The paper considers the Carleman bilinearization procedure for developing weakly nonlinear models of electric power systems by representing the original nonlinear equations of the system in a bilinear form. In particular, we describe a general algorithm for representing a system of nonlinear differential equations in a bilinear form and proposed techniques for optimizing the dimension of the obtained bilinear system. An example of a single-machine infinite bus system demonstrates the bilinearization procedure and model order reduction.