Solving under-determined system of equations and constraints with Sympy and Scipy

Question

I'm working on a python script that will help with annual portfolio rebalancing. For those of you unfamiliar with the concept, it basically means annually buying/selling assets that have gained/lost value so that your portfolio composition actually matches what your desired asset allocation is. If on January 1st, you have a portfolio with a 50/50 stock/bond split, and the stock market has a great year, by Dec 31st of that same year you'll be overweight stocks in your portfolio, which means you'll need to sell some stocks to buy bonds in order to get that split back to 50/50. Figuring out how to do this is trivially easy if all the assets are in one account, but it can get complicated if you have investments in multiple accounts, especially if they're retirement accounts that you're not supposed to withdraw from until retirement age.

For the rebalancing step, I set up a system of 8 simple linear equations. They look roughly like this (see augmented matrix in 2nd code block):

Asset1 + Asset2 + ... = us_bonds_target 
Asset1 + Asset2+ ... = total_401k_value
...and so on

I then use Sympy's solve_linear_system function in the following manner:

from sympy import Matrix, symbols, solve_linear_system

Assets_Matrix = Matrix([[1, 0, 0, 0, 1, 0, 0, 1, 9571165],
                        [0, 1, 0, 0, 0, 0, 0, 0, 7298011],
                        [0, 0, 0, 1, 0, 0, 0, 0, 4665941],
                        [0, 0, 1, 0, 0, 1, 0, 0, 7178371],
                        [0, 0, 0, 0, 0, 0, 1, 0, 7178371],
                        [1, 1, 1, 1, 0, 0, 0, 0, 22550494],
                        [0, 0, 0, 0, 1, 0, 0, 0, 7200311],
                        [0, 0, 0, 0, 0, 1, 1, 1, 6141054]])

Asset1, Asset2, Asset3, Asset4, Asset5, Asset6, Asset7, Asset8 = symbols('Asset1, Asset2, Asset3, Asset4, Asset5, Asset6, Asset7, Asset8', nonnegative = True)
solution = solve_linear_system(Assets_Matrix, Asset1, Asset2, Asset3, Asset4, Asset5, Asset6, Asset7, Asset8)

Running this code produces output like this:

{Asset5: 7200311,
 Asset6: -Asset8 - 1037317,
 Asset7: 7178371,
 Asset3: Asset8 + 8215688,
 Asset4: 4665941,
 Asset2: 7298011,
 Asset1: -Asset8 + 2370854}

This is pretty close to what I want, which is a list of assets and their target values for rebalancing. However, there are a few limitations:

1. Sometimes it will give me negative solutions, which in this context isn't useful. The asset values can't be less than 0 in real life.

2. This is a bigger issue: I can't specify constraints on any of the variables. Most of these assets are mutual funds, and it's advantageous to maintain certain minimum values for each asset, because some mutual funds have different "classes" that have lower expense ratios (how much it costs to hold the investment annually) if you meet an investment minimum. For most of these assets, I'd like to have no less than $10k invested in each. I realize that sometimes this will cause the systemto be unsolvable, but I'd at least like to attempt solving it with these constraints first, then relaxing them if the solver fails.

After poking around on stack overflow and google, I learned that this solve of problem should be solvable using linear programming. So I set up the problem as follows (note that I haven't incorporated minimum values requirements for assets into the code yet--other than them needing to be positive--I'm just trying to prove this approach will produce useful solutions):

from scipy.optimize import linprog

A = [[1, 0, 0, 0, 1, 0, 0, 1],
     [0, 1, 0, 0, 0, 0, 0, 0],
     [0, 0, 0, 1, 0, 0, 0, 0],
     [0, 0, 1, 0, 0, 1, 0, 0],
     [0, 0, 0, 0, 0, 0, 1, 0],
     [1, 1, 1, 1, 0, 0, 0, 0],
     [0, 0, 0, 0, 1, 0, 0, 0],
     [0, 0, 0, 0, 0, 1, 1, 1]]
B = [9571165, 7298011, 4665941, 7178371, 7178371, 22550494, 7200311, 6141054]
C = [0, 0, 0, 0, 0, 0, 0, 0]

linprog(c = C, A_eq = A, b_eq = B, bounds = (0, None), method = 'interior-point')

However, when I run this code, I get the following output:

     con: array([  2370852.29007765,         0.        ,         0.        ,
         7178369.8786112 ,         0.        ,  10586541.29564287,
               0.        ,  -1037319.12695402])
     fun: 0.0
 message: 'The algorithm terminated successfully and determined that the problem is infeasible.'
     nit: 4
   slack: array([], dtype=float64)
  status: 2
 success: False
       x: array([  3.97865052e-02,   7.29801100e+06,   6.64570623e-01,
         4.66594100e+06,   7.20031100e+06,   4.56818174e-01,
         7.17837100e+06,   1.67013585e+00])

It seems like linprog doesn't like my equations for some reason. Is my problem well-suited for the linprog function? And if so, what am I doing wrong? Or should I be approaching this problem in a different manner?

Answer 1

The algorithm terminated successfully and determined that the problem is infeasible.

This means that the problem is impossible to solve within the constraints. To understand why, let's take a closer look at the sympy results:

Asset5: 7200311
Asset6: -Asset8 - 1037317
Asset7: 7178371
Asset3: Asset8 + 8215688
Asset4: 4665941
Asset2: 7298011
Asset1: -Asset8 + 2370854

• Assets 2, 4, 5, 7 are uniquely defined.
• Assets 1, 3, 6, 8 share one degree of freedom

This means we can choose any value for one of assets 1, 3, 6, or 8 and the others are derived from that. Now apply the constraint that assets can't be negative.

• From Asset1 >= 0 follows that Asset8 <= 2370854
• From Asset3 >= 0 follows that Asset8 >= -8215688
• From Asset6 >= 0 follows that Asset8 <= -1037317 ... this is incompatible with Asset8 >= 0.

In summary, the problem can't be solved if all assets must be positive.