To help resolve the issues of non-realizability and restriction to homogeneity faced by the analytical theories of turbulence, we explore three-dimensional homogeneous shear turbulence of an incompressible Newtonian fluid within the context of optimal control and convex optimization. The framework is composed of multi-point spatial correlations of velocity and pressure fluctuations up to the degenerate fourth order, their evolution equations, and the constraints. The integral of the trace of the second order correlations is argued as the objective functional to be maximized. The sources of the constraints are discussed, such as the Cauchy-Schwarz inequality and the nonnegativity of variance of products. Two models are introduced: the second-order model uses the contracted and degenerate third order correlations as control variables; the third-order model takes the degenerate fourth order correlations as control variables. Both model are second-order cone programs when discretized. The nature of large-scale and huge-scale computations and the link to big data are commented on. The asymptotic steady state of the second-order model is solved numerically. The predicted values of the anisotropy tensor are consistent with experimental data qualitatively (in regard to the relative numerical order pattern of the diagonal components), albeit with significant quantitative differences. Such differences are attributed to the non-enforceability of the non-negativity of variance of products within the model. The third-order model is expected to improve predictions, because of its ability to include constraints generated by this non-negativity requirement. The issue of how to solve this huge-scale problem is yet open.