For the green vehicle routing problem with simultaneous pickup and delivery, a mixed integer programming model is developed with the optimization objective of minimizing the distribution cost with carbon emission cost, and a three-stage Lagrangian heuristic algorithm combining mathematical programming methods and heuristics is proposed for solving the problem. In the first stage, the Lagrangian dual model of the problem is obtained using Lagrangian relaxation technique. In the second stage, an improved sub-gradient algorithm is designed to solve the dual model iteratively, and a repair mechanism is introduced to repair the solution corresponding to the lower bound obtained in each iteration to a high-quality feasible solution of the original problem. The feasible solution is then used to update the sub-gradient direction and step size in the next iteration. In the third stage, a heuristic local search algorithm is designed to optimize the feasible solution obtained in the second stage and further improve the quality of the solution to obtain the approximate optimal solution of the original problem. Experiments show that the proposed algorithm can obtain a high-quality solution to the problem and provide a compact lower bound to evaluate the quality of the solution.