Optimum well spacing is an essential element for the economic development of shale gas reservoirs. We present an integrated assisted history matching (AHM) and embedded discrete fracture model (EDFM) workflow for well spacing optimization by considering multiple uncertainty realizations and economic analysis. This workflow is applied in shale gas reservoirs of the Sichuan Basin in China. Firstly, we applied the AHM to calibrate ten matrix and fracture uncertain parameters using a real shale-gas well, including matrix permeability, matrix porosity, three relative permeability parameters, fracture height, fracture half-length, fracture width, fracture conductivity, and fracture water saturation. There are 71 history matching solutions obtained to quantify their posterior distributions. Integrating these uncertainty realizations with five well spacing scenarios, which are 517 ft, 620 ft, 775 ft, 1030 ft, and 1550 ft, we generated 355 cases to perform production simulations using the EDFM method coupled with a reservoir simulator. Then, P10, P50, and P90 values of gas estimated ultimate recovery (EUR) for different well spacing scenarios were determined. Additionally, the degradation of EUR with and without well interference was analyzed. Next, we calculated the NPVs of all simulation cases and trained the K-Nearest Neighbors (KNN) proxy, which describes the relationship between the NPV and all uncertain matrix and fracture parameters and varying well spacing. After that, the KNN proxy was used to maximize the NPV under the current operation cost and natural gas price. Finally, the maximum NPV of 3 million USD with well spacing of 766 ft was determined.