Harmonic distortion in distribution networks is becoming increasingly severe. However, existing studies on harmonic mitigation typically focus on harmonic suppression from individual devices, neglecting coordinated optimization among multiple devices. This results in limited effectiveness and insufficient economic performance of harmonic mitigation strategies. To address this issue, a multi-objective optimal harmonic power flow model for three-phase unbalanced distribution systems based on collaborative scheduling of photovoltaic (PV) systems and electric vehicles (EV) is proposed. The model reduces power losses and voltage unbalance while simultaneously lowering harmonic distortion. However, the introduction of harmonic variables increases the convergence difficulty of the non-convex model. To overcome this challenge, the original optimal harmonic power flow model is decomposed into an optimal power flow sub-model and a quadratic programming sub-model, and an alternating optimization approach is employed to achieve fast and efficient solution. Finally, the efficacy of the proposed methodology is verified using the IEEE 34-bus and IEEE 123-bus systems.