Although a multi-stage hydraulically fractured horizontal well in a shale reservoir initially produces gas at a high production rate, this production rate declines rapidly within a short period and the cumulative gas production is only a small fraction (20-30%) of the estimated gas in place. In order to maximize the gas recovery rate (GRR), this study proposes a multi-parameter optimization model for a typical multi-stage hydraulically fractured shale gas horizontal well. This is achieved by combining the response surface methodology (RSM) for the optimization of objective function with a fully coupled hydro-mechanical FEC-DPM for forward computation. The objective function is constructed with seven uncertain parameters ranging from matrix to hydraulic fracture. These parameters are optimized to achieve the GRR maximization in short-term and long-term gas productions, respectively. The key influential factors among these parameters are identified. It is established that the gas recovery rate can be enhanced by 10% in the short-term production and by 60% in the long-term production if the optimized parameters are used. Therefore, combining hydraulic fracturing with an auxiliary method to enhance the gas diffusion in matrix may be an effective alternative method for the economic development of shale gas.