Otimização SciPy com limites agrupados

Eu estou tentando executar uma otimização de portfólio que retorna os pesos que maximizam minha função de utilidade. Eu posso fazer essa parte muito bem, incluindo a restrição que pesos soma a um e que os pesos também me dão um risco de destino. Eu também incluí os limites para [0 <= pesos <= 1]. Este código é o seguinte:

def rebalance(PortValue, port_rets, risk_tgt):
    #convert continuously compounded returns to simple returns
    Rt = np.exp(port_rets) - 1 
    covar = Rt.cov()

    def fitness(W):
        port_Rt = np.dot(Rt, W)
        port_rt = np.log(1 + port_Rt)
        q95 = Series(port_rt).quantile(.05)
        cVaR = (port_rt[port_rt < q95] * sqrt(20)).mean() * PortValue
        mean_cVaR = (PortValue * (port_rt.mean() * 20)) / cVaR
        return -1 * mean_cVaR

    def solve_weights(W):
        import scipy.optimize as opt
        b_ = [(0.0, 1.0) for i in Rt.columns]
        c_ = ({'type':'eq', 'fun': lambda W: sum(W) - 1},
              {'type':'eq', 'fun': lambda W: sqrt(np.dot(W, np.dot(covar, W))\
                                                          * 252) - risk_tgt})
        optimized = opt.minimize(fitness, W, method='SLSQP', constraints=c_, bounds=b_)  

        if not optimized.success: 
           raise BaseException(optimized.message)
        return optimized.x  # Return optimized weights


    init_weights = Rt.ix[1].copy()
    init_weights.ix[:] = np.ones(len(Rt.columns)) / len(Rt.columns)

    return solve_weights(init_weights)

Agora eu posso me aprofundar no problema, eu tenho meus pesos armazenados em uma série de pandas MultIndex de tal forma que cada ativo é um ETF correspondente a uma classe de ativos. Quando um portfólio com pesos iguais é impresso, fica assim:

Saída [263]:
equity       CZA     0.045455
             IWM     0.045455
             SPY     0.045455
intl_equity  EWA     0.045455
             EWO     0.045455
             IEV     0.045455
bond         IEF     0.045455
             SHY     0.045455
             TLT     0.045455
intl_bond    BWX     0.045455
             BWZ     0.045455
             IGOV    0.045455
commodity    DBA     0.045455
             DBB     0.045455
             DBE     0.045455
pe           ARCC    0.045455
             BX      0.045455
             PSP     0.045455
hf           DXJ     0.045455
             SRV     0.045455
cash         BIL     0.045455
             GSY     0.045455
Name: 2009-05-15 00:00:00, dtype: float64

Como posso incluir um requisito de limites adicionais de forma que, quando eu agrupar esses dados juntos, a soma do peso caia entre os intervalos de alocação que eu predeterminei para essa classe de ativo?

Então, concretamente, quero incluir um limite adicional tal que

init_weights.groupby(level=0, axis=0).sum()
Saída [264]:
equity         0.136364
intl_equity    0.136364
bond           0.136364
intl_bond      0.136364
commodity      0.136364
pe             0.136364
hf             0.090909
cash           0.090909
dtype: float64

está dentro desses limites

[(.08,.51), (.05,.21), (.05,.41), (.05,.41), (.2,.66), (0,.16), (0,.76), (0,.11)]

[UPDATE] Eu imaginei que mostraria meu progresso com uma desajeitada solução de psuedo com a qual não estou muito feliz. Ou seja, porque não resolve os pesos usando o conjunto de dados inteiro, mas sim a classe de ativo por classe de ativo. A outra questão é que, em vez disso, ele retorna a série, em vez dos pesos, mas tenho certeza de que alguém mais apto do que eu poderia oferecer algum insight em relação à função groupby.

Então, com um ligeiro ajuste no meu código inicial, tenho:

PortValue = 100000
model = DataFrame(np.array([.08,.12,.05,.05,.65,0,0,.05]), index= port_idx, columns = ['strategic'])
model['tactical'] = [(.08,.51), (.05,.21),(.05,.41),(.05,.41), (.2,.66), (0,.16), (0,.76), (0,.11)]


def fitness(W, Rt):
    port_Rt = np.dot(Rt, W)
    port_rt = np.log(1 + port_Rt)
    q95 = Series(port_rt).quantile(.05)
    cVaR = (port_rt[port_rt < q95] * sqrt(20)).mean() * PortValue
    mean_cVaR = (PortValue * (port_rt.mean() * 20)) / cVaR
    return -1 * mean_cVaR  

def solve_weights(Rt, b_= None):
    import scipy.optimize as opt
    if b_ is None:
       b_ = [(0.0, 1.0) for i in Rt.columns]
    W = np.ones(len(Rt.columns))/len(Rt.columns)
    c_ = ({'type':'eq', 'fun': lambda W: sum(W) - 1})
    optimized = opt.minimize(fitness, W, args=[Rt], method='SLSQP', constraints=c_, bounds=b_)

    if not optimized.success: 
        raise ValueError(optimized.message)
    return optimized.x  # Return optimized weights

O seguinte one-liner retornará a série um pouco otimizada

port = np.dot(port_rets.groupby(level=0, axis=1).agg(lambda x: np.dot(x,solve_weights(x))),\ 
solve_weights(port_rets.groupby(level=0, axis=1).agg(lambda x: np.dot(x,solve_weights(x))), \
list(model['tactical'].values)))

Series(port, name='portfolio').cumsum().plot()

[Atualização 2]

As alterações a seguir retornarão os pesos restritos, embora ainda não sejam ideais, pois são divididos e otimizados nas classes de ativos constituintes, portanto, quando a restrição para o risco alvo é considerada apenas uma versão reduzida da matriz de covariância inicial, está disponível

def solve_weights(Rt, b_ = None):

    W = np.ones(len(Rt.columns)) / len(Rt.columns)
    if b_ is None:
        b_ = [(0.01, 1.0) for i in Rt.columns]
        c_ = ({'type':'eq', 'fun': lambda W: sum(W) - 1})
    else:
        covar = Rt.cov()
        c_ = ({'type':'eq', 'fun': lambda W: sum(W) - 1},
              {'type':'eq', 'fun': lambda W: sqrt(np.dot(W, np.dot(covar, W)) * 252) - risk_tgt})

    optimized = opt.minimize(fitness, W, args = [Rt], method='SLSQP', constraints=c_, bounds=b_)  

    if not optimized.success: 
        raise ValueError(optimized.message)

    return optimized.x  # Return optimized weights

class_cont = Rt.ix[0].copy()
class_cont.ix[:] = np.around(np.hstack(Rt.groupby(axis=1, level=0).apply(solve_weights).values),3)
scalars = class_cont.groupby(level=0).sum()
scalars.ix[:] = np.around(solve_weights((class_cont * port_rets).groupby(level=0, axis=1).sum(), list(model['tactical'].values)),3)

return class_cont.groupby(level=0).transform(lambda x: x * scalars[x.name])