A coupled binary linear programming-differential evolution (BLP-DE) approach is proposed in this paper to optimize the design of water distribution systems (WDSs). Three stages are involved in the proposed BLP-DE optimization method. In the first stage, the WDS that is being optimized is decomposed into trees and the core using a graph algorithm. Binary linear programming (BLP) is then used to optimize the design of the trees during the second stage. In the third stage, a differential evolution (DE) algorithm is utilized to deal with the core design while incorporating the optimal solutions for the trees obtained in the second stage, thereby yielding near-optimal solutions for the original whole WDS. The proposed method takes advantage of both BLP and DE algorithms: BLP is capable of providing global optimal solution for the trees (no loops involved) with great efficiency, while a DE is able to efficiently generate good quality solutions for the core (loops involved) with a reduced search space compared to the original full network. Two benchmark WDS case studies and one real-world case study (with multiple demand loading cases) with a number of decision variables ranging from 21 to 96 are used to verify the effectiveness of the proposed BLP-DE optimization