The chaos of flow turbulence causes most quantity-of-interests J to become significantly non-convex over time, so that gradients of J soon fail to guide numerical simulations toward a useful local optimum. This challenge is investigated, and a multi-step penalty method is introduced to circumvent it, for the common case in which a relatively deterministic, and hence controllable, component of the flow is masked by the chaos. To do this the simulation time is divided into intervals, with periods approximating the chaos time scales, and solution discontinuity ∆q are introduced between the intervals. Model systems show how this temporarily enlarges search scale by making a related JA more convex. The discontinuities are then increasingly penalized to approach the ∆q → 0 physical limit. The challenge and method are illustrated with a logistics map, the Lorenz system, and an advection augmented Kuramoto–Sivashinsky Equation, before application to turbulent Kolmogorov flow, for which it far outperforms standard adjoint-based gradient search. The utility of such an optimized chaotic solution is discussed