We develop novel numerical methods for optimization problems subject to constraints given by nonlinear hyperbolic systems of conservation and balance laws in one space dimension. These types of control problems arise in a variety of applications, in which inverse problems for the corresponding initial value problems are to be solved. The optimization method can be seen as a block Gauss-Seidel iteration. The optimization requires one to numerically solve the hyperbolic system forward in time and the corresponding linear adjoint system backward in time. We test the optimization method on a number of control problems constrained by nonlinear hyperbolic systems of PDEs with both smooth and discontinuous prescribed terminal states. The theoretical foundation of the introduced scheme is provided in the case of scalar hyperbolic conservation laws on an unbounded domain with a strictly convex flux. In addition, we empirically demonstrate that using a higher-order temporal discretization helps to substantially improve both the efficiency and accuracy of the overall numerical method.