In recent years, increased attention has been shifted to the design of cooling systems in the injection molding process, as it becomes clear that the cooling systems affect significantly both productivity and the part quality. In order to systematically improve the performance of a cooling system the mold designer may need a computer aided optimal design system for designing the injection mold cooling systems and determining the process conditions during the cooling stage. In this chapter, an efficient optimization procedure for this problem is proposed utilizing (i) the special boundary element analysis, (ii) the corresponding design sensitivity analysis using the direct differentiation approach, and (iii) the optimization algorithm. For this optimal design, an objective function is proposed to minimize a weighted combination of the cooling time and the temperature nonuniformity over the part surface. The former has to do with the warpage in the final part, while the latter is directly related to the overall productivity of the injection molding process. In this optimization program, various design variables are considered as follows: (i) (design variables related to processing conditions) the inlet coolant bulk temperature and inlet coolant volumetric flow rate of each cooling channel and (ii) (design variables related to mold cooling system design) the radius and location of each cooling channel. Each step of the proposed optimization procedure will be briefly explained below. First, in thermal analysis, mold heat transfer is considered as a cyclic-steady, three-dimensional conduction; heat transfer within the melt region is treated as a transient, one-dimensional conduction; heat exchange between the cooling channel surfaces and the coolant is considered to be steady; heat exchange between ambient air and mold exterior surfaces is also considered steady. Numerical implementation includes the application of a hybrid scheme consisting of a modified three-dimensional boundary element method for the mold region and a finite difference method with a variable mesh for the melt region. However, it was found that seemingly negligible inaccuracy in the thermal analysis result sometimes lead to a meaningless sensitivity analysis result. In this study, the thermal analysis system based on the above-mentioned modified boundary element method has been improved and rigorous treatments of boundary conditions appropriate for sensitivity analysis have been developed by considering the following issues: (i) numerical convergency, (ii) the series solution in part thermal analysis, (iii) the treatment of tip surface of line elements, (iv) the treatment of coolant, and (v) the treatment of mold exterior surface. Using an example, the importance of these issue is amply demonstrated. Next, the sensitivity analysis program developed in the present studies utilizes the implicit differentiation of the boundary integral equations and the boundary conditions presented in thermal analysis with respect to all design variables to obtain the sensitivity equations. A sample problem is solved to demonstrate the accuracy and the efficiency of the present sensitivity analysis formulation as well as to discuss the characteristics of each design variable. Finally, the CONMIN algorithm is applied for the optimization program with the help of the above thermal analysis and corresponding design sensitivity analysis. In this optimization program, the proper constraints imposed upon the design variables are considered to maintain design reality. The CONMIN algorithm employs the augmented Lagrangian multiplier method to deal with the equality constraints and the Davidon–Fletcher–Powell method for the unconstrained minimization during the successive unconstrained minimization procedure. Two sample problems were solved to demonstrate the efficiency and the usefulness of the objective function. The developed computer aided optimal design system would be very useful for injection mold designers in obtaining an optimal configuration of an injection mold cooling system in terms of radii and locations of cooling channels, as well as determining the optimal processing conditions of the cooling stage in terms of the inlet coolant bulk temperature and the inlet coolant volumetric flow rate of each cooling channel by minimizing certain objective functions related to the part quality and/or the productivity in the injection molding processes.