Quadratic constraints

[Midnighter]

Within a short amount of time, I have heard from three different developers that they require quadratic constraints. @KristianJensen you probably have the best idea why those were not implemented yet. Can you give a rough idea of the scope of adding this to optlang?

[erling-d-andersen]

I assume you convex quadratic constraints. If yes, then we might be interested in contributing a Mosek interface if is possible,

Also maybe you should look at conic optimization because that is likely easier to add to your language and you can support a much richer set of convex nonlinearities e.g. ln() and exp().

[KristianJensen]

@Midnighter The main reason it was not implemented was that there was no need for it at the time.
Overall I don't think it would be a lot of work to implement it. The parsing of quadratic expression is already implemented (used for quadratic objectives), and could easily be used to parse quadratic constraints.

Optlang itself imposes no restrictions on the nature of the constraints, so it is only a question of whether solvers support a given type of constraints, and then implementing the functions to translate the expression into whatever format the solver needs.

[matthiaskoenig]

I also need quadratic terms. Please support this in optlang.

[A2P2]

I'm in the process of adding quadratic constraints to the cplex interface. The minimal version is shown below. So far it's possible to pass constraints like (x1^2 + x2^2) and (x1*x2 + x2*x3), but the second type is usually nonconvex and can't be solved.
Parsing of a constraint, containing both quadratic and linear parts (x1^2 + x2) doesn't work. I have to change expression_parsing.py for that to work. I'll do some work more and submit a pull request.

Any suggestions on what else to do or comments regarding the code?

            elif constraint.is_Quadratic:
                offset, linear_coefficients, quadratic_coefficients = parse_optimization_expression(constraint, quadratic=True, expression=constraint.expression)
                
                sense, rhs, _ = _constraint_lb_and_ub_to_cplex_sense_rhs_and_range_value(
                    constraint.lb,
                    constraint.ub
                )
                lin_indices = [var.name for var in linear_coefficients]
                lin_values = [float(val) for val in linear_coefficients.values()]
                
                #Unpack quadratic variables
                quad_indices_1 = []
                quad_indices_2 = []
                #quadratic_coefficients.keys() is a dict_keys() of frozentsets, so we have to unpack it.
                for var_set in quadratic_coefficients.keys():
                    var_list = list(var_set)
                    # If x1*x2, a bit of overkill at the moment, since such problems usually nonconvex
                    if len(var_list) == 2: 
                        quad_indices_1.append(var_list[0].name)
                        quad_indices_2.append(var_list[1].name)
                    # If x1**2
                    elif len(var_list) == 1: 
                        quad_indices_1.append(var_list[0].name)
                        quad_indices_2.append(var_list[0].name)
                    #Which error to raise?
                    # else:
                    #     RaiseError 
                quad_values = [float(val) for val in quadratic_coefficients.values()]
                
                quadratic_constraints['lin_expr'].append(cplex.SparsePair(ind=lin_indices, val=lin_values))
                quadratic_constraints['quad_expr'].append(cplex.SparseTriple(ind1=quad_indices_1, ind2=quad_indices_2, val=quad_values))
                quadratic_constraints['sense'].append(sense)
                quadratic_constraints['rhs'].append(rhs)
                quadratic_constraints['name'].append(constraint.name)
            else:
                raise ValueError(
                    "CPLEX only supports linear or quadratic constraints. %s is neither linear nor quadratic." % constraint)
            constraint.problem = self
        self.problem.linear_constraints.add(**linear_constraints)
        
        #Cplex can not unpack quadratic_constraints the same way as linear_constraints
        #So we unpack constraints ourselves here
        for i in range(len(quadratic_constraints['quad_expr']) ):
            quadratic_constraint = {key:value[i] for key, value in quadratic_constraints.items()}
            self.problem.quadratic_constraints.add(**quadratic_constraint)

[asaldivar93]

@A2P2 I've been working with a similar workaround for gurobi. So far I've found that you need to modify the _get_expression method for constraints, since its build for linear expression only. If you are working with cobrapy, you wiĺl also find a bug with summary methods but i have not dig enough as to find the source of it.

the code i use for gurobi is this:

`
else:

            offset, linear_coefficients, quadratic_coefficients = parse_optimization_expression(
                constraint, quadratic=True,
            )
            
            grb_terms = []
            for var, coef in linear_coefficients.items():
                var = self.problem.getVarByName(var.name)
                grb_terms.append(coef * var)
                
            for key, coef in quadratic_coefficients.items():
                if len(key) == 1:
                    var = six.next(iter(key))
                    var = self.problem.getVarByName(var.name)
                    grb_terms.append(coef * var * var)
                else:
                    var1, var2 = key
                    var1 = self.problem.getVarByName(var1.name)
                    var2 = self.problem.getVarByName(var2.name)
                    grb_terms.append(coef * var1 * var2)

            lhs = gurobipy.quicksum(grb_terms)
            sense, rhs, range_value = _constraint_lb_and_ub_to_gurobi_sense_rhs_and_range_value(constraint.lb,
                                                                                                constraint.ub)
            if range_value != 0:
                aux_var = self.problem.addVar(name=constraint.name + '_aux', lb=0, ub=range_value)
                self.problem.update()
                lhs = lhs - aux_var

            self.problem.addQConstr(lhs, sense, rhs, name=constraint.name)
            self.problem.params.NonConvex = 2

`

[mrivas]

@A2P2 Your code looks very promising. I tried it but I get the following error:

CPLEX Error 1031: Not available for QCP.

I would greatly appreciate if you could share how you avoided this error?