7. A guided tour


7.1. Introduction

One of the big advantages of MATLAB is that it makes it very easy to do experiments and try out things without doing a lot of programming. The MOSEK optimization toolbox has been designed with this in mind. Hence, it should be very easy to solve optimization problems using MOSEK.

Moreover, a guided tour to the optimization toolbox has been designed to introduce the toolbox by examples. After having studyied these examples, the reader should be able to solve his or her own optimization problems without much further effort. Nevertheless, for the user interested in exploiting the toolbox to the limits, a detailed discussion and command reference are provided in the following chapters.

7.2. The tour starts

The MOSEK optimization toolbox consists of two layers of functions. The procedures in the top layer are application specific functions which have an easy-to-use interface. Currently, there are five procedures in the top layer:

msklpopt

Performs linear optimization.

mskqpopt

Performs quadratic optimization.

mskenopt

Performs entropy optimization.

mskgpopt

Performs geometric optimization (posynomial case).

mskscopt

Performs separable convex optimization.

The bottom layer of the MOSEK optimization toolbox consists of one procedure named mosekopt. This procedure provides a very flexible and powerful interface to the MOSEK optimization package. However, the price for this flexibility is a more complicated calling procedure.

For compatibility with the MATLAB optimization toolbox MOSEK also provides an implementation of linprog, quadprog and so forth. For details about these functions we refer the reader to Chapter 8.

In the following sections usage of the MOSEK optimization toolbox is demonstrated using examples. Most of these examples are available in

mosek\6\toolbox\examp\

7.3. The MOSEK terminolgy

First, some MOSEK terminology is introduced which will make the following sections easy to understand.

The MOSEK optimization toolbox can solve different classes of optimization problems such as linear, quadratic, conic, and mixed-integer optimization problems. Each of these problems is solved by one of the optimizers in MOSEK. Indeed MOSEK includes the following optimizers:

Depending on the optimizer different solution types may be produced, e.g. the interior-point optimizers produce a general interior-point solution whereas the simplex optimizer produces a basic solution.

7.4. Linear optimization

The first example is the linear optimization problem

\begin{math}\nonumber{}\begin{array}{lccccl}\nonumber{}\mbox{minimize} &  &  & x_{1}+2x_{2} &  & \\\nonumber{}\mbox{subject to} & 4 & \leq{} & x_{1}+x_{3} & \leq{} & 6,\\\nonumber{} & 1 & \leq{} & x_{1}+x_{2}, &  & \\\nonumber{} &  &  & 0\leq{}x_{1},x_{2},x_{3.} &  &\end{array}\end{math} (7.4.1)

7.4.1. Using msklpopt

A linear optimization problem such as (7.4.1) can be solved using the msklpopt function which is designed for solving the problem

\begin{math}\nonumber{}\begin{array}{llcccc}\nonumber{}\mbox{minimize} &  &  & c^{T}x &  & \\\nonumber{}\mbox{subject to} & l^{c} & \leq{} & Ax & \leq{} & u^{c},\\\nonumber{} & l^{x} & \leq{} & x & \leq{} & u^{x}.\end{array}\end{math} (7.4.2)

[[MathCmd 5]] and [[MathCmd 6]] are called constraint bounds whereas [[MathCmd 7]] and [[MathCmd 8]] are variable bounds.

The first step in solving the example (7.4.1) is to setup the data for problem (7.4.2) i.e. the c, A, etc. Afterwards the problem is solved using an appropriate call to msklpopt.

% lo1.m 

c     = [1 2 0]';
a     = [[1 0 1];[1 1 0]];
blc   = [4 1]';
buc   = [6 inf]';
blx   = sparse(3,1);
bux   = [];
[res] = msklpopt(c,a,blc,buc,blx,bux);
sol   = res.sol;

% Interior-point solution.

sol.itr.xx'      % x solution.
sol.itr.sux'     % Dual variables corresponding to buc.
sol.itr.slx'     % Dual variables corresponding to blx.

% Basic solution.

sol.bas.xx'      % x solution in basic solution.

Please note that

  • Infinite bounds are specified using -inf and inf. Moreover, the bux = [] means that all upper bounds [[MathCmd 8]] are plus infinite.
  • The [res] = msklpopt(c,a,blc,buc) call implies that the lower and upper bounds on x are minus and plus infinity respectively.
  • The lines after the msklpopt call can be omitted, but the purpose of those lines is to display different parts of the solutions. The res.sol field contains one or more solutions. In this case both the interior-point solution (sol.itr) and the basic solution (sol.bas) are defined.

7.4.2. Using mosekopt

The msklpopt function is in fact just a wrapper around the real optimization routine mosekopt. Therefore, an alternative to using the msklpopt is to call mosekopt directly. In general, the syntax for a mosekopt call is

[rcode,res] = mosekopt(cmd,prob,param) 

The arguments prob and param are optional. The purpose of the arguments are as follows:

cmd

A string telling mosekopt what to do, e.g. 'minimize info' tells mosekopt that the objective should be minimized and information about the optimization should be returned.

prob

A MATLAB structure specifying the problem that should be optimized.

param

A MATLAB structure specifying parameters controlling the behavior of the MOSEK optimizer. However, in general it should not be necessary to change the parameters.

The following MATLAB commands demonstrate how to set up the prob structure for the example (7.4.1) and solve the problem using mosekopt:

% lo2.m

clear prob;

% Specify the c vector.
prob.c = [ 1 2 0]';

% Specify a in sparse format.
subi   = [1 2 2 1];
subj   = [1 1 2 3];
valij  = [1.0 1.0 1.0 1.0];

prob.a = sparse(subi,subj,valij);

% Specify lower bounds of the constraints.
prob.blc  = [4.0 1.0]';

% Specify  upper bounds of the constraints.
prob.buc  = [6.0 inf]';

% Specify lower bounds of the variables.
prob.blx  = sparse(3,1);

% Specify upper bounds of the variables.
prob.bux = [];   % There are no bounds.

% Perform the optimization.
[r,res] = mosekopt('minimize',prob); 

% Show the optimal x solution.
res.sol.bas.xx

Please note that

  • A MATLAB structure named prob containing all the relevant problem data is defined.
  • All fields of this structure are optional except prob.a which is required to be a sparse matrix.
  • Different parts of the solution can be viewed by inspecting the solution field res.sol.

7.5. Convex quadratic optimization

A frequently occurring problem type is the quadratic optimization problem which consists of minimizing a quadratic objective function subject to linear constraints. One example of such a problem is:

\begin{math}\nonumber{}\begin{array}{lcccl}\nonumber{}\mbox{minimize} &  &  & x_{1}^{2}+0.1x_{2}^{2}+x_{3}^{2}-x_{1}x_{3}-x_{2} & \\\nonumber{}\mbox{subject to} & 1 & \leq{} & x_{1}+x_{2}+x_{3} & \\\nonumber{} &  &  & x\geq{}0. &\end{array}\end{math} (7.5.1)

In general, a quadratic optimization problem has the form

\begin{math}\nonumber{}\begin{array}{lccccl}\nonumber{}\mbox{minimize} &  &  & \frac{1}{2}x^{T}Qx+c^{T}x &  & \\\nonumber{}\mbox{subject to} & l^{c} & \leq{} & Ax, & \leq{} & u^{c},\\\nonumber{} & l^{x} & \leq{} & x & \leq{} & u^{x},\end{array}\end{math} (7.5.2)

which for the example (7.5.1) implies that

\begin{math}\nonumber{}Q=\left[\begin{array}{ccc}\nonumber{}2 & 0 & -1\\\nonumber{}0 & 0.2 & 0\\\nonumber{}-1 & 0 & 2\end{array}\right],\quad{}c=\left[\begin{array}{c}\nonumber{}0\\\nonumber{}-1\\\nonumber{}0\end{array}\right],\quad{}A=\left[\begin{array}{ccc}\nonumber{}1 & 1 & 1\end{array}\right],\end{math} (7.5.3)

and that

\begin{displaymath}\nonumber{}l^{c}=1,\quad{}u^{c}=\infty ,\quad{}l^{x}=\left[\begin{array}{c}\nonumber{}0\\\nonumber{}0\\\nonumber{}0\end{array}\right]\mbox{ and }u^{x}=\left[\begin{array}{c}\nonumber{}\infty \\\nonumber{}\infty\end{array}\right]\end{displaymath}

Please note the explicit [[MathCmd 14]] in the objective function of (7.5.2) which implies that diagonal elements must be doubled in Q, i.e. [[MathCmd 15]], whereas the coefficient in (7.5.1) is 1 in front of [[MathCmd 16]].

7.5.1. Two important assumptions

MOSEK assumes that the Q matrix is symmetric, i.e.

\begin{displaymath}\nonumber{}Q=Q^{T}\end{displaymath}

and that Q is positive semi-definite . A matrix is positive semi-definite if the smallest eigenvalue of the matrix is nonnegative. An alternative statement of the positive semi-definite requirement is

\begin{displaymath}\nonumber{}x^{T}Qx\geq{}0,~\forall x.\end{displaymath}

If Q is not positive semi-definite, then MOSEK will not produce reliable results or work at all.

One way of checking whether Q is positive semi-definite is to check whether all the eigenvalues of Q are nonnegative. The MATLAB command eig computes all eigenvalues of a matrix.

7.5.2. Using mskqpopt

The subsequent MATLAB statements solve the problem (7.5.1) using the mskqpopt MOSEK function

% qo1.m

% Set up Q.
q     = [[2 0 -1];[0 0.2 0];[-1 0 2]];

% Set up the linear part of the problem.
c     = [0 -1 0]';
a     = ones(1,3);
blc   = [1.0];
buc   = [inf];
blx   = sparse(3,1);
bux   = [];

% Optimize the problem.
[res] = mskqpopt(q,c,a,blc,buc,blx,bux);

% Show the primal solution.
res.sol.itr.xx

It should be clear that the format for calling mskqpopt is very similar to calling msklpopt except that the Q matrix is included as the first argument of the call. Similarly, the solution can be inspected by viewing the res.sol field.

7.5.3. Using mosekopt

The following sequence of MATLAB commands solves the quadratic optimization example by calling mosekopt directly.

% qo2.m

clear prob;

% c vector.
prob.c = [0 -1 0]';

% Define the data.

% First the lower triangular part of q in the objective 
% is specified in a sparse format. The format is:
%
%   Q(prob.qosubi(t),prob.qosubj(t)) = prob.qoval(t), t=1,...,4

prob.qosubi = [ 1  3 2   3]';
prob.qosubj = [ 1  1 2   3]';
prob.qoval  = [ 2 -1 0.2 2]';

% a, the constraint matrix
subi  = ones(3,1);
subj  = 1:3;
valij = ones(3,1);

prob.a = sparse(subi,subj,valij);

% Lower bounds of constraints.
prob.blc  = [1.0]';

% Upper bounds of constraints.
prob.buc  = [inf]';

% Lower bounds of variables.
prob.blx  = sparse(3,1);

% Upper bounds of variables.
prob.bux = [];   % There are no bounds.

[r,res] = mosekopt('minimize',prob);

% Display return code.
fprintf('Return code: %d\n',r);

% Display primal solution for the constraints.
res.sol.itr.xc'

% Display primal solution for the variables.
res.sol.itr.xx'

This sequence of commands looks much like the one that was used to solve the linear optimization example using mosekopt except that the definition of the Q matrix in prob. mosekopt requires that Q is specified in a sparse format. Indeed the vectors qosubi, qosubj, and qoval are used to specify the coefficients of Q in the objective using the principle

\begin{displaymath}\nonumber{}Q_{{\mbox{{qosubi(t)}},\mbox{{qosubj(t)}}}}=\mbox{{qoval(t)}},\quad{}\mbox{for}~t=1,\ldots ,\mbox{length}(\mbox{{qosubi}}).\end{displaymath}

An important observation is that due to Q being symmetric, only the lower triangular part of Q should be specified.

7.6. Conic optimization

One way of generalizing a linear optimization problem is to include a constraint of the form

\begin{displaymath}\nonumber{}x\in{}\mathcal{C}\end{displaymath}

in the problem definition where [[MathCmd 21]] is required to be a convex cone. The resulting class of problems is known as conic optimization.

MOSEK can solve a subset of all conic problems and subsequently it is demonstrated how to solve this subset using the mosekopt toolbox function.

7.6.1. The conic optimization problem

A conic optimization problem has the following form

\begin{math}\nonumber{}\begin{array}{lcccccl}\nonumber{}\mbox{minimize} &  &  & c^{T}x+c^{f} &  &  & \\\nonumber{}\mbox{subject to} & l^{c} & \leq{} & Ax & \leq{} & u^{c}, & \\\nonumber{} & l^{x} & \leq{} & x & \leq{} & u^{x}, & \\\nonumber{} &  &  & x\in{}\mathcal{C}, &  &  &\end{array}\end{math} (7.6.1)

where [[MathCmd 21]] must satisfy the following requirements. Let

\begin{displaymath}\nonumber{}x^{t}\in{}\mathbb{R}^{{n^{t}}},~t=1,\ldots ,k\end{displaymath}

be vectors comprised of parts of the decision variable vector x such that each decision variable is a member of exactly one [[MathCmd 25]] vector, e.g.:

\begin{displaymath}\nonumber{}x^{1}=\left[\begin{array}{c}\nonumber{}x_{1}\\\nonumber{}x_{4}\\\nonumber{}x_{7}\end{array}\right]\mbox{ and }x^{2}=\left[\begin{array}{c}\nonumber{}x_{6}\\\nonumber{}x_{5}\\\nonumber{}x_{{3}}\\\nonumber{}x_{2}\end{array}\right].\end{displaymath}

Next, define

\begin{displaymath}\nonumber{}\mathcal{C}:=\left\lbrace{}x\in{}\mathbb{R}^{n}:~x^{t}\in{}\mathcal{C}_{t},~t=1,2,\ldots ,k\right\rbrace{}\end{displaymath}

where [[MathCmd 28]] must have one of the following forms.

  • [[MathCmd 29]] set:

    \begin{displaymath}\nonumber{}\mathcal{C}_{t}=\lbrace{}x\in{}\mathbb{R}^{{n^{t}}}\rbrace{}.\end{displaymath}
  • Quadratic cone:

    \begin{displaymath}\nonumber{}\mathcal{C}_{t}=\left\lbrace{}x\in{}\mathbb{R}^{{n^{t}}}:x_{1}\geq{}\sqrt{\sum \limits _{{j=2}}^{{n^{t}}}x_{j}^{2}}\right\rbrace{}.\end{displaymath}
  • Rotated quadratic cone:

    \begin{displaymath}\nonumber{}\mathcal{C}_{t}=\left\lbrace{}x\in{}\mathbb{R}^{{n^{t}}}:2x_{1}x_{2}\geq{}\sum \limits _{{j=3}}^{{n^{t}}}x_{j}^{2},~x_{1},x_{2}\geq{}0\right\rbrace{}.\end{displaymath}

A variable is by default members of the [[MathCmd 29]] set unless it explicitly belongs to a specific cone.

Although the cones MOSEK can handle give rise to a limited class of conic problems it includes linear, quadratic, quadratically constrained optimization, and other classes of nonlinear convex optimization problems. See Section 10.3 for a discussion.

7.6.2. Solving an example

The problem

\begin{math}\nonumber{}\begin{array}{lccccc}\nonumber{}\mbox{minimize} & x_{5}+x_{6} &  & \\\nonumber{}\mbox{subject to} & x_{1}+x_{2}+x_{3}+x_{4} & = & 1,\\\nonumber{} & x_{1},x_{2},x_{3},x_{4} & \geq{} & 0,\\\nonumber{} & x_{5}\geq{}\sqrt{x_{1}^{2} + x_{3}^{2}}, &  & \\\nonumber{} & x_{6}\geq{}\sqrt{x_{2}^{2} + x_{4}^{2}} &  &\end{array}\end{math} (7.6.2)

is an example of a conic quadratic optimization problem. The problem involves some linear constraints and two quadratic cones. The linear constraints are specified as if the problem was a linear problem whereas the cones are specified using a MATLAB cell array named cones. cones must contain one cell per cone, where a cell must contain the two fields type and sub. type is used to specify the type of the cone and sub is used to specify the member variables of the cone.

The following MATLAB code demonstrates how to solve the example (7.6.2) using MOSEK.

% cqo1.m

clear prob;

% Specify the non-confic part of the problem.

prob.c   = [0 0 0 0 1 1];
prob.a   = sparse([1 1 1 1 0 0]);
prob.blc = 1;
prob.buc = 1;
prob.blx = [0 0 0 0 -inf -inf];
prob.bux = inf*ones(6,1);

% Specify the cones.
% Define an empty cell array names 'cones' containing one cell per cone.
 
prob.cones   = cell(2,1);

% The first cone is specified.

prob.cones{1}.type = 'MSK_CT_QUAD';
prob.cones{1}.sub  = [5 3 1];

% The second cone is specified.

prob.cones{2}.type = 'MSK_CT_QUAD';
prob.cones{2}.sub  = [6 2 4];

% The field 'type' specifies the cone type, i.e. whether it is quadratic cone
% or rotated quadratic cone. The keys for the two cone types are MSK_CT_QUAD
% MSK_CT_RQUAD respectively.
%
% The field 'sub' specifies the members of the cone, i.e. the above definitions
% imply that x(5) >= sqrt(x(3)^2+x(1)^2) and x(6) * x(2) >= x(4)^2.

% Optimize the problem. 

[r,res]=mosekopt('minimize',prob);

% Display the primal solution.

res.sol.itr.xx'

Note in partiucular that:

  • No variable can be member of more than one cone. This is not serious restriction — see the following section.
  • The [[MathCmd 29]] set is not specified explicitly.

7.6.3. Quadratic and conic optimization

The example

\begin{math}\nonumber{}\begin{array}{llcl}\nonumber{}\mbox{minimize} & x_{1}+x_{2}+x_{3} &  & \\\nonumber{}\mbox{subject to} & x_{1}^{2}+x_{2}^{2}+x_{3}^{2} & \leq{} & 1,\\\nonumber{} & x_{1}+0.5x_{2}^{2}+x_{3} & \leq{} & 0.5\end{array}\end{math} (7.6.3)

is not a conic quadratic optimization problem but can easily be reformulated as such.

Indeed the first constraint is equivalent to

\begin{math}\nonumber{}\begin{array}{rcl}\nonumber{}x_{4} & \geq{} & \sqrt{x_{1}^{2} + x_{2}^{2} + x_{3}^{2}},\\\nonumber{}x_{4} & = & 1\end{array}\end{math} (7.6.4)

where [[MathCmd 38]] is a new variable. This is a quadratic cone and a linear constraint. The second constraint in (7.6.3) is equivalent to

\begin{displaymath}\nonumber{}\begin{array}{rcl}\nonumber{}x_{1}+x_{3}+x_{5} & = & 0.5,\\\nonumber{}x_{2}-x_{7} & = & 0,\\\nonumber{}x_{5} & \geq{} & 0,\\\nonumber{}x_{6} & = & 1,\\\nonumber{}x_{7}^{2} & \leq{} & 2x_{5}x_{6},\end{array}\end{displaymath}

because this implies that

\begin{displaymath}\nonumber{}x_{5}\geq{}0.5x_{7}^{2}=0.5x_{2}^{2.}\end{displaymath}

and that

\begin{displaymath}\nonumber{}x_{1}+0.5x_{2}^{2}+x_{3}\leq{}x_{1}+x_{3}+x_{5}=0.5.\end{displaymath}

Please note that no variable can occur in more than one cone and therefore the additional constraint

\begin{displaymath}\nonumber{}x_{2}=x_{7}\end{displaymath}

is introduced and [[MathCmd 43]] is included in the second conic constraint instead of [[MathCmd 44]]. Using this “trick” it is always possible to obtain a formulation where no variable occurs in more than one cone.

Therefore, the example (7.6.3) is equivalent to the conic quadratic optimization problem

\begin{math}\nonumber{}\begin{array}{lccl}\nonumber{}\mbox{minimize} & x_{1}+x_{2}+x_{3} &  & \\\nonumber{}\mbox{subject to} & x_{1}+x_{3}+x_{5} & = & 0.5,\\\nonumber{} & x_{2}-x_{7} & = & 0,\\\nonumber{} & x_{4} & = & 1,\\\nonumber{} & x_{5} & \geq{} & 0,\\\nonumber{} & x_{6} & = & 1,\\\nonumber{} & x_{4}\geq{}\sqrt{x_{1}^{2} + x_{2}^{2} + x_{3}^{2}}, &  & \\\nonumber{} & 2x_{5}x_{6}\geq{}x_{7}^{2.} &  &\end{array}\end{math} (7.6.5)

This problem can be solved using MOSEK as follows:

% cqo2.m

% Set up the non-conic part of the problem.

prob          = [];
prob.c        = [1 1 1 0 0 0 0]';
prob.a        = sparse([[1 0 1 0 1 0 0];...
                        [0 1 0 0 0 0 -1]]);
prob.blc      = [0.5 0];
prob.buc      = [0.5 0];
prob.blx      = [-inf -inf -inf 1 -inf 1 -inf];
prob.bux      = [inf   inf  inf 1  inf 1  inf];

% Set up the cone information.

prob.cones         = cell(2,1);
prob.cones{1}.type = 'MSK_CT_QUAD';
prob.cones{1}.sub  = [4 1 2 3];
prob.cones{2}.type = 'MSK_CT_RQUAD';
prob.cones{2}.sub  = [5 6 7];

[r,res]       = mosekopt('minimize',prob);

% Display the solution.
res.sol.itr.xx'

7.6.4. Conic duality and the dual solution

The dual problem corresponding to the conic optimization problem (7.6.1) is given by

\begin{math}\nonumber{}\begin{array}{lcccccl}\nonumber{}\mbox{maximize} & (l^{c})^{T}s_{l}^{c}-(u^{c})^{T}s_{u}^{c} &  & \\\nonumber{} & +(l^{x})^{T}s_{l}^{x}-(u^{x})^{T}s_{u}^{x}+c^{f} &  & \\\nonumber{}\mbox{subject to} & -y+s_{l}^{c}-s_{u}^{c} & = & 0,\\\nonumber{} & A^{T}y+s_{l}^{x}-s_{u}^{x}+s_{n}^{x} & = & c,\\\nonumber{} & s_{l}^{c},s_{u}^{c},s_{l}^{x},s_{u}^{x} & \geq{} & 0,\\\nonumber{} & s_{n}^{x}\in{}\mathcal{C}^{*} &  &\end{array}\end{math} (7.6.6)

where the dual cone [[MathCmd 47]] is defined as follows. Let [[MathCmd 48]] be partitioned similar to x, i.e. if [[MathCmd 49]] is a member of [[MathCmd 25]], then [[MathCmd 51]] is a member of [[MathCmd 52]] as well. Now, the dual cone is defined by

\begin{displaymath}\nonumber{}\mathcal{C}^{*}:=\left\lbrace{}s_{n}^{x}\in{}\mathbb{R}^{{n^{t}}}:~(s_{n}^{x})^{t}\in{}\mathcal{C}^{*}_{t},~t=1,\ldots ,k\right\rbrace{}\end{displaymath}

where the type of [[MathCmd 54]] is dependent on the type of [[MathCmd 28]]. For the cone types MOSEK can handle the relation between the primal and dual cones is given as follows:

  • [[MathCmd 29]] set:

    \begin{displaymath}\nonumber{}\mathcal{C}_{t}=\left\lbrace{}x\in{}\mathbb{R}^{{n^{t}}}\right\rbrace{}\quad{}\Leftrightarrow \quad{}\mathcal{C}^{*}_{t}:=\left\lbrace{}s\in{}\mathbb{R}^{{n^{t}}}:~s=0\right\rbrace{}.\end{displaymath}
  • Quadratic cone:

    \begin{displaymath}\nonumber{}\mathcal{C}_{t}:=\left\lbrace{}x\in{}\mathbb{R}^{{n^{t}}}:x_{1}\geq{}\sqrt{\sum \limits _{{j=2}}^{{n^{t}}}x_{j}^{2}}\right\rbrace{}\quad{}\Leftrightarrow \quad{}\mathcal{C}^{*}_{t}=\mathcal{C}_{t}.\end{displaymath}
  • Rotated quadratic cone:

    \begin{displaymath}\nonumber{}\mathcal{C}_{t}:=\left\lbrace{}x\in{}\mathbb{R}^{{n^{t}}}:2x_{1}x_{2}\geq{}\sum \limits _{{j=3}}^{{n^{t}}}x_{j}^{2},~x_{1},x_{2}\geq{}0\right\rbrace{}\quad{}\Leftrightarrow \quad{}\mathcal{C}^{*}_{t}=\mathcal{C}_{t}.\end{displaymath}

For a more detailed discussion about conic duality see Section 10.3.

7.6.4.1. How to obtain the dual solution

When solving a conic optimization problem using MOSEK, the dual solution is available. The following MATLAB code fragment shows where the dual solution is stored.

% cqo3.m

[r,res]=mosekopt('minimize',prob);

% Solution record.
res.sol

% Dual variables for lower
% bounds of constraints.
res.sol.itr.slc'

% Dual variables for upper
% bounds of constraints.
res.sol.itr.suc'

% Dual variables for lower
% bounds on variables.
res.sol.itr.slx'

% Dual variables for upper
% bounds on variables.
res.sol.itr.sux'

% Dual variables with respect
% to the conic constraints.
res.sol.itr.snx'

7.6.5. Setting accuracy parameters for the conic optimizer

Three parameters control the accuracy of the solution obtained by the conic interior-point optimizer. The following example demonstrates which parameters should be reduced to obtain a more accurate solution, if required.

% How to change the parameters that controls
% the accuracy of a solution computed by the conic
% optimizer.

param = [];

% Primal feasibility tolerance for the primal solution
param.MSK_DPAR_INTPNT_CO_TOL_PFEAS = 1.0e-8;

% Dual feasibility tolerance for the dual solution
param.MSK_DPAR_INTPNT_CO_TOL_DFEAS = 1.0e-8;

% Relative primal-dual gap tolerance.
param.MSK_DPAR_INTPNT_CO_TOL_REL_GAP = 1.0e-8;

[r,res]=mosekopt('minimize',prob,param);

7.7. Quadratically constrained optimization

In the previous section a quadratically constrained optimization problem was solved using the conic optimizer. It is also possible to solve such a problem directly. An example of such an optimization problem is

\begin{math}\nonumber{}\begin{array}{llcl}\nonumber{}\mbox{minimize} & x_{1}+x_{2}+x_{3} &  & \\\nonumber{}\mbox{subject to} & x_{1}^{2}+x_{2}^{2}+x_{3}^{2} & \leq{} & 1,\\\nonumber{} & x_{1}+0.1x_{2}^{2}+x_{3} & \leq{} & 0.5.\end{array}\end{math} (7.7.1)

Please note that there are quadratic terms in both constraints. This problem can be solved using mosekopt as follows:

% qco1.m

clear prob;

% Specify the linear objective terms.
prob.c      = ones(3,1);

% Specify the quadratic objective terms.
prob.qcsubk = [1   1   1   2  ]';
prob.qcsubi = [1   2   3   2  ]';
prob.qcsubj = [1   2   3   2  ]';
prob.qcval  = [2.0 2.0 2.0 0.2]';

% Specify the linear constraint matrix
prob.a      = [sparse(1,3);sparse([1 0 1])];

prob.buc    = [1 0.5]';

[r,res]     = mosekopt('minimize',prob);

% Display the solution.
fprintf('\nx:');
fprintf(' %-.4e',res.sol.itr.xx');
fprintf('\n||x||: %-.4e',norm(res.sol.itr.xx));

Note that the quadratic terms in the constraints are specified using the fields prob.qcsubk, prob.qcsubi, prob.qcsubj, and prob.qcval as follows

\begin{displaymath}\nonumber{}Q_{{\mbox{{qcsubi(t)}},\mbox{{qcsubj(t)}}}}^{{\mbox{{qcsubk(t)}}}}=\mbox{{qcval(t)}},\quad{}\mbox{for}\quad{}t=1,\ldots ,\mbox{length}(\mbox{{qcsubk}})\end{displaymath}

where [[MathCmd 62]] is the quadratic term in the kth constraint. Also note that only the lower triangular part of the [[MathCmd 63]]s should be specified.

7.8. Linear least squares and related norm minimization problems

A frequently occurring problem in statistics and in many other areas of science is the problem

\begin{math}\nonumber{}\mbox{minimize}\quad{}\left\|Fx-b\right\|\end{math} (7.8.1)

where F and b are a matrix and vector of appropriate dimensions. x is the vector decision variables.

Typically, the norm used is the 1-norm, the 2-norm, or the infinity norm.

7.8.1. The case of the 2 norm

Initially let us focus on the 2 norm. In this case (7.8.1) is identical to the quadratic optimization problem

\begin{math}\nonumber{}\mbox{minimize}\quad{}1/2x^{T}F^{T}Fx+1/2b^{T}b-b^{T}Fx\end{math} (7.8.2)

in the sense that the set of optimal solutions for the two problems coincides. This fact follows from

\begin{displaymath}\nonumber{}\begin{array}{rcl}\nonumber{}\left\|Fx-b\right\|^{2} & = & (Fx-b)^{T}(Fx-b)\\\nonumber{} &  & x^{T}F^{T}Fx+b^{T}b+2b^{T}Fx.\end{array}\end{displaymath}

Subsequently, it is demonstrated how the quadratic optimization problem (7.8.2) is solved using mosekopt. In the example the problem data is read from a file, then data for the problem (7.8.2) is constructed and finally the problem is solved.

% nrm1.m

% Read data from 'afiro.mps'.
[r,res] = mosekopt('read(afiro.mps)');

% Get data for the problem
%             minimize ||f x - b||_2
f = res.prob.a';
b = res.prob.c;

% Solve the problem
%             minimize 0.5 xf'fx+0.5*b'*b-(f'*b)'*x

% Clear prob
clear prob;

% Compute the fixed term in the objective.
prob.cfix = 0.5*b'*b

% Create the linear objective terms
prob.c = -f'*b;

% Create the quadratic terms. Please note that only the lower triangular
% part of f'*f is used.
[prob.qosubi,prob.qosubj,prob.qoval] = find(sparse(tril(f'*f)))

% Obtain the matrix dimensions.
[m,n]   = size(f);

% Specify a.
prob.a  = sparse(0,n);

[r,res] = mosekopt('minimize',prob);

% The optimality conditions are f'*(f x - b) = 0.
% Check if they are satisfied:

fprintf('\nnorm(f^T(fx-b)): %e',norm(f'*(f*res.sol.itr.xx-b)));

Often the x variables must be within some bounds or satisfy some additional linear constraints. These requirements can easily be incorporated into the problem (7.8.2). E.g. the constraint [[MathCmd 67]] can be modeled as follows

% nrm2.m. Continuation of nrm1.m.

% Assume that the same objective should be
% minimized subject to -1 <= x <= 1

prob.blx = -ones(n,1);
prob.bux = ones(n,1);

[r,res] = mosekopt('minimize',prob);

% Check if the solution is feasible.
norm(res.sol.itr.xx,inf)

7.8.2. The case of the infinity norm

In some applications of the norm minimization problem (7.8.1) it is better to use the infinity norm than the 2 norm. However, the problem (7.8.1) stated as an infinity norm problem is equivalent to the linear optimization problem

\begin{math}\nonumber{}\begin{array}{lccl}\nonumber{}\mbox{minimize} & \tau  &  & \\\nonumber{}\mbox{subject to} & Fx+\tau e-b & \geq{} & 0,\\\nonumber{} & Fx-\tau e-b & \leq{} & 0,\end{array}\end{math} (7.8.3)

where e is the vector of ones of appropriate dimension. This implies that

\begin{displaymath}\nonumber{}\begin{array}{rcl}\nonumber{}\tau e & \geq{} & Fx-b\\\nonumber{}\tau e & \geq{} & -(Fx-b)\end{array}\end{displaymath}

and hence at optimum

\begin{displaymath}\nonumber{}\tau ^{*}=\left\|Fx^{*}-b\right\|_{\infty }\end{displaymath}

holds.

The problem (7.8.3) is straightforward to solve.

% nrm3.m. Continuation of nrm1.m.

% Let x(n+1) play the role as tau, then the problem is
% solved as follows.

clear prob;

prob.c   = sparse(n+1,1,1.0,n+1,1);
prob.a   = [[f,ones(m,1)];[f,-ones(m,1)]];
prob.blc = [b            ; -inf*ones(m,1)];
prob.buc = [inf*ones(m,1); b             ];

[r,res]  = mosekopt('minimize',prob);

% The optimal objective value is given by:
norm(f*res.sol.itr.xx(1:n)-b,inf) 

7.8.3. The case of the 1-norm

By definition, for the 1-norm we have that

\begin{displaymath}\nonumber{}\left\|Fx-b\right\|_{1}=\sum _{{i=1}}^{m}|f_{{i:}}x-b_{i}|.\end{displaymath}

Therefore, the norm minimization problem can be formulated as follows

\begin{math}\nonumber{}\begin{array}{lccll}\nonumber{}\mbox{minimize} & \sum \limits _{{i=1}}^{m}t_{i} &  &  & \\\nonumber{}\mbox{subject to} & |f_{{i:}}x-b_{i}| & = & t_{i}, & i=1,\ldots ,m,\end{array}\end{math} (7.8.4)

which in turn is equivalent to

\begin{math}\nonumber{}\begin{array}{lccll}\nonumber{}\mbox{minimize} & \sum \limits _{{i=1}}^{m}t_{i} &  &  & \\\nonumber{}\mbox{subject to} & f_{{i:}}x-b_{i} & \leq{} & t_{i}, & i=1,\ldots ,m,\\\nonumber{} & -(f_{{i:}}x-b_{i}) & \leq{} & t_{i}, & i=1,\ldots ,m.\end{array}\end{math} (7.8.5)

The reader should verify that this is really the case.

In matrix notation this problem can be expressed as follows

\begin{math}\nonumber{}\begin{array}{lccl}\nonumber{}\mbox{minimize} & e^{T}t &  & \\\nonumber{}\mbox{subject to} & Fx-te & \leq{} & b,\\\nonumber{} & Fx+te & \geq{} & b,\end{array}\end{math} (7.8.6)

where [[MathCmd 75]]. Next, this problem is solved.

% nrm4.m. Continuation of nrm1.m.

% Let x(n:(m+n)) play the role as t. Now,
% the problem can be solved as follows

clear prob;

prob.c   = [sparse(n,1)   ; ones(m,1)];
prob.a   = [[f,-speye(m)] ; [f,speye(m)]];
prob.blc = [-inf*ones(m,1); b];
prob.buc = [b             ; inf*ones(m,1)];

[r,res]  = mosekopt('minimize',prob);

% The optimal objective value is given by:
norm(f*res.sol.itr.xx(1:n)-b,1) 

7.8.3.1. A better formulation

It is possible to improve upon the formulation of the problem (7.8.5). Indeed problem (7.8.5) is equivalent to

\begin{math}\nonumber{}\begin{array}{lccll}\nonumber{}\mbox{minimize} & \sum \limits _{{i=1}}^{m}t_{i} &  &  & \\\nonumber{}\mbox{subject to} & f_{{i:}}x-b_{i}-t_{i}+v_{i} & = & 0, & i=1,\ldots ,m,\\\nonumber{} & -(f_{{i:}}x-b_{i})-t_{i} & \leq{} & 0, & i=1,\ldots ,m,\\\nonumber{} & v_{i}\geq{}0, &  &  & i=1,\ldots ,m.\end{array}\end{math} (7.8.7)

After eliminating the t variables then this problem is equivalent to

\begin{math}\nonumber{}\begin{array}{lccll}\nonumber{}\mbox{minimize} & \sum \limits _{{i=1}}^{m}(f_{{i:}}x-b_{i}+v_{i}) &  &  & \\\nonumber{}\mbox{subject to} & -2(f_{{i:}}x-b_{i})-v_{i} & \leq{} & 0, & i=1,\ldots ,m,\\\nonumber{} & v_{i}\geq{}0, &  &  & i=1,\ldots ,m.\end{array}\end{math} (7.8.8)

Please note that this problem has only half the number of general constraints than problem (7.8.5) since we have replaced constraints of the general form

\begin{displaymath}\nonumber{}f_{{i:}}x\leq{}b_{i}\end{displaymath}

with simpler constraints

\begin{displaymath}\nonumber{}v_{i}\geq{}0\end{displaymath}

which MOSEK treats in a special and highly efficient way. Furthermore MOSEK stores only the non-zeros in the coefficient matrix of the constraints. This implies that the problem (7.8.8) is likely to require much less space than the problem (7.8.7).

It is left as an exercise for the reader to implement this formulation in MATLAB.

7.9. More about solving linear least squares problems

Linear least squares problems with and without linear side constraints appear very frequently in practice and it is therefore important to know how such problems are solved efficiently using MOSEK.

Now, assume that the problem of interest is the linear least squares problem

\begin{math}\nonumber{}\begin{array}{lccccc}\nonumber{}\mbox{minimize} & \frac{1}{2}\left\|Fx-f\right\|_{2}^{2} &  & \\\nonumber{}\mbox{subject to} & Ax & = & b,\\\nonumber{} & l^{x}\leq{}x\leq{}u^{x}, &  &\end{array}\end{math} (7.9.1)

where F and A are matrices and the remaining quantities are vectors. x is the vector of decision variables. The problem (7.9.1) as stated is a convex quadratic optimization problem and can be solved as such.

However, if F has much fewer rows than columns then it will usually be more efficient to solve the equivalent problem

\begin{math}\nonumber{}\begin{array}{lccccc}\nonumber{}\mbox{minimize} & \frac{1}{2}\left\|z\right\|_{2}^{2} &  & \\\nonumber{}\mbox{subject to} & Ax & = & b,\\\nonumber{} & Fx-z & = & f,\\\nonumber{} & l^{x}\leq{}x\leq{}u^{x}. &  &\end{array}\end{math} (7.9.2)

Please note that a number of new constraints and variables has been introduced which of course seems to be disadvantageous but on the other hand the Hessian of the objective in problem (7.9.2) is much sparser than in problem (7.9.1). Frequently this turns out to be more important for the computational efficiency and therefore the latter formulation is usually the better one.

If F has many more rows than columns, then formulation (7.9.2) is not attractive whereas the corresponding dual problem is. Using the duality theory outlined in Section 10.4.1 we obtain the dual problem

\begin{math}\nonumber{}\begin{array}{lccccc}\nonumber{}\mbox{maximize} & b^{T}y+f^{T}\bar{y} &  & \\\nonumber{} & +(l^{x})^{T}s_{l}^{x}+(u^{x})^{T}s_{u}^{x} &  & \\\nonumber{} & -\frac{1}{2}\left\|z\right\|_{2}^{2} &  & \\\nonumber{}\mbox{subject to} & A^{T}y+F^{T}\bar{y}+s_{l}^{x}-s_{u}^{x} & = & 0,\\\nonumber{} & z-\bar{y} & = & 0,\\\nonumber{} & s_{l}^{x},s_{u}^{x}\geq{}0 &  &\end{array}\end{math} (7.9.3)

which can be simplified to

\begin{math}\nonumber{}\begin{array}{lccccc}\nonumber{}\mbox{maximize} & b^{T}y+f^{T}z &  & \\\nonumber{} & +(l^{x})^{T}s_{l}^{x}+(u^{x})^{T}s_{u}^{x} &  & \\\nonumber{} & -\frac{1}{2}\left\|z\right\|_{2}^{2} &  & \\\nonumber{}\mbox{subject to} & A^{T}y+F^{T}z+s_{l}^{x}-s_{u}^{x} & = & 0,\\\nonumber{} & s_{l}^{x},s_{u}^{x}\geq{}0 &  &\end{array}\end{math} (7.9.4)

after eliminating the [[MathCmd 84]] variables. Here we use the convention that

\begin{displaymath}\nonumber{}l_{j}^{x}=-\infty ~\Rightarrow (s_{l}^{x})_{j}=0\quad{}\mbox{and}\quad{}u_{j}^{x}=\infty \Rightarrow (s_{u}^{x})_{j}=0.\end{displaymath}

In practice such fixed variables in [[MathCmd 86]] and [[MathCmd 87]] should be removed from the problem.

Given our assumptions the dual problem (7.9.4) will have much fewer constraints than the primal problem (7.9.2); in general, the fewer constraints a problem contains, the more efficient MOSEK tends to be. A question is: If the dual problem (7.9.4) is solved instead of the primal problem (7.9.2), how is the optimal x solution obtained? It turns out that the dual variables corresponding to the constraint

\begin{displaymath}\nonumber{}A^{T}y+F^{T}z+s_{l}^{x}-s_{u}^{x}=0\end{displaymath}

are the optimal x solution. Therefore, due to the fact that MOSEK always reports this information as the

res.sol.itr.y

vector, the optimal x solution can easily be obtained.

In the following code fragment it is investigated whether it is attractive to solve the dual rather than the primal problem for a concrete numerical example. This example has no linear equalities and F is a 2000 by 400 matrix.

% nrm5.m

% Read data from a file.
[rcode,res] = mosekopt('read(lsqpd.mps) echo(0)');

% Define the problem data.
F           = res.prob.a;
f           = res.prob.blc;
blx         = res.prob.blx;
bux         = [];

% In this case there are no linear constraints
% First we solve the primal problem:
%
% minimize   0.5|| z ||^2
% subject to F x - z = f
%            l <= x <= u

% Note that m»n
[m,n]       = size(F);

prob        = [];

prob.qosubi = n+(1:m);
prob.qosubj = n+(1:m);
prob.qoval  = ones(m,1);
prob.a      = [F,-speye(m,m)];
prob.blc    = f;
prob.buc    = f;
prob.blx    = [blx;-inf*ones(m,1)];
prob.bux    = bux;


fprintf('m=%d  n=%d\n',m,n);

fprintf('First try\n');

tic
[rcode,res] = mosekopt('minimize echo(0)',prob);

% Display the solution time.
fprintf('Time           : %-.2f\n',toc);

try 
  % x solution:
  x = res.sol.itr.xx;

  % objective value:
  fprintf('Objective value: %-6e\n',norm(F*x(1:n)-f)^2);

  % Check feasibility.
  fprintf('Feasibility    : %-6e\n',min(x(1:n)-blx(1:n)));
catch
  fprintf('MSKERROR: Could not get solution')
end

% Clear prob.
prob=[];

%
% Next, we solve the dual problem.

% Index of lower bounds that are finite:
lfin        = find(blx>-inf);

% Index of upper bounds that are finite:
ufin        = find(bux<inf);

prob.qosubi = 1:m;
prob.qosubj = 1:m;
prob.qoval  = -ones(m,1);
prob.c      = [f;blx(lfin);-bux(ufin)];
prob.a      = [F',...
               sparse(lfin,(1:length(lfin))',...
                      ones(length(lfin),1),...
                      n,length(lfin)),...
               sparse(ufin,(1:length(ufin))',...
                      -ones(length(ufin),1),...
                      n,length(ufin))];
prob.blc    = sparse(n,1);
prob.buc    = sparse(n,1);
prob.blx    = [-inf*ones(m,1);...
               sparse(length(lfin)+length(ufin),1)];
prob.bux    = [];

fprintf('\n\nSecond try\n');
tic
[rcode,res] = mosekopt('maximize echo(0)',prob);

% Display the solution time.
fprintf('Time           : %-.2f\n',toc);

try
  % x solution:
  x = res.sol.itr.y;

  % objective value:
  fprintf('Objective value: %-6e\n',...
          norm(F*x(1:n)-f)^2);

  % Check feasibility.
  fprintf('Feasibility    : %-6e\n',...
          min(x(1:n)-blx(1:n)));
catch
  fprintf('MSKERROR: Could not get solution')
end

Here is the output produced:

m=2000  n=400
First try
Time           : 2.07
Objective value: 2.257945e+001
Feasibility    : 1.466434e-009


Second try
Time           : 0.47
Objective value: 2.257945e+001
Feasibility    : 2.379134e-009 

Both formulations produced a strictly feasible solution having the same objective value. Moreover, using the dual formulation leads to a reduction in the solution time by about a factor 5: In this case we can conclude that the dual formulation is far superior to the primal formulation of the problem.

7.9.1. Using conic optimization on linear least squares problems

Linear least squares problems can also be solved using conic optimization because the linear least squares problem

\begin{math}\nonumber{}\begin{array}{lccccc}\nonumber{}\mbox{minimize} & \left\|Fx-f\right\|_{2} &  & \\\nonumber{}\mbox{subject to} & Ax & = & b,\\\nonumber{} & l^{x}\leq{}x\leq{}u^{x} &  &\end{array}\end{math} (7.9.5)

is equivalent to

\begin{math}\nonumber{}\begin{array}{lccccc}\nonumber{}\mbox{minimize} &  &  & t &  & \\\nonumber{}\mbox{subject to} &  &  & Ax & = & b,\\\nonumber{} &  &  & Fx-z & = & f,\\\nonumber{} & l^{x} & \leq{} & x & \leq{} & u^{x},\\\nonumber{} &  &  & \left\|z\right\|_{2}\leq{}t. &  &\end{array}\end{math} (7.9.6)

This problem is a conic quadratic optimization problem having one quadratic cone and the corresponding dual problem is

\begin{math}\nonumber{}\begin{array}{lccl}\nonumber{}\mbox{maximize} & b^{T}y+f^{T}\bar{y}+(l^{x})^{T}s_{l}^{x}-(u^{x})^{T}s_{u}^{x} &  & \\\nonumber{}\mbox{subject to} & A^{T}y+F^{T}\bar{y}+s_{l}^{x}-s_{u}^{x} & = & 0,\\\nonumber{} & -\bar{y}+s_{z} & = & 0,\\\nonumber{} & s_{t} & = & 1,\\\nonumber{} & \left\|s_{z}\right\|\leq{}s_{t}, &  & \\\nonumber{} & s_{l}^{x},s_{u}^{x}\geq{}0 &  &\end{array}\end{math} (7.9.7)

which can be reduced to

\begin{math}\nonumber{}\begin{array}{lccl}\nonumber{}\mbox{maximize} & b^{T}y+f^{T}s_{z}+(l^{x})^{T}s_{l}^{x}-(u^{x})^{T}s_{u}^{x} &  & \\\nonumber{}\mbox{subject to} & A^{T}y-F^{T}\bar{s}_{z}+s_{l}^{x}-s_{u}^{x} & = & 0,\\\nonumber{} & s_{t} & = & 1,\\\nonumber{} & \left\|s_{z}\right\|\leq{}s_{t}, &  & \\\nonumber{} & s_{l}^{x},s_{u}^{x}\geq{}0. &  &\end{array}\end{math} (7.9.8)

Often the dual problem has much fewer constraints than the primal problem. In such cases it will be more efficient to solve the dual problem and obtain the primal solution x as the dual solution of the dual problem.

7.10. Entropy optimization

7.10.1. Using mskenopt

An entropy optimization problem has the following form

\begin{math}\nonumber{}\begin{array}{lccccl}\nonumber{}\mbox{minimize} &  &  & \sum \limits _{{j=1}}^{n}d_{j}x_{j}\ln (x_{j})+c^{T}x &  & \\\nonumber{}\mbox{subject to} & l^{c} & \leq{} & Ax & \leq{} & u^{c},\\\nonumber{} &  &  & 0\leq{}x, &  &\end{array}\end{math} (7.10.1)

where all the components of d must be nonnegative, i.e. [[MathCmd 94]]. An example of an entropy optimization problem is

\begin{math}\nonumber{}\begin{array}{lccccl}\nonumber{}\mbox{minimize} &  &  & x_{1}\ln (x_{1})-x_{1}+x_{2}\ln (x_{2}) &  & \\\nonumber{}\mbox{subject to} & 1 & \leq{} & x_{1}+x_{2} & \leq{} & 1,\\\nonumber{} &  &  & 0\leq{}x_{1},x_{2.} &  &\end{array}\end{math} (7.10.2)

This problem can be solved using the mskenopt command as follows

d     = [1 1]';
c     = [-1 0]';
a     = [1 1];
blc   = 1;
buc   = 1;
[res] = mskenopt(d,c,a,blc,buc);
res.sol.itr.xx;

7.11. Geometric optimization

A so-called geometric optimization problem can be stated as follows

\begin{math}\nonumber{}\begin{array}{lccll}\nonumber{}\mbox{minimize} & \sum \limits _{{k\in{}J_{0}}}c_{k}\prod \limits _{{j=1}}^{n}t_{j}^{{a_{{kj}}}} &  &  & \\\nonumber{}\mbox{subject to} & \sum \limits _{{k\in{}J_{i}}}c_{k}\prod \limits _{{j=1}}^{n}t_{j}^{{a_{{kj}}}} & \leq{} & 1, & i=1,\ldots ,m,\\\nonumber{} & t>0, &  &  &\end{array}\end{math} (7.11.1)

where it is assumed that

\begin{displaymath}\nonumber{}\cup _{{k=0}}^{m}J_{k}=\lbrace{}1,\ldots ,T\rbrace{}\end{displaymath}

and if [[MathCmd 98]], then

\begin{displaymath}\nonumber{}J_{i}\cap J_{j}=\emptyset .\end{displaymath}

Hence, A is a [[MathCmd 100]] matrix and c is a vector of length t. In general, the problem (7.11.1) is very hard to solve, but the posynomial case where

\begin{displaymath}\nonumber{}c>0\end{displaymath}

is relatively easy. Using the variable transformation

\begin{math}\nonumber{}t_{j}=e^{{x_{j}}}\end{math} (7.11.2)

we obtain the problem

\begin{math}\nonumber{}\begin{array}{lccll}\nonumber{}\mbox{minimize} & \sum \limits _{{k\in{}J_{0}}}c_{k}e^{{a_{{k:}}x}} &  &  & \\\nonumber{}\mbox{subject to} & \sum \limits _{{k\in{}J_{i}}}c_{k}e^{{a_{{k:}}x}} & \leq{} & 1, & i=1,\ldots ,m,\end{array}\end{math} (7.11.3)

which is convex in x for c>0. We apply the [[MathCmd 104]] function to obtain the equivalent problem

\begin{math}\nonumber{}\begin{array}{lccll}\nonumber{}\mbox{minimize} & \log (\sum \limits _{{k\in{}J_{0}}}c_{k}e^{{a_{{k:}}x}}) &  &  & \\\nonumber{}\mbox{subject to} & \log (\sum \limits _{{k\in{}J_{i}}}c_{k}e^{{a_{{k:}}x}}) & \leq{} & \log (1), & i=1,\ldots ,m,\end{array}\end{math} (7.11.4)

which is also a convex optimization problem since [[MathCmd 104]] is strictly increasing. Hence, the problem (7.11.4) can be solved by MOSEK.

For further details about geometric optimization we refer the reader to [17, pp. 531-538].

7.11.1. Using mskgpopt

MOSEK cannot handle a geometric optimization problem directly, but the transformation (7.11.4) can be solved using the MOSEK optimization toolbox function mskgpopt. Please note that the solution to the transformed problem can easily be converted into a solution to the original geometric optimization problem using relation (7.11.2).

Subsequently, we will use the example

\begin{math}\nonumber{}\begin{array}{lccccl}\nonumber{}\mbox{minimize} & 40t_{1}^{{-1}}t_{2}^{{-1/2}}t_{3}^{{-1}}+20t_{1}t_{3}+40t_{1}t_{2}t_{3} &  & \\\nonumber{}\mbox{subject to} & \frac{1}{3}t_{1}^{{-2}}t_{2}^{{-2}}+\frac{4}{3}t_{2}^{{1/2}}t_{3}^{{-1}} & \leq{} & 1,\\\nonumber{} & 0<t_{1},t_{2},t_{3} &  &\end{array}\end{math} (7.11.5)

to demonstrate how a geometric optimization problem is solved using mskgpopt. Please note that both the objective and the constraint functions consist of a sum of simple terms. These terms can be specified completely using the matrix

\begin{displaymath}\nonumber{}A=\left[\begin{array}{ccc}\nonumber{}-1 & -0.5 & -1\\\nonumber{}1 & 0 & 1\\\nonumber{}1 & 1 & 1\\\nonumber{}-2 & -2 & 0\\\nonumber{}0 & 0.5 & -1\end{array}\right],\end{displaymath}

and the vectors

\begin{displaymath}\nonumber{}c=\left[\begin{array}{c}\nonumber{}40\\\nonumber{}20\\\nonumber{}40\\\nonumber{}\frac{1}{3}\end{array}\right]~\mbox{and}~\mathrm{map}=\left[\begin{array}{c}\nonumber{}0\\\nonumber{}0\\\nonumber{}0\\\nonumber{}1\\\nonumber{}1\end{array}\right].\end{displaymath}

The interpretation is this: Each row of A,c describes one term, e.g. the first row of A and the first element of c describe the first term in the objective function. The vector [[MathCmd 110]] indicated whether a term belongs to the objective or to a constraint. If [[MathCmd 111]] equals zero, the kth term belongs to the objective function, otherwise it belongs to the [[MathCmd 111]]th constraint.

The following MATLAB code demonstrates how the example is solved using mskgpopt.

% go1.m

c     = [40 20 40 1/3 4/3]';
a     = sparse([[-1  -0.5  -1];[1 0 1];...
                [1 1 1];[-2 -2 0];[0 0.5 -1]]);
map   = [0 0 0 1 1]';
[res] = mskgpopt(c,a,map);

fprintf('\nPrimal optimal solution to original gp:');
fprintf(' %e',exp(res.sol.itr.xx));
fprintf('\n\n');

% Compute the optimal objective value and
% the constraint activities.
v = c.*exp(a*res.sol.itr.xx);

% Add appropriate terms together.
f = sparse(map+1,1:5,ones(size(map)))*v;

% First objective value. Then constraint values.
fprintf('Objective value: %e\n',log(f(1)));
fprintf('Constraint values:');
fprintf(' %e',log(f(2:end)));
fprintf('\n\n');

% Dual multipliers (should be negative)
fprintf('Dual variables (should be negative):');
fprintf(' %e',res.sol.itr.y);
fprintf('\n\n');

The code also computes the objective value and the constraint values at the optimal solution. Moreover, the optimal dual Lagrange multipliers for the constraints are shown and the gradient of the Lagrange function at the optimal point is computed.

7.11.2. Comments

7.11.2.1. Solving large scale problems

If you want to solve a large problem, i.e. a problem where A has large dimensions, then A must be sparse or you will run out of space. Recall that a sparse matrix contains few non-zero elements, so if A is a sparse matrix, you should construct it using MATLAB's sparse sparse as follows

A = sparse(subi,subj,valij);

where

\begin{displaymath}\nonumber{}a_{{\mathtt{subi}[k],\mathtt{subj}[k]}}=\mathtt{valij}[k].\end{displaymath}

For further details on the sparse function, please enter

help sparse

in MATLAB.

7.11.2.2. Preprocessing tip

Before solving a geometric optimization problem it is worthwhile to check if a column of the A matrix inputted to mskgpopt contains only positive elements. If this is the case, the corresponding variable [[MathCmd 114]] can take the value zero in the optimal solution: This may cause problems for MOSEK so it is better to remove such variables from the problem — doing so will have no influence on the optimal solution.

7.12. Separable convex optimization

This section discusses separable convex nonlinear optimization problems. A general separable nonlinear optimization problem can be specified as follows:

\begin{math}\nonumber{}\begin{array}{lccccc}\nonumber{}\mbox{minimize} &  &  & f(x)+c^{T}x &  & \\\nonumber{}\mbox{subject to } &  &  & g(x)+Ax-x^{c} & = & 0,\\\nonumber{} & l^{c} & \leq{} & x^{c} & \leq{} & u^{c},\\\nonumber{} & l^{x} & \leq{} & x & \leq{} & u^{x},\end{array}\end{math} (7.12.1)

where

This implies that the ith constraint essentially has the form

\begin{displaymath}\nonumber{}l_{i}^{c}\leq{}g_{i}(x)+\sum \limits _{{j=1}}^{n}a_{{ij}}x_{j}\leq{}u_{i}^{c}\end{displaymath}

when the [[MathCmd 127]] variable has been eliminated.

The problem (7.12.1) must satisfy the three important requirements:

  1. Separability: This requirement implies that all nonlinear functions can be written on the form

    \begin{displaymath}\nonumber{}f(x)=\sum \limits _{{j=1}}^{n}f^{j}(x_{j})\end{displaymath}

    and

    \begin{displaymath}\nonumber{}g_{i}(x)=\sum \limits _{{j=1}}^{n}g_{i}^{j}(x_{j}).\end{displaymath}

    Hence, the nonlinear functions can be written as a sum of functions which depends on only one variable.

  2. Differentiability: All functions should be twice differentiable for all [[MathCmd 49]] satisfying

    \begin{displaymath}\nonumber{}l_{j}^{x}<x<u_{j}^{x}\end{displaymath}

    if [[MathCmd 49]] occurs in at least one nonlinear function. Hence, if [[MathCmd 133]] appears in the problem, then the lower bound on [[MathCmd 44]] should be 0.

  3. Convexity: The problem should be a convex optimization problem. See Section 10.4 for a discussion of this requirement.

7.12.1. Using mskscopt

Subsequently, we will use the following example

\begin{math}\nonumber{}\begin{array}{lccl}\nonumber{}\mbox{minimize} & x_{1}-\ln (x_{1}+2x_{2}) &  & \\\nonumber{}\mbox{subject to} & x_{1}^{2}+x_{2}^{2} & \leq{} & 1\end{array}\end{math} (7.12.2)

to demonstrate solving a convex separable optimization problem using the MOSEK optimization toolbox function mskscopt.

First, note that the problem (7.12.2) is not a separable optimization problem due to the fact that the logarithmic term in objective is not a function of a single variable. However, by introducing one additional constraint and variable the problem can be made separable as follows

\begin{math}\nonumber{}\begin{array}{lccl}\nonumber{}\mbox{minimize} & x_{1}-\ln (x_{3}) &  & \\\nonumber{}\mbox{subject to} & x_{1}^{2}+x_{2}^{2} & \leq{} & 1,\\\nonumber{} & x_{1}+2x_{2}-x_{3} & = & 0,\\\nonumber{} & x_{3}\geq{}0. &  &\end{array}\end{math} (7.12.3)

This problem is separable and equivalent to the previous problem. Moreover, note that all nonlinear functions are well defined for values of x satisfying the variable bounds strictly, i.e.

\begin{displaymath}\nonumber{}x_{3}>0.\end{displaymath}

This (almost) makes sure that function evaluation errors will not occur during the optimization process since MOSEK will only evaluate [[MathCmd 138]] for [[MathCmd 139]].

When using the mskscopt function to solve problem (7.12.3), the linear part of the problem, such as a c and A, is specified as usual using MATLAB vectors and matrices. However, the nonlinear functions must be specified using five arrays which in the case of problem (7.12.3) can have the form

opr  = ['log'; 'pow'; 'pow'];
opri = [0;     1;     1    ];
oprj = [3;     1;     2    ];
oprf = [-1;    1;     1    ];
oprg = [0;     2;     2;   ];

Hence, opr(k,:) specifies the type of a nonlinear function, opri(k) specifies in which constraint the nonlinear function should be added (zero means objective), and oprj(k) means that the nonlinear function should be applied to [[MathCmd 49]]. Finally, oprf(k) and oprg(k) are parameters used by the mskscopt function according to the table:

opr(k,:) opri(k) oprj(k) oprf(k) oprg(k) function
ent i j f (not used) [[MathCmd 141]]
exp i j f g [[MathCmd 142]]
log i j f (not used) [[MathCmd 143]]
pow i j f g [[MathCmd 144]]

The i value indicates which constraint the nonlinear function belongs to. However, if i is identical to zero, then the function belongs to the objective. Using this notation a separable convex optimization problem can be solved with the function:

mskscopt(opr,
         opri,
         oprj,
         oprf,
         oprg,
         c,
         a,
         blc,
         buc,
         blx,
         bux)

All the elements for solving a separable convex nonlinear optimization problem have now been discussed and therefore we will conclude this section by showing the MATLAB code that will solve the example problem (7.12.3).

% sco1.m

% Specify the linear part of the problem.

c           = [1;0;0];
a           = sparse([[0 0 0];[1 2 -1]]);
blc         = [-inf; 0];
buc         = [1;0];
blx         = [-inf;-inf;0];

% Specify the nonlinear part.

opr         = ['log'; 'pow'; 'pow'];
opri        = [0;     1;     1    ];
oprj        = [3;     1;     2    ];
oprf        = [-1;    1;     1    ];
oprg        = [0;     2;     2    ];

% Call the optimizer.
% Note that bux is an optional parameter which should be added if the variables
% have an upper bound. 

[res]       = mskscopt(opr,opri,oprj,oprf,oprg,c,a,blc,buc,blx); 
                                                                 

% Print the solution.
res.sol.itr.xx

7.13. Mixed-integer optimization

Up until now it has been assumed that the variables in an optimization problem are continuous. Hence, it has been assumed that any value between the bounds on a variable is feasible. In many cases this is not a valid assumption because some variables are integer-constrained. E.g. a variable may denote the number of persons assigned to a given job and it may not be possible to assign a fractional person.

Using a mixed-integer optimizer MOSEK is capable of solving linear and quadratic optimization problems where one or more of the variables are integer-constrained.

7.13.1. Solving an example

Using the example

\begin{math}\nonumber{}\begin{array}{lccl}\nonumber{}\mbox{minimize} & -2x_{1}-3x_{2} &  & \\\nonumber{}\mbox{subject to } & 195x_{1}+273x_{2} & \leq{} & 1365,\\\nonumber{} & 4x_{1}+40x_{2} & \leq{} & 140,\\\nonumber{} & x_{1}\leq{}4, &  & \\\nonumber{} & x_{1},x_{2}\geq{}0, &  & \mbox{ and integer}\end{array}\end{math} (7.13.1)

we will demonstrate how to solve an integer optimization problem using MOSEK.

% milo1.m

% Specify the linear problem data as if
% the problem is a linear optimization
% problem.

clear prob       
prob.c        = [-2 -3];
prob.a        = sparse([[195 273];[4 40]]);
prob.blc      = -[inf inf];
prob.buc      = [1365 140];
prob.blx      = [0 0];
prob.bux      = [4 inf];

% Specify indexes of variables that are integer
% constrained.

prob.ints.sub = [1 2];

% Optimize the problem.
[r,res] = mosekopt('minimize',prob);

try 
  % Display the optimal solution.
  res.sol.int
  res.sol.int.xx'
catch
  fprintf('MSKERROR: Could not get solution')
end

Please note that compared to a linear optimization problem with no integer-constrained variables:

  • The prob.ints.sub field is used to specify the indexes of the variables that are integer-constrained.
  • The optimal integer solution is returned in the res.sol.int MATLAB structure.

7.13.2. Speeding up the solution of a mixed-integer problem

In general, a mixed-integer optimization problem can be very difficult to solve. Therefore, in some cases it may be necessary to improve upon the problem formulation and “assist” the mixed-integer optimizer.

How to obtain a good problem formulation is beyond the scope of this section and the reader is referred to [18]. However, two methods for assisting the mixed-integer optimizer are discussed subsequently.

7.13.2.1. Specifying an initial feasible solution

In many cases a good feasible integer solution to the optimization problem may be known. If this is the case, it is worthwhile to inform the mixed-integer optimizer since this will reduce the solution space searched by the optimizer.

Consider the problem:

\begin{math}\nonumber{}\begin{array}{lccl}\nonumber{}\mbox{maximize} & 7x_{0}+10x_{1}+x_{2}+5x_{3} &  & \\\nonumber{}\mbox{subject to} & x_{0}+x_{1}+x_{2}+x_{3} & \leq{} & 2.5\\\nonumber{} & x_{3}\geq{}0 &  & \\\nonumber{} & x_{0},x_{1},x_{2}\geq{}0 &  & \mbox{and integer},\end{array}\end{math} (7.13.2)

where only some of the variables are integer and the remaining are continuous. A feasible solution to this problem is:

\begin{math}\nonumber{}x_{0}=0,x_{1}=2,x_{2}=0,x_{3}=0.5\end{math} (7.13.3)

The following example demonstrates how to input this initial solution to MOSEK.

% milo2.m

clear prob
clear param
[r,res]         = mosekopt('symbcon');
sc              = res.symbcon;


prob.c          = [7 10 1 5];
prob.a          = sparse([1 1 1 1 ]);
prob.blc        = -[inf];
prob.buc        = [2.5];
prob.blx        = [0 0 0 0];
prob.bux        = [inf inf inf inf];
prob.ints.sub   = [1 2 3];

prob.sol.int.xx = [0 2 0 0.5]';

% Optionally set status keys too. 
% prob.sol.int.skx = [sc.MSK_SK_SUPBAS;sc.MSK_SK_SUPBAS;...
%                     sc.MSK_SK_SUPBAS;sc.MSK_SK_BAS] 
% prob.sol.int.skc = [sc.MSK_SK_UPR]

[r,res] = mosekopt('maximize',prob);

try
  % Display the optimal solution.
  res.sol.int.xx'
catch
  fprintf('MSKERROR: Could not get solution')
end

It is also possible to specify only the values of the integer variables and then let MOSEK compute values for the remaining continuous variables in order to obtain a feasible solution. If the MSK_IPAR_MIO_CONSTRUCT_SOL parameter is set to MSK_ON then MOSEK triesl to compute a feasible solution from the specified values of the integer variables. MOSEK generates the feasible solution by temporarily fixing all integer variables to the specified values and then optimizing the resulting continuous linear optimization problem. Hence, using this feature it is necessary to specify only the values of prob.sol.int.xx corresponding to the integer-constrained variables.

Suppose it is known that [[MathCmd 148]] are candidates for good integer values to our problem, then the following example demonstrates how to optimize the problem () using a feasible starting solution generated from the integer values as [[MathCmd 148]].

% milo3.m

[r,res]       = mosekopt('symbcon');
sc            = res.symbcon;

clear prob

prob.c        = [7 10 1 5];
prob.a        = sparse([1 1 1 1 ]);
prob.blc      = -[inf];
prob.buc      = [2.5];
prob.blx      = [0 0 0 0];
prob.bux      = [inf inf inf inf];
prob.ints.sub = [1 2 3];

% Values for the integer variables are specified.
prob.sol.int.xx  = [0 2 0 0]';

% Tell Mosek to construct a feasible solution from a given integer
% value. 
param.MSK_IPAR_MIO_CONSTRUCT_SOL = sc.MSK_ON;

[r,res] = mosekopt('maximize',prob,param);

try
  % Display the optimal solution.
  res.sol.int.xx'
catch
  fprintf('MSKERROR: Could not get solution')
end

7.13.2.2. Using branching priorities

The mixed-integer optimizer in MOSEK employs the so-called branch-and-bound algorithm to search for the optimal solution. See [18, pp. 91-112] for details about the branch-and-bound algorithm. The branch-and-bound algorithm can benefit from knowing about priorities of the integer variables.

E.g. in an optimization model some integer variables may denote which factories to build and other variables which products to make in the factories. It seems natural to decide upon which factories to build first and then decide upon which products to make in which factories. Hence, some integer variables are more important than others.

In MOSEK it is possible to assign priorities to all the integer variables. The higher priority assigned to a variable the more important the variable is considered by the branch-and-bound algorithm. Priorities are specified using the prob.ints.pri field as follows:

prob.ints.sub = [4 1 2 3];  % Integer variables.
prob.ints.pri = [5 10 2 4]; % Priorities.

This implies that variable 4 has priority 5, variable 1 has priority 10 and so forth.

An example of the usage of priorities can be seen in [18, pp. 232-235].

7.14. Sensitivity analysis

Given an optimization problem it is often useful to obtain information about how the optimal objective value changes when a problem parameter is perturbed. E.g. the objective function may reflect the price of a raw material such as oil which may not be known with certainty. Therefore, it is interesting to know how the optimal objective value changes as the oil price changes.

Analyzing how the optimal objective value changes when the problem data is changed is called sensitivity analysis.

Consider the problem:

minimize

\begin{math}\nonumber{}\begin{array}{ccccccccccccccl}\nonumber{}1x_{{11}} & + & 2x_{{12}} & + & 5x_{{23}} & + & 2x_{{24}} & + & 1x_{{31}} & + & 2x_{{33}} & + & 1x_{{34}}\end{array}\end{math} (7.14.1)

subject to

\begin{math}\nonumber{}\begin{array}{ccccccccccccccl}\nonumber{}x_{{11}} & + & x_{{12}} &  &  &  &  &  &  &  &  &  &  & \leq{} & 400,\\\nonumber{} &  &  &  & x_{{23}} & + & x_{{24}} &  &  &  &  &  &  & \leq{} & 1200,\\\nonumber{} &  &  &  &  &  &  &  & x_{{31}} & + & x_{{33}} & + & x_{{34}} & \leq{} & 1000,\\\nonumber{}x_{{11}} &  &  &  &  &  &  & + & x_{{31}} &  &  &  &  & = & 800,\\\nonumber{} &  & x_{{12}} &  &  &  &  &  &  &  &  &  &  & = & 100,\\\nonumber{} &  &  &  & x_{{23}} & + &  &  &  &  & x_{{33}} &  &  & = & 500,\\\nonumber{} &  &  &  &  &  & x_{{24}} & + &  &  &  &  & x_{{34}} & = & 500,\\\nonumber{}x_{{11}}, &  & x_{{12}}, &  & x_{{23}}, &  & x_{{24}}, &  & x_{{31}}, &  & x_{{33}}, &  & x_{{34}} & \geq{} & 0.\end{array}\end{math} (7.14.2)

The example below demonstrate how sensitivity analysis can answer questions of the type: What happens to the optimal solution if we decrease the upper bound of the first constraint with 1? For more information on sensitivity analysis see Chapter 14.

% sensitivity2.m

% Setup problem data.
clear prob
prob.a = sparse([1,    1,    0,    0,    0,    0,    0;
                 0,    0,    1,    1,    0,    0,    0;
                 0,    0,    0,    0,    1,    1,    1;
                 1,    0,    0,    0,    1,    0,    0;
                 0,    1,    0,    0,    0,    0,    0;
                 0,    0,    1,    0,    0,    1,    0;
                 0,    0,    0,    1,    0,    0,    1]);

prob.c =  [1,2,5,2,1,2,1];
prob.blc = [-Inf,-Inf,-Inf,800,100,500, 500];
prob.buc =[400,1200,1000,800,100,500,500];
prob.bux(1:7) = Inf; 
prob.blx(1:7) = 0;

% Analyze upper bound of constraint 1.
prob.prisen.cons.subu = [1];  

[r,res] = mosekopt('minimize echo(0)',prob); 
fprintf ('Optimal objective value: %e\n',prob.c * res.sol.bas.xx  );
fprintf('Sensitivity results for constraint 1:');
res.prisen.cons

% If we change the upper bound of constraint 1 with a
% value v in [res.prisen.cons.lr_bu(1),res.prisen.cons.rr_bu(1)] 
% then the optimal objective changes with - v * ls_bu(0) 
% e.g. changing prob.buc(1) with -1
prob.buc(1) =  prob.buc(1) - 1;
new_sol_predicted = prob.c * res.sol.bas.xx  + 1 * res.prisen.cons.ls_bu(1);
fprintf ('New optimal objective after changing bound predicted to:%e\n', ...
         new_sol_predicted);
[r,res] = mosekopt('minimize echo(0)',prob); 
fprintf ('New optimal objective value: %e\n',prob.c * res.sol.bas.xx  );

The output from running the example is given below:

Optimal objective value: 3.000000e+03
Sensitivity results for constraint 1:
ans =

    lr_bl: []
    rr_bl: []
    ls_bl: []
    rs_bl: []
    lr_bu: -300
    rr_bu: 0
    ls_bu: 3
    rs_bu: 3

New optimal objective after changing bound predicted to:3.003000e+03
New optimal objective value: 3.003000e+03

7.15. Inspecting a problem

The problem analyzer (discussed in detail in Sec. 13.1) provides useful diagnostics about an optimization problem, and is quick way to verify that a model has been specified correctly. For example, executing the command

mosekopt('anapro',prob)

will generate a report looking like

Constraints               Bounds                    Variables
 upper bd:        19       lower bd: all             cont: all
 fixed   :         8                                          

-------------------------------------------------------------------------------

Objective, cx
   range: min |c|: 0.00000   min |c|>0: 0.320000     max |c|: 10.0000
 distrib:        |c|        vars                                     
                   0          27                                     
           [0.32, 1)           4                                     
             [1, 10]           1                                     

-------------------------------------------------------------------------------

Constraint matrix A has
        27 rows (constraints)
        32 columns (variables)
        83 (9.60648%) nonzero entries (coefficients)

Row nonzeros, A_i
   range: min A_i: 1 (3.125%)    max A_i: 9 (28.125%)
 distrib:        A_i        rows       rows%        acc%
                   1           2        7.41        7.41
                   2          16       59.26       66.67
              [3, 7]           8       29.63       96.30
              [8, 9]           1        3.70      100.00

Column nonzeros, A|j
   range: min A|j: 1 (3.7037%)    max A|j: 4 (14.8148%)
 distrib:        A|j        cols       cols%        acc%
                   1           1        3.12        3.12
                   2          21       65.62       68.75
              [3, 4]          10       31.25      100.00

A nonzeros, A(ij)
   range: min |A(ij)|: 0.107000     max |A(ij)|: 2.42900
 distrib:      A(ij)      coeffs
          [0.107, 1)          17
           [1, 2.43]          66

-------------------------------------------------------------------------------

Constraint bounds, lb <= Ax <= ub
 distrib:        |b|             lbs             ubs
                   0               7              20
           [10, 100)               1               3
         [100, 1000]                               4

Variable bounds, lb <= x <= ub
 distrib:        |b|             lbs             ubs
                   0              32

The report provides an overview of the objective function, the number of constraints and bounds, as well as sparsity information and distributions of nonzero elements.

7.16. The solutions

Whenever an optimization problem is solved using MOSEK, one or more optimal solutions are reported depending on which optimizer is used. These solutions are available in the

res.sol

structure, which has one or more of the subfields

res.sol.itr  % Interior solution.
res.sol.bas  % Basic    solution.
res.sol.int  % Integer  solution.

The interior (point) solution is an arbitrary optimal solution which is computed using the interior-point optimizer. The basic solution is available only for linear problems and is produced by the simplex optimizer or the basis identification process which is an add-on to the interior-point optimizer. Finally, the integer solution is available only for problems having integer-constrained variables and is computed using the integer optimizer.

Each of the three solutions may contain one or more of the following subfields:

.prosta

Problem status. See Appendix F.38.

.solsta

Solution status. See Appendix F.51.

.skc

Constraint status keys. See Tablen 7.1 below.

.skx

Variable status keys. See Table 7.1 below.

.xc

Constraint activities.

.xx

Variable activities.

.y

Identical to -.slc+.suc.

.slc

Dual variables corresponding to lower constraint bounds.

.suc

Dual variables corresponding to upper constraint bounds.

.slx

Dual variables corresponding to lower variable bounds.

.sux

Dual variables corresponding to upper variable bounds.

.snx

Dual variables corresponding to the conic constraints.

7.16.1. The constraint and variable status keys

In a solution both constraints and variables are assigned a status key which indicates whether the constraint or variable is at its lower limit, its upper limit, is super basic and so forth in the optimal solution. For interior-point solutions these status keys are only indicators which the optimizer produces.

In Table 7.1 the possible values for the status keys are shown accompanied with an interpretation of the key.

Symbolic Numeric String Interpretation
constant constant code
MSK_SK_UNK 0 UN Unknown status
MSK_SK_BAS 1 BS Is basic
MSK_SK_SUPBAS 2 SB Is superbasic
MSK_SK_LOW 3 LL Is at the lower limit (bound)
MSK_SK_UPR 4 UL Is at the upper limit (bound)
MSK_SK_FIX 5 EQ Lower limit is identical to upper limit
MSK_SK_INF 6 ** Is infeasible i.e. the lower limit is
  greater than the upper limit.
Table 7.1: Constraint and variable status keys.

By default the constraint and variable status keys are reported using string codes but it is easy to have MOSEK report the numeric codes instead. Indeed in the example

% Status keys in string format.
[rcode,res]=mosekopt('minimize statuskeys(0)',prob);
res.sol.skc(1)
res.sol.prosta 

the status keys are represented using string codes whereas in the example

% Status keys in string format.
[rcode,res]=mosekopt('minimize statuskeys(1)',prob);
res.sol.skc(1)
res.sol.prosta 

the status keys are represented using numeric codes.

7.17. Viewing the task information

In MOSEK the optimization problem and the related instructions with respect to the optimization process are called an optimization task or for short a task. Whenever MOSEK performs operations on a task it stores information in the task information database. Examples of information that is stored are the number of interior-point iterations performed to solve the problem and time spent doing the optimization.

All the items stored in the task information database are listed in Appendixes F.13 and F.17. It is possible to see the whole or part of the task information database from within MATLAB.

% Solve a problem and obtain
% the task information database.
[r,res]=mosekopt('minimize info',prob);

% View one item
res.info.MSK_IINF_INTPNT_ITER

% View the whole database
res.info

7.18. Inspecting and setting parameters

A large number of parameters controls the behavior of MOSEK, e.g. there is a parameter controlling which optimizer is used, one that limits the maximum number of iterations allowed, and several parameters specifying the termination tolerance. All these parameters are stored in a database internally in MOSEK. The complete parameter database can be obtained and viewed using the commands:

[r,res]=mosekopt('param');
res.param 

We will not describe the purpose of each parameter here but instead refer the reader to Appendix E where all the parameters are presented in detail.

In general, it should not be necessary to change any of the parameters but if required, it is easily done. In the following example code it is demonstrated how to modify a few parameters and afterwards performing the optimization using these parameters.

% Obtain all symbolic constants
% defined by MOSEK.

[r,res]  = mosekopt('symbcon');
sc       = res.symbcon;

param    = [];

% Basis identification is unnecessary.
param.MSK_IPAR_INTPNT_BASIS   = sc.MSK_OFF;

% Alternatively you can use
%
%  param.MSK_IPAR_INTPNT_BASIS   = 'MSK_OFF';
%

% Use another termination tolerance.
param.MSK_DPAR_INTPNT_TOLRGAP = 1.0e-9;

% Perform optimization using the
% modified parameters.

[r,res] =  mosekopt('minimize',prob,param);

7.19. Advanced start (hot-start)

In practice it frequently occurs that when an optimization problem has been solved, then the same problem slightly modified should be reoptimized. Moreover, if it is just a small the modification, it can be expected that the optimal solution to the original problem is a good approximation to the modified problem. Therefore, it should be efficient to start the optimization of the modified problem from the previous optimal solution.

Currently, the interior-point optimizer in MOSEK cannot take advantage of a previous optimal solution, however, the simplex optimizer can exploit any basic solution.

7.19.1. Some examples using hot-start

Using the example

\begin{math}\nonumber{}\begin{array}{lccccl}\nonumber{}\mbox{minimize} &  &  & x_{1}+2x_{2} &  & \\\nonumber{}\mbox{subject to} & 4 & \leq{} & x_{1}+x_{3} & \leq{} & 6,\\\nonumber{} & 1 & \leq{} & x_{1}+x_{2}, &  & \\\nonumber{} &  &  & 0\leq{}x_{1},x_{2},x_{3} &  &\end{array}\end{math} (7.19.1)

the hot-start facility using the simplex optimizer will be demonstrated. A quick inspection of the problem indicates that [[MathCmd 153]] is an optimal solution. Hence, it seems to be a good idea to let the initial basis consist of [[MathCmd 154]] and [[MathCmd 155]] and all the other variables be at their lower bounds. This idea is used in the example code:

% advs1.m

clear prob param bas

% Specify an initial basic solution.
bas.skc      = ['LL';'LL'];
bas.skx      = ['BS';'LL';'BS'];
bas.xc       = [4 1]';
bas.xx       = [1 3 0]';

prob.sol.bas = bas;

% Specify the problem data.
prob.c       = [ 1 2 0]';
subi         = [1 2 2 1];
subj         = [1 1 2 3];
valij        = [1.0 1.0 1.0 1.0];
prob.a       = sparse(subi,subj,valij);
prob.blc     = [4.0 1.0]';
prob.buc     = [6.0 inf]';
prob.blx     = sparse(3,1);
prob.bux     = [];

% Use the primal simplex optimizer.
param.MSK_IPAR_OPTIMIZER = 'MSK_OPTIMIZER_PRIMAL_SIMPLEX';
[r,res] = mosekopt('minimize',prob,param)

Some comments:

  • In the example the dual solution is defined. This is acceptable because the primal simplex optimizer is used for the reoptimization and it does not exploit a dual solution. In the future MOSEK will also contain a dual simplex optimizer and if that optimizer is used, it will be important that a “good” dual solution is specified.
  • The status keys bas.skc and bas.skx must contain only the entries BS, EQ, LL, UL, and SB. Moreover, e.g. EQ must be specified only for a fixed constraint or variable. LL and UL can be used only for a variable that has a finite lower or upper bound respectively.
  • The number of constraints and variables defined to be basic must correspond exactly to the number of constraints, i.e. the row dimension of A.

7.19.2. Adding a new variable

Next, assume that the problem

\begin{math}\nonumber{}\begin{array}{lccccl}\nonumber{}\mbox{minimize} &  &  & x_{1}+2x_{2}-x_{4} &  & \\\nonumber{}\mbox{subject to} & 4 & \leq{} & x_{1}+x_{3}+x_{4} & \leq{} & 6,\\\nonumber{} & 1 & \leq{} & x_{1}+x_{2}, &  & \\\nonumber{} &  &  & 0\leq{}x_{1},x_{2},x_{3},x_{4.} &  &\end{array}\end{math} (7.19.2)

should be solved. It is identical to the problem (7.19.1) except that a new variable [[MathCmd 38]] has been added. In continuation of the previous example this problem can be solved as follows (using hot-start):

% advs2.m. Continuation of advs1.m.

prob.c       = [prob.c;-1.0];
prob.a       = [prob.a,sparse([1.0 0.0]')];
prob.blx     = sparse(4,1);

% Reuse the old optimal basic solution.
bas          = res.sol.bas;

% Add to the status key.
bas.skx      = [res.sol.bas.skx;'LL'];

% The new variable is at it lower bound.
bas.xx       = [res.sol.bas.xx;0.0];
bas.slx      = [res.sol.bas.slx;0.0];
bas.sux      = [res.sol.bas.sux;0.0];

prob.sol.bas = bas;

[rcode,res]  = mosekopt('minimize',prob,param);

% The new primal optimal solution
res.sol.bas.xx'

7.19.3. Fixing a variable

In e.g. branch-and-bound methods for integer programming problems it is necessary to reoptimize the problem after a variable has been fixed to a value. This can easily be achieved as follows:

% advs3.m. Continuation of advs2.m.

prob.blx(4)  = 1;
prob.bux     = [inf inf inf 1]';

% Reuse the basis.
prob.sol.bas = res.sol.bas;

[rcode,res]  = mosekopt('minimize',prob,param);

% Display the optimal solution.
res.sol.bas.xx'

The [[MathCmd 38]] variable is simply fixed at the value 1 and the problem is reoptimized. Please note that the basis from the previous optimization can immediately be reused.

7.19.4. Adding a new constraint

Now, assume that the constraint

\begin{math}\nonumber{}x_{1}+x_{2}\geq{}2\end{math} (7.19.3)

should be added to the problem and the problem should be reoptimized. The following example demonstrates how to do this.

% advs4.m. A continuation of advs3.m. 

% Modify the problem.
prob.a       = [prob.a;sparse([1.0 1.0 0.0 0.0])];
prob.blc     = [prob.blc;2.0];
prob.buc     = [prob.buc;inf];

% Obtain the previous optimal basis.
bas          = res.sol.bas;

% Set the solution to the modified problem.
bas.skc      = [bas.skc;'BS'];
bas.xc       = [bas.xc;bas.xx(1)+bas.xx(2)];
bas.y        = [bas.y;0.0];
bas.slc      = [bas.slc;0.0];
bas.suc      = [bas.suc;0.0];

% Reuse the basis.
prob.sol.bas = bas;

% Reoptimize.
[rcode,res]  = mosekopt('minimize',prob,param);

res.sol.bas.xx' 

Please note that the slack variable corresponding to the new constraint are declared basic. This implies that the new basis is nonsingular and can be reused.

7.19.5. Using numeric values to represent status key codes

In the previous examples the constraint and variable status keys are represented using string codes. Although the status keys are easy to read they are sometimes difficult to work with in a program. Therefore, the status keys can also be represented using numeric values as demonstrated in the example:

% sk1.m

% Obtain all symbolic constants
% defined in MOSEK.

clear prob bas;

[r,res]  = mosekopt('symbcon');
sc       = res.symbcon;

% Specify an initial basic solution.
% Please note that symbolic constants are used.
% I.e. sc.MSK_SK_LOW instead of 4.
bas.skc      = [sc.MSK_SK_LOW;sc.MSK_SK_LOW];
bas.skx      = [sc.MSK_SK_BAS;sc.MSK_SK_LOW;sc.MSK_SK_BAS];
bas.xc       = [4 1]';
bas.xx       = [1 3 0]';
prob.sol.bas = bas;

% Specify the problem data.
prob.c   = [ 1 2 0]';
subi     = [1 2 2 1];
subj     = [1 1 2 3];
valij    = [1.0 1.0 1.0 1.0];
prob.a   = sparse(subi,subj,valij);
prob.blc = [4.0 1.0]';
prob.buc = [6.0 inf]';
prob.blx = sparse(3,1);
prob.bux = [];

% Use the primal simplex optimizer.
clear param;
param.MSK_IPAR_OPTIMIZER = sc.MSK_OPTIMIZER_PRIMAL_SIMPLEX;

[r,res] = mosekopt('minimize statuskeys(1)',prob,param)

% Status keys will be numeric now i.e.

res.sol.bas.skc'

% is a vector of numeric values.

Please note that using the commands

[r,res]  = mosekopt('symbcon');
sc       = res.symbcon;

all the symbolic constants defined within MOSEK are obtained and used in the lines

bas.skc  = [sc.MSK_SK_LOW;sc.MSK_SK_LOW];
bas.skx  = [sc.MSK_SK_BAS;sc.MSK_SK_LOW;sc.MSK_SK_BAS]; 

These two lines are in fact equivalent to

bas.skc  = [1;1];
bas.skx  = [3;1;3];

However, it is not recommended to specify the constraint and variable status keys this way because it is less readable and portable. Indeed if e.g. MOSEK later changes the definition that 1 is equivalent to `LL', all programs using numerical keys will be incorrect whereas using the symbolic constants the programs remain correct.

7.20. Using names

In MOSEK it is possible to give the objective, each constraint, each variable, and each cone a name. In generalm such names are not really needed except in connection with reading and writing MPS files. See Section 7.21 for details.

All the names are specified in the prob.names structure.

% The problem is named.
prob.names.name   = 'CQO example';

% Objective name.
prob.names.obj    = 'cost';

% The two constraints are named.
prob.names.con{1} = 'constraint_1';
prob.names.con{2} = 'constraint_2';

% The six variables are named.
prob.names.var    = cell(6,1);
for j=1:6
  prob.names.var{j} = sprintf('x%d',j);
end

% Finally the two cones are named.
prob.names.cone{1} = 'cone_a';
prob.names.cone{2} = 'cone_b';

7.20.1. Blanks in names

Although it is allowed to use blanks (spaces) in names it is not recommended to do so except for the problem name. In general, avoid names like “x 1” or “con 1”.

7.21. MPS files

An industry standard format for storing linear optimization problems in an ASCII file is the so-called MPS format. For readers not familiar with the MPS format a specification of the MPS format supported by MOSEK can be seen in Appendix A.

The advantage of the MPS format is that problems stored in this format can be read by any commercial optimization software, so it facilitates communication of optimization problems.

7.21.1. Reading an MPS file

It is possible to use mosekopt to read an MPS file containing the problem data. In this case mosekopt reads data from an MPS file and returns both the problem data and the optimal solution, if required. Assume that afiro.mps is the MPS file from which mosekopt should read the problem data, then this task is performed using the command

[r,res] = mosekopt('read(afiro.mps'));

In this case res.prob will contain several fields with the problem data. E.g.

res.prob.c'

will display the c-vector.

The names used in the MPS file is also available in the prob.names structure.

% All names.
prob.names

% Constraint names.
prob.names.con

The quadratic terms of a problem can be accessed and displayed in a similar manner:

% mpsrd.m

% Read data from the file wp12-20.mps.

[r,res] = mosekopt('read(wp12-20.mps)');

% Looking at the problem data
prob = res.prob;
clear res;

% Form the quadratic term in the objective.
q = sparse(prob.qosubi,prob.qosubj,prob.qoval);

% Get a graphical picture.
spy(q) % Notice that only the lower triangular part is defined.

7.21.2. Writing a MPS files

It is possible to write an MPS file using MOSEK. To write a problem contained in a MATLAB structure prob to the file “datafile.mps”, use the command:

% Write the data defined by prob to an MPS file
% named datafile.mps
mosekopt('write(datafile.mps)',prob);

If the prob.names field is defined, MOSEK will use those names when writing the MPS file, otherwise MOSEK will use generic (automatically generated) names.

7.22. User call-back functions

A call-back function is a user-defined MATLAB function to be called by MOSEK on a given event. The optimization toolbox supports two types of call-back functions which are presented below.

7.22.1. Log printing via call-back function

When using mosekopt it is possible to control the amount of information that mosekopt prints to the screen, e.g.

[r,res] = mosekopt('minimize echo(0)',prob)

forces mosekopt to not print log information — the string echo(0) indicates that no output should be printed during optimization. A high value in the echo(n) command, e.g. echo(3), forces MOSEK to display more log information.

It is possible to redirect the MOSEK log printing almost anywhere using a user-defined log call-back function. It works as follows. Create an m-file to handle the log output, similar to:

function myprint(handle,str)
% handle: Is user defined data structure
% str   : Is a log string.
%

fprintf(handle,'%s',str); 

The name and actions of the function are not important, but its argument list must be identical to the example: It must accept two arguments. The first argument, handle, is a user-defined MATLAB structure and the second argument, str, is a text string. In the example above myprint prints the string to a file defined by handle.

The following code fragment shows how to tell MOSEK to send log output to the myprint function.

%
% In this example the MOSEK log info
% should be printed to the screen and to a file named
% mosek.log.
%

fid                = fopen('mosek.log','wt');
callback.log       = 'myprint';
callback.loghandle = fid;

%
% The 'handle' argument in myprint() will be identical to
% callback.loghandle when called.
%

mosekopt('minimize',prob,[],callback);

7.22.2. The iteration call-back function

It is possible to specify a function to be called frequently during the optimization. Typically this call-back function is used to display information about the optimization process or to terminate it.

The iteration call-back function has the following form:

function [r] =  myiter(handle,where,info)
% handle: Is a user-defined data structure.
% where : Is an integer indicating from where in the optimization
%         process the callback was invoked.
% info  : A MATLAB structure containing information about the state of the
%         optimization.

r = 0;   % r should always be assigned a value.

if handle.symbcon.MSK_CALLBACK_BEGIN_INTPNT==where
    fprintf('Interior point optimizer started\n');
end


% Print primal objective
fprintf('Interior-point primal obj.: %e\n',info.MSK_DINF_INTPNT_PRIMAL_OBJ);

% Terminate when cputime > handle.maxtime
if info.MSK_DINF_INTPNT_TIME > handle.maxtime
      r = 1;
else
  r = 0;
end
     

if handle.symbcon. MSK_CALLBACK_END_INTPNT==where
    fprintf('Interior-point optimizer terminated\n');
end

The function accepts three arguments: The first argument, handle, is a user-defined MATLAB structure, the second argument, where, indicates from where in the optimization process the call-back was invoked and the third argument, info, is a structure containing information about the process. For details about info see Section 8.1.7g If the function returns a non-zero value, MOSEK will terminate the optimization process immediately.

In order to inform MOSEK about the iteration call-back function the fields iter and iterhandle are initialized as shown in the following example.

[r,res]             = mosekopt('symbcon');
data.maxtime        = 100.0;
data.symbcon        = res.symbcon;

callback.iter       = 'myiter';
callback.iterhandle = data;

mosekopt('minimize',prob,[],callback);

7.23. The license system

By default a license token remains checked out for the duration of the matlab session. This can be changed such that the license is returned after each call to mosek by setting the parameter MSK_IPAR_CACHE_LICENSE.

param.MSK_IPAR_CACHE_LICENSE = 'MSK_OFF'; %set parameter.
[r,res] = mosekopt('minimize',prob,param); %call mosek.

By default an error will be returned if no license token is available. By setting the parameter MSK_IPAR_LICENSE_WAIT mosek can be instructed to wait until a license token is available.

param.MSK_IPAR_LICENSE_WAIT = 'MSK_ON'; %set parameter.
[r,res] = mosekopt('minimize',prob,param); %call mosek.
Wed Feb 29 16:16:54 2012