Structural gradients in metallic materials can give rise to substantial extra strengths compared to their non-gradient counterparts. This strengthening effect originates from the plastic strain gradients arising in plastically deformed gradient structures. Here we develop a gradient theory of plasticity by incorporating the strengthening effect of plastic strain gradients into the classical J₂ flow theory. This gradient theory is numerically implemented to study the strain gradient plasticity in gradient nanotwinned (GNT) Cu with different twin thickness gradients and accordingly different built-in strength gradients. Numerical simulations reveal the primary features of gradient plasticity in GNT Cu under uniaxial tension, including progressive yielding, gradient distributions of plastic strain and extra plastic flow resistance. The combined simulation and experimental results show that the extra strength of GNT Cu depends on the square root of the built-in strength gradient and also on the square root of the saturated plastic strain gradient in fully yielded GNT samples. This finding enables the predictions of optimal gradient structures and associated gradient strength distributions, which suggest possible routes for achieving the maximum strength of GNT Cu.