The aim of this study was to design and optimize a transdermal liposomes gel formulation for paeonol (PAE). A three-factor, three-level Box-Behnken design was used to derive a second-order polynomial equation to construct three-dimensional (3-D) contour plots for prediction of responses. Independent variables studied were the DC-Chol concentration (X-1), molar ratio of lipid/drug (X-2), and the polymer concentration (X-3), and the levels of each factor were low, medium, and high. The dependent variables studied were the encapsulation efficiency (%EE) of PAE (Y-1), flux of PAE (Y-2), and viscosity of the gels (Y-3). Response surface plots were drawn and statistical validity of the polynomials was established to find the compositions of optimized formulation, which was evaluated using the Franz diffusion cell. The %EE of PAE increased proportionally with the molar ratio of lipid/drug, but decreased with polymer concentration, whereas the flux of PAE increased proportionally with polymer concentration and the DC-Chol concentration. The viscosity of gels increased with the polymer concentration. Gels showed a non-Fickian diffusion release mechanism for PAE, and the in vitro release profiles were fit for Higuchi's order model. The design demonstrated the role of the derived polynomial equation and 3-D contour plots in predicting the values of dependent variables for the preparation and optimization of gel formulation for transdermal drug release.