Shale gas reservoir is characterized by complex pore structure, ultra-low permeability and a large number of natural fractures. Therefore, it is necessary to apply horizontal well drilling with hydraulic fracturing to stimulate the shale gas reservoirs by forming complex fractures network composed of main fractures (i.e., hydraulic fractures) and branch fractures (i.e., activated natural fractures). During shale gas fracturing, hydraulic fractures will initiate simultaneously, deflect under the stress interference with each other, activate adjacent natural fractures and connect with them. These complex behaviors of fractures propagation would greatly complicate the production prediction after shale gas fracturing. Currently, the majority of the production prediction models neglect the complexity of fractures network. Therefore, based on Fick's diffusion law, Langmuir isothermal adsorption equation and dual-medium seepage theory, this paper derived a productivity prediction model by considering complexity of fractures network, stress-sensitive effect and ad-desorption by using point source method, Duhamel principle and Laplace transform, and solved the model with perturbation theory, discrete superposition and Stehfest numerical inversion. Then, the production prediction model was verified by a field example analysis, and the influence of the morphology of main fractures and branch fractures on productivity is explored. It indicates that the main fracture shape has a great influence on the production, and the trapezoidal mode is closer to the production practice. The large deflection of main fractures, the high conductivity and large approaching angle of branch fractures are beneficial to improving the production of fractured shale gas horizontal wells.