Divisor must be a constant
Answered
Hello,
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
-
Official comment
This post is more than three years old. Some information may not be up to date. For current information, please check the Gurobi Documentation or Knowledge Base. If you need more help, please create a new post in the community forum. Or why not try our AI Gurobot?. -
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 risk-adjusted 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 the following model, which we refer to as \( {\bf(3)} \):
$$\begin{alignat*}{2} \min_y\ && y^\top Q y \enspace & \\ \textrm{s.t.}\ && \mu^\top y = {} &1 \\ && y \geq {} &0.\end{alignat*}$$
\( {\bf(3)} \) 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 \).
3 -
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, self-contained 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 non-convex by the quadratic equality constraints):
model.Params.NonConvex = 2Note that solving non-convex quadratic problems is only supported in Gurobi 9.0 onward.
0 -
Thanks a lot. It works
0 -
Hi Eli Towle - in your first solution above you state, "A solution y_bar to the above problem can be mapped to the original problem using the same transformation"... Essentially, divide the optimized values by the sum of the values and it'll give you the correct values. This works for me, but I can't find a way in Python to update the Gurobi variable values.
Here are the values I get after optimizing:
model.getVars()
Out[445]:
[<gurobi.Var Stock[adj_close,CVX] (value 0.2238213647992848)>,
<gurobi.Var Stock[adj_close,GLD] (value 2.89692620697449)>,
<gurobi.Var Stock[adj_close,JNJ] (value 0.8208641232459375)>,
<gurobi.Var Stock[adj_close,JPM] (value 0.360227339970831)>,
<gurobi.Var Stock[adj_close,KO] (value 5.721567870427095e-06)>,
<gurobi.Var Stock[adj_close,MSFT] (value 2.066410404873706)>,
<gurobi.Var Stock[adj_close,VZ] (value 2.4384480927725286e-08)>]I then run this:
portfolio_weights = stocks.gppd.X/sum(stocks.gppd.X)
and it gives me the values below which are the correct values:
adj_close CVX 0.0351
GLD 0.4549
JNJ 0.1289
JPM 0.0566
KO 0.0000
MSFT 0.3245
VZ 0.0000I just want to divide by the sum, or do the transformation inside of the Gurobi variable so I can use it elsewhere.
Thanks, Erik
1 -
The values of the X variable attributes come from the model solution and cannot be manually edited. That said, why do you want to edit the variable values in the first place? With the code
portfolio_weights = stocks.gppd.X/sum(stocks.gppd.X)
you create a Python object that contains the portfolio weights. You can use these portfolio weights wherever needed.
0 -
hi Eli Towle, Thanks for the response. I'm trying to develop an efficient frontier (Efficient Frontier: What It Is and How Investors Use It (investopedia.com)) by modeling many scenarios with the same objective function. I'm using this example (portfolio_selection_optimization_gcl.ipynb - Colaboratory (google.com))
When it iterates through multiple scenarios, it pulls out a different portfolio weight for each iteration. My thought is that if I could update the variables in the model, I wouldn't need to edit every solution of the iteration.
Thanks - Erik
0 -
What do you mean by "edit every solution of the iteration"? Are you referring to mapping a solution to the convex QP \( {\bf(3)} \) to the original (non-convex) problem \( {\bf(1)} \)? If so, this is a very fast operation, especially relative to the amount of work required to solve a model.
To compute the efficient frontier, we solve a separate instance of model \( {\bf(3)} \) for different target risk/return values. To compute the "correct" portfolio weights, we transform our solution to the space of model \( {\bf(1)} \) with the mapping \( \bar{x}_i := \bar{y}_i / \sum_{j=1}^n \bar{y}_j \) for \( i = 1, \ldots, n \). The corresponding Sharpe ratio is calculated as \( \mu^\top \bar{x} / \sqrt{\bar{x}^\top Q \bar{x}} \). In each iteration, you can simply store the computed portfolio weights and Sharpe ratio for later. It's not necessary to fix/inject a solution into model \( {\bf(3)} \). In general, a solution to \( {\bf(1)} \) will not even satisfy the constraints of model \( {\bf(3)} \).
0 -
Hallo everyone,
I do have a similar problem with a variable as a divisor, but i could not find a workaround based on the previous answers. Can anybody held me based on the following minimal example? The error is caused in Constraint 2.
#----------import packages-----------------------------------------------------
import numpy as np
from scipy.stats import norm
import gurobipy as gp
from gurobipy import GRB#----------parameters----------------------------------------------------------
p = 10 #price for selling a product
c1 = 5 #costs for pre-order
c2 = 7 #costs for re-order
T = 7 #number of periods
beta = 0.9 #target service level required
times, d, std = gp.multidict({-2: [0,0], -1: [0,0], 0: [0,0],
1: [100,10], 2: [200,20], 3: [300,30], 4: [400,40], 5: [500,50], 6: [600,60], 7: [700,70],
8: [0,0], 9: [0,0]}) #periods, forecasted average and standard deviation
time = [*range(0,T+1)]
M = 100000000 #large number
#----------model---------------------------------------------------------------m = gp.Model("Production planning")
Q1 = m.addVars(time, vtype=gp.GRB.INTEGER, lb=0, ub=M, name="Pre-eorder")
Q2 = m.addVars(time, vtype=gp.GRB.INTEGER, lb=0, ub=M, name="Re-order")def L(q, mu, sigma): #beta service level from expected lost sales depending on the inventory, and demand distribution
x = (q-mu)/sigma
y = np.exp(-(x**2)/2)/m.sqrt(2*m.pi) - x * (1 - norm.cdf(x))
z = 1-(y*sigma/mu)
return(z)Con01 = m.addConstrs(Q2[t] >= 0.1*Q1[t] for t in time)
Con02 = m.addConstrs(L(Q1[t]+Q2[t], d[t], std[t]) >= beta for t in time) #this line causes the ERROR, respectively the arguments for my defined function L.
obj = gp.quicksum(p*d[t]-(c1*Q1[t]+c2*Q2[t]) for t in time)m.setObjective(obj, GRB.MAXIMIZE) #Objective minimize holding costs, variable costs
m.optimize()
#----------outputs-------------------------------------------------------------
print(Q1[1],Q1[2],Q1[3],Q1[4],Q1[5],Q1[6],Q1[7],Q2[1],Q2[2],Q2[3],Q2[4],Q2[5],Q2[6],Q2[7])0 -
In general, you can't build a model's constraints or objective using arbitrary functions. The following two articles contain more information about the types of models you can solve with Gurobi:
Perhaps you could build a piecewise-linear approximation of your \(\texttt{L}\) function using Model.addGenConstrPWL().
0
Post is closed for comments.
Comments
14 comments