Parametrized QPs

QP_Param stands for Parametrized Quadratic Program, a mathematical program in the following form:

\[\begin{split}\min_y \frac{1}{2}y^TQy + c^Ty + (Cx)^T y + d^T x \\ \text{s.t.} \quad Ax + By \le b\\ \quad \quad y \ge 0\end{split}\]

Where \(y\) are the decision variables for the program, and \(x\) are parameters. You can find the API information in MathOpt::QP_Param. This class is an inheritor of MathOpt::MP_Param.

Modeling the problem

Consider the following quadratic program.

\[ \begin{align}\begin{aligned}\min_{y_1, y_2, y_3} (y_1 + y_2 - 2y_3)^2 + 2 x_1y_1 + 2 x_2y_1 + 3 x_1y_3 + y_1-y_2+y_3\\\text{s.t.}\;\;\;\;\; y_1, y_2, y_3 &\ge 0\\\;\;\;\;\;\;\;\; y_1 + y_2 + y_3 &\le 10\\\;\;\;\;\;\;\;-y_1 +y_2 -2y_3 &\le -1 + x_1 + x_2\end{aligned}\end{align} \]

We can model the problem above as follows:

#include "zero.h"
//[...] Your other code here
unsigned int numParams = 2, numVars = 3, numConstr = 2;

//Create Q (from dense to sparse)
arma::mat Qd(3, 3);
Qd << 1 << 1 << -2 << arma::endr
<< 1 << 1 << -2 << arma::endr
<< -2 << -2 << 4 << arma::endr;
//Convert from dense to sparse
arma::sp_mat Q = sp_mat(2 * Qd);

// The matrix for x and y interaction
arma::sp_mat C(3, 2);
arma::vec d(3);
C.zeros(); C(0, 0) = 2; C(0, 1) = 2; C(2, 0) = 3;
// The vector for linear terms in y
arma::vec c(3);
c << 1 << arma::endr << -1 << arma::endr << 1 << arma::endr;
// Constraint matrix for x terms
arma::sp_mat A(2, 2);
A.zeros(); A(1, 0) = -1; A(1, 1) = -1;
// Constraint matrix for y terms
arma::mat Bd(2, 3); Bd << 1 << 1 << 1 << arma::endr << -1 << 1 << -2 << arma::endr;
arma::sp_mat B = sp_mat(Bd);
//RHSs
arma::vec b(2);
b(0) = 10;
b(1) = -1;
// Create a Gurobi environment to handle any solving related calls.
GRBEnv env = GRBEnv();

Now the required object can be constructed in multiple ways.

// Method 1: Make a call to the constructor
MathOpt::QP_Param q1(Q, C, A, B, c, b, d, &env);

// Method 2: Using QP_Param::set member function
MathOpt::QP_Param q2(&env);
q2.set(Q, C, A, B, c, b, d);

// Method 3: Reading from a file. This requires that such an object is saved to a file at first.
q1.save("dat/q1dat.dat"); // Saving the file so it can be retrieved.
MathOpt::QP_Param q3(&env);
q3.load("dat/q1dat.dat");

// Checking they are the same
assert(q1==q2);
assert(q2==q3);

We can now feed some values for the parameters and compute the corresponding optimal solution. Assume \((x_1, x_2) = (-1, 0.5)\). Then, we the problem becomes a standard QP as:

\[ \begin{align}\begin{aligned}\min_{y_1, y_2, y_3} (y_1 + y_2 - 2y_3)^2 -y_2 -2y_3\\\text{s.t.}\;\;\;\; y_1, y_2, y_3 &\ge 0\\\;\;\;\;\;\;\;\;y_1 + y_2 + y_3 &\le 10\\\;\;\;\;\;\;\;\;-y_1 +y_2 -2y_3 &\le -1.5\end{aligned}\end{align} \]

Correspondingly, we have the following code:

// Enter the value of x in an arma::vec
arma::vec x(2);
x(0) = -1;
x(1) = 0.5;

// Uses Gurobi to solve the model, returns a unique_ptr to GRBModel.
// With the second parameters, we require a model which has already been solved
auto FixedModel = q2.solveFixed(x,true);

Computing solutions

FixedModel holds the GRBModel object, and all operations native to GRBModel, like accessing the value of a variable, a dual multiplier, saving the problem to an .lp file or a .mps file. In particular, we can compare the solution with a hand-calculated one as below.

arma::vec sol(3);
// Hard-coding the solution as calculated outside
sol << 0.5417 << arma::endr << 5.9861 << arma::endr << 3.4722;
for (unsigned int i = 0; i < numVars; i++)
    assert(abs(sol(i)- FixedModel->getVar(i).get(GRB_DoubleAttr_X)) <= 1e-5);
cout<<FixedModel->get(GRB_DoubleAttr_ObjVal<<endl; // Will print -12.757

In many cases, one might need the “KKT” conditions of a convex quadratic program. We can then use the function MathOpt::QP_Param::KKT():

The function returns M, N and q, where the KKT conditions can be written as \(0 \leq y \perp Mx + Ny + q \geq 0\).

arma::sp_mat M, N; arma::vec q;
q1.KKT(M, N, q);
M.print("M");
N.print("N");
q.print("q");