Time-delays are frequently occurred in various fields of science and engineering, such as communication networks, chemical processes, transmission lines, biological systems, population growth and power systems. Up to now, many research works have been devoted to the theoretical aspects and the numerical treatments of time-delay systems. Optimal control of time-delay systems is one of the most challenging mathematical problems in control theory. It is well-known that except for simple cases, it is either extremely difficult or impossible to obtain a closed-form solution for the above mentioned problems. So far, several computational techniques have been employed to obtain an approximate solution for delayed optimal control problems. Most of the existing well-established numerical algorithms can efficiently solve optimal control problems whose solutions are infinitely smooth or well-behaved. Owing to the inherent behavior of time-delay systems, the analytical solutions of these systems are different functions over the distinct subintervals. In such situations, neither the continuous basis functions (CBFs) including Chebyshev polynomials nor piecewise constant basis functions (PCBFs) such as block-pulse functions taken alone would form an efficient basis in the representation of piecewise smooth functions. The Chebyshev polynomials have been successfully used for solving various problems whose solutions are infinitely smooth and well-behaved. An extremely important property of Chebyshev polynomials is the good representation of smooth functions provided that the desired function is infinitely differentiable. The coefficients in Chebyshev expansion approach zero faster than any power of 1/m as m goes to infinity. Although, the Chebyshev polynomials enjoy some excellent properties, they also suffer from many drawbacks. For instance, they cannot properly model the inherent behavior of piecewise smooth optimal control problems. In additions, the convergence rate of the Chebyshev polynomials depends highly on the smoothness of the solution. So, we expect a slow convergence rate and the results obtained by Chebyshev polynomials are not in good agreement with the optimal solutions. Indeed, these drawbacks are due to the lack of smoothness in the exact response of optimal control. As mentioned before, the analytical solutions of optimal control problems whose solutions are piecewise smooth cannot be produced solely either by continuous basis functions or by piecewise constant basis functions. In fact, PCBFs are appropriate to deal with discontinuous functions, whereas CBFs are better for continuous functions. To overcome these drawbacks, it behoves us to seek a suitable hybrid basis inherently possessing the required features of the solutions corresponding to non-smooth optimal control problems. Combining Chebyshev polynomials and block-pulse functions constitutes a flexible framework to efficiently represent piecewise functions. In this thesis, a direct approach is successfully implemented for solving nonlinear optimal control problems involving multiple time delays. The method is based on a hybrid of block-pulse functions and Chebyshev polynomials using the well-known Chebyshev-Gauss-Lobatto points. The associated operational matrices of integration, derivative, product and delay are used to convert the delayed optimal control problem into a parameter optimization problem whose solution is much easier than the original one. Various types of multi-delay systems are investigated to show the efficiency and accuracy of the proposed approximation method.