In the management of groundwater in coastal areas, saltwater intrusion associated with extensive groundwater pumping, is an important problem. The groundwater optimization model is an advanced method to study the aquifer and decide the optimal pumping rates or optimal well locations. Cheng and Park gave the analytical solutions to the optimization problems basing on Strack's analytical solution. However, the analytical solutions have some limitations of the property of aquifer, boundary conditions, and so on. A simulation-optimization numerical method presented in this study can deal with non-homogenous aquifers and various complex boundary conditions. This simulation-optimization model includes the sharp interface solution which solves the same governing equation with Strack's analytical solution, therefore, the freshwater head and saltwater thickness should be in the same conditions, that can lead to the comparable results in optimal pumping rates and optimal well locations for both of the solutions. It is noticed that the analytical solutions can only be applied on the infinite domain aquifer, while it is impossible to get a numerical model with infinite domain. To compare the numerical model with the analytical solutions, calculation of the equivalent boundary flux was planted into the numerical model so that the numerical model can have the same conditions in steady state with analytical solutions.