avoid infinity from logarithmic function in the objective

Answered

November 10, 2021 22:45

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!

Comments

5 comments

Official comment
Simranjit Kaur

Gurobi Staff

June 23, 2025 02:44 Edited
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?.

Maliheh Aramon

Gurobi Staff

November 13, 2021 08:01

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.326901

def 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

Jiazhen Zhou

November 15, 2021 15:24 Edited

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.025850929940457

Do 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

Maliheh Aramon

Gurobi Staff

November 15, 2021 19:11

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.02585092993819

Maybe your print function prints up to a certain decimal value.

Best regards,

Maliheh

Jiazhen Zhou

November 16, 2021 17:43
I changed the tolerance of the problem and it is showing the same results as yours.

Thank you, Maliheh!

Jiazhen
0

Post is closed for comments.

avoid infinity from logarithmic function in the objective

Comments

Didn't find what you were looking for?