Soft robot designs based on auxetic (i.e., negative Poisson's ratio) lattices could offer superior dexterity, tunable local kinematics, and morphological intelligence. However, the design and control of these structures for robotic tasks, requiring multiple states and motions, remains a challenging problem. Finite element models (FEMs) offer a promising way of predicting robot behaviour that might be used for design and control optimization. Yet, these physics-based models often have high computational cost and can not provide explicit gradient information to guide the search for optimal designs. In this paper, we abstract the physical predictions of FEMs through differentiable surrogate models and demonstrate design and trajectory optimization using a gradient-based optimizer. We compare the performance of convolutional neural networks (CNNs) and graph neural networks (GNNs) as surrogate models. We then demonstrate the use of a gradient-based optimizer to find optimal designs for a specified deformation and optimal pairs of designs and actuation inputs for a trajectory specified by waypoints. In each case, the differentiable surrogate model enables the gradient-based optimizer to discover novel designs lying outside of the training data that achieve the required motions (with an relative error ≤10% for trajectories).