6 Conic Modeling

6.1 The model

A model built using Fusion is always a conic optimization problem and it is convex by definition. These problems can be succinctly characterized as

(1)\[\begin{split}\begin{array}{lll} \minimize &c^T x&\\ \st & Ax+b \in \K & \end{array}\end{split}\]

where \(\K\) is a product of the following basic types of cones:

  • linear: \(\real^n\), \(\{ x \in \R^n ~:~ x \leq u\}\), \(\{ x \in \R^n ~:~ l\leq x \}\), \(\{ x \in \R^n ~:~ l\leq x \leq u\}\), a point,
  • quadratic: \(\Q^n = \{x\in\real^n~:~x_1\geq\sqrt{x_2^2+\cdots+x_n^2}\}\),
  • rotated quadratic: \(\Q_r^n = \{x\in\real^n~:~2x_1x_2\geq x_3^2+\cdots+x_n^2,\ x_1,x_2\geq 0\}\),
  • semidefinite:: \(\PSD^n=\{X\in\real^{n\times n}~:~ X\ \textrm{is symmmetric positive semidefinite}\}\).

The main thing about a Fusion model is that it can be specified in a convenient way without explicitly constructing the representation (1). Instead the user has access to variables which are used to construct linear operators that appear in constraints. The cone types described above are the domains of those constraints. A Fusion model can potentially contain many different building blocks of that kind. To facilitate manipulations with a large number of variables Fusion defines various logical views of parts of the model.

This section briefly summarizes the constructions and techniques available in Fusion. See Section 7 for a basic tutorial and Section 10 for more advanced case studies. This section is only an introduction: detailed specification of the methods and classes mentioned here can be found in the API reference.

A Fusion model is represented by the class Model and created by a simple construction

            using (Model M = new Model()) 
            { 

The model object is the user’s interface to the optimization problem, used in particular for

  • formulating the problem by defining variables, constraints and objective,
  • solving the problem and retrieving the solution status and solutions,
  • interacting with the solver: setting up parameters, registering for callbacks, performing I/O, obtaining detailed information from the optimizer etc.
  • memory management.

Almost all elements of the model: variables, constraints and the model itself can be constructed with or without names. If used, the names for each type of object must be unique. Choosing a good naming convention can make the problem more readable when dumped to a file. Most Fusion components also support some degree of pretty printing (toString method).

6.2 Variables

Continuous variables can be scalars, vectors or higher-dimensional arrays. The are added to the model with the method Model.Variable which returns a representing object of type Variable. The shape of a variable (number of dimensions and length in each dimension) has to be specified at creation. Optionally a variable may be created in a restricted domain (by default variables are unbounded, that is in \(\real\)). For instance, to declare a variable \(x\in\real_+^n\) we could write

              Variable x = M.Variable("x",n, Domain.GreaterThan(0.0));

A multi-dimensional variable is declared by specifying an array with all dimension sizes. Here is an \(n\times n\) variable:

              Variable x = M.Variable( new int[]{n,n}, Domain.Unbounded() ) ;

The specification of dimensions can also be part of the domain, as in this declaration of a symmetric positive semidefinite variable of dimension \(n\):

              Variable v = M.Variable(Domain.InPSDCone(n));

Integer variables are specified with an additional domain modifier. To add an integer variable \(z\in [1,10]\) we write

              Variable z = M.Variable("z", Domain.Integral(Domain.InRange(1.0,10.0)) );

The function Domain.Binary is a shorthand for binary variables often appearing in combinatorial problems:

              Variable y = M.Variable("y", Domain.Binary());

Integrality requirement can be switched on and off using the methods Variable.MakeInteger and Variable.MakeContinuous.

The Variable object provides the primal (Variable.Level) and dual (Variable.Dual) solution values of the variable after optimization, and it enters in the construction of linear expressions involving the variable.

6.3 Linear algebra

Linear expressions are constructed combining variables and matrices by linear operators. The result is an object that represents the linear expression itself. Fusion only allows for those combinations of operators and arguments that yield linear functions of the variables. Expressions have shapes and dimensions in the same fashion as variables. For instance, if \(x\in \real^n\) and \(A \in \real^{m\times n}\), then \(Ax\) is a vector expression of length \(m\). Note, however, that the internal size of \(Ax\) is \(mn\), because each entry is a linear combination for which \(m\) coefficients have to be stored.

Expressions are concrete implementations of the virtual interface Expression. In typical situations, however, all operations on expressions can be performed using the static methods and factory methods of the class Expr.

Table 6.1 Linear Operators
Method Description
Expr.Add Element-wise addition of two matrices
Expr.Sub Element-wise subtraction of two matrices
Expr.Mul Matrix or matrix-scalar multiplication
Expr.Neg Sign inversion
Expr.Outer Vector outer-product
Expr.Dot Dot product
Expr.Sum Sum over a given dimension
Expr.MulDiag Sum over the diagonal of a matrix which is the result of a matrix multiplication
Expr.ConstTerm Return a constant term

Operations on expressions must adhere to the rules of matrix algebra regarding dimensions; otherwise a DimensionError exception will be thrown.

Expression can be composed, nested and used as building blocks in new expressions. For instance \(Ax + By\) can be implemented as:

              Expr.Add( Expr.Mul(A,x), Expr.Mul(B,y) );

For operations involving multiple variables and expressions the users should consider list-based methods. For instance, a clean way to write \(x+y+z+w\) would be:

              Expr.Add( new Variable[]{x, y, z, w});

Note that a single variable (object of class Variable) can also be used as an expression. Once constructed, expressions are immutable.

6.4 Constraints and objective

Constraints are declared within an optimization model using the method Model.Constraint. Every constraint in Fusion has the form

Expression belongs to a Domain.

Objects of type Domain correspond roughly to the types of convex cones \(\mathcal{K}\) mentioned at the beginning of this section. For instance, the following set of linear constraints

(2)\[\begin{split}\begin{array}{lccccr} x_1 & + & 2 x_2 & & &= 0\\ & + & x_2 & + & x_3 &= 0\\ x_1 & & & & &= 0 \end{array}\end{split}\]

could be declared as

              double [,]A = new double[,] { { 1.0, 2.0, 0.0}, 
                                             { 0.0, 1.0, 1.0},
                                             { 1.0, 0.0, 0.0} };

              Variable x = M.Variable("x",3,Domain.Unbounded());
              Constraint c = M.Constraint( Expr.Mul(A,x), Domain.EqualsTo(0.0));

Note that the scalar domain Domain.EqualsTo consisting of a single point \(0\) scales up to the dimension of the expression and applies to all its elements. This allows many constraints to be comfortably expressed in a vectorized form. See also Section 6.7.

The Constraint object provides the dual (Constraint.Dual) value of the constraint after optimization and the primal value of the constraint expression (Constraint.Level).

The typical domains used to specify constraints are listed below. Note that they can also be used directly at variable creation, whenever that makes sense.

  Type Domain
Linear equality Domain.EqualsTo
inequality \(\leq\) Domain.LessThan
inequality \(\geq\) Domain.GreaterThan
two-sided bound Domain.InRange
Conic Quadratic quadratic cone Domain.InQCone
rotated quadratic cone Domain.InRotatedQCone
Semidefinite PSD matrix Domain.InPSDCone
Integral Integers in domain D Domain.Integral(D)
\(\{0,1\}\) Domain.Binary

Having discussed variables and constraints we can finish by defining the optimization objective with Model.Objective. The objective function is a scalar expression and the objective sense is specified by the enumeration ObjectiveSense as either minimize or maximize. The typical linear objective function \(c^T x\) can be declared as

              M.Objective( ObjectiveSense.Minimize,  Expr.Mul(c,x) );         

6.5 Matrices

At some point it becomes necessary to specify linear expressions such as \(Ax\) where \(A\) is a (large) constant data matrix. Such coefficient matrices can be represented in dense or sparse format. Dense matrices can always be represented using the standard data structures for arrays and two-dimensional arrays built into the language. Alternatively, or when sparsity can be exploited, matrices can be constructed as objects of the class Matrix. This can have some advantages: a more generic code that can be ported across platforms and can be used with both dense and sparse matrices without modifications.

Dense matrices are constructed with a variant of the static factory method Matrix.Dense. The values of all entries must be specified all at once and the resulting matrix is immutable. For example the matrix

\[\begin{split}A = \left[ \begin{array}{cccc} 1 & 2 & 3 & 4 \\ 5& 6 & 7 & 8 \end{array} \right]\end{split}\]

can be defined with:

              double[,] A= new double[,]{ {1.0,2.0,3.0,4.0}, {5.0,6.0,7.0,8.0} };
              Matrix Ad= Matrix.Dense(A);

or from a flattened representation:

              double[] Af={ 1,2,3,4,5,6,7,8 };
              Matrix Aff= Matrix.Dense(2,4,Af);

Sparse matrices are constructed with a variant of the static factory method Matrix.Sparse. This is both speed- and memory-efficient when the matrix has few nonzero entries. A matrix \(A\) in sparse format is given by a list of triples \((i, j, v)\), each defining one entry: \(A_{i,j}=v\). The order does not matter. The entries not in the list are assumed to be \(0\). For example, take the matrix

\[\begin{split}A = \left[ \begin{array}{cccc} 1.0 & 0.0 & 0.0 & 2.0\\ 0.0 & 3.0 & 0.0 & 4.0 \end{array} \right].\end{split}\]

Assuming we number rows and columns from \(0\), the corresponding list of triplets is:

\[A=\{(0,0,1.0),(0,3,2.0),(1,1,3.0),(1,3,4.0)\}\]

The Fusion definition would be:

              int[]    rows =  { 0,   0,   1,   1 };
              int[]    cols =  { 0,   3,   1,   3 };
              double[] values= { 1.0, 2.0, 3.0, 4.0 };
             
              Matrix ms = Matrix.Sparse(rows.Length, cols.Length, rows, cols, values);

The Matrix class provides more standard constructions such as the identity matrix, a constant value matrix, block diagonal matrices etc.

6.6 Stacking and views

Fusion provides a way to construct logical views of parts of existing expressions or combinations of existing expressions. They are still represented by objects of type Variable or Expression that refer to the original ones. This can be useful in some scenarios:

  • retrieving only the values of a few variables, and ignoring the remaining auxiliary ones,
  • stacking vectors or matrices to perform various matrix operations,
  • bundling a number of similar constraints into one (see Vectorization),
  • adding constraints between parts of the same variable, etc.

All these operations do not require new variables or expressions, but just lightweight logical views. In what follows we will concentrate on expressions; the same techniques are available for variables. These techniques will be familiar to the users of numerical tools such as Matlab or NumPy.

Picking and slicing

Expression.Pick picks a subset of entries from a variable or expression. Special cases of picking are Expression.Index, which picks just one scalar entry and Expression.Slice which picks a slice, that is restricts each dimension to a subinterval. Slicing is a frequently used operation.

_images/item-mapping-2.png

Fig. 6.1 Two dimensional slicing.

Both displayed regions are slices of the two-dimensional \(4\times 4\) expression, which can be selected as follows:

              Expression Axs1 = Ax.Slice(new int[]{0,0}, new int[]{2,2});
              Expression Axs2 = Ax.Index(new int[]{3,3});

Reshaping

Expressions can be reshaped creating a view with the same number of coordinates arranged in a different way. A particular example of this operation if flattening, which converts any multi-dimensional expression into a one-dimensional vector.

Stacking

Stacking refers to the concatenation of expressions to form a new larger one. For example, the next figure depicts the vertical stacking of two vectors of shape \(1\times 3\) resulting in a matrix of shape \(2\times 3\).

_images/stacking-2.png
              Expression c = Expr.Vstack(new Expression[]{a, b});

Vertical stacking (Expr.Vstack) of expressions of shapes \(d_1\times d_2\) and \(d_1'\times d_2\) has shape \((d_1+d_1')\times d_2\). Similarly, horizontal stacking (Expr.Hstack) of expressions of shapes \(d_1\times d_2\) and \(d_1\times d_2'\) has shape \(d_1\times (d_2+d_2')\). Fusion supports also more general versions of stacking for multi-dimensional variables, as described in Expr.Stack. A special case of stacking is repetition (Expr.Repeat), equivalent to stacking copies of the same expression.

6.7 Vectorization

Using Fusion one can compactly express sequences of similar constraints. For example, if we want to express

\[Ax_i = b_i,\ i=1,\ldots,n\]

we can think of \(x_i\in\real^m, b_i\in\real^k\) as the columns of two matrices \(X=[x_1,\ldots,x_n]\in\real^{m\times n}\), \(B=[b_1,\ldots,b_n]\in\real^{k\times n}\), and write simply

\[AX-B = 0.\]
              Variable X = Var.Hstack( new Variable[]{ xi[0], xi[1], xi[2], xi[3] } );
              Expression B = Expr.Hstack( new Expression[]{ bi[0], bi[1], bi[2], bi[3] } );

              M.Constraint(Expr.Sub(Expr.Mul(A, X), B), Domain.EqualsTo(0.0));

In this example the domain Domain.EqualsTo scales to apply to all the entries of the expression.

Another powerful case of vectorization and scaling domains is the ability to define a sequence of conic constraints in one go. Suppose we want to find an upper bound on the 2-norm of a sequence of vectors, that is we want to express

\[t \geq \|y_i\|,\ i=1,\ldots,n\]

Suppose that the vectors \(y_i\) are arranged in the rows of a matrix \(Y\). Then we can simply write:

              Variable t = M.Variable();
              
              M.Constraint(Expr.Hstack(Var.Vrepeat(t, n), Y), Domain.InQCone());

Here, again, the conic domain Domain.InQCone is by default applied to each row of the matrix separately, yielding the desired constraints in a loop-free way (the \(i\)-th row is \((t,y_i)\)). The direction along which conic constraints are created within multi-dimensional expressions can be changed with Domain.Axis.

We recommend vectorizing the code whenever possible. It is not only more elegant and portable but also more efficient — loops are eliminated and the number of Fusion API calls is reduced.

6.8 Reoptimization

Between optimizations the user can modify the model in two ways:

  • Add new constraints with Model.Constraint. This is useful for solving a sequence of optimization problems with more and more restrictions on the feasible set. See for example Section 10.9.

  • Replace the objective with a new one. This is particularly useful when solving a sequence of problems with the same data but different objectives, for instance in multi-objective optimization. For simplicity, suppose we want to minimize \(f(x) = \gamma x + \beta y\), for varying choices of \(\gamma>0\). Then we could write:

                  double[] gamma= new double[] {0.0, 0.5, 1.0};     // Choices for gamma
                  double beta = 2.0;
                  Variable x= M.Variable("x", 1, Domain.GreaterThan(0.0));
                  Variable y= M.Variable("y", 1, Domain.GreaterThan(0.0));
                  Expression beta_y = Expr.Mul(beta,y);
    
                  for(int i=0;i<3;i++)
                  {
                    M.Objective( ObjectiveSense.Minimize, Expr.Add(Expr.Mul(gamma[i],x), beta_y) );
                    M.Solve();
                  }
    
  • Add a new expression to an existing constraint (Constraint.Add).

Otherwise all Fusion objects are immutable.