•Dependence of harmonic spectrum of PV current on solar irradiance is considered.•Update of PV current harmonic spectrum is included in the optimization problem.•Optimization problem is solved by both probabilistic and deterministic methods.•Optimal PV penetration levels are obtained in unbalanced distribution systems.•PV and load uncertainties are considered in the optimization problem. Photovoltaic (PV) system is one of the most important technologies among the renewable energy resources. With the increasing penetration of PV systems in distribution networks, harmonic distortions also rise due to the inevitable effects of PV inverters. Therefore, determination of the optimal PV penetration level is a vital issue in terms of preventing power quality problems. In this paper, probabilistic and deterministic approaches are proposed to determine the optimal penetration levels of PV systems in unbalanced distorted distribution networks by taking into account the uncertainty of load profile and the intermittent characteristic of PV system output power due to changes in solar irradiance. The Interior Point (IP) method is applied to solve the optimization problem in conjunction with the Monte Carlo Simulation (MCS) and K – means clustering, respectively. The methodology is based on the dependence of the harmonic spectrum of PV injected current on solar irradiance. The harmonic power quality parameters are calculated by using loop frame of reference based three phase harmonic power flow method in the unbalanced distribution network, where the nonlinear loads and PV systems are interfaced. The allowable PV penetration levels are determined based on power quality parameters comprising of total harmonic voltage and individual harmonic voltage distortions, and RMS bus voltage limits.