Divisor must be a constant
AnsweredHello,
I've been trying to formulate a basic portfolio optimization problem in python where I am looking to maximize the expected portfolio return divided by the portfolio variance by varying the weights of the component stocks.
X is a list of the weights (decision variables A through E), V is the covariance matrix
num = simple sumproduct of X and the expected return of the corresponding stocks
denom = numpy.matmul(numpy.matmul(X,V),X)
When i set the objective to
Obj = num/denom
prob.setObjective(Obj, GRB.MAXIMIZE)
I get an error "Divisor must be a constant".
I have tried to reformulate this by setting up a constraint with a dummy variable and different variations of 1/denom to get around this but nothing seems to work.
Any ideas would be...much appreciated.
Thanks,
Marek

As the error suggests, Gurobi doesn't support dividing by variables. For more information, see What types of models can Gurobi solve?.
You could model this problem by introducing a few auxiliary variables and bilinear constraints. However, your problem looks very similar to problem of maximizing the Sharpe ratio (though you use variance instead of standard deviation), for which there is a nice convex reformulation. Let \( \mu \in \mathbb{R}^n \) be the vector of expected returns, and let \( Q \in \mathbb{S}_{+}^{n \times n} \) be the (symmetric, positive semidefinite) covariance matrix. I'll assume you are trying to maximize the Sharpe ratio \( \mu^\top x / \sqrt{x^\top Q x} \); I'm not aware of a common riskadjusted return measure that uses variance instead of standard deviation.
Your problem, which we refer to as \( {\bf(1)} \), can be written as
$$\begin{alignat*}{2} \max_x\ && \frac{\mu^\top x}{\sqrt{x^\top Q x}} & \\ \textrm{s.t.}\ && \sum_{i=1}^n x_i ={} &1 \\ && x \geq {} &0.\end{alignat*}$$
Consider the following problem, which we refer to as \( {\bf(2)} \):
$$\begin{alignat*}{2} \max_y\ && \frac{1}{\sqrt{y^\top Q y}} & \\ \textrm{s.t.}\ && \mu^\top y = {} &1 \\ && y \geq {} &0.\end{alignat*}$$
We'll show that \( {\bf(1)} \) and \( {\bf(2)} \) are "equivalent", in the sense that given a solution to either problem, we can construct a solution to the other of equal or better objective value. We assume there exists some \( \hat{x} \) for which \( \mu^\top \hat{x} > 0 \), which should hold in any reasonable portfolio optimization problem. It follows from this assumption that the optimal solution to \( {\bf(1)} \) has a positive expected return.
First, let \( \bar{x} \) be a solution to \( {\bf(1)} \). We construct a solution to \( {\bf(2)} \) of equal objective value. For all \( i = 1, \ldots, n \), define \( \bar{y}_i := \bar{x}_i / \mu^\top \bar{x} \). It holds that \( \bar{y} \geq 0 \), because \( \bar{x} \) is nonnegative and \( \mu^\top \bar{x} > 0 \) by our assumption. Additionally, \( \mu^\top \bar{y} = \mu^\top \bar{x} / \mu^\top \bar{x} = 1 \). Thus, \( \bar{y} \) satisfies the constraints of \( {\bf(2)} \). Finally, the objective of \( {\bf(2)} \) with respect to the solution \( \bar{y} \) is \( 1 / \sqrt{\bar{y}^\top Q \bar{y}} = 1 / \sqrt{ \bar{x}^\top Q \bar{x} / (\mu^\top \bar{x})^2} = \mu^\top \bar{x} / \sqrt{\bar{x}^\top Q \bar{x}} \).
Next, let \( \bar{y} \) be a solution to \( {\bf(2)} \). We construct a solution to \( {\bf(1)} \) of equal objective value. Let \( \bar{x}_i := \bar{y}_i / \sum_{j=1}^n \bar{y}_j \) for \( i = 1, \ldots, n \). Using the same ideas as above, it's easy to verify that \( \bar{x} \) satisfies the constraints of \( {\bf(1)} \) and has objective value equal to \( 1 / \sqrt{\bar{y}^\top Q \bar{y}} \).
So, instead of solving your original problem \( {\bf(1)} \), we can just solve \( {\bf(2)} \) and map the optimal solution \( \bar{y} \) back to the original problem using the translation \( \bar{x}_i := \bar{y}_i / \sum_{j=1}^n \bar{y}_j \) for \( i = 1, \ldots, n \). However, because \( Q \) is positive semidefinite, the solutions of \( {\bf(2)} \) are the same as the solutions to
$$\begin{alignat*}{2} \min_y\ && y^\top Q y \enspace & \\ \textrm{s.t.}\ && \mu^\top y = {} &1 \\ && y \geq {} &0.\end{alignat*}$$
This is a convex QP that Gurobi can solve directly. A solution \( \bar{y} \) to the above problem can be mapped to the original problem using the same transformation \( \bar{x}_i := \bar{y}_i / \sum_{j=1}^n \bar{y}_j \) for \( i = 1, \ldots, n \).
If you have other side constraints to \( {\bf(1)} \), see Erwin Kalvelagen's blog post on maximizing the Sharpe ratio. Basically, additional constraints \( Ax \leq b \) in \( {\bf(1)} \) translate to constraints \( Ay \leq b \kappa \) in \( {\bf(2)} \), where \( \kappa = \sum_{j=1}^n y_j \).
2 
My goodness Eli, this was a fantastic explanation. Thank you for taking the time to work through this and write this detailed of a reply. I am happy to say that I implemented this into my code without issue and got the exact result I was expecting. And yes, this is exactly the problem of maximizing the Sharpe ratio, but I had left out the square root due to other issues in an attempt to solve it in piecemeal fashion.
Thanks again!!
Marek
0 
Hi together, I have a similar problem, I understand Elis answer but am unable to implement the solution. How do I implement the described approach?
I have:
for tk in TK:
model.addConstr(R[tk]*A[tk]==E[tk])R, A and E are all variables in the model.
I would be glad if somebody could help here.
Kind regards
Moritz
0 
Which part of the implementation do you need help with? Your code looks okay, depending on how you defined \( \texttt{TK} \), \( \texttt{R} \), \( \texttt{A} \), and \( \texttt{E} \). For example, the following code runs without error:
import gurobipy as gp
TK = [1, 2, 3]
model = gp.Model()
R = model.addVars(TK, name='R')
A = model.addVars(TK, name='A')
E = model.addVars(TK, name='E')
for tk in TK:
model.addConstr(R[tk] * A[tk] == E[tk])If you encounter an error, it would be helpful to see a minimal, selfcontained code snippet that reproduces the error.
0 
Hi again,
below I send the snippet. It took a while to reduce the code to the problem, sorry for coming back late. I marked the line which makes trouble with a comment. Its about R[t]. Can you please support here? I am not aware of a a reasonable solution.
Thanks in advance and kind regards.
Moritz
from gurobipy import *
# =============================================================================
model = Model("Smart Factory Strategie")
model.modelSense = GRB.MAXIMIZE# =============================================================================
A=['1','2','3']
D={'1':1.5,'2':1.8,'3':1.5}
KWS={'1':[10842,1355,1355,1355,1355,1355,1355,1355,1355,1355,1355,1355,0,0,0,0],'2':[5421,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],'3':[10842,1355,1355,1355,1355,1355,1355,1355,1355,1355,1355,1355,0,0,0,0]}
KWE={'1':[35090,35090,35090,35090,35090,35090,35090,35090,35090,35090,35090,35090,0,0,0,0],'2':[22562,22562,22562,22562,22562,22562,22562,22562,22562,22562,22562,22562,0,0,0,0],'3':[18492,18492,18492,18492,18492,18492,18492,18492,18492,18492,18492,18492,0,0,0,0]}
Tmax= 10
T= [i for i in range(0,Tmax+1)]
# =============================================================================Y={}
for a in A:
for t in T:
Y[a,t]=model.addVar(name="Y_%s_%s" % (a,t), vtype=GRB.BINARY)Au={}
for t in T:
Au[t]=model.addVar(name="Au_%s" % (t), vtype=GRB.CONTINUOUS, lb=1*10**30)
BW={}
for t in T:
BW[t]=model.addVar(name="BW_%s" % (t), vtype =GRB.CONTINUOUS, lb=1*10**30)
BWX= model.addVar(vtype =GRB.CONTINUOUS , obj = 1 , name='BWX')
Er={}
for t in T:
Er[t]=model.addVar(name="Er_%s" % (t), vtype=GRB.CONTINUOUS, lb=1*10**30)
R={}
for t in T:
R[t]=model.addVar(name="R_%s" % (t), vtype =GRB.CONTINUOUS, lb=1*10**30)
# =============================================================================
model.addConstr(BWX==BW[Tmax])for t in T:
model.addConstr(Au[t] == quicksum(KWS[a][t] * Y[a,t] for a in A))
model.addConstr(Er[t] == quicksum(KWE[a][t] * Y[a,t] for a in A))
for t in T:
model.addConstr(BW[t]==Er[t]+Au[t])
for t in T:
model.addConstr(R[t]*Au[t]==Er[t]) #This line of code makes trouble# =============================================================================
model.update()
model.optimize()
# =============================================================================0 
It looks like you just need to set the NonConvex parameter to 2 to direct Gurobi to solve the model (which is made nonconvex by the quadratic equality constraints):
model.Params.NonConvex = 2
Note that solving nonconvex quadratic problems is only supported in Gurobi 9.0 onward.
0 
Thanks a lot. It works
0
Please sign in to leave a comment.
Comments
7 comments