A mathematical algorithm is presented for the computation of optimum doping profiles that maximize the breakdown voltage and on-state current in insulated gate bipolar transistors (IGBT). The algorithm is based on the evaluation of doping sensitivity functions, which are defined as the functional derivatives of the breakdown voltage and on-state current with respect to doping concentration. These functions are computed using an adjoint method and are used in combination with a gradient-based technique to search the optimization space of possible doping configurations efficiently. The mathematical algorithm is implemented numerically to optimize semiconductor devices that are simulated using finite element models and, then, applied to punch-through IGBTs with planar structure. In order to optimize the breakdown voltage it is shown that it is optimum to decrease the doping concentration in the drift region, particularly near the p-type junction on the emitter side and introduce p-type layers with low doping concentration in the drift region. In the case of the on-state current it is optimum to increase the n-type concentration the drift region, near the emitter junction. Depending on the initial structure and criteria imposed during the optimization, it is possible to increase the breakdown voltage by at least 5-10% and decrease the on-state voltage by at least 200 mV. The algorithm presented in this article can be easily extended to the optimization of three-dimensional doping profiles and to the optimization of other power devices, such that power p-n junctions and power metal-oxide-semiconductor field-effect-transistors.