Shape optimization in unsteady flow problems enables the consideration of dynamic effects on design. The ability to treat unsteady effects is attractive, as it can provide performance gains when compared to steady-state design methods for a variety of applications in which time-varying flows are of paramount importance. This is the case, for example, in turbomachinery or rotorcraft design. Given the high computational cost involved in time-accurate design problems, adjoint-based shape optimization is a promising option. However, efficient sensitivity analysis should also be accompanied by a significant decrease in computational cost for the primal flow solution, as well. Reduced-order models, like those based on the harmonic balance concept, in combination with the calculation of gradients via adjoint methods, are proposed for the efficient solution of a certain class of aerodynamics optimization problems. The harmonic balance method is applicable if the flow is characterized by discrete finite dominant flow frequencies that do not need to be integer multiples of a fundamental harmonic. A fully-turbulent harmonic balance discrete adjoint formulation based on a duality-preserving approach is proposed. The method is implemented by leveraging algorithmic differentiation and is applied to two test cases: the constrained shape optimization of both a pitching airfoil and a turbine cascade. A key advantage of the current approach is the accurate computation of gradients as compared to second order finite differences without any approximation in the linearization of the turbulent viscosity. The shape optimization results show significant improvements for the selected time-dependent objective functions, demonstrating that design problems involving almost-periodic unsteady flows can be tackled with manageable computational effort.