In order to account for big uncertainties such as well interferences, hydraulic and natural fractures’ properties and matrix properties in shale gas reservoirs, it is paramount to develop a robust and efficient approach for well spacing optimization. In this study, a novel well spacing optimization workflow is proposed and applied to a real shale gas reservoir with two-phase flow, incorporating the systematic analysis of uncertainty reservoir and fracture parameters. One hundred combinations of these uncertainties, considering their interactions, were gathered from assisted history matching solutions, which were calibrated by the actual field production history from the well in the Sichuan Basin. These combinations were used as direct input to the well spacing optimization workflow, and five “wells per section” spacing scenarios were considered, with spacing ranging from 157 m (517 ft) to 472 m (1550 ft). An embedded discrete fracture model was used to efficiently model both hydraulic fractures and complex natural fractures non-intrusively, along with a commercial compositional reservoir simulator. Economic analysis after production simulation was then carried out, by collecting cumulative gas and water production after 20 years. The net present value (NPV) distributions of the different well spacing scenarios were calculated and presented as box-plots with a NPV ranging from 15 to 35 million dollars. It was found that the well spacing that maximizes the project NPV for this study is 236 m (775 ft), with the project NPV ranging from 15 to 35 million dollars and a 50th percentile (P50) value of 25.9 million dollars. In addition, spacings of 189 m (620 ft) and 315 m (1033 ft) can also produce substantial project profits, but are relatively less satisfactory than the 236 m (775 ft) case when comparing the P25, P50 and P75 values. The results obtained from this study provide key insights into the field pilot design of well spacing in shale gas reservoirs with complex natural fractures.