Shale gas production has achieved great success in the United States and Canada. It ’s role in reshaping the global energy landscape is becoming increasingly key. Multi-stage hydraulic fracturing in horizontal wells is widely accepted as the most economic and effective technique to unlock shale gas reservoirs. As shale gas production demand grows, efforts to improve the long-term productivity after hydraulic fracturing are needed.Eight published groups of gas and water relative permeability data were employed in this work. They come from Wood Ford, Green River, Wolf Camp and Eagle Ford shales. Descriptive statistics were utilized to understand the context of the gas and water relative permeability features of each data set used. The standard errors were within the range of 0.2-0.28 and the central tendencies concentrated in the range of 0.05-0.2. The methods of one-way ANOVA, two-way ANOVA and Factorial ANOVA were used to compare eight means of data sets. The results show that, all the groups of experimental data of gas and water relative permeability must be governed by a similar percolation law.Based on the abovementioned descriptive statistics, the relationship between gas and water relative permeability, shale’s absolute permeability, pore structure, lithology and gas properties was identified by deriving the mathematical equations. From these equations, the universal exponents in the equations were determined and then the mathematical model for gas and water relative permeability was developed. In comparing real and simulation data, the shale rocks approach Permeability Jail when the water saturations range from 0.557 to 0.663. This important phenomenon partially explains the rapid production decline and low water production rates in typical shale formations.By introducing the concept of dual-porosity to theories of mass conservation and volume equilibrium, a dynamic material balance equation was created that considers gas desorption and rock expansion. The formula for a pressure drop in two-phase (gas + water) flow were derived while considering slippage effect and friction loss along the horizontal path. The effective wellbore radius and fluid channeling from matrix to fracture were introduced and the differential equations of percolation with gas and water phases in a dual-porosity model were approximately solved while considering the interference between matrix blocks and fractures. The gas and water production over time in the drainage area of the stimulated reservoir region associated with a hydraulic fracture was obtained. This model was involved with heterogeneous fracture lengths and different production rates from each fracture.This research focused on the development of a novel simulator using mathematical models to anticipate the gas and water production from dual-porosity shale reservoirs over time, and optimize the post-fracturing deliverability with fewest fracturing stages. After investigating the workflow of shale gas production, the interface and functions of the novel simulator were designed. The interface was realized with QT and performed by coding with C++ language. The Production Optimization section was completed. The real data from two shale gas wells of the Shuangyang (Well 1) and Weiyuan (Well 2) reservoirs in China were utilized to validate the developed simulator. The flowing pressure, gas and water production rates, cumulative gas and water productions are output in terms of time steps. The real and simulation data fit each other well, demonstrating its good reliability. The validation results show that eight stages are the optimal for Well 1 and twelve stages are the best for Well 2. Finally, conclusions and recommendations for future work are provided.