def prodmix(I,K,a,p,epsilon,LB):
    """prodmix:  robust production planning using soco
    Parameters:
        I - set of materials
        K - set of components
        a[i][k] -  coef. matrix
        p[i] - price of material i
        LB[k] - amount needed for k
    Returns a model, ready to be solved.
    """

    model = Model("robust product mix")

    x,rhs = {},{}
    for i in I:
        x[i] = model.addVar(vtype="C", name="x(%s)"%i)
    for k in K:
        rhs[k] = model.addVar(vtype="C", name="rhs(%s)"%k)

    model.addCons(quicksum(x[i] for i in I) == 1)
    for k in K:
        model.addCons(rhs[k] == -LB[k]+ quicksum(a[i,k]*x[i] for i in I) )
        model.addCons(quicksum(epsilon*epsilon*x[i]*x[i] for i in I) <= rhs[k]*rhs[k])

    model.setObjective(quicksum(p[i]*x[i] for i in I), "minimize")

    model.data = x,rhs
    return model