avoid infinity from logarithmic function in the objective
回答済みIs it possible to avoid infinity from the logarithmic function in objective?
In particular, I am trying to implement the log barrier function method for constrained optimization problem with the Gurobi addGenConstrLog() function (I know I can simply use the addConstr() function, but this is for the research purpose). What I did is something like following:
M = gp.Model()
x = M.addVars(n,m, lb = 0, ub = 1e5, vtype=GRB.CONTINUOUS)
% add some equality constraits
for i in someEqualityIndexSet:
M.addConstrs(x.sum(i,'*') == h[i])
%set objective funciton
objExp = gp.QuadExpr()
% add quadratic objective function
for i in range(n):
for j in range(m):
objExp.add(x[i,j], c[i,j]*x[i,j]])
% add barrier function for inequality constraints
gamma = 0.001
for i in someInequalityIdxSet:
mu = M.addVar(lb = 0, ub = 1e20, vtype = GRB.CONTINUOUS)
M.addConstr(mu == g[i] - x.sum(i,'*'))
y = M.addVar(lb = -1e20, ub = 1e20, vtype = GRB.CONTINUOUS)
M.addGenConstrLog(mu , y,"FuncPieces = -2 FuncPieceError=0.0001")
objExp.add(-gamma*y)
M.setObjective(objExp, gp.MINIMIZE)
M.update()
M.optimize()
After the model is solved, there are some active inequality constraints, i.e some mu is 0, which should make the barrier function goes to infinity. I wonder if there is some way I can change the model such that the log barrier function will work as intended, that is some variable mu close to 0 but not exactly 0.
I hope I stated my question clear. Thank you in advance!
-
正式なコメント
-
Gurobi Staff
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?. -
-
-
Hi Jiazhen,
Gurobi solves both convex and non-convex QCPs and uses the barrier algorithm for solving the convex QCP. Therefore, we highly recommend using Gurobi directly to solve your problem. That being said, let us address your question.
If the barrier function is represented correctly, it should be impossible for the values of the \(\mu\) variables to equal zero because \(- \gamma y = -\gamma log(\mu)\) would then equal \(\infty\). This cannot happen in a minimization problem. There is likely another issue.
Let us look at the qcp.py problem as a concrete example below. The example problem has the optimal objective function of 0.326992. Note that this is a maximization problem.
\[\begin{align} \max & ~~x \\ \text{s.t.} & ~~x+y+z = 1 \\ & x^2 + y^2 \leq z^2 \\ & x^2 \leq yz \\ & x, y, z \geq 0\end{align}\]
The logarithmic barrier function associated with the above optimization problem is:
\[\begin{align} B(\gamma) : \max & ~~x + \gamma \big[\ln(z^2 - x^2 -y^2) + \ln(yz - x^2)+ \ln(x) + \ln(y) + \ln(z)\big] \\ \text{s.t. } & ~~x+y+z = 1 \end{align}\]
Here \(\gamma\) is a positive scalar and as it approaches 0, the solution of \(B(\gamma)\) should converge to the solution of the qcp.py example problem.
Running the code snippet (copied below) shows that the value of the objective function \(B(\gamma)\) approaches the objective function of the qcp.py example problem as the \(\gamma\) decreases.
gamma:1 -- obj: -91.776410
gamma:0.1 -- obj: -8.883347
gamma:0.01 -- obj: -0.594040
gamma:0.001 -- obj: 0.234890
gamma:0.0001 -- obj: 0.317783
gamma:1e-05 -- obj: 0.326073
gamma:1e-06 -- obj: 0.326901def solve_barrier_function():
def solve_per_mu(gamma):
model = gurobipy.Model("")
x = model.addVar(name="x")
y = model.addVar(name="y")
z = model.addVar(name="z")
logx = model.addVar(lb=-GRB.INFINITY, name="logx")
logy = model.addVar(lb=-GRB.INFINITY, name="logy")
logz = model.addVar(lb=-GRB.INFINITY, name="logz")
# Auxiliary variables
w = model.addVar(name="w")
v = model.addVar(name="v")
logw = model.addVar(lb=-GRB.INFINITY, name="logw")
logv = model.addVar(lb=-GRB.INFINITY, name="logv")
# Set barrier objective
obj = 1.0 * x + gamma * (logw + logv + logx + logy + logz)
model.setObjective(obj, GRB.MAXIMIZE)
# Add constraint: x+y+z=1
model.addConstr(x + y + z == 1, "c0")
# Add barrier functions associated with non-negativity constraints
model.addGenConstrLog(x, logx)
model.addGenConstrLog(y, logy)
model.addGenConstrLog(z, logz)
# Add w=z^2-x^2-y^2 (this quadratic equality constraint makes the
# problem non-convex)
model.addConstr(w == z ** 2 - x ** 2 - y ** 2, "w=z^2-x^2-y^2")
# Add barrier associated with z^2-x^2-y^2 >= 0
model.addGenConstrLog(w, logw)
# Define v=yz-x^2 (this quadratic equality constraint makes the
# problem non-convex)
model.addConstr(v == y * z - x ** 2, "v=yz-x^2")
# Add barrier associated with yz-x^2 >= 0
model.addGenConstrLog(v, logv)
model.params.NonConvex = 2
model.params.OutputFlag = 0
model.optimize()
print(f"gamma:{gamma} -- obj: {model.objVal:.6f}")
for gamma in [1, 0.1, 0.01, 0.001, 0.0001, 0.00001, 0.000001]:
solve_per_mu(gamma)I would like to emphasize again that it is much more efficient and numerically stable to use Gurobi directly.
Best regards,
Maliheh
0 -
-
Hi Maliheh,
Thank you for this illustrative example. In running this example, however, there are still active constraints. In particular, variable \(v\) is zero for all \(\gamma\). Following is the printed results for \(v\) and \(\log(v)\):
gamma:1 -- obj: -91.776408
v:0.0 logv:-23.025850929940457
Changed value of parameter NonConvex to 2
Prev: -1 Min: -1 Max: 2 Default: -1
gamma:0.1 -- obj: -8.883347
v:0.0 logv:-23.025850929940457
Changed value of parameter NonConvex to 2
Prev: -1 Min: -1 Max: 2 Default: -1
gamma:0.01 -- obj: -0.594041
v:0.0 logv:-23.02585092994046
Changed value of parameter NonConvex to 2
Prev: -1 Min: -1 Max: 2 Default: -1
gamma:0.001 -- obj: 0.234891
v:0.0 logv:-23.025850929940457
Changed value of parameter NonConvex to 2
Prev: -1 Min: -1 Max: 2 Default: -1
gamma:0.0001 -- obj: 0.317784
v:0.0 logv:-23.025850929940457
Changed value of parameter NonConvex to 2
Prev: -1 Min: -1 Max: 2 Default: -1
gamma:1e-05 -- obj: 0.326073
v:0.0 logv:-23.025850929940457
Changed value of parameter NonConvex to 2
Prev: -1 Min: -1 Max: 2 Default: -1
gamma:1e-06 -- obj: 0.326902
v:0.0 logv:-23.025850929940457Do you think it is due to some precision setting for the log function or some other general setting for the model?
Thank you again for the help!
Best,
Jiazhen
0 -
-
-
Hi Jiazhen,
Running the snippet code posted before results into \(\texttt{v}\) being equal to \(10^{-10}\) and not zero. The log value that you posted is the \(log(10^{-10})\). When running the code, I get the following:
gamma:1 -- obj: -91.776410
v: 1e-10 -- logv: -23.02585092993819
gamma:0.1 -- obj: -8.883347
v: 1e-10 -- logv: -23.02585092993819
gamma:0.01 -- obj: -0.594040
v: 1e-10 -- logv: -23.02585092993819
gamma:0.001 -- obj: 0.234890
v: 1e-10 -- logv: -23.02585092994046
gamma:0.0001 -- obj: 0.317783
v: 1e-10 -- logv: -23.02585092994046
gamma:1e-05 -- obj: 0.326073
v: 1e-10 -- logv: -23.02585092993819
gamma:1e-06 -- obj: 0.326901
v: 1e-10 -- logv: -23.02585092993819Maybe your print function prints up to a certain decimal value.
Best regards,
Maliheh
0 -
-
I changed the tolerance of the problem and it is showing the same results as yours.
Thank you, Maliheh!
Jiazhen
0 -
投稿コメントは受け付けていません。
コメント
5件のコメント