This paper extends the recently developed hybrid method to find the optimal designs of systems with correlated non-gaussian random parameters. A double-bounded density function is used to approximate marginal distribution and a Frank copula is used to define dependence (a more general concept than correlation) among the random parameters. We use a Piecewise Ellipsoidal method to approximate the constraint region by a set of quadratic functions. The yield is estimated by a joint cumulative density function over a portion of the tolerance body contained in the feasible region. Yield maximization is done for positive and negative correlations and non-symmetrical marginal distributions, and tested on an example using Monte-Carlo simulation.