[Request] Implementation of a solver for the quadratic programming using the CPLEX Python API

[Yu-ki-Koyama]

Hi. I usually work with @mmmarinho. This time, I had trouble solving the quadratic programming and implemented a solver using the CPLEX Python API. @mmmarinho suggested me to share this code here, so I would like to do that.

Trouble

I'm using your DQ Robotics to develop our robotic system. And, I've utilized DQ_QuadprogSolver to solve the quadratic programming. Recently, I've been trying to solve tasks with priority. In this process, I need to solve two serial quadratic programmings, and I add the result of the first one to the second one as an equality constraint to ensure the first cost function remains minimal. Though I don't comment on the details, the second quadratic programming must have at least one answer theoretically. However, my program stops while solving the second one with an error like 'constraints are inconsistent, no solution.'

What I did

I chose to use another solver and utilized the CPLEX Python API. Then, my program never stops with errors. I'm not sure the error was because of DQ_QuadprogSolver, but I wrote a class file of the CPLEX Python API, so I would like to share it.

Cplex_QuadraticSolver.py

try:
    import cplex
    import numpy as np

    class Cplex_QuadraticSolver():
        def __init__(self):
            # Make and set a solver instance
            self.P = cplex.Cplex()
            self.P.objective.set_sense(self.P.objective.sense.minimize)
            self.P.set_problem_type(self.P.problem_type.QP)
            self.P.set_results_stream(results_file=None)

        def solve_quadratic_program(self, H, f, A, b, Aeq, beq):
            problem_size = H.shape[0]
            inequality_constraint_size = b.shape[0]
            equality_constraint_size = beq.shape[0]

            # Define the number of variables by naming them
            q_dot = ["q{}".format(i) for i in range(problem_size)]
            if self.P.variables.get_num() != problem_size:
                self.P.variables.delete()
                self.P.variables.add(names=q_dot)

            # Decide the range of variables
            for i in range(problem_size):
                self.P.variables.set_lower_bounds(i, -1*cplex.infinity)
                self.P.variables.set_upper_bounds(i, cplex.infinity)

            # Set f'x
            f = f.tolist()
            for i in range(problem_size):
                self.P.objective.set_linear(i, f[i][0])

            # Set x^T H x
            H = H.tolist()
            for i in range(problem_size):
                for j in range(i, problem_size):
                    self.P.objective.set_quadratic_coefficients(i, j, H[i][j])

            # Set linear inequality constraints and equality constraints
            # Define the number of total constraints by naming them
            c = ["c{}".format(i) for i in range(inequality_constraint_size + equality_constraint_size)]
            if self.P.linear_constraints.get_num() != inequality_constraint_size + equality_constraint_size:
                self.P.linear_constraints.delete()
                self.P.linear_constraints.add(names=c)

            # Inequalities
            A = A.tolist()
            b = b.tolist()
            for i in range(inequality_constraint_size):
                self.P.linear_constraints.set_linear_components(c[i], [q_dot, A[i]])
                self.P.linear_constraints.set_senses(c[i], 'L')
                self.P.linear_constraints.set_rhs(c[i], b[i][0])

            # Equalities
            Aeq = Aeq.tolist()
            beq = beq.tolist()
            for i in range(equality_constraint_size):
                self.P.linear_constraints.set_linear_components(c[inequality_constraint_size + i], [q_dot, Aeq[i]])
                self.P.linear_constraints.set_senses(c[inequality_constraint_size + i], 'E')
                self.P.linear_constraints.set_rhs(c[inequality_constraint_size + i], beq[i][0])

            # Solve the problem
            self.P.solve()
            delta_thetas = np.array(self.P.solution.get_values())

            return delta_thetas

except:
    pass

How to use

I use this class file as following.

qp_solver = Cplex_QuadraticSolver()

# When we don't have any equality constraints
delta_thetas = qp_solver.solve_quadratic_program(2 * (H + damping * np.eye(joint_number)), c, W, w, np.zeros([1, joint_number]), np.zeros([1, 1]))

# When we have some equality constraints
delta_thetas = qp_solver.solve_quadratic_program(2 * (H + damping * np.eye(joint_number)), c, W, w, Aeq, beq)

Thank you.

[mmmarinho]

Thank you for your contribution, @Yu-ki-Koyama.

I added it to DQ_Robotics Python as DQ_CPLEXSolver. Could you test it and let me know how it goes?

[Yu-ki-Koyama]

Hi, @mmmarinho. Thank you. I tested it, and it worked well!

[mmmarinho]

@Yu-ki-Koyama Great! Thanks again for the contribution.