Quadratic Programming

Quadratic programs use a specialized kind of decision variable, LinearProgramVariable. The objective function has a linear component. The coefficient of each variable in this linear component is usually called the cost. It can be set or retrieved through the Cost property.

By convention, quadratic program variables have a lower bound of zero by default.

The LowerBound and UpperBound properties define constraints on the values of a variable. Standard variables are positive. They have a lower bound of zero and no upper bound. When a variable is not bounded from above, its upper bound is DoublePositiveInfinity. When a variable is not bounded from below, its lower bound is DoubleNegativeInfinity.

The LinearConstraint class is used to represent constraints in a quadratic program. Every constraint must have a unique name, which can be accessed through the Name property.

A constraint is a linear combination of variables that is required to lie between specified boundaries. The LowerBound and UpperBound properties define these boundaries. When a constraint is not bounded from above, its upper bound is DoublePositiveInfinity. When a constraint is not bounded from below, its lower bound is DoubleNegativeInfinity.

Variables and constraints can be accessed through the quadratic program's Variables and Constraints collections. These collections are indexed by position and by name.

A quadratic program is represented by the QuadraticProgram class. There are three ways to construct QuadraticProgram objects.

The quickest way to define a quadratic program is by passing the vectors and matrices that appear in the definition above to a constructor. There are three constructors that can be used in this way.

The first constructor takes four arguments. It defines a quadratic program in standard form. This is a quadratic program where variables are constrained to be positive, and all constraints are inequalities with an upper bound. The first argument is a VectorT containing the coefficients of the linear portion of the objective function. It corresponds to the vector c in the definition. The second argument is a MatrixT that represents the Hessian and contains the coefficients of the quadratic terms in the objective function. This matrix must be symmetrical. It corresponds to the vector H in the definition.

The third argument is a MatrixT that contains the coefficients of the inequalities. It corresponds to the matrix A in the definition. The final argument is a VectorT containing the right-hand-sides of the inequalities. The code below creates the following a quadratic program:

Vector<double> c = Vector.CreateConstant(4, 0.0);
Matrix<double> R = Matrix.CreateSymmetric(4,
    new double[] {
         0.08,-0.05,-0.05,-0.05,
        -0.05, 0.16,-0.02,-0.02,
        -0.05,-0.02, 0.35, 0.06,
        -0.05,-0.02, 0.06, 0.35
    }, MatrixTriangle.Upper, MatrixElementOrder.RowMajor);
Matrix<double> A = Matrix.Create(2, 4, new double[]
{
    1, 1, 1, 1,
    -0.05, 0.2, -0.15, -0.30
}, MatrixElementOrder.RowMajor);
Vector<double> b = Vector.Create(10000.0, -1000.0);
QuadraticProgram qp1 = new QuadraticProgram(c, R, A, b, 0);

Dim c As Vector(Of Double) = Vector.CreateConstant(4, 0.0)
Dim R As Matrix(Of Double) = Matrix.CreateSymmetric(4, New Double() {
         0.08, -0.05, -0.05, -0.05,
        -0.05, 0.16, -0.02, -0.02,
        -0.05, -0.02, 0.35, 0.06,
        -0.05, -0.02, 0.06, 0.35
    }, MatrixTriangle.Upper, MatrixElementOrder.ColumnMajor)
Dim A As Matrix(Of Double) = Matrix.Create(2, 4, New Double() {
    1, 1, 1, 1,
    -0.05, 0.2, -0.15, -0.3
}, MatrixElementOrder.RowMajor)
Dim b As Vector(Of Double) = Vector.Create(10000.0, -1000.0)
Dim qp1 As QuadraticProgram = New QuadraticProgram(c, R, A, b, 0)

let c = Vector.CreateConstant(4, 0.0)
let R = Matrix.CreateSymmetric(4, 
            [|
                 0.08;-0.05;-0.05;-0.05;
                -0.05; 0.16;-0.02;-0.02;
                -0.05;-0.02; 0.35; 0.06;
                -0.05;-0.02; 0.06; 0.35
            |], MatrixTriangle.Upper, MatrixElementOrder.ColumnMajor)
let A = Matrix.Create(2, 4, 
            [|
                1.0; 1.0; 1.0; 1.0;
                -0.05; 0.2; -0.15; -0.30
            |], MatrixElementOrder.RowMajor)
let b = Vector.Create(10000.0, -1000.0)
let qp1 = QuadraticProgram(c, R, A, b, 0)

The third constructor is the most flexible. It takes seven arguments in all. The first two are once again the linear and quadratic parts of the objective function. The third argument is the coefficient matrix. The fourth and fifth parameters are vectors that contain the lower and upper bounds of the constraints. Each type of constraint can be expressed using a lower and upper bound, as shown in the table below.

The sixth and seventh parameters are vectors that contain the lower and upper bound for the variables. The different types of constraints on variables can be expressed the same way as for constraints. The following code sample builds the same quadratic program:

Vector<double> c = Vector.CreateConstant(4, 0.0);
Matrix<double> R = Matrix.CreateSymmetric(4,
    new double[] {
         0.08,-0.05,-0.05,-0.05,
        -0.05, 0.16,-0.02,-0.02,
        -0.05,-0.02, 0.35, 0.06,
        -0.05,-0.02, 0.06, 0.35
    }, MatrixTriangle.Upper, MatrixElementOrder.RowMajor);
Matrix<double> A = Matrix.Create(2, 4, new double[]
{
    1, 1, 1, 1,
    -0.05, 0.2, -0.15, -0.30
}, MatrixElementOrder.RowMajor);
Vector<double> cl = Vector.CreateConstant(2, double.PositiveInfinity);
Vector<double> cu = Vector.Create(10000.0, -1000.0);
Vector<double> vl = Vector.CreateConstant(2, 0.0);
Vector<double> vu = Vector.CreateConstant(2, double.PositiveInfinity);
QuadraticProgram qp1 = new QuadraticProgram(c, R, A, cl, cu, vl, vu);

Dim c As Vector(Of Double) = Vector.CreateConstant(4, 0.0)
Dim R As Matrix(Of Double) = Matrix.CreateSymmetric(4,
    New Double() {
         0.08, -0.05, -0.05, -0.05,
        -0.05, 0.16, -0.02, -0.02,
        -0.05, -0.02, 0.35, 0.06,
        -0.05, -0.02, 0.06, 0.35
    }, MatrixTriangle.Upper, MatrixElementOrder.ColumnMajor)
Dim A As Matrix(Of Double) = Matrix.Create(2, 4, New Double() {
    1, 1, 1, 1,
    -0.05, 0.2, -0.15, -0.3
}, MatrixElementOrder.RowMajor)
Dim cl As Vector(Of Double) = Vector.CreateConstant(2, Double.PositiveInfinity)
Dim cu As Vector(Of Double) = Vector.Create(10000.0, -1000.0)
Dim vl As Vector(Of Double) = Vector.CreateConstant(2, 0.0)
Dim vu As Vector(Of Double) = Vector.CreateConstant(2, Double.PositiveInfinity)
Dim qp1 As QuadraticProgram = New QuadraticProgram(c, R, A, cl, cu, vl, vu)

let c = Vector.CreateConstant(4, 0.0)
let R = Matrix.CreateSymmetric(4,
            [|
                 0.08;-0.05;-0.05;-0.05;
                -0.05; 0.16;-0.02;-0.02;
                -0.05;-0.02; 0.35; 0.06;
                -0.05;-0.02; 0.06; 0.35
            |], MatrixTriangle.Upper, MatrixElementOrder.ColumnMajor)
let A = Matrix.Create(2, 4, 
            [|
                1.0; 1.0; 1.0; 1.0;
                -0.05; 0.2; -0.15; -0.30
            |], MatrixElementOrder.RowMajor)
let cl = Vector.CreateConstant(2, infinity)
let cu = Vector.Create(10000.0, -1000.0)
let vl = Vector.CreateConstant(2, 0.0)
let vu = Vector.CreateConstant(2, infinity)
let qp1 = QuadraticProgram(c, R, A, cl, cu, vl, vu)

The QuadraticProgram class also has a default constructor that creates an empty quadratic program. You can then call the quadratic program's AddVariable and AddLinearConstraint methods.

The AddVariable method has three overloads. The first has one parameter: the name of the variable. The second overload takes two arguments. The first is again the name of the variable. The second is the cost associated with the variable. The third overload has two additional parameters that specify the lower and upper bound of the variable. If no values are specified, a cost of zero, a lower bound of zero and an upper bound of infinity are assumed.

The AddLinearConstraint method has four overloads. The first takes four parameters. The first is the name of the constraint. The second is a VectorT that contains the coefficients of the variables. The third parameter is a ConstraintType value that specifies the type of constraint. The possible values are listed in the table below. The last parameter is the right-hand side of the constraint.

The second overload also takes four parameters. The first two are the same as before. The third and fourth parameters are the lower and upper bound of the constraint. The different types of constraints can be specified according to the examples in Table 1.

The third and fourth overloads are similar to the first two. The only difference is that the coefficients are specified using a Double array instead of a vector.

Finally, the SetQuadraticCoefficient method lets you specify the coefficients of a quadratic term. This method has two overloads, which both take 3 parameters. In the first overload, the first two parameters are the variables for which the coefficient is being set. The third parameter contains this value. In the second overload, the first two parameters are the names of the variables.

Note that when a coefficient is specified in this way for a cross term with two different variables, its value is usually twice its corresponding value in the Hessian matrix. This is because each combination appears twice in the matrix.

The following example creates a quadratic program that is identical to the last example:

QuadraticProgram qp2 = new QuadraticProgram();
qp2.AddVariable("X1", 0.0);
qp2.AddVariable("X2", 0.0);
qp2.AddVariable("X3", 0.0);
qp2.AddVariable("X4", 0.0);
qp2.AddLinearConstraint("C1", Vector.Create(1.0, 1.0, 1.0, 1.0), ConstraintType.LessThanOrEqual, 10000);
qp2.AddLinearConstraint("C2", Vector.Create(0.05, -0.2, 0.15, 0.3), ConstraintType.GreaterThanOrEqual, 1000);
qp2.SetQuadraticCoefficient("X1", "X1", 0.08);
qp2.SetQuadraticCoefficient("X1", "X2", -0.05 * 2);
qp2.SetQuadraticCoefficient("X1", "X3", -0.05 * 2);
qp2.SetQuadraticCoefficient("X1", "X4", -0.05 * 2);
qp2.SetQuadraticCoefficient("X2", "X2", 0.16);
qp2.SetQuadraticCoefficient("X2", "X3", -0.02 * 2);
qp2.SetQuadraticCoefficient("X2", "X4", -0.02 * 2);
qp2.SetQuadraticCoefficient("X3", "X3", 0.35);
qp2.SetQuadraticCoefficient("X3", "X4", 0.06 * 2);
qp2.SetQuadraticCoefficient("X4", "X4", 0.35);

Dim qp2 As QuadraticProgram = New QuadraticProgram()
qp2.AddVariable("X1", 0.0)
qp2.AddVariable("X2", 0.0)
qp2.AddVariable("X3", 0.0)
qp2.AddVariable("X4", 0.0)
qp2.AddLinearConstraint("C1", Vector.Create(1.0, 1.0, 1.0, 1.0), ConstraintType.LessThanOrEqual, 10000)
qp2.AddLinearConstraint("C2", Vector.Create(0.05, -0.2, 0.15, 0.3), ConstraintType.GreaterThanOrEqual, 1000)
qp2.SetQuadraticCoefficient("X1", "X1", 0.08)
qp2.SetQuadraticCoefficient("X1", "X2", -0.05 * 2)
qp2.SetQuadraticCoefficient("X1", "X3", -0.05 * 2)
qp2.SetQuadraticCoefficient("X1", "X4", -0.05 * 2)
qp2.SetQuadraticCoefficient("X2", "X2", 0.16)
qp2.SetQuadraticCoefficient("X2", "X3", -0.02 * 2)
qp2.SetQuadraticCoefficient("X2", "X4", -0.02 * 2)
qp2.SetQuadraticCoefficient("X3", "X3", 0.35)
qp2.SetQuadraticCoefficient("X3", "X4", 0.06 * 2)
qp2.SetQuadraticCoefficient("X4", "X4", 0.35)

let qp2 = QuadraticProgram()
qp2.AddVariable("X1", 0.0) |> ignore
qp2.AddVariable("X2", 0.0) |> ignore
qp2.AddVariable("X3", 0.0) |> ignore
qp2.AddVariable("X4", 0.0) |> ignore
qp2.AddLinearConstraint("C1", 
    Vector.Create(1.0, 1.0, 1.0, 1.0), 
    ConstraintType.LessThanOrEqual, 10000.0) |> ignore
qp2.AddLinearConstraint("C2", 
    Vector.Create(0.05, -0.2, 0.15, 0.3), 
    ConstraintType.GreaterThanOrEqual, 1000.0) |> ignore
qp2.SetQuadraticCoefficient("X1", "X1", 0.08)
qp2.SetQuadraticCoefficient("X1", "X2", -0.05 * 2.0)
qp2.SetQuadraticCoefficient("X1", "X3", -0.05 * 2.0)
qp2.SetQuadraticCoefficient("X1", "X4", -0.05 * 2.0)
qp2.SetQuadraticCoefficient("X2", "X2", 0.16)
qp2.SetQuadraticCoefficient("X2", "X3", -0.02 * 2.0)
qp2.SetQuadraticCoefficient("X2", "X4", -0.02 * 2.0)
qp2.SetQuadraticCoefficient("X3", "X3", 0.35)
qp2.SetQuadraticCoefficient("X3", "X4", 0.06 * 2.0)
qp2.SetQuadraticCoefficient("X4", "X4", 0.35)

Although not strictly necessary, it is advisable to create the variables first, followed by the constraints. If the vector containing the coefficients for a constraint is longer than the number of variables that have been defined, additional variables are added automatically. The kth variable is given the name x<g>.

Whenever possible, constraints on variables should be preferred over explicit quadratic program constraints. The time to solve a quadratic program depends most heavily on the number of constraints.

Use the MpsReader class to read a QuadraticProgram from a file in MPS format. The format was named after an early linear programming system from IBM. It has since become a de facto standard for defining linear and quadratic programs in text format in commercial optimization systems. There are many online resources that describe the format in detail.

The MpsReader class has two method, ReadQuadraticProgram, which is overloaded. The first method takes a string containing the path to an MPS file and returns a QuadraticProgram object that represents the quadratic program in the file. The second method does the same but gets its input from a System.IOStreamReader.

Several extensions of the MPS format have been defined over the years. The current implementation of MpsReader doesn't support these extensions. A SystemFormatException is thrown when an unrecognized extension is found.

Once a quadratic program has been defined, solving it is very straightforward. The Solve method does all the work.

The Solve method returns a VectorT containing the optimal solution. You should inspect the Status property to make sure an optimal solution was indeed found. This property is of type OptimizationModelStatus and can take on the following values:

Finally, the OptimalValue property returns the value of the objective function at the solution.

The code below solves the quadratic program defined earlier:

Vector<double> x = qp2.Solve();
Console.WriteLine("Status: {0:F1}", qp2.Status);
Console.WriteLine("Solution: {0:F1}", x);
Console.WriteLine("Optimal value:   {0:F1}", qp2.OptimalValue);

Dim x As Vector(Of Double) = qp2.Solve()
Console.WriteLine("Status: 0:F1", qp2.Status)
Console.WriteLine("Solution: 0:F1", x)
Console.WriteLine("Optimal value:   0:F1", qp2.OptimalValue)

let x = qp2.Solve()
Console.WriteLine("Status: 0:F1", qp2.Status)
Console.WriteLine("Solution: 0:F1", x)
Console.WriteLine("Optimal value:   0:F1", qp2.OptimalValue)