This article presents constrained numerical optimization of high-order linearly and algebraically stable diagonally-implicit Runge–Kutta methods. After satisfying the desired order conditions, undetermined coefficients are optimized with respect to objective functions which consider accuracy, stability, and computational cost. Constraints are applied during the optimization to enforce stability properties, to ensure a well-conditioned method, and to limit the domain of the abscissa. Two promising third-order methods are derived using this approach, labelled SDIRK[3,(1,2,2)](3)L_14 and SDIRK[3,1](4)L_SA_5. Both optimized schemes have a good balance of properties. The relative error norm of the latter, the L2-norm scaled by a function of the number of implicit stages, is a factor of two smaller than comparable methods found in the literature. Variations on these methods are discussed relative to trade-offs in their accuracy and stability properties. A novel fifth-order scheme SDIRK[5,1](5)L_02 is derived with a significantly lower relative error norm than the comparable fifth-order A-stable reference method. In addition, the optimized scheme is L-stable. The accuracy and relative efficiency of the Runge–Kutta methods are verified through numerical simulation of van der Pol's equation, as well as numerical simulation of vortex shedding in the laminar wake of a circular cylinder, and in the turbulent wake of a NACA 0012 airfoil. These results demonstrate the value of numerical optimization for selecting undetermined coefficients in the construction of high-order Runge–Kutta methods with a balance between competing objectives.