Handling Cardinality Constraints in Bond Portfolio Optimization

# Bond portfolio optimization with Python reaching iteration limit


# Bond portfolio optimization with Python reaching iteration limit












I am using scipy to optimize a hypothetical bond portfolio for maximum yield by choosing from a list of bonds in the portfolio's investable universe while adhering to portfolio constraints such as minimum and maximum % of corporate bonds, duration constraints, minimum and maximum % in a single issuer, etc. I previously had no constraint for the total number of bonds in the optimized portfolio, and the code was obviously returning an optimized portfolio that was highly concentrated in a small number of bonds (say 20 to 30) with the highest yield. When I add a constraint for a range of the number of bonds the portfolio must be invested in (currently set as 70 securities as the lower bound and 150 securities as the upper bound, and for a security to be counted in this total number, it must have a weight of at least 0.5%) the code hits the maximum iteration limit and the optimization does not converge. The code for the lower and upper limit of the number of securities to be invested in is this line:

`def securities_count_constraint(weights): invested_securities_count = np.sum(weights > min_investment_weight) return [invested_securities_count - min_securities, max_securities - invested_securities_count]`

The specific values for this constraint is just pulled from an excel sheet where I set the various constraints.

Does anyone have any suggestions on what I should add/remove to fix this issue? Here is the code below:

```
import numpy as np
import pandas as pd
from scipy.optimize import minimize
from openpyxl import load_workbook

np.set_printoptions(suppress=True, precision=6)

file_path = ''
sheet_name = 'Optimizer'

df = pd.read_excel(file_path, sheet_name=sheet_name, header=None, engine='openpyxl')

df.dropna(how='all', inplace=True)

start_row = 4
start_col = 10
end_col = 22

last_row = df[10].last_valid_index()

bond_data = df.loc[start_row:last_row, start_col:end_col]

bond_data.columns = ['Cusip', 'Public or Private', 'Issuer', 'Coupon', 'Maturity', 'Rating',
                     'Industry', 'Country', 'Price', 'Yield', 'Duration', 'Duration Bucket','Current Weight']

#Loading constraints
bounds_df = pd.read_excel(file_path, sheet_name=sheet_name, usecols="G:H",header=None, engine='openpyxl')

#Extract bounds
min_securities = bounds_df.iloc[4,0]
max_securities = bounds_df.iloc[4,1]

Single_bond_lower_bound = bounds_df.iloc[5,0]
Single_bond_upper_bound = bounds_df.iloc[5,1]

min_investment_weight = bounds_df.iloc[6,0]

cash_lower_bound = bounds_df.iloc[8,0]
cash_upper_bound = bounds_df.iloc[8,1]

duration_lower_bound = bounds_df.iloc[9,0]
duration_upper_bound = bounds_df.iloc[9,1]

Government_lower_bound = bounds_df.iloc[10,0]
Government_upper_bound = bounds_df.iloc[10,1]

Credit_lower_bound = bounds_df.iloc[11,0]
Credit_upper_bound = bounds_df.iloc[11,1]

Us_lower_bound = bounds_df.iloc[12,0]
US_upper_bound = bounds_df.iloc[12,1]

#By rating
AAA_lower_bound = bounds_df.iloc[16,0]
AAA_upper_bound = bounds_df.iloc[16,1]

AA_lower_bound = bounds_df.iloc[17,0]
AA_upper_bound = bounds_df.iloc[17,1]

A_lower_bound = bounds_df.iloc[18,0]
A_upper_bound = bounds_df.iloc[18,1]

BBB_lower_bound = bounds_df.iloc[19,0]
BBB_upper_bound = bounds_df.iloc[19,1]

#By sector
Treasury_lower_bound = bounds_df.iloc[26,0]
Treasury_upper_bound = bounds_df.iloc[26,1]

Communication_lower_bound = bounds_df.iloc[30,0]
Communication_upper_bound = bounds_df.iloc[30,1]

Energy_lower_bound = bounds_df.iloc[31,0]
Energy_upper_bound = bounds_df.iloc[31,1]

Financial_lower_bound = bounds_df.iloc[32,0]
Financial_upper_bound = bounds_df.iloc[32,1]

Industrial_lower_bound = bounds_df.iloc[33,0]
Industrial_upper_bound = bounds_df.iloc[33,1]

Infrastructure_lower_bound = bounds_df.iloc[34,0]
Infrastructure_upper_bound = bounds_df.iloc[34,1]

Realestate_lower_bound = bounds_df.iloc[35,0]
Realestate_upper_bound = bounds_df.iloc[35,1]

#By duration buckets
A_bucket_lower_bound = bounds_df.iloc[42,0]
A_bucket_upper_bound = bounds_df.iloc[42,1]

B_bucket_lower_bound = bounds_df.iloc[43,0]
B_bucket_upper_bound = bounds_df.iloc[43,1]

C_bucket_lower_bound = bounds_df.iloc[44,0]
C_bucket_upper_bound = bounds_df.iloc[44,1]

D_bucket_lower_bound = bounds_df.iloc[45,0]
D_bucket_upper_bound = bounds_df.iloc[45,1]

E_bucket_lower_bound = bounds_df.iloc[46,0]
E_bucket_upper_bound = bounds_df.iloc[46,1]

F_bucket_lower_bound = bounds_df.iloc[47,0]
F_bucket_upper_bound = bounds_df.iloc[47,1]

G_bucket_lower_bound = bounds_df.iloc[48,0]
G_bucket_upper_bound = bounds_df.iloc[48,1]

H_bucket_lower_bound = bounds_df.iloc[49,0]
H_bucket_upper_bound = bounds_df.iloc[49,1]

I_bucket_lower_bound = bounds_df.iloc[50,0]
I_bucket_upper_bound = bounds_df.iloc[50,1]

J_bucket_lower_bound = bounds_df.iloc[51,0]
J_bucket_upper_bound = bounds_df.iloc[51,1]

#By issuer
issuer_df = pd.read_excel(file_path, sheet_name=sheet_name, usecols="B, G:H", skiprows=56, header=None, engine='openpyxl')
issuer_df.columns = ['Issuer', 'Lower Bound', 'Upper Bound']
issuer_df = issuer_df.dropna(subset=['Issuer'])

#Functions

yields = bond_data['Yield'].values
durations = bond_data['Duration'].values
num_bonds = len(bond_data)

initial_weights = bond_data['Current Weight'].values

print("Initial weights:")
print(initial_weights)
print("Sum of initial weights:", np.sum(initial_weights))

#Objective function to maximize
def portfolio_yield(weights,yields):
    return np.dot(weights,yields)

def neg_portfolio_yield(weights, yields):
    return -portfolio_yield(weights, yields)

#Constraints:
def create_bounds(bond_data):
    bounds = []
    for idx, row in bond_data.iterrows():
        if row['Industry'] in ['Cash', 'Treasury']:
            bounds.append((0,1))
        else:
            bounds.append((Single_bond_lower_bound, Single_bond_upper_bound))
    return bounds

bounds = create_bounds(bond_data)

def cash_constraint(weights):
    cash_weights = weights[bond_data['Industry'] == 'Cash']
    total_cash_weight = np.sum(cash_weights)
    return [total_cash_weight - cash_lower_bound, cash_upper_bound - total_cash_weight]

def duration_constraint(weights):
    weighted_avg_duration = np.dot(weights, durations)
    return [weighted_avg_duration - duration_lower_bound, duration_upper_bound - weighted_avg_duration]

weighted_avg_duration = np.dot(initial_weights, durations)
print(f"Initial Weighted Average Duration: {weighted_avg_duration}")

def government_constraint(weights):
    government_bond_weights = weights[bond_data['Industry'].str.contains('Federal|Treasury', regex=True)]
    total_government_bond_weight = np.sum(government_bond_weights)
    return [total_government_bond_weight - Government_lower_bound, Government_upper_bound - total_government_bond_weight]

def credit_constraint(weights):
    credit_weights = weights[~bond_data['Industry'].str.contains('Federal|Treasury|Cash', regex=True)]
    total_credit_bond_weights = np.sum(credit_weights)
    return [total_credit_bond_weights - Credit_lower_bound, Credit_upper_bound - total_credit_bond_weights]

def us_constraint(weights):
    us_weights = weights[bond_data['Country'] == 'US']
    total_us_weight = np.sum(us_weights)
    return [total_us_weight - Us_lower_bound, US_upper_bound - total_us_weight]

#Def by rating
def aaa_constraint(weights):
    aaa_weights = weights[bond_data['Rating'] == 'AAA']
    total_aaa_weight = np.sum(aaa_weights)
    return [total_aaa_weight - AAA_lower_bound, AAA_upper_bound - total_aaa_weight]

def aa_constraint(weights):
    aa_weights = weights[bond_data['Rating'] == 'AA']
    total_aa_weight = np.sum(aa_weights)
    return [total_aa_weight - AA_lower_bound, AA_upper_bound - total_aa_weight]

def a_constraint(weights):
    a_weights = weights[bond_data['Rating'] == 'A']
    total_a_weight = np.sum(a_weights)
    return [total_a_weight - A_lower_bound, A_upper_bound - total_a_weight]

def bbb_constraint(weights):
    bbb_weights = weights[bond_data['Rating'] == 'BBB']
    total_bbb_weight = np.sum(bbb_weights)
    return [total_bbb_weight - BBB_lower_bound, BBB_upper_bound - total_bbb_weight]

#Def by sector
def treasury_constraint(weights):
    treasury_weights = weights[bond_data['Industry'] == 'Treasury']
    total_treasury_weight = np.sum(treasury_weights)
    return [total_treasury_weight - Treasury_lower_bound, Treasury_upper_bound - total_treasury_weight]

def communication_constraint(weights):
    communication_weights = weights[bond_data['Industry'] == 'Communication']
    total_communication_weight = np.sum(communication_weights)
    return [total_communication_weight - Communication_lower_bound, Communication_upper_bound - total_communication_weight]

def energy_constraint(weights):
    energy_weights = weights[bond_data['Industry'] == 'Energy']
    total_energy_weight = np.sum(energy_weights)
    return [total_energy_weight - Energy_lower_bound, Energy_upper_bound - total_energy_weight]

def financial_constraint(weights):
    financial_weights = weights[bond_data['Industry'] == 'Financial']
    total_financial_weight = np.sum(financial_weights)
    return [total_financial_weight - Financial_lower_bound, Financial_upper_bound - total_financial_weight]

def industrial_constraint(weights):
    industrial_weights = weights[bond_data['Industry'] == 'Industrial']
    total_industrial_weight = np.sum(industrial_weights)
    return [total_industrial_weight - Industrial_lower_bound, Industrial_upper_bound - total_industrial_weight]

def infrastructure_constraint(weights):
    infrastructure_weights = weights[bond_data['Industry'] == 'Infrastructure']
    total_infrastructure_weight = np.sum(infrastructure_weights)
    return [total_infrastructure_weight - Infrastructure_lower_bound, Infrastructure_upper_bound - total_infrastructure_weight]

def realestate_constraint(weights):
    realestate_weights = weights[bond_data['Industry'] == 'Real Estate']
    total_realestate_weight = np.sum(realestate_weights)
    return [total_realestate_weight - Realestate_lower_bound, Realestate_upper_bound - total_realestate_weight]

def a_bucket_constraint(weights):
    a_bucket_weights = weights[bond_data['Duration Bucket'] == 'a']
    total_a_bucket_weight = np.sum(a_bucket_weights)
    return [total_a_bucket_weight - A_bucket_lower_bound, A_bucket_upper_bound - total_a_bucket_weight]

def b_bucket_constraint(weights):
    b_bucket_weights = weights[bond_data['Duration Bucket'] == 'b']
    total_b_bucket_weight = np.sum(b_bucket_weights)
    return [total_b_bucket_weight - B_bucket_lower_bound, B_bucket_upper_bound - total_b_bucket_weight]

def c_bucket_constraint(weights):
    c_bucket_weights = weights[bond_data['Duration Bucket'] == 'c']
    total_c_bucket_weight = np.sum(c_bucket_weights)
    return [total_c_bucket_weight - C_bucket_lower_bound, C_bucket_upper_bound - total_c_bucket_weight]

def d_bucket_constraint(weights):
    d_bucket_weights = weights[bond_data['Duration Bucket'] == 'd']
    total_d_bucket_weight = np.sum(d_bucket_weights)
    return [total_d_bucket_weight - D_bucket_lower_bound, D_bucket_upper_bound - total_d_bucket_weight]

def e_bucket_constraint(weights):
    e_bucket_weights = weights[bond_data['Duration Bucket'] == 'e']
    total_e_bucket_weight = np.sum(e_bucket_weights)
    return [total_e_bucket_weight - E_bucket_lower_bound, E_bucket_upper_bound - total_e_bucket_weight]

def f_bucket_constraint(weights):
    f_bucket_weights = weights[bond_data['Duration Bucket'] == 'f']
    total_f_bucket_weight = np.sum(f_bucket_weights)
    return [total_f_bucket_weight - F_bucket_lower_bound, F_bucket_upper_bound - total_f_bucket_weight]

def g_bucket_constraint(weights):
    g_bucket_weights = weights[bond_data['Duration Bucket'] == 'g']
    total_g_bucket_weight = np.sum(g_bucket_weights)
    return [total_g_bucket_weight - G_bucket_lower_bound, G_bucket_upper_bound - total_g_bucket_weight]

def h_bucket_constraint(weights):
    h_bucket_weights = weights[bond_data['Duration Bucket'] == 'h']
    total_h_bucket_weight = np.sum(h_bucket_weights)
    return [total_h_bucket_weight - H_bucket_lower_bound, H_bucket_upper_bound - total_h_bucket_weight]

def i_bucket_constraint(weights):
    i_bucket_weights = weights[bond_data['Duration Bucket'] == 'i']
    total_i_bucket_weight = np.sum(i_bucket_weights)
    return [total_i_bucket_weight - I_bucket_lower_bound, I_bucket_upper_bound - total_i_bucket_weight]

def j_bucket_constraint(weights):
    j_bucket_weights = weights[bond_data['Duration Bucket'] == 'j']
    total_j_bucket_weight = np.sum(j_bucket_weights)
    return [total_j_bucket_weight - J_bucket_lower_bound, J_bucket_upper_bound - total_j_bucket_weight]

#Issuer constraints
def issuer_constraints(weights):
    constraints = []
    for idx, row in issuer_df.iterrows():
        issuer = row['Issuer']
        #If no bounds are specified in columns G and H, use default 0% to 10%
        lower_bound = row['Lower Bound'] if not pd.isna(row['Lower Bound']) else 0.0
        upper_bound = row['Upper Bound'] if not pd.isna(row['Upper Bound']) else 0.10

        issuer_weights = weights[bond_data['Issuer'] == issuer]
        total_issuer_weight = np.sum(issuer_weights)
        #Add the constraints for the issuer weight
        constraints.append(total_issuer_weight - lower_bound)
        constraints.append(upper_bound - total_issuer_weight)
    return constraints

def securities_count_constraint(weights):
    invested_securities_count = np.sum(weights > min_investment_weight)
    return [invested_securities_count - min_securities, max_securities - invested_securities_count]

#Optimized weight constraints
def weight_sum_constraint(weights):
    return np.sum(weights) - 1

#Constraints dictionary
constraints = [
    {'type': 'eq', 'fun': weight_sum_constraint},
    {'type': 'ineq', 'fun': cash_constraint},
    {'type': 'ineq', 'fun': duration_constraint},
    {'type': 'ineq', 'fun': government_constraint},
    {'type': 'ineq', 'fun': credit_constraint},
    {'type': 'ineq', 'fun': us_constraint},
    {'type': 'ineq', 'fun': aaa_constraint},
    {'type': 'ineq', 'fun': aa_constraint},
    {'type': 'ineq', 'fun': a_constraint},
    {'type': 'ineq', 'fun': bbb_constraint},
    {'type': 'ineq', 'fun': treasury_constraint},
    {'type': 'ineq', 'fun': communication_constraint},
    {'type': 'ineq', 'fun': energy_constraint},
    {'type': 'ineq', 'fun': financial_constraint},
    {'type': 'ineq', 'fun': industrial_constraint},
    {'type': 'ineq', 'fun': infrastructure_constraint},
    {'type': 'ineq', 'fun': realestate_constraint},
    {'type': 'ineq', 'fun': a_bucket_constraint},
    {'type': 'ineq', 'fun': b_bucket_constraint},
    {'type': 'ineq', 'fun': c_bucket_constraint},
    {'type': 'ineq', 'fun': d_bucket_constraint},
    {'type': 'ineq', 'fun': e_bucket_constraint},
    {'type': 'ineq', 'fun': f_bucket_constraint},
    {'type': 'ineq', 'fun': g_bucket_constraint},
    {'type': 'ineq', 'fun': h_bucket_constraint},
    {'type': 'ineq', 'fun': i_bucket_constraint},
    {'type': 'ineq', 'fun': j_bucket_constraint},
    {'type': 'ineq', 'fun': issuer_constraints},
    {'type': 'ineq', 'fun': securities_count_constraint},

]

#Perform optimization
result = minimize(neg_portfolio_yield, initial_weights, args=(yields,),
    method='SLSQP', bounds=bounds, constraints=constraints, tol=1e-9,
    #Delete this if no longer works
    options={'maxiter': 1000}
)

if not result.success:
    print("Optimization did not converge. Check constraints or objective function.")
    print(result.message)

# Print optimization details
print("Optimization Result:")
print(result)

# Print optimized weights for easy copy-paste into Excel
print("Optimized Weights (Ready for Copy-Paste):")
for weight in result.x:
    print(f"{weight:.6f}")
```

## Answer by Attack68 (score 2, accepted)

https://quant.stackexchange.com/a/80701

This is a wild guess but this function is not like your others:

```
def securities_count_constraint(weights):
    invested_securities_count = np.sum(weights > min_investment_weight)
    return [invested_securities_count - min_securities, max_securities - invested_securities_count]
```

The partial derivative of this with respect to weights simply doesn't exist, making a gradient based optimiser essentially oscillate, and hence you get iteration limit breaches. SLSQP is not appropriate for all problems and this is certainly an element of the problem not suitable.

Since the number of bonds you are targeting seems to be fairly arbitrary anyway, what I would do is to try to add a penalty function to your objective function, which forces weights to be more distributed in more bonds than a few concentrated bonds.

Noting that:

$$ min_{\textbf{w}} \mathbf{w^Tw}, \quad s.t. \quad \mathbf{w^T1} = 1 \quad \implies \quad w_i = 1/n $$ $$ max_{\textbf{w}} \mathbf{w^Tw}, \quad s.t. \quad \mathbf{w^T1} = 1 \quad \implies \quad w_1 = 1, w_{i,i\ne1}=0 $$

provides a basis for a some parametrised penalty function: $\gamma \mathbf{w^Tw}$, which you must tune to arrive at an appropriate number of included bonds. Start with gamma very low and progressively increase it until you have a reasonable number of included bonds.