This dissertation addresses optimization of cubic polynomial problems. Heuristics for finding good quality feasible solutions and for improving on existing feasible solutions for a complex industrial problem, involving cubic and pooling constraints among other complicating constraints, have been developed. The heuristics for finding feasible solutions are developed based on linear approximations to the original problem that enforce a subset of the original problem constraints while it tries to provide good approximations for the remaining constraints, obtaining in this way nearly feasible solutions. The performance of these heuristics has been tested by using industrial case studies that are of appropriate size, scale and structure. Furthermore, the quality of the solutions can be quantified by comparing the obtained feasible solutions against upper bounds on the value of the problem. In order to obtain these upper bounds we have extended efficient existing techniques for bilinear problems for this class of cubic polynomial problems. Despite the efficiency of the upper bound techniques good upper bounds for the industrial case problem could not be computed efficiently within a reasonable time limit (one hour). We have applied the same techniques to subproblems with the same structure but about one fifth of the size and in this case, on average, the gap between the obtained solutions and the computed upper bounds is about 3%. In the remaining part of the thesis we look at global optimization of cubic polynomial problems with non-negative bounded variables via branch and bound. A theoretical study on the properties of convex underestimators for non-linear terms which are quadratic in one of the variables and linear on the other variable is presented. A new underestimator is introduced for this class of terms. The final part of the thesis describes the numerical testing of the previously mentioned underestimators together with approximations obtained by considering lifted approximations of the convex hull of the (x x y) terms. Two sets of instances are generated for this test and the descriptions of the procedures to generate the instances are detailed here. By analyzing the numerical results we can conclude that our proposed underestimator has the best behavior in the family of instances where the only non-linear terms present are of the form (x x y). Problems originating from least squares are much harder to solve than the other class of problems. In this class of problems the efficiency of linear programming solvers plays a big role and on average the methods that use these solvers perform better than the others.