We propose new mathematical optimization models for generating sparse dynamical graphs, or networks, that can achieve synchronization. The synchronization phenomenon is studied using the Kuramoto model, defined in terms of the adjacency matrix of the graph and the coupling strength of the network, modelling the so-called coupled oscillators. Besides sparsity, we aim to obtain graphs which have good connectivity properties, resulting in small coupling strength for synchronization. We formulate three mathematical optimization models for this purpose. Our first model is a mixed integer optimization problem, subject to ODE constraints, reminiscent of an optimal control problem. As expected, this problem is computationally very challenging, if not impossible, to solve, not only because it involves binary variables but also some of its variables are functions. The second model is a continuous relaxation of the first one, and the third is a discretization of the second, which is computationally tractable by employing standard optimization software. We design dynamical graphs that synchronize, by solving the relaxed problem and applying a practical algorithm for various graph sizes, with randomly generated intrinsic natural frequencies and initial phase variables. We test robustness of these graphs by carrying out numerical simulations with random data and constructing the expected value of the network's order parameter and its variance under this random data, as a guide for assessment.