We consider $k$ square integrable random variables $Y_1,...,Y_k$ and $k$ random (row) vectors of length $p$, $X_1,...,X_k$ such that $X_i(l)$ is square integrable for $1\le i\le k$ and $1\le l\le p$. No assumptions whatsoever are made of any relationship between the $X_i$:s and $Y_i$:s. We shall refer to each pairing of $X_i$ and $Y_i$ as an environment. We form the square risk functions $R_i(β)=\mathbb{E}\left[(Y_i-βX_i)^2\right]$ for every environment and consider $m$ affine combinations of these $k$ risk functions. Next, we define a parameter space $Θ$ where we associate each point with a subset of the unique elements of the covariance matrix of $(X_i,Y_i)$ for an environment. Then we study estimation of the $\arg\min$-solution set of the maximum of a the $m$ affine combinations the of quadratic risk functions. We provide a constructive method for estimating the entire $\arg\min$-solution set which is consistent almost surely outside a zero set in $Θ^k$. This method is computationally expensive, since it involves solving polynomials of general degree. To overcome this, we define another approximate estimator that also provides a consistent estimation of the solution set based on the bisection method, which is computationally much more efficient. We apply the method to worst risk minimization in the setting of structural equation models.