Abstract: In order to imporved such disadvatages as the matrix, applied in the embedding method, requires large memory for storage, and as the global optimum solution can not be guaranteed in response matrix method, here an integral unsteady response matrix of the whole region, derived by the analytical FEM, is proposed as a constraint condition, embeddling into the linear programming model. Since the later can express the behavior of the drawdown in the whole domain, it can be found that the method, given in this paper, will facilitate optimal development and management of ground water resources.